Skip to content

Commit 7bfe65a

Browse files
stevengjKristofferC
authored andcommitted
clarify pinv documentation (#29793)
* clarify pinv documentation The `pinv` documentation falsely implied that it discarded singular values `σ > tol`, when in fact it discards `σ > tol max(σ)`. This PR corrects the docstring. (`σ > tol max(σ)` is a good behavior! `tol` should be a dimensionless quantity that doesn't depend on the overall scaling/units of the matrix.) * rename tol -> rtol (cherry picked from commit 2996521)
1 parent 7d1875b commit 7bfe65a

File tree

1 file changed

+17
-17
lines changed

1 file changed

+17
-17
lines changed

stdlib/LinearAlgebra/src/dense.jl

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1219,20 +1219,20 @@ factorize(A::Transpose) = transpose(factorize(parent(A)))
12191219
## Moore-Penrose pseudoinverse
12201220

12211221
"""
1222-
pinv(M[, tol::Real])
1222+
pinv(M[, rtol::Real])
12231223
12241224
Computes the Moore-Penrose pseudoinverse.
12251225
12261226
For matrices `M` with floating point elements, it is convenient to compute
1227-
the pseudoinverse by inverting only singular values above a given threshold,
1228-
`tol`.
1227+
the pseudoinverse by inverting only singular values greater than
1228+
`rtol * maximum(svdvals(M))`.
12291229
1230-
The optimal choice of `tol` varies both with the value of `M` and the intended application
1231-
of the pseudoinverse. The default value of `tol` is
1230+
The optimal choice of `rtol` varies both with the value of `M` and the intended application
1231+
of the pseudoinverse. The default value of `rtol` is
12321232
`eps(real(float(one(eltype(M)))))*minimum(size(M))`, which is essentially machine epsilon
12331233
for the real part of a matrix element multiplied by the larger matrix dimension. For
12341234
inverting dense ill-conditioned matrices in a least-squares sense,
1235-
`tol = sqrt(eps(real(float(one(eltype(M))))))` is recommended.
1235+
`rtol = sqrt(eps(real(float(one(eltype(M))))))` is recommended.
12361236
12371237
For more information, see [^issue8859], [^B96], [^S84], [^KY88].
12381238
@@ -1262,7 +1262,7 @@ julia> M * N
12621262
12631263
[^KY88]: Konstantinos Konstantinides and Kung Yao, "Statistical analysis of effective singular values in matrix rank determination", IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. [doi:10.1109/29.1585](https://doi.org/10.1109/29.1585)
12641264
"""
1265-
function pinv(A::StridedMatrix{T}, tol::Real) where T
1265+
function pinv(A::StridedMatrix{T}, rtol::Real) where T
12661266
m, n = size(A)
12671267
Tout = typeof(zero(T)/sqrt(one(T) + one(T)))
12681268
if m == 0 || n == 0
@@ -1273,7 +1273,7 @@ function pinv(A::StridedMatrix{T}, tol::Real) where T
12731273
maxabsA = maximum(abs.(diag(A)))
12741274
B = zeros(Tout, n, m)
12751275
for i = 1:min(m, n)
1276-
if abs(A[i,i]) > tol*maxabsA
1276+
if abs(A[i,i]) > rtol*maxabsA
12771277
Aii = inv(A[i,i])
12781278
if isfinite(Aii)
12791279
B[i,i] = Aii
@@ -1286,14 +1286,14 @@ function pinv(A::StridedMatrix{T}, tol::Real) where T
12861286
SVD = svd(A, full = false)
12871287
Stype = eltype(SVD.S)
12881288
Sinv = zeros(Stype, length(SVD.S))
1289-
index = SVD.S .> tol*maximum(SVD.S)
1289+
index = SVD.S .> rtol*maximum(SVD.S)
12901290
Sinv[index] = one(Stype) ./ SVD.S[index]
12911291
Sinv[findall(.!isfinite.(Sinv))] .= zero(Stype)
12921292
return SVD.Vt' * (Diagonal(Sinv) * SVD.U')
12931293
end
12941294
function pinv(A::StridedMatrix{T}) where T
1295-
tol = eps(real(float(one(T))))*min(size(A)...)
1296-
return pinv(A, tol)
1295+
rtol = eps(real(float(one(T))))*min(size(A)...)
1296+
return pinv(A, rtol)
12971297
end
12981298
function pinv(x::Number)
12991299
xi = inv(x)
@@ -1303,12 +1303,12 @@ end
13031303
## Basis for null space
13041304

13051305
"""
1306-
nullspace(M[, tol::Real])
1306+
nullspace(M[, rtol::Real])
13071307
13081308
Computes a basis for the nullspace of `M` by including the singular
1309-
vectors of A whose singular have magnitude are greater than `tol*σ₁`,
1309+
vectors of A whose singular have magnitude are greater than `rtol*σ₁`,
13101310
where `σ₁` is `A`'s largest singular values. By default, the value of
1311-
`tol` is the smallest dimension of `A` multiplied by the [`eps`](@ref)
1311+
`rtol` is the smallest dimension of `A` multiplied by the [`eps`](@ref)
13121312
of the [`eltype`](@ref) of `A`.
13131313
13141314
# Examples
@@ -1332,14 +1332,14 @@ julia> nullspace(M, 2)
13321332
0.0 0.0 1.0
13331333
```
13341334
"""
1335-
function nullspace(A::AbstractMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
1335+
function nullspace(A::AbstractMatrix, rtol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
13361336
m, n = size(A)
13371337
(m == 0 || n == 0) && return Matrix{eltype(A)}(I, n, n)
13381338
SVD = svd(A, full=true)
1339-
indstart = sum(s -> s .> SVD.S[1]*tol, SVD.S) + 1
1339+
indstart = sum(s -> s .> SVD.S[1]*rtol, SVD.S) + 1
13401340
return copy(SVD.Vt[indstart:end,:]')
13411341
end
1342-
nullspace(a::AbstractVector, tol::Real = min(size(a)...)*eps(real(float(one(eltype(a)))))) = nullspace(reshape(a, length(a), 1), tol)
1342+
nullspace(a::AbstractVector, rtol::Real = min(size(a)...)*eps(real(float(one(eltype(a)))))) = nullspace(reshape(a, length(a), 1), rtol)
13431343

13441344
"""
13451345
cond(M, p::Real=2)

0 commit comments

Comments
 (0)