High-order algorithms for solving eigenproblems over discrete surfaces

# High-order algorithms for solving eigenproblems over discrete surfaces

Sheng-Gwo Chen Mei-Hsiu Chi Jyh-Yang Wu Department of Applied Mathematics, National Chiayi University, Chia-Yi 600, Taiwan. Department of Mathematics, National Chung Cheng University, Chia-Yi 621, Taiwan.
###### Abstract

The eigenvalue problem of the Laplace-Beltrami operators on curved surfaces plays an essential role in the convergence analysis of the numerical simulations of some important geometric partial differential equations which involve this operator. In this note we shall combine the local tangential lifting (LTL) method with the configuration equation to develop a new effective and convergent algorithm to solve the eigenvalue problems of the Laplace-Beltrami operators acting on functions over discrete surfaces. The convergence rates of our algorithms of discrete Laplace-Beltrami operators over surfaces is , , where represents the size of the mesh of discretization of the surface. The problem of high-order accuracies will also be discussed and used to compute geometric invariants of the underlying surfaces. Some convergence tests and eigenvalue computations on the sphere, tori and a dumbbell are presented.

###### keywords:
Eigenproblem, Local tangential lifting method, Configuration equation, Discrete Laplace-Beltrami operator.
journal:

## 1 Introduction

Let be a smooth regular surface in the 3D space. The Laplace-Beltrami (LB) operator is a natural generalization of the classical Laplacian from the Euclidean space to a curved space. To understand the LB operators on curved surfaces, it is natural to investigate their associated eigenvalue problems:

 ΔΣϕ=λϕ. (1)

Or, more generally,

 ∇Σ⋅(h∇Σϕ)=λϕ (2)

where is real function on .

The eigenvalue problem of the LB operator plays important roles not just in the study of geometric properties of curved spaces, but also in many applications in the fields of physics, engineering and computer science. The LB operator has recently many applications in a variety of different areas, such as surface processingClarenz (); Sapiro (), signal processingRomeny () and geometric partial differential equationsDesbrun ().

Since the objective underlying surfaces to be considered are usually represented as discrete meshes in these applications, it is useful in practice to discretize the LB operators and solving the eigenproblems over discrete surfaces. There are many approaches for estimating Laplace-Beltrami operator and solving the Laplace-Beltrami eigenproblems Shi1 (); Shi2 (). In 2011, MacdonaldMacdonald () proposed an elegant method to solve the Laplace-Bletrami eigenproblems for Equations (1) and (2) by the closest point methodRuuth (). In this paper we shall describe simple and effective methods with high-order accuracies to define the discrete LB operator on functions on a triangular mesh.

In 2012, Ray et al.Ray () used the method of least square to obtain high-order approximations of derivatives and integrations. In this paper, we shall use ideas developed in Chen, Chi and WuChen3 (); Chen4 () where we try to estimate the discrete partial derivatives of functions on 2D scattered data points. Indeed, the ideas that we shall use to develop our algorithms are divided into two main steps: first we lift the 1-neighborhood points to the approximating tangent space and obtain a local tangential polygon. Second, we use some geometric idea to lift functions to the tangent space. We call this a local tangential lifting (LTL) methodChen5 (); Wu (). Then we present a new algorithm, the configuration method, to compute their Laplacians in the 2D tangent space. This means that the LTL process allows us to reduce the 2D curved surface problem to the 2D Euclidean problem.

In other words, we shall combine the local tangential lifting (LTL) method with the configuration equation to develop a new effective and convergent algorithm to solve the eigenpair problems of the Laplace-Beltrami operators acting on functions over curved surfaces. We shall also present a mathematical proof of the convergence of our algorithm. Our algorithm is not only conceptually simple, but also easy to implement. Indeed, the convergence rate of our new algorithms of discrete Laplace-Beltrami operators over surfaces is , , where represents the size of the mesh of discretization of the surface. In section 2, we introduce the gradient of a function, the divergence of a vector field and the Laplace-Beltrami operator on regular surfaces. Our -LTL configuration method is discussed in section 3. In section 4, we discuss how to improve the methods to have high-order accuracies. We also give some numerical simulations to support these results in section 5.

## 2 The gradient, divergence, and Laplace-Beltrami operator

In order to describe the gradient, divergence and the LB operator on functions or vector fields in a regular surface in the 3D Euclidean space , we consider a parameterization at a point , where is an open subset of the 2D Euclidean space . We can choose, at each point of , a unit normal vector . The map is the local Gauss map from an open subset of the regular surface to the unit sphere in the 3D Euclidean space . Denote the tangent space of at the point by . The tangent space is a linear space spanned by where are coordinates for .

The gradient of a smooth function on can be computed from

 ∇Σg=guG−gvFEG−F2xu+gvE−guFEG−F2xv (3)

where , and are the coefficients of the first fundamental form and

 gu=∂g(x(u,v))∂u and gv=∂g(x(u,v))∂v. (4)

See do CarmoDocarmo (); Docarmo2 () for the details.

Let be a local vector field on . The divergence, , of is defined as a function given by the trace of the linear mapping for . A direct computation gives

 ∇Σ⋅X=1√EG−F2(∂∂u(A√EG−F2)+∂∂v(B√EG−F2)). (5)

The LB operator acting on the function is defined by

 ΔΣg=∇Σ⋅∇Σg. (6)

for all smooth function on . A direct computation yields the following local representation for the LB operator on a smooth function :

 ΔΣg=1√EG−F2[∂∂u(G√EG−F2∂g∂u)−∂∂u(F√EG−F2∂g∂v)]+1√EG−F2[∂∂v(E√EG−F2∂g∂v)−∂∂v(F√EG−F2∂g∂u)]. (7)

## 3 An O(r)-LTL configuration method

In this section, we shall introduce a new algorithm to solve the eigenpair problems, Equations (1) and (2), by the LTL configuration method.

Consider a triangular surface mesh , where is the list of vertices and is the list of triangles.

To describe the local tangential lifting (LTL) method, we introduce the approximating tangent plane and the local tangential polygon at the vertex of as follows:

1. The normal vector at the vertex in is given by

 NA(v)=∑T∈T(v)ωTNT∥∑T∈T(v)ωTNT∥ (8)

where is the set of triangles that contain the vertex , is the unit normal to a triangle face with for all and the centroid weight is given in Chen (); Chen2 () by

 ωT=1∥GT−v∥2∑~T∈T(v)1∥G~T−v∥2 (9)

where is the centroid of the triangle face determined by

 GT=v+vi+vj3. (10)

Note that the letter in the notation stands for the word ”Approximation”.

2. The approximating tangent plane of at is now determined by .

3. The local tangential polygon of in is formed by the vertices which is the lifting vertex of adjacent to in :

 ¯vi=(vi−v)−NA(v) (11)

as in figure 1.

4. We can choose an orthonormal basis for the tangent plane of at and obtain an orthonormal coordinates for vectors by . We set with respect to the orthonormal basis .

Now we explain how to lift locally a function defined on to the local tangential polygon . Consider a function on . We will lift locally the function to a function of two variables , denoted by , on the vertices in by simply setting

 ¯ϕ(xi,yi)=ϕ(vi) (12)

and where is the origin of . Then one can extend the function to be a piecewise linear function on the whole polygon in a natural and obvious way.

Next, we introduce the configuration matrices for and in Equations (1) and (2).

### 3.1 The Configuration matrix for ΔΣϕ

According to the LTL method, the differential quantities at a point on curved surfaces correspond to planar differential quantities in . Hence, we only need to estimate the Laplace-Beltrami operator on planar triangular meshes. Given a function on an open domain with the origin , Taylor’s expansion for two variables and gives

 f(x,y)=f(0,0)+xfx(0,0)+yfy(0,0)+x22fxx(0,0)+xyfxy(0,0)+y22fyy(0,0)+O(r3) (13)

when is small.

Consider a family of neighboring points , , of the origin . Take some constants , , with . Then one has

 n∑j=1αj(f(xj,yj)−f(0,0))=(n∑j=1αjxj)fx(0,0)+(n∑j=1αjyj)fy(0,0)\par+12(n∑j=1αjx2j)fxx(0,0)+(n∑j=1αjxjyj)fxy(0,0)\par+12(n∑j=1αjy2j)fyy(0,0)+O(r3), (14)

where . To estimate the Laplacian, , at , we choose the constants , with , so that they satisfy the following equations:

and

 n∑j=1αjx2j=n∑j=1αjy2j

or equivalently

 n∑j=1αj(x2j−y2j)=0

One can rewrite these equations in a matrix form with the condition and obtain the following equation:

 (15)

The solutions of this equation allow us to obtain a formula for the Laplacian :

 Δf(0,0)=fxx(0,0)+fyy(0,0)=2n∑j=1αj(f(xj,yj)−f(0,0))n∑j=1αjx2j+O(r) (16)
###### Remark 1.
1. Equation (15) is called the configuration equation of the Laplace-Beltrami operator. We call the matrix in Equation (15) the configuration matrix of Laplace-Beltrami operator and the solution in Equation (15) the configuration coefficients of the Laplace-Beltrami operator, respectively.

2. For the reason of symmetry, Equation (16) also gives

 Δf(0,0)=4n∑j=1αj(f(xj,yj)−f(0,0))n∑j=1αj(x2j+y2j)+O(r) (17)

since we have .

3. For simplicity, the scalar in Equation (15) can be chosen to be .

4. It is worth to point out that Equation (17) is an generalization of the well-known 5-point Laplacian formula. In the 5-point Laplacian case, we have the origin along with 4 neighboring points , ( and for sufficiently small positive number . One can find a solution for , in this case.

Now, let be a regular surface with a triangular surface mesh of . The set is the list of vertices on , is the list of triangles on and is the number of 1-neighbors of on . Suppose that is a function on . For each vertex on , is the tangential polygon of the neighbors of with coordinates and is the lifting function of . By the configuration equation of the Laplacian in Equation (15), the Laplace-Beltrami operator, on , is defined by

 ΔSϕ(vi)=2∑N(i)j=1αi,j(ϕ(xi,j,yi,j)−ϕ(0,0))∑N(i)j=1(αi,jx2i,j). (18)

Then one can prove

###### Theorem 1.

Given a smooth function on a closed regular surface and a triangular surface mesh with mesh size , one has

 ΔΣϕ(v)=ΔSϕ(v)+O(r) (19)

where the discrete LB operator is given in Equation (18).

We will prove Theorem 1 by the following Lemmas. Indeed, Theorems 2 and 3 in this section can also be proved in a similar way.

Given a smooth function on a regular surface , we can lift via the exponential map locally to obtain a smooth function defined on by setting

 ^h(w)=h(expp(w)) (20)

for . Fix an orthonormal basis for the tangent space . This gives us a coordinate system on . Namely, for we have for two constants and . Without loss of ambiguity, we can identify the vector with the vector with respect to the orthonormal basis . In this way, the function can also give us a smooth function of two variables and by defining

 ~h(x,y)=^h(w) (21)

for . Using these notations, we will prove

###### Lemma 1.

One has

 ΔΣh(p)=Δ^h(0)=Δ~h(0,0). (22)
###### Proof.

It is well-known that the LB operator acting on a smooth function at a point can be computed from the second derivatives of along any two perpendicular geodesics with unit speed. See do Carmo Docarmo2 () for details. Indeed, we consider the following two perpendicular geodesics with unit speed in by using the orthonormal vectors :

 ci(t)=expp(t~ei),  i=1,2 (23)

with and . One has

 ΔΣh(p)=d2dt2h(c1(t))|t=0+d2dt2h(c2(t))|t=0=d2dt2^h(t~e1)|t=0+d2dt2^h(t~e2)|t=0=Δ^h(→0)=∂2~h∂x2(0,0)+∂2~h∂y2(0,0)=Δ~h(0,0). (24)

Next we consider a triangular surface mesh for the regular surface , where is the list of vertices and is the list of triangles and the mesh size is less than . Fix a vertex in . For each face containing , we have

 NΣ(v)=NT+O(r) (25)

where is the unit normal vector of the true tangent plane of at and is the unit normal vector of the face . Since the approximating normal vector , defined in section 3 is a weighted sum of these neighboring face normals , we have

###### Lemma 2.

One has

 NΣ(v)=NAv+O(r). (26)

Due to this lemma, the orthonormal basis for the tangent plane will give us an orthonormal basis for the approximating tangent space by the Gram-Schmidt process in linear algebra:

 e1=~e1−<~e1,NA(v)>NA(v)∥~e1−<~e1,NA(v)>NA(v)∥,

and

 e2=~e2−<~e2,NA(v)>NA(v)−<~e2,e1>e1∥~e2−<~e2,NA(v)>NA(v)−<~e2,e1>e1∥.

Logically speaking, one can first choose an orthonormal basis for the approximating tangent space and then apply the Gram-Schmidt process to obtain an orthonormal basis for the tangent plane . In either way, we always have by Lemma 2 the following relations.

###### Lemma 3.

One has

 ~ei=ei+O(r),  i=1,2. (27)

Consider a neighboring vertex of in . For small enough, we can use the inverse of the exponential map to lift the vertex up to the tangent plane and obtain

 ~vi=exp−1v(vi)∈TΣ(v)

and

 ~vi=~xi~e1+~yi~e2

for some constants. As discussed in section 3, we can also lift the vertex up to the approximating tangent plane and get

 ¯vi=(vi−v)−NA(v)

and for some constants . Then Lemmas 2 and 3 yield

###### Lemma 4.

One has

 {~xi=xi+O(r2)~yi=yi+O(r2). (28)

Using these relations, one can solve the configuration equation (15) for and respectively and obtain their corresponding solutions and with the relation

 ~αi=αi+O(r). (29)

Note that the lifting function is a smooth function of two variables and . Equation (17) now gives an approximation of the Laplacian :

 Δ~h(0,0)=4n∑i=1~αi(~h(xi,yi)−~h(0,0))n∑i=1~αi(~xi2+~yi2)+O(r). (30)

The relations (21), (26) and (27) imply

 Δ~h(0,0)=4n∑i=1αi(h(vi)−h(v))n∑i=1αi(x2i+y2i)+O(r). (31)

This along with Lemma 1 proves Theorem 1.

For each vertex , we have

 ΔΣϕ(vi)=2∑N(i)j=1(αi,j(ϕ(vi,j)−ϕ(vi)))∑N(i)j=1(αi,jx2i,j)+O(r). (32)

Denote . Since , Equation (32) can be rewritten as

 ΔΣϕ(vi)=ωi⎡⎣⎛⎝N(i)∑j=1(αi,jϕ(vi,j)⎞⎠−ϕ(vi)⎤⎦+O(r). (33)

Furthermore the vector can be easily extended to a vector by

 ai,k=⎧⎪⎨⎪⎩αi,j if there exists j∈{1,2,⋯,N(i)} such that vk=vi,j,−1 if i=k,0 otherwise. (34)

Obviously,

 ΔΣϕ(vi)=ωi[(ai,1,ai,2,⋯,ai,nV)(ϕ(v1),ϕ(v2),⋯,ϕ(vnV))T]+O(r). (35)
###### Remark 2.

is the set of indices of all vertices in and denotes the ith vertex in . For each , is the set of indices of one-neighbors of and denotes the jth one-neighbor of in . Obviously, every one-neighbor of is corresponding to a unique vertex in while and may be not equal.

This implies that

 (36)

where and are two matrices. For simplicity, we rewrite Equation (36) as

 ΔΣϕ(V)=(WA)ϕ(V)+O(r). (37)

Hence, we have an eigenvalue approximation result by the method discussed in Strange1 ().

###### Theorem 2.

Let be a closed regular surface, be a triangular mesh of with mesh width . If is the ith eigenvalue of the Laplace-Beltrami operator on and is the ith eigenvalue of the matrix ., then we have, for sufficiently small ,

 λi=¯λi+O(r).
###### Remark 3.

The matrix equation

 (38)

is called the configuration equation of the Laplace-Beltrami operator at on . The constants are called the configuration coefficients of Laplace-Beltrami operator at on . The matrix defined in Equation (34) is called the configuration matrix of the Laplace-Beltrami operator on .

### 3.2 The Configuration matrix for ∇Σ⋅(h∇Σϕ)

Let be a bounded smooth function defined on a regular surface . We introduce the configuration matrix of the quantity by a similar method as in subsection 3.1. First, let us consider two smooth functions and defined on an open domain in with the original point . Since , we need to estimate the quantity in . Taylor expansions of and are given by

 f(x,y)−f(0,0)=xfx(0,0)+yfy(0,0)+x22fxx(0,0)+xyfxy(0,0)+y22fyy(0,0)+O(r3) (39)

and

 ⎧⎪⎨⎪⎩f(x,y)−f(0,0)=xfx(0,0)+yfy(0,0)+O(r2)g(x,y)−g(0,0)=xgx(0,0)+ygy(0,0)+O(r2), (40)

when is small. From Equation (40), one has

 (f(x,y)−f(0,0))(g(x,y)−g(0,0))=(xfx(0,0)+yfy(0,0))(xgx(0,0)+ygy(0,0))+O(r3)=x2fx(0,0)gx(0,0)+xy(fx(0,0)gy(0,0)+fy(0,0)gx(0,0))  +y2fy(0,0)gy(0,0)+O(r3). (41)

These imply

 (f(x,y)−f(0,0))(g(x,y)−g(0,0))+2g(0,0)(f(x,y)−f(0,0))=2xg(0,0)fx(0,0)+2yg(0,0)fy(0,0)+xy(2g(0,0)fxy(0,0)+fx(0,0)gy(0,0)  +fy(0,0)gx(0,0))+x2(g(0,0)fxx(0,0)+fx(0,0)gx(0,0))  +y2(g(0,0)fyy(0,0)+fy(0,0)gy(0,0))+O(r3). (42)

Consider a family of neighboring points ,¸ ¶, of the origin . Take some constants ¿ with . We have

 (43)

To compute at , we choose the constants , , so that they satisfy the following equations:

 (44)

The solutions of this equation gives a formula for at :

 (∇g⋅∇f+gΔf)(0,0)=gx(0,0)fx(0,0)+gy(0,0)fy(0,0)+g(0,0)[fxx(0,0)+fyy(0,0)]=1n∑j=1βjx2jn∑j=1[(f(xj,yj)−f(0,0))(g(xj,yj)−g(0,0))+2g(0,0)(f(xj,yj)−f(0,0))]+O(r)=n∑j=1(f(xj,yj)−f(0,0))(g(xj,yj)+g(0,0))n∑j=1βjx2j+O(r) (45)

Using the notations of subsection 3.1, the quantity at on is given by

 ∇S(h∇Sϕ)(vi)=N(i)∑j=1(βi,j(f(xi,yi)−f(0,0))(g(xi,yi)+g(0,0)))N(i)∑j=1βi,jx2i,j. (46)

Then one can prove

###### Theorem 3.

Given two smooth functions on a closed regular surface with a triangular surface mesh , one has

 ∇Σ(h∇Σ)(vi)=∇S(h∇Sϕ)(vi)+O(r) (47)

where the quantity is given in Equation (46) and is the mesh size of .

###### Remark 4.

The error terms in Theorems 1, 2 and 3 depend only on some geometric invariants of and the function since is a closed regular surface.

Similarly, We extend these scalars for each and for each to a matrix with

 bi,k=⎧⎪⎨⎪⎩βi,j if there exists j∈{1,2,⋯,N(i)} such that vk=vi,j,−1 if i=k,0 otherwise.. (48)

And, we have

 (49)

where , .

## 4 High-order approximations

In this section we shall discuss how to use the LTL method to obtain high-order approximations of these differential operators. The ideas are very simple. First, we shall propose an algorithm to construct a high-order approximation of the underlying surface . Second, we also give a method to obtain a high-order approximation of smooth functions on .Third, using these approximations, we can compute the differential quantities under consideration with high-order accuracies.

As before, we consider a triangular surface mesh , of the smooth surface where with mesh size is the list of vertices and is the list of triangles. To obtain a high-order approximation of the underlying surface around a vertex , we will try to construct a local parametrization by representing the smooth surface as locally a graph surface around the vertex . Let be the approximating normal vector at the vertex in as in (8). The approximating tangent plane of at is given by .

We can choose an orthonormal basis for the approximating tangent plane of at and obtain an orthonormal coordinates for vectors by . The approximating tangent plane is nearly tangential to the surface and the -coordinate is orthogonal to the -plane, the approximating tangent plane