LAPACK exception with gels!

Specific Domains Numerics
,
1

If I run the following program then Julia quits with an exception. I suspect this has something to do with rounding in the factorization. The matrix B below can be replaced with a matrix of zeros of the same size and the problem persists, but if A is replaced with a matrix of zeros then we get a solution. What is the best fix and what is the easiest fix for problems like this? Thanks!

To be clear: I understand that there is a discontinuity at zeros versus close to zero. But even if I set the remaining elements to exact zeros, I’m not getting something like a generalized inverse of A times B.

using LinearAlgebra
import LinearAlgebra.LAPACK.gels!

A = [0.06044523175910898 0.05417624043156454 -2.6133655153109947e-9 -2.6133655153109947e-9; 0.05417624043156454 0.1355075110293397 -6.536641412174804e-9 -6.536641412174804e-9; -2.6133655153109947e-9 -6.5366414139095275e-9 -4.0441353880288355e-8 -4.0441353880288355e-8; -2.6133655153109947e-9 -6.5366414139095275e-9 -4.0441353880288355e-8 -4.0441353880288355e-8]

B = [-0.08117075255109502 -0.04961506097021058 0.018683685925334124 0.018683688132780544 0.018683675632683132 0.01868368133744242 0.01868368337844292 0.018683710556390808 0.01868369486710419; -6.93889390390725e-18 -0.06474891160657195 0.009249843663413143 0.00924984475626679 0.00924983856776912 0.009249841392058273 0.009249842402508547 0.00924985585765707 0.009249848090268504; 1.1153959005309217e-25 5.84905020402288e-25 8.170386227774643e-9 8.170387188473825e-9 8.170381712421108e-9 8.170384213190477e-9 8.17038512532588e-9 8.170396991872857e-9 8.170390125098772e-9; 1.1479437019748901e-41 0.0 8.271806125530277e-25 8.271806125530277e-25 8.271806125530277e-25 8.271806125530277e-25 0.0 8.271806125530277e-25 8.271806125530277e-25]

gels!('N',A,B)
2

The matrix A is numerically singular

julia> cond(A)
2.2306564779458504e38

so you probably want to use a method that determines the rank as part of the solve and gives a solution to the rank deficient problem. In LAPACK that is xgelsy whereas xgels assumes that A has full rank. You can also get this with Julia’s API, i.e.

julia> ldiv!(qr!(copy(A), Val(true)), copy(B), eps())[1] ≈ gelsy!(copy(A), copy(B))[1]
true

where the eps() is the tolerance used in the determination of the rank.

1 Like