Metric-Free Natural Gradient for Joint-Training of Boltzmann Machines
Authors: Guillaume Desjardins, Razvan Pascanu, Aaron Courville, Yoshua Bengio
ICLR 2013

Abstract: This paper introduces the Metric-Free Natural Gradient (MFNG) algorithm for training Boltzmann Machines. Similar in spirit to the Hessian-Free method of Martens [8], our algorithm belongs to the family of truncated Newton methods and exploits an efficient matrix-vector product to avoid explicitely storing the natural gradient metric $L$. This metric is shown to be the expected second derivative of the log-partition function (under the model distribution), or equivalently, the variance of the vector of partial derivatives of the energy function. We evaluate our method on the task of joint-training a 3-layer Deep Boltzmann Machine and show that MFNG does indeed have faster per-epoch convergence compared to Stochastic Maximum Likelihood with centering, though wall-clock performance is currently not competitive.

Metric-Free Natural Gradient for Joint-Training of
Boltzmann Machines
Guillaume Desjardins, Razvan Pascanu, Aaron Courville and Yoshua Bengio
D´epartement d’informatique et de recherche op´erationnelle
Universit´e de Montr´eal
Abstract
This paper introduces the Metric-Free Natural Gradient (MFNG) algorithm for
training Boltzmann Machines. Similar in spirit to the Hessian-Free method of
Martens [8], our algorithm belongs to the family of truncated Newton methods and
exploits an efﬁcient matrix-vector product to avoid explicitly storing the natural
gradient metric L. This metric is shown to be the expected second derivative
of the log-partition function (under the model distribution), or equivalently, the
covariance of the vector of partial derivatives of the energy function. We evaluate
our method on the task of joint-training a 3-layer Deep Boltzmann Machine and
show that MFNG does indeed have faster per-epoch convergence compared to
Stochastic Maximum Likelihood with centering, though wall-clock performance
is currently not competitive.
1 Introduction
Boltzmann Machines (BM) have become a popular method in Deep Learning for performing fea-
ture extraction and probability modeling. The emergence of these models as practical learning
algorithms stems from the development of efﬁcient training algorithms, which estimate the negative
log-likelihood gradient by either contrastive [4] or stochastic [18, 19] approximations. However, the
success of these models has for the most part been limited to the Restricted Boltzmann Machine
(RBM) [6], whose architecture allows for efﬁcient exact inference. Unfortunately, this comes at the
cost of the model’s representational capacity, which is limited to a single layer of latent variables.
The Deep Boltzmann Machine (DBM) [15] addresses this by deﬁning a joint energy function over
multiple disjoint layers of latent variables, where interactions within a layer are prohibited. While
this affords the model a rich inference scheme incorporating top-down feedback, it also makes train-
ing much more difﬁcult, requiring until recently an initial greedy layer-wise pretraining scheme.
Since, Montavon and Muller [9] have shown that this difﬁculty stems from an ill-conditioning of
the Hessian matrix, which can be addressed by a simple reparameterization of the DBM energy
function, a trick called centering (an analogue to centering and skip-connections found in the de-
terministic neural network literature [17, 14]). As the barrier to joint-training 1 is overcoming a
challenging optimization problem, it is apparent that second-order gradient methods might prove to
be more effective than simple stochastic gradient methods. This should prove especially important
as we consider models with increasingly complex posteriors or higher-order interactions between
latent variables.
To this end, we explore the use of the Natural Gradient [2], which seems ideally suited to the stochas-
tic nature of Boltzmann Machines. Our paper is structured as follows. Section 2 provides a detailed
derivation of the natural gradient, including its speciﬁc form for BMs. While most of these equations
1Joint-training refers to the act of jointly optimizingθ(the concatenation of all model parameters, across all
layers of the DBM) through maximum likelihood. This is in contrast to [15], where joint-training is preceded
by a greedy layer-wise pretraining strategy.
1
arXiv:1301.3545v2  [cs.LG]  16 Mar 2013
have previously appeared in [3], our derivation aims to be more accessible as it attempts to derive the
natural gradient from basic principles, while minimizing references to Information Geometry. Sec-
tion 3 represents the true contribution of the paper: a practical natural gradient algorithm for BMs
which exploits the persistent Markov chains of Stochastic Maximum Likelihood (SML) [18], with
a Hessian-Free (HF) like algorithm [8]. The method, named Metric-Free Natural Gradient (MFNG)
(in recognition of the similarities of our method to HF), avoids explicitly storing the natural gradient
metric Land uses a linear solver to perform the required matrix-vector product L−1Eq[∇log pθ].
Preliminary experimental results on DBMs are presented in Section 4, with the discussion appearing
in Section 5.
2 The Natural Gradient
2.1 Motivation and Derivation
The main insight behind the natural gradient is that the space of all probability distributions P=
{pθ(x); θ ∈Θ,x ∈χ}forms a Riemannian manifold. Learning, which typically proceeds by
iteratively adapting the parameters θ to ﬁt an empirical distribution q, thus traces out a path along
this manifold. An immediate consequence is that following the direction of steepest descent in the
original Euclidean parameter space does not correspond to the direction of steepest descent alongP.
To do so, one needs to account for the metric describing the local geometry of the manifold, which
is given by the Fisher Information matrix [1], shown in Equation 4. While this metric is typically
derived from Information Geometry, a derivation more accessible to a machine learning audience
can be obtained as follows.
The natural gradient aims to ﬁnd the search direction ∆θwhich minimizes a given objective func-
tion, such that the Kullback–Leibler divergence KL(pθ ∥ pθ+∆θ) remains constant throughout
optimization. This constraint ensures that we make constant progress regardless of the curvature
of the manifold Pand enforces an invariance to the parameterization of the model . The natural
gradient for maximum likelihood can thus be formalized as:
∇N := ∆θ∗ ← argmin∆θ Eq[−log pθ+∆θ(x)] (1)
s.t. KL(pθ ∥pθ+∆θ) = const.
In order to derive a useful parameter update rule, we will consider the KL divergence under the
assumption ∆θ →0. We also assume we have a discrete and bounded domain χover which we
deﬁne the probability mass function2 pθ. Taking the Taylor series expansion oflog pθ+∆θ around θ,
and denoting ∇f as the column vector of partial derivatives with ∂f
∂θi
as the i-th entry, and ∇2f the
Hessian matrix with ∂2f
∂θi∂θj
in position (i,j), we have:
KL(pθ ∥pθ+∆θ) ≈
∑
χ
pθlog pθ −
∑
χ
pθ
[
log pθ + (∇log pθ)T ∆θ+ 1
2∆θT (
∇2 log pθ
)
∆θ
]
= 1
2∆θTEpθ
[
−∇2 log pθ
]
∆θ (2)
with the transition stemming from the fact that ∑
χpθ
∂log pθ
∂θi
= ∂
∂θi
∑
x∈χpθ(x) = 0 . Replacing
the objective function of Equation 1 by its ﬁrst-order Taylor expansion and rewriting the constraint
as a Lagrangian, we arrive at the following formulation for L(θ,∆θ), the loss function which the
natural gradient seeks to minimize.
L(θ,∆θ) = Eq[−log pθ] + Eq[−∇log pθ]T ∆θ+ λ
2 ∆θTEpθ
[
−∇2 log pθ
]
∆θ.
Setting ∂L
∂∆θ to zero yields the natural gradient direction ∇N:
∇N = L−1Eq[∇log pθ] with L = Epθ
[
−∇2 log pθ
]
(3)
or equivalently L = Epθ
[
∇log pθ∇T log pθ
]
(4)
2When clear from context, we will drop the argument of pθ to save space.
2
While its form is reminiscent of the Newton direction, the natural gradient multiplies the estimated
gradient by the inverse of the expected Hessian of log pθ (Equation 3) or equivalently by the Fisher
Information matrix (FIM, Equation 4). The equivalence between both expressions can be shown
trivially, with the details appearing in the Appendix. We stress that both of these expectations
are computed with respect to the model distribution, and thus computing the metric L does not
involve the empirical distribution in any way. The FIM for Boltzmann Machines is thus not equal
to the uncentered covariance of the maximum likelihood gradients. In the following, we pursue our
derivation from the form given in Equation 4.
2.2 Natural Gradient for Boltzmann Machines
Derivation. Boltzmann machines deﬁne a joint distribution over a vector of binary random vari-
ables x ∈ {0,1}N by way of an energy function E(x) = −∑
k<lWklxkxl −∑
kbkxk, with
weight matrix W ∈RN×N and bias vector b ∈RN. Energy and probability are related by the
Boltzmann distribution, such that p(x) = 1
Z exp (−E(x)), with Z the partition function deﬁned by
Z = ∑
xexp (−E(x)).
Starting from the expression of Lfound in Equation 3, we can derive the natural gradient metric for
Boltzmann Machines.
L(BM) = Epθ
[
∇2E(x) + ∇2 log Z
]
= Epθ
[
∇2 log Z
]
The natural gradient metric for ﬁrst-order BMs takes on a surprisingly simple form: it is the expected
Hessian of the log-partition function. With a few lines of algebra (whose details are presented in the
Appendix), we can rewrite it as follows:
L(BM) = Epθ
[
(∇E(x) −Epθ [∇E(x)])T (∇E(x) −Epθ [∇E(x)])
]
(5)
L(BM) is thus given by the covariance of∇E, measured under the model distributionpθ. Concretely,
if we denote Wkl and Wmn as the i and j-th parameters of the model respectively, the entryLij will
take on the value −E[xkxlxmxn] + E[xkxl] E[xmxn].
Discussion. When computing the Taylor expansion of the KL divergence in Equation 2, we
glossed over an important detail. Namely, how to handle latent variables in pθ(x), a topic ﬁrst dis-
cussed in [3]. If x= [v,h], we could just as easily have derived the natural gradient by considering
the constraint KL(∑
hpθ(v,h) ∥∑
hpθ+∆θ(v,h)) = const. Alternatively, since the distinction
between visible and hidden units is entirely artiﬁcial (since the KL divergence does not involve the
empirical distribution), we may simply wish to consider the distribution obtained by analytically
integrating out a maximal number of random variables. In a DBM, this would entail marginalizing
over all odd or even layers, a strategy employed with great success in the context of AIS [15]. In this
work however, we only consider the metric obtained by considering theKLdivergence between the
full joint distributions pθ and pθ+∆θ.
3 Metric-Free Natural Gradient Implementation
We can compute the natural gradient∇N by ﬁrst replacing the expectations of Equation 5 by a ﬁnite
sample approximation. We can do this efﬁciently by reusing the model samples generated by the
persistent Markov chains of SML. Given the size of the matrix being estimated however, we expect
this method to require a larger number of chains than is typically used. The rest of the method is
similar to the Hessian-Free (HF) algorithm of Martens [8]: we exploit an efﬁcient matrix-vector
implementation combined with a linear-solver, such as Conjugate Gradient or MinRes[12], to solve
the system Ly = Eq[∇log pθ] for y ∈RN. Additionally, we replace the expectation on the rhs.
of this previous equation by an average computed over a mini-batch of training examples (sampled
from the empirical distribution q), as is typically done in the stochastic learning setting.
For Boltzmann Machines, the matrix-vector product Lycan be computed in a straightforward man-
ner, without recourse to Pearlmutter’s R-operator [13]. Starting from a sampling approximation to
Equation 5, we simply push the dot product inside of the expectation as follows:
3
Algorithm 1MFNG iteration(θ,X+,Z−
old)
θ: parameters of the model. N :=|θ|.
X+: mini-batch of training examples, with X+ = {xm; m∈[1,M]}.
Z−
old: previous state of Markov chains, with Z= {zm := (vm,h(1)
m ,h(2)
m ); m∈[1,M]}
•Generate positive phase samples Z+ = {z+
m := (xm,h(1)+
m ,h(2)+
m ); m∈[1,M]}
•Initializing M Markov chains from state Z−
old, generate negative phase samples Z−
new.
•Compute the vectors s+
m = ∂E(z+
m)
∂θ and s−
m = ∂E(z−
m)
∂θ , ∀m.
•Compute negative log-likelihood gradient as g= 1
M
∑
m(s+
m −s−
m).
•Denote S ∈RM×N as the matrix with rows s−
m and ¯S = 1
M
∑
ms−
m.
# Solve the system “Ly= g” for y, given L= (S−¯S)T(S−¯S) and an initial zero vector.
# computeLy is a function which performs equation 6, without instantiating L.
•∇Nθ←CGSolve(computeLy,S,g, zeros(N))
L(BM)y ≈
(
S−¯S
)T [(
S−¯S
)
y
]
(6)
with S ∈RM×N, the matrix with entries smj = ∂E(xm)
∂θj
and ¯S ∈RN, the vector with entries sj = 1
M
∑
m
smj
and xm ∼pθ(x),m ∈[1,M].
By ﬁrst computing the matrix-vector product (S −¯S)y, we can easily avoid computing the full
N ×N matrix L. Indeed, the result of this operation is a vector of length M, which is then left-
multiplied by a matrix of dimension N ×M, yielding the matrix-vector product Ly ∈RN. A
single iteration of the MFNG is presented in Algorithm 1. A full open-source implementation is
also available online.3.
4 Experiments
We performed a proof-of-concept experiment to determine whether our Metric-Free Natural Gra-
dient (MFNG) algorithm is suitable for joint-training of complex Boltzmann Machines. To this
end, we compared our method to Stochastic Maximum Likelihood and a diagonal approximation of
MFNG on a 3-layer Deep Boltzmann Machine trained on MNIST [7]. All algorithms were run in
conjunction with the centering strategy of Montavon and Muller [9], which proved crucial to suc-
cessfully joint-train all layers of the DBM (even when using MFNG) 4. We chose a small 3-layer
DBM with 784-400-100 units at the ﬁrst, second and third layers respectively, to be comparable to
[10]. Hyper-parameters were varied as follows. For inference, we ran 5 iterations of either mean-
ﬁeld as implemented in [15] or Gibbs sampling. The learning rate was kept ﬁxed during training and
chosen from the set {5 ·10−3,10−3,10−4}. For MinRes, we set the damping coefﬁcient to 0.1 and
used a ﬁxed tolerance of 10−5 (used to determine convergence). Finally, we tested all algorithms
on minibatch sizes of either 25, 128 or 256 elements 5. Finally, since we are comparing optimiza-
tion algorithms, hyper-parameters were chosen based on the training set likelihood (though we still
report the associated test errors). All experiments used the MinRes linear solver, both for its speed
and its ability to return pseudo-inverses when faced with ill-conditioning.
3https://github.com/gdesjardins/MFNG
4The centering coefﬁcients were initialized as in [9], but were otherwise held ﬁxed during training.
5We expect larger minibatch sizes to be preferable, however simulating this number of Markov chains in
parallel (on top of all other memory requirements) was sufﬁcient to hit the memory bottlenecks of GPUs.
4
Figure 1: Estimated model likelihood as a function of (left) epochs and (right) CPU-time for MFNG,
its diagonal approximation (MFNG-diag) and SML. All methods were run in conjunction with the
DBM centering trick [9], with centering coefﬁcients held ﬁxed during training. Our grid search
yielded the following hyper-parameters: batch size of 256/128 for MFNG(-diag)/SGD; 5 steps of
mean-ﬁeld / sampling-based inference for MFNG(-diag)/SGD and a learning rate of 5 ·10−3.
Figure 1 (left) shows the likelihood as estimated by Annealed Importance Sampling [15, 11] as a
function of the number of epochs 6. Under this metric, MFNG achieves the fastest convergence,
obtaining a training/test set likelihood of −71.26/−72.84 nats after 94 epochs. In comparison,
MFNG-diag obtains −73.22/−74.05 nats and SML −80.12/−79.71 nats in 100 epochs. The picture
changes however when plotting likelihood as a function of CPU-time, as shown in Figure 1 (right).
Given a wall-time of 8000s for MFNG and SML, and 5000s for MFNG-diag 7, SML is able to
perform upwards of 1550 epochs, resulting in an impressive likelihood score of −64.94 / −67.73.
Note that these results were obtained on the binary-version of MNIST (thresholded at 0.5) in order
to compare to [10]. These results are therefore not directly comparable to [15], which binarizes
the dataset through sampling (by treating each pixel activation as the probability pof a Bernouilli
distribution).
Figure 2 shows a breakdown of the algorithm runtime, for various components of the algorithm.
These statistics were collected in the early stages of training, but are generally representative of the
bigger picture. While the linear solver clearly dominates the runtime, there are a few interesting
observations to make. For small models and batch sizes greater than 256, a single evaluation of Ly
appears to be of the same order of magnitude as a gradient evaluation. In all cases, this cost is smaller
than that of sampling, which represents a non-negligible part of the total computational budget.
This suggests that MFNG could become especially attractive for models which are expensive to
sample from. Overall however, restricting the number of CG/MinRes iterations appears key to
computational performance, which can be achieved by increasing the damping factor α. How this
affects convergence in terms of likelihood is left for future work.
5 Discussion and Future Work
While the wall-clock performance of MFNG is not currently competitive with SML, we believe there
are still many avenues to explore to improve computational efﬁciency. Firstly, we performed almost
no optimization of the various MinRes hyper-parameters. In particular, we ran the algorithm to
convergence with a ﬁxed tolerance of 10−5. While this typically resulted in relatively few iterations
(around 15), this level of precision might not be required (especially given the stochastic nature of
6While we do not report error margins for AIS likelihood estimates, the numbers proved robust to changes in
the number of particles and temperatures being simulated. To obtain such robust estimates, we implemented all
the tricks described in Salakhutdinov and Hinton [15] and [16]:pAa zero-weight base-rate model whose biases
are set by maximum likelihood; interpolating distributions pi ∝ p(1−βi)
A p(βi)
B , with pB the target distribution;
and ﬁnally analytical integration of all odd-layers.
7This discrepancy will be resolved in the next revision.
5
Figure 2: Breakdown of algorithm runtime, when we vary (left) the batch size (with ﬁxed model
architecture 784 −400 −100) and (right) the model size (with ﬁxed batch size of 256). Runtime
is additive in the order given by the labels (top to bottom). Dotted lines denote intermediate steps,
while continuous lines denote full steps. Data was collected on a Nvidia GTX 480 card.
the algorithm). Additionally, it could be worth exploiting the same strategy as HF where the linear
solver is initialized by the solution found in the previous iteration. This may prove much more
efﬁcient than the current approach of initializing the solver with a zero vector. Pre-conditioning
is also a well-known method for accelerating the convergence speed of linear solvers [5]. Our
implementation used a simple diagonal regularization of L. The Jacobi preconditioner could be
implemented easily however by computing the diagonal of Lin a ﬁrst-pass.
Finally, while our single experiment offers little evidence in support of either conclusion, it may very
well be possible that MFNG is simply not computationally efﬁcient for DBMs, compared to SML
with centering. In this case, it would be worth applying the method to either (i) models with known
ill-conditioning, such as factored 3-rd order Boltzmann Machines or (ii) models and distributions
exhibiting complex posterior distributions. In such scenarios, we may wish to maximize the use
of the positive phase statistics (which were obtained at a high computational cost) by performing
larger jumps in parameter space. It remains to be seen how this would interact with SML, where the
burn-in period of the persistent chains is directly tied to the magnitude of ∆θ.
Appendix
We include the following derivations for completeness.
5.1 Expected Hessian oflog Zand Fisher Information.
Epθ
[
−∂2 log p(x)
∂θi∂θj
]
= Epθ
[ 1
p(x)2
∂p(x)
∂θj
∂p(x)
∂θi
− 1
p(x)
∂2p(x)
∂θi∂θj
]
= Epθ
[( 1
p(x)
∂p(x)
∂θi
)( 1
p(x)
∂p(x)
∂θj
)
− 1
p(x)
∂2p(x)
∂θi∂θj
]
= Epθ
[∂log p(x)
∂θi
∂log p(x)
∂θj
]
−
∑
x
p(x) 1
p(x)
∂2p(x)
∂θi∂θj
= Epθ
[∂log p(x)
∂θi
∂log p(x)
∂θj
]
−
∑
x
∂2p(x)
∂θi∂θj
= Epθ
[∂log p(x)
∂θi
∂log p(x)
∂θj
]
−∂2 ∑
xp(x)
∂θi∂θj
= Epθ
[∂log p(x)
∂θi
∂log p(x)
∂θj
]
6
5.2 Derivation of Equation 5
log p(x) = −E(x) −log Z
∂log p(x)
∂θi
= −∂E(x)
∂θi
−1
Z
∑
x
∂
∂θi
[exp (−E(x))]
= −∂E(x)
∂θi
+ Epθ
[∂E(x)
∂θi
]
Epθ
[
−∂2 log p(x)
∂θi∂θj
]
= Epθ
[(∂E(x)
∂θi
−Epθ
[∂E(x)
∂θi
])(∂E(x)
∂θj
−Epθ
[∂E(x)
∂θj
])]
= Epθ
[∂E(x)
∂θi
∂E(x)
∂θj
]
−Epθ
[∂E(x)
∂θi
]
Epθ
[∂E(x)
∂θj
]
References
[1] Amari, S. (1985). Differential geometrical methods in statistics. Lecture notes in statistics, 28.
[2] Amari, S. (1998). Natural gradient works efﬁciently in learning. Neural Computation, 10(2), 251–276.
[3] Amari, S., Kurata, K., and Nagaoka, H. (1992). Information geometry of Boltzmann machines. IEEE
Trans. on Neural Networks, 3, 260–271.
[4] Carreira-Perpi ˜nan, M. A. and Hinton, G. E. (2005). On contrastive divergence learning. In AISTATS’2005,
pages 33–40.
[5] Chapelle, O. and Erhan, D. (2011). Improved preconditioner for hessian free optimization. NIPS Workshop
on Deep Learning and Unsupervised Feature Learning.
[6] Freund, Y . and Haussler, D. (1992). A fast and exact learning rule for a restricted class of Boltzmann
machines. pages 912–919, Denver, CO. Morgan Kaufmann, San Mateo.
[7] LeCun, Y ., Bottou, L., Bengio, Y ., and Haffner, P. (1998). Gradient based learning applied to document
recognition. Proc. IEEE.
[8] Martens, J. (2010). Deep learning via Hessian-free optimization. pages 735–742.
[9] Montavon, G. and Muller, K.-R. (2012). Deep Boltzmann machines and the centering trick. In G. Mon-
tavon, G. Orr, and K.-R. Muller, editors, Neural Networks: Tricks of the Trade , volume 7700 of Lecture
Notes in Computer Science, pages 621–637.
[10] Montavon, G. and M ¨uller, K.-R. (2012). Learning feature hierarchies with centered deep Boltzmann
machines. CoRR, abs/1203.4416.
[11] Neal, R. M. (2001). Annealed importance sampling. Statistics and Computing, 11(2), 125–139.
[12] Paige, C. C. and Saunders, M. A. (1975). Solution of Sparse Indeﬁnite Systems of Linear Equations.
SIAM Journal on Numerical Analysis, 12(4), 617–629.
[13] Pearlmutter, B. (1994). Fast exact multiplication by the Hessian. Neural Computation, 6(1), 147–160.
[14] Raiko, T., Valpola, H., and LeCun, Y . (2012). Deep learning made easier by linear transformations in
perceptrons. In AISTATS’2012.
[15] Salakhutdinov, R. and Hinton, G. (2009). Deep Boltzmann machines. In Proceedings of the Twelfth
International Conference on Artiﬁcial Intelligence and Statistics (AISTATS 2009), volume 8.
[16] Salakhutdinov, R. and Murray, I. (2008). On the quantitative analysis of deep belief networks. volume 25,
pages 872–879.
[17] Schraudolph, N. N. (1998). Centering neural network gradient factors. In G. B. Orr and K.-R. Muller,
editors, Neural Networks: Tricks of he Trade, pages 548–548. Springer.
[18] Tieleman, T. (2008). Training restricted Boltzmann machines using approximations to the likelihood
gradient. pages 1064–1071.
[19] Younes, L. (1999). On the convergence of Markovian stochastic algorithms with rapidly decreasing
ergodicity rates. Stochastics and Stochastic Reports, 65(3), 177–228.
7


Peer review:
[Review]:
This paper presents a natural gradient algorithm for deep Boltzmann machines. The authors must be commended for their extremely clear and succinct description of the natural gradient method in Section 2. This presentation is particularly useful because, indeed, many of the papers on information geometry are hard to follow. The derivations are also correct and sound. The derivations in the appendix are classical statistics results, but their addition is likely to improve readability of the paper.

The experiments show that the natural gradient approach does better than stochastic maximum likelihood when plotting estimated likelihood against epochs. However, per unit computation, the stochastic maximum likelihood method still does better. 

I was not able to understand remark 4 about mini-batches. Why are more parallel chains needed? Why not simply use a single chain but have longer memory. I strongly think this part of the paper could be improved if the authors write down the pseudo-code for their algorithm. Another suggestion is to use automatic algorithm configuration to find the optimal hyper-parameters for each method, given that they are so close.

The trade-offs of second order versus first order optimization methods are well known in the deterministic case. There is is also some theoretical guidance for the stochastic case. I encourage the authors to look at the following papers for this:

A Stochastic Gradient Method with an Exponential Convergence Rate for Finite Training Sets. N. Le Roux, M. Schmidt, F. Bach. NIPS, 2012. 

Hybrid Deterministic-Stochastic Methods for Data Fitting.
M. Friedlander, M. Schmidt. SISC, 2012. 

'On the Use of Stochastic Hessian Information in Optimization Methods for Machine Learning' R. Byrd, G. Chin and W. Neveitt, J. Nocedal.
SIAM J. on Optimization, vol 21, issue 3, pages 977-995 (2011).

'Sample Size Selection in Optimization Methods for Machine Learning'
R. Byrd, G. Chin, J. Nocedal and Y. Wu. to appear in Mathematical Programming B (2012).

In practical terms, given that the methods are so close, how does the choice of implementation (GPUs, multi-cores, single machine) affect the comparison? Also, how data dependent are the results. I would be nice to gain a deeper understanding of the conditions under which the natural gradient might or might not work better than stochastic maximum likelihood when training Boltzmann machines.

Finally, I would like to point out a few typos to assist in improving the paper:

Page 1: litterature should be literature
Section 2.2 cte should be const for consistency.
Section 3: Avoid using x instead of grad_N in the linear equation for Lx=E(.)  This causes overloading. For consistency with the previous section, please use grad_N instead.
Section 4: Add a space between MNIST and [7].
Appendix 5.1: State that the expectation is with respect to p_{	heta}(x).
Appendix 5.2: The expectation with respect to q_	heta should be with respect to p_{	heta}(x) to ensure consistency of notation, and correctness in this case.
References: References [8] and [9] appear to be duplicates of the same paper by J. Martens.