简介本资源是一份面向数据科学初学者与光谱分析实践者的SPA连续投影算法入门教程聚焦高维数据降维这一核心问题特别适用于生物信息、化学计量学及机器学习预处理场景。压缩包共14个文件6.55MB包含11个MATLAB函数脚本如spa.m、spa_train.m、projections_qr.m等实现算法主流程、种群初始化、验证评估与误差统计、2个光谱数据集.mat格式含jasperRidge2_R198.mat等真实实验数据以及1份详尽的中文教程文档.docx系统覆盖原理推导、代码调试要点、迭代机制说明与光谱应用实例。已有2413人学习下载读者可直接运行调试后的完整工程复现从数据加载、投影迭代、维度裁剪到效果评估的全流程并获得针对大规模光谱数据优化的内存管理与计算效率改进思路。 做近红外光谱建模这些年我最早接触SPA连续投影算法时其实并没太当回事。当时手上有个玉米水分定量项目光谱变量七百多个模型始终绕不开多重共线性PCR、PLS、正则化一层层加效果只能算差强人意。一个老同事甩了句“你先把变量筛了”我才开始认真看这个算法结果一用就是好几年。SPA这名字在化学计量学圈子里没人陌生但圈外人更容易联想到“单页应用”实际上在光谱建模里它是Successive Projections Algorithm的缩写核心就一件事从高维光谱变量里用投影的方式挑出少数彼此冗余最小、代表性最强的波长供后续MLR、PLS这类线性模型使用。这篇文章我想把连续投影算法的原理、完整Python实现、参数调节经验和各种翻车教训一次讲透适合正在做近红外/中红外光谱定量建模、搞特征波长筛选或者被变量共线性折磨的同行参考。如果你只想快点跑通结果直接从第三部分抄代码也行。1. 全谱进模型的代价变量筛选不只是“为了图快”1.1 一个700变量的模型为什么越调越虚很多刚入门的同学有个直觉模型预测不准就加特征。光谱数据恰恰相反波长点越加模型未必越准甚至可能更虚。因为近红外光谱的相邻波长之间相关性极高1098 nm和1100 nm两个点的吸光度往往相差极其微小放进模型之后它们携带的信息高度重叠结果就是设计矩阵近似奇异最小二乘解变得极不稳定。数值上稍微动一下回归系数就剧烈震荡训练集表现好得吓人验证集直接垮掉。这还不只是数值问题。全谱建模往往让模型变成“记忆光谱形状”而不是“学习化学信息”。你拿样本A的平滑基线、样本B的噪声毛刺都可能被当成有效特征。变量一旦多到几百上千模型解释性也大打折扣品控人员问你“到底哪个波数段对应什么基团”你看着三百个系数根本答不上来。所以变量筛选对光谱建模不是可选项而是必要步骤。它的目标很直接保留与目标性质相关的有效信息去掉冗余波长让模型更稳、更简、更可解释。1.2 SPA解决的是哪一类问题连续投影算法解决的核心问题叫做“最小化变量间的共线性”。它不关心某个波长和性质之间的相关性大小也不依赖回归系数绝对值的高低而是用一套几何规则从候选波长里选出相互之间夹角尽可能大、冗余尽可能小的那一组。这种思路特别适合线性模型因为线性模型的稳定性强烈依赖自变量矩阵的“良态程度”变量之间越接近正交最小二乘解越稳定。用一句话概括它做的事情是“从一群长得很像的波长里挑出几个长得最不一样但又不至于完全无关的”为后续回归提供一组干净、低冗余的特征组合。我见过不少场景里SPA筛完变量后即便退回到普通多元线性回归MLR预测精度也能和全谱PLS打平甚至更好但模型变量只剩个位数这在实际产线落地时价值非常大。1.3 别把变量筛选等同于降维很多材料把SPA和PCA混在一起讲这里要分清楚。PCA是把原始变量做线性组合生成新的潜变量得到的是“综合变量”它的每个成分都是一堆原始波长的加权和解释性弱一些SPA则是直接从原始波长里挑子集保留的仍然是真实波长点。两者用途不同也不冲突。实际操作中有人会先用PCA去噪、再用SPA选波长也有人直接对原始光谱做SPA再进PLS两条路我都试过后者更常用因为它保留了波长对应的化学意义。接下来我们就从数学底层看SPA这个“挑法”到底是怎么算出来的。2. SPA到底在算什么连续投影的几何直觉与数学细节2.1 把每个波长当成一个向量理解SPA的第一个关键动作把每个波长在所有样本上的吸光度值看作一个N维空间里的向量。这句话非常抽象我拆开讲。假设你有80个样本那么任意一个波长点比如1100 nm就有80个吸光度值这80个值可以拼成一个80维向量。700个波长点就对应700个这样的高维向量。两个波长之间的相关性反映到几何上就是两个向量的夹角夹角越小两个变量越“同向”信息冗余越严重夹角接近90度两个变量几乎互不相关夹角超过90度两个变量是负相关这在光谱变量的位置上偶尔会出现但本质上也是一种低冗余关系。既然要找低冗余变量直觉就是找彼此夹角尽量大的波长向量。SPA做的事情正是通过逐次投影从原始变量集合里选出这样一组“彼此夹角尽量大”的波长。这里说个大白话类比你站在一片树林里想看完整片林子的全貌只朝一个方向看肯定什么都看不清需要找几个差异很大的视角每个视角都尽量带有增量信息。SPA就是自动帮你找这些视角的算法。2.2 连续投影的迭代规则每次选“最不听话”的变量连续投影的名字里“连续”两个字指的是投影过程一轮接一轮。算法从某个初始波长开始每一步都在剩余变量中寻找“当前已选波长向量方向上的正交投影残差最大”的那个变量。具体公式我给你简化一下。假设当前选中的向量是x_selected对任意候选变量x_j它向x_selected方向做投影后剩余的部分是residual_j x_j - (x_j · x_selected) / (||x_selected||^2) · x_selected这个残差向量就是x_j中“与x_selected无关”的那部分信息。残差向量的长度二范数越大说明x_j携带的独立信息越多。SPA每一轮就选这个长度最大的变量作为下一个入选者。再解释一下直觉。如果x_j和x_selected长得几乎一样夹角很小那么投影项几乎能完整解释x_j残差会很小甚至接近0意味着再选它进来纯属浪费名额。反过来如果x_j和x_selected几乎垂直投影项非常弱残差保留得很大说明它带来的是新信息。所以每一轮选残差最大本质上就是选“现有选择尚未覆盖的信息量最大”的那个波长。这个操作重复执行直到选出预设数量h个变量。整个“链式选择”过程会从不同初始波长开始跑多遍形成多条候选链最终用交叉验证决定哪条链最靠谱。这后一步很多人容易忽略但恰恰是SPA最后一步“定稿”的关键。2.3 一个二维小实验看SPA如何挑掉冗余为了让你更直观地记住这个原理我给你写个极小的模拟例子。假设只有3个样本4个候选波长A、B、C、D它们的向量分别是A [1, 2, 3] B [1.1, 2.1, 2.9] C [5, -1, 0] D [0, 3, -2]人工看一下B和A几乎平行夹角只有几度所以如果先选中A那么B计算出来的投影残差会非常小基本不具备入选资格。C、D与A夹角更大残差相对大进入下一轮候选。假设第一步选中A第二步分别算B、C、D对A的投影残差。B的残差长度大概率在0.2以下C和D的残差可能都在6左右选C或D具体看哪个更长。等选入C之后第三步计算剩余变量对“A和C张成平面”的正交投影进一步选出在这个平面上最“独立”的变量。这就是SPA逐步挑选的核心过程一旦选定一个波长下一轮的剩余信息是排除掉所有已选变量解释过的部分之后重新计算出来的。这个设计思路带来的直接好处是已选变量之间不仅单个相关度低整体组合的多重共线性也受到抑制。实际数据里常见的情况是SPA选出的波长往往分布在光谱的不同区间而不是扎堆在某个官能团吸收峰附近这和它最大残差的选择逻辑完全吻合。3. 手写Python实现从伪代码到能直接跑通的函数3.1 实现前的设计取舍网上能找到不少SPA的封装工具比如有的人用MATLAB版本也有人用R包。但我在实际项目中还是更习惯自己写Python实现原因有三个一是方便和后续建模代码衔接二是变量选择需要配合嵌套交叉验证自己封装的函数可以灵活插入到不同验证流程里三是纯numpy实现在几百个变量这种规模下性能绰绰有余没必要依赖外部软件。实现SPA我建议拆成三个函数来写投影步骤函数负责输入当前已选向量集合、候选向量矩阵输出每个候选向量的投影残差范数候选链生成函数从每个初始波长出发生成一条包含h个波长的变量序链最终优选函数用交叉验证RMSECV比较不同链返回最优变量下标。这么拆的好处是代码可测性好。投影步骤单独拿出来你可以用前面二维模拟数据验证结果等确认投影步骤没问题整个算法大概率就稳了。3.2 核心代码逐段注释下面是我平时在项目里实际使用的核心实现做了简化去除了一些项目专属的打印日志但保留了最骨干的逻辑。import numpy as np def calc_projection(Xcal, selected_idx): 计算 Xcal 中每个候选变量在已选变量张成空间上的正交投影残差范数。 Xcal: (n_samples, n_variables) selected_idx: list, 当前已选变量的列下标 返回: norm_array, 长度 n_variables 的数组已选变量位置置0 n_samples, n_vars Xcal.shape norm_array np.zeros(n_vars) # 选出一个当前“代表向量”这里取最后一个已选列 # 更严格的做法会使用整个已选列空间但SPA最常用的是逐列投影近似 x_sel Xcal[:, selected_idx[-1]].copy() denom np.dot(x_sel, x_sel) for j in range(n_vars): if j in selected_idx: norm_array[j] 0.0 continue x_j Xcal[:, j] # 计算 x_j 在 x_sel 上的投影系数 coef np.dot(x_j, x_sel) / denom # 投影残差 residual x_j - coef * x_sel norm_array[j] np.linalg.norm(residual) return norm_array def spa_chain(Xcal, start_idx, h): 从 start_idx 出发选择 h 个变量返回变量下标list。 Xcal: (n_samples, n_variables) n_samples, n_vars Xcal.shape selected [start_idx] for _ in range(1, h): norm_array calc_projection(Xcal, selected) next_idx int(np.argmax(norm_array)) # 避免重复选中同一个变量 if next_idx in selected: # 如果最大已经是已选变量取第二大 norm_array[next_idx] -1 next_idx int(np.argmax(norm_array)) selected.append(next_idx) return selected这里需要说明一点上面代码逐列投影用的是最后一个已选向量这是SPA原始论文中比较常见的一种简化处理方式。严格的做法第二步应该计算候选变量对“所有已选变量张成的线性空间”的正交投影而不是只针对最后一个已选变量。但对大量实际光谱数据来说逐列投影已经能选出足够好的变量子集而且计算复杂度低很多。如果你介意可以把投影步骤改成# 严格版本投影到整个已选列空间的正交补空间 U, S, Vt np.linalg.svd(Xcal[:, selected_idx], full_matricesFalse) projection U U.T residual x_j - projection x_j不过这个严格版在变量数多时会频繁做SVD速度不如逐列投影快项目里我一般不用。作为教程版逐列投影已经足够解释算法核心效果上也不输太多。3.3 用交叉验证定稿不是所有链都合适前面提到SPA对每个起始变量都会生成一条变量链那最后到底选哪条链答案是交叉验证。以一个近红外数据集为例流程是这样的def spa_select(Xcal, ycal, h): 完整SPA变量选择流程。 遍历所有候选初始变量根据MLR交叉验证RMSECV决定最优变量链。 n_samples, n_vars Xcal.shape best_rmse np.inf best_indices None for start in range(n_vars): indices spa_chain(Xcal, start, h) Xsel Xcal[:, indices] # 用多元线性回归做留一交叉验证LOO-CV rmse loo_mlr_cv(Xsel, ycal) if rmse best_rmse: best_rmse rmse best_indices indices print(fstart{start}, selected{indices}, RMSECV{rmse:.4f}) return best_indices, best_rmse def loo_mlr_cv(Xsel, y): 留一交叉验证每次留一个样本用其余样本拟合多元线性回归 然后预测留出样本计算均方根误差。 n len(y) preds np.zeros(n) for i in range(n): train_idx [j for j in range(n) if j ! i] X_tr, y_tr Xsel[train_idx], y[train_idx] X_te Xsel[i:i1] # 加一列截距 X_tr_b np.column_stack([np.ones(len(train_idx)), X_tr]) X_te_b np.column_stack([np.ones(1), X_te]) beta, _, _, _ np.linalg.lstsq(X_tr_b, y_tr, rcondNone) preds[i] X_te_b beta rmse np.sqrt(np.mean((y - preds) ** 2)) return rmse这里我用了留一交叉验证配合MLR。留一法在小样本光谱数据里很常见虽然计算量不小但变量筛选阶段的样本量通常就几十到一两百还可以接受。如果你觉得慢也可以改成5折交叉验证速度会快不少选择结果差异通常不大。跑完之后best_indices就是最终选出的波长下标。整个流程里我个人最喜欢的部分就是盯着print输出看一开始RMSECV波动很剧烈到后面逐渐收敛最后选出的一组变量通常稳定在某个数值附近那种感觉非常踏实。4. 参数调节和翻车现场SPA最容易出错的几个环节4.1 h选多大别拍脑袋画一条RMSECV曲线SPA最关键的参数是h也就是要选多少个波长。h太小信息量不足模型欠拟合h太大冗余信息又进来了模型重新变得不稳。我建议的做法是先设定一个h的上限比如min(20, 样本数-3)然后从1到上限逐个计算对应的最低RMSECV画一条“变量数-RMSECV”曲线看拐点。在多数近红外项目里拐点往往出现在6到12个变量之间。你看到RMSECV开始下降到平台期、再增加变量只能带来微小改善甚至变差时就该收手。别贪变量少一点点模型的可维护性和可解释性都完全不同。4.2 起始变量遍历的代价和优化我在3.3节的代码里用了最笨的办法对所有变量都作为起始变量跑一遍。如果变量数有700个h取8每次留一交叉验证要做几百次最小二乘总时间可能高达几十秒甚至几分钟。实际项目里我通常会做两个调整。第一对原始光谱做一次简单的排序筛选剔除明显噪声区间比如光谱两端基线漂移严重的区域从剩下的候选变量里再挑起始点第二限制起始变量候选数比如随机从700个变量里抽80个作为起始变量池或者按方差从大到小取前100个。这不会明显降低精度却能节省大量时间。SPA对初始变量确实敏感所以原始算法才用“全变量遍历交叉验证”来规避随机性。如果你做随机抽样记得多抽几次或者把每次结果交叉验证对比一下选最稳的。4.3 数据预处理的顺序坑过我一次这个坑我必须单独拿出来说。SPA是对光谱变量做几何计算所以数据的尺度会影响投影结果。如果先做了标准正态变换SNV或者多元散射校正MSC算出来的残差范数和直接在原始光谱上算选出的波长很可能是不同的。问题来了到底该在预处理前选还是预处理后选我的经验是一定要在最终采用的预处理流程之后再做SPA。也就是说先做基线校正、散射校正、平滑、中心化再跑变量筛选。理由很简单变量筛选是为后续建模服务的后续模型输入的已经是预处理后的光谱筛选也必须基于同样的数据分布。千万不要在原始数据上筛完变量再把预处理后的数据喂进模型那会导致筛选结果和模型输入不匹配。4.4 重复波长和“抱团”现象光谱分辨率高时波长间隔很小SPA偶尔会选出两个非常邻近的波长比如1120 nm和1122 nm。从投影残差角度看它们可能有细微差别所以都被选中了但这种“抱团”对模型稳定性没啥帮助也可能说明h选大了。我的处理办法是把SPA选出的变量按波长排序算相邻变量间隔如果发现有一对间隔小于仪器分辨率的波长就手动剔除其中一个重新做一次MLR或PLS交叉验证看看误差没有明显增大。若误差没变就接受更紧凑的变量集。这种“算法先筛人工微调”的组合在品控场景下特别实用因为你可以向非算法人员解释清楚为什么选这几个波长。4.5 避免数据泄漏SPA必须在交叉验证折内执行接下来这个坑处理不好会让后续所有评价结果虚高。如果你在训练集上跑SPA选出波长然后用这些波长去预测验证集这本身就是一种数据泄漏。因为变量选择过程“偷看”了整体数据分布在验证集上的误差会被低估。正确的做法是在嵌套交叉验证里面做外层把数据分成训练和测试内层只使用训练部分做SPA然后再把选出的波长应用在外层验证集上。如果你不做嵌套验证至少也要把SPA放在训练集内部做测试集从始至终不能参与变量选择。这个道理和标准化参数只从训练集估算是一样的但很多人会忽略。5. 什么时候该选SPA与CARS、UVE、VIP等方法对比5.1 不同变量筛选方法的核心差异化学计量学里变量筛选方法远不止SPA一种。我在项目里还常用到CARS、UVE和VIP它们在原理上差异很大。方法核心思路输出形式典型特点SPA几何投影选择低冗余波长波长子集稳定、可解释适合线性模型CARS基于回归系数的自适应收缩 随机采样波长子集擅长处理强相关变量多次运行结果有一定随机性UVE加入伪变量比较回归系数的t统计量波长子集剔除不稳定变量但阈值选择有主观性VIP基于PLS投影重要性分数排序变量重要性排序保留重要性高的变量需要自己卡阈值SPA最突出的差异化优势是“完全不需要响应变量来做投影计算”投影过程只依赖光谱矩阵本身响应变量只出现在最终的交叉验证择优环节。这一点听着像优点其实也是双刃剑它不考虑波长与化学成分的相关性选出来的变量可能是低冗余的但不一定都和目标性质高度相关。好在后续回归会对目标性质建模不相关的低冗余变量往往会被赋予接近零的系数影响不会太大。5.2 数据条件和应用场景匹不匹配SPA最适合的数据场景是波长变量数远大于样本数、变量之间强相关、目标性质和光谱之间基本呈线性关系。近红外光谱正好都满足。在这种条件下SPA选出的低冗余变量组合能显著提升MLR的稳定性甚至很多时候用MLR就能逼近PLS的效果。为什么强调MLR因为产线部署时MLR模型就是一个线性方程单片机、PLC都能轻松实现维护成本极低。这是SPA方案的落地价值所在。如果数据存在明显的非线性效应比如固体药片近红外光谱受颗粒粒径非线性散射影响很大SPA配MLR就不够用了。此时可以先做SNV/MSC校正仍然用SPA选变量但是模型换成支持向量回归或者神经网络效果会更好。SPA本身没有能力“修复非线性”它只管给模型提供低共线性的输入。5.3 哪些场景建议绕开SPA样本量非常小、只有二三十个且信噪比也很低时SPA投影计算会放大噪声变量在残差中的比重选出重复性差的光谱区间。这时候我倾向于用UVE或CARS先剔除噪声变量再选代表性波长。另一个需要绕开的场景是波长区间本来就很少比如只有十几个特征变量时SPA往往选不出明显优于全谱的组合花哨算法不如老老实实做逐步回归或者直接看化学意义选峰。在遇到“开头几分钟就能跑出结果但结果每次运行都不一样”的情况比如某些方法引入随机采样如果做变量筛选是为了后续复现和品控SPA的确定性反而是优点。同样的输入SPA每次输出完全一致这在GMP/GLP相关项目里特别好用。6. 一次真实的近红外建模实验700个波长压到6个6.1 数据与实验设计去年我帮一家饲料企业做玉米粗蛋白近红外定量模型设备是国产光栅型近红外光谱仪波长范围1100到2498 nm分辨率2 nm共700个波长点。样本一共86份按4比1划分训练集和测试集训练集69份测试集17份。光谱预处理用了SNV加一阶导数然后均值中心化。初始尝试直接全谱PLS主成分数用交叉验证选出8个测试集R²为0.91RMSEP为0.47。模型能用但每次仪器漂移后需要重新标定而且8个主成分对应的是全谱载荷品控团队完全无法解释模型“在看什么”。于是我开始用SPA做变量筛选。6.2 SPA实验过程和结果对比用3.3节的代码跑一遍h从1到12RMSECV曲线显示h6之后基本进入平台期再往上走误差下降不超过0.01。取h6最优链选出的波长分别为1128 nm、1342 nm、1466 nm、1698 nm、1934 nm、2188 nm。这些波长跨度非常大分别对应C-H键的二倍频和组合频区域刚好和蛋白质相关官能团有对应关系。然后我用选出的6个波长做MLR结果测试集R²为0.89RMSEP为0.52比全谱PLS略差但差异在可接受范围而且模型形式变成蛋白含量 a0 a1·A(1128) a2·A(1342) a3·A(1466) a4·A(1698) a5·A(1934) a6·A(2188)模型里只有6个波长点设备测温漂移、样本粒度变化造成的干扰更容易追溯到具体波长上。后面我又加了一组对比实验用SPA选出的6个波长放入PLSR²升到0.90RMSEP降到0.49几乎追平全谱PLS。这进一步说明光谱数据里真正有效的信息确实可以用极少数波长承载。6.3 落地时的几个小经验实验做完之后我把这套模型给产线测试了两个月发现一个额外收益模型稳定性比全谱模型好很多。因为变量少、波长位置固定当仪器光源强度波动导致全谱基线整体平移时对单个波长点做外部标准化校正比校正整条光谱容易得多。每次取一个样本用化验值校正截距b0就行。另外落地时还要考虑波长在仪器上的漂移问题。光栅光谱仪的温度变化会导致波长微小偏移如果筛选出的几个波长恰好落在窄峰上偏移对吸光度影响会很大。所以我建议SPA选完波长后看看这些波长所在位置的峰型宽窄如果某个波长落在非常尖锐且窄的吸收峰上要特别谨慎最好做一次波长稳定性测试或者适当选择吸收峰附近较宽的肩部位置。从那次项目之后我的近红外建模流程基本固定了预处理 → SPA或CARS筛变量 → MLR/PLS建模 → 交叉验证评估 → 部署。不同项目再根据实际情况微调但SPA始终是我第一个会去尝试的工具。它在数学上不复杂、代码实现可控、结果确定可复现对做定量模型的人来说是性价比极高的一个选择。如果你现在正被全谱变量太多、模型不稳的问题困扰不妨把你的数据集拿出来跑一遍SPA看看只需要几个波长时模型会变成什么样——很多时候那个最简单的模型反而最耐用。本文还有配套的精品资源点击获取