前言

机器学习的东西真是又多又杂啊，之前刷吴恩达就见过主成分分析，即PCA(Principle Component Analysis)，没见过线性判别分析，即LDA(Linear Discriminant Analysis)。才知道前者是对无标签数据作压缩的，而后者可以对有标签的二分类、多分类任务的数据作压缩。

本文主要讲述了LDA算法的推导，以及展示了LDA算法的代码实现和应用。

参考资料：

周志华《机器学习》，人称西瓜书

更详细的算法推导：机器学习算法推导&手写实现03——线性判别分析 - 知乎 (zhihu.com)

一、算法推导

由前言可知，LDA相比于PCA，是用来处理有标签数据的，为了使压缩后的数据仍然能表示分类的情况，我们应该让同类的样本尽可能靠近，异类的样本尽可能远离。

数学表示为：给定数据集$D = \{x_i, y_i\}_{i=1}^m$，该数据集共有$m$个样本，每个样本包含一个数据，即一个$n$维的向量$x_i$，以及一个标签，即一个数$y_i$，代表该数据点所在的类别。

当$y_i \in \{0, 1\}$，说明这是一个二分类任务，我们先从简单的二分类开始讨论。

1. 二分类

1.1 优化目标

接下来，我们要用数学语言来表示同类的样本尽可能靠近，异类的样本尽可能远离这两个任务。

同类的样本尽可能靠近，也就是同类之间的方差要尽可能小

设压缩前，两个类的协方差矩阵分别为$\Sigma_0$和$\Sigma_1$

那么，压缩后的两个类的协方差矩阵为$\omega^T \Sigma_0 \omega$和$\omega^T \Sigma_1 \omega$

定义“类内散度矩阵”（with-class scatter matrix）为

$S_{\omega} = \Sigma_0 + \Sigma_1 = \sum \limits_{x \in X_0} (x - \mu_0)^T + \sum \limits_{x \in X_1} (x - \mu_1)^T$

那么，同类之间的方差要尽可能小，就是要最小化$\omega^T S_{\omega} \omega$

异类的样本尽可能远离，也就是异类样本的均值的差值要尽可能大

设压缩前，两个类的均值分别为$\mu_0$和$\mu_1$

那么，压缩后的两个类的均值为$\omega^T \mu_0$和$\omega^T \mu_1$

定义“类间散度矩阵”（between-class scatter matrix）为

$S_b = (\mu_0 - \mu_1)(\mu_0 - \mu_1)^T$

那么，异类的样本尽可能远离，就是要最大化$\omega^T S_b\omega$

将上述两条件统一到一个式子，我们的优化目标就变为了

最大化 $J = \frac{\omega^T S_b\omega}{\omega^T S_{\omega} \omega}$

1.2 数学推导

由于$J = \frac{\omega^T S_b\omega}{\omega^T S_{\omega} \omega}$分子分母都是关于$\omega$的二次项，因此，求解上述问题与$\omega$的大小无关，只与$\omega$的的方向有关。

我们不妨假设$\omega^T S_{\omega} \omega = 1$，问题转化为

$min_{\omega} \space -\omega^T S_b\omega$

$st. \space \omega^T S_{\omega} \omega = 1$

由拉格朗日乘子法，得

$L(\omega) = -\omega^T S_b\omega + \lambda (\omega^T S_{\omega} \omega - 1)$

令 $\frac{\partial L(\omega)}{\partial \omega} = -2 S_b \omega + 2 \lambda S_{\omega} \omega = 0$

得 $S_b \omega = \lambda S_{\omega} \omega \space \space \textcircled{1}$

又 $S_b \omega = (\mu_0 - \mu_1)(\mu_0 - \mu_1)^T \omega$，其中$(\mu_0 - \mu_1)^T \omega$可看作是个标量，因此$S_b \omega$的方向恒与$\mu_0 - \mu_1$相同，不妨令

$S_b \omega = \lambda (\mu_0 - \mu_1) \space \space \textcircled{2}$

由$\textcircled{1} \textcircled{2}$得

$S_{\omega} \omega = \mu_0 - \mu_1$

$\omega = S_{\omega}^{-1} (\mu_0 - \mu_1)$

在实践中，常常对$S_{\omega}$做奇异值分解：$S_{\omega} = U \Sigma V^T$，故有$S_{\omega}^{-1} = V \Sigma^{-1} U^T$，从而得到$S_{\omega}^{-1}$

2. 多分类

2.1 问题求解

将问题扩展到多分类，定义所有样本的均值为$\mu$，每个类对应的样本数量为$m_i$，那么重新定义“类间散度矩阵”为

$S_b = \sum\limits_{i=1}^N m_i (\mu_i - \mu) ( \mu_i - \mu )^T$

设总共由$N$个类，重新定义“类内散度矩阵”为

$S_{\omega} = \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x - \mu_i) (x - \mu_i )^T$

在上述二分类中，我们将样本投影到一维空间，也就是一条直线上，而在多分类中，我们将样本投影到$N - 1$维空间，记投影矩阵 $W \in R^{d \times (N -1)}$，其中$d$表示原样本的维度数。

同理可得

$S_b W = \lambda S_{\omega} W$

$S_{\omega}^{-1} S_b W = \lambda W$

因此$W$的解是$S_{\omega}^{-1} S_b$的最大$N - 1$个特征值所对应特征向量所组成的矩阵

2.2 全局散度矩阵

除了使用$S_b$和$S_{\omega}$作为度量的标准，多分类LDA还可以有多种不同的实现，引入“全局散度矩阵”

$S_t = \sum \limits_{i=1}^m (x_i - \mu)(x_i - \mu)^T$

可以证明 $S_t = S_b + S_{\omega}$，如下：

$\begin{aligned} S_b &= \sum\limits_{i=1}^N m_i (\mu_i - \mu) ( \mu_i - \mu )^T \\ &= \sum\limits_{i=1}^N m_i (\mu_i \mu_i^T - \mu \mu_i^T - \mu_i \mu^T + \mu \mu^T) \end{aligned}$

其中，$m_i \mu \mu_i^T == \sum\limits_{x \in X_i} \mu x^T$，$m_i \mu_i \mu^T = \sum\limits_{x \in X_i} X \mu^T$，那么有

$S_b = \sum\limits_{i=1}^N \sum\limits_{x \in X_i} (\mu_i \mu_i^T -\mu x^T - x \mu^T + \mu \mu^T)$

又有

$\begin{aligned} S_{\omega} &= \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x - \mu_i) (x - \mu_i )^T \\ &= \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x x^T - \mu_i x^T - x \mu_i^T + \mu_i \mu_i^T) \end{aligned}$

那么，

$S_b + S_{\omega} = \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x x^T -\mu x^T - x \mu^T + \mu \mu^T - \mu_i x^T - x \mu_i^T + 2 \mu_i \mu_i^T)$

其中，$\sum\limits_{i=1}^N \sum \limits_{x \in X_i} (- \mu_i x^T - x \mu_i^T + 2 \mu_i \mu_i^T) = 0$

则有

$\begin{aligned} S_b + S_{\omega} &= \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x x^T -\mu x^T - x \mu^T + \mu \mu^T) \\ &= \sum\limits_{i=1}^N \sum \limits_{x \in X_i} (x - \mu) (x - \mu)^T \\ &= \sum\limits_{i=1}^m (x_i - \mu) (x_i - \mu)^T = S_t \end{aligned}$

证毕

3. 等价模型表示

以上的算法推导是基于$max_{\omega} \frac{\omega^T S_b\omega}{\omega^T S_{\omega} \omega}$，或者说，基于$max_W \frac{ tr(W^T S_b W) }{ tr(W^T S_{\omega} W) }$，我们称之为LDA的除法模型。

实际上，LDA还有多种等价的模型表示，比如

(1) 减法模型：

$max_W tr(W^T S_b W) - tr(W^T S_{\omega} W)$

$st. \space \space tr(W^T W) = 1$

求解得到

$(S_{\omega} - S_b) W = \lambda W$

即$W$是$S_{\omega} - S_b$前$N-1$个最大特征值对应的特征向量组成的矩阵

(2) 除法正则模型：

$max_W \frac{ tr(W^T S_b W) }{ tr(W^T S_{\omega} W) } + \lambda tr(W^T W)$

$st. \space \space tr(W^T S_{\omega} W) = 1$

由拉格朗日乘子法，得

$L = tr(W^T S_b W) + \lambda tr(W^T W) + \beta (tr(W^T S_{\omega} W) - 1)$

$\frac{\partial L}{\partial W} = 2 S_b W + 2 \lambda W + 2 \beta S_{\omega} W = 0$

得

$S_{\omega}^{-1}(S_b + \lambda I) W = -\beta W$

即$W$是$S_{\omega}^{-1}(S_b + \lambda I)$前$N-1$个最大特征值对应的特征向量组成的矩阵

(3) 减法正则模型：

$max_W tr(W^T S_b W) - tr(W^T S_{\omega} W) - \lambda tr(W^T W)$

$st. \space \space tr(W^T W) = 1$

求得问题的解仍然是$S_{\omega} - S_b$前$N-1$个最大特征值对应的特征向量组成的矩阵

二、代码示例

大半都是gpt写的（

X, y = read_data(PATH)

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, random_state=42)

accuracy_pca = []
accuracy_lda = []

n_components_range = range(1, 35, 2)
for n_components in n_components_range:
    # PCA
    pca = PCA(n_components=n_components)
    X_train_pca = pca.fit_transform(X_train)
    X_test_pca = pca.transform(X_test)
    knn_pca = KNeighborsClassifier()
    knn_pca.fit(X_train_pca, y_train)
    y_pred_pca = knn_pca.predict(X_test_pca)
    accuracy_pca.append(accuracy_score(y_test, y_pred_pca))

    # LDA
    lda = LinearDiscriminantAnalysis(n_components=n_components)
    X_train_lda = lda.fit_transform(X_train, y_train)
    X_test_lda = lda.transform(X_test)
    knn_lda = KNeighborsClassifier()
    knn_lda.fit(X_train_lda, y_train)
    y_pred_lda = knn_lda.predict(X_test_lda)
    accuracy_lda.append(accuracy_score(y_test, y_pred_lda))

fig, ax = plt.subplots()
ax.plot(n_components_range, accuracy_pca, label='PCA')
ax.plot(n_components_range, accuracy_lda, label='LDA')
ax.set_xlabel('n_components')
ax.set_ylabel('Accuracy')
ax.legend()
plt.show()

在ORL数据集的测试结果如下：

其实可以再探讨一下减法模型还有正则模型的性能，但是sklearn貌似不支持(?)，得手搓，急着交作业就没再搞了。。