Kernel Multigrid on Manifolds

Consider a compact $d$-dimensional Riemannian manifold without boundary $\mathbb{M}$. Here we list some useful tools and their properties which hold in this setting. We direct the reader to [19, 20] for relevant background. Results about covariant derivatives and Sobolev spaces can be found in [21].

The tangent bundle is $T\mathbb{M}$ and the cotangent bundle is $T^{*}\mathbb{M}$. We denote by $T^{k}_{r}\mathbb{M}$ the vector bundle of tensors with contravariant order $k$ and covariant order $r$; we say type $(k,r)$ for short. Thus $T\mathbb{M}=T^{1}_{0}\mathbb{M}$ and $T^{*}\mathbb{M}=T^{0}_{1}\mathbb{M}$. We will denote the fiber at $x\in\mathbb{M}$ by $T^{k}_{r}\mathbb{M}_{x}$. The space of tensor fields of type $(k,r)$ (known also as sections; i.e., maps $\mathbf{S}:\mathbb{M}\to T_{r}^{k}\mathbb{M}$ with $\mathbf{S}(x)\in T_{r}^{k}\mathbb{M}_{x}$ for every $x\in\mathbb{M}$) is denoted $\mathcal{T}_{r}^{k}\mathbb{M}$.

In this article, we are concerned primarily with covariant tensors (i.e., tensors of type $(0,k)$), so we use the short hand notation $T_{k}\mathbb{M}=T_{k}^{0}\mathbb{M}$ and $\mathcal{T}_{k}\mathbb{M}=\mathcal{T}_{k}^{0}\mathbb{M}$.

For a chart $(U,\phi)$ for $\mathbb{M}$ from which we get the usual vector fields $\frac{\partial}{\partial x^{j}}$ and forms $dx^{j}$ ($1\leq j\leq d$), which act as local bases for $T\mathbb{M}$ and $T^{*}\mathbb{M}$ over $U$. These can be used to generate bases for tensor fields. In particular, for given covariant rank $k$ and $\boldsymbol{\mathbf{i}}=(i_{1},\dots,i_{k})\in\{1,\dots,d\}^{k}$ we have basis element $dx^{\boldsymbol{\mathbf{i}}}:=dx^{i_{1}}\cdots dx^{i_{k}}$ . This allows us to write $\mathbf{S}\in\mathcal{T}_{k}\mathbb{M}$ in coordinates as $\mathbf{S}(x)=\sum_{\boldsymbol{\mathbf{i}}\in\{1,\dots,d\}^{k}}(S(x))_{\boldsymbol{\mathbf{i}}}dx^{\boldsymbol{\mathbf{i}}}$.

Because $\mathbb{M}$ is a Riemannian manifold, at each $x\in\mathbb{M}$, $T\mathbb{M}_{x}$ has an inner product $\langle\cdot,\cdot\rangle_{x}$ and induced norm $\|\cdot\|_{x}$. This means that there is by $g\in\mathcal{T}_{2}\mathbb{M}$ (the metric tensor), so that for tangent vectors in $V,W\in T\mathbb{M}_{x}$ we have $g(x)(V,W)=\langle V,W\rangle_{x}$. The inner product extends to the dual: for cotangent vectors $\mu,\nu\in T^{*}\mathbb{M}_{x}$ we have $\langle\mu,\nu\rangle_{x}=\sum\mu_{j}\nu_{k}g^{j,k}$, where $\sum g_{j,k}g^{k,\ell}=\delta_{j,\ell}$. From this, it naturally extends to tensors. For $\mathbf{T},\mathbf{S}\in{T}_{k}\mathbb{M}_{x}$, written in coordinates as $\mathbf{T}=\sum_{\boldsymbol{\mathbf{j}}\in\{1,\dots,d\}^{k}}T_{\boldsymbol{\mathbf{j}}}\,d{x^{j_{1}}}\dots d{x^{j_{k}}}$ and $\mathbf{S}=\sum_{\boldsymbol{\mathbf{i}}\in\{1,\dots,d\}^{k}}S_{\boldsymbol{\mathbf{i}}}\,d{x^{i_{1}}}\dots d{x^{i_{k}}}$, we have

We denote the Riemannian distance on $\mathbb{M}$ by $\mathrm{dist}:\mathbb{M}\times\mathbb{M}\to[0,\infty)$; it is given by the formula $\mathrm{dist}(x,y)=\inf\{\int_{a}^{b}\|\gamma^{\prime}(t)\|_{\gamma(t)}\mathrm{d}t\mid\gamma(a)=x,\ \gamma(b)=y\}$ where the infimum is taken over piecewise smooth curves connecting $x$ and $y$. The Riemannian metric gives rise to a volume form $\mathrm{d}\mu=\sqrt{\det(g(x))}\mathrm{d}x$. By compactness, there exist constants $0<\alpha_{\mathbb{M}}\leq\beta_{\mathbb{M}}$ so that

or any ball $B(x,r):=\{y\in\mathbb{M}\mid\mathrm{dist}(x,y)<r\}$ centered at $x\in\mathbb{M}$ and having radius $0<r<\mathrm{diam}(\mathbb{M})$.

For a finite subset $\Xi\subset\mathbb{M}$, we can define the following useful quantities: the separation distance, $q$ of $\Xi$ and the fill distance, $h$, of $\Xi$ in $\mathbb{M}$. They are given by

The finite sets considered throughout this paper will be quasiuniform, with mesh ratio $\rho:=h/q$ bounded by a fixed constant. For this reason, quantities which are controlled above or below by a power of $q$ can be likewise controlled by a power of $h$ – in short, whenever possible, we express estimates in terms of the fill distance $h$, allowing constants to depend on $\rho$. For instance, the cardinality $|\Xi|$ is bounded above and below by

Similarly, if $f:[0,\infty)\to\mathbb{R}$ is continuous then for any $x\in\mathbb{M}$,

Moreover, we use the following notation: $\mathbb{R}^{A}=\{f:A\to\mathbb{R}\}$ denotes the functions from the set $A$ to $\mathbb{R}$. Of course, if $A$ is a discrete finite set, we can identity $\mathbb{R}^{A}\cong\mathbb{R}^{|A|}$.

The covariant derivative $\nabla$ maps tensor fields of type $(r,s)$ to fields of type $(r,s+1)$. Its adjoint (with respect to the $L_{2}$ inner products on the space of sections of $T_{r}^{s}(\mathbb{M})$) is denoted $\nabla^{*}$. For functions, this is fairly elementary. The covariant derivative of a (scalar) function $f:\mathbb{M}\to\mathbb{R}$ equals its exterior derivative; in coordinates, we have

For a 1-form $\omega=\sum_{j=1}^{d}\omega_{j}dx^{j}$, we have

A direct calculation shows $\int_{\mathbb{M}}f(x)\nabla^{*}\omega(x)\mathrm{d}\mu(x)=\int_{\mathbb{M}}\langle\omega(x),\nabla f(x)\rangle_{x}\mathrm{d}\mu(x)$.

For $\Omega\subset\mathbb{M}$ Sobolev space $W_{p}^{k}(\Omega)$ is defined to be the set of functions $f:\Omega\to\mathbb{R}$ which satisfy

where the $\ell$-th order covariant derivatives can be found in [21, Section 2.2].

For a scalar function $f:\mathbb{M}\to\mathbb{R}$, the Laplace-Beltrami operator is given in local coordinates as $\Delta f=\sum_{j=1}^{d}\sum_{k=1}^{d}\frac{1}{\sqrt{|\det g|}}\frac{\partial}{\partial x^{j}}\bigl{(}\sqrt{|\det g|}g^{jk}\frac{\partial}{\partial x^{k}}f\bigr{)}$. Thus for scalar functions $\Delta f=-\nabla^{*}\nabla f$. For any integer $k\in\mathbb{N}$ and $p\in(1,\infty)$, the Bessel-potential norm $\|(1-\Delta)^{k/2}f\|_{L_{p}(\mathbb{M})}$ is equivalent to $\|f\|_{W_{p}^{k}(\mathbb{M})}$, as demonstrated in [22, Theorem 4(ii)] (although when $k=1$ and $p=2$, the two norms are equal; this can be observed directly).

Lemma 3.2 from [21] applies, so there are uniform constants $\Gamma_{1},\Gamma_{2}$ and $r_{\mathbb{M}}>0$ so that the family of exponential maps $\{\mathrm{Exp}_{x}:B(0,r_{\mathbb{M}})\to\mathbb{M}\mid x\in\mathbb{M}\}$ (which are diffeomorphisms taking $0$ to $x$) provides local metric equivalences: for any open set $\Omega\subset B(0,r_{\mathbb{M}})\subset\mathbb{R}^{d}$, we have

This shows that $W_{p}^{k}(\mathbb{M})$ can be endowed with equivalent norms using a partition of unity $(\varphi_{j})_{j\leq N}$ subordinate to a cover $\{\mathcal{O}_{j}\}_{j\leq N}$ with associated charts $\psi_{j}:\mathcal{O}_{j}\to\mathbb{R}^{d}$ to obtain $\|u\|_{W_{p}^{k}(\mathbb{M})}^{p}\sim\sum_{j=1}^{N}\|(\varphi_{j}u)\circ\psi_{j}^{-1}\|_{W_{p}^{k}(\mathbb{R}^{d})}^{p}$. Here constants of equivalence depend on the partition of unity and charts.

A useful result in this setting, which we will use explicitly in this article, but is also behind a number of background results in section 2.5, is the following zeros estimate [23, Corollary A.13 ], which holds for Sobolev spaces: If $u\in W_{2}^{m}(\mathbb{M})$ satisfies $u\left|{}_{X}\right.=0$, then for any $k\leq m$ we have

Here $C_{\mathrm{zeros}}$ depends on $m$ and $\mathbb{M}$. The result can also be obtained on bounded, Lipschitz domain $\Omega\subset\mathbb{M}$ that satisfies a uniform cone condition, with the cone having radius $R_{\Omega}\leq r_{\mathbb{M}}/3$, although the constant $C_{\mathrm{zeros}}$ will depend on the aperture of the cone in that case. See [23, Theorem A.11] for a precise statement and definitions of the involved quantities.

We write the operator $\mathcal{L}$ in divergence form:

Here $\textsc{a}_{0}$ is a smooth function, $\textsc{a}_{1}$ is a smooth tensor field of type $(1,0)$ and $\textsc{a}_{2}$ is a smooth tensor field of type $(2,0)$ which generates the field $\textsc{a}_{2}^{\flat}$ of type (1,1). In coordinates, $\textsc{a}_{2}$ has the form $\sum_{j=1}^{d}\sum_{k=1}^{d}a^{jk}\frac{\partial}{\partial x_{j}}\frac{\partial}{\partial x_{k}}$, and $\textsc{a}_{2}^{\flat}$ is $\sum_{j=1}^{d}\sum_{k=1}^{d}\tensor{a}{{}_{k}^{j}}dx_{k}\frac{\partial}{\partial x_{j}}$ with $\tensor{a}{{}_{k}^{j}}=\sum_{\ell=1}^{d}g_{k\ell}a^{\ell j}$. Here we have used the index lowering operator ^♭ which ensures, for any $\mu,\nu\in T_{1}\mathbb{M}_{x}$, that $\langle\textsc{a}_{2}^{\flat}(\mu,\cdot),\nu\rangle_{x}=\textsc{a}_{2}(\mu,\nu)$. This follows by a direct calculation from the above expressions in coordinates (index lowering and the ^♭ operator are discussed in [19, p. 341] or [20, p. 27]).

Furthermore, we require $c_{0}>0$ to be a constant so that

As a basic example, we consider $\textsc{a}_{2}(x)=\langle\cdot,\cdot\rangle_{x}$, $\textsc{a}_{1}=0$ and $\textsc{a}_{0}=1$; then (6) holds with $c_{0}=1$. In this case, $a^{jk}=g^{jk}$, $\tensor{a}{{}_{k}^{j}}=\delta_{j,k}$ and $\mathcal{L}=1+\nabla^{*}\nabla=1-\Delta$.

With the identity $\int_{\mathbb{M}}u(\textsc{a}_{1}\nabla u)=\frac{1}{2}\int_{\mathbb{M}}(\textsc{a}_{1}\nabla(u^{2}))=\frac{1}{2}\int_{\mathbb{M}}(\nabla^{*}\textsc{a}_{1})u^{2}$, the second part of (6) ensures that

We have also that $\int_{\mathbb{M}}u(x)\nabla^{*}(\textsc{a}_{2}^{\flat}\nabla u)(x)\mathrm{d}\mu(x)=\int_{\mathbb{M}}\langle\textsc{a}_{2}^{\flat}\nabla u,\nabla u\rangle_{x}\mathrm{d}\mu(x)$. The definition of $\textsc{a}_{2}^{\flat}$ ensures that this equals $\int_{\mathbb{M}}\textsc{a}_{2}(x)\bigl{(}\nabla u(x),\nabla u(x)\bigr{)}\mathrm{d}\mu(x)$, and the first part of (6) then ensures that

Thus, (6) guarantees that the bilinear form $a(u,v):=\int_{\mathbb{M}}v\mathcal{L}u$, defined initially for smooth functions, is bounded on $W_{2}^{1}(\mathbb{M})$ and is coercive. Thus, we have that the energy quasi-norm $[u]^{2}_{\mathcal{L}}:=a(u,u)$ satisfies the metric equivalence

(if $a$ is symmetric, this is a norm, and we write $\|u\|_{\mathcal{L}}=[u]_{\mathcal{L}}$.)

Please note that this excludes differential operators with a null space at the moment. Having in mind the time dependent problems, this assumption is justified. The technical more challenging analysis for those more general operators is left to future research.

We fix $f\in L_{2}(\mathbb{M})$ and consider $u\in W^{1}_{2}(\mathbb{M})$ as solution to

Regularity estimates ([24, Chapter 5.11, Theorem 11.1]) yield that $u\in W^{2}_{2}(\mathbb{M})$ and $\|u\|_{W^{2}_{2}(\mathbb{M})}\leq C\|f\|_{L_{2}(\mathbb{M})}$. Consider now a family of finite dimensional subspaces $(V_{h})$ with $V_{h}\subset W_{2}^{1}(\mathbb{M})$ associated to a parameter²²2In the sequel, these will be kernel spaces generated by a subset $\Xi\subset\mathbb{M}$, and $h=h(\Xi,\mathbb{M})$ will be the fill distance. For now, we consider a more abstract setting, with $h>0$ only playing a role in establishing the approximation property $\mathrm{dist}_{W_{2}^{1}}(g,V_{h})\leq Ch\|g\|_{W_{2}^{2}(\mathbb{M})}$ below. $h>0$. Define $P_{V_{h}}:W^{1}_{2}(\mathbb{M})\to V_{h}$ so that

Because $a\left(u-P_{V_{h}}u,P_{V_{h}}u-v\right)=0$, the classical Céa Lemma holds:

And thus we get

This can be improved by a Nitsche-type argument to get the following result, whose proof can be found in many textbooks on numerical methods for partial differential equations.

We point out, that Lemma 1 does not depend on the choice of a specific basis in the finite dimensional space $V_{h}$.

We consider a continuous function $\phi:\mathbb{M}\times\mathbb{M}\to\mathbb{R}$, the kernel, which satisfies a number of analytic properties which we explain in this section.

Most important is that $\phi$ is conditionally positive definite with respect to some (possibly trivial) finite dimensional subspace³³3generally a space spanned by some eigenfunctions of the Laplace-Beltrami operator $\Pi\subset C^{\infty}(\mathbb{M})$. Conditional positive definiteness with respect to $\Pi$ means that for any $\Xi\subset\mathbb{M}$, the collocation matrix $\bigl{(}\phi(\xi,\zeta)\bigr{)}_{\xi,\zeta}$ is positive definite on the vector space $\{a\in\mathbb{R}^{\Xi}\mid(\forall p\in\Pi)\,\sum_{\xi\in\Xi}a_{\xi}p(\xi)=0\}$. As a result, if $\Xi$ separates elements of $\Pi$, then the space

has dimension $\Xi$. In case $\Pi=\{0\}$, the kernel is positive definite, and the collocation matrix is strictly positive definite on $\mathbb{R}^{\Xi}$, and $V_{\Xi}=\mathrm{span}_{\xi\in\Xi}\phi(\cdot,\xi)$.

For a conditionally positive definite kernel, there is an associated reproducing kernel semi-Hilbert space $\mathcal{N}(\phi)\subset C(\mathbb{M})$ with the property that $\Pi=\mathrm{null}(\|\cdot\|_{\mathcal{N}(\phi)})$ and that for any $a\in\mathbb{R}^{\Xi}$ for which $\sum_{\xi\in\Xi}a_{\xi}\delta_{\xi}\perp\Pi$, the following identity

holds for all $f\in\mathcal{N}(\phi)$. It follows that if $\Xi$ separates elements of $\Pi$, interpolation with $V_{\Xi}$ is well defined on $\Xi$ and the projection $I_{\Xi}:\mathcal{N}(\phi)\to V_{\Xi}$ is orthogonal with respect to the semi-norm on $\mathcal{N}(\phi)$. Of special interest are (conditionally) positive definite kernels that have Sobolev native spaces.

As a consequence, the kernel Galerkin solution $u_{\Xi}\in V_{\Xi}$ to $\mathcal{L}u=f$ with $f\in W_{2}^{j+1}(\mathbb{M})$ satisfies $\|u-u_{\Xi}\|_{W_{2}^{1}(\mathbb{M})}\leq Ch^{j-1}\|f\|_{W_{2}^{j+1}(\mathbb{M})}$. More importantly, for our purposes, the hypothesis of Lemma 1 is satisfied by the space $V_{\Xi}$. Indeed, by (10), we have the following approximation property:

which, together with a smoothing property, forms the backbone of the convergence theory for the multigrid method.

We now discuss the stiffness matrix and some of its properties. Most of these have appeared in [26], with earlier versions for the sphere appearing in [21] and [8].

A consequence of the results of this section is that the problem of calculating the Galerkin solution to $\mathcal{L}u=f$ from $V_{\Xi}$ involves treating a problem whose condition number grows like $\mathcal{O}(h^{-2})$ – this is the fundamental issue that the multigrid method seeks to overcome.

The analysis map for $\bigl{(}\chi_{\xi}\bigr{)}_{\xi\in\Xi_{\ell}}$ with respect to the bilinear form $a$ defined in (8) is

The analysis map is a surjection.

The synthesis map is

The range of the synthesis map is clearly $V_{\Xi}$; in other words, it is the natural isomorphism between Euclidean space and the finite dimensional kernel space; indeed, (15) shows that it is bounded above and below between $L_{2}(\mathbb{M})$ and $\ell_{2}(\Xi)$. By abusing notation slightly, we write $\bigl{(}\sigma^{\ast}_{\Xi}\bigr{)}^{-1}:V_{\Xi}\to\mathbb{R}^{\Xi}$. This permits a direct matrix representation of linear operators on $V_{\Xi}$ via conjugation: $S\mapsto\mathbf{S}:=(\sigma_{\Xi}^{*})^{-1}S\sigma_{\Xi}^{*}\in\mathbb{R}^{\Xi\times\Xi}$. Furthermore, by the Riesz property (15), we have

A simple calculation shows that $a\left(\sigma^{\ast}_{\Xi}(\boldsymbol{\mathbf{w}}),v\right)=\left\langle\boldsymbol{\mathbf{w}},\sigma_{\Xi}(v)\right\rangle_{2}$, so $\sigma_{\Xi}^{\ast}$ is the $a$-adjoint of $\sigma_{\Xi}$. Of course, when $a$ is symmetric, we also have $\left\langle\sigma_{\ell}(v),\boldsymbol{\mathbf{w}}\right\rangle_{2}=a\left(v,\sigma^{\ast}_{\ell}(\boldsymbol{\mathbf{w}})\right)$.

The stiffness matrix is defined as

It represents the operator $\mathcal{L}$ on the finite dimensional space $V_{\Xi}$. Using the analysis and synthesis maps, $\boldsymbol{\mathbf{A}}_{\Xi}=\sigma_{\Xi}\circ\sigma^{\ast}_{\Xi}:\left(\mathbb{R}^{\Xi},(\cdot,\cdot)_{2}\right)\to\left(\mathbb{R}^{\Xi},(\cdot,\cdot)_{2}\right):\boldsymbol{\mathbf{c}}\mapsto\left(a(\chi_{\xi},\chi_{\zeta})\right)_{\xi,\zeta\in\Xi}\boldsymbol{\mathbf{c}}$, and we have that the Galerkin projector $P_{\Xi}:W_{2}^{1}(\mathbb{M})\to V_{\Xi}$ can be expressed as $P_{\Xi}=\sigma^{\ast}_{\Xi}\left(\sigma_{\Xi}\circ\sigma^{\ast}_{\Xi}\right)^{-1}\sigma_{\Xi}$.