linalg: Moore-Penrose pseudo-inverse (pinv)

[perazz]

Compute pseudo-inverse of a square or rectangular matrix using singular value decomposition (SVD).
Array $A$, real or complex, has rank n==2.

The pseudo-inverse $A^{+}$ is a generalization of the matrix inverse and satisfies the following properties:

• $A \cdot A^{+} \cdot A = A$
• $A^{+} \cdot A \cdot A^{+} = A^{+}$
• $(A \cdot A^{+})^H = A \cdot A^{+}$
• $(A^{+} \cdot A)^H = A^{+} \cdot A$

The user may specify an eigenvalue cutoff relative tolerance, rtol.

Proposed interface

• ap = pinv(a [, rtol] [, err]): function interface
• call pseudoinvert(a, ap [, rtol] [, err]): subroutine interface (preallocated ap)
• .pinv.a: operator interface

Key facts

• The implementation is made to match that of inv, invert and .inv.
• Fortran Discourse thread
• Use SVD for performance: PERF: scipy.linalg.pinv is very slow compared to numpy.linalg.pinv scipy/scipy#12196
• if rtol is not provided, just use $\text{max}(m,n)\cdot\varepsilon$

Progress

• base implementation
• tests
• documentation
• submodule
• examples

Prior art

• Numpy: linalg.pinv(a, rcond=None, hermitian=False, *, rtol=<no value>)
• Scipy: pinv(a, *, atol=None, rtol=None, return_rank=False, check_finite=True)
• Matlab: B = pinv(a [,tol])

cc: @jvdp1 @jalvesz @Beliavsky @CRquantum @fortran-lang/stdlib

[jvdp1]

Thank you @perazz . It looks really good to me!

[perazz]

Thank you for the reviews and the approvals @jalvesz @jvdp1, I believe all fixes are addressed and no comments for a while, so I will merge this.