Grid independent convergence using multilevel circulant preconditioning: Poisson’s equation

We consider the iterative solution of the discrete Poisson’s equation with Dirichlet boundary conditions. The discrete domain is embedded into an extended domain and the resulting system of linear equations is solved using a fixed point iteration combined with a multilevel circulant preconditioner. Our numerical results show that the rate of convergence is independent of the grid’s step sizes and of the number of spatial dimensions, despite the fact that the iteration operator is not bounded as the grid is refined. The embedding technique and the preconditioner is derived with inspiration from theory of boundary integral equations. The same theory is used to explain the behaviour of the preconditioned iterative method.

For large-scale systems, it is desirable that the solution algorithm is of optimal order, both in terms of computational cost and computer memory. Therefore: 1. The rate of convergence should not slow down when the grid is refined. This is not a trivial task since differential operators are unbounded and the condition number of A consequently increases rapidly as the grid is refined. 2. Both the number of arithmetic operations and the amount of computer memory required for performing one iteration should be proportional to the number of unknowns. This requires a fast and memory lean algorithm for applying M −1 .
For Poisson's equation, geometric and algebraic multigrid methods are of optimal order and among the most effective strategies available [6,16]. However, these techniques contain many parameters and may be difficult to tune. Achieving scalability on parallell architectures may also be a challenge. Other approaches are therefore worth considering.
One alternative line of research originates from Strang [15]. He proposed an algorithm for Toeplitz systems based on the conjugate gradient method with a circulant preconditioner: If M is a circulant approximation of A, the Fast Fourier Transform (FFT) can be used when applying M −1 , and the resulting algorithm is memory lean and has a computational cost of near optimal order. Strang also provided examples where all but a few eigenvalues of M −1 A clustered around 1, yielding fast convergence. This turned out to be true for a wide class of Toeplitz problems [4].
The technique can be extended to multilevel Toeplitz systems, but the rate of convergence is in general degraded [14]. Furthermore, for a discretisation of a differential equation that is ill-posed with periodic boundary conditions, a uni-or multilevel circulant approximation of A is singular, and the technique fails. One possible remedy is to use a semicirculant approximation, as suggested by Holmgren and Otto [9]. Semicirculant preconditioners have circulant structures on all levels but one and an FFT-based algorithm is still available. This turns out to be an effective remedy for first order partial differential equations, and grid independent convergence has been demonstrated for convective flow problems, see e.g. [2].
For diffusive flow problems, the rate of convergence is significantly improved by a semicirculant preconditioner, but the required number of iterations is not bounded as the grid is refined, see e.g. [8]. This holds in particular for elliptic problems [3]. More generally, Manteuffel and Parter [12] have proved that if A and M are discrete elliptic operators, then the condition number of M −1 A is bounded as the grid is refined if and only if A and M has the same boundary conditions. A multilevel circulant or semicirculant approximation M of a discrete elliptic operator A is elliptic if A is elliptic, and circulant and semicirculant preconditioning therefore fits into this framework. Consequently, since circulant approximations imposes periodic boundary conditions, the condition number grows as the grid is refined if A has non-periodic boundary conditions.
In this paper, the pursuit for grid independent convergence continues. We show that Strang's idea can be extended in such a way that grid independent convergence is achieved for a multilevel elliptic problem. Our inspiration comes from Neumann [13] and the research about linear integral equations and potential theory that Neumann contributed to. In 1877, Neumann gave the first rigorous proof for the existence of a solution to Laplace's equation by transforming the elliptic partial differential equation into a boundary integral equation on fixed point form, and by showing that the corresponding fixed point iteration converges. We show that Neumann's fixed point iteration can be viewed as a preconditioned simple iteration for the differential equation and our numerical experiments indicate that the corresponding procedure yields grid independent convergence in the discrete case. We also improve the convergence rate further by using a slightly more advanced iterative method.
The resulting preconditioning procedure involves an embedding of the discrete domain into an extended domain where the preconditioner has a multilevel circulant structure and can be applied using FFT. Applying our preconditioner to a

Problem
We consider Poisson's equation with Dirichlet boundary conditions where is a bounded domain in R d with boundary of class C 2 , and f and g are given continuous functions.
A fundamental solution E of the operator − is given by where γ (d) is the volume of the unit sphere in R d , and |x| is the Euclidean norm of x. Let G(x, y) = E(x − y) be the free space Green's function. Then, (2.1) can be transformed into (see Sect. 4) where n(y) denotes the normal unit vector at y ∈ , directed into the exterior of . The Green's function G(x, y) is singular at x = y, but the integrals on the right hand side of (2.2) exist for points x not lying on the boundary (see e.g. [11]). Also, from the so called jump conditions in potential theory (see e.g. [5]) follows that the integrals can be continuously extended from to¯ = ∪ , if f , g and the normal derivative of u on the boundary is continuous. In this sense, (2.2) holds on¯ . Let F(u) denote the right hand side of (2.2), continuously extended to¯ . Then, given some initial guess u (0) ∈ C(¯ ), the operator F can be used to form a sequence of approximations converges uniformly to a unique solution of (2.2) since the normal derivative of u (k) on the boundary converges uniformly to a unique solution (see Sect. 3). Two tasks are of interest to us in this paper. First, we demonstrate that it may be possible to improve the rate of convergence by using an iterative method that is slightly more advanced than the fixed point iteration (2.3). This is the topic of Sect. 3. The analysis in Sect. 3 guides us when we construct an iterative method for the discrete problem in Sect. 5.
Secondly, we show that it is possible to frame the transformation of (2.1) into (2.2) as a preconditioning operation that can be performed using Fourier technique. This is the topic of Sect. 4. This framing serves as a blueprint for the derivation of an FFT-based algorithm for the discrete problem in Sect. 6.

Iterative method
Given some smooth initial guess u (0) , we consider an iterative method of the form u (k+1) = u (k) − α k+1 r (u (k) ), k = 0, 1, 2, . . . , (3.1) where α k are some chosen constants, and r (u) = u − F(u) is the residual. If α k = 1 for all k, then (3.1) is the fixed point iteration (2.3). Let u * denote the exact solution to (2.2). Then, from (2.2) and (3.1) it follows that the error e (k) = u (k) − u * satisfies for k = 0, 1, 2, . . . , and x ∈¯ . From differentiation of (3.2), it follows that the normal derivative of the error satisfies for k = 0, 1, 2, . . ., x ∈ , and some small and positive parameter h. Letting h → 0 + uniformly on and using the jump conditions in potential theory (see e.g. [5]) yields for k = 0, 1, 2, . . . , and x ∈ . Here, the integral on the right hand side of (3.3) is singular, but exists as an improper integral if the normal derivative of the error on the boundary is continuous. For ease of notation, we introduce Then, (3.3) can be written for k = 0, 1, 2, . . . , and x ∈ . If (3.5) is a convergent iteration and the normal derivative of the error on the boundary tends to the zero function with a certain rate, then (3.2) can be used to estimate the rate of convergence of the error in the entire domain. This suggests that a convergence analysis of (3.1) should start with the analysis of (3.5), and, in particular, with the analysis of the integral operator K given by (3.4). Denote K 's spectrum by σ (K ). Then (see e.g. [11]): -K is compact on C( ), and, as for any other compact operator, λ = 0 belongs to σ (K ) and σ (K ) \ {0} consists of at most a countable number of eigenvalues. -The nullspace of I + K , where I is the identity operator, has dimension one and is spanned by a function that is not identical to zero, so λ = −1 belongs to σ (K ). -The eigenvalues of K are real and σ (K ) ⊂ [−1, 1).
In general, fixed point iteration requires an infinite number of iterations. However, as the following theorem indicates, if K has a finite number of eigenvalues and the eigenfunctions form an orthonormal basis in L 2 ( ), it is possible to choose the coefficients α k in such a way that a solution can be found in a finite number of steps. Theorem 3.1 If K has a finite number of eigenvalues, λ 1 , . . . , λ s , and the eigenfunctions μ m (x) of K form an orthonormal basis in L 2 ( ), then, letting for some constants β m , it follows from (3.5) that for some constants c m . For terms where K μ m = λ 1 μ m , it follows that and, consequently, that ψ (1) is orthogonal to all eigenfunctions corresponding to the eigenvalue λ 1 . By repeating the argument k times, k < s, it follows that ψ (k+1) is orthogonal to all eigenfunctions corresponding to eigenvalues λ 1 , . . . , λ k+1 . From the completeness of the orthonormal set of eigenfunctions then follows that ψ (s) = 0, and from (3.2) that e (s+1) = 0, so u (s+1) = u * .
The result in Theorem 3.1 is not surprising, since a discrete counterpart of (3.5) is a Krylov subspace method. For Krylov subspace methods, the rate of convergence is determined by the eigenvalues of the coefficient matrix, as long as the coefficient matrix has a complete set of orthonormal eigenvectors, see e.g. [7]. Consequently, results simliar to the one in Theorem 3.1 can be derived for Krylov subspace methods. For example, one can use error estimates for the GMRES method to show that GMRES converges in s + 1 iterations if there are only s distinct eigenvalues. A discrete counterpart of the following example is examined in Sect. 7.
Example 3.1 Let be the unit disc |x| < 1 in R 2 and introduce polar coordinates x 1 = r cos θ, x 2 = r sin θ, y 1 = ρ cos ϕ, and y 2 = ρ sin ϕ. Then, it follows from a simple calculation that and, consequently, that Here, σ (K ) consists of two eigenvalues (see [1]): λ = −1 is an eigenvalue with corresponding eigenfunction 1 and λ = 0 is an eigenvalue with corresponding eigenfunctions cos mθ and sin mθ for m ∈ N .
In the case of fixed point iteration, i.e. α k = 1 for all k, relation (3.5) reads and σ (I /2 + K /2) = {0, 1/2}. In fact, Thus, with the exception for the first iteration, the normal derivative of the error on the boundary is halved in each iteration. Now, let α 1 = 1, α 2 = 2, and α 3 = 1, as suggested by Theorem 3.1. Then, it follows from (3.5) that ψ (1) However, as illustrated by the following example, the unit disc is a degenerate case. In general, the operator K has an infinite number of eigenvalues.

Example 3.2
Introduce the elliptic coordinates x 1 = cosh r cos θ , x 2 = sinh r sin θ , y 1 = cosh ρ cos ϕ, and y 2 = sinh ρ sin ϕ, and let be the ellipse r < 1 in R 2 . Then, a not so simple calculation [1] shows that −1, −e −2m , and e −2m are eigenvalues of K with corresponding eigenfunctions 1, cos mθ , and sin mθ for m ∈ N. Here, λ = 0 belongs to the spectrum of K , but as a limit point as m tends to infinity and consequently as a member of the continuous spectrum rather than being an eigenvalue.

Transformation
Equation (2.1) can be transformed into (2.2) by the use of Green's second identity: Using (− G)u dV = u for x ∈ , and substituting u with g on the boundary gives the result. In this section, we frame this transformation as a preconditioning operation that can be performed using Fourier technique. The framing serves as a blueprint for the derivation of an FFT-based algorithm for the discrete problem in Sect. 5. Let be the characteristic function of the domain and define the distributions for distributions v and continuous functions u, as linear functionals acting on infinitely differentiable test functions w with compact support. It then follows from Green's second identity that for every test function w and an arbitrary continuous function u. If u satisfies (2.1), Equation (4.1) is our blue print for the system of linear equations Au = b in the discrete case. When using this formulation, rather than (2.1), the preconditioning procedure is simply a convolution with a fundamental solution: Define the convolution u * v, where u is a function and v is a distribution with compact support, from for functions u and w. Convolving both sides of (4.1) with a fundamental solution E yields Since E * (− v) = v for every distribution v with compact support (see e.g. [10]), it can be shown that (4.2) holds for every test function w if and only if which is the same as (2.2) when x ∈ . Equation (4.2) is our blue print for the preconditioned system of linear equations M −1 Au = M −1 b in the discrete case. Applying the preconditioner M −1 is then simply a convolution with a discrete fundamental solution over a larger, but bounded domain where we impose periodic boundary conditions. Thus, viewed as a matrix, the preconditioner has a multilevel circulant structure and can be applied using FFT.
Let i = (i 1 , . . . , i d ) ∈ Z d be a multi-index and introduce a uniform grid {x i } i∈Z d in R d with step size h k in dimension k, k = 1, . . . , d. Also, letê k denote unit vectors in Z d , and define the discrete Laplace operator h for grid functions u according to

Consider the discrete Poisson's equation with Dirichlet boundary conditions
where the domain is a finite subset of Z d , is the set of points outside of that have at least one neighbour in , and f and g are given grid functions. Furthermore, let + be the set of points outside of¯ = ∪ that have at least one neighbour in¯ , be the characteristic function of the set . Applying − h to χ¯ u yields Comparing (5.2) with (4.1) suggests a splitting R = R 1 + R 2 , where R 1 u and R 2 u are discrete counterparts of ∂u ∂n δ and ∂(uδ ) ∂n .
Such a splitting can be done in several different ways. In this paper, we let If u is a solution to (5.1), then (5.2) turns into which is a discrete counterpart of (4.1). Let E be a fundamental solution of the operator − h (see Sect. 6) and define the convolution of two grid functions u and v according to Convolving both sides of (5.3) with E yields a discrete counterpart of (4.2): Since E * (− v) = v for every grid function v with finite support, it follows that Let F h (u) denote the right hand side of (5.4) and introduce the discrete residual Also, let u * denote the exact solution to (5.5). Then, given some initial guess u (0) , the iteration can be used to form a sequence of approximations u (k) , k = 1, 2, . . ., where the error for k = 0, 1, 2, . . ., and i ∈ Z d . Applying R 1 to both sides of (5.7) and restricting the domain from Z d to yields for k = 0, 1, 2, . . ., and i ∈ , where For ease of notation, introduce for k = 0, 1, 2, . . ., and i ∈ , and for i ∈ . Then, (5.8) can be written for k = 0, 1, 2, . . ., and i ∈ , which is a discrete counterpart of (3.5). This suggests that the convergence properties of (5.6) depends on the convergence properties of (5.10), and, in particular, on the properties of the operator K h , given by (5.9). The spectrum σ (K h ) and the convergence of (5.6) are investigated numerically in Section 7. The results show that despite the fact that F h is not bounded as the grid is refined, u (k) converges to the solution of (5.3) with a rate that is independent of h.

Preconditioning
Let B¯ + ⊂ Z d be a box of size N 1 × · · · × N d such that¯ + = ∪ ∪ + ⊂ B¯ + . Applying the preconditioner requires the computation of where E is a fundamental solution of the operator − h and u is a grid function with supp(u) ⊂¯ + . The sum in (6.1) can be written as a matrix-vector multiplication with a matrix T , where T has a multilevel Toeplitz structure. In this section, we rewrite the sum in (6.1) into a sum that can be written as a matrix-vector multiplication with a matrix C, where C has a multilevel circulant structure, and where the matrix-vector multiplication therefore can be computed using FFT. We also derive an expression for a fundamental solution E.
Now, let u c be a periodic grid function with period M = 2N such that u c i = u i for i ∈ B E . Then, since E c i− j = E i− j when i, j ∈ B¯ + , u c j = u j when j ∈ B E , and supp(u) ⊂ B¯ + . Hence, (6.1) can be computed by computing (E c * c u c ) i for all i ∈ B E and then sampling the result for i ∈ B¯ + . Also, since where u c denotes the discrete Fourier transform of u c , it follows that E c * c u c can be computed by transforming u c , multiplying with E c and taking the inverse transform of the result. Here, the discrete Fourier transform u c of an M-periodic grid function u c is given by whereî is the imaginary unit, and is the inverse discrete Fourier transform. Both the transform and its inverse can be computed by the FFT-algorithm in O(M 1 · · · · · M d · log(M 1 · · · · · M d )) arithmetic operations.
To derive an expression for E c , recall that We require (6.2) to hold also for i k = −N k + 1 and i k = N k − 1, k = 1, . . . , d, but allow for a more general right hand side when i k = N k , k = 1, . . . , d.
More specifically, we look for a grid function E c that satisfies 3) for some constants a k , k = 1, . . . , d. Applying the discrete Fourier transform to both sides of (6.3) and using the orthogonality of the Fourier basis yields a transformed equation 4) where the symbol P j is given by Note that P j = 0 if and only if j = (m 1 M 1 , . . . , m d M d ) for some m ∈ Z d , and, therefore, that (6.4) has a solution only if the right hand side of (6.4) is zero for such a j: Letting

Numerical experiments
In this section, we assume that N k is even for all k, and use a grid {x i } i∈Z d with and h k = 2/(N k − 4) for k = 1, . . . , d. We consider four different domains. First, let This d-dimensional discrete ball is a d-dimensional discrete counterpart of the unit disc in Example 3.1. The second domain is given by where | · | ∞ denotes the maximum norm. This domain is a discrete d-dimensional cube, and, consequently, a discrete counterpart of a domain with edges and corners. The third and fourth domain, B L and C L , are L-shaped modifications of B and C , where points i satisfying (x i ) k > 0 for all k have been removed. Figure 1 shows the four sets and their boundary sets for the case d = 2 and N 1 = N 2 = 16.
Domains with non-smooth boundaries are not covered by the standard results for boundary value problems of potential theory, since the integral operator K for such a boundary is not compact [11]. Fixed point theory can still be used to prove the existence of a unique solution if the norm of the iteration operator (I + K )/2 is less than 1. However, this criterion is only a sufficient condition for convergence and it may be that a fixed point iteration converges even though the norm is larger than or equal to 1. Hence, the theory of boundary integral equations becomes less helpful for understanding the behaviour of our preconditioner on domains with non-smooth boundaries. Still, since structured grids are often used on rectangular patches, nonsmooth cases are of interest.
We first study the spectrum of the summation operator K h given by (5.9). A grid refinement study for the case d = 2 is presented in Figs. 2 and 3. The spectral properties of K h thus resembles the spectral properties of K when = B . The introduction of corners perturbs the spectrum of K h , but those perturbations seem to be relatively small. Table 1 shows the results of a grid refinement study of the maximum norm of the iteration operator (I + K h )/2 for the case d = 2. These results indicate that the maximum norm cannot be used to prove convergence of fixed point iteration, since the norm of the iteration operator is larger than 1 already for small grid sizes.
Next, we solve the discrete Poisson's equation (5.1) on the four domains using two different versions of the iterative method (5.6). The first version is fixed point iteration, i.e. α k = 1 for all k. This method is denoted M1. The second version is denoted M2 and uses α k = 1 for k odd and α k = 2 for k even. This is the method introduced in Example 3.1 in restarted form. We construct the problem data f and g from a known solution u * i = exp(−2(x i ) 2 1 − · · · − 2(x i ) 2 d ) of the difference equation, use u (0) = 0 as initial guess, and iterate until the maximum norm of the error is decreased by a factor of 10 6 .  Tables 2 and 3 show the number of iterations required for convergence for the case d = 2 and d = 3 respectively. The convergence rate seems to be grid independent for all four domains.  Table 4 shows the number of iterations required for convergence for different number of dimensions d. The convergence rate seems to be independent of the number of dimensions for all four domains. Table 1 The maximum norm of (I + K h )/2 when d = 2

Summary and conclusions
In this paper, the iterative solution of the discrete Poisson's equation in d dimensions is considered. The discrete domain is embedded into an extended domain and the resulting system of linear equations is solved using fixed point iteration combined with a multilevel circulant preconditioner. The main objective of the paper is to show that Strang's idea of using circulant preconditioning for Toeplitz systems [15] can be successfully extended to the multilevel case for Poisson's equation.
Our inspiration comes from Neumann [13] and the research about linear integral equations and potential theory that Neumann contributed to. He transformed the elliptic partial differential equation into a boundary integral equation on fixed point form, u = F(u), and showed that the corresponding fixed point iteration, u (k+1) = F(u (k) ), converges. The operator F is not bounded, but the normal derivative of F(u) on the boundary is bounded if the normal derivative of u on the boundary is bounded, and convergence of the normal derivative of u (k) on the boundary turns out to be a sufficient condition for the convergence of u (k) .
We show that Neumann's findings holds also for the discrete Poisson's equation: The elliptic partial difference equation can be transformed into a boundary summation equation on fixed point form, u = F h (u), and our numerical results show that the corresponding fixed point iteration, u (k+1) = F h (u (k) ), converges with a rate that is independent of the grid's step sizes h, despite the fact that the operator F h is not bounded as the grid is refined.
We also show that Neumann's transformation of Poisson's equation can be viewed as a convolution with a fundamental solution to both sides of a modified partial differential equation on an extended domain, and that a corresponding transformation of the discrete Poisson's equation can be viewed as a multilevel circulant preconditioning of both sides of a modified partial difference equation on an extended domain.
Our preconditioned iterative solver is consequently of near optimal order: Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.