Does iterative method work for singular consistent linear system Ax=b?

Question

Recently, I have been studied iterative methods for large sparse linear system Ax=b, where A is nonsingular, so there is a unique solution x. And the stopping criterion is usually chosen with norm(b-Ax_k)/norm(b)< tol (may 1e-6).

But some situations may happen that matrix A is singular, for example, when solving Poisson equation with Neumann boundary condition. My question is that for those singular linear system, can we still use iterative method? Also, in those singular system, the right hand side $b\in range(A)$, i.e., b is in the range of matrix A such that the system Ax=b is a consistent one.

Iterative methods like Jacobi, SOR, CG, minres or gmres still can be used to compute the solution? but in this case, there are infinite solution. how do the iterative methods compute a solution? thanks very much.

Answer 1

When using Krylov solvers to solve $\mathbf A \mathbf x = \mathbf b$ systems, singularity in $\mathbf A$ can be pretty benign as long as $\mathbf b$ is consistent. Consider applying the conjugate gradient (CG) algorithm to a symmetric positive semidefinite $\mathbf A$ (singular, but all non-zero eigenvalues still positive). Further assume $\mathbf b$ is in the range of $\mathbf A$, and we pick the initial guess to be $\mathbf x = \mathbf 0$. Then, the search vector/residual $\mathbf r = \mathbf b - \mathbf A \mathbf x$ will always be orthogonal to the nullspace of $\mathbf A$ (because in CG, search vectors/residuals are always drawn by applying $\mathbf A$ to previous residuals, and this multiplication by $\mathbf A$ can't generate any component in the nullspace). The algorithm can never "sense" the singularity of $\mathbf A$, and will converge to the (unique!) $\mathbf x$ such that $\mathbf A \mathbf x = \mathbf b$ and $\mathbf x \perp \mathrm {null}(\mathbf A)$

A short demo in matlab:

clear all
close all

% Form positive semidefinite A.
N = 10;
[Q,~] = qr(rand(N));
D1 = linspace(1,N/2,N/2);
D2 = zeros(1,N/2);
D = diag([D1 D2]); 
A = Q*D*Q';
kappa = cond(A) % infinite/singular

% Form random b in range(A).
b = A * rand(N,1);

% Solve A*x=b using CG.
x = pcg(A,b);

% Check properties.
residual = norm(A*x-b)
nullity = norm(x'*null(A))

Similar arguments can be made for other Krylov solvers, but I can't really speak to "classical" iterative solvers (Gauss-Seidel etc) .. I doubt they are as well-behaved.

There are preconditioning strategies (deflation) based entirely on this idea: find a basis for the eigenvectors of $\mathbf A$ that are somehow problematic to the solver, then "deflate"/project them out (change the operator such that the corresponding eigenvalues are mapped to zero), then apply the same projection to $\mathbf b$, and iterate away. For more information about deflation preconditioning (which involves convergence analysis for singular systems), I recommend:

Tang, JM, et al. "A comparison of two-level preconditioners based on multigrid and deflation."

Erlangga, Yogi A., and Reinhard Nabben. "Deflation and balancing preconditioners for Krylov subspace methods applied to nonsymmetric matrices."

Answer 2

Yes, many iterative methods still converge, or at least break down when they have exhausted all search directions corresponding to the span of the eigenvectors corresponding to non-zero eigenvalues. For example, you can use CG for the pure Neumann problem.

Because all of these methods are deterministic, if they converge, they will converge to a particular element of the solution set. For example, for the CG method applied to the pure Neumann problem, I suspect that one could show that the mean value of the elements of the solution vector equal the mean value of the elements of the starting guess because one only ever adds vectors to the iterates that have mean value zero. In other words, the vector you get out at the end will be one particular solution among infinitely many; whether it's the one you want is a different question, of course.