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.


2 Answers 2


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."

  • $\begingroup$ Thanks Prof, thie matlab code is so clear that I understand that Krylov subspace methods still work for consistent singular systems. Thanks very much for your careful reply. $\endgroup$
    – Happy
    Oct 26, 2019 at 10:43

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.

  • $\begingroup$ thanks Prof, for Krylov subspace methods, like CG,GMRES,..., if we get a solution for the singular consitent system Ax=b, the solution is a least square solution because of infinite many solutions, right? $\endgroup$
    – Happy
    Oct 24, 2019 at 23:04
  • $\begingroup$ No, it's one out of the infinitely many solutions, but it is not necessarily the least squares solution. $\endgroup$ Oct 25, 2019 at 1:32
  • $\begingroup$ thanks Prof, in the monography, most authors do not discuss the singular linear system. Is this because of less theories for singular system? I find many monographies assume that matrix A is nonsingular before they discuss so many iterative methods. So, do we need to consider the singular system? or don't care the singular system? because in some matlab codes, I find the author still use MINRES to compute the iterative solution of the Stokes equations [A,B;B',O]*[v;p]=b, where the coefficient matrix is singular in enclosed fluid case. So, I am a little confused whether it works. $\endgroup$
    – Happy
    Oct 25, 2019 at 2:33
  • $\begingroup$ Generally, if you are required to solve a linear system that is singular, you've made a mistake in deriving your model or in implementing it. That's because every physical system should have a unique solution. $\endgroup$ Oct 26, 2019 at 14:51
  • $\begingroup$ thanks Prof, but for Poisson with pure Neumann boundary conditions has many solutions not a unique solution. I donot understand the 'unique solution' you said. and stokes with enclosed fluid also has many solution because the matrix is singular. $\endgroup$
    – Happy
    Oct 27, 2019 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.