How to obtain a pseudoinverse of an outer product of a matrix with itself given an inverse of its inner product?

Question

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

Answer 1

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!