机器学习——数据降维

数据降维

数据降维又称维数约简，就是降低数据的维度。其方法有很多种，从不同角度入手可以有不同的分类，主要分类方法有：根据数据的特性可以划分为线性降维和非线性降维，根据是否考虑和利用数据的监督信息可以划分为无监督降维、有监督降维和半监督降维，根据保持数据的结构可以分为全局保持降维、局部保持降维和全局与局部保持一致降维等。需要根据特定的问题选择合适的数据降维方法。

数据降维一方面可以解决“维数灾难”，缓解信息丰富、知识贫乏的现状，降低复杂度，另一方面可以更好地认识和理解数据。

1 维度灾难与降维

维度灾难用来描述当空间维度增加时，分析和组织高维空间，因体积指数增加而遇到各种问题场景。例如在高维情况下，数据样本稀疏的问题。例如k近邻算法中，任意样本附近任意小的距离内总能找到一个训练样本，即训练样本的采样密度足够大才能保证分类性能，当特征维度很大时，满足密采样的样本数量会呈现指数级增大，大到几乎无法达到。所以在高维度情况下，涉及距离、内积的计算变得困难。

缓解维度灾难的一个重要途径就是降维。一般来说，想获得低维子空间，最简单的方法就是对原始高维度空间进行线性变换：给定d维空间中的样本X=(x1,x2,···,xm)变换之后得到的d’<=d维空间中的样本Z=W^TX。

变换矩阵W可视为d’个d维基向量，新空间中的属性是原空间属性的线性组合，基于线性变换来进行降维的方法都成线性维方法，主要区别于对低维子空间的性质有所不同，相当于对W施加了不同的约束。

2 主成分分析

主成分分析是一种最常用的无监督降维方法，通过降维技术吧多个变量化为少数几个主成分的统计分析方法。这些主成分能够反映原始变量的绝大部分信息，它们通常表示为原始变量的某种线性组合。

2.1 PCA原理

为了便于维度变换有如下假设：

• 假设样本数据是n维的。
• 假设原始坐标系为由标准正交基向量{i1,i2,···,in}张成的空间。
• 假设经过线性变换后的新坐标系为由标准正交基向量{j1,j2,···,jn}张成的空间。

根据定义有：
记ws=(js·i1,js·i2,···,js·in)，其各个分量就是基向量js在原始坐标系中的投影。

根据定义，有：
令坐标变换矩阵W为：W=(w1,w2,···,wn)，则有：

这样W的第s列就是js在原始坐标系中的投影。

假设样本点xi在原始坐标系中的表示为：
令其中：

假设样本点xi在新坐标系中的表示为：
令其中：
于是可以推导出：

丢弃其中的部分坐标，将维度降低到d<n，则样本点xi在低维坐标系中的坐标：

那么，现在出现一个问题，即丢弃哪些坐标？思想是：基于降低之后的坐标重构样本时，尽量要与原始样本接近。

2.2 PCA算法

算法的输入为样本集D与低维空间维数d，输出为所求的投影矩阵W。

算法的步骤表现在：

• 对所有样本进行中心化操作：
  
• 计算样本的协方差矩阵XX^T。
• 对协方差矩阵作特征值分解。
• 取最大的d个特征值对应的特征向量构造投影矩阵W。

通常低维空间维数d的选取有两种方法：

• 通过交叉验证法选取较好的d。
• 从算法原理角度设置一个阈值，例如t=95%，然后选取使得下式成立的最小的d的值：
  

PCA降维有两个准则：

• 最近重构性：样本集中所有点，重构后的点距离原来点的误差之和最小。
• 最大可分性：样本点在低维空间的投影尽可能分开。

2.3 PCA算法实验

加载PCA算法包，并加载PCA算法设置降维后的主成分数目为2。

from sklearn.decomposition import PCA

pca=PCA(n_components=2)     
reduced_x=pca.fit_transform(x)
1
2
3
4

样本的数据结构如下图所示：

数据集包含150个数据集，分为三类，每类50个数据，每个数据包含4个特征。可以通过花萼长度，花萼宽度，花瓣长度，花瓣宽度4个特征预测鸢尾花卉属于（Setosa，Versicolour，Virginica）三个种类中的哪一类。

数据的结构化表示为：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets
from sklearn.decomposition import PCA
import pandas as pd
import seaborn
 
data = datasets.load_iris()
X =data['data']
y =data['target']
a = pd.DataFrame(X, columns=data.feature_names)

#seaborn.boxplot(data= a)
plt.plot(a)
plt.legend(data.feature_names)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

首先选取三个特征查看数据的分布情况：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets
  
data = datasets.load_iris()
X =data['data']
y =data['target']
ax = Axes3D(plt.figure())
for c, i, target_name in zip('rgb', [0, 1, 2], data.target_names):
    ax.scatter(X[y==i, 0], X[y==i, 2], c=c, label=target_name)
 
ax.set_xlabel(data.feature_names[0])
ax.set_xlabel(data.feature_names[1])
ax.set_xlabel(data.feature_names[2])
ax.set_title('Iris')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

可以发现红色数据点，也就是Setosa距离其他数据较远。

然后选择两个特征子集查看数据分布情况。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets
 
data = datasets.load_iris()
X =data['data']
y =data['target']
ax = Axes3D(plt.figure())
for c, i, target_name in zip('rgb', [0, 1, 2], data.target_names):
    ax.scatter(X[y==i, 0], X[y==i, 1], c=c, label=target_name)
 
ax.set_xlabel(data.feature_names[0])
ax.set_xlabel(data.feature_names[1])
ax.set_title('Iris')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

可以发现红色的Setosa仍然是线性可分的。

利用PCA将数据降到二维：

import matplotlib.pyplot as plt                
from sklearn.decomposition import PCA          
from sklearn.datasets import load_iris

data=load_iris()
y=data.target
x=data.data
pca=PCA(n_components=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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

最终得到的主成分降维效果如下图所示：

可以看到用sklearn库实现降维后，将原本的四维数据变成二维，以实现在平面中画出样本点的分布。

3 SVD降维

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

其中，
根据正交矩阵的性质，有：

若令M为：

则可以推导出：

那么，λi(i=1,2,···,r)就是XX^ T的特征值，其对应的特征向量组成正交矩阵V。因此SVD奇异值分解等价于PCA主成分分析，核心都是求解XX^ T的特征值以及对应的特征向量。

下面的实验利用SVD对给定的数据进行降维处理。

用户数据的SVD奇异值分解：

def _svd(self):
        self.U, self.S, self.VT = np.linalg.svd(self.data)
        return self.U, self.S, self.VT 
1
2
3

前k个奇异值的平方和占比要大于等于percentage，求出满足此条件的最小的k值。

def _calc_k(self, percentge):
        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
1
2
3
4
5
6
7
8
9
10
11

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

def _buildSD(self, k):
        self.SD = np.eye(self.k) * self.S[:self.k] 

        e = np.eye(self.k)
        for i in range(self.k):
            e[i,i] = self.S[i] 
        return self.SD
1
2
3
4
5
6
7

SVD降维。

def DimReduce(self, percentage):
        self._svd()

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

        self._buildSD(self.k)
        k,U,SD,VT = self.k,self.U, self.SD, self.VT
        
        a = U[:len(U), :k]
        b = np.dot(SD, VT[:k, :len(VT)])
        newData = np.dot(a,b)
        return newData
1
2
3
4
5
6
7
8
9
10
11
12
13

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

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('====================================================') 
1
2
3
4
5
6
7
8
9
10
11
12
13
14

最后运行得到的结果如下图所示：

4 核主成分分析降维

PCA方法假设从高维空间到低维空间的函数映射是线性的，但是在很多现实任务中，可能需要非线性的映射才能找到合适的降维空间来降维。

非线性降维的一种常用方法是基于核技巧对线性降维方法进行核化。

考虑一个线性不可分的数据集，如果我们用一个映射 ，将其映射到三维空间，比如可取：

可以看到，映射到高维空间后数据变得线性可分。这意味着，非线性映射使得我们能够处理非线性的数据。

假设要将高维特征空间中的数据投影到低维空间中，投影矩阵W根据PCA求解的结果需要求解方程：

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

定义核矩阵：

定义：


可以推导出：

最终通过矩阵的化简可以得到KA=λA。同样该问题也是一个特征值分解问题，取K的最大的d个特征对应的特征向量组成W即可。

对于新样本x，其投影后的坐标为：

可以看到，为了获取投影后的坐标，KPCA需要对所有样本求和，因此它的计算开销更大。

下面实现对非线性映射的降维。

def rbf_kernel_pca(X,gama,n_components):
    sq_dists = pdist (X, 'sqeuclidean')     

    mat_sq_dists=squareform(sq_dists)     

    K=exp(-gama * mat_sq_dists)     
     
    N=K.shape[0]
    one_n = np.ones((N,N))/N 
    K=K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)    

    eigvals,eigvecs = eigh(K)

    X_pc = np.column_stack((eigvecs[:,-i] for i in range(1,n_components+1))) 
    return X_pc 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

首先计算样本对欧几里得距离，并生成核矩阵k(x,y)=exp(-gama * ||x-y||^ 2)，x和y表示样本，构建一个NXN的核矩阵，矩阵值是样本间的欧氏距离值。然后，聚集核矩阵K’=K-LK-KL + LKL，其中L是一个nXn的矩阵(和核矩阵K的维数相同，所有的值都是1/n。

聚集核矩阵的必要性是：样本经过标准化处理后，当在生成协方差矩阵并以非线性特征的组合替代点积时，所有特征的均值为0；但用低维点积计算时并没有精确计算新的高维特征空间，也无法确定新特征空间的中心在零点。

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

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()
1
2
3
4

分离半月形数据，生成二维线性不可分数据。

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()
1
2
3
4
5
6
7
8
9
10
11
12
13

显然，这两个半月形不是线性可分的，而我们的目标是通过核 PCA 将这两个半月形数据展开，使得数据集成为适用于某一线性分类器的输入数据。 PCA降维，映射到主成分，但仍不能很好地进行线性分类。

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'))
1
2
3
4
5
6
7
8
9
10
11
12
13

利用基于RBF核的KPCA来实现线性可分。可以看到，两个类别（圆形和三角形）此时是线性可分的，这使得转换后的数据适合作为线性分类器的训练数据集。

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()
1
2
3
4

分离同心圆，生成同心圆数据。再进行标准的PCA映射，最后得到KPCA的降维效果图。

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()

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

通过实验可以发现，仔细观察PCA其实只对原始数据进行了旋转操作，这是由于其寻找的是数据的“主要分布方向”。从图中看出，在重构后的数据清晰的分为两类，且在主成分方向上分割明显，便于使用线性分类方法进行处理。KPCA可以将原始数据投影至线性可分情况，打破了数据的原有结构。但是KPCA不能将数据投影至完全线性可分的程度，这说明KPCA只是个无监督的降维算法，它不管样本的类别属性，只是降维而已。

总结

主成分分析是一种降维方法。在PCA中，数据从原来的坐标系转换到了新的坐标系，新坐标系由数据本身决定。在新坐标系中，第一个坐标轴选择的是原始数据中方差最大的方向，第二个坐标轴选择的是和第一个坐标轴正交且具有最大方差的方向。该过程一直重复，重复次数为原始数据中特征的数目。我们会发现，大部分方差都包含在最前面的几个新坐标轴中。因此，我们可以忽略余下的坐标轴，即对数据进行了降维处理。其优点是降低数据的复杂性，识别最重要的特征，但其缺点在于不一定需要，且可能损失有用信息。

奇异值分解，是一种矩阵因式分解。通过计算奇异值个数和奇异向量，生成一个可以代替原矩阵的近似矩阵，将数据集的奇异值表征按重要性排列，舍弃不重要的特征向量。可用来达到降维的目的，从而找出数据中的主成分。其优点在于并行化，缺点是计算量大，分解后的矩阵解释性较弱。我们说到PCA降维时，需要找到样本协方差矩阵X^ TX最大的d个特征向量，然后用着最大的d个特征向量组成的矩阵来做低维投影降维。可以看出，在这个过程中需要先求出协方差矩阵X^ TX，当样本数多、样本特征数也多的时候，比如10000*10000的矩阵，这个计算量是很大的。注意到SVD也可以求出协方差矩阵X^ TX最大的d个特征向量组成的矩阵，但是SVD有个好处，就是可以不求出协方差矩阵X^TX，也能通过某些算法求出右奇异矩阵V。也就是说，PCA算法可以不用做特征分解，而是用SVD来进行完成。

核主成分分析降维，为了更好地处理非线性数据，引入非线性映射函数，将原空间中的数据映射到高维空间，注意，这个是隐性的，利用核函数进行处理，无需知道映射函数的具体形式，它让降维变得更有意义