Bug in LinearAlgebra pinv()

General Usage
1

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
A2A
[1.0 0.0;0.0 1.0] correct

1 Like
2

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
3

Can you spot it? https://github.com/JuliaLang/julia/blob/bf534986350a991e4a1b29126de0342ffd76205e/stdlib/LinearAlgebra/src/dense.jl#L1427

4

OK, a hint: is A diagonal?

5

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.

6

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

7

then call factorize(A) and reuse that.

8

https://github.com/JuliaLang/julia/issues/44234

2 Likes
9

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.

10

Thank you.