Generalized svd throws

[andreasnoack]

Refiling JuliaLinearAlgebra/LAPACK.jl#6 here

julia> A = [-0.33872753963694624 1.124096715384297 -0.6293570718176809; 0.03919190688122216 -0.1300617417823436 0.07281871376668783];

julia> B = [-1.5303758632785613 5.136068273894432 -2.9372584484394606; 0.5364872797265587 -2.4543618264129545 2.0986693466314685];

julia> U,V,Q,C,S,R = svd(A,B)
ERROR: LAPACKException(1)
Stacktrace:
 [1] chklapackerror at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lapack.jl:38 [inlined]
 [2] ggsvd3!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lapack.jl:1868
 [3] svd!(::Array{Float64,2}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/svd.jl:277
 [4] svd(::Array{Float64,2}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/svd.jl:281
 [5] top-level scope at none:0

I think this might be an issue with the LAPACK routine and not our wrappers. The error means

the Jacobi-type procedure failed to converge. For further details, see subroutine DTGSJA.

and the error is DTGSJA means

the procedure does not converge after MAXIT cycles.

cc: @alanedelman

[alanedelman]

Here is code that produces LAPACK bugs, also has a curiosity that maybe someone would like
to investigate which is why C is 2x2 82% of the time , and 2x3 18% of the time?
subtle rank decisions perhaps?

using LinearAlgebra
#print all the digits
Base.show(io::IO, x::Union{Float64,Float32}) = Base.Grisu._show(io, x, Base.Grisu.SHORTEST, 0, true, false)
m1,m2,n,r,rA,rB =2,2,3,2,1,2
v=[0,0,0]
trials = 200_000
for i = 1:trials
   
 E = randn(r,n)
 A = randn(m1,rA)*randn(rA,r)*E
 B = randn(m2,rB)*randn(rB,r)*E
 try
       U,V,Q,C,S,R =svd(A,B)
       v[size(C,2)] +=1
   catch
       println("Found an LAPACK Bug on trial #$i")
       show(A);show(B)
       println(v*100 ./trials)
   end
  
end
println("% of time, C is 2x2: ", v[2]*100 ./ trials)
println("% of time, C is 2x3: ", v[3]*100 ./ trials);

Sample run:

Found an LAPACK Bug on trial #9127
[-0.04648273093474292 0.07525280330253786 0.08835704902460387; 0.16474011316682655 -0.2667045391456321 -0.31314748429551553][1.146724235077151 -1.8763720982311796 -2.1938877941189032; 1.0614108598749987 -1.5736190848804763 -1.9147923724497076][0.0, 3.762, 0.801]
Found an LAPACK Bug on trial #185514
[0.04922786153958221 -0.04475784056737846 -0.0271191928324028; -1.5459064845767538 1.4055340574387678 0.8516261877918397][-1.1647247355858317 1.045582474734974 0.5750783254692098; -4.986924429583816 4.571096370718112 2.9312676204742663][0.0, 76.494, 16.262]
% of time, C is 2x2: 82.49
% of time, C is 2x3: 17.509

[ianzwaan]

Both of you are almost certainly correct. What I suspect happens is the following.

According to the documentation of xGGSVD3 size(C, 2) should equal rank([A; B]), which would be min(m1+m2,r) == 2 in exact arithmetic. However, in the preprocessing step (xGGSVP3) B is first changed to B = V * [0 B12; 0 0] * Z' using QR with pivoting, a rank decision, and an RQ decomposition. Then A is updated and partitioned such that [A11 A12] = A * Z', and the whole QR+pivoting, rank decision, and RQ thing is repeated for A11. This roughly results in rank([A; B]) being computed as rank(A11) + rank(B12), which may be 3 because QR with pivoting isn't as good for rank decisions as the SVD and because rank decisions are ill-posed anyway.

The problem with xTGSJA also has to do with floating-point woes, but is unrelated otherwise. For context, xTGSJA works on n-by-n upper-triangular matrices A and B, where B is additionally nonsingular, and iterates over a bunch of 2x2 subproblems. The routine stops iterating when the rows of A and B are approximately parallel (so that A * inv(B) is approximately diagonal). The tolerance used here is quite small, something like 2 * eps() * min(norm(A), norm(B)) for the 2x2 problem from this issue. This is a problem because xTGSJA uses xLAGS2 to solve the 2x2 subproblems, but the rows of the output of xLAGS2 aren't guaranteed to be THAT close to being parallel; see Bai and Demmel '93 for details.

TL;DR: these are LAPACK issues, not Julia issues.

[ViralBShah]

Should we file this upstream?

[ianzwaan]

I'm not sure. Rank decisions are ill-posed problems and can be notoriously difficult, especially, in my experience, when rank decisions of submatrices are involved that depend on other rank decisions. The second problem could be addressed with a more relaxed tolerance; however, then you'd have to decide what convergence means in the face of round-off errors and if it is acceptable to stop iterating when the columns of A and B are not within machine precision of being parallel.

It could be interesting to talk to the LAPACK people about this, but at the same time I think it's important to realize that methods like this can fail.