第 17 章 二阶优化Chapter 17 Second-order optimization

我们在式 16.2 中看到，梯度下降可以从变分视角来解释：在该视角下，我们借助由当前梯度得到的目标函数在当前迭代值附近的线性近似。类似地，若我们不仅能取得目标函数的梯度，还能取得其 Hessian 矩阵，便可使用目标函数在当前迭代值附近的二次近似。更准确地说，给定函数 $L({\bm{w}})$，我们可以考虑极小化 $L({\bm{w}})$ 在当前迭代值 ${\bm{w}}^{t}$ 附近的二阶泰勒近似：We saw in Eq. 16.2 that gradient descent can be motivated from a variational perspective, in which we use a linear approximation of the objective around the current iterate, obtained from the current gradient. Similarly, if we have access not only to the gradient but also to the Hessian of the objective, we can use a quadratic approximation of the objective around the current iterate. More precisely, given a function $L({\bm{w}})$, we may consider minimizing the second-order Taylor approximation of $L({\bm{w}})$ around the current iterate ${\bm{w}}^{t}$,

$$L({\bm{w}})\approx L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}-{\bm{w}}^{t}\rangle+\frac{1}{2}\langle{\bm{w}}-{\bm{w}}^{t},\nabla^{2}L({\bm{w}}^{t})({\bm{w}}-{\bm{w}}^{t})\rangle.$$

牛顿法所做的，正是围绕当前迭代值 ${\bm{w}}^{t}$ 反复极小化这一二次近似，即Newton’s method simply iteratively minimizes this quadratic approximation around the current iteration ${\bm{w}}^{t}$, namely,

$$\displaystyle{\bm{w}}^{t+1}=\argmin_{{\bm{w}}\in\mathcal{W}}L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}{-}{\bm{w}}^{t}\rangle+\frac{1}{2}\langle{\bm{w}}{-}{\bm{w}}^{t},\nabla^{2}L({\bm{w}}^{t})({\bm{w}}{-}{\bm{w}}^{t})\rangle.$$

(17.1)

若Hessian 矩阵在 ${\bm{w}}^{t}$ 处正定（记作 $\nabla^{2}L({\bm{w}}^{t})\succ 0$），则该极小值点良定义且唯一（例如，当 $L$ 为严格凸时即属此情形）。此时，迭代值可解析地写为If the Hessian is positive definite at ${\bm{w}}^{t}$, which we denote by $\nabla^{2}L({\bm{w}}^{t})\succ 0$, then the minimum is well-defined and unique (this is for example the case if $L$ is strictly convex). The iterates can then be written analytically as

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\nabla^{2}L({\bm{w}}^{t})^{-1}\nabla L({\bm{w}}^{t}).$$

若Hessian 矩阵并非正定，极小值就可能无定义。此时若无视这一问题而径直采用该解析表达式，则可能带来危险：例如，当该二次函数为严格凹（即 $\nabla^{2}L({\bm{w}})\prec 0$）时，这样做所算出的可能反而成了该二次函数的极大值。If the Hessian is not positive definite, the minimum may not be defined. Ignoring this issue and taking the analytical formulation could be dangerous, as it could amount to computing the maximum of the quadratic instead if, for example, the quadratic was strictly concave (i.e., $\nabla^{2}L({\bm{w}})\prec 0$).

规避这一问题的一种简单做法，是给 Hessian 矩阵添加一个正则化项。也就是说，从变分视角出发，我们可以加入邻近项 $\frac{1}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2}$，促使迭代值保持在当前 ${\bm{w}}^{t}$ 附近。于是，这种正则化牛顿法的迭代值具有如下形式：A simple technique to circumvent this issue consists in adding a regularization term to the Hessian. Namely, from a variational viewpoint, we can add a proximity term $\frac{1}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2}$, encouraging it to stay close to the current ${\bm{w}}^{t}$. The iterates of this regularized Newton method then take the form

$$\displaystyle{\bm{w}}^{t+1}=\argmin_{{\bm{w}}\in\mathcal{W}} \$$	$$\displaystyle L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}{-}{\bm{w}}^{t}\rangle+\frac{1}{2}\langle{\bm{w}}{-}{\bm{w}}^{t},\nabla^{2}L({\bm{w}}^{t})({\bm{w}}{-}{\bm{w}}^{t})\rangle$$
	$$\displaystyle+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2},$$

其中 $\eta^{t}$ 控制正则化强度。假设 $\eta^{t}>0$ 足够强，足以使 $\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I}$ 正定，则有where $\eta^{t}$ controls the regularization strength. Assuming $\eta^{t}>0$ is strong enough to make $\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I}$ positive-definite, we have

$${\bm{w}}^{t+1}={\bm{w}}^{t}-{\bm{d}}^{t},$$

其中我们定义了方向where we defined the direction

$${\bm{d}}^{t}\coloneqq(\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t}).$$

(17.2)

规避这一问题的其他技术还包括采用三次正则化，以及修改 Hessian 矩阵的谱分解——或对特征值作阈值处理，或取其绝对值。感兴趣的读者可参阅例如 [146, 211] 以了解详情。Other techniques to circumvent this issue include using cubic regularization and modifying the spectral decomposition of the Hessian, by thresholding the eigenvalues or taking their absolute values. We refer the interested reader to, e.g., [146, 211] for more details.

我们注意到，实现这种二阶优化算法的一个主要障碍在于：即便我们能免费获取目标函数的 Hessian 矩阵，且该矩阵正定，要精确计算式 17.2 中的方向 ${\bm{d}}^{t}$，仍需计算梯度 $\nabla L({\bm{w}}^{t})$ 的逆 Hessian—向量乘积（IHVP）。而精确计算该乘积，需要求解一个线性方程组We observe a main impediment for implementing such a second-order optimization algorithm: even if we had access to the Hessian of the objective for free and this Hessian was positive definite, computing the exact direction ${\bm{d}}^{t}$ in Eq. 17.2 requires computing an inverse-Hessian vector product (IHVP) with the gradient $\nabla L({\bm{w}}^{t})$. Doing so exactly requires solving a linear system

$$(\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I}){\bm{d}}^{t}=\nabla L({\bm{w}}^{t}),$$

这一求解原则上需要 $O(P^{3})$ 时间。不过在实践中，我们可以近似计算 IHVP，具体如 9.4 节所述。which a priori takes $O(P^{3})$ time. In practice, however, we can compute IHVPs approximately, as explained in Section 9.4.

实现牛顿法虽需付出更高的计算代价，却也能换来更快的收敛速率。简言之，若牛顿法从一个足够接近 $\mu$-强凸函数极小值点 ${\bm{w}}^{\star}$ 的点 ${\bm{w}}^{0}\in\mathcal{W}$ 出发，且该函数的Hessian 矩阵满足 $M$-Lipschitz 连续性（即 $\|{\bm{w}}^{0}-{\bm{w}}^{*}\|_{2}\leq\frac{2\mu}{3M}$），则牛顿法以二次速率收敛 [146]，即 $R(t)\leq O(\exp(-2^{t}))$（性能保证的简要介绍见 15.5 节）。这远胜于梯度下降。正是这种高效性推动了内点法的发展；借助对数障碍罚函数，内点法成为约束优化领域的一项突破 [146]。While implementing Newton’s method comes at a higher computational cost, it can also benefit from faster convergence rates. Briefly, if Newton’s method is initialized at a point ${\bm{w}}^{0}\in\mathcal{W}$ close enough to the minimizer ${\bm{w}}^{\star}$ of a $\mu$-strongly convex function with $M$-Lipschitz continuous Hessian (namely $\|{\bm{w}}^{0}-{\bm{w}}^{*}\|_{2}\leq\frac{2\mu}{3M}$), then Newton’s method converges at a quadratic rate [146], that is, $R(t)\leq O(\exp(-2^{t}))$ (see Section 15.5 for a brief introduction to performance guarantees). This is far superior to gradient descent. Such an efficiency motivated the development of interior point methods, that have been a breakthrough in constrained optimization, thanks to the use of log-barrier penalties [146].

实践中，我们未必能获得足够接近极小值点的初始点。此时，即便对于严格凸函数——牛顿步对其是良定义的——取 ${\bm{w}}^{t+1}={\bm{w}}^{t}-{\bm{d}}^{t}$ 也未必能保证目标值下降。不过，方向 ${\bm{d}}^{t}$ 仍可能按下文定义构成一个下降方向。In practice, we may not have access to an initial point close enough to the minimizer. In that case, even for strictly convex functions for which Newton’s steps are well-defined, taking ${\bm{w}}^{t+1}={\bm{w}}^{t}-{\bm{d}}^{t}$ may not ensure a decrease of the objective values. Nevertheless, the direction ${\bm{d}}^{t}$ may define a descent direction as defined below.

对于未加正则化的牛顿法，只要 $\nabla L({\bm{w}}^{t})\neq 0$ 且 $\nabla^{2}L({\bm{w}}^{t})\succ 0$，则 ${\bm{d}}^{t}=\nabla^{2}L({\bm{w}}^{t})^{-1}\nabla L({\bm{w}}^{t})$ 就是 ${\bm{w}}^{t}$ 处的一个下降方向。若 $\nabla^{2}L({\bm{w}}^{t})\not\succ 0$，则选取 $\eta^{t}>0$ 使 $\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I}\succ 0$，同样能保证 ${\bm{d}}^{t}=-(\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t})$ 是下降方向（只要 $\nabla L({\bm{w}}^{t})\neq 0$）。于是，牛顿法通常会配合一种线搜索方法，该方法试图取如下形式的步长For Newton’s method without regularization, ${\bm{d}}^{t}=\nabla^{2}L({\bm{w}}^{t})^{-1}\nabla L({\bm{w}}^{t})$ is then a descent direction at ${\bm{w}}^{t}$, as long as $\nabla L({\bm{w}}^{t})\neq 0$ and $\nabla^{2}L({\bm{w}}^{t})\succ 0$. If $\nabla^{2}L({\bm{w}}^{t})\not\succ 0$, choosing $\eta^{t}>0$ such that $\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I}\succ 0$, also ensures that ${\bm{d}}^{t}=-(\nabla^{2}L({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t})$ is a descent direction (as long as $\nabla L({\bm{w}}^{t})\neq 0$). Newton’s method is then generally equipped with a linesearch method that attempts to take steps of the form

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\gamma^{t}{\bm{d}}^{t}$$

其中 $\gamma^{t}$ 取为 $\{\rho^{\tau},\tau\in\mathbb{N}\}$（对 $\rho\in(0,1)$）中使目标函数满足充分下降条件（例如对 $c\in(0,1)$ 有）的最大步长，with $\gamma^{t}$ chosen as the largest stepsize among $\{\rho^{\tau},\tau\in\mathbb{N}\}$ for $\rho\in(0,1)$ until a sufficient decrease of the objective is satisfied such as, for $c\in(0,1)$,

$$L({\bm{w}}^{t}-\gamma^{t}{\bm{d}}^{t})\leq L({\bm{w}}^{t})-c\gamma^{t}\langle\nabla L({\bm{w}}^{t}),\nabla^{2}L({\bm{w}}^{t})^{-1}\nabla L({\bm{w}}^{t})\rangle.$$

对于强凸函数，这种实现会分为两个阶段：第一阶段用步长 $\gamma^{t}<1$ 对牛顿步进行“阻尼”；第二阶段则超快收敛，此时步长取 $\gamma^{t}=1$，目标函数迅速下降。即便离最优点很远，牛顿方向也能适应目标函数的局部几何，从而比常规梯度下降更快收敛，详见下文。For strongly convex functions, such an implementation exhibits two phases: a first phase during which Newton’s steps are “damped” by using a stepsize $\gamma^{t}<1$ and a second phase of super-fast convergence during which stepsizes $\gamma^{t}=1$ are taken, and the objective decreases very fast. Even far from the optimum, Newton directions can advantageously adapt to the local geometry of the objective to speed up convergence compared to a regular gradient descent as explained below.

为理解牛顿法相对于梯度下降的高效性，考虑一个简单二次函数的最小化问题To understand the efficiency of Newton’s method compared to gradient descent, consider the minimization of a simple quadratic

$$L({\bm{w}})=\frac{1}{2}aw_{1}^{2}+\frac{1}{2}bw_{2}^{2}$$

对于 $a\gg b\geq 0$，如图 17.1 所示。梯度下降沿 $\nabla L({\bm{w}})=(aw_{1},bw_{2})^{\top}$ 的方向移动，其步长受第一个坐标变化的制约，因而产生振荡。若仅用 $(a,b)$ 对梯度作缩放，即取如下形式的步长for $a\gg b\geq 0$, as illustrated in Fig. 17.1. A gradient descent moves along the directions $\nabla L({\bm{w}})=(aw_{1},bw_{2})^{\top}$ and its stepsize is limited by the variations in the first coordinate leading to some oscillations. If we were simply rescaling the gradient by $(a,b)$, i.e., taking steps of the form

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\gamma\operatorname{\bf diag}(a^{-1},b^{-1})\nabla L({\bm{w}}^{t}),$$

则两个坐标上的变化都归一为一，步长可直接取 $\gamma=1$，从而直接得到 ${\bm{w}}^{\star}$。换言之，让方向的几何结构与目标函数诱导的几何结构相匹配，便能绕开振荡。the variations in both coordinates would be normalized to one and the stepsize could simply be chosen to $\gamma=1$ to directly get ${\bm{w}}^{\star}$. In other words, by adapting the geometry of the directions with the geometry induced by the objective, we can circumvent the oscillations.

牛顿法正是如此：它借助 Hessian 矩阵的逆来修正梯度方向。形式上，在第 $t$ 次迭代时，考虑修改后的目标函数That’s exactly what Newton’s method does by modifying the gradient direction using the inverse of the Hessian. Formally, at iteration $t$, consider the modified objective

$$\tilde{L}({\bm{v}})=L(A{\bm{v}})\ \text{for}\ A=\nabla^{2}L({\bm{w}}^{t})^{-1/2},$$

其中 $L$ 严格凸，$A$ 为 Hessian 矩阵的逆矩阵平方根。不难验证，牛顿步等价于在 $\tilde{L}$ 上做梯度步，即with $L$ strictly convex and $A$ the inverse matrix square root of the Hessian. One easily verifies that a Newton step is equivalent to a gradient step on $\tilde{L}$, that is,

$${\bm{v}}^{t+1}={\bm{v}}^{t}-\nabla\tilde{L}({\bm{v}}^{t})\ \iff\ {\bm{w}}^{t+1}={\bm{w}}^{t}-(\nabla^{2}L({\bm{w}}^{t}))^{-1}\nabla L({\bm{w}}^{t})$$

其中where

$${\bm{w}}^{t}=A{\bm{v}}^{t}=\nabla^{2}L({\bm{w}}^{t})^{-1/2}{\bm{v}}^{t}.$$

在 $A$ 诱导的几何结构下，目标函数通常具有更好的条件数，如图 17.1 所示。这正是牛顿法高效的原因。特别地，对任意强凸二次函数，一次牛顿步即可在一次迭代中抵达最优点，而梯度步则可能需要多得多的迭代。In the geometry induced by $A$, the objective is generally better conditioned as illustrated in Fig. 17.1. This explains the efficiency of Newton’s method. In particular for any strongly convex quadratic, a Newton step reaches the optimum in one iteration, while a gradient step can take many more iterations.

现在考虑一个期望损失Consider now an expected loss

$$\min_{{\bm{w}}\in\mathcal{W}}\ L({\bm{w}})\coloneqq\mathbb{E}_{S\sim\rho}\left[L({\bm{w}};S)\right].$$

此时，可以仿照构造梯度估计的方式构造Hessian 矩阵的估计，依据是In that case, an estimate of the Hessian can be constructed just like for the gradient using that

$$\mathbb{E}_{S\sim\rho}\left[\nabla^{2}L({\bm{w}};S)\right]=\nabla^{2}L({\bm{w}}).$$

于是记Les us then denote

$$g({\bm{w}};S)\approx\nabla L({\bm{w}}),\quad H({\bm{w}};S^{\prime})\approx\nabla^{2}L({\bm{w}})$$

它们分别是梯度与Hessian 矩阵的某种随机估计，其中 $S,S^{\prime}$ 独立地采自 $\rho$，或来自小批量大小各异的小批量近似。随机牛顿法的一种实现如下：some stochastic estimates of the gradient and the Hessian, respectively, with $S,S^{\prime}$ independently drawn from $\rho$ or from mini-batch approximations with varying mini-batch sizes. One implementation of a stochastic Newton method can then be

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\gamma^{t}(H({\bm{w}}^{t};S^{\prime})+\eta^{t}\operatorname{I})^{-1}g({\bm{w}}^{t};S),$$

其中 $\eta^{t}\geq 0$ 满足 $(H({\bm{w}}^{t};S^{\prime})+\eta^{t}\operatorname{I})^{-1}\succ 0$，且 $\gamma^{t}$ 固定不变，或按某个充分下降条件选取。更多细节与变体，请感兴趣的读者参阅例如 [212]。for $\eta^{t}\geq 0$ such that $(H({\bm{w}}^{t};S^{\prime})+\eta^{t}\operatorname{I})^{-1}\succ 0$ and $\gamma^{t}$ fixed or chosen to satisfy some sufficient decrease condition. We refer the interested reader to, e.g., [212], for more details and variants.

考虑如下形式的复合目标函数Consider a composite objective of the form

$$L({\bm{w}})\coloneqq\ell(f({\bm{w}})),$$

其中 $\ell:\mathcal{M}\rightarrow\mathbb{R}$ 是凸函数，例如作用在给定样本上的凸损失函数；$f:\mathcal{W}\rightarrow\mathcal{M}$ 是非线性函数，例如参数为 ${\bm{w}}\in\mathcal{W}$ 的神经网络，并在同一样本上求值。我们已经看到，梯度下降与牛顿法分别相当于在当前迭代值 ${\bm{w}}^{t}$ 附近对 $L({\bm{w}})$ 使用线性近似与二次近似。作为两者之间的折中，Gauss-Newton 方法使用 $f$ 在 ${\bm{w}}^{t}$ 附近的线性化where $\ell:\mathcal{M}\rightarrow\mathbb{R}$ is a convex function, such as a convex loss function applied on a given sample, and $f:\mathcal{W}\rightarrow\mathcal{M}$ is a nonlinear function, such as a neural network with parameters ${\bm{w}}\in\mathcal{W}$, evaluated on the same sample. We saw that gradient descent and Newton’s method amount to using linear and quadratic approximations of $L({\bm{w}})$ around the current iterate ${\bm{w}}^{t}$, respectively. As a middle ground between the two, the Gauss-Newton method uses the linearization of $f$ around ${\bm{w}}^{t}$

$$f({\bm{w}})\approx f({\bm{w}}^{t})+\partial f({\bm{w}}^{t})({\bm{w}}-{\bm{w}}^{t})$$

但保持 $\ell$ 不变，从而得到目标函数but keeps $\ell$ as is to obtain the objective

$$\displaystyle{\bm{w}}^{t+1}$$	$$\displaystyle\coloneqq\argmin_{{\bm{w}}\in\mathcal{W}}\ell(f({\bm{w}}^{t})+\partial f({\bm{w}}^{t})({\bm{w}}-{\bm{w}}^{t}))$$
	$$\displaystyle=\argmin_{{\bm{w}}\in\mathcal{W}}\ell(\partial f({\bm{w}}^{t}){\bm{w}}+f({\bm{w}}^{t})-\partial f({\bm{w}}^{t}){\bm{w}}^{t})$$
	$$\displaystyle=\argmin_{{\bm{w}}\in\mathcal{W}}\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t}),$$

其中我们定义了简写 $J^{t}\coloneqq\partial f({\bm{w}}^{t})$ 与 ${\bm{\delta}}^{t}\coloneqq f({\bm{w}}^{t})-\partial f({\bm{w}}^{t}){\bm{w}}^{t}$。我们把 $\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$ 称为 $L=\ell\circ f$ 在 ${\bm{w}}^{t}$ 处的部分线性化，以区别于梯度下降中所使用的 $L$ 的完全线性化。where we defined the shorthands $J^{t}\coloneqq\partial f({\bm{w}}^{t})$ and ${\bm{\delta}}^{t}\coloneqq f({\bm{w}}^{t})-\partial f({\bm{w}}^{t}){\bm{w}}^{t}$. We call $\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$ the partial linearization of $L=\ell\circ f$ at ${\bm{w}}^{t}$, as opposed to the full linearization of $L$ used in gradient descent.

由于凸函数与线性函数的复合仍是凸函数，因此即便 $L({\bm{w}})$ 是非凸的，该目标函数也是凸的。实践中，我们常常加入一个邻近项作为正则化，从而定义Since the composition of a convex function and of linear function is convex, this objective is convex even if $L({\bm{w}})$ is nonconvex. In practice, we often add a proximity term as regularization to define

$${\bm{w}}^{t+1}\coloneqq\argmin_{{\bm{w}}\in\mathcal{W}}\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2}.$$

(17.3)

我们可以把这一更新视为对邻近点更新的近似We can see this update as an approximation of the proximal point update

$$\argmin_{{\bm{w}}\in\mathcal{W}}L({\bm{w}})+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2},$$

其中 $L({\bm{w}})$ 已被替换为其部分线性化。用基于梯度的求解器求解式 17.3，需要计算 ${\bm{w}}\mapsto\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$ 的梯度，即 ${\bm{w}}\mapsto(J^{t})^{*}\nabla\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$。因此，用自动微分计算该梯度需要执行一次前向传播来计算 JVP $J^{t}{\bm{w}}$，再执行一次反向传播来计算 VJP $(J^{t})^{*}\nabla\ell({\bm{z}})$。关于这些算子的介绍见 2.3 节，关于自动微分的介绍见第 8 章。where $L({\bm{w}})$ has been replaced by its partial linearization. Solving Eq. 17.3 using gradient-based solvers requires computing the gradient of ${\bm{w}}\mapsto\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$, which is ${\bm{w}}\mapsto(J^{t})^{*}\nabla\ell(J^{t}{\bm{w}}+{\bm{\delta}}^{t})$. Computing this gradient by autodiff therefore requires performing a forward pass to compute the JVP $J^{t}{\bm{w}}$ and a backward pass to compute the VJP $(J^{t})^{*}\nabla\ell({\bm{z}})$. See Section 2.3 for an introduction to these operators and Chapter 8 for an introduction to autodiff.

带任意凸外层损失的 Gauss-Newton 方法，常被称为修正 Gauss-Newton 方法 [145] 或 prox-linear 方法 [57]。经典的 Gauss-Newton 方法与 Levenberg-Marquardt 方法 [118, 126] 源自非线性最小二乘，当 $\ell({\bm{z}})$ 为二次函数时，便可还原出这两种方法 [105]，例如 $\ell({\bm{z}})\coloneqq\frac{1}{2}\|{\bm{z}}-{\bm{y}}\|^{2}_{2}$，其中 ${\bm{y}}$ 为某个参考目标。在经典意义上，Gauss-Newton 方法对应于不使用正则化的情形（即 $\eta^{t}=0$），而 Levenberg-Marquardt 法则使用了正则化（通常称为阻尼，并可能在不同迭代中改变 $\eta^{t}$）。关于各种变体，可参见例如 [133] 的综述。The Gauss-Newton method with arbitrary convex outer loss is often called modified Gauss-Newton [145] or prox-linear [57]. The classical Gauss-Newton and Levenberg-Marquardt [118, 126] methods originate from nonlinear least-squares and are recovered when $\ell({\bm{z}})$ is quadratic [105], such as $\ell({\bm{z}})\coloneqq\frac{1}{2}\|{\bm{z}}-{\bm{y}}\|^{2}_{2}$, for ${\bm{y}}$ some reference target. The Gauss-Newton method corresponds classically to not using regularization (i.e., $\eta^{t}=0$) and the Levenberg-Marquardt method uses regularization (usually called damping, potentially changing $\eta^{t}$ across iterations). See e.g., [133], for a survey of different variants.

Gauss-Newton 方法的另一种变体，是把凸损失 $\ell$ 替换为它在 ${\bm{z}}^{t}\coloneqq f({\bm{w}}^{t})$ 附近的二次近似，Another variant of the Gauss-Newton method consists in replacing the convex loss $\ell$ with its quadratic approximation around ${\bm{z}}^{t}\coloneqq f({\bm{w}}^{t})$,

$$q^{t}({\bm{z}})\coloneqq\ell({\bm{z}}^{t})+\langle\nabla\ell({\bm{z}}^{t}),{\bm{z}}-{\bm{z}}^{t}\rangle+\frac{1}{2}\langle{\bm{z}}-{\bm{z}}^{t},\nabla^{2}\ell({\bm{z}}^{t})({\bm{z}}-{\bm{z}}^{t})\rangle\approx\ell({\bm{z}})$$

由此定义如下更新to define the update

$${\bm{w}}^{t+1}\coloneqq\argmin_{{\bm{w}}\in\mathcal{W}}q^{t}(J^{t}{\bm{w}}+{\bm{\delta}}^{t})+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|_{2}^{2}.$$

注意，此处已将 $\ell$ 替换为其二次近似 $q^{t}$。该目标始终为凸二次函数；相比之下，式 17.1 中牛顿法的目标函数在 $f$ 为非线性时，先验地是一个非凸二次函数。简单的计算表明Notice that $\ell$ has been replaced by its quadratic approximation $q^{t}$. This objective is always a convex quadratic, unlike the objective of the Newton method in Eq. 17.1, which is a priori a nonconvex quadratic, if $f$ is nonlinear. Simple calculations show that

$$\displaystyle{\bm{w}}^{t+1}=\argmin_{{\bm{w}}\in\mathcal{W}} \$$	$$\displaystyle L({\bm{w}}^{t})+\langle{\bm{q}}^{t},J^{t}({\bm{w}}{-}{\bm{w}}^{t})\rangle$$
	$$\displaystyle+\frac{1}{2}\langle{\bm{w}}{-}{\bm{w}}^{t},(J^{t})^{*}Q^{t}J^{t}({\bm{w}}{-}{\bm{w}}^{t})\rangle+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|^{2}_{2}$$

其中 ${\bm{q}}^{t}\coloneqq\nabla\ell(f({\bm{w}}^{t}))\in\mathcal{M}=\mathbb{R}^{Z},Q^{t}\coloneqq\nabla^{2}\ell(f({\bm{w}}^{t}))\in\mathbb{R}^{Z\times Z}.$。其闭式解为where ${\bm{q}}^{t}\coloneqq\nabla\ell(f({\bm{w}}^{t}))\in\mathcal{M}=\mathbb{R}^{Z},Q^{t}\coloneqq\nabla^{2}\ell(f({\bm{w}}^{t}))\in\mathbb{R}^{Z\times Z}.$ The closed form solution is

$$\displaystyle{\bm{w}}^{t+1}$$	$$\displaystyle={\bm{w}}^{t}-((J^{t})^{*}Q^{t}J^{t}+\eta^{t}\operatorname{I})^{-1}(J^{t})^{*}{\bm{q}}^{t}$$
	$$\displaystyle={\bm{w}}^{t}-(\nabla^{2}_{\mathrm{GN}}(\ell\circ f)({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t}),$$

这里我们使用了 $L=\ell\circ f$ 的（广义）Gauss-Newton 矩阵，该矩阵的定义见 9.2 节。where we used the (generalized) Gauss-Newton matrix of $L=\ell\circ f$, defined in Section 9.2.

与牛顿法类似，Gauss-Newton 方法若单独使用，其迭代值也可能发散。不过，方向 $-(\nabla^{2}_{\mathrm{GN}}(\ell\circ f)({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t})$ 对任意 $\eta^{t}>0$ 都构成下降方向，并且可以与步长 $\gamma^{t}$（通常通过线搜索选取）结合，得到如下形式的迭代值Similarly to Newton’s method, the iterates of a Gauss-Newton method may diverge when used alone. However, the direction $-(\nabla^{2}_{\mathrm{GN}}(\ell\circ f)({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t})$ defines a descent direction for any $\eta^{t}>0$ and can be combined with a stepsize $\gamma^{t}$ (typically chosen using a linesearch) to obtain iterates of the form

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\gamma^{t}(\nabla^{2}_{\mathrm{GN}}(\ell\circ f)({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t}).$$

在深度学习中，目标函数通常由损失函数与网络函数的复合关于样本的期望构成：In deep learning, the objective generally consists in an expectation over samples of the composition between a loss function and a network function:

$$L({\bm{w}})=\mathbb{E}_{S\sim\rho}[L({\bm{w}};S)]=\mathbb{E}_{(X,Y)\sim\rho}[\ell(f({\bm{w}};X);Y)]$$

其中 $S=(X,Y)$ 表示由输入 $X$ 及其对应标签 $Y$ 构成的样本对。此时，正如 9.2 节已讨论过的，Gauss-Newton 矩阵 $\nabla^{2}_{\mathrm{GN}}L$ 是各个 Gauss-Newton 矩阵的期望where $S=(X,Y)$ denotes a sample pair of input $X$ with associated label $Y$. In that case, as already studied in Section 9.2, the Gauss-Newton matrix $\nabla^{2}_{\mathrm{GN}}L$ is the expectation of the individual Gauss-Newton matrices

$$\displaystyle\nabla^{2}_{\mathrm{GN}}L({\bm{w}};{\bm{x}},{\bm{y}})$$	$$\displaystyle\coloneqq\partial f({\bm{w}};{\bm{x}})^{\top}\nabla^{2}\ell(f({\bm{w}};{\bm{x}}))\partial f({\bm{w}};{\bm{x}}),$$
$$\displaystyle\nabla^{2}_{\mathrm{GN}}L({\bm{w}})$$	$$\displaystyle\coloneqq\mathbb{E}_{(X,Y)\sim\rho}[\nabla^{2}_{\mathrm{GN}}L({\bm{w}};{\bm{x}},{\bm{y}})].$$

我们可以分别用 ${\bm{g}}({\bm{w}};S)\approx\nabla L({\bm{w}})$ 和 $G({\bm{w}};S^{\prime})\approx\nabla^{2}_{\mathrm{GN}}L({\bm{w}})$ 来估计梯度与 Gauss-Newton 矩阵，其中 $S,S^{\prime}\sim\rho$，或改用小批量近似。因此，随机 Gauss-Newton 方法执行的迭代为We can estimate the gradient and the Gauss-Newton matrix by, respectively, ${\bm{g}}({\bm{w}};S)\approx\nabla L({\bm{w}})$, and $G({\bm{w}};S^{\prime})\approx\nabla^{2}_{\mathrm{GN}}L({\bm{w}})$ for $S,S^{\prime}\sim\rho$ or using mini-batch approximations. A stochastic Gauss-Newton method therefore performs iterates

$${\bm{w}}^{t+1}\coloneqq{\bm{w}}^{t}-\gamma^{t}(G({\bm{w}}^{t};S^{\prime})+\eta^{t}\operatorname{I})^{-1}{\bm{g}}({\bm{w}}^{t};S),$$

其中 $\eta^{t}\geq 0$ 与 $\gamma^{t}$ 保持固定，或按某种准则加以选取。for $\eta^{t}\geq 0$ and $\gamma^{t}$ fixed or selected to satisfy some criterion.

自然梯度下降的出发点，是如下形式的更新：Natural gradient descent is motivated by updates of the form

$${\bm{w}}^{t+1}=\argmin_{{\bm{w}}\in\mathcal{W}}L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}-{\bm{w}}^{t}\rangle+\mathrm{KL}(q_{{\bm{w}}^{t}},q_{\bm{w}}),$$

其中 $\mathrm{KL}(p,q)\coloneqq\int p({\bm{z}})\log\frac{p({\bm{z}})}{q({\bm{z}})}d{\bm{z}}$ 为 Kullback-Leibler（KL）散度。因此，与梯度下降不同，这里的邻近项衡量的是当前分布 $q_{{\bm{w}}^{t}}$ 与候选概率分布 $q_{{\bm{w}}}$ 之间的差异。上述问题在一般情况下难以求解，因为 KL 散度可能没有闭式形式。不过可以证明[5]，其二次近似具有如下简洁形式：where $\mathrm{KL}(p,q)\coloneqq\int p({\bm{z}})\log\frac{p({\bm{z}})}{q({\bm{z}})}d{\bm{z}}$ is the Kullback-Leibler (KL) divergence. Unlike gradient descent, the proximity term is therefore between the current distribution $q_{{\bm{w}}^{t}}$ and a candidate probability distribution $q_{{\bm{w}}}$. The above problem is intractable in general, as the KL may not have a closed form. Nevertheless, its quadratic approximation can be shown [5] to admit a simple form,

$$\mathrm{KL}(q_{{\bm{w}}^{t}},q_{\bm{w}})\approx\frac{1}{2}\langle{\bm{w}}-{\bm{w}}^{t},\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})({\bm{w}}-{\bm{w}}^{t})\rangle$$

其中用到了 9.3 节所研究的 Fisher 信息矩阵 $\nabla^{2}_{\mathrm{F}}L({\bm{w}})$。有了 KL 散度的这一二次近似，自然梯度下降便相当于按如下方式计算迭代值：where we used the Fisher information matrix $\nabla^{2}_{\mathrm{F}}L({\bm{w}})$, studied in Section 9.3. Equipped with this quadratic approximation of the KL divergence, natural gradient descent amounts to computing iterates as

$$\displaystyle{\bm{w}}^{t+1}\coloneqq\argmin_{{\bm{w}}\in\mathcal{W}} \$$	$$\displaystyle L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}-{\bm{w}}^{t}\rangle$$
	$$\displaystyle+\frac{1}{2}\langle{\bm{w}}-{\bm{w}}^{t},\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})({\bm{w}}-{\bm{w}}^{t})\rangle+\frac{\eta^{t}}{2}\|{\bm{w}}-{\bm{w}}^{t}\|^{2}_{2},$$

其中加入了二次邻近项，以保证解唯一。由于 $\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$ 半正定，因此这是一个严格凸问题。其闭式解为where a quadratic proximity term was added to ensure a unique solution. This is a strictly convex problem as $\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$ is positive semi-definite. The closed-form solution is

$${\bm{w}}^{t+1}={\bm{w}}^{t}-(\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})+\eta^{t}\operatorname{I})^{-1}\nabla L({\bm{w}}^{t}).$$

由于当 $p_{\bm{\theta}}$ 为指数族分布时，Gauss-Newton 矩阵与 Fisher 信息矩阵等价（命题 9.3.5），因此在这种情况下，Gauss-Newton 方法与自然梯度方法彼此一致。Because the Gauss-Newton and Fisher information matrices are equivalent when $p_{\bm{\theta}}$ is an exponential family distribution (Proposition 9.3.5), the Gauss-Newton and natural gradient methods coincide in this case.

在实践中，$\nabla L({\bm{w}}^{t})$ 是关于 $\rho$ 的期望，我们未必能获得其闭式形式；此外，$\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$ 也可能无法以闭式形式计算。为估计 Fisher 信息矩阵，我们可以利用（见 9.3 节）如下简写 ${\bm{\theta}}\coloneqq f({\bm{w}},X)$，In practice, we may not have access to $\nabla L({\bm{w}}^{t})$ in closed form as it is an expectation over $\rho$. Moreover, $\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$ may not be computable in closed form either. To estimate the Fisher information matrix, we can use that (see Section 9.3) using the shorthand ${\bm{\theta}}\coloneqq f({\bm{w}},X)$,

$$\nabla^{2}_{\mathrm{F}}L({\bm{w}})=\mathbb{E}_{X\sim\rho_{\mathcal{X}}}\mathbb{E}_{Y\sim p_{\bm{\theta}}}[\nabla L({\bm{w}};X,Y)\otimes\nabla L({\bm{w}};X,Y)].$$

然后，我们可以基于从 $\rho$ 中采样得到的 $S$，以及从 $q_{{\bm{w}}^{t}}({\bm{x}},{\bm{y}})=p_{\mathcal{X}}({\bm{x}})\rho_{{\bm{\theta}}}({\bm{y}})$ 中采样得到的 $S^{\prime}$，构建估计量 ${\bm{g}}({\bm{w}}^{t};S)\approx\nabla L({\bm{w}}^{t},S)$ 与 $F({\bm{w}}^{t};S^{\prime})\approx\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$。随机自然梯度下降便可按如下方式实现：We can then build estimates ${\bm{g}}({\bm{w}}^{t};S)\approx\nabla L({\bm{w}}^{t},S)$ and $F({\bm{w}}^{t};S^{\prime})\approx\nabla^{2}_{\mathrm{F}}L({\bm{w}}^{t})$ for $S$ sampled from $\rho$ and $S^{\prime}$ sampled from $q_{{\bm{w}}^{t}}({\bm{x}},{\bm{y}})=p_{\mathcal{X}}({\bm{x}})\rho_{{\bm{\theta}}}({\bm{y}})$. A stochastic natural gradient descent can then be implemented as

$${\bm{w}}^{t+1}={\bm{w}}^{t}-\gamma^{t}(F({\bm{w}}^{t};S^{\prime})+\eta^{t}\operatorname{I})^{-1}{\bm{g}}({\bm{w}}^{t};S),$$

其中 $\gamma^{t}$ 是步长，可由线搜索选取。where $\gamma^{t}$ is a stepsize, possibly chosen by linesearch.

在深度学习中，与 Fisher 矩阵或 Gauss-Newton 矩阵之逆相乘的计算代价可能仍然很高。为此，人们已经提出了若干近似方法，例如 KFAC [128, 31]，它利用计算效率高的结构近似来逼近这些矩阵。In deep learning, the product with the inverse Fisher or Gauss-Newton matrices can remain costly to compute. Several approximations have been proposed, such as KFAC [128, 31], which uses a computationally efficient structural approximation to these matrices.

拟牛顿法的一个著名例子是 BFGS 方法 [37, 67, 80, 183]，其缩写来自其作者的名字。BFGS 更新的依据再次源于变分观点。我们希望为目标 $h^{t}({\bm{w}})=L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}-{\bm{w}}^{t}\rangle+\frac{1}{2}\langle{\bm{w}}-{\bm{w}}^{t},Q^{t}({\bm{w}}-{\bm{w}}^{t})\rangle$ 建立一个简单的二次模型，所用的某个 $Q^{t}$ 是沿迭代过程构造出来的，而不是取作 $\nabla^{2}L({\bm{w}}^{t})$。这种二次模型的一个理想性质是，其在相邻迭代点处的梯度与原函数的梯度相匹配，即 $\nabla h^{t}({\bm{w}}^{t})=\nabla L({\bm{w}}^{t})$ 和 $\nabla h^{t}({\bm{w}}^{t-1})=\nabla L({\bm{w}}^{t-1})$。一个更简单的条件称为割线条件，它在于考虑这些向量的差，也就是确保A celebrated example of quasi-Newton method is the BFGS method [37, 67, 80, 183], whose acronym follows from its authors’ names. The rationale of the BFGS update stems once again from a variational viewpoint. We wish to build a simple quadratic model of the objective $h^{t}({\bm{w}})=L({\bm{w}}^{t})+\langle\nabla L({\bm{w}}^{t}),{\bm{w}}-{\bm{w}}^{t}\rangle+\frac{1}{2}\langle{\bm{w}}-{\bm{w}}^{t},Q^{t}({\bm{w}}-{\bm{w}}^{t})\rangle$ for some $Q^{t}$ built along the iterations rather than taken as $\nabla^{2}L({\bm{w}}^{t})$. One desirable property of such a quadratic model would be that its gradients at consecutive iterates match the gradients of the original function, i.e., $\nabla h^{t}({\bm{w}}^{t})=\nabla L({\bm{w}}^{t})$ and $\nabla h^{t}({\bm{w}}^{t-1})=\nabla L({\bm{w}}^{t-1})$. A simpler condition, called the secant condition consists in considering the differences of these vectors, that is, ensuring that

	$$\displaystyle\nabla h^{t}({\bm{w}}^{t})-\nabla h^{t}({\bm{w}}^{t-1})=\nabla L({\bm{w}}^{t})-\nabla L({\bm{w}}^{t-1})$$
$$\displaystyle\iff$$	$$\displaystyle Q^{t}({\bm{w}}^{t}-{\bm{w}}^{t-1})=\nabla L({\bm{w}}^{t})-\nabla L({\bm{w}}^{t-1})$$
$$\displaystyle\iff$$	$$\displaystyle{\bm{w}}^{t}-{\bm{w}}^{t-1}=B^{t}(\nabla L({\bm{w}}^{t})-\nabla L({\bm{w}}^{t-1})),$$

对于 $B^{t}=(Q^{t})^{-1}$。接着，可按如下方式构建满足割线方程的 $B^{t}$，即 Hessian 矩阵之逆的替代量：for $B^{t}=(Q^{t})^{-1}$. Building $B^{t}$, a surrogate of the inverse of the Hessian satisfying the secant equation, can then be done as

$$B^{t+1}\coloneqq\left(\operatorname{I}-\rho^{t}{\bm{s}}^{t}({\bm{y}}^{t})^{\top}\right)B^{t}\left(\operatorname{I}-\rho^{t}{\bm{s}}^{t}({\bm{y}}^{t})^{\top}\right)+\rho^{t}{\bm{s}}^{t}({\bm{s}}^{t})^{\top}$$

其中where

$$\displaystyle{\bm{s}}^{t}$$	$$\displaystyle\coloneqq{\bm{w}}^{t+1}-{\bm{w}}^{t}$$
$$\displaystyle{\bm{y}}^{t}$$	$$\displaystyle\coloneqq\nabla L({\bm{w}}^{t+1})-\nabla L({\bm{w}}^{t})$$
$$\displaystyle\rho^{t}$$	$$\displaystyle\coloneqq\frac{1}{\langle{\bm{s}}^{t},{\bm{y}}^{t}\rangle}.$$

BFGS 的一种典型实现会在内存中存储 $B^{t}\in\mathbb{R}^{P\times P}$；当 $P$ 很大时，这种做法的内存开销会大得令人却步。A typical implementation of BFGS stores $B^{t}\in\mathbb{R}^{P\times P}$ in memory, which is prohibitive when $P$ is large.

在实践中，人们往往更青睐 BFGS 的有限内存版本，即 LBFGS [119]。LBFGS 的关键观察是，我们并不需要在内存中具体化 $B^{t}$：我们只需要将其与梯度 $\nabla L({\bm{w}}^{t})$ 相乘即可，也就是说，可以把 $B^{t}$ 视为一个线性映射。幸运的是，只要在内存中存储 $({\bm{s}}^{1},{\bm{y}}^{1},\rho^{1}),\dots,({\bm{s}}^{t},{\bm{y}}^{t},\rho^{t})$，就能高效地计算 $B^{t}$ 与任意向量 ${\bm{v}}$ 的乘积。实践中，人们会利用过去取值的一小段历史来降低内存与计算代价。由于 LBFGS 兼具类二阶方法的好处而代价大大降低，它已成为事实上的标准算法，对于没有特殊结构的中等规模问题，其表现优于大多数其他算法 [119]。In practice, the limited-memory counterpart of BFGS, called LBFGS [119], is often preferred. The key observation of LBFGS is that we do not need to materialize $B^{t}$ in memory: we only need to multiply it with the gradient $\nabla L({\bm{w}}^{t})$. That is, we can see $B^{t}$ as a linear map. Fortunately, the product between $B^{t}$ and any vector ${\bm{v}}$ can be computed efficiently if we store $({\bm{s}}^{1},{\bm{y}}^{1},\rho^{1}),\dots,({\bm{s}}^{t},{\bm{y}}^{t},\rho^{t})$ in memory. In practice, a small history of past values is used to reduce memory and computational cost. Because LBFGS has the benefits of second-order-like methods with much reduced cost, it has become a de-facto algorithm, outperforming most other algorithms for medium-scale problems without particular structure [119].