上一篇寫過線性判別分析處理二分類問題https://blog.csdn.net/z962013489/article/details/79871789，當使用LDA處理多分類問題時，通常是作為一個降維工具來使用的。若我們有一個D維的樣本集，該樣本集包含C個類別共n個樣本，希望將D維降維成K維。之前在二分類問題中，我們定義的類間散度矩陣為：

Sb=(μ1−μ2)(μ1−μ2)T$$S_b=({\mu }_1-{\mu }_2)({\mu }_1-{\mu }_2{)}^T$$

當類別為3時就已經不再適用，在這裡引出全域性散度矩陣的概念：

St=Sw+Sb=∑i=1n(xi−μ)(xi−μ)T$$S_t=S_w+S_b=\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^T$$

其中 μ$\mu$ 為整個樣本集的均值向量，其中的類內散度矩陣

Sw$S_w$定義為所有類別的散度矩陣之和，與二分類類似，基於全域性散度矩陣，我們就可以求得類間散度矩陣Sb$S_b$定義為：

Sb=St−Sw=∑i=1nmi(μi−μ)(μi−μ)T$$S_b=S_t-S_w=\sum _{i=1}^n{m}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T$$
其中mi$m_i$為第i個類別樣本總數，μi${\mu }_i$為第i個類別樣本的均值向量。

上圖的例子可以作為參考便於理解。
我們要計算一個W矩陣，使樣本向量在該投影矩陣的作用下能夠實現其與同類樣本的中心距離儘量小，而與異類樣本的中心距離儘量大。常見的一種優化目標為：

maxWtr(WTSbW)tr(WTSwW)$$\underset{W}{max}\frac{tr(W^T{S}_{b}W)}{tr(W^T{S}_{w}W)}$$
其中的tr()為矩陣的跡，一個n×n的對角矩陣A的主對角線（從左上方至右下方的對角線）上各個元素的總和被稱為矩陣A的跡（或跡數），一般記作tr(A)。
這個優化目標實際上等價於求解多個w組合成W，那麼該問題就等價於求解多個上一章的優化目標，使用相同的方法，可以求得下式：

SbW=λSwW$$S_{b}W=\lambda S_{w}W$$
即S−1wSbW=λW$$S_w^{-1}{S}_{b}W=\lambda W$$
W的閉式解為S−1wSb$S_w^{-1}{S}_b$的k個最大非零廣義特徵值對應的特徵向量組成的矩陣，LDA降維最多可以降至C-1，C為樣本類別數，與原始特徵維數n無關。
使用python實現LDA降維程式碼，對Iris資料集從四維降至二維，繪圖如下：
其中上圖為自己實現的方法降維表現，下圖為sklearn自帶方法降維表現。

python3.6實現如下：

import numpy as np
from sklearn.datasets import load_iris
import matplotlib.pyplot as
 
 plt
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis


def LDA_dimensionality(X, y, k):
    '''
    X為資料集，y為label，k為目標維數
    '''
    label_ = list(set(y))

    X_classify = {}

    for label in label_:
        X1 = np.array([X[i] for i in range(len(X)) if y[i] == label])
        X_classify[label] = X1

    mju = np.mean(X, axis=0)
    mju_classify = {}

    for label in label_:
        mju1 = np.mean(X_classify[label], axis=0)
        mju_classify[label] = mju1

    #St = np.dot((X - mju).T, X - mju)

    Sw = np.zeros((len(mju), len(mju)))  # 計算類內散度矩陣
    for i in label_:
        Sw += np.dot((X_classify[i] - mju_classify[i]).T,
                     X_classify[i] - mju_classify[i])

    # Sb=St-Sw

    Sb = np.zeros((len(mju), len(mju)))  # 計算類內散度矩陣
    for i in label_:
        Sb += len(X_classify[i]) * np.dot((mju_classify[i] - mju).reshape(
            (len(mju), 1)), (mju_classify[i] - mju).reshape((1, len(mju))))

    eig_vals, eig_vecs = np.linalg.eig(
        np.linalg.inv(Sw).dot(Sb))  # 計算Sw-1*Sb的特徵值和特徵矩陣

    sorted_indices = np.argsort(eig_vals)
    topk_eig_vecs = eig_vecs[:, sorted_indices[:-k - 1:-1]]  # 提取前k個特徵向量
    return topk_eig_vecs


if '__main__' == __name__:

    iris = load_iris()
    X = iris.data
    y = iris.target

    W = LDA_dimensionality(X, y, 2)
    X_new = np.dot((X), W)
    plt.figure(1)
    plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)

    # 與sklearn中的LDA函式對比
    lda = LinearDiscriminantAnalysis(n_components=2)
    lda.fit(X, y)
    X_new = lda.transform(X)
    print(X_new)
    plt.figure(2)
    plt.scatter(X_new[:, 0], X_new[:, 1], marker='o', c=y)

    plt.show()