# Fastest algorithm for pseudoinverse of skinny matrices

For a performance-sensitive problem, I need to compute the pseudoinverse of a skinny matrix (#rows = 1000–10000, #cols= 10–20).

I already employ the traditional SVD econ method. For some problem sizes, this takes up most of my computational time. Is there any faster method? The matrix is dense, and usually ($$A^TA$$) is well-conditioned

• Do you perhaps mean $\mathbf A^T \mathbf A$ is well-conditioned? (For tall/skinny $\mathbf A$, the "outer product" $\mathbf A \mathbf A^T$ would be singular). – rchilton1980 Sep 5 at 14:23
• I corrected the typo. – Enzo Ferrazzano Sep 6 at 15:31

If $$A$$ is of full column rank and $$A^{T}A$$ is non-singular and well-conditioned, then you can compute the pseudoinverse as:

$$A^{\dagger}=(A^{T}A)^{-1}A^{T}$$.

This will be faster than computing a QR or SVD factorization of $$A$$ but be careful about the conditioning of $$A^{T}A$$.

• At least using IntelMKL & C++, they seems to be of the same speed. Perhaps under the hood they all use SVD. – Enzo Ferrazzano Sep 6 at 15:30
• @EnzoFerrazzano if you're calling MKL routines to do the inverse of $A^{T}A$, then you're not computing an SVD. My testing in MATLAB found this formula to be about 4 times faster than MATLAB's SVD based pinv() function. I'd be interested in seeing your C++ code. Note that your matrices are so small, that this will run very fast and it can be difficult to time it. – Brian Borchers Sep 6 at 15:48
• Note that you have to be a bit careful with how you code it. Are you using syrk rather than gemm to compute $A^TA$? Are you using sysv rather than getri for the system solve? – Federico Poloni Sep 8 at 14:45
• I was finally able to get the speedup. Problem was that I was mixing Armadillo (which choses the algorithm, and seems to use SVD in many places) and LAPACK, so performance was inconsistent. @Federico: will try. My matrix are very skinny (10-20 cols), so the inversion is not a big problem. I will test them nonetheless. For sysv, I though of it, but since the RHS is A^T, gemm takes care it. Is there a way to achieve that with sysv? – Enzo Ferrazzano Sep 13 at 9:51
• Yes, just compute each column of the product by solving a symmetric linear system instead of computing the inverse. But I agree that probably it won't matter if the matrix is very thin. Exploiting symmetry in the product (syrk), on the other hand, reduces the total complexity by a factor 2, so make syre you don't miss it. – Federico Poloni Sep 13 at 10:01

You can "frontend" the SVD of a tall/skinny matrix using a QR decomposition, then just SVD the remaining small/square matrix R. Here's a snippet of matlab code that computes the pseudoinverse based on this idea:

clear all
close all

% Make input.
I = 10000;
J = 20;
A = rand(I,J);

% Pseudoinverse B := pinv(A)
tic
B = pinv(A);
toc
norm(A - A*B*A,'fro')
norm(B - B*A*B,'fro')

% Pseudoinverse C := V*inv(S)*(Q*U)', based on qr(A) and svd(R).
tic
[Q,R] = qr(A,0);
[U,S,V] = svd(R);
C = V*inv(S)*(Q*U)';
toc
norm(A - A*C*A,'fro')
norm(C - C*A*C,'fro')


Example output:

Elapsed time is 0.0120101 seconds.
ans =   6.7672e-013
ans =   2.1392e-016
Elapsed time is 0.00999999 seconds.
ans =   6.7672e-013
ans =   2.1392e-016


Maybe it was a little faster? Though to be honest, a high quality implementation of SVD would already incorporate this idea. This might be what you meant already by "traditional SVD econ method", just thought I'd point it out.