3
$\begingroup$

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: enter image description here

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?

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – wim Dec 28 '18 at 10:16
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Chen Guevara Dec 30 '18 at 16:50
  • $\begingroup$ 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. $\endgroup$ – Chen Guevara Dec 30 '18 at 16:58
  • $\begingroup$ 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$. $\endgroup$ – wim Dec 30 '18 at 23:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.