SciML/LinearSolve.jl

Factorization based solvers returning trivial result for least square solve

Opened this issue · 2 comments

using LinearAlgebra, LinearSolve

A = rand(10, 4) # tall and skinny system
A[:, 1] .= 0    # A has a column with all zeros
f = rand(10)    # RHS

prob1 = LinearProblem(A, f)           # tall and skinny system
prob2 = LinearProblem(A' * A, A' * f) # square singular system (has a row and a column of all zeros)

sol_df = solve(prob1)                    # [0.0, 0.0, 0.0, 0.0]
sol_lu = solve(prob1, LUFactorization()) # [0.0, 0.0, 0.0, 0.0]
sol_qr = solve(prob1, QRFactorization()) # [0.0, 0.0, 0.0, 0.0]
sol_gm = solve(prob2, SimpleGMRES())     # [0.0, -0.4916139636428349, 0.44349181143744504, 0.7631772710700259]
A \ f                                    # [0.0, -0.4916139636428351, 0.44349181143744565, 0.7631772710700254]

sol_df.retcode # ReturnCode.Failure = 9
sol_lu.retcode # ReturnCode.Failure = 9
sol_qr.retcode # ReturnCode.Failure = 9
sol_gm.retcode # ReturnCode.Success = 1

Come to think about it, I shouldn't be surprised that factorizations fail for rank-deficient/singular matrices. Better to use Krylov methods and get a non-unique solution

QR you probably just need to turn on pivoting.