線性判別分析（LDA）數據降維及案例實戰

一、LDA是什麼

二、計算散佈矩陣

三、線性判別式及特徵選擇

四、樣本數據降維投影

五、完整程式碼

結語

線性判別分析（LDA）數據降維及案例實戰

一、LDA是什麼

LDA概念及與PCA區別

LDA線性判別分析（Linear Discriminant Analysis）也是一種特徵提取、數據壓縮技術。在模型訓練時候進行LDA數據處理可以提高計算效率以及避免過擬合。它是一種有監督學習演算法。

與PCA主成分分析（Principal Component Analysis）相比，LDA是有監督數據壓縮方法，而PCA是有監督數據壓縮及特徵提取方法。PCA目標是尋找數據集最大方差方向作爲主成分，LDA目標是尋找和優化具有可分性特徵子空間。其實兩者各有優勢，更深入詳細的區分和應用等待之後的學習，這裏我仍然以葡萄酒數據集分類爲案例記錄原理知識的學習和具體實現步驟。

對比我之前記錄的PCA請看：PCA數據降維原理及python應用（葡萄酒案例分析）

LDA內部邏輯實現步驟

標準化d維數據集。
計算每個類別的d維均值向量。
計算跨類散佈矩陣 $\large S_{B}$ 和類內散佈矩陣 $\large S_{W}$ .
線性判別式及特徵計算。
按特徵值降序排列，與對應的特徵向量成對排序。
選擇最具線性判別性的前k個特徵，構建變換矩陣 $\large W_{d\times k}$ .
通過變換矩陣將原數據投影至k維子空間。

二、計算散佈矩陣

1、數據集下載

下載葡萄酒數據集wine.data到本地，或者到時在載入數據程式碼是從遠端伺服器獲取，爲了避免載入超時推薦下載本地數據集。

下載之後用記事本開啓wine.data可見得，第一列爲葡萄酒數據類別標籤，共有3類，往後的13列爲特徵值。

數據載入以及標準化數據處理與PCA技術一樣，具體可以翻看PCA數據降維原理及python應用（葡萄酒案例分析），或者本文第五部分完整程式碼有具體實現程式碼。

2、計算散佈矩陣第一步，先計算每個類別每個樣本的均值向量。

公式： $\large m_{i}=\frac{1}{n_{i}}\sum_{x \in D}^{c}x_{m}$ , i =1,2,3 表示類別，每個特徵取平均值。

得到三個均值向量爲： $\large m_{i}=\begin{bmatrix} \\ m_{i,1} \\ m_{i,2} \\ m_{i,3} \\ . \\ . \\ . \\ \\ m_{i,13} \end{bmatrix}$

程式碼實現：

# 計算均值向量
np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1, 4):
    mean_vecs.append(np.mean(x_train_std[y_train == label], axis=0))

列印檢視結果：

3、計算類內散佈矩陣。

每個樣本 i 的散佈矩陣：

$\large S_{i}=\sum_{x\in D}(x-m_i)(x-m_i)^{T}$

類內散佈矩陣即每個樣本的累加：

$\large S_{W}=\sum_{i=1}^{13}S_i$

程式碼實現：

# 計算類內散佈矩陣
k = 13
Sw = np.zeros((k, k))
for label, mv in zip(range(1, 4), mean_vecs):
    Si = np.zeros((k, k))
    Si = np.cov(x_train_std[y_train == label].T)  
    Sw += Si
print("類內散佈矩陣：",Sw.shape[0],"*",Sw.shape[1])

矩陣規模：

4、計算跨類散佈矩陣。

公式： $\large S_{B}=\sum_{1}^{13}n_i(m_i-m)(m_i-m)^{T}$

公式中，m是所有樣本總均值向量，也就是不分類的情況下計算特徵平均值。

程式碼實現：

# 計算跨類散佈矩陣
mean_all = np.mean(x_train_std, axis=0)
Sb = np.zeros((k, k))
for i, col_mv in enumerate(mean_vecs):
    n = x_train[y_train == i + 1, :].shape[0]
    col_mv = col_mv.reshape(k, 1)  # 列均值向量
    mean_all = mean_all.reshape(k, 1)
    Sb += n * (col_mv - mean_all).dot((col_mv - mean_all).T)

三、線性判別式及特徵選擇

LDA其他步驟與PCA相似，但是，PCA是分解協方差矩陣提取特徵值，LDA則是求解矩陣 $\large S_{W}^{-1}S_B$ 得到廣義特徵值，實現：

# 計算廣義特徵值
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(Sw).dot(Sb))
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)

對特徵值降序排列之後，列印看看：

print("特徵值降序排列：")
for eigen_val in eigen_pairs:
    print(eigen_val[0])

從捕捉到的特徵值發現，前兩個可以佔據大部分數據集特徵了，接下來視覺化表示更加直觀地觀察：

# 線性判別捕捉，計算辨識力
tot = sum(eigen_vals.real)
discr = []
# discr=[(i/tot) for i in sorted(eigen_vals.real,reverse=True)]
for i in sorted(eigen_vals.real, reverse=True):
    discr.append(i / tot)
# print(discr)
cum_discr = np.cumsum(discr)  # 計算累加方差
plt.rcParams['font.sans-serif'] = ['SimHei']  # 顯示中文
plt.bar(range(1,14),discr,alpha=0.5,align='center',label='獨立辨識力')
plt.step(range(1,14),cum_discr,where='mid',label='累加辨識力')
plt.ylabel('"辨識力"比')
plt.xlabel('線性判別')
plt.ylim([-0.1,1.1])
plt.legend(loc='best')
plt.show()

很明顯，最具線性判別的前兩個特徵捕捉了100%的資訊，下面下麪以此構建變換矩陣 W.

四、樣本數據降維投影

構建變換矩陣：

# 變換矩陣
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
print(w)

來瞅瞅，這就是前兩個特徵向量的矩陣表示。

現在有了變換矩陣，就可以將樣本訓練數據投影到降維特徵空間了： $\large X^{'}=XW$ . 並展示分類結果：

# 樣本數據投影到低維空間
x_train_lda = x_train_std.dot(w)
colors = ['r', 'g', 'b']
marks = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, marks):
    plt.scatter(x_train_lda[y_train == l, 0],
                x_train_lda[y_train == l, 1] * -1,
                c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.show()

很明顯，三個類別線性可分，效果也不錯，相比較PCA方法，我感覺LDA分類結果更好，我們知道，LDA是有監督的方法，有利用到數據集的標籤。

五、完整程式碼

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# load data
df_wine = pd.read_csv('D:\\PyCharm_Project\\maching_learning\\wine_data\\wine.data', header=None)  # 本地載入
# df_wine=pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',header=None)#伺服器載入

# split the data，train：test=7:3
x, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, stratify=y, random_state=0)

# standardize the feature 標準化單位方差
sc = StandardScaler()
x_train_std = sc.fit_transform(x_train)
x_test_std = sc.fit_transform(x_test)

# 計算均值向量
np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1, 4):
    mean_vecs.append(np.mean(x_train_std[y_train == label], axis=0))
    # print("Mean Vectors %s:" % label,mean_vecs[label-1])

# 計算類內散佈矩陣
k = 13
Sw = np.zeros((k, k))
for label, mv in zip(range(1, 4), mean_vecs):
    Si = np.zeros((k, k))
    # for row in x_train_std[y_train==label]:
    #     row,mv=row.reshape(n,1),mv.reshape(n,1)
    #     Si+=(row-mv).dot((row-mv).T)
    Si = np.cov(x_train_std[y_train == label].T)
    Sw += Si
# print("類內散佈矩陣：",Sw.shape[0],"*",Sw.shape[1])
# print("類內標籤分佈：",np.bincount(y_train)[1:])

# 計算跨類散佈矩陣
mean_all = np.mean(x_train_std, axis=0)
Sb = np.zeros((k, k))
for i, col_mv in enumerate(mean_vecs):
    n = x_train[y_train == i + 1, :].shape[0]
    col_mv = col_mv.reshape(k, 1)  # 列均值向量
    mean_all = mean_all.reshape(k, 1)
    Sb += n * (col_mv - mean_all).dot((col_mv - mean_all).T)
# print("跨類散佈矩陣：", Sb.shape[0], "*", Sb.shape[1])

# 計算廣義特徵值
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(Sw).dot(Sb))
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
# print(eigen_pairs[0][1][:,np.newaxis].real) # 第一特徵向量
# print("特徵值降序排列：")
# for eigen_val in eigen_pairs:
#     print(eigen_val[0])

# 線性判別捕捉，計算辨識力
tot = sum(eigen_vals.real)
discr = []
# discr=[(i/tot) for i in sorted(eigen_vals.real,reverse=True)]
for i in sorted(eigen_vals.real, reverse=True):
    discr.append(i / tot)
# print(discr)
cum_discr = np.cumsum(discr)  # 計算累加方差
# plt.rcParams['font.sans-serif'] = ['SimHei']  # 顯示中文
# plt.bar(range(1,14),discr,alpha=0.5,align='center',label='獨立辨識力')
# plt.step(range(1,14),cum_discr,where='mid',label='累加辨識力')
# plt.ylabel('"辨識力"比')
# plt.xlabel('線性判別')
# plt.ylim([-0.1,1.1])
# plt.legend(loc='best')
# plt.show()

# 轉換矩陣
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
# print(w)

# 樣本數據投影到低維空間
x_train_lda = x_train_std.dot(w)
colors = ['r', 'g', 'b']
marks = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, marks):
    plt.scatter(x_train_lda[y_train == l, 0],
                x_train_lda[y_train == l, 1] * -1,
                c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.show()

結語

這篇記錄了這幾天學習的LDA實現數據降維的方法，仍然以葡萄酒數據集爲案例，在上面一步步的拆分中，我們更加清楚線性判別分析LDA方法的內部實現，在這個過程，對於初步學習的我感覺能夠認識和理解更深刻，當然以後數據處理使用LDA方法時候會用到一些第三方庫的類，實現起來更加方便，加油學習，期待下一篇LDA實現更簡潔的方法！

我的部落格園：

我的CSDN部落格：https://blog.csdn.net/Charzous/article/details/108007441

LDA數據壓縮原理及python應用（葡萄酒案例分析）