# Golub-Kahan-Lanczos Bidiagonalization Procedure implementation doesn't produce bidiagonal matrix

I'm trying to implement the aforementioned procedure using this website as a reference. At the end of the page the algorithm is described as follows:

I think I've mapped the given algorithm to code correctly. However when I use the U and V obtained in the procedure to compute $$U^{T}AV$$ I don't get a bidiagonal matrix. The values along the diagonals indeed are the ones computed within the procedure but the rest of the values are not zeros. Here's my code:

def golub_kahan(a):
n = a.shape[1]
v = np.ones(n, dtype="float32") / np.sqrt(n)
u = np.zeros(a.shape[0], dtype="float32")
beta = 0
U, V = np.zeros_like(a, dtype="float32"), np.zeros((n,n), dtype="float32")

for i in range(n):
V[:, i] = v
u = a @ v - beta * u
alpha = np.linalg.norm(u)
u /= alpha
U[:, i] = u
v = a.T @ u - alpha * v
beta = np.linalg.norm(v)
v /= beta

return U, V


I'm not sure where my mistake is, did I misunderstand the algorithm somehow?

• Your algorithm is correct, but due to rounding errors you may see small entries outside the bidiagonal, instead of exact zeros. These values are about 1e-16 times as large as the other values.
– wim
Dec 28, 2018 at 10:16

This algorithm suffers from similar numerical stability problems as the symmetric Lanczos tridiagonalization algorithm, see here

In exact arithmetic, after $$k$$ steps the following holds: $$U_k^T U_k =I$$, $$V_k^T V_k =I$$, and $$U_k^T A V_k = B_k$$, where $$B_k$$ is the bidiagonal matrix with diagonal elements $$\alpha_i$$ and superdiagonal elements $$\beta_i$$. However, in floating point arithmetic, even for modest $$k$$, the matrices $$U_k$$ and $$V_k$$ may become far from orthogonal. Moreover, $$U_k^* A V_k$$ may have relatively large entries not only on the bidiagonals, but also at positions far away from the diagonal. Nevertheless, the largest singular value of $$B_k$$ is often a quite good approximation of the largest singular value of $$A$$, even for relatively small $$k$$.

A much more stable bidiagonalization algorithm is the Householder bidiagonalization algorithm, which is described in Paragraph 5.4.8 of Matrix computations by Golub and Van Loan, 4th edition. This method uses Householder reflections to bidiagonalize matrix $$A$$, which is relatively expensive but also very stable.

It is possible to improve the Lanczos bidiagonalization algorithm, and to maintain the orthogonality, by reothogonalization See also: Some remarks on bidiagonalization and its implementation. Note that Lanczos bidiagonalization only needs matrix vector products with $$A$$ and $$A^T$$, while the Householder bidiagonalization algorithm needs explicit access to the matrix entries of $$A$$, which is a disadvantage

• Thanks for the answer! I've modified the algorithm to use double precision and now it seems to be good enough for my purposes. Now I want to use QR algorithm applied to a symmetric matrix the way it's described here but when I use my implementation of QR algorithm the results are all close to 0 (around 1e-14). When i use np.eigvals instead of my implementation the positive eigenvalues indeed match the singular values computed using np.svd on the original matrix. Dec 30, 2018 at 16:50
• Here's how I've implemented the algorithm pastebin.com/4dtdfurm . It seems that this approach is too naive, I was unable to find an alternative implementation that was understandable to me. Dec 30, 2018 at 16:58
• Probably you are mixing up the QR method and the QR factorization. The QR factorization (scipy.linalg.qr) is a factorization of the form $A=QR$, with $Q$ orthogonal and $R$ upper triangular. The $QR$ factorization does not compute the eigenvalues. The QR method (scipy.linalg.schur) computes a factorization of the form $A=U T U^H$, with a unitary $U$ (or orthogonal $U$ if $A$ is real symmetric), and $T$ upper triangular. The eigenvalues of $A$ are on the diagonal of $T$.
– wim
Dec 30, 2018 at 23:31