From 9cf5df31f327b0e2c631f3b71cfe1cee0fbd3c3f Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 00:44:07 -0600 Subject: [PATCH 1/3] 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 aec01351161234ab3ef6da2211025c74ca9bf679 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Fri, 29 Nov 2024 18:51:02 -0600 Subject: [PATCH 2/3] simplified DOCSTRING for rank(::SVD) --- src/svd.jl | 25 ++++++------------------- 1 file changed, 6 insertions(+), 19 deletions(-) diff --git a/src/svd.jl b/src/svd.jl index b1e2976d..9e80a2bc 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -248,28 +248,15 @@ 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. +Compute the numerical rank of given an SVD object by counting how many singular values are greater +than `max(atol, rtol*σ₁)` where `σ₁` is the 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. +The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of UΣV' +and `ϵ` is the [`eps`](@ref) of the element type of `S`. +!!! compat "Julia 1.12" + The `rank(::SVD)` method requires at least Julia 1.12. """ - 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]) From 743b94dcca9b4f539755898ff5a782960affa307 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 30 Nov 2024 19:38:38 -0600 Subject: [PATCH 3/3] simplified rank(::SVD) implementation --- src/svd.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/svd.jl b/src/svd.jl index 9e80a2bc..ea9f279b 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -248,8 +248,9 @@ 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 given an SVD object by counting how many singular values are greater +Compute the numerical rank of the SVD object `S` by counting how many singular values are greater than `max(atol, rtol*σ₁)` where `σ₁` is the largest calculated singular value. +This is equivalent to the default [`rank(::AbstractMatrix)`](@ref) method except that it re-uses an existing SVD factorization. `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 UΣV' and `ϵ` is the [`eps`](@ref) of the element type of `S`. @@ -257,10 +258,9 @@ and `ϵ` is the [`eps`](@ref) of the element type of `S`. !!! compat "Julia 1.12" The `rank(::SVD)` method requires at least Julia 1.12. """ -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) +function rank(S::SVD; atol::Real = 0.0, rtol::Real = (min(size(S)...)*eps(real(float(eltype(S)))))) + tol = max(atol, rtol*S.S[1]) + count(>(tol), S.S) end ### SVD least squares ###