sklearn與機器學習系列專題之降維(二)一文弄懂LDA特徵篩選&降維

在上一篇推文中，我們詳解了PCA算法。這是機器學習中最為常用的降維方法，能降低算法的計算開銷，使得數據集更容易處理，且完全無參數的限制。但是，如果用戶對觀測對象有一定的先驗知識，掌握了數據的一些特徵，卻很難按照預想的方法對處理過程進行幹預，可能達不到預期的效果，在非高斯分布的情況下，PCA方法得出的主元可能也並不是最優的。

這時候，就要線性判別分析降維（Linear discriminant Analysis，LDA）隆重登場啦！線性判別分析(也有叫做Fisher Linear Discriminant)，是一種有監督的線性降維算法，與PCA算法尋找數據集中方差最大的方向作為主成分分量的軸不同，LDA是為了使得降維後的數據點儘可能地容易被區分。

線性判別分析的中心思想是，最大化類間距離，最小化類內距離。

如下圖，紅色的點代表class1類別的數據，藍色代表class2的數據，根據PCA算法，數據應該映射到方差最大的方向，即Y軸，但是class1和class2兩個不同類別的數據就會完全的混合在一起，很難區分開。所以使用PCA算法進行降維後再進行分類的效果會非常差，這時候就需要我們使用LDA算法，將數據映射到X軸上。

由於理論推導可以通過網絡、書籍等渠道了解到很多，此處以二分類的LDA為例，簡要介紹原理。

上一節說道，LDA投影后，希望同一類別數據的投影點儘可能地靠近，而不同類數據的類別中心之間儘可能地大，因此，我們在推導過程中，將這一中心思想量化，如下，分子為不同類別之間中心點距離的差的平方，分母為同類樣本之間的離散程度。最終目標是將J(W)最大化。

先將每一類中的散列值展開：

其中，x是投影前的特徵，w為投影矩陣。

再對分子展開：

代入原式，目標函數最終可簡化為：

其中SB代表類間散布矩陣，Sw代表類內散布矩陣。利用拉格朗日乘子法對上式變換，可得：

觀察一下式子，可將它視為線性代數中的特徵向量。因此，在線性判別分析中，只需要得到類內和類間散布矩陣，再求解其特徵向量，即可得到投影方向，再對數據執行相應的矩陣變化，即可完成全部的降維過程。

LDA在sklearn中有自帶的包，導入方法為

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
使用自帶庫時，將上一篇推文代碼中的PCA換成LDA即可，此處不再贅述代碼過程。為了更清楚地向讀者展示LDA算法的原理，以下將以「萬能」的鳶尾花數據集為例，手撕LDA算法（選自《跟著迪哥學python》）！
# 自己來定義列名feature_dict = {i:label for i,label in zip(                range(4),                  ('sepal length in cm','sepal width in cm','petal length in cm','petal width in cm', ))}
label_dict = {i:label for i,label in zip(                range(1,4),                  ('Setosa','Versicolor','Virginica'                 ))}import pandas as pd# 數據讀取，大家也可以先下載下來直接讀取df = pd.io.parsers.read_csv(    filepath_or_buffer='https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data',    header=None,    sep=',',    )# 指定列名df.columns = [l for i,l in sorted(feature_dict.items())] + ['class label']
df.head()
數據集共有150條數據，每條數據有4個特徵，現在需要將四維特徵降成二維。觀察輸出結果可以發現，其特徵已經是數值數據，不需要做額外處理，但是需要轉換一下標籤：
from sklearn.preprocessing import LabelEncoder
X = df[['sepal length in cm','sepal width in cm','petal length in cm','petal width in cm']].valuesy = df['class label'].values
# 製作標籤{1: 'Setosa', 2: 'Versicolor', 3:'Virginica'}enc = LabelEncoder()label_encoder = enc.fit(y)y = label_encoder.transform(y) + 1
上述代碼使用了sklearn工具包中的「LabelEncoder」用於快速完成標籤轉換，可以發現基本上所有sklearn中的數據處理操作都是分兩步走，先fit再transform。
在計算過程中需要基於均值來判斷距離，因此先要對數據中各個特徵求均值，但是只求4個特徵的均值能滿足要求嗎？不要忘記任務中還有3種花，相當於3個類別，所以也要對每種花分別求其各個特徵的均值。
import numpy as np#設置小數點的位數np.set_printoptions(precision=4)#這裡會保存所有的均值mean_vectors = []# 要計算3個類別for cl in range(1,4):# 求當前類別各個特徵均值    mean_vectors.append(np.mean(X[y==cl], axis=0))    print('均值類別 %s: %s\n' %(cl, mean_vectors[cl-1]))
輸出如下：
接下來計算類內散布矩陣
# 原始數據中有4個特徵S_W = np.zeros((4,4))# 要考慮不同類別，自己算自己的for cl,mv in zip(range(1,4), mean_vectors):    class_sc_mat = np.zeros((4,4))# 選中屬於當前類別的數據for row in X[y == cl]:# 這裡相當於對各個特徵分別進行計算，用矩陣的形式        row, mv = row.reshape(4,1), mv.reshape(4,1)# 跟公式一樣        class_sc_mat += (row-mv).dot((row-mv).T)    S_W += class_sc_matprint('類內散布矩陣:\n', S_W)
輸出為：
繼續計算類間散布矩陣：
# 全局均值overall_mean = np.mean(X, axis=0)# 構建類間散布矩陣S_B = np.zeros((4,4))# 對各個類別進行計算for i,mean_vec in enumerate(mean_vectors):    #當前類別的樣本數n = X[y==i+1,:].shape[0]mean_vec = mean_vec.reshape(4,1)overall_mean = overall_mean.reshape(4,1)    # 如上述公式進行計算S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)
print('類間散布矩陣:\n', S_B)
輸出為：
得到類內和類間散布矩陣後，還需將它們組合在一起，然後求解矩陣的特徵向量：
#求解矩陣特徵值，特徵向量eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))# 拿到每一個特徵值和其所對應的特徵向量for i in range(len(eig_vals)):    eigvec_sc = eig_vecs[:,i].reshape(4,1)    print('\n特徵向量 {}: \n{}'.format(i+1, eigvec_sc.real))    print('特徵值 {:}: {:.2e}'.format(i+1, eig_vals[i].real))
輸出特徵向量：
輸出結果得到4個特徵值和其所對應的特徵向量。特徵向量直接觀察起來比較麻煩，因為投影方向在高維上很難理解；特徵值還是比較直觀的，這裡可以認為特徵值代表的是其所對應特徵向量的重要程度，也就是特徵值越大，其所對應的特徵向量就越重要，所以接下來可以對特徵值按大小進行排序，排在前面的越重要，排在後面的就沒那麼重要了。
#特徵值和特徵向量配對eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
# 按特徵值大小進行排序eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
print('特徵值排序結果:\n')for i in eig_pairs:    print(i[0])
特徵值排序結果為：
32.27195779972981
0.27756686384003953
1.1483362279322388e-14
3.422458920849769e-15
print('特徵值佔總體百分比:\n')eigv_sum = sum(eig_vals)for i,j in enumerate(eig_pairs):print('特徵值 {0:}: {1:.2%}'.format(i+1, (j[0]/eigv_sum).real))
特徵值佔總體百分比：
特徵值 1: 99.15%
特徵值 2: 0.85%
特徵值 3: 0.00%
特徵值 4: 0.00%
可以看出，列印出來的結果差異很大，第一個特徵值佔據總體的99.15%，第二個特徵值只佔0.85%，第三和第四個特徵值，看起來微不足道。這表示對鳶尾花數據進行降維時，可以把特徵數據降到二維甚至一維，但沒必要降到三維。
既然已經有結論，選擇把數據降到二維，只需選擇特徵值1、特徵值2所對應的特徵向量即可：
W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1)))print('矩陣W:\n', W.real)
矩陣W:
[[-0.2049 -0.009 ]
 [-0.3871 -0.589 ]
 [ 0.5465  0.2543]
 [ 0.7138 -0.767 ]]
這也是最終所需的投影方向，只需和原始數據組合，就可以得到降維結果：
# 執行降維操作X_lda = X.dot(W)X_lda.shape
現在可以看到數據維度從原始的（150,4）降到（150,2），到此就完成全部的降維工作。接下來對比分析一下降維後結果，為了方便可視化展示，在原始四維數據集中隨機選擇兩維進行繪圖展示。
from matplotlib import pyplot as plt# 可視化展示def plot_step_lda():
    ax = plt.subplot(111)    for label,marker,color in zip(        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):
        plt.scatter(x=X[:,0].real[y == label],                y=X[:,1].real[y == label],                marker=marker,                color=color,                alpha=0.5,                label=label_dict[label]                )
    plt.xlabel('X[0]')    plt.ylabel('X[1]')
    leg = plt.legend(loc='upper right', fancybox=True)    leg.get_frame().set_alpha(0.5)    plt.title('Original data')
    # 把邊邊角角隱藏起來    plt.tick_params(axis="both", which="both", bottom="off", top="off",            labelbottom="on", left="off", right="off", labelleft="on")
    # 為了看的清晰些，儘量簡潔    ax.spines["top"].set_visible(False)    ax.spines["right"].set_visible(False)    ax.spines["bottom"].set_visible(False)    ax.spines["left"].set_visible(False)
    plt.grid()    plt.tight_layout    plt.show()
plot_step_lda()
從上述輸出結果可以發現，如果對原始數據集隨機取兩維數據，數據集並不能按類別劃分開，很多數據都堆疊在一起（尤其是圖中方塊和圓形數據點）。再來看看降維後的數據點分布，繪圖代碼保持不變，只需要傳入降維後的兩維數據即可。
from matplotlib import pyplot as plt# 可視化展示def plot_step_lda():
    ax = plt.subplot(111)    for label,marker,color in zip(        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):
        plt.scatter(x=X_lda[:,0].real[y == label],                y=X_lda[:,1].real[y == label],                marker=marker,                color=color,                alpha=0.5,                label=label_dict[label]                )
    plt.xlabel('LD1')    plt.ylabel('LD2')
    leg = plt.legend(loc='upper right', fancybox=True)    leg.get_frame().set_alpha(0.5)    plt.title('LDA on iris')
    # 把邊邊角角隱藏起來    plt.tick_params(axis="both", which="both", bottom="off", top="off",            labelbottom="on", left="off", right="off", labelleft="on")
    # 為了看的清晰些，儘量簡潔    ax.spines["top"].set_visible(False)    ax.spines["right"].set_visible(False)    ax.spines["bottom"].set_visible(False)    ax.spines["left"].set_visible(False)
    plt.grid()    plt.tight_layout    plt.show()
plot_step_lda()
可以明顯看到，坐標軸變成LD1與LD2，這就是降維後的結果，從數據點的分布來看，混雜在一起的數據不多，劃分起來就更容易。這就是經過一步步計算得到的最終降維結果。
在實驗過程中，我們也可以發現，當拿到一份規模較大的數據集時，一方面可以通過觀察特徵值排序結果來決定，另一方面還需要通過實驗交叉驗證。
小編食言了，上一篇結尾說好「t檢驗」沒有發出來，因此這次就不「預告」了，下一篇依然是以降維為主題，「t檢驗」也會在後期結合醫學數據分析實戰慢慢補上。