机器学习强基计划8-2：详细推导多维缩放MDS算法(附Python实现)

目录

• 0 写在前面
• 1 线性降维技术
• 2 多维缩放算法推导
• 3 Python实现

0 写在前面

机器学习强基计划聚焦深度和广度，加深对机器学习模型的理解与应用。“深”在详细推导算法模型背后的数学原理；“广”在分析多个机器学习模型：决策树、支持向量机、贝叶斯与马尔科夫决策、强化学习等。强基计划实现从理论到实践的全面覆盖，由本人亲自从底层编写、测试与文章配套的各个经典算法，不依赖于现有库，可以大大加深对算法的理解。

🚀详情：机器学习强基计划(附几十种经典模型源码)

---

1 线性降维技术

降维(dimension reduction)是缓解维数灾难的一个重要途径，因为样本数据往往以某种与学习任务密切相关的低维分布的形式，嵌入在高维空间内，如图所示。关于为什么要引入降维的详细解释请回顾：机器学习强基计划8-1：图解主成分分析PCA算法(附Python实现)

在降维技术中有一类经典的线性降维方法

给定 d d d维空间的样本集 X = [ x 1 x 2 ⋯ x m ] ∈ R d × m \boldsymbol{X}=\left[ \begin{matrix} \boldsymbol{x}_1& \boldsymbol{x}_2& \cdots& \boldsymbol{x}_m\\\end{matrix} \right] \in \mathbb{R} ^{d\times m} X=[x1​​x2​​⋯​xm​​]∈Rd×m，在线性降维技术中采用线性映射 W = [ w 1 w 2 ⋯ w d ′ ] ∈ R d × d ′ ( d ′ ⩽ d ) \boldsymbol{W}=\left[ \begin{matrix} \boldsymbol{w}_1& \boldsymbol{w}_2& \cdots& \boldsymbol{w}_{d'}\\\end{matrix} \right] \in \mathbb{R} ^{d\times d'}\left( d'\leqslant d \right) W=[w1​​w2​​⋯​wd′​​]∈Rd×d′(d′⩽d)将原始高维空间变换到低维空间

Z = W T X \boldsymbol{Z}=\boldsymbol{W}^T\boldsymbol{X} Z=WTX

其中 { w 1 , w 2 , ⋯   , w d ′ } \left\{ \boldsymbol{w}_1, \boldsymbol{w}_2, \cdots , \boldsymbol{w}_{d'} \right\} {w1​,w2​,⋯,wd′​}是 d ′ d' d′维空间的标准正交基在 d d d维空间的表示； Z ∈ R d ′ × m \boldsymbol{Z}\in \mathbb{R} ^{d'\times m} Z∈Rd′×m是样本集 X \boldsymbol{X} X在 d ′ d' d′维空间的表示。 z i = W T x i \boldsymbol{z}_i=\boldsymbol{W}^T\boldsymbol{x}_i zi​=WTxi​是高维样本 x i \boldsymbol{x}_i xi​与低维基向量做内积的降维过程，相反地， x i = ∑ j = 1 d ′ z i j w j \boldsymbol{x}_i=\sum\nolimits_{j=1}^{d'}{z_{ij}\boldsymbol{w}_j} xi​=∑j=1d′​zij​wj​是低维向高维扩展的重构过程。

不同的线性降维技术区别在于对降维映射 W \boldsymbol{W} W施加不同的约束，使降维后的空间符合某种特性或假设。

2 多维缩放算法推导

多维缩放算法(Multiple Dimensional Scaling, MDS)限制样本在经过降维映射 W \boldsymbol{W} W得到的低维空间中的欧式距离，等价于原始空间，即 ∥ z i − z j ∥ 2 = ∥ x i − x j ∥ 2 \left\| \boldsymbol{z}_i-\boldsymbol{z}_j \right\| _2=\left\| \boldsymbol{x}_i-\boldsymbol{x}_j \right\| _2 ∥zi​−zj​∥2​=∥xi​−xj​∥2​。

2.1 中心化约束

设样本 X = [ x 1 x 2 ⋯ x m ] \boldsymbol{X}=\left[ \begin{matrix} \boldsymbol{x}_1& \boldsymbol{x}_2& \cdots& \boldsymbol{x}_m\\\end{matrix} \right] X=[x1​​x2​​⋯​xm​​]在 d d d维空间中的距离矩阵为 D ∈ R m × m \boldsymbol{D}\in \mathbb{R} ^{m\times m} D∈Rm×m，其中的元素 d i j d_{ij} dij​表示样本 x i \boldsymbol{x}_i xi​与 x j \boldsymbol{x}_j xj​间的欧式距离，则根据MDS算法的设计，有

d i j 2 = ∥ x i − x j ∥ 2 2 = ∥ z i − z j ∥ 2 2 = ∥ z i ∥ 2 2 + ∥ z j ∥ 2 2 − 2 z i T z j d_{ij}^{2}=\left\| \boldsymbol{x}_i-\boldsymbol{x}_j \right\| _{2}^{2}=\left\| \boldsymbol{z}_i-\boldsymbol{z}_j \right\| _{2}^{2}=\left\| \boldsymbol{z}_i \right\| _{2}^{2}+\left\| \boldsymbol{z}_j \right\| _{2}^{2}-2\boldsymbol{z}_{i}^{T}\boldsymbol{z}_j dij2​=∥xi​−xj​∥22​=∥zi​−zj​∥22​=∥zi​∥22​+∥zj​∥22​−2ziT​zj​

必须指出，若对样本 Z \boldsymbol{Z} Z进行平移，则

∥ z i − z j ∥ 2 2 = ∥ ( z i + λ ) − ( z j + λ ) ∥ 2 2 \left\| \boldsymbol{z}_i-\boldsymbol{z}_j \right\| _{2}^{2}=\left\| \left( \boldsymbol{z}_i+\lambda \right) -\left( \boldsymbol{z}_j+\lambda \right) \right\| _{2}^{2} ∥zi​−zj​∥22​=∥(zi​+λ)−(zj​+λ)∥22​

换言之，降维样本有冗余自由度，使其平移样本簇都能满足与原始空间的距离等价性，因此，通常约束 ∑ i = 1 m z i = 0 \sum\nolimits_{i=1}^m{\boldsymbol{z}_i=\mathbf{0}} ∑i=1m​zi​=0，既便于推导算法，又不破坏假设条件。

2.2 内积矩阵特征值分解

设降维空间的内积矩阵 B = Z T Z ∈ R m × m \boldsymbol{B}=\boldsymbol{Z}^T\boldsymbol{Z}\in \mathbb{R} ^{m\times m} B=ZTZ∈Rm×m，其中的元素 b i j = z i T z j b_{ij}=\boldsymbol{z}_{i}^{T}\boldsymbol{z}_j bij​=ziT​zj​，则

d i j 2 = b i i + b j j − 2 b i j ⇔ b i j = 1 2 ( b i i + b j j − d i j 2 ) d_{ij}^{2}=b_{ii}+b_{jj}-2b_{ij}\Leftrightarrow b_{ij}=\frac{1}{2}\left( b_{ii}+b_{jj}-d_{ij}^{2} \right) dij2​=bii​+bjj​−2bij​⇔bij​=21​(bii​+bjj​−dij2​)

考虑到 Z \boldsymbol{Z} Z被中心化，因此 ∑ i = 1 m b i j = ∑ j = 1 m b i j = ( ∑ i = 1 m z i T ) z j = 0 \sum\nolimits_{i=1}^m{b_{ij}}=\sum\nolimits_{j=1}^m{b_{ij}}=\left( \sum\nolimits_{i=1}^m{\boldsymbol{z}_{i}^{T}} \right) \boldsymbol{z}_j=\mathbf{0} ∑i=1m​bij​=∑j=1m​bij​=(∑i=1m​ziT​)zj​=0。现建立 D \boldsymbol{D} D与 B \boldsymbol{B} B的联系

{ ∑ i = 1 m d i j 2 = t r ( B ) + m b j j ∑ j = 1 m d i j 2 = t r ( B ) + m b i i ∑ i = 1 m ∑ j = 1 m d i j 2 = 2 m t r ( B ) \begin{cases} \sum_{i=1}^m{d_{ij}^{2}}=\mathrm{tr}\left( \boldsymbol{B} \right) +mb_{jj}\\ \sum_{j=1}^m{d_{ij}^{2}}=\mathrm{tr}\left( \boldsymbol{B} \right) +mb_{ii}\\ \sum_{i=1}^m{\sum_{j=1}^m{d_{ij}^{2}}}=2m\mathrm{tr}\left( \boldsymbol{B} \right)\\\end{cases} ⎩ ⎨ ⎧​∑i=1m​dij2​=tr(B)+mbjj​∑j=1m​dij2​=tr(B)+mbii​∑i=1m​∑j=1m​dij2​=2mtr(B)​

联立上式可得

b i j = 1 2 ( b i i + b j j − d i j 2 ) = 1 2 ( 1 m ( ∑ i = 1 m d i j 2 − t r ( B ) ) + 1 m ( ∑ j = 1 m d i j 2 − t r ( B ) ) − d i j 2 ) = 1 2 ( 1 m ∑ i = 1 m d i j 2 + 1 m ∑ j = 1 m d i j 2 − 1 m 2 ∑ i = 1 m ∑ j = 1 m d i j 2 − d i j 2 ) \begin{aligned}b_{ij}&=\frac{1}{2}\left( b_{ii}+b_{jj}-d_{ij}^{2} \right) \\&=\frac{1}{2}\left( \frac{1}{m}\left( \sum_{i=1}^m{d_{ij}^{2}}-\mathrm{tr}\left( \boldsymbol{B} \right) \right) +\frac{1}{m}\left( \sum_{j=1}^m{d_{ij}^{2}}-\mathrm{tr}\left( \boldsymbol{B} \right) \right) -d_{ij}^{2} \right) \\&=\frac{1}{2}\left( \frac{1}{m}\sum_{i=1}^m{d_{ij}^{2}}+\frac{1}{m}\sum_{j=1}^m{d_{ij}^{2}}-\frac{1}{m^2}\sum_{i=1}^m{\sum_{j=1}^m{d_{ij}^{2}}}-d_{ij}^{2} \right)\end{aligned} bij​​=21​(bii​+bjj​−dij2​)=21​(m1​(i=1∑m​dij2​−tr(B))+m1​(j=1∑m​dij2​−tr(B))−dij2​)=21​(m1​i=1∑m​dij2​+m1​j=1∑m​dij2​−m21​i=1∑m​j=1∑m​dij2​−dij2​)​

由此即可根据已知的距离矩阵 D \boldsymbol{D} D求出在距离不变形约束下的低维内积矩阵 B \boldsymbol{B} B。现对 B \boldsymbol{B} B进行特征值分解可得

B = U Λ U T \boldsymbol{B}=\boldsymbol{U\varLambda U}^T B=UΛUT
其中 Λ = d i a g ( λ 1 , λ 2 , ⋯   , λ m ) \boldsymbol{\varLambda }=\mathrm{diag}\left( \lambda _1,\lambda _2,\cdots ,\lambda _m \right) Λ=diag(λ1​,λ2​,⋯,λm​)为特征值构成的对角矩阵，假设 λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ λ m \lambda _1\geqslant \lambda _2\geqslant \cdots \geqslant \lambda _m λ1​⩾λ2​⩾⋯⩾λm​； U \boldsymbol{U} U是特征值对应的特征向量组成的正交矩阵。

实际应用中为有效降维，往往仅需降维后与原始空间样本距离近似而非严格相等，因此通常取 d ′ ≪ d d'\ll d d′≪d个最大特征值构成 Λ ~ ∈ R d ′ × d ′ \boldsymbol{\tilde{\varLambda}}\in \mathbb{R} ^{d'\times d'} Λ~∈Rd′×d′、 U ~ ∈ R m × d ′ \boldsymbol{\tilde{U}}\in \mathbb{R} ^{m\times d'} U~∈Rm×d′。考虑到

B = U ~ Λ ~ 1 / 2 Λ ~ 1 / 2 U ~ T = ( U ~ Λ ~ 1 / 2 ) ( U ~ Λ ~ 1 / 2 ) T \boldsymbol{B}=\boldsymbol{\tilde{U}\tilde{\varLambda}}^{{{1}/{2}}}\boldsymbol{\tilde{\varLambda}}^{{{1}/{2}}}\boldsymbol{\tilde{U}}^T=\left( \boldsymbol{\tilde{U}\tilde{\varLambda}}^{{{1}/{2}}} \right) \left( \boldsymbol{\tilde{U}\tilde{\varLambda}}^{{{1}/{2}}} \right) ^T B=U~Λ~1/2Λ~1/2U~T=(U~Λ~1/2)(U~Λ~1/2)T

结合 B = Z T Z \boldsymbol{B}=\boldsymbol{Z}^T\boldsymbol{Z} B=ZTZ可知

Z = Λ ~ 1 / 2 U ~ T ∈ R d ′ × m \boldsymbol{Z}=\boldsymbol{\tilde{\varLambda}}^{{{1}/{2}}}\boldsymbol{\tilde{U}}^T\in \mathbb{R} ^{d'\times m} Z=Λ~1/2U~T∈Rd′×m

降维完成。

3 Python实现

3.1 算法流程

3.2 可视化

核心代码如下

def mds(self, D, outDim):
	# 计算内积矩阵
	B = self.__calB(D)
	# 特征值分解
	eigVal, eigVec = np.linalg.eig(B)
	# 获取最大的d'个特征值对应的索引, np.argsort是按从小到大排序, 所以对特征值取负号
	index = np.argsort(-eigVal)[0:outDim]
	eigVal_ = eigVal[index].real
	eigVec_ = eigVec[:, index]
	# 计算低维样本
	Z = np.dot(np.diag(eigVal_)**0.5, eigVec_.T)
	return Z

以鸢尾花数据集为例执行降维，效果如下

---

🔥 更多精彩专栏：

• 《ROS从入门到精通》
• 《Pytorch深度学习实战》
• 《机器学习强基计划》
• 《运动规划实战精讲》
• …