企业官网建设流程全解析

说实话，第一次接触偏最小二乘回归（PLS）的时候，我整个人是有点懵的。网上一搜全是“潜变量”“载荷矩阵”“交叉验证”这些词，配套的公式密密麻麻，看起来就像某个数学系博士生用来折磨人的毕业设计。但后来因为工作里要部署一个跨平台的小工具，没法把整套Python环境搬过去，只能手写核心推理逻辑，我才被迫把这层窗户纸捅破。等真把NIPALS算法一步步拆开之后，我只有一种感受：PLS的核心逻辑，拆开了就是矩阵来回折腾——这正是这篇文章想和你一起做的事。

这篇文章不依赖任何机器学习的黑盒库，从头到尾用numpy手写PLS的核心实现，重点包括数据标准化的必要性、潜变量提取的迭代细节、残差矩阵更新方式，以及回归系数从标准化空间还原到原始尺度时的换算逻辑。整个过程我会围绕PLS的矩阵运算展开，把代码逐段讲明白、把每一步动作的理由说清楚。不管你是刚接触多元统计分析的学生，还是需要在嵌入式环境或纯计算场景下手推算法的工程师，这篇文章都适合你跟着敲一遍。

1. 为什么“手动实现”这件事值得做

1.1 调库和懂原理之间隔着一条认知断层

我用sklearn用得很早，PLSRegression(n_components=5).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 手写实现的基本路线图

整篇文章的实现路线很清晰，一共四步：

1. 数据标准化：把X和y的均值、方差统一到同一量纲，这是所有协方差计算的前提。
2. 迭代提取潜变量：用NIPALS（非线性迭代偏最小二乘）算法，一次提一个成分，每次提取后更新残差矩阵。
3. 拟合回归系数：在标准化空间里利用权重W、载荷P和y载荷Q计算回归系数B。
4. 还原尺度：把标准化空间的系数映射回原始数据尺度，算出截距，得到最终预测公式。

下面的内容，我就按照这条路线一步步展开。

2. 数据标准化：PLS的第一步，也是很多人忽略的定海神针

2.1 标准化到底在做什么

很多人在入门时觉得标准化就是“减均值除以标准差”，好像很机械，但这一步在PLS里不是例行公事，而是直接参与算法核心计算的。实现上就是两个动作，代码很轻量：

X_mean = X.mean(axis=0) X_std = X.std(axis=0, ddof=0) # 总体标准差，和sklearn默认行为一致 Xs = (X - X_mean) / X_std y_mean = y.mean() y_std = y.std(ddof=0) ys = (y - y_mean) / y_std

这里特别提醒一个细节：ddof参数的取值。ddof=0算的是总体标准差，分母是n；ddof=1算的是样本标准差，分母是n-1。sklearn.preprocessing.StandardScaler默认使用ddof=0，sklearn的PLSRegression在scale=True时也按这种方式标准化。如果你想和现有工具链的结果对得上，这里必须统一，否则预测值可能差一个小数点后几位的数值。

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_size=0.3, random_state=42) # 在训练集上计算参数 X_mean = X_train.mean(axis=0) X_std = X_train.std(axis=0, ddof=0) y_mean = y_train.mean() y_std = y_train.std(ddof=0) # 用训练集参数变换训练集和测试集 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的原始变量“解释”了多少方差。
• q（y载荷）：y方向上的载荷标量，q = y^T t / (t^T t)，可以理解为t与y的最佳回归斜率。

整个NIPALS循环，本质上就是在反复求解“X的权重”和“y的载荷”，直到得到收敛的一组潜变量。

3.2 单个潜变量提取的完整循环

假设现在只提取一个成分，具体步骤如下。这个循环也是整个PLS迭代的核心：

1. 初始化u为y（在PLS1中，u直接取残差y向量即可；在PLS2中，一般取Y矩阵中方差最大的一列）。
2. 计算X的权重：w = X^T u / ||X^T u||。这一步算的是X每一列与u的相关性，归一化保证方向一致。
3. 计算X的得分：t = Xw。这是所有样本在那个新方向上的投影。
4. 计算y的载荷：q = (y^T t) / (t^T t)。这是t对y做最小二乘回归得到的系数。
5. 更新u：u = yq。
6. 检查收敛：如果新的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 多成分提取的关键：残差矩阵的更新（deflation）

PLS和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, tol=1e-8, max_iter=500): # 1. 保存原始尺度参数并标准化 X_mean = X.mean(axis=0) X_std = X.std(axis=0, ddof=0) y_mean = y.mean() y_std = y.std(ddof=0) 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, Q

4.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_size=0.3, random_state=7) # 手写实现 B_manual, intercept_manual, W, T, P, Q = pls_manual(X_train, y_train, n_components=3) 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_components=3, scale=True) 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(f"RMSE 手写实现: {rmse_manual:.6f}") print(f"RMSE 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，先按下面顺序排查：

• 数据是否包含NaN？NaN在矩阵运算里会像病毒一样传染，导致整个矩阵失效。
• 标准化是否成功？如果某个变量标准差为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了”，而是以后再遇到奇异矩阵、共线性诊断、变量重要性这类问题时，脑子里会自动浮现出矩阵迭代的过程，而不是一个黑盒的输入输出。你在实现中如果真的遇到问题，先别急着怀疑算法公式——回到数据标准化那里看一眼，十次里有八次是那里出了问题。