Fisher投影加速分类

线性判别分析（LDA）详解：从理论到代码实现

一、引言

线性判别分析（Linear Discriminant Analysis, LDA）是一种经典的监督学习降维与分类方法。其核心思想是通过最大化类间散度、最小化类内散度，将高维数据投影到低维空间并实现高效分类。本文以鸢尾花数据集为例，详解LDA的数学原理，并逐行解析代码实现。

二、LDA的核心原理

1. 核心目标

LDA的目标是找到一个最优投影方向 $ W $，使得以下准则成立：

类内散度最小：同一类样本在投影后的空间尽可能紧凑。
类间散度最大：不同类样本在投影后的空间尽可能分离。

数学表达式为：

$$
J(W) = \frac{W^T S_B W}{W^T S_W W}
$$

其中：

$ S_W $：类内散度矩阵（Within-class Scatter Matrix）
$ S_B $：类间散度矩阵（Between-class Scatter Matrix）

2. 数学推导

(1) 类均值向量

对每个类别 $ c $，计算其特征均值向量 $ \mu_c $：

$$
\mu_c = \frac{1}{N_c} \sum_{x \in C_c} x
$$

其中 $ N_c $ 是类别 $ c $ 的样本数。

(2) 类内散度矩阵 $ S_W $

衡量同一类样本的离散程度：

$$
S_W = \sum_{c} \sum_{x \in C_c} (x - \mu_c)(x - \mu_c)^T
$$

(3) 类间散度矩阵 $ S_B $

衡量不同类中心之间的距离：

$$
S_B = \sum_{c} N_c (\mu_c - \mu)(\mu_c - \mu)^T
$$

其中 $ \mu $ 是全局均值向量。

(4) 投影矩阵计算

通过广义特征值问题求解最优投影方向：

$$
S_W^{-1} S_B w = \lambda w
$$

选择最大特征值对应的特征向量 $ w $ 作为投影方向。

三、代码详解与实现

1. 数据准备与预处理

from sklearn.datasets import load_iris
iris = load_iris()
X = iris.data[y != 2]  # 筛选前两类（0和1）
y = y[y != 2]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

功能：加载鸢尾花数据集并筛选前两类，构造二分类任务。
原理：LDA适用于多分类，但二分类更简单直观。

2. 类均值计算

def class_means(X, y):
    classes = np.unique(y)
    means = [np.mean(X[y == c], axis=0) for c in classes]
    return np.array(means)

功能：计算每个类别的均值向量。
数学对应：实现公式 $ \mu_c $ 的计算。

3. 类内散度矩阵

def within_class_scatter(X, y):
    means = class_means(X, y)
    num_features = X.shape[1]
    S_w = np.zeros((num_features, num_features))
    for c, mean in zip(np.unique(y), means):
        for x in X[y == c]:
            x = x.reshape(-1, 1)
            mean = mean.reshape(-1, 1)
            S_w += (x - mean).dot((x - mean).T)
    return S_w

功能：计算类内散度矩阵 $ S_W $。
关键点：通过外积累加每个样本与类均值的偏差。

4. 类间散度矩阵

def between_class_scatter(X, y):
    means = class_means(X, y)
    mean1, mean2 = means[0].reshape(-1, 1), means[1].reshape(-1, 1)
    return (mean1 - mean2).dot((mean1 - mean2).T)

功能：计算类间散度矩阵 $ S_B $。
局限性：仅适用于二分类，多分类需扩展。

5. 投影矩阵求解

def projection_matrix(X, y):
    S_w = within_class_scatter(X, y)
    S_b = between_class_scatter(X, y)
    eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_w).dot(S_b))
    return eig_vecs[:, np.argsort(eig_vals)[::-1][:1]]

功能：求解广义特征值问题，取最大特征值对应的特征向量作为投影方向。
数学对应：求解 $ S_W^{-1} S_B $ 的最大特征值对应的特征向量。

6. 数据投影与分类

W = projection_matrix(X_train, y_train)
X_train_lda = X_train.dot(W)  # 降维后的训练数据
X_test_lda = X_test.dot(W)    # 降维后的测试数据

# 计算阈值并分类
threshold = (X_train_lda[y_train == 0].mean() + X_train_lda[y_train == 1].mean()) / 2
y_pred = (X_test_lda >= threshold).astype(int)

决策规则：使用两类均值的中点作为阈值进行分类。

四、实验结果与可视化

1. 可视化设计

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6))
# 绘制训练集、测试集及合并结果
for ax, X_lda, y_true, title in zip([ax1, ax2, ax3], 
                                   [X_train_lda, X_test_lda, X_train_lda], 
                                   [y_train, y_test, y_train], 
                                   ["Train", "Test", "Combined"]):
    for cls, color in zip(np.unique(y), ['r', 'b']):
        ax.scatter(X_lda[y_true == cls], np.zeros(np.sum(y_true == cls)), 
                   c=color, label=f'Class {cls}')
    ax.axvline(x=threshold, c='g', linestyle='--')
    ax.set_title(f"{title} (Acc: {accuracy_score(y_true, y_pred):.2%})")
plt.legend()
plt.show()