Return Failure ReturnCode if Linear Solve failed

[avik-pal]

This came up in discourse, our LS solvers are not robust and throw LAPACK exception. Changing it to QRPivoted, but this will probably slow down some end-user code by around 3x (if I remember that number correctly).

An alternative would be to introduce a special version of QR which does the following:

1. Is wide rectangular system --> always use QRPivoted
2. For long rectangular system, use QR by default, but if that fails (the linear solve not QR) refactorize using pivoted. This might seem expensive but the assumption is that the matrix (jacobian/hessian) is singular only at limited places.

[codecov]

Codecov Report

Attention: Patch coverage is 81.03448% with 11 lines in your changes are missing coverage. Please review.

Project coverage is 62.71%. Comparing base (c3f4c4f) to head (9e49a76).

Files	Patch %	Lines
src/LinearSolve.jl	81.03%	11 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #492      +/-   ##
==========================================
+ Coverage   62.55%   62.71%   +0.16%     
==========================================
  Files          29       29              
  Lines        2214     2221       +7     
==========================================
+ Hits         1385     1393       +8     
+ Misses        829      828       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ChrisRackauckas]

@oscardssmith @YingboMa thoughts?

[oscardssmith]

I think this is probably a step in the right direction since without pivoting, QR will often just fail, but we probably want to go further eventually. Specifically, we should definitely try out something like https://ieeexplore.ieee.org/document/8287754 and for more rectangular matrices, we should possibly be using one of the super magic randomized algorithms that use matrix sketching.

[YingboMa]

I am not sure about defaulting to pivoted QR. Pivoted QR and QR solve different classes of problems. QR is fine as long it's a full-rank least squares problem. SVD and pivoted QR can solve rank-deficient systems where the solution are not unique. Mathematically they are different classes of problems. Just like we won't use pivoted QR or SVD as defaults for square matrices, I am not sure we should use pivoted QR as a default for least squares problems. Should we add a hint on the rank of least squares problems?

That said, I think MATLAB defaults to QR with column pivoting though.

[ChrisRackauckas]

The problem with ideas of trying non-pivoting first and the doing pivoting is that if we do this with mutation then we're destroying our original matrix. I haven't looked into it but I presume that would be a non-recoverable issue here, so we'd have to make a copy?

[oscardssmith]

QR is fine as long it's a full-rank

The issue with this is that determining whether or not you have a full rank system is not easy to do (without doing a factorization in the first place)

[avik-pal]

Just like we won't use pivoted QR or SVD as defaults for square matrices, I am not sure we should use pivoted QR as a default for least squares problems.

I agree here, the problem which came up in discourse is not LS specific, if the jacobian was constructed to be not full rank even regular root finding would fail, that doesn't mean we default to PivotedQR even in the square case.

We can probably just document what to do when such failures happen.

[YingboMa]

It's common for the outer level algorithm to handle factorization failures. That's what people do for most square systems already. Not reporting singularity and just solve a different problem silently isn't a good behavior.

[avik-pal]

The issue here is it throws an exception instead of giving a return code (I think I do have an internal linear solve failure for some of the cases I detect) unless there is a way to avoid that in linear algebra.

[YingboMa]

Do you have an example that throws an exception? Fixing that is probably the right solution but not changing the default.

[YingboMa]

julia> sol = solve!(cache)
ERROR: LAPACKException(2)
Stacktrace:
  [1] chklapackerror
    @ LinearAlgebra.LAPACK ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:40 [inlined]
  [2] trtrs!(uplo::Char, trans::Char, diag::Char, A::SubArray{…}, B::SubArray{…})
    @ LinearAlgebra.LAPACK ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:3564
  [3] generic_trimatdiv!
    @ LinearAlgebra ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:830 [inlined]
  [4] _ldiv!
    @ LinearAlgebra ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:752 [inlined]
  [5] ldiv!
    @ LinearAlgebra ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:757 [inlined]
  [6] ldiv!(A::LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}, b::Vector{Float64})
    @ LinearAlgebra ~/.bin/julia-1.10.0/share/julia/stdlib/v1.10/LinearAlgebra/src/qr.jl:520

The [5]-call points to the ldiv! with a triangular matrix, so you can just check the factorization by

julia> F = qr([2.0 1.0
                         0.0 0.0
                         0.0 0.0]);

julia> any(iszero, F.factors[diagind(F.factors)])
true

Note that for a serious implementation, we should check from the end of the diagonal where zeros are more likely to appear.

[YingboMa]

The failure isn't from QR factorization, and it's just from ldiv!. We can use this fake one liner to check if a QR factorization is singular.

julia> using LinearAlgebra

julia> issingular_qr(F::LinearAlgebra.QRCompactWY) = ((m, n) = size(F); U = view(F.factors, 1:min(m,n), 1:n); any(iszero, Iterators.reverse(@view U[diagind(U)])))
issingular_qr (generic function with 1 method)

julia> A = zeros(4, 4);

julia> issingular_qr(qr(A))
true

[avik-pal]

Okay based on this discussion, the current state is we give a failure code instead of erroring.

[avik-pal]

I have no clue why the diff is so weird

[ChrisRackauckas]

And just make the NonlinearSolve default alg try with a different linear solve if the first fails? That's fine.