New insights and perspectives on the natural gradient method

Introduction and Overview

The natural gradient descent approach, pioneered by Amari and collaborators (e.g. Amari, 1998), is a popular alternative to traditional gradient descent methods which has received a lot of attention over the past several decades, motivating many new and related approaches. It has been successfully applied to a variety of problems such as blind source separation (Amari and Cichocki, 1998), reinforcement learning (Peters and Schaal, 2008), and neural network training (e.g. Park et al., 2000; Martens and Grosse, 2015; Desjardins et al., 2015).

Natural gradient descent is generally applicable to the optimization of probabilistic modelsThis includes neural networks, which can be cast as conditional models., and involves the use of the so-called “natural gradient”, which is defined as the gradient times the inverse of the model’s Fisher information matrix (aka “the Fisher” ; see Section 5), in place of the standard gradient. In many applications, natural gradient descent seems to require far fewer iterations than gradient descent, making it a potentially attractive alternative method. Unfortunately, for models with very many parameters such as large neural networks, computing the natural gradient is impractical due to the extreme size of the Fisher matrix. This problem can be addressed through the use of various approximations to the Fisher (e.g Le Roux et al., 2008; Ollivier, 2015; Grosse and Salakhudinov, 2015; Martens and Grosse, 2015) that are designed to be easier to compute, to store and finally to invert, than the exact Fisher.

Natural gradient descent is classically motivated as a way of implementing steepest descent“Steepest descent” is a common synonym for gradient descent, which emphasizes the interpretation of the gradient as being the direction that descends down the loss surface most “steeply”. Here “steepness” is measured as the amount loss reduction per unit of distance traveled, where distance is measured according to some given metric. For standard gradient descent this metric is the Euclidean distance on the default parameter space. in the space of realizable distributions“Realizable distributions” means distributions which correspond to some setting of the model’s parameters. instead of the space of parameters, where distance in the distribution space is measured with a special “Riemannian metric” (Amari and Nagaoka, 2007). This metric depends only on the properties of the distributions themselves and not their parameters, and in particular is defined so that it approximates the square root of the Kullback–Leibler divergence within small neighborhoods. Under this interpretation (discussed in detail in Section 6), natural gradient descent is invariant to any smooth and invertible reparameterization of the model, putting it in stark contrast to gradient descent, whose performance is parameterization dependent.

In practice however, natural gradient descent still operates within the default parameter space, and works by computing directions in the space of distributions and then translating them back to the default space before taking a step. Because of this, the above discussed interpretation breaks down unless the step-size becomes arbitrarily small, and as discussed in Section 10, this breakdown has important implications for designing a natural gradient method that can work well in practice. Another problem with this interpretation is that it doesn’t provide any obvious reason why a step of natural gradient descent should make more progress reducing the objective than a step of standard gradient descent (assuming well-chosen step-sizes for both). Moreover, given a large step-size one also loses the parameterization invariance property of the natural gradient method, although it will still hold approximately under certain conditions which are described in Section 12.

In Section 10 we argue for an alternative view of natural gradient descent: as a type of 2nd-order methodBy “2nd-order method” we mean any iterative optimization method which generates updates as the (possibly approximate) solution of a non-trivial local quadratic model of the objective function. This extends well beyond the classical Newton’s method, and includes approaches like (L-)BFGS, and methods based on the Gauss-Newton matrix. which utilizes the Fisher as an alternative to the Hessian. As discussed in Section 7, 2nd-order methods work by forming a local quadratic approximation to the objective around the current iterate, and produce the next iterate by optimizing this approximation within some region where the approximation is believed to be accurate. According to this view, natural gradient descent ought to make more progress per step than gradient descent because it uses a local quadratic model/approximation of the objective function which is more detailed (which allows it to be less conservative) than the one implicitly used by gradient descent.

In support of this view is the fact that the Fisher can be cast as an approximation of the Hessian in at least two different ways (provided the objective has the form discussed in Section 4). First, as discussed in Section 5, it corresponds to the expected Hessian of the loss under the model’s distribution over predicted outputs (instead of the usual empirical one used to compute the exact Hessian). Second, as we establish in Section 9, it is very often equivalent to the so-called “Generalized Gauss-Newton matrix” (GGN) (discussed in Section 8), a generalization of the classical Gauss-Newton matrix (e.g. Dennis Jr and Schnabel, 1996; Ortega and Rheinboldt, 2000; Nocedal and Wright, 2006), which is a popular alternative/approximation to the Hessian that has been used in various practical 2nd-order optimization methods designed specifically for neural networks (e.g Schraudolph, 2002; Martens, 2010; Vinyals and Povey, 2012), and which may actually be a better choice than the Hessian in the context of neural net training (see Section 8.1).

Viewing natural gradient descent as a 2nd-order method is also prescriptive, since it suggests the use of various damping/regularization techniques often used in the optimization literature to account for the limited accuracy of local quadratic approximations (especially over long distances). Indeed, such techniques have been successfully applied in 2nd-order methods designed for neural networks (e.g. Martens, 2010; Martens and Grosse, 2015), where they proved crucial in achieving fast and robust performance in practice. And before that have had a long history of application in the context of practical non-linear regression procedures (Tikhonov, 1943; Levenberg, 1944; Marquardt, 1963; Moré, 1978).

The “empirical Fisher”, which is discussed in Section 11, is an approximation to the Fisher whose computation is easier to implement in practice using standard automatic-differentiation libraries. The empirical Fisher differs from the usual Fisher in subtle but important ways, which as we show in Section 11.1, make it considerably less useful as an approximation to the Fisher, or as a curvature matrix to be used in 2nd-order methods. Using the empirical Fisher also breaks some of the theory justifying natural gradient descent, although it nonetheless preserves its (approximate) parameterization invariance (as we show in Section 12). Despite these objections, the empirical Fisher has been used in many approaches, such as TONGA (Le Roux et al., 2008), and the recent spate of methods that use the diagonal of this matrix such as RMSprop (Tieleman and Hinton, 2012) and Adam (Ba and Kingma, 2015) (which we examine in Section 11.2).

A well-known and oft quoted result about stochastic natural gradient descent is that it is asymptotically “Fisher efficient” (Amari, 1998). Roughly speaking, this means that it provides an asymptotically unbiased estimate of the parameters with the lowest possible variance among all unbiased estimators (that see the same amount of data), thus achieving the best possible expected objective function value. Unfortunately, as discussed in Section 14.1, this result comes with several important caveats which significantly limit its applicability. Moreover, even when it is applicable, it only provides an asymptotically accurate characterization of the method, which may not usefully describe its behavior given a finite number of iterations.

To address these issues we build on the work of Murata (1998) in Section 14.2 and Section 14.3 to develop a more powerful convergence theory for stochastic 2nd-order methods (including natural gradient descent) as applied to convex quadratic objectives. Our results provide a more precise expression for the convergence speed of such methods than existing results doAlternate versions of several of our results have appeared in the literature before (e.g. Polyak and Juditsky, 1992; Bordes et al., 2009; Moulines and Bach, 2011). These formulations tend to be more general than ours (we restrict to the quadratic case), but also less precise and interpretable. In particular, previous results tend to omit asymptotically negligible terms that are nonetheless important pre-asymptotically. Or when they do include these terms, their bounds are simultaneously looser and more complicated than our own, perhaps owing to their increased generality. We discuss connections to some prior work on convergence bounds in Sections 14.2.2 and 14.3.2., and properly account for the effect of the starting point. And as we discuss in Section 14.2.1 and Section 14.3.1 they imply various interesting consequences about the relative performance of various 1st and 2nd-order stochastic optimization methods. Perhaps the most interesting conclusion of this analysis is that, while stochastic gradient descent with Polyak-style parameter averaging achieves the same asymptotic convergence speed as stochastic natural gradient descent (and is thus also “Fisher efficient”, as was first shown by Polyak and Juditsky (1992)), stochastic 2nd-order methods can possess a much more favorable dependence on the starting point, which means that they can make much more progress given a limited iteration budget. Another interesting observation made in our analysis is that stochastic 2nd-order methods that use a decaying learning rate (of the form $\alpha_{k}=1/(k+a+1)$) can, for certain problems, achieve an asymptotic objective function value that is better than that achieved in the same number of iterations by stochastic gradient descent (with a similar decaying learning rate), by a large constant factor.

Unfortunately, these convergence theory results fail to explain why 2nd-order optimization with the GGN/Fisher works so much better than classical 2nd-order schemes based on the Hessian for neural network training (Schraudolph, 2002; Martens, 2010; Vinyals and Povey, 2012). In Section 15 we propose several important open questions in this direction that we leave for future work.

Table of Notation

Neural Networks

where $a_{0}=x$. Here, $a_{i}$ is the vector of values (“activities”) of the network’s $i$-th layer, and $\phi_{i}(\cdot)$ is the vector-valued non-linear “activation function” computed at layer $i$, and is often given by some simple scalar function applied coordinate-wise.

Note that most of the results discussed in this document will apply to the more general setting where $f(x,\theta)$ is an arbitrary differentiable function (in both $x$ and $\theta$).

Supervised Learning Framework

The goal of optimization in the context of supervised learning is to find some setting of $\theta$ so that, for each input $x$ in the training set, the output of the network (which we will sometimes call its “prediction”) matches the given target outputs as closely as possible, as measured by some loss. In particular, given a training set $S$ consisting of pairs $(x,y)$, we wish to minimize the objective function

where $L(y,z)$ is a “loss function” which measures the disagreement between $y$ and $z$.

The prediction $f(x,\theta)$ may be a guess for $y$, in which case $L$ might measure the inaccuracy of this guess (e.g. using the familiar squared error $\frac{1}{2}\|y-z\|^{2}$). Or $f(x,\theta)$ could encode the parameters of some simple predictive distribution. For example, $f(x,\theta)$ could be the set of probabilities which parameterize a multinomial distribution over the possible discrete values of $y$, with $L(y,f(x,\theta))$ being the negative log probability of $y$ under this distribution.

KL Divergence Objectives

The natural gradient method of Amari (1998) can potentially be applied to any objective function which measures the performance of some statistical model. However, it enjoys richer theoretical properties when applied to objective functions based on the KL divergence between the model’s distribution and the target distribution, or certain approximations/surrogates of these. In this section we will establish the basic notation and properties of these objective functions, and discuss the various ways in which they can be formulated. Each of these formulations will be analogous to a particular formulation of the Fisher information matrix and natural gradient (as defined in Section 5), which will differ in subtle but important ways.

In the idealized setting, input vectors $x$ are drawn independently from a target distribution $Q_{x}$ with density function $q(x)$, and the corresponding (target) outputs $y$ from a conditional target distribution $Q_{y|x}$ with density function $q(y|x)$.

We define the goal of learning as the minimization of the KL divergence from the target joint distribution $Q_{x,y}$, whose density is $q(y,x)=q(y|x)q(x)$, to the learned distribution $P_{x,y}(\theta)$, whose density is $p(x,y|\theta)=p(y|x,\theta)q(x)$. (Note that the second $q(x)$ is not a typo here, since we are not learning the distribution over $x$, only the conditional distribution of $y$ given $x$.) Our objective function is thus

This is equivalent to the expected KL divergence $\operatorname{E}_{Q_{x}}[\operatorname{KL}(Q_{y|x}\|P_{y|x}(\theta))]$ as can be seen by

It is often the case that we only have samples from $Q_{x}$ and no direct knowledge of its density function. Or the expectation w.r.t. $Q_{x}$ in eqn. 2 may be too difficult to compute. In such cases, we can substitute an empirical training distribution $\hat{Q}_{x}$ in for $Q_{x}$, which is given by a set $S_{x}$ of samples from $Q_{x}$. This results in the objective

Provided that $q(y|x)$ is known for each $x$ in $S_{x}$, and that $\operatorname{KL}(Q_{y|x}\|P_{y|x}(\theta))$ can be efficiently computed, we can use the above expression as our objective. Otherwise, as is often the case, we might only have access to a single sample $y$ from $Q_{y|x}$ for each $x\in S_{x}$, giving an empirical training distribution $\hat{Q}_{y|x}$. Substituting this in for $Q_{y|x}$ gives the objective function

where we have extended $S_{x}$ to a set $S$ of the $(x,y)$ pairs (which agrees with how $S$ was defined in Section 3). Here, the proportionality is with respect to $\theta$, and it hides an additive constant which is technically infinityThe constant corresponds to the differential entropy of the Dirac delta distribution centered at $y$. One can think of this as approaching infinity under the limit-based definition of the Dirac.. This is effectively the same objective that is minimized in standard maximum likelihood learning.

This kind of objective function fits into the general supervised learning framework described in Section 3 as follows. We define the learned conditional distribution $P_{y|x}(\theta)$ to be the composition of the deterministic prediction function $f(x,\theta)$ (which may be a neural network), and an “output” conditional distribution $R_{y|z}$ (with associated density function $r(y|z)$), so that

We then define the loss function as $L(y,z)=-\log r(y|z)$.

Given a loss function $L$ which is not explicitly defined this way one can typically still find a corresponding $R$ to make the definition apply. In particular, if $\exp(-L(y,z))$ has the same finite integral w.r.t. $y$ for each $z$, then one can define $R$ by taking $r(y|z)\propto\exp(-L(y,z))$, where the proportion is w.r.t. both $y$ and $z$.

Various Definitions of the Natural Gradient and the Fisher Information Matrix

The Fisher information matrix $F$ of $P_{x,y}(\theta)$ w.r.t. $\theta$ (aka the “Fisher”) is given by

where gradients and Hessians are taken w.r.t. $\theta$. It can be immediately seen from the first of these expressions for $F$ that it is positive semi-definite (PSD) (since it’s the expectation of something which is trivially PSD, a vector outer-product). And from the second expression we can see that it also has the interpretation of being the negative expected Hessian of $\log p(x,y|\theta)$.

The usual definition of the natural gradient (Amari, 1998) which appears in the literature is

where $F$ is the Fisher and $h$ is the objective function.

Because $p(x,y|\theta)=p(y|x,\theta)q(x)$, where $q(x)$ doesn’t depend on $\theta$, we have

and so $F$ can also be written as the expectation (w.r.t. $Q_{x}$) of the Fisher information matrix of $P_{y|x}(\theta)$ as follows:

In Amari (1998), this version of $F$ is computed explicitly for a basic perceptron model (basically a neural network with 0 hidden layers) in the case where $Q_{x}=N(0,I)$. However, in practice the real $q(x)$ may not be directly available, or it may be difficult to integrate $H_{\log p(y|x,\theta)}$ over $Q_{x}$. For example, the conditional Hessian $H_{\log p(y|x,\theta)}$ corresponding to a multi-layer neural network may be far too complicated to be analytically integrated, even for a very simple $Q_{x}$. In such situations $Q_{x}$ may be replaced with its empirical version $\hat{Q}_{x}$, giving

This is the version of $F$ considered in Park et al. (2000).

From these expressions we can see that when $L(y,z)=-\log r(y|z)$ (as in Section 4), the Fisher has the interpretation of being the expectation under $P_{x,y}$ of the Hessian of $L(y,f(x,\theta))$:

Meanwhile, the Hessian $H$ of $h$ is also given by the expected value of the Hessian of $L(y,f(x,\theta))$, except under the distribution $\hat{Q}_{x,y}$ instead of $P_{x,y}$ (where $\hat{Q}_{x,y}$ is given by the density function $\hat{q}(x,y)=\hat{q}(y|x)\hat{q}(x)$). In other words

Thus $F$ and $H$ can be seen as approximations of each other in some sense.

Geometric Interpretation

The negative gradient $-\nabla h$ can be interpreted as the steepest descent direction for $h$ in the sense that it yields the greatest instantaneous rate of reduction in $h$ per unit of change in $\theta$, where change in $\theta$ is measured using the standard Euclidean norm $\|\cdot\|$. More formally we have

This interpretation highlights the strong dependence of the gradient on the Euclidean geometry of the parameter space (as defined by the norm $\|\cdot\|$).

One way to motivate the natural gradient is to show that it (or more precisely its negation) can be viewed as a steepest descent direction, much like the negative gradient can be, except with respect to a metric that is intrinsic to the distributions being modeled, as opposed to the default Euclidean metric which is tied to the given parameterization. In particular, the natural gradient can be derived by adapting the steepest descent formulation to use an alternative definition of (local) distance based on the “information geometry” (Amari and Nagaoka, 2000) of the space of probability distributions. The particular distance functionNote that this is not a formal “distance” function in the usual sense since it is not symmetric. which gives rise to the natural gradient turns out to be

To formalize this, one can use the well-known connection between the KL divergence and the Fisher, given by the Taylor series approximation

where “$O(d^{3})$” is short-hand to mean terms that are order 3 or higher in the entries of $d$. Thus, $F$ defines the local quadratic approximation of this distance, and so gives the mechanism of local translation between the geometry of the space of distributions, and that of the original parameter space with its default Euclidean geometry.

To make use of this connection, Arnold et al. (2011) proves for general PSD matrices $A$ that

where the notation $\|v\|_{B}$ is defined by $\|v\|_{B}=\sqrt{v^{\top}Bv}$. Taking $A=\frac{1}{2}F$ and using the above Taylor series approximation to establish that

as $\epsilon\to 0$, (Arnold et al., 2011) then proceed to show that

Thus the negative natural gradient is indeed the steepest descent direction in the space of distributions where distance is measured in small local neighborhoods by the KL divergence.

By using the smoothly varying PSD matrix $F$ to locally define a metric tensor at every point in parameter space, a Riemannian manifold can be generated over the space of distributions. Note that the associated metric of this space won’t be the square root of the KL divergence (this isn’t even a valid metric), although it will be “locally equivalent” to it in the sense that the two functions will approximate each other within a small local neighborhood.

2nd-order Optimization

Gradient descent, the canonical 1st-order method, can be viewed in the framework of 2nd-order methods as making the choice $B_{k}=\beta I$ for some $\beta$, resulting in the update $\delta_{k}^{*}=-\frac{1}{\beta}\nabla h(\theta_{k})$. In the case where $h$ is convex and Lipschitz-smoothBy this we mean that $\|\nabla h(\theta)-\nabla h(\theta^{\prime})\|\leq\mathcal{L}\|\theta-\theta^{\prime}\|$ for all $\theta$ and $\theta^{\prime}$. with constant $\mathcal{L}$, a safe/conservative choice that will ensure convergence with $\alpha_{k}=1$ is $\beta=\mathcal{L}$ (e.g. Nesterov, 2013). The intuition behind this choice is that $B$ will act as a global upper bound on the curvature of $h$, in the sense that $B_{k}=\mathcal{L}I\succeq H(\theta)$Here we define $A\succeq C$ to mean that $A-C$ is PSD. for all $\theta$, so that $\delta_{k}^{*}$ never extends past the point that would be safe in the worst-case scenario where the curvature is at its upper bound $\mathcal{L}$ the entire way along $\delta^{*}$. More concretely, one can show that given this choice of $\beta$, $M_{k}(\delta)$ upper bounds $h(\theta_{k}+\delta)$, and will therefore never predict a reduction in $h(\theta_{k}+\delta)$ where there is actually a sharp increase (e.g. due to $h$ curving unexpectedly upward on the path from $\theta_{k}$ to $\theta_{k}+\delta$). Minimizing $M_{k}(\delta)$ is therefore guaranteed not to increase $h(\theta_{k}+\delta)$ beyond the current value $h(\theta_{k})$ since $M_{k}(0)=h(\theta_{k})$. But despite these nice properties, this choice will almost always overestimate the curvature in most directions, leading to updates that move unnecessarily slowly along directions of consistent low curvature.

While neural networks haven’t been closely studied by optimization researchers until somewhat recently, many of the local optimization issues related to neural network learning can be seen as special cases of problems which arise more generally in continuous optimization. For example, tightly coupled parameters with strong local dependencies, and large variations in scale along different directions in parameter space (which may arise due to the “vanishing gradient” phenomenon (Hochreiter et al., 2000)), are precisely the sorts of issues for which 2nd-order optimization is well suited. Gradient descent on the other hand is well known to be very sensitive to such issues, and in order to avoid large oscillations and instability must use a learning rate which is inversely proportional to $\mathcal{L}$. 2nd-order optimization methods provide a much more powerful and elegant solution to the problem of variations in scale/curvature along different directions by selectively re-scaling the gradient along different eigen-directions of the curvature matrix $B_{k}$ according to their associated curvature (eigenvalue), instead of employing a one-size-fits-all curvature estimate.

In the classical Newton’s method we take $B_{k}=H(\theta_{k})$, in which case $M_{k}(\delta)$ becomes the 2nd-order Taylor-series approximation of $h$ centered at $\theta_{k}$. This choice gives us the most accurate local model of the curvature possible, and allows for rapid exploration of low-curvature directions and thus faster convergence.

Unfortunately, naive implementations of Newton’s method can run into numerous problems when applied to neural network training objectives, such as $H$ being sometimes indefinite (and thus $M_{k}(\delta)$ being unbounded below in directions of negative curvature) and related issues of “model trust”, where the method implicitly trusts its own local quadratic model of the objective too much, causing it to propose very large updates that may actually increase the $h$. These problems are usually not encountered with first order methods, but only because they use a very conservative local model that is intrinsically incapable of generating large updates. Fortunately, using the Gauss-Newton approximation to the Hessian (as discussed in Section 8), and/or applying various update damping/trust-region techniques (as discussed in Section 10), the issue of model trust issue in 2nd-order methods can be mostly overcome.

Another important issue preventing the naive application of 2nd-order methods to neural networks is the typically very high dimensionality of the parameter space ($n$), which prohibits the calculation/storage/inversion of the $n^{2}$-entry curvature matrix $B_{k}$. To address this, various approximate Newton methods have been developed within the optimization and machine learning communities. These methods work by approximating $B_{k}$ with something easier to compute/store/invert such as a low-rank or diagonal matrix, or by performing only approximate/incomplete optimization of $M_{k}(\delta)$. A survey of such methods is outside the scope of this report, but many good references and reviews are available (e.g. Nocedal and Wright, 2006; Fletcher, 2013). Martens (2016) reviews these approaches specifically in the context of neural networks.

Finally, it is worth observing that the local optimality of the Hessian-based 2nd-order Taylor series approximation to $h$ won’t necessarily yield the fastest possible optimization procedure, as it is possible to imagine quadratic models that take a “longer view” of the objective. (As an extreme example, given knowledge of a global minimizer $\theta^{*}$ of $h$, one could construct a quadratic model whose minimizer is exactly $\theta^{*}$ but which is a very poor local approximation to $h(\theta)$.) It is possible that the Fisher might give rise to such quadratic models, which would help explain its observed superiority to the Hessian in neural network optimization (Schraudolph, 2002; Martens, 2010; Vinyals and Povey, 2012). We elaborate more on this speculative theory in Section 8.1.

The Generalized Gauss-Newton Matrix

This section discusses the Generalized Gauss-Newton matrix of Schraudolph (2002), and justifies its use as an alternative to the Hessian. Its relevance to our discussion of natural gradient methods will be made clear later in Section 9, where we establish a correspondence between this matrix and the Fisher.

The classical Gauss-Newton matrix (or more simply the Gauss-Newton matrix) is the curvature matrix $G$ which arises in the Gauss-Newton method for non-linear least squares problems (e.g. Dennis Jr and Schnabel, 1996; Ortega and Rheinboldt, 2000; Nocedal and Wright, 2006). It is applicable to our standard neural network training objective $h$ in the case where $L(y,z)=\frac{1}{2}\|y-z\|^{2}$, and is given by

where $J_{f}$ is the Jacobian of $f(x,\theta)$ w.r.t. the parameters $\theta$. It is usually defined as a modified version of the Hessian $H$ of $h$ (w.r.t. $\theta$), obtained by dropping the second term inside the sum in the following expression for $H$:

where $H_{[f]_{j}}$ is the Hessian (w.r.t. $\theta$) of the $j$-th component of $f(x,\theta)$. We can see from this expression that $G=H$ when $y=f(x,\theta)$. And more generally, if the $y$’s are well-described by the model $f(x,\theta)+\epsilon$ for i.i.d. noise $\epsilon$ then $G=H$ will hold approximately.

Beyond the fact that the resulting matrix is PSD and has other nice properties discussed below, there doesn’t seem to be any obvious justification for linearizing $f$ (or equivalently, dropping the corresponding term from the Hessian). It’s likely that the reasonableness of doing this depends on problem-specific details about $L$ and $f$, and how the curvature matrix will be used by the optimizer. In Subsection 8.1.3 we discuss how linearizing $f$ might be justified, specifically for wide neural networks, by some recent theoretical analyses.

Schraudolph (2002) showed how the idea of the Gauss-Newton matrix can be generalized to the situation where $L(y,z)$ is any loss function which is convex in $z$. The generalized formula for $G$ is

where $H_{L}$ is the Hessian of $L(y,z)$ w.r.t. $z$, evaluated at $z=f(x,\theta)$. Because $L(y,z)$ is convex in $z$, $H_{L}$ will be PSD for each $(x,y)$, and thus so will $G$. We will call this $G$ the Generalized Gauss-Newton matrix (GGN). (Note that this definition is sensitive to where we draw the dividing line between the loss function $L$ and the network itself (i.e. $z$), in contrast to the definition of the Fisher, which is invariant to this choice.)

Analogously to the case of the classical Gauss-Newton matrix (which assumed $L(y,z)=\frac{1}{2}\|y-z\|^{2}$), the GGN can be obtained by dropping the second term inside the sum of the following expression for the Hessian $H$:

Here $\left.\nabla_{z}L(y,z)\right|_{z=f(x,\theta)}$ is the gradient of $L(y,z)$ w.r.t. $z$, evaluated at $z=f(x,\theta)$. Note if we have for some local optimum $\theta^{*}$ that $\left[\left.\nabla_{z}L(y,z)\right|_{z=f(x,\theta^{*})}\right]_{j}\approx 0$ for each $(x,y)$ and $j$, which corresponds to the network making an optimal prediction for each training case over each dimension, then $G(\theta^{*})=H(\theta^{*})$. In such a case, the behavior of a 2nd-order optimizer using $G$ will approach the behavior of standard Newton’s method as it converges to $\theta^{*}$. A weaker condition implying equivalence is that $\frac{1}{S_{y}(x)}\sum_{y\in S_{y}(x)}\left[\left.\nabla_{z}L(y,z)\right|_{z=f(x,\theta^{*})}\right]_{j}\approx 0$ for all $x\in S_{x}$ and $j$, where $S_{y}(x)$ denotes the set of $y$’s s.t. $(x,y)\in S$, which corresponds to the network making an optimal prediction for each $x$ in the presence of intrinsic uncertainty about the target $y$. (This can be seen by noting that $H_{[f]_{j}}$ doesn’t depend on $y$.)

Like the Hessian, the GGN can be used to define a local quadratic model of $h$, as given by:

In 2nd-order methods based on the GGN, parameter updates are computed by minimizing $M_{k}(\delta)$ w.r.t. $\delta$. The exact minimizerThis formula assumes $G(\theta_{k})$ is invertible. If it’s not, the problem will either be unbounded, or the solution can be computed using the pseudo-inverse instead. $\delta^{*}=-G(\theta_{k})^{-1}\nabla h(\theta_{k})$ is often too difficult to compute, and so practical methods will often only approximately minimize $M_{k}(\delta)$ (e.g Dembo et al., 1982; Steihaug, 1983; Dennis Jr and Schnabel, 1996; Martens, 2010; Vinyals and Povey, 2012).

Since computing the whole matrix explicitly is usually too expensive, the GGN is typically accessed via matrix-vector products. To compute such products efficiently one can use the method of Schraudolph (2002), which is a generalization of the well-known method for computing such products with the classical Gauss-Newton (and is also related to the TangentProp method of Simard et al. (1992)). The method is similar in cost and structure to standard backpropagation, although it can sometimes be tricky to implement (see Martens and Sutskever (2012)).

Unlike the Hessian, the GGN is positive semi-definite (PSD). This means that it never models the curvature as negative in any direction. The most obvious problem with negative curvature is that the quadratic model will predict an unbounded improvement in the objective for moving in the associated directions. Indeed, without the use of some kind of trust-region or damping technique (as discussed in Section 10), or pruning/modification of negative curvature directions (Vinyals and Povey, 2012; Dauphin et al., 2014), or self-terminating Newton-CG scheme (Steihaug, 1983), the update produced by minimizing the quadratic model will be infinitely large in such directions.

However, attempts to use such methods in combination with the Hessian have yielded lackluster results for neural network optimization compared to methods based on the GGN (Martens, 2010; Martens and Sutskever, 2012; Vinyals and Povey, 2012). So what might be going on here? While the true curvature of $h(\theta)$ can indeed be negative in a local neighborhood (as measured by the Hessian), we know it must quickly become non-negative as we travel along any particular direction, given that our loss $L(y,z)$ is convex in $z$ and bounded below. Meanwhile, positive curvature predicts a quadratic penalty, and in the worst case merely underestimates how badly the objective will eventually increase along a particular direction. We can thus say that negative curvature is somewhat less “trustworthy” than positive curvature for this reason, and speculate that a 2nd-order method based on the GGN won’t have to rely as much on trust-regions etc (which restrict the size of the update and slow down performance) to produce reliable updates.

There is also the issue of estimation from limited data. Because contributions made to the GGN for each training case and each individual component of $f(x,\theta)$ are PSD, there can be no cancellation between positive and negative/indefinite contributions. This means that the GGN can be more robustly estimated from subsets of the training data than the Hessian. (By analogy, consider how much harder it is to estimate the scale of the mean value of a variable when that variable can take on both positive and negative values, and has a mean close to .) This property also means that positive curvature from one case or component will never be cancelled out by negative curvature from another case or component. And if we believe that negative curvature is less trustworthy than positive curvature over larger distances, this is probably a good thing.

Despite these nice properties, the GGN is notably not an upper bound on the Hessian (in the PSD sense), as it fails to model all of the positive curvature contained in the latter. But crucially, it only fails to model the (positive or negative) curvature coming from the network function $f(x,\theta)$, as opposed to the curvature coming from the loss function $L(y,z)$. (To see this, recall the decomposition of the Hessian from eqn. 6, noting that the term dropped from the Hessian depends only on the gradients of $L$ and the Hessian of components of $f$.) Curvature coming from $f$, whether it is positive or negative, is arguably less trustworthy/stable across long distance than curvature coming from $L$, as argued below.

1.2 A More Detailed View of the Hessian vs the GGN

Consider the following decomposition of the Hessian, which is a generalization of the one given in eqn. 6:

Here, $a_{i}$, $\phi_{i}$ and $s_{i}$ are defined as in Section 2, $\nabla_{a_{i}}L(y,f)$ is the gradient of $L(y,f)$ w.r.t. $a_{i}$, $H_{[\phi_{i}(s_{i})]_{j}}$ is the Hessian of $\phi_{i}(s_{i})$ (i.e. the function which computes $a_{i}$) w.r.t. $s_{i}$, $J_{s_{i}}$ is the Jacobian of $s_{i}$ (viewed as a function of $\theta$ and $x$) w.r.t. $\theta$, and $C$ is given by

where $\otimes$ denotes the Kronecker product.

Aside from $C+C^{\top}$, the contribution to the Hessian made by $f$ comes from a sum of terms of the form $\left[\nabla_{a_{i}}L(y,f)\right]_{j}J_{s_{i}}^{\top}H_{[\phi_{i}(s_{i})]_{j}}J_{s_{i}}$. These terms represent the curvature of the activation functions, and will be zero in a linear network (since we would have $H_{[\phi_{i}(s_{i})]_{j}}=0$). It seems reasonable to suspect that the sign of these terms will be subject to rapid and unpredictable change, resulting from sign changes in both $H_{[\phi_{i}(s_{i})]_{j}}$ and $\left[\nabla_{a_{i}}L(y,f)\right]_{j}$. The former is the “local Hessian” of $\phi_{i}$, and will change signs as the function $\phi_{i}$ enters its different convex and concave regions ($\phi_{i}$ is typically non-convex). $\left[\nabla_{a_{i}}L(y,f)\right]_{j}$ meanwhile is the loss derivative w.r.t. that unit’s output, and depends on the behavior of all of the layers above $a_{i}$, and on which “side” of the training target the network’s current prediction is (which may flip back and forth at each iteration).

This is to be contrasted with the term $J_{f}^{\top}H_{L}J_{f}$, which represents the curvature of the loss function $L$, and which remains PSD everywhere (and for each individual training case). Arguably, this term will be more stable w.r.t. changes in the parameters, especially when averaged over the training set.

1.3 Some Insights From Other Works

In the case of the squared error loss $L(y,z)=\frac{1}{2}\|y-z\|^{2}$ (which means that the GGN reduces to the standard Gauss-Newton matrix) with $m=1$, Chen (2011) established that the GGN is the unique matrix which gives rise to a local quadratic approximation of $L(y,f(x,\theta))$ which is both non-negative (as $L$ itself is), and vanishing on a subspace of dimension $n-1$ (which is the dimension of the manifold on which $L$ itself vanishes). Notably, the quadratic approximation produced using the Hessian need not have either of these properties. By summing over output components and averaging over the training set $S$, one should be able to generalize this result to the entire objective $h(\theta)$ with $m\geq 1$. Thus, we see that the GGN gives rise to a quadratic approximation which shares certain global characteristics with the true $h(\theta)$ that the 2nd-order Taylor series doesn’t, despite being a less precise approximation to $h(\theta)$ in a strictly local sense.

Botev et al. (2017) observed that for networks with piece-wise linear activation functions, such as the popular RELUs (given by $[\phi_{i}(s_{i})]_{j}=\max([s_{i}]_{j},0)$), the GGN and the Hessian will coincide on the diagonal blocks whenever the latter is well-defined. This can be seen from the above decomposition of $H$ by noting that $C+C^{\top}$ is zero on the diagonal blocks, and that for piece-wise linear activation functions we have $H_{[\phi_{i}(s_{i})]_{j}}=0$ everywhere that this quantity exists (i.e. everywhere except the “kinks” in the activation functions).

Finally, under certain realistic assumptions on the network architecture and initialization point, and a lower bound on the width of the layers, recent results have shown that the $f$ function for a neural network behaves very similarly to its local linear approximation (taken at the initial parameters) throughout the entirety of optimization. This happens both for gradient descent (Du et al., 2018; Jacot et al., 2018; Lee et al., 2019), and natural gradient descent / GGN-based methods (Zhang et al., 2019; Cai et al., 2019), applied to certain choices for $L$. Not only does this allow one to prove strong global convergence guarantees for these algorithms, it lends support to the idea that modeling the curvature in $f$ (which is precisely the part of the Hessian that the GGN throws out) may be pointless for the purposes of optimization in neural networks, and perhaps even counter-productive.

Computational Aspects of the Natural Gradient and Connections to the Generalized Gauss-Newton Matrix

where $J_{f}$ is the Jacobian of $f(x,\theta)$ w.r.t. $\theta$, and $\nabla_{z}\log r(y|z)$ is the gradient of $\log r(y|z)$ w.r.t. $z$, evaluated at $z=f(x,\theta)$ (with $r$ defined as near the end of Section 4).

As was first shown by Park et al. (2000), the Fisher information matrix is thus given by

where $F_{R}$ is the Fisher information matrix of the predictive distribution $R_{y|z}$ at $z=f(x,\theta)$. $F_{R}$ is itself given by

where $H_{\log r(y|z)}$ is the Hessian of $\log r(y|z)$ w.r.t. $z$, evaluated at $z=f(x,\theta)$.

Note that even if $Q_{x}$’s density function $q(x)$ is known, and is relatively simple, only for certain choices of $R_{y|z}$ and $f(x,\theta$) will it be possible to analytically evaluate the expectation w.r.t. $Q_{x}$ in the above expression for $F$. For example, if we take $Q_{x}=\mathcal{N}(0,I)$, $R_{y|z}=\mathcal{N}(z,\sigma^{2})$, and $f$ to be a simple neural network with no hidden units and a single tan-sigmoid output unit, then both $F$ and its inverse can be computed efficiently (Amari, 1998). This situation is exceptional however, and for even slightly more complex models, such as neural networks with one or more hidden layers, it has never been demonstrated how to make such computations feasible in high dimensions.

Fortunately the situation improves significantly if $Q_{x}$ is replaced by $\hat{Q}_{x}$, as this gives

which is easy to compute assuming $F_{R}$ is easy to compute. Moreover, this is essentially equivalent to the expression in eqn. 5 for the generalized Gauss-Newton matrix (GGN), except that we have the Fisher $F_{R}$ of the predictive distribution ($R_{y|z}$) instead of Hessian $H_{L}$ of the loss ($L$) as the “inner” matrix.

Eqn. 7 also suggests a straightforward and efficient way of computing matrix-vector products with $F$, using an approach similar to the one in Schraudolph (2002) for computing matrix-vector products with the GGN. In particular, one can multiply by $J_{f}$ using a linearized forward pass (aka forward-mode automatic differentiation), then multiply by $F_{R}$ (which will be easy if $R_{y|z}$ is sufficiently simple), and then finally multiply by $J_{f}^{\top}$ using standard backprop.

2 Qualified Equivalence of the GNN and the Fisher

As we shall see in this subsection, the connections between the GGN and Fisher run deeper than just similar expressions and algorithms for computing matrix-vector products.

Heskes (2000) showed that the Fisher and the classical Gauss-Newton matrix are equivalent in the case of the squared error loss, and proposed using the Fisher as an alternative to the Hessian in more general contexts. Concurrently with this work, Pascanu and Bengio (2014) showed that for several common loss functions like cross-entropy and squared error, the GGN and Fisher are equivalent.

We will show that in fact there is a much more general equivalence between the two matrices, starting from the observation that the expressions for the GGN in eqn. 5 and Fisher in eqn. 7 are identical up to the equivalence of $H_{L}$ and $F_{R}$.

First, note that $L(y,z)$ might not even be convex in $z$, so that it wouldn’t define a valid GGN matrix. But even if $L(y,z)$ is convex in $z$, it won’t be true in general that $F_{R}=H_{L}$, and so the GGN and Fisher will differ. However, there is an important class of $R_{y|z}$’s for which $F_{R}=H_{L}$ will hold, provided that we have $L(y,z)=-\log r(y|z)$ (putting us in the framework of Section 4).

Notice that $F_{R}=-\operatorname{E}_{R_{y|f(x,\theta)}}[H_{\log r(y|z)}]$, and $H_{L}=-H_{\log r(y|z)}$ (which follows from $L(y,z)=-\log r(y|z)$). Thus, the two matrices being equal is equivalent to the condition

While this condition may seem arbitrary, it is actually very natural and holds in the important case where $R_{y|z}$ corresponds to an exponential family model with “natural” parameters given by $z$. Stated in terms of equations this condition is

for some function $T(y)$, where $Z(z)$ is the normalizing constant/partition function. In this case we have $H_{\log r(y|z)}=-H_{\log Z}$ (which doesn’t depend on $y$), and so eqn. 8 holds trivially.

multivariate normal distributions where $z$ parameterizes only the mean $\mu$

multivariate normal distributions where $z$ is the concatenation of $\Sigma^{-1}\mu$ and the vectorization of $\Sigma^{-1}$

multinomial distributions where the softmax of $z$ is the vector of probabilities for each class

Note that the loss function $L$ corresponding to the multivariate normal is the familiar squared error, and the loss corresponding to the multinomial distribution is the familiar cross-entropy.

Interestingly, the relationship observed by Ollivier et al. (2018) between natural gradient descent and methods based on the extended Kalman filter for neural network training relies on precisely the same condition on $R_{y|z}$. This makes intuitive sense, since the extended Kalman filter is derived by approximating $f$ as affine and then applying the standard Kalman filter for linear/Gaussian systems (which implicitly involves computing a Hessian of a linear model under a squared loss), which is the same approximation that can be used to derive the GGN from the Hessian (see Section 8).

As discussed in Section 8, when constructing the GGN one must pay attention to how $f$ and $L$ are defined with regards to what parts of the neural network’s computation are performed by each function (this choice is irrelevant to the Fisher). For example, the softmax computation performed at the final layer of a classification network is usually considered to be part of the network itself and hence to be part of $f$. The output $f(x,\theta)$ of this computation are normalized probabilities, which are then fed into a cross-entropy loss of the form $L(y,z)=-\sum_{j}y_{j}\log z_{j}$. But the other way of doing it, which Schraudolph (2002) recommends, is to have the softmax function be part of $L$ instead of $f$, which results in a GGN which is slightly closer to the Hessian due to “less” of the computational pipeline being linearized before taking the 2nd-order Taylor series approximation. The corresponding loss function is $L(y,z)=-\sum_{j}y_{j}z_{j}+\log(\sum_{j}\exp(z_{j}))$ in this case. As we have established above, doing it this way also has the nice side effect of making the GGN equivalent to the Fisher, provided that $R_{y|z}$ is an exponential family model with $z$ as its natural parameters.

This (qualified) equivalence between the Fisher and the GGN suggests how the GGN can be generalized to cases where it might not otherwise be well-defined. In particular, it suggests formulating the loss as the negative log density for some distribution, and then taking the Fisher of this distribution. Sometimes, this might be as simple as defining $r(y|z)\propto\exp(-L(y,z))$ as per the discussion at the end of Section 4.

For example, suppose our loss is defined as the negative log probability of a multi-variate normal distribution $R_{y|z}=N(\mu,\sigma^{2})$ parameterized by $\mu$ and $\gamma=\log\sigma^{2}$ (so that $z=\begin{bmatrix}\mu\\ \gamma\end{bmatrix}$). In other words, suppose that

In this case the loss Hessian is equal to

It is not hard to verify that this matrix is indefinite for certain settings of $x$ and $z$ (e.g. $x=2$, $\mu=\gamma=0$). Therefore, $L$ is not convex in $z$ and we cannot define a valid GGN matrix from it.

To resolve this problem we can use the Fisher $F_{R}$ in place of $H_{L}$ in the formula for the GGN, which by eqn. 7 yields $F$. Alternatively, we can insert reparameterization operations into our network to transform $\mu$ and $\gamma$ into the natural parameters $\frac{\mu}{\sigma^{2}}=\frac{\mu}{\exp(\gamma)}$ and $-\frac{1}{2\sigma^{2}}=-\frac{1}{2\exp(\gamma)}$, and then proceed to compute the GGN as usual, noting that $H_{L}=F_{R}$ in this case, so that $H_{L}$ will be PSD. Either way will yield the same curvature matrix, due to the above discussed equivalence of the Fisher and GGN matrix for natural parameterizations.

Constructing Practical Natural Gradient Methods, and the Critical Role of Damping

Assuming that it is easy to compute, the simplest way to use the natural gradient in optimization is to substitute it in place of the standard gradient within a basic gradient descent approach. This gives the iteration

where $\{\alpha_{k}\}_{k}$ is a schedule of step-sizes/learning-rates.

Choosing the step-size schedule can be difficult. There are adaptive schemes which are largely heuristic in nature (Amari, 1998) and some non-adaptive prescriptions such as $\alpha_{k}=\rho/k$ for some constant $\rho$, which have certain theoretical convergence guarantees in the stochastic setting, but which won’t necessarily work well in practice.

In principle, we could apply the natural gradient method with infinitesimally small steps and produce a smooth idealized path through the space of realizable distributions. But since this is usually impossible in practice, and we don’t have access to any other simple description of the class of distributions parameterized by $\theta$ that we could work with more directly, our only option is to take non-negligible discrete steps in the given parameter spaceIn principle, we could move to a much more general class of distributions, such as those given by some non-parametric formulation, where we could work directly with the distributions themselves. But even assuming such an approach would be practical from a computational efficiency standpoint, we would lose the various advantages that we get from working with powerful parametric models like neural networks. In particular, we would lose their ability to generalize to unseen data by modeling the “computational process” which explains the data, instead of merely using smoothness and locality to generalize..

The fundamental problem with simple schemes such as the one in eqn. 9 is that they implicitly assume that the natural gradient is a good direction to follow over non-negligible distances in the original parameter space, which will not be true in general. Traveling along a straight line in the original parameter space will not yield a straight line in distribution space, and so the resulting path may instead veer far away from the target that the natural gradient originally pointed towards. This is illustrated in Figure 1.

Fortunately, we can exploit the (qualified) equivalence between the Fisher and the GGN in order to produce natural gradient-like updates which will often be appropriate to take with $\alpha_{k}=1$. In particular, we know from the discussion in Section 8 that the GGN matrix $G$ can serve as a reasonable proxy for the Hessian $H$ of $h$, and may even be superior in certain contexts. Meanwhile, the update $\delta$ produced by minimizing the GGN-based local quadratic model $M_{k}(\delta)=\frac{1}{2}\delta^{\top}G(\theta_{k})\delta+\nabla h(\theta_{k})^{\top}\delta+h(\theta_{k})$ is given by $-G(\theta_{k})^{-1}\nabla h(\theta_{k})$, which will be equal to the negative natural gradient when $F=G$. Thus, the (negative) natural gradient, with scaling factor $\alpha=1$, can be seen as the optimal update according to a particular local 2nd-order approximation of $h$. And just as in the case of other 2nd-order methods, the break-down in the accuracy of this quadratic approximation over long distances, combined with the potential for the natural gradient to be very large (e.g. when $F$ contains some very small eigenvalues), can often lead to very large and very poor update proposals. Simply re-scaling the update by reducing $\alpha$ may be too crude a mechanism to deal with this subtle problem, as it will affect all eigen-directions (of $F$) equally, including those in which the natural gradient is already sensible, or even overly conservative.

Instead, the connection between natural gradient descent and 2nd-order methods motivates the use of “update damping” techniques that have been developed for the latter, which work by constraining or penalizing the solution for $\delta$ in various ways during the optimization of $M_{k}(\delta)$. Examples include Tikhonov regularization/damping and the closely related trust-region method (e.g. Tikhonov, 1943; Moré and Sorensen, 1983; Conn et al., 2000; Nocedal and Wright, 2006), and other ones such as the “structural damping” approach of Martens and Sutskever (2011), or the approach present in Krylov Subspace Descent (Vinyals and Povey, 2012). See Martens and Sutskever (2012) for an in-depth discussion of these and other damping techniques in the context of neural network optimization.

This idea is supported by practical experience in neural network optimization. For example, the Hessian-free optimization approach of Martens (2010) generates its updates using a Tikhonov damping scheme applied to the exact GGN matrix (which was equivalent to the Fisher in that work). These updates, which can be applied with a step-size of 1, make a lot more progress optimizing the objective than updates computed without any damping (which must instead rely on a carefully chosen step-size to even be feasible).

It is worth pointing out that other interpretations of natural gradient descent can also motivate the use of damping/regularization terms. In particular, Ollivier et al. (2018) has shown that online natural gradient descent, with a particular flavor of Tikhonov regularization, closely resembles a certain type of extended Kalman filter-based training algorithm for neural networks (Singhal and Wu, 1989; Ruck et al., 1992), where $\theta$ is treated as an evolving hidden state that is estimated by the filter (using training targets as noisy observations and inputs as control signals).

The Empirical Fisher

An approximation of the Fisher known as the “empirical Fisher” (Schraudolph, 2002), which we denote by $\bar{F}$, is commonly used in practical natural gradient methods. It is obtained by taking the inner expectation of eqn. 3 over the target distribution $Q_{x,y}$ (or its empirical surrogate $\hat{Q}_{x,y}$) instead of the model’s distribution $P_{x,y}$.

In the case where one uses $\hat{Q}_{x,y}$, this yields the following simple form:

This matrix is often incorrectly referred to as the Fisher, or even the Gauss-Newton, even though it is not equivalent to either of these matrices in general.

Like the Fisher $F$, the empirical Fisher $\bar{F}$ is PSD. But unlike $F$, it is essentially free to compute, provided that one is already computing the gradient of $h$. And it can also be applied to objective functions which might not involve a probabilistic model in any obvious way.

Compared to $F$, which is of rank $\leq|S|\operatorname{rank}(F_{R})$, $\bar{F}$ has a rank of $\leq|S|$, which can make it easier to work with in practice. For example, the problem of computing the diagonal (or various blocks) is easier for the empirical Fisher than it is for higher rank matrices like the standard Fisher (Martens et al., 2012). This has motivated its use in optimization methods such as TONGA (Le Roux et al., 2008), and as the diagonal preconditioner of choice in the Hessian-free optimization method (Martens, 2010). Interestingly however, there are stochastic estimation methods (Chapelle and Erhan, 2011; Martens et al., 2012) which can be used to efficiently estimate the diagonal (or various blocks) of the standard Fisher $F$, and these work quite well in practice. (These include the obvious method of sampling $y$’s from the model’s conditional distribution and computing gradients from them, but also includes methods based on matrix factorization and random signs. See Martens et al. (2012) for comparative analysis of the variance of these methods.)

Despite the various practical advantages of using $\bar{F}$, there are good reasons to use true Fisher $F$ instead of $\bar{F}$ whenever possible. In addition to Amari’s extensive theory developed for the exact natural gradient (which uses $F$), perhaps the best reason for using $F$ over $\bar{F}$ is that $F$ turns out to be a reasonable approximation/substitute to the Hessian $H$ of $h$ in certain important special cases, which is a property that $\bar{F}$ lacks in general.

For example, as discussed in Section 5, when the loss is given by $-\log p(y|x)$ (as in Section 4), $F$ can be seen as an approximation of $H$, because both matrices have the interpretation of being the expected Hessian of the loss under some distribution. Due to the similarity of the expression for $F$ in eqn. 3 and the one above for $\bar{F}$, it might be tempting to think that $\bar{F}$ is given by the expected Hessian of the loss under $\hat{Q}_{x,y}$ (which is actually the formula for $H$) in the same way that $F$ is given by eqn. 4. But this is not the case in general.

And as we saw in Section 9, given certain assumptions about how the GGN is computed, and some additional assumptions about the form of the loss function $L$, $F$ turns out to be equivalent to the GGN. This is very useful since the GGN can be used to define a local quadratic approximation of $h$, whereas $F$ normally doesn’t have such an interpretation. Moreover, Schraudolph (2002) and later Martens (2010) compared $\bar{F}$ to the GGN and observed that the latter performed much better as a curvature matrix within various neural network optimization methods.

As concrete evidence for why the empirical Fisher is, at best, a questionable choice for the curvature matrix, we will consider the following example. Set $n=1$, $f(x,\theta)=\theta$, $R_{y|z}=\mathcal{N}(z,1)$, and $S=\{(0,0)\}$, so that $h(\theta)$ is a simple convex quadratic function of $\theta$, given by $h(\theta)=\frac{1}{2}\theta^{2}$. In this example we have that $\nabla h=\theta$, $\bar{F}=\theta^{2}$, while $F=1$. If we use $\bar{F}^{\xi}$ as our curvature matrix for some exponent $\frac{1}{2}\leq\xi\leq 1$, then it is easy to see that an iteration of the form

will fail to converge to the minimizer (at $\theta=0$) unless $\xi<1$ and the step-size $\alpha_{k}$ goes to sufficiently fast. And even when it does converge, it will only be at a rate comparable to the speed at which $\alpha_{k}$ goes to , which in typical situations will be either $\mathcal{O}(1/k)$ or $\mathcal{O}(1/\sqrt{k})$. Meanwhile, a similar iteration of the form

which uses the exact Fisher $F$ as the curvature matrix, will experience very fast linear convergenceHere we mean “linear” in the classical sense that $|\theta_{k}-0|\leq|\theta_{0}-0||1-\alpha|^{k}$. with rate $|1-\alpha|$, for any fixed step-size $\alpha_{k}=\alpha$ satisfying $0<\alpha<2$.

It is important to note that this example uses a noise-free version of the gradient, and that this kind of linear convergence is (provably) impossible in most realistic stochastic/online settings. Nevertheless, we would argue that a highly desirable property of any stochastic optimization method should be that it can, in principle, revert to an optimal (or nearly optimal) behavior in the deterministic setting. This might matter a lot in practice, since the gradient may end up being sufficiently well estimated in earlier stages of optimization from only a small amount of data (which is a common occurrence in our experience), or in later stages provided that larger mini-batches or other variance-reducing procedures are employed (e.g. Le Roux et al., 2012; Johnson and Zhang, 2013). More concretely, the pre-asymptotic convergence rate of stochastic 2nd-order optimizers can still depend strongly on the choice of the curvature matrix, as we will show in Section 14.

2 A Discussion of Recent Diagonal Methods Based on the Empirical Fisher

Recently, a spate of stochastic optimization methods have been proposed that are all based on diagonal approximations of the empirical Fisher $\bar{F}$. These include the diagonal version of AdaGrad (Duchi et al., 2011), RMSProp (Tieleman and Hinton, 2012), Adam (Ba and Kingma, 2015), etc. Such methods use iterations of the following form (possibly with some slight modifications):

where the curvature matrix $B_{k}$ is taken to be a diagonal matrix $\operatorname{diag}(u_{k})$ with $u_{k}$ adapted to maintain some kind of estimate of the diagonal of $\bar{F}$ (possibly using information from previous iterates/mini-batches), $g_{k}(\theta_{k})$ is an estimate of $\nabla h(\theta_{k})$ produced from the current mini-batch, ${\alpha_{k}}_{k}$ is a schedule of step-sizes, and $0<\lambda$ and $0<\xi\leq 1$ are hyperparameters (discussed later in this section).

There are also slightly more sophisticated methods (Schaul et al., 2013; Zeiler, 2013) which use preconditioners that combine the diagonal of $\bar{F}$ with other quantities (such as an approximation of the diagonal of the Gauss-Newton/Fisher in the case of Schaul et al. (2013)) in order to correct for how the empirical Fisher doesn’t have the right “scale” (which is ultimately the reason why it does poorly in the example given at the end of Section 11.1).

A diagonal preconditioner (Nash, 1985) of the form used in eqn. 10 was also used by (Martens, 2010) to accelerate the conjugate gradient (CG) sub-optimizations performed within a truncated-Newton method (using the GGN matrix). In the context of CG, the improper scale of $\bar{F}$ is not as serious an issue due to the fact that CG is invariant to the overall scale of its preconditioner (since it computes an optimal “step-size” at each step which automatically adjusts for the scale). However, it still makes more sense to use the diagonal of the true Fisher $F$ as a preconditioner, and thanks to the method proposed by Chapelle and Erhan (2011), this can be estimated efficiently and accurately.

The idea of using the diagonal of $F$, $\bar{F}$, or the Gauss-Newton as a preconditioner for stochastic gradient descent (SGD) and was likely first applied to neural networks with the work of Lecun and collaborators (Becker and LeCun, 1989; LeCun et al., 1998), who proposed an iteration of the form in eqn. 10 with $\xi=1$ where $u_{k}$ approximates the diagonal of the Hessian or the Gauss-Newton matrix (which as shown in Section 9, is actually equivalent to $F$ for the common squared-error loss). Following this work, various neural network optimization methods have been developed over the last couple of decades that use diagonal, block-diagonal, low-rank, or Krylov-subspace based approximations of $F$ or $\bar{F}$ as a curvature matrix/preconditioner. In addition to methods based on diagonal approximations already mentioned, some methods based on non-diagonal approximations include the method of Park et al. (2000), TONGA (Le Roux et al., 2008), Natural Newton (Le Roux and Fitzgibbon, 2010), HF (Martens, 2010), KSD (Vinyals and Povey, 2012) and many more.

The idea of computing an estimate of the (empirical) Fisher using a history of previous iterates/mini-batches also appeared in various early works. The particular way of doing this proposed Duchi et al. (2011), which is to use an equally weighted average of all past gradients, was motivated from a regret-based asymptotic convergence analysis and tends not to work well in practice (Tieleman and Hinton, 2012). The traditional and more intuitive approach of using an exponentially decayed running average (e.g. LeCun et al., 1998; Park et al., 2000) works better, at least pre-asymptotically, as it is able to naturally “forget” very old contributions to the estimate (which are based on stale parameter values).

It is important to observe that the way $\bar{F}$ is estimated can affect the convergence characteristics of an iteration like eqn. 10 in subtle but important ways. For example, if $\bar{F}$ is estimated using gradients from previous iterations, and especially if it is the average of all past gradients (as in AdaGrad), it may shrink sufficiently slowly that the convergence issues seen in the example at the end of Section 11.1 are avoided. Moreover, for reasons related to this phenomenon, it seems likely that the proofs of regret bounds in Duchi et al. (2011) and the related work of Hazan et al. (2007) could not be modified to work if the exact $\bar{F}$, computed only at the current $\theta$, were used. Developing a better understanding of this issue, and the relationship between methods developed in the online learning literature (such as AdaGrad), and classical stochastic 2nd-order methods based on notions of curvature, remains an interesting direction for future research.

3 The Constants λ𝜆\lambda and ξ𝜉\xi

The constants $\lambda$ and $\xi$ present in eqn. 10 are often thought of as fudge factors designed to correct for the “poor conditioning” (Becker and LeCun, 1989) of the curvature matrix, or to guarantee boundedness of the updates and prevent the optimizer from “blowing up” (LeCun et al., 1998). However, these explanations are oversimplifications that reference the symptoms instead of the cause. A more compelling and functional explanation, at least in the case of $\lambda$, comes from viewing the update in eqn. 10 as being the minimizer of a local quadratic approximation $M_{k}(\delta)=\frac{1}{2}\delta^{\top}B_{k}\delta+\nabla h(\theta_{k})^{\top}\delta+h(\theta_{k})$ to $h(\theta_{k}+\delta)$, as discussed in Section 10. In this view, $\lambda$ plays the role of a Tikhonov damping parameter (Tikhonov, 1943; Conn et al., 2000; Nocedal and Wright, 2006; Martens and Sutskever, 2012) which is added to $B_{k}$ in order to ensure that the proposed update stays within a certain region around zero in which $M_{k}(\delta)$ remains a reasonable approximation to $h(\theta_{k}+\delta)$. Note that this explanation implies that no single fixed value of $\lambda$ will be appropriate throughout the entire course of optimization, since the local properties of the objective will change, and so an adaptive adjustment scheme, such as the one present in HF (Martens, 2010) (which is based on the Levenberg-Marquardt method), should be used.

The use of the exponent $\xi=3/4$ first appeared in HF as part of its diagonal preconditioner for CG, and was justified as a way of making the curvature estimate “more conservative” by making it closer to a multiple of the identity, to compensate for the diagonal approximation being made (among other things). Around the same time, Duchi et al. (2011) proposed to use $\xi=1/2$ within an update of the form of eqn. 10, which was required in order to prove certain regret bounds for non-strongly-convex objectives.

However, it is important to note that Hazan et al. (2007) also proves a $\mathcal{O}(\log(k))$ bound on the regret for a basic version of SGD, and that what actually differentiates the various methods they analyze is the constant hidden in the big-O notation, which is much larger for the version of SGD they consider than for their approximate Newton method. In particular, the former depends on a quantity which grows with the condition number of the Hessian $H$ at $\theta^{*}$ while the latter does not, in a way that echos the various analyses performed on stochastic gradient descent and stochastic approximations of Newton’s method in the more classical “local-convergence” setting (e.g. Murata, 1998; Bottou and LeCun, 2005).

A Critical Analysis of Parameterization Invariance

One of the main selling points of the natural gradient method is its invariance to reparameterizations of the model. In particular, the smooth path through the space of distributions generated by the idealized natural gradient method with infinitesimally small steps will be invariant to any smooth invertible reparameterization of the $f$.

In this section we will examine this “smooth path parameterization invariance” property more closely in order to answer the following questions:

How can we characterize it using only basic properties of the curvature matrix?

Is there an elementary proof that can be applied in a variety of settings?

What other kinds of curvature matrices give rise to it, and is the Hessian included among these?

Will this invariance property imply that practical optimization algorithms based on the natural gradient (i.e. those that use large steps) will behave in a way that is invariant to the parameterization?

Let $\zeta$ be as above, and let $d_{\theta}$ and $d_{\gamma}$ be updates given in $\theta$-space and $\gamma$-space (resp.). Additively updating $\gamma$ by $d_{\gamma}$ and translating it back to $\theta$-space via $\zeta$ gives $\zeta(\gamma+d_{\gamma})$. Measured by some non-specific norm $\|\cdot\|$, this differs from $\theta+d_{\theta}$ by:

where $J_{\zeta}$ is the Jacobian of $\zeta$, and we have used $\theta=\zeta(\gamma)$.

The first term on the RHS of eqn. 11 measures the extent to which $\zeta(\gamma+d_{\gamma})$ fails to be predicted by the first-order Taylor series approximation of $\zeta$ centered at $\gamma$ (i.e. the local affine approximation of $\zeta$ at $\gamma$). This quantity will depend on the size of $d_{\gamma}$, and the amount of curvature in $\gamma$. In the case where $\zeta$ is affine, it will be exactly . We can further bound it by applying Taylor’s theorem for each component of $\zeta$, which gives

for some $c_{i}\in(0,1)$. If we assume that there is some $C>0$ so that for all $i$ and $\gamma$, $\|H_{[\zeta]_{i}}(\gamma)\|_{2}\leq C$, then using the fact that $|d_{\gamma}^{\top}H_{[\zeta]_{i}}(\gamma+c_{n}d_{\gamma})d_{\gamma}|\leq\frac{1}{2}\|H_{[\zeta]_{i}}(\gamma+c_{i}d_{\gamma})\|_{2}\|d_{\gamma}\|^{2}$, we can further upper bound this by $\frac{1}{2}C\sqrt{n}\|d_{\gamma}\|^{2}$.

The second term on the RHS of eqn. 11 will be zero when

which (as we will see) is a condition that is satisfied in certain natural situations. A slightly weakened version of this condition is that $J_{\zeta}d_{\gamma}\propto d_{\theta}$. Because we have

this condition can thus be interpreted as saying that $d_{\gamma}$, when translated appropriately via $\zeta$, points in the same direction away from $\theta$ that $d_{\theta}$ does. In the smooth path case, where the optimizer only moves an infinitesimally small distance in the direction of $d_{\gamma}$ (or $d_{\theta}$) at each iteration before recomputing it at the new $\gamma$ (or $\theta$), this condition is sufficient to establish that the path in $\gamma$ space, when mapped back to $\theta$ space via the $\zeta$ function, will be the same as the path which would have been taken if the optimizer had worked directly in $\theta$ space.

However, for a practical update scheme where we move the entire distance of $d_{\gamma}$ or $d_{\theta}$ before recomputing the update vector, such as the one in eqn. 9, this kind of invariance will not strictly hold even when $J_{\zeta}d_{\gamma}=d_{\theta}$. But given that $J_{\zeta}d_{\gamma}=d_{\theta}$, the per-iteration error will be bounded by the first term on the RHS of eqn. 11, and will thus be small provided that $d_{\gamma}$ is sufficiently small and $\zeta$ is sufficiently smooth (as shown above).

Now, suppose we generate the updates $d_{\theta}$ and $d_{\gamma}$ from curvature matrices $B_{\theta}$ and $B_{\gamma}$ according to $d_{\theta}=-\alpha B_{\theta}^{-1}\nabla h$ and $d_{\gamma}=-\alpha B_{\gamma}^{-1}\nabla_{\gamma}h$, where $\nabla_{\gamma}h$ is the gradient of $h(\zeta(\gamma))$ w.r.t. $\gamma$. Then noting that $\nabla_{\gamma}h=J_{\zeta}^{\top}\nabla h$, the condition in eqn. 13 becomes equivalent to

For this to hold, a sufficient condition is that $B_{\theta}^{-1}=J_{\zeta}B_{\gamma}^{-1}J_{\zeta}^{\top}$. Since $J_{\zeta}$ is invertible (because $\zeta$ is) an equivalent condition is

The following theorem summarizes our results so far.

Suppose that $\theta=\zeta(\gamma)$ and $B_{\theta}$ and $B_{\gamma}$ are invertible matrices satisfying

Then we have that additively updating $\theta$ by $d_{\theta}=-\alpha B_{\theta}^{-1}\nabla h$ is approximately equivalent to additively updating $\gamma$ by $d_{\gamma}=-\alpha B_{\gamma}^{-1}\nabla_{\gamma}h$, in the sense that $\zeta(\gamma+d_{\gamma})\approx\theta+d_{\theta}$, with error bounded according to

Moreover, this error can be further bounded as in eqn. 12, and will be exactly if $\zeta$ is affine. And if there is a $C\geq 0$ such that $\|H_{[\zeta]_{i}}(\gamma)\|_{2}\leq C$ for all $i$ and $\gamma$, then we can even further bound this as $\frac{1}{2}C\sqrt{n}\|d_{\gamma}\|^{2}$.

Because the error bound is zero when $\zeta$ is affine, this result will trivially extend to entire sequences of arbitrary number of steps for such $\zeta$’s. And in the more general case, since the error scales as $\alpha^{2}$, we can obtain equivalence of sequences of $T/\alpha$ steps in the limit as $\alpha\to 0$. Because the length of the updates scale as $\alpha$, and we have $T/\alpha$ of them, the sequences converges to smooth paths of fixed length in the limit. The following corollary establishes this result, under a few additional (mild) hypotheses. Its proof is in Appendix E.

Suppose that $B_{\theta}$ and $B_{\gamma}$ are invertible matrices satisfying

for all values of $\theta$. Then the path followed by an iterative optimizer working in $\theta$-space and using additive updates of the form $d_{\theta}=-\alpha B_{\theta}^{-1}\nabla h$ is the same as the path followed by an iterative optimizer working in $\gamma$-space and using additive updates of the form $d_{\gamma}=-\alpha B_{\gamma}^{-1}\nabla_{\gamma}h$, provided that the optimizers use equivalent starting points (i.e. $\theta_{0}=\zeta(\gamma_{0})$), and that either

or $d_{\theta}/\alpha$ is uniformly continuous as a function of $\theta$, $d_{\gamma}/\alpha$ is uniformly bounded (in norm), there is a $C$ as in the statement of Theorem 1, and $\alpha\to 0$.

Note that in the second case we allow the number of steps in the sequences to grow proportionally to $1/\alpha$ so that the continuous paths they converge to have non-zero length as $\alpha\to 0$.

So from these results we see that natural gradient-based methods that take finite steps will not be invariant to smooth invertible reparameterizations $\zeta$, although they will be approximately invariant, and in a way that depends on the degree of curvature of $\zeta$ and the size $\alpha$ of the step-size.

Suppose the curvature matrix $B_{\theta}$ has the form

To obtain the analogous curvature matrix $B_{\gamma}$ for the $\gamma$ parameterization we replace $f$ by $f\circ\zeta$ which gives

Then noting that $J_{f\circ\zeta}=J_{f}J_{\zeta}$, where $J_{\zeta}$ is the Jacobian of $\zeta$, we have

(Here we have used the fact that the reparameterization function $\zeta$ is independent of $x$ and $y$.) Thus, this type of curvature matrix satisfies the sufficient condition in eqn. 14.

The Hessian on the other hand does not satisfy this sufficient condition, except in certain special cases. To see this, note that taking the curvature matrix to be the Hessian gives

where $H=B_{\theta}$ is the Hessian of $h$ w.r.t. $\theta$. Thus, when the curvature matrix is the Hessian, the sufficient condition $J_{\zeta}^{\top}B_{\theta}J_{\zeta}=J_{\zeta}^{\top}HJ_{\zeta}\propto B_{\gamma}$ holds if and only if

where $\nabla L$ is the gradient of $L(y,z)$ w.r.t. $z$ (evaluated at $z=f(x,\theta)$), and we allow a proportionality constant of . Rearranging this gives

This relation is unlikely to be satisfied unless the left hand side is equal to . One situation where this will occur is when $H_{[\zeta]_{j}}=0$ for each $j$, which holds when $[\zeta]_{j}$ is an affine function of $\gamma$. Another situation is where we have $\nabla h=0$ for each $(x,y)\in S$.

A New Interpretation of the Natural Gradient

As discussed in Section 10, the negative natural gradient is given by the minimizer of a local quadratic approximation $M(\delta)$ to $h$ whose curvature matrix is the Fisher $F$. And if we have that the gradient $\nabla h$ and $F$ are computed on the same set $S$ of data points, $M(\delta$) can be written as

where $F_{R}$ is the Fisher of the predictive distribution $R_{y|z}$ (as originally defined in Section 9), $\|v\|_{F_{R}}=\sqrt{v^{\top}F_{R}v}$, and $c=h(\theta)-\frac{1}{2}(\sum_{(x,y)\in S}\nabla_{z}\log r(y|z)^{\top}F_{R}^{-1}\nabla_{z}\log r(y|z))/|S|$ is a constant (independent of $\delta$).

Note that for a given $(x,y)\in S$, $F_{R}^{-1}\nabla_{z}\log r(y|z)$ can be interpreted as the natural gradient direction in $z$-space for an objective corresponding to the KL divergence between the predictive distribution $R_{y|z}$ and a delta distribution on the given $y$. In other words, it points in the direction which moves $R_{y|z}$ most quickly towards to said delta distribution, as measured by the KL divergence (see Section 6). And assuming that the GGN interpretation of $F$ holds (as discussed in Section 9), we know that it also corresponds to the optimal change in $z$ according to the 2nd-order Taylor series approximation of the loss function $L(y,z)$.

Thus, $M(\delta)$ can be interpreted as the sum of squared distances (as measured using the Fisher metric tensor) between these “optimal” changes in the $z$’s, and the changes in the $z$’s which result from adding $\delta$ to $\theta$, as predicted using 1st-order Taylor-series approximations to $f$.

In addition to giving us a new interpretation for the natural gradient, this expression also gives us an easy-to-compute bound on the largest possible improvement to $h$ (as predicted by $M(\delta)$). In particular, since the squared error terms are non-negative, we have

Given $F_{R}=H_{L}$, this quantity has the simple interpretation of being the optimal improvement in $h$ (as predicted by a 2nd-order order model of $L(y,z)$ for each case in $S$) achieved in the hypothetical scenario where we can change the $z$’s independently for each case.

The existence of this bound shows that the natural gradient can be meaningfully defined even when $F^{-1}$ may not exist, provided that we compute $F$ and $\nabla h$ on the same data, and that each $F_{R}$ is invertible. In particular, it can be defined as the minimizer of $M(\delta)$ that has minimum norm (which must exist since $M(\delta)$ is bounded below), which in practice could be computed by using the pseudo-inverse of $F$ in place of $F^{-1}$. Other choices are possible, although care would have to be taken to ensure invariance of the choice with respect to parameterization.

Asymptotic Convergence Speed

A property of natural gradient descent which is frequently referenced in the literature is that it is “Fisher efficient”. In particular, Amari (1998) showed that an iteration of the form

which is (asymptotically) the smallestWith the usual definition of $\preceq$ for matrices: $A\preceq C$ iff $C-A$ is PSD. possible variance matrix that any unbiased estimator computed from $k$ training cases can have, according to the Cramér-Rao lower boundNote that to apply the Cramér-Rao lower bound in this context one must assume that the training data set, on which we compute the objective (and which determines $\theta^{*}$), is infinitely large, or more precisely that its conditional distribution over $y$ has a density function. For finite training sets one can easily obtain an estimator with exactly zero error for a sufficiently large $k$ (assuming a rich enough model class), and so these requirements are not surprising. If we believe that there is a true underlying distribution of the data that has a density function, and from which the training set is just a finite collection of samples, then Cramér-Rao can be thought of as applying to the problem of estimating the true parameters of this distribution from said samples (with $F$ computed using the true distribution), and will accurately bound the rate of convergence to the true parameters until we start to see samples repeat. After that point, convergence to the true parameters will slow down and eventually stop, while convergence on the training objective may start to beat the bound. This of course implies that any convergence on the training set that happens faster than the Cramér-Rao bound will necessarily correspond to over-fitting (since this faster convergence cannot happen for the test loss)..

This result can also be straightforwardly extended to handle the case where $g_{k}(\theta_{k})$ is computed using a mini-batch of size $m$ (which uses $m$ independently sampled cases at each iteration), in which case the above asymptotic variance bound becomes

which again matches the Cramér-Rao lower bound.

This result applies to the version of natural gradient descent where $F$ is computed using the training distribution $\hat{Q}_{x}$ and the model’s conditional distribution $P_{y|x}$ (see Section 5). If we instead consider the version where $F$ is computed using the true data distribution $Q_{x}$, then a similar result will still apply, provided that we sample $x$ from $Q_{x}$ and $y$ from $Q_{y|x}$ when computing the stochastic gradient $g_{k}(\theta_{k})$, and that $\theta^{*}$ is defined as the minimum of the idealized objective $\operatorname{KL}(Q_{x,y}\|P_{x,y}(\theta))$ (see Section 4).

While this Fisher efficiency result would seem to suggest that natural gradient descent is the best possible optimization method in the stochastic setting, it unfortunately comes with several important caveats and conditions, which we will discuss. (Moreover, as we will later, it is also possessed by much simpler methods, and so isn’t a great justification for the use of natural gradient descent by itself.)

Firstly, the proof assumes that the iteration in eqn. 15 eventually converges to the global optimum $\theta^{*}$ (at an unspecified speed). While this assumption can be justified when the objective $h$ is convex (provided that $\alpha_{k}$ is chosen appropriately), it won’t be true in general for non-convex objectives, such as those encountered in neural network training. In practice however, a reasonable local optimum $\theta^{*}$ might be a good enough surrogate for the global optimum, in which case a property analogous to Fisher efficiency may still hold, at least approximately.