做了十来年数据分析回头一看真正让我在工作里反复受益的并不是那些花哨的新模型反而是最基础的最小二乘法。帮传感器做标定曲线、给实验数据找趋势线、评估投放效果、做量化因子分析底层都是这一套。很多人一听到“最小二乘法公式”就觉得高深其实说白了它就是在回答一个非常实际的问题面对一堆散点怎么找到一条最合适的曲线让每个点到曲线的误差平方和最小。这篇内容适合刚接触回归的初学者也适合那些公式背过、但还没真正理解它为什么能行的朋友。我见过不少人都把最小二乘当成“会用np.polyfit就够了”的黑盒结果遇到数据稍微脏一点就翻车。这篇文章我会从思想、公式推导、代码实现、再到各种坑完整讲一遍我自己的理解。看完你不仅能手推最小二乘法公式还能在实际项目里用得更稳。1. 最小二乘法的核心思想找一条让误差平方和最小的线1.1 从散点拟合说起它到底在解决什么问题先假设一个常见场景。我手里有一组实验数据横轴是温度纵轴是传感器的输出电压。理论上它们应该是线性的但实际测量总会有噪声画出来就是一堆大致呈直线趋势的散点。这时候我想找一条直线用它来描述“温度和电压的关系”这样以后只要测到电压我就能反推温度。问题来了能用这一堆点确定无数条直线哪一条才是“最合适”的最小二乘法给出的答案是让所有数据点到这条直线的垂直距离的平方和最小。垂直距离就是每个点的真实值减去直线预测值也就是残差。把所有残差平方一下再加起来得到一个数字这个数字越小说明直线整体上离所有点越近。为什么要“平方”而不是直接加绝对值这个问题特别关键。如果只取误差的绝对值再求和确实也能衡量“整体偏差”但绝对值函数在零点不可导后面做最优化很麻烦。平方之后的函数处处光滑可导能用求导解方程的方式直接算出最优参数这条路太顺了。更深一层说“误差平方和最小”背后有一个统计假设。如果数据的噪声服从独立同分布的高斯分布那么最小化误差平方和等价于最大化高斯分布下的似然函数。换句话说当你觉得测量误差是随机扰动、没有系统性偏差时最小二乘法给出的参数就是最可能产生这组观测值的参数。这也是它名字里“法”字的底气所在。1.2 为什么偏偏是“平方”换绝对值行不行这个问题我当年也纠结过很久。后来在项目里处理过一批带尖峰噪声的数据才真正体会到平方损失的脾气它对大误差的惩罚是二次的一个偏离很远的点它的影响会随着距离变大迅速增加。这既是优点也是缺点。优点是在误差近似服从正态分布的时候平方损失能充分利用所有数据的信息估计出来的参数方差小、稳定性好。统计里有个高斯-马尔可夫定理说的就是在线性回归的框架下如果误差满足零均值、同方差、互不相关这三个条件最小二乘估计在所有线性无偏估计里方差最小。简单说就是它在这个前提下是最优的。缺点是它对异常值极度敏感。一个离群点如果离主趋势特别远它的残差平方可能比其余几十个点的残差总和还大拟合出来的直线会被这个点强行拉过去。实际工程数据里传感器偶尔跳变、人工录入错一位数都可能导致这种情况。这时候绝对值损失反而更稳所以后来又发展出Huber损失这种折中方案——小误差时用平方大误差时用绝对值。这个我们后面还会提。所以“平方”不是随手选的它是在正态噪声假设下从最大似然原理自然导出的结果。理解这一点你才能真正明白为什么有时候最小二乘不好用——不是公式错了而是你的噪声分布不符合它的假设。1.3 适用场景与前置假设根据我自己的经验最小二乘法在下面这些场景里最可靠数据特点典型场景最小二乘是否合适近似线性关系噪声均匀传感器标定、电压-温度曲线非常适合关系明显非线性但可变量变换指数衰减、对数增长变换后适合特征高度相关房价回归、经济指标分析不直接适合需正则化离群点较多实验数据含粗大误差不直接适合需鲁棒方法样本量小于特征数基因表达、文本分类严重不适合欠定问题在正式上手之前我会习惯性地做两个检查第一画散点图看关系是不是大致线性。第二算一下残差的标准差判断噪声是不是近似恒定。如果噪声随着x增大而明显变大说明方差不齐直接做普通最小二乘就会让大噪声区域的点过度影响结果。这时候要么对y做加权要么先做变换。记住一个总原则最小二乘法不是“万能拟合器”它是一个带前置假设的参数估计方法。前置假设满足得越好结果越接近真实关系假设被破坏得越厉害结果越离谱。2. 从公式到代码一元线性回归的完整推导2.1 建立目标函数现在正式来推一次最小二乘法公式。以最简单的一元线性回归为例假设拟合直线是y wx b其中w是斜率b是截距。对于第i个样本点(x_i, y_i)直线预测值是 wx_i b真实值和预测值的差就是e_i y_i - (wx_i b)把所有点的误差平方加起来得到了目标函数L(w, b) Σ(y_i - wx_i - b)²我们要求的就是让L取最小值时的w和b。只要能把目标函数写出来剩下的全是数学操作。这一步特别重要因为无论后面变成多元回归、多项式回归还是岭回归都是在这个目标函数上加加减减。2.2 求偏导、解方程L是关于w和b的二次函数这种函数要么有唯一最小值要么没有最大值形状像一个碗。要求碗底的位置就对w和b分别求偏导并令偏导等于0。对b求偏导∂L/∂b -2Σ(y_i - wx_i - b) 0化简得Σy_i - wΣx_i - nb 0也就是 b ȳ - wx̄这里 x̄ 是x的均值ȳ 是y的均值。这个式子的含义很直观拟合直线一定经过点(x̄, ȳ)也就是所有样本点的重心。对w求偏导∂L/∂w -2Σx_i(y_i - wx_i - b) 0把b的表达式代进去整理之后得到w Σ(x_i - x̄)(y_i - ȳ) / Σ(x_i - x̄)²这就是最小二乘法公式的核心形态。分子是x和y的协方差乘上常数项分母是x的方差乘上常数项所以斜率的本质是“y随x的变化量占x自身变化量的比例”。2.3 一个手算实例光推公式容易晕我拿一组真实数据完整算一遍。假设我有5个样本点ixy112.6224.3336.1448.25510.4先算几个求和项Σx 12345 15Σy 2.64.36.18.210.4 31.6Σx² 1491625 55Σxy 2.68.618.332.852 114.3n 5x均值 3y均值 6.32。用前面推导的公式w (Σxy - n·x̄·ȳ) / (Σx² - n·x̄²)w (114.3 - 5×3×6.32) / (55 - 5×9)w (114.3 - 94.8) / (55 - 45) 19.5 / 10 1.95b ȳ - wx̄ 6.32 - 1.95×3 6.32 - 5.85 0.47所以拟合直线是 y 1.95x 0.47。实际生成这组数据时我用的真实关系是 y 2x 0.5加了少量高斯噪声。可以看到最小二乘估计和真实参数非常接近。2.4 为什么要写成协方差和方差的形式我在项目里经常用另一种等价形式w Cov(x, y) / Var(x)。如果你用numpy算过相关性一定见过这个结构。协方差度量的是两个变量共同变化的趋势方差度量的是x自身的离散程度。当y和x有强正向关系时协方差为正斜率就为正协方差为负斜率就为负。这个形式还有个好处它天然把量纲问题暴露出来了。如果x的单位是米y的单位是元那么斜率的单位就是元/米这和物理直觉一致。遇到特征缩放问题时记住这个形式也有助于理解为什么数据标准化之后系数的含义会发生变化。纯对比公式和实际计算我觉得数学上并不难难的是理解每一步在做什么。你只要记住三个关键量目标函数、偏导为0、重心点必定在线上后面所有线性回归的推导都能自洽。3. 升级到多元矩阵形式与正规方程3.1 设计矩阵与正规方程实际项目中不太可能只有一个自变量。预测房价既看面积又看地段还要看楼层每个自变量都对应一个系数。这时候用一元线性回归的标量公式去推不仅写法冗长逻辑也会变乱。更优雅的方式是把所有数据组织成矩阵用矩阵运算一步到位。把n个样本、d个特征写成设计矩阵X形状是n×(d1)每行是一个样本每列是一个特征另外在最前面加一列全1作为截距项。系数向量θ维度是d1包括截距和每个特征的系数。目标函数写成L(θ) ‖y - Xθ‖²展开就是 (y - Xθ)ᵀ(y - Xθ)对θ求导并令导数为0得到XᵀXθ Xᵀy这就是正规方程。只要XᵀX可逆解就是θ (XᵀX)⁻¹Xᵀy很多同学一开始看到这个公式会懵其实它就是一元场景的推广。在一元里我们求的是 w Cov(x,y)/Var(x)在多元里则是对整个协方差矩阵求逆本质上做了同样的事用自变量之间的结构去消除特征之间的互相影响。3.2 正规方程怎么来的向量求导我补充一下这个公式的推导过程能帮你彻底放下“为什么是这样”的疑问。先展开LL(θ) yᵀy - 2θᵀXᵀy θᵀXᵀXθ对θ求梯度第一项yᵀy不含θ直接消失第二项梯度是-2Xᵀy第三项用二次型求导公式梯度是2XᵀXθ。令梯度为0-2Xᵀy 2XᵀXθ 0于是XᵀXθ Xᵀy。这就是正规方程。这一段推导特别值得手写一遍。理解了它后面看岭回归的解析解 (XᵀX λI)⁻¹Xᵀy就会发现只是在对角线上加了一个λ原理完全一样。所有线性模型的解析解都是一脉相承的。3.3 实际构造设计矩阵的注意事项构造设计矩阵有几个容易踩坑的地方我吃过的亏不少第一截距项别漏。如果忘了加全1那一列拟合出来的直线会被强制穿过原点。某些场景下过原点确实合理比如x0时物理上y必须为0但如果你没想清楚就默认不过原点那模型会带一个系统性偏差。用sklearn的LinearRegression默认会拟合截距但用numpy.linalg.lstsq的时候必须自己在X上拼一列1。第二特征量纲差太大会导致数值问题。XᵀX中的元素会带上特征的单位如果面积单位是平方米数量级可能上千而楼层的数量级只有个位数两者乘起来矩阵的条件数会变得很大。条件数大意味着矩阵接近奇异求逆的结果对微小扰动极其敏感。解决办法是先做中心化或标准化再送入正规方程。第三一元里用np.polyfit很顺手但多元时不要手动用np.linalg.inv(X.TX)X.Ty更推荐用np.linalg.lstsq或np.linalg.pinv。原因后面会详细讲可以理解为数值稳定性的差异。4. 实操用Python从零实现最小二乘拟合4.1 直接按公式手写一元拟合先把前面推导的公式翻译成代码。我用numpy几行就能写完import numpy as np x np.array([1, 2, 3, 4, 5], dtypefloat) y np.array([2.6, 4.3, 6.1, 8.2, 10.4], dtypefloat) n len(x) w (n * np.sum(x * y) - np.sum(x) * np.sum(y)) / (n * np.sum(x**2) - np.sum(x)**2) b (np.sum(y) - w * np.sum(x)) / n print(fw {w:.4f}, b {b:.4f}) # w 1.9500, b 0.4700我建议新手至少手写一次这个公式而不是直接调库。手写能让你确认自己对公式的每一步理解都是对的排查异常结果时心里也更有底。写出这种代码之后再看任何框架里的LinearRegression都不过是这个过程的封装。4.2 多元线性回归与多项式拟合多元场景直接用正规方程的矩阵形式。这里我会用np.linalg.lstsq它内部通过SVD分解求解比直接求逆稳定得多x np.array([1, 2, 3, 4, 5], dtypefloat) y np.array([2.6, 4.3, 6.1, 8.2, 10.4], dtypefloat) X np.vstack([x, np.ones_like(x)]).T # 设计矩阵第一列是特征第二列是截距 theta, residuals, rank, s np.linalg.lstsq(X, y, rcondNone) print(fw {theta[0]:.4f}, b {theta[1]:.4f}) # w 1.9500, b 0.4700如果要做多项式拟合比如二次曲线只需要把设计矩阵从[x, 1]扩展成[x², x, 1]x np.array([1, 2, 3, 4, 5], dtypefloat) y np.array([3.2, 4.8, 7.1, 10.8, 15.9], dtypefloat) X_design np.column_stack([x**2, x, np.ones_like(x)]) theta np.linalg.lstsq(X_design, y, rcondNone)[0] print(fy {theta[0]:.4f}x^2 {theta[1]:.4f}x {theta[2]:.4f})当然实际项目中不用每次都手动拼设计矩阵np.polyfit和np.polyval两个函数就是专门干这个的。polyfit本质上也走的是最小二乘路线只是帮你处理了设计矩阵的构造。4.3 拟合效果怎么评估拟合做完不等于任务结束关键是要知道这条线到底拟合得好不好。最常用的指标是R²和RMSE。R²的计算公式是R² 1 - SS_res / SS_tot其中SS_res是残差平方和SS_tot是y与均值之差的平方和。R²越接近1表示拟合效果越好它衡量的是模型解释了多少y的方差。如果R²变成了负数说明你的模型比直接用y的均值还差基本等于没有拟合。代码也很直接y_hat w * x b ss_res np.sum((y - y_hat)**2) ss_tot np.sum((y - np.mean(y))**2) r2 1 - ss_res / ss_tot rmse np.sqrt(np.mean((y - y_hat)**2)) print(fR^2 {r2:.4f}, RMSE {rmse:.4f})我在实际项目里通常两个指标一起看R²判断整体解释力RMSE看预测误差的绝对值水平。如果RMSE只有几块钱哪怕R²只有0.8业务上可能也够用反之R²达到0.99但RMSE对应几十万还是要看是什么业务场景。5. 最小二乘法的局限与翻车现场5.1 异常值一个离群点就能拉偏整条线之前提过平方损失对大误差的惩罚是二次的这在异常值存在时是最致命的问题。我举个例子10个点在x1到10范围内大致线性分布但其中一个点由于记录错误y值比正常情况高了20。这个点的残差平方可能占全部残差平方的80%以上优化器为了压低这个点的误差会把整条直线往它的方向旋转。处理方式有几种第一干活前先画箱线图或散点图排查异常值第二用RANSAC这类随机抽样一致性算法它先随机抽一部分点拟合再找支持这条线的内点集合能有效抵御外点干扰第三把损失函数换成Huber损失让大误差处的梯度增长从线性变成饱和状态。我自己的习惯是数据量小的时候逐点检查数据量大的时候至少看一眼标准化残差分布。5.2 多重共线性与病态矩阵多元回归里一个非常隐蔽的坑是多重共线性。简单说就是两个自变量之间存在强相关关系比如分析房价时同时放入“房屋面积”和“房间数量”这两个特征强相关矩阵XᵀX的行列式会接近0求逆结果变得极不稳定。结果就是系数估计值方差巨大可能今天拟合出面积系数是1万明天换一批相似的数据就变成负2万。模型在训练集上R²还不错但系数完全无法解释泛化能力也很差。判断方法看特征的相关系数矩阵或者计算XᵀX的条件数。条件数超过100就要警惕超过1000基本属于病态。解决办法是加正则化。岭回归就是在XᵀX对角线上加一个λI使得矩阵更稳定代价是系数略有偏倚但换来的是方差大幅降低。Lasso则是加L1惩罚能把部分系数直接压缩成0起到特征选择的作用。在这些场景里最小二乘只是退后一步让位给它的正则化变体。5.3 过拟合与欠拟合最小二乘法只保证在给定函数形式上找到最优参数它不负责帮你选“函数形式”。如果你用最高次数50的多项式去拟合5个点误差平方和可以降到几乎为0但曲线会在点与点之间剧烈震荡离样本稍远的地方预测值完全失控。这就是过拟合。欠拟合的典型则是关系明明是非线性的你却硬用直线去拟合。比如数据呈指数增长最小二乘拟合出一条直线残差图会呈现明显的弯曲模式说明模型没有捕捉到数据的主要结构。一个通用经验用多项式拟合时从1次开始逐次增加每次比较验证集上的误差。误差先降后升的拐点就是最合适的次数。只看训练集误差的话次数越高越好那是骗人的。5.4 什么时候该换方法聊局限不是为了否定最小二乘而是为了知道它的边界。当数据量达到百万级、特征维度高解析解需要算(n×d)矩阵的乘法甚至求逆计算量会成倍增长这时候梯度下降法更合适。当特征间有复杂的非线性交互树模型和神经网络更合适。当样本量小于特征数普通最小二乘得不到唯一解岭回归和Lasso才是常规选择。我的原则是能用简单模型解释清楚的绝不上复杂模型。最小二乘永远是第一个尝试对象。它给出的结果既是基准线也是判断复杂模型是否值得的信号。6. 避开这些坑实战中的经验与建议6.1 参数求解优先用数值稳定的库函数我见过太多人用np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)来求解这是不推荐的。原因在于直接求逆会把矩阵的数值误差放大很多倍尤其当XᵀX接近奇异时误差是灾难级别的。np.linalg.lstsq内部用的是SVD或QR分解不直接求逆数值稳定性好很多。还有一个选择是np.linalg.pinv也就是伪逆。它通过SVD分解得到即使矩阵不满秩也能返回最小范数解。我自己的习惯是工程脚本里一律用lstsq如果是快速验证一两个特征也可以用polyfit但基本不用伪逆因为可读性略差容易让队友困惑。6.2 数据预处理比你想象中更重要有次处理一个回归任务特征里既有年龄30-50又有收入5万-50万还有是否已婚0或1。直接做最小二乘系数估计结果乱七八糟原因不是模型不好而是量纲差异导致数值优化和矩阵求逆都出现了不稳定。后来把所有连续特征标准化到零均值、单位方差结果立刻稳定了。中心化的作用不只是数值上的。当特征经过中心化后截距项的含义会更清晰它变成了所有特征取均值时预测值的期望。这个解释在业务汇报中非常好用。另外如果特征量纲差异极大正规方程解对特征缩放敏感先做标准化再拟合几乎总是更稳妥。6.3 检查残差别只盯着R²我踩过最大的坑就是只看R²就收工。R²高只能说明线性关系强不能说明模型假设正确。正确的做法是拟合后立刻画残差图横轴是预测值纵轴是残差。如果残差随机分布在零点附近呈无规则的带状分布说明模型基本合理。如果残差呈现明显的漏斗形随着预测值增大而发散说明方差不齐。如果残差呈现弯曲的弧线说明模型漏掉了非线性项。这些模式都不是R²能反映的。我在实际项目中每次拟合完都会花两分钟看残差图这已经成了条件反射。6.4 评估指标的选择R²、RMSE、MAE这三个指标各有适用场景。R²适合对比不同模型在同一个数据集上的解释力但它受异常值影响很大而且当数据集不同时不能直接对比。RMSE在业务上最直观单位与y一致但它对大误差更敏感如果你想更稳健地衡量典型误差MAE更好。我通常结合着看R²判断模型的相对解释力RMSE判断真实误差水平MAE判断典型误差水平。如果RMSE远大于MAE基本说明数据里存在小幅比例的较大误差需要回去查异常值。这种交叉验证的方法比单看任何一个指标都要靠谱得多。最后分享一点个人习惯我现在拿到一组数据第一件事永远是画散点图然后用最小二乘拟合一遍再画残差图。这个过程只需要几分钟但能过滤掉一半以上后续可能翻车的方向。最小二乘法最大的价值不在于它本身多么复杂而在于它给你提供了一面照妖镜数据关系是否线性、噪声是否均匀、有没有异常点在这个最简单的方法面前全都藏不住。把这一套逻辑想透了后面再学岭回归、逻辑回归、神经网络你都会发现它们只是站在同一个思想原点上的不同分支。