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
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
Can you spot it? https://github.com/JuliaLang/julia/blob/bf534986350a991e4a1b29126de0342ffd76205e/stdlib/LinearAlgebra/src/dense.jl#L1427
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.
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.