Pseudo-inverse and the backslash operator

Specific Domains Optimization (Mathematical)
1

So, I have a sparse matrix A and a vector b and I want to compute pinv(A)*b.

After reading this post I tried the following:

    using LinearAlgebra

    julia>  A = [1 1 1; 2 2 2; 1 0 1]
    3×3 Array{Int64,2}:
     1  1  1
     2  2  2
     1  0  1
    
    julia> b = [1,0,0]
    3-element Array{Int64,1}:
     1
     0
     0
    
    julia> pinv(A)
    3×3 Array{Float64,2}:
     -2.22045e-16  2.77556e-17   0.5
      0.2          0.4          -1.0
     -1.52656e-16  1.66533e-16   0.5
    
    julia> pinv(A)*b
    3-element Array{Float64,1}:
     -2.220446049250313e-16
      0.20000000000000048
     -1.5265566588595902e-16

Everything is fine up to this point, but when I try the backslash operator i get:

    julia> qr(A)\b
    3-element Array{Float64,1}:
     -1.8385868488076965e15
     -1.8242825833464351e-16
      1.838586848807696e15

I’m not sure what happened here. I don’t understand what the result is supposed to be, but clearly it is not pinv(A)*b. Am I using this wrongly or is it a bug?

2

Using a pivoted QR factorization will give you the result you are looking for:

julia> qr(A, Val(true)) \ b
3-element Array{Float64,1}:
 -2.379544922143407e-16 
  0.20000000000000057   
 -2.1619865292617238e-16

The standard QR algorithm only gives the same result as a multiplication with the pseudo-inverse if the matrix has full column rank.

4 Likes
3

That works, thanks! However, I’m a bit perplexed by this. What does qr(A)\b actually do? What is the result?

4

I think it solve the linear system Rx=Q^Tb but assumes that R has full rank. In the algorithm for the pivoted QR is an explicit check on the rank of R (which is also easier to do there) and then computes a least square solution for the system.

3 Likes