1. 这不是“调包跑通就完事”的 Iris 实验——它是一把解剖机器学习底层逻辑的手术刀你可能在 dozens 个教程里见过这个标题用 Jupyter 做 Iris 数据集的 Fisher 线性分类。但绝大多数内容止步于from sklearn.discriminant_analysis import LinearDiscriminantAnalysis一行代码然后画个散点图说“看准确率98%”。这就像教人修车只让你拧紧螺丝却不告诉你为什么这个扭矩值是25N·m、螺纹牙型如何影响预紧力、垫片失效会怎样引发缸体漏油。Fisher 线性判别FLD从来不是 sklearn 里一个黑箱函数它是统计学与几何直觉的完美结晶是理解“分类本质”的第一道分水岭。我带过十几期机器学习实训发现一个惊人现象能手推 Fisher 准则函数 $J(w) \frac{w^T S_B w}{w^T S_W w}$ 的学员后续理解 SVM 的最大间隔、PCA 的方差最大化、甚至神经网络的特征解耦速度比只背代码的快3倍以上。因为 FLD 强迫你直面三个核心问题类内离散度怎么量化类间分离度怎么建模投影方向为何必须同时压制噪声、拉大距离这篇文章不提供“一键运行”的 notebook 链接而是带你用纯 NumPy 从零手写 Fisher 投影矩阵用 Matplotlib 动态可视化决策边界如何随 $S_W$ 和 $S_B$ 的数值变化而呼吸起伏。你会看到当 setosa 类样本被误标为 versicolor 时$S_W$ 矩阵的 (0,0) 元素如何像血压计一样飙升当 virginica 的花瓣宽度整体偏移0.3cm$S_B$ 的迹trace又怎样悄然改变投影轴的倾斜角度。这不是一次数据集练习而是一次对“判别”二字的临床解剖。适合刚学完线性代数想验证理论、或已用过 sklearn 但总感觉“隔层纱”的实践者。接下来所有代码、公式、图表都服务于一个目标让你合上电脑后能在白纸上徒手推导出二维 Iris 的最优投影向量。2. 核心设计思路为什么放弃 sklearn坚持手写三重计算链2.1 拒绝黑箱从“调用函数”到“重建数学契约”很多初学者卡在第一步为什么 Fisher 判别叫“线性”却和线性回归的“线性”不是一回事线性回归的线性指输出 $y$ 是输入 $x$ 的线性组合而 Fisher 的线性指决策边界是超平面即分类规则 $w^T x b 0$。这个根本差异只有亲手构建 $w$ 才能刻进肌肉记忆。sklearn 的LinearDiscriminantAnalysis内部做了太多优化自动处理多类扩展LDA、正则化防止 $S_W$ 奇异、协方差矩阵收缩估计……这些对工程落地极重要但对理解原始 Fisher 思想却是干扰项。我们坚持手写就是要守住1936年 Fisher 论文里的纯粹契约仅针对两类问题仅用原始样本均值与协方差仅求解使类间散度与类内散度比值最大的单位向量 $w$。这个契约有三重不可妥协的约束约束一严格二类限定。Iris 数据集有三类但我们只取 setosa标签0和 versicolor标签1彻底规避多类 LDA 的复杂加权。这是为了聚焦核心矛盾——如何让两类在投影线上“拉开距离各自抱团”。约束二协方差矩阵显式计算。不用np.cov()的默认参数而是手动计算 $S_W \sum_{i0}^{1} \sum_{x \in C_i} (x - \mu_i)(x - \mu_i)^T$其中 $\mu_i$ 是第 $i$ 类样本均值。这个过程会暴露一个关键事实$S_W$ 不是所有样本的协方差而是每类内部协方差的加权和权重为各类样本数。当你看到 versicolor 类因样本数多50 vs 50其协方差对 $S_W$ 的贡献占比达50%就会明白为何类别不平衡会扭曲 Fisher 方向。约束三广义特征值问题显式求解。Fisher 解等价于求解 $S_B w \lambda S_W w$。我们不用scipy.linalg.eig(S_B, S_W)一步到位而是先计算 $S_W^{-1} S_B$再对其做标准特征值分解。这样做有两个硬收益第一当 $S_W$ 接近奇异时如某维度方差极小你能亲眼看到 $S_W^{-1}$ 的巨大元素如何放大数值误差从而理解为何实际工程中必须加正则项第二特征向量排序逻辑变得透明——最大特征值对应的向量就是使 $J(w)$ 最大的 $w$无需查文档确认排序规则。2.2 可视化驱动让抽象公式在二维平面上“活过来”Iris 数据集的妙处在于它有4个特征萼片长、萼片宽、花瓣长、花瓣宽但前两个主成分能解释97%的方差。这意味着我们可以将4D Fisher 投影压缩到2D平面可视化这是其他高维数据集无法提供的教学红利。我们的可视化不是静态截图而是构建三重联动视图左图原始4D空间的类中心连线。用向量箭头画出 $\mu_1 - \mu_0$这是 Fisher 方向的“直觉起点”——类中心连线方向天然有分离倾向但未考虑类内散布。中图$S_W$ 主成分椭圆。对 $S_W$ 做特征分解画出其特征向量构成的椭圆长轴最大特征值方向直观显示“类内噪声”的主导方向。你会发现花瓣长度维度的椭圆半径远大于萼片宽度说明该维度类内变异更强Fisher 必须在此方向施加更大压制。右图Fisher 投影线动态叠加。将计算出的 $w$ 向量归一化后以红色粗线画在中图上并在其上绘制两类投影后的直方图。当 $w$ 完美对齐类中心连线时直方图分离度最大一旦 $w$ 偏离直方图开始重叠——这个重叠面积就是 $J(w)$ 公式中分母 $w^T S_W w$ 的几何意义。这种三图联动把 $J(w) \frac{w^T (\mu_1 - \mu_0)(\mu_1 - \mu_0)^T w}{w^T S_W w}$ 这个分式转化成了可触摸的空间关系分子是 $w$ 在类中心差向量上的投影长度平方分母是 $w$ 在类内噪声椭圆上的“穿透深度”。你不再背诵公式而是用眼睛丈量它。2.3 工程鲁棒性手写代码如何比封装库更“抗造”有人质疑手写难道不怕 bug恰恰相反正是通过手写我们植入了三层防御机制让代码比直接调用 sklearn 更健壮第一层协方差矩阵条件数监控。在计算 $S_W^{-1}$ 前先用np.linalg.cond(S_W)检查条件数。Iris 数据中若错误地将所有4个特征投入计算$S_W$ 的条件数常达 $10^5$ 量级此时直接求逆会导致 $w$ 向量出现数量级为 $10^3$ 的异常值。我们的代码会触发警告“S_W condition number {cond:.2e} 1e4, consider feature scaling or regularization”并自动切换到伪逆np.linalg.pinv(S_W)。这个判断逻辑是 sklearn 默认关闭的调试开关。第二层投影向量归一化强制校验。Fisher 要求 $||w||1$但数值计算中累积误差会使 $w^T w 1.0000000002$ 或 $0.9999999997$。我们的代码在每次更新 $w$ 后执行w w / np.linalg.norm(w)并断言abs(np.linalg.norm(w) - 1) 1e-10。这个看似琐碎的步骤能避免后续计算中因范数漂移导致的决策边界偏移。第三层决策阈值自适应计算。sklearn 的predict()直接返回类别但 Fisher 的本质是找投影线上的最佳分割点。我们手写find_optimal_threshold(X_proj, y)函数遍历投影值范围对每个候选阈值 $t$计算误分类样本数取最小值对应的 $t$。这个 $t$ 就是 $(\mu_0 \mu_1)/2$ 在 $w$ 方向的投影它让两类在投影线上“公平切割”。当数据分布偏斜时如 versicolor 投影值明显右偏这个自适应阈值比固定中点更准——这正是工业场景中模型微调的核心技巧。这三重设计让手写代码不再是教学玩具而成为可嵌入真实 pipeline 的验证模块。某次为某高校实验室做图像特征筛选他们用我们的手写 Fisher 模块快速定位到某通道特征因光照变化导致 $S_W$ 崩溃及时剔除了该维度将下游分类器准确率从82%提升至89%。3. 核心细节解析从数据加载到 Fisher 向量的七步精解3.1 数据准备为什么必须手动切分而非用 train_test_splitIris 数据集虽小但数据切分方式直接影响 Fisher 结果的可解释性。sklearn 的train_test_split默认随机打乱这会破坏 Fisher 的统计假设——Fisher 要求训练集能代表总体的类内/类间分布。我们坚持手动切分# 加载原始数据使用 sklearn 内置但绝不调用其 LDA from sklearn.datasets import load_iris iris load_iris() X, y iris.data, iris.target # 严格按顺序取前100个样本setosa(0-49), versicolor(50-99) # 这确保两类样本在原始采集顺序中保持时空连续性 X_binary X[y 2] # 取前两类 y_binary y[y 2] # 手动切分前40个为训练集每类20个后10个为测试集每类5个 X_train np.vstack([X_binary[0:20], X_binary[50:70]]) y_train np.hstack([y_binary[0:20], y_binary[50:70]]) X_test np.vstack([X_binary[20:25], X_binary[70:75]]) y_test np.hstack([y_binary[20:25], y_binary[70:75]])这个切分法有三重深意意义一保留采集批次效应。Iris 样本按种类分批采集顺序切分让训练集包含同一批次的完整 setosa 和 versicolor避免随机切分将同一采集批次的样本拆散到训练/测试集导致 $S_W$ 低估类内变异。意义二控制样本量平衡。每类训练样本严格20个消除类别不平衡对 $S_W$ 权重的影响。若用train_test_split(test_size0.2)因随机性可能导致 setosa 训练样本18个、versicolor 22个此时 $S_W$ 中 versicolor 协方差权重为22/4055%偏离理论50%。意义三测试集设计为“压力测试”。测试集每类仅5个样本且取自训练集之后的连续样本模拟真实场景中“新批次数据”的分布偏移。当 Fisher 在此测试集上准确率骤降时你能立刻追溯是 $S_W$ 对新批次变异估计不足而非泛化能力问题。提示切分后务必检查np.unique(y_train, return_countsTrue)确认(array([0, 1]), array([20, 20]))。任何偏差都会污染 Fisher 方向的物理意义。3.2 类内散度矩阵 $S_W$手算协方差的“去中心化”陷阱$S_W$ 的定义是 $S_W \sum_{i0}^{1} \sum_{x \in C_i} (x - \mu_i)(x - \mu_i)^T$其中 $\mu_i$ 是第 $i$ 类样本均值。新手常犯的致命错误是用全体样本均值 $\mu$ 替代 $\mu_i$。我们用具体数字揭示后果设 setosa 类20个样本的萼片长度均值 $\mu_0^{(0)} 5.01$ cmversicolor 类20个样本的 $\mu_1^{(0)} 5.94$ cm。若错误使用全体均值 $\mu^{(0)} (5.015.94)/2 5.475$ cm 计算 $S_W$则 setosa 类每个样本的残差 $(x^{(0)} - 5.475)$ 平均被放大 $|5.475-5.01|0.465$ cmversicolor 类被缩小 $|5.94-5.475|0.465$ cm。这导致 $S_W$ 的 (0,0) 元素萼片长度方差被系统性高估Fisher 会过度压制萼片长度维度使 $w$ 向量偏向花瓣维度——而实际上萼片长度对区分这两类至关重要。正确做法是分步计算# 分别计算两类均值 mu0 np.mean(X_train[y_train 0], axis0) # shape (4,) mu1 np.mean(X_train[y_train 1], axis0) # shape (4,) # 初始化 S_W 为 4x4 零矩阵 S_W np.zeros((4, 4)) # 对 setosa 类y0累加 (x - mu0)(x - mu0)^T for x in X_train[y_train 0]: diff x - mu0 S_W np.outer(diff, diff) # outer product, not dot! # 对 versicolor 类y1累加 (x - mu1)(x - mu1)^T for x in X_train[y_train 1]: diff x - mu1 S_W np.outer(diff, diff)注意np.outer(diff, diff)与np.dot(diff.reshape(-1,1), diff.reshape(1,-1))等价但前者更直观体现“外积”本质。np.outer的结果是一个4x4矩阵其 (i,j) 元素是 $diff_i \times diff_j$这正是协方差矩阵的构成单元。若误用np.dot(diff, diff)得到的是标量残差平方和完全丢失维度间相关性信息。3.3 类间散度矩阵 $S_B$为什么是 $(\mu_1 - \mu_0)(\mu_1 - \mu_0)^T$ 而非 $\mu_1 \mu_1^T - \mu_0 \mu_0^T$$S_B$ 的标准形式是 $S_B (\mu_1 - \mu_0)(\mu_1 - \mu_0)^T$这是一个秩为1的矩阵。新手易混淆为 $\mu_1 \mu_1^T - \mu_0 \mu_0^T$后者在数学上不等于前者。我们用 Iris 的实际数值演示$\mu_0 [5.01, 3.42, 1.46, 0.24]$setosa 均值$\mu_1 [5.94, 2.77, 4.26, 1.33]$versicolor 均值$\mu_1 - \mu_0 [0.93, -0.65, 2.80, 1.09]$计算 $(\mu_1 - \mu_0)(\mu_1 - \mu_0)^T$ 的 (0,2) 元素萼片长与花瓣长交叉项 $$ (0.93) \times (2.80) 2.604 $$而 $\mu_1 \mu_1^T - \mu_0 \mu_0^T$ 的 (0,2) 元素是 $$ (5.94 \times 4.26) - (5.01 \times 1.46) 25.3044 - 7.3146 17.9898 $$二者相差近7倍根本原因在于$S_B$ 描述的是类中心之间的相对位移其几何意义是连接 $\mu_0$ 和 $\mu_1$ 的向量张成的方向。$(\mu_1 - \mu_0)(\mu_1 - \mu_0)^T$ 正是该向量的外积它唯一确定了类间分离的主方向。而 $\mu_1 \mu_1^T - \mu_0 \mu_0^T$ 混合了绝对位置信息当 $\mu_0$ 和 $\mu_1$ 同时平移时如所有样本萼片长度1cm$S_B$ 应不变但后者会剧烈变化违背 Fisher 的平移不变性假设。因此$S_B$ 的计算必须严格mu_diff mu1 - mu0 S_B np.outer(mu_diff, mu_diff) # 秩为1shape (4,4)3.4 Fisher 准则函数 $J(w)$梯度上升法的手动实现与收敛监控Fisher 的目标是最大化 $J(w) \frac{w^T S_B w}{w^T S_W w}$。虽然理论上可通过广义特征值求解但手动实现梯度上升能深化理解。$J(w)$ 对 $w$ 的梯度为 $$ \nabla_w J(w) \frac{2 S_B w}{w^T S_W w} - \frac{2 w^T S_B w}{(w^T S_W w)^2} S_W w $$ 我们用此梯度进行迭代更新# 初始化 w 为随机单位向量 w np.random.randn(4) w w / np.linalg.norm(w) # 存储 J(w) 历史值用于收敛判断 J_history [] for iter in range(100): numerator w.T S_B w denominator w.T S_W w J_w numerator / denominator J_history.append(J_w) # 计算梯度 grad (2 * S_B w) / denominator - (2 * numerator * S_W w) / (denominator ** 2) # 更新 w步长 eta0.1 w_new w 0.1 * grad # 强制单位化 w_new w_new / np.linalg.norm(w_new) # 检查收敛J(w) 变化小于 1e-6 if iter 0 and abs(J_history[-1] - J_history[-2]) 1e-6: print(fConverged at iteration {iter}) break w w_new这个手动梯度法的价值在于暴露数值不稳定性当 $S_W$ 接近奇异时denominator极小梯度grad会爆炸。此时你会看到J_history剧烈震荡从而理解为何必须监控条件数。验证理论最优性最终收敛的 $w$应与广义特征值解完全一致。我们用np.allclose(w_manual, w_eig, atol1e-5)断言失败则说明某步计算有误。理解学习率影响若将eta0.1改为eta1.0$w$ 会发散改为eta0.01收敛变慢。这让你亲身体会优化算法的敏感性。3.5 投影与决策从4D向量到1D判决的完整映射链得到最优 $w$ 后分类流程是投影对任意样本 $x$计算标量投影值 $x_{proj} w^T x$阈值判决若 $x_{proj} t$判为 versicolor否则为 setosa阈值 $t$ 的选择是关键。理论最优是两类投影均值的中点 $$ t \frac{1}{2} (w^T \mu_0 w^T \mu_1) w^T \left( \frac{\mu_0 \mu_1}{2} \right) $$ 但实际中我们采用更鲁棒的“最小误分类”法def find_optimal_threshold(X_proj, y): # X_proj: 1D array of projection values # y: binary labels [0,0,...,1,1...] thresholds np.linspace(X_proj.min(), X_proj.max(), 100) errors [] for t in thresholds: pred (X_proj t).astype(int) error np.mean(pred ! y) # 误分类率 errors.append(error) best_idx np.argmin(errors) return thresholds[best_idx], errors[best_idx] # 对训练集投影 X_train_proj X_train w t_opt, min_error find_optimal_threshold(X_train_proj, y_train)这个函数返回的t_opt是使训练集误分类最少的阈值。它自动适应数据分布若 versicolor 投影值整体右偏t_opt会右移避免将大量 versicolor 误判为 setosa。这比理论中点更贴近工程实践。3.6 可视化验证三图联动的 Matplotlib 实现细节三图联动可视化是理解 Fisher 的核心。关键代码如下import matplotlib.pyplot as plt from matplotlib.patches import Ellipse fig, axes plt.subplots(1, 3, figsize(15, 5)) # 左图原始4D空间的类中心连线降维到前2主成分 pca PCA(n_components2) X_pca pca.fit_transform(X_train) mu0_pca pca.transform(mu0.reshape(1,-1)).flatten() mu1_pca pca.transform(mu1.reshape(1,-1)).flatten() axes[0].scatter(X_pca[y_train0,0], X_pca[y_train0,1], cred, labelsetosa, alpha0.6) axes[0].scatter(X_pca[y_train1,0], X_pca[y_train1,1], cblue, labelversicolor, alpha0.6) axes[0].arrow(mu0_pca[0], mu0_pca[1], mu1_pca[0]-mu0_pca[0], mu1_pca[1]-mu0_pca[1], head_width0.05, length_includes_headTrue, colorblack, lw2) axes[0].set_title(Class Centers in PCA Space) # 中图S_W 主成分椭圆 Fisher 投影线 # 计算 S_W 在PCA空间的投影2x2矩阵 S_W_pca pca.components_ S_W pca.components_.T eigvals_W, eigvecs_W np.linalg.eigh(S_W_pca) # 绘制椭圆长轴最大特征值方向 ell Ellipse(xy(0,0), width2*np.sqrt(eigvals_W[1]), height2*np.sqrt(eigvals_W[0]), anglenp.degrees(np.arctan2(eigvecs_W[1,1], eigvecs_W[0,1])), edgecolorgreen, facecolornone, lw2, linestyle--) axes[1].add_patch(ell) # Fisher w 在PCA空间的投影 w_pca pca.transform(w.reshape(1,-1)).flatten() w_pca w_pca / np.linalg.norm(w_pca) # unit vector axes[1].arrow(0, 0, w_pca[0]*3, w_pca[1]*3, head_width0.1, colorred, lw3) axes[1].set_xlim(-3,3); axes[1].set_ylim(-3,3) axes[1].set_title(S_W Ellipse Fisher Direction) # 右图投影直方图 X_proj X_train w axes[2].hist(X_proj[y_train0], bins15, alpha0.6, labelsetosa, colorred) axes[2].hist(X_proj[y_train1], bins15, alpha0.6, labelversicolor, colorblue) axes[2].axvline(t_opt, colorblack, linestyle--, labelfOptimal t{t_opt:.3f}) axes[2].legend() axes[2].set_title(Projection Histograms)这段代码的精妙之处在于PCA 降维保真pca.components_是2x4矩阵pca.components_ S_W pca.components_.T将4D的 $S_W$ 精确投影到2D PCA 空间确保椭圆形状反映真实类内散布。椭圆角度计算np.arctan2(eigvecs_W[1,1], eigvecs_W[0,1])计算第二特征向量对应大特征值与x轴夹角这是椭圆长轴方向直观显示“噪声最强”的维度。Fisher 红线长度归一化w_pca / np.linalg.norm(w_pca)确保红线长度为1在图中清晰可见其方向而非被缩放失真。3.7 测试集评估超越准确率的三重诊断在测试集上我们不只报告准确率而是进行三重诊断X_test_proj X_test w y_pred (X_test_proj t_opt).astype(int) acc np.mean(y_pred y_test) # 1. 混淆矩阵 from sklearn.metrics import confusion_matrix cm confusion_matrix(y_test, y_pred) print(Confusion Matrix:) print(cm) # [[5,0],[0,5]] 表示完美 # 2. 投影分布对比训练vs测试 plt.figure(figsize(10,4)) plt.subplot(1,2,1) plt.hist(X_train_proj[y_train0], bins10, alpha0.5, labeltrain setosa) plt.hist(X_train_proj[y_train1], bins10, alpha0.5, labeltrain versicolor) plt.title(Training Set Projections) plt.subplot(1,2,2) plt.hist(X_test_proj[y_test0], bins10, alpha0.5, labeltest setosa) plt.hist(X_test_proj[y_test1], bins10, alpha0.5, labeltest versicolor) plt.title(Test Set Projections) plt.legend()混淆矩阵揭示错误模式。若出现[4,1]setosa 误判1个说明 Fisher 方向对 setosa 类的某维度压制不足。投影分布对比图若测试集直方图相比训练集整体右移表明新批次 versicolor 花瓣更长$S_W$ 未能捕捉此变异需在 $S_W$ 中加入正则项。4. 实操过程从零开始的完整代码与逐行注释4.1 环境准备与依赖声明我们使用最简依赖避免版本冲突# 创建干净环境推荐 conda create -n fisher-iris python3.9 conda activate fisher-iris pip install numpy matplotlib scikit-learn关键点不安装 tensorflow/pytorch。Fisher 是纯线性代数引入深度学习框架是认知污染。所有计算用numpy原生操作确保每一步都可控。4.2 完整可运行代码含详细注释# -*- coding: utf-8 -*- Fisher Linear Discriminant from Scratch on Iris Dataset Author: A practitioner with 12 years in ML deployment Date: 2023-10-15 import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import load_iris from sklearn.decomposition import PCA from sklearn.metrics import confusion_matrix # STEP 1: DATA LOADING MANUAL SPLIT print( STEP 1: Data Loading and Manual Split ) iris load_iris() X, y iris.data, iris.target # Extract only first two classes: setosa (0) and versicolor (1) # This is CRITICAL for pure Fisher (not multi-class LDA) X_binary X[y 2] y_binary y[y 2] # Manual split: first 20 of each class for training, next 5 for test # Preserves temporal order and class balance X_train np.vstack([X_binary[0:20], X_binary[50:70]]) y_train np.hstack([y_binary[0:20], y_binary[50:70]]) X_test np.vstack([X_binary[20:25], X_binary[70:75]]) y_test np.hstack([y_binary[20:25], y_binary[70:75]]) print(fTraining set: {X_train.shape[0]} samples ({np.sum(y_train0)} setosa, {np.sum(y_train1)} versicolor)) print(fTest set: {X_test.shape[0]} samples ({np.sum(y_test0)} setosa, {np.sum(y_test1)} versicolor)) # STEP 2: COMPUTE CLASS MEANS print(\n STEP 2: Computing Class Means ) mu0 np.mean(X_train[y_train 0], axis0) # shape (4,) mu1 np.mean(X_train[y_train 1], axis0) # shape (4,) print(fmu0 (setosa): {mu0.round(3)}) print(fmu1 (versicolor): {mu1.round(3)}) # STEP 3: COMPUTE WITHIN-CLASS SCATTER S_W print(\n STEP 3: Computing Within-Class Scatter S_W ) S_W np.zeros((4, 4)) # Accumulate for setosa for x in X_train[y_train 0]: diff x - mu0 S_W np.outer(diff, diff) # Outer product: (4,) x (4,) - (4,4) # Accumulate for versicolor for x in X_train[y_train 1]: diff x - mu1 S_W np.outer(diff, diff) # Check condition number - critical for numerical stability cond_S_W np.linalg.cond(S_W) print(fS_W condition number: {cond_S_W:.2e}) if cond_S_W 1e4: print(WARNING: S_W is ill-conditioned. Using pseudo-inverse.) S_W_inv np.linalg.pinv(S_W) else: S_W_inv np.linalg.inv(S_W) # STEP 4: COMPUTE BETWEEN-CLASS SCATTER S_B print(\n STEP 4: Computing Between-Class Scatter S_B ) mu_diff mu1 - mu0 S_B np.outer(mu_diff, mu_diff) # Rank-1 matrix, shape (4,4) print(fS_B shape: {S_B.shape}, rank: {np.linalg.matrix_rank(S_B)}) # STEP 5: SOLVE GENERALIZED EIGENVALUE PROBLEM print(\n STEP 5: Solving S_B w lambda S_W w ) # Compute S_W^{-1} S_B if cond_S_W 1e4: M S_W_inv S_B else: M np.linalg.inv(S_W) S_B # Standard eigenvalue decomposition eigvals, eigvecs np.linalg.eig(M) # Sort eigenvalues descending, get corresponding eigenvectors idx np.argsort(eigvals)[::-1] eigvals eigvals[idx] eigvecs eigvecs[:, idx] # Fishers optimal w is the eigenvector for largest eigenvalue w eigvecs[:, 0] # Force unit norm w w / np.linalg.norm(w) print(fOptimal w (unit vector): {w.round(4)}) print(fCorresponding eigenvalue: {eigvals[0]:.4f}) # STEP 6: FIND OPTIMAL THRESHOLD print(\n STEP 6: Finding Optimal Decision Threshold ) X_train_proj X_train w def find_optimal_threshold(X_proj, y): thresholds np.linspace(X_proj.min(), X_proj.max(), 100) errors [] for t in thresholds: pred (X_proj t).astype(int) errors.append(np.mean(pred ! y)) best_idx np.argmin(errors) return thresholds[best_idx], errors[best_idx] t_opt, min_error find_optimal_threshold(X_train_proj, y_train) print(fOptimal threshold t*: {t_opt:.4f} (min error: {min_error:.4f})) # STEP 7: EVALUATE ON TEST SET print(\n STEP 7: Test Set Evaluation ) X_test_proj X_test w y_pred (X_test_proj t_opt).astype(int) acc np.mean(y_pred