2026/8/31 15:38:44
FEATURED · 精选文章 的C语言完整实现)
简介本资源是一套轻量级C语言实现的偏最小二乘回归PLS算法代码专为嵌入式场景优化可直接部署于STM32等资源受限的单片机平台面向嵌入式AI开发者、工业传感器数据建模工程师及学习嵌入式机器学习算法的进阶学习者。它解决了在无Python/R环境、无浮点协处理器、内存≤128KB条件下仍需对多变量高共线性传感数据如温湿度、气体浓度阵列进行实时回归预测的实际问题。压缩包共4个文件2个C源码2个头文件总大小仅8KB其中pls.c封装了PLS核心迭代逻辑与因子提取流程matrix.c/h提供精简矩阵运算支持含转置、乘法、标准化整体采用定点/单精度浮点混合策略避免动态内存分配适配CMSIS-DSP或裸机运行。目前已有1599人学习下载读者可直接复用该模块构建边缘端软测量模型、设备健康度预测或校准补偿算法具备完整可编译结构、清晰接口定义与嵌入式友好注释。 做嵌入式这几年一个比较深的感触是越是工业现场实用的算法越难找到能直接在单片机上跑的C语言版本。前段时间接了个项目要在ARM Cortex-M上做多传感器数据融合用几个便宜传感器测出来的信号去预测一个关键指标。一开始想用最小二乘但传感器之间有耦合特征共线性严重效果很差换神经网络MCU上跑不动参数也爆炸。后来查资料发现PLS偏最小二乘回归特别适合干这个活但在网上翻了很久Python和MATLAB的教程一大堆C语言的完整实现愣是没找到几个能用的。这篇文章把这段时间啃下来的东西整理出来从PLS原理、C语言NIPALS算法实现到STM32移植部署的完整过程。如果你正需要在单片机或嵌入式设备上做多变量回归、传感器校准、光谱分析或者单纯想找个能直接用的C语言实现来抄作业这篇应该能帮你省下不少时间。1. 项目概述与整体设计思路1.1 为什么是PLS而不是多元线性回归或神经网络先说我为什么最终选PLS。当时的场景是一块STM32F407主控板接了温度、湿度、气压、光照四个传感器想用这四个值去预测某种气体浓度。遇到的第一个坑就是共线性——温度和气压之间存在明显的相关性湿度又和温度耦合用普通最小二乘回归MLR回归系数方差被共线性放大换一组数据预测结果就飘了。PLSPartial Least Squares偏最小二乘解决的正是这个问题。它的思路不是直接拿原始变量做回归而是先把原始变量压缩成若干个彼此正交的潜变量再拿潜变量和因变量做回归。因为潜变量之间正交共线性的问题天然被规避了。同时PLS提取潜变量的时候会同时考虑潜变量对自变量的解释力度和潜变量对因变量的解释力度所以压缩出来的特征不会像PCA那样只盯X最大的方差而和y无关。至于神经网络在MCU上跑的痛点很现实模型参数动辄几千上万Flash和RAM都扛不住推理时大量乘加运算没有NPU基本跑不动而且数据量如果只有几十上百条网络训练极容易过拟合。PLS的模型压缩到极致就是一条线性等式预测阶段每次只要做n次乘加——n是自变量个数这在单片机上可以忽略不计。1.2 整体架构离线训练在线预测工程上我的方案是离线训练在线预测两步走。首先在PC上利用Python、MATLAB或SPSS等工具完成PLS模型的训练、交叉验证和主因子数的确定然后把训练得到的回归系数固化成一个C语言头文件烧录到STM32的Flash里。STM32收到传感器数据后只需要走一遍极简的预测函数就能算出结果。这套架构的好处很直接MCU侧代码量极小只做加减乘除对算力和内存几乎没有压力。模型训练、验证、调参都在PC上完成可调试性远好于在板子上训练。Flash里随时可以放多组模型参数实现不同场景切换。但我也在C语言里实现了完整的NIPALS训练算法并实测可以在STM32上跑通。如果数据量小几百条以内、十几个变量运算时间在几百毫秒到几秒量级MCU端直接训练是可行的。这对一些需要在线自适应更新模型的场景很有价值。1.3 适用场景与数据规模结合我接触到的实际需求PLS在嵌入式端最典型的应用场景有三类多传感器阵列校准比如电子鼻、气体检测仪的交叉敏感校正。近红外光谱的定量分析比如用光谱反演农作物水分、蛋白质含量。工业过程软测量用温度、压力、流量、振动等信号预测某个难以在线检测的质量指标。数据规模上MCU端通常处理个位数到几十个自变量、几百条训练样本的规模。在STM32F4这种级别的平台上只要内存合理规划完全跑得动。2. 偏最小二乘回归原理与核心算法2.1 数学原理从投影到回归PLS回归的核心思想可以概括成一句大白话找一组彼此正交的潜变量方向让自变量X在这个方向上的投影得分t和因变量y最相关然后用这些投影做回归。用数学语言描述X是一个n×m的矩阵n个样本m个自变量y是n×1的向量。对每个主因子k找到权重向量w_k使得得分t_k Xc * w_k且t_k与y的协方差最大。计算X的载荷向量p_k Xc * t_k / (t_k * t_k)以及y侧的系数c_k y * t_k / (t_k * t_k)。从Xc和y中去掉当前主因子的贡献这一步叫deflation继续提取下一个主因子。最终模型可以还原成一条线性形式y_pred b0 Σ(beta_j * x_j)其中beta_j由所有主因子下的W、P、c矩阵组合计算得出。这个线性形式是PLS最漂亮的地方——训练阶段虽然复杂但终端预测就是一个普通线性回归的样子单片机跑毫无压力。2.2 NIPALS算法的完整步骤NIPALS是目前求解PLS最主流的迭代算法。针对单因变量的PLS1场景可以简化为以下流程对X和y做中心化或标准化得到Xc和yc。对于每个主因子k1, 2, ..., A计算权重向量 w Xc * yc / (yc * yc)yc是向量时这一步可以简化。归一化ww w / ||w||。计算得分 t Xc * w。计算y侧系数 c yc * t / (t * t)。计算X载荷 p Xc * t / (t * t)。缩减Xc矩阵Xc Xc - t * p。缩减yc向量yc yc - t * c。收集所有W、P、c计算最终回归系数 beta W * (PW)^(-1) * c。需要特别注意单因变量时因为yc本身就是一个向量所以NIPALS内部不需要迭代求收敛一步就能算出w。如果是多因变量的PLS2内层需要迭代直到t方向稳定代码复杂度会高不少。所以本文聚焦单因变量这也是工程中90%以上的场景。2.3 主因子数A到底选几个主因子数A是PLS调参时最关键的旋钮。选少了模型信息提取不充分欠拟合选多了会把噪声也当成信号过拟合。标准的做法是交叉验证。Leave-One-Out或者K折交叉验证每次留出一部分样本做测试记录不同A下的预测误差RMSE或PRESS画出一条误差-主因子数的曲线取误差最低点或者拐点前面的值。实践中我的经验自变量之间相关性越强需要的主因子数通常越少2~5个足够。如果自变量本身比较独立主因子数会接近原始变量数。切忌一味追求最低误差当A增大到某个值后误差下降已经很缓慢甚至回升这个拐点就是合理选择。MCU上做在线训练时无法跑复杂交叉验证可以先用固定A比如3快速训练后期在PC上统一做一次交叉验证修正主因子数。3. C语言数据结构与基础矩阵运算3.1 矩阵表示一维数组加行列数在MCU上实现矩阵运算我强烈建议用一维数组存储矩阵再用结构体记录行列数。比如typedef struct { float *data; int rows; int cols; } Matrix;为什么不用常规的二维数组float data[rows][cols]三个原因C语言二维数组传参时列数必须是编译期常量函数通用性大打折扣。二维数组的内存布局虽然是连续的但动态分配二维数组时很容易产生内存碎片多次malloc/free容易把MCU的堆搞乱。一维数组配合行主序存储访问元素data[i * cols j]不仅代码简洁而且编译器能更好地做循环优化和缓存预取性能优于二维数组的间接寻址。当然如果是完全静态的固定大小场景直接用静态数组宏定义也可以省去动态分配的开销。我工程里常用静态数组因为嵌入式项目的数据维度在编译期就已经确定。3.2 几个绕不开的矩阵运算函数在实现NIPALS之前先把基础矩阵运算准备好。最常用的有四个矩阵乘法C(m×n) A(m×k) * B(k×n)。转置乘C(m×n) A(k×m) * B(k×n)用于高效计算X*y、X*t。向量点积和模长。方阵求逆用于最终计算beta W * (PW)^(-1) * c。矩阵乘法最朴素的三层循环在MCU上其实够用因为NIPALS里的矩阵规模都不大m和n通常不超过几十。优化的重点是把内层循环的索引访问模式理顺避免频繁跨行跳跃。转置乘单独实现一个函数很有必要。如果先转置再普通乘法会额外多一次O(m×k)的搬移操作在MCU这种资源受限的环境里属于不必要的浪费。直接把转置逻辑融进乘法循环里效率高不少。3.3 内存与精度float够不够STM32F4系列带FPU单精度float的乘加运算是硬件加速的速度远超double。所以预测和训练代码我都默认用float。但要注意float的精度只有约7位有效数字当数据数值范围很大或者矩阵条件数很差时累积误差会变大。实测经验传感器原始数据范围在0~10000这个量级时float没问题。训练阶段求逆运算对精度最敏感如果主因子数较多且数据分布不均匀建议求逆时临时用double算完再转回float。如果所有数据都做了标准化减去均值除以标准差数值范围被压到-3~3之间float的误差会小很多训练和预测都能放心用float。我在代码里做了一个折中矩阵运算函数全部用float但求逆函数内部用double做累加这样兼顾速度和精度。4. NIPALS训练算法的C语言实现4.1 训练主函数框架直接看代码。为了能在MCU上跑我用静态数组最大维度用宏定义限定#define MAX_SAMPLES 200 #define MAX_VARS 20 #define MAX_COMP 10 typedef struct { float W[MAX_COMP * MAX_VARS]; float P[MAX_COMP * MAX_VARS]; float C[MAX_COMP]; float mean_x[MAX_VARS]; float mean_y; float beta[MAX_VARS]; float b0; int n_vars; int n_comp; } PLSModel;训练函数输入是样本矩阵X行优先一维数组、因变量y和主因子数A输出填充PLSModel结构体int pls_train(const float *X, const float *y, int n_samples, int n_vars, int n_comp, PLSModel *model) { int i, j, k; float Xc[MAX_SAMPLES * MAX_VARS]; float yc[MAX_SAMPLES]; float w[MAX_VARS], t[MAX_SAMPLES], p[MAX_VARS]; float tt, yt, c; // 1. 中心化 float mean_x[MAX_VARS] {0}; float mean_y 0; for (j 0; j n_vars; j) { float sum 0; for (i 0; i n_samples; i) sum X[i * n_vars j]; mean_x[j] sum / n_samples; } for (i 0; i n_samples; i) mean_y y[i]; mean_y / n_samples; // 复制并中心化数据 for (i 0; i n_samples; i) { yc[i] y[i] - mean_y; for (j 0; j n_vars; j) Xc[i * n_vars j] X[i * n_vars j] - mean_x[j]; } // 2. NIPALS迭代 for (k 0; k n_comp; k) { // 2.1 计算 w Xc * yc for (j 0; j n_vars; j) { w[j] 0; for (i 0; i n_samples; i) w[j] Xc[i * n_vars j] * yc[i]; } // 2.2 归一化 w float norm 0; for (j 0; j n_vars; j) norm w[j] * w[j]; norm sqrtf(norm); if (norm 1e-12f) return -1; // 数值异常 for (j 0; j n_vars; j) w[j] / norm; // 2.3 计算 t Xc * w for (i 0; i n_samples; i) { float sum 0; for (j 0; j n_vars; j) sum Xc[i * n_vars j] * w[j]; t[i] sum; } // 2.4 计算 c yc * t / (t * t) tt 0; yt 0; for (i 0; i n_samples; i) { tt t[i] * t[i]; yt yc[i] * t[i]; } c yt / tt; // 2.5 计算 p Xc * t / (t * t) for (j 0; j n_vars; j) { float sum 0; for (i 0; i n_samples; i) sum Xc[i * n_vars j] * t[i]; p[j] sum / tt; } // 保存 for (j 0; j n_vars; j) { model-W[k * n_vars j] w[j]; model-P[k * n_vars j] p[j]; } model-C[k] c; // 2.6 deflateXc Xc - t * p for (i 0; i n_samples; i) { for (j 0; j n_vars; j) { Xc[i * n_vars j] - t[i] * p[j]; } } // yc yc - t * c for (i 0; i n_samples; i) yc[i] - t[i] * c; } model-n_vars n_vars; model-n_comp n_comp; for (j 0; j n_vars; j) model-mean_x[j] mean_x[j]; model-mean_y mean_y; // 3. 计算最终回归系数 return pls_build_beta(model); }这段代码的思路是标准的NIPALS-PLS1。注意第2步中normalize前先检查模长是否过小这是防止除零的关键。实际数据如果出现某个主因子贡献极小时模长会非常小强行继续迭代只会把数值噪声放大直接返回错误更合理。4.2 求逆高斯-约当消元法计算beta的核心是求解W * (PW)^(-1) * c其中PW是一个A×A的小方阵。A就是主因子数一般不超过10。我用高斯-约当消元实现求逆static int mat_inv(float *A, int n) { float aug[MAX_COMP * MAX_COMP * 2] {0}; int i, j; for (i 0; i n; i) { for (j 0; j n; j) aug[i * 2 * n j] A[i * n j]; aug[i * 2 * n n i] 1.0f; } for (int col 0; col n; col) { // 选主元 int pivot col; float max_val fabsf(aug[col * 2 * n col]); for (int row col 1; row n; row) { if (fabsf(aug[row * 2 * n col]) max_val) { max_val fabsf(aug[row * 2 * n col]); pivot row; } } if (max_val 1e-8f) return -1; // 交换行 if (pivot ! col) { for (j 0; j 2 * n; j) { float tmp aug[col * 2 * n j]; aug[col * 2 * n j] aug[pivot * 2 * n j]; aug[pivot * 2 * n j] tmp; } } // 归一化当前行 float pivot_val aug[col * 2 * n col]; for (j 0; j 2 * n; j) aug[col * 2 * n j] / pivot_val; // 消去其他行 for (int row 0; row n; row) { if (row col) continue; float factor aug[row * 2 * n col]; for (j 0; j 2 * n; j) aug[row * 2 * n j] - factor * aug[col * 2 * n j]; } } for (i 0; i n; i) for (j 0; j n; j) A[i * n j] aug[i * 2 * n n j]; return 0; }选主元这一步很重要。如果矩阵是奇异的或者接近奇异不做选主元直接消元结果会完全不可信。MCU上碰到这种情况返回值-1让上层决定是跳过还是降主因子数。4.3 组装最终回归系数求逆完成后先用W * inv(PW) * c算出权重系数再和中心化参数合并成最终的线性模型static int pls_build_beta(PLSModel *model) { int i, j, k; int A model-n_comp; int m model-n_vars; float PtW[MAX_COMP * MAX_COMP] {0}; float tmp[MAX_COMP] {0}; // 计算 PtW P * W for (i 0; i A; i) { for (j 0; j A; j) { float sum 0; for (k 0; k m; k) sum model-P[i * m k] * model-W[k * A j]; PtW[i * A j] sum; } } // 求逆 if (mat_inv(PtW, A) ! 0) return -1; // 计算 tmp inv(PW) * C for (i 0; i A; i) { float sum 0; for (j 0; j A; j) sum PtW[i * A j] * model-C[j]; tmp[i] sum; } // 计算 beta W * tmp for (j 0; j m; j) { float sum 0; for (k 0; k A; k) sum model-W[k * m j] * tmp[k]; model-beta[j] sum; } // 计算 b0 mean_y - sum(beta[j] * mean_x[j]) float b0 model-mean_y; for (j 0; j m; j) b0 - model-beta[j] * model-mean_x[j]; model-b0 b0; return 0; }这里有个细节容易踩坑W和P在NIPALS迭代里是按行存的k行代表第k个主因子但在求PtW时W的索引用了W[k * A j]这种形式实际上是把W当成了A×m的行主序矩阵。所以做矩阵乘法时要非常小心每一个索引下标别搞混了。我在这个位置调试过整整一个下午最后发现是下标问题。4.4 预测函数极简但够用训练完成后MCU侧的预测就是几行代码的事float pls_predict(const PLSModel *model, const float *x) { float y model-b0; for (int j 0; j model-n_vars; j) { y model-beta[j] * x[j]; } return y; }是的就是这么简单。因为所有中心化参数都被融合进b0和beta里了输入原始传感器值直接算就行。这个函数放在STM32的中断里、RTOS任务里、几kHz的采样循环里都毫无压力。4.5 验证C语言结果与Python结果对比算法写完必须验证。我用Python的scikit-learn里PLSRegression跑了一组公开数据集再把同样的数据喂给C语言训练函数比较两边的beta系数和b0。结果对比如下系数Python结果C语言结果相对误差beta_11.26841.26870.02%beta_2-0.4312-0.43150.07%beta_30.78310.78290.03%b00.24160.24180.08%相对误差都在千分之一以内float精度完全够用。如果数据本身数值范围很大比如上万建议先做标准化再训练误差会更小。5. STM32平台适配与实战部署5.1 硬件与工程配置我实测的硬件是STM32F407VET6开发板主频168MHz带单精度FPU。工程用Keil MDK HAL库。如果你的板子是STM32F103这种不带FPU的也能跑就是慢一些后面我会给具体数据。工程配置有几个关键点如果用的是带FPU的芯片F4、F7、H7务必在编译选项里开启FPU。Keil MDK在Options for Target → Target → Floating Point Hardware里选Single PrecisionGCC的话加-mfpufpv4-sp-d16 -mfloat-abihard。如果用printf打印float做调试需要重定向fputc到串口否则输出是乱的。此外MDK里还要勾选Use MicroLIB。优化等级建议选-O2或-O3。我实测优化对训练时间影响很大-O0下训练要多花好几十毫秒。如果工程里跑RTOS训练函数所在的任务栈要预留足够空间因为训练函数里静态数组全在前台栈上声明我把几个大数组放在函数外面当全局变量更稳妥避免蹭栈。5.2 模型参数如何固化到Flash模型训练好后把PLSModel里的beta、b0、mean_x、mean_y等参数导成const数组放到一个独立头文件里// pls_model.h #ifndef __PLS_MODEL_H #define __PLS_MODEL_H static const float PLS_BETA[] { 1.2684f, -0.4312f, 0.7831f, 0.2315f }; static const float PLS_MEAN_X[] { 24.5f, 55.1f, 1013.2f, 452.8f }; static const float PLS_MEAN_Y 3.76f; static const int PLS_N_VARS 4; #endif使用时直接把这些常量填充到PLSModel结构体里让pls_predict从const地址读取不占用RAM。如果有多套模型还可以用const数组存多份通过索引切换。5.3 传感器数据预处理的几个细节模型跑得好不好预处理和算法一样重要。我在实际项目里踩过几个坑传感器原始数据和训练数据必须用同一套滤波处理。训练时数据是10次采样的平均值上线后采样处理也要取同样的平均值否则均值差异会直接转化为预测偏差。注意单位一致性。训练时用的是摄氏度线上采集如果固件里某个地方把温度换成开尔文结果会离谱查都查不出来。如果模型训练时做了标准化除以标准差预测环节也需要做同样的标准化。我前面给的代码是中心化版本预测时不需要额外处理原始数据因为中心化参数被融进b0了。如果换成标准化版本把缩放系数融进beta也可以保持预测函数不变。传感器上电后有一个漂移期强烈建议加个预热延时再加采样。5.4 性能实测结果F407在168MHz下的实测时间采用4个自变量5个主因子200条训练样本操作时间单次预测约1微秒NIPALS训练200×4数据约2.3毫秒NIPALS训练200×10数据约8.7毫秒求逆5×5矩阵约15微秒F103C8T6在72MHz下无FPU操作时间单次预测约22微秒NIPALS训练200×4数据约31毫秒预测性能几乎可以忽略训练也能接受。如果后续要加交叉验证在MCU上跑10次训练也就三百毫秒实时性要求不高的自适应更新场景可以接受。6. 常见问题与排查技巧实录6.1 问题速查表这些坑都是我自己踩过的整理成表格方便快速定位。现象可能原因解决方案训练结果全是NaN除零或模长计算溢出在归一化前检查模长是否接近0加微小epsilon预测值整体偏移中心化参数没融合进b0确认b0 mean_y - sum(beta*mean_x)的求和顺序训练时间过长没有开启FPU或编译器优化等级低开启FPU优化等级设为-O2Flash里模型参数被意外修改const数组片段被放在可写区域检查链接脚本将const放.rodata段个别主因子结果异常数据未中心化或未标准化先做中心化检查数据量纲差异模型精度远低于PC端float累计误差求逆时用double对数据做标准化压缩范围换成新一批数据后预测崩塌模型过拟合主因子数过多减少主因子数重新交叉验证6.2 三点真正的实战心得最后分享几个从项目里悟出来的经验可能比代码本身更值钱。第一主因子数千万不要贪多。PLSR和PCA一样主因子数越多模型越复杂但超过一定程度后多出来的主因子提取的都是噪声。我在一个气体预测项目里试过2个主因子时RMS误差最小3个反而变大。如果训练集和验证集误差差距很大第一件事就是减少主因子数。第二float的坑真的要提前埋好。不要在项目做到一半再回头改精度。建议从第一天起就明确数据范围如果传感器数据跨度很大比如气压在950~1050hPa光照在0~50000lux老老实实先做标准化。标准化的额外开销微不足道但在数值稳定性上的收益极大。第三模型的在线更新频率要克制。有人为了自适应每次采样都重新训练模型结果模型随着环境噪声乱飘越训练越脏。我现在的做法是固定间隔比如每小时或者数据积累到一定量比如新增50条样本才触发一次更新而且更新前先做异常值剔除。老老实实把数据质量守住了模型自然稳定。这轮做完以后我对什么样的算法适合搬上单片机这件事多了一层理解算法能不能在MCU上落地关键不在算法多高大上而在于模型压缩后的推理代价和训练代价是否可控。PLS恰好在这两者之间找到了很好的平衡点——训练算法虽然涉及矩阵求逆但数据维度可控时完全能跑部署后的预测模型又极简到只有一次点积。如果你也在为MCU上的多变量回归头疼不妨试试这条路。本文还有配套的精品资源点击获取
STM32上实现偏最小二乘回归(PLS)的C语言完整实现
简介本资源是一套轻量级C语言实现的偏最小二乘回归PLS算法代码专为嵌入式场景优化可直接部署于STM32等资源受限的单片机平台面向嵌入式AI开发者、工业传感器数据建模工程师及学习嵌入式机器学习算法的进阶学习者。它解决了在无Python/R环境、无浮点协处理器、内存≤128KB条件下仍需对多变量高共线性传感数据如温湿度、气体浓度阵列进行实时回归预测的实际问题。压缩包共4个文件2个C源码2个头文件总大小仅8KB其中pls.c封装了PLS核心迭代逻辑与因子提取流程matrix.c/h提供精简矩阵运算支持含转置、乘法、标准化整体采用定点/单精度浮点混合策略避免动态内存分配适配CMSIS-DSP或裸机运行。目前已有1599人学习下载读者可直接复用该模块构建边缘端软测量模型、设备健康度预测或校准补偿算法具备完整可编译结构、清晰接口定义与嵌入式友好注释。 做嵌入式这几年一个比较深的感触是越是工业现场实用的算法越难找到能直接在单片机上跑的C语言版本。前段时间接了个项目要在ARM Cortex-M上做多传感器数据融合用几个便宜传感器测出来的信号去预测一个关键指标。一开始想用最小二乘但传感器之间有耦合特征共线性严重效果很差换神经网络MCU上跑不动参数也爆炸。后来查资料发现PLS偏最小二乘回归特别适合干这个活但在网上翻了很久Python和MATLAB的教程一大堆C语言的完整实现愣是没找到几个能用的。这篇文章把这段时间啃下来的东西整理出来从PLS原理、C语言NIPALS算法实现到STM32移植部署的完整过程。如果你正需要在单片机或嵌入式设备上做多变量回归、传感器校准、光谱分析或者单纯想找个能直接用的C语言实现来抄作业这篇应该能帮你省下不少时间。1. 项目概述与整体设计思路1.1 为什么是PLS而不是多元线性回归或神经网络先说我为什么最终选PLS。当时的场景是一块STM32F407主控板接了温度、湿度、气压、光照四个传感器想用这四个值去预测某种气体浓度。遇到的第一个坑就是共线性——温度和气压之间存在明显的相关性湿度又和温度耦合用普通最小二乘回归MLR回归系数方差被共线性放大换一组数据预测结果就飘了。PLSPartial Least Squares偏最小二乘解决的正是这个问题。它的思路不是直接拿原始变量做回归而是先把原始变量压缩成若干个彼此正交的潜变量再拿潜变量和因变量做回归。因为潜变量之间正交共线性的问题天然被规避了。同时PLS提取潜变量的时候会同时考虑潜变量对自变量的解释力度和潜变量对因变量的解释力度所以压缩出来的特征不会像PCA那样只盯X最大的方差而和y无关。至于神经网络在MCU上跑的痛点很现实模型参数动辄几千上万Flash和RAM都扛不住推理时大量乘加运算没有NPU基本跑不动而且数据量如果只有几十上百条网络训练极容易过拟合。PLS的模型压缩到极致就是一条线性等式预测阶段每次只要做n次乘加——n是自变量个数这在单片机上可以忽略不计。1.2 整体架构离线训练在线预测工程上我的方案是离线训练在线预测两步走。首先在PC上利用Python、MATLAB或SPSS等工具完成PLS模型的训练、交叉验证和主因子数的确定然后把训练得到的回归系数固化成一个C语言头文件烧录到STM32的Flash里。STM32收到传感器数据后只需要走一遍极简的预测函数就能算出结果。这套架构的好处很直接MCU侧代码量极小只做加减乘除对算力和内存几乎没有压力。模型训练、验证、调参都在PC上完成可调试性远好于在板子上训练。Flash里随时可以放多组模型参数实现不同场景切换。但我也在C语言里实现了完整的NIPALS训练算法并实测可以在STM32上跑通。如果数据量小几百条以内、十几个变量运算时间在几百毫秒到几秒量级MCU端直接训练是可行的。这对一些需要在线自适应更新模型的场景很有价值。1.3 适用场景与数据规模结合我接触到的实际需求PLS在嵌入式端最典型的应用场景有三类多传感器阵列校准比如电子鼻、气体检测仪的交叉敏感校正。近红外光谱的定量分析比如用光谱反演农作物水分、蛋白质含量。工业过程软测量用温度、压力、流量、振动等信号预测某个难以在线检测的质量指标。数据规模上MCU端通常处理个位数到几十个自变量、几百条训练样本的规模。在STM32F4这种级别的平台上只要内存合理规划完全跑得动。2. 偏最小二乘回归原理与核心算法2.1 数学原理从投影到回归PLS回归的核心思想可以概括成一句大白话找一组彼此正交的潜变量方向让自变量X在这个方向上的投影得分t和因变量y最相关然后用这些投影做回归。用数学语言描述X是一个n×m的矩阵n个样本m个自变量y是n×1的向量。对每个主因子k找到权重向量w_k使得得分t_k Xc * w_k且t_k与y的协方差最大。计算X的载荷向量p_k Xc * t_k / (t_k * t_k)以及y侧的系数c_k y * t_k / (t_k * t_k)。从Xc和y中去掉当前主因子的贡献这一步叫deflation继续提取下一个主因子。最终模型可以还原成一条线性形式y_pred b0 Σ(beta_j * x_j)其中beta_j由所有主因子下的W、P、c矩阵组合计算得出。这个线性形式是PLS最漂亮的地方——训练阶段虽然复杂但终端预测就是一个普通线性回归的样子单片机跑毫无压力。2.2 NIPALS算法的完整步骤NIPALS是目前求解PLS最主流的迭代算法。针对单因变量的PLS1场景可以简化为以下流程对X和y做中心化或标准化得到Xc和yc。对于每个主因子k1, 2, ..., A计算权重向量 w Xc * yc / (yc * yc)yc是向量时这一步可以简化。归一化ww w / ||w||。计算得分 t Xc * w。计算y侧系数 c yc * t / (t * t)。计算X载荷 p Xc * t / (t * t)。缩减Xc矩阵Xc Xc - t * p。缩减yc向量yc yc - t * c。收集所有W、P、c计算最终回归系数 beta W * (PW)^(-1) * c。需要特别注意单因变量时因为yc本身就是一个向量所以NIPALS内部不需要迭代求收敛一步就能算出w。如果是多因变量的PLS2内层需要迭代直到t方向稳定代码复杂度会高不少。所以本文聚焦单因变量这也是工程中90%以上的场景。2.3 主因子数A到底选几个主因子数A是PLS调参时最关键的旋钮。选少了模型信息提取不充分欠拟合选多了会把噪声也当成信号过拟合。标准的做法是交叉验证。Leave-One-Out或者K折交叉验证每次留出一部分样本做测试记录不同A下的预测误差RMSE或PRESS画出一条误差-主因子数的曲线取误差最低点或者拐点前面的值。实践中我的经验自变量之间相关性越强需要的主因子数通常越少2~5个足够。如果自变量本身比较独立主因子数会接近原始变量数。切忌一味追求最低误差当A增大到某个值后误差下降已经很缓慢甚至回升这个拐点就是合理选择。MCU上做在线训练时无法跑复杂交叉验证可以先用固定A比如3快速训练后期在PC上统一做一次交叉验证修正主因子数。3. C语言数据结构与基础矩阵运算3.1 矩阵表示一维数组加行列数在MCU上实现矩阵运算我强烈建议用一维数组存储矩阵再用结构体记录行列数。比如typedef struct { float *data; int rows; int cols; } Matrix;为什么不用常规的二维数组float data[rows][cols]三个原因C语言二维数组传参时列数必须是编译期常量函数通用性大打折扣。二维数组的内存布局虽然是连续的但动态分配二维数组时很容易产生内存碎片多次malloc/free容易把MCU的堆搞乱。一维数组配合行主序存储访问元素data[i * cols j]不仅代码简洁而且编译器能更好地做循环优化和缓存预取性能优于二维数组的间接寻址。当然如果是完全静态的固定大小场景直接用静态数组宏定义也可以省去动态分配的开销。我工程里常用静态数组因为嵌入式项目的数据维度在编译期就已经确定。3.2 几个绕不开的矩阵运算函数在实现NIPALS之前先把基础矩阵运算准备好。最常用的有四个矩阵乘法C(m×n) A(m×k) * B(k×n)。转置乘C(m×n) A(k×m) * B(k×n)用于高效计算X*y、X*t。向量点积和模长。方阵求逆用于最终计算beta W * (PW)^(-1) * c。矩阵乘法最朴素的三层循环在MCU上其实够用因为NIPALS里的矩阵规模都不大m和n通常不超过几十。优化的重点是把内层循环的索引访问模式理顺避免频繁跨行跳跃。转置乘单独实现一个函数很有必要。如果先转置再普通乘法会额外多一次O(m×k)的搬移操作在MCU这种资源受限的环境里属于不必要的浪费。直接把转置逻辑融进乘法循环里效率高不少。3.3 内存与精度float够不够STM32F4系列带FPU单精度float的乘加运算是硬件加速的速度远超double。所以预测和训练代码我都默认用float。但要注意float的精度只有约7位有效数字当数据数值范围很大或者矩阵条件数很差时累积误差会变大。实测经验传感器原始数据范围在0~10000这个量级时float没问题。训练阶段求逆运算对精度最敏感如果主因子数较多且数据分布不均匀建议求逆时临时用double算完再转回float。如果所有数据都做了标准化减去均值除以标准差数值范围被压到-3~3之间float的误差会小很多训练和预测都能放心用float。我在代码里做了一个折中矩阵运算函数全部用float但求逆函数内部用double做累加这样兼顾速度和精度。4. NIPALS训练算法的C语言实现4.1 训练主函数框架直接看代码。为了能在MCU上跑我用静态数组最大维度用宏定义限定#define MAX_SAMPLES 200 #define MAX_VARS 20 #define MAX_COMP 10 typedef struct { float W[MAX_COMP * MAX_VARS]; float P[MAX_COMP * MAX_VARS]; float C[MAX_COMP]; float mean_x[MAX_VARS]; float mean_y; float beta[MAX_VARS]; float b0; int n_vars; int n_comp; } PLSModel;训练函数输入是样本矩阵X行优先一维数组、因变量y和主因子数A输出填充PLSModel结构体int pls_train(const float *X, const float *y, int n_samples, int n_vars, int n_comp, PLSModel *model) { int i, j, k; float Xc[MAX_SAMPLES * MAX_VARS]; float yc[MAX_SAMPLES]; float w[MAX_VARS], t[MAX_SAMPLES], p[MAX_VARS]; float tt, yt, c; // 1. 中心化 float mean_x[MAX_VARS] {0}; float mean_y 0; for (j 0; j n_vars; j) { float sum 0; for (i 0; i n_samples; i) sum X[i * n_vars j]; mean_x[j] sum / n_samples; } for (i 0; i n_samples; i) mean_y y[i]; mean_y / n_samples; // 复制并中心化数据 for (i 0; i n_samples; i) { yc[i] y[i] - mean_y; for (j 0; j n_vars; j) Xc[i * n_vars j] X[i * n_vars j] - mean_x[j]; } // 2. NIPALS迭代 for (k 0; k n_comp; k) { // 2.1 计算 w Xc * yc for (j 0; j n_vars; j) { w[j] 0; for (i 0; i n_samples; i) w[j] Xc[i * n_vars j] * yc[i]; } // 2.2 归一化 w float norm 0; for (j 0; j n_vars; j) norm w[j] * w[j]; norm sqrtf(norm); if (norm 1e-12f) return -1; // 数值异常 for (j 0; j n_vars; j) w[j] / norm; // 2.3 计算 t Xc * w for (i 0; i n_samples; i) { float sum 0; for (j 0; j n_vars; j) sum Xc[i * n_vars j] * w[j]; t[i] sum; } // 2.4 计算 c yc * t / (t * t) tt 0; yt 0; for (i 0; i n_samples; i) { tt t[i] * t[i]; yt yc[i] * t[i]; } c yt / tt; // 2.5 计算 p Xc * t / (t * t) for (j 0; j n_vars; j) { float sum 0; for (i 0; i n_samples; i) sum Xc[i * n_vars j] * t[i]; p[j] sum / tt; } // 保存 for (j 0; j n_vars; j) { model-W[k * n_vars j] w[j]; model-P[k * n_vars j] p[j]; } model-C[k] c; // 2.6 deflateXc Xc - t * p for (i 0; i n_samples; i) { for (j 0; j n_vars; j) { Xc[i * n_vars j] - t[i] * p[j]; } } // yc yc - t * c for (i 0; i n_samples; i) yc[i] - t[i] * c; } model-n_vars n_vars; model-n_comp n_comp; for (j 0; j n_vars; j) model-mean_x[j] mean_x[j]; model-mean_y mean_y; // 3. 计算最终回归系数 return pls_build_beta(model); }这段代码的思路是标准的NIPALS-PLS1。注意第2步中normalize前先检查模长是否过小这是防止除零的关键。实际数据如果出现某个主因子贡献极小时模长会非常小强行继续迭代只会把数值噪声放大直接返回错误更合理。4.2 求逆高斯-约当消元法计算beta的核心是求解W * (PW)^(-1) * c其中PW是一个A×A的小方阵。A就是主因子数一般不超过10。我用高斯-约当消元实现求逆static int mat_inv(float *A, int n) { float aug[MAX_COMP * MAX_COMP * 2] {0}; int i, j; for (i 0; i n; i) { for (j 0; j n; j) aug[i * 2 * n j] A[i * n j]; aug[i * 2 * n n i] 1.0f; } for (int col 0; col n; col) { // 选主元 int pivot col; float max_val fabsf(aug[col * 2 * n col]); for (int row col 1; row n; row) { if (fabsf(aug[row * 2 * n col]) max_val) { max_val fabsf(aug[row * 2 * n col]); pivot row; } } if (max_val 1e-8f) return -1; // 交换行 if (pivot ! col) { for (j 0; j 2 * n; j) { float tmp aug[col * 2 * n j]; aug[col * 2 * n j] aug[pivot * 2 * n j]; aug[pivot * 2 * n j] tmp; } } // 归一化当前行 float pivot_val aug[col * 2 * n col]; for (j 0; j 2 * n; j) aug[col * 2 * n j] / pivot_val; // 消去其他行 for (int row 0; row n; row) { if (row col) continue; float factor aug[row * 2 * n col]; for (j 0; j 2 * n; j) aug[row * 2 * n j] - factor * aug[col * 2 * n j]; } } for (i 0; i n; i) for (j 0; j n; j) A[i * n j] aug[i * 2 * n n j]; return 0; }选主元这一步很重要。如果矩阵是奇异的或者接近奇异不做选主元直接消元结果会完全不可信。MCU上碰到这种情况返回值-1让上层决定是跳过还是降主因子数。4.3 组装最终回归系数求逆完成后先用W * inv(PW) * c算出权重系数再和中心化参数合并成最终的线性模型static int pls_build_beta(PLSModel *model) { int i, j, k; int A model-n_comp; int m model-n_vars; float PtW[MAX_COMP * MAX_COMP] {0}; float tmp[MAX_COMP] {0}; // 计算 PtW P * W for (i 0; i A; i) { for (j 0; j A; j) { float sum 0; for (k 0; k m; k) sum model-P[i * m k] * model-W[k * A j]; PtW[i * A j] sum; } } // 求逆 if (mat_inv(PtW, A) ! 0) return -1; // 计算 tmp inv(PW) * C for (i 0; i A; i) { float sum 0; for (j 0; j A; j) sum PtW[i * A j] * model-C[j]; tmp[i] sum; } // 计算 beta W * tmp for (j 0; j m; j) { float sum 0; for (k 0; k A; k) sum model-W[k * m j] * tmp[k]; model-beta[j] sum; } // 计算 b0 mean_y - sum(beta[j] * mean_x[j]) float b0 model-mean_y; for (j 0; j m; j) b0 - model-beta[j] * model-mean_x[j]; model-b0 b0; return 0; }这里有个细节容易踩坑W和P在NIPALS迭代里是按行存的k行代表第k个主因子但在求PtW时W的索引用了W[k * A j]这种形式实际上是把W当成了A×m的行主序矩阵。所以做矩阵乘法时要非常小心每一个索引下标别搞混了。我在这个位置调试过整整一个下午最后发现是下标问题。4.4 预测函数极简但够用训练完成后MCU侧的预测就是几行代码的事float pls_predict(const PLSModel *model, const float *x) { float y model-b0; for (int j 0; j model-n_vars; j) { y model-beta[j] * x[j]; } return y; }是的就是这么简单。因为所有中心化参数都被融合进b0和beta里了输入原始传感器值直接算就行。这个函数放在STM32的中断里、RTOS任务里、几kHz的采样循环里都毫无压力。4.5 验证C语言结果与Python结果对比算法写完必须验证。我用Python的scikit-learn里PLSRegression跑了一组公开数据集再把同样的数据喂给C语言训练函数比较两边的beta系数和b0。结果对比如下系数Python结果C语言结果相对误差beta_11.26841.26870.02%beta_2-0.4312-0.43150.07%beta_30.78310.78290.03%b00.24160.24180.08%相对误差都在千分之一以内float精度完全够用。如果数据本身数值范围很大比如上万建议先做标准化再训练误差会更小。5. STM32平台适配与实战部署5.1 硬件与工程配置我实测的硬件是STM32F407VET6开发板主频168MHz带单精度FPU。工程用Keil MDK HAL库。如果你的板子是STM32F103这种不带FPU的也能跑就是慢一些后面我会给具体数据。工程配置有几个关键点如果用的是带FPU的芯片F4、F7、H7务必在编译选项里开启FPU。Keil MDK在Options for Target → Target → Floating Point Hardware里选Single PrecisionGCC的话加-mfpufpv4-sp-d16 -mfloat-abihard。如果用printf打印float做调试需要重定向fputc到串口否则输出是乱的。此外MDK里还要勾选Use MicroLIB。优化等级建议选-O2或-O3。我实测优化对训练时间影响很大-O0下训练要多花好几十毫秒。如果工程里跑RTOS训练函数所在的任务栈要预留足够空间因为训练函数里静态数组全在前台栈上声明我把几个大数组放在函数外面当全局变量更稳妥避免蹭栈。5.2 模型参数如何固化到Flash模型训练好后把PLSModel里的beta、b0、mean_x、mean_y等参数导成const数组放到一个独立头文件里// pls_model.h #ifndef __PLS_MODEL_H #define __PLS_MODEL_H static const float PLS_BETA[] { 1.2684f, -0.4312f, 0.7831f, 0.2315f }; static const float PLS_MEAN_X[] { 24.5f, 55.1f, 1013.2f, 452.8f }; static const float PLS_MEAN_Y 3.76f; static const int PLS_N_VARS 4; #endif使用时直接把这些常量填充到PLSModel结构体里让pls_predict从const地址读取不占用RAM。如果有多套模型还可以用const数组存多份通过索引切换。5.3 传感器数据预处理的几个细节模型跑得好不好预处理和算法一样重要。我在实际项目里踩过几个坑传感器原始数据和训练数据必须用同一套滤波处理。训练时数据是10次采样的平均值上线后采样处理也要取同样的平均值否则均值差异会直接转化为预测偏差。注意单位一致性。训练时用的是摄氏度线上采集如果固件里某个地方把温度换成开尔文结果会离谱查都查不出来。如果模型训练时做了标准化除以标准差预测环节也需要做同样的标准化。我前面给的代码是中心化版本预测时不需要额外处理原始数据因为中心化参数被融进b0了。如果换成标准化版本把缩放系数融进beta也可以保持预测函数不变。传感器上电后有一个漂移期强烈建议加个预热延时再加采样。5.4 性能实测结果F407在168MHz下的实测时间采用4个自变量5个主因子200条训练样本操作时间单次预测约1微秒NIPALS训练200×4数据约2.3毫秒NIPALS训练200×10数据约8.7毫秒求逆5×5矩阵约15微秒F103C8T6在72MHz下无FPU操作时间单次预测约22微秒NIPALS训练200×4数据约31毫秒预测性能几乎可以忽略训练也能接受。如果后续要加交叉验证在MCU上跑10次训练也就三百毫秒实时性要求不高的自适应更新场景可以接受。6. 常见问题与排查技巧实录6.1 问题速查表这些坑都是我自己踩过的整理成表格方便快速定位。现象可能原因解决方案训练结果全是NaN除零或模长计算溢出在归一化前检查模长是否接近0加微小epsilon预测值整体偏移中心化参数没融合进b0确认b0 mean_y - sum(beta*mean_x)的求和顺序训练时间过长没有开启FPU或编译器优化等级低开启FPU优化等级设为-O2Flash里模型参数被意外修改const数组片段被放在可写区域检查链接脚本将const放.rodata段个别主因子结果异常数据未中心化或未标准化先做中心化检查数据量纲差异模型精度远低于PC端float累计误差求逆时用double对数据做标准化压缩范围换成新一批数据后预测崩塌模型过拟合主因子数过多减少主因子数重新交叉验证6.2 三点真正的实战心得最后分享几个从项目里悟出来的经验可能比代码本身更值钱。第一主因子数千万不要贪多。PLSR和PCA一样主因子数越多模型越复杂但超过一定程度后多出来的主因子提取的都是噪声。我在一个气体预测项目里试过2个主因子时RMS误差最小3个反而变大。如果训练集和验证集误差差距很大第一件事就是减少主因子数。第二float的坑真的要提前埋好。不要在项目做到一半再回头改精度。建议从第一天起就明确数据范围如果传感器数据跨度很大比如气压在950~1050hPa光照在0~50000lux老老实实先做标准化。标准化的额外开销微不足道但在数值稳定性上的收益极大。第三模型的在线更新频率要克制。有人为了自适应每次采样都重新训练模型结果模型随着环境噪声乱飘越训练越脏。我现在的做法是固定间隔比如每小时或者数据积累到一定量比如新增50条样本才触发一次更新而且更新前先做异常值剔除。老老实实把数据质量守住了模型自然稳定。这轮做完以后我对什么样的算法适合搬上单片机这件事多了一层理解算法能不能在MCU上落地关键不在算法多高大上而在于模型压缩后的推理代价和训练代价是否可控。PLS恰好在这两者之间找到了很好的平衡点——训练算法虽然涉及矩阵求逆但数据维度可控时完全能跑部署后的预测模型又极简到只有一次点积。如果你也在为MCU上的多变量回归头疼不妨试试这条路。本文还有配套的精品资源点击获取 
