第3章
CHAPTER 3


Python数据降维




伴随ICT(通信与信息技术)和互联网技术的不断发展，人们收集和获得数据的能力越来越强。而这些数据已呈现出维数高、规模大和结构复杂等特点。
人们想利用这些大数据(维数大、规模大、复杂大)，挖掘其中有意义的知识和内容以指导实际生产和具体应用，数据的降维就显得尤为重要了。数据降维又称为维数约简。顾名思义，就是降低数据的维数。为什么要降低数据的维数？如何有效地降低数据的维数？由此问题引发了广泛的研究和应用。
数据降维，一方面可以解决“维数灾难”，缓解“信息丰富、知识贫乏”现状，降低复杂度； 另一方面可以更好地认识和理解数据。
截止到目前，数据降维的方法很多。从不同的角度入手可以有着不同的分类，主要分类方法有： 根据数据的特性可以划分为线性降维和非线性降维，根据是否考虑和利用数据的监督信息可以划分为无监督降维、有监督降维和半监督降维，根据保持数据的结构可以划分为全局保持降维、局部保持降维和全局与局部保持一致降维等。
总之，数据降维意义重大，数据降维方法众多，很多时候需要根据特定问题选用合适的数据降维方法。数据降维是机器学习领域中非常重要的内容。
3.1维度灾难与降维
1.  维度灾难
维度灾难(curse of dimensionality)用来描述当(数学)空间维度增加时，分析和组织高维空间(通常有成百上千维)，因体积指数增加而遇到各种问题场景。在机器学习中，维度灾难常指以下问题： 
在高维情况下，数据样本稀疏。
例如，k近邻法的讨论中经常涉及维度灾难，是因为k近邻法基于一个重要的基本假设： 任意样本附近任意小的距离内总能找到一个训练样本，即训练样本的采样密度足够大，也称为“密采样”，才能保证分类性能； 当特征维度很大时，满足密采样的样本数量会呈指数级增长，大到几乎无法达到。 
在高维情况下，涉及距离、内积的计算变得困难。
其实，不仅是k近邻，其他机器学习算法几乎都会遇到维度灾难的问题。
2.  降维
缓解维度灾难的一个重要途径就是降维。
1) 为什么能够进行降维
这是因为很多时候，数据是高维的，但是与学习任务(分类、回归等)密切相关的仅是某个低维分布，即高维空间中的某个低维难嵌入。因此，很多情况下，高维空间中的样本点，在低维嵌入子空间中更容易学习。
2) 线性降维
一般来说，想获得低维子空间，最简单的方法是对原始高维空间进行线性变换： 
给定d维空间中的样本X=(x1,x2,…,xm)∈Rd×m，变换之后得到d′≤d维空间中的样本

Z=WTX

其中，W∈Rd×d′是变换矩阵，Z∈Rd′×m是样本在新空间中的表达。
变换矩阵W可视为d′个d维基向量，新空间中的属性是原空间属性的线性组合，基于线性变换来进行降维的方法称为线性维方法，都符合上面的式子，主要区别在于对低维子空间的性质有所不同，相当于对W施加了不同的约束。
3) 降维效果的评估
通常通过比较降维前后学习器性能，若性能有提高，则认为降维起到了作用。针对已经降到二维或者三维的情况，可以利用可视化技术直观地判断降维效果。
3.2主成分分析
主成分分析(Principal Component Analysis,PCA)是一种最常用的无监督降维方法，通过降维技术把多个变量化为少数几个主成分(综合变量)的统计分析方法。这些主成分能够反映原始变量的绝大部分信息，它们通常表示为原始变量的某种线性组合。
3.2.1PCA原理
为了便于维度变换，有如下假设。
 假设样本数据是n维的。
 假设原始坐标系为： 由标准正交基向量{i→1,i→2,…,i→n}张成的空间，其中‖i→s‖=1； i→s·i→t=0，s≠t。
 假设经过线性变换后的新坐标系为： 由标准正交基向量j→1,j→2,…,j→n张成的空间，其中‖j→s‖=1； j→s·j→t=0，s≠t。
根据定义，有： 

j→s=(i→1,i→2,…,i→n)j→s·i→1
︙
j→s·i→n，s=1,2,…,n

记ws=(j→s·i→1,j→s·i→2,…,j→s·i→n)T(它是一个列向量，但是为了与基向量做区分，这里没有给出向量的箭头符号)，其各分量就是基向量j→s在原始坐标系{i→1,i→2,…,i→n}中的投影。即： j→s=(i→1,i→2,…,i→n)·ws。根据标准正交基的性质，有： 
 ‖ws‖=1，s=1,2,…,n； 
 ws·wt=0，s≠t。
根据定义，有： 

(j→1,j→2,…,j→n)=(i→1,i→2,…,i→n)(w1,w2,…,wn)

令坐标变换矩阵W为： 

W=(w1,w2,…,wn)=j→1·i→1j→2·i→1…j→n·i→1
j→1·i→2j→2·i→2…j→n·i→2
︙︙︙
j→1·i→nj→2·i→n…j→n·i→n

则有

(j→1,j→2,…,j→n)=(i→1,i→2,…,i→n)W

W的第s列就是j→s在原始坐标系{i→1,i→2,…,i→n}中的投影，且有W=WT，WWT=I(即它的逆矩阵就是它的转置)。
假设样本点x→i在原始坐标系中的表示为： 

x→i=(i→1,i→2,…,i→n)x(1)i
x(2)i
︙
x(n)i

令xi=(x(1)i,x(2)i，…,x(n)i)T，则x→i=(i→1,i→2,…,i→n)xi。
假设样本点x→i在新坐标系中的表示为： 

x→i=(j→1,j→2,…,j→n)z(1)i
z(2)i
︙
z(n)i

令zi=(z(1)i,z(2)i，…,z(n)i)T，则x→i=(j→1,j→2,…,j→n)zi。根据x→i=x→i，则有： 

(j→1,j→2,…,j→n)zi=(i→1,i→2,…,i→n)Wzi=(i→1,i→2,…,i→n)xi

于是有： 

zi=W-1xi=WTxi

丢弃其中的部分坐标，将维度降低到d<n，则样本点x→i在低维坐标系中的坐标为z′i=(z(1)i,z(2)i，…,z(n)i)T。现在的问题是： 最好丢弃哪些坐标？想法是： 基于降低之后的坐标重构样本时，尽量要与原始样本相近。
如果基于降维后的坐标z′i来重构x→i： 

xi=(j→1,j→2,…,j→d)z(1)i
z(2)i
︙
z(d)i=(i→1,i→2,…,i→n)(w1,w2,…,wd)z(1)i
z(2)i
︙
z(d)i
=(i→1,i→2,…,i→n)(w1,w2,…,wd)wT1·xi
wT2·xi
︙
wTd·xi
=(i→1,i→2,…,i→n)(w1,w2,…,wd)wT1
wT2
︙
wTd·xi

令Wd=(w1,w2,…,wd)，即它是坐标变换矩阵W的前d列，则： 

xi=(i→1,i→2,…,i→n)WdWTdxi

考虑整个训练集，原样本点x→i和基于投影重构的样本点xi之间的距离为(即所有重构的样本点与原样本点的整体误差)： 

∑Ni=1‖xi-x→i‖22=∑Ni=1xi-WdWTdxi‖22

考虑： 

WdWTdxi=(w1,w2,…,wd)wT1
wT2
︙
wTdxi=∑ds=1ws(wTsxi)

由于wTsxi是标量，所以有： 

WdWTdxi=∑ds=1(wTsxi)ws

由于wTsxi是标量，所以它的转置等于它本身，所以有： 

WdWTdxi=∑ds=1(xTiws)ws

于是有： 

∑Ni=1‖xi-x→i‖22=∑Ni=1‖xi-WdWTdxi‖22=∑Ni=1‖xi-∑ds=1(xTiws)ws‖22

定义矩阵X=(x1,x2,…,xN)，即矩阵X的第i列就是xi。即可以证明： 

‖XT-XTWdWTd‖2F=∑Ni=1‖xi-∑ds=1(xTiws)ws‖22

其中，‖·‖F为矩阵的Frobenius范数（简称F范数）。接下来的证明过程中，要用到矩阵的F范数和矩阵的迹的性质： 
(1) 矩阵A的F范数定义为： ‖A‖F=∑i∑ja2ij，即矩阵所有元素的平方和的开方。F范数的性质有： 
 ‖A‖F=‖AT‖F。
 ‖A‖F=tr(ATA)，tr为矩阵的迹。
(2) 对于方阵，矩阵的迹定义为： tr(A)=∑iaii，即矩阵对角线元素之和。矩阵的迹的性质有： 
 tr(A)=tr(AT)。
 tr(A±B)=tr(A)±tr(B)。
 如果A为m×n阶矩阵，B为n×m阶矩阵，则tr(AB)=tr(BA)。
 矩阵的迹等于矩阵的特征值之和： tr(A)=λ1+λ2+…+λn。
 对任何正整数k有： tr(Ak)=λk1+λk2+…+λkn。
要求解最优化问题： 

W*d=argminWd‖xi-x→i‖22
=argminWd‖XT-XTWdWTd‖2F
=argminWdtr［(XT-XTWdWTd)T(XT-XTWdWTd)］
=argminWdtr［(X-WdWTdX)(XT-XTWdWTd)］
=argminWdtr［XXT-XXTWdWTd-WdWTdXXT+WdWTdXXTWdWTd］
=argminWdtr［tr(XXT)-tr(XXTWdWTd)-tr(WdWTdXXT)+tr(WdWTdXXTWdWTd)］

因为矩阵及其转置的迹相等，因此tr(XXTWdWTd)=tr(WdWTdXXT)。由于可以在tr(·)中调整矩阵的顺序，则tr(WdWTdXXTWdWTd)=tr(XXTWdWTdWdWTd)。
考虑到： 

WdWTd=wT1
wT2
︙
wTd(w1,w2,…,wd)=Id×d

代入上式有： 

tr(WdWTdXXTWdWTd)=tr(XXTWdWTd)

于是： 

W*d=argminWd［tr(XXT)-2tr(XXTWdWTd)+tr(XXTWdWTd)］
=argminWd［tr(XXT)-tr(XXTWdWTd)］

由于tr(XXT)与Wd无关，因此，

W*d=argminWd-tr(XXTWdWTd)=argminWdtr(XXTWdWTd)

调整矩阵顺序，有： 

W*d=argminWdtr(WTdXXTWd)

该最优化问题的求解就是求解XXT的特征值。因此只需要对矩阵XXT(也称为样本的协方差矩阵，它是一个n阶方阵)进行特征值分解，将求得的特征值排序： λ1≥λ2≥…≥λn，然后取前d个特征值对应的特征向量构成W=(w1,w2,…,wd)。
3.2.2PCA算法
下面介绍PCA算法。
(1) 输入： 样本集D={x→1,x→2,…,x→N}； 低维空间维数d。
(2) 输出： 投影矩阵W=(w→1,w→2,…,w→d)。
(3) 算法步骤表现在： 
 对所有样本进行中心化操作

x→i←x→i-1N∑Nj=1x→j

 计算样本的协方差矩阵XXT； 
 对协方差矩阵XXT做特征值分解； 
 取最大的d个特征值对应的特征向量w→1,w→2,…,w→d，构造投影矩阵W=(w→1,w→2,…,w→d)。
通常低维空间维数d的选取有两种方法： 
 通过交叉验证法选取较好的d(在降维后的学习器的性能比较好)。
 从算法原理的角度设置一个阈值，例如t=95%，然后选取使得下式成立的最小的d的值： 

∑di=1λi∑ni=1λi≥t

其中，λi从大到小排列。
3.2.3PCA降维的两个准则
PCA降维的准则有以下两个： 
 最近重构性——就是前面介绍的样本集中所有点，重构后的点距离原来的点的误差之和最小。
 最大可分性——样本点在低维空间的投影尽可能分开。
可以证明，最近重构性就等价于最大可分性。证明如下： 对于样本点x→i，它在降维后空间中的投影是z→i。根据： 

xi=(w→1,w→2,…,w→d)z(1)i
z(2)i
︙
z(d)i=Wz→i

由投影矩阵的性质，以及xi与x→i的关系，有z→i=WTx→i。
由于样本数据进行了中心化： 即∑ix→i=(0,0,…,0)T，因此投影后，样本点的方差为： 

∑Ni=1WTx→ix→TiW

令X=(x→1,x→2,…,x→N)为n×N维矩阵，于是根据样本点的方差最大，优化目标可为： 

maxWtr(WTXXTW)

s.t.WTW=I

这就是前面最后重构性推导的结果。
【例31】通过Python的sklearn库来实现鸢尾花数据进行降维，数据本身是四维的降维后变成二维，可以在平面中画出样本点的分布。样本数据结构如图31所示。 


图31鸢尾花数据


其中样本总数为150，鸢尾花的类别有3种，分别标记为0、1、2。

import matplotlib.pyplot as plt #加载matplotlib用于数据的可视化

from sklearn.decomposition import PCA #加载PCA算法包

from sklearn.datasets import load_iris



data=load_iris()

y=data.target

x=data.data

pca=PCA(n_components=2) #加载PCA算法，设置降维后主成分数目为2

reduced_x=pca.fit_transform(x)#对样本进行降维



red_x,red_y=［］,［］

blue_x,blue_y=［］,［］

green_x,green_y=［］,［］



for i in range(len(reduced_x)):

if y［i］ ==0:

red_x.append(reduced_x［i］［0］)

red_y.append(reduced_x［i］［1］)

elif y［i］==1:

blue_x.append(reduced_x［i］［0］)

blue_y.append(reduced_x［i］［1］)

else:

green_x.append(reduced_x［i］［0］)

green_y.append(reduced_x［i］［1］)

#可视化

plt.scatter(red_x,red_y,c='r',marker='x')

plt.scatter(blue_x,blue_y,c='b',marker='D')

plt.scatter(green_x,green_y,c='g',marker='.')

plt.show()

运行程序，得到如图32所示的效果。


图32主成分降维效果


【例32】利用PCA对给定的数据data2.txt进行降维处理。
其实现步骤为： 
(1) 首先引入numpy，由于测试中用到了pandas和matplotlib，所以这里一并加载。

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

(2) 定义一个均值函数。

#计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征

def meanX(dataX):

return np.mean(dataX,axis=0)#axis=0表示按照列来求均值，如果输入list，则axis=1

(3) 编写pca方法，具体解释参考注释。

"""

参数： 
 XMat——传入的是一个numpy的矩阵格式，行表示样本数，列表示特征。
 k——表示取前k个特征值对应的特征向量。
返回值： 
 finalData——指的是返回的低维矩阵，对应于输入reconData。
 reconData——对应的是移动坐标轴后的矩阵。

"""

def pca(XMat, k):

average = meanX(XMat) 

m, n = np.shape(XMat)

data_adjust = ［］

avgs = np.tile(average, (m, 1))

data_adjust = XMat - avgs

covX = np.cov(data_adjust.T) #计算协方差矩阵

featValue, featVec=np.linalg.eig(covX)#求解协方差矩阵的特征值和特征向量

index = np.argsort(-featValue) #按照featValue进行从大到小排序

finalData = ［］

if k > n:

print("k must lower than feature number")

return

else:

#注意特征向量时列向量，而numpy的二维矩阵(数组)a［m］［n］中，a［1］表示第1行值

selectVec = np.matrix(featVec.T［index［:k］］) #所以这里需要进行转置

finalData = data_adjust * selectVec.T 

reconData = (finalData * selectVec) + average

return finalData, reconData

(4) 编写一个加载数据集的函数。

#输入文件的每行数据都以＼t隔开

def loaddata(datafile):

return np.array(pd.read_csv(datafile,sep="＼t",header=-1)).astype(np.float)

(5) 可视化结果。将维数k指定为2，所以可以使用下面的函数将其绘制出来： 

def plotBestFit(data1, data2):

dataArr1 = np.array(data1)

dataArr2 = np.array(data2)



m = np.shape(dataArr1)［0］

axis_x1 = ［］

axis_y1 = ［］

axis_x2 = ［］

axis_y2 = ［］

for i in range(m):

axis_x1.append(dataArr1［i,0］)

axis_y1.append(dataArr1［i,1］)

axis_x2.append(dataArr2［i,0］) 

axis_y2.append(dataArr2［i,1］)

fig = plt.figure()

ax = fig.add_subplot(111)

ax.scatter(axis_x1, axis_y1, s=50, c='red', marker='s')

ax.scatter(axis_x2, axis_y2, s=50, c='blue')

plt.xlabel('x1'); plt.ylabel('x2');

plt.savefig("outfile.png")

plt.show()

(6) 测试方法。将测试方法写入main函数中，然后直接执行main函数即可： 

#根据数据集data.txt

def main():

datafile = "data2.txt"

XMat = loaddata(datafile)

k = 2

return pca(XMat, k)

if __name__ == "__main__":

finalData, reconMat = main()

plotBestFit(finalData, reconMat)

运行程序，效果如图33所示。


图33PCA降维处理效果

3.3SVD降维
奇异值分解(SVD)： 设X为n×N阶矩阵，且rank(X)=r，则n阶正交矩阵V和N阶正交矩阵U，使得： 

VTXU=Σ0
00n×N

其中，

Σ=σ10…0
0σ2…0
︙︙︙
00…σr

其中，σ1≥σ2≥…≥σr>0。
根据正交矩阵的性质，VVT=I，UUT=I，有： 

X=VΣ0
00n×NUTXT=UΣ0
00n×NVT

则有XXT=VMVT，其中M是n阶对角矩阵： 

M=Σ0
00n×NΣ0
00N×n=λ100…0
0λ20…0
︙︙︙︙
000…λnn×n

λi=σ2ii=1,2,…,r

λi=0i=r+1,r+2,…,n

于是有： XXTV=VM。根据M是对角矩阵的性质，有VM=MV，则有： 

XXTV=MV

则λi（i=1,2,…,r）就是XXT的特征值，其对应的特征向量组成正交矩阵V。因此SVD奇异值分解等价于PCA主成分分析，核心都是求解XXT的特征值以及对应的特征向量。
【例33】利用SVD对给定的数据进行降维处理。

import numpy as np 

class CSVD(object):

'''

实现SVD分解降维应用示例的Python代码

'''

def __init__(self, data):

self.data = data #用户数据

self.S = ［］#用户数据矩阵的奇异值序列

self.U = ［］#svd后的单位正交向量

self.VT = ［］#svd后的单位正交向量

self.k = 0 #满足self.p的最小k值(k表示奇异值的个数)

self.SD = ［］ #对角矩阵，对角线上元素是奇异值

def _svd(self):

'''

用户数据矩阵的SVD奇异值分解

'''

self.U, self.S, self.VT = np.linalg.svd(self.data)

return self.U, self.S, self.VT

def _calc_k(self, percentge):

'''确定k值：前k个奇异值的平方和占比 >=percentage, 求满足此条件的最小k值

:param percentage, 奇异值平方和的占比的阈值

:return 满足阈值percentage的最小k值

'''

self.k = 0

#用户数据矩阵的奇异值序列的平方和

total = sum(np.square(self.S))

svss = 0 #奇异值平方和

for i in range(np.shape(self.S)［0］):

svss += np.square(self.S［i］)

if (svss/total) >= percentge:

self.k = i+1

break

return self.k



def _buildSD(self, k):

'''构建由奇异值组成的对角矩阵

:param k,根据奇异值开放和的占比阈值计算出来的k值

:return 由k个前奇异值组成的对角矩阵

'''

#方法1：用数组乘方法

self.SD = np.eye(self.k) * self.S［:self.k］ 

#方法2：用自定义方法

e = np.eye(self.k)

for i in range(self.k):

e［i,i］ = self.S［i］ 

return self.SD



def DimReduce(self, percentage):

'''

SVD降维

:param percentage, 奇异值开方和的占比阈值

:return 降维后的用户数据矩阵

'''

#Step1:svd奇异值分解

self._svd()

#Step2:计算k值

self._calc_k(percentage)

print('＼n按照奇异值开方和占比阈值percentage=%d, 求得降维的k=%d'%(percentage, self.k))

#Step3:构建由奇异值组成的对角矩阵

self._buildSD(self.k)

k,U,SD,VT = self.k,self.U, self.SD, self.VT

#Step4:按照svd分解公式对用户数据矩阵进行降维，得到降维压缩后的数据矩阵

a = U［:len(U), :k］

b = np.dot(SD, VT［:k, :len(VT)］)

newData = np.dot(a,b)

return newData



def CSVD_manual():

#训练数据集，用户对商品的评分矩阵，行为多个用户对单个商品的评分，列为用户对每个

#商品的评分

data = np.array(［［5, 5, 0, 5］,

［5, 0, 3, 4］,

［3, 4, 0, 3］,

［0, 0, 5, 3］,

［5, 4, 4, 5］,

［5, 4, 5, 5］］)

percentage = 0.9

svdor = CSVD(data)

ret = svdor.DimReduce(percentage)

print('====================================================')

print('原始用户数据矩阵:＼n', data)

print('降维后的数据矩阵:＼n', ret)

print('====================================================')

if __name__=='__main__':

CSVD_manual()

运行程序，输出如下：

按照奇异值开方和占比阈值percentage=0, 求得降维的k=2

====================================================

原始用户数据矩阵:

［［5 5 0 5］

［5 0 3 4］

［3 4 0 3］

［0 0 5 3］

［5 4 4 5］

［5 4 5 5］］

降维后的数据矩阵:

［［ 5.288493595.162728120.214912374.45908018］

［ 3.276809941.902085433.740019723.80580978］

［ 3.532418273.54790444 -0.133168882.89840405］

［ 1.14752376 -0.641713684.947235862.3845504 ］

［ 5.072687063.663995353.788689655.31300375］

［ 5.108565953.401879054.6166049 5.58222363］］

====================================================

3.4核主成分分析降维
PCA方法假设从高维空间到低维空间的函数映射是线性的，但是在很多现实任务中，可能需要非线性映射才能找到合适的降维空间来降维。非线性降维的一种常用方法是基于核技巧对线性降维方法进行核化(kernelized)。核主成分分析(Kernelized PCA，KPCA)是对PCA的一种推广。

假定原始属性空间中的样本点x→i通过将映射到高维特征空间的坐标为x→i,，即x→i,=(x→i)。且假设高维特征空间是n维的，即： x→i,∈Rn。
假定要将高维特征空间中的数据投影到低维空间中，投影矩阵W为n×d维矩阵，根据PCA推导的结果，要求解方程： 

XXTW=λW

其中，X=(x→1,,x→2,,…,x→N,)为n×N维矩阵，于是有： 

∑Ni=1(x→i)(x→i)TW=λW

通常并不清楚的解析表达式，于是引入核函数： 

k(x→i,x→j)=(x→i)T(x→j)

定义核矩阵： 

K=k(x→1,x→1)k(x→1,x→2)…k(x→1,x→N)
k(x→2,x→1)k(x→2,x→2)…k(x→2,x→N)
︙︙︙
k(x→N,x→1)k(x→N,x→2)…k(x→N,x→N)

则有XTX=K。
定义： 

α→i=x→Ti,Wλ

则α→i为1×d维行向量。
定义： A=(α→1,α→2,…,α→N)T为N×d维矩阵，则有： 

W=1λ∑Ni=1x→i,x→Ti,W=∑Ni=1x→i,x→Ti,Wλ=∑Ni=1x→i,α→i=XA

将W=XA代入

XXTW=λW

得到： 

XXTXA=λXA

两边同时左乘以XT，再代入XTX=K，有： 

KKA=λKA

如果要求核矩阵可逆，则上式两边同时左乘以K-1，则有： 

KA=λA

同样该问题也是一个特征值分解问题，取K最大的d个特征值对应的特征向量组成W即可。
对于新样本x→，其投影后第j(j=1,2,…,d)维的坐标为： 

zj=wTj(x→)=∑Ni=1α(j)i(x→i)T(x→)=∑Ni=1α(j)ik(x→i,x→)

其中，α(j)i为行向量α→i的第j个分量。可以看到，为了获取投影后的坐标，KPCA需要对所有样本求和，因此它的计算开销更大。
【例34】对数据实现非线性映射降维(KPCA方法)。

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from scipy.spatial.distance import pdist,squareform

from scipy import exp

from scipy.linalg import eigh

from sklearn.datasets import make_moons

from sklearn.datasets import make_circles

from sklearn.decomposition import PCA

from matplotlib.ticker import FormatStrFormatter

def rbf_kernel_pca(X,gama,n_components):

#1：计算样本对欧几里得距离，并生成核矩阵k(x,y)=exp(-gama *||x-y||^2)，x和y表示#样本，构建一个NXN的核矩阵，矩阵值是样本间的欧氏距离值。计算两两样本间欧几里得距离

sq_dists = pdist (X, 'sqeuclidean') 

##距离平方

mat_sq_dists=squareform(sq_dists) 

##计算对称核矩阵

K=exp(-gama * mat_sq_dists) 

#2:聚集核矩阵K'=K-L*K-K*L + L*K*L，其中L是一个n×n的矩阵(和核矩阵K的维数相#同，所有的值都是1/n。聚集核矩阵的必要性是：样本经过标准化处理后，当在生成协方差矩阵#并以非线性特征的组合替代点积时，所有特征的均值为0；但用低维点积计算时并没有精确计#算新的高维特征空间，也无法确定新特征空间的中心在零点

N=K.shape［0］

one_n = np.ones((N,N))/N #N×N单位矩阵

K=K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

#3：对聚集后的核矩阵求取特征值和特征向量 

eigvals,eigvecs = eigh(K)



#4：选择前K个特征值所对应的特征向量，和PCA不同，KPCA得到的K个特征，不是主成分轴，而#是高维映射到低维后的低维特征数量核化过程是低维映射到高维，PCA是降维，经过核化后的#维度已经不是原来的特征空间。核化是低维映射到高维，但并不是在高维空间计算(非线性特#征组合)而是在低维空间计算(点积)，做到这点关键是核函数，核函数通过两个向量点积来度#量向量间相似度，能在低维空间内近似计算出高维空间的非线性特征空间

X_pc = np.column_stack((eigvecs［:,-i］ for i in range(1,n_components+1))) 

return X_pc 

###分离半月形数据

##生成二维线性不可分数据

X,y=make_moons(n_samples=100,random_state=123)

plt.scatter(X［y==0,0］,X［y==0,1］,color='red',marker='^',alpha=0.5)

plt.scatter(X［y==1,0］,X［y==1,1］,color='blue',marker='o',alpha=0.5)

plt.show()

##PCA降维，映射到主成分，仍不能很好地进行线性分类

sk_pca = PCA(n_components=2)

X_spca=sk_pca.fit_transform(X)

fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))

ax［0］.scatter(X_spca［y==0,0］,X_spca［y==0,1］,color='red',marker='^',alpha=0.5)

ax［0］.scatter(X_spca［y==1,0］,X_spca［y==1,1］,color='blue',marker='o',alpha=0.5)

ax［1］.scatter(X_spca［y==0,0］,np.zeros((50,1))+0.02,color='red',marker='^',alpha=0.5)

ax［1］.scatter(X_spca［y==1,0］,np.zeros((50,1))-0.02,color='blue',marker='^',alpha=0.5)

ax［0］.set_xlabel('PC1')

ax［0］.set_ylabel('PC2')

ax［1］.set_ylim(［-1,1］)

ax［1］.set_yticks(［］)

ax［1］.set_xlabel('PC1')

plt.show()

##利用基于RBF核的KPCA来实现线性可分

X_kpca=rbf_kernel_pca(X, gama=15, n_components=2)

fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))

ax［0］.scatter(X_kpca［y==0,0］,X_kpca［y==0,1］,color='red',marker='^',alpha=0.5)

ax［0］.scatter(X_kpca［y==1,0］,X_kpca［y==1,1］,color='blue',marker='o',alpha=0.5)

ax［1］.scatter(X_kpca［y==0,0］,np.zeros((50,1))+0.02,color='red',marker='^',alpha=0.5)

ax［1］.scatter(X_kpca［y==1,0］,np.zeros((50,1))-0.02,color='blue',marker='^',alpha=0.5)

ax［0］.set_xlabel('PC1')

ax［0］.set_ylabel('PC2')

ax［1］.set_ylim(［-1,1］)

ax［1］.set_yticks(［］)

ax［1］.set_xlabel('PC1')

ax［0］.xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))

ax［1］.xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))

plt.show()



###分离同心圆

##生成同心圆数据

X,y=make_circles(n_samples=1000,random_state=123,noise=0.1,factor=0.2)

plt.scatter(X［y==0,0］,X［y==0,1］,color='red',marker='^',alpha=0.5)

plt.scatter(X［y==1,0］,X［y==1,1］,color='blue',marker='o',alpha=0.5)

plt.show()

##标准PCA映射

sk_pca = PCA(n_components=2)

X_spca=sk_pca.fit_transform(X)

fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))

ax［0］.scatter(X_spca［y==0,0］,X_spca［y==0,1］,color='red',marker='^',alpha=0.5)

ax［0］.scatter(X_spca［y==1,0］,X_spca［y==1,1］,color='blue',marker='o',alpha=0.5)

ax［1］.scatter(X_spca［y==0,0］,np.zeros((500,1))+0.02,color='red',marker='^',alpha=0.5)

ax［1］.scatter(X_spca［y==1,0］,np.zeros((500,1))-0.02,color='blue',marker='^',alpha=0.5)

ax［0］.set_xlabel('PC1')

ax［0］.set_ylabel('PC2')

ax［1］.set_ylim(［-1,1］)

ax［1］.set_yticks(［］)

ax［1］.set_xlabel('PC1')

plt.show()

##RBF-KPCA映射

X_kpca=rbf_kernel_pca(X, gama=15, n_components=2)

fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(7,3))

ax［0］.scatter(X_kpca［y==0,0］,X_kpca［y==0,1］,color='red',marker='^',alpha=0.5)

ax［0］.scatter(X_kpca［y==1,0］,X_kpca［y==1,1］,color='blue',marker='o',alpha=0.5)

ax［1］.scatter(X_kpca［y==0,0］,np.zeros((500,1))+0.02,color='red',marker='^',alpha=0.5)

ax［1］.scatter(X_kpca［y==1,0］,np.zeros((500,1))-0.02,color='blue',marker='^',alpha=0.5)

ax［0］.set_xlabel('PC1')

ax［0］.set_ylabel('PC2')

ax［1］.set_ylim(［-1,1］)

ax［1］.set_yticks(［］)

ax［1］.set_xlabel('PC1')

ax［0］.xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))

ax［1］.xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))

plt.show()

运行程序，效果如图34所示。


图34KPCA降维效果

3.5流形学习降维
流形学习(Manifold Learning)是一类借鉴了拓扑流形概念的降维方法，被认为属于非线性降维的一个分支。流形学习假设所处理的数据点分布在嵌入于欧氏空间的一个潜在的流形体上，或者说这些数据点可以构成这样一个潜在的流形体。流形学习的方法有很多，并具有一些共同的特征： 首先需要构造流形上样本点的局部邻域结构，然后用这些局部邻域结构将样本点全局地映射到一个降维空间。这些方法之间的不同之处主要在于构造的局部邻域结构不同，以及利用这些局部邻域结构来构造全局的低维嵌入方法的不同。
3.6多维缩放降维
3.6.1原理

多维缩放(Multiple Dimensional Scaling，MDS)要求原始空间中样本之间的距离在低维空间中得到保持。
假设N个样本在原始空间中的距离矩阵为D=(di，j)N×N： 

D=d1,1d1,2…d1,N
d2,1d2,2…d2,N
︙︙︙
dN,1dN,2…dN,N

其中，di，j=‖x→i-x→j‖为样本x→i到样本x→j的距离。
假设原始样本是在n维空间，我们的目标是在n,n′<n维空间里获取样本，欧氏距离保持不变。
假设样本集在原空间的表示X=(x→1,x→2,…,x→N)为n×N维矩阵，样本集在降维后空间的坐标Z=(z→1,z→2,…,z→N)为n′×N维矩阵。所求的正是Z矩阵。
令B=ZTZ为N×N维矩阵，即

B=b1,1b1,2…b1,N
b2,1b2,2…b2,N
︙︙︙
bN,1bN,2…bN,N

其中，bi,j=z→i·z→j为降维后样本的内积。
根据降维前后样本的欧氏距离保持不变，有： 

d2i，j=‖z→i-z→j‖2=‖z→i‖2+‖z→j‖2-2z→Tiz→j=bi,i+bj,j-2bi,j

假设降维后的样本集Z被中心化，即∑Ni=1z→i=0→，则矩阵B的每行之和均为零，每列之和均为零，即： 

∑Ni=1bi,j=0,j=1,2,…,N

∑Nj=1bi,j=0,i=1,2,…,N

于是有： 

∑Ni=1d2i，j=∑Ni=1bi,i+Nbj,j=tr(B)+Nbj,j

∑Nj=1d2i，j=∑Ni=1bj,j+Nbi,i=tr(B)+Nbi,i

∑Ni=1∑Ni=1d2i，j=∑Ni=1tr(B)+Nbi,i=2Ntr(B)

其中，tr(B)表示矩阵B的迹。
令

d2i,·=1N∑Ni=1d2ij=tr(B)N+bi,i

d2j,·=1N∑Nj=1d2ij=tr(B)N+bj,j

d2·,·=1N2∑Ni=1∑Nj=1d2ij=2tr(B)N

代入d2i，j=bi,i+bj,j-2bi,j，有

bi,j=bi,i+bj,j-d2i，j2=d2i,·+d2j,·-d2·,·-d2i，j2

上面右式根据di，j给出了bi,j，因此可以根据原始空间中的距离矩阵D求出在降维后空间的内积矩阵B。现在的问题是： 已知内积矩阵B=ZTZ，如何求得矩阵Z。
对矩阵B做特征值分解，设B=VΛVT，其中Λ=diag(λ1,λ2,…,λN)为特征值构成的对角矩阵，其中λ1≥λ2≥…≥λN，V为特征向量矩阵。
假定特征中有n*个非零特征值，它们构成对角矩阵Λ*=diag(λ1,λ2,…,λn*)。令V*为对应的特征向量矩阵，则

Z=Λ1/2*VT*

其中，Z为n*×N阶矩阵，此时有n′=n*。
在现实应用中，为了有效降维，往往仅需要降维后的距离与原始空间中的距离尽可能相等，而不必严格相等。此时可以取n′n*≤n个最大特征值构成的对角矩阵为： 

Λ～=diag(λ1,λ2,…,λn′)

令V～表示对应的特征向量矩阵，则

Z=Λ～1/2VT∈Rn′×N

3.5.2MDS算法
多维缩放(MDS)算法如下。
(1) 输入： 距离矩阵D∈RN×N； 低维空间维数n′。
(2) 输出： 样本集在低维空间中的矩阵Z。
(3) 算法步骤为： 
 根据下列式子计算d2i,·、d2j,·、d2·,·。

d2i,·=1N∑Ni=1d2i，j=tr(B)N+bi,i

d2j,·=1N∑Nj=1d2i，j=tr(B)N+bj,j

d2·,·=1N2∑Ni=1∑Nj=1d2i，j=2tr(B)N

 根据正式计算矩阵B： 

bi,j=bi,i+bj,j-d2i，j2=d2i,·+d2j,·-d2·,·-d2i，j2

 对矩阵B进行特征值分解。
 Λ～为n′个最大特征值所构成的对角矩阵，V～表示对应的特征向量矩阵，则： 

Z=Λ～1/2VT∈Rn′×N

【例35】利用MDS算法实现数据的降维。

"""

MDS降维

"""

import numpy as np

import matplotlib.pyplot as plt

from sklearn import datasets,manifold

def load_data():

'''

加载用于降维的数据

:return: 一个元组，依次为训练样本集和样本集的标记

'''

iris=datasets.load_iris()#使用 scikit-learn 自带的 iris 数据集

returniris.data,iris.target

def test_MDS(*data):

'''

测试 MDS 的用法

:param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、训练样本的标记

'''

X,y=data

for n in ［4,3,2,1］: #依次考查降维目标为四维、三维、二维、一维

mds=manifold.MDS(n_components=n)

mds.fit(X)

print('stress(n_components=%d) : %s'% (n, str(mds.stress_)))

def plot_MDS(*data):

'''

绘制经过 使用 MDS 降维到二维之后的样本点

:param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、训练样本的标记 

'''

X,y=data

mds=manifold.MDS(n_components=2)

X_r=mds.fit_transform(X) #原始数据集转换到二维

###绘制二维图形

fig=plt.figure()

ax=fig.add_subplot(1,1,1)

colors=((1,0,0),(0,1,0),(0,0,1),(0.5,0.5,0),(0,0.5,0.5),(0.5,0,0.5),

(0.4,0.6,0),(0.6,0.4,0),(0,0.6,0.4),(0.5,0.3,0.2),)#颜色集合，不同标记的样本
#染不同的颜色

for label ,color in zip( np.unique(y),colors):

position=y==label

ax.scatter(X_r［position,0］,X_r［position,1］,label="target= %d"%label,color=color)

ax.set_xlabel("X［0］")

ax.set_ylabel("X［1］")

ax.legend(loc="best")

ax.set_title("MDS")

plt.show()

if __name__=='__main__':

X,y=load_data() #产生用于降维的数据集

test_MDS(X,y) #调用 test_MDS

plot_MDS(X,y) #调用 plot_MDS

运行程序，输出如下，效果如图35所示。

stress(n_components=4) : 11.887221490372065

stress(n_components=3) : 27.590871404910008

stress(n_components=2) : 113.30127128702554

stress(n_components=1) : 28321.42085807291



图35MDS降维效果

3.7等度量映射降维
等度量映射(Isometric Mapping，Isomap)原理如下： 
(1) 首先建立近邻连接图： 利用流形在局部上与欧氏空间同胚这个性质，基于欧氏距离，对每个点找出它在低维流形上的近邻点，建立近邻连接图。
(2) 计算任意两点之间的距离： 计算近邻连接图上任意两点之间的最短路径问题，作为两点之间的距离。
(3) 在得到任意两点的距离之后，就可以通过MDS算法来获得样本点在低维空间中的坐标。
Isomap算法如下。
(1) 输入： 样本集D={x→1,x→2,…,x→N}； 近邻参数k； 低维空间维数n′。
(2) 输出： 样本集在低维空间中的矩阵Z。
(3) 算法步骤为： 
 对每个样本点x→i，计算它的k近邻； 同时将x→i与它的k近邻的距离设置为欧氏距离，与其他点的距离设置为无穷大。
 调用最短路径算法计算任意两个样本点之间的距离，获得距离矩阵D∈RN×N。
 调用多维缩放MDS算法，获得样本集在低维空间中的矩阵Z。
Isomap算法有个很大的问题： 对于新样本，难以将其映射到低维空间。理论上可以将新样本添加到样本集中，重新调用Isomap算法，这种方案计算量太大。一般的解决方法是： 训练一个回归学习器来对新样本的低维空间进行预测。
近邻图有如下两种类型。
(1) k近邻图： 指定近邻点个数，如指定距离最近的k个点为近邻点。
(2) ε近邻图： 指定距离阈值ε，距离小于ε的点被认为是近邻点。
在建立近邻图的时候要注意控制近邻图的范围，否则容易出现“短路”或者“断路”问题。
(1) “短路”问题： 近邻范围指定过大，距离很远的点也被误认为是近邻。
(2) “断路”问题： 近邻范围指定过小，本应该相连的区域被认为是断开的。
【例36】利用等度量映射(Isomap)算法对数据进行降维。
(1) 导入必要的编程库。

import numpy as np

import matplotlib.pyplot as plt

from sklearn import datasets,decomposition,manifold

(2) 加载数据。

def load_data():

iris=datasets.load_iris()

return iris.data,iris.target

(3) 使用somap od。

def test_Isomap(*data):

X,y=data

for n in ［4,3,2,1］:

isomap=manifold.Isomap(n_components=n)

isomap.fit(X)

print('reconstruction_error(n_components=%d):%s'%(n,

isomap.reconstruction_error()))

X,y=load_data()

test_Isomap(X,y)

(4) 降维后的样本分布图。

def plot_Isomap(*data):

X,y=data

Ks=［1,5,25,y.size-1］

fig=plt.figure()

for i,k in enumerate(Ks):

isomap=manifold.Isomap(n_components=2,n_neighbors=k)

X_r=isomap.fit_transform(X)

ax=fig.add_subplot(2,2,i+1)

colors=((1,0,0),(0,1,0),(0,0,1),(0.5,0.5,0),(0,0.5,0.5),(0.5,0,0.5),

 (0.4,0.6,0),(0.6,0.4,0),(0,0.6,0.4),(0.5,0.3,0.2),)

for label,color in zip(np.unique(y),colors):

position=y==label



ax.scatter(X_r［position,0］,X_r［position,1］,label='target=%d'%label,color=color)

ax.set_xlabel('X［0］')

ax.set_ylabel('X［1］')

ax.legend(loc='best')

ax.set_title("k=%d"%k)

plt.suptitle('Isomap')

plt.show()

plot_Isomap(X,y)

(5) 将原始数据的特征直接压缩到一维。

def plot_Isomap_k_d1(*data):

X,y=data

Ks=［1,5,25,y.size-1］

fig=plt.figure()

for i,k in enumerate(Ks):

isomap=manifold.Isomap(n_components=2,n_neighbors=k)

X_r=isomap.fit_transform(X)

ax=fig.add_subplot(2,2,i+1)

colors=((1,0,0),(0,1,0),(0,0,1),(0.5,0.5,0),(0,0.5,0.5),(0.5,0,0.5),

 (0.4,0.6,0),(0.6,0.4,0),(0,0.6,0.4),(0.5,0.3,0.2),)

for label,color in zip(np.unique(y),colors):

position=y==label

ax.scatter(X_r［position］,np.zeros_like(X_r［position］),label='target=%d'%label,color=color)

ax.set_xlabel('X［0］')

ax.set_ylabel('Y')

ax.legend(loc='best')

ax.set_title("k=%d"%k)

plt.suptitle('Isomap')

plt.show()

plot_Isomap_k_d1(X,y)

运行程序，输出如下，效果如图36和图37所示。

reconstruction_error(n_components=4):1.0097180068081741

reconstruction_error(n_components=3):1.0182845146289834

reconstruction_error(n_components=2):1.0276983764330463

reconstruction_error(n_components=1):1.0716642763207656



图36降维后的样本分布图




图37将原始数据的特征直接压缩到一维效果


3.8局部线性嵌入
3.8.1原理

局部线性嵌入(Locally Linear Embedding，LLE)的目标是： 保持邻域内样本之间的线性关系。
对每个样本x→i，首先寻找其近邻点，假设这些近邻点的下标集合为Qi。然后需要计算基于x→i的近邻点对x→j进行线性重构的系数w→i。定义样本集重构误差为： 

err=∑Ni=1x→i-∑j∈Qiwi,jx→j22

其中，wi,j为w→i的分量。我们的目标是使样本集重构误差最小，即： 

minw→1,w→2,…,w→N∑Ni=1x→i-∑j∈Qiwi,jx→j22

这样的解有无数个，对权重进行归一化处理，即： 

∑j∈Qiwi,j=1,i=1,2,…,N

这样一来，就是求解最优化问题： 

minw→1,w→2,…,w→N ∑Ni=1x→i-∑j∈Qiwi,jx→j22
s.t.∑j∈Qiwi,j=1,i=1,2,…,N

该最优化问题有解析解。令Cj,k=(x→i-x→j)T(x→i-x→k)，则可以解出： 

wi,j=∑k∈QiC-1j,k∑l,s∈QiC-1l,s,j∈Qi

求出了线性重构的系数w→i之后，LLE在低维空间中保持w→i不变。设z→i对应的低维坐标为z→j，已知线性重构的系数w→i，定义样本集在低维空间中重构误差为： 

err′=∑Ni=1z→i-∑j∈Qiwi,jz→j22

现在的问题是要求出z→i，从而使上式最小。即求解： 

minz→1,z→2,…,z→N∑Ni=1z→i-∑j∈Qiwi,jz→j22

令Z=(z→1,z→2,…,z→N)∈Rn′×N，其中n′为低维空间的维数(n为原始样本所在的高维空间的维数)。令： 

W=w1,1w1,2…w1,N
w2,1w2,2…w2,N
︙︙︙
wN,1wN,2…wN,N

定义M=(I-W)T(I-W)，于是最优化问题可重写为： 

minZtr(ZMZT)

该最优化问题有无数个解。添加约束ZZT=I，于是最优化问题为： 

minZtr(ZMZT)

s.t.ZZT=I

该最优化问题可以通过特征值分解求解： 选取M最小的n′个特征值对应的特征向量组成的矩阵即为ZT。
3.8.2LLE算法
LLE算法如下。
(1) 输入： 样本集D={x→1,x→2,…,x→N}： 近邻参数k； 低维空间维数n′。
(2) 输出： 样本集在低维空间中的矩阵Z。
(3) 算法步骤： 
 对于样本集中的每个点x→i,i=1,2,…,N，执行下列操作。

◎ 确定x→i的k近邻，获得其近邻下标集合Qi。

◎ 对于j∈Qi，根据下式计算wi,j。

wi,j=∑k∈QiC-1j,k∑l,s∈QiC-1l,s

Cj,k=(x→i-x→j)T(x→i-x→k)

◎ 对于jQi，wi,j=0。
 根据wi,j构建矩阵W。
 计算M=(I-W)T(I-W)。
 对M进行特征值分解，取其最小的n′个特征值对应的特征向量，即得到样本集在低维空间中的矩阵Z。
【例37】利用局部线性嵌入对数据进行降维。

import numpy as np

import matplotlib.pyplot as plt

from sklearn import datasets,manifold



def load_data():

'''

加载用于降维的数据

:return: 一个元组，依次为训练样本集和样本集的标记

'''

iris=datasets.load_iris()#使用 scikit-learn 自带的 iris 数据集

returniris.data,iris.target

def test_LocallyLinearEmbedding(*data):

'''

测试 LocallyLinearEmbedding 的用法

:param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、训练样本的标记

'''

X,y=data

for n in ［4,3,2,1］:#依次考查降维目标为四维、三维、二维、一维

lle=manifold.LocallyLinearEmbedding(n_components=n)

lle.fit(X)

print('reconstruction_error(n_components=%d) : %s'%

(n, lle.reconstruction_error_))

def plot_LocallyLinearEmbedding_k(*data):

'''

测试 LocallyLinearEmbedding 中 n_neighbors 参数的影响，其中降维至二维

:param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、训练样本的标记

:return: None

'''

X,y=data

Ks=［1,5,25,y.size-1］#n_neighbors参数的候选值的集合

fig=plt.figure()

for i, k in enumerate(Ks):

lle=manifold.LocallyLinearEmbedding(n_components=2,n_neighbors=k)

X_r=lle.fit_transform(X)#原始数据集转换到二维

ax=fig.add_subplot(2,2,i+1)##两行两列，每个单元显示不同 n_neighbors 参数的 LocallyLinearEmbedding 的效果图

colors=((1,0,0),(0,1,0),(0,0,1),(0.5,0.5,0),(0,0.5,0.5),(0.5,0,0.5),

(0.4,0.6,0),(0.6,0.4,0),(0,0.6,0.4),(0.5,0.3,0.2),)#颜色集合，不同标记的样
#本染不同的颜色

for label ,color in zip( np.unique(y),colors):

position=y==label

ax.scatter(X_r［position,0］,X_r［position,1］,label="target= %d"

%label,color=color)

ax.set_xlabel("X［0］")

ax.set_ylabel("X［1］")

ax.legend(loc="best")

ax.set_title("k=%d"%k)

plt.suptitle("LocallyLinearEmbedding")

plt.show()

def plot_LocallyLinearEmbedding_k_d1(*data):

'''

测试 LocallyLinearEmbedding中n_neighbors 参数的影响，其中降维至一维

:param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、训练样本的标记 

'''

X,y=data

Ks=［1,5,25,y.size-1］#n_neighbors参数的候选值的集合

fig=plt.figure()

for i, k in enumerate(Ks):

lle=manifold.LocallyLinearEmbedding(n_components=1,n_neighbors=k)

X_r=lle.fit_transform(X)#原始数据集转换到一维

ax=fig.add_subplot(2,2,i+1)#两行两列，每个单元显示不同n_neighbors参数

#的 LocallyLinearEmbedding 的效果图

#颜色集合，不同标记的样本染不同的颜色

colors=((1,0,0),(0,1,0),(0,0,1),(0.5,0.5,0),(0,0.5,0.5),(0.5,0,0.5),

(0.4,0.6,0),(0.6,0.4,0),(0,0.6,0.4),(0.5,0.3,0.2),)

for label ,color in zip( np.unique(y),colors):

position=y==label

ax.scatter(X_r［position］,np.zeros_like(X_r［position］),

label="target= %d"%label,color=color)#FFFFFF

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.legend(loc="best")

ax.set_title("k=%d"%k)

plt.suptitle("LocallyLinearEmbedding")

plt.show()

if __name__=='__main__':

X,y=load_data() #产生用于降维的数据集

test_LocallyLinearEmbedding(X,y) #调用 test_LocallyLinearEmbedding

plot_LocallyLinearEmbedding_k(X,y) #调用 plot_LocallyLinearEmbedding_k

plot_LocallyLinearEmbedding_k_d1(X,y)#调用 plot_LocallyLinearEmbedding_k_d1

运行程序，输出如下，效果如图38和图39所示。

reconstruction_error(n_components=4) : 7.199368860901911e-07

reconstruction_error(n_components=3) : 3.870605052928055e-07

reconstruction_error(n_components=2) : 6.641420116916785e-08

reconstruction_error(n_components=1) : 1.6515558016659176e-16



图38数据LLE降维效果




图39将数据降到一维的效果


3.9非负矩阵分解
非负矩阵分解(Nonnegative Matrix Factorization，NMF)是在矩阵中所有元素均为非负数约束条件之下的矩阵分解方法。其基本思想： 给定一个非负矩阵V，NMF能够找到一个非负矩阵W和一个非负矩阵H，使得矩阵W和H的乘积近似等于矩阵V中的值。

Vn×m=Wn×k×Hk×m

其中，W为基础图像矩阵，相当于从原始矩阵V中抽取出来的特征； H矩阵为系数矩阵。
NMF广泛应用于图像分析、文本挖掘和语音处理等领域。
最小化W矩阵H矩阵的乘积和原始矩阵之间的差别，其目标函数为： 

argmin12‖H-WH‖2=12∑i,j(xij-WHij)2

基于KL散度的优化目标，其损失函数为： 

argminJ(W,H)=∑i,jxijlnxijWHij-xij+WHij

【例38】在sklearn封装了NMF的实现，可以非常方便地使用，其实现基本和理论部分的实现是一致的，但应注意sklearn中输入数据的格式是(samples, features)。

#导入必要的编程库

from sklearn.decomposition import NMF

from sklearn.datasets import load_iris

#载入数据

X, _ = load_iris(True)

#最重要的参数是n_components，alpha，l1_ratio，solver

nmf = NMF(n_components=2,#k value,默认会保留全部特征

init=None,#W H 的初始化方法，包括'random' | 'nndsvd'(默认) |'nndsvda' | 'nndsvdar' | 'custom'.

solver='cd',#'cd' | 'mu'

#{'frobenius', 'kullback-leibler', 'itakura-saito'}，一般默认就好

beta_loss='frobenius', 

tol=1e-4, #停止迭代的极限条件

max_iter=200, #最大迭代次数

random_state=None,

alpha=0., #正则化参数

l1_ratio=0., #正则化参数

verbose=0, #冗长模式

shuffle=False #针对"cd solver"

)

#-----------------函数------------------------

print('params:', nmf.get_params()) #获取构造函数参数的值，也可以通过

#nmf.attr得到，所以下面会省略这些属性

#下面的4个函数很简单，也最核心

nmf.fit(X)

W = nmf.fit_transform(X)

W = nmf.transform(X)

nmf.inverse_transform(W)

#-----------------属性------------------------

H = nmf.components_#H矩阵

print('reconstruction_err_', nmf.reconstruction_err_)#损失函数值

print('n_iter_', nmf.n_iter_)#实际迭代次数

运行程序，输出如下：

params: {'alpha': 0.0, 'beta_loss': 'frobenius', 'init': None, 'l1_ratio': 0.0, 'max_iter': 200, 'n_components': 2, 'random_state': None, 'shuffle': False, 'solver': 'cd', 'tol': 0.0001, 'verbose': 0}

reconstruction_err_ 3.9480195652425465

n_iter_ 199

在以上代码中，各参数的含义为： 
 init参数中，nndsvd(默认)更适用于sparse factorization，其变体则适用于dense factorization。
 solver参数中，如果初始化中产生很多零值，Multiplicative Update (mu) 不能很好地更新。所以mu一般不和nndsvd一起使用，而和其变体nndsvda、nndsvdar一起使用。
 solver参数中，cd只能优化Frobenius norm函数，而mu可以更新所有损失函数。
【例39】一个NMF在图像特征提取的应用例子。

#导入必要的编程库

from time import time

from numpy.random import RandomState

import matplotlib.pyplot as plt

from sklearn.datasets import fetch_olivetti_faces

from sklearn import decomposition

#设置参数

n_row, n_col = 2, 3

n_components = n_row * n_col

image_shape = (64, 64)

rng = RandomState(0)

#载入face数据

dataset = fetch_olivetti_faces('./', True, random_state=rng)

faces = dataset.data

n_samples, n_features = faces.shape

print("Dataset consists of %d faces, features is %s" % (n_samples, n_features))

#显示原始图像

def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):

plt.figure(figsize=(2. * n_col, 2.26 * n_row))

plt.suptitle(title, size=16)

for i, comp in enumerate(images):

plt.subplot(n_row, n_col, i + 1)

vmax = max(comp.max(), -comp.min())

#显示压缩后的图像

plt.imshow(comp.reshape(image_shape), cmap=cmap,

interpolation='nearest',

vmin=-vmax, vmax=vmax)

plt.xticks(())

plt.yticks(())

plt.subplots_adjust(0.01, 0.05, 0.99, 0.93, 0.04, 0.)



estimators = ［

('Non-negative components - NMF',

decomposition.NMF(n_components=n_components, init='nndsvda', tol=5e-3))

］

#绘制输入数据的示例

plot_gallery("First centered Olivetti faces", faces［:n_components］)



#估算并绘制它

for name, estimator in estimators:

print("Extracting the top %d %s..." % (n_components, name))

t0 = time()

data = faces

estimator.fit(data)

train_time = (time() - t0)

print("done in %0.3fs" % train_time)

components_ = estimator.components_

print('components_:', components_.shape, '＼n**＼n', components_)

plot_gallery('%s - Train time %.1fs' % (name, train_time),

components_)

plt.show()

运行程序，输出如下，效果如图310和图311所示。

downloading Olivetti faces from https://ndownloader.figshare.com/files/5976027 to ./

Dataset consists of 400 faces, features is 4096

Extracting the top 6 Non-negative components - NMF...

done in 0.500s

components_: (6, 4096)

**

［［0.  0.   0. ... 1.89640523 1.78331733 1.68142998］

［0.95390665 1.02885565 1.09771352 ... 0.25895402 0.31447183 0.32589285］

［0.02238854 0.01136337 0.06086569 ... 0.13581616 0.14109341 0.15405469］

［0.18907826 0.26982826 0.38558988 ... 0. 0. 0.］

［0. 0. 0. ... 0. 0. 0.］

［0.51045833 0.52450599 0.50316775 ... 0.03719877 0.02762018 0.04051619］］




图310降维前的图像




图311降维后的图像



3.10小结
数据降维基本原理是将样本点从输入空间通过线性或非线性变换映射到一个低维空间，从而获得一个关于原数据集紧致的低维表示。本章从维度灾难与降维、主成分分析、SVD降维、核主成分分析(KPCA)降维、多维缩放(MDS)降维、局部线性嵌入(LLE)、非负矩阵分解等多个方面介绍了数据降维相关内容，每个小节都是通过理论、图文、实例相结合进行数据降维介绍，让读者快速上手利用Python解决实际降维问题。
3.11习题
1. 数据降维，一方面可以解决，缓解、现状，降低复杂度； 另一方面可以更好地和数据。
2. 根据是否考虑和利用数据的监督信息可以划分为、和。
3. 缓解维度灾难的一个重要途径是什么？
4. 什么是PCA？
5.  PCA降维的准则有几个？分别是什么？