search

SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
244 stars 52 forks source link

Factorization based solvers returning trivial result for least square solve #531

Open vpuri3 opened 3 weeks ago

vpuri3 commented 3 weeks ago 
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 commented 3 weeks ago

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

ChrisRackauckas commented 2 weeks ago

QR you probably just need to turn on pivoting.