More Accurate Bidiagonal Reduction for Computing the Singular Value Decomposition
Jesse L. Barlow · SIAM Journal on Matrix Analysis and Applications · 2002
Bidiagonal reduction is the preliminary stage for the fastest stable algorithms for computing the singular value decomposition (SVD) now available. However, the best-known error bounds on bidiagonal reduction methods on any matrix are of the form \[ A + \delta A = UBV^T , \] \[ \|\delta A\|_2 \leq \varepsilon_M f(m,n) \|A\|_2, \] where B is bidiagonal, U and V are orthogonal, $\varepsilon_M$ is machine precision, and f(m,n) is a modestly growing function of the dimensions of A. A preprocessing technique analyzed by Higham [Linear Algebra Appl., 309 (2000), pp. 153--174] uses orthogonal factorization with column pivoting to obtain the factorization \[ A=Q \left( \begin{array}{c} C^T \\ 0 \end{array} \right) P^T, \] where Q is orthogonal, C is lower triangular, and P is permutation matrix. Bidiagonal reduction is applied to the resulting matrix C. To do that reduction, a new Givens-based bidiagonalization algorithm is proposed that produces a bidiagonal matrix B that satisfies $C + \delta C = U (B + \delta B ) V^T$ where $\delta B$ is bounded componentwise and $\delta C$ satisfies a columnwise bound (based upon the growth of the lower right corner of C) with U and V orthogonal to nearly working precision. Once we have that reduction, there is a good menu of algorithms that obtain the singular values of the bidiagonal matrix B to relative accuracy, thus obtaining an SVD of C that can be much more accurate than that obtained from standard bidiagonal reduction procedures. The additional operations required over the standard bidiagonal reduction algorithm of Golub and Kahan [J. Soc. Indust. Appl. Math. Ser. B Numer. Anal., 2 (1965), pp. 205--224] are those for using Givens rotations instead of Householder transformations to compute the matrix V, and 2n 3/3 flops to compute column norms.