在下文中,我们将提供一些关于数值优化算法的笔记。由于存在众多方法 ,我们将只关注所谓的梯度方法 。我们之所以将此类方法作为考虑数值优化算法的自然起点,主要有以下两点原因。一方面,这些方法在该领域确实是主力军,因此它们在实践中的频繁使用证明了它们在这里的覆盖范围。另一方面,这种方法非常直观,因为它在某种程度上从众所周知的最优值性质 自然地推导出。我们将特别关注该类方法的三个例子:牛顿方法 、最速下降法 以及可变度量方法 ,其中包括拟牛顿方法 。
在开始之前,我们要强调,似乎不存在“唯一”算法,而特定算法的性能总是取决于待解决的特定问题。因此,经验和“试错”在应用工作中非常重要。为了阐明这一点,我们将提供几个应用程序,其中可以以图形方式比较不同算法的性能。此外,最后还提供了一个关于最大似然估计 的具体示例。尤其是对于统计学家和计量经济学家 来说,最大似然估计可能是在实践中必须依赖数值优化算法的最重要示例。
任何数值优化算法都需要解决找到函数的“可观察”性质的问题,以便计算机程序知道何时达到解。由于我们处理的是优化问题,因此两个众所周知的结论似乎是进行这种性质推导的合理起点。
如果 f 可微,并且 x ⋆ {\displaystyle x^{\star }} 是(局部)最小值,那么
( 1 a ) D f ( x ⋆ ) = 0 {\displaystyle (1a)\quad Df(x^{\star })=0}
即雅可比矩阵 D f ( x ) {\displaystyle Df(x)} 等于零
以及
如果 f 二次可微,并且 x ⋆ {\displaystyle x^{\star }} 是(局部)最小值,那么
( 1 b ) x T D 2 f ( x ⋆ ) x ≥ 0 {\displaystyle (1b)\quad x^{T}D^{2}f(x^{\star })x\geq 0}
即海森矩阵 D 2 f ( x ) {\displaystyle D^{2}f(x)} 是正半定 的。
在下文中,我们将始终用 x ⋆ {\displaystyle x^{\star }} 表示最小值。尽管这两个条件似乎代表了有助于找到最优值 x ⋆ {\displaystyle x^{\star }} 的陈述,但它们给出了 x ⋆ {\displaystyle x^{\star }} 作为函数 f {\displaystyle f} 的最优值的含义,存在一个小问题。但为了我们的目的,我们需要相反的含义,即我们最终想要得到一个形式为:"如果某个条件 g ( f ( x ⋆ ) ) {\displaystyle g(f(x^{\star }))} 为真,则 x ⋆ {\displaystyle x^{\star }} 是一个最小值" 的陈述。但上面的两个条件显然不足以实现这一点(例如考虑 f ( x ) = x 3 {\displaystyle f(x)=x^{3}} 的情况,其中 D f ( 0 ) = D 2 f ( 0 ) = 0 {\displaystyle Df(0)=D^{2}f(0)=0} ,但 x ⋆ ≠ 0 {\displaystyle x^{\star }\neq 0} )。因此,我们必须查看 x ⋆ {\displaystyle x^{\star }} 的整个邻域,如以下检测最优值的充分条件中所述
如果 D f ( x ⋆ ) = 0 {\displaystyle Df(x^{\star })=0} 且 x T D 2 f ( z ) x ≥ 0 , ∀ x ∈ R n {\displaystyle x^{T}D^{2}f(z)x\geq 0,\forall x\in \mathbb {R} ^{n}} 且 z ∈ B ( x ⋆ , δ ) {\displaystyle z\in {\mathcal {B}}(x^{\star },\delta )} ,则: x ⋆ {\displaystyle x^{\star }} 是一个局部最小值。
证明:对于 x ∈ B ( x ⋆ , δ ) {\displaystyle x\in {\mathcal {B}}(x^{\star },\delta )} 令 z = x ⋆ + t ( x − x ⋆ ) ∈ B {\displaystyle z=x^{\star }+t(x-x^{\star })\in {\mathcal {B}}} 。泰勒近似 得到: f ( x ) − f ( x ⋆ ) = 0 + 1 2 ( x − x ⋆ ) T D 2 f ( z ) ( x − x ⋆ ) ≥ 0 {\displaystyle f(x)-f(x^{\star })=0+{\frac {1}{2}}(x-x^{\star })^{T}D^{2}f(z)(x-x^{\star })\geq 0} ,其中 B ( x ⋆ , δ ) {\displaystyle {\mathcal {B}}(x^{\star },\delta )} 表示以 x ⋆ {\displaystyle x^{\star }} 为中心的开球,即 B ( x ⋆ , δ ) = { x : | | x − x ⋆ | | ≤ δ } {\displaystyle {\mathcal {B}}(x^{\star },\delta )=\{x:||x-x^{\star }||\leq \delta \}} 对于 δ > 0 {\displaystyle \delta >0} 。
与上述两个条件相反,此条件足以检测最优值——考虑两个简单的例子
f ( x ) = x 3 {\displaystyle f(x)=x^{3}} ,其中 D f ( x ⋆ = 0 ) = 0 {\displaystyle Df(x^{\star }=0)=0} 但 x T D 2 f ( z ) x = 6 z x 2 ≱ 0 ( e . g . z = − δ 2 ) {\displaystyle x^{T}D^{2}f(z)x=6zx^{2}\not \geq 0\quad (e.g.\,\,z=-{\frac {\delta }{2}})}
以及
f ( x ) = x 4 {\displaystyle f(x)=x^{4}} ,其中 D f ( x ⋆ = 0 ) = 0 {\displaystyle Df(x^{\star }=0)=0} 并且 x T D 2 f ( z ) x = 12 z 2 x 2 ≥ 0 ∀ z {\displaystyle x^{T}D^{2}f(z)x=12z^{2}x^{2}\geq 0\quad \forall z} 。
牢记这一小点,我们现在可以转向数值优化程序。
以下所有算法都将依赖于以下假设
(A1) 集合 N ( f , f ( x ( 0 ) ) = { x ∈ R n | f ( x ) ≤ f ( x ( 0 ) ) } {\displaystyle N(f,f(x^{(0)})=\{x\in \mathbb {R} ^{n}|f(x)\leq f(x^{(0)})\}} 是 紧凑 的。
其中 x ( 0 ) {\displaystyle x^{(0)}} 是算法给定的某个初始值。这个假设的意义在于 *魏尔斯特拉斯定理*,它指出每个紧凑集都包含其 *上确界 * 和其 *下确界 *。因此,*(A1)* 保证了在 N ( f , f ( x ( 0 ) ) {\displaystyle N(f,f(x^{(0)})} 中存在某个解。在这个全局最小值 x ⋆ {\displaystyle x^{\star }} 处,当然有 D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 。因此,考虑到上述讨论,优化问题基本上可以归结为求解方程组 D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 。
这种方法的问题当然很普遍,因为 D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 对于 极大值和鞍点 也是成立的。因此,好的算法应该确保极大值和鞍点都被排除为可能的解。极大值可以通过要求 f ( x ( k + 1 ) ) < f ( x ( k ) ) {\displaystyle f(x^{(k+1)}) ,即我们限制自己于一个 序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} ,使得函数值在每一步都减小。问题当然是在这是否总是可能的。幸运的是,这是可能的。其基本见解如下。在构建映射 x ( k + 1 ) = φ ( x ( k ) ) {\displaystyle x^{(k+1)}=\varphi (x^{(k)})} (即,从 x ( k ) {\displaystyle x^{(k)}} 到 x ( k + 1 ) {\displaystyle x^{(k+1)}} 的规则)时,我们有两个自由度,即
方向 d ( k ) {\displaystyle d^{(k)}} 和
步长 σ ( k ) {\displaystyle \sigma ^{(k)}} 。
因此,我们可以选择我们想要移动的方向以到达 x ( k + 1 ) {\displaystyle x^{(k+1)}} 以及 这种移动需要多远。所以,如果我们选择 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 以“正确的方式”,我们可以有效地确保函数值下降。以下将提供此推理的正式表示。
引理:如果 d ( k ) ∈ R n {\displaystyle d^{(k)}\in \mathbb {R} ^{n}} 并且 D f ( x ) T d ( k ) < 0 {\displaystyle Df(x)^{T}d^{(k)}<0} ,那么: ∃ σ ¯ > 0 {\displaystyle \exists {\bar {\sigma }}>0} 使得
f ( x + σ ( k ) d ( k ) ) < f ( x ) ∀ σ ∈ ( 0 , σ ¯ ) {\displaystyle f(x+\sigma ^{(k)}d^{(k)})<f(x)\quad \forall \sigma \in (0,{\bar {\sigma }})}
证明:因为 D f ( x ) T d ( k ) < 0 {\displaystyle Df(x)^{T}d^{(k)}<0} 并且 D f ( x ) T d ( k ) = lim σ → 0 f ( x + σ ( k ) d ( k ) ) − f ( x ) σ ( k ) {\displaystyle Df(x)^{T}d^{(k)}=\lim _{\sigma \rightarrow 0}{\frac {f(x+\sigma ^{(k)}d^{(k)})-f(x)}{\sigma ^{(k)}}}} ,因此得出 f ( x + σ ( k ) d ( k ) ) < f ( x ) {\displaystyle f(x+\sigma ^{(k)}d^{(k)})<f(x)} 对于足够小的 σ ( k ) {\displaystyle \sigma ^{(k)}} 。
满足此条件的方向向量 d ( k ) {\displaystyle d^{(k)}} 被称为下降方向 。在实践中,此引理允许我们使用以下程序来数值求解优化问题。
1. 定义序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} ,通过 x ( k + 1 ) = x ( k ) + σ ( k ) d ( k ) {\displaystyle x^{(k+1)}=x^{(k)}+\sigma ^{(k)}d^{(k)}} 递归地定义。
2. 从点 x ( k ) {\displaystyle x^{(k)}} 的局部信息中选择方向 d ( k ) {\displaystyle d^{(k)}} 。
3. 选择一个步长 σ ( k ) {\displaystyle \sigma ^{(k)}} ,确保算法的收敛 。
4. 如果 | f ( x ( k + 1 ) ) − f ( x ( k ) ) | < ϵ {\displaystyle |f(x^{(k+1)})-f(x^{(k)})|<\epsilon } ,其中 ϵ > 0 {\displaystyle \epsilon >0} 是某个选定的最小值容差,则停止迭代。
这个过程已经暗示了 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 的选择不是独立的,而是相互依赖的。特别要注意,即使该方法是下降方法(即 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 是根据引理 1 选择的),也不能保证收敛到最小值。乍一看,这似乎有点令人费解。如果我们找到一个序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} ,使得函数值在每一步都下降,人们可能会认为,在某个阶段,即在k 趋于无穷大的极限中,我们应该找到解决方案。为什么情况并非如此,可以从借鉴于 W. Alt (2002, p. 76) 的以下例子中看出。
考虑以下例子,它虽然明显是下降的,但并不收敛。令准则函数由下式给出
f ( x ) = x 2 {\displaystyle f(x)=x^{2}} ,令初始值为 x ( 0 ) = 1 {\displaystyle x^{(0)}=1} ,考虑一个(常数)方向向量 d ( k ) = − 1 {\displaystyle d^{(k)}=-1} ,并选择步长为 σ ( k ) = ( 1 2 ) k + 2 {\displaystyle \sigma ^{(k)}=({\frac {1}{2}})^{k+2}} 。因此,序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 的递归定义如下:
( 2 ) x ( k + 1 ) = x ( k ) + ( 1 2 ) k + 2 ( − 1 ) = x ( k − 1 ) − ( 1 2 ) k + 1 − ( 1 2 ) k + 2 = x ( 0 ) − ∑ j = 0 k ( 1 2 ) j + 2 . {\displaystyle (2)\quad x^{(k+1)}=x^{(k)}+({\frac {1}{2}})^{k+2}(-1)=x^{(k-1)}-\left({\frac {1}{2}}\right)^{k+1}-\left({\frac {1}{2}}\right)^{k+2}=x^{(0)}-\sum _{j=0}^{k}\left({\frac {1}{2}}\right)^{j+2}.}
注意 x ( k ) > 0 ∀ k {\displaystyle x^{(k)}>0\,\,\forall \,\,k} ,因此 f ( x ( k + 1 ) ) < f ( x ( k ) ) ∀ k {\displaystyle f(x^{(k+1)})<f(x^{(k)})\,\,\forall \,\,k} ,因此该方法显然是下降法。然而,我们发现
( 3 ) lim k → ∞ x ( k ) = lim k → ∞ x ( 0 ) − ∑ j = 0 k − 1 ( 1 2 ) j + 2 = lim k → ∞ 1 − 1 4 ( 1 − ( 1 2 ) k 1 2 ) = lim k → ∞ 1 2 + ( 1 2 ) k + 1 = 1 2 ≠ 0 = x ⋆ . {\displaystyle (3)\quad \lim _{k\rightarrow \infty }x^{(k)}=\lim _{k\rightarrow \infty }x^{(0)}-\sum _{j=0}^{k-1}\left({\frac {1}{2}}\right)^{j+2}=\lim _{k\rightarrow \infty }1-{\frac {1}{4}}\left({\frac {1-\left({\frac {1}{2}}\right)^{k}}{\frac {1}{2}}}\right)=\lim _{k\rightarrow \infty }{\frac {1}{2}}+\left({\frac {1}{2}}\right)^{k+1}={\frac {1}{2}}\neq 0=x^{\star }.}
这种不收敛的原因要归咎于步长 σ ( k ) {\displaystyle \sigma ^{(k)}} 下降过快。对于较大的 k ,步长 x ( k + 1 ) − x ( k ) {\displaystyle x^{(k+1)}-x^{(k)}} 会变得非常小,从而导致收敛无法实现。因此,我们必须将步长与下降方向 d ( k ) {\displaystyle d^{(k)}} 联系起来。
这种联系的显而易见的思路是要求实际下降与一阶近似成正比,即选择 σ ( k ) {\displaystyle \sigma ^{(k)}} 使得存在一个常数 c 1 > 0 {\displaystyle c_{1}>0} ,满足
( 4 ) f ( x ( k ) + σ ( k ) d ( k ) ) − f ( x ( k ) ) ≤ c 1 σ ( k ) D ( f ( x ( k ) ) ) d ( k ) < 0. {\displaystyle (4)\quad f(x^{(k)}+\sigma ^{(k)}d^{(k)})-f(x^{(k)})\leq c_{1}\sigma ^{(k)}D(f(x^{(k)}))d^{(k)}<0.}
请注意,我们仍然只关注下降方向,因此 D f ( x ( k ) ) T d ( k ) < 0 {\displaystyle Df(x^{(k)})^{T}d^{(k)}<0} ,如上文引理 1 中所述。因此, N ( f , f ( x ( k ) ) ) {\displaystyle N(f,f(x^{(k)}))} 的紧致性意味着 LHS 的 收敛 ,并且由 (4) 得出
( 5 ) lim k → ∞ σ ( k ) D ( f ( x ( k ) ) ) d ( k ) = 0. {\displaystyle (5)\quad \lim _{k\rightarrow \infty }\sigma ^{(k)}D(f(x^{(k)}))d^{(k)}=0.}
最后,我们希望选择一个序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 使得 lim k → ∞ D ( f ( x ( k ) ) ) = 0 {\displaystyle \lim _{k\rightarrow \infty }D(f(x^{(k)}))=0} ,因为这正是我们想要解决的必要的一阶条件。在什么条件下,(5) 实际上意味着 lim k → ∞ D ( f ( x ( k ) ) ) = 0 {\displaystyle \lim _{k\rightarrow \infty }D(f(x^{(k)}))=0} ?首先,步长 σ ( k ) {\displaystyle \sigma ^{(k)}} 不能太快地趋于零。这正是我们在上面示例中遇到的情况。因此,要求步长从下方有界似乎是合理的,即:
( 6 ) σ ( k ) ≥ − c 2 D f ( x ( k ) ) T d ( k ) | | d ( k ) | | 2 > 0 {\displaystyle (6)\quad \sigma ^{(k)}\geq -c_{2}{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||^{2}}}>0}
对于某个常数 c 2 > 0 {\displaystyle c_{2}>0} 。将 (6) 代入 (5) 最终得到
( 7 ) f ( x ( k ) + σ ( k ) d ( k ) ) − f ( x ( k ) ) ≤ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 , c = c 1 c 2 {\displaystyle (7)\quad f(x^{(k)}+\sigma ^{(k)}d^{(k)})-f(x^{(k)})\leq -c\left({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}\right)^{2},\quad c=c_{1}c_{2}}
这里, N ( f , f ( x ( k ) ) ) {\displaystyle N(f,f(x^{(k)}))} 的 紧致性 保证了左侧的 收敛 ,因此
( 8 ) lim k → ∞ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 = lim k → ∞ D f ( x ( k ) ) T d ( k ) | | d ( k ) | | = 0 {\displaystyle (8)\quad \lim _{k\rightarrow \infty }-c({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}})^{2}=\lim _{k\rightarrow \infty }{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}=0}
满足 (4) 和 (6) 的步长被称为 *有效步长*,并用 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 表示。条件 (6) 的重要性将在以下示例 1 的续篇中说明。
请注意,正是 (6) 的失效导致了示例 1 未收敛。将示例中的步长代入 (6) 得到
( 6.1 ) σ ( k ) = ( 1 2 ) ( k + 2 ) ≥ − c 2 2 x ( k ) ( − 1 ) 1 = c 2 ⋅ 2 ( 1 2 + ( 1 2 ) k + 1 ) ⇔ 1 4 ( 1 + 2 ( k ) ) ≥ c 2 > 0 {\displaystyle (6.1)\quad \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{(k+2)}\geq -c_{2}{\frac {2x^{(k)}(-1)}{1}}=c_{2}\cdot 2({\frac {1}{2}}+\left({\frac {1}{2}})^{k+1}\right)\Leftrightarrow {\frac {1}{4(1+2^{(k)})}}\geq c_{2}>0}
因此,对于所有k ,没有 常数 c 2 > 0 {\displaystyle c_{2}>0} 满足此不等式,如 (6) 所要求的那样。因此,步长没有从下界约束,并且下降得太快。为了真正认识 (6) 的重要性,让我们稍微修改一下示例,并假设 σ ( k ) = ( 1 2 ) k + 1 {\displaystyle \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{k+1}} 。然后我们发现
( 6.2 ) lim k → ∞ x ( k + 1 ) = lim k → ∞ x ( 0 ) − 1 2 ∑ i ( 1 2 ) i = lim k → ∞ ( 1 2 ) k + 1 = 0 = x ⋆ , {\displaystyle (6.2)\quad \lim _{k\rightarrow \infty }x^{(k+1)}=\lim _{k\rightarrow \infty }x^{(0)}-{\frac {1}{2}}\sum _{i}\left({\frac {1}{2}}\right)^{i}=\lim _{k\rightarrow \infty }\left({\frac {1}{2}}\right)^{k+1}=0=x^{\star },}
即,收敛 确实发生了。此外,认识到此示例实际上满足条件 (6),因为
( 6.3 ) σ ( k ) = ( 1 2 ) ( k + 1 ) ≥ − c 2 2 x ( k ) ( − 1 ) 1 = c 2 ⋅ 2 ( 1 2 ) k ⇔ 1 4 ≥ c 2 > 0. {\displaystyle (6.3)\quad \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{(k+1)}\geq -c_{2}{\frac {2x^{(k)}(-1)}{1}}=c_{2}\cdot 2\left({\frac {1}{2}}\right)^{k}\Leftrightarrow {\frac {1}{4}}\geq c_{2}>0.}
我们已经论证过, σ ( k ) {\displaystyle \sigma ^{(k)}} 和 d ( k ) {\displaystyle d^{(k)}} 的选择是相互交织的。因此,选择“正确”的 d ( k ) {\displaystyle d^{(k)}} 总是取决于相应的步长 σ ( k ) {\displaystyle \sigma ^{(k)}} 。那么在这个语境中,“正确”意味着什么?在上面,我们在方程 (8) 中展示了,选择一个有效步长 意味着
( 8 ′ ) lim k → ∞ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 = lim k → ∞ D f ( x ( k ) ) T d ( k ) | | d ( k ) | | = 0. {\displaystyle (8')\quad \lim _{k\rightarrow \infty }-c({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}})^{2}=\lim _{k\rightarrow \infty }{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}=0.}
因此,“正确”的方向向量将保证 (8') 意味着
( 9 ) lim k → ∞ D f ( x ( k ) ) = 0 {\displaystyle (9)\quad \lim _{k\rightarrow \infty }Df(x^{(k)})=0}
因为 (9) 是所选序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 收敛的条件。因此,让我们探索哪些方向可以选择以产生 (9)。假设步长 σ k {\displaystyle \sigma _{k}} 是 _有效的_ 并定义
( 10 ) β ( k ) = D f ( x ( k ) ) T d ( k ) | | D f ( x ( k ) ) | | | | d ( k ) | | ⇔ β ( k ) | | D f ( x ( k ) ) | | = D f ( x ( k ) ) T d ( k ) | | d ( k ) | | {\displaystyle (10)\quad \beta ^{(k)}={\frac {Df(x^{(k)})^{T}d^{(k)}}{||Df(x^{(k)})||||d^{(k)}||}}\quad \Leftrightarrow \quad \beta ^{(k)}||Df(x^{(k)})||={\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}}
根据 (8') 和 (10),我们知道
( 11 ) lim k → ∞ β ( k ) | | D f ( x ( k ) ) | | = 0. {\displaystyle (11)\quad \lim _{k\rightarrow \infty }\beta ^{(k)}||Df(x^{(k)})||=0.}
因此,如果我们从下方对 β ( k ) {\displaystyle \beta ^{(k)}} 进行约束(即 β ( k ) ≤ − δ < 0 {\displaystyle \beta ^{(k)}\leq -\delta <0} ),(11) 意味着
( 12 ) lim k → ∞ β ( k ) | | D f ( x ( k ) ) | | = lim k → ∞ | | D f ( x ( k ) ) | | = lim k → ∞ D f ( x ( k ) ) = 0 , {\displaystyle (12)\quad \lim _{k\rightarrow \infty }\beta ^{(k)}||Df(x^{(k)})||=\lim _{k\rightarrow \infty }||Df(x^{(k)})||=\lim _{k\rightarrow \infty }Df(x^{(k)})=0,}
其中 (12) 给出了序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 收敛到解 x ⋆ {\displaystyle x^{\star }} 的条件。由于 (10) 通过 β ( k ) {\displaystyle \beta ^{(k)}} 隐式定义了方向向量 d ( k ) {\displaystyle d^{(k)}} ,对 β ( k ) {\displaystyle \beta ^{(k)}} 的要求直接转化为对 d ( k ) {\displaystyle d^{(k)}} 的要求。
当考虑对 β ( k ) {\displaystyle \beta ^{(k)}} 的条件时,很明显“梯度方法”一词从何而来。当 β ( k ) {\displaystyle \beta ^{(k)}} 由以下给出时
β k = D ( f ( x ) ) d ( k ) | | D f ( x ( k ) ) | | | | d ( k ) | | = cos ( D f ( x ( k ) ) , d ( k ) ) {\displaystyle \beta _{k}={\frac {D(f(x))d^{(k)}}{||Df(x^{(k)})||||d^{(k)}||}}=\cos(Df(x^{(k)}),d^{(k)})}
我们有以下结果
鉴于 σ ( k ) {\displaystyle \sigma ^{(k)}} 是有效选择的,并且 d ( k ) {\displaystyle d^{(k)}} 满足
( 13 ) cos ( D f ( x ( k ) ) , d ( k ) ) = β k ≤ − δ < 0 {\displaystyle (13)\quad \cos(Df(x^{(k)}),d^{(k)})=\beta _{k}\leq -\delta <0}
我们有
( 14 ) lim k → ∞ D f ( x ( k ) ) → 0 {\displaystyle (14)\quad \lim _{k\rightarrow \infty }Df(x^{(k)})\rightarrow 0}
因此:如果 x ( k ) {\displaystyle x^{(k)}} 处的负梯度和方向 d ( k ) {\displaystyle d^{(k)}} 之间的角度始终小于直角,则会发生收敛。依赖于满足 (13) 的 d ( k ) {\displaystyle d^{(k)}} 的方法被称为梯度方法。
换句话说:只要不朝与梯度垂直的方向移动,并且如果步长选择有效,梯度方法 就可以保证收敛到解 x ⋆ {\displaystyle x^{\star }} 。
现在,让我们探索这一类别中的三种特定算法,它们在各自选择的 d ( k ) {\displaystyle d^{(k)}} 上有所不同。
牛顿法 是该领域中最受欢迎的方法。它是一种众所周知的方法,用于求解所有类型方程的根 ,因此很容易应用于优化问题。牛顿法的主要思想是线性化方程组,得到
( 15 ) g ( x ) = g ( x ^ ) + D g ( x ^ ) T ( x − x ^ ) = 0. {\displaystyle (15)\quad g(x)=g({\hat {x}})+Dg({\hat {x}})^{T}(x-{\hat {x}})=0.}
(15) 可以很容易地求解出 x ,因为解只是(假设 D g ( x ^ ) T {\displaystyle Dg({\hat {x}})^{T}} 是非奇异 的)
( 16 ) x = x ^ − [ D g ( x ^ ) T ] − 1 g ( x ^ ) . {\displaystyle (16)\quad x={\hat {x}}-[Dg({\hat {x}})^{T}]^{-1}g({\hat {x}}).}
为了我们的目的,我们选择 g ( x ) {\displaystyle g(x)} 为梯度 D f ( x ) {\displaystyle Df(x)} ,并得出
( 17 ) d N ( k ) = x ( k + 1 ) − x ( k ) = − [ D 2 f ( x ( k ) ) ] − 1 D f ( x ( k ) ) {\displaystyle (17)\quad d_{N}^{(k)}=x^{(k+1)}-x^{(k)}=-[D^{2}f(x^{(k)})]^{-1}Df(x^{(k)})}
其中 d N ( k ) {\displaystyle d_{N}^{(k)}} 是所谓的牛顿方向 。
分析(17)揭示了牛顿法的主要性质
如果 D 2 f ( x ( k ) ) {\displaystyle D^{2}f(x^{(k)})} 是正定 ,则 d N k {\displaystyle d_{N}^{k}} 是按照引理1 意义上的下降方向。
牛顿法 使用一阶和 二阶导数的局部信息来计算 d N k {\displaystyle d_{N}^{k}} 。
( 18 ) x ( k + 1 ) = x ( k ) + d N ( k ) {\displaystyle (18)\quad x^{(k+1)}=x^{(k)}+d_{N}^{(k)}}
牛顿法 使用 σ ( k ) = 1 {\displaystyle \sigma ^{(k)}=1} 的固定步长。因此,牛顿法不 一定是一种按照引理1 意义上的下降方法。原因是固定步长 σ ( k ) = 1 {\displaystyle \sigma ^{(k)}=1} 可能大于引理1 中给出的临界步长 σ k ¯ {\displaystyle {\bar {\sigma _{k}}}} 。下面我们以Rosenbrock函数为例说明牛顿法不是下降方法。
该方法可能很耗时,因为在每一步k 中计算 [ D 2 f ( x ( k ) ) ] − 1 {\displaystyle [D^{2}f(x^{(k)})]^{-1}} 很麻烦。在应用工作中,可以考虑使用近似方法。例如,可以每s 步更新一次Hessian,或者可以依赖于局部近似 。这被称为拟牛顿法 ,将在关于可变度量法 的部分中讨论。
为了确保该方法为下降方法,可以使用一个有效的步长 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} ,并设置
( 19 ) x ( k + 1 ) = x ( k ) − σ E ( k ) d N ( k ) = x ( k ) − σ E ( k ) [ D 2 f ( x k ) ] − 1 D f ( x ( k ) ) {\displaystyle (19)\quad x^{(k+1)}=x^{(k)}-\sigma _{E}^{(k)}d_{N}^{(k)}=x^{(k)}-\sigma _{E}^{(k)}[D^{2}f(x^{k})]^{-1}Df(x^{(k)})}
另一种常用的方法是最速下降法 。这种方法的思想是选择方向 d ( k ) {\displaystyle d^{(k)}} ,使得函数值f 的下降最大。尽管这个过程乍一看很合理,但它实际上使用的信息比牛顿法 少,因为它忽略了Hessian关于函数曲率的信息。特别是在下面的应用中,我们将看到几个关于这个问题的例子。
最速下降法 的方向向量由下式给出:
( 20 ) d S D ( k ) = a r g m a x d : | | d | | = r { − D f ( x ( k ) ) T d } = a r g m i n d : | | d | | = r { D f ( x ( k ) ) T d } = − r D f ( x ) | | D f ( x ) | | {\displaystyle (20)\quad d_{SD}^{(k)}=argmax_{d:||d||=r}\{-Df(x^{(k)})^{T}d\}=argmin_{d:||d||=r}\{Df(x^{(k)})^{T}d\}=-r{\frac {Df(x)}{||Df(x)||}}}
证明:根据柯西-施瓦茨不等式 ,可以得到
( 21 ) D f ( x ) T d | | D f ( x ) | | | | d | | ≥ − 1 ⇔ D f ( x ) T d ≥ − r | | D f ( x ) | | . {\displaystyle (21)\quad {\frac {Df(x)^{T}d}{||Df(x)||||d||}}\geq -1\quad \Leftrightarrow \quad Df(x)^{T}d\geq -r||Df(x)||.}
显然 (21) 在 d ( k ) = d S D ( k ) {\displaystyle d^{(k)}=d_{SD}^{(k)}} 下成立,如 (20) 所示。
特别要注意的是,对于 r = | | D f ( x ) | | {\displaystyle r=||Df(x)||} ,我们有 d S D ( k ) = − D f ( x ( k ) ) {\displaystyle d_{SD}^{(k)}=-Df(x^{(k)})} ,即我们只是沿着负梯度的方向“行走”。与牛顿法 相比,最速下降法 没有使用固定的步长,而是选择了一个有效的步长 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 。因此,最速下降法 定义了序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 为
( 22 ) x ( k + 1 ) = x ( k ) + σ E ( k ) d S D ( k ) , {\displaystyle (22)\quad x^{(k+1)}=x^{(k)}+\sigma _{E}^{(k)}d_{SD}^{(k)},}
其中 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 是一个有效的步长,而 d S D ( k ) {\displaystyle d_{SD}^{(k)}} 是在 (20) 中给出的最速下降方向。
对于 d S D ( k ) = − r D f ( x ) | | D f ( x ) | | {\displaystyle d_{SD}^{(k)}=-r{\frac {Df(x)}{||Df(x)||}}} ,最速下降法定义了一个下降方向,正如引理 1 所述,因为
D f ( x ) T d S D ( k ) = D f ( x ) T ( − r D f ( x ) | | D f ( x ) | | ) = − r | | D f ( x ) | | D f ( x ) T D f ( x ) < 0. {\displaystyle Df(x)^{T}d_{SD}^{(k)}=Df(x)^{T}\left(-r{\frac {Df(x)}{||Df(x)||}}\right)=-{\frac {r}{||Df(x)||}}Df(x)^{T}Df(x)<0.}
最速下降法 仅在局部上是合理的,因为它忽略了二阶信息。
特别是当目标函数很平坦(即解 x ⋆ {\displaystyle x^{\star }} 位于一个“山谷”中)时,由最速下降法定义的序列会剧烈波动(见下面的应用,尤其是 Rosenbrock 函数的例子)。
由于它不需要海森矩阵,因此最速下降法 的计算和实现很容易且快速。
牛顿法 和最速下降法 都属于可变度量法 类的一般方法。此类方法依赖于更新公式
( 23 ) x k + 1 = x k − σ E ( k ) [ A k ] − 1 D f ( x k ) . {\displaystyle (23)\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}[A^{k}]^{-1}Df(x^{k}).}
如果 A k {\displaystyle A^{k}} 是一个 对称 且 正定 矩阵,(23) 定义了一个下降方法,因为 [ A k ] − 1 {\displaystyle [A^{k}]^{-1}} 正定当且仅当 A k {\displaystyle A^{k}} 也正定。要看到这一点,只需考虑 谱分解
( 24 ) A k = Γ Λ Γ T {\displaystyle (24)\quad A^{k}=\Gamma \Lambda \Gamma ^{T}}
其中 Γ {\displaystyle \Gamma } 和 Λ {\displaystyle \Lambda } 分别是具有 特征向量 和 特征值 的矩阵。如果 A k {\displaystyle A^{k}} 正定,所有特征值 λ i {\displaystyle \lambda _{i}} 严格为正。因此它们的逆 λ i − 1 {\displaystyle \lambda _{i}^{-1}} 也为正,因此 [ A k ] − 1 = Γ Λ − 1 Γ T {\displaystyle [A^{k}]^{-1}=\Gamma \Lambda ^{-1}\Gamma ^{T}} 显然是正定的。但是然后,用 d ( k ) = [ A k ] − 1 D f ( x k ) {\displaystyle d^{(k)}=[A^{k}]^{-1}Df(x^{k})} 代入得到
( 25 ) D f ( x k ) T d ( k ) = − D f ( x k ) T [ A k ] − 1 D f ( x k ) ≡ − v T [ A k ] − 1 v ≤ 0 , {\displaystyle (25)\quad Df(x^{k})^{T}d^{(k)}=-Df(x^{k})^{T}[A^{k}]^{-1}Df(x^{k})\equiv -v^{T}[A^{k}]^{-1}v\leq 0,}
即该方法确实是下降的。到目前为止,我们还没有指定矩阵 A k {\displaystyle A^{k}} ,但很容易看出,对于两种特定的选择,可变度量法 恰好与最速下降法 和牛顿法 分别一致。
对于 A k = I {\displaystyle A^{k}={\mathcal {I}}} (其中 I {\displaystyle {\mathcal {I}}} 是 单位矩阵 ),可以得出
( 22 ′ ) x k + 1 = x k − σ E ( k ) D f ( x k ) {\displaystyle (22')\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}Df(x^{k})}
这仅仅是 *最速下降法*。
对于 A k = D 2 f ( x k ) {\displaystyle A^{k}=D^{2}f(x^{k})} ,可以得出
( 19 ′ ) x k + 1 = x k − σ E ( k ) [ D 2 f ( x k ) ] − 1 D f ( x k ) {\displaystyle (19')\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}[D^{2}f(x^{k})]^{-1}Df(x^{k})}
这仅仅是使用步长 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 的 *牛顿法*。
对于 *变量度量法* 而言,另一个自然的候选方法是 *拟牛顿法*。与标准 *牛顿法* 相比,它使用了一种有效的步长,使其成为一种下降方法,并且与 *最速下降法* 相比,它没有完全忽略有关函数曲率的局部信息。因此,*拟牛顿法* 由矩阵 A k {\displaystyle A^{k}} 上的两个要求定义
A k {\displaystyle A^{k}} 应该近似于海森矩阵 D 2 f ( x k ) {\displaystyle D^{2}f(x^{k})} 以利用有关曲率的信息,以及
更新 A k → A k + 1 {\displaystyle A^{k}\rightarrow A^{k+1}} 应该很容易,这样算法仍然相对快(即使在高维情况下)。
为了确保第一个要求, A k + 1 {\displaystyle A^{k+1}} 应该满足所谓的 *拟牛顿方程*
( 26 ) A k + 1 ( x ( k + 1 ) − x ( k ) ) = D f ( x ( k + 1 ) ) − D f ( x ( k ) ) {\displaystyle (26)\quad A^{k+1}(x^{(k+1)}-x^{(k)})=Df(x^{(k+1)})-Df(x^{(k)})}
因为所有满足 (26) 的 A k {\displaystyle A^{k}} 反映了有关海森矩阵的信息。为了看到这一点,请考虑定义为以下函数的 g ( x ) {\displaystyle g(x)}
( 27 ) g ( x ) = f ( x k + 1 ) + D f ( x k + 1 ) T ( x − x k + 1 ) + 1 2 ( x − x k + 1 ) T A k + 1 ( x − x k + 1 ) . {\displaystyle (27)\quad g(x)=f(x^{k+1})+Df(x^{k+1})^{T}(x-x^{k+1})+{\frac {1}{2}}(x-x^{k+1})^{T}A^{k+1}(x-x^{k+1}).}
那么很明显, g ( x k + 1 ) = f ( x k + 1 ) {\displaystyle g(x^{k+1})=f(x^{k+1})} 并且 D g ( x k + 1 ) = D f ( x k + 1 ) {\displaystyle Dg(x^{k+1})=Df(x^{k+1})} 。所以 g ( x ) {\displaystyle g(x)} 和 f ( x ) {\displaystyle f(x)} 在 x ( k + 1 ) {\displaystyle x^{(k+1)}} 的邻域内是相当相似的。为了确保 g ( x ) {\displaystyle g(x)} 在 x ( k ) {\displaystyle x^{(k)}} 处也是一个好的近似,我们希望选择 A k + 1 {\displaystyle A^{k+1}} 使得在 x ( k ) {\displaystyle x^{(k)}} 处的梯度相同。通过
( 28 ) D g ( x k ) = D f ( x k + 1 ) − A k + 1 ( x k + 1 − x k ) {\displaystyle (28)\quad Dg(x^{k})=Df(x^{k+1})-A^{k+1}(x^{k+1}-x^{k})}
很明显, D g ( x k ) = D f ( x k ) {\displaystyle Dg(x^{k})=Df(x^{k})} 如果 A k + 1 {\displaystyle A^{k+1}} 满足 (26) 中给出的拟牛顿方程 。但随后可以得出
( 29 ) A k + 1 ( x k + 1 − x k ) = D f ( x k + 1 ) − D g ( x k ) = D f ( x k + 1 ) − D f ( x k ) = D 2 f ( λ x ( k ) + ( 1 − λ ) x ( k + 1 ) ) ( x k + 1 − x k ) . {\displaystyle (29)\quad A^{k+1}(x^{k+1}-x^{k})=Df(x^{k+1})-Dg(x^{k})=Df(x^{k+1})-Df(x^{k})=D^{2}f(\lambda x^{(k)}+(1-\lambda )x^{(k+1)})(x^{k+1}-x^{k}).}
因此,只要 x ( k + 1 ) {\displaystyle x^{(k+1)}} 和 x ( k ) {\displaystyle x^{(k)}} 之间相差不大, A k + 1 {\displaystyle A^{k+1}} 满足 (26) 是 D 2 f ( x ( k ) ) {\displaystyle D^{2}f(x^{(k)})} 的一个很好的近似。
现在让我们来谈谈第二个要求,即 A k {\displaystyle A^{k}} 的更新应该很容易。一种特定的算法是所谓的 BFGS 算法 。该算法的主要优点是它只使用已经计算出的元素 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 和 { D f ( x ( k ) ) } k {\displaystyle \{Df(x^{(k)})\}_{k}} 来构建更新 A ( k + 1 ) {\displaystyle A^{(k+1)}} 。因此,不需要计算新的实体,只需要跟踪 x 序列和梯度序列即可。作为 BFGS 算法的起点,可以提供任何正定矩阵(例如,单位矩阵或在 x ( 0 ) {\displaystyle x^{(0)}} 处的 Hessian)。然后,BFGS 更新公式 由下式给出
( 30 ) A k = A k − 1 − ( A k − 1 ) T γ k − 1 T γ k − 1 A k − 1 γ k − 1 T A k − 1 γ k − 1 + Δ k − 1 Δ k − 1 T Δ k − 1 T γ k − 1 {\displaystyle (30)\quad A^{k}=A^{k-1}-{\frac {(A^{k-1})^{T}\gamma _{k-1}^{T}\gamma _{k-1}A^{k-1}}{\gamma _{k-1}^{T}A^{k-1}\gamma _{k-1}}}+{\frac {\Delta _{k-1}\Delta _{k-1}^{T}}{\Delta _{k-1}^{T}\gamma _{k-1}}}}
其中 Δ k − 1 = D f ( x ( k ) ) − D f ( x ( k − 1 ) ) {\displaystyle \Delta _{k-1}=Df(x^{(k)})-Df(x^{(k-1)})} 以及 γ k − 1 = x ( k ) − x ( k − 1 ) {\displaystyle \gamma _{k-1}=x^{(k)}-x^{(k-1)}} 。此外,(30) 确保所有 A k {\displaystyle A^{k}} 都是正定的,正如变尺度法 所要求的那样,是下降的。
它使用关于 f ( x ) {\displaystyle f(x)} 曲率的二阶信息,矩阵 A k {\displaystyle A^{k}} 与海森矩阵 D 2 f ( x ) {\displaystyle D^{2}f(x)} 有关。
尽管如此,它确保了简单快速的更新(例如,通过 BFGS 算法),因此它比标准牛顿法更快。
它是一种下降方法,因为 A k {\displaystyle A^{k}} 是正定的。
它相对容易实现,因为 BFGS 算法在大多数数值或统计软件包中都可用。
为了比较这些方法并说明算法之间的差异,我们现在将评估最速下降法 、标准牛顿法 和拟牛顿法 在有效步长下的性能。我们使用该领域中的两个经典函数,即 Himmelblau 函数和 Rosenbrock 函数。
Himmelblau 函数由下式给出
( 31 ) f ( x , y ) = ( x 2 + y − 11 ) 2 + ( x + y 2 − 7 ) 2 {\displaystyle (31)\quad f(x,y)=(x^{2}+y-11)^{2}+(x+y^{2}-7)^{2}}
这个四阶多项式有四个极小值,四个鞍点和一个极大值,因此算法有足够的失败可能性。在下图中,我们显示了 等高线图 和函数的 3D 图,用于不同的起始值。
在图 1 中,我们显示了函数以及所有三种方法在起始值为 ( 2 , − 4 ) {\displaystyle (2,-4)} 时的路径。显然,三种方法没有找到相同的极小值。原因当然是由最速下降法 的不同方向向量造成的——通过忽略关于曲率的信息,它选择的方向与两个牛顿法 完全不同(尤其是在图 1 的右面板中)。
图 1:两个牛顿法收敛到相同的极小值,而最速下降法收敛到不同的极小值。
现在考虑起始值为 ( 4.5 , − 0.5 ) {\displaystyle (4.5,-0.5)} ,如图 2 所示。最重要的是,现在所有 方法都找到了不同的解。最速下降法 找到的解与两个牛顿法 不同,这并不令人惊讶。但是,两个牛顿法 收敛到不同的解,这表明了步长 σ {\displaystyle \sigma } 的重要性。使用拟牛顿法 在第一次迭代中选择一个有效的步长,两种方法在第一次迭代后的所有迭代中都有不同的步长和 方向向量。从图中可以看到:结果可能非常重要。
图 2:即使所有方法都找到了不同的解。
Rosenbrock 函数由下式给出
( 32 ) f ( x , y ) = 100 ( y − x 2 ) 2 + ( 1 − x ) 2 {\displaystyle (32)\quad f(x,y)=100(y-x^{2})^{2}+(1-x)^{2}\ }
虽然这个函数只有一个极小值,但它对于优化问题来说是一个有趣的函数。原因是这个 U 形函数的谷非常平坦(参见图 3 和 4 的右面板)。对于 计量经济学家 来说,这个函数可能很有趣,因为在最大似然估计的情况下,平坦的准则函数经常出现。因此,在下面图 3 和 4 中显示的结果对于具有此问题的函数来说似乎是相当通用的。
我使用这个函数和所采用的算法的经验是,图 3(给定一个初始值为 ( 2 , − 5 ) {\displaystyle (2,-5)} )似乎很有代表性。与上面的 Himmelblau 函数相比,所有算法都找到了相同的解,鉴于只有一个最小值,这是可以预期的。更重要的是,不同方法选择的路徑,反映了各自方法的不同属性。可以看出,**最速下降法**波动相当剧烈。这是由于它不使用曲率信息,而是来回跳跃“山丘”之间的“山谷”。两种牛顿法选择更直接的路徑,因为它们使用二阶信息。两种牛顿法的主要区别当然在于步长。图 3 显示,**拟牛顿法**在穿过山谷时使用非常小的步长。相反,**牛顿法**的步长是固定的,因此它直接跳向解的方向。虽然人们可能会认为这是**拟牛顿法**的缺点,但需要注意的是,通常这些较小的步长伴随着更高的稳定性,即算法不太可能跳跃到不同的解。这可以在图 4 中看到。
图 3:所有方法都找到了相同的解,但最速下降法波动很大。
图 4 考虑了初始值为 ( − 2 , − 2 ) {\displaystyle (-2,-2)} ,显示了使用固定步长的**牛顿法**的主要问题 - 该方法可能会“越界”,因为它没有下降。在第一步中,**牛顿法**(在图中显示为紫色线)跳出了山谷,只是在下一步迭代中反弹回来。在这种情况下,仍然会收敛到最小值,因为每侧的梯度都指向中心的单个山谷,但可以很容易地想象这种情况不成立的函数。这种跳跃的原因是二阶导数非常小,因此步长 [ D f ( x ( k ) ) ] − 1 D f ( x ( k ) ) ) {\displaystyle [Df(x^{(k)})]^{-1}Df(x^{(k)}))} 由于 Hessian 的逆而变得非常大。根据我的经验,我建议使用有效的步长来更好地控制各自方法选择的路徑。
图 2:由于固定步长,牛顿法越界。
对于计量经济学家和统计学家,**最大似然估计器**可能是数值优化算法最重要的应用。因此,我们将简要介绍估计程序如何适应上述框架。
像往常一样,令
( 33 ) f ( Y | X ; θ ) {\displaystyle (33)\quad f(Y|X;\theta )}
是给定 X 时的 Y 的**条件密度**,参数为 θ {\displaystyle \theta } ,以及
( 34 ) l ( θ ; Y | X ) {\displaystyle (34)\quad l(\theta ;Y|X)}
是参数 θ {\displaystyle \theta } 的**条件似然函数**。
如果我们假设数据是**独立同分布 (iid)**,那么样本对数似然函数为
( 35 ) L ( θ ; Y 1 , … , Y N ) = ∑ i N L ( θ ; Y i ) = ∑ i N log ( l ( θ ; Y i ) ) . {\displaystyle (35)\quad {\mathcal {L}}(\theta ;Y_{1},\dots ,Y_{N})=\sum _{i}^{N}{\mathcal {L}}(\theta ;Y_{i})=\sum _{i}^{N}\log(l(\theta ;Y_{i})).}
因此,最大似然估计归结为对参数 θ {\displaystyle \theta } 最大化 (35)。为了简单起见,我们只决定使用**牛顿法**来解决这个问题,序列 { θ ( k ) } k {\displaystyle \{\theta ^{(k)}\}_{k}} 由递归定义为
( 36 ) D θ L ( θ ( k + 1 ) ) = D θ L ( θ ( k ) ) + D θ θ L ( θ ( k ) ) ( θ ( k + 1 ) − θ ( k ) ) = 0 ⇔ θ ( k + 1 ) = θ ( k ) − [ D θ θ L ( θ ( k ) ) ] − 1 D θ L ( θ ( k ) ) {\displaystyle (36)\quad D_{\theta }{\mathcal {L}}(\theta ^{(k+1)})=D_{\theta }{\mathcal {L}}(\theta ^{(k)})+D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})(\theta ^{(k+1)}-\theta ^{(k)})=0\Leftrightarrow \theta ^{(k+1)}=\theta ^{(k)}-[D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})]^{-1}D_{\theta }{\mathcal {L}}(\theta ^{(k)})}
其中 D θ L {\displaystyle D_{\theta }{\mathcal {L}}} 和 D θ θ L {\displaystyle D_{\theta \theta }{\mathcal {L}}} 分别表示关于参数向量 θ {\displaystyle \theta } 的一阶导数和二阶导数,而 [ D θ θ L ( θ ( k ) ) ] − 1 D θ L ( θ ( k ) ) {\displaystyle [D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})]^{-1}D_{\theta }{\mathcal {L}}(\theta ^{(k)})} 定义了 (17) 中给出的牛顿方向。由于最大似然估计总是假设条件密度(即误差项的分布)在参数 θ {\displaystyle \theta } 范围内是已知的,因此上述方法可以很容易地应用。
假设一个简单的线性模型
( 37 a ) Y i = β 1 + β x X i + U i {\displaystyle (37a)\quad Y_{i}=\beta _{1}+\beta _{x}X_{i}+U_{i}}
其中 θ = ( β 1 , β 2 ) ′ {\displaystyle \theta =(\beta _{1},\beta _{2})'} 。那么,Y 的条件分布由U 的分布决定,即
( 37 b ) p ( Y i − β 1 − β x X i ) ≡ p ∣ X i ( Y i ) = p ( U i ) , {\displaystyle (37b)\quad p(Y_{i}-\beta _{1}-\beta _{x}X_{i})\equiv p_{\mid X_{i}}(Y_{i})=p(U_{i}),}
其中 p {\displaystyle p} 表示 密度函数 。一般来说,最大化 (35) 没有闭合形式解(至少如果U 恰好不是 正态分布 ),因此必须使用数值方法。因此,假设U 服从 学生 t 分布 ,自由度为m 自由度 ,因此 (35) 由下式给出
( 38 ) L ( θ ; Y ∣ X ) = ∑ log ( Γ ( m + 1 2 ) π m Γ ( m 2 ) ( 1 + ( y i − x i T β ) 2 m ) − m + 1 2 ) {\displaystyle (38)\quad {\mathcal {L}}(\theta ;Y_{\mid X})=\sum \log \left({\frac {\Gamma ({\frac {m+1}{2}})}{{\sqrt {\pi m}}\Gamma ({\frac {m}{2}})}}(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}})^{-{\frac {m+1}{2}}}\right)}
我们只是使用了t分布密度函数的定义。(38)可以简化为
( 39 ) L ( θ ; Y ∣ X ) = N [ log ( Γ ( m + 1 2 ) ) − log ( π m Γ ( m 2 ) ) ] − m + 1 2 ∑ log ( 1 + ( y i − x i T β ) 2 m ) {\displaystyle (39)\quad {\mathcal {L}}(\theta ;Y_{\mid X})=N\left[\log \left(\Gamma \left({\frac {m+1}{2}}\right)\right)-\log \left({\sqrt {\pi m}}\Gamma \left({\frac {m}{2}}\right)\right)\right]-{\frac {m+1}{2}}\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)}
所以(如果我们假设自由度m 是已知的)
( 40 ) argmax { L ( θ ; Y ∣ X ) } = argmax { − m + 1 2 ∑ log ( 1 + ( y i − x i T β ) 2 m ) } = argmin { ∑ log ( 1 + ( y i − x i T β ) 2 m ) } . {\displaystyle (40)\quad {\textit {argmax}}\left\{{\mathcal {L}}(\theta ;Y_{\mid X})\right\}={\textit {argmax}}\left\{-{\frac {m+1}{2}}\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)\right\}={\textit {argmin}}\left\{\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)\right\}.}
使用准则函数
( 41 ) f ( β 1 , β 2 ) = ∑ log ( 1 + ( y i − β 1 − β 2 x i ) 2 m ) {\displaystyle (41)\quad f(\beta _{1},\beta _{2})=\sum \log \left(1+{\frac {(y_{i}-\beta _{1}-\beta _{2}x_{i})^{2}}{m}}\right)}
上述方法可以很容易地应用于计算最大似然估计量 ( β ^ 1 , M L , β ^ 2 , M L ) {\displaystyle ({\hat {\beta }}_{1,ML},{\hat {\beta }}_{2,ML})} ,最大化(41)。
Alt, W. (2002): "Nichtlineare Optimierung", Vieweg: Braunschweig/Wiesbaden
Härdle, W. and Simar, L. (2003): "Applied Multivariate Statistical Analysis", Springer: Berlin Heidelberg
Königsberger, K. (2004): "分析 I", 施普林格出版社: 柏林 海德堡
Ruud, P. (2000): "经典计量经济学理论", 牛津大学出版社: 纽约