Shape Analysis课程笔记02-Linear and Variational Problems
这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要复习课程所需的数学知识。
Linear Algebra
Vector
线性代数是本课程最重要的工具之一，因此我们首先来复习线性代数中的相关内容。向量(vector)是线性代数中最基本的概念，在本课程中我们默认向量都是列向量，这样向量内积还有二次型可以使用向量转置来表示。
更严格来说，向量的坐标是向量本身和基向量\(\mathbf{e}_k\)的内积。
Two Roles for Matrices
矩阵(matrix)是另一个非常重要的概念。尽管矩阵更常见的形式是数的方阵，但在本课程中我们会把矩阵更多地理解为满足如下关系的线性算子(linear operator)。更进一步还可以证明任意满足要求的线性算子和二次型都对应着一个矩阵。
因此矩阵既可以表示线性变换也可以表示向量内积，我们在使用时需要加以注意。
Einstein Notation
在本课程中我们还会大量使用Einstein记号来表示求和的过程。
使用Einstein记号的一个好处在于我们可以显式区分线性映射以及内积：
\[w^j = L_k^j v^k\]
从线性算子的角度来认识矩阵还可以帮助我们把线性代数的概念推广到无限维空间中。比如说类似于矩阵特征向量的概念，我们可以定义算子的特征函数。
Linear System of Equations
线性代数的另一大应用是求解线性系统\(A\mathbf{x} = \mathbf{b}\)。当矩阵\(A\)是可逆矩阵时，我们可以通过高斯消元等方法来求解线性方程组。这里需要额外注意的是在任何情况下我们都应该避免直接计算矩阵的逆阵。
除了线性方程组外，很多线性算子也可以表示为线性系统。此时同样可以利用相关的技术来进行求解。
在求解线性系统时，如果我们已知系统的一些特性还可以利用这些特性来加速求解。
目前常用的求解器可以分为直接求解或是迭代求解两类。
在几何问题中最常见的线性系统是稀疏系统。利用系统的稀疏性我们可以极大地减少存储空间，同时还可以加速求解过程。
当然在本课程中我们不需要亲自去实现一个线性系统求解器，像matlab等数值计算软件中都有非常成熟稳定的求解器以供调用。
Variational Calculus
Optimization Terminology
在正式介绍variational calculus前我们先来复习一下优化问题。一般来说一个优化问题由优化目标、等式约束以及不等式约束三部分组成。
Differential
在本课程中我们会大量使用微分(differential)的技术来求解优化问题。我们定义函数\(f\)在\(\mathbf{x}_0\)处沿方向\(\mathbf{v}\)的微分为：
\[d f_{\mathbf{x}_0} (\mathbf{v}) = \lim_{h \rightarrow 0} \frac{f(\mathbf{x}_0 + h \mathbf{v}) - f(\mathbf{x}_0)}{h}\]利用泰勒展开可以证明：
\[d f_{\mathbf{x}_0} (\mathbf{v}) = \sum_k \frac{\partial f}{\partial x^k} \cdot v^k = \nabla f(\mathbf{x}_0) \cdot \mathbf{v}\]因此我们可以把\(df\)看作是关于向量\(\mathbf{v}\)的一个线性算子。
Notions from Calculus
除此之外，函数的梯度、Jacobian矩阵以及Hessian矩阵也都可以表示为矩阵的形式。
Optimization to Root-Finding
对于无约束优化问题，我们可以利用梯度的概念将它转换为对梯度求根的问题。不过需要注意的是梯度的根只是优化问题的极值点，想要确认具体是极大值还是极小值还需要结合其他的条件。
Encapsulates Many Problems
实际上很多问题都可以转换为一个优化问题。
不过需要注意的是很多情况下使用优化工具来求解这些问题并不一定比直接求解原始问题来得高效。以二次型的优化问题为例，假设优化目标和约束为：
\[\begin{aligned} \min \ \ &\frac{1}{2} \mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x} + c \\ \text{s.t.} \ \ &M \mathbf{x} = \mathbf{v} \end{aligned}\]利用Lagrangian乘子法可以将其转换为一个无约束优化问题：
\[\Lambda (\mathbf{x}, \lambda) = \frac{1}{2} \mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x} + c + \lambda^T (M \mathbf{x} - \mathbf{v})\]分别对\(\mathbf{x}\)和\(\lambda\)求导可以得到：
\[\nabla_{\mathbf{x}} = A \mathbf{x} - \mathbf{b} + M^T \lambda = 0\] \[\nabla_{\lambda} = M \mathbf{x} - \mathbf{v} = 0\]整理后得到一个新的线性系统：
\[\begin{bmatrix} A & M^T \\ M & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x} \\ \lambda \end{bmatrix} = \begin{bmatrix} \mathbf{b} \\ \mathbf{v} \end{bmatrix}\]显然对于这种问题直接求解线性系统更加稳定和高效。类似地，最小二乘法这样的问题也建议使用构造线性系统的方式来进行求解。
Mesh Embedding
在后面的课程里我们还会介绍网格嵌入(mesh embedding)的问题，它的目标是把一个三维空间中的网格展开为一个二维平面。
实际上这样的问题可以转换为一个带约束的二次型优化问题，我们可以通过前面介绍过的方法来进行求解。
更进一步，边界点的位置对于这样的参数化问题起着重要的作用。当我们没有显式设置边界点的位置时则需要引入一些额外的边界条件防止解的退化。
Unconstrained Optimization
当优化问题没有各种约束时就得到了无约束优化问题(unconstrained optimization)。求解无约束优化问题的基本方法是梯度下降法，它在机器学习中有着非常广泛的应用。
当然梯度下降法并不是求解无约束优化问题最快的方法。如果我们可以计算优化目标的Hessian矩阵则可以使用Newton法来进行加速。
我们后面会介绍的形状插值就是使用这种方法来进行求解的。
对于这些优化问题也无需自己编写求解算法，直接调用现成的求解器即可。
Lagrangian Multiplier
Lagrangian乘子法(Lagrangian multiplier)是我们求解约束优化的重要方法。对于等式约束优化问题，优化目标\(f(\mathbf{x})\)仅在\(\nabla f(\mathbf{x})\)和\(\nabla g(\mathbf{x})\)共线时取最小值。
以二次型为例，利用Lagrangian乘子法不难发现带约束的二次型优化等价于求解特征值问题。
Lagrangian乘子法最大的价值在于它可以把约束优化问题转化为无约束优化问题(函数求根问题)，这样我们就可以使用各种求解器来进行计算。
在后面的课程中我们会陆续介绍约束优化以及Lagrangian乘子法在几何处理上的应用。
Gateaux Derivative
本节最后我们就可以介绍variational calculus了。在几何问题中我们处理的对象往往不是数字而是函数，在这种情况下我们不能直接通过求导的方法来进行优化。
为了对函数进行求导我们需要引入Gateaux导数(Gateaux derivative)的概念。Gateaux导数和方向导数的概念非常类似，它们的主要区别在于Gateaux导数是在函数空间中进行求导：
\[d \mathcal{F}_u (\psi) = \frac{d}{d h} \mathcal{F} [u + h \ \psi] \bigg\vert_{h=0} = \lim_{h \rightarrow 0} \frac{\mathcal{F} [u + h \ \psi] - \mathcal{F} [u]}{h}\]这里\(u\)是满足优化目标\(\mathcal{F}\)约束的函数，而\(\psi\)则可以理解为对\(u\)进行的扰动。
基于Gateaux导数我们还可以推导出各种常用的微分方程。假设在空间\(\Omega\)上存在向量场\(\mathbf{v}(\mathbf{x})\)，我们的目标是寻找一个标量场\(f(\mathbf{x})\)使得它的梯度\(\nabla f(\mathbf{x})\)尽可能接近\(\mathbf{v}(\mathbf{x})\)。因此我们可以定义一个优化问题：
\[\min_f \int_\Omega \Vert \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \Vert_2^2 \ d \mathbf{x}\]为了求解这个优化问题我们可以定义\(\mathcal{F} [f]\)为优化目标：
\[\mathcal{F} [f] = \int_\Omega \Vert \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \Vert_2^2 \ d \mathbf{x}\]假设\(f\)即为所求的最优函数，那么对于任意函数\(g\)我们可以定义函数\(p(h)\)：
\[p(h) = \mathcal{F} [f + h \cdot g] = \int_\Omega \Vert \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) - h \cdot \nabla g(\mathbf{x}) \Vert_2^2 \ d \mathbf{x}\]显然\(p(h)\)仅在\(h=0\)处取最小值，因此有：
\[\begin{aligned} p'(0) &= \frac{d}{d h} \mathcal{F} [f + h \cdot g] \bigg\vert_{h=0} \\ &= \frac{d}{d h} \int_\Omega \Vert \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) - h \cdot \nabla g(\mathbf{x}) \Vert_2^2 \ d \mathbf{x} \ \bigg\vert_{h=0} \\ &= \int_\Omega \frac{d}{d h} \Vert \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) - h \cdot \nabla g(\mathbf{x}) \Vert_2^2 \ \bigg\vert_{h=0} \ d \mathbf{x} \\ &= \int_\Omega \bigg( \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) - h \cdot \nabla g(\mathbf{x}) \bigg) \cdot \nabla g(\mathbf{x}) \ \bigg\vert_{h=0} \ d \mathbf{x} \\ &= \int_\Omega \bigg( \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \bigg) \cdot \nabla g(\mathbf{x}) \ d \mathbf{x} \\ &= 0 \end{aligned}\]利用分部积分可以得到：
\[\begin{aligned} \int_\Omega \bigg( \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \bigg) \cdot \nabla g(\mathbf{x}) \ d \mathbf{x} &= \oint_{\partial \Omega} g(\mathbf{x}) \ \bigg( \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \bigg) \cdot \hat{\mathbf{n}} \ d l(\mathbf{x}) - \int_\Omega g(\mathbf{x}) \ \nabla \cdot \bigg( \mathbf{v}(\mathbf{x}) - \nabla f(\mathbf{x}) \bigg) \ d\mathbf{x} \\ &= 0 \end{aligned}\]上式对任意函数\(g\)都成立。因此在\(\Omega\)内部有
\[\nabla \cdot \mathbf{v} = \nabla \cdot \nabla f\]实际上，上式即为著名的Poisson方程(Poisson’s equation)。
Variational Splines
以曲线为例，我们可以更好地理解variational calculus的意义和价值。样条曲线(splines)是在工程中广泛应用的一类曲线。样条曲线的一个经典问题是根据端点的位置和导数来确定曲线的表达式，同时我们希望曲线在满足端点约束的同时还要尽可能地光滑。因此我们可以把曲线二阶导数平方的积分作为优化目标：
\[E[f] = \int_0^1 \vert f''(t) \vert^2 \ dt\]而整个优化问题可以表示为：
\[\begin{aligned} \min_{f:[0,1] \rightarrow \mathbb{R}} \ &E[f] \\ \text{s.t.} \ &f(0) = a \\ &f(1) = b \\ &f'(0) = c \\ &f'(1) = d \end{aligned}\]
这个优化问题的难点在于满足约束的曲线是无穷多的，不过在继续讨论如何求解优化问题之前我们先来回顾一下有限维空间上的优化问题。对于\(\mathbb{R}^n\)上的优化问题，极值点需要满足梯度为0的条件，更进一步我们可以结合方向导数的概念得到等价的表达形式：
\[\nabla g(\mathbf{x}) = \mathbf{0} \iff \mathbf{v} \cdot \nabla g(\mathbf{x}) = 0, \ \forall \mathbf{v} \in \mathbb{R}^n\] \[\frac{d}{d h} g(\mathbf{x} + h \cdot \mathbf{v}) \bigg\vert_{h=0} = 0, \ \forall \mathbf{v} \in \mathbb{R}^n\]上面的推导说明极值点处\(g(\mathbf{x})\)任意方向上的方向导数均为0。更重要的是我们可以把原始优化问题转换为一个关于\(h\)的单变量优化问题，这在很多情况下会方便我们进行处理。
回到样条曲线问题上，根据优化目标的定义我们在\(f(t)\)上施加任意扰动\(p(t)\)都会增加能量项\(E\)的值。结合Gateaux导数的概念可以得到：
\[\begin{aligned} dE_f [p] &= \frac{d}{d h} E[f + h \cdot p] \bigg\vert_{h=0} \\ &= \frac{d}{d h} \int_0^1 \big( f''(t) + h \cdot p''(t) \big)^2 \ dt \ \bigg\vert_{h=0} \\ &= \int_0^1 2 \ ( f''(t) + h \cdot p''(t) ) \cdot p''(t) \ dt \ \bigg\vert_{h=0} \\ &= 2 \int_0^1 f''(t) \ p''(t) \ dt\\ &= 0 \end{aligned}\]利用分部积分可以得到：
\[\begin{aligned} 2\int_0^1 f''(t) \ p''(t) \ dt &= \bigg[ 2f''(t) \ p'(t) \bigg]_0^1 - 2 \int_0^1 f'''(t) \ p'(t) \ dt \\ &= \bigg[ 2f''(t) \ p'(t) \bigg]_0^1 - \bigg[ 2f''(t) \ p(t) \bigg]_0^1 + 2 \int_0^1 f^{(4)}(t) \ p(t) \ dt \end{aligned}\]
考虑到扰动后的函数\(f(t) + h \cdot p(t)\)仍然需要满足边界条件，我们可以得到扰动\(p(t)\)需要满足：
\[p(0) = p(1) = p'(0) = p'(1) = 0\]代回分部积分中得到：
\[\bigg[ 2f''(t) \ p'(t) \bigg]_0^1 - \bigg[ 2f''(t) \ p(t) \bigg]_0^1 + 2 \int_0^1 f^{(4)}(t) \ p(t) \ dt = 2 \int_0^1 f^{(4)}(t) \ p(t) \ dt = 0\]上式说明对于任意扰动\(p(t)\)，当能量取最小值时曲线\(f^{(4)}(t)\)需要恒为0：
\[f^{(4)} (t) \equiv 0\]换句话说我们可以使用一个三次多项式来表达样条曲线，这样定义的曲线是最光滑的。
