diff --git a/src/svd.jl b/src/svd.jl index 7a88c4a6..ea9f279b 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -245,6 +245,24 @@ 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 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`. + +!!! compat "Julia 1.12" + The `rank(::SVD)` method requires at least Julia 1.12. +""" +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 ### 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