DDG课程笔记3-Laplace算子

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

Overview

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

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

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

\[\Delta u = \sum_{i=1}^n \frac{\partial^2 u}{\partial x_i^2}\]

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

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

\[(L u)_i = \big( \frac{1}{\text{deg}(i)} \sum_{ij \in E} u_j \big) - u_i\]

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

Definitions

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

\[\Delta u = \sum_{i=1}^n \sum_{j=1}^n \frac{1}{\sqrt{\det (g)}} \frac{\partial}{\partial x_i} \Big( \sqrt{\det (g)} (g^{-1})_{ij} \frac{\partial}{\partial x_j} u \Big)\]

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

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

\[u(x_0 + w) \approx u(x_0) + \langle \nabla u(x_0), w \rangle + \frac{1}{2} \langle \nabla^2 u(x_0) w, w \rangle\]

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

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

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

\[\nabla_{X, Y}^2 u (X, Y) = (\nabla_X du) (Y)\]

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

\[\Delta u = \text{tr}(\nabla^2 u) = \text{tr}(\nabla du)\]

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

\[\Delta u = \nabla \cdot \nabla u\]

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

\[\left\{ \begin{aligned} &\frac{d}{dt} u = \Delta u \\ &u \vert_{t=0} = \phi \end{aligned} \right. \Rightarrow u(t, x) = \int_\Omega k_t(x, y) \phi(y) dy = \mathbb{E}[\phi (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 \vert \nabla u \vert^2 dA = - \int_M u \Delta u dA = - \langle \langle \Delta u, u \rangle \rangle\]

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

\[u = \arg \min E(u) \Leftrightarrow \Delta u = 0\]

Properties

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

线性算子：\(\Delta (f + \alpha g) = \Delta f + \alpha \Delta g\)
常值函数满足\(\Delta \alpha = 0\)
对刚体变换保持不变
可交换性：对等距映射\(\eta\)满足\(\eta \circ \Delta = \Delta \circ \eta\)
自伴随：\(\int_M f \Delta g dA = \int_M g \Delta f dA\)
(半)负定：\(\int_M f \Delta f dA \leq 0\)

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

\[\Delta \phi_i = \lambda_i \phi_i\] \[\int_M \Vert \phi \Vert dA = 1\]

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

\[f = \sum_i \alpha_i \phi_i\]

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

\[E(f) = \int_M \vert \nabla f \vert^2 dA = -\int_M f \Delta f dA = \sum_i \alpha_i^2 (-\lambda_i)\]

Discretization via DEC

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

\[\Delta u = \star d \star d u\]

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

\[\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{\vert e_{ij}^\star \vert}{\vert e_{ij} \vert} (u_j - u_i)\]

其中，

\[\frac{\vert e_{ij}^\star \vert}{\vert e_{ij} \vert} = \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{\vert e_{ij}^\star \vert}{\vert e_{ij} \vert} (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算子为：

\[(\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 = \Delta u\]

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

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

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

\[u^{k+1} = (I + h L) u^k\]

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

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

整理得到：

\[(I - h L) u^{k+1} = u^k\]

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

Poisson Equation

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

\[\left\{ \begin{aligned} & \Delta u = -f, & & \text{on} \ \Omega \\ & u = g, & & \text{on} \ \partial \Omega \end{aligned} \right.\]

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

\[\left\{ \begin{aligned} & \frac{d}{dt} u = \Delta u + f, & & \text{on} \ \Omega \\ & u = g, & & \text{on} \ \partial \Omega \end{aligned} \right.\]

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

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

\[E(u) = \int_M \vert \nabla u - X \vert^2 dA\]

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

\[\Delta u = \nabla \cdot X\]

Boundary Conditions

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

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

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

\[\int_\Omega 0 \ dA = \int_\Omega \Delta u \ dA = \int_\Omega \nabla \cdot \nabla u \ dA = \int_{\partial \Omega} n \cdot \nabla u \ ds\]

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

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

\[L u = M f\]

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

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

\[\begin{bmatrix} L_{II} & L_{IB} \\ L_{BI} & L_{BB} \\ \end{bmatrix} \begin{bmatrix} u_I \\ u_B \end{bmatrix} = \begin{bmatrix} M_I & 0 \\ 0 & M_B \\ \end{bmatrix} \begin{bmatrix} f_I \\ f_B \end{bmatrix}\]

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

\[L_{II} u_I = M_{II} f_I - L_{IB} u_B\]

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

\[(\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边界的泊松方程：

\[\begin{bmatrix} L_{II} & L_{IB} \\ L_{BI} & L_{BB} \\ \end{bmatrix} \begin{bmatrix} u_I \\ u_B \end{bmatrix} = \begin{bmatrix} M_I f_I \\ M_B f_B - g \end{bmatrix}\]

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