Null Space Projection for Singular Systems

Question

Let $A$ be a general symmetric matrix operator and $P$ be the unique orthogonal projection onto $\operatorname{Range}(A) = \operatorname{Null}(A)^\perp$.

Analytically, the system $$Ax = Pb$$ should have a solution for any vector $b$. If we further stipulate that this solution lie in $\operatorname{Range}(A)$, the solution should be unique. (And in the positive definite case, this is precisely the solution that the conjugate gradient algorithm will find in exact arithmetic.)

The above is a common tactic when solving the Poisson equation with either Neumann or periodic boundary conditions; if $A$ is taken to be the Neumann or periodic Laplacian, then $$ P = I - \frac{1}{N}(1_N)(1_N)^T,$$ where $1_N = [1, \cdots, 1]^T$, and we can solve $Ax = Pb$ using a conjugate gradient approach.

But what if the projection operator $P$ doesn't have a simple analytic expression? Can it be calculated numerically from $A$? Is there a way to solve this type of singular problem in the case that $A$ is a black-box symmetric operator?

I realize one approach would be to solve the normal equations $A^2x = Ab$, which avoids using a projection operator entirely, but this also severely worsens the conditioning of the system, so I'm interested to know if there's another way.

Answer 1

Computing the null space without any further knowledge is expensive in general. There are some iterative methods that converge to minimum-norm solutions even when presented with inconsistent right hand sides. Choi, Paige, and Saunders' MINRES-QLP is a nice example of such a method. For non-symmetric problems, see Reichel and Ye's Breakdown-free GMRES.

In practice, usually some characterization of the null space is important for effective preconditioning. Since most practical problems require preconditioning, the purely iterative methods have seen limited adoption. Note that in case of very large null space, preconditioners will often be used in an auxiliary space where the null space has been removed. See the "auxiliary-space Maxwell" methods for more details.

Answer 2

Jed has already answered your question, implicitly assuming that your matrix is large and sparse. If this should not be the case, you can do an eigenvalue decomposition and write $A=X\Lambda X^T$ where the columns of $X$ correspond to the eigenvectors and the entries of the diagonal matrix $\Lambda$ are the eigenvalues. (Note that this decomposition is not easily computable if $A$ is large, as Jed mentioned.) Then, the null space consists of exactly those vectors spanned by the eigenvectors corresponding to zero (or very small) eigenvalues, i.e., the projector onto the null space can be written as $$ 1_N = X \tilde\Lambda X^T $$ where $(\tilde\Lambda)_{ii} = 0$ if $|\Lambda_{ii}|\ge \varepsilon$ and $(\tilde\Lambda)_{ii} = \Lambda_{ii}$ if $|\Lambda_{ii}|<\varepsilon$.