Should \(SVD,...) truncate small singular values?

[andreasnoack]

This is a follow up on the discussion in JuliaLang/julia#32126. Currently, \ and inv truncate the smaller singular values, i.e. it actually computes a pseudo-inverse (or applies regularization if that is your preferred language). That is sometimes beneficial and the SVD is used for this purpose but I'm wondering if we should make the truncation/regularization explicit, i.e. let inv be the untruncated inverse and use pinv for the truncated version. Similarly, we could add a relative tolerance keyword argument to ldiv! which would default to zero, which would mean no truncation. The issue is \ which is an operator so we can't really control the behavior with an extra argument (when used as infix). Thoughts?

[clason]

Nitpick: "regularization" is actually the correct language -- the truncation level corresponds to a regularization parameter. The pseudoinverse corresponds to inverting all non-zero singular values (and is therefore not numerically stable in the presence of very small singular values -- the only benefit of the pseudoinverse over the actual inverse is that it always exists even if the matrix is not invertible, in which case it corresponds to the least-squares solution of minimal norm).

As a mathematician, Matlab's silent use of (regularized) pseudoinverses when asked to solve a linear system has always bugged me to no end. My personal vote would therefore be to be fully explicit:

• Make inv invert all singular values and either warn or error if zeros are encountered.
• Make pinv invert all non-zero singular values (possibly warning if very small values are encountered).
• Add a keyword to pinv (but not inv) that allows passing a truncation level (defaulting to zero).
• Have \ behave as inv, but always error in the presence of zero values.
• I have no strong opinion about ldiv! either way, but it might be useful to discuss what the tolerance keyword should control (and what "relative" is relative to -- does it involve the specific right-hand side?)

The error messages could (and should) then point towards the regularized version. This way users will know what they are doing and not rely on silent "do-what-I-mean" (where "I" is the software) behavior.

If you want to keep the default truncation (replacing "non-zero" above by "numerically non-zero"), I would advocate changing the default value to be absolute and related to the floating-point accuracy but not the dimension.