From 26f6d394926857e09a1301073b260849228e29d5 Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 19:38:52 +0900 Subject: [PATCH 1/9] fix bug in LinearAlgebra.pinv --- stdlib/LinearAlgebra/src/dense.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index e5909a459395c..d23dca5e6488e 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -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) From c0a1562ef3d626a85b41bb405bfd55c7beb69b28 Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 22:55:46 +0900 Subject: [PATCH 2/9] update test_pinv function --- stdlib/LinearAlgebra/test/pinv.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index d3eafb26797a9..2ef53b92baa42 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -88,12 +88,11 @@ end function test_pinv(a,m,n,tol1,tol2,tol3) apinv = @inferred pinv(a) - @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)))))) + apinv = @inferred pinv(a,sqrt(eps(real(one(eltype(a)))))) @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 From 9106e4d8205238b0f736f9c58750a4d75c43623f Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:16:44 +0900 Subject: [PATCH 3/9] remove unused argument `tol3` in test_pinv --- stdlib/LinearAlgebra/test/pinv.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 2ef53b92baa42..4ef47a9b5a9af 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -86,7 +86,7 @@ function randn_float32(m::Integer, n::Integer) end -function test_pinv(a,m,n,tol1,tol2,tol3) +function test_pinv(a,m,n,tol1,tol2) apinv = @inferred pinv(a) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 x0 = randn(n); b = a*x0; x = apinv*b @@ -103,28 +103,26 @@ end 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, m, n, 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, m, n, 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, m, n, 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, m, m, default_tol, default_tol) end @testset "Vector" begin a = rand(eltya, m) From d6021aa237fa13e7dc9930c4b031575cca3c378b Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:23:10 +0900 Subject: [PATCH 4/9] update variable definitions in test_pinv --- stdlib/LinearAlgebra/test/pinv.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 4ef47a9b5a9af..39d563b9d02d7 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -89,13 +89,19 @@ end function test_pinv(a,m,n,tol1,tol2) apinv = @inferred pinv(a) @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 + for _ in 1:100 + b = a*randn(n) + x = apinv*b + @test norm(a*x-b)/norm(b) ≈ 0 atol=tol1 + end apinv = @inferred pinv(a,sqrt(eps(real(one(eltype(a)))))) @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 + 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) From d871977919f8ba36accaa221b71a5c9baa881983 Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:25:17 +0900 Subject: [PATCH 5/9] add more tests in test_pinv --- stdlib/LinearAlgebra/test/pinv.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 39d563b9d02d7..8b107e45e23c5 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -89,6 +89,7 @@ end function test_pinv(a,m,n,tol1,tol2) apinv = @inferred pinv(a) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 + @test norm(apinv*a*apinv-apinv)/norm(apinv) ≈ 0 atol=tol1 for _ in 1:100 b = a*randn(n) x = apinv*b @@ -97,6 +98,7 @@ function test_pinv(a,m,n,tol1,tol2) apinv = @inferred pinv(a,sqrt(eps(real(one(eltype(a)))))) @test norm(a*apinv*a-a)/norm(a) ≈ 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 From 54fedceeb09330d4593f7203d3a7ab0798820c37 Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:32:18 +0900 Subject: [PATCH 6/9] remove arguments (m,n) from test_pinv --- stdlib/LinearAlgebra/test/pinv.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 8b107e45e23c5..fef702c38158a 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -86,8 +86,11 @@ function randn_float32(m::Integer, n::Integer) end -function test_pinv(a,m,n,tol1,tol2) +function test_pinv(a,tol1,tol2) + m,n = size(a) + apinv = @inferred pinv(a) + @test size(apinv) == (n,m) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 @test norm(apinv*a*apinv-apinv)/norm(apinv) ≈ 0 atol=tol1 for _ in 1:100 @@ -97,6 +100,7 @@ function test_pinv(a,m,n,tol1,tol2) end 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 @test norm(apinv*a*apinv-apinv)/norm(apinv) ≈ 0 atol=tol2 for _ in 1:100 @@ -118,19 +122,19 @@ 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) + 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) + 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) + 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) + test_pinv(a, default_tol, default_tol) end @testset "Vector" begin a = rand(eltya, m) From 2087cd3003ac0095718c3c9390dc31134b864dff Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:35:47 +0900 Subject: [PATCH 7/9] remove unused functions randn_float64 and randn_float32 --- stdlib/LinearAlgebra/test/pinv.jl | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index fef702c38158a..c15d380fbda11 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -63,29 +63,6 @@ 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) @@ -120,7 +97,6 @@ end tol2 = 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, tol1, tol2) end From ec56620c8e22efd548a86dd9cccb73ad52cb64cf Mon Sep 17 00:00:00 2001 From: hyrodium Date: Sun, 17 Apr 2022 23:43:52 +0900 Subject: [PATCH 8/9] add pinv tests for non-square diagonal matrices --- stdlib/LinearAlgebra/test/pinv.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index c15d380fbda11..8975e4014a9e6 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -149,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) From 61611f11dc59062853630c23711987cd37735fe6 Mon Sep 17 00:00:00 2001 From: hyrodium Date: Mon, 18 Apr 2022 21:59:56 +0900 Subject: [PATCH 9/9] remove for _ in 1:100 --- stdlib/LinearAlgebra/test/pinv.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 8975e4014a9e6..c7268865a0505 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -70,21 +70,17 @@ function test_pinv(a,tol1,tol2) @test size(apinv) == (n,m) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 @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 + b = a*randn(n) + x = apinv*b + @test norm(a*x-b)/norm(b) ≈ 0 atol=tol1 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 @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 + b = a*randn(n) + x = apinv*b + @test norm(a*x-b)/norm(b) ≈ 0 atol=tol2 end @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)