Shape Analysis课程笔记10-Distance Metrics and Embeddings

这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍度量空间的相关概念。

在上一章节中我们介绍了测地线和测地距离的相关概念，而在这一章节我们会考虑它的逆问题，即已知样本点之间的距离然后恢复它们本身的几何。这一问题也称为inverse distance problem。

实际上这样的问题不仅在几何处理中有大量的应用，在机器学习等相关领域中也同样有着重要的意义。

具体来说，我们的基本目标是在给定样本之间距离的基础上计算样本的嵌入(embedding)。

Metric Space

我们先来介绍一下度量空间(metric space)的概念。在集合\(M\)上赋予一个距离函数\(d\)就得到了度量空间，当然这个距离函数\(d\)需要满足一些性质：

非负性(non-negativity)：\(d(x, y) \geq 0\)
同一性：\(d(x, y) = 0 \Leftrightarrow x = y\)
对称性(symmetry)：\(d(x, y) = d(y, x)\)
三角不等式(triangle inequality)：\(d(x, z) \leq d(x, y) + d(y, z)\)

除了我们常见的欧式空间外，配备了测地距离的流形以及\(C^\infty(\mathbb{R})\)的函数空间都是度量空间。

Embedding a Metric Space

由于度量空间是非常丰富的，在实践中我们更关心如何把一个广义上的度量空间映射到一个方便进行处理的空间(如欧式空间)上。如果存在一个映射能够将沟通两个度量空间，而且保证映射前后样本之间的距离是不变的，则称这样的一个映射是等距(isometry)的。

显然我们希望可以把一个任意的度量空间通过等距映射变换到欧式空间上。不过遗憾的是这基本是不可能的，即使是有限度量空间也不能保证一定可以通过等距映射转换为欧式空间。

但对于有限度量空间可以证明存在等距映射将样本嵌入到某个由无穷范数定义的欧式空间\(l_\infty(\mathbb{R}^n)\)上。

这里我们简单证明一下上面的结论。假设存在有限度量空间\((M, d)\)，其中\(M = \{ a_1, a_2, ..., a_n \}\)。接下来定义映射\(\phi: M \rightarrow \mathbb{R}^n\)为：

\[\phi(a_k) = \big( d(a_1, a_k), d(a_2, a_k), ..., d(a_n, a_k) \big)\]

不难证明\(\phi\)是一个等距映射。首先考虑映射后\(\phi(a_i)\)和\(\phi(a_j)\)之间的距离：

\[\Vert \phi(a_i) - \phi(a_j) \Vert_\infty = \max_k \vert d(a_i, a_k)- d(a_j, a_k) \vert\]

对于右边有：

\[\begin{aligned} \max_k \vert d(a_i, a_k)- d(a_j, a_k) \vert &\geq \vert d(a_i, a_i)- d(a_j, a_i) \vert \\ &= d(a_i, a_j) \end{aligned}\]

同时根据三角不等式，对于任意\(1 \leq k \leq n\)满足：

\[d(a_i, a_k) \leq d(a_i, a_j) + d(a_j, a_k) \Leftrightarrow d(a_i, a_k) - d(a_j, a_k) \leq d(a_i, a_j)\] \[d(a_j, a_k) \leq d(a_j, a_i) + d(a_i, a_k) \Leftrightarrow d(a_j, a_k) - d(a_i, a_k) \leq d(a_i, a_j)\]

即

\[\vert d(a_i, a_k) - d(a_j, a_k) \vert \leq d(a_i, a_j)\] \[\max_k \vert d(a_i, a_k) - d(a_j, a_k) \vert \leq d(a_i, a_j)\]

综合上面的推导可以得到：

\[d(a_i, a_j) \leq \max_k \vert d(a_i, a_k) - d(a_j, a_k) \vert \leq d(a_i, a_j)\] \[\max_k \vert d(a_i, a_k) - d(a_j, a_k) \vert = d(a_i, a_j)\]

即原始度量空间\(M\)可以通过等距映射\(\phi\)变换到度量空间\(l_\infty(\mathbb{R}^n)\)上。

Approximate Embedding

Approximate Isometry

尽管上面的推导说明我们可以把有限度量空间映射到欧式空间中，但映射后的度量空间\(l_\infty(\mathbb{R}^n)\)其维度依赖于原始度量空间中的样本数量。在大多数情况下上面的方法会构造出一个非常高维的空间，这既不利于计算机表示也不利于我们进行分析和可视化。

因此我们可以放松等距映射的要求，只要两个度量空间是近似等距(approximate isometry)即可。对于原始度量空间上的距离函数\(\rho\)和映射后度量空间的距离函数\(\mu\)，我们使用distortion来描述映射\(f\)偏离等距映射的程度。显然对于等距映射distortion严格为1；而随着distortion的增长，\(f\)会不断地偏离等距映射。

Fréchet Embedding

具体地，Fréchet embedding就是这样的一种映射方法。我们只需要在原始空间中选择一系列子集\(S_1, ..., S_r \subseteq M\)，然后构造出\(\mathbb{R}^r\)中的向量即可。

Bourgain’s Theorem

同时，Bourgain定理(Bourgain’s theorem)指出我们可以把有限度量空间近似等距地映射到一个相对低维的欧式空间中。

Embedding Metrics into Euclidean Space

Euclidean Problem

接下来考虑如何从度量来恢复几何的问题。假设在欧式空间\(\mathbb{R}^m\)中存在样本\(\{ \mathbf{x}_1, ..., \mathbf{x}_n \}\)，我们已知任意两个样本之间的距离为\(\mathbf{P}_{ij} = \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2^2\)，我们的目标是重建\(\mathbb{R}^m\)中样本的坐标。

实际上这个问题并不简单。比如说我们可以对样本进行平移和旋转，此时样本之间的距离是保持不变的。这说明对于样本重建问题可能有多个可行解。

Gram Matrix

根据样本矩阵我们可以定义Gram矩阵(Gram matrix)：

\[\mathbf{G} = \mathbf{X}^T \mathbf{X}\]

Gram矩阵的每一个元素表示对应样本的内积：

\[\mathbf{G}_{ij} = \mathbf{x}_i \cdot \mathbf{x}_j\]

Classical MDS

利用Gram矩阵就可以把样本的重建问题分解为2步：首先利用\(\mathbf{P}\)矩阵恢复Gram矩阵，然后再利用Gram矩阵来重建样本坐标。

Recover Gram Matrix

假设样本的中心为原点，即\(\sum_i \mathbf{x}_i = \mathbf{0}\)，利用欧式空间下内积的定义有：

\[\begin{aligned} \mathbf{P}_{ij} &= \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2^2 \\ &= \Vert \mathbf{x}_i \Vert_2^2 + \Vert \mathbf{x}_j \Vert_2^2 - 2 \mathbf{x}_i \cdot \mathbf{x}_j \\ &= \mathbf{G}_{ii} + \mathbf{G}_{jj} - 2 \mathbf{G}_{ij} \end{aligned}\]

注意到\(\mathbf{G}_{ii}\)和\(\mathbf{G}_{jj}\)是Gram矩阵的对角元素，我们可以把\(P\)矩阵使用Gram矩阵来表示：

\[\mathbf{P} = -2 \mathbf{G} + \text{diag}(\mathbf{G}) \ \mathbb{1}^T + \mathbb{1} \ \text{diag}(\mathbf{G})^T\]

根据中心化条件\(\sum_i \mathbf{x}_i = \mathbf{0}\)，可以得到：

\[\mathbf{P} \ \mathbf{1} = \sum_i \mathbf{x}_i = \mathbf{0}\] \[\mathbf{G} \ \mathbf{1} = \begin{bmatrix} \mathbf{x}_1 \cdot \mathbf{x}_1 + \mathbf{x}_1 \cdot \mathbf{x}_2 + \dots \mathbf{x}_1 \cdot \mathbf{x}_n \\ \vdots \\ \mathbf{x}_n \cdot \mathbf{x}_1 + \mathbf{x}_n \cdot \mathbf{x}_2 + \dots \mathbf{x}_n \cdot \mathbf{x}_n \end{bmatrix} = \begin{bmatrix} \mathbf{x}_1 \cdot (\sum_i \mathbf{x}_i) \\ \vdots \\ \mathbf{x}_n \cdot (\sum_i \mathbf{x}_i) \end{bmatrix} = \mathbf{0}\]

因此我们可以定义一个辅助矩阵\(\mathbf{J} = \mathbf{I} - \frac{1}{n} \mathbf{1} \mathbf{1}^T\)，此时

\[\begin{aligned} -\frac{1}{2} \mathbf{J}^T \mathbf{P} \mathbf{J} &= -\frac{1}{2} (\mathbf{I} - \frac{1}{n} \mathbf{1} \mathbf{1}^T)^T \ (-2 \mathbf{G} + \text{diag}(\mathbf{G}) \ \mathbb{1}^T + \mathbb{1} \ \text{diag}(\mathbf{G})^T) \ (\mathbf{I} - \frac{1}{n} \mathbf{1} \mathbf{1}^T) \\ &= -\frac{1}{2} (-2 \mathbf{G} + \text{diag}(\mathbf{G}) \ \mathbb{1}^T + \mathbb{1} \ \text{diag}(\mathbf{G})^T + \frac{2}{n} \mathbf{1} \mathbf{1}^T \mathbf{G} - \frac{1}{n} \mathbf{1} \mathbf{1}^T \text{diag}(\mathbf{G}) \ \mathbb{1}^T - \frac{1}{n} \mathbf{1} \mathbf{1}^T \mathbb{1} \ \text{diag}(\mathbf{G})^T ) (\mathbf{I} - \frac{1}{n} \mathbf{1} \mathbf{1}^T) \\ &= -\frac{1}{2} (-2 \mathbf{G} + \text{diag}(\mathbf{G}) \ \mathbb{1}^T - \frac{1}{n} \text{tr}(\mathbf{G}) \mathbf{1} \mathbf{1}^T) (\mathbf{I} - \frac{1}{n} \mathbf{1} \mathbf{1}^T) \\ &= -\frac{1}{2} (-2 \mathbf{G} + \text{diag}(\mathbf{G}) \ \mathbb{1}^T - \frac{1}{n} \text{tr}(\mathbf{G}) \mathbf{1} \mathbf{1}^T -\frac{1}{n} \text{diag}(\mathbf{G}) \mathbf{1}^T \mathbf{1} \mathbb{1}^T + \frac{1}{n^2} \text{tr}(\mathbf{G}) \mathbf{1} \mathbf{1}^T \mathbf{1} \mathbf{1}^T) \\ &= \mathbf{G} \end{aligned}\]

这样我们就从\(\mathbf{P}\)矩阵恢复了Gram矩阵。

Recover \(\mathbf{X}\) from Gram Matrix

在继续推导如何重建样本坐标前，我们先来看一下Gram矩阵的性质。首先Gram矩阵的元素都是非负的，\(\mathbf{G} \succeq 0\)，而且它的秩是有界的\(\text{rank}(\mathbf{G}) \leq \min \{ m, n \}\)。根据特征值分解有：

\[\mathbf{G} = \mathbf{Q} \ \mathbf{\Lambda} \ \mathbf{Q}^T\]

其中矩阵\(\mathbf{Q}\)不一定是一个方阵，且\(\mathbf{Q}^T \mathbf{Q} = \mathbf{I}\)。这样就可以重建出样本的坐标：

\[\mathbf{\bar{X}} = \sqrt{\mathbf{\Lambda}} \ \mathbf{Q}^T\]

实际上前面推导的算法称为classical multidimensional scaling(classical MDS)，它在机器学习以及数据可视化中都有重要的应用。

Landmark MDS

classical MDS算法的一个变体是landmark MDS。注意到：

\[\begin{aligned} \mathbf{\bar{X}} \mathbf{P} &= (\sqrt{\mathbf{\Lambda}} \ \mathbf{Q}^T) \ (-2 \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T + \text{diag}(\mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T) \ \mathbb{1}^T + \mathbb{1} \ \text{diag}(\mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T)^T) \\ &= -2 \mathbf{\Lambda}^{3/2} \mathbf{Q}^T + \mathbf{\bar{X}} \ \text{diag}(\mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T) \ \mathbb{1}^T \end{aligned}\] \[\mathbf{\bar{X}} (\mathbf{P} - \text{diag}(\mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T) \ \mathbb{1}^T ) = -2 \mathbf{\Lambda}^{3/2} \mathbf{Q}^T = -2 \mathbf{\Lambda} \mathbf{\bar{X}}\]

整理后得到：

\[\mathbf{\bar{X}} = -\frac{1}{2} \mathbf{\Lambda}^{-1} \mathbf{\bar{X}} (\mathbf{P} - \text{diag}(\mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T) \ \mathbb{1}^T )\]

从列的角度看，上式等价于：

\[\mathbf{\bar{x}}_i = \frac{1}{2} \mathbf{\Lambda}^{-1} \mathbf{\bar{X}} \ (\mathbf{p}_i - \mathbf{g})\]

其中\(\mathbf{p}_i\)是\(\mathbf{P}\)矩阵的第\(i\)列，而\(\mathbf{g}=\text{diag}(\mathbf{\bar{x}}^T \mathbf{\bar{x}})\)。上式的意义在于当我们得到新的数据时无需重新计算特征值分解，而是可以利用已知的信息来计算新样本点的嵌入。

SMACOF and Variants

除了上面介绍的MDS相关算法外我们还可以从优化的角度来认识样本重建问题。此时我们希望嵌入后样本之间的距离矩阵可以尽可能地接近已知的距离矩阵，这样可以得到优化问题：

\[\min_{\mathbf{X}} \sum_{ij} (\mathbf{D}_{0ij} - \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2)^2\]

其中\(\mathbf{D}_0\)为已知的样本距离矩阵。遗憾的是上面定义的优化问题实际上是一个非凸优化问题，我们很难找到它的全局最优。尽管如此，我们仍然可以使用SMACOF算法(scaling by majorizing a complicated function)来进行求解。

首先对目标函数进行展开：

\[\begin{aligned} \sum_{ij} (\mathbf{D}_{0ij} - \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2)^2 &= \sum_{ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2^2 - 2 \sum_{ij} \mathbf{D}_{0ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2 + C \end{aligned}\]

对第一项进行展开可以得到：

\[\begin{aligned} \sum_{ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2^2 &= \sum_{ij} (\mathbf{G}_{ii} + \mathbf{G}_{jj} - 2 \mathbf{G}_{ij}) \\ &= 2n \cdot \text{tr} (\mathbf{G}) - 2 \sum_{ij} \mathbf{G}_{ij} \\ &= 2n \cdot \text{tr} (\mathbf{G}) - 2 \cdot \text{tr} (\mathbf{1}^T \mathbf{G} \mathbf{1}) \\ &= 2n \cdot \text{tr} (\mathbf{G}) - 2 \cdot \text{tr} (\mathbf{G} \mathbf{1} \mathbf{1}^T) \\ &= \text{tr} \big( \mathbf{G} (2n \ \mathbf{I} - 2 \mathbf{1} \mathbf{1}^T) \big) \\ &= \text{tr} (\mathbf{X}^T \mathbf{X} \mathbf{V}) \\ &= \text{tr} (\mathbf{X} \mathbf{V} \mathbf{X}^T) \end{aligned}\]

其中\(\mathbf{V} = 2n \ (\mathbf{I} - \mathbf{1} \mathbf{1}^T) = 2n \mathbf{J}\)。

接下来考虑第二项\(-2\sum_{ij} \mathbf{D}_{0ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2\)，实际上整个优化问题的非凸性都来源于它。对其进行展开可以得到：

\[\begin{aligned} -2\sum_{ij} \mathbf{D}_{0ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2 &= -2 \sum_{ij, \mathbf{x}_i \neq \mathbf{x}_j} \frac{\mathbf{D}_{0ij}}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2^2 \\ &= -2 \ \text{tr} (\mathbf{X} \mathbf{B} \mathbf{X}^T) \end{aligned}\] \[\mathbf{B}_{ij} (\mathbf{X}) = \begin{cases} -\frac{\mathbf{D}_{0ij}}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2} & \text{if}\ \mathbf{x}_i \neq \mathbf{x}_j \ \text{and}\ i \neq j \\ 0 & \text{if}\ \mathbf{x}_i = \mathbf{x}_j \ \text{and}\ i \neq j \\ -\sum_{j \neq i} \mathbf{B}_{ij} & \text{if}\ i=j \end{cases}\]

这样整个优化问题的非凸性都位于矩阵\(\mathbf{B}\)中。接下来我们定义一个新的目标函数\(\tau\)：

\[\tau (\mathbf{X}, \mathbf{Z}) = \text{tr} (\mathbf{X} \mathbf{V} \mathbf{X}^T) - 2 \ \text{tr} (\mathbf{X} \ \mathbf{B}(\mathbf{Z}) \ \mathbf{Z}^T) + C\]

显然原始的优化目标等于\(\tau (\mathbf{X}, \mathbf{X})\)，而且任意矩阵\(\mathbf{Z}\)都是原始目标函数的一个上界。这里简单证明一下，首先根据Cauchy-Schwarz不等式有：

\[(\mathbf{x}_i - \mathbf{x}_j) \cdot (\mathbf{z}_i - \mathbf{x}_j) \leq \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2 \ \Vert \mathbf{z}_i - \mathbf{z}_j \Vert_2\] \[\Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2 \geq \frac{(\mathbf{x}_i - \mathbf{x}_j) \cdot (\mathbf{z}_i - \mathbf{x}_j)}{\Vert \mathbf{z}_i - \mathbf{z}_j \Vert_2}\]

然后考虑\(\tau\)的第二项，

\[\begin{aligned} -2 \ \text{tr} (\mathbf{X} \ \mathbf{B}(\mathbf{X}) \ \mathbf{X}^T) &= -2\sum_{ij} \mathbf{D}_{0ij} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert_2 \\ &\leq -2\sum_{ij} \mathbf{D}_{0ij} \frac{(\mathbf{x}_i - \mathbf{x}_j) \cdot (\mathbf{z}_i - \mathbf{x}_j)}{\Vert \mathbf{z}_i - \mathbf{z}_j \Vert_2} \\ &= -2 \ \text{tr} (\mathbf{X} \ \mathbf{B}(\mathbf{Z}) \ \mathbf{Z}^T) \end{aligned}\]

由于\(\tau\)的剩余两项都与\(\mathbf{Z}\)无关，上式说明\(\tau (\mathbf{X}, \mathbf{Z})\)是\(\tau (\mathbf{X}, \mathbf{X})\)的一个上界：

\[\tau (\mathbf{X}, \mathbf{X}) \leq \tau (\mathbf{X}, \mathbf{Z})\]

这样我们就可以通过不断优化上界\(\tau (\mathbf{X}, \mathbf{Z})\)的方式来最小化原始目标函数。具体地，我们固定\(\tau\)的一个参数\(\mathbf{X}\)然后通过迭代的方式来最小化目标函数：

\[\mathbf{X}^{k+1} \leftarrow \arg \min_{\mathbf{X}} \ \tau (\mathbf{X}, \mathbf{X}^k)\]

可以证明这样的迭代方式是单调不增的：

\[\begin{aligned} f(\mathbf{X}^{k+1}) = \tau (\mathbf{X}^{k+1}, \mathbf{X}^{k+1}) &\leq \tau (\mathbf{X}^{k+1}, \mathbf{X}^k) \\ &= \tau (\mathbf{X}^k, \mathbf{X}^k) \\ &= f(\mathbf{X}^{k}) \end{aligned}\]

这种迭代式算法的另一个好处在于当我们固定了\(\tau\)的参数\(\mathbf{X}^k\)时，\(\tau\)关于\(\mathbf{X}\)是一个凸函数而且我们可以显式计算它的极值点：

\[\nabla_\mathbf{X} \ \tau(\mathbf{X}, \mathbf{X}^k) = 2 \mathbf{X} \mathbf{V} - 2 \mathbf{X}^k \mathbf{B} (\mathbf{X^k}) = 0\]

即迭代格式可以表示为：

\[\mathbf{X}^{k+1} \leftarrow \frac{1}{2n} \mathbf{X}^k \ \mathbf{B} (\mathbf{X^k}) \ \bigg( \mathbf{I}_{n \times n} - \frac{\mathbf{1} \mathbf{1}^T}{n} \bigg)\]