别再死记公式了！用Python手撸一个LDA分类器，从鸢尾花数据集开始

从零实现LDA分类器：用Python代码解锁鸢尾花数据集的分类奥秘

当我在第一次接触线性判别分析(LDA)时，那些复杂的数学公式让我望而生畏。直到有一天，我决定抛开教科书，直接用Python代码一步步实现这个算法，才真正理解了它的精妙之处。本文将带你体验这段从理论到实践的旅程，我们不仅会用NumPy手写LDA的核心计算，还会在经典的鸢尾花数据集上进行可视化展示，让你获得"亲手实现"的成就感。

1. 环境准备与数据探索

在开始编码前，我们需要准备好Python环境和必要的数据集。推荐使用Anaconda创建虚拟环境，这能避免包依赖冲突：

conda create -n lda_demo python=3.8 conda activate lda_demo pip install numpy matplotlib scikit-learn

鸢尾花数据集是机器学习领域的"Hello World"，它包含三种鸢尾花的四个特征：萼片长度、萼片宽度、花瓣长度和花瓣宽度。我们先加载并观察数据：

from sklearn.datasets import load_iris import numpy as np iris = load_iris() X = iris.data y = iris.target feature_names = iris.feature_names print(f"特征矩阵形状: {X.shape}") print(f"类别分布: {np.bincount(y)}") print("特征名称:", feature_names)

输出结果会显示我们有150个样本，均匀分布在三个类别中。为了简化问题，我们暂时只考虑二分类场景，选择前两个类别：

# 提取前两类数据 X = X[y != 2] y = y[y != 2]

2. LDA核心算法实现

LDA的核心思想是找到最佳投影方向，使得同类样本尽可能聚集，不同类样本尽可能分离。这需要计算两个关键矩阵：

2.1 计算类内散度矩阵

类内散度矩阵(S_w)衡量同类样本的离散程度。我们需要先计算每个类别的协方差矩阵，然后求和：

def compute_within_class_scatter(X, y): classes = np.unique(y) n_features = X.shape[1] S_w = np.zeros((n_features, n_features)) for c in classes: X_c = X[y == c] # 计算类内协方差矩阵，注意ddof=1表示无偏估计 cov_matrix = np.cov(X_c, rowvar=False, ddof=1) S_w += cov_matrix return S_w S_w = compute_within_class_scatter(X, y) print("类内散度矩阵:\n", S_w)

2.2 计算类间散度矩阵

类间散度矩阵(S_b)衡量不同类别中心点的离散程度：

def compute_between_class_scatter(X, y): overall_mean = np.mean(X, axis=0) classes = np.unique(y) n_features = X.shape[1] S_b = np.zeros((n_features, n_features)) for c in classes: n_c = np.sum(y == c) mean_c = np.mean(X[y == c], axis=0) # 计算每个类别的贡献 diff = (mean_c - overall_mean).reshape(-1, 1) S_b += n_c * np.dot(diff, diff.T) return S_b S_b = compute_between_class_scatter(X, y) print("类间散度矩阵:\n", S_b)

2.3 求解最优投影方向

LDA的最优投影方向是最大化广义瑞利商的方向，这可以通过求解广义特征值问题得到：

def lda(X, y): S_w = compute_within_class_scatter(X, y) S_b = compute_between_class_scatter(X, y) # 数值稳定的求逆方法 eigenvalues, eigenvectors = np.linalg.eig(np.linalg.pinv(S_w) @ S_b) # 选择最大特征值对应的特征向量 idx = eigenvalues.argsort()[::-1] eigenvectors = eigenvectors[:, idx] w = eigenvectors[:, 0].real # 归一化投影向量 w = w / np.linalg.norm(w) return w w = lda(X, y) print("最优投影方向:", w)

在实际项目中，我们可能会遇到矩阵奇异的问题。这时可以添加一个小的正则化项来保证数值稳定性：

S_w_reg = S_w + 1e-6 * np.eye(S_w.shape[0])

3. 结果可视化与分析

理解算法最好的方式就是看到它的实际效果。让我们将原始数据和投影后的结果可视化：

import matplotlib.pyplot as plt # 原始数据散点图 plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) for c in np.unique(y): plt.scatter(X[y == c, 0], X[y == c, 1], label=f'Class {c}') plt.xlabel('Sepal length') plt.ylabel('Sepal width') plt.title('Original Data') plt.legend() # 计算投影后的数据 projected = X @ w.reshape(-1, 1) # 投影结果可视化 plt.subplot(1, 2, 2) for c in np.unique(y): plt.hist(projected[y == c], alpha=0.7, label=f'Class {c}', bins=20) plt.xlabel('Projection Value') plt.ylabel('Frequency') plt.title('LDA Projection Results') plt.legend() plt.tight_layout() plt.show()

这段代码会生成两个子图：左边显示原始特征空间中的数据分布，右边显示投影到最佳方向后的分布情况。你会看到投影后的两类数据明显更易区分。

4. 性能评估与边界绘制

为了量化我们的LDA实现效果，我们可以计算分类准确率，并绘制决策边界：

from sklearn.metrics import accuracy_score # 计算投影后的类别中心 mean_0 = np.mean(projected[y == 0]) mean_1 = np.mean(projected[y == 1]) # 计算决策阈值 threshold = (mean_0 + mean_1) / 2 # 预测类别 y_pred = (projected > threshold).astype(int).flatten() # 计算准确率 acc = accuracy_score(y, y_pred) print(f"分类准确率: {acc:.2%}") # 绘制决策边界 plt.figure(figsize=(8, 6)) for c in np.unique(y): plt.scatter(X[y == c, 0], X[y == c, 1], label=f'Class {c}') # 生成网格点用于绘制决策边界 x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1 y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1 xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.1), np.arange(y_min, y_max, 0.1)) # 计算网格点的投影值 Z = np.c_[xx.ravel(), yy.ravel()] @ w.reshape(-1, 1) Z = (Z > threshold).astype(int).reshape(xx.shape) # 绘制决策边界 plt.contourf(xx, yy, Z, alpha=0.3, levels=[-0.5, 0.5, 1.5]) plt.xlabel('Sepal length') plt.ylabel('Sepal width') plt.title('LDA Decision Boundary') plt.legend() plt.show()

5. 扩展到多类别场景

虽然我们演示的是二分类场景，但LDA天然支持多类别分类。关键在于计算投影矩阵时选择多个特征向量：

def multiclass_lda(X, y, n_components=None): S_w = compute_within_class_scatter(X, y) S_b = compute_between_class_scatter(X, y) # 解决数值稳定性问题 S_w_reg = S_w + 1e-6 * np.eye(S_w.shape[0]) eigenvalues, eigenvectors = np.linalg.eig(np.linalg.pinv(S_w_reg) @ S_b) # 按特征值降序排序 idx = eigenvalues.argsort()[::-1] eigenvectors = eigenvectors[:, idx] eigenvalues = eigenvalues[idx] # 选择前n_components个特征向量 if n_components is not None: eigenvectors = eigenvectors[:, :n_components] return eigenvectors.real # 使用全部三类数据 X_full = iris.data y_full = iris.target W = multiclass_lda(X_full, y_full, n_components=2) print("投影矩阵:\n", W) # 投影数据 X_lda = X_full @ W # 可视化投影结果 plt.figure(figsize=(8, 6)) for c in np.unique(y_full): plt.scatter(X_lda[y_full == c, 0], X_lda[y_full == c, 1], label=f'Class {c}') plt.xlabel('LD1') plt.ylabel('LD2') plt.title('Multiclass LDA Projection') plt.legend() plt.show()

在多类别场景中，LDA实际上是将数据投影到一个最多有C-1维的子空间（C是类别数）。对于三分类问题，我们最多可以得到两个有意义的线性判别方向。

6. 实际应用中的注意事项

在真实项目中应用LDA时，有几个关键点需要特别注意：

• 特征缩放：虽然LDA不受特征尺度影响（因为涉及协方差计算），但为了数值稳定性，建议标准化数据：

from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X)

• 维度灾难：当特征维度很高而样本数较少时，S_w可能奇异。解决方法包括：

  • 增加正则化项
  • 先使用PCA降维
  • 使用伪逆代替常规逆

• 类别不平衡：原始LDA假设各类别分布均匀。对于不平衡数据，可以调整类间散度矩阵的计算方式：

# 考虑类别权重的S_b计算 def weighted_between_class_scatter(X, y): overall_mean = np.mean(X, axis=0) classes, counts = np.unique(y, return_counts=True) n_features = X.shape[1] S_b = np.zeros((n_features, n_features)) total_samples = len(y) for c, n_c in zip(classes, counts): mean_c = np.mean(X[y == c], axis=0) diff = (mean_c - overall_mean).reshape(-1, 1) weight = n_c / total_samples # 使用类别比例作为权重 S_b += weight * np.dot(diff, diff.T) return S_b

• 与PCA的区别：虽然PCA和LDA都是线性降维方法，但它们的优化目标不同：

  特性	PCA	LDA
  优化目标	最大化方差	最大化类间分离度
  监督/无监督	无监督	有监督
  投影方向数	由用户指定	最多C-1个（C为类别数）
  适用场景	探索性分析、特征提取	分类任务、可解释性降维

在完成这个项目后，我最大的收获是理解了数学公式背后的实际意义。比如，类内散度矩阵实际上是在量化每个类别的"紧凑程度"，而类间散度矩阵则衡量不同类别的"分离程度"。这种从代码实现反推理论理解的方式，比单纯学习数学推导要直观得多。