KNOWECON WIKI

ARTICLE

迭代重加权最小二乘法

迭代重加权最小二乘法 (Iteratively Reweighted Least Squares, IRLS) 迭代重加权最小二乘法（IRLS）是一种求解广义线性模型（GLM）和稳健回归中参数估计的数值优化算法。其核心思想是将一个非线性参数估计问题转化为一系列加权最小二乘法（Weighted Least Squares, WLS）的迭代求解，每步利用当前参数

浏览 5 更新 2025-10-26

迭代重加权最小二乘法 (Iteratively Reweighted Least Squares, IRLS)

迭代重加权最小二乘法（IRLS）是一种求解广义线性模型（GLM）和稳健回归中参数估计的数值优化算法。其核心思想是将一个非线性参数估计问题转化为一系列加权最小二乘法（Weighted Least Squares, WLS）的迭代求解，每步利用当前参数估计值重新计算各观测的权重，继而求解新的加权最小二乘问题。IRLS 在计算上等价于费舍尔得分法（Fisher Scoring），因而在 GLM 框架下具有自然的统计学基础。

核心思想

设因变量 YiY_i 与自变量 xi\mathbf{x}_i 之间的关系由参数 β\boldsymbol{\beta} 决定，且模型通过连接函数建立线性预测值与期望之间的联系。IRLS 的基本策略如下：在第 tt 次迭代中，根据当前估计 β(t)\boldsymbol{\beta}^{(t)} 构造一个"工作响应变量"（working response）zi(t)z_i^{(t)} 和一组权重 wi(t)w_i^{(t)}，使得对 (zi(t),xi,wi(t))(z_i^{(t)}, \mathbf{x}_i, w_i^{(t)}) 求解加权最小二乘即可得到更新后的 β(t+1)\boldsymbol{\beta}^{(t+1)}。这一做法避免了直接优化非二次的似然函数，将问题规约为反复求解线性方程组。

算法步骤

给定初始值 β(0)\boldsymbol{\beta}^{(0)}（通常由最小二乘法或零向量获得），IRLS 重复以下步骤直至收敛：

  1. 计算当前线性预测值：ηi(t)=xiTβ(t)\eta_i^{(t)} = \mathbf{x}_i^T \boldsymbol{\beta}^{(t)}。
  2. 通过连接函数的反函数得到拟合值：μi(t)=g−1(ηi(t))\mu_i^{(t)} = g^{-1}(\eta_i^{(t)})。
  3. 计算工作响应变量： \[ z_i^{(t)} = \eta_i^{(t)} + \frac{Y_i - \mu_i^{(t)}}{g'(\mu_i^{(t)})} \]
  4. 计算权重： \[ w_i^{(t)} = \frac{1}{\operatorname{Var}(Y_i) \cdot [g'(\mu_i^{(t)})]^2} \] 其中 Var⁡(Yi)\operatorname{Var}(Y_i) 是当前估计下的方差函数。
  5. 求解加权最小二乘问题： \[ \boldsymbol{\beta}^{(t+1)} = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{n} w_i^{(t)} \left(z_i^{(t)} - \mathbf{x}_i^T \boldsymbol{\beta}\right)^2 \] 其解为 β(t+1)=(XTW(t)X)−1XTW(t)z(t)\boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}，其中 W(t)=diag⁡(w1(t),…,wn(t))\mathbf{W}^{(t)} = \operatorname{diag}(w_1^{(t)}, \ldots, w_n^{(t)})。
  6. 检查收敛性：若 ∥β(t+1)−β(t)∥\|\boldsymbol{\beta}^{(t+1)} - \boldsymbol{\beta}^{(t)}\| 或对数似然变化量小于预设阈值，则停止。

数学推导

IRLS 可从极大似然估计的得分方程出发加以推导。设对数似然函数为 ℓ(β)\ell(\boldsymbol{\beta})，其一阶条件（得分方程）为：

∂ℓ∂β=∑i=1nYi−μiVar⁡(Yi)⋅∂μi∂ηi⋅xi=0\frac{\partial \ell}{\partial \boldsymbol{\beta}} = \sum_{i=1}^{n} \frac{Y_i - \mu_i}{\operatorname{Var}(Y_i)} \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot \mathbf{x}_i = \mathbf{0}

记 u(β)=∇βℓu(\boldsymbol{\beta}) = \nabla_{\boldsymbol{\beta}} \ell。对得分方程应用牛顿-拉夫森迭代：

β(t+1)=β(t)−[H(t)]−1u(β(t))\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - [\mathbf{H}^{(t)}]^{-1} u(\boldsymbol{\beta}^{(t)})

其中 H\mathbf{H} 为 Hessian 矩阵。在 GLM 中，若将 Hessian 替换为其期望——即费舍尔信息矩阵的负值 I(β)=−E[H]\mathcal{I}(\boldsymbol{\beta}) = -E[\mathbf{H}]——则得到费舍尔得分法。经矩阵代数整理后，费舍尔得分法的迭代公式恰好等同于上述 IRLS 的加权最小二乘形式。因此，对于典型 GLM（如逻辑回归、泊松回归），IRLS 与费舍尔得分法完全等价。

与广义线性模型的关系

在广义线性模型框架下，IRLS 是参数估计的默认算法，它与 GLM 的各组成部分紧密对应：

  • 线性预测器 ηi=xiTβ\eta_i = \mathbf{x}_i^T \boldsymbol{\beta}：定义系统化成分。
  • 连接函数 g(⋅)g(\cdot)：将均值 μi=E[Yi]\mu_i = E[Y_i] 映射到线性尺度，即 g(μi)=ηig(\mu_i) = \eta_i。
  • 方差函数 V(μi)V(\mu_i)：刻画 Var⁡(Yi)\operatorname{Var}(Y_i) 如何依赖于 μi\mu_i（如泊松分布中 V(μ)=μV(\mu) = \mu，二项分布中 V(μ)=μ(1−μ)V(\mu) = \mu(1-\mu)）。

不同的分布族和连接函数组合决定了工作响应变量和权重的具体形式。例如：

  • 逻辑回归（二项分布 + logit 连接）：g(μ)=log⁡μ1−μg(\mu) = \log\frac{\mu}{1-\mu}，V(μ)=μ(1−μ)V(\mu) = \mu(1-\mu)，权重 wi=μi(1−μi)w_i = \mu_i(1-\mu_i)。
  • 泊松回归（泊松分布 + log 连接）：g(μ)=log⁡μg(\mu) = \log\mu，V(μ)=μV(\mu) = \mu，权重 wi=μiw_i = \mu_i。
  • 线性回归（正态分布 + identity 连接）：V(μ)=1V(\mu) = 1，g′(μ)=1g'(\mu) = 1，所有权重恒为 wi=1/σ2w_i = 1/\sigma^2，此时 IRLS 一步收敛即为 OLS。

稳健回归中的应用

IRLS 在稳健回归中发挥关键作用。当数据存在异常值时，标准最小二乘法会受到严重干扰。稳健回归通过赋予异常观测较小的权重来缓解这一问题。在稳健 M 估计框架下，目标函数为：

∑i=1nρ(Yi−xiTβs)\sum_{i=1}^{n} \rho\left(\frac{Y_i - \mathbf{x}_i^T \boldsymbol{\beta}}{s}\right)

其中 ρ(⋅)\rho(\cdot) 是稳健损失函数（如Huber 损失或Tukey 双权重），ss 为尺度参数。令 ψ=ρ′\psi = \rho'，一阶条件为 ∑ψ(ri/s)xi=0\sum \psi(r_i/s) \mathbf{x}_i = \mathbf{0}。定义权重 wi=ψ(ri/s)/(ri/s)w_i = \psi(r_i/s) / (r_i/s)，则该一阶条件可重新表达为加权最小二乘的正规方程，从而直接套用 IRLS 框架迭代求解。这一过程被称为迭代重加权最小二乘稳健估计。

收敛性与收敛速度

IRLS 的收敛性取决于具体模型和数据：

  • 对于典型的 GLM（逻辑回归、泊松回归等），由于对数似然函数是凹函数，IRLS（即费舍尔得分法）通常具有全局收敛性，从任意合理的初始值出发均可收敛到全局极大值。
  • 收敛速度为二次收敛（quadratic convergence），与牛顿-拉夫森法相当。在接近最大值时，每步迭代的有效位数翻倍。
  • 对于稳健回归中的非凸损失函数（如 Tukey 双权重），IRLS 可能陷入局部极值，需要仔细选择初始值（通常先以 Huber 损失做初步估计）。
  • 实践中迭代次数通常很少——逻辑回归和泊松回归在 5-10 次以内即可收敛到机器精度。

优缺点总结

IRLS 的主要优点在于概念直观、实现简便，将复杂的非线性优化转化为反复求解熟悉的加权最小二乘问题；每次迭代只需处理一个对角权重矩阵，计算代价低；在 GLM 框架下具有严格的统计理论基础。其局限性包括：对高度共线性的数据，加权设计矩阵 XTWX\mathbf{X}^T \mathbf{W} \mathbf{X} 可能接近奇异；对于某些非标准连接函数或分布族，费舍尔得分法与牛顿-拉夫森法不再等价，IRLS 可能收敛变慢；稳健回归中需要额外估计尺度参数 ss。

IRLS 是现代统计软件的标配算法。R 语言中 \texttt{glm()} 函数、Python 中 \texttt{statsmodels} 库的 GLM 模块以及 Stata 的 \texttt{glm} 命令均默认或可选地使用 IRLS 完成参数估计。

历史背景与发展

IRLS 的思想渊源可追溯至二十世纪中叶。丹麦统计学家乔治·拉什（Georg Rasch）在 1960 年代研究项目反应模型时首次系统阐述了迭代加权的思路。随后，英国统计学家约翰·内尔德（John Nelder）和罗伯特·韦德伯恩（Robert Wedderburn）在 1972 年发表的广义线性模型奠基性论文中，将 IRLS 确立为 GLM 参数估计的标准算法，并将其与费舍尔得分法统一于同一框架之下。

在稳健统计领域，彼得·胡伯（Peter Huber）在 1964 年提出 M 估计理论时，建议用 IRLS 作为计算工具。此后，各种稳健权重函数——如Huber 权重、Tukey 双权重函数（bisquare）、Hampel 三部分递减权重等——均通过 IRLS 框架得以高效实现。

数值示例：逻辑回归

以最简单的二分类逻辑回归为例，说明 IRLS 的具体计算过程。假设有 nn 个观测，因变量 yi∈{0,1}y_i \in \{0, 1\}，自变量 xi\mathbf{x}_i。逻辑回归的连接函数为 logit：g(μ)=log⁡(μ/(1−μ))g(\mu) = \log(\mu/(1-\mu))，其导数为 g′(μ)=1/[μ(1−μ)]g'(\mu) = 1/[\mu(1-\mu)]。

在第 tt 次迭代中：

  1. 计算概率：p^i(t)=1/(1+e−xiTβ(t))\hat{p}_i^{(t)} = 1/(1 + e^{-\mathbf{x}_i^T \boldsymbol{\beta}^{(t)}})。
  2. 工作响应：zi(t)=xiTβ(t)+yi−p^i(t)p^i(t)(1−p^i(t))z_i^{(t)} = \mathbf{x}_i^T \boldsymbol{\beta}^{(t)} + \frac{y_i - \hat{p}_i^{(t)}}{\hat{p}_i^{(t)}(1-\hat{p}_i^{(t)})}。
  3. 权重：wi(t)=p^i(t)(1−p^i(t))w_i^{(t)} = \hat{p}_i^{(t)}(1-\hat{p}_i^{(t)})。
  4. 更新：β(t+1)=(XTW(t)X)−1XTW(t)z(t)\boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}。

该过程不断重复。当模型预测概率 p^i\hat{p}_i 接近 0 或 1 时，权重趋近于 0，意味着这些"几乎确定"的观测对参数更新的贡献极小——这符合直觉：分类边界附近的观测才是区分性的关键。

与其他优化算法的比较

| 方法 | 需要二阶导数 | 每步代价 | 收敛速度 | 适用场景 | |---|---|---|---|---| | IRLS / 费舍尔得分法 | 否（仅需期望信息阵） | 中等（矩阵求逆） | 二次收敛 | GLM、稳健回归 | | 牛顿-拉夫森法 | 是（Hessian 矩阵） | 较高 | 二次收敛 | 任意光滑似然 | | 梯度下降法 | 否 | 低 | 线性收敛 | 大规模、在线学习 | | BFGS / L-BFGS | 否（近似 Hessian） | 中等 | 超线性收敛 | 通用非线性优化 |

IRLS 在 GLM 场景中有独特优势：由于费舍尔信息矩阵在典型 GLM 中恒为正定，IRLS 的每次迭代均保证下降方向正确；而牛顿-拉夫森法在远离极值时可能因 Hessian 非正定而发散。但当数据维度极高（如 p>104p > 10^4）时，矩阵求逆代价过高，梯度类方法或坐标下降法（如 \texttt{glmnet} 使用的算法）更为可取。

统计推断的连接

IRLS 不仅提供了参数估计，其收敛时的权重矩阵 W\mathbf{W} 直接给出了参数估计的渐近协方差矩阵：

Cov⁡^(β^)=(XTW^X)−1\widehat{\operatorname{Cov}}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1}

这是费舍尔信息矩阵逆的估计，可用于构造Wald 检验和置信区间。此外，IRLS 的工作响应变量 ziz_i 的残差分析可用于诊断模型拟合质量：偏差残差（deviance residuals）和皮尔逊残差（Pearson residuals）均可在 IRLS 框架下自然计算。