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

Question

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?

Answer 1

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