简介面向希望用 C# 理解并复现线性拟合的开发者与数据分析初学者这一压缩包围绕最小二乘法求解直线 ymxb 展开重点解决斜率 m、截距 b 的计算及预测函数封装问题并可与 Excel 线性拟合结果对照验证。包内共 45 个文件以 14 个 cs 源码为主配合 csproj、sln 工程文件以及 dll、exe、config、resx、txt 等运行与配置资料另有少量 pdb、cache、suo 等构建缓存整体约 87KB目录保留了较完整项目结构。已有 413 人学习适合据此梳理均值—偏差—斜率—截距的计算流程阅读拟合函数实现理解数据读取、结果输出与误差处理的组织方式并在此基础上扩展到多项式等曲线拟合场景提升 C# 数据处理能力。1. 线性拟合公式求解一行代码背后还有三道判断把一组散点丢进np.polyfit斜率和截距立刻就出来了这大概是很多人对线性拟合公式求解的全部印象。可真到用的时候麻烦往往不在公式本身明明 R² 有 0.98换个样本系数就翻了符号明明只是加了一列特征程序直接抛Singular matrix明明两组数据长得差不多一组求出来的截距是 3.2另一组是 -47。这些问题都不是“公式记错了”而是求解路径、矩阵形态和数据量纲在背后主导。这篇内容面向的是需要把线性拟合真正用起来的场景传感器标定、实验数据建模、报表里的趋势线、特征工程里的基线模型。读完你该能做到三件事从误差平方和独立推出正规方程不再死记知道正规方程、QR、SVD 三条路各自什么时候会出错能对着一份跑歪的结果按现象倒推是哪一环出了问题。2. 最小二乘闭式解是怎么推出来的从误差平方和到正规方程2.1 把“拟合得好”翻译成一个能求导的目标函数“拟合得好”这句话本身没法计算得先翻译成一个标量。最常用的翻译方式是残差平方和对每个样本预测值与真实值之差叫残差把所有残差平方后加起来得到一个关于参数 w 的函数。平方的作用有两个一是让正负残差不会互相抵消二是让这个函数处处可导后面能直接求极值。写成矩阵形式设设计矩阵 X 形状为 m×nm 个样本、n 个待求参数第一列通常是全 1 的常数列对应截距响应向量 y 形状为 m×1参数向量 w 形状为 n×1那么残差向量是 y - Xw平方和可以展开成S(w) (y - Xw)ᵀ(y - Xw) yᵀy - 2wᵀXᵀy wᵀXᵀXw这一步展开没有做任何近似纯粹是转置和乘法的分配律。展开之后能看清三件事yᵀy 与 w 无关是个常数中间项是 w 的一次项最后一项是 w 的二次型系数矩阵 XᵀX 一定是对称的。也就是说不管数据长什么样S(w) 都是一个开口向上的二次曲面半正定最小值要么在唯一一点取到要么在一整片区域上取到后者就对应 X 列不满秩的情况。这里有个容易被含糊过去的概念“线性”指的是对参数 w 线性不是对输入 x 线性。所以 x 取平方、取对数、取两个特征的乘积放进 X模型依然叫线性回归闭式解照样成立。很多人以为多项式拟合是另一套东西其实只是把设计矩阵换成了范德蒙矩阵求解过程一模一样。2.2 求导置零正规方程怎么来的以及它默认了什么对 S(w) 关于 w 求梯度用到两条矩阵求导规则一次项 wᵀa 的梯度是 a二次型 wᵀAwA 对称的梯度是 2Aw。于是∂S/∂w -2Xᵀy 2XᵀXw令梯度为零两边同除 2得到 XᵀXw Xᵀy这就是正规方程。因为 XᵀX 半正定这个驻点必然是全局最小值点不需要再验证二阶条件。当 X 列满秩时XᵀX 正定可逆解唯一写出来是 w (XᵀX)⁻¹Xᵀy。正规方程还有一个几何解释比代数推导更直观把等式写成 Xᵀ(y - Xw) 0意思是残差向量与 X 的每一列都正交。而 X 的列张成的空间就是所有可能的预测值构成的空间所以最小二乘解等价于“把 y 正交投影到 X 的列空间上”。这也解释了为什么残差之和恒为零常数列在其中以及为什么加了常数列之后残差均值必然为 0。推导干净但它默认了三件事这三件事在真实数据里经常不成立。第一误差零均值、同方差、互不相关违反这一条就该上加权最小二乘第二X 列满秩违反这一条正规方程直接无解或有无穷多解第三数值上 XᵀX 可逆且良态这一条最隐蔽因为数学上可逆不等于浮点运算里算得准。前两条是建模问题第三条是数值问题后面两章分别处理。2.3 20 行 numpy 复现闭式解从设计矩阵到系数把上面的推导落地最短的实现只需要一个矩阵乘法加一次线性方程求解。关键点不在代码量而在两个容易写错的地方常数列有没有加、解方程用的是solve还是inv。import numpy as np def design_matrix(x, add_interceptTrue): 把一维特征扩成设计矩阵多特征时直接传二维数组 x np.asarray(x, dtypenp.float64) if x.ndim 1: x x.reshape(-1, 1) if add_intercept: # 常数列放第一列对应参数 w[0] 即截距 x np.column_stack([np.ones(x.shape[0]), x]) return x def fit_normal_equation(X, y): 正规方程闭式解 w (X^T X)^-1 X^T y X np.asarray(X, dtypenp.float64) y np.asarray(y, dtypenp.float64).reshape(-1) if X.shape[0] ! y.shape[0]: raise ValueError(f样本数不一致: X {X.shape}, y {y.shape}) if not (np.isfinite(X).all() and np.isfinite(y).all()): raise ValueError(输入含 NaN/Inf先清洗再拟合) A X.T X # n x n 对称矩阵 b X.T y # n 维右端向量 w np.linalg.solve(A, b) # 用 solve不要用 inv(A) b return w def residuals(X, y, w): return np.asarray(y, dtypenp.float64) - X w if __name__ __main__: rng np.random.default_rng(7) x np.linspace(0, 10, 200) y 2.5 * x - 1.8 rng.normal(0, 1.2, sizex.size) # 真值 k2.5, b-1.8 X design_matrix(x) w fit_normal_equation(X, y) print(截距, 斜率 , w) # 应接近 [-1.8, 2.5] r residuals(X, y, w) print(残差均值 , r.mean()) # 含常数列时应接近 0design_matrix里np.column_stack的顺序决定了系数的顺序截距在前、特征在后后面所有解释都按这个约定来别中途换顺序否则系数对不上号。fit_normal_equation里用np.linalg.solve(A, b)而不是np.linalg.inv(A) b原因是前者走 LU 分解只解一次方程后者显式求逆浮点误差更大复杂度也更高而且求出来的逆矩阵大部分元素根本用不上。这是我在代码评审里见得最多的一处“看起来对但没必要”的写法。residuals单独拎出来是个习惯问题但它值得单独成函数后面做残差图、算残差标准误、检查自相关都要反复用它写成一个函数比到处复制y - X w靠谱。最后那段自测用固定随机种子生成带噪声的直线真值 k2.5、b-1.8跑出来的系数应该贴着真值残差均值应该在 1e-15 量级——如果残差均值明显不为 0八成是常数列丢了。3. 正规方程、QR 与 SVD 怎么选同一份数据三种解法的实测差异3.1 正规方程的代价条件数被平方正规方程的复杂度很讨喜构造 XᵀX 需要约 mn² 次乘加解 n×n 方程组约 n³/3整体 O(mn² n³)n 小的时候非常快。问题出在数值上。矩阵的条件数 κ 衡量的是“输入微小扰动被放大多少倍”而 κ(XᵀX) κ(X)²也就是说构造 XᵀX 这一步把病态程度平方了。float64 有大约 16 位十进制有效数字。如果 κ(X) 1e8那么 κ(XᵀX) 1e16正规方程求解过程中 16 位精度基本被吃干净算出来的系数可能只剩两三位可信符号翻转都不奇怪。经验阈值我一般这么定κ(X) 超过 1e6 就别信正规方程的系数了超过 1e8 基本只能看预测值不能看系数。这不是玄学是浮点数的硬约束。病态在真实数据里怎么来的三种最常见量纲差异巨大一个特征单位是毫米、另一个是千牛数值差 1e9、特征之间高度相关两个传感器测同一个物理量、多项式阶数偏高x 从 0 到 100x⁵ 就是 1e10 量级各列差着好几个数量级。前两种可以靠中心化、标准化缓解第三种靠换基或换解法。3.2 QR 与 SVD多花的浮点运算买到什么QR 分解把 X 拆成 X QRQ 是 m×n 列正交矩阵R 是 n×n 上三角矩阵。代进正规方程RᵀQᵀQRw RᵀQᵀy因为 QᵀQ I化简成 Rw Qᵀy而 R 是上三角回代一下就解完了。整个过程不构造 XᵀX所以条件数没有被平方稳定性比正规方程高一个量级。代价是运算量Householder QR 约 2mn² - 2n³/3大概是正规方程的两倍但换来的是可信的系数。SVD 把 X 拆成 USVᵀ解写成 w VS⁺Uᵀy其中 S⁺ 是把奇异值取倒数再转置。它最稳而且额外给了两样东西奇异值谱能直接看出矩阵的秩和病态程度截断小奇异值把小于阈值的 1/s 直接置零就能得到最小范数解这对列不满秩的情况是唯一的正确处理方式。代价是运算量最大约 4mn² - 4n³/3是 QR 的两倍左右。实际选型的经验是n 小于 30、数据量纲统一、κ(X) 在 1e4 以内正规方程够用且最快数据来源杂、量纲不统一或者需要看系数是否可信直接上 QR列可能不满秩比如加了独热编码但某类只有一个样本、需要明确处理共线性用 SVD 或直接调np.linalg.lstsq底层走 SVD 的 gelsd。lstsq是我在业务代码里的默认选择除非性能真的卡住否则没必要为了省那点浮点运算去手写正规方程。方法典型复杂度数值稳定性适用场景正规方程mn² n³/3低条件数平方n 小、数据已标准化、快速试算QR 分解2mn² - 2n³/3中高一般业务回归、需可信系数SVD4mn² - 4n³/3最高可截断共线严重、可能不满秩、需诊断3.3 同一份数据跑三种解法把差异打出来光看表格没感觉构造一个刻意病态的例子跑一遍更直观两列特征几乎共线再加一列量纲巨大的特征让 κ(X) 上到 1e7 左右然后看三种解法给出的系数差多少。import numpy as np def solve_three_ways(X, y): X np.asarray(X, dtypenp.float64) y np.asarray(y, dtypenp.float64).reshape(-1) # 1) 正规方程 w_ne np.linalg.solve(X.T X, X.T y) # 2) QRX QR, R w Q^T y Q, R np.linalg.qr(X, modereduced) w_qr np.linalg.solve(R, Q.T y) # 3) SVD 伪逆按机器精度截断小奇异值 U, s, Vt np.linalg.svd(X, full_matricesFalse) tol s.max() * max(X.shape) * np.finfo(np.float64).eps s_inv np.where(s tol, 1.0 / s, 0.0) w_svd Vt.T (s_inv * (U.T y)) return w_ne, w_qr, w_svd rng np.random.default_rng(42) m 300 x1 rng.normal(0, 1, m) x2 x1 rng.normal(0, 1e-4, m) # 与 x1 几乎共线 x3 rng.normal(0, 1, m) * 1e6 # 量纲差 6 个数量级 X np.column_stack([np.ones(m), x1, x2, x3]) y 3.0 * x1 5.0 * x3 rng.normal(0, 1, m) print(cond(X) , np.linalg.cond(X)) for name, w in zip([normal, qr, svd], solve_three_ways(X, y)): print(f{name:7}: {np.array2string(w, precision4)})跑出来会看到cond(X)在 1e7 上下正规方程给出的 x1、x2 系数可能一个是 800、一个是 -800量级离谱且符号相对QR 和 SVD 的结果彼此接近但 x1、x2 的系数依然不稳定——因为数据本身就分不清这两列各自的贡献任何解法都救不了只能说明共线性是数据问题不是算法问题。这一个实验同时说明了两件事换解法能解决“算不准”解决不了“分不清”。参数上有几个点值得说。np.linalg.qr的modereduced返回 m×n 的 Q 和 n×n 的 R在 m 远大于 n 时内存友好如果只需要解方程scipy.linalg.lstsq指定lapack_drivergelsy走 QR 路线比 numpy 的 SVD 版本快。SVD 里的截断阈值tol用了 numpy 内部同款公式最大奇异值 × 矩阵最大维度 × 机器精度自己写伪逆时直接把这段抄过去别硬编码 1e-10 这种魔数——数据量级一变硬编码的阈值立刻失效。4. 把公式改成能用的形态加权、岭回归与多项式拟合的三种变体4.1 加权最小二乘异方差数据先开根号再回归普通最小二乘假设每个样本的误差方差相同。真实数据经常不满足高浓度区间的测量误差大、后期的数据点更可信、聚合数据里样本量大的点更稳。这时候给每个样本一个权重 wᵢ目标函数改成 Σwᵢ(yᵢ - xᵢᵀw)²求导置零得到 XᵀWXw XᵀWy解是 w (XᵀWX)⁻¹XᵀWy。直接按这个式子写也能跑但更省事、数值上也更稳的做法是做变量替换令 √W 为权重的平方根构成的对角矩阵定义 X √W X、y √W y那么加权问题就变成了普通最小二乘问题直接调lstsq就行。这个变换的价值在于复用了已经验证过的稳定求解器不用为加权单独写一套。权重的取法有讲究。如果已知每个观测的测量标准差 σᵢ权重取 1/σᵢ²这是统计上最优的如果不知道但知道每个点是 k 次重复测量的均值权重取 k如果是凭经验给可信度那就只能是工程近似这时候要接受结果不如前两种严谨。我一般会先画一张“残差绝对值 vs 拟合值”的散点图如果呈喇叭形拟合值越大残差越散说明异方差确实存在加权是必要的不是可选项。import numpy as np def fit_wls(X, y, w): 加权最小二乘先开根号做变量替换再走普通最小二乘 X np.asarray(X, dtypenp.float64) y np.asarray(y, dtypenp.float64).reshape(-1) w np.asarray(w, dtypenp.float64).reshape(-1) if (w 0).any(): raise ValueError(权重必须为正0 权重请直接删样本) sw np.sqrt(w).reshape(-1, 1) Xw X * sw # 逐行缩放 yw y * np.sqrt(w) return np.linalg.lstsq(Xw, yw, rcondNone)[0]w 0的检查不能省。权重为 0 意味着这个样本不参与拟合与其留着不如直接从 X 和 y 里删掉否则 √0 0 会让该行全零虽然数学上不报错但会悄悄改变矩阵的秩后面算条件数时会看到一个莫名其妙的零奇异值。负权重则直接是数据录入错误。4.2 岭回归给 XᵀX 加一条对角线到底在压什么当特征共线严重时XᵀX 的最小特征值接近 0它的倒数巨大导致系数的方差爆炸——换个样本系数可能从 500 跳到 -480。岭回归在目标函数里加一个 L2 惩罚项 λ‖w‖²正规方程变成 (XᵀX λI)w Xᵀy最小特征值被抬到至少 λ条件数从 σmax²/σmin² 降到 (σmax² λ)/(σmin² λ)病态被强行压住。代价是解变成有偏的λ 越大系数被压缩得越厉害偏差越大但方差越小。λ 的选择就是在这对矛盾里找平衡点常用留一交叉验证岭回归有个闭式解可以免去重训 n 次的开销也可以简单网格搜索。λ 0 退化回普通最小二乘λ 趋于无穷时所有系数趋于 0。一个必须注意的细节截距不参与惩罚。因为惩罚项是给“特征的贡献”加的截距反映的是数据的整体水平把它一起压缩会让预测整体偏移。实现上的标准做法是先把 X 和 y 都中心化减去各自的均值中心化之后截距项自动消失惩罚只作用在特征系数上最后再用 b ȳ - x̄ᵀw 把截距还原回来。import numpy as np def fit_ridge(X, y, lam): 岭回归中心化后惩罚特征系数截距单独还原 X np.asarray(X, dtypenp.float64) y np.asarray(y, dtypenp.float64).reshape(-1) x_mean X.mean(axis0) y_mean y.mean() Xc X - x_mean yc y - y_mean n Xc.shape[1] A Xc.T Xc lam * np.eye(n) w np.linalg.solve(A, Xc.T yc) b y_mean - x_mean w return w, b def ridge_path(X, y, lams): return [(lam, *fit_ridge(X, y, lam)) for lam in lams]lam的量级取决于 XᵀX 的尺度如果特征没标准化XᵀX 的对角元素可能是 1e8 也可能是 1e-2同一个 λ 在两个数据集上完全不可比。所以岭回归之前必须先标准化特征减均值除标准差让 XᵀX 的对角线接近样本量 mλ 才有可迁移的意义。这是我在复现别人的岭回归结果时最常发现的差异来源。4.3 多项式拟合范德蒙矩阵必须先中心化缩放多项式拟合的设计矩阵是范德蒙矩阵第一列全 1第二列 x第三列 x²一直到 x^d。它的条件数随阶数指数增长x 范围在 [0, 10] 时5 阶的 κ 就能到 1e8 以上10 阶直接爆掉。所以多项式拟合的第一原则不是选阶数而是先把 x 变换到 [-1, 1] 或者做 z-score 标准化。变换之后同样阶数的范德蒙矩阵条件数能降好几个数量级正规方程都能算得动。但这里有个必须守住的一致性变换参数均值、标准差或者区间的中点、半宽是从训练集算出来的预测时必须用同一组参数不能用新数据的均值和标准差重新算一遍否则模型相当于被偷偷改掉了。import numpy as np def poly_design(x, deg, statsNone): 范德蒙设计矩阵stats(center, scale) 用于复用训练集的变换参数 x np.asarray(x, dtypenp.float64).reshape(-1) if stats is None: center x.mean() scale x.std(ddof0) scale scale if scale 0 else 1.0 # 常数输入的保护 stats (center, scale) z (x - stats[0]) / stats[1] # increasingTrue 让第一列是 z^0常数列系数顺序从低阶到高阶 return np.vander(z, deg 1, increasingTrue), stats def fit_poly(x, y, deg): X, stats poly_design(x, deg) w np.linalg.lstsq(X, np.asarray(y, dtypenp.float64), rcondNone)[0] return w, stats def predict_poly(x, w, stats, deg): X, _ poly_design(x, deg, stats) # 必须传 stats不能重新计算 return X wnp.vander默认是降幂排列第一列是 x^d加increasingTrue改成升幂系数顺序才和直觉一致w[0] 是常数项。这个参数不写的话很容易出现“拟合看起来没问题但常数项解读错了”的情况尤其是把系数拿去和业务含义对照的时候。还有个绕不开的现象阶数高了之后拟合区间内很贴合稍微外推一点点就飞出去残差图两端翘得像张弓。这是高次多项式的固有性质不是求解器的问题换 SVD 也救不了。可行的做法是控制阶数一般不超过 35 阶、在数据边界处不做外推、或者改用分段拟合、样条这类局部方法。我自己在工艺曲线拟合上踩过这个坑为了把 R² 从 0.97 提到 0.995 硬上到 7 阶结果超出标定范围 5% 的预测偏差到了 30%。5. 线性拟合求解的五个踩坑现场从现象倒推到根因5.1 数值层的翻车NaN、奇异矩阵和精度不够现象一程序抛LinAlgError: Singular matrix或者结果里全是 NaN。原因通常有两类X 里存在完全重复的列比如两次导出了同一个特征、独热编码没有去掉冗余类别或者数据里有 NaN/Inf 悄悄参与了运算。重复列让 XᵀX 秩亏solve直接报错而lstsq遇到 NaN 不会报错它会一路算下去返回 NaN这才是最危险的。解决办法是先做一次体检np.isfinite(X).all()查 NaN/Infnp.linalg.matrix_rank(X)和X.shape[1]对比查秩亏独热编码用drop_firstTrue去掉一列。现象二同一份数据numpy 算出来的系数和别人用表格软件算的对不上差在小数点后好几位。原因大概率是精度如果 X 或 y 是 float32从某些传感器或图像流程里带出来的有效数字只有 7 位左右条件数一上千就明显偏差。另一个常见原因是某个环节用了np.linalg.inv显式求逆比solve多一轮误差累积。解决办法很简单入口处强制dtypenp.float64全链路不用inv。血泪经验是.astype(np.float32)在某些流程里是为了省内存但绝不能出现在回归求解这一步。现象三系数偶尔正常、偶尔正常得离谱换个随机种子就变天。原因是条件数处在临界区1e51e7 之间正规方程的误差是不稳定的有时候舍入误差恰好抵消有时候恰好放大。这种“看起来能用”的状态最坑。解决办法是改用 QR 或 SVD同时在日志里打一行np.linalg.cond(X)一旦超过 1e6 就告警。别指望靠调参蒙对这是个确定性的数值问题。5.2 建模层的翻车截距、共线性和高阶外推现象四拟合线被强行拉过原点残差均值明显不为零。原因是设计矩阵没加常数列模型退化成 y kx。这时候残差和不为零正规方程的解等价于“在过原点的直线里找最贴近的”如果数据整体有个明显截距整条线会被拽歪。判断方法很直接算一下残差均值含截距的模型残差均值应该在 1e-12 以下相对量级明显不为 0 就是截距丢了。解决办法是构造设计矩阵时统一加常数列或者直接调带fit_interceptTrue的封装别手写裸矩阵。现象五R² 高达 0.99但换一批新数据预测就明显偏或者两列特征系数一个大正一个大负且量级夸张。前者的原因是过拟合典型触发场景就是高阶多项式和特征数量接近样本量后者的原因是共线性模型分不清两列各自的贡献只能靠一正一负的夸张组合去凑。这两种情况的共同点是训练集指标好看、泛化差所以不能用训练集 R² 判断模型好坏。解决办法多项式阶数控制在 35 阶以内共线时先看 VIF 或相关矩阵相关系数超过 0.95 就删一列或合并实在要保留就上岭回归验证一律留出独立测试集看测试集 RMSE 而不是训练集 R²。这五个坑里前三个属于“算得对不对”后两个属于“模型对不对”。排查顺序我习惯先从数值层查起因为数值问题会伪装成模型问题——系数符号翻转看起来像业务结论其实只是精度丢了。6. 拟合结果的验收与进阶把系数的不确定性量化出来系数算出来不等于结论成立还得知道它有多准。最小二乘框架下有个现成的答案残差方差 σ² RSS/(m - n)其中 RSS 是残差平方和m - n 是自由度减去参数个数。系数的协方差矩阵是 σ²(XᵀX)⁻¹对角线开根号就是每个系数的标准误系数除标准误得到 t 值t 的绝对值小于 2 左右的时候就别急着宣称这个特征有影响了。这个式子有个前提XᵀX 得算得准。前面说过条件数大的时候 (XᵀX)⁻¹ 本身就不可信标准误跟着失真。所以我的习惯是先看条件数条件数大的时候改用 bootstrap对样本做有放回重采样重跑几百次求解看系数的分布。bootstrap 不在乎矩阵病不病态它只关心“换一批样本系数会怎么变”这个视角反而更贴近真实的不确定性。诊断项计算方式判断标准残差均值r.mean()含截距时应接近 0残差标准差r.std(ddofn_params)与业务可接受误差对比条件数np.linalg.cond(X)超过 1e6 系数不可信系数 t 值w / se绝对值小于 2 视为不显著留出集 RMSE在独立测试集上算与训练集 RMSE 差距不应过大还有一个我常用来验证实现是否正确的技巧造数据反解。先设一组已知的真值系数用 X w_true 生成无噪声的 y跑一遍求解看能不能精确还原到 1e-10 以内再加微小噪声看系数是否在真值附近小幅波动。这个自检能揪出绝大部分实现错误——列顺序搞反、常数列漏加、幂次排列颠倒全都会在这一步暴露。比拿真实数据反复调参快得多也是我每次重写拟合模块时的第一件事。最后留一个习惯任何一次线性拟合求解日志里至少记下样本数、特征数、条件数、残差标准差四个数字。出事的时候这四个数字就是后悔药——条件数大说明是数值问题残差标准差大说明模型漏了变量样本数和特征数接近说明过拟合风险。我以前图省事只记 R²结果线上预测偏了一个量级回头查连当时用了哪些特征都对不上只能整段重跑。希望帮到你。本文还有配套的精品资源点击获取