julia
julia copied to clipboard
JuliaLang
`\` for square singular matrices
From this Discourse discussion.
using LinearAlgebra
julia> ones(3,2) \ ones(3)
2-element Vector{Float64}:
0.5
0.4999999999999999
julia> ones(3,3) \ ones(3)
ERROR: SingularException(2)
julia> ones(3,4) \ ones(3)
4-element Vector{Float64}:
0.2500000000000001
0.25
0.25
0.25
The issue here is that \ applies the LU decomposition to a square input but a pivoted QR to a rectangular input. The LU solver cannot handle singular inputs while QR can.
julia> qr(ones(3,3),ColumnNorm()) \ ones(3)
3-element Vector{Float64}:
0.33333333333333337
0.33333333333333337
0.33333333333333337
A naive fix might just be to check the lu(A) call within \ and re-try with qr(A,ColumnNorm()) if it fails via singularity.
I would much rather have an error thrown at me when the matrix is singular. If you want a least-squares solution you should, IMO, use qr explicitly.
I agree that a SingularException is often useful - and that sidestepping it could cause plenty of chaos. However, I don't think it's terribly consistent that this works for all systems EXCEPT those that happen to be square.
At the very least, this might be a good target for an error hint indicating that a manual qr call can provide a minimum-norm least-squares solution.
However, I don't think it's terribly consistent that this works for all systems EXCEPT those that happen to be square.
Agreed, but I would prefer \ to only apply to square matrices instead (but obviously that ship has sailed).
I don't think it is feasible to change the behavior of \ at this point. A more helpful error message is very possible, however.
However, we might need to add an additional field to SingularException to add more context, as this error can be thrown in other contexts than ldiv!/rdiv! (e.g. in inv or cond) where a suggestion to use qr would not be appropriate, similar to what we currently do with throw_complex_domainerror.
