Shape Analysis课程笔记12-The Laplacian Operator
这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍Laplace算子的相关概念。
Laplace算子(Laplacian operator)在数学和物理等相关领域中有着大量的应用。在几何处理中，我们可以利用Laplace算子将几何形状变换到频域然后从频谱的角度进行分析，这样的方法称为spectral geometry。
在本节课中我们会从简单的情况入手逐渐揭示Laplace算子在不同几何对象上的意义和应用。
不过在具体介绍Laplace算子之前我们先回忆一下线性算子的概念。当算子满足分配律和数乘时就称其是一个线性算子，比如说矩阵就是一个典型的线性算子。
使用线性算子的一个优势在于我们可以不再受限于空间的维度，线性算子可以自然地推广到无穷维的空间上。
同时矩阵的谱定理对于线性算子也有相应的推广。
Line Segments
1D Spring Network
首先考虑一维的情况。假设有一维质点弹簧系统，其中每个质点的质量为\(m\)，质点数为\(n+1\)，弹簧刚度为\(k\)。我们固定两端质点使其无法移动，
根据牛顿运动定律，每个质点的运动方程可以表示为：
\[\begin{aligned} F &= m \cdot a = m \cdot \frac{d^2 u_i}{d t^2} \\ &= k (u_{i+1} - u_i) + k (u_{i-1} - u_i) \\ &= k (u_{i+1} + u_{i-1} - 2 u_i) \end{aligned}\]整理之后可以得到线性系统：
\[U'' (t) = -k L \ U(t)\] \[L = \begin{bmatrix} 2 &-1 \\ -1 & 2 &-1 \\ & & \ddots & \\ & & -1 & 2 &-1 \\ & & &-1 & 2 \end{bmatrix}\]对于微分项我们使用差分来进行近似：
\[\frac{d^2 u}{d x^2} \approx \frac{\frac{u_{i+1} - u_i}{h} - \frac{u_i - u_{i-1}}{h}}{h} = \frac{u_{i+1} + u_{i-1} - 2u_i}{h^2}\]这样线性算子\(L\)可以理解为求二阶导数：
\[L \approx - \frac{d^2}{d x^2}\]此时求解运动方程相当于求解偏微分方程(PDE)：
\[\frac{\partial^2 U}{\partial t^2} = c^2 \ \frac{d^2 U}{d x^2}\]上式即为一维波动方程(wave equation)，它描述了质点弹簧系统的振动。
波动方程更常见的形式为：
\[\frac{\partial^2 u}{\partial t^2} - \frac{\partial^2 u}{\partial x^2} = 0\]
Minus Second Derivative Operator
我们定义算子\(\mathcal{L}[\cdot]: u \rightarrow - \frac{\partial^2 u}{\partial x^2}\)，可以证明它是一个(半)正定算子。
对于离散的情况，\(\mathcal{L}\)即为上面定义的矩阵\(L\)。对于任意向量\(u\)有：
\[\begin{aligned} u^T L u &= \sum_i u_i (Lu)_i \\ &= \sum_i u_i (-u_{i-1} + 2 u_i - u_{i+1}) \\ &= 2 \sum_i u_i^2 - 2 \sum_i u_{i} u_{i-1} \\ &= \sum_i u_i^2 + \sum_i u_{i-1}^2 - 2 \sum_i u_{i} u_{i-1} \\ &= \sum_i (u_i - u_{i-1})^2 \\ &\geq 0 \end{aligned}\]而对于连续的情况则有：
\[\begin{aligned} \langle u, \mathcal{L}[u] \rangle &= \int_a^b u(x) \cdot - \bigg( \frac{\partial^2 u}{\partial x^2} \bigg) \ dx \\ &= \bigg[ - \frac{\partial u}{\partial x} \cdot u \bigg]_a^b + \int_a^b \frac{\partial u}{\partial x} \cdot \frac{\partial u}{\partial x} \ dx \\ &= \int_a^b \bigg( \frac{\partial u}{\partial x} \bigg)^2 \ dx \\ &\geq 0 \end{aligned}\]
Eigenfuctions of Second Derivative Operator
除了半正定外，在连续情况下算子\(\mathcal{L}\)还有类似于矩阵特征向量的概念，称为特征函数(eigenfunction)。容易验证\(\mathcal{L}\)的特征函数和对应特征值为：
\[\phi_k (x) = \sqrt{\frac{2}{l}} \sin \bigg( \frac{\pi k x}{l} \bigg)\] \[\lambda_k = \bigg( \frac{\pi k}{l} \bigg)^2\]而且不同特征函数之间是相互正交的：
\[\langle \phi_k, \phi_m \rangle = \delta_{km}\]
利用特征向量(函数)的正交性，我们可以对位移向量(函数)进行分解。对于离散情况，根据谱定理有：
\[U(t) = \sum_i U^i(t) \phi_i\]其中\(U^i(t)\)为随时间变化的系数。带入微分方程可以得到：
\[\begin{aligned} \sum_i (U^i)'' (t) \ \phi_i &= -k \mathcal{L} \ \sum_i U^i(t) \ \phi_i \\ &= -k \ \sum_i U^i(t) \ \lambda_i \phi_i \end{aligned}\]进一步有：
\[(U^i)'' (t) = -k \ U^i(t) \ \lambda_i\]因此我们只需要求解\(n\)个常微分方程(ODE)就可以得到波动方程的解。具体地，每个常微分方程的解为：
\[U^i (t) = a^i \sin (\sqrt{k \lambda_i} t) + b^i \cos (\sqrt{k \lambda_i} t)\]这表示我们可以通过\(L\)矩阵的特征值与特征向量来计算波动方程的闭式解，而对于连续的情况也有类似的结论。
Regions in \(\mathbb{R}^n\)
Typical Notation
在二维平面上同样存在波动方程，它描述了平面上的振动。其中算子\(\Delta\)即为二维平面上的Laplace算子。
\[\frac{\partial^2 u}{\partial t^2} = - \Delta u\]
Laplace算子有很多种定义方式，其中最经典的是使用梯度的散度来定义它。
Intrinsic Operator
Laplace算子的一个重要性质是它与坐标变换无关。记坐标变换前的函数为\(f(x)\)，变换后的函数为\(g(x) = f(R x + t)\)，我们可以证明刚体变换后Laplace算子的效果是一样的：
\[\frac{\partial}{\partial x^i} = \sum_j \frac{\partial y^j}{\partial x^i} \frac{\partial}{\partial y^j} = \sum_j R_{ji} \ \frac{\partial}{\partial y^j}\] \[\nabla_x = R^T \nabla_y\]因此，
\[\Delta_x = -\nabla_x^T \ \nabla_x = -\nabla_y R \ R^T \nabla_y = \Delta_y\]
Dirichlet Energy
Laplace算子可以描述函数的光滑性。我们定义函数的Dirichlet能量(Dirichlet energy)为梯度的大小(范数)在区域\(\Omega\)上的积分：
\[E[u] = \frac{1}{2} \int_\Omega \Vert \nabla u (\mathbf{x}) \Vert_2^2 \ d A(\mathbf{x})\]直观来看，Dirichlet能量越小表示函数在区域中的变化越平缓，函数也就越光滑。能够证明在给定边界条件下最小化Dirichlet能量等价于求解Laplace方程：
\[\min_{u} E[u] \Leftrightarrow \Delta u = 0\]
这里我们简单证明一下：
\[\begin{aligned} E[u] &= \frac{1}{2} \int_\Omega \Vert \nabla u (\mathbf{x}) \Vert_2^2 \ d A(\mathbf{x}) \\ &= \frac{1}{2} \int_\Omega \nabla u \cdot \nabla u \ d A(\mathbf{x}) \\ &= \frac{1}{2} \oint_{\partial \Omega} \dots \ dl - \frac{1}{2} \int_\Omega u (\mathbf{x}) \nabla \cdot \nabla u (\mathbf{x}) \ d A(\mathbf{x}) \\ &= \frac{1}{2} \int_\Omega u(\mathbf{x}) \Delta u(\mathbf{x}) \ d A(\mathbf{x}) + \text{boundary term} \\ &= \frac{1}{2} \langle u, \Delta u \rangle + \text{boundary term} \end{aligned}\]接下来利用variational calculus来处理内积项的极值问题。假设函数\(u\)即为所求，我们在\(u\)上添加扰动\(w\)得到新的Dirichlet能量\(E[u + hw]\)，其中\(h\)为步长且扰动\(w\)满足边界条件\(w_{\partial \Omega} = 0\)。由于\(u\)为最优函数，扰动后的Dirichlet能量需要满足：
\[\frac{d}{d h} E[u + hw] \bigg\vert_{h=0} = \frac{d}{d h} \bigg[ \frac{1}{2} \int_\Omega (u+hw) \ \Delta (u+hw) \ d A \bigg] = 0\]对Dirichlet能量进行展开可以得到：
\[\begin{aligned} \frac{d}{d h} E[u + hw] \bigg\vert_{h=0} &= \frac{d}{d h} \bigg[ \frac{1}{2} \int_\Omega (u+hw) \ \Delta (u+hw) \ d A \bigg] \\ &= \frac{1}{2} \int_\Omega \bigg[ w \ \Delta u + u \ \Delta w \bigg] \ d A \\ &= \int_\Omega w \ \Delta u \ d A \end{aligned}\]这表示对于任意扰动\(w\)，\(\Delta u\)需要恒为0。即在\(\Omega\)内部\(u\)需要满足Laplace方程：
\[\Delta u \equiv 0\]
实际上满足Laplace方程的函数也称为调和函数(harmonic function)，它们有很多优雅的几何性质。
Positivity, Self-Adjointness
Laplace算子还满足正定性和自伴随性质，因此可以看成正定矩阵以及自伴随矩阵在函数空间中的推广：
\[\begin{aligned} \langle f, \Delta f \rangle &= - \int_\Omega f(\mathbf{x}) \ \nabla \cdot \nabla f(\mathbf{x}) \ dA \\ &= \int_\Omega \Vert \nabla f(\mathbf{x}) \Vert_2^2 \ dA \\ &\geq 0 \end{aligned}\] \[\begin{aligned} \langle f, \Delta g \rangle &= - \int_\Omega f(\mathbf{x}) \ \nabla \cdot \nabla g(\mathbf{x}) \ dA \\ &= \int_\Omega \nabla f(\mathbf{x}) \cdot \nabla g(\mathbf{x}) \ dA \\ &= \langle \Delta f, g \rangle \end{aligned}\]
Laplacian Eigenfunctions
Laplace算子的特征方程与Dirichlet能量也有一定的联系，实际上小特征值对应的特征函数具有比较小的Dirichlet能量。
有时我们需要固定\(\Omega\)内的一些点然后求解Laplace方程，这种情况下可能会导致一些病态的问题。
Graphs
Basic Setup
在图这种拓扑结构上同样可以定义Laplace算子。首先我们定义图上的函数是将顶点到实数的映射。
这样图上的Dirichlet能量可以定义为相邻顶点上函数值之差的平方和。
Graph Laplacian
要形式化Dirichlet能量需要在图上定义差分算子\(D\)。
在此基础上就可以将图上的Dirichlet能量表示为差分的范数。
对于无向图，我们可以把差分算子\(D\)进行展开从而得到Laplace算子\(L=D^T D\)。可以证明这样定义的Laplace算子是对称且半正定的。
Properties
图上Laplace算子的第二小特征值称为Fiedler向量(Fiedler vector)，它描述了图的连通性和聚类性质。
满足\(Lx=0\)的节点表示该位置上的函数值等于相邻节点上值的平均。
使用Laplace算子来研究图性质的方法称为spectral graph theory。
不过需要注意我们可能无法使用Laplace算子来对图进行区分，一些图有着不同的拓扑形式但是相同的特征向量。
Submanifold Laplacians
Gradient Vector Field
在几何处理中我们希望可以利用Laplace算子来研究曲面或者说是流形的性质。我们定义曲面上的标量函数\(f\)是点到实数的映射，而\(f\)的微分\(df\)是切平面上的线性算子，它将切平面上的向量\(\mathbf{v}\)映射为\(f\)在该方向上的导数。
这样就可以定义曲面上的梯度场\(\nabla f(\mathbf{p})\)，它满足任意切平面上的向量\(\mathbf{v}\)有：
\[df_\mathbf{p} (\mathbf{v}) = \mathbf{v} \cdot \nabla f(\mathbf{p})\]
这里我们证明一下梯度场的存在性。首先在\(\mathbf{p}\)点的切平面上\(T_\mathbf{p} \mathcal{M}\)取一组基\(\{ \mathbf{b}_1, ..., \mathbf{b}_m \}\)，而它们在\(df_\mathbf{p}\)作用下可以得到：
\[a_i = df_\mathbf{p} (\mathbf{b}_i)\]同时定义一个矩阵\(\mathbf{g}\)，它的每个元素为基\(\mathbf{b}_i\)和\(\mathbf{b}_j\)的内积：
\[\mathbf{g}_{ij} = \mathbf{b}_i \cdot \mathbf{b}_j\]而它的逆阵\(\mathbf{g}^{-1}\)中的每个元素记为\(\mathbf{g}^{ij}\)。
接下来我们引入一个量\(\mathbf{x}\)，它满足：
\[\mathbf{x} = \sum_{ij} a_i \ \mathbf{g}^{ij} \ \mathbf{b}_j\]这样对于切平面上的向量\(\forall \mathbf{v} \in T_\mathbf{p} \mathcal{M}\)可以使用基向量进行分解：
\[\mathbf{v} = \sum_i v^i \ \mathbf{b}_i\]对\(\mathbf{v}\)和\(\mathbf{x}\)进行内积可以得到：
\[\begin{aligned} \mathbf{v} \cdot \mathbf{x} &= (v^i \ \mathbf{b}_i) \cdot (a_k \ \mathbf{g}^{kl} \ \mathbf{b}_l) \\ &= v^i \ a_k \ \mathbf{g}^{kl} \ \mathbf{g}_{il} \\ &= v^i \ a_k \ (\mathbf{g}^{-1} \mathbf{g})_i^k \\ &= v^i \ a_k \ \delta_k^i \\ &= v^i \ a_i \\ &= v^i \ df_\mathbf{p} (\mathbf{b}_i) \\ &= df_\mathbf{p} (\mathbf{v}) \end{aligned}\]注意这里我们使用了Einstein记号来化简推导过程。上式说明我们定义的量\(\mathbf{x}\)实际上就是\(f\)的梯度场\(\mathbf{x} = \nabla f(\mathbf{p})\)，而且可以证明这样定义的梯度场是唯一的。
Dirichlet Energy
有了梯度场后就可以定义曲面上的Dirichlet能量。和前面介绍过的定义一样，Dirichlet能量是梯度范数的积分。
Laplace–Beltrami Operator
除了使用Dirichlet能量来推导Laplace算子外，我们也可以使用内积来推导Laplace算子的定义。这样得到的Laplace算子也称为Laplace–Beltrami算子(Laplace–Beltrami operator)。
Divergence
或者我们也可以使用梯度的散度来推导Laplace算子。
Other Definitions
实际上还有一些其它的方法来推导Laplace算子，需要注意这些Laplace算子之间是存在一定差异的。
Applications
有了曲面上的Laplace算子就可以使用前面介绍过的特征函数等工具来分析曲面的各种性质。
另外Laplace算子在物理上有大量的应用，各种常见的PDE都与Laplace算子密切相关。
在球面上满足Laplace方程的函数称为球面谐波函数(spherical harmonics)，它在分子物理中有着重要的意义。
回忆前面介绍过的平均曲率法向，不难发现它实际上就是Laplace算子在曲面上作用的结果。
Divergence
本节课的extra content讨论了流形上散度的定义。在很多资料中都使用了梯度的散度来定义Laplace算子：
\[\Delta = - \nabla \cdot \nabla\]这样的定义在\(\mathbb{R}^n\)上没有什么大的问题，但在曲面或者说是流形上则是不严谨的，这主要是因为我们没有严格地定义过流形上的散度。
对于n维空间中的m维流形\(\mathcal{M} \subset \mathbb{R}^n\)，记它的向量场为\(\mathbf{v}: \mathcal{M} \rightarrow \mathbb{R}^n\)。我们对欧式空间中散度的定义进行推广，散度可以定义为流形到实数域的映射：
\[\nabla \cdot \mathbf{v} : \mathcal{M} \rightarrow \mathbb{R}\]通过引入\(\mathbf{p} \in \mathcal{M}\)处的局部参数化\(\phi: \mathbb{R}^m \rightarrow \mathbb{R}^n\)我们就可以定义流形上的散度。记\(\mathbb{R}^n\)中\(\mathbf{p}\)点处切平面normal space中的基向量为\(\bar{\mathbf{n}}_{m+1}, \dots, \bar{\mathbf{n}}_{n}\)，我们定义\(\bar{\phi}: \mathbb{R}^n \rightarrow \mathbb{R}^n\)如下：
\[\bar{\phi} (x^1, \dots, x^n) = \phi(x^1, \dots, x^m) + \sum_{k=m+1}^n x^k \cdot \bar{\mathbf{n}}_{k}\]即在前\(m\)维上与参数化\(\phi\)相同，而在后\(n-m\)维上取normal space中基向量的线性组合。显然如果我们令后\(n-m\)维上的坐标为0，\(\bar{\phi}\)即为参数化\(\phi\)。
实际上这样定义的\(\bar{\phi}\)在局部是一个单射(injective)，在足够小的邻域内甚至是双射(bijective)。因此\(\mathbf{p}\)点附近的任意点\(\mathbf{x} \in \mathbb{R}^n\)都可以通过\(\bar{\phi}^{-1}\)映射回切平面上。
更重要的是我们可以利用\(\bar{\phi}\)对向量场\(\mathbf{v}\)进行推广：
\[\bar{\mathbf{v}} (\mathbf{x}) = \mathbf{v} (\bar{\phi}^{-1} (\mathbf{x}))\]上面定义的\(\bar{\mathbf{v}} (\mathbf{x})\)是\(\mathbb{R}^n\)上的向量场，这样我们就可以把欧式空间中散度的定义推广到流形上：
\[\begin{aligned} \nabla \cdot \bar{\mathbf{v}} &= \sum_{k=1}^m \mathbf{e}_k \cdot d \bar{\mathbf{v}} (\mathbf{e}_k) + \sum_{k=m+1}^n \bar{\mathbf{n}}_k \cdot d \bar{\mathbf{v}} (\bar{\mathbf{n}}_k) \\ &= \sum_{k=1}^m \mathbf{e}_k \cdot d \bar{\mathbf{v}} (\mathbf{e}_k) \\ &= \sum_{k=1}^m \mathbf{e}_k \cdot d \mathbf{v} (\mathbf{e}_k) \\ &= \nabla \cdot \mathbf{v} \end{aligned}\]其中，\(\mathbf{e}_1, \dots, \mathbf{e}_m \in T_\mathbf{p} \mathcal{M}\)为流形上的基。
回忆\(\mathbb{R}^n\)上的散度定理(divergence theorem)：
\[\int_C \nabla \cdot \bar{\mathbf{v}} \ d \text{vol} = \int_{\partial C} \bar{\mathbf{v}} \cdot \mathbf{n}_C \ d A\]在\(\mathbf{p}\)点附近的邻域上，等式左边可以近似为：
\[\int_C \nabla \cdot \bar{\mathbf{v}} \ d \text{vol} \approx \nabla \cdot \mathbf{v} (\mathbf{p}) \cdot \text{Area} (C_0) \cdot \varepsilon^{n-m}\]而等式右边则有：
\[\begin{aligned} \int_{\partial C} \bar{\mathbf{v}} \cdot \mathbf{n}_C \ d A &= \int_{\partial C_0} \mathbf{v} \cdot \mathbf{n}_C \ d A + \int_{\partial C_1} \bar{\mathbf{v}} \cdot \mathbf{n}_C \ d A + \int_{\partial B} \bar{\mathbf{v}} \cdot \mathbf{n}_C \ d A \\ &\approx \int_{\partial B} \bar{\mathbf{v}} \cdot \mathbf{n}_C \ d A \\ &\approx \varepsilon^{n-m} \oint_{\partial C_0} \mathbf{v} \cdot \mathbf{n}_C \ d l \end{aligned}\]
这样\(\mathbf{p}\)点的散度可以近似为：
\[\begin{aligned} \nabla \cdot \mathbf{v} (\mathbf{p}) &= \sum_{k=1}^m \mathbf{e}_k \cdot d \mathbf{v} (\mathbf{e}_k) \\ &= \lim_{r \to 0} \frac{1}{\text{Area } C_0 (r)} \oint_{\partial C_0(r)} \mathbf{v} \cdot \mathbf{n}_C \ d l \end{aligned}\]而流形上的散度定理则可以表示为：
\[\int_{\mathcal{M}} \nabla \cdot \mathbf{v} \ dA = \oint_{\partial \mathcal{M}} \mathbf{v} \cdot \mathbf{n}_{\partial \mathcal{M}} \ dl\]