# 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 '18 at 10:16

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. – Chen Guevara Dec 30 '18 at 16:50
• 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 '18 at 23:31