1. 为什么用C语言手搓多元线性回归很多人第一次接触机器学习都是从Python的scikit-learn或者R语言的lm()函数开始的。一行代码就能拟合出模型方便是方便但模型内部到底在算什么很多人其实说不清楚。我在带新人的时候经常发现一个现象调包调得很熟练但一旦遇到矩阵奇异、特征量纲差异巨大、迭代不收敛这些问题就完全不知道从哪里下手排查。用C语言从零实现多元线性回归最大的价值不在于用C更快——对于小规模数据Python的numpy底层也是C和Fortran速度差距没那么夸张。真正的价值在于你被迫把每一步数学运算都写出来被迫管理内存被迫处理数值稳定性问题。这个过程会让你对最小二乘法的理解从知道公式变成知道每一步在干什么。多元线性回归的模型形式很简洁y β₀ β₁x₁ β₂x₂ ... βₚxₚ ε其中y是因变量x₁到xₚ是p个自变量β₀是截距β₁到βₚ是各自变量的系数ε是误差项。所谓拟合就是找到一组β使得预测值和真实值之间的差距尽可能小。最常用的准则是最小化残差平方和也就是最小二乘法。用矩阵形式写出来就是Y Xβ ε其中X是n×(p1)的设计矩阵第一列全是1对应截距项Y是n×1的观测向量β是(p1)×1的系数向量。最小二乘解的正规方程是β (XᵀX)⁻¹XᵀY这个公式看起来简单但直接求逆在数值上很不稳定。实际实现中我们通常用高斯消元法或者QR分解来解这个线性方程组而不是显式地求逆。这一点在后面会详细展开。这篇文章适合谁看如果你已经学过线性代数基础知道矩阵乘法怎么算会写基本的C语言程序但对怎么把数学公式变成能跑的代码还没有完整经验那这篇内容就是为你准备的。我会从数据结构设计开始一步步把多元线性回归的核心计算写出来包括设计矩阵的构建、正规方程的组装、线性方程组的求解以及一些实际编码中容易踩的坑。需要说明的是这是系列的第一篇重点放在最小二乘法的直接求解上。梯度下降法、正则化、模型评估这些内容会放在后续篇章里展开。先把最核心的数学到代码的映射搞清楚后面的扩展才有根基。2. 数据结构设计与内存布局的取舍2.1 用一维数组还是二维数组存矩阵在C语言里表示矩阵最直觉的做法是二维数组double X[n][p1]。但这里有个问题n和p在编译时往往不确定而C99的变长数组虽然支持运行时确定大小但在一些编译器上支持并不完整而且大矩阵放在栈上容易爆栈。我的选择是用一维数组模拟二维矩阵配合手动索引计算。具体来说一个n行m列的矩阵用一个长度为n*m的double数组存储访问第i行第j列的元素用mat[i*m j]。这样做的好处是内存连续缓存友好而且动态分配方便。typedef struct { double *data; // 一维数组按行优先存储 int rows; // 行数 int cols; // 列数 } Matrix;行优先存储意味着同一行的元素在内存中是连续的。为什么选行优先而不是列优先因为在组装正规方程XᵀX的时候我们主要做的是行与行的点积运算行优先的布局能让每次点积访问的内存尽量连续减少缓存未命中。2.2 矩阵的分配、释放与访问封装动态分配一个矩阵需要两步先分配结构体本身再分配数据区。为了避免忘记释放导致内存泄漏我把分配和释放都封装成函数Matrix* matrix_create(int rows, int cols) { Matrix *m (Matrix*)malloc(sizeof(Matrix)); if (!m) return NULL; m-rows rows; m-cols cols; m-data (double*)calloc(rows * cols, sizeof(double)); if (!m-data) { free(m); return NULL; } return m; } void matrix_free(Matrix *m) { if (m) { free(m-data); free(m); } }这里用calloc而不是malloc是因为calloc会把内存清零。对于设计矩阵和正规方程矩阵来说清零能避免很多因为未初始化内存导致的诡异bug。我踩过这个坑有一次用malloc分配了XᵀX矩阵忘了初始化结果累加的时候把垃圾值也算进去了拟合出来的系数完全离谱排查了半天才发现是初始化的问题。访问元素的宏定义#define MAT(m, i, j) ((m)-data[(i) * (m)-cols (j)])用宏而不是函数是为了避免函数调用开销。虽然现代编译器大概率会内联但宏更直观地表达了这就是个索引计算的意图。2.3 设计矩阵的构建截距列的处理设计矩阵X的第一列全是1这是为了把截距β₀统一进矩阵运算。假设原始数据有n个样本每个样本有p个特征那么X就是n行(p1)列。Matrix* build_design_matrix(double *features, double *target, int n, int p) { Matrix *X matrix_create(n, p 1); if (!X) return NULL; for (int i 0; i n; i) { MAT(X, i, 0) 1.0; // 截距列 for (int j 0; j p; j) { MAT(X, i, j 1) features[i * p j]; } } return X; }原始特征数据我用一维数组按行优先存储第i个样本的第j个特征在features[i*p j]。这样从原始数据到设计矩阵的映射很直接不需要额外的拷贝。注意截距列放在第一列还是最后一列会影响后续系数向量的读取顺序。我习惯放在第一列这样β[0]就是截距β[1]到β[p]对应各个特征读起来比较自然。如果你选择放在最后一列记得在输出结果时对应调整。2.4 目标向量的存储选择目标向量Y是n×1的理论上也可以用Matrix结构体存。但为了简化我直接用double*数组。因为Y在正规方程里只参与XᵀY的计算不需要矩阵的行列信息。double* y (double*)malloc(n * sizeof(double));这样做减少了结构体嵌套的复杂度代码更清爽。代价是失去了统一的矩阵接口但对于Y这种一维数据来说直接操作数组反而更方便。3. 正规方程的组装XᵀX与XᵀY的计算3.1 为什么不用公式直接求逆前面提到β (XᵀX)⁻¹XᵀY理论上只要算出XᵀX求个逆再乘XᵀY就完事了。但实际编码中显式求逆有几个致命问题。第一求逆的运算量是O(p³)而解线性方程组的运算量也是O(p³)但常数因子更小。第二也是更关键的当XᵀX接近奇异时求逆会放大数值误差。举个例子如果两个特征高度相关XᵀX的某两列几乎线性相关行列式接近零求逆得到的矩阵元素会非常大微小的输入扰动会导致系数剧烈变化。正确的做法是把问题转化为解线性方程组Aβ b其中A XᵀXb XᵀY然后用高斯消元法求解。这样不需要显式计算逆矩阵数值稳定性好得多。3.2 XᵀX的对称性与计算优化XᵀX是一个(p1)×(p1)的方阵而且是对称矩阵。第i行第j列的元素等于X第i列和第j列的点积(XᵀX)ᵢⱼ Σₖ Xₖᵢ · Xₖⱼ因为对称我们只需要计算上三角部分然后镜像到下三角能省将近一半的计算量。对于p比较大的情况这个优化很可观。Matrix* compute_xtx(Matrix *X) { int n X-rows; int m X-cols; Matrix *XtX matrix_create(m, m); if (!XtX) return NULL; for (int i 0; i m; i) { for (int j i; j m; j) { double sum 0.0; for (int k 0; k n; k) { sum MAT(X, k, i) * MAT(X, k, j); } MAT(XtX, i, j) sum; MAT(XtX, j, i) sum; // 利用对称性 } } return XtX; }内层循环累加的时候我特意把MAT(X, k, i)和MAT(X, k, j)写在一起。因为X是行优先存储固定k时第i列和第j列的元素在内存中相隔cols个位置。如果p不大这两次访问大概率在同一个缓存行里效率不错。3.3 XᵀY的计算与内存访问模式XᵀY是一个(p1)×1的向量第i个元素是X第i列和Y的点积(XᵀY)ᵢ Σₖ Xₖᵢ · Yₖdouble* compute_xty(Matrix *X, double *y) { int n X-rows; int m X-cols; double *xty (double*)calloc(m, sizeof(double)); if (!xty) return NULL; for (int i 0; i m; i) { double sum 0.0; for (int k 0; k n; k) { sum MAT(X, k, i) * y[k]; } xty[i] sum; } return xty; }这里的内存访问模式是外层遍历列i内层遍历行k。对于固定的i访问的是X的第i列在行优先存储下这些元素之间间隔cols个double。这种跨步访问对缓存不太友好但因为只做一次影响可以接受。如果追求极致性能可以先转置X再算但那样需要额外的内存和转置时间对于中小规模数据不划算。3.4 数值量级差异带来的隐患在实际数据中不同特征的量纲可能差异巨大。比如一个特征是年龄20到60另一个特征是年收入几万到几百万。直接计算XᵀX时收入对应的对角线元素会是年龄对应元素的10¹⁰倍以上。这种量级差异会导致矩阵条件数极大高斯消元时出现大数吃小数的现象系数估计严重失真。解决办法有两个一是标准化特征把每个特征减去均值再除以标准差让所有特征在同一量级二是用带主元的高斯消元通过行交换把绝对值最大的元素换到主元位置。我建议两个都做标准化是预处理阶段的事主元选取是求解器内部的事双保险。标准化的代码很简单但要注意保存均值和标准差因为预测新样本时需要用同样的参数做变换void standardize(double *features, int n, int p, double *mean, double *std) { for (int j 0; j p; j) { double sum 0.0; for (int i 0; i n; i) sum features[i * p j]; mean[j] sum / n; double sq_sum 0.0; for (int i 0; i n; i) { double diff features[i * p j] - mean[j]; sq_sum diff * diff; } std[j] sqrt(sq_sum / n); if (std[j] 1e-12) std[j] 1.0; // 防止除零 for (int i 0; i n; i) { features[i * p j] (features[i * p j] - mean[j]) / std[j]; } } }提示标准化之后截距项的含义会变化。原始空间中截距是当所有特征为0时的预测值标准化空间中截距是当所有特征取均值时的预测值。如果你需要原始空间的截距可以在标准化之前先把目标向量也中心化或者拟合完之后再做逆变换。这个细节在后续篇章讲预测的时候会再展开。4. 高斯消元求解线性方程组4.1 朴素高斯消元的实现与局限有了A XᵀX和b XᵀY问题变成解Aβ b。最基础的高斯消元分两步前向消元把A变成上三角矩阵回代求解β。int gaussian_elimination(double *A, double *b, double *x, int n) { // 前向消元 for (int k 0; k n - 1; k) { if (fabs(A[k * n k]) 1e-12) return -1; // 主元过小 for (int i k 1; i n; i) { double factor A[i * n k] / A[k * n k]; for (int j k; j n; j) { A[i * n j] - factor * A[k * n j]; } b[i] - factor * b[k]; } } // 回代 for (int i n - 1; i 0; i--) { double sum b[i]; for (int j i 1; j n; j) { sum - A[i * n j] * x[j]; } if (fabs(A[i * n i]) 1e-12) return -1; x[i] sum / A[i * n i]; } return 0; }这个实现的问题在于如果主元A[k][k]恰好接近零即使矩阵本身非奇异消元也会失败或产生巨大误差。比如矩阵[[0, 1], [1, 0]]第一主元是0朴素消元直接卡住但交换两行就能正常求解。4.2 列主元选取让消元更稳健列主元的思想很简单在第k步消元之前从第k行到第n-1行中找到第k列绝对值最大的元素把它所在的行和第k行交换。这样主元的绝对值尽可能大消元时的乘数因子就不会爆炸。int gaussian_elimination_pivot(double *A, double *b, double *x, int n) { for (int k 0; k n - 1; k) { // 找列主元 int max_row k; double max_val fabs(A[k * n k]); for (int i k 1; i n; i) { if (fabs(A[i * n k]) max_val) { max_val fabs(A[i * n k]); max_row i; } } if (max_val 1e-12) return -1; // 矩阵奇异 // 交换行 if (max_row ! k) { for (int j 0; j n; j) { double tmp A[k * n j]; A[k * n j] A[max_row * n j]; A[max_row * n j] tmp; } double tmp b[k]; b[k] b[max_row]; b[max_row] tmp; } // 消元 for (int i k 1; i n; i) { double factor A[i * n k] / A[k * n k]; for (int j k; j n; j) { A[i * n j] - factor * A[k * n j]; } b[i] - factor * b[k]; } } // 回代 for (int i n - 1; i 0; i--) { double sum b[i]; for (int j i 1; j n; j) { sum - A[i * n j] * x[j]; } if (fabs(A[i * n i]) 1e-12) return -1; x[i] sum / A[i * n i]; } return 0; }列主元选取的额外开销是每步消元多了一次O(n)的扫描相对于消元本身的O(n²)来说可以忽略。但它带来的数值稳定性提升是巨大的。我实测过一个条件数在10⁸左右的矩阵朴素消元的结果完全不可用列主元消元还能保持3到4位有效数字。4.3 消元过程中的精度损失与应对即使有列主元高斯消元仍然可能损失精度。主要来源是消元过程中的减法操作当两个相近的大数相减时有效数字会大量丢失。这在正规方程方法中尤其严重因为XᵀX的条件数是X的条件数的平方。假设X的条件数是10⁴XᵀX的条件数就是10⁸double类型大约有15到16位有效数字消元后可能只剩7到8位。如果X的条件数达到10⁶XᵀX就是10¹²结果基本不可信了。应对策略有几个层次。最基础的是标准化特征把条件数降下来。进阶一点可以用迭代精化先用高斯消元得到一个近似解然后计算残差r b - Aβ再解Aδ r得到修正量δ把β δ作为更精确的解。迭代精化能显著改善病态方程组的解质量代价是每次迭代多一次消元。void iterative_refinement(double *A, double *b, double *x, int n, int max_iter) { double *r (double*)malloc(n * sizeof(double)); double *delta (double*)malloc(n * sizeof(double)); double *A_copy (double*)malloc(n * n * sizeof(double)); for (int iter 0; iter max_iter; iter) { // 计算残差 r b - A*x for (int i 0; i n; i) { double sum 0.0; for (int j 0; j n; j) { sum A[i * n j] * x[j]; } r[i] b[i] - sum; } // 解 A*delta r memcpy(A_copy, A, n * n * sizeof(double)); gaussian_elimination_pivot(A_copy, r, delta, n); // 更新解 double max_delta 0.0; for (int i 0; i n; i) { x[i] delta[i]; if (fabs(delta[i]) max_delta) max_delta fabs(delta[i]); } if (max_delta 1e-10) break; // 收敛 } free(r); free(delta); free(A_copy); }注意迭代精化需要保留原始的A矩阵因为每次迭代都要用原始A计算残差。上面的代码里我用了A_copy来保存消元前的副本实际工程中应该把原始A和用于消元的工作副本分开管理。4.4 求解结果的验证残差检查解出β之后怎么知道结果靠不靠谱最直接的检查是计算残差向量r Y - Xβ看它的范数是否足够小。如果残差范数相对于Y的范数很大说明拟合效果差可能是模型设定有问题也可能是数值求解出了问题。double compute_residual_norm(Matrix *X, double *y, double *beta, int n) { double sum_sq 0.0; for (int i 0; i n; i) { double pred 0.0; for (int j 0; j X-cols; j) { pred MAT(X, i, j) * beta[j]; } double resid y[i] - pred; sum_sq resid * resid; } return sqrt(sum_sq); }如果残差范数异常大可以进一步检查XᵀX的条件数。条件数的粗略估计可以用最大对角元素和最小对角元素的比值精确计算需要做奇异值分解那超出了这篇的范围。但粗略估计已经能发现大部分问题。我个人的经验是标准化之后如果XᵀX的对角元素最大值和最小值之比超过10⁶就要警惕了。这时候要么检查特征之间是否存在严重的多重共线性要么考虑用正则化方法岭回归来稳定求解。正则化的内容会在后续篇章展开。5. 完整流程串联与实测中的意外情况5.1 从原始数据到系数输出的完整调用链把前面各个模块串起来一个完整的多元线性回归求解流程是这样的int multiple_linear_regression( double *features, // n*p行优先 double *target, // n int n, int p, double *beta_out // p1输出系数 ) { // 1. 标准化特征原地修改 double *mean (double*)malloc(p * sizeof(double)); double *std (double*)malloc(p * sizeof(double)); standardize(features, n, p, mean, std); // 2. 构建设计矩阵 Matrix *X build_design_matrix(features, target, n, p); if (!X) { free(mean); free(std); return -1; } // 3. 组装正规方程 Matrix *XtX compute_xtx(X); double *Xty compute_xty(X, target); if (!XtX || !Xty) { /* 清理并返回 */ } // 4. 高斯消元求解 int ret gaussian_elimination_pivot(XtX-data, Xty, beta_out, p 1); // 5. 残差检查 if (ret 0) { double resid compute_residual_norm(X, target, beta_out, n); printf(Residual norm: %.6f\n, resid); } // 6. 清理 matrix_free(X); matrix_free(XtX); free(Xty); free(mean); free(std); return ret; }这个流程里标准化是在构建设计矩阵之前做的所以设计矩阵里的特征列已经是标准化后的值。截距列仍然是1对应标准化空间中的截距。5.2 一个容易忽略的坑标准化与截距的交互标准化之后截距β₀的含义变了。在标准化空间中当所有特征取值为0时也就是原始空间中所有特征取均值时预测值等于β₀。所以β₀实际上等于目标变量的均值如果目标没有中心化的话。如果你需要原始空间中的截距可以用下面的公式转换β₀_original β₀_standardized - Σⱼ (βⱼ_standardized * meanⱼ / stdⱼ)而原始空间中的斜率就是βⱼ_standardized / stdⱼ。这个转换在需要解释模型系数的时候很重要因为标准化空间的系数不能直接解读为特征每增加一个单位目标变化多少。我踩过这个坑有一次做完标准化拟合直接把系数拿去做业务解释结果发现截距完全对不上。后来才意识到标准化改变了截距的基准点。从那以后我养成了一个习惯拟合完成后先把系数转换回原始空间再做任何解释。5.3 多重共线性的识别与处理当两个或多个特征高度相关时XᵀX接近奇异系数估计会变得极不稳定。表现是系数的绝对值很大符号可能和业务直觉相反而且换一组样本系数就剧烈变化。识别多重共线性的简单方法是计算方差膨胀因子VIF。对每个特征j用其他所有特征对它做回归得到R²ⱼ则VIFⱼ 1 / (1 - R²ⱼ)。VIF大于10通常认为存在严重共线性。在C语言里计算VIF需要多次调用回归代码量不小。一个更简单的替代方案是检查XᵀX的对角元素。标准化之后XᵀX的对角元素应该都接近n因为每个特征的方差为1平方和为n。如果某个对角元素远大于n说明该特征和其他特征存在线性关联。处理多重共线性的方法有删除相关性高的特征之一、用主成分分析降维、或者用岭回归。这些方法各有适用场景没有万能方案。我的经验是如果业务上能解释清楚哪些特征冗余优先删除如果所有特征都重要考虑岭回归。5.4 数据量n和特征数p的关系正规方程方法要求n ≥ p1否则XᵀX是奇异的方程组有无穷多解。即使n略大于p1XᵀX也可能接近奇异系数估计不稳定。经验法则是n至少是p的5到10倍模型才比较可靠。如果n不够要么收集更多数据要么减少特征数量要么用正则化方法。正则化通过给XᵀX的对角线加上一个正数λ保证矩阵可逆即使n p也能得到唯一解。这是岭回归的核心思想后续篇章会详细讲。5.5 实测中的性能观察我用一个模拟数据集做了简单测试n10000p20在普通笔记本上编译运行gcc -O2从数据标准化到输出系数总耗时大约0.15秒。其中XᵀX的计算占了大约60%的时间高斯消元占30%其余是标准化和残差检查。如果p增大到100XᵀX的计算时间按O(n·p²)增长高斯消元按O(p³)增长。p100时总耗时大约2秒其中高斯消元占比上升到50%以上。如果p继续增大高斯消元会成为瓶颈这时候可以考虑用Cholesky分解因为XᵀX是对称正定的运算量能减半。不过对于大多数实际应用场景p很少超过几百这个性能已经足够了。真正需要担心的是数值稳定性而不是速度。5.6 编译与调试的实用建议编译的时候一定要开优化和警告gcc -O2 -Wall -Wextra -o regression regression.c -lm-Wall -Wextra能帮你发现很多潜在问题比如未使用的变量、隐式类型转换、格式化字符串不匹配等。-lm链接数学库因为用到了fabs和sqrt。调试数值代码的时候printf是最直接的工具。在关键步骤打印中间结果XᵀX的对角元素、消元后的上三角矩阵、残差范数。如果发现某个值异常比如NaN或Inf往回追溯是哪一步产生的。还有一个技巧用一个小规模的手算例子验证代码。比如n5p2自己用计算器算一遍XᵀX和XᵀY和程序输出对比。这能帮你快速定位是数学推导错了还是代码实现错了。我每次写新的数值算法都会先做这个验证比直接上大规模数据调试效率高得多。提示如果程序输出NaN最常见的原因是除零或者对负数开方。检查标准化时std是否为零检查高斯消元的主元是否为零。如果输出Inf通常是数值溢出检查是否有极大的数参与运算。5.7 后续可以扩展的方向这篇实现的是最基础的最小二乘直接求解。在实际项目中还有几个方向值得深入一是梯度下降法适合n和p都很大的场景不需要显式构造XᵀX二是岭回归和Lasso处理多重共线性和特征选择三是模型评估包括R²、调整R²、交叉验证等四是预测新样本需要保存标准化参数并做一致的变换。这些内容我会在后续篇章里逐一展开。先把这篇的直接求解方法吃透理解每一步的数学含义和数值特性后面的扩展就是水到渠成的事。