I can't find a single solid explanation of how to implement this -- whitepapers too detailed/confusing. Closest I came to an answer was this:


<http://www.hep.ucl.ac.uk/~bino/libbpm/doc/pro/html/gsl__linalg_8c-source.html>


see:


[1] gsl\_linalg\_bidiag\_decomp


[2] gsl\_linalg\_SV\_decomp (which calls [1])


Which perform it, but really, really obfusticate the process underneath.



Bidiagonalization and SVD of matrix
**Golub-Kahan Bidiagonalisation**


In this process householder reflectors are applied alternatively on the left and then the right. The $i^{\text{th}}$ left reflector introduces zeros below the diagonal in the $i^{\text{th}}$ column. The $i^{\text{th}}$ right reflector introduces zeros to the right of the first super-diagonal in the $i^{\text{th}}$ row.


In software packages I suspect they use a mixture of this Golub-Kahan bidiagonalisation and a process called Lawson-Hanson-Chan (LHC) bidiagonalisation depending on the size of the matrix.


**Computing the SVD**


The first phase of computing the SVD is bidiagonalising the matrix. Then the SVD of the bidiagonal matrix is determined by a process very similar to the QR algorithm. This process is described in Golub and Kahn, ["Calculating the singular values and pseudo-inverse of a matrix"](http://www.amath.washington.edu/~narc/win08/papers/golub-kahan.pdf) (1960).
Since this paper there have been some alterations to provide better accuracy when the singular values are small, see Demmel and Kahan: "Accurate singular values of bidiagonal matrices" (1990)


Having been lectured by N. Trefethen on this very subject he briefly mentioned a divide-and-conquer type algorithm was now state-of-the-art though I don't know much about the details.


