svds returns non-orthogonal singular vectors when there are degenerate singular values

[mhauru]

A colleague of mine spotted a bug with the following:

function bugdemo(D)
        M = diagm(D)
        U, S, V = svds(M, nsv=4)
        C = U' * M * V
        C[abs(C) .< 1e-12] = 0.
        println(diagm(S))  # These two should
        println(C)         # be the same.

        ide = U' * U
        ide[abs(ide) .< 1e-12] = 0
        println(ide)
end

D = rand(5)
D[5] = 0
bugdemo(D)

println()
D[2] = D[1]
bugdemo(D)

That prints out for instance

[0.9425868903202645 0.0 0.0 0.0
 0.0 0.7436506177105782 0.0 0.0
 0.0 0.0 0.3137218131235404 0.0
 0.0 0.0 0.0 0.1446005101113296]
[0.942586890320264 0.0 0.0 0.0
 0.0 0.7436506177105787 0.0 0.0
 0.0 0.0 0.3137218131235404 0.0
 0.0 0.0 0.0 0.14460051011132963]
[0.9999999999999993 0.0 0.0 0.0
 0.0 1.0000000000000009 0.0 0.0
 0.0 0.0 1.0 0.0
 0.0 0.0 0.0 0.9999999999999993]

[0.9425868903202645 0.0 0.0 0.0
 0.0 0.9425868903202641 0.0 0.0
 0.0 0.0 0.7436506177105784 0.0
 0.0 0.0 0.0 0.31372181312354036]
[0.942586890320264 -0.8015220638016495 0.0 0.0
 -0.8015220638016496 0.9425868903202647 0.0 0.0
 0.0 0.0 0.7436506177105782 0.0
 0.0 0.0 0.0 0.31372181312354047]
[0.9999999999999996 -0.8503428936183435 0.0 0.0
 -0.8503428936183435 1.0000000000000007 0.0 0.0
 0.0 0.0 0.9999999999999984 0.0
 0.0 0.0 0.0 1.0000000000000004]

So in the first case there are no degenerate singular values and everything goes fine: After doing an SVD of M, we can multiply M with U and S and get a diagonal matrix with the singular values on the diagonal. The columns of U are orthogonal.

However, if we set D[2] = D[1] so that there's degeneracy, the same thing fails. The singular values come out right, but the singular vectors don't: You can see that the block of C corresponding to the degenerate singular values is not diagonal and the columns of U are not orthogonal.

I'm guessing that this is related to the non-uniqueness of singular vectors corresponding to degenerate singular values, but I haven't looked at the code for svds. Note that if you run that code a few times, sometimes things work out correctly by luck even in the latter case.

[ViralBShah]

If the same thing happens in octave/python/matlab, then this is likely an upstream issue.

[jiahao]

The behavior described is not a bug, but an intrinsic limitation of the algorithm used by svds that is documented in all good textbooks on numerical linear algebra. Krylov subspace methods can only find at most one eigenvector (singular vector) of a multiple eigenvalue (singular value), since all the other eigenvectors for the multiple eigenvalue must be orthogonal to the one eigenvector found, and hence by construction can never be found in the Krylov subspace containing the one eigenvector found.

For example, Parlett writes in The Symmetric Eigenvalue Problem, Ch. 12:

The fact that svds still works occasionally is due to roundoff error and is also documented in textbooks.

The only workaround for sparse matrices is to either use a different algorithm, such as a sparse direct eigensolver, or a block Krylov method (c.f. Parlett 13.10, pp. 316 ff.), neither of which are provided by svds.

The only reasonable actions here are:

• close as invalid - not a bug, or
• change to documentation issue: document the limitation.

[andreasnoack]

@jiahao I'm not sure this is the complete story. In general, I find that eigs([0 A;A' 0]) finds orthogonal eigenvectors for matrices of the type mentioned in this issue. E.g.

julia> A = diagm(sort(rand(5), rev = true));

julia> A[2,2]=A[1,1];

julia> AA = [zeros(5,5) A; A' zeros(5,5)];

julia> VV = eigs(AA, nev = 8)[2]
10×8 Array{Float64,2}:
  0.396673     -0.301188     -0.639755      0.585363     -2.77556e-17  -2.77556e-17   0.0           2.25514e-17
 -0.585363     -0.639755      0.301188      0.396673      2.77556e-16   1.11022e-16   1.11022e-16   3.05311e-16
  9.71445e-17  -2.48711e-16   1.14518e-16   0.0          -0.707107      0.707107      7.91034e-16  -5.55112e-16
  0.0           5.81047e-16  -2.42155e-17  -5.55112e-17  -2.54198e-16   7.77156e-16  -0.707107      0.707107   
  2.77556e-17   0.0           0.0           0.0          -1.94289e-16  -5.55112e-17  -2.91434e-16  -1.66533e-16
 -0.396673     -0.301188     -0.639755     -0.585363      8.32667e-17   0.0           2.77556e-17   3.90313e-16
  0.585363     -0.639755      0.301188     -0.396673      3.88578e-16  -1.11022e-16   5.55112e-17   3.85109e-16
 -1.11022e-16  -4.77305e-16   1.2396e-16    4.16334e-17  -0.707107     -0.707107     -9.02056e-16  -3.05311e-16
  0.0           3.10657e-16   4.91347e-17   1.11022e-16  -5.57705e-16  -7.49401e-16   0.707107      0.707107   
  5.55112e-17  -3.33067e-16   0.0           5.55112e-17  -1.11022e-16  -2.22045e-16  -7.21645e-16  -3.33067e-16

Parlett is careful about mentioning "...exact arithmetic..." in the paragraphs after the first you list and, as I read him, he suggests that, in practice, rounding errors solve most of the problem.

To me, the problem seems to be that we'll have to be a little more careful about which of the columns of the eigenvectors of [0 A;A' 0] to extract, e.g. in this case, we should extract [1,4,5,7] and not 1:2:8=[1,3,5,7] as we do in

julia/base/linalg/arnoldi.jl

Line 382 in 15626c3

left_sv = sqrt(2) * ex[2][ 1:size(X,1), ind ] .* sign(ex[1][ind]')

[jiahao]

@andreasnoack When roundoff breaks the degeneracy in the subspace traversed and allows some component in other eigenvectors, what you suggest is reasonable. However, I'm skeptical that users can rely on floating point roundoff to make their computations work in general. I don't recall Parlett ever saying that "rounding errors solve most of the problem", rather only that he can't prove much in finite precision about what happens.

[andreasnoack]

The problem you suggest is not a problem for our eigs or MATLAB's svds so I've changed the title back to the original one.

[jiahao]

I've just checked MATLAB 2016a's implementation of svds and it uses the implicitly restarted Lanczos method of Baglama and Reichl (equivalent to IterativeSolvers.svdl but with Larsen's partial reorthogonalization strategy implemented), not ARPACK on an augmented matrix. So we are comparing different algorithms in MATLAB's svds and our own. (MATLAB's eigs still uses ARPACK)

(Apparently this is new in 2016a; @andreasnoack has seen that svds in MATLAB 2015b uses ARPACK on an augmented matrix.)