cond(A,1) fails to produce Inf for A exactly singular

[andreasvarga]

The following error occurs with Julia V1.1:

A = [1. 1.; 0. 0.]
cond(A,1)
ERROR: SingularException(2)
Stacktrace:
[1] checknonsingular at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\factorization.jl:12 [inlined]
[2] #lu!#103(::Bool, ::Function, ::Array{Float64,2}, ::Val{true}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lu.jl:41
[3] #lu#107 at .\none:0 [inlined]
[4] lu at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lu.jl:142 [inlined] (repeats 2 times)
[5] _cond1Inf(::Array{Float64,2}, ::Int64) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\dense.jl:1390
[6] cond(::Array{Float64,2}, ::Int64) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\dense.jl:1386
[7] top-level scope at none:0

The result should be Inf, as correctly computed by cond(A,2).

[andreasnoack]

Thanks for reporting this. We'd have to disable exception throwing from the factorization and return Inf if exactly singular. Would you be able to prepare a PR?

[andreasvarga]

I am afraid I will be not able to do a PR (actually I never did one). Just a comment: I am facing a similar case when estimating the reciprocal condition number of a singular linear operator (e.g., encountered when solving Sylvester and Lyapunov matrix equations, see https://github.com/andreasvarga/MatrixEquations.jl ). Since my solution is far from a perfect one, my hope is to learn how the exact singularity aspect can be handled in the Linear Algebra package of Julia. I hope there exists an elegant way which could serve as paradigm for similar problems (the Sylvester equation solver LAPACK.trsyl! also exits with an error in the singular case). I implemented a wrapper to the LAPACK routine DLACON to estimate 1-norm of an operator. My solution is to force DLACON to produce a large value for the estimated norm by providing large input vectors, when the solution process fails (I used a “try and catch” scheme for detecting singularity errors). The problem is that ANY error is caught in this way. A way to label errors (as done in MATLAB) would be a way for a more specific error control. Many thanks for your time, Andreas Varga Von: Andreas Noack [mailto:[email protected]] Gesendet: Dienstag, 17. September 2019 13:01 An: JuliaLang/julia Cc: andreasvarga; Author Betreff: Re: [JuliaLang/julia] cond(A,1) fails to produce Inf for A exactly singular (#33297) Thanks for reporting this. We'd have to disable exception throwing from the factorization and return Inf if exactly singular. Would you be able to prepare a PR? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#664?email_source=notifications&email_token=ALJDHEF3ZAVIPIQZTR3I4CDQKC2H3A5CNFSM4IXOC2XKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD64EUEI#issuecomment-532171281> , or mute the thread <https://github.com/notifications/unsubscribe-auth/ALJDHECCERXWAQBSEJSFWLLQKC2H3ANCNFSM4IXOC2XA> . <https://github.com/notifications/beacon/ALJDHEBLI5HRNO5HHYGFCFDQKC2H3A5CNFSM4IXOC2XKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD64EUEI.gif> [ { "@context": "http://schema.org", "@type": "EmailMessage", "potentialAction": { "@type": "ViewAction", "target": "#664?email_source=notifications\u0026email_token=ALJDHEF3ZAVIPIQZTR3I4CDQKC2H3A5CNFSM4IXOC2XKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD64EUEI#issuecomment-532171281", "url": "#664?email_source=notifications\u0026email_token=ALJDHEF3ZAVIPIQZTR3I4CDQKC2H3A5CNFSM4IXOC2XKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD64EUEI#issuecomment-532171281", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { "@type": "Organization", "name": "GitHub", "url": "https://github.com" } } ]

[antoine-levitt]

I am afraid I will be not able to do a PR (actually I never did one).

Your choice of course, but doing this kind of small fixes is a great way to contribute, and lowers the bar for future PR (to this project or others). https://www.digitalocean.com/community/tutorials/how-to-create-a-pull-request-on-github is a good tutorial on how to create a PR if you're already familiar with git; for this kind of small PR you can even do the whole thing online by just going into a file on the github web interface and clicking the edit button on top. The offending line here is
_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A), p, opnorm(A, p)) in stdlib/LinearAlgebra/dense.jl, which calls lu without checking it succeeds.

[andreasvarga]

Thanks for the hint. The fix for this bug is easy. Replace the offending line

_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A), p, opnorm(A, p))

in stdlib/LinearAlgebra/dense.jl, with

_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A;check=false), p, opnorm(A, p))

Not so easy is to fix the error in the next line

_cond1Inf(A::AbstractMatrix, p::Real) = opnorm(A, p)*opnorm(inv(A), p)

where inv(A) fails for singular A with

ERROR: LAPACKException(2).

This failure can be generated with

a = [1 1; 0 0]
cond(a,1)

[andreasvarga]

Digging a little bit more, the following errors came up:

a = [1 1;0 0]
cond(UpperTriangular(a),1)
ERROR: LAPACKException(2)
Stacktrace:
[1] chklapackerror at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lapack.jl:38 [inlined]
[2] trtrs!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\lapack.jl:3348
[3] ldiv! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\triangular.jl:583 [inlined]
[4] inv(::UpperTriangular{Int64,Array{Int64,2}}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\triangular.jl:630
[5] _cond1Inf(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\dense.jl:1391
[6] cond(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\dense.jl:1386
[7] top-level scope at none:0

but it works for

cond(UpperTriangular(float(a)),1)

Similar errors come up for LowerTriangular and Diagonal matrices. Apparently all these errors are related to the attempt to compute the inverses of singular matrices.

[antoine-levitt]

Two possibilities here: computing the inverse using an explicit lu factorization and check that that worked with something like

function _cond1Inf(A::AbstractMatrix, p::Real)
    fac = lu(A, check=true)
    issuccess(fac) || return convert(real(eltype(A)), Inf)
    opnorm(A, p)*opnorm(inv(fac), p)
end

or using a try/catch to catch the error in the inv. I think the explicit lu factorization is simpler here.

[andreasvarga]

The lu-based method seems to be the right way, because you have the full control in detecting exact singularity. The solution with try and catch can be unreliable, because there are many ways to produce errors or overlook them (for example, inv(rand(2,1)) or inv(0.)).

Many thanks for your help in addressing this issue.

[andreasvarga]

A small correction: I think the right code is:

function _cond1Inf(A::AbstractMatrix, p::Real)
    fac = lu(A, check=false)
    issuccess(fac) || return convert(real(eltype(A)), Inf)
    opnorm(A, p)*opnorm(inv(fac), p)
end

Still the warning "Failed factorization of type LU{Float64,Array{Float64,2}}" in the case of singularity is disturbing!

[goggle]

There is a problem with introducing lu to solve this issue: lu currently does not
work for all kind of matrices, e.g.

lu(Diagonal(ones(3,3)))

currently fails.
This means that

cond(Diagonal(ones(3,3)), 1)

would fail instead of resulting in

julia> cond(Diagonal(ones(3,3)), 1)
1.0

So what can we do here? Should we first make sure that lu accepts all kind of matrix types, or use another approach to solve this issue?

[andreasnoack]

What about a cond(Diagonal) method? Even disregarding the lu issue, it seems like a good idea.

[goggle]

It's not only the Diagonal type, but also Bidiagonal, SymTridiagonal, LowerTriangular, UpperTriangular, UnitLowerTriangular, UnitUpperTriangular and probably some more, for which the lu factorization fails.

[andreasnoack]

Maybe we can add a small helper function around the place where we call lu.