svd stuck in infinite loop for certain matrices

[cvsvensson]

Using svd on certain matrices leads to an infinite loop:

using LinearAlgebra, GenericLinearAlgebra
m = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
svd(m) # ok
svd(BigFloat.(m)) # stuck in infinite loop
svdvals!(Float64.(m), debug=true) # also stuck, so it's not BigFloat's fault.

The svd of some permutations of the matrix works:

perm = [2,3,1,4]
svd(BigFloat.(m[perm,perm])) # ok

For this particular matrix, any small perturbation (e.g. m + 1e-10*rand(4,4)) seem to work fine but there are other matrices, such as

m2 = [0.4 0 0 0 0 0 0 0; 0 0.8 0 0 0 0 1 0; 0 0 -0.3 0 0 0 0 0; 0 0 0 -0.4 0 0 0 0; 0 0 0 0 0.4 0 0 0; 0 0 0 0 0 -0.3 0 0; 0 1 0 0 0 0 0.8 0; 0 0 0 0 0 0 0 -0.4]

where svdvals!(m2+1e-2*rand(8,8), debug=true) often gets stuck.

The problem may be related to the code linked below. If I replace shift^2 < eps(B.dv[n1]^2) with true, the infinite loop is avoided, but I don't know if this is a proper thing to do.

GenericLinearAlgebra.jl/src/svd.jl

Lines 181 to 188 in 9f9b46f

	shift = svdvals2x2(d[n2 - 1], d[n2], e[n2 - 1])[1]
	if shift^2 < eps(B.dv[n1]^2)
	# Shift is too small to affect the iteration so just use
	# zero-shift algorithm anyway
	svdDemmelKahan!(B, n1, n2, U, Vᴴ)
	else
	svdIter!(B, n1, n2, shift, U, Vᴴ)
	end

[andreasnoack]

Thanks for reporting this. I can reproduce the issue. I'll post an update once I've figured out what is going on.

[andreasnoack]

#82 seems to fix the issue for the reported matrices here. Please let me know if you see any convergence problems for other matrices once the patch is released.

[cvsvensson]

Thanks for the fix!

However, even with the patch, here's a another matrix with problems:

using GenericLinearAlgebra
m = [
0.3   0.0   0.0  0.0  0.0  0.2  0.3   0.0;
0.0   0.0   0.0  0.0  0.1  0.0  0.0   0.0;
0.0  -0.2   0.0  0.0  0.0  0.0  0.0  -0.2;
0.3   0.0   0.0  0.0  0.0  0.2  0.4   0.0;
0.0   0.4  -0.2  0.0  0.0  0.0  0.0   0.3;
0.2   0.0   0.0  0.0  0.0  0.0  0.2   0.0;
0.0   0.0   0.0  0.1  0.0  0.0  0.0   0.0;
0.0   0.3  -0.2  0.0  0.0  0.0  0.0   0.3]
svdvals!(copy(m), debug=true) #stuck in infinite loop

Curiously, the convergence properties can be changed with a simple multiplication:

svdvals!(10*copy(m), debug=true) #converges!

[andreasnoack]

This is indeed still an issue. I believe we need a dedicated solver for the 2x2 case.

[shinaoka]

Just for your information, I've faced a similar problem when trying to decompose a (ill-conditioned) Hilbert matrix using Double64. svd gets stuck in an infinite loop if the matrix size is larger than 16.

using DoubleFloats
using GenericLinearAlgebra

# Hilbert matrix
# N <= 16 is OK
N = 17
A = zeros(Double64, N, N)
for j in 1:N, i in 1:N
    A[i, j] = 1/Double64(i+j-1)
end

svd = GenericLinearAlgebra.svd(A)

[rabbott999]

I ran into a similar issue to the one described here and I believe I've found the issue. It looks like there's a typo at

GenericLinearAlgebra.jl/src/svd.jl

Line 117 in f3e9d27

μ = abs(dv[j + 1])*(μ*(μ + abs(ev[j])))

the line

μ = abs(dv[j + 1])*(μ*(μ + abs(ev[j])))

should be

μ = abs(dv[j + 1])*(μ/(μ + abs(ev[j])))

This is what matches eq. 2.4 in http://www.netlib.org/lapack/lawnspdf/lawn03.pdf, and also the current version doesn't pass the check that scaling the matrix should scale μ by the same amount. If the bidiagonal matrix has small off-diagonal entries then μ can become very small (for the first matrix above matrix above I'm getting ~1e-41), which results in an impossible-to-meet stopping condition. Making the above change fixed all of the infinite loops I've found, including both of the matrices above.

[shinaoka]

I also have confirmed that this fixes the issue!