企业官网建设流程全解析

用Python手把手复现LDA二分类例题：从协方差矩阵到投影结果（附完整代码）

线性判别分析（LDA）作为经典的监督降维算法，在金融风控、生物特征识别等领域有着广泛应用。但对于初学者而言，理论推导与代码实现之间往往存在断层——你可能已经理解了最大化类间距离、最小化类内距离的数学原理，却不知道如何用NumPy实现协方差矩阵计算；或者能手动解出投影向量，但面对实际数据集时无从下手。本文将以一个二维特征空间的二分类问题为例，带你用Python完整复现LDA从数据预处理到结果可视化的全流程。

1. 环境准备与数据加载

首先确保你的Python环境已安装以下库：

import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import make_blobs # 用于生成演示数据

我们使用与原始例题相同的二维数据集，包含10个样本点（每个类别5个）：

# 负类样本（类别0） X_neg = np.array([[4,2], [2,4], [2,3], [3,6], [4,4]]) # 正类样本（类别1） X_pos = np.array([[9,10], [6,8], [9,5], [8,7], [10,8]]) X = np.vstack((X_neg, X_pos)) # 合并特征矩阵 y = np.array([0]*5 + [1]*5) # 标签向量

通过散点图直观观察数据分布：

plt.scatter(X_neg[:,0], X_neg[:,1], c='blue', label='Class 0') plt.scatter(X_pos[:,0], X_pos[:,1], c='red', label='Class 1') plt.xlabel('Feature 1'); plt.ylabel('Feature 2') plt.legend(); plt.grid() plt.show()

2. 核心矩阵计算实现

2.1 计算类均值向量

LDA的第一步是计算每个类别的均值向量。通过NumPy的mean函数可高效实现：

mu_0 = np.mean(X_neg, axis=0) # 负类均值 [3. , 3.8] mu_1 = np.mean(X_pos, axis=0) # 正类均值 [8.4, 7.6]

2.2 构建类间散度矩阵Sb

类间散度矩阵反映类别中心的差异，计算公式为：

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

对应的Python实现：

mu_diff = mu_0 - mu_1 # 均值差向量 Sb = np.outer(mu_diff, mu_diff) # 外积运算 print("Sb矩阵:\n", Sb)

输出结果应与手动计算一致：

[[29.16 22.68] [22.68 17.64]]

2.3 计算类内散度矩阵Sw

类内散度矩阵需要分别计算每个类的协方差矩阵。我们封装一个covariance_matrix函数：

def covariance_matrix(X): n_samples = X.shape[0] X_centered = X - np.mean(X, axis=0) return (X_centered.T @ X_centered) / (n_samples - 1) Sigma0 = covariance_matrix(X_neg) # 负类协方差 Sigma1 = covariance_matrix(X_pos) # 正类协方差 Sw = Sigma0 + Sigma1 # 类内散度矩阵

验证输出：

Sw矩阵: [[ 3.3 -0.3 ] [-0.3 5.5 ]]

3. 求解投影方向

3.1 矩阵求逆与特征分解

LDA的最优投影方向是广义特征方程$S_w^{-1}S_b\omega = \lambda\omega$的最大特征值对应特征向量：

Sw_inv = np.linalg.inv(Sw) # 矩阵求逆 A = Sw_inv @ Sb # 矩阵乘积 eigvals, eigvecs = np.linalg.eig(A) # 特征分解

3.2 选择最优投影向量

提取最大特征值对应的特征向量并进行归一化：

max_idx = np.argmax(eigvals) # 最大特征值索引 w = eigvecs[:, max_idx] # 最优投影方向 w = w / np.linalg.norm(w) # 单位向量化 print("投影向量:", w)

得到的投影方向约为[-0.83, -0.56]，与手动计算结果一致。

4. 数据投影与可视化

4.1 计算投影坐标

将原始数据投影到判别方向上：

projection = X @ w # 矩阵乘法计算投影值

4.2 可视化投影结果

绘制原始数据点及其在投影直线上的映射：

# 绘制原始数据 plt.scatter(X_neg[:,0], X_neg[:,1], c='blue', label='Class 0') plt.scatter(X_pos[:,0], X_pos[:,1], c='red', label='Class 1') # 绘制投影直线 x_line = np.linspace(0, 11, 100) y_line = (w[1]/w[0]) * x_line plt.plot(x_line, y_line, 'k--', label='Projection Line') # 绘制投影连线 for point, proj in zip(X, projection): plt.plot([point[0], proj*w[0]], [point[1], proj*w[1]], 'gray', alpha=0.3) plt.xlabel('Feature 1'); plt.ylabel('Feature 2') plt.legend(); plt.grid(); plt.axis('equal') plt.show()

4.3 一维分布直方图

观察投影后的类间分离效果：

plt.hist(projection[y==0], bins=5, alpha=0.5, color='blue', label='Class 0') plt.hist(projection[y==1], bins=5, alpha=0.5, color='red', label='Class 1') plt.xlabel('Projection Value'); plt.ylabel('Frequency') plt.legend(); plt.show()

5. 完整代码与扩展建议

将上述步骤整合为完整函数：

def lda_manual(X, y): # 分割类别 X0 = X[y==0]; X1 = X[y==1] # 计算均值 mu0 = np.mean(X0, axis=0) mu1 = np.mean(X1, axis=0) # 计算散度矩阵 Sb = np.outer((mu0 - mu1), (mu0 - mu1)) Sw = covariance_matrix(X0) + covariance_matrix(X1) # 特征分解 A = np.linalg.inv(Sw) @ Sb _, eigvecs = np.linalg.eig(A) w = eigvecs[:, np.argmax(np.abs(eigvals))] return w / np.linalg.norm(w)

实际应用中可考虑以下优化：

• 添加正则化项防止Sw矩阵奇异
• 扩展多分类场景（使用OvR或OvO策略）
• 与PCA结合进行分层降维