# Bootstrap Multigrid for the Shifted Laplace-Beltrami Eigenvalue Problem

## Abstract

This paper introduces bootstrap two-grid and multigrid finite element approximations to the Laplace-Beltrami (surface Laplacian) eigen-problem on a closed surface. The proposed multigrid method is suitable for recovering eigenvalues having large multiplicity, computing interior eigenvalues, and approximating the shifted indefinite eigen-problem. Convergence analysis is carried out for a simplified two-grid algorithm and numerical experiments are presented to illustrate the basic components and ideas behind the overall bootstrap multigrid approach.

ootstrap multigrid, multigrid eigensolver, Laplace-Beltrami eigenvalue problem, surface finite element method

58C40, 65N25, 65N30, 65N55

## 1 Introduction

We consider developing multigrid methods for the surface finite element (SFEM) approximation to the Laplace-Beltrami eigen-problem

(1.1) |

where denotes the Laplace-Beltrami operator on a 2-dimensional, smooth, orientable, and closed surface , is the eigenvalue to the continuous eigenproblem, and denotes the associated eigenfunction. Letting

(1.2) |

the variational formulation of (1.1) is as follows: Fnd and such that

(1.3) |

where is the Sobolev space defined on :

(1.4) |

equipped with the -norm and -seminorm:

(1.5) |

For a closed surface, it is known that the first eigenvalue of is always 0 with a constant eigenfunction (e.g. see [8] Chapter 1), and the integral of eigenfunction associated with a nonzero eigenvalue on the whole surface is 0. Instead of working on recovering the eigenpairs with this “zero average” compatibility condition, it is convenient to compute approximate eigenfunctions in as in (1.3) and set all eigenfunctions associated with a nonzero eigenvalue perpendicular to the first eigenfunction with 0 eigenvalue, which is the approach that we use in the proposed algorithm.

Throughout this paper, we assume that a surface finite element (SFEM) discretization is used to approximate problem (1.3), an approach first introduced by Dziuk in [14]. A summary of existing works on SFEM is found in [15]. The linear algebra aspects of these problems, e.g., the condition number of their associated stiffness and mass matrices on certain triangulations, are discussed in [23]. Recently in [28] Chen et. al. studied patch recovery techniques for the SFEM.

For the Laplace eigen-problem on the plane, a two-grid eigensolver was proposed by Xu [29], and later improved using a Newton type iteration by Hu and Zhou in [17]. A similar approach was designed for the Maxwell eigen-problem in [30]. These two-grid methods involve a coarse mesh and a fine mesh and the finite element spaces defined on these meshes. In addition, they use a direct solve, e.g., eig in Matlab, to solve the coarse space eigenvalue problem, and then Newton’s method is applied on the fine mesh in order to solve the nonlinear eigen-problem, i.e., solving an appropriately chosen linear source problem. To the author’s best knowledge, there is no known two-grid (or multigrid) eigensolver for the Laplace-Beltrami operator on surfaces in the finite element setting, which is the focus of this paper. The extension we consider involves a suitable geometric projection from coarse spaces to fine spaces defined on a sequence of refined and non-nested meshes and the use of a bootstrap procedure to iteratively enrich the coarse spaces until the desired approximation is computed to sufficient accuracy.

Though two-grid methods do provide improvements when compared to single-grid methods, such as the Arnoldi algorithm in terms of their computational complexity, they have two main drawbacks in practice. First, two-grid methods are generally not optimal as the mesh spacing goes to zero since the coarse eigen-problem needs to be solved to a high accuracy. Second, to resolve eigenpairs corresponding to large eigenvalues in the discrete spectra, the coarse mesh used in the two-grid method must be fine enough. We observe a “loss of spectra” phenomena for the two-grid methods when the coarse mesh is not fine enough. This issue is overcome in the proposed algorithm by using a bootstrap multilevel approach from [4]. We note in addition that, these two-grid eigensolvers require solving a linear source problem on the fine mesh that is indefinite so that using optimal solvers such as multigrid can become problematic. As we show numerically in this paper, it is not necessary to solve this indefinite problem directly and a few sweeps of an iterative solver suffices to obtain an optimal multigrid algorithm. In fact, for certain cases, we show that the shift can be moved to the right hand side using an interpolated coarse approximation to the eigenfunction of interest.

In [6], multilevel analogues of the two-grid solvers noted above are developed. Specifically, the paper develops multilevel approaches for nearly singular elliptic problems and eigen-problems. It should be noted that these methods are able to approximate the components in the eigenspace with small eigenvalues of (1.3) and as presented can’t be used to approximate larger eigenpairs.

In a recent paper by Lin and Xie ([20]), another multilevel approach was developed. The main ingredient in this method is to solve the coarse eigen-problem in an enriched space. This enrichment is achieved by including a single extra function in the coarse space that is obtained by solving a positive definite source problem on a finer mesh. Then, this two-grid correction scheme is used repeatedly to span multiple levels, resembling the bootstrapping procedure developed in [4].

It is known that the Laplace-Beltrami eigenvalue has very high multiplicity for closed surfaces (e.g. see [27] Chapter 3). For example, for the Laplacian-Beltrami operator on the 2-sphere, the number of linearly independent eigenfunctions associated with -th distinct eigenvalue is . Thus, an approach that enriches the coarsest space with a subspace of linearly independent functions is needed for this model problem. The approach that we consider here is the bootstrap multilevel eigensolver proposed in [4]. Therein, a bootstrap algebraic multigrid (BAMG) process is proposed that can be used for computing multiple (smooth and oscillatory) eigenpairs for symmetric and positive definite matrices and, as shown in [5], the approach can also be extended to non-Hermitian (or Hermitian and indefinite) systems. The main component of the BAMG approach is its enrichment of the coarsest space with multiple functions obtained from approximating the source problems on finer meshes with appropriately chosen right hand sides.

We note that the idea to enrich the coarse space in designing eigensolvers goes back to [18] and [21]. The authors also analyze an iterative method for computing the smallest eigenpair under the somewhat restrictive condition that the initial guess of the eigenfunction is sufficiently close to the smallest one, namely that its Rayleigh quotient lies between the smallest and second smallest eigenvalues. In [7], the method from [18] is extended to both two-grid and multigrid methods and an algorithm for computing a given number of the smallest eigenpairs is presented. The paper also presents a convergence theory with less restrictive assumptions on the initial guess.

In this paper, we develop a geometric bootstrap multigrid solver (BMG) for the surface finite element discretization of the shifted Laplace-Beltrami eigen-problem. The approach we propose can be seen as a generalization of the approaches proposed in [18, 21, 7, 20] in that the coarse space is enriched with a subspace, instead of a single eigenfunction and we consider computing interior eigenvalues directly. In addition, we consider solving the shifted and indefinite Helmholtz type eigen-problem. Alternatively, our proposed approach can be viewed as a simplification of the BAMG algorithm in that we use the finite element spaces to explicitly define the components of the multilevel method, including interpolation and restriction operators among different levels, and the enriched coarsest space eigen-problem.

This paper is organized as follows. In Section 2, we provide some preliminary notations and extend the classical a priori estimate of the finite element method for elliptic eigen-problems on the plane to the Laplace-Beltrami eigen-problem on a surface. In Section 3, we introduce the standard two-grid method for Laplace-Beltrami eigen-problems on surfaces mimicking the approaches developed for elliptic eigen-problem introduced in [29] and further developed in [17, 30]. In addition, we prove the convergence of this method for the case of a closed and orientable surface. In Section 4, we introduce a finite element bootstrap multigrid method for the same problem and give details on the approach. Section 5 contains results of numerical experiments for both the two-grid and multigrid methods applied to the model problem on the 2-sphere . Note that by fixing the geometry we are able to study the algorithm in a detailed and systematic way, our future research will focus on developing robust error estimators that will allow us to optimize the proposed algorithm.

## 2 Notation and Preliminary Results

In this section, the finite element approximation, together with its a priori error estimate, to the eigen-problem (1.3) is established. Approximating (1.3) with the finite element method involves two discrete approximations. First, a polygonal surface with a finite set of vertices is generated to approximate the original smooth surface . This discrete surface consists of triangles, i.e., a triangulation is constructed. Second, a finite element discretization is constructed to approximate the continuous eigen-problem on this discrete polygonal surface.

The surface gradient on a 2-dimensional smooth orientable surface that can be embedded into is defined as follows: {definition} For any , define the surface gradient operator as:

(2.1) |

where is a smooth extension of to a 3-dimensional neighborhood of , is the weak gradient operator in , and is the unit normal pointing to the outside of this closed surface at point . The Laplace-Beltrami operator is defined in a distributional sense:

(2.2) |

For a more detailed definition and the technicalities that arise when defining a differential operator on surfaces, we refer to [11, 15, 28].

### 2.1 The eigen-problem on the discrete surface

With being the triangulation mentioned in the beginning of Section 2, let be that polygonal surface approximating the continuous surface , where stands for the “flat” triangular element. is assumed to be quasi-uniform and regular. The mesh size is then defined as the maximum of the diameter of all the triangles: . Furthermore, the set of all vertices is denoted by . For any , it is assumed that , i.e., any vertex in the triangulation lies on the original continuous surface .

Note that the surface gradient on a smooth surface in Definition 2 carries over naturally to a discrete surface : the unit normal is now a constant vector for each point . For ease of notation, the surface gradient on and on will both be denoted by , where the definition should be clear from the context. With these definitions the bilinear forms on is as follows

(2.3) |

The fact that this discrete surface is piecewise linear affine, which is a -surface, implies that the Sobolev space is well-defined on this surface (see [14]).

The continuous variantional formulation for the eigen-problem on the discrete surface is now given by: find and such that

(2.4) |

Using the Poincaré inequality ([11] Lemma 2.2) or the compact embedding of when is a piecewise linear affine manifold ([1] Chapter 2), and the geometric error estimate between (2.2) and (2.3) (e.g. see [15] Section 4), it follows that if the mesh is sufficiently fine (required for the coercivity), then for any

(2.5) | ||||

and |

Throughout this article, and are used for convenience to represent and respectively, where and are two constants independent of the mesh size and eigenvalues. The constants in these inequalities may in certain cases depend on specific eigenvalue(s) and when such dependence exists, it will be stated explicitly.

If , , the coercivity and continuity of (2.5) implies that induces a bounded, compact, and self-adjoint operator, which is exactly the Laplacian-Beltrami operator defined in the distributional sense (2.2). By the Hilbert-Schmidt theory, and the spectrum theory of the Laplacian-Beltrami operator on compact surfaces (every closed surface being compact, see [8]), problem (2.4) is a well-posed self-adjoint eigen-problem. The eigenvalues for problem (2.4) form a discrete sequence, starting from 0, with no accumulation point:

Moreover, the eigenfunctions associated with are orthogonal in the sense that .

### 2.2 Finite element approximation

In this subsection, the surface finite element discretization (2.9) of the eigen-problem (2.4) is established. Here, if the geometric error introduced by the discrete surface is sufficiently small, then the surface finite element approximates the eigen-problem (1.3) on the original smooth surface. In the last part of this subsection, a priori error estimation for the surface finite element eigen-problem using a direct eigensolve is proved, giving Lemma 2.2. Note that, the orders of the approximation errors for the computed eigenpairs given in this lemma are useful in determining the effectiveness of an iterative procedure to obtain approximate eigenpairs, namely, the two-grid or multigrid method needs to compute approximations with the same order of approximation error.

The finite element approximation to problem (2.4) uses piecewise affine linear polynomials which are continuous across the inter-element edges as the test function and trial function spaces denoted by :

(2.8) |

The finite element approximation to problem (2.4) is then: find and such that

(2.9) |

Note that the finite element approximation problem (2.9) serves as a straightforward conforming discretization to (2.4) on the discrete polygonal surface, but not directly to the original eigen-problem (1.3). The connection between the finite element solution on and its continuous counterpart on the surface is established through a bijective lifting operator (see [11]).

It is assumed that there is a bijective mapping between any triangle to a curvilinear triangle on . Then, for any , its lifting to the continuous surface can be defined as follows: for any point , there is a unique point , such that

(2.10) |

where is the signed distantce to at point and is positive when is outside of the closed surface , with . When all the vertices of lie on the continuous surface with -smoothness (), it is known that (see [15] Lemma 4.1).

Conversely, through this lifting bijection, for any , its restriction on the discrete surface can be defined as such that for any

(2.11) |

In [14], the following lemma comparing the -seminorms between the original function lying on the discrete surface and the lifting to the continuous surface is proved. {lemma} If , then for any

(2.12) | |||

Let measure the approximation, under the -seminorm, of the discrete space to the continuous eigenspace on the discrete triangulated surface :

(2.13) |

Using the estimate from the quasi-interpolation of Clément-type as introduced in [10] and extended to surface finite elements in [11], we have . Notice that a standard global Bramble-Hilbert estimate on a flat domain , like , cannot be applied. The reason is that is not well-defined for a polygonal surface , due to the fact that is not continuous where is the co-normal of (See remarks in Section 2 [14]). Here, to be well-defined stands for the Sobolev space containing functions with second weak derivatives -integrable:

(2.14) |

And using the lifting as defined in (2.10), we define the following analogous measure:

(2.15) |

Given these definitions we now present the main result of this subsection in Lemma 2.2. Back in [12], the author establishes the well-known cotangent formulation (for summary and history please refer to [22]) approximating the Laplace-Beltrami operator on a surface that can be embedded in . Here, in our setting the orientability of guarantees that it can be embedded into . The following a priori estimate can be proved for the eigenvalues of the discrete Laplace-Beltrami operator using cotangent formula, namely that the error is bounded by with a factor related to the magnitude of the eigenvalue .

The cotangent formulation of the Laplace-Beltrami operator is known to be equivalent to the linear finite element formulation on the triangulated surface, e.g. [25] points out this relation, and [14] uses it implicitly. Thanks to these bridging results, the a priori estimate from [12] can be used as a guideline to establish the error estimate for the surface finite element approximation to the eigen-problem (2.9).

Combining the a priori estimates from [2, 12, 14, 19], we have the following a priori estimates for both and the corresponding eigenfuntions from problem (2.9).

If the mesh size is small enough and all the vertices of lie on , then for an eigenvalue of problem (1.3) with multiplicity on , there exist ’s that are the eigenvalues of problem (2.9) on , and the following estimate holds

(2.16) |

Moreover, let , where ’s are the eigenfunctions associated with , then for any eigenfunction ,

(2.17) |

We present a brief proof as follows mainly bridging the results of eigen-problem approximations with results for surface finite elements.

Define the elliptic projection of an eigenfunction with respect to the inner product as follows:

(2.18) |

Moreover, . Then by a standard estimate from [2], the identity bridging with in Section 2.3 of [11], and Lemma 2.2 we have

(2.19) |

assuming the mesh size is small enough. Then by the estimate of the Clément-type interpolant in Lemma 2.2 of [11] with the -norm of the metric distortion tensor being 1 plus a higher order term, , the estimate (2.16) follows.

To prove estimate (2.17), assume that the true eigenvalue to be approximated has multiplicity 2. When has bigger multiplicity the same proof follows without essential modifications. Suppose the orthogonal finite element approxmations to the eigenpairs are and from problem (2.9). Let be the projection of onto the discrete eigenspace:

Then a standard estimate can be obtained as follows, with the modification using the lifting operator estimates in Lemma 2.2:

(2.20) | ||||

## 3 A Two-Grid Eigensolver

In this section, a two-grid algorithm to approximate problem (2.9) is presented in Algorithm 1. Algorithm 1 is the extension to the surface case of the two-grid methods introduced in [17, 30, 29].

Assume that a pair of hierarchical meshes, where the fine mesh is uniformly refined from the coarse mesh . The newly created vertices, which are the midpoints of the edges of , are projected onto the continuous surface . Denote the coarse finite element approximation space (2.8) by and the fine space by . The subscript is the solution by a direct eigensolve of eigen-problem (2.9) in . While the superscript in (3.4) of Algorithm 1 is a source problem approximation in the fine space . In the two-grid method, approximates the direct solve solution . The procedures involves a direct eigensolve solution in a coarser space (nonlinear problem), and a source problem approximation using as data in a finer space . This nomenclature, where a subscript stands for direct eigensolve and a superscript stands for source problem approximation, is adopted throughout the remainder of the paper.

An important difference in the surface case as considered in this paper, when being compared with previous works for the Laplacian eigen-problem on the plane, is that the projection operators need extra care. When the mesh is refined, the finite element space on the coarse mesh is not a subspace of the fine mesh. The natural inclusion does not hold even though all vertices from the coarse mesh are defined such that is a subset of when refining.

Consider the geometric projection operator such that

(3.1) |

The definition of is then given by: for any vertex on the fine mesh

(3.2) |

where in the second case, is the projected midpoint between the coarse mesh vertices and (see Figure 1 (c)).

The prolongation operator (or natural inclusion) can then be defined as

To define this operator in matrix notation, suppose the finite element approximation spaces have the following nodal basis set

If , where , then the prolongation operator has the following matrix form:

where is the matrix representation of the prolongation operator. Note that here the geometric projection is implicitly imposed.

For example, by (3.2), , and if , where is a newly created vertex in by projecting the midpoint of and onto the continuous surface (see Figure 1 (c)).

Similarly, the restriction operator, that restricts a finite element function on the fine mesh to the coarse mesh is defined as follows:

Here is opted to be the transpose of the geometric projection defined as a mapping from the fine space to the coarse space. Its matrix representation is

where is the matrix representation of the restriction operator.

(3.3) |

(3.4) |

Given these definitions, the two-grid method approximating the exact solution of problem (2.9) is given by Algorithm 1.

###### Remark \thetheorem (Natural extension to a multilevel method)

When multiple levels of meshes are available ( for with ), Algorithm 1 can be naturally extended to be a multilevel method by being applied in a cascading fashion between two adjacent levels. For example, starting from , when a two-grid eigenpair approximation is obtained, we set . Then step 2 and step 3 in Algorithm 1 are repeated for level 3 to level , where the shift is only added into the Rayleigh quotient on the -th level.

###### Remark \thetheorem (Approximation accuracy of the source problem)

In Algorithm 1, the source problem (3.4) can be approximated by a direct or multilevel method. We note that if a multilevel hierarchy exists and the two-grid method is applied in a cascading fashion, then the accuracy that the source problem needs to be approximated in on the -th level () does not necessarily require accuracy, and often we observe to be sufficient accuracy. This implies that a smoother (relaxation method) can be applied to problem (3.4) in Step 2 of Algorithm 1, with the approximation from the previous level as initial guess, instead of a solve. Below, in the Algorithms 2 and 3, when the term “approximate” is used for the source problems, the user can choose a direct/multigrid solver or smoother. We illustrate this numerically below in Section 5.3.

### 3.1 Convergence analysis

The main ingredient in proving convergence of the two-grid method for the Laplace-Beltrami finite element eigen-problem is to bridge the connection between the geometric projection on the surface with existing two-grid convergence results for the Laplace eigen-problem on the plane. The error introduced by the projection between non-hierarchical finite element spaces that arises in this setting is accounted for in the final estimate we derive.

We use the following lemmas to obtain the convergence estimate of the two-grid method for the Laplace-Beltrami eigen-problem. The first lemma, Lemma 3.1 (e.g. see [2, 17]), gives the stability estimate for the discrete shifted problem.

[Discrete Shift Inf-Sup Condition] If is not an exact eigenvalue to problem (2.9), then there exists a constant such that

The next lemma, Lemma 3.1, is an important identity used to prove the rate of convergence for the approximation of a certain eigenvalue (e.g. see [2]). {lemma} Let be an eigenpair for problem (1.3), then for any

(3.5) |

Aside from in (2.15) which measures of the approximability of the discrete space to the eigenspace , set

(3.6) |

where the operators acts on an functionals to get an Riesz representation:

(3.7) |

Using Lemma 3.1 together with Lemma 3.1, the following lemma is proved for the convergence of the two-grid method for the finite element approximation of Laplace eigen-problem on the plane (see [17]).

Assume that are obtained from the two-grid Algorithm 1 with , then for some eigenfunction , the error estimates are

(3.8) | ||||

This is the known general estimate for the two-grid approximation if the hierarchical coarse and fine finite element spaces are used, i.e. . For finite element eigen-problems on a plane domain , one can normally assume that the eigenspace has certain regularity, e.g., one assumes that is at least twice weakly differentiable. As a result, the infima in the definition of and can be bounded by the canonical finite element interpolation estimate using a Bramble-Hilbert argument:

Then, it is easy to show that the following estimates hold:

(3.9) | ||||

However, in the surface case we are interested in here, Lemma 3.1 cannot be directly applied due to the facts that (a) only is well-defined on a piecewise linear triangulation, but not (see remarks above (2.14)), (b) when refining, the finite element spaces are not hierarchical (see Figure 1).

In the rest of this subsection, a modified two-grid convergence proof is presented in Theorem 3.1 following the idea from [30], and similar bounds are obtained as the standard two-grid convergence results in (3.9).

[Convergence of the two-grid method] Let be an approximation to the eigenvalue of problem (1.3) satisfying the a priori estimate in Lemma 2.2, and assume the eigenpair is obtained from the two-grid method given in Algorithm 1 with . Then there exists an eigenfunction such that the following estimate holds

(3.10) | ||||

and |

Assume the coarse approximation is not an eigenvalue of the discrete eigen-problem (2.9) on the fine mesh. Let an auxiliary solution , where solves problem (3.4) with , then it can be verified that this satisfies:

Now let be a true eigenpair that is obtained from the direct solve for problem (2.9). Then for any test function . Taking the difference of these two equations yields the error equation for the two-grid method as follows: for any

(3.11) |

Applying the discrete inf-sup stability estimate in Lemma 3.1, we have

(3.12) | ||||

By the triangle inequality and the a priori estimate from Lemma 2.2, we have

(3.13) |

The proof of the second estimate is done by restricting the true eigenfunction with its continuous mapping to the discrete surface or , and using the geometric error estimate in Lemma 2.2:

(3.14) | ||||

Now for and , since on each element , the eigenfunction is smooth, then by the interpolation estimate in Lemma 2.2 from [11], we have

(3.15) |

where only depends on the geometry. Using this result and the a priori estimate for the eigenfunction from Lemma 2.2, the following estimate holds:

(3.16) |

Then using the triangle inequality and again using the geometric error estimates for both and ,

(3.17) |

where is assumed to be small enough such that the geometric error, which is , can be omitted comparing with other two terms. Lastly, using the fact that where is the two-grid approximation, we have

(3.18) |

By Lemma 3.1 and assuming the eigenfunction is normalized to have unit norm, and the fact that the geometric error is of higher order, we can get the estimate for the two-grid approximation to the true eigenvalue:

(3.19) | ||||

###### Remark \thetheorem

Theorem 3.1 implies that if the mesh sizes are chosen such that between neighboring levels, then the optimal linear rate of convergence for the eigenfunction in and the quadratic convergence for the eigenvalue follow. In our setting, assuming multiple levels of meshes (obtained by uniformly refining the mesh from previous level) and that we project the vertices onto the surface, holds and the estimate follows. Assume Algorithm 1 is applied in a cascading fashion spanning multiple levels, then the optimal convergence rates of these two algorithms depend on the assumption that the coarse mesh is fine enough, that is, they depend on the assumption that the geometric error is sufficiently small.

## 4 The Bootstrap Multigrid Method

In this section, we propose a finite element BMG (bootstrap multigrid) eigensolver based on the bootstrap algebraic multigrid (BAMG) framework [4]. The essence in the bootstrap approach proposed in this section is to continuously enrich the coarse space with computed eigenfunction approximations coming from finer meshes.

The motivation of the eigensolve in an enriched coarse space is to overcome the drawbacks of the standard two-grid method. Essentially, all the two-grid methods, when using the original geometrically defined coarse space, accelerate the eigensolve on the finest mesh. However, the number of correctly approximated eigenpairs by these two-grid method depend on the dimension of the coarse space (see the numerical example in Section 5.1). With the BMG eigensolver, the entire spectrum of the original Laplace-Beltrami operator can be approximated with same order of accuracy as the direct eigensolve achieves on the finest mesh, assuming certain mesh size relations are satisfied between consecutive levels.

### 4.1 The two-grid bootstrap algorithm

A two-grid bootstrap algorithm is outlined in Algorithm 2 and illustrated briefly in Figure a. The algorithm takes as input the original geometry from the finite element formulation (represented on a coarse mesh), the shift , and the tolerance tol.

In the two-grid bootstrap method, the bilinear forms and are defined as follows. Let be any function such that , where , and . Then, for any test function