Efficient projection onto the kernel of a matrix

Question

Suppose I have a positive semidefinite matrix $M = \sum_i^N A_i^T A_i$ where each $A_i$ is a fat matrix of shape (m,n) and $m << n$, we can also assume that $A_i$ is full rank(but the stacked matrix might not be). I can compute the projection onto the kernel of $M$ by computing the individual projections onto the kernel of $M$, given by $P_M = I - M^T(M M^T)^{-1} M$.

How can I improve this algorithm by exploiting the fact that the $m \times m$ system is easy to work with while keeping the memory cost low? Specifically, I would like a more stable and possibly more memory-efficient algorithm for this. Also, I only care about projection vector products, not the projection matrix itself which I can't store in memory.

Answer 1

If $n$ is sufficiently small you can compute $M= \sum_i A_i^T A_i$ which is $n\times n$. The projection on its kernel can be written as $P_M = I - MM^{+}$. So say you have a vector $v$ and you want to compute $P_Mv$ then you can first solve $M^TM x = M^Tv$ to get $x=M^{+}v$. If you initialize the conjugate gradient solver for the normal equations (CGNR/CGLS) with an intial guess of zero then you would get the desired solution. After that just compute $P_Mv = v - Mx$. If $n$ is sifficiently large and $M$ is too dense you can instead implement matrix-vector products with $M$ and $M^T$ through the $A_i$ without explicitly forming $M$.

Answer 2

The complement of the null space, $Z^\perp$, is of course spanned by the rows of the matrices $A_i$. If $Nm$ is relatively small, you can compute an orthonormal basis of $Z^\perp$ by running the Gram-Schmidt process on the vectors that correspond to the rows of the $A_i$, and once you have a basis of $Z^\perp$, you can easily project onto $Z^\perp$ (let's call this projection operator $\Pi_{Z^\perp}$). The projection $\Pi_Z$ onto the null space $Z$ of $M$ is then of course $I-\Pi_{Z^\perp}$.