Golub-Kahan bidiagonalization
/ GOH-loob KAH-han /
How does a computer actually compute the SVD A = U Sigma V^T? The naive idea — form A^T A, eigendecompose it to get V and the squared singular values, then recover U — is a known accuracy trap, because squaring A also squares its condition number and shreds the small singular values. Golub-Kahan bidiagonalization is the clever first stage of the right way: it reduces A to a simple bidiagonal shape using orthogonal transformations that work directly on A, never forming A^T A, so accuracy is preserved.
The reduction sandwiches A between Householder reflections applied alternately from the left and the right: A = U_B B V_B^T, where B is BIDIAGONAL (nonzero only on the main diagonal and the one diagonal just above it) and U_B, V_B are orthogonal. A left reflection zeros a column below the diagonal; a right reflection zeros a row to the right of the superdiagonal; alternating them leaves exactly the bidiagonal skeleton. This costs a fixed O(m n^2) for an m-by-n matrix and is done once. The crucial fact: B has exactly the same singular values as A (orthogonal transformations preserve singular values), so the hard problem has been reduced to finding the SVD of a tiny-bandwidth bidiagonal matrix.
The second stage then computes the SVD of the bidiagonal B using an iterative process that is, in disguise, the implicitly-shifted symmetric QR algorithm applied to B^T B WITHOUT ever forming B^T B (the Golub-Reinsch algorithm, or for high relative accuracy on the small singular values, the dqds algorithm). This two-stage structure — direct orthogonal bidiagonalization, then implicit symmetric-QR on the bidiagonal — is exactly how LAPACK's SVD routines work, delivering backward-stable singular values and vectors. It is the same architectural idea as the dense eigensolver (reduce to a simple band form once, then iterate cheaply), specialized to the SVD.
For a tall 1000x50 data matrix A, Golub-Kahan first reduces it to a 50x50 bidiagonal B (just 50 + 49 = 99 nonzero numbers) via alternating left/right Householder reflections, accumulating U_B and V_B. Then a few sweeps of implicit QR on B turn the superdiagonal entries to zero, leaving the 50 singular values on the diagonal — all without ever computing the 50x50 matrix A^T A, whose condition number would be the square of A's.
Stage 1: orthogonal reduction of A to a bidiagonal matrix with the same singular values. Stage 2: implicit symmetric-QR sweeps on the bidiagonal.
The whole point of bidiagonalization is to AVOID forming A^T A, whose condition number is the square of A's — forming it would lose roughly half your digits on the small singular values. The reflections are applied alternately from both sides (not just one), which is what produces the bidiagonal rather than triangular form.