)
用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], cblue, labelClass 0) plt.scatter(X_pos[:,0], X_pos[:,1], cred, labelClass 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, axis0) # 负类均值 [3. , 3.8] mu_1 np.mean(X_pos, axis0) # 正类均值 [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, axis0) 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], cblue, labelClass 0) plt.scatter(X_pos[:,0], X_pos[:,1], cred, labelClass 1) # 绘制投影直线 x_line np.linspace(0, 11, 100) y_line (w[1]/w[0]) * x_line plt.plot(x_line, y_line, k--, labelProjection Line) # 绘制投影连线 for point, proj in zip(X, projection): plt.plot([point[0], proj*w[0]], [point[1], proj*w[1]], gray, alpha0.3) plt.xlabel(Feature 1); plt.ylabel(Feature 2) plt.legend(); plt.grid(); plt.axis(equal) plt.show()4.3 一维分布直方图观察投影后的类间分离效果plt.hist(projection[y0], bins5, alpha0.5, colorblue, labelClass 0) plt.hist(projection[y1], bins5, alpha0.5, colorred, labelClass 1) plt.xlabel(Projection Value); plt.ylabel(Frequency) plt.legend(); plt.show()5. 完整代码与扩展建议将上述步骤整合为完整函数def lda_manual(X, y): # 分割类别 X0 X[y0]; X1 X[y1] # 计算均值 mu0 np.mean(X0, axis0) mu1 np.mean(X1, axis0) # 计算散度矩阵 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结合进行分层降维