Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

pinv bug: matrix misdiagnosed as diagonal #911

Closed
PetrKryslUCSD opened this issue Feb 18, 2022 · 3 comments · Fixed by JuliaLang/julia#45009
Closed

pinv bug: matrix misdiagnosed as diagonal #911

PetrKryslUCSD opened this issue Feb 18, 2022 · 3 comments · Fixed by JuliaLang/julia#45009
Labels
bug Something isn't working

Comments

@PetrKryslUCSD
Copy link

PetrKryslUCSD commented Feb 18, 2022

The code

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

produces

[1.0 0.0 1.0;0.0 0.0 0.0]

An incorrect answer, which is due to the identification of A as diagonal on line
https://github.com/JuliaLang/julia/blob/bf534986350a991e4a1b29126de0342ffd76205e/stdlib/LinearAlgebra/src/dense.jl#L1433

I think the solution could be to check that the matrix is square and isdiag returns true.

Or, perhaps isdiag does not work correctly, and the solution is to fix that function to identify matrix as diagonal only if it is square?

The svd code below the if block would have worked correctly.
https://discourse.julialang.org/t/bug-in-linearalgebra-pinv/76663/2

@PetrKryslUCSD PetrKryslUCSD added the bug Something isn't working label Feb 18, 2022
@martinholters
Copy link
Member

martinholters commented Feb 18, 2022

I think the problem is the indexing operation in
https://github.com/JuliaLang/julia/blob/bf534986350a991e4a1b29126de0342ffd76205e/stdlib/LinearAlgebra/src/dense.jl#L1439
Using PermutedDimsArray(B, (2,1))[ind] on the LHS might do the trick.
EDIT: On second thought, B[diagind(B)] should work just as well and is a bit simpler.

@PetrKryslUCSD
Copy link
Author

That might work.

@hyrodium
Copy link
Contributor

Note that this bug doesn't exist on Julia v1.6.6.

v1.6.6

julia> A = [1.0 0.0;0.0 1.0;0.0 0.0]
3×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  0.0

julia> A1 = pinv(A)
2×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0

v1.7.2

julia> A = [1.0 0.0;0.0 1.0;0.0 0.0]
3×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  0.0

julia> A1 = pinv(A)
2×3 Matrix{Float64}:
 1.0  0.0  1.0
 0.0  0.0  0.0

Maybe JuliaLang/julia#39756 is related?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants