# Bug in LinearAlgebra pinv()

A = [1.0 0.0;0.0 1.0;0.0 0.0]

A1 = pinv(A):
[1.0 0.0 1.0;0.0 0.0 0.0] wrong
A1*A
[1.0 0.0;0.0 0.0] wrong

A2 = inv(A’*A)A’:
[1.0 0.0 0.0;0.0 1.0 0.0] correct
A2
A
[1.0 0.0;0.0 1.0] correct

1 Like

Agreed, it looks like a bug:

``````julia> SVD = svd(A)
LinearAlgebra.SVD{Float64, Float64, Matrix{Float64}}
U factor:
3×2 Matrix{Float64}:
1.00000e+00  0.00000e+00
0.00000e+00  1.00000e+00
0.00000e+00  0.00000e+00
singular values:
2-element Vector{Float64}:
1.00000e+00
1.00000e+00
Vt factor:
2×2 Matrix{Float64}:
1.00000e+00  0.00000e+00
0.00000e+00  1.00000e+00

julia> (SVD.Vt * Diagonal(SVD.S)) * SVD.U' * A
2×2 Matrix{Float64}:
1.00000e+00  0.00000e+00
0.00000e+00  1.00000e+00
``````

OK, a hint: is `A` diagonal?

We should fix this bug, but it’s worth noting that you should pretty much never use `pinv`. Doing linear solves directly is pretty much always better.

Right, but I am actually computing a projection matrix which will be used many times.

then call `factorize(A)` and reuse that.

2 Likes

The internal part of `pinv` is just the SVD, as the above posts show.
So you might want to look at using `svd(A)` and using the SVD factors to construct your projection.

Thank you.