说实话第一次接触偏最小二乘回归PLS的时候我整个人是有点懵的。网上一搜全是“潜变量”“载荷矩阵”“交叉验证”这些词配套的公式密密麻麻看起来就像某个数学系博士生用来折磨人的毕业设计。但后来因为工作里要部署一个跨平台的小工具没法把整套Python环境搬过去只能手写核心推理逻辑我才被迫把这层窗户纸捅破。等真把NIPALS算法一步步拆开之后我只有一种感受PLS的核心逻辑拆开了就是矩阵来回折腾——这正是这篇文章想和你一起做的事。这篇文章不依赖任何机器学习的黑盒库从头到尾用numpy手写PLS的核心实现重点包括数据标准化的必要性、潜变量提取的迭代细节、残差矩阵更新方式以及回归系数从标准化空间还原到原始尺度时的换算逻辑。整个过程我会围绕PLS的矩阵运算展开把代码逐段讲明白、把每一步动作的理由说清楚。不管你是刚接触多元统计分析的学生还是需要在嵌入式环境或纯计算场景下手推算法的工程师这篇文章都适合你跟着敲一遍。1. 为什么“手动实现”这件事值得做1.1 调库和懂原理之间隔着一条认知断层我用sklearn用得很早PLSRegression(n_components5).fit(X, y)这种代码写过不下几十次训练、预测、画图一条龙轻松得像在点外卖。直到有次组会导师指着我的结果问了一句“你选的5个主成分依据是什么权重的迭代过程你能给我画出来吗”我当场愣住支支吾吾半天也没说出个所以然。那之后我才意识到一个模型如果只能跑通不能讲清那在真正需要它做决策的场合心里是发虚的。手动实现的意义不在于“不用调库显得很厉害”而在于把黑盒里的每一个矩阵变换摊开之后你对算法边界、异常行为和结果可靠性的判断会完全不同。比如你会在代码里亲自看到为什么PLS能处理共线性数据为什么成分数不能选太多为什么标准化这一步不能省这些问题如果不进到公式和迭代细节里是永远得不到直观答案的。1.2 PLS到底在解决什么问题PLS全称Partial Least Squares在化学计量学、光谱建模、经济预测等领域非常常见。它最典型的应用场景有三个特征自变量维度高、自变量之间强共线性、样本量可能小于变量数。这种情况下普通最小二乘的X^T X矩阵不可逆或条件数极差直接算回归系数会得到一堆离谱的值。PLS的做法是不在原始变量空间里硬回归而是在X中提取少数几个“潜变量”得分向量t并且让这些潜变量尽可能与y相关。换句话说它同时兼顾了“对X的解释”和“对y的预测”。有个形象的类比PCA有点像选一个各科成绩都好的“全才”而PLS选的是一个对“你关心的那门课”特别有促进作用的“偏才”。这个“偏才”就是潜变量它是原始变量按照权重w组合出来的新特征。1.3 手写实现的基本路线图整篇文章的实现路线很清晰一共四步数据标准化把X和y的均值、方差统一到同一量纲这是所有协方差计算的前提。迭代提取潜变量用NIPALS非线性迭代偏最小二乘算法一次提一个成分每次提取后更新残差矩阵。拟合回归系数在标准化空间里利用权重W、载荷P和y载荷Q计算回归系数B。还原尺度把标准化空间的系数映射回原始数据尺度算出截距得到最终预测公式。下面的内容我就按照这条路线一步步展开。2. 数据标准化PLS的第一步也是很多人忽略的定海神针2.1 标准化到底在做什么很多人在入门时觉得标准化就是“减均值除以标准差”好像很机械但这一步在PLS里不是例行公事而是直接参与算法核心计算的。实现上就是两个动作代码很轻量X_mean X.mean(axis0) X_std X.std(axis0, ddof0) # 总体标准差和sklearn默认行为一致 Xs (X - X_mean) / X_std y_mean y.mean() y_std y.std(ddof0) ys (y - y_mean) / y_std这里特别提醒一个细节ddof参数的取值。ddof0算的是总体标准差分母是nddof1算的是样本标准差分母是n-1。sklearn.preprocessing.StandardScaler默认使用ddof0sklearn的PLSRegression在scaleTrue时也按这种方式标准化。如果你想和现有工具链的结果对得上这里必须统一否则预测值可能差一个小数点后几位的数值。2.2 为什么PLS对标准化这么敏感回到算法内部。PLS在计算权重w时要算X^T u这本质上是X每一列与u的协方差在中心化后。如果X的第1列是范围在0.1到0.2的微小量第2列是范围在1000到10000的大数量那么X^T u的结果必然被第2列主导第1列的信息几乎被淹没。这会导致提取出的潜变量严重偏向量纲大的变量模型给出的“重要性”完全失真。你可以这样理解不标准化相当于让一个体重用“斤”计、一个身高用“纳米”计的人去参加体质测试还没开始跑胜负已经被单位决定了一大半。PLS最终解出的权重w是变量间的相对比较只有把所有变量都放在以“标准差为单位”的同一把尺子下这个比较才公平。2.3 标准化里最容易被忽略的坑参数泄漏标准化看似简单但实际使用中一个非常隐蔽的坑就是“参数泄漏”。具体表现是有些人先对全部数据做标准化然后再划分训练集和测试集。这个做法是错的——因为测试集的均值、标准差参与了对训练数据的变换相当于把测试集的信息倒灌给了训练阶段。举个例子假设训练集和测试集的整体均值是100但训练集本身的均值是80。如果先整体标准化再划分训练集用的就是100这个均值而不是80。这会让训练过程“偷看”到测试集的分布交叉验证时的误差估计会过于乐观模型上线后表现往往明显下降。正确的做法是# 先划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42) # 在训练集上计算参数 X_mean X_train.mean(axis0) X_std X_train.std(axis0, ddof0) y_mean y_train.mean() y_std y_train.std(ddof0) # 用训练集参数变换训练集和测试集 X_train_s (X_train - X_mean) / X_std X_test_s (X_test - X_mean) / X_std y_train_s (y_train - y_mean) / y_std后面实现中保存这些标准化参数还有一个作用在线预测时必须用同一套参数训练得到的均值和标准差来变换新数据而不是重新算新数据的均值。这一点在部署时尤其关键我见过不止一个同学在部署环节里对单个新样本做标准化结果预测结果完全跑偏。3. NIPALS算法PLS核心迭代逻辑的逐步拆解3.1 先理清权重、得分、载荷这些术语NIPALS算法是PLS最经典的求解方式。在进入代码之前先把几个核心矩阵向量的角色弄清楚后面才不会绕晕。X标准化后的自变量矩阵形状为(n_samples, n_features)。y标准化后的响应变量向量这里为了方便讲解先做单输出PLS1多输出的PLS2后面单独说形状为(n_samples,)。w权重X的列组合系数长度为p且被归一化为单位长度。它表示潜变量是由哪些原始变量如何组合出来的。t得分样本在潜变量方向上的投影值长度为n计算公式是t Xw。它是“新特征”的具体取值。p载荷把得分t映射回X空间的系数p X^T t / (t^T t)。它表示这个潜变量对X的原始变量“解释”了多少方差。qy载荷y方向上的载荷标量q y^T t / (t^T t)可以理解为t与y的最佳回归斜率。整个NIPALS循环本质上就是在反复求解“X的权重”和“y的载荷”直到得到收敛的一组潜变量。3.2 单个潜变量提取的完整循环假设现在只提取一个成分具体步骤如下。这个循环也是整个PLS迭代的核心初始化u为y在PLS1中u直接取残差y向量即可在PLS2中一般取Y矩阵中方差最大的一列。计算X的权重w X^T u / ||X^T u||。这一步算的是X每一列与u的相关性归一化保证方向一致。计算X的得分t Xw。这是所有样本在那个新方向上的投影。计算y的载荷q (y^T t) / (t^T t)。这是t对y做最小二乘回归得到的系数。更新uu yq。检查收敛如果新的u与旧的u变化量小于阈值就认为迭代收敛否则回到第2步重新算。这段迭代逻辑我直接用代码写出来u y_res.copy() for _ in range(max_iter): w X_res.T u w w / np.linalg.norm(w) t X_res w q (y_res t) / (t t) u_new y_res * q if np.linalg.norm(u_new - u) tol: u u_new break u u_new有朋友第一次看这段代码会问既然u最终只是y的q倍方向并不会发生根本变化那这个内层迭代到底在收敛什么实际上在PLS1里因为u的方向从一开始就由y决定这个循环通常一两轮就稳定了。但它并非无用——它会校正t和q的数值尺度而且当扩展到PLS2的多输出场景时这个循环才真正体现“交替迭代”的意义。3.3 多成分提取的关键残差矩阵的更新deflationPLS和PCA有点类似都不会只提取一个成分而是提取多个成分来逐步逼近X和y的关系。但这里有一个关键问题如果我们对原始X反复执行同一个提取过程得到的永远是第一个权重方向。怎么得到第二个、第三个不同的方向呢答案就是deflation也就是“矩阵收缩”。每次提取完一个成分后要把它对X和y的贡献从当前残差矩阵中减掉p (X_res.T t) / (t t) X_res X_res - np.outer(t, p) y_res y_res - t * q这里np.outer(t, p)计算的是t与p的外积代表这个潜变量对X的“重构贡献”。减掉之后X_res中剩下的就是第一个成分解释不了的部分。下一次迭代就在这个残差矩阵上进行提取出与第一个方向正交的第二个方向。这个过程就像剥洋葱一层层把信息结构剥出来。每次deflation后残差矩阵的方差逐渐下降当剩余成分基本是噪声时继续提取就没有意义了。3.4 收敛判据与迭代细节动手实现NIPALS时有几个数值细节直接决定代码能不能稳定跑出结果。收敛阈值我习惯设为1e-8或1e-10太大会提前收敛导致精度损失太小会白白增加迭代次数。最大迭代次数必须加一个上限比如max_iter 500防止极端情况下死循环。维度问题计算X_res.T u时要确保u是二维列向量或一维数组时numpy广播正确。推荐统一用一维数组加运算简单直观。分母保护q (y_res t) / (t t)可能遇到极小的t t这种情况通常发生在数据被严重噪声主导时。如果遇到数值异常可以给分母加一个很小的eps值比如1e-12。我在实际实现中还发现残差矩阵在数学上应该保持对称或半正定特性但浮点累积误差会有轻微破坏。如果后续要做更精细的诊断可以用对称化处理或定期重构不过对于一般回归预测这点误差可以忽略。4. 用numpy从零实现PLS完整代码与逐段讲解4.1 函数整体设计现在把前面的逻辑串起来实现一个完整的PLS1函数。函数输入是训练数据X、y需要指定的成分数n_components以及可选的收敛阈值和最大迭代次数。函数返回核心结果原始尺度下的回归系数B、截距intercept另外把标准化空间的W、T、P、Q也一并返回方便后续做诊断和变量重要性分析。import numpy as np def pls_manual(X, y, n_components, tol1e-8, max_iter500): # 1. 保存原始尺度参数并标准化 X_mean X.mean(axis0) X_std X.std(axis0, ddof0) y_mean y.mean() y_std y.std(ddof0) Xs (X - X_mean) / X_std ys (y - y_mean) / y_std n, p Xs.shape W np.zeros((p, n_components)) # 权重矩阵 T np.zeros((n, n_components)) # 得分矩阵 P np.zeros((p, n_components)) # 载荷矩阵 Q np.zeros(n_components) # y载荷 X_res Xs.copy() y_res ys.copy() # 2. 主循环逐个提取潜变量 for k in range(n_components): u y_res.copy() for _ in range(max_iter): w X_res.T u w w / np.linalg.norm(w) t X_res w q (y_res t) / (t t) u_new y_res * q if np.linalg.norm(u_new - u) tol: u u_new break u u_new p (X_res.T t) / (t t) # 保存当前成分 W[:, k] w T[:, k] t P[:, k] p Q[k] q # deflation X_res X_res - np.outer(t, p) y_res y_res - t * q # 3. 标准化空间中的回归系数 # B_std W inv(P.T W) Q WTP_inv np.linalg.pinv(P.T W) B_std W WTP_inv Q # 4. 还原到原始尺度 B B_std * (y_std / X_std) intercept y_mean - X_mean B return B, intercept, W, T, P, Q4.2 主循环里每个矩阵的维度变化新手最容易卡住的地方是维度不匹配。我把每步的维度写一下方便对照排查X_res(n, p)u(n,)w(p,)t(n,)q标量p(p,)np.outer(t, p)(n, p)正好和X_res形状一致整个循环里没有任何一处需要出现(n, 1)这种多余的维度如果你在自己的实现中遇到广播错误先检查是否多写了一个维度的reshape。这是NIPALS实现里最常见的低级错误报错信息还特别迷惑。4.3 为什么回归系数要还原到原始尺度在标准化空间里我们拟合的关系是y_std ≈ X_std B_std注意这里的B_std是在“标准化后的X和y”下的系数它默认所有变量的单位都是一倍的“标准差”。但实际预测时用户拿到的是原始量纲的数据不能直接往这个公式里套。把公式展开一下就很清楚了原始关系y_hat X B intercept标准化空间的预测(y_hat - y_mean) / y_std ((X - X_mean) / X_std) B_std整理得y_hat y_mean y_std * (((X - X_mean) / X_std) B_std)进一步展开y_hat (y_std / X_std) * B_std X (y_mean - X_mean ((y_std / X_std) * B_std))所以B B_std * (y_std / X_std)intercept y_mean - X_mean B。这一段推导建议自己动手在纸上写一遍很多部署问题都出在对这步换算的理解上。我在代码里使用了np.linalg.pinv伪逆而不是np.linalg.inv这是因为P.T W在成分数较多时可能接近奇异直接用逆容易得到异常大的系数。伪逆的数值稳定性更好代价是计算时间略高但对一般的PLS场景来说完全可接受。5. 用模拟数据和sklearn对比验证5.1 构建一份带共线性的实验数据为了验证手写实现是否正确我准备了一份带共线性的模拟数据。场景是这样的有10个自变量但真正影响y的只有其中3个剩下几个变量与它们高度相关形成典型的共线性结构。np.random.seed(42) n, p 300, 10 X np.random.randn(n, p) # 真实系数只有第0、3、7个变量影响y true_beta np.zeros(p) true_beta[[0, 3, 7]] [2.0, -1.5, 0.8] y X true_beta np.random.randn(n) * 0.3 # 人为制造共线性第4、5列与第0、3列强相关 X[:, 4] X[:, 0] * 0.9 np.random.randn(n) * 0.1 X[:, 5] -X[:, 3] * 0.8 np.random.randn(n) * 0.1这种数据放到普通线性回归里X^T X的条件数会非常大回归系数极不稳定。但在PLS里由于我们只提取少数几个潜变量就可以绕开共线性问题。5.2 完整对比预测值、RMSE、回归系数然后用训练集拟合在测试集上比较手写实现和sklearn的结果。from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state7) # 手写实现 B_manual, intercept_manual, W, T, P, Q pls_manual(X_train, y_train, n_components3) y_pred_manual X_test B_manual intercept_manual rmse_manual np.sqrt(mean_squared_error(y_test, y_pred_manual)) # sklearn实现 pls_sk PLSRegression(n_components3, scaleTrue) pls_sk.fit(X_train, y_train) y_pred_sk pls_sk.predict(X_test).flatten() rmse_sk np.sqrt(mean_squared_error(y_test, y_pred_sk)) print(fRMSE 手写实现: {rmse_manual:.6f}) print(fRMSE sklearn: {rmse_sk:.6f}) print(f回归系数差值最大绝对值: {np.max(np.abs(B_manual - pls_sk.coef_.flatten())):.2e})在我本机跑出来的结果中两个RMSE几乎一致回归系数差值的最大绝对值在1e-8到1e-7量级。这说明手写实现和sklearn在算法层面是等价的。5.3 为什么会有细微差异虽然结果一致但如果你把系数和预测值逐位对比还是能发现极小差异。这些差异不是谁错了而是来自几个方面sklearn内部对收敛阈值和最大迭代次数的默认设置与手写代码不完全一样sklearn的PLSRegression在求解过程中可能启用了SVD或其他数值优化路径而不是纯NIPALS的幂迭代浮点运算顺序不同会带来微小的舍入误差。这些差异在1e-7量级对实际预测结果的影响可以忽略。但有个用途如果在自己的项目里需要确认另一个实现是否可靠可以用这种“和sklearn对比差值是否在1e-6以内”的办法做回归测试非常有效。6. 实操中那些文档不会告诉你的经验6.1 成分数量选择不能只看累计解释方差很多人用PLS时会参考PCA的习惯根据X的累计解释方差来选成分数。但PLS是监督方法潜变量不仅要解释X还要和y相关。只看X方差容易选出一些对y毫无解释能力的“大方差噪声成分”。我自己的做法是做交叉验证对不同成分数分别计算验证集RMSE画出一条RMSE随成分数变化的曲线。通常情况下RMSE会先快速下降然后进入平台期最后因过拟合略微上升。选那个进入平台期对应的最小成分数——如果再增加成分RMSE下降不显著就没必要多算。需要特别注意的是交叉验证的每一折内部都要重新计算标准化的均值、标准差不能直接用全量数据的标准化结果。否则又会出现前面说的参数泄漏。6.2 收敛失败时先检查标准化和初值如果NIPALS迭代一直不收敛不要急着调大max_iter先按下面顺序排查数据是否包含NaNNaN在矩阵运算里会像病毒一样传染导致整个矩阵失效。标准化是否成功如果某个变量标准差为0比如工程上常见的常量列会出现除以0的问题。要么删除该列要么给它加一个极小扰动。初值u的选择是否合适在PLS2中u一般取Y矩阵中方差最大的一列。若所有y都接近常量迭代方向就会随机漂移。我调试时还有一个习惯在每次内层迭代打印一下np.linalg.norm(u_new - u)如果这个值振荡不降或者跳到NaN往往不是算法公式错而是前面的标准化环节埋了雷。6.3 从PLS1扩展到PLS2的改动现实中常常需要同时预测多个响应变量比如用光谱数据同时预测水分含量、蛋白质含量和脂肪含量。这就是PLS2多输出PLS。代码需要改动的地方其实很小u初值从Y矩阵中选方差最大的一列q从标量变成向量q (Y_res.T t) / (t t)形状为(m,)m是输出维度u_new Y_res q / (q q)deflation时Y_res Y_res - np.outer(t, q)最终回归系数B_std用(p, n_components) (n_components, n_components) (n_components, m)的维度计算还原尺度时要注意对每个输出列分别乘对应的y_std。只要理解了单输出版本的循环逻辑扩展到多输出其实只是把几个标量换成向量矩阵没有本质难度。但有几个易错点值得留意Q矩阵的维度、Y标准化的方式、还原时intercept需要按列计算。我建议在扩展之前先把单输出版本的推导过程完整写一遍。6.4 数值稳定性和性能提升思路手写实现最大的好处就是能针对具体计算环境调优。如果矩阵规模较大可以把np.linalg.pinv(P.T W)替换成np.linalg.solve(P.T W, np.eye(n_components))这样更快。如果数据维度极高比如光谱数据上万列NIPALS的每次迭代都会做一次大型矩阵乘法计算压力不小。此时可以考虑Kernel PLS用核矩阵代替原始X的运算但这是另一个主题了暂时不展开。对于普通规模的数据几千行、几百列纯numpy实现已经足够快。因为numpy底层的矩阵乘法是高度优化的BLAS实现手写循环最忌讳在Python层逐元素遍历矩阵只要用向量化矩阵运算性能就完全在可接受范围内。最后说点实在的。我在写完这套实现之后最大的收获不是“会手写PLS了”而是以后再遇到奇异矩阵、共线性诊断、变量重要性这类问题时脑子里会自动浮现出矩阵迭代的过程而不是一个黑盒的输入输出。你在实现中如果真的遇到问题先别急着怀疑算法公式——回到数据标准化那里看一眼十次里有八次是那里出了问题。
用numpy手写偏最小二乘回归：NIPALS算法与矩阵运算详解
说实话第一次接触偏最小二乘回归PLS的时候我整个人是有点懵的。网上一搜全是“潜变量”“载荷矩阵”“交叉验证”这些词配套的公式密密麻麻看起来就像某个数学系博士生用来折磨人的毕业设计。但后来因为工作里要部署一个跨平台的小工具没法把整套Python环境搬过去只能手写核心推理逻辑我才被迫把这层窗户纸捅破。等真把NIPALS算法一步步拆开之后我只有一种感受PLS的核心逻辑拆开了就是矩阵来回折腾——这正是这篇文章想和你一起做的事。这篇文章不依赖任何机器学习的黑盒库从头到尾用numpy手写PLS的核心实现重点包括数据标准化的必要性、潜变量提取的迭代细节、残差矩阵更新方式以及回归系数从标准化空间还原到原始尺度时的换算逻辑。整个过程我会围绕PLS的矩阵运算展开把代码逐段讲明白、把每一步动作的理由说清楚。不管你是刚接触多元统计分析的学生还是需要在嵌入式环境或纯计算场景下手推算法的工程师这篇文章都适合你跟着敲一遍。1. 为什么“手动实现”这件事值得做1.1 调库和懂原理之间隔着一条认知断层我用sklearn用得很早PLSRegression(n_components5).fit(X, y)这种代码写过不下几十次训练、预测、画图一条龙轻松得像在点外卖。直到有次组会导师指着我的结果问了一句“你选的5个主成分依据是什么权重的迭代过程你能给我画出来吗”我当场愣住支支吾吾半天也没说出个所以然。那之后我才意识到一个模型如果只能跑通不能讲清那在真正需要它做决策的场合心里是发虚的。手动实现的意义不在于“不用调库显得很厉害”而在于把黑盒里的每一个矩阵变换摊开之后你对算法边界、异常行为和结果可靠性的判断会完全不同。比如你会在代码里亲自看到为什么PLS能处理共线性数据为什么成分数不能选太多为什么标准化这一步不能省这些问题如果不进到公式和迭代细节里是永远得不到直观答案的。1.2 PLS到底在解决什么问题PLS全称Partial Least Squares在化学计量学、光谱建模、经济预测等领域非常常见。它最典型的应用场景有三个特征自变量维度高、自变量之间强共线性、样本量可能小于变量数。这种情况下普通最小二乘的X^T X矩阵不可逆或条件数极差直接算回归系数会得到一堆离谱的值。PLS的做法是不在原始变量空间里硬回归而是在X中提取少数几个“潜变量”得分向量t并且让这些潜变量尽可能与y相关。换句话说它同时兼顾了“对X的解释”和“对y的预测”。有个形象的类比PCA有点像选一个各科成绩都好的“全才”而PLS选的是一个对“你关心的那门课”特别有促进作用的“偏才”。这个“偏才”就是潜变量它是原始变量按照权重w组合出来的新特征。1.3 手写实现的基本路线图整篇文章的实现路线很清晰一共四步数据标准化把X和y的均值、方差统一到同一量纲这是所有协方差计算的前提。迭代提取潜变量用NIPALS非线性迭代偏最小二乘算法一次提一个成分每次提取后更新残差矩阵。拟合回归系数在标准化空间里利用权重W、载荷P和y载荷Q计算回归系数B。还原尺度把标准化空间的系数映射回原始数据尺度算出截距得到最终预测公式。下面的内容我就按照这条路线一步步展开。2. 数据标准化PLS的第一步也是很多人忽略的定海神针2.1 标准化到底在做什么很多人在入门时觉得标准化就是“减均值除以标准差”好像很机械但这一步在PLS里不是例行公事而是直接参与算法核心计算的。实现上就是两个动作代码很轻量X_mean X.mean(axis0) X_std X.std(axis0, ddof0) # 总体标准差和sklearn默认行为一致 Xs (X - X_mean) / X_std y_mean y.mean() y_std y.std(ddof0) ys (y - y_mean) / y_std这里特别提醒一个细节ddof参数的取值。ddof0算的是总体标准差分母是nddof1算的是样本标准差分母是n-1。sklearn.preprocessing.StandardScaler默认使用ddof0sklearn的PLSRegression在scaleTrue时也按这种方式标准化。如果你想和现有工具链的结果对得上这里必须统一否则预测值可能差一个小数点后几位的数值。2.2 为什么PLS对标准化这么敏感回到算法内部。PLS在计算权重w时要算X^T u这本质上是X每一列与u的协方差在中心化后。如果X的第1列是范围在0.1到0.2的微小量第2列是范围在1000到10000的大数量那么X^T u的结果必然被第2列主导第1列的信息几乎被淹没。这会导致提取出的潜变量严重偏向量纲大的变量模型给出的“重要性”完全失真。你可以这样理解不标准化相当于让一个体重用“斤”计、一个身高用“纳米”计的人去参加体质测试还没开始跑胜负已经被单位决定了一大半。PLS最终解出的权重w是变量间的相对比较只有把所有变量都放在以“标准差为单位”的同一把尺子下这个比较才公平。2.3 标准化里最容易被忽略的坑参数泄漏标准化看似简单但实际使用中一个非常隐蔽的坑就是“参数泄漏”。具体表现是有些人先对全部数据做标准化然后再划分训练集和测试集。这个做法是错的——因为测试集的均值、标准差参与了对训练数据的变换相当于把测试集的信息倒灌给了训练阶段。举个例子假设训练集和测试集的整体均值是100但训练集本身的均值是80。如果先整体标准化再划分训练集用的就是100这个均值而不是80。这会让训练过程“偷看”到测试集的分布交叉验证时的误差估计会过于乐观模型上线后表现往往明显下降。正确的做法是# 先划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42) # 在训练集上计算参数 X_mean X_train.mean(axis0) X_std X_train.std(axis0, ddof0) y_mean y_train.mean() y_std y_train.std(ddof0) # 用训练集参数变换训练集和测试集 X_train_s (X_train - X_mean) / X_std X_test_s (X_test - X_mean) / X_std y_train_s (y_train - y_mean) / y_std后面实现中保存这些标准化参数还有一个作用在线预测时必须用同一套参数训练得到的均值和标准差来变换新数据而不是重新算新数据的均值。这一点在部署时尤其关键我见过不止一个同学在部署环节里对单个新样本做标准化结果预测结果完全跑偏。3. NIPALS算法PLS核心迭代逻辑的逐步拆解3.1 先理清权重、得分、载荷这些术语NIPALS算法是PLS最经典的求解方式。在进入代码之前先把几个核心矩阵向量的角色弄清楚后面才不会绕晕。X标准化后的自变量矩阵形状为(n_samples, n_features)。y标准化后的响应变量向量这里为了方便讲解先做单输出PLS1多输出的PLS2后面单独说形状为(n_samples,)。w权重X的列组合系数长度为p且被归一化为单位长度。它表示潜变量是由哪些原始变量如何组合出来的。t得分样本在潜变量方向上的投影值长度为n计算公式是t Xw。它是“新特征”的具体取值。p载荷把得分t映射回X空间的系数p X^T t / (t^T t)。它表示这个潜变量对X的原始变量“解释”了多少方差。qy载荷y方向上的载荷标量q y^T t / (t^T t)可以理解为t与y的最佳回归斜率。整个NIPALS循环本质上就是在反复求解“X的权重”和“y的载荷”直到得到收敛的一组潜变量。3.2 单个潜变量提取的完整循环假设现在只提取一个成分具体步骤如下。这个循环也是整个PLS迭代的核心初始化u为y在PLS1中u直接取残差y向量即可在PLS2中一般取Y矩阵中方差最大的一列。计算X的权重w X^T u / ||X^T u||。这一步算的是X每一列与u的相关性归一化保证方向一致。计算X的得分t Xw。这是所有样本在那个新方向上的投影。计算y的载荷q (y^T t) / (t^T t)。这是t对y做最小二乘回归得到的系数。更新uu yq。检查收敛如果新的u与旧的u变化量小于阈值就认为迭代收敛否则回到第2步重新算。这段迭代逻辑我直接用代码写出来u y_res.copy() for _ in range(max_iter): w X_res.T u w w / np.linalg.norm(w) t X_res w q (y_res t) / (t t) u_new y_res * q if np.linalg.norm(u_new - u) tol: u u_new break u u_new有朋友第一次看这段代码会问既然u最终只是y的q倍方向并不会发生根本变化那这个内层迭代到底在收敛什么实际上在PLS1里因为u的方向从一开始就由y决定这个循环通常一两轮就稳定了。但它并非无用——它会校正t和q的数值尺度而且当扩展到PLS2的多输出场景时这个循环才真正体现“交替迭代”的意义。3.3 多成分提取的关键残差矩阵的更新deflationPLS和PCA有点类似都不会只提取一个成分而是提取多个成分来逐步逼近X和y的关系。但这里有一个关键问题如果我们对原始X反复执行同一个提取过程得到的永远是第一个权重方向。怎么得到第二个、第三个不同的方向呢答案就是deflation也就是“矩阵收缩”。每次提取完一个成分后要把它对X和y的贡献从当前残差矩阵中减掉p (X_res.T t) / (t t) X_res X_res - np.outer(t, p) y_res y_res - t * q这里np.outer(t, p)计算的是t与p的外积代表这个潜变量对X的“重构贡献”。减掉之后X_res中剩下的就是第一个成分解释不了的部分。下一次迭代就在这个残差矩阵上进行提取出与第一个方向正交的第二个方向。这个过程就像剥洋葱一层层把信息结构剥出来。每次deflation后残差矩阵的方差逐渐下降当剩余成分基本是噪声时继续提取就没有意义了。3.4 收敛判据与迭代细节动手实现NIPALS时有几个数值细节直接决定代码能不能稳定跑出结果。收敛阈值我习惯设为1e-8或1e-10太大会提前收敛导致精度损失太小会白白增加迭代次数。最大迭代次数必须加一个上限比如max_iter 500防止极端情况下死循环。维度问题计算X_res.T u时要确保u是二维列向量或一维数组时numpy广播正确。推荐统一用一维数组加运算简单直观。分母保护q (y_res t) / (t t)可能遇到极小的t t这种情况通常发生在数据被严重噪声主导时。如果遇到数值异常可以给分母加一个很小的eps值比如1e-12。我在实际实现中还发现残差矩阵在数学上应该保持对称或半正定特性但浮点累积误差会有轻微破坏。如果后续要做更精细的诊断可以用对称化处理或定期重构不过对于一般回归预测这点误差可以忽略。4. 用numpy从零实现PLS完整代码与逐段讲解4.1 函数整体设计现在把前面的逻辑串起来实现一个完整的PLS1函数。函数输入是训练数据X、y需要指定的成分数n_components以及可选的收敛阈值和最大迭代次数。函数返回核心结果原始尺度下的回归系数B、截距intercept另外把标准化空间的W、T、P、Q也一并返回方便后续做诊断和变量重要性分析。import numpy as np def pls_manual(X, y, n_components, tol1e-8, max_iter500): # 1. 保存原始尺度参数并标准化 X_mean X.mean(axis0) X_std X.std(axis0, ddof0) y_mean y.mean() y_std y.std(ddof0) Xs (X - X_mean) / X_std ys (y - y_mean) / y_std n, p Xs.shape W np.zeros((p, n_components)) # 权重矩阵 T np.zeros((n, n_components)) # 得分矩阵 P np.zeros((p, n_components)) # 载荷矩阵 Q np.zeros(n_components) # y载荷 X_res Xs.copy() y_res ys.copy() # 2. 主循环逐个提取潜变量 for k in range(n_components): u y_res.copy() for _ in range(max_iter): w X_res.T u w w / np.linalg.norm(w) t X_res w q (y_res t) / (t t) u_new y_res * q if np.linalg.norm(u_new - u) tol: u u_new break u u_new p (X_res.T t) / (t t) # 保存当前成分 W[:, k] w T[:, k] t P[:, k] p Q[k] q # deflation X_res X_res - np.outer(t, p) y_res y_res - t * q # 3. 标准化空间中的回归系数 # B_std W inv(P.T W) Q WTP_inv np.linalg.pinv(P.T W) B_std W WTP_inv Q # 4. 还原到原始尺度 B B_std * (y_std / X_std) intercept y_mean - X_mean B return B, intercept, W, T, P, Q4.2 主循环里每个矩阵的维度变化新手最容易卡住的地方是维度不匹配。我把每步的维度写一下方便对照排查X_res(n, p)u(n,)w(p,)t(n,)q标量p(p,)np.outer(t, p)(n, p)正好和X_res形状一致整个循环里没有任何一处需要出现(n, 1)这种多余的维度如果你在自己的实现中遇到广播错误先检查是否多写了一个维度的reshape。这是NIPALS实现里最常见的低级错误报错信息还特别迷惑。4.3 为什么回归系数要还原到原始尺度在标准化空间里我们拟合的关系是y_std ≈ X_std B_std注意这里的B_std是在“标准化后的X和y”下的系数它默认所有变量的单位都是一倍的“标准差”。但实际预测时用户拿到的是原始量纲的数据不能直接往这个公式里套。把公式展开一下就很清楚了原始关系y_hat X B intercept标准化空间的预测(y_hat - y_mean) / y_std ((X - X_mean) / X_std) B_std整理得y_hat y_mean y_std * (((X - X_mean) / X_std) B_std)进一步展开y_hat (y_std / X_std) * B_std X (y_mean - X_mean ((y_std / X_std) * B_std))所以B B_std * (y_std / X_std)intercept y_mean - X_mean B。这一段推导建议自己动手在纸上写一遍很多部署问题都出在对这步换算的理解上。我在代码里使用了np.linalg.pinv伪逆而不是np.linalg.inv这是因为P.T W在成分数较多时可能接近奇异直接用逆容易得到异常大的系数。伪逆的数值稳定性更好代价是计算时间略高但对一般的PLS场景来说完全可接受。5. 用模拟数据和sklearn对比验证5.1 构建一份带共线性的实验数据为了验证手写实现是否正确我准备了一份带共线性的模拟数据。场景是这样的有10个自变量但真正影响y的只有其中3个剩下几个变量与它们高度相关形成典型的共线性结构。np.random.seed(42) n, p 300, 10 X np.random.randn(n, p) # 真实系数只有第0、3、7个变量影响y true_beta np.zeros(p) true_beta[[0, 3, 7]] [2.0, -1.5, 0.8] y X true_beta np.random.randn(n) * 0.3 # 人为制造共线性第4、5列与第0、3列强相关 X[:, 4] X[:, 0] * 0.9 np.random.randn(n) * 0.1 X[:, 5] -X[:, 3] * 0.8 np.random.randn(n) * 0.1这种数据放到普通线性回归里X^T X的条件数会非常大回归系数极不稳定。但在PLS里由于我们只提取少数几个潜变量就可以绕开共线性问题。5.2 完整对比预测值、RMSE、回归系数然后用训练集拟合在测试集上比较手写实现和sklearn的结果。from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state7) # 手写实现 B_manual, intercept_manual, W, T, P, Q pls_manual(X_train, y_train, n_components3) y_pred_manual X_test B_manual intercept_manual rmse_manual np.sqrt(mean_squared_error(y_test, y_pred_manual)) # sklearn实现 pls_sk PLSRegression(n_components3, scaleTrue) pls_sk.fit(X_train, y_train) y_pred_sk pls_sk.predict(X_test).flatten() rmse_sk np.sqrt(mean_squared_error(y_test, y_pred_sk)) print(fRMSE 手写实现: {rmse_manual:.6f}) print(fRMSE sklearn: {rmse_sk:.6f}) print(f回归系数差值最大绝对值: {np.max(np.abs(B_manual - pls_sk.coef_.flatten())):.2e})在我本机跑出来的结果中两个RMSE几乎一致回归系数差值的最大绝对值在1e-8到1e-7量级。这说明手写实现和sklearn在算法层面是等价的。5.3 为什么会有细微差异虽然结果一致但如果你把系数和预测值逐位对比还是能发现极小差异。这些差异不是谁错了而是来自几个方面sklearn内部对收敛阈值和最大迭代次数的默认设置与手写代码不完全一样sklearn的PLSRegression在求解过程中可能启用了SVD或其他数值优化路径而不是纯NIPALS的幂迭代浮点运算顺序不同会带来微小的舍入误差。这些差异在1e-7量级对实际预测结果的影响可以忽略。但有个用途如果在自己的项目里需要确认另一个实现是否可靠可以用这种“和sklearn对比差值是否在1e-6以内”的办法做回归测试非常有效。6. 实操中那些文档不会告诉你的经验6.1 成分数量选择不能只看累计解释方差很多人用PLS时会参考PCA的习惯根据X的累计解释方差来选成分数。但PLS是监督方法潜变量不仅要解释X还要和y相关。只看X方差容易选出一些对y毫无解释能力的“大方差噪声成分”。我自己的做法是做交叉验证对不同成分数分别计算验证集RMSE画出一条RMSE随成分数变化的曲线。通常情况下RMSE会先快速下降然后进入平台期最后因过拟合略微上升。选那个进入平台期对应的最小成分数——如果再增加成分RMSE下降不显著就没必要多算。需要特别注意的是交叉验证的每一折内部都要重新计算标准化的均值、标准差不能直接用全量数据的标准化结果。否则又会出现前面说的参数泄漏。6.2 收敛失败时先检查标准化和初值如果NIPALS迭代一直不收敛不要急着调大max_iter先按下面顺序排查数据是否包含NaNNaN在矩阵运算里会像病毒一样传染导致整个矩阵失效。标准化是否成功如果某个变量标准差为0比如工程上常见的常量列会出现除以0的问题。要么删除该列要么给它加一个极小扰动。初值u的选择是否合适在PLS2中u一般取Y矩阵中方差最大的一列。若所有y都接近常量迭代方向就会随机漂移。我调试时还有一个习惯在每次内层迭代打印一下np.linalg.norm(u_new - u)如果这个值振荡不降或者跳到NaN往往不是算法公式错而是前面的标准化环节埋了雷。6.3 从PLS1扩展到PLS2的改动现实中常常需要同时预测多个响应变量比如用光谱数据同时预测水分含量、蛋白质含量和脂肪含量。这就是PLS2多输出PLS。代码需要改动的地方其实很小u初值从Y矩阵中选方差最大的一列q从标量变成向量q (Y_res.T t) / (t t)形状为(m,)m是输出维度u_new Y_res q / (q q)deflation时Y_res Y_res - np.outer(t, q)最终回归系数B_std用(p, n_components) (n_components, n_components) (n_components, m)的维度计算还原尺度时要注意对每个输出列分别乘对应的y_std。只要理解了单输出版本的循环逻辑扩展到多输出其实只是把几个标量换成向量矩阵没有本质难度。但有几个易错点值得留意Q矩阵的维度、Y标准化的方式、还原时intercept需要按列计算。我建议在扩展之前先把单输出版本的推导过程完整写一遍。6.4 数值稳定性和性能提升思路手写实现最大的好处就是能针对具体计算环境调优。如果矩阵规模较大可以把np.linalg.pinv(P.T W)替换成np.linalg.solve(P.T W, np.eye(n_components))这样更快。如果数据维度极高比如光谱数据上万列NIPALS的每次迭代都会做一次大型矩阵乘法计算压力不小。此时可以考虑Kernel PLS用核矩阵代替原始X的运算但这是另一个主题了暂时不展开。对于普通规模的数据几千行、几百列纯numpy实现已经足够快。因为numpy底层的矩阵乘法是高度优化的BLAS实现手写循环最忌讳在Python层逐元素遍历矩阵只要用向量化矩阵运算性能就完全在可接受范围内。最后说点实在的。我在写完这套实现之后最大的收获不是“会手写PLS了”而是以后再遇到奇异矩阵、共线性诊断、变量重要性这类问题时脑子里会自动浮现出矩阵迭代的过程而不是一个黑盒的输入输出。你在实现中如果真的遇到问题先别急着怀疑算法公式——回到数据标准化那里看一眼十次里有八次是那里出了问题。