非线性方程数值解及优化方法

对于 logistic 回归的广义线性模型, 它是由属于指数分布族的 Binomial 分布而得到的. 此 时响应的分布为 yi|xi ∼ BIN(πi, mi), i = 1, . . . , n, 且相互独立. 则对于第 i 个观测,
f (πi) = log{πi/(1 − πi)} = xTi β ≡ αi,
§4.3 无约束优化问题
45
我们希望设计一种新方法, 尽可能保持牛顿法收敛快的优点, 但又不必计算 Hessian 阵及 逆阵. 本节介绍的变尺度法 (variable metric), 其找到一个适当的 R, 能达到上述目的. 由于其 亦可看作是牛顿法的推广, 故也称之为拟牛顿法. 下面介绍其基本思想. 考虑基于 (4.20) 的迭 代公式
算法 4.3.5. 牛顿法:
(0) 令 t = 0, 给定一个初值 x0. 设定两个精度 1 和 2;
(1) 计算 ∇gt, 若 D(∇gt, 0) < 1, 则停止迭代, 取 x∗ = xt 为解. 否则, 求解线性方程组 Htd = −∇gt, 可得 dt = −H−t 1∇gt;
(2) 令 xt+1 = xt + dt;
当函数 g 是非二次的, Hessian 矩阵很可能是非正定的. 这时, 我们可能无法从 Htd = −∇gt 中解出 dt, 或者即便能求出来, 也不能保证 dt 是一个下降方向. 因此, 为了在梯度不为 零的点 xt 找到一个下降方向, 我们可以强迫 Htd = −∇gt 中的 Ht 恒取正定阵, 这样总能得 到下降方向. 其中一种常用的方法是采用 Gt = Ht + Et, 其中 Et 为一个对角阵. 这样的一个 矩阵可通过强迫矩阵的Cholesky分解而得到.
42
第4章 非线性方程数值解及优化方法
即
g;
1 2
dT
Hg
(x)d
作为
g(x
+
d)
的二次近似.
当 Hg
为正定时, dT ∇g(x) +
1 2
dT
Hg (x)d
有最小值满足
∇g(x) + Hg(x)d = 0.
则取 d = −H−g 1(x)∇g(x) 为搜索方向, 得 xt+1 = xt − H−g 1(xt)∇g(xt). 该式即为牛顿法的迭 代公式. 不难看出, 牛顿法是在 (4.20) 中 R 取为 Hessian 矩阵的下降法.
(4.23)
其中 W 为第 i 个对角元等于 miπi(1 − πi) 的对角阵. 于是, Newton法的更新方程为
βt+1 = βt − l (βt)−1l (βt) = βt + XT WtX −1 XT (y − µt) ,
(4.24) (4.25)
其中 µt 为对应着 βt 的 µ 的值, Wt 为在 µt 处取值的对角权重阵. 注意到 Hessian 阵 l (βt) 不依赖于 y. 于是, Fisher信息阵等于观测的信息量, 即 I(β) =
(3) 若达到收敛条件 (如 D(xt+1, xt) < 2 或 |g(xt+1) − g(xt)| < 2), 则停止迭代, 返回 x∗ = xt+1; 否则置 t = t + 1 转 (1).
当初始点选择的不好时, 也就是说离极小点比较远, 牛顿法则不能够保证收敛或者收敛 速度很慢, 其主要原因是可能在其中的几次迭代中未能保证下降. 为了克服这一缺点, 可考 虑改进 xt+1 的迭代公式, 即只把增量的方向设定为搜素方向 dt = −H−t 1∇gt, 引入步长因 子 λt, 并由一维搜索方法确定 λt. 即, 令 xt+1 = xt + λtdt, 其中 λt 是一维搜索得到的步长 因子, λt = arg min g(xt + λdt), 这样修改的算法称为阻尼牛顿法 (damp Newton method). 当 取 λt = 1, 即是牛顿法. 通常情况下, 阻尼牛顿法每次迭代比原始的牛顿法下降得更多, 所以 通常收敛会更快. 当然, 有时根据计算量的考虑, 我们也不一定必须求得上述一维的优化的最 优值, 只需要保证下降即可. 这时, 可使用类似前一节中的进退算法找到一个适当的值即可.
et = y − µt, zt = Xβt + Wt−1et.
则 Fisher 得分法的更新方程可以写成
βt+1 = βt + XT WtX −1 XT et = XT WtX −1 XT WtXβt + XT WtWt−1et = XT WtX −1 XT Wtzt.
(4.26)
从(4.26)可以看出, 由于 βt+1 是 zt 关于 X 的加权最小二乘的回归系数, 且其权重为 Wt (我 们通常称其为工作权矩阵working weighted matrix) 的对角元, 而称 zt 为工作响应(working response). 在每一步迭代, 都要重新计算一个新的工作响应和权向量, 且更新方程可由一个 加权最小二乘拟合得到. 对于广义线性模型, IWLS 是下将介绍的处理非线性最小二乘问题 的 Gauss-Newton 法的一种特殊情况.
n
l (β) = XT y − miπixi = XT (y − µ),
i=1
(4.22)
44
第4章 非线性方程数值解及优化方法
其中 µ = E(y) = (m1π1, . . . , mnπn)T . 其 Hessian 矩阵为
l
(β) =
∂ ∂β
(XT
(y
−
µ))
=
−
∂µ ∂β
T
X = −XT WX,
(Gt + Et) gt = xt,
(4.28)
因而有 Et gt = xt − Gt gt.
若能选取向量 ut 和 vt 使其满足规范化条件
uTt gt = vtT gt = 1,
则满足 (4.28) 式的 Et 可以简单地取以下形式:
Et = xtvtT − Gt gtuTt ,
故我们有
Gt+1 = Gt + xtvtT − Gt gtuTt .
yi ∼ BIN(πi, mi).
y1, . . . , yn 的联合似然函数可表示为
n
L(π; y) = Cmyii πiyi (1 − πi)mi−yi ,
i=1
其中 π = (π1, . . . , πn)T , y = (y1, . . . , yn)T .
将 πi =
eαi 1+eαi
代入上式, 取对数并将与参数 β 无
(3) 计算 g(xt+1) 和 ∇gt+1. 若 D(∇gt+1, 0) < 1, 或 |g(xt+1)−g(xt)| < 2, 或 D(xt+1, xt) < 3, 则停止迭代, 返回 x∗ = xt+1; 否则继续执行 (4).
(4) 判别 g(xt+1) ≥ g(xt) 是否成立, 如果成立令 x0 = xt, 转 (1); 否则继续执行 (5).
§4.3 无约束优化问题
43
例 4.3.2. (Logistic 回归模型的极大似然估计)
logistic model 是一著名的广义线性模型, 现考虑其参数的MLE. 在广义线性模型中, 响 应变量 yi 独立来自某参数为 θi 的分布 (i = 1, 2, . . . , n). 虽然不同类型的响应用不同的分布 来拟合, 但其分布通常均属于某指数分布族, 此分布族的形式为 f (y|θ) = exp{[yθ−b(θ)]/a(φ)+ c(y, φ)}, 其中 θ 为自然或典则参数, 而 φ 为散度参数. 此分布族的两个最有用的性质为：E{Y } = b (θ) 和 var{Y } = b (θ)a(φ).
关项舍去, 我们得到如下的联合对数似然函数
n
l(β) =
yi log
eαi 1 + eαi
+ (mi
− yi) log
1 1 + eαi
i=1
n
= [yiαi − mi log(1 + eαi )]
i=1
n
= yT Xβ − mi log(1 + exp{xTi β}),
i=1
其中 X = (x1, . . . , xn)T 是 n × (p + 1) 矩阵. 现考虑利用 Newton 法求最大化此似然的 β, 此时的得分函数为
• 由于
∇gt+1 − ∇gt ≈ Ht+1(xt+1 − xt),
(4.27)
所以我们自然也要求
Gt+1 gt = xt,
其中我们定义 gt = ∇gt+1 − ∇gt, xt = xt+1 − xt. 注意到这个条件有些类似于我们正割法 中的基于差分的正割公式, 我们一般称此条件为拟牛顿条件.
根据 Gt 应满足的基本条件, 特别是拟牛顿条件, 可以构造出 Gt 的各种各样的方法, 从 而形成各种变尺度算法. 下面介绍最具代表性的 DFP 算法中 Gt 的构造方法. 要求 Gt+1 = Gt + Et 满足拟牛顿条件, 得
E{−l (β)} = E{XT WX} = −l (β). 因此, 对于本例, Fisher 得分法等同于 Newton 法.
利用 Fisher 得分法来求广义线性模型的极大似然估计亦可看作是 迭代加权最小二乘(iterative weighted least quare; IWLS) 方法的应用. 记
(5) 计算 xt, gt, 及矩阵
Gt+1 = Gt +