If we have $A=U\Sigma V^T$ as the SVD of $A$, what is the SVD of the Moore-Penrose Pseudoinverse $M=(A^TA)^{-1}A^T$? I've started to simplify it but eventually I get stuck.
$(A^TA)^{-1}A^T$
$((U\Sigma V^T)^TU\Sigma V^T)^{-1}(U\Sigma V^T)^T$
$((U\Sigma V^T)^TU\Sigma V^T)^{-1}(U\Sigma V^T)^T$
$(V\Sigma^T U^TU\Sigma V^T)^{-1}V\Sigma^TU^T$
$(V\Sigma^T \Sigma V^T)^{-1}V\Sigma^TU^T$
I don't know where to go from here. I need the expression to have a $N\times N$ matrix ($U$) followed by an $N \times d$ matrix ($\Sigma$) followed by a $d \times d$ matrix ($V^T$) but I can't make this expression into that.
If you look up the properties of the Moore-Penrose pseudo inverse $A^+$, you will find that, for example, if $A$ has orthonormal columns, or $B$ has orthonormal rows, then $(AB)^+ = B^+A^+$, just like the regular inverse. So, if $A = U\Sigma V^T$, you should easily be able to prove that $A^+ = V\Sigma^+U^T$ (recall that $U$ and $V$ are unitary so they are invertible with inverse equal to there conjugate transpose).
Now this could be enough for your purposes, but if you want to go deeper, $\Sigma^+$ is actually very easy to compute given $\Sigma$. Recall that $\Sigma$ is "diagonal" with the singular values along the diagonal. It turns out, and you should be able to prove that, to get the pseudo-inverse of such a matrix, simply transpose it, and then invert all the non-zero singular values.
*As extra intuition, verify that this answer makes geometric sense to you. Indeed, if you pretend for a second that $A$ is actually invertible, then this result is telling you that, if $A$ is a "rotation", followed by a scaling, followed by another rotation, then $A^{-1}$ should logically be the undoing of the last rotation, followed by inverting the scaling, followed by undoing the first rotation. Makes sense. Try to see if you can fit in zero singular values in that picture.