Solving Projection onto Linear Equality for Multiple Input Vectors

I am after solving the projection onto a Linear Equation problem:

\arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} \; \text{ subject to } \boldsymbol{A} \boldsymbol{x} = \boldsymbol{b}

with \boldsymbol{A} \in \mathbb{R}^{m \times n} with m \ll n.

The closed form solution is:

\begin{align*} \boldsymbol{x} & = \boldsymbol{y} - \boldsymbol{A}^{\top} {\left( \boldsymbol{A} \boldsymbol{A}^{\top} \right)}^{-1} \left( \boldsymbol{A} \boldsymbol{y} - \boldsymbol{b} \right) \\ & = \boldsymbol{y} - \boldsymbol{A}^{\top} {\left( \boldsymbol{A} \boldsymbol{A}^{\top} \right)}^{-1} \boldsymbol{A} \boldsymbol{y} + \boldsymbol{A}^{\top} {\left( \boldsymbol{A} \boldsymbol{A}^{\top} \right)}^{-1} \boldsymbol{b} \end{align*}

I wonder what the optimal pre processing needed to solve the above repeatedly for many cases of \boldsymbol{y} (The Matrix \boldsymbol{A} and vector \boldsymbol{b} are the same).

Specifically, how would you handle for the cases:

  • \boldsymbol{A} is a dense matrix.
  • \boldsymbol{A} is a sparse matrix.

For the case \boldsymbol{A} is dense, a simple solution is using the SVD of the matrix:

\begin{align*} \boldsymbol{A}^{\top} {\left( \boldsymbol{A} \boldsymbol{A}^{\top} \right)}^{-1} \boldsymbol{A} & = \boldsymbol{V} \boldsymbol{S} \boldsymbol{U}^{\top} {\left( \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^{\top} \boldsymbol{V} \boldsymbol{S} \boldsymbol{U}^{\top} \right)}^{-1} \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^{\top} \\ & = \boldsymbol{V} \boldsymbol{V}^{\top} \end{align*}
\begin{align*} \boldsymbol{A}^{\top} {\left( \boldsymbol{A} \boldsymbol{A}^{\top} \right)}^{-1} & = \boldsymbol{V} \boldsymbol{S} \boldsymbol{U}^{\top} {\left( \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^{\top} \boldsymbol{V} \boldsymbol{S} \boldsymbol{U}^{\top} \right)}^{-1} \\ & = \boldsymbol{V} \boldsymbol{S}^{\dagger} \boldsymbol{U}^{\top} \end{align*}