LAPACK Archives

[Lapack] psuedo inverse algorithm

Hi Adam,

Few comments.

- Please see my post on the LAPACK forum:
on how to compute the pseudo inverse with LAPACK. Essentially you can do this 
through the (1) QR factorization
of the matrix A (DGEQRF), (2) you can also use partial pivoting (DGEQP3), (3) 
you can also use SVD. Of course
(4) you can always go straight at the formula: ( A^T A )^(-1) A^T, but that 
will not be numerically stable if A is

- LAPACK does not have a routine to compute the pseudo-inverse of a matrix. 
From time to time we wonder if we should have such a routine,
No user has ever asked for a routine to compute the pseudo-inverse, (besides 
for the reason: "I want to compute the pseudo inverse, because
I want to"), and I am not aware of any application which requires it.

- The main advantage of (2, QR with pivoting) and (3, SVD) is that these 
methods enable filtering. So that you can discard small singular values
of A before inverting. This feature is needed for most applications.

- If you go with a QR based method (1 or 2), you number of flops would be: 
(assuming m>=n)
QR factorization: 2 mn^2 - 2/3n^3
Construct Q: 2 mn^2 - 2/3n^3
R \ Q^T: mn^2
So cost is 5 mn^2 - 4/3 n^3.
This counts multiplications and additions. You may divide by two for 
multiplications only.

- Now if we go with ( A^T A )^(-1) A^T, (not stable but let us do this anyway,)
A^TA : (SYRK ): mn^2
Cholesky factorization of A^T A : ( POTRF ): n^3 / 3
Cholesky solve ( A^T A )^(-1) A^T: 2mn^2
So cost is 3 mn^2 + 1/3 n^3.
This counts multiplications and additions. You may divide by two for 
multiplications only.

Now looking at the algorithm you are giving in your post:

- The algorithm does not look stable. I think it is as stable as computing the 
pseudo-inverse through the formula ( A^T A )^(-1) A^T. Essentially if you 
|| A^dagger * A - I || you will see a dependency in K( A )^2 where K( A ) is 
the condition number of the matrix A. A "stable" algorithm should get a decency 
in kappa(A).

- The algorithm has less flops when using elementary transformation, I believe 
we end up in the same number of flops for a QR factorization based method when
the algorithm uses Householder reflectors. LAPACK uses Householder reflectors 
for pretty much anything except LU factorization, Cholesky factorization, etc. 
For example,
we could very well use elementary operations for reduction to Hessenberg form, 
and this would pretty much valid, and halves the number of flops. We use
Householder reflectors for most of these operations (despite the double number 
of flops) because these transformations are always stable.

So all in all, this is a quick review of the algorithm. It is fun. I was not 
aware about it.
I am not sure we are interested in including it in LAPACK, the main reason is 
that (1) not sure there is a user base for such a functionality, (2) we can
obtain the pseudo-inverse from functionality in LAPACK, (3) the proposed 
algorithm is not the most stable one.


On Mar 2, 2013, at 9:17 AM, adam W 
<aswmac@Domain.Removed<mailto:aswmac@Domain.Removed>> wrote:


I have recently published an algorithm that computes the Moore-Penrose 
pseudoinverse, and does so in exact arithmetic. It accomplishes this with 2mn^2 
multiplications (dimensions with m>n), thus is comparable to a single matrix 
multiply. I believe your function makes use of singular value decomposition. 
Mine does not achieve the singular values, which is why it works quicker for 
just the pseudoinverse.

Here is the (thus far only) link to my algorithm:

I believe that this could enhance your wonderful software, how may I help this 
happen? Do I need to submit publication into a journal, or maybe write some 
code? (I have my code in python currently, but the algorithm is quite simple 
really and could easily be ported-- at least it would be easy for me to do.)


Lapack mailing list

-------------- next part --------------
An HTML attachment was scrubbed...

<Prev in Thread] Current Thread [Next in Thread>

For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or