降维算法之PCA及其实战

      主成分分析（Principal components analysis，以下简称PCA）是最重要的降维方法之一。在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用。为什么说PCA应用广泛呢？PCA是一种无监督算法，也就是我们不需要标签也能对数据做降维，这就使得其应用范围更加广泛了。一般我们提到降维最容易想到的算法就是PCA，下面我们就对PCA的原理做一个总结。
1.PCA的思想
      PCA就是找出数据里最主要的方面，用数据里最主要的方面来代替原始数据。具体可以这么理解，假如我们的数据集是n维的，共有m个数据。我们希望这个m个数据的维度从n维降到n’维，希望这m个n’维的数据集尽可能地代表原始数据。当然，毕竟降低了维度，所以数据多少都会产生损失，我们能做的就是把损失降得越低越好。但是PCA也有一个问题，原来的数据中比如包括了年龄，性别，身高等指标降维后的数据既然维度变小了，那么每一维都是什么含义呢？这个就很难解释了，所以PCA本质来说是无法解释降维后的数据的物理含义，换句话说就是降维之后，计算机能更好地认识这些数据，但是我们就很难理解它们了。
      我们先看看最简单的情况，也就是n=2，n’=1,也就是将数据从二维降维到一维。数据如下图。我们希望找到某一个维度方向，它可以代表这两个维度的数据。图中列了两个向量方向，u1和u2，那么哪个向量可以更好的代表原始数据集呢？从直观上也可以看出，u1比u2好。

      为什么u1比u2好呢？可以有两种解释，第一种解释是样本点到这个直线的距离足够近，第二种解释是样本点在这个直线上的投影能尽可能的分开。所以我们引出降维标准：样本点到这个超平面的距离足够近（最近重构性：基于最小投影距离）或者说样本点在这个超平面上的投影能尽可能的分开（最大可分性：基于最大投影方差）。
2.PCA数学原理推导（最大可分性）
      内积： $\left( {{a}_{1}},{{a}_{2}},...,{{a}_{n}} \right)\cdot {{\left( {{b}_{1}},{{b}_{2}},...,{{b}_{n}} \right)}^{T}}={{a}_{1}}{{b}_{1}}+{{a}_{2}}{{b}_{2}}+...+{{a}_{n}}{{b}_{n}}$(a1​,a2​,...,an​)⋅(b1​,b2​,...,bn​)T=a1​b1​+a2​b2​+...+an​bn​
$A\cdot B=\left| A \right|\left| B \right|\cos \left( a \right)$A⋅B=∣A∣∣B∣cos(a)
      设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度。

      基是正交的（即内积为0，或者说相互垂直）。要求基线性无关。如下图蓝线所示:

      基变换：数据与一个基做内积运算，结果作为第一个新的坐标分量，然后与第二个基做内积运算，结果作为第二个新坐标的分量。
$\left( \begin{matrix} {{p}_{1}} \\ {{p}_{2}} \\ ... \\ {{p}_{n'}} \\ \end{matrix} \right)\left( \begin{matrix} {{a}_{1}} & {{a}_{2}} & ... & {{a}_{m}} \\ \end{matrix} \right)=\left( \begin{matrix} {{p}_{1}}{{a}_{1}} & {{p}_{1}}{{a}_{2}} & ... & {{p}_{1}}{{a}_{m}} \\ {{p}_{2}}{{a}_{1}} & {{p}_{2}}{{a}_{2}} & ... & {{p}_{2}}{{a}_{m}} \\ ... & ... & ... & ... \\ {{p}_{n'}}{{a}_{1}} & {{p}_{n'}}{{a}_{2}} & ... & {{p}_{n'}}{{a}_{m}} \\ \end{matrix} \right)$⎝⎜⎜⎛​p1​p2​...pn′​​⎠⎟⎟⎞​(a1​​a2​​...​am​​)=⎝⎜⎜⎛​p1​a1​p2​a1​...pn′​a1​​p1​a2​p2​a2​...pn′​a2​​............​p1​am​p2​am​...pn′​am​​⎠⎟⎟⎞​
      将a矩阵中的每一列列向量变换到p矩阵中每一行行向量为基所表示的空间中去。
      例如：数据(3,2)映射到基中坐标： $\left( \begin{matrix} \begin{matrix} 1/\sqrt{2} & 1/\sqrt{2} \\ \end{matrix} \\ \begin{matrix} -1/\sqrt{2} & -1/\sqrt{2} \\ \end{matrix} \\ \end{matrix} \right)\left( \begin{matrix} 3 \\ 2 \\ \end{matrix} \right)=\left( \begin{matrix} 5/\sqrt{2} \\ -1/\sqrt{2} \\ \end{matrix} \right)$(1/2​​1/2​​−1/2​​−1/2​​​)(32​)=(5/2​−1/2​​)
      首先将数据进行中心化处理（均值为0，方差为1）。寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。
      协方差： $\operatorname{cov}\left( a,b \right)=\frac{1}{m-1}\sum\limits_{i=1}^{m}{\left( {{a}_{i}}-0 \right)\left( {{b}_{i}}-0 \right)}$cov(a,b)=m−11​i=1∑m​(ai​−0)(bi​−0)
      如果单纯只选择方差最大的方向，后续方向应该会和方差最大的方向接近重合。为了让后续方向尽量分开，我们希望所有方向线性无关，从而引入协方差来进行计算。当协方差为0时，表示两个字段完全独立。为了让协方差为0，选择第二个基时，只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。
      优化目标： $\underset{W}{\mathop{\max }}\,{{W}^{T}}X{{X}^{T}}W$Wmax​WTXXTW
$s.t.{{W}^{T}}W=I$s.t.WTW=I
      求最大化，我们是不是又想到了拉格朗日乘子法？建立新的式子为：$c\left( W \right)={{W}^{T}}X{{X}^{T}}W-\lambda \left( {{W}^{T}}W-1 \right)$c(W)=WTXXTW−λ(WTW−1)
      对W进行求导： $\frac{dc}{dW}=X{{X}^{T}}W-\lambda W\Rightarrow X{{X}^{T}}W=\lambda W$dWdc​=XXTW−λW⇒XXTW=λW
      从而只需要求出原始数据的协方差即可矩阵即可。
      将一组n维向量降为n’维（n’大于0且小于n）目标选择n’个单位正交基，使原始数据变换到这组基上后，各个字段两两之间协方差为0，字段的方差则尽可能大。
      协方差矩阵： $X=\left( \begin{matrix} \begin{matrix} {{a}_{1}} & {{a}_{2}} & ... & {{a}_{n'}} \\ \end{matrix} \\ \begin{matrix} {{b}_{1}} & {{b}_{1}} & ... & {{b}_{n'}} \\ \end{matrix} \\ \end{matrix} \right),\frac{1}{n'}X{{X}^{T}}=\left( \begin{matrix} \begin{matrix} \frac{1}{n'}\sum\limits_{i=1}^{n'}{a_{i}^{2}} & \frac{1}{n'}\sum\limits_{i=1}^{n'}{{{a}_{i}}{{b}_{i}}} \\ \end{matrix} \\ \begin{matrix} \frac{1}{n'}\sum\limits_{i=1}^{n'}{{{a}_{i}}{{b}_{i}}} & \frac{1}{n'}\sum\limits_{i=1}^{n'}{b_{i}^{2}} \\ \end{matrix} \\ \end{matrix} \right)$X=(a1​​a2​​...​an′​​b1​​b1​​...​bn′​​​),n′1​XXT=⎝⎜⎜⎛​n′1​i=1∑n′​ai2​​n′1​i=1∑n′​ai​bi​​n′1​i=1∑n′​ai​bi​​n′1​i=1∑n′​bi2​​​⎠⎟⎟⎞​
      矩阵对角线上的两个元素分别是两个字段的方差，而其他元素是a和b的协方差。
      协方差矩阵对角化：即除对角线外的其他元素化为0，并且在对角线上将元素按大小从上到下排列。
$PC{{P}^{T}}=\Lambda =\left( \begin{matrix} {{\lambda }_{1}} & {} & {} & {} \\ {} & {{\lambda }_{2}} & {} & {} \\ {} & {} & ... & {} \\ {} & {} & {} & {{\lambda }_{n}} \\ \end{matrix} \right)$PCPT=Λ=⎝⎜⎜⎛​λ1​​λ2​​...​λn​​⎠⎟⎟⎞​
      实对称矩阵：一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量。
      根据特征值的绝对值从大到小，将特征向量从上到下排列，则用前n’行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。
      实例如下图所示：

3.实战
3.1获取鸢尾花数据集

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
df = pd.read_csv('iris.data')
# 修改列属性名称
df.columns=['sepal_len', 'sepal_wid', 'petal_len', 'petal_wid', 'class']
# 将数据进行分离，特征数据X，标签数据y
X = df.iloc[:,0:4].values
y = df.iloc[:,4].values

3.2四个特征对应与三个标签，观察一下分布情况

label_dict = {1: 'Iris-Setosa',
              2: 'Iris-Versicolor',
              3: 'Iris-Virgnica'}

feature_dict = {0: 'sepal length [cm]',
                1: 'sepal width [cm]',
                2: 'petal length [cm]',
                3: 'petal width [cm]'}
plt.figure(figsize=(8,6))
for cnt in range(4):
    plt.subplot(2,2,cnt+1)
    for lab in ('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'):
        plt.hist(X[y == lab,cnt],label=lab,bins=10,alpha=0.3)
        plt.xlabel(feature_dict[cnt])
        plt.legend(loc='upper right',fancybox=True,fontsize=8)

plt.tight_layout()
plt.show()

3.3数据标准化之后，计算协方差、特征值和特征向量
（1）标准化数据

from sklearn.preprocessing import StandardScaler
X_std = StandardScaler().fit_transform(X)
# 求每列特征的均值
mean_vec = np.mean(X_std,axis=0)

（2）求协方差

cov_mat = (X_std - mean_vec).T.dot(X_std - mean_vec) / (X_std.shape[0]-1)
print('Covariance matrix \n%s' %cov_mat)

输出结果：

Covariance matrix 
[[ 1.00675676 -0.10448539  0.87716999  0.82249094]
 [-0.10448539  1.00675676 -0.41802325 -0.35310295]
 [ 0.87716999 -0.41802325  1.00675676  0.96881642]
 [ 0.82249094 -0.35310295  0.96881642  1.00675676]]

Numpy科学计算包求协方差

print('NumPy covariance matrix: \n%s' %np.cov(X_std.T))

输出结果：

NumPy covariance matrix: 
[[ 1.00675676 -0.10448539  0.87716999  0.82249094]
 [-0.10448539  1.00675676 -0.41802325 -0.35310295]
 [ 0.87716999 -0.41802325  1.00675676  0.96881642]
 [ 0.82249094 -0.35310295  0.96881642  1.00675676]]

（3）求特征值与特征向量

cov_mat = np.cov(X_std.T)
eig_vals,eig_vecs = np.linalg.eig(cov_mat)
# 构造特征值与特征向量的元组
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
# 计算特征值之和
tot = sum(eig_vals)
# 计算每个特征值的重要程度
var_exp = [(i / tot)*100 for i in sorted(eig_vals,reverse = True)]
print(var_exp)
cum_var_exp = np.cumsum(var_exp)
cum_var_exp

输出结果：

[72.6200333269203, 23.14740685864416, 3.711515564584526, 0.521044249851025]
array([ 72.62003333,  95.76744019,  99.47895575, 100.        ])

画出特征值重要程度图

plt.figure(figsize=(6,4))
plt.bar(range(4),var_exp,align='center',label='individual explained variance')
plt.step(range(4),cum_var_exp,where='mid',label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

# 选择最大的两个特征值所对应的特征向量
matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1),eig_pairs[1][1].reshape(4,1)))

3．4降维

Y = X_std.dot(matrix_w)

（1）未进行降维的原始据对应的散点图

plt.figure(figsize=(6, 4))
for lab, col in zip(('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'),
                        ('blue', 'red', 'green')):
     plt.scatter(X[y==lab, 0],
                X[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('sepal_len')
plt.ylabel('sepal_wid')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

（2）降维之后的数据对应的散点图

plt.figure(figsize=(6, 4))
for lab, col in zip(('Iris-setosa', 'Iris-versicolor', 'Iris-virginica'),
                        ('blue', 'red', 'green')):
     plt.scatter(Y[y==lab, 0],
                Y[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend(loc='lower center')
plt.tight_layout()
plt.show()

      通过上面的两个图可以看出，降维在一定程度上可以使得类的界线更加明显。

参考文档：
周志华《机器学习》
https://www.cnblogs.com/pinard/p/6239403.html
https://blog.csdn.net/tangyudi/article/details/80188302