机器学习笔记之（6）——数据降维（PCA与SVD）

数据降维（数据压缩）是属于非监督学习的一种，但是其实它也属于一种数据处理的手段。也就是说，通过数据降维，对输入的数据进行降维处理，由此剔除数据中的噪声并通过机器学习算法的性能，用于数据预处理。主要有：主成分分析（PCA）和奇异值分解（SVD）。PCA——通过按照数据方差最大方向调整数据的逼近来达到降维的目的。SVD——是矩阵分解技术中的一种，通过对原始数据的逼近来达到降维的目的。（之前介绍过的监督学习Fisher分类器也可以实现数据的降维）

在低维下，数据更容易进行处理。另外，其相关特征可能在数据中明确地显示出来。也可以起到去除噪声的作用。

使用特征抽取来减少数据集中特征的数量。特征抽取算法会将数据转换或映射到一个新的特征空间。特征抽取可以理解为：在尽可能多地保持相关信息的情况下，对数据进行压缩的一种方法。特征抽取通常用于提高计算效率，同样也可以帮助我们降低“维数灾难”——尤其当模型不适用于正则化处理时。

所谓的降维就是指采用某种映射方法，将原高维空间中的数据点映射到低维度的空间中。在原始的高维空间中，包含冗余信息以及噪声信息，在实际应用中例如图像识别中会造成误差，降低准确率；而通过降维，可以减少信息的冗余、提高识别等的精度，还可以寻求数据内部的本质结构特征。

如果特征维度成千上万，则需要的训练样本的数量几乎不可能满足。而且高维空间距离计算也比较麻烦（大部分机器学习算法都是基于特征距离的计算）。在高维度情形下出现的数据样本稀疏、距离计算困难等问题是所有机器学习方法共同面临的难题，称为“维度灾难”。可以通过数据降维来解决。

PCA：

主成分分析（principal component analysis）。在PCA中，数据从原来的坐标系转换到了新的坐标系，新坐标系的选择是由数据本身决定的。第一个新坐标轴选择的是原始数据中方差最大的方向，第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向。该过程一直重复，重复的次数为原始数据中特征的数目。我们会发现，大部分方差都包含在最前面的几个新坐标轴中。因此，可以忽略余下的坐标轴，即对数据进行了降维处理。

第一个主成分就是从数据差异性最大（即方差最大）的方向提取出来的，第二个主成分则来自于数据差异性次大的方向，并且该方向与第一个主成分方向正交。通过对数据集的协方差矩阵及其特征值分析，可以求得这些主成分的值。一旦得到了协方差矩阵的特征向量，我们就可以保留最大的N个值。这些特征向量也给出了N个最重要特征的真实结构。可以通过将数据乘上这N个特征向量而将它转换到新的空间。

PCA的目标是在高维数据中找到最大方差的方向，并将数据映射到一个维度不大于原始数据的新的子空间上。新特征的坐标是相互正交的，且新的子空间上正交的坐标轴（主成分）可以被解释为方差最大的方向。

对于d维原始特征空间，通过构建一个维度为d*k维的转换矩阵，将样本映射到新的k维特征空间上去（k<d）。

PCA算法步骤如下：

1、对原始d维数据做预处理。包括划分数据集与标准化处理。由于主成分的方向对数据值的范围高度敏感，因此要先对特征值做标准化处理。

2、构造样本的协方差矩阵。d*d维协方差矩阵是沿主对角线对称的，其中d为数据集的维度，此矩阵成对地存储了不同特征之间的协方差。如两个特征和，可通过如下公式计算它们之间的协方差：

其中，和分别为特征和的均值。若对数据集进行了标准化处理，则样本的均值会为0.两个特征之间的协方差如果为正，则它们会同时增减，而一个负的协方差则表示两个特征会朝相反的方向变动。一个包含三个特征的协方差矩阵可以记为：

3、计算协方差矩阵的特征值和相应的特征向量。协方差矩阵的特征向量代表主成分（最大方差方向），而对应的特征值大小就决定了特征向量的重要性。协方差矩阵的特征对计算如下：

特征值为一个标量。

4、选择与前k个最大特征值对应的特征向量，其中k为新特征空间的维度（k<=d）。由于要将数据集压缩到一个新的特征子空间上来实现数据降维，所以只选择那些包含最多信息（方差最大）的特征向量（主成分）组成子集。由于特征值的大小决定了特征向量的重要性，因此需要将特征值按降序排列，感兴趣的是排序在前k个的特征值所对应的特征向量。

特征值的方差贡献率：特征值与所有特征值和的比值：

5、通过前k个特征向量构建映射矩阵W。

6、通过映射矩阵W将d维的输入数据集X转换到新的k维特征子空间。将数据集中的信息转换到新的主成分轴上。

下面给出上述过程的python代码：

##############################数据的读入、划分、标准化##########################
import pandas as pd
#Python Data Analysis Library 或 pandas 是基于NumPy 的一种工具，该工具是为了解决数据分析任务而创建的。
#Pandas 纳入了大量库和一些标准的数据模型，提供了高效地操作大型数据集所需的工具。
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',header=None)
#读入csv文件，形成一个数据框
#使用葡萄酒数据集


from sklearn.cross_validation import train_test_split
x, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
x_train, x_test, y_train, y_test =train_test_split(x, y, test_size=0.3,stratify=y,random_state=0)
#葡萄酒数据集划分为训练集和测试集，分别占据数据集的70%和30%

#使用单位方差对数据集进行标准化
from sklearn.preprocessing import StandardScaler
sc=StandardScaler()
#注意，当改为from sklearn.preprocessing import StandardScaler as sc时，会报错
# fit_transform() missing 1 required positional argument: 'X'。这就是由于sc没有StandardScaler()而仅仅是StandardScaler
x_train_std=sc.fit_transform(x_train)
x_test_std=sc.fit_transform(x_test)



##############################构造协方差矩阵、获得特征值与特征向量###############
import numpy as np
cov_mat=np.cov(x_train_std.T)#求数据集的协方差矩阵
#进行特征分解
eigen_vals,eigen_vecs=np.linalg.eig(cov_mat)#用numpy中的linalg.eig函数来计算数据集协方差矩阵的特征对
print('\n特征值 \n%s' % eigen_vals)
#对于此处导入的葡萄酒数据集，可以获得13*13的协方差矩阵和13个特征向量及其对应的特征值。
#eigen_vals为存放了13个特征值的向量。而特征向量以列的方式存储于一个13*13维的矩阵中eigen_vecs。


#########################绘制特征值的方差贡献率图像##############################
tot=sum(eigen_vals)
var_exp=[(i/tot) for i in sorted(eigen_vals,reverse=True)]#用sorted函数进行排序，并用for循环求出各个特征值的方差贡献率，组成一个列表
#使用numpy的cumsum函数计算累计方差，并通过matplotlib的step函数绘制
cum_var_exp=np.cumsum(var_exp)

#画图
import matplotlib.pyplot as plt
###########下面两句使得图片内可以出现中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
###########
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',label= "单个方差贡献")
plt.step(range(1, 14), cum_var_exp, where='mid',label='累计方差贡献')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.tight_layout()
plt.show()


##########################################特征转换#############################
#按特征值的降序排列特征对：
eigen_pairs=[(np.abs(eigen_vals[i]),eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)
#接下来选取两个对应的特征值最大的特征向量，这两个值之和占据数据集总体方差的60%（可见特征值的方差贡献率图像）
w=np.hstack((eigen_pairs[0][1][:, np.newaxis],eigen_pairs[1][1][:, np.newaxis]))
#np.hstack用于把几个小数组合并成一个大数组,hstack表示轴1合并。
#hstack的字母h来自于horizontal，表示两个数组是水平的，hstack((a,b))将把b排在a的右边的意思。
#np.newaxis 在使用和功能上等价于 None，其实就是 None 的一个别名。
print('\nW矩阵 \n%s' % w)

print('\n将一个样本转换的PCA的子空间上，得到： \n%s' % x_train_std[0].dot(w))
print('\n样本转换前为： \n%s' % x_train_std[0])

##将数据集进行特征转换
x_train_pca=x_train_std.dot(w)
#画图
X_train_pca = x_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(x_train_pca[y_train == l, 0],X_train_pca[y_train == l, 1],c=c, label=l, marker=m)

plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

结果如下图所示：

上面是给出PCA过程的Python代码，接下来看一下scikit-learn中的PCA类（建议与上面程序结合着理解）

#####################################定义一个画图函数###########################
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], 
                    y=X[y == cl, 1],
                    alpha=0.6, 
                    c=cmap(idx),
                    edgecolor='black',
                    marker=markers[idx], 
                    label=cl)


##############################数据的读入、划分、标准化###########################
import pandas as pd
#Python Data Analysis Library 或 pandas 是基于NumPy 的一种工具，该工具是为了解决数据分析任务而创建的。
#Pandas 纳入了大量库和一些标准的数据模型，提供了高效地操作大型数据集所需的工具。
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',header=None)
#读入csv文件，形成一个数据框
#使用葡萄酒数据集


from sklearn.cross_validation import train_test_split
x, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
x_train, x_test, y_train, y_test =train_test_split(x, y, test_size=0.3,stratify=y,random_state=0)
#葡萄酒数据集划分为训练集和测试集，分别占据数据集的70%和30%

#使用单位方差对数据集进行标准化
from sklearn.preprocessing import StandardScaler
sc=StandardScaler()
#注意，当改为from sklearn.preprocessing import StandardScaler as sc时，会报错
# fit_transform() missing 1 required positional argument: 'X'。这就是由于sc没有StandardScaler()而仅仅是StandardScaler
x_train_std=sc.fit_transform(x_train)
x_test_std=sc.fit_transform(x_test)


#####################################使用sklearn中的PCA类#######################
from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
pca=PCA(n_components=2)#类似于上面的程序，选取两个最大的特征值作为主成分
#做降维处理
x_train_pca=pca.fit_transform(x_train_std)
x_test_pca=pca.transform(x_test_std)

#用逻辑斯蒂回归进行分类(对训练数据进行处理)
lr=LogisticRegression()
lr.fit(x_train_pca,y_train)
plot_decision_regions(x_train_pca,y_train,classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

#用逻辑斯蒂回归进行分类(对测试数据进行处理)
lr=LogisticRegression()
lr.fit(x_train_pca,y_train)
plot_decision_regions(x_test_pca,y_test,classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

结果如下图所示：

比较PCA类与自行实现的PCA分析结果，可以发现，左图为右图沿Y轴对称反转获得。这仅仅是由于特征分析方法不一样而导致，本质上没有区别。

当将主成分由2改为None时，可得到：

当将n_components设置为None时，将按照方差贡献率递减的顺序返回所有主成分。

核主成分分析：

许多机器学习算法都是假定输入数据是线性可分的。降维处理也一样，PCA和LDA都是针对线性可分数据降维。而对于非线性可分的数据要进行降维处理可以使用核技巧的PCA，或称为核PCA。通过使用核PCA可以将非线性可分的数据转换线性可分的低维子空间上（在机器学习笔记5中的核SVM介绍过类似的概念）

线性与非线性可分问题图示如下：

由于本人对于核函数与核技巧理解还不是很深入，下面直接给出书本上的照片（后面有新的体会会及时更新本博文）：

下面给出Python代码：

################################基于RBF核的PCA##################################
from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):#RBF核函数结合着理论部分来看
    """
    RBF kernel PCA implementation.

    Parameters
    ------------
    X: {NumPy ndarray}, shape = [n_samples, n_features]
        
    gamma: float
      Tuning parameter of the RBF kernel
        
    n_components: int
      Number of principal components to return

    Returns
    ------------
     X_pc: {NumPy ndarray}, shape = [n_samples, k_features]
       Projected dataset   

    """
    # Calculate pairwise squared Euclidean distances
    # in the MxN dimensional dataset.
    sq_dists = pdist(X, 'sqeuclidean')

    # Convert pairwise distances into a square matrix.
    mat_sq_dists = squareform(sq_dists)

    # Compute the symmetric kernel matrix.
    K = exp(-gamma * mat_sq_dists)

    # Center the kernel matrix.
    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)

    # Obtaining eigenpairs from the centered kernel matrix
    # scipy.linalg.eigh returns them in ascending order
    eigvals, eigvecs = eigh(K)
    eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]

    # Collect the top k eigenvectors (projected samples)
    X_pc = np.column_stack((eigvecs[:, i]
                            for i in range(n_components)))

    return X_pc


################################分离半月形数据##################################
from sklearn.datasets import make_moons
x, y = make_moons(n_samples=100, random_state=123)

#画出半月形数据
import matplotlib.pyplot as plt
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.tight_layout()
plt.show()   


#########################先用标准的PCA来进行降维处理#############################
from sklearn.decomposition import PCA
scikit_pca = PCA(n_components=2)
x_spca = scikit_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='o', 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.tight_layout()
plt.show()


#########################使用核PCA来进行降维处理#################################
x_kpca = rbf_kernel_pca(x, gamma=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='o', 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.tight_layout()
plt.show()


################################分离同心圆数据##################################
from sklearn.datasets import make_circles
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.tight_layout()
plt.show()


#########################先用标准的PCA来进行降维处理#############################
scikit_pca = PCA(n_components=2)
X_spca = scikit_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='o', 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.tight_layout()
plt.show()


#########################使用核PCA来进行降维处理#################################
X_kpca = rbf_kernel_pca(X, gamma=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='o', 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.tight_layout()
plt.show()

结果如下图所示：



下面给出使用sklearn中核PCA类来进行处理的Python代码：

################################分离半月形数据##################################
from sklearn.datasets import make_moons
x, y = make_moons(n_samples=100, random_state=123)

#画出半月形数据
import matplotlib.pyplot as plt
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.tight_layout()
plt.show()   


#########################核PCA来进行降维处理#############################
from sklearn.decomposition import KernelPCA
KPCA=KernelPCA(n_components=2,kernel='rbf',gamma=15)
x_KPCA=KPCA.fit_transform(x)


#画出结果
import matplotlib.pyplot as plt
import numpy as np

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='o', 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.tight_layout()
plt.show()

结果如下图所示：

SVD：

奇异值分解（singular value decomposition，SVD）是一种矩阵分解的技术，是强大的信息提取的工具。SVD的应用场景：隐性语义索引、推荐系统。SVD的数学推导见下图片：

SVD将原始的数据集矩阵Data分解成为如下式所示的形式：

上述分解会构建出一个矩阵，该矩阵只有对角元素，其他元素均为0。另一个惯例就是，该矩阵的对角元素是从大到小排列的。这些对角元素称为奇异值。对应的就是原始数据集矩阵的奇异值。

奇异值就是矩阵特征值的平方根。

由于奇异值分解后的矩阵只有从大到小排列的对角元素。在科学和工程中，一直认为：在某个奇异值的数目（r个）之后，其他的奇异值都置为0。这就意味着数据集中仅有r个重要的特征，而其余特征则都是噪声或冗余特征。

在奇异值跳变的地方选取奇异值（或者，保留矩阵中90%的能量信息。为了计算总能量，可以将所有奇异值求其平方和。于是可以将奇异值的平方和累加到总值的90%为止），那么原始的矩阵就可以用以下结果来近似：

Python中实现奇异值分解的函数如下图所示

此处的Sigma以行向量的形式返回，但实际上它是一个矩阵，只是为了节省空间而已~

给出Python代码如下：

#######################################定义一个数据矩阵#########################
def loadData():
    return[[0, 0, 0, 2, 2],
           [0, 0, 0, 3, 3],
           [0, 0, 0, 1, 1],
           [1, 1, 1, 0, 0],
           [2, 2, 2, 0, 0],
           [5, 5, 5, 0, 0],
           [1, 1, 1, 0, 0]]


##########################################对上面矩阵进行奇异值分解###############
Data=loadData()
import numpy as np
U,Sigma,VT=np.linalg.svd(Data)
print("奇异值")
print(Sigma)#通过打印的结果可以看出，前两个奇异值比后三个奇异值大很多，故此把后三个去掉（在奇异值跳变的地方选）

######################################重构原始矩阵##############################
Sig=np.mat([[Sigma[0],0],[0,Sigma[1]]])
Data1=U[:,:2]*Sig*VT[:2,:]
print("奇异值分解降维后矩阵")
print(Data1)

结果如下图所示：

感觉对于奇异值分解降维理解还不够深入，接下来我更深入的体会会继续更新本博文

当然~除了以上降维方法以外，还有流形学习降维、多维缩放降维、等度量映射降维、局部线性嵌入等。

主要参考资料如下：

• 《机器学习实战》
• 《Python机器学习》
• 《Python大战机器学习》