Solve $X^TX b = a$ for $b$ using $XX^T$ for a short and wide matrix $X$ I have a matrix $X$ of dimensions $n \times p$ and a fixed $p$-dimensional vector $a$, with $p \gg n$. How can I efficiently solve a problem of the following form? $$X^TXb = a$$
Perhaps using the $n \times n$ matrix $XX^T$?
 A: Expanding on @Analyst's approach:
Let $X = U \Sigma V^T$ be the compact SVD (where $U$ is $n \times r$, $\Sigma$ is $r \times r$ diagonal, $V$ is $p \times r$, with $r \le n$; probably you have $r = n$).
Then $X^T X = V \Sigma U^T U \Sigma V = V \Sigma^2 V^T$, and $(X^T X)^\dagger a = V \Sigma^{-2} V^T a$ is the minimum-norm solution to your problem. Given $V$ and $\Sigma$, this multiplication takes $O(n p + n^2)$ time.
Note that you can get the SVD of $X$ as follows:


*

*Compute $X X^T$, an $n \times n$ matix, in $O(n^2 p)$ time.

*Eigendecompose $X X^T = U \Sigma^2 U^T$, in $O(n^3)$ time.

*Compute $V$ by $X^T U \Sigma^{-1} = V \Sigma U^T U \Sigma^{-1} = V$, in $O(n^2 p)$ time.


Thus this method takes $O(n^3 + n^2 p)$ time, compared to $O(n p^2 + p^3)$ for multiplying out $X^T X$ and then Cholesky-factoring it: potentially enormous savings when $n \ll p$.
If that's still not fast enough, you could look into whether any iterative methods (e.g. Gauss-Siedel or conjugate gradient) can run efficiently in this case....
A: You could use an generalized inverse via minimun norm solution.  
https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse#Applications
