LinearSolve.jl icon indicating copy to clipboard operation
LinearSolve.jl copied to clipboard

Factorization based solvers returning trivial result for least square solve

Open vpuri3 opened this issue 1 year ago • 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

vpuri3 avatar Aug 22 '24 14:08 vpuri3

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

vpuri3 avatar Aug 22 '24 18:08 vpuri3

QR you probably just need to turn on pivoting.

ChrisRackauckas avatar Aug 26 '24 11:08 ChrisRackauckas