bidiagonal_decomposition Struct Template Reference

Golub-Kahan decomposition of a bidiagonal matrix B. More...

#include <svd.h>

List of all members.

Public Types

typedef FloatType scalar_t

Public Member Functions

 bidiagonal_decomposition (af::ref< scalar_t > const &diagonal, af::ref< scalar_t > const &off_diagonal, int kind, af::ref< scalar_t, af::mat_grid > const &u0, bool accumulate_u, af::ref< scalar_t, af::mat_grid > const &v0, bool accumulate_v, scalar_t epsilon=std::numeric_limits< scalar_t >::epsilon(), int max_iteration_multiplier=6)
 Prepare for the decomposition of the given bidiagonal matrix.
void compute ()
 Compute the decomposition.
void test_downward_iteration_convergence ()
void test_upward_iteration_convergence ()
void solve_2x2_case ()
void compute_trailing_wilkinson_shift ()
void compute_leading_wilkinson_shift ()
void do_implicit_shift_qr_iteration_downward (bool compute_shift=false)
 Do one Golub-Kahan iteration chasing nonzeroes from top-left to bottom-right.
void do_implicit_shift_qr_iteration_upward (bool compute_shift=false)
 Do one Golub-Kahan iteration chasing nonzeroes from bottom-right to top-left.
void do_implicit_zero_shift_qr_iteration_downward (bool compute_shift=false)
void do_implicit_zero_shift_qr_iteration_upward (bool compute_shift=false)
void sort ()
 Sort singular values in descending order.
std::size_t numerical_rank (scalar_t delta)
 Numerical rank.

Public Attributes

af::ref< scalar_t > d
af::ref< scalar_t > f
scalar_t eps
scalar_t tol
scalar_t thresh
af::ref< scalar_t, af::mat_grid > u
af::ref< scalar_t, af::mat_grid > v
givens::product< scalar_t > q_u
givens::product< scalar_t > q_v
int n_iterations
int n_max_iterations
bool has_iteration_converged
bool has_converged
bool sorted
int r
int s
int r0
int s0
scalar_t s_lower
scalar_t s_upper
bool shall_chase_nonzero_down
scalar_t shift

Detailed Description

template<typename FloatType>
struct scitbx::matrix::svd::bidiagonal_decomposition< FloatType >

Golub-Kahan decomposition of a bidiagonal matrix B.

See [1] for an introduction. Our code is based on [2,3,4].

References: [1] Golub and Van Loan, Matrix Computations, 3rd edition, section 8.6.2 [2] LAPACK source code (version 3.0) [3] Ake Bjorck, Numerical Methods for Least-Squares Problems [4] Accurate Singular Values of Bidiagonal Matrices, Demmel, J. and Kahan, W. SIAM J. Sci. Stat. Comput. 11, 873-912, 1990

Constructor & Destructor Documentation

bidiagonal_decomposition ( af::ref< scalar_t > const &  diagonal,
af::ref< scalar_t > const &  off_diagonal,
int  kind,
af::ref< scalar_t, af::mat_grid > const &  u0,
bool  accumulate_u,
af::ref< scalar_t, af::mat_grid > const &  v0,
bool  accumulate_v,
scalar_t  epsilon = std::numeric_limits<scalar_t>::epsilon(),
int  max_iteration_multiplier = 6 
) [inline]

Prepare for the decomposition of the given bidiagonal matrix.

The argument kind can be either upper_diagonal or lower_diagonal.

u0 and v0 shall be the U and V matrices from the bidiagonalisation. Each time a Givens rotation is performed, it may be accumulated into either u0 and v0 to eventually produce the final U and V of the SVD. Whether accumulation takes place depends on the arguments accumulate_u and accumulate_v.

References rotation::chase_nonzero_from_x1_to_y0(), and SCITBX_ASSERT.

Member Function Documentation

void compute (  )  [inline]

Compute the decomposition.

Reference: [2, DBDSQR] and [3, 2.6.4] A few important points:

  • we handle zeroes on the diagonal with Givens rotation theta = pi/2 as in [2, DBDSQR] (that alleviates the need of the decoupling stages described in [1]);
  • the Wilkinson shift is computed as a singular value of the trailing 2x2 block of B(r:s, r:s), as advised in [3, 2.6.4] and implemented in [2, DLAS2] (instead of an eigenvalue of the trailing 2x2 block of B(r:s, r:s)^T B(r:s, r:s) in the original Golub-Kahan algorithm described in [1]).
  • the shift of the leading 2x2 block of B(r:s, r:s)^T B(r:s, r:s) is carefully done to avoid overflow and underflow [2, DBDSQR, line 577]
  • when B(r:s, r:s) is a 2x2 matrix, we compute its SVD decomposition analytically with bidiagonal_2x2_decomposition (as in [2, DBDSQR] with DLASV2)
  • we check the direction of grading on B(r:s, r:s) diagonal and apply the algorithm from top to bottom or from bottom to top, depending on which end is the biggest ([2] and [3, 2.6.4]). This is essential for efficient convergence (or convergence at all when there are zero entries on the diagonal).

References bidiagonal_decomposition::do_implicit_shift_qr_iteration_downward(), bidiagonal_decomposition::do_implicit_shift_qr_iteration_upward(), bidiagonal_decomposition::do_implicit_zero_shift_qr_iteration_downward(), and bidiagonal_decomposition::do_implicit_zero_shift_qr_iteration_upward().

void do_implicit_zero_shift_qr_iteration_downward ( bool  compute_shift = false  )  [inline]

Do one Demmel-Kahan iteration, chasing nonzeroes from top-left to bottom-right

References rotation::apply_assuming_null_x0(), demmel_kahan_rotations::chase_non_zero_from_z_to_t(), rotation::chase_nonzero_from_x1_to_y0(), demmel_kahan_rotations::g1, demmel_kahan_rotations::g2, and rotation::zero_x1().

Referenced by bidiagonal_decomposition::compute().

void do_implicit_zero_shift_qr_iteration_upward ( bool  compute_shift = false  )  [inline]

Do one Demmel-Kahan iteration, chasing nonzeroes from bottom-right to top-left

References rotation::apply_assuming_null_x0(), demmel_kahan_rotations::chase_non_zero_from_z_to_t(), rotation::chase_nonzero_from_x1_to_y0(), demmel_kahan_rotations::g1, demmel_kahan_rotations::g2, and rotation::zero_x1().

Referenced by bidiagonal_decomposition::compute().

void sort (  )  [inline]

Sort singular values in descending order.

Implementation note: we use selection sort if U or V is accumulated.

Any sort algorithm needs to swap values which are out-of-order. Here it means swapping not only two singular values but also the associated singular vectors. Thus the cost of one swap is 2n ops, compared to the cost of one comparison which is 1 op. Since selection sort reaches the theoretical minimum number of swaps (n) while performing n(n+1)/2 comparisons, it is the algorithm of choice here.

Referenced by bidiagonal_decomposition::numerical_rank().

std::size_t numerical_rank ( scalar_t  delta  )  [inline]

Numerical rank.

A matrix A is said to have numerical delta-rank equal to k if

k = min { rank(B) | ||A-B||_2 <= delta }

After ordering the singular values in decreasing order, we can use the following result sigma[1] >= ... >= sigma[k] >= delta >= sigma[k+1] >= ... >= sigma[n]

References bidiagonal_decomposition::sort().

The documentation for this struct was generated from the following file:
Generated on Tue Sep 1 17:12:36 2009 for cctbx by  doxygen 1.5.6