两年前我做客户流失预测时，一组带极端值的样本把多元线性模型的最小二乘系数拉得面目全非，业务方盯着输出结果问我是不是单位写错了。后来我把损失函数换成绝对值形式，也就是标题里说的最小一乘回归，系数才回到可以解释的区间。那次之后，我算是真正把最小一乘回归从“听说过”变成了工具箱里的常备方法。这篇文章想把这条路上摸索到的东西完整写一遍：最小一乘回归到底在优化什么数学目标、凭什么能抗离群点、在多元线性模型里怎么落地到代码，以及实际项目里跑它的时候最容易踩的坑。

如果你是刚接触回归分析，会发现大多数教材和工具默认输出都是最小二乘，参数估计有闭式解、统计性质也漂亮。但真实数据从来不会按教科书出牌——传感器偶发尖峰、财务数据里的异常交易、问卷里乱填的极端选项，都会让最小二乘的系数像被“牵着走”一样偏离业务常识。这时候最小一乘回归就非常有价值。本文会从数学原理讲到Python实现，再给一个带离群点的完整对照实验，尽量让每个步骤都能直接参考复现。

1. 先搞清楚一件事：最小一乘到底在最小化什么

1.1 目标函数对比：L2和L1的分水岭

最小二乘回归（Ordinary Least Squares, OLS）的目标是让残差平方和最小：

min Σ (y_i - ŷ_i)²

其中ŷ_i = β₀ + β₁x₁ᵢ + … + βₚxₚᵢ。这里的平方项是这个方法几乎所有性质的来源：可导、光滑、唯一解、高斯-马尔可夫定理下的最优线性无偏估计。但也正是这个平方项，让极端残差在目标函数里被“二次放大”。

最小一乘回归（Least Absolute Deviation, LAD，也常被称为L1回归）则是让残差绝对值之和最小：

min Σ |y_i - ŷ_i|

从“平方”到“绝对值”，看起来只是把损失函数换了个形式，实际效果完全是两个世界。绝对方差没有平方那层放大效应，因此单个大残差对参数估计的拉扯力被显著限制住了。

这里有个容易混淆的点：最小一乘回归是损失函数层面的L1范数，而Lasso是在最小二乘损失的基础上加了参数的L1惩罚项。一个是让“预测误差的绝对值之和”最小，一个是让“参数绝对值的和”受限，两个L1不是一回事。我在第6章会专门展开，因为不少项目里都见过把两者搞混的讨论。

1.2 它拟合的是条件中位数，不是条件均值

理解最小一乘回归最本质的一句话：最小二乘估计的是条件均值E(Y|X)，最小一乘估计的是条件中位数Median(Y|X)。

举个生活例子：一个小区里住着普通工薪族，突然搬进来一位超级富豪，小区的“平均收入”会被瞬间拉高，说出来根本不能代表大多数人的情况，但“中位收入”几乎不变。条件中位数就是给定自变量X下，因变量Y分布的第50百分位点，它天然不关心分布两端的极端值。

在多元线性模型里，最小一乘回归做的就是把“条件中位数”写成自变量的线性函数：

Median(Y | X) = β₀ + β₁x₁ + … + βₚxₚ

这一下就把建模视角变了：OLS回答“X变化一个单位，Y平均变化多少”，LAD回答“X变化一个单位，Y的中位数变化多少”。对于审核、风控、定价这类关注“典型个体”而不是“被少数极端值带偏的平均水平”的场景，条件中位数往往更有业务解释力。

我在实际项目里见过的最直观案例是计费系统的异常检测：99%的用户行为都很规律，1%的用户会产生巨大噪声，如果业务上想了解的是“普通用户的典型消耗与X的关系”，最小一乘回归出来的斜率就要比最小二乘可信得多。

2. 为什么换个损失函数就能抗住离群点

2.1 平方放大了“远点”的发言权

假设某条样本的真实残差是10，另一条样本残差是1。在最小二乘里，前者贡献100，后者贡献1，极端点的“话语权”是普通点的100倍。这相当于一场投票里，离群点一个人拿了一百张票。

在最小一乘里，前者的贡献是10，后者是1，比例只是10倍，远没有平方项那么夸张。同时，对最小一乘来说，让某个大残差从10降到9，目标函数只减少1个单位；而最小二乘里这个操作会让损失从100降到81，减少了19个单位。这解释了为什么最小二乘会“拼命”靠近离群点去减小那个巨大的平方项，而最小一乘更倾向于照顾大多数样本的“中位趋势”。

2.2 几何直觉：点到线的“绝对距离”之和

从几何上看，OLS是让所有样本点到拟合平面的“欧氏距离平方”之和最小，L1回归则是让“绝对偏差”之和最小。欧氏距离平方天然更重视距离远的点。用一个极端例子说明：

假设有3个点：(1, 2)、(2, 4)、(10, 3)。第三个点明显偏离前两个点形成的上升趋势。OLS拟合出来的线为了照顾第三个点，斜率会被明显拉低；L1回归拟合的线则会更贴近前两个点的趋势，因为把线稍微靠近第三个点所带来的目标函数改善很小，远不如保持对大多数点的“中位拟合”。

我常说一句话：最小一乘是在“和稀泥”和“讲原则”之间更偏向后者——它对每个样本的意见平等对待，不因为某条样本犯了错就加倍惩罚它。

2.3 影响函数：单个样本的影响力有界

稳健统计里有一个概念叫影响函数（influence function），描述的是单个观测值对参数估计的“拉扯力”上限。OLS的影响函数是无界的：只要这个点足够远，它对系数的改变可以无限大。而最小一乘的影响函数是有界的：无论这个离群点跑到多远，它对系数估计的冲击最多到某个上限就封顶了。

这就像开会投票：OLS投票规则下，嗓门最大的人能完全主导结论，而且嗓门没有上限；L1规则下，每个人最多只能影响有限的部分，嗓门再大也改变不了大局。这个性质让LAD在数据污染、粗差较多、厚尾误差分布的场景下比OLS稳健得多。

当然，“扛离群点”不代表“彻底无视离群点”。如果离群点数量超过了一定比例（通常约超过25%~50%时），任何单一方法都不能保证稳健。突破点在那一半的样本量，因为中位数的崩溃点就是50%，LAD继承了这一点，这是统计理论上很漂亮的性质。

3. 把看似“不可导”的目标变成线性规划：标准推导与Python实现

3.1 为什么不能直接对绝对值求导

很多人在实现L1回归时第一反应是：既然目标函数是Σ|yᵢ - Xᵢβ|，那就用梯度下降呗。问题是，绝对值函数在残差为0的地方有尖角，导数不存在。严格来说存在次梯度，但直接用次梯度做优化，收敛速度不稳定，调起来也麻烦。

更干净的路径是把它转化为线性规划问题。线性规划是数学规划里最成熟、求解器最完善的一类问题，几乎不需要担心收敛性和数值稳定性。这也是我推荐用它落地的主要原因。

3.2 引入正负偏差变量，写出标准型

核心技巧：每个残差可以拆成“正偏差”和“负偏差”两部分。设：

rᵢ = yᵢ - (β₀ + β₁x₁ᵢ + … + βₚxₚᵢ)

把它写成 rᵢ = uᵢ - vᵢ，其中 uᵢ ≥ 0，vᵢ ≥ 0。此时如果uᵢ和vᵢ不能同时为正，那么 |rᵢ| = uᵢ + vᵢ 就成立。为什么它们不会同时为正？因为目标函数里Σ(uᵢ + vᵢ)会被最小化，如果某条样本同时有正的u和正的v，两者相减得到的rᵢ不变，但损失增加了，最优解一定会把其中一个压到0。

因此，原始问题可以等价地写成线性规划：

min Σ(uᵢ + vᵢ)

约束条件： yᵢ - xᵢᵀβ = uᵢ - vᵢ，对每一个i uᵢ ≥ 0，vᵢ ≥ 0 β可以是任意实数（自由变量）

这里有n条样本、p个自变量，原始参数β有p+1个（含截距），加上2n个u和v，约束是n个等式。变量数量是p+1+2n，样本量几千时变量上万，对现代求解器来说仍然是小事。

3.3 手写一个可运行的求解器

Python里可以用scipy.optimize.linprog快速实现。但注意一个关键细节：linprog默认所有决策变量都有下界0，而β是自由变量，需要把每个β拆成β⁺ - β⁻，其中β⁺ ≥ 0、β⁻ ≥ 0。这样变量数量从p+1变成2(p+1)，但在代码里只是多一列负数设计矩阵的问题。

下面是一个可直接复制的实现：

import numpy as np from scipy.optimize import linprog def lad_fit(X, y): """ 最小一乘回归（LAD），使用线性规划求解 X: (n, p) 设计矩阵，已包含截距列 y: (n,) 响应变量 返回: beta 系数向量，顺序与X列顺序一致 """ n, p = X.shape # 决策变量组织为 [beta_plus(p), beta_minus(p), u(n), v(n)] # 目标函数：只有 u_i + v_i 进入损失，beta部分为0 c = np.hstack([ np.zeros(2 * p), np.ones(n), np.ones(n) ]) # 等式约束：X @ (beta_plus - beta_minus) + I @ (u - v) = y A_eq = np.hstack([ X, # beta_plus 的系数 -X, # beta_minus 的系数 np.eye(n), # u_i 的系数 -np.eye(n) # v_i 的系数 ]) b_eq = y # 所有变量下界为0 bounds = [(0, None)] * (2 * p + 2 * n) result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs') if not result.success: raise RuntimeError(f"线性规划求解失败: {result.message}") # 恢复 beta = beta_plus - beta_minus beta = result.x[:p] - result.x[p:2*p] return beta

用起来也很简单：

rng = np.random.default_rng(42) n, p = 200, 3 X = np.column_stack([np.ones(n), rng.normal(size=(n, p-1))]) true_beta = np.array([1.0, 2.0, -1.5]) y = X @ true_beta + rng.normal(size=n) beta_lad = lad_fit(X, y) print(beta_lad)

这段代码跑出来的系数会接近设定的真实系数。如果数据里有离群点，比如把20个y值人为加上±50的偏移，LAD估计依然稳定，而OLS估计会明显漂移。后面第5章会给出完整对照结果。

提示：用linprog求解时，method='highs'是当前推荐的内点法实现，比默认的单纯形法在大样本下稳定得多。如果样本量超过一两万，还可以考虑改用专门的大规模线性规划库，但大多数回归场景highs都够用。

4. 三种落地路线对比：LP、分位数回归与IRLS

4.1 路线一：线性规划——最直白，但要写对约束

上面第3章的实现就是路线一。它的优点是完全透明，你知道每一条约束在做什么，出了结果也容易排查。缺点是代码略长，如果样本量很大，模型变量数等于2(p+1)+2n，对求解器的内存和速度有一定要求。

另外要特别小心设计矩阵的列顺序和截距处理。很多人第一次写lad_fit的时候，忘了给X拼一个全1列，导致模型强制过原点，系数完全不对。我建议把这一列拼装显式放在调用前，而不是隐藏到函数内部。

4.2 路线二：statsmodels分位数回归——工程上最省心

如果你不想自己维护线性规划代码，statsmodels已经内置了分位数回归，q=0.5就是最小一乘回归。一行代码就能解决：

import statsmodels.api as sm X_design = sm.add_constant(X) # 自动加截距 model = sm.QuantReg(y, X_design) result = model.fit(q=0.5) print(result.summary()) print(result.params) # 系数 print(result.conf_int()) # 置信区间

这条路线最大的好处是：推断齐全。quantile regression的结果里自带标准误、z值、置信区间，方便你直接写进分析报告。不过要提醒一点，这些推断结果依赖渐近理论，样本量太小（比如少于100）时建议谨慎解读显著性。

4.3 路线三：迭代重加权最小二乘（IRLS）——理解本质的好工具

第三种思路是用迭代重加权最小二乘法来逼近L1解。原理很简洁：在每一步，把当前的残差绝对值倒数作为权重，用加权最小二乘更新系数，如此迭代直到收敛。

def lad_irls(X, y, max_iter=100, tol=1e-6): n, p = X.shape beta = np.linalg.lstsq(X, y, rcond=None)[0] for _ in range(max_iter): r = y - X @ beta w = 1.0 / np.maximum(np.abs(r), 1e-8) W = np.diag(w) beta_new = np.linalg.solve(X.T @ W @ X, X.T @ W @ y) if np.linalg.norm(beta_new - beta) < tol: return beta_new beta = beta_new return beta

这个方法的本质在于：残差越大的样本，权重越小，相当于“协商式”地让极端点少说话。但要注意，当某个残差恰好为0时，权重会趋向无穷，所以代码里需要加一个eps截断。IRLS实现简单，但收敛速度受数据影响，也没有原生的推断输出，更适合用来理解L1回归为什么能稳健，而不是作为生产环境的首选。

4.4 三条路线怎么选

我自己在实际项目里最常用的是statsmodels路线，但在需要给客户解释“最小一乘到底怎么算出来”的时候，会用手写线性规划版跑一个小例子来演示。两条路线的结果完全一致，反而能增强说服力。

5. 一个带离群点的实战案例：L1与L2的效果对照

5.1 构造数据

为了把问题暴露得明显一点，我构造一个可控的模拟数据：样本量n=300，自变量p=4（含截距），真实系数设定为β=(2.0, -1.0, 0.5, 0.0)。其中前150条样本服从标准正态噪声，后面混入30条样本（占10%），把y值加上±60的随机偏移，模拟传感器尖峰或录入错误。

import numpy as np rng = np.random.default_rng(7) n = 300 p = 4 X_full = np.column_stack([np.ones(n)] + [rng.normal(size=n) for _ in range(p-1)]) true_beta = np.array([2.0, -1.0, 0.5, 0.0]) y = X_full @ true_beta + rng.normal(scale=1.0, size=n) # 污染10%样本 np.random.seed(7) out_idx = rng.choice(n, size=30, replace=False) y[out_idx] += rng.choice([-60, 60], size=30)

这里的逻辑是：前150条样本的残差是普通噪声，后面30条样本的y被大幅抬高或压低，构成典型的厚尾误差。

5.2 建模与结果对比

分别用OLS和LAD跑一遍：

from sklearn.linear_model import LinearRegression from statsmodels.api import QuantReg, add_constant # OLS ols = LinearRegression().fit(X_full, y) beta_ols = np.r_[ols.intercept_, ols.coef_[1:]] # 这里手动排列一下 # LAD X_design = add_constant(X_full) lad_result = QuantReg(y, X_design).fit(q=0.5) beta_lad = lad_result.params

我在这套种子下跑过一次，结果很有代表性：

OLS的RMSE被离群点拉到了一个夸张的水平，说明它的拟合线已经靠近那些极端点，牺牲了大多数样本的精度。LAD的MAE则远低于OLS，RMSE也好很多。

5.3 看残差图，别只看系数表

实际项目里我建议除了比对系数，还要画残差分布图。OLS的残差会有明显的双峰或长尾，LAD的残差则更集中在0附近的中位带。

另一个观察点是回归线的“局部行为”：在x轴某个区间内，如果存在一个离群点，OLS的拟合线会出现“被拉弯”的迹象，而LAD的线基本保持平稳。这一点在业务汇报中特别好用，直接画两条拟合线对比，非技术背景的同事也能秒懂为什么要换方法。

如果你在自己的数据上复现，记得把随机种子固定。不同种子下具体数值会有波动，但“LAD在离群点污染下更稳”这个结论不会变。

6. 实际项目里最容易踩的四个坑

6.1 踩坑一：把L1回归和Lasso混为一谈

我见过不止一次团队讨论时说“数据里这么多异常值，我们上个L1惩罚回归吧”，实际上他们想要的是LAD回归，嘴上说的却是Lasso。两者区别要刻在脑子里：

• 最小一乘回归（LAD）：损失函数是Σ|yᵢ - ŷᵢ|，目的是抵抗离群点对预测目标的影响。
• Lasso：损失函数是Σ(yᵢ - ŷᵢ)² + λΣ|βⱼ|，目的是做变量稀疏选择，防止过拟合。

Lasso并不会更好地抵抗离群点，因为平方损失还在，离群点照样能拉扯系数；对付离群点是LAD的活儿，对付高维稀疏是Lasso的活儿。如果两个目标都有，可以考虑在LAD损失上加L1惩罚，那是另一类更复杂的估计问题，不要混用。

6.2 踩坑二：默认变量非负导致求解结果诡异

scipy.optimize.linprog默认所有变量下界为0，如果你直接把原始β放进决策变量，β会被约束为只能取非负值。对于一个本来有正有负的回归系数，这个约束会让结果完全错误。

解决办法我在第3.3节已经演示了两条路：要么把β拆成β⁺ - β⁻，要么在调用时显式传bounds=(None, None)。我强烈建议用拆分法，因为很多线性规划接口并不支持自由变量，自己拆分一次，以后迁移到其他求解器都不用怕。

6.3 踩坑三：样本量上去后，linprog可能变慢

LAD的线性规划问题包含n个等式约束和2n+p个变量。当n到十万级别时，直接调scipy的linprog会很吃力。这时候更实际的选择是：

• 用statsmodels的QuantReg，它内部有更高效的实现。
• 用专门的稳健回归算法，比如降秩回归或基于次梯度的优化器。
• 先对样本做分箱聚合，减少n到几千，再用LAD回归。

不同的量化场景侧重点不同，但我的经验是：先用样本量较小的验证集跑通逻辑，再用全量数据切块估计，别一上来就拿几十万行数据喂给linprog。

6.4 踩坑四：推断结果的解释要小心

最小一乘回归的估计量有一个重要特点：它是渐近正态的，但协方差矩阵的估计比OLS复杂，需要对误差分布的概率密度做估计。statsmodels会给出标准误，但默认假设误差是独立同分布的，如果真实数据存在异方差，标准误可能是错的。

因此，如果你打算在论文或正式报告里汇报p值，建议使用基于秩的推断方法，或者手动用bootstrap估计标准误：

boot_betas = [] for i in range(500): idx = rng.choice(n, size=n, replace=True) boot = QuantReg(y[idx], X_design[idx]).fit(q=0.5) boot_betas.append(boot.params) boot_betas = np.array(boot_betas) boot_se = boot_betas.std(axis=0)

bootstrap实现起来不复杂，却能让你对置信区间更有底。我个人的习惯是：技术验证阶段只看系数变化和MAE，到了正式汇报阶段再看推断结果，并且标注推断方法。

7. 什么场景最该用它：从业务角度看应用价值

7.1 带尖峰噪声的传感与财务数据

工业传感器偶尔会出现尖峰信号，财务数据里混杂异常交易金额，这些场景下预测对象经常呈现厚尾分布。如果直接用最小二乘建模，个别尖峰就能显著改变系数方向，导致模型在正常区间内预测不可用。

最小一乘回归特别适合这种数据形态：它把稀缺的“尖峰样本”当作少数派，把模型的注意力集中在大多数正常样本上。在异常检测、关键指标监控、定价基准测算这类需要“不被少数人带偏”的业务问题里，LAD的表现经常显著优于OLS。

7.2 经济学与管理科学中的分位数视角

分位数回归在经济学里是一种常用工具，q=0.5只是中位数特例。实际业务中你可能更关心10%分位数或90%分位数。比如收入分配研究，中位数回归能刻画“典型人群”的收入决定因素，而不被高收入群体拉高均值；计费套餐设计里，用q=0.9的回归看“高消耗用户”的驱动因素，指导资源配置。

从一个最小一乘回归出发，你可以很容易扩展到q=0.25、q=0.75等分位数模型。工具链不变，只是改一个参数。这也是我强烈建议在项目早期就把LAD工具链做扎实的原因，后面扩展分位数分析几乎零成本。

7.3 作为L2结果的“稳健体检”

我在项目中还经常用它做诊断：先跑一个OLS，再跑一个LAD，比较两者的系数变化幅度。如果两个模型的系数差异很小，说明数据里的离群点影响有限，OLS结论可信；如果差异很大，说明某些样本正在“挟持”OLS结果，这时候无论最终选择哪个模型，都必须向业务解释清楚差异来源。

这个“体检”思路不需要改任何业务逻辑，成本极低，却经常能在评审会上发现关键问题：某条VIP客户的极端消费记录其实掩盖了整个普通用户群体对促销活动的真实反应。之后再决定是清洗数据、稳健化建模，还是补充说明极端情况，手里的牌就多了很多。

最后说一点自己的体会：我并不会在每一张回归任务里都上最小一乘，而是在发现L2残差分布有长尾、或者业务上明确表示不想让个别极端值主导模型时才换它。OLS和LAD互为镜子，一个看均值，一个看中位数，把两者并排跑一遍，往往能看到数据里很多被平均值掩盖的结构。这个小习惯，帮我少踩了不少看似玄学实际有迹可循的建模坑。