Currently, svd defaults to using LinearAlgebra.LAPACK.gesdd! under the hood. I vaguely mentioned some time ago that gesdd! gives incorrect results in some cases which can be circumvented by explicitly calling LinearAlgebra.LAPACK.gesvd! instead. However, for large enough matrices, gesvd! is much slower than gesdd!, which uses a divide-and-conquer approach. Whether the inaccuracies (which can be huge) are due to a bug in the implementation or a shortcoming of this different method I don’t know.
In any case, it naturally raises the question which method we should default to. It seems to me that we should prefer accuracy over speed and default to gesvd! in svd. Of course, we could introduce a keyword flag divide or similar that allows one to switch to the faster method. Just getting wrong results without a warning or anything, as it is the case right now, seems bad.
Note that Octave, which was also defaulting to gesdd has recently made precisely this switch as discussed in this issue.
Sort of taking the example from this thread, here is a MWE that demonstrates the issue:
using LinearAlgebra
function compan(c::AbstractVector{T}) where T
N = length(c)
A = diagm(-1 => ones(T, N-1))
@inbounds for i in 1:N
A[1,i] = - c[i]/c[1]
end
A
end
function svdvaldiff(N)
A = compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) ));
svals_gesdd = LinearAlgebra.LAPACK.gesdd!('A', copy(A))[2];
svals_gesvd = LinearAlgebra.LAPACK.gesvd!('A', 'A', copy(A))[2];
return maximum(svals_gesdd - svals_gesvd)
end
function gesdd_gesvd_comparison(Nmax=40)
for N in 1:Nmax
println("N = $(N), max diff svdvals: ", svdvaldiff(N))
end
end
Which gives the following output:
julia> gesdd_gesvd_comparison()
N = 1, max diff svdvals: 0.0
N = 2, max diff svdvals: 0.0
N = 3, max diff svdvals: 0.0
N = 4, max diff svdvals: 0.0
N = 5, max diff svdvals: 0.0
N = 6, max diff svdvals: 0.0
N = 7, max diff svdvals: 0.0
N = 8, max diff svdvals: 0.0
N = 9, max diff svdvals: 0.0
N = 10, max diff svdvals: 0.0
N = 11, max diff svdvals: 0.0
N = 12, max diff svdvals: 0.0
N = 13, max diff svdvals: 0.0
N = 14, max diff svdvals: 0.0
N = 15, max diff svdvals: 0.0
N = 16, max diff svdvals: 0.0
N = 17, max diff svdvals: 0.0
N = 18, max diff svdvals: 0.0
N = 19, max diff svdvals: 0.0
N = 20, max diff svdvals: 0.0
N = 21, max diff svdvals: 0.0
N = 22, max diff svdvals: 0.0
N = 23, max diff svdvals: 0.0
N = 24, max diff svdvals: 0.0
N = 25, max diff svdvals: 749.7516745517077
N = 26, max diff svdvals: 166.767736950586
N = 27, max diff svdvals: 584.3097624800283
N = 28, max diff svdvals: 635.9955873150014
N = 29, max diff svdvals: 752.0372910883041
N = 30, max diff svdvals: 934.3056077045787
N = 31, max diff svdvals: 532.6321090336445
N = 32, max diff svdvals: 644.2861698543643
N = 33, max diff svdvals: 362.0563181932318
N = 34, max diff svdvals: 524.8519471908377
N = 35, max diff svdvals: 681.9577012165146
N = 36, max diff svdvals: 959.8179240905338
N = 37, max diff svdvals: 116.35110997021529
N = 38, max diff svdvals: 510.2091542800206
N = 39, max diff svdvals: 245.12683184534578
N = 40, max diff svdvals: 6.869915826281117
