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)0
  • 同一性d(x,y)=0x=y
  • 对称性(symmetry)d(x,y)=d(y,x)
  • 三角不等式(triangle inequality)d(x,z)d(x,y)+d(y,z)

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

Embedding a Metric Space

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

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

但对于有限度量空间可以证明存在等距映射将样本嵌入到某个由无穷范数定义的欧式空间l(Rn)上。

这里我们简单证明一下上面的结论。假设存在有限度量空间(M,d),其中M={a1,a2,...,an}。接下来定义映射ϕ:MRn为:

ϕ(ak)=(d(a1,ak),d(a2,ak),...,d(an,ak))

不难证明ϕ是一个等距映射。首先考虑映射后ϕ(ai)ϕ(aj)之间的距离:

ϕ(ai)ϕ(aj)=maxk|d(ai,ak)d(aj,ak)|

对于右边有:

maxk|d(ai,ak)d(aj,ak)||d(ai,ai)d(aj,ai)|=d(ai,aj)

同时根据三角不等式,对于任意1kn满足:

d(ai,ak)d(ai,aj)+d(aj,ak)d(ai,ak)d(aj,ak)d(ai,aj) d(aj,ak)d(aj,ai)+d(ai,ak)d(aj,ak)d(ai,ak)d(ai,aj)

|d(ai,ak)d(aj,ak)|d(ai,aj) maxk|d(ai,ak)d(aj,ak)|d(ai,aj)

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

d(ai,aj)maxk|d(ai,ak)d(aj,ak)|d(ai,aj) maxk|d(ai,ak)d(aj,ak)|=d(ai,aj)

即原始度量空间M可以通过等距映射ϕ变换到度量空间l(Rn)上。

Approximate Embedding

Approximate Isometry

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

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

Fréchet Embedding

具体地,Fréchet embedding就是这样的一种映射方法。我们只需要在原始空间中选择一系列子集S1,...,SrM,然后构造出Rr中的向量即可。

Bourgain’s Theorem

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

Embedding Metrics into Euclidean Space

Euclidean Problem

接下来考虑如何从度量来恢复几何的问题。假设在欧式空间Rm中存在样本{x1,...,xn},我们已知任意两个样本之间的距离为Pij=xixj22,我们的目标是重建Rm中样本的坐标。

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

Gram Matrix

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

G=XTX

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

Gij=xixj

Classical MDS

利用Gram矩阵就可以把样本的重建问题分解为2步:首先利用P矩阵恢复Gram矩阵,然后再利用Gram矩阵来重建样本坐标。

Recover Gram Matrix

假设样本的中心为原点,即ixi=0,利用欧式空间下内积的定义有:

Pij=xixj22=xi22+xj222xixj=Gii+Gjj2Gij

注意到GiiGjj是Gram矩阵的对角元素,我们可以把P矩阵使用Gram矩阵来表示:

P=2G+diag(G) 1T+1 diag(G)T

根据中心化条件ixi=0,可以得到:

P 1=ixi=0 G 1=[x1x1+x1x2+x1xnxnx1+xnx2+xnxn]=[x1(ixi)xn(ixi)]=0

因此我们可以定义一个辅助矩阵J=I1n11T,此时

12JTPJ=12(I1n11T)T (2G+diag(G) 1T+1 diag(G)T) (I1n11T)=12(2G+diag(G) 1T+1 diag(G)T+2n11TG1n11Tdiag(G) 1T1n11T1 diag(G)T)(I1n11T)=12(2G+diag(G) 1T1ntr(G)11T)(I1n11T)=12(2G+diag(G) 1T1ntr(G)11T1ndiag(G)1T11T+1n2tr(G)11T11T)=G

这样我们就从P矩阵恢复了Gram矩阵。

Recover X from Gram Matrix

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

G=Q Λ QT

其中矩阵Q不一定是一个方阵,且QTQ=I。这样就可以重建出样本的坐标:

X¯=Λ QT

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

Landmark MDS

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

X¯P=(Λ QT) (2QΛQT+diag(QΛQT) 1T+1 diag(QΛQT)T)=2Λ3/2QT+X¯ diag(QΛQT) 1T X¯(Pdiag(QΛQT) 1T)=2Λ3/2QT=2ΛX¯

整理后得到:

X¯=12Λ1X¯(Pdiag(QΛQT) 1T)

从列的角度看,上式等价于:

x¯i=12Λ1X¯ (pig)

其中piP矩阵的第i列,而g=diag(x¯Tx¯)。上式的意义在于当我们得到新的数据时无需重新计算特征值分解,而是可以利用已知的信息来计算新样本点的嵌入。

SMACOF and Variants

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

minXij(D0ijxixj2)2

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

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

ij(D0ijxixj2)2=ijxixj222ijD0ijxixj2+C

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

ijxixj22=ij(Gii+Gjj2Gij)=2ntr(G)2ijGij=2ntr(G)2tr(1TG1)=2ntr(G)2tr(G11T)=tr(G(2n I211T))=tr(XTXV)=tr(XVXT)

其中V=2n (I11T)=2nJ

接下来考虑第二项2ijD0ijxixj2,实际上整个优化问题的非凸性都来源于它。对其进行展开可以得到:

2ijD0ijxixj2=2ij,xixjD0ijxixj2xixj22=2 tr(XBXT) Bij(X)={D0ijxixj2if xixj and ij0if xi=xj and ijjiBijif i=j

这样整个优化问题的非凸性都位于矩阵B中。接下来我们定义一个新的目标函数τ

τ(X,Z)=tr(XVXT)2 tr(X B(Z) ZT)+C

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

(xixj)(zixj)xixj2 zizj2 xixj2(xixj)(zixj)zizj2

然后考虑τ的第二项,

2 tr(X B(X) XT)=2ijD0ijxixj22ijD0ij(xixj)(zixj)zizj2=2 tr(X B(Z) ZT)

由于τ的剩余两项都与Z无关,上式说明τ(X,Z)τ(X,X)的一个上界:

τ(X,X)τ(X,Z)

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

Xk+1argminX τ(X,Xk)

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

f(Xk+1)=τ(Xk+1,Xk+1)τ(Xk+1,Xk)=τ(Xk,Xk)=f(Xk)

这种迭代式算法的另一个好处在于当我们固定了τ的参数Xk时,τ关于X是一个凸函数而且我们可以显式计算它的极值点:

X τ(X,Xk)=2XV2XkB(Xk)=0

即迭代格式可以表示为:

Xk+112nXk B(Xk) (In×n11Tn)

因此SMACOF算法的本质是优化原始目标的一个凸的上界来逐步最小化目标函数。尽管这种方法理论上并不能保证收敛到全局最优值,但大量的实践表明它仍然可以收敛到一个足够好的解上。

SMACOF算法有很多的应用,比如说可以使用它来可视化图或者对网格进行变形等。

限于课时有限本节课只能对multidimensional scaling这类方法进行简单的介绍,更多的细节可以参考专业教材。

Reference