Why does julia systematically fails when doing \ operation?

General Usage
,
9

It appears that the convention adopted by Julia is that the backslash symbol used in x = A\b with a square A matrix really stands exclusively for the inverse of a matrix (computed perhaps using LU decomposition). If A is singular (or detected close to singular), the call of x= A\b fails. In such situations, SVD is your friend:

julia> G = svd(A)
SVD{Float64, Float64, Matrix{Float64}}
U factor:
3×3 Matrix{Float64}:
 -0.332289   0.850246   0.408248
 -0.549449   0.177311  -0.816497
 -0.766609  -0.495624   0.408248
singular values:
3-element Vector{Float64}:
 11.218599757513006
  0.3781791648532636
  1.9743360888694275e-16
Vt factor:
3×3 Matrix{Float64}:
 -0.332574  -0.479504  -0.812078
 -0.745695   0.660865  -0.08483
  0.57735    0.57735   -0.57735

julia> xₛ = G\b
3-element Vector{Float64}:
 0.9999999999999953
 2.000000000000006
 2.999999999999999

Obviously, when the backslash is used on the SVD object (a variable of type SVD), it does the standard routine of removing all those small enough singular values and the corresponding columns in U and V matrices and then doing the multiplication of b by the inverse of the reduced SVD of A

In fact, the same ideas are behind the concept of a pseudoinverse. And Julia can compute it too:

julia> x = pinv(A)*b
3-element Vector{Float64}:
 1.0
 2.0
 2.9999999999999996

It is just that the computation is a bit less efficient (you need to alocate memory for the pseudoinverse first, then compute it, and then multiply with it the b vector).

A bit cheaper (compared to the SVD-based solution) is a solution based on a QR decomposition:

julia> F = qr(A)
LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:
 -0.267261   0.872872   0.408248
 -0.534522   0.218218  -0.816497
 -0.801784  -0.436436   0.408248
R factor:
3×3 Matrix{Float64}:
 -3.74166  -5.34522   -9.08688
  0.0       0.654654   0.654654
  0.0       0.0       -1.11022e-16

julia> xᵣ = F\b
3-element Vector{Float64}:
  4.000000000000012
  4.99999999999999
 -0.0
3 Likes