Natural Gradient Revisited
Authors: Razvan Pascanu, Yoshua Bengio
ICLR 2013

Abstract: The aim of this paper is two-folded. First we intend to show that Hessian-Free optimization (Martens, 2010) and Krylov Subspace Descent (Vinyals and Povey, 2012) can be described as implementations of Natural Gradient Descent due to their use of the extended Gauss-Newton approximation of the Hessian. Secondly we re-derive Natural Gradient from basic principles, contrasting the difference between the two version of the algorithm that are in the literature.

Revisiting Natural Gradient for Deep Networks
Razvan Pascanu and Yoshua Bengio
Dept. IRO
University of Montreal
Montreal, QC
Abstract
The aim of this paper is three-fold. First we show that Hessian-Free (Martens,
2010) and Krylov Subspace Descent (Vinyals and Povey, 2012) can be described
as implementations of natural gradient descent due to their use of the extended
Gauss-Newton approximation of the Hessian. Secondly we re-derive natural gra-
dient from basic principles, contrasting the difference between two versions of
the algorithm found in the neural network literature, as well as highlighting a few
differences between natural gradient and typical second order methods. Lastly we
show empirically that natural gradient can be robust to overﬁtting and particularly
it can be robust to the order in which the training data is presented to the model.
1 Introduction
Several recent papers tried to address the issue of using better optimization techniques for machine
learning, especially for training deep architectures or neural networks of various kinds. Hessian-Free
optimization (Martens, 2010; Sutskever et al., 2011; Chapelle and Erhan, 2011), Krylov Subspace
Descent (Vinyals and Povey, 2012), natural gradient descent (Amari, 1997; Park et al., 2000; Le
Roux et al., 2008; Le Roux et al., 2011) are just a few of such recently proposed algorithms. They
usually can be split in two different categories: those which make use of second order information
and those which use the geometry of the underlying parameter manifold (natural gradient).
One particularly interesting pipeline to scale up such algorithms was originally proposed in Pearl-
mutter (1994), ﬁnetuned in Schraudolph (2001) and represents the backbone behind both Hessian-
Free optimization (Martens, 2010) and Krylov Subspace Descent (Vinyals and Povey, 2012). The
core idea behind it is to make use of the forward (renamed toR-operator in Pearlmutter (1994)) and
backward pass of automatic differentiation to compute efﬁcient products between Jacobian or Hes-
sian matrices and vectors. These products are used within a truncated-Newton approach (Nocedal
and Wright, 2000) which considers the exact Hessian and only inverts it approximately without the
need for explicitly storing the matrix in memory, as opposed to other approaches which perform a
more crude approximation of the Hessian (or Fisher) matrix (either diagonal or block-diagonal).
The contributions of this paper to the study of the natural gradient are as follows. We provide a de-
tailed derivation of the natural gradient, avoiding elements of information geometry. We distinguish
natural gradient descent from TONGA and provide arguments suggesting that natural gradient may
also beneﬁts from a form of robustness that should yield better generalization. The arguments for
this robustness are different from those invoked for TONGA. We show experimentally the effects
of this robustness when we increase the accuracy of the metric using extra unlabeled data. We also
provide evidence that the natural gradient is robust to the order of training examples, resulting in
lower variance as we change the order. The ﬁnal contribution of the paper is to show that Martens’
Hessian-Free approach of Martens (2010) (and implicitly Krylov Subspace Descent (KSD) algo-
rithm) can be cast into the framework of the natural gradient, showing how these methods can be
seen as doing natural gradient rather then second order optimization.
1
arXiv:1301.3584v4  [cs.LG]  13 Mar 2013
2 Natural Gradient
Natural gradient can be traced back to Amari’s work on information geometry (Amari, 1985) and
its application to various neural networks (Amariet al., 1992; Amari, 1997), though a more in depth
introduction can be found in Amari (1998); Park et al. (2000); Arnold et al. (2011). The algorithm
has also been successfully applied in the reinforcement learning community (Kakade, 2001; Peters
and Schaal, 2008) and for stochastic search (Sun et al., 2009). Le Roux et al. (2007) introduces a
different formulation of the algorithm for deep models. Although similar in name, the algorithm is
motivated differently and is not equivalent to Amari’s version, as will be shown in section 4.1.
Let us consider a family of density functions F: RP →(B →[0,1]), where for every θ ∈RP,
F(θ) deﬁnes a density function from B →[0,1] over the random variable z ∈B, where B is some
suitable numeric set of values, for e.g. B = RN. We also deﬁne a loss function that we want to
minimize L: RP →R. Any choice of θ ∈RP deﬁnes a particular density function pθ(z) = F(θ)
and by considering all possible θ values, we explore the set F, which is our functional manifold.
Because we can deﬁne a similarity measures between nearby density functions, given by the KL-
divergence which in its inﬁnitesimal form behaves like a distance measure, we are dealing with a
Riemannian manifold whose metric is given by the Fisher Information matrix. Natural gradient
attempts to move along the manifold by correcting the gradient ofLaccording to the local curvature
of the KL-divergence surface 1:
∇NL(θ) = ∂L(θ)
∂θ Ez
[(∂log pθ(z)
∂θ
)T (∂log pθ(z)
∂θ
)]−1
= ∂L(θ)
∂θ G−1 (1)
We can derive this resultwithout relying on information geometry. We consider the natural gradient
to be deﬁned as the algorithm which, at each step, picks a descent direction such that the KL-
divergence between pθ and pθ+∆θ is constant. At each step, we need to ﬁnd ∆θsuch that:
arg min∆θL(θ+ ∆θ)
s. t. KL(pθ||pθ+∆θ) = constant (2)
Using this constraint we ensure that we move along the functional manifold with constant speed,
without being slowed down by its curvature. This also makes learning robust to re-parametrizations
of the model, as the functional behaviour of pdoes not depend on how it is parametrized.
Assuming ∆θ→0, we can approximate the KL divergence by its Taylor series:
KL(pθ(z) ∥pθ+∆θ(z)) ≈ (Ez [log pθ] −Ez [log pθ]) −Ez
[∂log pθ
∂θ
]
∆θ−1
2∆θTEz
[∂2 log pθ
∂θ2
]
∆θ
= 1
2∆θTEz
[
−∂2 log pθ(z)
∂θ2
]
∆θ (3)
= 1
2∆θTEz
[(∂log pθ(z)
∂θ
)T (∂log pθ(z)
∂θ
)]
∆θ (4)
The ﬁrst term cancels out and because Ez
[
∂log pθ(z)
∂θ
]
= 0, 2 we are left with only the last term.
The Fisher Information Matrix form can be obtain from the expected value of the Hessian through
algebraic manipulations (see the Appendix).
We now express equation (2) as a Lagrangian, where the KL divergence is approximated by (4) and
L(θ+ ∆θ) by its ﬁrst order Taylor series L(θ) + ∂L(θ)
∂θ ∆θ:
1 Throughout this paper we use the mathematical convention that a partial derivative∂log pθ
∂θ is a row-vector
2Proof: Ez
[
∂log pθ(z)
∂θ
]
= ∑
z
(
pθ(z) 1
pθ(z)
∂pθ(z)
∂θ
)
= ∂
∂θ
(∑
θpθ(z)
)
= ∂1
∂θ = 0. The proof holds for
the continuous case as well, replacing sums for integrals.
2
L(θ) + ∂L(θ)
∂θ ∆θ+ 1
2λ∆θTEz
[(∂log pθ(z)
∂θ
)T (∂log pθ(z)
∂θ
)]
∆θ= 0 (5)
Solving equation (5) for ∆θ gives us the natural gradient formula (1). Note that we get a scalar
factor of 2 1
λ times the natural gradient. We fold this scalar into the learning rate, and hence the
learning rate also controls the difference between pθ and pθ+∆θ that we impose at each step. Also
the approximations we make are meaningful only around θ. Schaul (2012) suggests that using a
large step size might be harmful for convergence. We deal with such issues both by using damping
(i.e. setting a trust region around θ) and by properly selecting a learning rate.
3 Natural Gradient for Neural Networks
The natural gradient for neural networks relies on their probabilistic interpretation (which induces a
similarity measure between different parametrization of the model) given in the form of conditional
probabilities pθ(t|x), with x representing the input and t the target.
We make use of the following notation. q(x) describes the data generating distribution of x and
q(t|x) is the distribution we want to learn. y is the output of the model, and by an abuse of no-
tation, it will refer to either the function mapping inputs to outputs, or the vector of output acti-
vations. r is the output of the model before applying the output activation function σ. t(i) and
x(i) are the i-th target and input samples of the training set. Jy stands for the Jacobian matrix
Jy =


∂y1
∂θ1
.. ∂y1
∂θP
.. .. ..
∂yo
∂θ1
.. ∂y0
∂θP

. Finally, lower indices such as yi denote the i-th element of a vector.
We deﬁne the neural network loss as follows:
L(θ) = 1
n
n∑
i
[
log pθ(t(i)|y(x(i)))
]
= 1
n
n∑
i
[
log pθ(t(i)|σ(r(x(i))))
]
(6)
Because we have a conditional density function pθ(t|x) the formulation for the natural gradient
changes slightly. Each value of x now deﬁnes a different family of density functions pθ(t|x), and
hence a different manifold. In order to measure the functional behaviour of pθ(t|x) for different
values of x, we use the expected value (with respect to x ∼˜q(x)) of the KL-divergence between
pθ(t|x) and pθ+∆θ(t|x).
arg min∆θL(θ+ ∆θ)
s. t. Ex∼˜q(x) [KL(pθ(t|x)||pθ+∆θ(t|x))] = constant (7)
The metric G is now an expectation over ˜q(x) of an expectation over p(t|x). The former aver-
ages over possible manifolds generated by different choices of x, while the latter comes from the
deﬁnition of the Fisher Information Matrix.
∇NL(θ) = ∂L(θ)
∂θ Ex∼˜q(x)
[
Et∼p(t|x)
[(∂log pθ(t|x)
∂θ
)T (∂log pθ(t|x)
∂θ
)]]−1
= ∂L(θ)
∂θ G−1
(8)
Note that we use the distribution ˜qinstead of the empirical q. This is done in order to emphesis that
the theory does not force us to use the empirical distribution. However, in practice, we do want ˜qto
be as close as possible to q such that the curvature of the KL-divergence matches (in some sense)
the curvature of the error surface. To clarify the effects of ˜q let us consider an example. Assume
that ˜qis unbalanced with respect to q. Namely it contains twice the amount of elements of a classA
versus the other B. This means that a change inθthat affects elements of classA is seen as having a
larger impact on p(in the KL sense) than a change that affects the prediction of elements in B. Due
to the formulation of natural gradient, we will move slower along any direction that affectsA at the
expense of B landing with higher probability on solutions θ∗that favour predicting A. In practice
we approximate this expectation over ˜qby a sample average over minibatches.
3
In what follows we consider typical output activation functions and the metrics G they induce. In
the Appendix we provide a detailed description of how these matrices were obtained starting from
equation (8). Similar derivations were done in Park et al. (2000), which we repeat for convenience.
The formulas we get for the linear, sigmoid and softmax activation functions are:
Glinear = β2Ex∼˜q
[
∂y
∂θ
T ∂y
∂θ
]
= β2Ex∼˜q
[
JT
y Jy
]
(9)
Gsigmoid = Ex∼˜q
[
JT
y diag( 1
y(1 −y))Jy
]
(10)
Gsoftmax = Ex∼˜q
[ o∑
i
1
yi
(∂yi
∂θ
)T ∂yi
∂θ
]
(11)
To efﬁciently implement the natural gradient, we use a truncated Newton approach following the
same pipeline as Hessian-Free (Martens, 2010) (more details are provided in the Appendix). We
rely on Theano (Bergstra et al., 2010) for both ﬂexibility and in order to use GPUs to speed up.
The advantages of this pipeline are two-fold: (1) it uses the full-rank matrix, without the need for
explicitely storing it in memory and (2) it does not rely on a smoothness assumption of the metric.
Unlike other algorithms such as nonlinear conjugate gradient or BFGS, it does not assume that the
curvature changes slowly as we change θ. This seems to be important for recurrent neural networks
(as well as probably for deep models) where the curvature can change quickly (Pascanuet al., 2013).
4 Insights into natural gradient
Figure 1 considers a one hidden unit auto-encoder, where we minimize the error
(x−w·sigmoid(wx+ b) + b)2 and shows the path taken by Newton method (blue), natural gradi-
ent (gray), Le Roux’s version of natural gradient (orange) and gradient descent (purple). On the left
(larger plot) we show the error surface as a contour plot, where the x-axis represents band y-axis is
w. We consider two different starting points ([1] and [3]) and draw the ﬁrst 100 steps taken by each
algorithm towards a local minima. The length of every other step is depicted by a different shade
of color. Hyper-parameters, like learning rate and damping constant, were chosen such to improve
convergence speed while maintaining stability (i.e. we looked for a smooth path). Values are pro-
vided in the Appendix. At every point θduring optimization, natural gradient considers a different
KL divergence surface, KL(pθ||pθ+∆θ) parametrized by ∆θ, which has a minima at origin. On the
right we have contour plots of four different KL surfaces. They correspond to locations indicated by
black arrows on the path of natural gradient. The x-axis is ∆band y-axis is ∆wfor the KL surfaces
subplots. On top of these contour plots we show the direction and length of the steps proposed by
each of the four considered algorithms.
The point of this plot is to illustrate that each algorithm can take a different path in the parame-
ter space towards local minima. In a regime where we have a non-convex problem, with limited
resources, these path can result in qualitatively different kinds of minima. We can not draw any
general conclusions about what kind of minima each algorithm ﬁnds based on this toy example,
however we make two observations. First, as showed on the KL-surface plot for [3] the step taken
by natural gradient can be smaller than gradient descent (i.e. the KL curvature is high) even though
the error surface curvature is not high (i.e. Newton’s method step is larger than gradient descent
step). Secondly, the direction chosen by natural gradient can be quite different from that of gradient
descent (see for example in [3] and [4]), which can result in ﬁnding a different local minima than
gradient descent (for e.g. when the model starts at [3]).
4.1 A comparison between Amari’s and Le Roux’s natural gradient
In Le Roux et al. (2007) a different approach is taken to derive natural gradient. Speciﬁcally one
assumes that the gradients computed over different minibatches are distributed according to a Gaus-
sian centered around the true gradient with some covariance matrix C. By using the uncertainty
provided by C we can correct the step that we are taking to maximize the probability of a downward
move in generalization error (expected negative log-likelihood), resulting in a formula similar to that
4
Figure 1: Path taken by four different learning algorithms towards a local minima. Newton method
(blue), natural gradient (gray), Le Roux’s natural gradient (orange) and gradient descent (purple).
See text for details.
of natural gradient. If g = ∂L
∂θ is the gradient, then Le Roux et al. (2007) proposes following the
direction ˜g= ∂L(θ)
∂θ C−1 where C is:
C = 1
n
∑
i(g −⟨g⟩)T (g −⟨g⟩) (12)
While the probabilistic derivation requires the use of the centered covariance, equation (12), in
Le Roux et al. (2007) it is argued that using the uncentered covariance U is equivalent up to a
constant resulting in a simpliﬁed formula which is sometimes confused with the metric derived by
Amari.
U = 1
n
∑
igTg ≈E(x,t)∼q
[(
∂log p(t|x)
∂θ
)T (
∂log p(t|x)
∂θ
)]
(13)
The misunderstanding comes from the fact that the equation has the form of an expectation, though
the expectation is over the empirical distribution q(x,t). It is therefore not clear if U tells us how
pθ would change, whereas it is clear that G does. The two methods are just different, and one can
not straightforwardly borrow the interpretation of one for the other. However, we believe that there
is an argument strongly suggesting that the protection against drops in generalization error afforded
by Le Roux’sU is also a property shared by the natural gradient’sG.
If KL(p∥q) is small, than U can be seen as an approximation to G. Speciﬁcally we approximate
the second expectation from equation (8), i.e. the expectation overt ∼pθ(t|x), by a single point, the
corresponding t(i). This approximation makes sense when t(i) is a highly probable sample under p
which happens when we converge. Note that at convergence,U,G and the Hessian are very similar,
hence both versions of natural gradient and most second order methods would behave similarly.
An interesting question is if these different paths taken by each algorithm represent qualitatively dif-
ferent kinds of solutions. We will address this question indirectly by enumerating what implications
each choice has.
The ﬁrst observation has to do with numerical stability. One can express G as a sum of n×oouter
products (where nis the size of the minibatch over which we estimate the matrix andois the number
5
of output units) while U is a sum of only nouter products. Since the number of terms in these sums
provides an upper bound on the rank of each matrix, it follows that one could expect that U will be
lower rank than G for the same size of the minibatch n. This is also pointed out by Schraudolph
(2002) to motivate the extended Gauss-Newton matrix as middle ground between natural gradient
and the true Hessian3
A second difference regards plateaus of the error surface. Given the formulation of our error function
in equation (6) (which sums the log of pθ(t|x) for speciﬁc values of t and x), ﬂat regions of the
objective function are intrinsically ﬂat regions of the functional manifold4. Moving at constant speed
in the functional space means we should not get stalled near such plateaus. In Parket al. (2000) such
plateaus are found near singularities of the functional manifold, providing a nice framework to study
them (as is done for example in Rattray et al. (1998) where they hypothesize that such singularities
behave like repellors for the dynamics of natural gradient descent). An argument can also be made
in favour of U at plateaus. If a plateau at θ exists for most possible inputs x, than the covariance
matrix will have a small norm (because the vectors in each outer product will be small in value). The
inverse of U consequentially will be large, meaning that we will take a large step, possibly out of
the plateau region. This suggest both methods should be able to escape from some plateaus, though
the reasoning behind the functional manifold approach more clearly motivates this advantage.
Another observation that is usually made regarding the functional manifold interpretation is that it
is parametrization-independent. That means that regardless of how we parametrize our model we
should move at the same speed, property assured by the constraint on the KL-divergence betweenpθ
and pθ+∆θ. In Sohl-Dickstein (2012), following this idea, a link is made between natural gradient
and whitening in parameter space. This property does not transfer directly to the covariance matrix.
On the other hand Le Roux’s method is designed to obtain better generalization errors by moving
mostly in the directions agreed upon by the gradients on most examples. We will argue that the
functional manifold approach can also provide a similar property.
One argument relies on large detrimental changes of the expected log-likelihood, which is what Le
Roux’s natural gradient step protects us from with higher probability. The metric of Amari’s natural
gradient measures the expected (over x) KL-divergence curvature. We argue that if ∆θ induces a
large change in log-likelihood computed over x ∈D, where D corresponds to some minibatch, then
it produces a large change in pθ (in the KL sense), i.e. it results in a high KL-curvature. Because
we move at constant speed on the manifold, we slow down in these high KL-curvature regions,
and hence we do not allow large detrimental changes to happen. This intuition becomes even more
suitable when D is larger than the training set, for example by incorporating unlabeled data, and
hence providing a more accurate measure of how pθ changes (in the KL sense). This increase in
accuracy should allow for better predictions of large changes in the generalization error as opposed
to only the training error.
A second argument comes from looking at the Fisher Information matrix
which has the form of an uncentered weighted covariance matrix of gradients,
Ex
[∑
t pθ(t|x)
(
∂log pθ(t|x)
∂θ
)T (
∂log pθ(t|x)
∂θ
)]
. Note that these are not the gradients ∂L
∂θ
that we follow towards a local minima. By using this matrix natural gradient moves in the expected
direction of low variance for pθ. As the cost Ljust evaluates pθ at certain points t(i) for a given
x(i), we argue that with high probability expected directions of low variance for pθ correspond to
directions of low variance for L. Note that directions of high variance for Lindicates direction in
which pθ changes quickly which should be reﬂected in large changes of the KL. Therefore, in the
same sense as TONGA, natural gradient avoids directions of high variance that can lead to drops in
generalization error.
3Note that Schraudolph (2002) assumes a form of natural gradient that uses U as a metric, similar to the
Le Roux et al. (2007) proposed, an assumption made in Martens (2010) as well.
4One could argue that pmight be such that L has a low curvature while the curvature of KL is much larger.
This would happen for example if pis not very sensitive to θfor the values x(i),t(i) provided in the training
set, but it is for other pairings of xand t. However we believe that in such a scenario is more useful to move
slowly, as the other parings of xand tmight be relevant for the generalization error
6
4.2 Natural gradient descent versus second order methods
Even though natural gradient is usually assumed to be a second order method it is more useful, and
arguably more correct to think of it as a ﬁrst order method. While it makes use of curvature, it is
the curvature of the functional manifold and not that of the error function we are trying to minimize.
The two quantities are different. For example the manifold curvature matrix is positive semi-deﬁnite
by construction while for the Hessian we can have negative curvature.
To make this distinction clear we can try to see what information carries the metric that we invert
(as it was done in Roux and Fitzgibbon (2010) for Newton’s and Le Roux’s methods).
The functional manifold metric can be written as either the expectation of the Hessian ∂2pθ
∂θ2 or the
expectation of the Fisher Information Matrix
[(
∂pθ
∂θ
)T ∂pθ
∂θ
]
(see (3) and (4)). The ﬁrst form tells us
the that the matrix measures how a change inθaffects the gradients ∂pθ
∂θ of pθ (as the Hessian would
do for the error). The second form tells us how the change in the input affects the gradients ∂pθ
∂θ , as
the covariance matrix would do for Le Roux’s TONGA. However, while the matrix measures both
the effects of a change in the input and θit does so on the functional behaviour of pθ who acts as a
surrogate for the training error. As a consequence we need to look for density functions pθ which
are correlated with the training error, as we do in the examples discussed here.
Lastly, compared to second order methods, natural gradient lends itself very well to the online
optimization regime. In principle, in order to apply natural gradient we need an estimate of the
gradient, which can be the stochastic gradient over a single sample and some reliable measure of
how our model, through pθ, changes with θ (in the KL sense), which is given by the metric. For
e.g. in case of probabilistic models like DBMs, the metric relies only on negative samples obtained
from pθ and does not depend on the empirical distribution q at al Desjardins et al. (2013), while
for a second order method the Hessian would depend on q. For conditional distributions (as is the
case for neural networks), one good choice is to compute the metric on a held out subset of input
samples, offering this way an unbiased estimate of how p(t|x) changes with θ. This can easily be
done in an online regime. Given that we do not even need to have targets for the data over which we
compute the metric, as G integrates out the random variable t, we could even use unlabeled data to
improve the accuracy as long as it comes from the same distribution q, which can not be done for
second order methods.
5 Natural gradient robustness to overﬁtting
We explore the robustness hypothesis from section 4.1 empirically. The results of all experiments
carried out are summarized in table 1 present in the Appendix. Firstly we consider the effects of
using extra unlabeled data to improve the accuracy of the metric . A similar idea was proposed
in Sun et al. (2009). The idea is that for G to do a good job in this robustness sense, it has to
accurately predict the change in KL divergence in every direction. If G is estimated from too little
data (e.g., a small labeled set) and that data happens to be the training set, then it might “overﬁt” and
underestimate the effect of a change in some directions where the training data would tend to push
us. To protect us against this, what we propose here is the use alarge unlabeled set to obtain a more
generalization-friendly metric G.
Figure 2 describes the results on the Toronto Face Dataset (TFD), where using unlabeled data results
in 83.04% accuracy vs 81.13% without. State of the art is 85%Rifaiet al. (2012), though this result is
obtained by a larger model that is pre-trained. Hyper-parameters were validated using a grid-search
(more details in the Appendix).
As you can see from the plot, it suggests that using unlabeled data helps to obtain better testing error,
as predicted by our argument in Sec. 4.1. This comes at a price. Convergence (on the training error)
is slower than when we use the same training batch.
Additionally we explore the effect of using different batches of training data to compute the metric.
The full results as well as experimental setup are provided in the Appendix. It shows that, as most
second order methods, natural gradient has a tendency to overﬁt the current minibatch if both the
metric and the gradient are computed on it. However, as suggested in Vinyals and Povey (2012)
using different minibatches for the metric helps as we tend not to ignore directions relevant for other
minibatches.
7
Figure 2: (left) train error (cross entropy over the entire training set) on a log scale and (right) test
error (percentage of misclassiﬁed examples) as a function of number of updates for the Toronto
Faces Dataset. ‘kl, unlabeled‘ stands for the functional manifold version of natural gradient, where
the metric is computed over unlabeled data. for ’KL, different training minibatch’ we compute
the metric on a different minibatch from the training set, while ’KL, same minibatch’ we compute
the metric over the same minibatch we computed the gradient hence matching the standard use of
hessian-free. ’covariance’ stands for tonga that uses the covariance matrix as a metric, while msgd
is minibatch stochastic gradient descent. note that the x axis was interrupted, in order to improve the
visibility of how the natural gradient methods behave.
Figure 3: The plot describes how much the model is inﬂuenced by different parts of an online
training set, for the two learning strategies compared (minibatch stochastic gradient descent and
natural gradient descent). The x-axis indicates which part (1st 10th, 2nd 10th, etc.) of the ﬁrst half
of the data was randomly resampled, while the y-axis measures the resulting variance of the output
due to the change in training data.
6 Natural gradient is robust to the order of the training set
We explore the regularization effects of natural gradient descent by looking at the variance of the
trained model as a function of training samples that it sees. To achieve this we repeat the experiment
described in Erhan et al. (2010) which looks at how resampling different fraction of the training set
affects the variance of the model and focuses speciﬁcally to the relative higher variance of the early
examples. Our intuition is that by forbidding large jumps in the KL divergence of pθ and following
the direction of low variance natural gradient will try to limit the amount of overﬁtting that occurs
at any stage of learning.
We repeat the experiment from Erhan et al. (2010), using the NISTP dataset introduced in Bengio
et al. (2011) (which is just the NIST dataset plus deformations) and use 32.7M samples of this data.
We divide the ﬁrst 16.3M data into 10 equal size segments. For each data point in the ﬁgure, we ﬁx
9 of the 10 data segments, and over 5 different runs we replace the 10th with 5 different random sets
8
of samples. This is repeated for each of the 10 segments to produce the down curves. By looking
at the variance of the model outputs on a held out dataset (of 100K samples) after the whole 32.7M
online training samples, we visualize the inﬂuence of each of the 10 segments on the function learnt
(i.e., at the end of online training). The curves can be seen in ﬁgure 3.
There are two observation to be made regarding this plot. Firstly, it seems that early examples have a
relative larger effect on the behaviour of the function than latter ones (phenomena sometimes called
early-overﬁtting). This happens for both methods, natural gradient and stochastic gradient descent.
The second observation regards the overall variance of the learnt model.
Note that the variance at each point on the curve depends on the speed with which we move in
functional space. For a ﬁxed number of examples one can artiﬁcially tweak the curves for e.g. by
decreasing the learning rate. With a smaller learning rate we move slower, and since the model,
from a functional point of view, does not change by much, the variance is lower. In the limit, with
a learning rate of 0, the model always stays the same. If we increase the number of steps we take
(i.e. measure the variance after ktimes more samples) the curve recovers some of its shape.This is
because we allow the model to move further away from the starting point.
In order to be fair to the two algorithms, we use the validation error as a measure of how much we
moved in the functional space. This helps us to chose hyper-parameters such that after 32.7M sam-
ples both methods achieve the same validation error of 49.8% (see Appendix for hyper-parameters).
The results are consistent with our hypothesis that natural gradient avoids making large steps in
function space during training, staying on the path that induces least variance. Such large steps
may be present with SGD, possibly yielding the model to overﬁt (e.g. getting forced into some
quadrant of parameter space based only on a few examples) resulting in different models at the end.
By reducing the variance overall the natural gradient becomes more invariant to the order in which
examples are presented. Note that the relative variance of early examples to the last re-sampled
fraction is about the same for both natural gradient and stochastic gradient descent. However, the
amount of variance induced in the learnt model by the early examples for natural gradient is on the
same magnitude as the variance induce by the last fraction of examples for MSGD (i.e. in a global
sense natural gradient is less sensitive the order of samples it sees).
7 The relationship between Hessian-Free and natural gradient
Hessian-Free as well as Krylov Subspace Descent rely on the extended Gauss-Newton approxima-
tion of the Hessian, GN instead of the actual Hessian (see Schraudolph (2002)).
GN = 1
n
∑
i
[(∂r
∂θ
)T ∂2 log p(t(i)|x(i))
∂r2
(∂r
∂θ
)]
= Ex∼˜q
[
JT
r
(
Et∼˜q(t|x) [HL◦r]
)
Jr
]
(14)
The reason is not computational, as computing both can be done equally fast, but rather better
behaviour during learning. This is usually assumed to be caused by the fact that the Gauss-Newton
is positive semi-deﬁnite by construction, so one needs not worry about negative curvature issues.
In this section we show that in fact the extended Gauss-Newton approximation matches perfectly
the natural gradient metric, and hence by choosing this speciﬁc approximation, one can view both
algorithms as being implementations of natural gradient rather than typical second order methods.
The last step of equation (14) is obtained by using the normal assumption that (x(i),t(i)) are i.i.d
samples. We will consider the three activation functions and corresponding errors for which the
extended Gauss-Newton is deﬁned and show it matches perfectly the natural gradient metric for the
same activation.
For the linear output units with square errors we can derive the matrix HL◦r as follows:
HL◦rij,i̸=j = ∂2 ∑
k(rk−tk)2
∂ri∂rj
= ∂2(ri−ti)
∂rj
= 0
HL◦rii = ∂2 ∑
k(rk−tk)2
∂ri∂ri
= ∂2(ri−ti)
∂ri
= 2
(15)
9
GN = 1
n
∑
x(i),t(i)
JT
r HL◦rJr = 1
n
∑
x(i),t(i)
JT
y HL◦yJy = 1
n
∑
x(i)
JT
y (2I) Jy = 2Ex∈q(x)
[
JT
y Jy
]
(16)
The result is summarized in equation 16, where we make use of the fact that r = y. It matches the
corresponding natural gradient metric, equation (24) from section 3, up to a constant.
In the case of sigmoid units with cross-entropy objective (σis the sigmoid function), HL◦r is
HL◦rij,i̸=j = ∂2 ∑
k(−tklog(σ(rk))−(1−tk) log(1−σ(rk)))
∂ri∂rj
=
∂
(
−ti 1
σ(ri) σ(ri)(1−σ(ri))+(1−ti) 1
1−σ(ri) σ(ri)(1−σ(ri))
)
∂rj
= ∂σ(ri)−ti
∂rj
= 0
HL◦rii = ...= ∂σ(ri)−ti
∂ri
= σ(ri)(1 −σ(ri))
(17)
If we insert this back into the Gauss-Newton approximation of the Hessian and re-write the equation
in terms of Jy instead of Jr, we get, again, the corresponding natural gradient metric, equation (10).
GN = 1
n
∑
x(i),t(i) JT
r HL◦rJr = 1
n
∑
x(i) JT
r diag (y(1 −y)) diag
(
1
y(1−y)
)
diag (y(1 −y)) Jr
= Ex∼˜q
[
JT
y diag
(
1
y(1−y)
)
Jy
]
(18)
The last matching activation and error function that we consider is the softmax with cross-entropy.
HL◦rij,i̸=j = ∂2 ∑
k(−tklog(φ(rk)))
∂ri∂rj
= ∂∑
k(tkφ(ri))−ti
∂rj
= −φ(ri)φ(rj)
HL◦rii = ...= ∂φ(ri)−ti
∂ri
= φ(ri) −φ(ri)φ(ri)
(19)
Equation (20) starts from the natural gradient metric and singles out a matrixM in the formula such
that the metric can be re-written as the productJT
r MJr (similar to the formula for the Gauss-Newton
approximation). In (21) we show that indeed M equals HL◦r and hence the natural gradient metric
is the same as the extended Gauss-Newton matrix for this case as well. Note that δis the Kronecker
delta, where δij,i̸=j = 0 and δii = 1.
G = Ex∼˜q
[∑o
k=1
1
yk
(
∂yk
∂θ
)T ∂yk
∂θ
]
= Ex∼˜q
[
JT
r
(∑o
k=1
1
yk
(
∂yk
∂r
)T (
∂yk
∂r
))
Jr
]
= 1
N
∑
x(i)
(
JT
r MJr
) (20)
Mij,i̸=j = ∑o
k=1
1
yk
∂yk
∂ri
∂yk
∂rj
= ∑o
k=1(δki −yi)yk(δkj −yj) = yiyj −yiyj −yiyj = −φ(ri)φ(rj)
Mii = ∑o
k=1
1
yk
∂yk
∂yi
∂yk
∂rj
= y2
i (∑o
k=1 yk) + yi −2y2
i = φ(ri) −φ(ri)φ(ri)
(21)
8 Conclusion
In this paper we re-derive natural gradient, by imposing that at each step we follow the direction
that minimizes the error function while resulting in a constant change in the KL-divergence of the
probability density function that represents the model. This approach minimizes the amount of
differential geometry needed, making the algorithm more accessible.
We show that natural gradient, as proposed by Amari, is not the same as the algorithm proposed by
Le Roux et al, even though it has the same name. We highlight a few differences of each algorithm
and hypothesis that Amari’s natural gradient should exhibit the same robustness against overﬁtting
that Le Roux’s algorithm has, but for different reasons.
10
We explore empirically this robustness hypothesis, by proving better test errors whenunlabeled data
is used to improve the accuracy of the metric. We also show that natural gradient may reduce the
worrisome early specialization effect previously observed with online stochastic gradient descent
applied to deep neural nets, and reducing the variance of the resulting learnt function (with respect
to the sampled training data).
By computing the speciﬁc metrics needed for standard output activation functions we showed that
the extended Gauss-Newton approximation of the Hessian coincides with the natural gradient metric
(provided that the metric is estimated over the same batch of data as the gradient). Given this identity
one can re-interpret the recently proposed Hessian-Free and Krylov Subspace Descent as natural
gradient.
Finally we point out a few differences between typical second order methods and natural gradient.
The latter seems more suitable for online or probabilistic models, and relies on a surrogate probabil-
ity density function pθ in place of the error function in case of deterministic models.
Acknowledgements
We would like to thank Guillaume Desjardens, Aaron Courville, Li Yao, David Warde-Farley and
Ian Goodfellow for the interesting discussion on the topic, or for any help provided during the
development of this work. Reviewers at ICLR were particularly helpful, and we want to thank
them, especially one of the reviewers that suggested several links with work from the reinforcement
learning community. Also special thanks goes to the Theano development team as well (particularly
to Frederic Bastien, Pascal Lamblin and James Bergstra) for their help.
We acknowledge NSERC, FQRNT, CIFAR, RQCHP and Compute Canada for the resources they
provided.
References
Amari, S. (1985). Differential geometrical methods in statistics. Lecture notes in statistics, 28.
Amari, S. (1997). Neural learning in structured parameter spaces - natural Riemannian gradient. InIn Advances
in Neural Information Processing Systems, pages 127–133. MIT Press.
Amari, S., Kurata, K., and Nagaoka, H. (1992). Information geometry of Boltzmann machines. IEEE Trans.
on Neural Networks, 3, 260–271.
Amari, S.-I. (1998). Natural gradient works efﬁciently in learning. Neural Comput., 10(2), 251–276.
Arnold, L., Auger, A., Hansen, N., and Olivier, Y . (2011). Information-geometric optimization algorithms: A
unifying picture via invariance principles. CoRR, abs/1106.3708.
Bengio, Y ., Bastien, F., Bergeron, A., Boulanger-Lewandowski, N., Breuel, T., Chherawala, Y ., Cisse, M., Cˆot´e,
M., Erhan, D., Eustache, J., Glorot, X., Muller, X., Pannetier Lebeuf, S., Pascanu, R., Rifai, S., Savard, F.,
and Sicard, G. (2011). Deep learners beneﬁt more from out-of-distribution examples. In JMLR W&CP:
Proc. AISTATS’2011.
Bergstra, J., Breuleux, O., Bastien, F., Lamblin, P., Pascanu, R., Desjardins, G., Turian, J., Warde-Farley, D.,
and Bengio, Y . (2010). Theano: a CPU and GPU math expression compiler. In Proceedings of the Python
for Scientiﬁc Computing Conference (SciPy).
Chapelle, O. and Erhan, D. (2011). Improved Preconditioner for Hessian Free Optimization. InNIPS Workshop
on Deep Learning and Unsupervised Feature Learning.
Choi, S.-C. T., Paige, C. C., and Saunders, M. A. (2011). MINRES-QLP: A Krylov subspace method for
indeﬁnite or singular symmetric systems. 33(4), 1810–1836.
Desjardins, G., Pascanu, R., Courville, A., and Bengio, Y . (2013). Metric-free natural gradient for joint-training
of boltzmann machines. CoRR, abs/1301.3545.
Erhan, D., Courville, A., Bengio, Y ., and Vincent, P. (2010). Why does unsupervised pre-training help deep
learning? In JMLR W&CP: Proc. AISTATS’2010, volume 9, pages 201–208.
Kakade, S. (2001). A natural policy gradient. In T. G. Dietterich, S. Becker, and Z. Ghahramani, editors,NIPS,
pages 1531–1538. MIT Press.
Le Roux, N., Manzagol, P.-A., and Bengio, Y . (2007). Topmoumoute online natural gradient algorithm. Tech-
nical Report 1299, D´epartement d’informatique et recherche op´erationnelle, Universit´e de Montr´eal.
Le Roux, N., Manzagol, P.-A., and Bengio, Y . (2008). Topmoumoute online natural gradient algorithm. In
NIPS’07.
11
Le Roux, N., Bengio, Y ., and Fitzgibbon, A. (2011). Improving ﬁrst and second-order methods by modeling
uncertainty. In Optimization for Machine Learning. MIT Press.
Martens, J. (2010). Deep learning via hessian-free optimization. In ICML, pages 735–742.
Nocedal, J. and Wright, S. J. (2000). Numerical Optimization. Springer.
Park, H., Amari, S.-I., and Fukumizu, K. (2000). Adaptive natural gradient learning algorithms for various
stochastic models. Neural Networks, 13(7), 755 – 764.
Pascanu, R., Mikolov, T., and Bengio, Y . (2013). On the difﬁculty of training recurrent neural networks.CoRR,
abs/1211.5063.
Pearlmutter, B. A. (1994). Fast exact multiplication by the hessian. Neural Computation, 6, 147–160.
Peters, J. and Schaal, S. (2008). Natural actor-critic. (7-9), 1180–1190.
Rattray, M., Saad, D., and Amari, S. I. (1998). Natural Gradient Descent for On-Line Learning. Physical
Review Letters, 81(24), 5461–5464.
Rifai, S., Bengio, Y ., Courville, A., Vincent, P., and Mirza, M. (2012). Disentangling factors of variation for
facial expression recognition. In Proceedings of the European Conference on Computer Vision (ECCV 6) ,
pages 808–822.
Roux, N. L. and Fitzgibbon, A. W. (2010). A fast natural newton method. In J. F ¨urnkranz and T. Joachims,
editors, ICML, pages 623–630. Omnipress.
Schaul, T. (2012). Natural evolution strategies converge on sphere functions. In Genetic and Evolutionary
Computation Conference (GECCO).
Schraudolph, N. N. (2001). Fast curvature matrix-vector products. In ICANN, pages 19–26.
Schraudolph, N. N. (2002). Fast curvature matrix-vector products for second-order gradient descent. Neural
Computation, 14(7), 1723–1738.
Sohl-Dickstein, J. (2012). The natural gradient by analogy to signal whitening, and recipes and tricks for its
use. CoRR, abs/1205.1828.
Sun, Y ., Wierstra, D., Schaul, T., and Schmidhuber, J. (2009). Stochastic search using the natural gradient. In
ICML, page 146.
Sutskever, I., Martens, J., and Hinton, G. E. (2011). Generating text with recurrent neural networks. In ICML,
pages 1017–1024.
Vinyals, O. and Povey, D. (2012). Krylov Subspace Descent for Deep Learning. In AISTATS.
Appendix
8.1 Expected Hessian to Fisher Information Matrix
The Fisher Information Matrix form can be obtain from the expected value of the Hessian :
Ez
[
−∂2 log pθ
∂θ
]
= Ez
[
−
∂ 1
pθ
∂pθ
∂θ
∂θ
]
= Ez
[
− 1
pθ(z)
∂2pθ
∂θ2 +
(1
pθ
∂pθ
∂θ
)T (1
pθ
∂pθ
∂θ
)]
= −∂2
∂θ2
(∑
z
pθ(z)
)
+ Ez
[(∂log pθ(z)
∂θ
)T (∂log pθ(z)
∂θ
)]
= Ez
[(∂log pθ(z)
∂θ
)T (∂log pθ(z)
∂θ
)]
(22)
8.2 Derivation of the natural gradient metrics
8.2.1 Linear activation function
In the case of linear outputs we assume that each entry of the vector t, ti comes from a Gaussian
distribution centered around yi(x) with some standard deviation β. From this it follows that:
pθ(t|x) =
o∏
i=1
N(ti|y(x,θ)i,β2) (23)
12
G = Ex∼˜q
[
Et∼N(t|y(x,θ),β2I)
[∑o
i=1
(
∂logθp(ti|y(x)i
∂θ
)T (
∂log pθ(ti|y(x)i
∂θ
)]]
= Ex∼˜q
[∑o
i=1
[
Et∼N(t|y(x,θ),β2I)
[(
∂(ti−yi)2
∂θ
)T (
∂(ti−yi)2
∂θ
)]]]
= Ex∼˜q
[∑o
i=1
[
Et∼N(t|y(x,θ),β2I)
[
(ti −yi)2
(
∂yi
∂θ
)T (
∂yi
∂θ
)]]]
= Ex∼˜q
[∑o
i=1
[
Et∼N(t|y(x,θ),βI)
[
(ti −yi)2](
∂yi
∂θ
)T (
∂yi
∂θ
)]]
= Ex∼˜q
[∑o
i=1
[
Et∼N(t|y(x,θ),β2I)
[
(ti −yi)2](
∂yi
∂θ
)T (
∂yi
∂θ
)]]
= Ex∼˜q
[∑o
i=1
[
Et∼N(t|y(x,θ),β2I)
[
(ti −yi)2](
∂yi
∂θ
)T (
∂yi
∂θ
)]]
= β2Ex∼˜q
[∑o
i=1
[(
∂yi
∂θ
)T (
∂yi
∂θ
)]]
= β2Ex∼˜q
[
JT
y Jy
]
(24)
8.2.2 Sigmoid activation function
In the case of the sigmoid units, i,e,y = sigmoid(r), we assume a binomial distribution which gives
us:
p(t|x) =
∏
i
yti
i (1 −yi)1−ti (25)
log p gives us the usual cross-entropy error used with sigmoid units. We can compute the Fisher
information matrix as follows:
G = Ex∼˜q
[
Et∼p(t|x)
[∑o
i=1
(ti−yi)2
y2
i(1−yi)2
(
∂yi
∂θ
)T ∂yi
∂θ
]]
= Ex∼˜q
[∑o
i=1
1
yi(1−yi)
(
∂yi
∂θ
)T ∂yi
∂θ
]
= Ex∼˜q
[
JT
y diag( 1
y(1−y) )Jy
]
(26)
8.2.3 Softmax activation function
For the softmax activation function, y = softmax(r), p(t|x) takes the form of a multinomial:
p(t|x) =
o∏
i
yti
i (27)
G = Ex∼˜q
[ o∑
i
1
yi
(∂yi
∂θ
)T ∂yi
∂θ
]
(28)
8.3 Implementation Details
We have implemented natural gradient descent using a truncated Newton approach similar to the
pipeline proposed by Pearlmutter (1994) and used by Martens (2010). In order to better deal with
singular and ill-conditioned matrices we use the MinRes-QLP algorithm (Choi et al., 2011) instead
of linear conjugate gradient. Both Minres-QLP as well as linear conjugate gradient can be found im-
plemented in Theano at https://github.com/pascanur/theano optimize. We used the Theano library
(Bergstra et al., 2010) which allows for a ﬂexible implementation of the pipeline, that can automat-
ically generate the computational graph of the metric times some vector for different models:
13
import theano.tensor as TT
# ‘params‘ is the list of Theano variables containing the parameters
# ‘vs‘ is the list of Theano variable representing the vector ‘v‘
# with whom we want to multiply the metric
# ‘Gvs‘ is the list of Theano expressions representing the product
# between the metric and ‘vs‘
# ‘out_smx‘ is the output of the model with softmax units
Gvs = TT.Lop(out_smx,params,
TT.Rop(out_smx,params,vs)/(out_smx*out_smx.shape[0]))
# ‘out_sig‘ is the output of the model with sigmoid units
Gvs = TT.Lop(out_sig,params,
TT.Rop(out_sig,params,vs)/(out_sig*
(1-out_sig)*
out_sig.shape[0]))
# ‘out‘ is the output of the model with linear units
Gvs = TT.Lop(out,params,TT.Rop(out,params,vs)/out.shape[0])
The full pseudo-code of the algorithm (which is very similar to the one for Hessian-Free) is given
below. The full Theano implementation can be retrieved from https://github.com/pascanur/natgrad.
Algorithm 1 Pseudocode for natural gradient algorithm
# ‘gfn‘ is a function that computes the metric times some vector
gfn ←(lambda v→Gv)
while not early stopping condition do
g ←∂L
∂θ
# linear cg solves the linear system Gx= ∂L
∂θ
ng←linear cg(gfn, g, max iters = 20, rtol=1e-4)
# γis the learning rate
θ←θ−γng
end while
Even though we are ensured that G is positive semi-deﬁnite by construction, and MinRes-QLP is
able to ﬁnd a suitable solutions in case of singular matrices, we still use a damping strategy for two
reasons. The ﬁrst one is that we want to take in consideration the inaccuracy of the metric (which
is approximated only over a small minibatch). The second reason is that natural gradient makes
sense only in the vicinity of θ as it is obtained by using a Taylor series approximation, hence (as
for ordinary second order methods) it is appropriate to enforce a trust region for the gradient. See
Schaul (2012), where the convergence properties of natural gradient (in a speciﬁc case) are studied.
Following the functional manifold interpretation of the algorithm, we can recover the Levenberg-
Marquardt heuristic used in Martens (2010) by considering a ﬁrst order Taylor approximation, where
for any function f,
f
(
θt −ηG−1 ∂f(θt)
∂θt
T)
≈f(θt) −η∂f(θt)
∂θt
G−1 ∂f(θt)
∂θt
T
(29)
This gives as the reduction ratio given by equation (30) which can be shown to behave identically
with the one in Martens (2010).
ρ=
f
(
θt −ηG−1 ∂f(θt)
∂θt
T)
−f(θt)
−η∂f(θt)
∂θt
G−1 ∂f(θt)
∂θt
T (30)
14
8.4 Additional experimental results
For the one hidden unit auto-encoder we selected hyper-parameters such to ensure stability of train-
ing, while converging as fast as possible to a minima. We compute the inverse of the metric or
Hessian exactly (as it is just a 2 by 2 matrix). The learning rate for SGD is set to .1, for Amari’s
natural gradient .5 and for the covariance of gradience 1. (Newton’s method usually does not use
a learning rate). We damped the Hessian and the covariance of gradients by adding I and Amari’s
metric using 0.01 ·I.
8.5 Restricted MNIST experiment
For the restricted MNIST, we train a one hidden layer MLP of 1500 hidden units. The hyper-
parameters where chosen based on a grid search over learning rate, damping factor and damping
strategy. Note that beside using unlabeled data, the regularization effect of natural gradient is
strongly connected to the damping factor which accounts for the uncertainty in the metric (in a
similar way to how it does in the uncentered covariance version of natural gradient). The minibatch
size was kept constant to 2500 samples for natural gradient methods and 250 for MSGD. We used a
constant learning rate and used a budged of 2000 iterations for natural gradient and 40000 iterations
for MSGD.
We used a learning rate of 1.0 for MSGD and 5.0 for the functional manifold NGD using unlabeled
data or the covariance based natural gradient. For the functional manifold NGD using either the
same training minibatch or a different batch from the training set for computing the metric we set
the learning rate to 0.1. We use a Levenberg-Marquardt heuristic only when using unlabeled data,
otherwise the damping factor was kept constant. Its initial value was 2.0 for when using unlabeled
data, and 0.01 for every case except when using the covariance of the gradients as the metric, when
is set to 0.1.
Figure 4: (left) train error (cross entropy over the entire training set) on a log scale in order to im-
prove visibility and (right) test error (percentage of misclassiﬁed examples) as a function of number
of updates for the restricted mnist dataset.
8.6 MNIST experiment
The model used has 3 layers, where the ﬁrst two are convolutional layers both with ﬁlters of size
5x5. We used 32 ﬁlters on the ﬁrst layer and 64 on the second. The last layer forms an MLP with 750
hidden units. We used minibatches of 10000 examples (for both the gradient and the metric), and a
1
t decaying learning rate strategy. The learning rate was kept constant for the ﬁrst 200 updates and
then it was computed based on the formula l0
1+ t−200
20
, where tis the number of the current update.
We used a budged ot 2000 update.
The learning rate was set to 0.5 for the functional manifold approach when using a different batch
for computing the metric and 1.0 when using the same batch for computing the metric, or for using
the covariance of gradients as metric. We use a Levenberg-Marquardt heuristic to adapt the damping
15
Table 1: Results on the three datasets considered (restricted MNIST, MNIST and TFD). Note that
different models are used for different datasets. The training error is given as cross-entropy error,
while the test error is percentage of miss-classiﬁed examples. The algorithms name are the same as
in the legend of ﬁgure 2
DATA SET DATA FOLD MSGD KL, UNLABELED KL, DIFFERENT KL, SAME COVARIANCE
BATCH BATCH
RESTRICTED TRAIN 0.0523 0.0017 0.0012 0.0023 0.0006
MNIST TEST 5.22% 4.63% 4.89% 4.91% 4.74%
MNIST TRAIN 0.00010 0.0011 0.024
TEST 0.78% 0.82% 1.07%
TFD TRAIN 0.054 0.098
TEST 16.96% 18.87%
factor which initially is 5.0 for the functional manifold approach, and a constant damping factor of
0.1 for using the covariance as metric. These values were validated by a grid search.
8.7 TFD experiment
The Toronto Face Dataset (TFD), has a large amount of unlabeled data of poorer quality than the
training set. To ensure that the noise in the unlabeled data does not affect the metric, we compute the
metric over the training batch plus unlabeled samples. We used a three hidden layer model, where
the ﬁrst layer is a convolutional layer of 300 ﬁlters of size 12x12. The second two layers from a 2
hidden layer MLP of 2048 and 1024 hidden units respectively.
For the TFD experiment we used the same decaying learning rate strategy introduced above, in
subsection 8.6, where we computed gradients over the minibatch of 960 examples. When using the
unlabeled data, we added 480 unlabeled examples to the 960 used to compute the gradient (therefore
the metric was computed over 1440 examples) otherwise we used the same 960 examples for the
metric. In both cases we used an initial damping factor of 8, and the Levenberg-Marquardt heuristic
to adapt this damping value. Initial learning rate l0 was set to 1 in both cases.
Note that we get only 83.04% accuracy on this dataset, when the state of the art is 85.0% Rifaiet al.
(2012), but our ﬁrst layer is roughly 3 times smaller (300 ﬁlters versus 1024).
8.8 NISTP exepriment (robustness to the order of training samples)
The model we experimented with was an MLP of only 500 hidden units. We compute the gradients
for both MSGD and natural gradient over minibatches of 512 examples. In case of natural gradient
we compute the metric over the same input batch of 512 examples. Additionally we use a constant
damping factor of 3 to account for the noise in the metric (and ill-conditioning since we only use
batches of 512 samples). The learning rates were kept constant, and we use .2 for the natural gradient
and .1 for MSGD.
16
Figure 5: Train and test error (cross entropy) on a log scale as a function of number of updates for
the MNIST dataset. The legend is similar to ﬁgure 2
17


Peer review:
[Review]:
Thank you for your comments. We will soon push a revision to fix all the grammar and language mistakes you pointed out.

Regarding equation (1) and equation (7), mathbf{G} represents the Fisher Information Matrix form of the metric resulting when you consider respectively p(x) vs p(y|x). Equation (1) is introduced in section 1, which presents the generic case of a family of distributions p_{	heta}(x). From section 3 onwards we adapt these equations specifically to neural networks, where, from a probabilistic point of view, we are dealing with conditional probabilities p(y|x).

Could you please be more specific regarding the elements of the paper that you found confusing?
We would like to reformulate the conclusion to make our contributions clearer. The novel points we are trying to make are:

(1) Hessian Free optimization and Krylov Subspace Descent, as long as they use the Gauss-Newton approximation of the Hessian, can be understood as Natural Gradient, because the Gauss-Newton matrix matches the metric of Natural Gradient (and the rest of the pipeline is the same).
(2) Possibly due to the regularization effect discussed in (6), we hypothesize and support with empirical results that Natural Gradient helps dealing with the early overfitting problem introduced by Erhan et al.  This early overfitting problem might be a serious issue when trying to scale neural networks to large models with very large datasets.
(3) We make the observation that since the targets get integrated out when computing the metric of Natural Gradient, one can use unlabeled data to improve the accuracy of this metric that dictates the speed with which we move in parameter space.
(4) Natural Gradient introduced by Nicolas Le Roux et al has a fundamental difference with Amari's. It is not just a different justification, but a different algorithm that might behave differently in practice.
(5) Natural Gradient is different from a second order method because while one uses second order information, it is not the second order information of the error function, but of the KL divergence (which is quite different). For e.g. it is always positive definite by construction, while the curvature is not. Also, when considering the curvature of the KL, is not the curvature of the same surface throughout learning. At each step we have a different KL divergence and hence a different surface, while for second order methods the error surface stays constant through out learning. The second distinction is that Natural Gradient is naturally suited for online learning, provided that we have sufficient statistics to estimate the KL divergence (the metric). Theoretically, second order methods are meant to be batch methods (because the Hessian is supposed over the whole dataset) where the Natural Gradient metric only depends on the model.
(6) The standard understanding of Natural Gradient is that by imposing the KL divergence between p_{theta}(y|x) and p_{theta+delta}(y|x) to be constant it ensures that some amount of progress is done at every step and hence it converges faster. We add that it also ensures that you do not move too far in some direction (which would make the KL change quickly), hence acting as a regularizer.

Regarding the paper not being formal enough we often find that a dry mathematical treatment of the problem does not help improving the understanding or eliminating confusions.  We believe that we were formal enough when showing the equivalence between the generalized Gauss-Newton and Amari's metric. Point (6) of our conclusion is a hypothesis which we validate empirically and we do not have a formal treatment for it.