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

1

Dear all,

Pseudo inverse (Moore-Penrose inverse), usually is called pinv in Matlab or some Python library, etc.
In Fortran, what is the best pseudo inverse (Moore-Penrose inverse) subroutine?

Thanks much in advance!

1 Like
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
3

Yes. Use one of the LAPACK suite of solvers.

If performance is an issue then make sure that you link to a BLAS implementation that is optimized for your platform and compiler.

1 Like
4

Dear all, sorry to revive this old thread -

I thought it’d be useful to link to my proposed implementation of the pseudo-inverse
at stdlib here for those interested.

Looking forward to your comments!

4 Likes
5

Thanks! Exciting!

I have a very basic and irrelevant question.
It looks like you write the code in fypp file. What IDE or code editor do you use to write code in fypp file?
I mean if I write f90 file I use visual studio, and it can easily compile and build the code and show errors messages if any, and I can easily locate to the lines which contain errors. For fypp, is there a IDE like visual studio? Thanks!

1 Like
6

Thank you @CRquantum, I use the Code::Blocks IDE. I have no experience with VS, but I think this trick applies to basically any editors: just set the syntax to Fortran for files with extension .fypp, that does the job for me (including linting, etc.). For preprocessing fypp files into .f90, you’ll have to script your IDE to call fypp behind the scenes though.

1 Like