前一篇提到的人脸识别中，我们在使用SVM支持向量机做人脸分类之前使用到PCA提取人脸数据中的主要成分，降低计算的维度，那么具体PCA是如何提取的呢？下文了解一下。

PCA is a method to project data in a higher dimensional space into lower dimensional space by maximizing the variance of each dimension --wiki

PCA is mostly used as a tool in exploratory data analysis and for making predictive models. It is often used to visualize genetic distance and relatedness between populations. PCA is either done by singular value decomposition of a design matrix or by doing the following 2 steps:

1. calculating the data covariance( or correlation) matrix of the original data

2. performing eigenvalue decomposition(特征值分解) on the covariance matrix(协方差矩阵).

--wiki

主成分分析，是一种统计方法，通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。对于这些成分，有些就显得没有那么必要。那么我们可以从这些成分中挑取比较有用的，而对于我们模型的训练作用不大的可以去除掉。这个过程其实也是一个降维的过程，对高维度的数据集进行降维，减少一些信息含量，保留主要的那些维度。主要的那些维度是指在这些维度上具有更大的差异，而对于那些差异化不大的特征，有时候反而会影响我们的模型学习，而成为一种噪音，也会影响我们training时候的速度。所以在机器学习的过程中，去掉那些无关紧要的变量会使得我们的学习模型效率提高。

因此我们要找到那些差异化比较大的特征，也就是说在表示它的特征的这一系列数中，它的variance（方差：是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量）是比较大的，比较离散的，而对于那些不明显的特征，用数据来衡量他们的这种特征程度的话，他们更容易聚集在一起，而不那么离散，也就是variance的值很小。下面使用代码画图来形象地了解一下。

%matplotlib inline 
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()

rng=np.random.RandomState(1)
X=np.dot(rng.rand(2,2),rng.rand(2,200)).T
plt.scatter(X[:,0],X[:,1])
plt.axis('equal')

先画出散点图，观察一下我们的samples：

Numpy.dot(a,b,out=None): do product of two arrays.

就是两个矩阵的乘积，A矩阵的第一行乘以B矩阵的第一列所得的元素的和作为新矩阵C的第一列中第一行的数，A矩阵第二行乘以B矩阵第一列所得元素的和作为C矩阵的第一列中第二行的数，以此类推。

Numpy.random.RandomState.rand(2，200) random values in a given shape, 产生一个2x200的矩阵。矩阵中的每一个element（元）都是随机产生的介于0到1之间的一个数。

X=np.dot(rng.rand(2,2),rng.rand(2,200)).T

这个代码的意思是先随机产生两个点，然后再用这两个点进行一个列变换（左乘行变换，右乘列变换）产生与左边矩阵具有相同行，与右乘的矩阵右相同列的shape为2x200的新矩阵，可以知道这个新的矩阵中的每个点的x或y坐标都不会大于2（1x1+1x1=2）。尾部的.T进行了矩阵的转置，使得矩阵从2x200的shape变成200x2。

使用PCA来解出上图中散点的主成分，图中每个点需要两个参数（x,y坐标）才能表达准确，那么对于该图上的点来说，number of components=2，有两个成分。

from sklearn.decomposition import PCA
pca=PCA(n_components=2)
pca.fit(X)

看起来只有三行代码，但是相关的计算几乎全在里面了。

Sklearn.decomposition.PCA(n_components=None, *, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', random_state=None)

PCA: principal component analysis(PCA)

上述代码中的主成分分析使用的是SVD（singular value decomposition,奇异值分解）方法。

Linear dimensionality reduction using Singular Value Decomposition of the data to project it to a lower dimensional space.

PCA skips less significant component. Obviously, we can use SVD to find PCA by truncating the less important basis vectors in the original SVD matrix.

PCA方法寻找主成分的过程，通过采用SVD找到所有的基础向量之后，再把不怎么重要的基础向量裁断掉。所以这里需要先使用SVD找到矩阵X的基础向量。

从Sklearn.decomposition.PCA()的原代码中抽出计算基础向量部分：

scipy.linalg.svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd')

Singular Value Decomposition. 奇异值分解。

A:(M,N) array_like, matrix to decompose，一个要分解的（m,n）矩阵。

SVD Return（返回值）:

U: unitary matrix having left singular vectors(左奇异矩阵) as columns. Of shape (M,M) or (M,K), depending on full_matrices.

S: ndarray, the singular values（奇异值）, sorted in non-increasing order. Of shape(K), with K=min(M,N).

Wh: Unitary matrix having right singular vectors(右奇异矩阵) as rows. Of shape (N,N) or (k,N).

任何一个矩阵A都可以拆分成如下的三个矩阵相乘。

（读作sigma）符号通常也用S来代替,下文中这两者通用哈。

上式中U和V分别是AAT和ATA的正交特征矩阵，S是拥有r个元素的对角矩阵，对角线上的元素的值等于AAT或ATA的正特征值的根，其中AAT和ATA拥有同样的正特征值。

S长下面这样：

加上维度完整地写出如下形式。

上面拆分过程中，A矩阵一次性可以做到的变换可以分为U、S、VT三个步骤完成，其每个步骤的作用如下图所示：

VT: orthogonal(rotation)， 正交化(旋转)；

Sigma: diagonal(scale)，对角化（调整比例）；

U：(orthogonal)(rotation)，正交化（旋转）。

其中U和V为正交矩阵，满足如下条件：

已经明了U,V,S大致的计算方法，下面使用具体的矩阵来计算体验一下这个过程，找点感觉。

假设A矩阵如下：

下面将A矩阵分解为U，V，S：

计算U时，因为U是AAT的正交特征矩阵，所以下面要先计算AAT，然后待求得特征值之后，再计算它的正交矩阵：

计算特征值：

得AAT的特征值为

代入特征值，并计算得特征向量和矩阵为：

由于V是ATA的正交特征矩阵，那么可以使用相同的方法计算V：

先求得ATA的特征值为

，

对应的特征向量和矩阵为：

由上面得到的U和V，可以通过下式求解得S:

把所求得的矩阵代回下式计算验证，最后等于A，没有问题。

至于为什么可以通过这样的方式计算得到左奇异矩阵，右奇异矩阵和奇异值，相关证明自己找文献哈。

下面对使用SVD方法找出的PCA结果画图查看一下。

def draw_vector(v1,v0,ax=None):
    ax=ax or plt.gca()
    arrowprops=dict(arrowstyle='<-',
                    linewidth=2,
                    shrinkA=0,shrinkB=0,color='Black')
    ax.annotate('',v1,v0,arrowprops=arrowprops)
    #  draw a FancyArrowPatch arrow between the positions xy and xytext(v1 and v0).

plt.scatter(X[:,0],X[:,1],alpha=0.3)
for length,vector in zip(pca.explained_variance_,pca.components_):
    v=vector*3*np.sqrt(length)
    draw_vector(pca.mean_,pca.mean_+v)
plt.axis('equal')

代码中相关参数解析：

Pca.components_:array,shape(n_components, n_features)

Components_=V

Principal axes in feature space, representing the direction of maximum variance in the data。

特征空间中的主轴（主要成分），展示数据中最大方差的方向，（要看你需要，需要两个主成分，n_components=2时就是前两大），确定的是主成分的方向。

Pca.explained_variance_:array,shape(n_components,)

Explained_variance_=(S**2)/(n_samples-1)

The amount of variance explained by each of the selected components.

选择的成分的方差的大小，确定的是对应的方向上的大小，也就是上图中两个黑色箭头的长短，排序是S中方差从大到小的排序。

Pca.mean_:array, shape(n_feature)

Per-feature empirical mean, estimated from the training set.Equal to X.mean,

计算出整个训练集中点群的中间点位置，即上图中两个黑色箭头的交叉点（箭头的起点）。

因为我们给定的n_components=2，所以上图中有两个向量，V确定了这两个矢量的方向，在计算时我们使用了sqrt(explained_variance_)，对explained_variance开了根号，而explained_variance是由S的平方除以（n_samples-1)得到的，所以这里的结果和S是线性的关系，从而确定了这两个矢量的大小，由此画出了上图中两个黑色的既有长度，又有大小的箭头。箭头越长说明图上的散点在这个方向上越分散，不同的点在这个方向上更能有较大的区别，而sigma值更小的，表现出来箭头也更短，所以长箭头表示的方向比短箭头表示的方向在这组数据中更为主要，图中两个黑色箭头互相垂直。

画箭头的function:

Ax.annotate(self, text, xy, *args, **kwargs)

Annotate the point xy with text text,

ax.annotate('',v1,v0,arrowprops=arrowprops)

直线是点xy 到xytext之间，即v1和v0这两个点连成一条直线，箭头的样式定义在arrowprops这个dictionary里面。

arrowprops=dict(arrowstyle='<-', linewidth=2, shrinkA=0, shrinkB=0, color='Black')

arrowstyle箭头的类型，arrowstyle变成’->’时，方向就会反过来。Linwidth:线条宽度，shrink：就是要不要缩短箭头，我们这里就不需要了，color：箭头颜色。

下面看看当我们把维度降下来，变成一个component的时候：

此时需要把n_components=1，之前的S里面有两个基础向量，现在把n_components设置为1后变成了一个。在S矩阵里面，从左上角到右下角，取出第一个sigma。

pca=PCA(n_components=1)
pca.fit(X)
X_pca=pca.transform(X)
print("original shape: ",X.shape)
print("transformed shape: ",X_pca.shape)

X_new=pca.inverse_transform(X_pca)
plt.scatter(X[:,0],X[:,1],alpha=0.3)
plt.scatter(X_new[:,0],X_new[:,1],alpha=0.8)
plt.axis('equal')

Pca.fit(X): fit the model.

Pca.transform(X) : Apply dimensionality reduction to X.

给原来的矩阵降维，这里设置的n_components=1，最终维度为1。

Pca.inverse_transform(X_pca): transform data back to its original space.

把降维后的数据还原到原来的空间中，这边是把二维降维到一维后的数据还原到二维空间中，还原到二维空间后，其中一个维度的信息被丢失就不能展示出来。

画出来的图像如上图，橘色部分的点，组成一条线条，和之前的较长边的向量（黑色箭头）具有相同的倾斜角度，并且蓝色的每一个点做垂线垂直于橘色点组成的直线的话，在垂足的地方有与蓝色点对应的橘色点，也就是蓝色点在这条直线上的投影。这里面在另外一个方向上，也就是较短的那个向量方向上的信息我们就丢掉了，inverse_transform之后的结果只有一个component inverse transform的结果，而在另外那个方向上，由于它的variance更小，我们选择丢掉它。

再来看一下更高维度的：

from sklearn.datasets import load_digits
digits=load_digits()
digits.data.shape

#shape 为: 1797,64，有64个维度。

def plot_digits(data):
    fig,axes=plt.subplots(4,10,figsize=(10,4),
                         subplot_kw={'xticks':[],'yticks':[]},
                         gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8,8),
                 cmap='binary',interpolation='nearest',
                 clim=(0,16))
plot_digits(digits.data)

画出来的图像如上图，现在人为给它加一点noise。

np.random.seed(42)
noisy=np.random.normal(digits.data,4)
plot_digits(noisy)

得到下面的图像：

然后只提取前面50%的variance，

pca=PCA(0.5).fit(noisy)
plt.plot(np.cumsum((pca.explained_variance_ratio_)))

画出上图所需要的相关计算如下：

所有components的方差求和：Total_var=explained_variance_.sum()

每个component的方差所占比例：Explained_variance_ratio_ = explained_variance_ / total_var

每个component的方差所占比例求和：Np.cumsum()。

上图可以看到对方差占比总和刚达50%时，components的数量为12个，也就是对方差按从大到小排名后，前面至少12个components的variance加起来才会占总的variance 50%以上，使用那12个components进行画图：

components=pca.transform(noisy)
filtered=pca.inverse_transform(components)
plot_digits(filtered)

上图看起来虽然没有加noise之前的图清晰，但是比增加了noise的那张图要清晰很多，这个提取主成分的过程过滤掉了很多不必要的成分，也就是一个过滤noise的过程，即使丢掉了一些细微的信息，但是图像的质量比有Noise那个清晰多了。