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.

Sign up for GitHub

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

Jump to bottom

Fix bug in pinv #45009

Merged
merged 9 commits into from
May 2, 2022
Merged

Fix bug in pinv #45009

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
26f6d39
fix bug in LinearAlgebra.pinv
hyrodium Apr 17, 2022
c0a1562
update test_pinv function
hyrodium Apr 17, 2022
9106e4d
remove unused argument `tol3` in test_pinv
hyrodium Apr 17, 2022
d6021aa
update variable definitions in test_pinv
hyrodium Apr 17, 2022
d871977
add more tests in test_pinv
hyrodium Apr 17, 2022
54fedce
remove arguments (m,n) from test_pinv
hyrodium Apr 17, 2022
2087cd3
remove unused functions randn_float64 and randn_float32
hyrodium Apr 17, 2022
ec56620
add pinv tests for non-square diagonal matrices
hyrodium Apr 17, 2022
61611f1
remove for _ in 1:100
hyrodium Apr 18, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions stdlib/LinearAlgebra/src/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1449,12 +1449,13 @@ function pinv(A::AbstractMatrix{T}; atol::Real = 0.0, rtol::Real = (eps(real(flo
return similar(A, Tout, (n, m))
end
if isdiag(A)
ind = diagind(A)
dA = view(A, ind)
indA = diagind(A)
dA = view(A, indA)
maxabsA = maximum(abs, dA)
tol = max(rtol * maxabsA, atol)
B = fill!(similar(A, Tout, (n, m)), 0)
B[ind] .= (x -> abs(x) > tol ? pinv(x) : zero(x)).(dA)
indB = diagind(B)
B[indB] .= (x -> abs(x) > tol ? pinv(x) : zero(x)).(dA)
return B
end
SVD = svd(A)
Expand Down
69 changes: 33 additions & 36 deletions stdlib/LinearAlgebra/test/pinv.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,69 +63,54 @@ function tridiag(T::Type, m::Integer, n::Integer)
end
tridiag(m::Integer, n::Integer) = tridiag(Float64, m::Integer, n::Integer)

function randn_float64(m::Integer, n::Integer)
a=randn(m,n)
b = Matrix{Float64}(undef, m, n)
for i=1:n
for j=1:m
b[j,i]=convert(Float64,a[j,i])
end
end
return b
end

function randn_float32(m::Integer, n::Integer)
a=randn(m,n)
b = Matrix{Float32}(undef, m, n)
for i=1:n
for j=1:m
b[j,i]=convert(Float32,a[j,i])
end
end
return b
end

function test_pinv(a,tol1,tol2)
m,n = size(a)

function test_pinv(a,m,n,tol1,tol2,tol3)
apinv = @inferred pinv(a)

@test size(apinv) == (n,m)
@test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1
x0 = randn(n); b = a*x0; x = apinv*b
@test norm(a*x-b)/norm(b) ≈ 0 atol=tol1
apinv = pinv(a,sqrt(eps(real(one(eltype(a))))))
@test norm(apinv*a*apinv-apinv)/norm(apinv) ≈ 0 atol=tol1
for _ in 1:100
b = a*randn(n)
x = apinv*b
@test norm(a*x-b)/norm(b) ≈ 0 atol=tol1
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why the loop here? Is this intended or something left over from editing? Same below in lines 83-87.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since randn(n) is included, I thought it would be better to verify multiple times rather than just one test.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think just one test should suffice. If we do this everywhere, the tests will take way too long.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, I have removed for _ in 1:100.


apinv = @inferred pinv(a,sqrt(eps(real(one(eltype(a))))))
@test size(apinv) == (n,m)
@test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol2
x0 = randn(n); b = a*x0; x = apinv*b
@test norm(a*x-b)/norm(b) ≈ 0 atol=tol2
@test norm(apinv*a*apinv-apinv)/norm(apinv) ≈ 0 atol=tol2
for _ in 1:100
b = a*randn(n)
x = apinv*b
@test norm(a*x-b)/norm(b) ≈ 0 atol=tol2
end
end

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
@testset for (m, n) in [(1000, 100), (100, 100), (100, 1000)]
default_tol = (real(one(eltya))) * max(m,n) * 10
tol1 = 1e-2
tol2 = 1e-5
tol3 = 1e-5
if real(eltya) == Float32
tol1 = 1e0
tol2 = 1e-2
tol3 = 1e-2
end
@testset "dense/ill-conditioned matrix" begin
### a = randn_float64(m,n) * hilb(eltya,n)
a = hilb(eltya, m, n)
test_pinv(a, m, n, tol1, tol2, tol3)
test_pinv(a, tol1, tol2)
end
@testset "dense/diagonal matrix" begin
a = onediag(eltya, m, n)
test_pinv(a, m, n, default_tol, default_tol, default_tol)
test_pinv(a, default_tol, default_tol)
end
@testset "dense/tri-diagonal matrix" begin
a = tridiag(eltya, m, n)
test_pinv(a, m, n, default_tol, tol2, default_tol)
test_pinv(a, default_tol, tol2)
end
@testset "Diagonal matrix" begin
a = onediag_sparse(eltya, m)
test_pinv(a, m, m, default_tol, default_tol, default_tol)
test_pinv(a, default_tol, default_tol)
end
@testset "Vector" begin
a = rand(eltya, m)
Expand Down Expand Up @@ -164,6 +149,18 @@ end
@test C ≈ ones(2,2)
end

@testset "non-square diagonal matrices" begin
A = eltya[1 0 ; 0 1 ; 0 0]
B = pinv(A)
@test A*B*A ≈ A
@test B*A*B ≈ B

A = eltya[1 0 0 ; 0 1 0]
B = pinv(A)
@test A*B*A ≈ A
@test B*A*B ≈ B
end

if eltya <: LinearAlgebra.BlasReal
@testset "sub-normal numbers/vectors/matrices" begin
a = pinv(floatmin(eltya)/100)
Expand Down