What is the best pseudo inverse (Moore-Penrose inverse) subroutine?

2

DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices.

If you really do need the inverse (see my reply in Why doesn't Fortran have a built-in matrix inverse function? - #4 by ivanpribec), then do your SVD / QR / LQ factorization first and then solve A X = I.

Note for example the Eigen C++ library has a complete orthogonal decomposition routine, which also provides a pseudo inverse method, however the documentation states explicitly:

Do not compute this->pseudoInverse()*rhs to solve a linear systems. It is more efficient and numerically stable to call this->solve(rhs) .

1 Like