跳转至

下降方法

约 1032 个字 预计阅读时间 3 分钟

梯度方法

所以就有了下面的方案:

\[\boxed{ \begin{array}{l}\Large\textbf{梯度方法}\\ \rule[4pt]{9.5cm}{0.05em}\\ \text{选择 }x_0 \in\mathbb{R}^n\\ \text{迭代 } x_{k+1}=x_{k}-h_k\nabla f(x_k),\; k=0,1,\dots \end{array}}\]

我们尝试用一种“普遍”的推理推出这种方案:取定点 \(x_0\),考虑函数 \(f\in C^1(\mathbb{R}^n)\) 与其下面的近似:

\[ \phi_2(x) = f(x_0) + \langle\nabla f(x_0),x-x_0\rangle + \dfrac{1}{2h}\lVert x-x_0\rVert^2 \]

其中参数 \(h>0\) 是一个正数,根据一阶最优性条件,该函数的无约束极小点 \(x_1^*\) 满足:

\[ \nabla\phi_2(x_1^*) = \nabla f(x_0)+\dfrac{1}{h}(x_1^*-x_0) = 0 \]

所以就有 \(x_1^* = x_0-h\nabla f(x_0)\),这就是我们梯度法的迭代方案。同时注意到如果 \(h\in \left( 0,\frac{1}{L} \right]\),那么函数 \(\phi_2\) 就是 \(f\) 的全局上近似,即:

\[ f(x) \leqslant \phi_2(x),\; \forall x\in\mathbb{R}^n \]

这就保证了梯度法的全局收敛性。记好这个 \(\phi_2\),我们在拟牛顿法的比较中还会用到(虽然也就是提了一嘴)。

学习率其实也是一种超参数,这是一种和神经网络的参数性质不同的参数,因为它不是通过训练得到的,而是需要人为设定的。如果学习率过小,有可能就会令学习花费更多时间,或者是没怎么更新就结束了;但是如果学习率过大,就有可能会导致学习发散而不能正确进行。

学习率的序列 \(\left\{h_k \right\}_{k = 1}^{\infty}\) 有这样几种选择:

  • 预先选择序列:
    • \(h_k = h > 0\) (固定步长)
    • \(h_k = \dfrac{h}{\sqrt{k+1}}\) (带有学习率衰减)
  • \(h_k = \arg\min\limits_{h\ge 0}f(x_k-h\nabla f(x_k))\) (精确线搜索)
  • Armijo-Goldstein 准则。

Armijo-Goldstein 准则的含义为:寻找一个 \(x_{k+1}=x_k-h\nabla f(x_k)\) ,使得:

\[ \begin{gather} \alpha\langle f(x_k),x_k-x_{k+1}\rangle \leqslant f(x_k)-f(x_{k+1})\\ \beta\langle f(x_k),x_k-x_{k+1}\rangle \geqslant f(x_k)-f(x_{k+1}) \end{gather} \]

其中 \(\alpha,\beta\in(0,1)\) 为固定的参数。我们可以对 Armijo-Goldstein 准则进行几何解释:固定 \(x\in\mathbb{R}^n\) ,考虑单变量函数 \(\phi(h) = f(x-h\nabla f(x)),\; h>0\) ,Armijo-Goldstein 准则可接受的步长值对应于 \(\phi\) 的图像的特定部分,该部分位于两个线性函数之间:

\[ \begin{gather} \phi_1(h)=f(x)-\alpha h\lVert\nabla f(x)\rVert^2\\ \phi_2(h)=f(x)-\beta h\lVert\nabla f(x)\rVert^2 \end{gather} \]

注意到 \(\phi(0)=\phi_1(0)=\phi_2(0)\) ,以及 \(\phi'(0)<\phi_2'(0)<\phi_1'(0)<0\) ,因此除非 \(\phi\) 没有下界,否则可接受的步长总是存在的,也就是 \(\phi\)\(\phi_1\)\(\phi_2\) 都一定会有第二个交点。有许多快速确定满足该条件的点的一维搜索算法。

梯度方法的性能估计

总的来说,我们考虑的是对不同类型的函数,梯度方法对下面的优化问题的性能:

\[ \underset{\boldsymbol{x}\in\mathbb{R}^n}{\operatorname{argmin}}f(\boldsymbol{x}) \]

连续函数类

凸函数类

强凸函数类

一般下降方法

我们已经介绍了梯度下降,也就是下降的方向是梯度的反方向,同样也证明了梯度下降确实下降最快(真的对吗?)。但是我们也希望考察是否有一种一般的下降方向:

我们讲的下降方向/Descent Direction 有如下定义:\(\boldsymbol{d}\)\(f\)\(x\) 处的一个下降方向,如果对于任意足够小的 \(t>0\)\(f(x+t\boldsymbol{d})<f(x)\)。并且也很轻松得到下面的命题:如果 \(f\)\(x\) 的一个邻域内连续可微,那么任意使得 \(\boldsymbol{d}^T\nabla f(x)<0\)\(\boldsymbol{d}\) 都是下降方向。

\(f\)\(x\) 处沿着方向 \(\boldsymbol{v}\in\mathbb{R}^n\) 的变化率可以用方向导数衡量:

\[ \nabla_{\boldsymbol{v}}f(x) = \lim_{\varepsilon\to 0}\dfrac{f(x+\varepsilon\boldsymbol{v})-f(x)}{\varepsilon} = \langle \boldsymbol{v},\nabla f(x)\rangle \]

\(f\)\(x\) 处关于范数 \(\lVert\cdot\rVert\) 的最速的下降方向为:

\[ \Delta_{\lVert\cdot\rVert}x = \operatorname*{argmin}_{\boldsymbol{v}}\langle \boldsymbol{v},\nabla f(x)\rangle \]

这里我们需要通过一个范数 \(\lVert\cdot\rVert\) 来约束 \(\boldsymbol{v}\) 的长度,这个范数很重要,因为选定不同的范数,可以得到不同的最速下降方向。假定 \(\nabla f(x)\neq 0\),有:

\[ \begin{align} \Delta_{\lVert\cdot\rVert_2}x &= -\lVert\nabla f(x)\rVert_2^{-1}\nabla f(x)\\ \Delta_{\lVert\cdot\rVert_1}x &= -\operatorname*{sign}\left(\dfrac{\partial f(x)}{\partial x_l}\right)\cdot\boldsymbol{e}_l,\; l = \operatorname*{argmax}_{i\in1,\dots,n}\left\lvert\dfrac{\partial f(x)}{\partial x_i}\right\rvert\\ \Delta_{\lVert\cdot\rVert_{\boldsymbol{A}}}x &= -\lVert\nabla f(x)\rVert_{\boldsymbol{A}^{-1}}^{-1}\boldsymbol{A}^{-1}\nabla f(x) \end{align} \]

其中 \(\boldsymbol{A}\) 是正定对称方阵;\(\boldsymbol{A}\) 范数 \(\lVert\cdot\rVert_{\boldsymbol{A}}\) 定义为 \(\lVert\boldsymbol{v}\rVert_{\boldsymbol{A}} = \sqrt{\boldsymbol{v}^T \boldsymbol{\boldsymbol{A}} \boldsymbol{v}}\)。这就发现了一个神奇的事实,在不同范数下,我们确实得到了不同的最速下降方向,回头看我们在上一节光滑优化基础提到的梯度下降是局部最速下降的陈述,现在我们就可以看到这个陈述更严谨的叙述了:

  • \(f\)\(x\) 处的 \(2\) 范数下的最速下降方向为 \(d_{gd} = -\nabla f(x)\),这就是我们所说的梯度下降;
  • \(f\)\(x\) 处的 \(1\) 范数下的最速下降方向为 \(d_{cd}=-\dfrac{\partial f(x)}{\partial x_l}\boldsymbol{e}_l,\; l = \operatorname*{argmax}_{i\in1,\dots,n}\left\lvert\dfrac{\partial f(x)}{\partial x_i}\right\rvert\),这叫做坐标下降/Coordinate Descent,即每次下降都依次沿坐标轴的方向最小化目标函数值;
  • \(f\)\(x\) 处的 \(\boldsymbol{A}\) 范数下的最速下降方向为 \(d_{\boldsymbol{A}}=-\boldsymbol{A}^{-1}\nabla f(x)\);如果我们选择 \(\boldsymbol{A}\) 为 Hessian 矩阵,那么得到的方法便是马上就要介绍的牛顿法。

此外,还有随机选择下降方向的策略,之需要满足下降方向的期望是负梯度方向即可。以坐标下降为例,\(f\)\(x\) 处的下降方向是随机选择的,由下式给出:

\[ d_{cd-rand} = -\nabla f(x)\boldsymbol{e}_{i_k} \]

这里的 \(i_k\) 每次从 \(\{1,2,\dots,n\}\) 中均匀随机选择;也就是说,在所有是下降方向的坐标轴方向之中,均匀随机选择一个方向。可以证明这种方法的期望正是负梯度方向。

总而言之,对于一种随机梯度方法,\(f\)\(x\) 处的梯度由下式给出:

\[ d_{sgc} = -g(x_k,\xi_k) \]

其中 \(\xi_k\) 是随机变量,使得 \(\mathbb{E}_{\xi_k}g(x_k,\xi_k)=\nabla f(x_k)\),即 \(g(x_k,\xi_k)\) 是对负梯度方向的无偏估计/Unbiased Estimation。

牛顿法

牛顿法基于线性逼近,其最出名的应用就是寻找但变量函数的零点。令 \(\phi(t):\mathbb{R}\mapsto\mathbb{R}\) ,考虑问题 \(\phi(t^*)=0\),我们找到了某个足够逼近 \(t^*\)\(t\),注意到:

\[ \begin{align} \phi(t+\Delta t) &= \phi(t)+\phi'(t)\Delta t+o(\lvert\Delta t\rvert) \end{align} \]

并且我们预估 \(t+\Delta t\)\(t^*\),那么由 \(\phi(t+\Delta t)=0\) 可以得到 \(\phi(t)+\phi'(t)\Delta t=0\),即 \(\Delta t = -\dfrac{\phi(t)}{\phi'(t)}\),我们可以认为 \(\Delta t\) 是对最优位移 \(t^*-t\) 的一个足够好的近似,也就是说在 \(t\) 处的下降步长应该是 \(-\dfrac{\phi(t)}{\phi'(t)}\)。将这个想法转化为算法的形式,就可以得到:

\[ t_{k+1} = t_k-\dfrac{\phi(t_k)}{\phi'(t_k)} \]

可以证明,如果函数是连续的,并且待求的零点是孤立的,那么在零点周围存在一个区域,只要初始值位于这个邻近区域内,那么牛顿法必定收敛,这保证了牛顿法的有效性。我们也可以将牛顿法推广到求解非线性方程组的问题,即 \(F(\boldsymbol{x})=\boldsymbol{0}\),其中 \(\boldsymbol{x}\in\mathbb{R}^n\)\(F:\mathbb{R}^n\mapsto\mathbb{R}^n\)。这种情况下,我们需要定义一个增量 \(\Delta\boldsymbol{x}\),作为如下线性方程组的解:

\[ F(\boldsymbol{x})+J_F(\boldsymbol{x})\Delta\boldsymbol{x} = \boldsymbol{0} \]

其中 \(J_F(\boldsymbol{x})\) 是雅各比矩阵/Jacobian Matrix,如果其是非退化的/Nondegenerate,即矩阵是满秩的,那么我们可以计算出位移 \(\Delta\boldsymbol{x} = -J_F(\boldsymbol{x})^{-1}F(\boldsymbol{x})\),迭代过程为:

\[ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k-J_F(\boldsymbol{x}_k)^{-1}F(\boldsymbol{x}_k) \]

考虑一个连续可微函数 \(f\) 的一阶条件,我们可以将无约束最小化问题转换成求解线性系统/Linear System \(\nabla f(\boldsymbol{x})=\boldsymbol{0}\) 的解,注意到梯度的雅各比矩阵正是 Hessian 矩阵,那么此时的牛顿系统为 \(\nabla f(\boldsymbol{x})+\nabla^2 f(\boldsymbol{x})\Delta\boldsymbol{x}=\boldsymbol{0}\),因此,应用牛顿法求解该优化问题,迭代形式应为:

\[ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k-\left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\nabla f(\boldsymbol{x}_k) \]

这样我们就得到了牛顿法的一般形式,下面我们将对牛顿法进行详细的讨论。

基本形式

换种方式重新推导:对于二阶连续可微函数 \(f\),考虑其在 \(\boldsymbol{x}_k\) 点处的二阶近似 \(\phi_1\)

\[ \phi_1(\boldsymbol{x}) = f(\boldsymbol{x}_k)+\langle\nabla f(\boldsymbol{x}_k),\boldsymbol{x}-\boldsymbol{x}_k\rangle+\dfrac{1}{2}\left\langle\nabla^2 f(\boldsymbol{x}_k)(\boldsymbol{x}-\boldsymbol{x}_k),\boldsymbol{x}-\boldsymbol{x}_k\right\rangle \]

\(\nabla^2 f(\boldsymbol{x})\succ 0\),也就是 \(f\) 为凸函数,我们可以选择二次函数 \(\phi_1\) 的极小点作为 \(\boldsymbol{x}_{k+1}\),即:

\[ \nabla\phi_1(\boldsymbol{x}_{k+1}) = \nabla f(\boldsymbol{x}_k)+\nabla^2 f(\boldsymbol{x})(\boldsymbol{x}_{k+1}-\boldsymbol{x}_k) = \boldsymbol{0} \]

此即牛顿法的迭代形式。在严格局部极小点的一个邻域内,牛顿法的收敛速度是非常快的,几次迭代就到了。但是牛顿法有一些严重的缺点:

  • 要求目标函数存在 Hessian 矩阵,且 Hessian 矩阵必须可逆;
  • 牛顿法有可能会发散,理论上只有在 \(\boldsymbol{x}^*\) 很近的位置才会收敛;
  • 虽然收敛更快,但牛顿法的计算较复杂,除需计算梯度而外,还需计算 Hessian 矩阵和其逆矩阵,对计算量、存储量都有要求,当维数很大时这个问题更加突出,尤其是计算逆矩阵的开销达到了 \(O(n^3)\)

下面有一个牛顿法发散的例子:考虑单变量函数 \(\phi(t) = \dfrac{t}{\sqrt{1+t^2}}\) 的零点 \(t^*\),显然 \(t^*=0\);由 \(\phi'(t) = \dfrac{1}{(1+t^2)^{\frac{3}{2}}}\),则牛顿法的迭代步骤为:

\[ t_{k+1} = t_k-\dfrac{\phi(t)}{\phi'(t)} = t_k-\dfrac{t_k}{\sqrt{1+t_k^2}}(1+t_k^2)^{\frac{3}{2}} = -t_k^3 \]

因此,如果 \(\lvert t_0\rvert<1\),那么牛顿法会迅速收敛;如果 \(\lvert t_0\rvert=1\),那么牛顿法会在 \(\pm 1\) 处反复震荡;如果 \(\lvert t_0\rvert>1\),那么牛顿法震荡且迅速发散。在实践中,为了防止发散,会采用阻尼牛顿法/Damped Newton Method,即为牛顿法添加一个步长 \(h_k > 0\)

\[ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k-h_k\left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\nabla f(\boldsymbol{x}_k) \]

在阻尼牛顿法的初始阶段可以使用与梯度方法相同的步长;后续阶段应当令 \(h_k=1\),即进入满足牛顿法的邻域之后,使用纯牛顿法快速下降。

收敛分析

简单说来,当起始点与局部极小点足够接近的时候,牛顿法的收敛速率是平方收敛/Qudratic Convergence,可以总结成以下定理:

定理:设 \(f(\boldsymbol{x})\) 满足 \(f\in C_M^{2,2}(\mathbb{R}^n)\),且 \(\nabla^2 f(\boldsymbol{x})\succeq l\boldsymbol{I}_n\),且初始点 \(\boldsymbol{x}_0\)\(\boldsymbol{x}^*\) 足够接近,即 \(\lVert\boldsymbol{x}_0-\boldsymbol{x}^*\rVert \leqslant \bar{r} = \dfrac{2l}{3M}\),那么对于 \(\lVert\boldsymbol{x}_k-\boldsymbol{x}^*\rVert<\bar{r}\),牛顿法以平方速率收敛:

\[ \lVert\boldsymbol{x}_{k+1}-\boldsymbol{x}^*\rVert \leqslant \dfrac{M\lVert\boldsymbol{x}_k-\boldsymbol{x}^*\rVert^2}{2(l-M\lVert\boldsymbol{x}_k-\boldsymbol{x}^*\rVert)} \]

忘掉这个定理,我们开始推导一下平方收敛是如何得出的,这里我们有两个任务:首先是得到 \(\lVert\boldsymbol{x}_{0} - \boldsymbol{x}^*\rVert\) 的上界,这规定了我们在什么情况下才适合使用牛顿法;其次是得到 \(\lVert\boldsymbol{x}_{k+1} - \boldsymbol{x}^*\rVert\)\(\lVert\boldsymbol{x}_{k} - \boldsymbol{x}^*\rVert\) 之间的关系,这表示了我们的收敛速率。

推导

我们考虑牛顿法的迭代过程:

\[ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k-\left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\nabla f(\boldsymbol{x}_k) \]

有:

\[ \begin{align} \boldsymbol{x}_{k+1}-\boldsymbol{x}^* &= \boldsymbol{x}_k-\boldsymbol{x}^*-\left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1} (\nabla f(\boldsymbol{x}_k) - 0)\\ &= \boldsymbol{x}_k-\boldsymbol{x}^*-\left(\nabla^2 f(\boldsymbol{x})\right)^{-1}\int_0^1\nabla^2 f(\boldsymbol{x}^*+\tau(\boldsymbol{x}_k-\boldsymbol{x}^*))(\boldsymbol{x}_k-\boldsymbol{x}^*)\,\mathrm{d}\tau\\ &= \left(\nabla^2 f(\boldsymbol{x_k})\right)^{-1} \cdot (\boldsymbol{x}_k - \boldsymbol{x}^*) \cdot \int_0^1 \left(\nabla^2 f(\boldsymbol{x}_k) - \nabla^2 f(\boldsymbol{x}^*+\tau(\boldsymbol{x}_k-\boldsymbol{x}^*))\right)\,\mathrm{d}\tau\\ &= \left(\nabla^2 f(\boldsymbol{x_k})\right)^{-1} G_k \cdot (\boldsymbol{x}_k - \boldsymbol{x}^*) \end{align} \]

其中我们希望 \(\boldsymbol{x}_k - \boldsymbol{x}^*\)\(\boldsymbol{x}_{k+1}-\boldsymbol{x}^*\) 的变换是一个压缩映射,因此需要考察作用在其上的算子 \(\left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}G_k\);这说明 \(\left(\nabla^2 f(\boldsymbol{x}_k\right)^{-1}\)\(G_k\) 都需要考察。

回忆我们在光滑优化基础在可微函数类一节的最后一个推论,有:

\[ \nabla^2 f(\boldsymbol{x}) - M\lVert y - x\rVert \boldsymbol{I}_n \preceq \nabla^2 f(\boldsymbol{y}) \preceq \nabla^2 f(\boldsymbol{x}) + M\lVert y - x\rVert \boldsymbol{I}_n \]

利用这个推论,记 \(r_k = \lVert\boldsymbol{x}_k-\boldsymbol{x}^*\rVert\),我们有:

\[ \begin{align} \lVert G_k\rVert &= \left\lVert\int_0^1\left(\nabla^2 f(\boldsymbol{x}_k)-\nabla^2 f(\boldsymbol{x}^*+\tau(\boldsymbol{x}_k-\boldsymbol{x}^*))\right)\,\mathrm{d}\tau\right\rVert\\ &\leqslant \int_0^1\left\lVert \nabla^2 f(\boldsymbol{x}_k)-\nabla^2 f(\boldsymbol{x}^*+\tau(\boldsymbol{x}_k - \boldsymbol{x}^*))\right\rVert\,\mathrm{d}\tau\\ &\leqslant \int_0^1 M(1-\tau)r_k\,\mathrm{d}\tau = \dfrac{r_k}{2}M \end{align} \]

另一方面,还是根据这个推论,有:

\[ \nabla^2 f(\boldsymbol{x}_k)\succeq \nabla^2 f(\boldsymbol{x}^*)-Mr_k\boldsymbol{I}_n\succeq (l-Mr_k)\boldsymbol{I}_n \]

因此,如果 \(r_k<\dfrac{l}{M}\),那么 \(\nabla^2 f(\boldsymbol{x}_k)\) 是正定的,有:

\[ \left\lVert \left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\right\rVert \leqslant (l-Mr_k)^{-1} \]

所以得到了 \(\left\lVert \left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\right\rVert \leqslant (l-Mr_k)^{-1}\),总结我们得到的条件:

  • \(r_{k+1} = \left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}G_kr_k\)
  • \(\lVert G_k\rVert \leqslant \dfrac{r_k}{2}M\)
  • \(\left\lVert \left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\right\rVert \leqslant (l-Mr_k)^{-1}\)

所以

\[ \begin{align} r_{k+1} = \lVert x_{k+1}-x^*\rVert &= \left\lVert \left(\nabla^2 f(\boldsymbol{x}_k)\right)^{-1}\right\rVert \cdot \lVert G_k\rVert \cdot \lVert x_k - x^*\rVert\\ &\leqslant \left(l - Mr_k\right)^{-1} \cdot \dfrac{r_k}{2}M \cdot r_k\\ &= \dfrac{M}{2(l-Mr_k)}r_k^2 \end{align} \]

我们还希望这个结论再紧凑一些,也就是说,我们希望估计 \(r_{k+1}\) 得到的上界是小于 \(r_k\) 的,所以就得到了

\[ r_{k+1} \leqslant \dfrac{M}{2(l-Mr_k)}r_k^2 < r_k \Rightarrow r_k < \dfrac{2l}{3M} \]

总结上述推导,就得到了牛顿法的平方收敛性质。

拟牛顿法

我们已经介绍过,牛顿法最吸引人的是迭代次数非常少,下降非常快,但是每次迭代的开销非常大,其中每次计算逆矩阵的复杂度就要达到大约 \(O(n^3)\)。我们希望能够找到一个 \(f(\boldsymbol{x})\) 的近似,其比梯度下降法使用到的逼近 \(\phi_2\) 要好,而又比牛顿法的二次逼近 \(\phi_1\) 开销更低。

思路很直接,既然求 Hessian 矩阵的复杂度是大头,那么我们可以用其他方式逼近 Hessian 矩阵的逆,如果我们的逼近复杂度不高且逼近效果很好,那么我们就可以“替代”牛顿法了。令 \(\boldsymbol{G}\) 是一个对称正定 \(n\times n\) 矩阵,记

\[ \phi_{\boldsymbol{G}}(\boldsymbol{x}) = f(\boldsymbol{x}_k)+\langle \nabla f(\boldsymbol{x}_k),\boldsymbol{x}-\boldsymbol{x}_k \rangle+\dfrac{1}{2}\langle \boldsymbol{G}(\boldsymbol{x}-\boldsymbol{x}_k),\boldsymbol{x}-\boldsymbol{x}_k \rangle \]

通过计算该式的极小值点,有更新公式:\(\boldsymbol{x}_{k+1} = \boldsymbol{x}_{\boldsymbol{G}}^* = \boldsymbol{x}_k-\boldsymbol{G}^{-1}\nabla f(\boldsymbol{x}_k)\)。如果我们能够有一系列逼近 \(\nabla^2 f(\boldsymbol{x})\) 的矩阵:\(\{\boldsymbol{G}_k\}:\;\boldsymbol{G}_k\to \nabla^2 f(\boldsymbol{x})\);或者一系列逼近 \(\left(\nabla^2 f(\boldsymbol{x})\right)^{-1}\) 的矩阵:\(\{\boldsymbol{H}_k\}:\;\boldsymbol{H}_k\equiv\boldsymbol{G}_k^{-1}\to \left(\nabla^2 f(\boldsymbol{x})\right)^{-1}\),那么就可以以更小的代价替代牛顿法了。

上述的方法叫做拟牛顿法/Quasi-Newton Method,也称为变尺度法/Variable Metric Method。拟牛顿法的关键是获取矩阵 \(\boldsymbol{G}_k\) 或者 \(\boldsymbol{H}_k\) 序列,不管这个问题,我们先介绍拟牛顿法的一般形式:

\[ \boxed{\begin{array}{l} \Large \textbf{拟牛顿法}\\ \rule[4pt]{12.5cm}{0.05em}\\ 0.\; 选择\;\boldsymbol{x}_0\in\mathbb{R}^n,\; 令\;\boldsymbol{H}_0 = \boldsymbol{I}_n,\; 计算\;f(\boldsymbol{x}_0)\;和\;\nabla f(\boldsymbol{x}_0)\\ 1.\; 第 \;k\; 次迭代:\;(k\ge 0)\\ \quad (a).\; 令\; \boldsymbol{p}_k = \boldsymbol{H}_k\nabla f(\boldsymbol{x}_k);\\ \quad (b).\;找到\; \boldsymbol{x}_{k+1} = \boldsymbol{x}_k - h_k\boldsymbol{p}_k;\\ \quad (c).\; 计算\;f(\boldsymbol{x}_{k+1})\;和\;\nabla f(\boldsymbol{x}_{k+1});\\ \quad (d).\; 由矩阵\;\boldsymbol{H}_k\;更新\; \boldsymbol{H}_{k+1}.\end{array}} \]

其中 \(\boldsymbol{p}_k\) 就是下降方向;\(h_k\) 是步长,其选择规则可以看第一节的内容。拟牛顿法与牛顿法最不一样的地方就在与步骤 \(1(d)\) 中:对序列 \(\{\boldsymbol{H}_k\}\) 进行更新需要用到步骤 \(1(c)\) 中的信息,更新 \(\boldsymbol{H}_{k}\) 的方法虽然有很多,但是需要满足一个统一的规则——拟牛顿规则/Quasi-Newton Rule

\[ \boldsymbol{H}_{k+1}(\boldsymbol{x}_{k+1}-\boldsymbol{x}_k) = \nabla f(\boldsymbol{x}_{k+1})-\nabla f(\boldsymbol{x}_k) \]

这个规则的想法来自于对于二次函数的讨论:对于二次函数 \(f(\boldsymbol{x}) = \alpha+\langle\boldsymbol{a},\boldsymbol{x}\rangle+\dfrac{1}{2}\langle\boldsymbol{A}\boldsymbol{x},\boldsymbol{x}\rangle\)\(\forall\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^n\),有 \(\nabla f(\boldsymbol{y})-\nabla f(\boldsymbol{x}) = \boldsymbol{A}(\boldsymbol{y}-\boldsymbol{x})\);我们也希望拟牛顿法构造出的 Hessian 逆近似满足这样的关系。并且我们可以证明:在严格局部极小点的邻域内,一般拟牛顿法局部收敛并且收敛速率为超线性,也就是

\[ \lVert \boldsymbol{x}_{k+1}-\boldsymbol{x}^*\rVert \leqslant C \cdot \lVert \boldsymbol{x}_k-\boldsymbol{x}^*\rVert \cdot \lVert \boldsymbol{x}_{k-n} - \boldsymbol{x}^*\rVert \]

往回看:我们想知道为什么拟牛顿法被称为变尺度法。简而言之,拟牛顿法是在范数 \(\lVert\cdot\rVert_{\boldsymbol{H}_k^{-1}}\) 下的最速下降法,我们下面仔细解释一下这个说法:

拟牛顿规则分析

共轭梯度法

共轭梯度法最早由 Hestenes 和 Stiefel 在 1952 年提出来作为解线性方程组的方法,由于解线性方程组等价于优化一个正定的二次函数,所以 1964 年 Fletcher 和 Reeves 提出了无约束极小化的共扼梯度法。下面针对二次函数的情况介绍共轭梯度法,并且给出一般情况下的共轭梯度法。

方法引入

简单来讲,先前的下降方法的所有的 \(\boldsymbol{x}_k\) 都来自 \(\operatorname*{Lin}\{\boldsymbol{x}_0,\nabla f(\boldsymbol{x}_i),i=1,\dots,k-1\}\),其中符号 \(\operatorname*{Lin}\) 表示线性张成/Linear Span 的空间,在其他材料中也经常记作 \(\operatorname*{span}\)。与其从 \(\mathbb{R}^n\) 里寻找一个 \(\boldsymbol{x}_k\),我们希望将搜索的范围缩小到一个简单求解的子空间之中,这个子空间我们将优化为 Krylov 子空间/Krylov Subspace。

对于二次函数 \(f(\boldsymbol{x}) = \alpha+\langle\boldsymbol{a},\boldsymbol{x}\rangle+\dfrac{1}{2}\langle\boldsymbol{A}\boldsymbol{x},\boldsymbol{x}\rangle\) 与最小化问题 \(\min\limits_{\boldsymbol{x}\in\mathbb{R}^n}f(\boldsymbol{x})\),设起始点为 \(\boldsymbol{x}_0\in\mathbb{R}^n\),考虑如下 Krylov 子空间:

\[ \mathcal{L}_k = \operatorname*{Lin}\{\boldsymbol{A}(\boldsymbol{x}_0-\boldsymbol{x}^*),\dots,\boldsymbol{A}^k(\boldsymbol{x}_0-\boldsymbol{x}^*)\},\; k\geqslant 1 \]

那么,共轭梯度法/Conjugate Gradient Method 生成如下序列 \(\{\boldsymbol{x}_k\}\)

\[ \boldsymbol{x}_k = \mathop{\arg\min}\{f(\boldsymbol{x})\;\vert\;\boldsymbol{x}\in\boldsymbol{x}_0+\mathcal{L}_k\},\; k\geqslant 1 \]

这个定义符合我们对共轭梯度法的期待,但是并不是很自然,也没有描述成下降算法的样子,并且下降方向也不是很明确。记每一步的下降方向为 \(\boldsymbol{p}_k\),我们希望共轭梯度法的样子是:

  • 要像一个下降算法,也就是对于 \(k=1,2,\dots\)\(\alpha_k>0\),有

    \[ \begin{align} &\boldsymbol{x}_k = \boldsymbol{x}_{k-1}+\alpha_{k-1}\boldsymbol{p}_{k-1} \\ \Rightarrow \;& (\boldsymbol{A}\boldsymbol{x}_k+\boldsymbol{a}) = (\boldsymbol{A}\boldsymbol{x}_{k-1}+\boldsymbol{a})+\alpha_{k-1}\boldsymbol{A}\boldsymbol{p}_{k-1} \\ \Rightarrow \;& \nabla f(\boldsymbol{x}_k) = \nabla f(\boldsymbol{x}_{k-1})+\alpha_{k-1}\boldsymbol{A}\boldsymbol{p}_{k-1}\qquad\qquad(*) \end{align} \]
  • 下降方向应当只根据梯度和最后一次下降方向进行更新;设 \(\beta_k>0\)\(\boldsymbol{p}_0=\nabla f(\boldsymbol{x}_0)\),对于 \(k=1,2,\dots\)

    \[ \boldsymbol{p}_k = \nabla f(\boldsymbol{x}_k)-\beta_{k-1}\boldsymbol{p}_{k-1} \]

基础理论

根据上面的想法,我们有下面的引理:

共轭梯度算法分析