你见过最全的主成分分析PAC与梯度上升法总结

主成分分析一个非监督学习算法，主要用于数据降维，通过降维可以发现数据更容易理解的特征，其他作用也有可视化、降噪等。
假设现有样本的分布如图。

样本有两个特征，如果对样本进行降维，首先可以考虑基于坐标轴进行降维。会有如下两种方式。

分别对应基于特征2降维和基于特征1降维。就降维之后的效果而言，右图会更好一些，因为样本之间的间距比较大，样本之间有更好的可区分度。但右图这种方案不一定是最好的方案。

如果将样本映射到这样的一根斜线上效果会更好，因此，降维的目标就是找到样本间距最大的轴。也就是样本映射到这个轴之后，使得样本间距最大。

在数学上，样本间距离反映了样本的离散程度，而方差可以很好的表示这种程度，因此我们定义样本间距离——方差。

主成分分析PCA

PCA算法是最经典的数据降维算法，步骤大致如下：

①将样本的均值归零（demean）
②求一个轴的方向w，使得映射后的方差值最大

均值归零

原有样本分布

归零之后

归零的作用是使得样本的均值为0，这样就可以在计算方差时，将原有公式中的均值置为0，即获得图中均值归零后的方差表达式。
注意这里的Xi表示的是映射之后样本的数据，也就是降维之后的数据。

映射轴和方差

映射的轴可以视为一个向量，w = (w1, w2, ..., wn)，n即为样本的特征数。
方差公式

我们的目标就是使得这个公式得到的方差值最大。
如果将样本的数据映射到坐标系中，那么每个样本就可以看做是个向量，而样本的均值也可以看做是一个向量。这样我们的方差公式可以变为

也就是求每个向量和均值向量之间距离的平方和再除以向量个数，得到的结果本质上也是方差。由于前一步中的均值归零操作，使得样本在各个维度的均值都是0，那么进一步化简公式为

假设在数据有两个维度，则可以得到一个二维平面。

Xi是样本点，有Xi,1和Xi,2两个特征，假设我们求得的最大样本间距的映射轴是w表示的向量，设w=(w1, w2)，那么Xi映射到w上的点如图中蓝色线表示。
高中数学中，两个向量相乘等于两个向量的模乘以两个向量的夹角余弦值。
通常在运算中，我们会将向量w设为单位向量（即模长为1的向量），因此等式右侧的w的模长可以省略。
又因为Xi的模长乘以夹角余弦值正好是蓝色部分的长度，因此我们得到最后一个等式。
因此，我们对之前求方差的算式进行修改得到

这里其实不应使用取模的符号，因为Xi样本有n个属性，可以理解为一个n个元素的向量，w也是有n个元素的向量，两个向量相乘的结果是一个实数。而放到矩阵来说，Xi是一个1 * n的矩阵，w理解为一个n * 1的矩阵，相乘之后也是一个实数，因此，最佳的表达式应该是

展开来看

合并一下w各个值

得到这个算式最大值的过程可以采用梯度上升法求解。

梯度上升法

回忆梯度下降法，通过求函数在一个点的导数值，得到函数在该点的斜率，之后使用减去斜率与学习率乘积的方式下降，逐渐靠近极小值点。
梯度上升法与之类似，上升不过就是将减法变为加法。第一步还是对方差函数f在每个特征维度求导。

对括号中的式子化简

提出每个式子的Xi * w

此时得到一个向量和X样本矩阵，由于向量w和X中每个样本都有相乘的操作，因此，可以将其写为2/m * (Xw)，左侧算式可以理解为Xi * w需要和每个样本的第i列做一次乘法，因此，可以认为最终是要乘以整个样本X的矩阵。因此得到最右侧的两个式子，下面的算式是为了方便运算而对上面的式子进行转换得到的。

主成分分析

单个主成分分析

先模拟一个数据集

import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10, size=100)

plt.scatter(X[:, 0], X[:, 1])
plt.show()

该数据集有两个属性，打印初始图像

接下来对数据进行主成分分析，第一步是均值归零，定义相应的函数

'''均值归零'''
def demean(X):
    # axis=0即按列计算均值，及每个属性的均值，1则是计算行的均值
    return (X - np.mean(X, axis=0))

X_demean = demean(X)
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.show()

计算均值调用的是numpy自带的函数，均值归零之后打印归零后的图像

分布没有变化，但从坐标轴可以看出均值归零后的效果。
接下来就是对方差函数和其导数函数的定义

'''方差函数'''
def f(w, X):
    return np.sum((X.dot(w)**2)) / len(X)

'''方差函数导数'''
def df_math(w, X):
    return X.T.dot(X.dot(w)) * 2. / len(X)

'''将向量化简为单位向量'''
def direction(w):
    return w / np.linalg.norm(w)

'''梯度上升法'''
def gradient_ascent(w, X, eta, n_iter=1e4, epsilon=0.0001):
    '''
    梯度上升法
    :param w:
    :param X:
    :param eta:
    :param n_iter:
    :param epsilon:
    :return:
    '''
    #先化简w为单位向量，方便运算
    w = direction(w)
    i_iter = 0
    while i_iter < n_iter:
        gradient = df_math(w, X)
        last_w = w
        w += gradient * eta
        #每次更新后将w化简为单位向量
        w = direction(w)
        if abs(f(w, X) - f(last_w, X)) < epsilon:
            break
        i_iter += 1
    return w

和线性回归梯度下降寻找极小值点类似，只不过这次的方向是上升。而且有需要注意限定循环次数，避免学习率过大造成陷入无限循环中。测试数据降维的结果

w_init = np.random.random(X_demean.shape[1])
eta = 0.01
w = gradient_ascent(w_init, X_demean, eta)
print(w)
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.plot([0, w[0] * 30], [0, w[1] * 30], color='r')
plt.show()

注意一点，线性回归中，通常将特征系数θ的值设为全部为0的向量，但在主成分分析中w的初始值不能为0！！！
得到映射向量如图中红线

向量w的值为

[0.75748163 0.65285648]

多个主成分分析

刚刚模拟了两个属性中选取一个主成分的过程，实际的降维过程可能会涉及到数据在多个维度的降维，需要依次求解多个主成分。
求解第一个主成分后，假设得到映射的轴为w所表示的向量，如果此时需要求解第二个主成分怎么做。
需要先将数据集在第一个主成分上的分量去掉，然后在没有第一个主成分的基础上再寻找第二个主成分。

由之前的推导所知，蓝色部分的模长是Xi向量和w向量的乘积，又因为w是单位向量，向量的模长乘以方向上的单位向量就可以得到这个向量，去掉Xi在w方向的分量得到新的数据Xi' = Xi - Xipro。
求第二主成分就是在新的数据Xi'上寻找第一主成分。以此类推，之后的主成分求法也是这个套路。

'''均值归零'''
def demean(X):
    # axis=0即按列计算均值，及每个属性的均值，1则是计算行的均值
    return (X - np.mean(X, axis=0))

'''方差函数'''
def f(w, X):
    return np.sum((X.dot(w)**2)) / len(X)

'''方差函数导数'''
def df_ascent(w, X):
    return X.T.dot(X.dot(w)) / 2 * len(X)

'''将向量化简为单位向量'''
def direction(w):
    return w / np.linalg.norm(w)

'''寻找第一主成分'''
def first_component(w_init, X, eta, n_iter=1e4, epsilon=0.0001):
    w = direction(w_init)
    i_iter = 0
    while i_iter < n_iter:
        last_w = w
        gradient = df_ascent(w, X)
        w += eta * gradient
        w = direction(w)
        if abs(f(w, X) - f(last_w, X)) < epsilon:
            break
        i_iter += 1
    return w

'''取前n个主成分'''
def first_n_component(n, X, eta, n_iter=1e4, epsilon=0.0001):
    '''先对数据均值归零'''
    X = demean(X)
    #res记录每个主成分
    res = []
    #进行n次
    for i in range(n):
        #每次初始化w向量，注意不能是0向量
        w = np.random.random(X.shape[1])
        #寻找当前数据的第一主成分并记录到res
        w = first_component(w, X, eta)
        res.append(w)
        #每次减去数据在主成分方向的分量获得新数据
        X = X - X.dot(w).reshape(-1, 1) * w
    return res

代码中的前5个函数和第一主成分分析中的一样，注意最后一个函数。每次求当前数据的第一主成分，并减去数据在第一主成分上的分量获得的新数据继续求第一主成分，知道完成n次为止。
出去数据在第一主成分的分量，如刚才的算式推导的，对于每个样本Xi减去主成分的过程可以写为

X_new[i] = X[i] - (X[i].dot(w)) * w

X[i]是第i个样本，可以看作是1 * n的向量，w是一个1 * n的向量，通过dot()方法可以得到一个常数，即为X[i]在w方向上的模长，模长再乘以w即为在w方向上的分量。
对于含有m个样本的数据集X，可以将整个过程视为

for i in range(len(X)):
    X_new[i] = X[i] - X[i].dot(w) * w

for循环方便理解但更好的方式是直接使用矩阵乘法的方式解决，上面的for循环等价于

X_new = X - X.dot(w).reshape(-1, 1) * w

X是m * n的矩阵，w是1 * n的矩阵，二者dot()之后得到m * 1的矩阵，这里需要使用reshape(-1, 1)使得X和w相乘之后的矩阵真正是一个m * 1的矩阵，之后再乘以w，得到的m * n的矩阵即为每个样本每个维度在w上的分量。再用原始数据X减去即可。
检验下第一次降维后的效果

'''模拟数据'''
X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10, size=100)
'''均值归零'''
X_demean = demean(X)
'''获得第一主成分'''
w_init = np.random.random(X_demean.shape[1])
eta = 0.01
w = first_component(w_init, X_demean, eta)
print(w)
'''画出第一主成分所在向量'''
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.plot([0, w[0] * 30], [0, w[1] * 30], color='r')
plt.show()
'''X减去X在第一主成分上的分量'''
X2 = X_demean - X_demean.dot(w).reshape(-1, 1) * w
plt.scatter(X2[:, 0], X2[:, 1])
plt.show()

第一主成分向量
减去第一主成分分量得到的新数据

下面在新的数据上获取第二主成分

'''获取第二主成分'''
w2_init = np.random.random(X2.shape[1])
w2 = first_component(w2_init, X2, eta)
plt.scatter(X2[:, 0], X2[:, 1])
plt.plot([0, w2[0] * 30], [0, w2[1] * 30], color='r')
plt.show()

得到图像

第二主成分向量

因为是二维的数据，因此前后所得的两个向量应该是垂直的关系，如上文图中所示，两个相互垂直的向量相乘值为0，实验下

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10, size=100)
eta = 0.01
res = first_n_component(2, X, eta)
print(res)
print(res[0].dot(res[1]))

调用之前定义的first_n_component函数，记录前后两次的w，打印res和相乘的结果

[array([0.77854679, 0.62758656]), array([ 0.62758656, -0.77854679])]
-3.885780586188048e-16

从两个向量和他们的乘积来看，乘积近乎为0，可以认为是垂直的。

主成分分析和线性回归的区别

这张图给出的是样本映射到方差最大的轴上的过程。很像线性回归中找一条拟合各个样本点的直线的图像。

这里有几点区别：

主成分分析中，坐标轴表示的是各个特征，而线性回归中，以一个特征的线性回归为例

y轴对应的是训练样本的输出标记。样本是垂直与特征所在的坐标轴，测量的是实际值和预测值的差距，我们的任务是使得这个距离尽可能的小。
主成分分析中坐标轴全部对应特征，映射的方向是垂直于方差最大的轴而非特征所在的坐标轴，且任务是使得样本间距尽可能大。

高维数据降维

主成分分析的作用就是选出能使样本方差最大的维度，选择完维度之后，进入对数据降维的操作。将高维数据映射为低维数据。
假设经过主成分分析之后，左侧X还是数据样本，一个m * n的矩阵，右侧是经过k轮分析之后所获得的k个主成分向量，形成一个k * n的矩阵。

自然想到使用线性代数将两个矩阵相乘，得到m * k的矩阵，将原本的n个特征降为k个特征，即达到了降维映射的目标。因此得到高维数据映射到低维数据的公式：

经过向低维映射得到新的m * k的矩阵Xk如图。低维数据也可以通过与主成分向量w相乘恢复成高维数据。Xk是m * k的矩阵与k * n的矩阵w相乘正好可以得到一个m * n的矩阵。但这个矩阵与原矩阵不一样！！！

低维数据映射回高维数据的公式：

代码实现PCA完整流程

import numpy as np
import matplotlib.pyplot as plt

class PCA:

    def __init__(self, n_component):
        assert n_component >= 1, 'n_component is invalidate'
        self.n_component = n_component
        self.components_ = None

    def __repr__(self):
        return 'PCA(n_component=%d)' % self.n_component

    def fit(self,X, eta, n_iter=1e4, epsilon=0.0001):
        '''
        主成分分析
        :param X: 
        :param eta: 
        :param n_iter: 
        :param epsilon: 
        :return: 
        '''
        assert X.shape[1] >= self.n_component, 'X is invalidate'

        '''均值归零'''
        def demean(X):
            return X - np.mean(X, axis=0)

        '''方差函数'''
        def f(w, X):
            return np.sum(X.dot(w)**2) / len(X)

        '''方差函数导数'''
        def df_ascent(w, X):
            return X.T.dot(X.dot(w)) * 2 / len(X)

        '''将向量化简为单位向量'''
        def direction(w):
            return w / np.linalg.norm(w)

        '''寻找第一主成分'''
        def first_component(w, X, eta, n_iter=1e4, epsilon=0.0001):
            i_iter = 0
            while i_iter < n_iter:
                last_w = w
                gradient = df_ascent(w, X)
                w += eta * gradient
                w = direction(w)
                if abs(f(w, X) - f(last_w, X)) < epsilon:
                    break
                i_iter += 1
            return w

        self.components_ = np.empty(shape=(self.n_component, X.shape[1]))
        X = demean(X)
        for i in range(self.n_component):
            w = np.random.random(X.shape[1])
            w = first_component(w, X, eta, n_iter, epsilon)
            X = X - (X.dot(w)).reshape(-1, 1) * w
            self.components_[i, :] = w
        return self

    def transform(self, X):
        '''
        将X映射到各个主成分中
        :param X:
        :return:
        '''
        assert X.shape[1] == self.components_.shape[1]
        return X.dot(self.components_.T)

    def inverse_transform(self, X):
        '''
        将低维数据转回高维
        :param X:
        :return:
        '''
        assert X.shape[1] == self.components_.shape[0]
        return X.dot(self.components_)

定义一个PCA类，有两个属性n_comnponent记录主成分个数，components_记录各个主成分向量。fit()函数实际就是通过梯度上升法选去n_component个主成分，并通过transform()函数从高维向低维映射，inverse_transform()函数实现了低维映射回高维。
使用测试数据进行测试

'''测试数据'''
X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10, size=100)
eta = 0.01

pca = PCA(1)
pca.fit(X, eta)
X_new = pca.transform(X)
print(pca.components_)
print(X_new.shape)
X_inverse = pca.inverse_transform(X_new)
plt.scatter(X[:, 0], X[:, 1], color='b', alpha=0.5)
plt.scatter(X_inverse[:, 0], X_inverse[:, 1], color='r', alpha=0.5)
plt.show()

得到主成分向量和降维后的数据集的大小

[[0.76233998 0.64717676]]
(100, 1)

同时绘制一下原数据集和低维数据映射回高维数据后的新数据集X_inverse，进行对比

可见，映射回高维的数据和原始数据不一样，说明PCA降维会丢失信息，低维的数据不能恢复成原先高维的数据。

sklearn中的PCA实现

sklearn的PCA在sklearn.decomposition的PCA类中

from sklearn.decomposition import PCA
import numpy as np
import matplotlib.pyplot as plt
import sklearn.datasets as dataset
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
import time

'''模拟数据集'''
X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10, size=100)
pca = PCA(n_components=1)

pca.fit(X)
print(pca.components_)

X_new = pca.transform(X)
print(X_new.shape)

X_inverse = pca.inverse_transform(X_new)
print(X_inverse.shape)
plt.scatter(X[:, 0], X[:, 1], color='b', alpha=0.5)
plt.scatter(X_inverse[:, 0], X_inverse[:, 1], color='r', alpha=0.5)
plt.show()

初始化一个PCA对象，主成分数为1，经过主成分分析后的主成分向量打印

[[0.75763785 0.65267518]]

和之前自定义的PCA类得到的结果近似。也可以通过transform()函数得到降维数据，也可以打印inverse_transform()函数得到映射回高维的数据。这里不给出了。
接下来使用真实数据测试PCA降维对性能的影响

'''真实数据'''
data = dataset.load_digits()
X = data.data
y = data.target
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

knn_clf = KNeighborsClassifier()
start = time.time()
knn_clf.fit(X_train, y_train)
end = time.time()
print("without PCA the KNN classifier's time cost is %s" % (end - start))
print("without PCA the KNN classifier's test score is %s" % knn_clf.score(X_test, y_test))

以KNN分类器为例，使用digits数据集，在不做降维的情况下进行分类的耗时和准确度为

without PCA the KNN classifier's time cost is 0.00571894645690918
without PCA the KNN classifier's test score is 0.9866666666666667

下面测试选2个主成分进行降维的情况

knn_clf = KNeighborsClassifier()
pca = PCA(n_components=2)
pca.fit(X_train)
X_train_new = pca.transform(X_train)
X_test_new = pca.transform(X_test)
start = time.time()
knn_clf.fit(X_train_new, y_train)
end = time.time()
print("with PCA the KNN classifier's time cost is %s" % (end - start))
print("with PCA the KNN classifier's test score is %s" % knn_clf.score(X_test_new, y_test))

得到结果

with PCA the KNN classifier's time cost is 0.0010619163513183594
with PCA the KNN classifier's test score is 0.6066666666666667

可见时间上降维的数据训练明显快了，但准确度降低了。这说明只选择2个主成分效果并不理想。此时可以用到PCA类中的explained_variance_ratio_属性。

print(pca.explained_variance_ratio_)

得到n_component=2时两个主成分对每个样本在主成分方向上的解释度。

[0.14566817 0.13735469]

第一个主成分能解释14%样本的方差，第二个主成分能解释13%样本的方差。即当前主成分量能够维持原数据集方差的百分比，经过2次主成分分析之后，新的数据集一共可以覆盖原始数据28%的方差。解释度越高，预测越准确，解释度之和相加的结果不会超过1，等于1是最理想的状态及可以100%预测准确。相当于n_component=2时，72%的原始信息丢失了，因此预测准确度不高。
如果对原始数据有多少个属性求多少个主成分向量会得到一下结果。

knn_clf = KNeighborsClassifier()
pca = PCA(n_components=X_train.shape[1])
start = time.time()
pca.fit(X_train)
X_train_new = pca.transform(X_train)
X_test_new = pca.transform(X_test)
end = time.time()
knn_clf.fit(X_train_new, y_train)
print("with 64 component PCA the KNN classifier's time cost is %s" % (end - start))
print("with 64 component PCA the KNN classifier's test score is %s" % knn_clf.score(X_test_new, y_test))
print(pca.explained_variance_ratio_)

输出

with 64 component PCA the KNN classifier's time cost is 0.007581233978271484
with 64 component PCA the KNN classifier's test score is 0.9866666666666667
[1.45668166e-01 1.37354688e-01 1.17777287e-01 8.49968861e-02
 5.86018996e-02 5.11542945e-02 4.26605279e-02 3.60119663e-02
 3.41105814e-02 3.05407804e-02 2.42337671e-02 2.28700570e-02
 1.80304649e-02 1.79346003e-02 1.45798298e-02 1.42044841e-02
 1.29961033e-02 1.26617002e-02 1.01728635e-02 9.09314698e-03
 8.85220461e-03 7.73828332e-03 7.60516219e-03 7.11864860e-03
 6.85977267e-03 5.76411920e-03 5.71688020e-03 5.08255707e-03
 4.89020776e-03 4.34888085e-03 3.72917505e-03 3.57755036e-03
 3.26989470e-03 3.14917937e-03 3.09269839e-03 2.87619649e-03
 2.50362666e-03 2.25417403e-03 2.20030857e-03 1.98028746e-03
 1.88195578e-03 1.52769283e-03 1.42823692e-03 1.38003340e-03
 1.17572392e-03 1.07377463e-03 9.55152460e-04 9.00017642e-04
 5.79162563e-04 3.82793717e-04 2.38328586e-04 8.40132221e-05
 5.60545588e-05 5.48538930e-05 1.08077650e-05 4.01354717e-06
 1.23186515e-06 1.05783059e-06 6.06659094e-07 5.86686040e-07
 9.18612290e-34 9.18612290e-34 9.18612290e-34 8.82949950e-34]

训练数据集是一个含有64个特征的数据集，经过主成分分析得到的explained_variance_ratio_也是一个含有64个元素的数据集，表示每个主成分对原始数据方差的解释度，是一个从大到小的排列。最后几个的解释度几乎为零，也就是说对原始方差基本没有作用。但这种方式虽然耗时增加了，但分类的准确度达到98%。在实际情况下，可能会忽略对原始方差影响小的成分，在时间和准确度之间做一个权衡。
例如要求降维达到解释度0.95即可。

knn_clf = KNeighborsClassifier()
pca = PCA(0.95)
start = time.time()
pca.fit(X_train)
X_train_new = pca.transform(X_train)
X_test_new = pca.transform(X_test)
end = time.time()
knn_clf.fit(X_train_new, y_train)
print("with 0.95 PCA the KNN classifier's time cost is %s" % (end - start))
print("with 0.95 component PCA the KNN classifier's test score is %s" % knn_clf.score(X_test_new, y_test))
print(pca.explained_variance_ratio_.shape)

得到

with 0.95 PCA the KNN classifier's time cost is 0.007498979568481445
with 0.95 component PCA the KNN classifier's test score is 0.98
(28,)

可见，当解释度0.95时，耗时缩短了一些，分类准确度还是98%，此时的主成分向量数为28，省去了其他36个特征维度。
n_component低可能造成信息丢失严重，但并非没有用。例如还是n_component=2进行降维。

pca = PCA(n_components=2)
pca.fit(X_train)
X_train_new = pca.transform(X_train)
for i in range(10):
    plt.scatter(X_train_new[y_train == i, 0], X_train_new[y_train == i, 1], alpha=0.8)
plt.show()

对降维后的数据进行分类别绘制得到

在二维空间中，类别见的区分度已经相对明显，如果仅仅是对蓝色和红色两种类别进行分析，则二维可能已经足够了。

PCA的其他用途

用途一：数据降噪

文章开头有提到PCA可以降噪，噪声在生产过程中可能因为各种原因出现，影响数据的准确性，PCA通过选取主成分将原有数据映射到低维数据再映射回高维数据的方式进行一定程度的降噪，以digits数据为例。

import numpy as np
import matplotlib.pyplot as plt
import sklearn.datasets as dataset
from sklearn.decomposition import PCA

'''加载数据'''
data = dataset.load_digits()
X = data.data
y = data.target
#有噪声的数据
noisy_digits = X + np.random.normal(0, 4, size=X.shape)
#取100个样例数据
example_digits = noisy_digits[y == 0, :][:10]

for i in range(1, 10):
    X_num = noisy_digits[y == i, :][:10]
    example_digits = np.vstack([example_digits, X_num])

'''绘图函数'''
def plot_digits(data):
    fig, axes = plt.subplots(10, 10, figsize=(10, 10),
                            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(example_digits)
plt.show()

在digits数据特征集中人工加入噪声，每类样本选去10个进行打印，打印的效果如

此时通过PCA进行主成分分析，获得映射后的低维数据，在映射回高维数据并以相同的方式打印

'''PCA数据降噪'''
#这里选用n_component=0.5进行主成分分析
pca = PCA(0.5)
pca.fit(example_digits)
#映射到低维数据example_digits_new
example_digits_new = pca.transform(example_digits)
#再从低维数据映射回高维数据并绘制
example_digits_no_noise = pca.inverse_transform(example_digits_new)
plot_digits(example_digits_no_noise)
plt.show()

得到的图像

明显的好于原始数据图像，降噪有一定效果。

用途二：人脸识别与特征脸

假设原始的人像数据是一个m * n的矩阵，在经过主成分分析之后，一个k * n的主成分矩阵。如果将主成分矩阵也看成是由k个样本组成的矩阵，那么可以理解为第一个样本是最重要的样本，第二个次之，依次往后。原始数据矩阵可以视为有m个人脸样本的集合，如果将主成分矩阵每一行也看做是一个样本的话，每一行相当于是一个由原始人脸数据矩阵经过主成分分析得到的特征脸矩阵，这个矩阵含有k个特征脸。每个特征脸表达了原有样本中人脸的部分特征。
使用sklearn的fetch_lfw_people数据集，选去36张人脸进行测试。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_lfw_people
from sklearn.decomposition import PCA
import ssl

ssl._create_default_https_context = ssl._create_unverified_context
'''加载数据'''
face = fetch_lfw_people()

#face数据集的特征集是13233 * 2914的矩阵，打乱顺序，取其中36个进行测试
random_index = np.random.permutation(len(face.data))
X = face.data[random_index]
example_face = X[:36, :]

#绘制原始的人脸
'''绘图函数'''
def plot_faces(data):
    fig, axes = plt.subplots(10, 10, figsize=(10, 10),
                            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(62, 47), cmap='binary', interpolation='nearest', clim=(0, 16))

plot_faces(example_face)
plt.show()

#数据集大，使用随机的方式求解pca，提升效率
pca = PCA(svd_solver='randomized')
pca.fit(example_face)
agent_faces = pca.components_[:36, :]
plot_faces(agent_faces)
plt.show()

选取36张人脸绘制原始图像，再根据这些数据进行PCA操作，再绘制特征脸。

特征脸