机器学习笔记——LDA线性判别分析

前言

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

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

参考资料：

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

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

一、算法推导

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

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

当 y i ∈ { 0 , 1 } y_i \in \{0, 1\} yi​∈{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 _{x\in X_0}(x-{\mu }_0{)}^T+\sum _{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$，问题转化为

$$mi{n}_{\omega }-{\omega }^T{S}_b\omega$$

$$st.{\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 }=-2S_b\omega +2\lambda S_{\omega }\omega =0$

得 $S_b\omega =\lambda S_{\omega }\omega \stackrel{◯}{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)\stackrel{◯}{2}$$

由$\stackrel{◯}{1}\stackrel{◯}{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 _{i=1}^N{m}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T$$

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

$$S_{\omega }=\sum _{i=1}^N\sum _{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 _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T$$

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

$\begin{array}{rl}{S}_b& =\sum _{i=1}^N{m}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T\\ & =\sum _{i=1}^N{m}_i({\mu }_i{\mu }_i^T-\mu {\mu }_i^T-{\mu }_i{\mu }^T+\mu {\mu }^T)\end{array}$

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

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

又有

$\begin{array}{rl}{S}_{\omega }& =\sum _{i=1}^N\sum _{x\in X_i}(x-{\mu }_i)(x-{\mu }_i{)}^T\\ & =\sum _{i=1}^N\sum _{x\in X_i}(x{x}^T-{\mu }_i{x}^T-x{\mu }_i^T+{\mu }_i{\mu }_i^T)\end{array}$

那么，

$S_b+S_{\omega }=\sum _{i=1}^N\sum _{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 _{i=1}^N\sum _{x\in X_i}(-{\mu }_i{x}^T-x{\mu }_i^T+2{\mu }_i{\mu }_i^T)=0$

则有

$\begin{array}{rl}{S}_b+S_{\omega }& =\sum _{i=1}^N\sum _{x\in X_i}(x{x}^T-\mu x^T-x{\mu }^T+\mu {\mu }^T)\\ & =\sum _{i=1}^N\sum _{x\in X_i}(x-\mu )(x-\mu {)}^T\\ & =\sum _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T=S_t\end{array}$

证毕

3. 等价模型表示

以上的算法推导是基于$ma{x}_{\omega }\frac{{\omega }^T{S}_b\omega }{{\omega }^T{S}_{\omega }\omega }$，或者说，基于$ma{x}_W\frac{tr(W^T{S}_{b}W)}{tr(W^T{S}_{\omega }W)}$，我们称之为LDA的除法模型。

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

(1) 减法模型：

$$ma{x}_{W}tr(W^T{S}_{b}W)-tr(W^T{S}_{\omega }W)$$

$$st.tr(W^{T}W)=1$$

求解得到

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

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

(2) 除法正则模型：

$$ma{x}_W\frac{tr(W^T{S}_{b}W)}{tr(W^T{S}_{\omega }W)}+\lambda tr(W^{T}W)$$

$$st.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}=2S_{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) 减法正则模型：

$$ma{x}_{W}tr(W^T{S}_{b}W)-tr(W^T{S}_{\omega }W)-\lambda tr(W^{T}W)$$

$$st.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貌似不支持(?)，得手搓，急着交作业就没再搞了。。