DDG课程笔记3-Laplace算子

这个系列是对CMU离散微分几何课程(CMU CS 15-458/858: Discrete Differential Geometry)的总结和回顾。本节我们介绍几何处理中最实用的工具-Laplace算子。

Overview

Laplace算子(Laplacian)在不同的领域里都有非常重要的应用。在物理学中我们可以通过Laplace算子来描述不同的物理现象：

在几何处理中Laplace算子描述了网格的局部信息，因此也称为微分坐标(differential coordinates)。同时，Laplace算子也和平均曲率以及网格的频率信息有着紧密的联系：

欧式空间中的Laplace算子可以定义为函数二阶导数的和：

$$\mathrm{\Delta }u=\sum _{i=1}^n\frac{{\partial }^{2}u}{\partial x_i^2}$$

因此Laplace算子描述了函数在局部邻域内的性质，我们可以通过Laplace算子来计算函数的极值点以及曲率：

在图论中，Laplace算子则描述了节点和它相邻节点平均值的相对关系：

$$(Lu{)}_i=(\frac{1}{\text{deg}(i)}\sum _{ij\in E}{u}_j)-u_i$$

在社交网络中这表示每一个人总是与和他相关的人具有相似的属性(观点)。

Definitions

尽管Laplace算子在不同的领域中有着不同的定义，我们可以看出Laplace算子描述了一个对象和它相邻对象平均值的关系。同时，针对Laplace算子的不同性质也可以给出很多等价的定义。考虑Laplace算子可以看做是偏导数的和，这里给出Laplace算子在流形上的定义：

$$\mathrm{\Delta }u=\sum _{i=1}^n\sum _{j=1}^n\frac{1}{\sqrt{det(g)}}\frac{\partial }{\partial x_i}(\sqrt{det(g)}(g^{-1}{)}_{ij}\frac{\partial }{\partial x_j}u)$$

可以验证在欧式空间下上式会退化为我们熟悉的二阶导数和的形式$\mathrm{\Delta }u=\sum _{i=1}^n\frac{{\partial }^{2}u}{\partial x_i^2}$。尽管这个式子十分复杂，在实际的计算中基本不会直接使用它。

对函数$u$进行泰勒展开并保留前两项得到：

$$u(x_0+w)\approx u(x_0)+⟨\mathrm{\nabla }u(x_0),w⟩+\frac{1}{2}⟨{\mathrm{\nabla }}^{2}u(x_0)w,w⟩$$

不难发现Laplace算子等于Hessian矩阵的迹：

$$\mathrm{\Delta }u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=\text{tr}({\mathrm{\nabla }}^{2}u)$$

在曲面上使用协变导数(covariant derivative)来代替欧式坐标中的导数，得到：

$${\mathrm{\nabla }}_{X,Y}^{2}u(X,Y)=({\mathrm{\nabla }}_{X}du)(Y)$$

此时Laplace算子仍然等于Hessian矩阵的迹：

$$\mathrm{\Delta }u=\text{tr}({\mathrm{\nabla }}^{2}u)=\text{tr}(\mathrm{\nabla }du)$$

Laplace算子还可以定义为梯度的散度，此时表达式为：

$$\mathrm{\Delta }u=\mathrm{\nabla }\cdot \mathrm{\nabla }u$$

此外，我们还可以通过随机游走(布朗运动(Brownian motion))来认识Laplace算子。对于单一粒子的随机游走其状态可由热传导方程进行刻画，带入初始条件$\varphi$得到：

$$\{\begin{array}{rl}& \frac{d}{dt}u=\mathrm{\Delta }u\\ & u{|}_{t=0}=\varphi \end{array}\Rightarrow u(t,x)={\int }_{\mathrm{\Omega }}{k}_t(x,y)\varphi (y)dy=\mathbb{E}[\varphi (X_t)]$$

其中$k_t(x,y)$称为热核(heat kernel)函数，表示从具有单位温度的热源$x$经时间$t$传播到$y$处时的温度，这里可以认为是粒子从$x$出发在$t$时刻位于$y$的概率密度；$X_t$表示随机游走的粒子在$t$时刻的位置。

最后来介绍Laplace算子与Dirichlet能量(Dirichlet energy)的关系。定义Dirichlet能量为：

$$E(u)={\int }_M|\mathrm{\nabla }u{|}^{2}dA=-{\int }_{M}u\mathrm{\Delta }udA=-⟨⟨\mathrm{\Delta }u,u⟩⟩$$

Dirichlet能量描述了函数的光滑程度：Dirichlet能量越小表示函数越光滑、变化越平稳；Dirichlet能量越大表示函数越粗糙、变化越剧烈。同时可以证明最小化Dirichlet能量相当于求解Laplace方程：

$$u=\arg minE(u)\iff \mathrm{\Delta }u=0$$

Properties

Laplace算子有非常丰富的性质，这里只介绍一些比较常用的性质：

• 线性算子：$\mathrm{\Delta }(f+\alpha g)=\mathrm{\Delta }f+\alpha \mathrm{\Delta }g$
• 常值函数满足$\mathrm{\Delta }\alpha =0$
• 对刚体变换保持不变
• 可交换性：对等距映射$\eta$满足$\eta \circ \mathrm{\Delta }=\mathrm{\Delta }\circ \eta$
• 自伴随：${\int }_{M}f\mathrm{\Delta }gdA={\int }_{M}g\mathrm{\Delta }fdA$
• (半)负定：${\int }_{M}f\mathrm{\Delta }fdA\le 0$

值得一提的是Laplace算子还具有谱性质(spectral properties)。类似于实对称矩阵的特征向量，Laplace算子存在正交的特征函数：

$$\mathrm{\Delta }{\varphi }_i={\lambda }_i{\varphi }_i$$ $${\int }_M\Vert \varphi \Vert dA=1$$

其中$0\ge {\lambda }_1\ge {\lambda }_2...$。上式表明可以把Laplace算子看做是无穷维的(半)负定矩阵，特征函数对应矩阵的特征向量。此外特征函数构成了函数空间中的基，因此可以对任意函数$f$进行分解：

$$f=\sum _i{\alpha }_i{\varphi }_i$$

利用特征函数分解还可以将Dirichlet能量写成级数的形式：

$$E(f)={\int }_M|\mathrm{\nabla }f{|}^{2}dA=-{\int }_{M}f\mathrm{\Delta }fdA=\sum _i{\alpha }_i^2(-{\lambda }_i)$$

Discretization via DEC

接下来推导离散Laplace算子的表达式，这里我们使用离散外微分来进行推导。首先给出光滑曲面上Laplace算子的表达式：

$$\mathrm{\Delta }u=\star d\star du$$

使用外微分形式的好处在于我们可以直接进行离散，得到网格上的Laplace算子：

$$\mathrm{\Delta }u={\star }_0^{-1}{d}_0^T{\star }_1{d}_{0}u$$

上式表明我们可以直接使用我们已经推导的离散微分算子来构造出Laplace算子。实际上我们还可以把它们写成更紧凑的形式：对于0-form $u$，它的离散外微分$du$等于顶点函数值的差：

$$(du{)}_{ij}={\int }_{e_{ij}}du={\int }_{\partial e_{ij}}u=u_j-u_i$$

在$du$上使用Hodge star算子等于将$du$缩放到dual边$e_{ij}^{\star }$上：

$$(\star du{)}_{ij}=\frac{|e_{ij}^{\star }|}{|e_{ij}|}(u_j-u_i)$$

其中，

$$\frac{|e_{ij}^{\star }|}{|e_{ij}|}=\frac{1}{2}(\cot {\alpha }_j+\cot {\beta }_j)$$

接着对dual 1-form $(\star du{)}_{ij}$进行微分，相当于把dual cell边界上的积分加起来，得到dual 2-form $(d\star du{)}_i$：

$$(d\star du{)}_i={\int }_{C_i}d\star du={\int }_{\partial C_i}\star du=\sum _j\frac{|e_{ij}^{\star }|}{|e_{ij}|}(u_j-u_i)=\frac{1}{2}\sum _j(\cot {\alpha }_j+\cot {\beta }_j)(u_j-u_i)$$

此时我们得到的dual 2-form $(d\star du{)}_i$ 与0-form $(\star d\star du{)}_i$ 只相差dual cell的面积，相当于是 $(\star d\star du{)}_i$ 在dual cell上面的积分。通常情况下我们会把dual cell的面积提出来单独处理，进而定义离散Laplace算子为：

$$(\mathrm{\Delta }u{)}_i=\frac{1}{2}\sum _j(\cot {\alpha }_j+\cot {\beta }_j)(u_j-u_i)$$

上式称为cotan-Laplace算子(cotan Laplacian)。注意到cotan-Laplace算子与平均曲率向量$(HN{)}_i$只相差一个符号，因此在有些资料中也会直接将平均曲率向量$(HN{)}_i$定义为Laplace算子，此时的Laplace算子为(半)正定。

Applications

Heat Equation

有了离散Laplace算子我们就可以在网格上解微分方程，比如说热传导方程(heat equation)：

$$\frac{d}{dt}u=\mathrm{\Delta }u$$

由于离散Laplace算子和平均曲率向量具有相同的形式，实际上平均曲率流的迭代算法也可以看做是求解热传导方程。取$u$为顶点坐标的函数，带入热传导方程并进行离散得到：

$$\frac{u^{k+1}-u^k}{h}=L{u}^k$$

其中$L$为Laplace算子的矩阵形式，$h$为时间步长。对方程进行整理可以得到迭代格式：

$$u^{k+1}=(I+hL)u^k$$

上式称为热传导方程的显式迭代(explicit update)。实际工程中更常用的是隐式迭代(implicit update)，此时需要将递推式修改为：

$$\frac{u^{k+1}-u^k}{h}=L{u}^{k+1}$$

整理得到：

$$(I-hL)u^{k+1}=u^k$$

和显式迭代相比隐式迭代在数值上更稳定，因此大部分情况下推荐使用隐式迭代算法。

Poisson Equation

在几何处理中另一类常见的微分方程是泊松方程(Poisson equation)：

$$\{\begin{array}{rlrl}& \mathrm{\Delta }u=-f,& & \text{on}\mathrm{\Omega }\\ & u=g,& & \text{on}\partial \mathrm{\Omega }\end{array}$$

实际上泊松方程可以看做是稳态传热情况下的热传导方程。假设$\mathrm{\Omega }$内存在热源$f$且边界条件为$g$，此时热传导方程可表示为：

$$\{\begin{array}{rlrl}& \frac{d}{dt}u=\mathrm{\Delta }u+f,& & \text{on}\mathrm{\Omega }\\ & u=g,& & \text{on}\partial \mathrm{\Omega }\end{array}$$

当$\mathrm{\Omega }$达到热平衡时$u$关于时间的导数为0，此时热传导方程即为泊松方程。

除此之外，泊松方程也可以理解为最小化Dirichlet能量。假设已知向量场$X$，我们希望寻找一个势场$u$使得$u$可以表示$X$，因此构造Dirichlet能量：

$$E(u)={\int }_M|\mathrm{\nabla }u-X{|}^{2}dA$$

可以证明最小化Dirichlet能量等价于求解泊松方程：

$$\mathrm{\Delta }u=\mathrm{\nabla }\cdot X$$

Boundary Conditions

我们知道边界条件对于微分方程起着至关重要的作用，因此求解泊松方程之前还要介绍一些常见边界条件。常见边界条件包括Dirichlet边界和Neumann边界，其中Dirichlet边界表示在边界上指定了函数值而Neumann边界则表示在边界上指定了函数梯度。以一维函数为例，Dirichlet边界指定了函数在端点的值而Neumann边界则指定了函数的导数：

同时需要注意对于Dirichlet边界方程一定有解，而对于Neumann边界则可能会出现无解的情况：

以二维Laplace方程为例，根据散度定理有：

$${\int }_{\mathrm{\Omega }}0dA={\int }_{\mathrm{\Omega }}\mathrm{\Delta }udA={\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot \mathrm{\nabla }udA={\int }_{\partial \mathrm{\Omega }}n\cdot \mathrm{\nabla }uds$$

也就是说Laplace方程有解需要保证方向导数 $n\cdot \mathrm{\nabla }u=\frac{\partial u}{\partial n}$ 沿边界的积分为0，因此对于Neumann边界需要小心处理。

最后我们来讨论一下如何求解泊松方程。对于无边界的情况，泊松方程退化为$\mathrm{\Delta }u=f$，对应的离散形式为：

$$Lu=Mf$$

这里$M$为顶点dual cell面积构成的对角矩阵，它来自于cotan-Laplace算子中忽略的${\star }_0^{-1}$项。此时微分方程转换为线性方程组，利用Laplace矩阵$L$的对称性可通过Cholesky分解来进行求解。

当泊松方程只存在Dirichlet边界时，我们可以将网格内部和边界分开来考虑。此时泊松方程可表示为：

$$\left[\begin{array}{cc}{L}_{II}& L_{IB}\\ L_{BI}& L_{BB}\end{array}\right]\left[\begin{array}{c}{u}_I\\ u_B\end{array}\right]=\left[\begin{array}{cc}{M}_I& 0\\ 0& M_B\end{array}\right]\left[\begin{array}{c}{f}_I\\ f_B\end{array}\right]$$

下标$I$和$B$分别表示网格的内部和边界。由于已知边界上的函数值$u_B$，我们只需要在网格内部进行求解即可。带入$u_B$整理得到：

$$L_{II}{u}_I=M_{II}{f}_I-L_{IB}{u}_B$$

如果存在Neumann边界则需要对泊松方程进行一定的修改。对于边界上的顶点，回忆Laplace算子是函数在dual cell上的积分，此时需要更新为：

$$(\mathrm{\Delta }u{)}_i=\frac{1}{2}(g_a+g_b)+\frac{1}{2}\sum _j(\cot {\alpha }_j+\cot {\beta }_j)(u_j-u_i)$$

其中$g_a$和$g_b$分别为给定方向导数在边界上的积分。

将$g_a$和$g_b$移到方程右侧，得到包含Neumann边界的泊松方程：

$$\left[\begin{array}{cc}{L}_{II}& L_{IB}\\ L_{BI}& L_{BB}\end{array}\right]\left[\begin{array}{c}{u}_I\\ u_B\end{array}\right]=\left[\begin{array}{c}{M}_I{f}_I\\ M_B{f}_B-g\end{array}\right]$$

此时按无边界的情况进行求解即可。

Reference

• CHAPTER 6: The Laplacian
• Laplace-Beltrami: The Swiss Army Knife of Geometry Processing
• Poisson Problem