More accurate bidiagonal reduction for computing the singular value decomposition

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 + δA = UBV^T, ∥δA∥₂ ≤ ε M f(m,n)∥A∥₂ where B is bidiagonal, U and V are orthogonal, ε_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 (C^T 0) 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 + δC = U(B + δB)V^T where δB is bounded componentwise and δ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 flops to compute column norms.