1
$\begingroup$

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

Improve this question
$\endgroup$
  • 2
    $\begingroup$ 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). $\endgroup$ – rchilton1980 Sep 5 '19 at 14:23
  • $\begingroup$ I corrected the typo. $\endgroup$ – Enzo Ferrazzano Sep 6 '19 at 15:31
5
$\begingroup$

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$.

Improve this answer
$\endgroup$
  • $\begingroup$ At least using IntelMKL & C++, they seems to be of the same speed. Perhaps under the hood they all use SVD. $\endgroup$ – Enzo Ferrazzano Sep 6 '19 at 15:30
  • 1
    $\begingroup$ @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. $\endgroup$ – Brian Borchers Sep 6 '19 at 15:48
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Sep 8 '19 at 14:45
  • $\begingroup$ 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? $\endgroup$ – Enzo Ferrazzano Sep 13 '19 at 9:51
  • $\begingroup$ 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. $\endgroup$ – Federico Poloni Sep 13 '19 at 10:01
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.