Riemannian metrics for neural networks I:
Feedforward networks

arXiv:1303.0818v5 [cs.NE] 3 Feb 2015

Yann Ollivier
February 4, 2015

Abstract
We describe four algorithms for neural network training, each adapted
to different scalability constraints. These algorithms are mathematically principled and invariant under a number of transformations in
data and network representation, from which performance is thus independent. These algorithms are obtained from the setting of differential
geometry, and are based on either the natural gradient using the Fisher
information matrix, or on Hessian methods, scaled down in a specific
way to allow for scalability while keeping some of their key mathematical properties.

The most standard way to train neural networks, backpropagation, has
several known shortcomings. Convergence can be quite slow. Backpropagation is sensitive to data representation: for instance, even such a simple
operation as exchanging 0’s and 1’s on the input layer will affect performance
(Figure 1), because this amounts to changing the parameters (weights and
biases) in a non-trivial way, resulting in different gradient directions in parameter space, and better performance with 1’s than with 0’s. (In the related context of restricted Boltzmann machines, the standard training technique by gradient ascent favors setting hidden units to 1, for similar reasons
[OAAH11, Section 5].) This specific phenomenon disappears if, instead of
the logistic function, the hyperbolic tangent is used as the activation function, or if the input is normalized. But this will not help if, for instance,
the activities of internal units in a multilayer network are not centered on
average. Scaling also has an effect on performance: for instance, a common recommendation [LBOM96] is to use 1.7159 tanh(2x/3) instead of just
tanh(x) as the activation function.
It would be interesting to have algorithms whose performance is insensitive to particular choices such as scaling factors in network construction,
parameter encoding or data representation. We call an algorithm invariant,
or intrinsic, if applying a change of variables to the parameters and activities
results in the same learning trajectory. This is not the case for backpropagation (even after changing the learning rate): for instance, changing from
1

50
BP trajectory on (x,x)
BP trajectory on (1-x,1-x)

45
40
Loss (bits)

35
30
25
20
15
10
5
0
0

20

40
60
Number of iterations

80

100

Figure 1: Backpropagation learns faster with 1’s. A neural network with
two layers of size 50 (no hidden layer) is trained to reproduce its input. A
random binary sequence x of length 50 with 75% of 1’s is generated. The
network is trained on the input-output pair (x, x) for 100 backpropagation
steps (learning rate 0.01). The experiment is repeated on the input-output
pair (1 − x, 1 − x). In both cases all initial weights are set to 0.

sigmoid to tanh activation amounts to dividing the connection weights by 4
and shifting the biases by half the weights (Eqs. (46)–(47)), which does not
preserve gradient directions.
Invariance of an algorithm means fewer arbitrary design choices, and
also more robustness: good performance on a particular problem indicates
good performance over a whole class of problems equivalent to the first one
by simple (e.g., affine) transformations.
Known invariant algorithms include Newton or quasi-Newton methods
[BL88], or the natural gradient [Ama98]. The latter, in particular, is invariant (and thus preserves performance) under a wide class of changes in
the representation of the data and of the network, while Newton-like methods are only invariant under affine transforms. However, these methods are
generally not scalable: the cost of maintaining the whole Hessian or Fisher
information matrix is quadratic in parameter dimension and prohibitive for
large networks [BL88, LBOM96]. The approximations made to ensure their
scalability, such as keeping only diagonal terms [LBOM96, SZL13], making small-rank approximations [RMB07], or using limited-memory BFGS
(e.g. [LNC+ 11] for a recent example), break their invariance properties.
Scalable Riemannian methods for neural networks. In this work
we introduce four invariant algorithms adapted to four different scalability
constraints. For this we develop a suitable theoretical framework in which
2

to build invariant algorithms, by treating neural networks in the context of
Riemannian geometry.
For a network with n units, nout output units and at most d incoming
connections per unit, processing each data sample with backpropagation has
an algorithmic cost O(nd). The most lightweight of our invariant algorithms
has O(nd) cost per sample as well. The heaviest one is more faithful to the
true natural gradient but has a cost O(nd2 + ndnout ) and thus requires that
the network is sparsely connected (the average number of units influencing
a given unit must not be too large) and that the output layer is not too
large. This latter condition is typically fulfilled in classification tasks.
The unitwise natural gradient is a scaled-down version of Amari’s natural
gradient [Ama98] in which the blocks of incoming parameters to each unit
are treated independently, thus dealing, at each unit, with a square matrix
indexed by the incoming parameters to this unit. This has been proposed
as far back as [Kur94] to train neural networks; however the algorithm in
[Kur94] is limited to networks with only one hidden layer, because it relies
on an explicit symbolic computation of entries of the Fisher matrix. Here
Proposition 27 allows for an efficient computation of the exact Fisher information matrix in arbitrary neural networks, by doing nout distinct backpropagations for each sample in the dataset. As a result, the unitwise natural
gradient is adapted to situations where both the connectivity of the network
and the output layer are reasonably small: the algorithmic cost of processing
each data sample is O(nd2 + ndnout ).
The backpropagated metric gradient can be described as a blockwise
quasi-Newton method in which several approximations (Gauss–Newton and
neglecting cross-unit terms) are used. However, we describe it in an intrinsic
way: it stems from a well-defined backpropagated metric on parameter space,
in which no approximations are involved. Invariance properties follow from
this viewpoint. It is adapted to networks with reasonably small connectivity
but output layers of arbitrary size: the cost of processing a data sample is
O(nd2 ).
The quasi-diagonal natural gradient and quasi-diagonal backpropagated
metric gradient apply a “quasi-diagonal reduction” to these two algorithms,
which removes the quadratic dependency on connectivity at each unit. This
is done in a specific way to keep some (but not all) of the invariance properties, such as insensitivity to using sigmoid or 1.7159 tanh(2x/3). The
quasi-diagonal natural gradient still requires that the output layer is not
too large, with a cost of O(ndnout ) per data sample, whereas the quasidiagonal backpropagated metric gradient has the same O(nd) complexity
as ordinary backpropagation. These quasi-diagonal methods have not been
described before, to the best of our knowledge.
In this context, we also clarify another method found in the literature
[APF00, RMB07]. It is related to, and sometimes confused with, the natural gradient (discussion in [PB13]). We call this method the outer product
3

metric gradient and introduce a scaled-down, invariant version. We prove
a novel interpretation of this method as the unique one that, at each step,
spreads the improvement most uniformly over all training samples (Proposition 15).
Organization of the text. In Section 1 we give the explicit form of the
algorithms, without justification, to serve as a reference for implementation.
In Section 2 we provide the mathematical principles behind the algorithms,
starting with the relationship between gradients, metrics, and choice of coordinates (Section 2.1), then using the tools of Riemannian geometry to build
several invariant metrics for neural networks (Section 2.2 and Appendix C)
together with a way of computing them. In Section 2.3 we introduce the
quasi-diagonal reduction of a metric. These metrics produce associated gradient descents (Section 2.4). In Section 2.5 we discuss in detail the case of the
Fisher metric for neural networks and various ways to compute or approximate it. In Section 3 we present some mathematical properties of these
algorithms, focusing on invariance by change of coordinates (Section 3.2)
and a “best-fit” interpretation (Section 3.3). Section 4 contains a set of
small-scale experiments as a proof of concept for the new algorithms.
A companion article [Oll13] develops related ideas for recurrent neural
networks and provides more in-depth experiments on complex symbolic data
sequences.
We now provide an introduction to how invariant algorithms are built,
and an overview of the experimental results.
Gradient descents and metrics. To build these invariant algorithms,
we use gradient descent in suitable invariant metrics.
Backpropagation is the simple gradient descent over parameter space.
Gradient descents follow the steepest direction of change in parameter space,
and implicitly rely on a norm (or quadratic form, or metric) to define the
steepest direction: the gradient step x ← x − η ∇f can also be rewritten
(for small enough η, up to O(η 2 ) and for regular enough functions f ) as


x ← arg min f (y) +
y

1
ky − xk2
2η



(1)

namely, the gradient descent moves into the direction yielding the smallest
values of f , penalized by the distance from the current point, measured
by the square norm ky − xk2 . For backpropagation this norm k·k2 is the
numerical change in the values of the parameters: backpropagation provides
the direction of largest improvement for a minimal change in these numerical
values. Hence simple changes in parametrization influence the behavior of
the algorithm. On the other hand, norms k·k2 based on what the network
4

does, rather than how it is represented as numbers, will lead to “intrinsic”
algorithms. This is one of the ideas behind Amari’s natural gradient.
In Section 2 we build several invariant norms, by placing neural networks
in the context of differential manifolds and Riemannian geometry. The gradient descent coming from an invariant norm (Riemannian metric) will itself
be invariant. Moreover, any gradient descent using any norm has the property that small enough learning rates ensure performance improvement at
each step.
The resulting algorithms are all invariant under a number of transformations, including affine reparametrization of the unit activities. Among
the invariance properties enjoyed by the unitwise natural gradient and the
backpropagated metric (but not their quasi-diagonal reductions) are linear
recombinations of the input received by a given unit in the network, so that
a unit receiving signals f and f + εg (as functions over the dataset) from
incoming units will learn an output correlated to g just as fast as a unit
receiving signals f and g (on the input layer this can be accounted for by
normalizing and de-correlating the inputs, but this could occur at internal
units as well). Thus these gradients have a “best-fit” interpretation (Section 3.3): at each unit they solve a least-square problem of matching input
signals and desired backpropagated output, an interpretation proposed in
[Kur94].
The quasi-diagonal reductions of these algorithms are based on the observation that there is a distinction between weights wik and wjk coming
to k from different units, but no intrinsic mathematical separation between
weights and biases. Intuitively, given that unit k receives a signal wik ai from
unit i, if we change wik to wik + δwik , the average signal to unit k will change
by δwik āi where āi is the average activation of i. Hence it might be a good
idea to automatically add −δwik āi to the bias of k, to compensate. The
quasi-diagonal algorithms we present are more sophisticated versions of this,
tuned for invariance and using weighted averages. The few added terms in
the update sometimes greatly improve performance (Fig 5 on page 50).
Arguably, none of these algorithms is second-order: the update on parameter θ takes the form θ ← θ − A(θ)−1 ∂θ f where A(θ) is a matrix depending
on the network but not on the objective function f . This matrix comes from
(Riemannian) metrics evaluating the magnitude of the effect on the output
of changes in a given direction, thus providing a suitable learning rate for
each direction, without estimating derivatives of gradients. Second-order
effects are emulated in the same way the Gauss–Newton algorithm emulates
the Newton method1 .
1

Actually, in the framework of differential geometry, without a metric, the Hessian is
only defined at local optima of the function [GHL87, paragraph 3.37], so one could say
that in such a setting the Newton method approximates the Gauss–Newton algorithm
rather than the other way around.

5

80
70

Loss (bits)

60
50

sigm backprop
tanh backprop
unitwise nat grad
backpropagated metric
quasi-diag nat grad
quasi-diag backpr metric

40
30
20
10
0
Computation time (arbitrary units)

Figure 2: Auto-encoding using a 100–30–10–30–100 deep sparsely connected
network. Comparison of backpropagation using sigmoid and tanh activation,
and the four algorithms described in Section 1, for a given computation time
budget.

Experimental proof of concept. While the focus of this article is mainly
the mathematics of neural network training, we quickly tested experimentally the implementability and impact of using the various methods. We
selected a very simple auto-encoding problem on which we expected that
any training method would perform well. A sparsely connected network
with 5 layers of size 100, 30, 10, 30, and 100 was built, and 16 random
length-100 binary strings were fed to the input layer, with the target equal
to the input. Ideally the network learns to encode each of the 16 samples
using 4 bits on the middle layer (thus with room to spare) and rewrites the
output from this. The details are given in Section 4.
Even then, backpropagation performs poorly: after 10, 000 batch passes
the average log-loss is about 36 bits per sample (out of 100) for sigmoid
backpropagation, and about 25 bits per sample for tanh backpropagation.
Note that 30 bits per sample would correspond to a method which learns
only the parameters of the output layer and can reproduce the output if
someone fills the last hidden layer with the correct 30 bits.
For the same total computation time equivalent to 10, 000 batch backpropagations2 , the quasi-diagonal algorithms have a log loss of about 1.5 to
2

as measured in CPU time on a personal computer, but this can depend a lot on
implementation details

6

3.5 bits per sample, and both the unitwise natural gradient and the backpropagated metric gradient have a log loss of 0.3 to 1.5 bit per sample, thus
essentially solving the problem. See Figure 2.
The impact of the few non-diagonal terms in the quasi-diagonal algorithms was tested by removing them. In this case the quasi-diagonal backpropagated metric gradient reduces to the diagonal Gauss–Newton method
([LBOM96, Section 7.4], also used in [SZL13]). This breaks the invariance
properties, thus the impact is different for sigmoid or tanh implementations.
The diagonal Gauss–Newton method in sigmoid implementation was found
to perform more poorly, with a final log-loss of about 12 bits per sample
(Figure 5 on page 50), while in tanh implementation it comes somewhat
close to the quasi-diagonal algorithms at about 3.5 bits per sample (presumably because in our problem, the activity of all units, not only input
units, stay perfectly centered during training). Thus the quasi-diagonal
backpropagated metric gradient can be seen as “the invariant way” to write
the diagonal Gauss–Newton method, while performance of the latter is not
at all invariant.
We also compared the exact unitwise natural gradient obtained thanks
to Proposition 27, to a variant of the natural gradient in which only the
gradient terms bb⊤ corresponding to the target for each sample are added
to the Fisher matrix ([APF00, RMB07] and Section 2.5 below). The latter,
when implemented unitwise, performs rather poorly on this auto-encoding
task, with a log loss of about 25 to 28 bits per sample. The reason, discussed
in Section 4, may be that quality of this approximation to the Fisher matrix
strongly depends on output dimensionality.
One lesson from the numerical experiments is that the regularization
term ε Id added to the matrices, needed to prevent bad behavior upon inversion, formally breaks the invariance properties: individual trajectories in
sigmoid or tanh implementations, initialized in the same way, start to differ
after a dozen iterations. Still, overall performance is not affected and is the
same in both implementations (Figure 4, p. 48).
Though the quasi-diagonal methods perform well, the only methods to
sometimes reach very small values of the loss function on this example (less
than 0.1 bit per sample) are the unitwise natural gradient and the backpropagated metric, which at each unit maintain a full matrix over the incoming
parameters and thus achieve invariance under affine recombination of incoming signals. These two methods are relevant only when network connectivity
is not too high. This highlights the interest of sparsely connected networks
from a theoretical viewpoint.

Notation for neural networks
Consider a directed neural network model: a set L of units together with
a set of directed edges i → j for i, j ∈ L, without cycle. Let Lout be the
7

output units, that is, the units with no outgoing edges3 , and similarly let Lin
be the set of units without incoming edges. Let s : R → R be an activation
function. Given an activation level for the input units, each unit j gets an
activation level
P

aj = s
(2)
i→j wij ai

depending on the activation levels of the units i pointing to j and on the
firing coefficients wij from4 i to j. Biases are treated as the weights w0j from
a special always-activated unit 0 (a0 ≡ 1) connected to every other unit. A
common choice for the activation function s is the logistic [RN03] function
1
eV
s(V ) = 1+e
V = 1+e−V , although for instance [LBOM96] recommends the
hyperbolic tangent s(V ) = tanh(V ), related to the logistic by tanh(V ) =
2( 1+e1−2V ) − 1. We refer to [RN03], which we mostly follow with minor
changes in notation.
For a given non-input unit j, we call the parameters w0j and wij for
i → j the set of incoming parameters to unit j.
The dataset for this network is a set D of inputs, where each input
x ∈ RLin is a real-valued vector over the input layer. For each input is given
a target y in an arbitrary space. We view the network as a probabilistic
generative model: given an input ai = xi on the input layer Lin , we assume
that the activations of the output layer are interpreted in a fixed way as a
probability distribution ω(y) over the target space. The goal is to maximize
the probability to output y on input x: we define the loss function
ℓ(ω, y) := − ln ω(y)

(3)

the sum of which over the dataset is to be minimized. For instance, interpreting the output layer activities (ak )k∈Lout as Gaussian variables with
mean ak and variance 1 leads to a quadratic loss function ℓ.
Example 1 (Square-loss, Bernoulli, and two classification interpretations). The square-loss interpretation of the output
layer sends the activations (ak )k∈Lout of the output layer to a random variable
Y = (Yk )k∈Lout of independent Gaussian variables, where each Yk ∼ N (ak , 1)
is a Gaussian of mean ak .
Assume that the activities ai of units in the network lie in [0; 1]. The
Bernoulli interpretation of the output layer is a Bernoulli distribution as
follows: given the activations (ak )k∈Lout of the output layer, the final output is a {0, 1}Lout -valued random variable Y = (Yk )k∈Lout of independent
Bernoulli variables, where the activations ak give the probability to have
Yk = 1, namely Pr(Yk = 1) = ak .
3

This restriction on output units is not necessary (any unit could be an output unit)
but simplifies notation.
4
What is wij for some authors is wji for others. Our convention is the same as [RN03]
but for instance [LBOM96] follows the opposite convention.

8

For classification, the interpretation must send the output activations
(ak )k∈Lout to a probability distribution over the indices k in the output
layer. In the softmax interpretation5 , the probability of class k ∈ Lout is
P
eak / k′ ∈Lout eak′ . In the spherical interpretation6 , the probability of class
P
k ∈ Lout is a2k /( k′ ∈Lout a2k′ ).

Remark 16 covers the case when the interpretation depends on extra paP
rameters ϑ, such as a softmax Pr(k) = eϑk ak / k′ ∈Lout eϑk′ ak′ with trainable
weights ϑk .
Backpropagation. A common way to train the network on a given target
value y is by backpropagation, which amounts to a gradient descent over the
parameters wij . For a given loss function, define the backpropagated value
bi at each unit i by
∂ℓ
(4)
bi := −
∂ai
so that bi indicates how we should modify the activities to decrease ℓ. Then
the value of bi satisfy the backpropagation equations [RN03] from the output
layer:
X
wij rj bj
for i 6∈ Lout
(5)
bi =
j, i→j

where the activation levels ai have first been computed by forward propagation, and where let rj stand for the value of the derivative of the activation
function at unit j:
rj := s



′

X

i→j



wij ai  =

(

aj (1 − aj ) for sigmoid activation function
1 − a2j
for tanh activation function

(6)

The backpropagated values on the output layer are computed directly by
(4), for instance,


yi − ai



yi −ai


 ai (1−ai )
bi = y − P eai
 i
eak


k∈Lout


2y
2a

i
 i −P
2
ai

k∈Lout

ak

(square-loss interpretation)
(Bernoulli interpretation)
(softmax interpretation)

(7)

(spherical interpretation)

for i ∈ Lout .
5

Usually combined with a linear activation function (s = Id) on the last layer
This latter example is motivated by a theoretical argument: the set of probability
distributions over a finite set, equipped with its Fisher metric, is isometric P
to the positive
quadrant in a sphere and so is naturally parametrized by numbers ak with
a2k = 1, and
these variables yield a slightly simpler expression for the Fisher matrix. Besides, taking
squares gives a boost to the most activated output unit, in a smooth way, as in the softmax
interpretation.
6

9

The backpropagated values are used to compute the gradient of the loss
function with respect to the parameters wij . Indeed we have
∂ℓ
= −ai rj bj
∂wij

(8)

for each edge (ij) in the network. (This includes the bias w0j using a0 ≡ 1.)
It is sometimes more convenient to work with the reduced variables
b̃i := ri bi

(9)

which satisfy the backpropagation equation
b̃i = ri

X

wij b̃j

for i 6∈ Lout

(10)

j, i→j

and

∂ℓ
= −ai b̃j
(11)
∂wij
The gradient descent with learning rate η > 0 is then defined as the
following update on the firing coefficients:
wij ← wij − η

1

∂ℓ
= wij + η ai rj bj = wij + η ai b̃j
∂wij

(12)

Four invariant gradient algorithms

We now describe four gradient algorithms for network training: the unitwise natural gradient, the quasi-diagonal natural gradient, the backpropagated metric gradient, and the quasi-diagonal backpropagated metric gradient. Each of these algorithms is adapted to a different scalability constraint.
The unitwise natural gradient requires low connectivity and a small output layer; the quasi-diagonal natural gradient requires a small output layer;
the backpropagated metric gradient requires low connectivity; the quasidiagonal backpropagated metric gradient has the same asymptotic complexity as backpropagation.
These algorithms are the implementation of the more general versions
described in Section 2. As they are designed for invariance properties, implementing them using either sigmoid or tanh activation function should result
in the same output, learning trajectories, and performance, provided the initialization is changed accordingly (Section 3.2). However, the Bernoulli and
classification interpretations of the output layer assumes that the activities
lie in [0; 1], as in sigmoid activation.
We first present the algorithms in a batch version. It is straightforward
to adapt them to use random mini-batches from the dataset. In Section 1.3
they are also adapted to an online setting: this can be done using standard
techniques because the main quantities involved take the form of averages
over the dataset, which can be updated online.
10

1.1

Unitwise natural gradient and quasi-diagonal natural gradient

The unitwise natural gradient has been proposed as far back as [Kur94]
to train neural networks; however the presentation in [Kur94] is limited to
networks with only one hidden layer, because it relies on an explicit symbolic
computation of entries of the Fisher matrix. Proposition 27 below allows
for an efficient computation of the exact Fisher information matrix by doing
nout distinct backpropagations for each sample in the dataset. This relies
on linearity of backpropagation, as follows.
Definition 2 (Backpropagation transfer rates). Fix an input
x for the network and compute the activities by forward propagation. Let
k be a unit in the network and kout be a unit in the output layer. Define
the backpropagation transfer rates Jkkout from kout to k by backpropagating
the value 1 at kout . Formally:
(

out
:= 1 ,
Jkkout := 0,
Jkkout
P
Jkkout := j, k→j wkj rj Jjkout

for k 6= kout in the output layer Lout
(13)
for non-output units k

where rj is the derivative of the activation function, given by (6).
These transfer rates have the property that if backpropagation values b
P
are set on the output layer, then bk = kout ∈Lout Jkkout bkout for any unit k
(see also [LBOM96, Section 7.2]).
Computation of the transfer rates can be done by nout distinct backpropagations. There are further simplifications, since the transfer rates for
k in the input layer are never used (as there are no incoming parameters),
and the transfer rates on the last hidden layer are readily computed as
Jkkout = wkkout rkout . Thus it is enough to backpropagate the transfer rates
from the last hidden layer to the first hidden layer. In particular, with only
one hidden layer (the case considered in [Kur94] for the Fisher matrix) no
backpropagation is needed.
Definition 3 (Fisher modulus). Fix an input x for the network
and compute the activities by forward propagation. For each unit k in the
network, define the Fisher modulus Φk (x) of unit k on input x as follows,
depending on output layer interpretation.
• For the Bernoulli interpretation, set
Φk (x) :=

X

kout ∈Lout

(Jkkout )2
akout (1 − akout )

(14)

• For the square-loss interpretation, set
Φk (x) :=

X

kout ∈Lout

11

(Jkkout )2

(15)

• For the softmax interpretation, set
1
Φk (x) :=
S
where S :=

X

kout ∈Lout

P



X
1
eakout (Jkkout )2 − 2 
S
k ∈L

kout ∈Lout

out

out

eakout .

2

eakout Jkkout 

(16)

• For the spherical interpretation, set
4
S

Φk (x) :=
where S :=

X

(Jkkout )2 −

kout ∈Lout

P

kout ∈Lout

a2kout .



4 
S2 k

X

out ∈Lout

2

akout Jkkout 

(17)

Definition 4 (Unitwise Fisher matrix). Let k be a unit in the
network. Let Ek be the set of incoming units to k (including the alwaysactivated unit 0). The unitwise Fisher matrix at unit k is the (#Ek )×(#Ek )
matrix F (k) defined by
(k)
Fij := Ex∈D ai aj rk2 Φk

(18)

for i and j in Ek (including the unit 0 with a0 ≡ 1), with Φk the Fisher
modulus of Definition 3. Here Ex∈D represents the average over samples
x in the dataset (all the terms ai , aj , rk , Φk depend on the input to the
network).
(k)

By Proposition 27 below, Fij is the block of the Fisher information
matrix associated with the incoming parameters to k, hence the name.
Definition 5 (Unitwise natural gradient). The unitwise natural gradient with learning rate η > 0 updates the parameters of the network
as follows.
For each unit k, define the vector G(k) by
(k)

Gi

:= −Ex∈D

∂ℓ(y)
= Ex∈D ai rk bk
∂wik

(19)

where bk is the backpropagated value at k obtained from the target y associated with x, and where i runs over the incoming units to k (including the
always-activated unit i = 0). Compute the vector
δw(k) := (F (k) )−1 G(k)

(20)

Then the parameters of the network are updated by
(k)

wik ← wik + ηδwi
12

(21)

The asymptotic algorithmic cost of the unitwise natural gradient is as
follows. Computing τ requires nout distinct backpropagations for each input
x. For a network with n units, nout output units, and at most d incoming
connections per unit, this costs O(nout nd) per data sample (this can be
reduced in some cases, as discussed above). Computing Φ takes O(nnout ).
(k)
Computing Fij for every k takes O(nd2 ). Computing the gradient requires
inverting these matrices, which takes O(nd3 ) but is done only once in each
(mini-)batch, thus if the size of batches is larger than d this cost is negligible;
if the size of batches is smaller than d or in an online setting, inversion can
be done recursively using the Sherman–Morrison formula at a cost of O(nd2 )
per data sample. So the overall cost of this gradient is O(nd2 + nout nd) per
data sample.
This algorithmic cost is fine if connectivity is low. We now define a more
light-weight version in case connectivity is large. Its computational cost is
equivalent to that of ordinary backpropagation provided the output layer is
small.
Definition 6 (Quasi-diagonal natural gradient). The quasidiagonal natural gradient with learning rate η > 0 updates the parameters
of the network as follows.
(k)
(k)
(k)
For each unit k, compute only the entries F00 , F0i , and Fii of the
(k)
unitwise Fisher matrix at k. Define the vector G as in (19) above. Define
the vector δw(k) by
(k)
δwi
(k)
δw0

(k)

(k)

(k)

(k)

=

Gi F00 − G0 F0i

=

G0

(k)

(k)

for i 6= 0

(k)

Fii F00 − (F0i )2
(k)

(k)
F00

−

X F (k)
0i

(k)
i6=0 F00

(k)

δwi

(22)
(23)

and update the parameters of the network by
(k)

wik ← wik + ηδwi

(24)

With respect to the unitwise natural gradient, the algorithmic cost does
not involve any O(nd2 ) terms because we only compute O(d) entries of the
matrices F and do not require a matrix inversion. The variables τ and Φ still
need to be computed. Thus the overall complexity is reduced to O(nout nd).
The quasi-diagonal formulas may seem arbitrary. If we remember that
(omitting the superscript (k)) F00 = Ex∈D rk2 Φk , F0i = Ex∈D ai rk2 Φk , and
Fii = Ex∈D a2i rk2 Φk , we can consider these sums as expectations over the
dataset with weights rk2 Φk . Then the weighted average of ai is Ai = F0i /F00
and its weighted variance is Vi = Fii /F00 − A2i so that we have
δwi =

Gi − G0 Ai
F00 Vi
13

(25)

and in particular the denominator is always positive unless the activity of
unit i is constant (in this case, the numerator vanishes too).
A possible interpretation is as follows: If the activity ai of unit i is not
centered on average over the dataset (with the weights above), increasing a
weight wik not only increases the influence of unit i over unit k, but also
shifts the average activity of unit k, which may not be desirable. Using the
method above, if ai is not centered, when we change wik a corresponding
term is automatically subtracted from the bias w0k so as not to shift the
average activity of unit k, as discussed in the introduction. On the other
hand if the activity ai is centered, then the update is diagonal, and scaled
by the inverse “variance” 1/Vi .

1.2

Backpropagated metric gradient and quasi-diagonal backpropagated metric gradient

Computing the Fisher matrix as above requires performing nout backpropagations for each sample. If one tries to compute the Fisher modulus Φ
directly by backpropagation, the backpropagation equation involves crossterms between different units. Neglecting these cross-terms results in a simpler version of the Fisher modulus which can be computed in one backward
pass; the corresponding backpropagation equation is well-known as an approximation of the Hessian [LBOM96, Section 7]. It turns out this quantity
and the associated metric are still intrinsic.
Definition 7 (Backpropagated modulus). Fix an input x for the
network and compute the activities by forward propagation. Define the
backpropagated modulus mk (x) for each unit k by
mk (x) :=

X

2 2
wkj
rj mj (x)

(26)

j, k→j

if k is not an output unit, and, depending on output interpretation,

mk (x) :=


1


 ak (1−ak )


1
a

(Bernoulli)
a

e k
e k


S (1 − S ),


2

ak
4
S (1

−

S

),

for k in the output layer.

S=

akout

kout ∈Lout e
P
2
kout ∈Lout akout

S=

P

(square-loss)
(softmax)

(27)

(spherical)

Definition 8 (Backpropagated metric). Let k be a unit in the
network. Let Ek be the set of incoming units to k (including the alwaysactivated unit 0). The backpropagated metric at unit k is the (#Ek )×(#Ek )
matrix M (k) defined by
(k)
Mij := Ex∈D ai aj rk2 mk

14

(28)

for i and j in Ek (including the unit 0 with a0 ≡ 1). Here Ex∈D represents
the average over samples x in the dataset (all the terms ai , aj , ak , mk depend
on the input to the network).
The backpropagated metric gradient can thus be described as an approximate, blockwise Hessian method in which the Hessian is approximated by
the Gauss–Newton technique with, in addition, cross-unit terms neglected.
Such a method turns out to be intrinsic.
Definition 9 (Backpropagated metric gradient). The backpropagated metric gradient with learning rate η > 0 updates the parameters
of the network as follows.
For each unit k, define the vector G(k) by
(k)

Gi

:= −Ex∈D

∂ℓ(y)
= Ex∈D ai rk bk
∂wik

(29)

where bk is the backpropagated value at k obtained from the target y associated with x, and where i runs over the incoming units to k (including the
always-activated unit i = 0). Compute the vector
δw(k) := (M (k) )−1 G(k)

(30)

Then the parameters of the network are updated by
(k)

wik ← wik + ηδwi

(31)

The algorithmic cost of the backpropagated metric gradient is O(nd2 )
per data sample, with notation as above. Indeed, computing m costs the
same as a backpropagation pass, namely O(nd). Computing the matrices M
costs O(nd2 ). Inverting the matrices has no impact on the overall complexity,
as explained after Definition 5. This cost is acceptable for small d (sparsely
connected networks). For large d we define the following.
Definition 10 (Quasi-diagonal backpropagated metric gradient). The quasi-diagonal backpropagated metric gradient with learning
rate η > 0 updates the parameters of the network as follows.
(k)
(k)
(k)
For each unit k, compute only the entries M00 , M0i , and Mii of
backpropagated metric at k. Define the vector G(k) as in (19) above. Define
the vector δw(k) by
(k)
δwi
(k)

(k)

=

δw0 =

(k)

(k)

(k)

Gi M00 − G0 M0i
(k)

(k)

(k)

Mii M00 − (M0i )2
(k)

G0

(k)
M00

−

X M (k)
0i
(k)
i6=0 M00

(k)

δwi

for i 6= 0

(32)
(33)

and update the parameters of the network by
(k)

wik ← wik + ηδwi
15

(34)

The same remarks as for the quasi-diagonal natural gradient apply for
(k)
(k)
(k)
interpreting the various terms. The denominator Mii M00 −(M0i )2 can be
seen as a weighted variance of the activity of unit i, and is positive unless ai
(k)
(k)
is constant over the dataset. The contribution of δwi to δw0 compensates
the change of average activity induced by a change of wik .
The asymptotic cost of this update is O(nd) per data sample, as for
backpropagation.
If, in the quasi-diagonal backpropagated metric gradient, the non-diagonal
(k)
terms are neglected (M0i is set to 0), then this reduces to the diagonal
Gauss–Newton method equations from [LBOM96, Section 7.4] (also used
for instance in [SZL13]).
Remark 11. On the incoming parameters to the output layer, the unitwise natural gradient and the backpropagated metric gradient coincide if
the Bernoulli or square-loss interpretation is used. (Actually, with learning
rate η = 1 they also both coincide with the Newton method restricted to
the output layer parameters.)
Remark 12. Since these algorithms rely on inverting matrices, regularization is an issue. In practice, terms ε Id have to be added to F and M
(k)
(k)
before inversion; terms ε have to be added to the diagonal terms F00 , Fii ,
(k)
(k)
M00 and Mii in the quasi-diagonal reduction. This formally breaks the
invariance properties. Section 3.3 elaborates on this. Still, this operation
preserves the guarantee of improvement for small enough learning rates.

1.3

Adaptation to an online setting

The unitwise natural gradient and unitwise backpropagated metric gradient
both update the weights by
δw = A−1 G

(35)

with G the gradient of the loss function over the dataset, and A a positivedefinite, symmetric matrix. A key feature here is that the matrix A takes
(k)
the form of an expectation over the dataset: Fij = Ex∈D ai aj rk2 Φk for the
(k)

Fisher matrix, and Mij = Ex∈D ai aj rk2 mk for the backpropagated metric.
Any such algorithm can be turned online using a standard construction
as follows (compare e.g. [RMB07]). Another possibility is, of course, to use
mini-batches.
In the following, A stands for either the unitwise Fisher matrix or the
backpropagated metric. Let A(x) be the corresponding contribution of each
(k)
input x in the expectation, namely, A(x)ij = ai aj rk2 Φk for the Fisher
(k)

metric and A(x)ij
A = Ex∈D A(x).

= ai aj rk2 mk for the backpropagated metric, so that

16

At each step t, we use one new sample in the dataset, update an estimate
A(t) of A, and follow a gradient step for this sample, as follows.
• Initialize the matrix A(0) by using a small subsample Dinit ⊂ D, for
instance the first ninit samples in the dataset:
A(0) := Ex∈Dinit A(x)

(36)

• Fix a discount factor 0 < γ < 1. For each new sample xt , compute its
contribution A(xt ) and update A by
A(t) := (1 − γ)A(t−1) + γA(xt )

(37)

• Compute the inverse of A(t) from the inverse of A(t−1) by the Sherman–
Morrison formula at each unit, using that A(xt ) is a rank-one matrix
at each unit. (This way matrix inversion is no more costly than the
rest of the step.)
• Compute the negative gradient G(xt ) of the loss function on input xt
by backpropagation.
• Update the parameters by
w ← w + ηt (A(t) )−1 G(xt )

(38)

where ηt is the learning rate.
• For the quasi-diagonal reductions of the algorithms, only the entries
A00 , Aii and A0i of the matrix A are updated at each step. No matrix
inversion is required for the update equations (22)–(23) and (32)–(33).
We could also initialize A(0) to a simple matrix like Id, but this breaks
the invariance properties of the algorithms.
The update rule for A(t) depends on the discount factor γ. It should
be large enough so that a large number of data points contribute to the
computation of A, but small enough to be reactive so that A evolves as
training gets along. In our setting, from the particular form of A(x) at a
unit k we see that each A(xt ) contributes a rank-one matrix. This means
that γ should be much smaller than 1/nk with nk the number of parameters
at unit k, because otherwise the estimated matrix A(t) will be close to a
low-rank matrix, and presumably a poor approximation of the true matrix
A, unreliable for numerical inversion.
The same remark applies to the number ninit of samples used for initialization: it should be somewhat larger than the number of parameters at
each unit, otherwise A(0) will be of low rank.

17

2

Constructing invariant algorithms: Riemannian
metrics for neural networks

2.1

Gradient descents and metrics, natural metrics

The gradient of a function f on Rd gives the direction of steepest ascent:
among all (very small) vectors with a given norm, it provides the greatest
variation of f . Formally, the gradient ∇f of a smooth function f is defined
by the property that
f (x + εv) = f (x) + εh∇f, vi + O(ε2 )

(39)

for any vector v, for small enough ε. This depends on the choice of a scalar
product h·, ·i. In an orthonormal basis, the coordinates of the gradient are
simply the partial derivatives ∂f /∂xi so that gradient descent is
xi ← xi − η ∂f /∂xi

(40)

in an orthonormal basis.
For a given norm of the vector v, the quantity h∇f, vi is maximal when v
is collinear with ∇f : so the gradient ∇f indeed gives the direction of steepest
ascent among all vectors with a given norm. The gradient step x ← x−η ∇f
can actually be rewritten (for small enough η, up to O(η 2 ) and for regular
enough functions f ) as
1
ky − xk2
x ← arg min f (y) +
2η
y




(41)

namely, the gradient descent moves into the direction yielding the smallest
values of f , penalized by the distance from the current point7 . This makes
it clear how the choice of the scalar product will influence the direction of
the gradient ∇f : indeed, another scalar product will define another norm
kvk2 = hv, vi for the vector v, so that the steepest direction among all vectors
with a given norm will not be the same. The norm thus defines directions v
in which it is “cheap” or ”expensive” to move; the gradient is the direction
of steepest ascent taking this into account.
If we happen to work in a non-orthonormal basis of vectors v1 , . . . , vd ,
the gradient is given by A−1 ∂f /∂x where ∂f /∂x is the vector of partial
derivatives with respect to the coordinates xi in the basis, and A is the
symmetric matrix made of the scalar products of the basis vectors with
themselves: Aij := h vi | vj i. Indeed, the change of variable x̃ := A1/2 x
provides an orthonormal basis, thus gradient descent for x̃ is x̃ ← x̃ −
7

This can be used to define or study the direction of the gradient in more general
metric spaces (e.g., [AGS05, Chapter 2]).

18

η ∂f /∂ x̃ = x̃ − ηA−1/2 ∂f /∂x. Thus, translating back on the variable x, the
gradient descent of f takes the form
x ← x − η A−1 ∂f /∂x

(42)

Conversely, we can start with a norm on Rd defined through a positivedefinite, symmetric matrix A (which thus defines “cheap” and ”expensive”
directions). The gradient descent using this norm will then be given by (42).
So any update of the form (42) above with A a symmetric, positivedefinite matrix can be seen as the gradient descent of f using some norm.
The matrix A may even depend on the current point x, defining a Riemannian metric in which the norm of a small change δx of x is
kδxk2 :=

X

δxi Aij (x)δxj = δx⊤A(x) δx

(43)

ij

An important feature of gradient descent, in any metric, is that for η
small enough, the step is guaranteed to decrease the value of f .
The choice of a metric A represents the choice of a norm for vectors in
parameter space. Conversely, choosing a set of parameters and using the
“naive” gradient ascent for these parameters amounts to implicitly deciding
that these parameters form an orthonormal basis.
∂ℓ
For the neural network above, the gradient ascent wij ← wij − η ∂w
ij
corresponds to the choice of A = Id on parameter space, namely, the norm
P
of a change of parameters δw = (δw ij ) is kδwk2 := | δw ij |2 . Thus, gradient
descent for wij gives the best δw for a given norm kδwk, that is, the best
change of f for a given change in the numerical value of the parameters.
Example: from sigmoid to tanh activation function. Neural networks using the sigmoid and tanh activation function are defined, respectively, by
X
ai wik )
(44)
ak = sigm(
i, i→k

and
a′k = tanh(

X

′
a′i wik
)

(45)

i, i→k
′ ). Since tanh(x) = 2 sigm(2x) − 1, the
(including the biases w0k and w0k
activities of the network correspond to each other via a′k = 2ak − 1 for all k
if we set
wik
′
(46)
wik
=
4
for i 6= 0, and
w0k
1X
′
w0k
=
wik
(47)
+
2
4 i6=0

19

for the biases.
Consequently, while the gradient for the sigmoid function will try to
improve performance while minimizing the change to the numerical values
of wik and w0k , the gradient for the tanh function will do the same for the
′ and w′ , obviously resulting in different updates. If
numerical values of wik
0k
we follow the tanh gradient and rewrite it back in terms of the variables wik ,
we see that the tanh update expressed in the variables wik is
1
wik ← wik + 16(δwik − δw0k )
2

(i 6= 0)

(48)

1
(δwik − δw0k )
2
i6=0

(49)

and
w0k ← w0k + 4δw0k − 8

X

where δwik is the update that would have been applied to wik if we were
following the standard sigmoid backpropagation. Indeed this takes the form
of a symmetric matrix applied to δwik (the cross-contributions of δw0k to
wik and of δwik to w0k are the same).
Apart from an obvious speedup factor, an important difference between
this update and ordinary (sigmoid) backpropagation on the wik is that each
time a weight wik is updated, there is an opposite, twice as small contribution
to w0k : in this sense, it is as if this update assumes that the activities ai are
centered around 1/2 so that when wik gets changed to wik + c, one “needs”
to add −c/2 to the bias so that things stay the same on average.
Newton’s method and gradient descent. To find the minimum of a
function f on R, one can use the Newton method to solve f ′ = 0, namely,
x ← x − f ′ (x)/f ′′ (x). In higher dimension this becomes
x ← x − (Hess f )−1 ∂f /∂x

(50)

where ∂f /∂x is the vector of partial derivatives, and (Hess f )ij := ∂ 2 f /∂xi ∂xj
is the Hessian matrix of f .
Around a non-degenerate minimum of f , the Hessian Hess f will be a
positive-definite matrix. So the Newton method can be seen as a gradient
descent with learning rate η = 1, in the metric A = Hess f , when one is
close enough to a minimum.
Intrinsic metrics. There could be a lot of arguing and counter-arguing
about the “right” way to write the parameters with respect to which the
gradient should be taken. The solution to avoid these choices is known:
use metrics that depend on what the system does, rather than on how the
parameters are decomposed as numbers.
20

The Fisher metric, which defines a natural gradient [AN00], is one such
metric. Namely: the size (norm) of a change of parameters is measured by
the change it induces on the probability distribution of the output of the
model. The symmetric matrix A used in the gradient update is then the
Fisher information matrix. We will use scaled-down versions of the Fisher
metric for better scalability.
We present another metric for neural networks, the backpropagated metric. The size of a change of parameters at a given unit is measured by the
effect it has on the units it directly influences, which is itself measured recursively in the same way up to the output layer. The matrix defining this
metric is obtained by well-known equations related to the Gauss–Newton
approximation of the Hessian.

2.2

Intrinsic metrics and their computation by backpropagation

Here we rewrite the definition of neural networks in the language of differential manifolds and Riemannian geometry; this allows us to define metrics
directly in an intrinsic way.
Consider a neural-like network made of units influencing each other. The
activity of each unit k takes values in a space Ak which we assume to be
a differentiable manifold (typically R without a preferred origin and scale,
but we allow room for multidimensional activations). Suppose that the
activation of the network follows
ak = fθkk (ai1 , . . . , aink )

(51)

where ai1 , . . . , aink are the units pointing to k, and where fθkk is a function
from Ai1 ×· · ·×Aink to Ak , depending on a parameter θk which itself belongs
to a manifold Θk . Here we have no special, always-activated unit coding for
biases: the biases are a part of the parameters θk .
We shall also assume that the output units in the network are interpreted
through a final decoding function to produce an object ω = ω((ak )k∈Lout )
relevant to the initial problem, also assumed to belong to a differentiable
manifold.
To implement any gradient ascent over the parameters θ, we first need a
(Riemannian) metric on the parameter space. Such a metric can be defined
by choosing a parametrization by Rd and deciding that the elementary vectors of Rd are orthogonal, but this is not intrinsic: different parametrizations
will lead to different learning trajectories.
In this setting, an object is said to be intrinsic if it does not depend
on a choice of parametrization of any of the manifolds involved (activities,
parameters, final output). Hopefully, casting the activities as elements in
an abstract manifold, and writing intrinsic algorithms that will not depend
21

on how this manifold is represented as numbers, allows the algorithms to be
agnostic as to any physical interpretation of these activities (activation levels,
activation frequencies, log-frequencies, synchronized activity of a group of
neurons...)
We assume that we are given a meaningful Riemannian metric on the
final output ω: that is, we know how to measure the size of a change in the
output. For instance, if ω describes a probability distribution over a target
variable y, we can use the Fisher metric over ω (usually simpler to work
with than the Fisher metric over the whole network parameter space). In
the case ω is a Gaussian of fixed variance centered on the output values, this
coincides with a Euclidean norm on these values.
Then there are several possibilities to define intrinsic Riemannian metrics
on parameter space. The most direct one is the Fisher metric: the output is
seen as a function of all parameters, and the norm of a change of parameter
δθ (over all parameters at once) is the norm of the change it induces on the
output. This is not scalable: for a neural network with n units, nout output
units, and d incoming edges per unit, processing each data sample takes
O(n2 d2 + n2 nout ), compared to O(nd) for backpropagation.
A more scalable version is to break down the change of parameter into
a sum of changes of incoming parameters to each unit and take the Fisher
metric at each unit independently. This is the unitwise Fisher metric. As
we will see, it scales well to sparsely connected networks if the output layer
is not too large: processing each data sample takes O(nd2 + nout nd).
An even simpler version is the backpropagated metric, defined by backwards induction from the output layer: the norm of a change of parameter
on the output layer is the norm of the change it induces on the final result,
and the norm of a change of parameter at an internal unit is the sum of the
norm of the resulting changes at the units it influences directly. Processing
each data sample takes O(nd2 ).
Quasi-diagonal reduction (Section 2.3) further produces simplified intrinsic metrics in which the O(nd2 ) terms reduce to O(nd).
Notation. In what follows, we use the standard objects of differential geometry but try to present them in an intuitive way; Appendix C gives a
fully formal treatment. The notation δa, δθ, δω denotes tangent vectors
on the corresponding manifolds (intuitively, differences between two very
∂ai
denotes the differential (total
close values of a or θ or ω). The notation ∂a
k
derivative) of the activity ai seen as a function of ak . In a basis it is represented as the Jacobian matrix of partial derivatives of the components of ai
w.r.t. those of ak (if activities are more than 1-dimensional). In particular,
for a numerical function f , ∂f
∂θ is represented as a row vector, and for an
infinitesimal change (tangent vector) δθ we have δf = ∂f
∂θ δθ. The various
metrics involved are (0, 2)-tensors, but we use standard matrix notation for
them. With this convention a metric gradient descent on θ takes the form
22

⊤

θ ← θ − M (θ)−1 ∂f
∂θ .
Definition 13 (Natural metric, unitwise natural metric,
P
backpropagated metric). Let kδωk2 = Ωij δωi δωj = δω⊤Ωδω be a
metric on the final output of the network, given by the symmetric, positivedefinite matrix Ω. We define three metrics on the parameter set.
• The natural metric on the parameter set θ = (θk ) is defined as follows.
Let x be an input in the dataset D and let ω(x) be the final output of
the network run with input x and parameter θ. Let δθ be a variation
of θ and let δω(x) be the resulting variation of ω(x). Let
kδθk2nat,x := kδω(x)k2

(52)

and then define the natural metric by
kδθk2nat := Ex∈D kδθk2nat,x

(53)

∂ω
In matrix form, we have δω(x) = ∂ω(x)
∂θ δθ where ∂θ is the Jacobian
matrix of ω(x) as a function of θ, so that the natural metric is given
by the matrix

∂ω(x)⊤ ∂ω(x)
δθ
Ω
kδθk2nat = Ex∼D δθ⊤
∂θ
∂θ
The natural metric is given by a matrix of size dim θ =
• The unitwise natural metric on the parameter set θ is
kδθk2u-nat :=

X

kδθk k2nat

(54)
P

k

dim θk .

(55)

k

where k runs over the units of the network and δθk is the variation of
the incoming parameters to unit k. This metric is given by keeping
only the block-diagonal terms incoming to each unit in the matrix
defining the natural metric.
In case ω is a probability distribution and the metric Ω on ω is the
Fisher metric, we also call kδθknat and kδθku-nat the Fisher metric and
unitwise Fisher metric.
• The backpropagated metric on θ is defined as follows. Let x be an
input in the data. We first define a metric on each of the activities ak ,
depending on the input x, working from the output layer backwards.
Given a change δakout in the activity at an output unit kout , let δω(x)
be the corresponding change in the final output and set
kδakout k2bp,x := kδω(x)k2
23

(56)

The metric on internal units k is defined as follows: Given a change
δak in the activity of unit k, let δai be the resulting changes in the
activities of units k → i directly influenced by k. Define by induction
from the output layer
kδak k2bp,x :=

X

kδai k2bp,x

(57)

i, k→i

Given a change δθk of the incoming parameters to unit k, let δak,x be
the change of activity of unit k resulting from the change δθk , when
the network is run on input x. Define the backpropagated metric by
kδθk k2bp := Ex∈D kδak,x k2bp,x
and

kδθk2bp :=

X

kδθk k2bp

(58)
(59)

k∈L

Another metric, the outer product (OP) metric, can be defined from
slightly different ingredients. It corresponds to an often-used variant of the
natural gradient (e.g., [APF00, RMB07]), in which the expectation under
the current probability distribution is replaced with a similar term involving
only the desired target y for each input x (more details in Section 2.5). It
can readily be computed by backpropagation from (62).
Whereas the metrics above depend on the actual output ω(x) for each
input x, together with a metric on ω, but not on any target value for x, the
OP metric depends on a loss function ℓ(ω(x), y(x)) encoding the deviation
of ω(x) from a desired target y(x) for x; but not on a choice of metric for ω.
Definition 14 (Outer product metric). For each input x in the
dataset D, let ω(x) be the final output of the network run with input x
and parameter θ. Let ℓ(ω(x), y(x)) be the loss function measuring how ω(x)
departs from the desired output y(x) for x.
The outer product metric is defined as follows. Let δθ be a variation of
θ and let δℓx be the resulting variation of ℓ(ω(x), y(x)). Define
kδθk2op := Ex∈D (δℓx )2

(60)

In matrix form, this metric is
kδθk2op = Ex∼D δθ⊤

∂ℓ ⊤∂ℓ
δθ
∂θ ∂θ

(61)

∂ℓ
is the row vector of partial derivatives (the differential) of the loss
where ∂θ
function. Thus this metric is given by the matrix

Ex∼D

∂ℓ⊤∂ℓ
∂θ ∂θ

24

(62)

hence its name.
The unitwise outer product metric is defined by
kδθk2u-op :=

X

kδθk k2op

(63)

k

where k runs over the units of the network and δθk is the variation of the
incoming parameters to unit k. This metric is given by keeping only the
block-diagonal terms incoming to each unit in the matrix defining the tensorsquare differential metric.
The OP metric has been used simply under the name “natural gradient”
in [APF00, RMB07], which can lead to some confusion because it is distinct
from the natural metric using the true Fisher information matrix (see the
discussion in [PB13]). Moreover the OP metric makes sense for optimization
situations more general than the natural gradient, in which the loss function
is not necessarily of the form ln pθ for a probabilistic model p. For these two
reasons we adopt a purely descriptive name8 .
The OP metric is characterized, among all possible metrics, by a unique
property: it provides a gradient step for which progress is most evenly distributed among all data samples.
Proposition 15 (OP gradient equalizes the gain over the
samples). Let L := Ex∈D ℓx be the average loss with ℓx the loss on input x.
The direction δθ given by the gradient of L computed in the outer product
metric (Def. 26) has the following property: Among all directions δθ yielding
the same infinitesimal increment δL at first order, it is the one for which
the increment is most evenly spread over the data samples x ∈ D, namely,
Varx∈D δℓx = Ex∈D (δℓx − δL)2 is minimal.
The proof is given in the Appendix. The unitwise OP metric does not,
in general, satisfy this property: instead, it minimizes the variance, over a
random data sample x ∈ D and a random unit k in the network, of the
contribution to δℓx of the change δθk at unit k, so that it tries to spread δL
uniformly both over data samples and units.
Remark 16 (Metric for output parameters). The case when
the decoding function ω((ak )k∈Lout ) depends on additional “output parameP
ters” ϑ (e.g., softmax output with variable coefficients Pr(k) = eϑk ak /( k′ eϑk′ ak′ ))
can be recovered by considering ω as an additional output unit to the network, so that ϑ becomes the parameter of the activation function of ω. In
8

The outer product metric is distinct from the “outer product (or Levenberg–
Marquardt) approximation” of the Hessian of the loss function [Bis06, 5.4.2]. The latter
can be obtained from the natural metric (54) in which the output metric Ω is replaced
with the Hessian of the loss function ℓx w.r.t. ω. It depends on the parametrization of
ω. For exponential families in the canonical parametrization it coincides with the Fisher
metric.

25

particular, applying the definitions above, the metric Ω on ω induces a metric on ϑ by
kδϑk2bp

=

kδϑk2nat

⊤

= δϑ

∂ω⊤ ∂ω
Ω
Ex∈D
∂ϑ
∂ϑ

!

δϑ

(64)

⊤

∂ω
given by the matrix Ex∈D ∂ω
∂ϑ Ω ∂ϑ . So ϑ can be trained by gradient descent
in this metric.

For the Fisher metric, the following is well-known.
Proposition 17 (Invariance). The natural metric, unitwise natural
metric, backpropagated metric, and plain and unitwise outer product metrics are intrinsic: kδθknat , kδθku-nat , kδθkbp , kδθkop , and kδθku-op do not
depend on a choice of parametrization for the activations ak and for the
parameter θk at each unit k.
Proof. These metrics have been defined without any reference to parametrizations, directly by defining the norm of a tangent vector δθ; see Appendix C
for a more formal treatment. Consequently the value of the norm kδθk is
the same expressed in any coordinate system [GHL87, 2.1].
The natural metric actually has stronger invariance properties than the
unitwise natural metric: it does not depend on a change of parametrization
of the whole parameter θ = (θk ) that would mix the various components.
As such, the unitwise natural metric depends on a choice of decomposition
of the network into units, while the natural metric is only a function of the
input-output relationship of the whole network. The same holds for the OP
and unitwise OP metrics.
Remark 18 (Unitwise metrics as change in activation profile). We saw above that the metric used to define a gradient represents
a “cost” of moving in certain directions. All three unitwise metrics (unitwise natural, backpropagated, and unitwise OP) share a common property:
these metrics decompose as a sum, over the units k, of terms of the form
kδθk k2 = Ex∈D cx,k kδak (x)k2 where δak (x) is the resulting change of activity at k on input x, and cx,k is a weight (different for these three metrics)
estimating the influence of k on the output. Thus, the “cost” of a change at
unit k according to these metrics, is an average square norm of the resulting
change in activation profile ak (x) over x in the dataset. This is related to
the best-fit interpretation of these metrics (Section 3.3).

26

Computing the metrics. These metrics can be explicitly computed as
follows.
∂ℓ
The outer product metric is the easiest to compute: the terms ∂θ
are
∂ℓ ∂ak
∂ℓ
directly computed by ordinary backpropagation, namely, ∂θk = ∂a
k ∂θk

∂ℓ
k
where ∂a
(= bk ) is computed by backpropagation and ∂a
∂θk is obtained from
k
the activation function at unit k. Then the matrix defining the metric is
∂ℓ⊤ ∂ℓ
∂ℓ ⊤ ∂ℓ
Ex∈D ∂θ
∂θ , or, for the unitwise version, Ex∈D ∂θk ∂θk at each unit k.
To compute the natural and unitwise natural metrics, it is enough to
compute the Jacobian matrix ∂ω
∂θ . This can be done by performing one
backpropagation for each component of the output layer, for each input
x ∈ D, as follows.

Definition 19 (Backpropagation transfer rates). Let kout be
an output unit and let k be any unit in the network. The backpropagation
transfer rate Jkkout from kout to k is the dim(akout ) × dim(ak ) matrix defined
by

kout


Jkout := Iddim(akout )

J kout := 0

k


kout ∂aj
J kout := P
j, k→j Jj
k
∂a

where

∂aj
∂ak

k

for k 6= kout in the output layer Lout
for non-output units k

(65)

is the Jacobian matrix of the activation function from unit k to

unit j. Then we have Jkkout =

∂akout
∂ak .

This depends on an input x: the activation state of the network has to be
computed by forward propagation before these quantities can be computed.
Typically the activities are one-dimensional, not multidimensional, so
that each Jkkout is just a number, not a matrix. In this case, all the transfer
rates Jkkout can be computed by performing nout distinct backpropagations
each initialized with a single 1 on the output layer.
Since the influence of the parameter θk on the output goes through the
activity of unit k, the unitwise natural metric at k can be computed from a
single number (if activities are one-dimensional) measuring the influence of
unit k on the output, the Fisher modulus.
Definition 20 (Fisher modulus). Let x be an input. Let k be a
unit in the network. Let Ω be the metric on the final output ω. The Fisher
modulus Φk (x) of k on input x is the dim(ak ) × dim(ak ) matrix given by


Φk (x) := 

X

kout ∈Lout

⊤



∂ω kout   X
J
Ω
∂akout k
k ∈L
out

out



∂ω kout 
J
∂akout k

(66)

For each input x, the Fisher modulus is an intrinsic metric on ak : for a given
input x, the norm
kδak k2F-mod := δa⊤kΦk δak
(67)
27

does not depend on any choice of parametrization.
Note that ∂a∂ω
depends on the output layer interpretation but not on
kout
any parameter θ. Thus, since the transfer rates J can be computed by
backpropagation, the Fisher modulus only involves known quantities.
Proposition 21 (Computation of the unitwise natural metric). The unitwise natural metric at unit k is given by
kδθk k2u-nat = Ex∈D kδak (x)k2F-mod

(68)

⊤

∂ak
∂ak
δθk
Φk
= Ex∈D δθ⊤
k
∂θk
∂θk

(69)

k
where δak (x) is the variation of ak (x) induced by δθ, and ∂a
∂θk is the Jacobian
matrix of the activation function at k. Thus the matrix defining the unitwise
natural metric at unit k is

F (k) = Ex∈D

∂ak ⊤ ∂ak
Φk
∂θk
∂θk

Proof. By definition of the transfer rates J we have Jkkout =
P

kout ∈Lout

kout
∂ω
∂akout Jk

=

∂ω
∂ak

∂akout
∂ak .

Thus

so that
Φk =

hence

(70)

∂ω ⊤ ∂ω
Ω
∂ak
∂ak

∂ak ⊤ ∂ω ⊤ ∂ω ∂ak
∂ω ⊤ ∂ω
∂ak ⊤ ∂ak
=
=
Φk
Ω
Ω
∂θk
∂θk
∂θk ∂ak
∂ak ∂θk
∂θk
∂θk

(71)

(72)

which, after averaging over the dataset, is the definition of the unitwise
natural metric at k.
An analogous formula can be defined for the full (rather than unitwise)
Fisher matrix, by defining a Fisher modulus Φkk′ indexed by two units, and
using unit k′ on the left and k on the right in (66). Then the block of entries
of the Fisher matrix corresponding to parameters θk and θk′ is
Ex∈D

∂ak′ ⊤
∂ak
Φkk′
∂θk′
∂θk

(73)

(see also Proposition 27).
The unitwise Fisher metric is costly to compute when the output layer
is large. We can define another intrinsic metric for the activation of unit k,
simply by backpropagating the metric of the output layer.

28

The changes in the output induced by a change of θk all transit through
the activation of unit k. So if we have an intrinsic metric kδak k2 for the
activation of unit k, we can immediately define an intrinsic metric for θk , by
k
looking at the resulting change δak = ∂a
∂θk δθk induced by a change δθk , and
defining the norm of δθk to be the norm of the resulting δak . If the metric
on ak is given, in some parametrization of ak , by kδak k2 = δa⊤kgk δak where
gk is a symmetric, positive-definite matrix of size (dim ak ) × (dim ak ), then
defining kδθk k to be the norm of this δak ,
kδθk k :=

∂ak
δθk
∂θk

yields
2

kδθk k =



∂ak
δθk
∂θk

⊤

gk



(74)

∂ak
δθk
∂θk



(75)

⊤

∂ak
k
in other words, the matrix defining this metric is ∂a
∂θk gk ∂θk .
The unitwise Fisher metric is obtained from the Fisher modulus by this
construction. We now define another intrinsic modulus playing the same
role for the backpropagated metric.

Proposition 22 (Backpropagated modulus and computation of the backpropagated metric). Let x be an input. Let
k be a unit in the network. Let Ω be the metric on the final output ω. The
backpropagated modulus mk (x) at k is the dim(ak ) × dim(ak ) matrix given
by

 ∂ω ⊤Ω ∂ω
for k in the output layer
∂ak
∂ak
(76)
mk (x) := P
∂aj ⊤
∂aj

for k an internal unit
j, k→j ∂ak mj ∂ak
Then, for each input x, the backpropagated metric on ak is given by the
backpropagated modulus:
kδak k2bp,x = δa⊤kmk δak
and so the backpropagated metric on θk is given by the matrix Ex∈D
namely,
∂ak ⊤ ∂ak
kδθk k2bp = Ex∈D δθ⊤
δθk
mk
k
∂θk
∂θk

(77)
∂ak
∂ak ⊤
∂θk mk ∂θk ,

(78)

Proof. Immediate from the definition of the backpropagated metric and
∂ai
k
δak and δak = ∂a
δai = ∂a
∂θk δθk .
k
Like the Fisher modulus, the backpropagated modulus is a single number
when activities are one-dimensional. The cost of its computation is the same
as one backpropagation pass.
The equation defining the backpropagated modulus is well-known: it is
related to the so-called Gauss–Newton approximation to the Newton method
29

(see for instance [LBOM96], Section 7), which consists in computing the
Hessian of the loss function and throwing away all terms involving the second derivative of the activation function (those could result in non–positivedefinite terms, in which case the Newton method is ill-behaved), with the
additional approximation that cross-terms between different units are also
thrown away. Here we see that no approximation is involved: we do not
throw away annoying terms, we simply define an intrinsic metric. There
is actually no meaningful notion of the Hessian of a function on manifolds
[GHL87, paragraph 3.37] unless additional structure (affine, Riemannian) is
given or we are at a critical point of the function; the annoying terms above
are precisely the terms that prevent such a notion to exist. So one could
even say, in the context of differential geometry, that the Newton method is
an approximation of the backpropagated metric rather than the other way
round.
The backpropagated modulus and the Fisher modulus are related: If one
tries to write a backpropagated equation to compute the Fisher modulus Φk
in terms of the Fisher modulus at units pointed by k, one finds a quadratic
(instead of linear) backpropagation equation with terms involving pairs of
units. Keeping only the terms involving a single unit yields the equation
defining the backpropagated modulus.

2.3

Quasi-diagonal reduction of a unitwise metric

The unitwise Fisher metric, backpropagated metric, and unitwise OP metric,
still involve a full matrix on the incoming parameter space at each unit, and
are thus not adapted if network connectivity is large. We now introduce two
metrics enjoying lesser invariance properties than the above, but quicker to
compute.
Given an intrinsic metric kδθk k on θk (such as the unitwise Fisher or
backpropagated metric), we are going to define a simpler one, kδθk kqd . The
inverse of the matrix defining this metric will be quasi-diagonal, with the
only non-zero diagonal terms being those between a weight and the bias.
This will allow for quick gradient steps costing no more than classical backpropagation.
This relies on the affine structure in neural networks: this simplification
makes sense in a somewhat more restricted setting than the general setting
above. Suppose that the activation function
ak = fθkk (ai1 , . . . , aink )

(79)

can be written as a composition of a fixed, non-linear activation function ϕ
and a quantity yk that is an affine function of ai1 , . . . , aink :
ak = ϕ(yk,θk (ai1 , . . . , aink ))

30

(80)

such that when θk ranges over its values, yk,θk ranges over all possible affine
functions of ai1 , . . . , aink . For this to make sense, we now have to assume
that the activities ak live in an affine space. So let us go back to activities
with values in R, but without any preferred origin and scale for activities:
we look for invariance under replacement of ai with αi ai + βi and likewise
for yi .
In any given parametrization (choice of origin and basis) of ai and yk we
can write
X
wik ai + w0k
(81)
yk =
i, i→k

for some values wik ; specifying the parameter θk is equivalent to specifying
these quantities.
But this decomposition will change if we change the affine parametrization of activities: if a′i = αi ai + βi and yk′ = γk yk + δk the relation becomes
P
P
′
(a′i − βi ) + γk w0k = i (γk wik α−1
yk′ = δk + i γk wik α−1
i
i )ai + (γk w0k + δk −
P
−1
−1
′ =γ w α
new weights are wik
k ik i and the new bias
i wik αi βi ) so that the
P
−1
′
is w0k = γk w0k + δk − i wik αi βi . In particular we see that there is no
intrinsic “separation” between the bias and the weights; but that there is a
separation between wik and wi′ k for different incoming units i and i′ . This
is formalized as follows.
Let δθk be a change of parameter θk . For i 6= 0, let δwik be the resulting
′ be the resulting
change of wik in a given parametrization, and let δwik
′ = 0
change in another parametrization. If δwik = 0, then we have δwik
−1
′ = γ w α
as well in any other affine parametrization, since wik
yields
k ik i
−1
′
δwik = γk δwik αi . Note that this does not depend on the input x either9 .
Thus, having δwik = 0 is a property of δθk that does not depend on the
chosen affine parametrization of activities: it is an intrinsic property of the
change of parameter δθk . Say that a change of parameter δθk does not
involve unit i if δwik vanishes.
′ = γ w +δ −
For the bias the situation is different: the expression w0k
k 0k
k
P
−1
w
α
β
giving
the
bias
in
a
parametrization
from
the
bias
in
another
i
i ik i
parametrization is more complex, and so the fact that δw0k = 0 does depend
on the parametrization. This is where the metric kδθk k we are trying to
simplify comes into play.
Say that a change of parameter δθk is pure bias if it does not involve
any unit i incoming to k, i.e., if δwik = 0 for all i 6= 0. This is an intrinsic
condition. Say that δθk is bias-free if it is orthogonal, in the metric kδθk k we
are trying to simplify, to all pure-bias vectors. Being bias-free is an intrinsic
condition, because by assumption the metric kδθk k is intrinsic. Being biasfree does not simply mean δw0k = 0; let us work it out in coordinates.
9
This relies on the affine form (81) of yk . If yk is not affine in (80), having a constant
∂yk /∂ai is not a well-defined notion and may depend on the input.

31

Let Aii′ be the symmetric matrix defining the metric kδθk k in a given
parametrization. The associated scalar product is
hδθk , δθk′ i =

XX

Aii′ δwik δwi′ ′ k +

i′

i

X

′
′
′
A0i (δw0k δwik
+δw0k
δwik )+A00 δw0k δw0k

i

(82)

with A0i = Ai0 .
In particular, if the only non-zero component of δθk is δwik , then its
′ will be A δw′ δw . On the other
scalar product with a pure bias δw0k
0i
ik
0k
hand, if to δθk we add a bias component δw0k = −A−1
00 A0i δwik , then the
scalar product with any pure bias will vanish. Such a δθk is thus bias-free.
In the case when the parameter θk allows to represent all affine functions of the incoming activations, we can decompose a variation δθk of
θk into components δθki each involving only one incoming unit i, and a
pure bias component δθk0 . This decomposition is unique if we impose
that each δθki is bias-free. Explicitly, if in some parametrization we have
δθk = (δw0k , δw1k , . . . , δwnk k ) this decomposition is
δθki = (−A−1
00 A0i δwik , 0, . . . , δwik , 0, . . . , 0)

(83)

and
δθk0 = (δw0k +

X

A−1
00 A0i δwik , 0, . . . , 0)

(84)

i

The decomposition δθk = δθk0 + i δθki is intrinsic.
We can then define a new intrinsic metric on δθk by setting
P

kδθk k2qd := kδθk0 k2 +

X

kδθki k2

(85)

i

which is readily computed:
kδθk k2qd

= A00 δw0k +
+

X

X

i
2
(Aii δwik

A−1
00 A0i δwik

!2

−1
2
− 2g0i (A−1
00 A0i δwik )δwik + A00 (A00 A0i δwik ) )

i

2
= A00 δw0k
+2

X

A0i δw0k δwik +

+

A−1
00 A0i A0i′ δwik δwi′ k

i,i′

i

X

X

2
2
(Aii − A−1
00 A0i )δwik

i

(86)
Thus, this metric is defined by a matrix Ã given by Ã00 = A00 , Ã0i = A0i
−1 2
and Ãii′ = A−1
00 A0i A0i′ + 1i=i′ (Aii − A00 A0i ).
Definition 23. Quasi-diagonal reduction is the process which, to an intrinsic metric defined by a matrix A in affine coordinates, associates the
metric defined by the matrix
−1
Ã := diag(A) + A−1
00 (v ⊗ v) − diag(A00 (v ⊗ v))

32

(87)

where
v = (A00 , . . . , A0i , . . .)

(88)

The quasi-diagonal backpropagated metric is the quasi-diagonal metric
obtained from the backpropagated metric. The quasi-diagonal Fisher metric
is the one obtained from the unitwise Fisher metric. The quasi-diagonal OP
metric is the one obtained from the unitwise OP metric.
The reasoning in this section can be summarized as follows.
Proposition 24. Assume that the activation function is a fixed nonlinear function composed with an affine function. Then the quasi-diagonal
reduction Ã of an intrinsic metric A is intrinsic.
−1
Importantly, the matrix Ã = diag(A) + A−1
00 (v ⊗ v) − diag(A00 (v ⊗ v))
is the sum of a diagonal matrix and a rank-1 matrix. This allows for easy
inversion, resulting in a quasi-diagonal inverse matrix.

Proposition 25 (Quasi-diagonal gradient step). Let Ã be the
quasi-diagonal reduction of A. Let b = (b0 , . . . , bi , . . .) and w = Ã−1 b. Then
w is given by
bi A00 − b0 A0i
Aii A00 − A20i
X A0i
b0
−
wi
w0 =
A00 i6=0 A00
wi =

for i 6= 0

(89)
(90)

Thus, only the entries A00 , Aii and A0i of the original matrix A need to
be known in order to implement gradient descent using the quasi-diagonal
metric defined by Ã.
Note that if the original matrix A is positive-definite, we have A00 > 0
and A00 Aii > A20i (by the Cauchy–Schwarz inequality applied to the first
and i-th basis vectors), so that the solution w above is well-defined and
unique.

2.4

Intrinsic gradients

Thanks to these intrinsic metrics we can define intrinsic gradient directions
in parameter space. Given a dataset D of inputs x and corresponding targets
y, the average loss function is
Lθ := Ex∈D ℓθ (y)

(91)

where we put a subscript θ to make explicit its dependency on the parameters
of the network. Given an intrinsic metric k·k, the differential
G=−
33

∂Lθ
∂θ

(92)

of the average loss with respect to the full parameter set θ, defines a gradient
direction ∇θ L by the usual definition: it is the only tangent vector such that
for any δθ we have
Lθ+δθ = Lθ + h∇θ L, δθi + O(kδθk2 )

(93)

where h·, ·i is the scalar product associated with the norm k·k. In a parametrization where k·k2 is given by a symmetric, positive definite matrix A, the
gradient is given by
∂L
∇θ L = A−1
= −A−1 G
(94)
∂θ
The gradient ∇θ L is an intrinsic tangent vector on the parameter set.
Definition 26. The natural gradient, unitwise natural gradient, backpropagated metric gradient, OP gradient, unitwise OP gradient, and their
quasi-diagonal reductions, respectively, are the following update rule for θ:
θ ← θ − η∇θ L

(95)

where ∇θ L = A−1 ∂L
∂θ is the gradient of the average loss function L computed with A the natural metric, unitwise natural metric, backpropagated
metric, OP metric, unitwise OP metric, and their quasi-diagonal reductions,
respectively.
The algorithms of Section 1 are the application of these updates to ordinary neural networks, written out with [0; 1]-valued activities and sigmoid
activation function. More details on how this works out are given below
(Section 2.5).
This update is intrinsic only under all affine reparametrizations of the
parameter θ. Indeed, even if the tangent vector ∇θ L giving the direction
of the gradient is fully intrinsic, adding a tangent vector to a parameter θ
is not an intrinsic operation (if two parametrizations differ by a non-affine
transformation, then the additions will not amount to the same).
On the other hand, the ideal limit when the learning rate η tends to 0
is intrinsic: the trajectories of the differential equation
dθ(t)
= −∇θ(t) L
dt

(96)

are intrinsic trajectories in parameter space for the unitwise natural gradient
and backpropagated metric.
For the quasi-diagonal algorithms, invariance is always restricted to affine
reparametrizations, since this is the setup in which they are well-defined.

34

2.5

The Fisher matrix for neural networks

The general definitions above depend on the choice of a metric on the output
ω of the network. When this metric is the Fisher metric on the output layer,
applying the general definitions above to ordinary neural networks leads to
the algorithms of Section 1. This is mostly by direct computation and we
do not reproduce it fully. Let us however discuss the Fisher metric in more
detail.
For each input x, the network defines a probability distribution ω on the
outputs y. This probability distribution depends on the parameters of the
network. Thus, for each input x, we can define a datum-wise Fisher matrix
on the parameter set:
F (x)wij wi′ j′ = Ey|x

∂ℓ(y) ∂ℓ(y)
∂wij ∂wi′ j ′

(97)

where as above ℓ(y) = ln ω(y) and where Ey|x denotes expectation for y
following the distribution ω defined by the input x.
The dataset together with the network define a probability distribution
on pairs (x, y), by first choosing at random an input x in the dataset, then
running the network on this input. The Fisher matrix associated with this
distribution on pairs (x, y) is the average of the datum-wise Fisher matrix
over the dataset
F = Ex∈D F (x)
(98)
(see [AN00], Section 8.2), or more explicitly
Fwij wi′ j′ = Ex∈D Ey|x

∂ℓ(y) ∂ℓ(y)
∂wij ∂wi′ j ′

(99)

Exact Fisher matrix versus one-sample Fisher matrix. One possible
way to train neural networks using the natural gradient is to estimate the
Fisher matrix by taking an input x in the dataset, taking a random output
∂ℓ(y|x)
y for this input, and add the term ∂ℓ(y|x)
∂wij ∂wi′ j ′ to the current estimate of
the Fisher matrix (with a discount factor for older contributions in an online
setting). This leads to the Monte Carlo natural gradient with K samples y
per input x:
K
∂ℓ(yk ) ∂ℓ(yk )
1 X
(100)
F̂wij wi′ j′ = Ex∈D
K k=1 ∂wij ∂wi′ j ′
where each yk is drawn according to the output probability distribution ω
defined by the output for input x. Even K = 1 can lead to reasonable results
(Section 4).
An important variant uses for y the actual target value for input x,
instead of taking y as a random sample given by the activations of the

35

output layer:
F̂wij wi′ j′ = Ex∈D

∂ℓ(y(x)) ∂ℓ(y(x))
∂wij
∂wi′ j ′

(101)

with y(x) the target for x: this is just the outer product (OP) metric of
Definition 14. It is not an unbiased estimate of the Fisher metric; still,
hopefully when the network converges towards the desired targets, the law
of y taken from the output distribution converges to the actual target and
the two variants get close. This variant has been present for a long time in
studies on natural gradient (as is clear, e.g., from Equation (14) in [APF00])
and is elaborated upon in [RMB07]. As pointed out in [PB13], the two
variants are often confused.
These two variants are both intrinsic. The OP variant, contrary to the
true natural gradient, depends on the targets and not only on the network
and inputs.
Both the true natural gradient, its Monte Carlo approximation, and
its “one-sample”/OP variant give rise to a unitwise version and to a quasidiagonal version (Section 2.3). For a network with n units and at most d
incoming connections per unit, the algorithmic cost of processing each data
sample is O(Kn2 d2 ) for the Monte Carlo natural gradient, O(Knd2 ) for its
unitwise version and O(Knd) for its quasi-diagonal reduction. Algorithmic
cost for the OP metric is the same with K = 1.
In Section 4 we compare performance of the unitwise natural gradient,
Monte Carlo unitwise natural gradient with K = 1, and unitwise OP natural
gradient. We will see that although the OP metric and the one-sample
(K = 1) Monte Carlo natural gradient look similar, the latter can perform
substantially better.
Exact Fisher matrix computation. It is possible compute the exact
Fisher matrix (rather than using a single value for y) by using the Fisher
modulus and backpropagation transfer rates. The latter can be computed
by doing nout backpropagations for each input. This is of course more convenient than computing the expectation Ey|x by summing over the (in the
Bernoulli case) 2nout possible outcomes for y.
The backpropagation transfer rates Jkkout from Definition 2 simply implement the general Definition 19 for ordinary neural networks. In Section 1,
the unitwise natural gradient was obtained from these transfer rates through
the Fisher modulus. Here we reproduce the corresponding formula for all
terms of the Fisher matrix, not only the terms of the unitwise Fisher matrix
incoming to a given unit, so we introduce a Fisher modulus indexed by pairs
of units.
Proposition 27 (Exact Fisher matrix for neural networks).
Let x be an input for the network. Compute the transfer rates Jkkout as in

36

Definition 2. Depending on output layer interpretation, set for each pair of
units k and k′ :

P

Φkk′ (x) := kout ∈Lout Jkkout Jkk′out



k
k

P
Jk out Jk′out



′ (x) :=
Φ
kk

k
∈L
out
out
a
(1−a

kout
kout )

P

kout kout

1
a
k

Φkk′ (x) := S kout
Jk′  

∈Lout e out Jk

P

P

1
akout kout

−

S2

kout ∈Lout

e

Jk

(square-loss)
(Bernoulli)
(softmax)
kout ∈Lout

akout

e

Jkk′out

P


where S := kout ∈Lout eakout



P

kout kout


Φkk′ (x) := S4 kout
Jk′

∈Lout Jk


 P



kout
kout
4 P

J
a
J
a
−

′
k
k
2
k
∈L
out
k
∈L
out
k
k

out
out
out
out
S


P

2
:=

where S

kout ∈Lout



(spherical)

akout

(102)
Then the entry of the datum-wise Fisher matrix F (x) associated with
parameters wik and wjk′ (including biases with i = 0 or j = 0) is
F (x)wik wjk′ = ai aj rk rk′ Φkk′

(103)

and thus the corresponding entry in the Fisher matrix is
Fwik wjk′ = Ex∈D ai aj rk rk′ Φkk′

(104)

The proof is given in the Appendix and is a more or less straightforward
application of the results of the previous section, together with an explicit
computation of the Fisher metric on the output in the Bernoulli, square-loss,
and classification interpretations.
So it is possible to compute the full Fisher matrix by performing nout
independent backpropagations for each sample input. The Fisher matrix F ,
being the average of F (x) over the dataset, may be approximated by the
standard online or small-batch techniques using samples from the dataset.
For a network with only one hidden layer, this simplifies and no additional backpropagations are needed. Indeed, the backpropagation transfer
rates of the input layer are never used, and on the hidden layer are given by
Jkkout = wkkout rkout

(105)

from which the Fisher modulus can be immediately computed. This is the
case treated in [Kur94] (for the Bernoulli interpretation).

3
3.1

Some properties of unitwise algorithms and their
quasi-diagonal approximations
Performance improvement at each step

A common feature of all gradient-based algorithms in any metric is that
the objective function improves at each step provided the learning rate is
37

small enough. Consequently this holds for the unitwise natural gradient,
backpropagated metric gradient, and their quasi-diagonal reductions.
Proposition 28. Suppose that training has not reached a local optimum, i.e., that the gradient vector G of Section 1 does not vanish. Suppose
that the metric considered is non-degenerate (i.e., respectively, that the
matrices F (k) or M (k) are invertible, or that the denominators in the quasidiagonal algorithms do not vanish), so that the algorithms considered are
well-defined. Suppose that the chosen interpretation ω of the output layer
depends smoothly on the output layer activities.
Then there exists a value ηC of the learning rate such that, for any
learning rate η < ηC , the value of the loss function strictly decreases after
one step of the unitwise natural gradient, backpropagated metric gradient,
or their quasi-diagonal reductions.
As usual, the value of ηC depends on the current state of the network
and thus may change over the course of training.

3.2

Invariance properties

The algorithms presented in Section 1 are the implementation of the gradients and metrics defined in Section 2, written out using [0; 1]-valued activities
and the logistic activation function. We could have written them out, for instance, using [−1; 1]-valued activities and the tanh activation function, and
the learning trajectory would be the same—provided, of course, that the
initialization was done so that both implementations of the network behave
the same at startup. We present a more precise formulation of this property.
Imagine that the inputs of the network are subjected to simple transformations such as scaling (ai ← αi ai for i in the input layer) or 0/1 inversion
(ai ← 1 − ai ). There is a simple way to change the parameters of subsequent
units so that the final activation of the network stays the same, namely,
wij ← wij /αi for scaling and wij ← −wij , w0j ← w0j + wij for 0/1 inversion. So clearly the expressivity of a neural network is not sensitive to such
changes.
However, training will behave differently. For instance, if we apply one
step of backpropagation training to the scaled inputs with the scaled network, the coefficients of units which have been scaled down (αi < 1) will
evolve more slowly and conversely for αi > 1. The final output of the
network after the update will be different. (Hence the common practice
of rescaling the activities of units.) The same goes for 0/1 inversion in a
slightly more complicated way: evolution of the bias depends on the activity
of input units, and the weights from input units with activity close to 0 will
evolve faster than those with activity close to 1, as seen on Figure 1.
We would like the following invariance for a training procedure: If we
start with two networks N and N ′ which are fed inputs x and x′ with x′
38

obtained from a simple transformation of x, and if the parameters of N ′
are set such that initially its output is the same as N , then we would like
the outputs of N and N ′ to stay the same after one step of the training
procedure.
This is not satisfied by backpropagation. However, for any affine transform of the activities of any unit, this is satisfied by the natural gradient,
unitwise natural gradient, unitwise outer product gradient, backpropagated
metric gradient, and their quasi-diagonal reductions.
The sigmoid and tanh networks correspond to each other by the following
P
rewriting, thanks to tanh(x) = 2 sigm(2x) − 1: if ak = sigm( i→k wik ai +
P
′ a′ + w′ ) (and interpretation of the output
w0k ) and a′k = tanh( i→k wik
i
0k
layer in the tanh case is done by putting back the activities in [0; 1] via
a′ 7→ 1/2 + a′ /2), then the two networks will behave the same if we set
P
′
′ (i 6= 0) and w
′
wik = 4wik
0k = 2w0k − 2
i6=0 wik .
Definition 29. Let k be an input or internal unit. Call (α, β, γ)-affine
reparametrization of unit k the following operation: Replace the activation
of unit k
(106)
ak = fθkk (ai1 , . . . , aink )
where θk = (w0k , (wik )i→k ), with
a′k = αfγk θ′ (ai1 , . . . , aink ) + β
k

(107)

′ , (w′ )
′
where θk′ = (w0k
ik i→k ). Send ak instead of ak to the next layer of the
network, with weights modified as follows:
′
wkj
= wkj /α,

′
w0j
= w0j − wkj β/α

(108)

′ = w /γ for all units i with i → k (infor all units j such that k → j, and wik
ik
cluding i = 0), so that the final outputs before and after the reparametrization are the same.

The passage from sigm to tanh consists in applying the (2, −1, 2)-reparametrization
to all units. We have restricted the definition to non-output units to simplify
notation; for output units a corresponding reparametrization of the output
interpretation has to be done.
The following result is an immediate consequence of the intrinsic definition of the algorithms. It is only part of the invariance properties of the
objects from Section 2. In particular, in the limit of small learning rates
(η → 0), the trajectories (96) of the unitwise natural gradient, backpropagated metric gradient, and unitwise OP gradient, are invariant under all
smooth one-to-one reparametrization and not only affine ones.

39

Proposition 30 (Invariance under affine reparametrization of activities). Consider a network obtained from an initial network by applying any number of (α, β, γ)-affine reparametrizations to any
number of units (where α, β and γ may depend on the unit).
Then, after one step of the unitwise natural gradient, backpropagated
metric gradient, Monte Carlo unitwise natural gradient, unitwise OP gradient, or their quasi-diagonal reductions, the final outputs of the non-reparametrized
and reparametrized networks are the same.
Consequently, the learning trajectories, and performance, of the two networks with these corresponding initializations are the same.
This may seem a simple thing, but we should keep in mind that this
property is not satisfied by backpropagation, or by quasi-Newton methods
if the latter use diagonal approximations of the Hessian.
In particular, the algorithms presented here are insensitive to shifting
and scaling of all units in the network. Traditionally, it is recommended to
normalize the activities on input units so that they average to 0 over the
dataset and have a prescribed variance: the algorithms here automatically do
the same in an implicit way, for all (not only input) units. As a consequence,
units with low activation levels get updated as fast as highly activated units.
(Note that as discussed after the definition of the quasi-diagonal algorithms,
these averages and variances are computed according to non-uniform weights
on the dataset given by the Fisher modulus or backpropagated modulus.)
Still the invariance above only applies if the two networks considered have
corresponding initializations. For instance, if the initial weights are random
with a variance set to 1 whatever the data, obviously the initial behavior
of the network will be sensitive to scaling of its input. So these methods
do not remove the need for traditional recommendations for initializing the
weights (either by normalizing the data and then taking weights of size 1,
or by taking initial weights depending on the variance or covariance matrix
of the data).
The unitwise gradients (natural, backpropagated metric, OP, and Monte
Carlo natural), but not their quasi-diagonal reductions, have a further, more
interesting invariance property: invariance under affine recombination of the
signals a unit receives from its various incoming units. For instance, if we
start with zero weights, an internal unit will evolve in the same way if it
receives f and f + εg (where f and g are seen as function of the input x)
as if it receives f and g. This is especially useful if g is correlated to the
desired output.
Proposition 31 (Invariance under affine recombination of
incoming signals). Consider a neural network and define a new one
in the following way. Let k be a non-input unit in the network, with nk

40

incoming units, and let ϕ : Rk → Rk be an invertible affine map. Define a
new network by replacing the activation function at unit k
ak = fθkk (ai1 , . . . , aink )

(109)

ak = fϕk∗ (θk ) (ϕ(ai1 , . . . , aink ))

(110)

with
still parametrized by θk , where ϕ∗ (θk ) results from applying the dual10 inverse affine transformation ϕ∗ to θk , so that initially the responses of the
original and reparametrized networks are the same.
Then, after one step of the unitwise natural gradient, backpropagated
metric gradient, Monte Carlo natural gradient, or unitwise OP gradient,
with respect to θk , the final outputs of the non-reparametrized and reparametrized
networks are the same.
Consequently, the learning trajectories, and performance, of the two networks with these corresponding initializations are the same.
Once more, this is not simply ϕ in one place cancelling out ϕ−1 in another: indeed, backpropagation or quasi-Hessian methods do not have this
property, and neither do the quasi-diagonally-reduced algorithms.
Proof. This comes as a consequence of the best-fit interpretation (Proposition 32) below.
It also follows from the intrinsic constructions by noting that, unlike the
quasi-diagonal reductions, the construction of these gradients never breaks
down the nk -tuple of incoming activities into its components from each incoming unit; thus, contrary to the quasi-diagonal reductions, we could have
written the unitwise natural and backpropagated metrics in a setting where
activation functions are given by ak = fθkk (g(ai1 , . . . , aink )) where g is a fixed,
parameterless map with values in a manifold.

3.3

Best-fit interpretation

The unitwise natural gradient, backpropagated metric gradient, and unitwise OP gradient (but not their quasi-diagonal reductions) share an interpretation as a least-squares regression problem at each unit. Namely, the
backpropagated value bk (x, y) on input x and target y indicates how the
activity of unit k should change on input x. Seeing bk as a function of the
input x, unit k has to use the activities of incoming units i (also seen as
functions of x) and combine them using the weights wik , to match bk (x, y)
as close as possible for each x. This idea is presented in [Kur94] in a more
10

θk is an affine form over the nk -tuple of incoming activities. ϕ∗ is defined, axiomatically, by the property that applying ϕ∗ (θk ) to the activities transformed by ϕ, is the same
as applying θk to the untransformed activities. Decomposing θk = (w0k , (wik )i→k ), the
affine matrix defining ϕ∗ is the transpose of the inverse of the affine matrix defining ϕ.

41

specific setting. This is also relevant to the behavior of the algorithms when
the matrices F and M defining the metrics are singular or close to singular,
as we will see.
Proposition 32 (Intrinsic gradients as best fit to b). Let k
be a non-input unit in the network. For x in the dataset D, let bk (x) be the
backpropagated value (4) obtained on input x and the corresponding target
y.
Consider the solution λ = (λi ) to the following weighted least-squares
problem:
λ = arg min


X


x∈D

X
i

bk (x)
λi ai (x) −
rk (x)Φk (x)

!2

Wx




(111)



where i runs over the incoming units to k (including i = 0 with a0 ≡ 1),
Φk (x) is the Fisher modulus (Definition 3), and the weights are
Wx := rk (x)2 Φk (x)

(112)

Then the unitwise natural gradient step (21) is given by λ, namely, the
update is wik ← wik + ηλi at each unit k.
The same holds for the backpropagated metric gradient using the backpropagated modulus mk (Definition 7) instead of the Fisher modulus Φk .
The same holds for the unitwise OP gradient using bk (x)2 instead of the
Fisher modulus Φk .
Thus, the gradient step depends on the linear span of the incoming activities (ai (x))i→k , seen as functions over the dataset (which, by the way, proves
Proposition 31 above). This is why the gradient step is the same whether
the unit receives signals f (x) and g(x) or f (x) and f (x)+εg(x). Thus, these
algorithms perform an implicit orthonormalization of the incoming signals
at each unit (not only input units) in the network.
Proof. A direct application of the well-known formula for the solution of the
weighted least-squares problem (111), with the choice of weight (112), yields
exactly the updates (20) and (30).
Non-invertibility and regularization of the matrices. In several situations the matrices F and M used to define the unitwise natural update
and backpropagated metric update can be singular.
This is the case, for instance, if one input unit is uniformly set to 0 over
all elements in the dataset: obviously such a unit is not informative, and
the corresponding term will vanish both in the metric and in the gradient.
This is also the case when, e.g., two units incoming to the same unit are

42

perfectly correlated. Correlation in the activation profiles happens systematically in case the size of the dataset is smaller than the number of incoming
parameters at a given unit.
The linear regression viewpoint limits, in theory, the seriousness of these
issues: this only means the linear regression problem has several solutions
(one can add any quantity to a non-informative weight), and any of them
will do as an update. Indeed, for instance, the matrix M in Definition 9
is of the form AA⊤, and M −1 is applied to the vector G which is of the
form AY , thus G always lies in the image of A and thus the linear system is
underdetermined, not overdetermined. From the gradient ascent viewpoint
this means the matrix M will be singular but the gradient term ∂L/∂w will
vanish in the corresponding directions.
Numerically, however, the issue must be dealt with. One can use the
Moore–Penrose pseudoinverse of M (or F in the Fisher matrix case), obtained by adding ε. Id to M or to F with very small ε. This is a standard
regularization technique. It has the advantage of producing a well-defined
update when ε → 0, asymptotically independent of ε.
Thus, if a formal definition is needed, one can decide to use the Moore–
Penrose pseudoinverse for M −1 and F −1 in the definition of the updates.
However, this formally breaks the invariance properties: the Moore–Penrose
pseudoinverse selects, among the several possible solutions (λi ) to an underP 2
determined least squares problem, the one with smallest norm
λi , and
this is not intrinsic.

4

A first experimental comparison

Although the main focus of this article is theoretical, we performed a light
set of experiments to ensure that the suggested methods are viable. The
companion article [Oll13] contains more in-depth experiments with recurrent
networks and complex symbolic sequences.
To test the influence of the different methods, we chose a very simple
problem in which a perfect solution is expected to be found. A sparsely
connected network with 5 layers of size 100, 30, 10, 30, and 100 was built,
and 16 random length-100 binary strings were fed to the input layer, with
the target equal to the input (auto-encoding). Ideally the network learns to
encode each of the 16 samples using 4 bits on the middle layer (thus with
room to spare) and uses the bottom layer parameters to rewrite the output
from this. This is purely an optimization task without a learning aspect, as
there is no generalization to be done and no underlying pattern. The focus
is on which methods are able to converge to a good solution.
The sparsely connected network is built at random in each instance as
follows. Each of the 100 units in the input layer is linked to 5 randomly
selected nodes in the first hidden layer. Each of the 30 units in the first

43

hidden layer is linked to 5 random nodes in the middle hidden layer. The
scheme is reversed for the bottom part of the model: each of the 100 output
units is linked to 5 random nodes in the last hidden layer, and each unit in
the last hidden layer is linked to 5 random nodes of the middle hidden layer.
For each instance, the dataset is made of 16 random binary strings of
length 100. The target for each input is identical to the input. We use
Bernoulli interpretation of the output.
Note that this setting is adverse for the unitwise and quasi-diagonal natural gradients, which require a small output layer; this must be remembered
in the comparisons below.
To test the influence of parametrization on non-invariant algorithms, and
to check invariance of the invariant ones, each algorithm was implemented
both using sigm(x) and tanh(x) as the activation function.
The methods tested are: backpropagation; unitwise natural gradient;
quasi-diagonal natural gradient; unitwise OP gradient; Monte Carlo unitwise or quasi-diagonal natural gradient with one sample (K = 1 in (100));
backpropagated metric gradient; quasi-diagonal backpropagated metric gradient; diagonal Gauss–Newton ([LBOM96, SZL13]; equivalent to keeping
only the diagonal terms in the quasi-diagonal backpropagated metric gradient); and a batch version of Adagrad/RMSprop [DHS11] in which the learning rate for each gradient component is divided by the root mean square of
this component over the samples.
Since the sample size is small, the algorithms were run in batch mode.
Regularization. The algorithms were taken directly from Section 1. To
all methods except backpropagation, we added a regularization term of
10−4 Id to the various matrices involved, to stabilize numerical inversion.
This value is not so small; values such as 10−7 seemed to affect performance.
This is probably due to the small sample size (16 samples): each sample contributes a rank-1 matrix to the various metrics. Larger sample sizes would
probably need less regularization.
Initialization. For the tanh activation function, all the weights werepinitialized to a centered Gaussian random variable of standard deviation 1/ dj
with dj the number of units pointing to unit j, and the biases set to 0. For
the sigmoid activation function, the initialization was the corresponding one
(using Eqs. 46 and 47) so that initially the responses of the networks are
the same: namely, each weightpwas set to a centered Gaussian random variable of standard deviation 4/ dj , and then the bias at unit k was set to
P
− i i→k wik /2. This initialization has the property that if the incoming
signals to a unit are independent, centered about 1/2 (sigmoid) or 0 (tanh)
and of variance σ with σ not too large, then the output of the unit is also
centered of variance ≈ σ. (The factor 4 in the sigmoid case compensates

44

for the derivative 1/4 of the sigmoid function at 0.) See the argument in
[GB10]11 .
Learning rate. A simple adaptive method was used for the learning rate.
All methods based on gradients in a metric have a guarantee of improvement
at each step if the learning rate is small enough. So in the implementation,
if a step was found to make the loss function worse (in a batch mode, thus
summed over all samples), the step was cancelled and the learning rate
was divided by 2. If the step improves the loss function, the learning rate
is increased by a factor 1.1. The initial learning rate was set to 0.01; in
practice the initial value of the learning rate is quickly forgotten and has
little influence.
Unfortunately this scheme only makes sense in batch mode, but it has the
advantage of automatically selecting learning rates that suit each method,
thus placing all methods on an equal footing.
Execution time and number of iterations. First, 10, 000 steps of backpropagation were performed on the whole dataset, in batch mode. The resulting running time was set aside and converted to an equivalent number
of iterations for all of the other algorithms. This is a very rough way to proceed, since the running times can depend on the implementation details12 ,
and vary from run to run (because floating point operations do not take the
same time depending on the numbers they are operating on, especially when
both very small and very large values are involved).
Most of all, the different methods scale in different ways with the network,
and so the network used here may not be representative of other situations.
In particular this auto-encoder setting with 100 output units puts the unitwise natural gradient and quasi-diagonal natural gradient at a disadvantage
(on the same time budget they must performe nout backpropagations per
sample), compared to, e.g., a classification setting.
Nevertheless, we give in Table 4 the numbers of iterations giving roughly
equal running time for each method.
The natural gradient was also tested (using the exact full Fisher matrix
as obtained from Proposition 27). The computational cost is very high and
only 10 iterations take place in the alloted time, too few for convergence.
Thus we do not report the associated results.
Results. In Table 2, we report the average loss per sample, in bits, and its
standard deviation, at the end of the allocated number of training iterations.
These values can be interpreted as representing the average number of bits
11

p

The other initialization suggested in [GB10], with weights of magnitude 6/(dj + d′j )
with d′j the number of edges from j, did not make any significant difference in our setup.
12
We tried to implement each method equally carefully.

45

Method
Backpropagation (sigmoid)
Backpropagation (tanh)
Natural gradient
Unitwise natural gradient
Quasi-diagonal natural gradient
Backpropagated metric
Quasi-diagonal backpropagated metric
Monte Carlo unitwise natural gradient (K = 1)
Monte Carlo quasi-diagonal natural gradient (K = 1)
Unitwise OP gradient
Diagonal Gauss–Newton
AdaGrad

Number of iterations
10,000
10,000
9 to 10
2,100 to 2,300
2,800 to 3,100
4,200 to 4,300
7,400 to 7,500
2,800 to 2,900
3,900 to 4,000
4,000 to 4,100
7,700 to 7,800
8,000 to 8,100

Table 1: Number of iterations resulting in approximately equal execution
times for our problem
of an output that the model did not learn to predict correctly (out of 100).
The results of the implementations using sigmoid and tanh activation are
reported separately.
Performance as a function of time is plotted in Figure 3.
The statistics were made using 20 independent runs for each method.
Interpretation. These results are mainly illustrative: the situation considered here may not be representative because of the small sample size and
network size involved.
Still, it is clear that for problems of this size, the more elaborate algorithms are very competitive. Only the tanh implementation of the diagonal
Gauss–Newton method comes close to the invariant algorithms (while its
performance in sigmoid implementation is not as good).
As can be expected, the invariant algorithms have similar performance
in sigmoid or tanh implementation: trajectories match each other closely
(Figure 4). The variations are caused, first, by random initialization of the
dataset and weights in each run, and second, by the inclusion of the regularization terms ε Id, which breaks invariance. If the effect of the latter is
isolated, by having the same initialization in tanh and sigmoid implementations, the trajectories coincide for the first few iterations but then start to
differ, without affecting overall performance.
In this setting, the natural gradient, in its unitwise and quasi-diagonal
versions, seems to perform slightly worse than the backpropagated metric
methods. This might be an effect of the large output layer size (which
directly affects their computational complexity) combined with a given com-

46

80
sigm backprop
tanh backprop
tanh unitwise nat grad
tanh backpropagated metric
tanh quasi-diag nat grad
tanh quasi-diag backpr metric
sigm diag Gauss-Newton
tanh diag Gauss-Newton
sigm UOP grad
tanh UOP grad
sigm adagrad
tanh adagrad

70

60

Loss (bits)

50

40

30

20

10

0
Computation time (arbitrary units)

Figure 3: Performance over time of all algorithms involved. For better
readability the trajectories of the invariant algorithms have been plotted only
in tanh implementation (Figs. 2 and 4 show them in sigmoid implementation
for completeness). Average over 20 runs.

47

80

70

60

60

50

50

Loss (bits)

Loss (bits)

70

80
sigm unitwise nat grad
tanh unitwise nat grad

40
30

40
30

20

20

10

10

0

0
Computation time (arbitrary units)

Computation time (arbitrary units)

80

80
sigm backpr metric
tanh backpr metric

70

60

60

50

50

Loss (bits)

Loss (bits)

70

sigm QD nat grad
tanh QD nat grad

40
30

40
30

20

20

10

10

0

0
Computation time (arbitrary units)

sigm QD backpr metric
tanh QD backpr metric

Computation time (arbitrary units)

Figure 4: Double-checking invariance: Comparison of the trajectories of the
invariant algorithms in tanh and sigmoid implementations

48

Average loss (bits) ± std-dev
sigmoid
tanh

Method
Non-invariant:
Backpropagation
Diagonal Gauss–Newton
AdaGrad
Invariant:
Unitwise natural gradient
Quasi-diagonal natural gradient
Backpropagated metric
Quasi-diagonal backpropagated metric
Monte Carlo unitwise natural gradient
Monte Carlo quasi-diagonal natural gradient
Unitwise OP gradient

35.9 ± 2.1
11.6 ± 2.5
51.1 ± 3.3

24.7 ± 2.2
3.5 ± 2.0
25.3 ± 2.0

0.9 ± 1
3.5 ± 1.2
0.8 ± 0.8
1.9 ± 1.2
12.9 ± 1.5
7.9 ± 2.3
24.7 ± 2.5

1.4 ± 1.8
3.4 ± 1.6
0.3 ± 0.5
1.5 ± 1.3
14.1 ± 2.2
10.1 ± 2.5
28.5 ± 3.4

Table 2: Average loss per sample (bits) after an execution time equivalent
to 10, 000 backpropagation passes, computed over 20 independent runs, together with standard deviation over the runs
putation time budget. Per iteration instead of computation time, the natural
gradient methods perform better than the backpropagated metric methods,
and we expect them to be more competitive for smaller output sizes, e.g.,
in classification tasks.
The two algorithms based on squared gradients, Adagrad and the unitwise OP gradient, both perform rather poorly in this setting. Adagrad
differs from the unitwise OP gradient by using a diagonal approximation
and introducing square roots, which breaks invariance, while the unitwise
OP gradient is invariant and is meant to approximate the natural gradient.
This is a surprise, as, for instance, methods close to the OP gradient have
been found to perform well, e.g., in [RMB07], or in [Oll13] for recurrent networks. The small size of the dataset in our setting is not enough to explain
this problem, as it does not seem to affect the other methods. This may
be related to the large dimensionality of the output layer compared to the
number of samples in our experiment (in contrast to [RMB07] or [Oll13]),
which damages the one-sample OP metric approximation of the natural gradient and could result in low-rank OP matrices. Indeed, reasoning on the
full (whole-network) metric, the OP gradient contributes a matrix of rank 1
for each data sample (see (101)); on the other hand, the exact Fisher matrix
contributes a sum of nout matrices of rank 1 for each data sample as can be
seen from (102)–(103). Thus from a theoretical viewpoint the quality of the
one-sample OP approximation of the natural gradient is likely to depend on
output dimensionality.

49

80
70

sigm diag Gauss-Newton
tanh diag Gauss-Newton
sigm quasi-diag BPM
tanh quasi-diag BPM

Loss (bits)

60
50
40
30
20
10
0
Computation time (arbitrary units)

Figure 5: Effect of introducing a few non-diagonal terms: Comparison of
the diagonal Gauss–Newton and the quasi-diagonal backpropagated metric
methods

Lack of invariance is striking for some algorithms, such as the diagonal
Gauss–Newton method: its performance is very different in the sigmoid and
tanh interpretations (Figure 5). The quasi-diagonal backpropagated metric
method only differs from diagonal Gauss–Newton by the inclusion of a small
number of non-diagonal terms in the update. This change brings the sigmoid
and tanh implementations in line with each other, improving performance
with respect to the best of the two diagonal Gauss–Newton implementations.
In settings where the activities of units (especially internal units, since the
input can always be centered) are not as well centered as here, we expect
the quasi-diagonal backpropagated metric method to outperform the tanh
diagonal Gauss–Newton implementation even more clearly. Thus the quasidiagonal backpropagated metric is arguably “the invariant way” to write the
diagonal Gauss–Newton algorithm.
In conclusion, although this experiment is a small-scale one, it clearly
emphasizes the interest of using invariant algorithms.

Conclusions
• It is possible to write invariant training algorithms for neural networks
that do not have the scalability issues of the natural gradient. For a
50

network with n units, nout output units, and at most d incoming connections per unit, we have provided four invariant algorithms for which
the cost of processing a data sample is respectively O(nd), O(nd2 ),
O(ndnout ) and O(nd2 + ndnout ). The slower methods are closer to the
natural gradient and have stronger invariance properties.
• All of these algorithms are mathematically motivated by building Riemannian metrics on the activity and parameter space of neural networks, treated as manifolds.
• The outer product metric encountered in the literature is also naturally
interpreted in this framework. It has a unique property of spreading
improvement most uniformly across the samples at each step.
• In a small-scale experiment involving an auto-encoding task with three
hidden layers, invariant methods substantially outperform non-invariant
methods. Lack of invariance is clear for some commonly used methods.
• The quasi-diagonal backpropagated metric is close to the diagonal
Gauss–Newton method [LBOM96, SZL13] but crucially differs by the
inclusion of a few well-chosen non-diagonal terms. While performance
of the Gauss–Newton method can change substantially from sigmoid
to tanh implementation, the quasi-diagonal backpropagated metric is
invariant and improves performance with respect to diagonal Gauss–
Newton in any implementation.
Acknowledgments. I would like to thank Youhei Akimoto, Jérémy Bensadon, Cécile Germain, Balázs Kégl, and Michèle Sebag for helpful conversations and their feedback on this article, as well as the anonymous referees
for careful reading and helpful suggestions.

Appendix
A

Proof of Proposition 27

The Fisher metric depends, of course, on the interpretation of the output
layer as a probability distribution ω. The final output ω is a probability
distribution over the values of y in the target space, parametrized by the
activities ak of the output units k. If the output activities change by δak ,
the probability distribution ω will change as well. The norm of this change
in Fisher metric is
kδωk2nat = Ey∼ω (δ ln ω(y))2
X
∂ ln ω(y) ∂ ln ω(y)
δak δak′
=
Ey∼ω
∂ak
∂ak′
k,k ′ ∈L
out

51

(113)
(114)

thus stemming from the matrix
Fkk′ := Ey∼ω

∂ ln ω(y) ∂ ln ω(y)
∂ak
∂ak′

(115)

over the output layer.
In the Bernoulli interpretation, for each component k of the output layer,
the random variable yk is a Bernoulli variable with parameter ak . In the
square-loss interpretation, each yk is a Gaussian random variable with mean
ak and variance 1. In the two classification interpretations, y is a discrete
P
random variable which takes value k with probability ωk = eak /( j∈Lout eaj )
P
or ωk = a2k /( j∈Lout a2j ).
Let us compute the Fisher metric in the space ω in each case. In the
Q
Bernoulli case, we have ω(y) = k∈Lout (1yk =1 ak + 1yk =0 (1 − ak )). Consequently
yi − ak
1y =1 1y =0
∂ ln ω(y)
= k − k
=
(116)
∂ak
ak
1 − ak
ak (1 − ak )
Since under the distribution ω we have Eyk = ak and Var yk = ak (1 − ak )
(with yk and yj independent for k 6= j) we find
Fkk′ =

1k=k′

(117)

ak (1 − ak )

for k and k′ in the output layer.

−(y −ω )2 /2

so that ∂ ln∂ωω(y)
In the Gaussian case we have ω(y) = k e i√2πk
=
k
yi − ωk . Since under the distribution ω we have Eyk = ωk and Var yk = 1
we find F = Id hence
(118)
Fkk′ = 1k=k′
Q

for k and k′ in the output layer.
In the softmax case the probability to have y = j is ω(y) = eaj /S with
P
S = i∈Lout eai . Thus ∂ ln ω(y)/∂ak = 1k=j −eak /S. Taking the expectation
over y (i.e., over j) we find
Fkk′ =

eak
1j=k −
S

X eaj 

j∈Lout

S

eak′
1j=k′ −
S







(119)


eak eak′
eak′ eak X eaj  eak eak′
eak
1k=k′ −
−
+
=
S
S S
S S
S
S S
j
=

eak
eak eak′
1k=k′ −
S
S2

(120)
(121)

Similarly, in the spherical case the probability to have y = j is ω(y) =
P
1
ak
with S = i∈Lout a2i . Thus ∂ ln ω(y)/∂ak = 2( k=j
aj − S ). Taking the

a2j /S

52

expectation over y we find
Fkk′ = 4

X a2j

j∈Lout

=4

1j=k

S

aj

a2k

ak
−
S

1j=k′

!

aj

a2k

ak′
−
S

!

(122)


a2k′

1 1
1 a
1 ak
−4
+4
1k=k′ − 4
S ak ak′
S ak S
S ak′ S
k′

X a2j
j

4
4ak ak′
= 1k=k′ −
S
S2

S




ak ak′
S S
(123)
(124)

These give the expression of the Fisher matrix Fkk′ for k and k′ in the
output layer. This is enough to compute the full Fisher matrix, as follows.
Let x be an input for the network. Given a variation δθ of the network
parameters θ, let δai be the resulting variation of unit i, and let δω be the
P
∂ω
resulting variation of the final output ω. We have δω = k∈Lout ∂a
δak .
k
The datum-wise Fisher metric on θ is
kδθk2nat = kδωk2nat
=

X

(125)
(126)

Fkk′ δak δak′

k,k ′ ∈Lout
k
For each k in the output layer, we have δak = i ∂a
∂θi δθi where the sum
∂ak ∂ai
k ∂ai
k
runs over all units i in the network. For each i we have ∂a
∂θi = ∂ai ∂θi = Ji ∂θi .
Plugging this into the above yields

P

kδθk2nat =

XX X
i

X

i′ k∈Lout k ′ ∈Lout

Fkk′ Jik Jik′

′

∂ai ∂ai′
∂θi ∂θi′

(127)

so that the term of the Fisher matrix corresponding to δθi and δθi′ is
P
k k ′ ∂ai ∂ai′
k ′ ∈Lout Fkk ′ Ji Ji′ ∂θi ∂θi′ .
k∈Lout
For standard neural networks we have δθi = (δwji )j, j→i and moreover
∂ai
∂wji = aj ri .
Plugging into this the expression for Fkk′ yields the results in Proposition 27.
P

B

Proof of Proposition 15

Let v be an infinitesimal variation of the parameter θ. Let vi be the coordinates of v in some coordinate system. At first order, the increment in the
P x
average loss function along v is Ex∈D i ∂ℓ
∂θi vi .
∂ℓx
Let us abbreviate ∂i ℓx = ∂θi . The matrix defining the OP metric is
Mij = Ex∈D ∂i ℓx ∂j ℓx . The corresponding gradient direction is M −1 Ex∈D ∂ℓx .
P
Let m = Ex∈D i ∂i ℓx vi be the change in loss function associated with v.
The variance, when x runs over the dataset, of the gain in the loss function
53

for x is Ex∈D (( i ∂i ℓx vi ) − m)2 = Ex∈D ( i ∂i ℓx vi )2 −m2 . For fixed average
P
gain m, this is minimal when Ex∈D ( i ∂i ℓx vi )2 is minimal.
This is a smooth convex function of v, whose minimum we have to find
P
over the hyperplane {v, Ex∈D i ∂i ℓx vi = m}. The minimum of a positiveP
P
definite quadratic functional ij Aij vi vj over a hyperplane i Bi vi = m, is
found at v = λA−1 B for some constant λ. Here we have Bi = Ex∈D ∂i ℓx , and
P
P
P
P
expanding Ex∈D ( i ∂i ℓx vi )2 = Ex∈D (( i ∂i ℓx vi )( j ∂j ℓx vj )) = ij Ex∈D ∂i ℓx ∂j ℓx vi vj
yields Aij = Ex∈D ∂i ℓx ∂j ℓx = Mij . Consequently, for any value of m, the
optimal v is a multiple of the OP gradient direction M −1 Ex∈D ∂ℓx .
P

P

C

Definition of the metrics in the formalism of differential
geometry

Let L be the neural network (directed acyclic finite graph of units); for
k ∈ L, let the activities of unit k belong to a manifold Ak . Let the actiQ
vation function for unit k be fk : (Θk × i→k Ai ) → Ak , (θk , (ai )i→k ) 7→
fθk ((ai )i→k ) where Θk is the manifold of the parameters of unit k. Let
Lin ⊂ L and Lout ⊂ L be the input and output layers, respectively; let X
Q
be the space to which the inputs belong, and let ι : X → k∈Lin Ak be the
input encoding. Let O be the manifold to which the outputs belong, and
Q
let ω : k∈Lout Ak → O be the output interpretation.
The values of ak and of the output ω can be seen as functions of the
Q
parameter θ = (θk ) ∈ k Θk and the input x, by using the induction relations defined by the network: ak (θ, x) := ι(x)k for k ∈ Lin , ak (θ, x) :=
fk (θk , (ai (θ, x))i→k ) for k 6∈ Lin , and by abuse of notation, ω(θ, x) := ω((ak (θ, x))k∈Lout ).
Let Ω be the output metric: a Riemannian metric on the output manifold
O, which to every vector δω tangent to O at point ω ∈ O, associates its
square norm Ωω (δω, δω) in a bilinear way. An important example is the
Fisher metric when O is a space of probability distributions.
Q
Let Θ := k Θk be the parameter manifold. We are going to define
the natural metric, unitwise natural metric, and backpropagated metric as
Riemannian metrics on Θ.
If ϕ : E → F is a linear map between vector spaces E and F , and g is a
bilinear form on F , we define the bilinear form g ◦ ϕ on E by
g ◦ ϕ : (e, e′ ) 7→ g(ϕ(e), ϕ(e′ ))

(128)

for any two vectors e, e′ ∈ E. If g is positive-semidefinite then so is g ◦ ϕ.
We denote by Tp M the tangent space to a manifold M at a point p ∈ M .
Recall [GHL87, 1.36] that if h : M → M ′ is a smooth map between manifolds,
′
its differential ∂h
∂p (p) at point p ∈ M is a linear map from Tp M to Th(p) M .
Let θ ∈ Θ. For an input x, let ∂ω
∂θ (θ, x) be the differential of the network
output ω(θ, x) with respect to θ: this is a linear map from Tθ Θ to Tω(θ,x) O.

54

Define the natural metric gnat as the bilinear form on Tθ Θ for each θ
given by


1 X
∂ω
gnat :=
(θ, x)
(129)
Ωω(θ,x) ◦
#D x∈D
∂θ

where x ranges over inputs in the dataset D. By construction, this metric
does not depend on any choice of parametrization and is thus intrinsic.
The unitwise natural metric is defined in a likewise manner except that
it first breaks down the tangent vector δθ ∈ Tθ Θ into its components along
L
Q
each unit k using that Θ = k Θk and thus Tθ Θ = k Tθk Θk . The effect is
to make the components δθk orthogonal. Namely:
gu-nat (δθ, δθ) :=

X

k
(δθk , δθk )
gnat

(130)

k

where δθ =

L

k

δθk , and where
k
gnat



1 X
∂ω
:=
(θ, x)
Ωω(θ,x) ◦
#D x∈D
∂θk

(131)

∂ω
(θ, x) the differential of the network
is the natural metric on Θk , with ∂θ
k
output with respect to θk , which is a linear map from Tθk Θk to Tω(θ,x) O.
The backpropagated metric is defined by backward induction in the directed acyclic graph L. First, for each input x and each unit k, let us define
Ak
a bilinear form gbp,x
on the tangent space Tak (θ,x) Ak to the activity at k.
On the output layer let us set
Ak
:= Ωω(θ,x) ◦
gbp,x

∂ω
((aj (θ, x))j∈Lout )
∂ak

for k ∈ Lout

(132)

∂ω
where ∂a
((aj (θ, x))j∈Lout ) is the differential of the output interpretation
k
Q
function ω : j∈Lout Aj → O with respect to ak , which is a linear map from
Tak (θ,x) Ak to Tω(θ,x) O. Then this is backpropagated through the network:
for each k we define a bilinear form on Tak (θ,x) Ak by
Ak
:=
gbp,x

X

Ai
gbp,x
◦

i, k→i

∂fi
(θi , (aj (θ, x))j→i )
∂ak

for k 6∈ Lout

(133)

with fi : Θi × j→i Aj → Ai the activation function of unit i. (If a unit
is both an output unit and influences some other units, we add the two
contributions.) This is transferred to a metric on Θk via
Q

Ak
Θk
:= gbp,x
◦
gbp,x

Finally, letting again δθ =
backpropagated metric by

L

gbp (δθ, δθ) :=

k

∂fk
(θk , (aj (θ, x))j→k )
∂θk

(134)

δθk be a tangent vector to Θ, define the

1 X X Θk
(δθk , δθk )
g
#D x∈D k bp,x
55

(135)

which is a metric on Θ.
Note that these metrics may be non-positive definite (e.g., if a parameter
has no influence on the output).
Since these metrics have been defined using only intrinsic objects without choosing a parametrization of any of the manifolds Θk , they are intrinsic
(for a given output metric Ω). Consequently, when working in explicit coordinates, the value of the norm of δθ is invariant with respect to any change
of variables for each θk (diffeomorphism of Θk ). The natural metric has the
additional property that it is invariant under changes of variables mixing the
Q
parameters of various units: its invariance group is Diff( k Θk ) whereas the
Q
invariance group is the smaller group k Diff(Θk ) for the unitwise natural
metric and backpropagated metric.

References
[AGS05]

Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient
flows in metric spaces and in the space of probability measures.
Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel,
2005.

[Ama98]

Shun-ichi Amari. Natural gradient works efficiently in learning.
Neural Comput., 10:251–276, February 1998.

[AN00]

Shun-ichi Amari and Hiroshi Nagaoka. Methods of information
geometry, volume 191 of Translations of Mathematical Monographs. American Mathematical Society, Providence, RI, 2000.
Translated from the 1993 Japanese original by Daishi Harada.

[APF00]

Shun-ichi Amari, Hyeyoung Park, and Kenji Fukumizu. Adaptive method of realizing natural gradient learning for multilayer
perceptrons. Neural Computation, 12(6):1399–1409, 2000.

[Bis06]

Christopher M. Bishop. Pattern recognition and machine learning. Springer, 2006.

[BL88]

Sue Becker and Yann LeCun. Improving the convergence of
back-propagation learning with second order methods. Technical Report CRG-TR-88-5, Department of Computer Science,
University of Toronto, 1988.

[DHS11]

John C. Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization.
Journal of Machine Learning Research, 12:2121–2159, 2011.

[GB10]

Xavier Glorot and Yoshua Bengio. Understanding the difficulty
of training deep feedforward neural networks. In International
56

Conference on Artificial Intelligence and Statistics, pages 249–
256, 2010.
[GHL87]

S. Gallot, D. Hulin, and J. Lafontaine. Riemannian geometry.
Universitext. Springer-Verlag, Berlin, 1987.

[Kur94]

Takio Kurita. Iterative weighted least squares algorithms for
neural networks classifiers. New Generation Comput., 12(4):375–
394, 1994.

[LBOM96] Yann LeCun, Léon Bottou, Genevieve B. Orr, and Klaus-Robert
Müller. Efficient backprop. In Genevieve B. Orr and KlausRobert Müller, editors, Neural Networks: Tricks of the Trade,
volume 1524 of Lecture Notes in Computer Science, pages 9–50.
Springer, 1996.
[LNC+ 11] Quoc V. Le, Jiquan Ngiam, Adam Coates, Ahbik Lahiri, Bobby
Prochnow, and Andrew Y. Ng. On optimization methods for
deep learning. In Lise Getoor and Tobias Scheffer, editors, ICML,
pages 265–272. Omnipress, 2011.
[OAAH11] Yann Ollivier, Ludovic Arnold, Anne Auger, and Nikolaus Hansen. Information-Geometric Optimization algorithms:
A unifying picture via invariance principles.
Preprint,
arXiv:1106.3708v2, 2011.
[Oll13]

Yann Ollivier. Riemannian metrics for neural networks II: recurrent networks and learning symbolic data sequences. Preprint,
arXiv:1306.0514, 2013.

[PB13]

Razvan Pascanu and Yoshua Bengio. Revisiting natural gradient for deep networks. Preprint, http://arxiv.org/abs/1301.3584,
2013.

[RMB07]

Nicolas Le Roux, Pierre-Antoine Manzagol, and Yoshua Bengio.
Topmoumoute online natural gradient algorithm. In John C.
Platt, Daphne Koller, Yoram Singer, and Sam T. Roweis, editors,
NIPS. Curran Associates, Inc., 2007.

[RN03]

Stuart Russell and Peter Norvig. Artificial Intelligence: A Modern Approach. Prentice-Hall, Englewood Cliffs, NJ, 2nd edition,
2003.

[SZL13]

Tom Schaul, Sixin Zhang, and Yann LeCun. No more pesky
learning rates. In Proc. International Conference on Machine
learning (ICML’13), 2013.

57

