Shape Analysis课程笔记17-Optimization on Manifolds
这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍流形上的优化问题。
Optimization with a Manifold Constraint
在机器学习和计算机视觉中很多问题都可以建模为在流形上进行处理的优化问题。

三维重建中对相机姿态的优化就是典型的例子:相机的旋转矩阵并不是任意的矩阵,它需要满足正交阵的约束。

冷冻电镜(Cryo-EM)在重建分子的结构时也会遇到类似的问题。

而对于如何优化传感器的布置,我们同样可以把它表示为一个优化问题,此时传感器的布置形式在刚体变换下是相互等价的。

Manifold Optimization Algorithms
Typical Approach
在过去,处理这些问题的经典方法是把流形上的优化变量嵌入到高维空间中,然后求解高维的约束优化问题。比如说三维空间的旋转构成李群SO(3),在优化旋转时可以先直接优化一个3×3的矩阵,然后再处理旋转矩阵的约束。

然而这种处理方法的缺陷在于它往往会具有非常复杂的表达形式,使得求解过程变得非常困难。

从几何上来看,求解这样的约束优化问题相当于在红色曲线上进行移动直到找到优化目标\(f\)的最优解。此时我们可以通过投影的方式进行处理,即先忽略约束直接求解无约束优化然后把找到的解投影到约束空间中。


当然投影过程往往也并不容易,尤其是在约束条件比较复杂的时候。

Instrinsic Perspective
因此更现代的处理方法是结合流形自身的结构,利用流形自身的性质来进行求解。


从更高的视角来看这样的方法相当于直接在流形上进行优化,从而无需考虑复杂的约束条件。

Gradient Descent
以梯度下降(gradient descent)为例,我们常用的梯度下降公式都假定了变量位于欧式空间中。


而在流形上,我们需要先把欧式空间上的梯度转换为切空间上的向量,然后按照切向量的方向沿着测地线进行移动从而实现优化过程。

当我们考虑了流形的结构时往往会得到更好的优化算法,此时的约束优化问题也一般会具有更优雅的数学性质。

需要说明的是流形上的优化是一门非常复杂的学科,想要更系统的学习可以参考专业的书籍或是优化软件手册。



Tangent Space
在继续介绍流形上的优化算法前我们先回忆一下方向导数以及切空间的概念。


在切空间的基础上,优化目标\(f\)的梯度实际上对应了切空间上的切向量。在整个流形上,这些切向量构成了一个梯度场。

Exponential Map
流形上的指数映射(exponential map)定义了从起点\(\mathbf{p}\)出发沿切向量\(\mathbf{v}\)方向的测地线移动的过程。

当然流形上测地线的计算往往是非常困难的,我们几乎无法给出一般流形上测地线的解析形式。

Retraction
不过我们可以利用retraction来放松测地线的约束,它是对指数映射的一种局部近似。

Riemannian Manifold
除此之外,我们还需要引入黎曼流形(Riemannian manifold)的概念。简单来说,当流形\(\mathcal{M}\)上配备了内积\(g_\mathbf{p}\)后就得到了黎曼流形。

在内积的基础上就可以定义长度、角度等几何量。

实际上很多几何处理任务都隐含了黎曼流形的假设。


Riemannian Gradient
在黎曼流形上我们可以定义黎曼梯度(Riemannian gradient)。和欧式空间中的梯度相比,黎曼梯度需要结合流形上的度量张量。

更具体来说,假设流形\(\mathcal{M}\)上\(\mathbf{p}\)点附近存在局部参数化\(\Phi: \mathbb{R}^n \rightarrow \mathcal{M}\),我们定义参数域上的内积为向量映射到\(\mathcal{M}\)切平面上的内积:
\[g_{\mathbf{x}}(\mathbf{v}, \mathbf{w}) = d \Phi_\mathbf{x} (\mathbf{v}) \cdot d \Phi_\mathbf{x} (\mathbf{w})\]进一步取\(\mathbf{v}\)和\(\mathbf{w}\)为参数域上的基向量,则内积可以表示为一个(半正定)矩阵\(g_{ij}\):
\[g_{ij} (\mathbf{x}) = g_{\mathbf{x}}(\mathbf{e}_i, \mathbf{e}_j)\]
更进一步可以证明流形\(\mathcal{M}\)上的梯度依赖于参数域\(\mathbb{R}^n\)上的内积定义:
\[\nabla_{\mathcal{M}} f = g^{-1} \nabla_{\mathbb{R}^n} f\]
结合上面的推导可以得到黎曼梯度下降(Riemannian gradient descent)算法,我们只需要把欧式空间中的梯度替换为黎曼梯度并且把沿梯度下降换成指数映射即可。

Extensions
对于线搜索(line search)这样的技术也可以进行类似的推广。

我们还可以把牛顿法这样的二阶方法也推广到黎曼流形上,不过这需要在黎曼流形上定义Hessian矩阵。

我们甚至可以在黎曼流形上重新定义凸性。

Example: Unit Sphere
接下来我们以\(\mathbb{R}^n\)上的单位球为例推导流形上的优化算法。显然此时的优化变量\(\mathbf{x}\)位于\(n-1\)维的流形上。

在\(\mathbf{p} \in S^{n-1}\)点上法向量由坐标直接给出,因此\(\mathbf{p}\)点的切空间可以表示为所有与法向量\(\mathbf{p}\)正交的向量集合:
\[T_\mathbf{p} S^{n-1} = \{ \mathbf{v} \in \mathbb{R}^n | \mathbf{v} \cdot \mathbf{p} = 0 \}\]
对于\(\mathbb{R}^n\)上的函数\(f\),我们需要把梯度投影到流形上:
\[\nabla_{S^{n-1}} f (\mathbf{p}) = (I_{n \times n} - \mathbf{p} \mathbf{p}^T) \nabla_{\mathbb{R}^{n}} f (\mathbf{p})\]
接下来我们设计retraction,在单位球上可以使用指数映射的解析形式或是先在切空间上进行移动再投影到流形上。当移动的步长比较小时,两种方法的结果是非常近似的。

Example: Stiefel Manifold
单位球构成的流形可能比较简单,接下来我们介绍Stiefel流形(Stiefel manifold)上的优化。Stiefel流形\(V_k (\mathbb{R}^n)\)是\(n\)维空间中具有\(k\)个正交列向量的矩阵构成的流形,实际上单位球可以理解为Stiefel流形在\(k=1\)情况下的特例。

考虑Stiefel流形上的切空间。根据定义,切空间可以表示为经过\(\mathbf{X}\)处曲线切向量构成的集合。带入\(V_k (\mathbb{R}^n)\)的约束可以得到:
\[T_\mathbf{X} V_k (\mathbb{R}^n) = \{ \xi: \mathbf{X}^T \xi + \xi^T \mathbf{X} = 0 \}\]
接下来定义Stiefel流形上的retraction:
\[R_\mathbf{X} (\xi) = (\mathbf{X} + \xi) (\mathbf{I} + \xi^T \xi)^{-\frac{1}{2}}\]可以验证上式满足retraction的要求。


Optimization Examples
到这里我们就可以求解流形上的优化问题。比如说\(\mathbb{R}^n\)上的瑞利商(Rayleigh quotient)问题就可以转换为单位球\(S^{n-1}\)上的优化问题:

类似地,regularized PCA算法可以等价地表示为Stiefel流形上的优化。我们稍后会介绍具体的推导过程。实际上很多PCA算法的变体都可以转化为Stiefel流形上的优化问题。


除此之外,很多复杂的约束优化问题都可以利用流形结构来进行处理。





Manifold Optimization for PCA Problems
PCA
PCA(principal component analysis)是机器学习中非常经典的方法。前面我们提到过PCA可以理解为Stiefel流形\(V_k (\mathbb{R}^n)\)上的优化问题:
\[\max_{\mathbf{X} \in V_k (\mathbb{R}^n)} \Vert \mathbf{X}^T \mathbf{A} \Vert_F^2\]其中矩阵\(\mathbf{A}\)是由已有数据构成的矩阵。求解PCA的经典方法是对数据矩阵进行SVD:
\[\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T\]利用Frobenius范数的性质可以证明PCA的最优解是\(\mathbf{A}\)矩阵最大的前\(k\)个奇异值所对应的左奇异向量:

Robust/Sparse PCA
在很多问题中我们还会对\(\mathbf{X}\)引入一些正则项,此时则可以使用流形上的优化进行求解。比如说我们希望\(\mathbf{X}\)具有稀疏性,此时需要在优化目标上添加L1正则:
\[\max_{\mathbf{X} \in V_k (\mathbb{R}^n)} \Vert \mathbf{X}^T \mathbf{A} \Vert_F^2 - \alpha \sum_{ij} \vert X_{ij} \vert\]
而为了提高算法的鲁棒性,我们可以通过L2范数重新定义优化目标:
\[\min_{\mathbf{X} \in V_k (\mathbb{R}^n)} \sum_i \Vert \mathbf{a}_i - \mathbf{X} \mathbf{X}^T \mathbf{a}_i \Vert_2\]
Manifold Optimization
在这些情况下我们可能无法像标准PCA那样直接求解,因此在Stiefel流形上进行迭代优化就是一种合适的处理方法。以标准PCA为例,我们只需要推导Stiefel流形上的梯度,然后带入流形上的retraction进行迭代即可。尽管对于标准PCA问题这种在流形上进行优化的方法比较复杂,但实际上我们只需要能够计算\(\mathbf{R}^n\)上优化目标的梯度,就可以处理各种PCA变体的优化问题:
