Skip to content

Commit

Permalink
Document numerical error in rank (#48127)
Browse files Browse the repository at this point in the history
Co-authored-by: Daniel Karrasch <@dkarrasch>
Co-authored-by: Fredrik Ekre <ekrefredrik@gmail.com>
Co-authored-by: Steven G. Johnson <stevenj@alum.mit.edu>
  • Loading branch information
3 people committed Feb 18, 2023
1 parent 12d329b commit f64463d
Showing 1 changed file with 13 additions and 4 deletions.
17 changes: 13 additions & 4 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -947,13 +947,22 @@ dot(x::AbstractVector, transA::Transpose{<:Real}, y::AbstractVector) = adjoint(d
rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
rank(A::AbstractMatrix, rtol::Real)
Compute the rank of a matrix by counting how many singular
values of `A` have magnitude greater than `max(atol, rtol*σ₁)` where `σ₁` is
`A`'s largest singular value. `atol` and `rtol` are the absolute and relative
Compute the numerical rank of a matrix by counting how many outputs of
`svdvals(A)` 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
Expand Down Expand Up @@ -981,7 +990,7 @@ function rank(A::AbstractMatrix; atol::Real = 0.0, rtol::Real = (min(size(A)...)
isempty(A) && return 0 # 0-dimensional case
s = svdvals(A)
tol = max(atol, rtol*s[1])
count(x -> x > tol, s)
count(>(tol), s)
end
rank(x::Union{Number,AbstractVector}) = iszero(x) ? 0 : 1

Expand Down

0 comments on commit f64463d

Please sign in to comment.