Projected Conjugate Gradient

Let $A\in {\mathbb{R}}^{n\times n}$ be symmetric positive semidefinite with null space $V=\ker A$ (orthonormal basis $V\in {\mathbb{R}}^{n\times k}$), and positive definite on $V^{\perp }=\mathrm{range}A$. Equivalently, $A$ is an SPD operator on the quotient space ${\mathbb{R}}^n/V$. For the pure Neumann Weighted Stiffness Matrix, $V=\mathrm{span}\{1\}$. Let

$$\Pi =I-V{V}^{\top }$$

be the orthogonal projector onto $V^{\perp }$.

Does plain CG work? In exact arithmetic, yes, as long as $\mathbf{b}\in \mathrm{range}A$ and ${\mathbf{x}}_0\in V^{\perp }$ (e.g. ${\mathbf{x}}_0=0$). All residuals and search directions then stay in $V^{\perp }$, where $A$ is SPD, and CG converges to the solution orthogonal to $V$ (the minimum-norm solution). The convergence rate depends on the effective condition number ${\lambda }_{\max }/{\lambda }_{\min }^{+}$ over the nonzero eigenvalues. Nothing in the CG recurrences is special to the nonsingular case. In practice there are three failure modes:

1. Inconsistent right-hand side. If $\mathbf{b}\notin \mathrm{range}A$, for example because the discrete load vector of a boundary current does not sum exactly to zero, then the component $V^{\top }\mathbf{b}$ can never be reduced. The residual stagnates at $\Vert V^{\top }\mathbf{b}\Vert$ and the iterate drifts along $V$ without bound. The tolerance is never reached.
2. Round-off. Even for consistent data, rounding errors introduce small null-space components into the residuals and directions, which accumulate over many iterations.
3. Preconditioners. $M^{-1}$ need not map $V^{\perp }$ into itself (Jacobi does not, AMG and incomplete factorisations generally do not). The preconditioned directions then acquire null-space components.

Projected CG removes all three by applying $\Pi$ where the iteration could leave $V^{\perp }$:

$$\begin{array}{rl}& {\mathbf{r}}_0=\Pi (\mathbf{b}-A{\mathbf{x}}_0),{\mathbf{z}}_0=\Pi M^{-1}{\mathbf{r}}_0,{\mathbf{p}}_0={\mathbf{z}}_0,\\ & {\alpha }_k=\frac{{\mathbf{r}}_k^{\top }{\mathbf{z}}_k}{{\mathbf{p}}_k^{\top }A{\mathbf{p}}_k},{\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k,{\mathbf{r}}_{k+1}={\mathbf{r}}_k-{\alpha }_{k}A{\mathbf{p}}_k,\\ & {\mathbf{z}}_{k+1}=\Pi M^{-1}{\mathbf{r}}_{k+1},{\beta }_k=\frac{{\mathbf{r}}_{k+1}^{\top }{\mathbf{z}}_{k+1}}{{\mathbf{r}}_k^{\top }{\mathbf{z}}_k},{\mathbf{p}}_{k+1}={\mathbf{z}}_{k+1}+{\beta }_k{\mathbf{p}}_k.\end{array}$$

• With inconsistent data it solves the nearest consistent problem $A\mathbf{x}=\Pi \mathbf{b}$, the least-squares solution. The removed part $\Vert (I-\Pi )\mathbf{b}\Vert /\Vert \mathbf{b}\Vert$ should be reported as a compatibility defect.
• The effective preconditioner is $\Pi M^{-1}\Pi$. It must be symmetric positive definite on $V^{\perp }$, which holds for a symmetric V-cycle or Jacobi.
• For $V=\mathrm{span}\{1\}$, applying $\Pi$ just subtracts the mean: $O(n)$ per application.
• Periodically recomputing the true residual $\Pi (\mathbf{b}-A{\mathbf{x}}_k)$ limits drift of the updated residual.

Choosing the representative. Projected CG returns the solution orthogonal to $V$. A different grounding, such as a zero sum over the boundary nodes, is obtained afterwards by an oblique projection along $V$. This is cheaper and better conditioned than building the constraint into the matrix, which is what pinning a node or adding $\epsilon I$ do (see Null Space of the Neumann Problem).

For many right-hand sides at once see Block Conjugate Gradient.

In ModularEIT.jl: pbcg, pbcg!.

References

1. E. F. Kaasschieter (1988). Preconditioned conjugate gradients for solving singular systems. J. Comput. Appl. Math. 24(1–2), 265–275. doi:10.1016/0377-0427(88)90358-5
2. P. Bochev, R. B. Lehoucq (2005). On the Finite Element Solution of the Pure Neumann Problem. SIAM Review 47(1), 50–66. doi:10.1137/S0036144503426074
3. M. R. Hestenes, E. Stiefel (1952). Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand. 49(6), 409–436. doi:10.6028/jres.049.044