Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
18 changes: 18 additions & 0 deletions src/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
12 changes: 12 additions & 0 deletions test/svd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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