From 9cf5df31f327b0e2c631f3b71cfe1cee0fbd3c3f Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 00:44:07 -0600 Subject: [PATCH 1/7] added rank method for SVD object alongwith unit tests for the method. issue #1126 --- src/svd.jl | 31 +++++++++++++++++++++++++++++++ test/svd.jl | 12 ++++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/svd.jl b/src/svd.jl index 7a88c4a6..b1e2976d 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -245,6 +245,37 @@ svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))] svdvals(x::Number) = abs(x) svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T} +""" + rank(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T} + +Compute the numerical rank of a `n × m` matrix by counting how many singular values are greater +than `max(atol, rtol*σ₁)` where `σ₁` is `A`'s largest calculated singular value. +`atol` and `rtol` are the absolute and relative tolerances, respectively. +The default relative tolerance is `n*ϵ`, +where `n` is the size of the smallest dimension of `A`, +and `ϵ` is the [`eps`](@ref) of the element type of `A`. +!!! note + Numerical rank can be a sensitive and imprecise characterization of + ill-conditioned matrices with singular values that are close to the threshold + tolerance `max(atol, rtol*σ₁)`. In such cases, slight perturbations to the + singular-value computation or to the matrix can change the result of `rank` + by pushing one or more singular values across the threshold. These variations + can even occur due to changes in floating-point errors between different Julia + versions, architectures, compilers, or operating systems. + +!!! compat "Julia 1.1" + The `atol` and `rtol` keyword arguments requires at least Julia 1.1. + In Julia 1.0 `rtol` is available as a positional argument, but this + will be deprecated in Julia 2.0. + +""" + +function rank(S::SVD{<:Any,T}; atol::Real = 0.0, rtol::Real = (min(size(S)...)*eps(real(float(one(eltype(S))))))) where {T} + svals = getfield(S, :S) + tol = max(atol, rtol*svals[1]) + count(>(tol), svals) +end + ### SVD least squares ### function ldiv!(A::SVD{T}, B::AbstractVecOrMat) where T m, n = size(A) diff --git a/test/svd.jl b/test/svd.jl index 9e8b5d5c..33c02517 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -294,4 +294,16 @@ end @test F.S ≈ F32.S end +@testset "rank svd" begin + # Test that the rank of an svd is computed correctly + @test rank(svd([1.0 0.0; 0.0 1.0])) == 2 + @test rank(svd([1.0 0.0; 0.0 0.9]), rtol=0.95) == 1 + @test rank(svd([1.0 0.0; 0.0 0.9]), atol=0.95) == 1 + @test rank(svd([1.0 0.0; 0.0 1.0]), rtol=1.01) == 0 + @test rank(svd([1.0 0.0; 0.0 1.0]), atol=1.01) == 0 + + @test rank(svd([1.0 2.0; 2.0 4.0])) == 1 + @test rank(svd([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0])) == 2 +end + end # module TestSVD From 77646012da5dde057cefe1163a8363b8beae4bc1 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 14:51:25 -0600 Subject: [PATCH 2/7] faster implementation of rank(::QRPivoted) --- src/qr.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/qr.jl b/src/qr.jl index 9a89e583..071555ec 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -538,8 +538,10 @@ end function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 - tol = max(atol, rtol*abs(A.R[1,1])) - return something(findfirst(i -> abs(A.R[i,i]) <= tol, 1:m), m+1) - 1 + rdiag = diag(getfield(A, :factors)) + tol = max(atol, rtol*abs(rdiag[1])) + + return something(findfirst(abs.(rdiag) .<= tol), m+1) - 1 end # Julia implementation similar to xgelsy From d7942fb82740adbfbf8aea7baa2b8c0ae156cf99 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 17:26:42 -0600 Subject: [PATCH 3/7] added stevengj's suggestion to qr. reversed changes in svd --- src/qr.jl | 6 +++--- src/svd.jl | 31 ------------------------------- test/svd.jl | 12 ------------ 3 files changed, 3 insertions(+), 46 deletions(-) diff --git a/src/qr.jl b/src/qr.jl index 071555ec..3b581c0d 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -538,10 +538,10 @@ end function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 - rdiag = diag(getfield(A, :factors)) - tol = max(atol, rtol*abs(rdiag[1])) + factors = getfield(A, :factors) + tol = max(atol, rtol*abs(factors[1,1])) - return something(findfirst(abs.(rdiag) .<= tol), m+1) - 1 + return something(findfirst(i -> abs(factors[i,i]) <= tol, 1:m), m+1) - 1 end # Julia implementation similar to xgelsy diff --git a/src/svd.jl b/src/svd.jl index b1e2976d..7a88c4a6 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -245,37 +245,6 @@ svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))] svdvals(x::Number) = abs(x) svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T} -""" - rank(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T} - -Compute the numerical rank of a `n × m` matrix by counting how many singular values are greater -than `max(atol, rtol*σ₁)` where `σ₁` is `A`'s largest calculated singular value. -`atol` and `rtol` are the absolute and relative tolerances, respectively. -The default relative tolerance is `n*ϵ`, -where `n` is the size of the smallest dimension of `A`, -and `ϵ` is the [`eps`](@ref) of the element type of `A`. -!!! note - Numerical rank can be a sensitive and imprecise characterization of - ill-conditioned matrices with singular values that are close to the threshold - tolerance `max(atol, rtol*σ₁)`. In such cases, slight perturbations to the - singular-value computation or to the matrix can change the result of `rank` - by pushing one or more singular values across the threshold. These variations - can even occur due to changes in floating-point errors between different Julia - versions, architectures, compilers, or operating systems. - -!!! compat "Julia 1.1" - The `atol` and `rtol` keyword arguments requires at least Julia 1.1. - In Julia 1.0 `rtol` is available as a positional argument, but this - will be deprecated in Julia 2.0. - -""" - -function rank(S::SVD{<:Any,T}; atol::Real = 0.0, rtol::Real = (min(size(S)...)*eps(real(float(one(eltype(S))))))) where {T} - svals = getfield(S, :S) - tol = max(atol, rtol*svals[1]) - count(>(tol), svals) -end - ### SVD least squares ### function ldiv!(A::SVD{T}, B::AbstractVecOrMat) where T m, n = size(A) diff --git a/test/svd.jl b/test/svd.jl index 33c02517..9e8b5d5c 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -294,16 +294,4 @@ end @test F.S ≈ F32.S end -@testset "rank svd" begin - # Test that the rank of an svd is computed correctly - @test rank(svd([1.0 0.0; 0.0 1.0])) == 2 - @test rank(svd([1.0 0.0; 0.0 0.9]), rtol=0.95) == 1 - @test rank(svd([1.0 0.0; 0.0 0.9]), atol=0.95) == 1 - @test rank(svd([1.0 0.0; 0.0 1.0]), rtol=1.01) == 0 - @test rank(svd([1.0 0.0; 0.0 1.0]), atol=1.01) == 0 - - @test rank(svd([1.0 2.0; 2.0 4.0])) == 1 - @test rank(svd([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0])) == 2 -end - end # module TestSVD From d270bda8cc66aadcd2180e9bacac969f522b9652 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 21:38:42 -0600 Subject: [PATCH 4/7] minor edit to qr(::QRPivoted) --- src/qr.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/qr.jl b/src/qr.jl index 3b581c0d..5def1eca 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -540,7 +540,6 @@ function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real( m == 0 && return 0 factors = getfield(A, :factors) tol = max(atol, rtol*abs(factors[1,1])) - return something(findfirst(i -> abs(factors[i,i]) <= tol, 1:m), m+1) - 1 end From c4aa5f3f4c2bdad39d31f236327e0164844e051b Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 30 Nov 2024 19:14:02 -0600 Subject: [PATCH 5/7] simplified rank(::QRPivoted) implementation --- src/qr.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/qr.jl b/src/qr.jl index 5def1eca..4be620bd 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -538,9 +538,8 @@ end function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 - factors = getfield(A, :factors) - tol = max(atol, rtol*abs(factors[1,1])) - return something(findfirst(i -> abs(factors[i,i]) <= tol, 1:m), m+1) - 1 + tol = max(atol, rtol*abs(A.factors[1,1])) + return something(findfirst(i -> abs(A.factors[i,i]) <= tol, 1:m), m+1) - 1 end # Julia implementation similar to xgelsy From d1c35e69de5ade948c3afd50cd279568b81a06ba Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 30 Nov 2024 19:17:11 -0600 Subject: [PATCH 6/7] simplified the rtol default eps call in rank(::QRPivoted) --- src/qr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qr.jl b/src/qr.jl index 4be620bd..eb23b638 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -535,7 +535,7 @@ function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T} return B end -function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(one(eltype(A.Q))))) * iszero(atol)) +function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(eltype(A.Q)))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 tol = max(atol, rtol*abs(A.factors[1,1])) From 3ee69e6034b673dcca8a79e8dd5458a7df546292 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 30 Nov 2024 19:18:59 -0600 Subject: [PATCH 7/7] using eltype(A) instead of eltype(A.Q) in the default val for rtol in rank(::QRPivoted) --- src/qr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qr.jl b/src/qr.jl index eb23b638..5020a2c8 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -535,7 +535,7 @@ function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T} return B end -function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(eltype(A.Q)))) * iszero(atol)) +function rank(A::QRPivoted; atol::Real=0, rtol::Real=min(size(A)...) * eps(real(float(eltype(A)))) * iszero(atol)) m = min(size(A)...) m == 0 && return 0 tol = max(atol, rtol*abs(A.factors[1,1]))