Documentation of `pinv` is wrong; also, methods behave differently

[Clemens-H]

(1) The documentation of pinv does not apply to all methods. In the general, the SVD is used to calculate the pseudoinverse as documented. There is a tolerance when elements in S are considered zero.

julia> pinv([1.0 0.0; 0.0 1.0e-17])
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  0.0

However, there are also special methods for diagonal matrices. The method pinv(D::Diagonal{T}) doesn't use any tolerance, which is different from the behavior of the general case.

julia> pinv(Diagonal([1.0, 1.0e-17]))
2×2 Diagonal{Float64}:
 1.0   ⋅    
  ⋅   1.0e17

Wouldn't it be better / more consistent to use optional arguments with default values for all methods?

(2) The documentation of pinv says the default value of the tolerance is eps(real(float(one(eltype(M)))))*maximum(size(A)). However, in the method pinv(A::StridedMatrix{T}, tol::Real), the values actually used are tol*maxabsA and tol*maximum(SVD.S) (where tol = eps(real(float(one(T))))*maximum(size(A))), which is good.

[stevengj]

+1 for using the same tolerance keywords for the Diagonal case.

We should also correct the docs, of course. Note that the Matlab pinv function uses a default tolerance of maximum(size(A))*eps(norm(A)), while the NumPy pinv function (also SciPy) uses roughly eps(norm(A)). I'm not sure what the rationale is for our current choice? At first glance, NumPy's choice makes more sense to me … I'm not sure why maximum(size(A)) comes in, but probably I'm missing something? Maybe it comes from some forward error bound on the svdvals?

[stevengj]

@alanedelman, do you know the story behind the default tolerance of pinv?

[jiahao]

Previous work on pinv tolerance and related least squares solvers: #59 JuliaLang/julia#8859 JuliaLang/julia#9008 JuliaLang/julia#9101

[stevengj]

@jiahao, the question is, why does the default rcond threshold go like eps(T)*maximum(size(A))? In NumPy and SciPy it is basically just eps(T) — where does the size of the matrix come in? I did a quick scan of the issues you linked and didn't see an explanation; did I miss something?

@Clemens-H, the case where tol*maxabsA is used is that of a diagonal matrix, where maxabsA is the same as norm(A). But I agree that the pinv documentation is misleading: the documentation says "inverting only singular values above a given threshold,
tol" which implies that tol is the threshold for the singular values, rather than for the singular values as a fraction of the largest singular value, i.e. the rcond (reciprocal conditional number).

[jiahao]

We had a long discussion in JuliaLang/julia#8859, leading to the current compromise default choice in JuliaLang/julia#9008. We've already had Per-Christian Hansen weigh in in JuliaLang/julia#8859 to confirm that there really is no good default choice and we're providing one purely for convenience based on the cases we studied in that issue.