主成分分析一個(gè)非監(jiān)督學(xué)習(xí)算法，主要用于數(shù)據(jù)降維霍狰，通過(guò)降維可以發(fā)現(xiàn)數(shù)據(jù)更容易理解的特征，其他作用也有可視化缓待、降噪等蚓耽。
假設(shè)現(xiàn)有樣本的分布如圖。

樣本有兩個(gè)特征旋炒，如果對(duì)樣本進(jìn)行降維，首先可以考慮基于坐標(biāo)軸進(jìn)行降維签杈。會(huì)有如下兩種方式瘫镇。

在數(shù)學(xué)上泽铛，樣本間距離反映了樣本的離散程度，而方差可以很好的表示這種程度辑鲤，因此我們定義樣本間距離——方差。

主成分分析PCA

PCA算法是最經(jīng)典的數(shù)據(jù)降維算法，步驟大致如下：

①將樣本的均值歸零（demean）
②求一個(gè)軸的方向w寸齐，使得映射后的方差值最大

均值歸零

原有樣本分布

歸零之后

映射軸和方差

映射的軸可以視為一個(gè)向量盐杂，w = (w1, w2, ..., wn)逗载，n即為樣本的特征數(shù)。
方差公式

我們的目標(biāo)就是使得這個(gè)公式得到的方差值最大链烈。
如果將樣本的數(shù)據(jù)映射到坐標(biāo)系中厉斟，那么每個(gè)樣本就可以看做是個(gè)向量，而樣本的均值也可以看做是一個(gè)向量强衡。這樣我們的方差公式可以變?yōu)?/p>

Xi是樣本點(diǎn)触幼，有Xi,1和Xi,2兩個(gè)特征，假設(shè)我們求得的最大樣本間距的映射軸是w表示的向量究飞，設(shè)w=(w1, w2)置谦，那么Xi映射到w上的點(diǎn)如圖中藍(lán)色線表示。
高中數(shù)學(xué)中亿傅，兩個(gè)向量相乘等于兩個(gè)向量的模乘以兩個(gè)向量的夾角余弦值媒峡。
通常在運(yùn)算中，我們會(huì)將向量w設(shè)為單位向量（即模長(zhǎng)為1的向量）葵擎，因此等式右側(cè)的w的模長(zhǎng)可以省略谅阿。
又因?yàn)閄i的模長(zhǎng)乘以?shī)A角余弦值正好是藍(lán)色部分的長(zhǎng)度，因此我們得到最后一個(gè)等式酬滤。
因此签餐，我們對(duì)之前求方差的算式進(jìn)行修改得到

梯度上升法

回憶梯度下降法乔询，通過(guò)求函數(shù)在一個(gè)點(diǎn)的導(dǎo)數(shù)值，得到函數(shù)在該點(diǎn)的斜率韵洋，之后使用減去斜率與學(xué)習(xí)率乘積的方式下降竿刁，逐漸靠近極小值點(diǎn)黄锤。
梯度上升法與之類似，上升不過(guò)就是將減法變?yōu)榧臃ㄊ嘲荨５谝徊竭€是對(duì)方差函數(shù)f在每個(gè)特征維度求導(dǎo)鸵熟。

對(duì)括號(hào)中的式子化簡(jiǎn)

主成分分析

單個(gè)主成分分析

先模擬一個(gè)數(shù)據(jù)集

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

該數(shù)據(jù)集有兩個(gè)屬性，打印初始圖像

接下來(lái)對(duì)數(shù)據(jù)進(jìn)行主成分分析鱼冀，第一步是均值歸零报破，定義相應(yīng)的函數(shù)

'''均值歸零'''
def demean(X):
    # axis=0即按列計(jì)算均值，及每個(gè)屬性的均值千绪，1則是計(jì)算行的均值
    return (X - np.mean(X, axis=0))

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

計(jì)算均值調(diào)用的是numpy自帶的函數(shù)，均值歸零之后打印歸零后的圖像

分布沒(méi)有變化梗脾，但從坐標(biāo)軸可以看出均值歸零后的效果荸型。
接下來(lái)就是對(duì)方差函數(shù)和其導(dǎo)數(shù)函數(shù)的定義

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

'''方差函數(shù)導(dǎo)數(shù)'''
def df_math(w, X):
    return X.T.dot(X.dot(w)) * 2. / len(X)

'''將向量化簡(jiǎn)為單位向量'''
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:
    '''
    #先化簡(jiǎn)w為單位向量，方便運(yùn)算
    w = direction(w)
    i_iter = 0
    while i_iter < n_iter:
        gradient = df_math(w, X)
        last_w = w
        w += gradient * eta
        #每次更新后將w化簡(jiǎn)為單位向量
        w = direction(w)
        if abs(f(w, X) - f(last_w, X)) < epsilon:
            break
        i_iter += 1
    return w

和線性回歸梯度下降尋找極小值點(diǎn)類似炸茧，只不過(guò)這次的方向是上升瑞妇。而且有需要注意限定循環(huán)次數(shù)，避免學(xué)習(xí)率過(guò)大造成陷入無(wú)限循環(huán)中梭冠。測(cè)試數(shù)據(jù)降維的結(jié)果

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

注意一點(diǎn)辕狰，線性回歸中，通常將特征系數(shù)θ的值設(shè)為全部為0的向量控漠，但在主成分分析中w的初始值不能為0Ｂ丁！盐捷！
得到映射向量如圖中紅線

向量w的值為

[0.75748163 0.65285648]

多個(gè)主成分分析

剛剛模擬了兩個(gè)屬性中選取一個(gè)主成分的過(guò)程偶翅，實(shí)際的降維過(guò)程可能會(huì)涉及到數(shù)據(jù)在多個(gè)維度的降維，需要依次求解多個(gè)主成分碉渡。
求解第一個(gè)主成分后聚谁，假設(shè)得到映射的軸為w所表示的向量，如果此時(shí)需要求解第二個(gè)主成分怎么做滞诺。
需要先將數(shù)據(jù)集在第一個(gè)主成分上的分量去掉形导，然后在沒(méi)有第一個(gè)主成分的基礎(chǔ)上再尋找第二個(gè)主成分环疼。

由之前的推導(dǎo)所知，藍(lán)色部分的模長(zhǎng)是Xi向量和w向量的乘積朵耕，又因?yàn)閣是單位向量炫隶，向量的模長(zhǎng)乘以方向上的單位向量就可以得到這個(gè)向量，去掉Xi在w方向的分量得到新的數(shù)據(jù)Xi' = Xi - Xipro憔披。
求第二主成分就是在新的數(shù)據(jù)Xi'上尋找第一主成分等限。以此類推，之后的主成分求法也是這個(gè)套路芬膝。

'''均值歸零'''
def demean(X):
    # axis=0即按列計(jì)算均值望门，及每個(gè)屬性的均值，1則是計(jì)算行的均值
    return (X - np.mean(X, axis=0))

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

'''方差函數(shù)導(dǎo)數(shù)'''
def df_ascent(w, X):
    return X.T.dot(X.dot(w)) / 2 * len(X)

'''將向量化簡(jiǎn)為單位向量'''
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個(gè)主成分'''
def first_n_component(n, X, eta, n_iter=1e4, epsilon=0.0001):
    '''先對(duì)數(shù)據(jù)均值歸零'''
    X = demean(X)
    #res記錄每個(gè)主成分
    res = []
    #進(jìn)行n次
    for i in range(n):
        #每次初始化w向量锰霜，注意不能是0向量
        w = np.random.random(X.shape[1])
        #尋找當(dāng)前數(shù)據(jù)的第一主成分并記錄到res
        w = first_component(w, X, eta)
        res.append(w)
        #每次減去數(shù)據(jù)在主成分方向的分量獲得新數(shù)據(jù)
        X = X - X.dot(w).reshape(-1, 1) * w
    return res

代碼中的前5個(gè)函數(shù)和第一主成分分析中的一樣筹误，注意最后一個(gè)函數(shù)。每次求當(dāng)前數(shù)據(jù)的第一主成分癣缅，并減去數(shù)據(jù)在第一主成分上的分量獲得的新數(shù)據(jù)繼續(xù)求第一主成分厨剪，知道完成n次為止。
出去數(shù)據(jù)在第一主成分的分量友存，如剛才的算式推導(dǎo)的祷膳，對(duì)于每個(gè)樣本Xi減去主成分的過(guò)程可以寫為

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

X[i]是第i個(gè)樣本，可以看作是1 * n的向量屡立，w是一個(gè)1 * n的向量直晨，通過(guò)dot()方法可以得到一個(gè)常數(shù)，即為X[i]在w方向上的模長(zhǎng)膨俐，模長(zhǎng)再乘以w即為在w方向上的分量勇皇。
對(duì)于含有m個(gè)樣本的數(shù)據(jù)集X，可以將整個(gè)過(guò)程視為

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

for循環(huán)方便理解但更好的方式是直接使用矩陣乘法的方式解決焚刺，上面的for循環(huán)等價(jià)于

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

X是m * n的矩陣敛摘，w是1 * n的矩陣，二者dot()之后得到m * 1的矩陣乳愉，這里需要使用reshape(-1, 1)使得X和w相乘之后的矩陣真正是一個(gè)m * 1的矩陣兄淫，之后再乘以w，得到的m * n的矩陣即為每個(gè)樣本每個(gè)維度在w上的分量匾委。再用原始數(shù)據(jù)X減去即可拖叙。
檢驗(yàn)下第一次降維后的效果

'''模擬數(shù)據(jù)'''
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()

第一主成分向量

減去第一主成分分量得到的新數(shù)據(jù)

下面在新的數(shù)據(jù)上獲取第二主成分

'''獲取第二主成分'''
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()

得到圖像

第二主成分向量

因?yàn)槭嵌S的數(shù)據(jù)，因此前后所得的兩個(gè)向量應(yīng)該是垂直的關(guān)系赂乐，如上文圖中所示薯鳍，兩個(gè)相互垂直的向量相乘值為0，實(shí)驗(yàn)下

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

調(diào)用之前定義的first_n_component函數(shù)，記錄前后兩次的w挖滤，打印res和相乘的結(jié)果

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

從兩個(gè)向量和他們的乘積來(lái)看崩溪，乘積近乎為0，可以認(rèn)為是垂直的斩松。

主成分分析和線性回歸的區(qū)別

這張圖給出的是樣本映射到方差最大的軸上的過(guò)程伶唯。很像線性回歸中找一條擬合各個(gè)樣本點(diǎn)的直線的圖像。

這里有幾點(diǎn)區(qū)別：

主成分分析中惧盹，坐標(biāo)軸表示的是各個(gè)特征乳幸，而線性回歸中，以一個(gè)特征的線性回歸為例

y軸對(duì)應(yīng)的是訓(xùn)練樣本的輸出標(biāo)記钧椰。樣本是垂直與特征所在的坐標(biāo)軸粹断，測(cè)量的是實(shí)際值和預(yù)測(cè)值的差距，我們的任務(wù)是使得這個(gè)距離盡可能的小嫡霞。
主成分分析中坐標(biāo)軸全部對(duì)應(yīng)特征瓶埋，映射的方向是垂直于方差最大的軸而非特征所在的坐標(biāo)軸，且任務(wù)是使得樣本間距盡可能大诊沪。

高維數(shù)據(jù)降維

主成分分析的作用就是選出能使樣本方差最大的維度养筒，選擇完維度之后，進(jìn)入對(duì)數(shù)據(jù)降維的操作端姚。將高維數(shù)據(jù)映射為低維數(shù)據(jù)晕粪。
假設(shè)經(jīng)過(guò)主成分分析之后，左側(cè)X還是數(shù)據(jù)樣本渐裸，一個(gè)m * n的矩陣兵多，右側(cè)是經(jīng)過(guò)k輪分析之后所獲得的k個(gè)主成分向量，形成一個(gè)k * n的矩陣橄仆。

自然想到使用線性代數(shù)將兩個(gè)矩陣相乘，得到m * k的矩陣衅斩，將原本的n個(gè)特征降為k個(gè)特征盆顾，即達(dá)到了降維映射的目標(biāo)。因此得到高維數(shù)據(jù)映射到低維數(shù)據(jù)的公式：

代碼實(shí)現(xià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)

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

        '''方差函數(shù)導(dǎo)數(shù)'''
        def df_ascent(w, X):
            return X.T.dot(X.dot(w)) * 2 / len(X)

        '''將向量化簡(jiǎn)為單位向量'''
        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映射到各個(gè)主成分中
        :param X:
        :return:
        '''
        assert X.shape[1] == self.components_.shape[1]
        return X.dot(self.components_.T)

    def inverse_transform(self, X):
        '''
        將低維數(shù)據(jù)轉(zhuǎn)回高維
        :param X:
        :return:
        '''
        assert X.shape[1] == self.components_.shape[0]
        return X.dot(self.components_)

定義一個(gè)PCA類溜畅，有兩個(gè)屬性n_comnponent記錄主成分個(gè)數(shù)捏卓，components_記錄各個(gè)主成分向量。fit()函數(shù)實(shí)際就是通過(guò)梯度上升法選去n_component個(gè)主成分慈格，并通過(guò)transform()函數(shù)從高維向低維映射怠晴，inverse_transform()函數(shù)實(shí)現(xiàn)了低維映射回高維遥金。
使用測(cè)試數(shù)據(jù)進(jìn)行測(cè)試

'''測(cè)試數(shù)據(jù)'''
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()

得到主成分向量和降維后的數(shù)據(jù)集的大小

[[0.76233998 0.64717676]]
(100, 1)

同時(shí)繪制一下原數(shù)據(jù)集和低維數(shù)據(jù)映射回高維數(shù)據(jù)后的新數(shù)據(jù)集X_inverse，進(jìn)行對(duì)比

可見(jiàn)蒜田，映射回高維的數(shù)據(jù)和原始數(shù)據(jù)不一樣稿械，說(shuō)明PCA降維會(huì)丟失信息，低維的數(shù)據(jù)不能恢復(fù)成原先高維的數(shù)據(jù)冲粤。

sklearn中的PCA實(shí)現(xiàn)

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

'''模擬數(shù)據(jù)集'''
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()

初始化一個(gè)PCA對(duì)象美莫，主成分?jǐn)?shù)為1，經(jīng)過(guò)主成分分析后的主成分向量打印

[[0.75763785 0.65267518]]

和之前自定義的PCA類得到的結(jié)果近似梯捕。也可以通過(guò)transform()函數(shù)得到降維數(shù)據(jù)厢呵，也可以打印inverse_transform()函數(shù)得到映射回高維的數(shù)據(jù)。這里不給出了科阎。
接下來(lái)使用真實(shí)數(shù)據(jù)測(cè)試PCA降維對(duì)性能的影響

'''真實(shí)數(shù)據(jù)'''
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數(shù)據(jù)集，在不做降維的情況下進(jìn)行分類的耗時(shí)和準(zhǔn)確度為

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

下面測(cè)試選2個(gè)主成分進(jìn)行降維的情況

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

得到結(jié)果

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

可見(jiàn)時(shí)間上降維的數(shù)據(jù)訓(xùn)練明顯快了锣笨，但準(zhǔn)確度降低了蝌矛。這說(shuō)明只選擇2個(gè)主成分效果并不理想。此時(shí)可以用到PCA類中的explained_variance_ratio_屬性错英。

print(pca.explained_variance_ratio_)

得到n_component=2時(shí)兩個(gè)主成分對(duì)每個(gè)樣本在主成分方向上的解釋度入撒。

[0.14566817 0.13735469]

第一個(gè)主成分能解釋14%樣本的方差，第二個(gè)主成分能解釋13%樣本的方差椭岩。即當(dāng)前主成分量能夠維持原數(shù)據(jù)集方差的百分比茅逮，經(jīng)過(guò)2次主成分分析之后，新的數(shù)據(jù)集一共可以覆蓋原始數(shù)據(jù)28%的方差判哥。解釋度越高献雅，預(yù)測(cè)越準(zhǔn)確，解釋度之和相加的結(jié)果不會(huì)超過(guò)1塌计，等于1是最理想的狀態(tài)及可以100%預(yù)測(cè)準(zhǔn)確挺身。相當(dāng)于n_component=2時(shí)，72%的原始信息丟失了锌仅，因此預(yù)測(cè)準(zhǔn)確度不高章钾。
如果對(duì)原始數(shù)據(jù)有多少個(gè)屬性求多少個(gè)主成分向量會(huì)得到一下結(jié)果。

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]

訓(xùn)練數(shù)據(jù)集是一個(gè)含有64個(gè)特征的數(shù)據(jù)集热芹，經(jīng)過(guò)主成分分析得到的explained_variance_ratio_也是一個(gè)含有64個(gè)元素的數(shù)據(jù)集贱傀，表示每個(gè)主成分對(duì)原始數(shù)據(jù)方差的解釋度，是一個(gè)從大到小的排列伊脓。最后幾個(gè)的解釋度幾乎為零府寒，也就是說(shuō)對(duì)原始方差基本沒(méi)有作用。但這種方式雖然耗時(shí)增加了，但分類的準(zhǔn)確度達(dá)到98%椰棘。在實(shí)際情況下纺棺，可能會(huì)忽略對(duì)原始方差影響小的成分，在時(shí)間和準(zhǔn)確度之間做一個(gè)權(quán)衡邪狞。
例如要求降維達(dá)到解釋度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,)

可見(jiàn)，當(dāng)解釋度0.95時(shí)帆卓，耗時(shí)縮短了一些巨朦，分類準(zhǔn)確度還是98%，此時(shí)的主成分向量數(shù)為28剑令，省去了其他36個(gè)特征維度糊啡。
n_component低可能造成信息丟失嚴(yán)重，但并非沒(méi)有用吁津。例如還是n_component=2進(jìn)行降維棚蓄。

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

對(duì)降維后的數(shù)據(jù)進(jìn)行分類別繪制得到

在二維空間中，類別見(jiàn)的區(qū)分度已經(jīng)相對(duì)明顯碍脏，如果僅僅是對(duì)藍(lán)色和紅色兩種類別進(jìn)行分析梭依，則二維可能已經(jīng)足夠了。

PCA的其他用途

用途一：數(shù)據(jù)降噪

文章開(kāi)頭有提到PCA可以降噪典尾，噪聲在生產(chǎn)過(guò)程中可能因?yàn)楦鞣N原因出現(xiàn)役拴，影響數(shù)據(jù)的準(zhǔn)確性，PCA通過(guò)選取主成分將原有數(shù)據(jù)映射到低維數(shù)據(jù)再映射回高維數(shù)據(jù)的方式進(jìn)行一定程度的降噪钾埂，以digits數(shù)據(jù)為例河闰。

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

'''加載數(shù)據(jù)'''
data = dataset.load_digits()
X = data.data
y = data.target
#有噪聲的數(shù)據(jù)
noisy_digits = X + np.random.normal(0, 4, size=X.shape)
#取100個(gè)樣例數(shù)據(jù)
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])

'''繪圖函數(shù)'''
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))

'''繪制帶有噪聲的數(shù)據(jù)'''
plot_digits(example_digits)
plt.show()

在digits數(shù)據(jù)特征集中人工加入噪聲，每類樣本選去10個(gè)進(jìn)行打印褥紫，打印的效果如

此時(shí)通過(guò)PCA進(jìn)行主成分分析姜性，獲得映射后的低維數(shù)據(jù)，在映射回高維數(shù)據(jù)并以相同的方式打印

'''PCA數(shù)據(jù)降噪'''
#這里選用n_component=0.5進(jìn)行主成分分析
pca = PCA(0.5)
pca.fit(example_digits)
#映射到低維數(shù)據(jù)example_digits_new
example_digits_new = pca.transform(example_digits)
#再?gòu)牡途S數(shù)據(jù)映射回高維數(shù)據(jù)并繪制
example_digits_no_noise = pca.inverse_transform(example_digits_new)
plot_digits(example_digits_no_noise)
plt.show()

得到的圖像

明顯的好于原始數(shù)據(jù)圖像髓考，降噪有一定效果污抬。

用途二：人臉識(shí)別與特征臉

假設(shè)原始的人像數(shù)據(jù)是一個(gè)m * n的矩陣，在經(jīng)過(guò)主成分分析之后绳军，一個(gè)k * n的主成分矩陣。如果將主成分矩陣也看成是由k個(gè)樣本組成的矩陣矢腻，那么可以理解為第一個(gè)樣本是最重要的樣本门驾，第二個(gè)次之，依次往后多柑。原始數(shù)據(jù)矩陣可以視為有m個(gè)人臉樣本的集合奶是，如果將主成分矩陣每一行也看做是一個(gè)樣本的話，每一行相當(dāng)于是一個(gè)由原始人臉數(shù)據(jù)矩陣經(jīng)過(guò)主成分分析得到的特征臉矩陣，這個(gè)矩陣含有k個(gè)特征臉聂沙。每個(gè)特征臉表達(dá)了原有樣本中人臉的部分特征秆麸。
使用sklearn的fetch_lfw_people數(shù)據(jù)集，選去36張人臉進(jìn)行測(cè)試及汉。

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
'''加載數(shù)據(jù)'''
face = fetch_lfw_people()

#face數(shù)據(jù)集的特征集是13233 * 2914的矩陣沮趣，打亂順序，取其中36個(gè)進(jìn)行測(cè)試
random_index = np.random.permutation(len(face.data))
X = face.data[random_index]
example_face = X[:36, :]

#繪制原始的人臉
'''繪圖函數(shù)'''
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()

#數(shù)據(jù)集大坷随，使用隨機(jī)的方式求解pca房铭，提升效率
pca = PCA(svd_solver='randomized')
pca.fit(example_face)
agent_faces = pca.components_[:36, :]
plot_faces(agent_faces)
plt.show()

選取36張人臉繪制原始圖像，再根據(jù)這些數(shù)據(jù)進(jìn)行PCA操作温眉，再繪制特征臉缸匪。

特征臉