# preconditioner for Laplace "without" boundary values

I'm looking at solving systems with the FEM discretization $$-\int_\Omega (\Delta u) v = \int_\Omega \nabla u \cdot \nabla v - \int_\Gamma (n\cdot\nabla u) v.$$ without applying Dirichlet- or Neumann-type boundary conditions. The resulting matrix is generally not self-adjoint (except in for 1D meshes).

The kernel consists of all linear functions over the domain, so there are always $$n+1$$ eigenvalues 0 (where $$n$$ is the dimensionality of the mesh).

The right-hand side is such that the system is consistent and I would like to find a solution of the system. The idea is to use GMRES, starting off with an all-zero initial guess.

Is there a good preconditioner for this known in literature? I played around with Dirichlet- and Neumann-Laplace without much success.

To get a feel, here's how to create the matrix and plot the spectrum in sckit-fem: import matplotlib.pyplot as plt
import meshzoo
import numpy as np
import skfem

@skfem.BilinearForm
def laplace(u, v, _):

@skfem.BilinearForm
def flux(u, v, w):

points, cells = meshzoo.disk(6, 20)

pT = np.ascontiguousarray(points.T)
cT = np.ascontiguousarray(cells.T)
mesh = skfem.MeshTri(pT, cT)

element = skfem.ElementTriP1()

basis = skfem.CellBasis(mesh, element)
facet_basis = skfem.FacetBasis(basis.mesh, basis.elem)

lap = skfem.asm(laplace, basis)
boundary_terms = skfem.asm(flux, facet_basis)

A = lap - boundary_terms

out = np.linalg.eigvals(A.toarray())
plt.plot(out.real, out.imag, "o")
plt.savefig("out.png", bbox_inches="tight")
plt.show()

• Have you tried multi-grid? Mar 29 at 2:41
• The question is as ill-posed as your problem. If you say you want to solve a linear system with this matrix, then this will not work: You have a singular matrix, so there will be infinitely many solutions. If you have a specific solution in mind, you have to add this information either to your system, or at least to the question here. Mar 29 at 3:32
• @AloneProgrammer The matrix isn't even symmetric so that's going to be tough. Unless you have a hint? Mar 29 at 7:43
• @WolfgangBangerth I clarified the question: I'm just looking for a solution, and I usually do that by starting off with Krylov and using an all-0 initial guess. Mar 29 at 7:44
• @NicoSchlömer OK, got it. I suggest you to use PETSC with a multi-grid preconditioner. The fact that your matrix is not symmetric is not a big problem. At least giving it a shot worth to try. By the way, if you realized multi-grid is not a good choice here there are plenty of other choices available in PETSC such as SOR. Mar 29 at 14:24

Here is at least an idea, whether it works is a different question.

Let's say you sort unknowns so that you have the ones in the interior of the domain first, and then all those at the boundary. Then the matrix that corresponds to your problem decomposes in the following way: $$A = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ C^{\partial\circ} & D^{\partial\partial} \end{pmatrix}$$ where $$\circ$$ indicates shape functions in the interior and $$\partial$$ on the boundary of the domain. $$A^{\circ\circ}$$ is simply the Dirichlet-boundary condition matrix for the Laplace operator and we know how to invert it efficiently. One could then think about using a Schur complement approach where you first solve for the boundary unknowns and then the interior unknowns, or maybe the other way around.

If you follow the arguments of the Silvester-Wathen preconditioner for the Stokes system, you will also be able to construct a good preconditioner for the whole matrix based on the Schur complement approach; it will probably look something like this: $$P^{-1} = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ 0 & S^{\partial\partial} \end{pmatrix}^{-1}$$ where $$S^{\partial\partial}=D^{\partial\partial}-D^{\partial\circ}[A^{\circ\circ}]^{-1}B^{\circ\partial}$$ is the Schur complement. It would be worthwhile thinking about whether one can approximate $$S^{\partial\partial}$$ by a simpler matrix, in the same way as for the Stokes equations, one approximates the Schur complement by the mass matrix.

As others have pointed out, (algebraic) multigrid can actually be a good preconditioner in this scenario. Below is a proof-of-concept implementation with scikit-fem and pyamg. It shows that pyamg's preconditioner makes the number of GMRES iterations more or less independent of the number of unknowns.

I only tested this with a mesh of a few thousand unknowns. The main computational cost is computing the left nullspace of A, used for making the right-hand side consistent. (A computation that isn't necessary in the actual application where consistency is guaranteed.) import matplotlib.pyplot as plt
import meshzoo
import numpy as np
import pyamg
import scipy.linalg
import scipyx
import skfem as fem
from skfem.helpers import dot
from skfem.models.poisson import laplace

rng = np.random.default_rng(0)

tol = 1.0e-8

for n in range(5, 41, 5):
# points, cells = meshzoo.rectangle_tri((0.0, 0.0), (1.0, 1.0), n)
points, cells = meshzoo.disk(6, n)
print(f"{n = }, {len(points) = }")

@fem.BilinearForm
def flux(u, v, w):

mesh = fem.MeshTri(points.T.copy(), cells.T.copy())
basis = fem.InteriorBasis(mesh, fem.ElementTriP1())
facet_basis = fem.FacetBasis(basis.mesh, basis.elem)

lap = fem.asm(laplace, basis)
boundary_terms = fem.asm(flux, facet_basis)

A = lap - boundary_terms

b = rng.random(A.shape)
# make the system consistent by removing A's left nullspace components from the
# right-hand side
lns = scipy.linalg.null_space(A.T.toarray()).T
for n in lns:
b -= np.dot(b, n) / np.dot(n, n) * n

ml = pyamg.smoothed_aggregation_solver(
A,
coarse_solver="jacobi",
symmetry="nonsymmetric",
max_coarse=100,
)
M = ml.aspreconditioner(cycle="V")

_, info = scipyx.gmres(A, b, tol=tol, M=M, maxiter=20)
# res = b - A @ info.xk

num_unknowns = A.shape
plt.semilogy(
np.arange(len(info.resnorms)), info.resnorms, label=f"N={num_unknowns}"
)

plt.xlabel("step")
plt.ylabel("residual")
plt.legend()
plt.savefig("out.png", bbox_inches="tight")
plt.show()


If you are not enforcing any boundary conditions, you should remove the nullspace of the problem (to make sure that there is a unique solution, MatSetNullSpace can be used to achieve this in PETSc) and then use multigrid on the problem (for example, PCMG of PETSc).

What surprises me is that there are n+1 zero eigenvalues. I would expect only one. Are you solving the Laplace's equation on surfaces (like a 2D surface embedded in 3D space) and the boundary is the boundary of that surface (like 1D boundary embedded in 3D space)? That is the only scenario I can imagine that the nullspace is that populous.

• The kernel consists of all linear functions over the domain, that's how you can see that dimensionality. (I've added this to the original question.) Q: Why should you remove the null space before applying multigrid? And what makes you think multigrid is a good choice here? The system has complex-valued eigenvalues after all. Mar 30 at 7:00
• Heuristically, by removing the nullspace you are preventing false divergence which results from the accumulation of coefficients corresponding to the basis vectors of the kernel. More mathematically, but still in a handwavy way, what you need for multigrid convergence is ellipticity on the kernel of the operator -which you have- (I couldn't find a direct citation for this, but check Boffi, Brezzi and Fortin's Mixed FEM and Applications). Removing the nullspace allows for effective restriction and prolongation. Lastly, I don't know whether multigrid is the best choice. But it is a must-try. Mar 30 at 15:48