How to obtain a pseudoinverse of an outer product of a matrix with itself given an inverse of its inner product? Given $A\in\mathbb{R}^{m\times n}$ where $m > n$, a tall skinny matrix, and the inverse $B = (A^TA)^{-1}$, is there a good way to obtain a pseudoinverse for $AA^T$?  Certainly, $AA^T$ is rank deficient, but it has the same eigenvalues as $A^TA$, so I'm hoping to cluster the eigenvalues of $AA^T$ by constructing some kind of pseudoinverse from $A^TA$.
 A: Evidently, I asked too soon.  From the Matrix Cookbook from Petersen and Pedersen, $(AA^T)^+=(A^T)^+A^+$.  Then, we have that $A^+=(A^TA)^{-1}A^T$ and $(A^T)^+=A(A^TA)^{-1}$.  Putting the two together, we have $(AA^T)^+=A(A^TA)^{-2}A^T$, which is what I wanted.

Edit 1
Per @TommasoF 's request, the solution above requires $A$ to have full rank.  The proof when it's full rank is actually easier and doesn't really require the pseudo inverse
$$
(A(A^TA)^{-2}A^T) (AA^T) = A(A^TA)^{-2}(A^TA)A^T = A(A^TA)^{-1}A^T
$$
but this is an orthogonal projector onto the range of $A$, so it has eigenvalues of 0 and 1, which is what I originally wanted.
Now, if $A$ is rank deficient, there's a trick using a modified QR factorization to pull out the right information.  Basically, a QR factorization with column pivoting will order the diagonal elements of $R$ in decreasing order.  When $R$ is rank deficient, they eventually become zero.  There's probably some nice way to do this with a rank revealing QR, but the code I posted in the question here seems to work.  Anyway, if we ignore pivoting, we basically take a modified QR factorization $QR=A^T$ where $A\in\mathbb{R}^{m\times n}$, $m>n$, $Q\in\mathbb{R}^{n\times r}$, $R\in\mathbb{R}^{r\times m}$, and $r$ is the rank of $A$.  Then
\begin{align*}
(R^T(RR^T)^{-2}R)(AA^T)
=& (R^T(RR^T)^{-2}R)(R^TQ^TQR)\\
=& (R^T(RR^T)^{-2}R)(R^TR)\\
=& R^T(RR^T)^{-2}(RR^T)R\\
=& R^T(RR^T)^{-1}R
\end{align*}
which is a projector and hence has eigenvalues of either 0 or 1.  Note, the inverse here is guaranteed to exist since $R$ is full rank.  Using the code that I linked to above, the following should demonstrate this
% Create a tall and skinny, rank deficient matrix
m = 10;
n = 5;
r = 3;
A = randn(m,r)*randn(r,n);

% Create a preconditioner that makes the rank deficient outter product
% spectrally simple
[q r p] = qrrd(A',1e-12);
P = r'*inv(r*r')^2*r;
[v d]=eig(P * A(p,:)*A(p,:)');
d=real(diag(d));
fprintf('Eigenvalues of Rt inv(R Rt)^2 R A(p,:) A(p,:)t: %d (0), %d (1)\n', ...
    length(find(abs(d)<1e-12)), length(find(abs(d-1)<1e-12)));

Anyway, there may be a better way to do it, but this seems to work.  If someone else has a better idea, please let me know!
