What Regularized Auto-Encoders Learn from the Data Generating Distribution
Authors: Guillaume Alain, Yoshua Bengio
ICLR 2013

Abstract: What do auto-encoders learn about the underlying data generating distribution? Recent work suggests that some auto-encoder variants do a good job of capturing the local manifold structure of data. This paper clarifies some of these previous intuitive observations by showing that minimizing a particular form of regularized reconstruction error yields a reconstruction function that locally characterizes the shape of the data generating density. We show that the auto-encoder captures the score (derivative of the log-density with respect to the input), along with the second derivative of the density and the local mean associated with the unknown data-generating density. This is the second result linking denoising auto-encoders and score matching, but in way that is different from previous work, and can be applied to the case when the auto-encoder reconstruction function does not necessarily correspond to the derivative of an energy function. The theorems provided here are completely generic and do not depend on the parametrization of the auto-encoder: they show what the auto-encoder would tend to if given enough capacity and examples. These results are for a contractive training criterion we show to be similar to the denoising auto-encoder training criterion with small corruption noise, but with contraction applied on the whole reconstruction function rather than just encoder. Similarly to score matching, one can consider the proposed training criterion as a convenient alternative to maximum likelihood, i.e., one not involving a partition function.

What Regularized Auto-Encoders Learn from the Data
Generating Distribution
Guillaume Alain and Yoshua Bengio
guillaume.alain@umontreal.ca, yoshua.bengio@umontreal.ca
Department of Computer Science and Operations Research
University of Montreal
Montreal, H3C 3J7, Quebec, Canada
Abstract
What do auto-encoders learn about the underlying data generating distribution? Recent work
suggests that some auto-encoder variants do a good job of capturing the local manifold structure
of data. This paper clariﬁes some of these previous observations by showing that minimizing
a particular form of regularized reconstruction error yields a reconstruction function that lo-
cally characterizes the shape of the data generating density. We show that the auto-encoder
captures the score (derivative of the log-density with respect to the input). It contradicts pre-
vious interpretations of reconstruction error as an energy function. Unlike previous results, the
theorems provided here are completely generic and do not depend on the parametrization of
the auto-encoder: they show what the auto-encoder would tend to if given enough capacity
and examples. These results are for a contractive training criterion we show to be similar to
the denoising auto-encoder training criterion with small corruption noise, but with contrac-
tion applied on the whole reconstruction function rather than just encoder. Similarly to score
matching, one can consider the proposed training criterion as a convenient alternative to max-
imum likelihood because it does not involve a partition function. Finally, we show how an
approximate Metropolis-Hastings MCMC can be setup to recover samples from the estimated
distribution, and this is conﬁrmed in sampling experiments.
1. Introduction
Machine learning is about capturing aspects of the unknown distribution from which the ob-
served data are sampled (the data-generating distribution). For many learning algorithms and
in particular in manifold learning , the focus is on identifying the regions (sets of points) in
the space of examples where this distribution concentrates, i.e., which conﬁgurations of the
observed variables are plausible.
Unsupervised representation-learningalgorithms attempt to characterize the data-generating
distribution through the discovery of a set of features or latent variables whose variations cap-
ture most of the structure of the data-generating distribution. In recent years, a number of
unsupervised feature learning algorithms have been proposed that are based on minimizing
some form of reconstruction error, such as auto-encoder and sparse coding variants (Olshausen
and Field, 1997; Bengio et al., 2007; Ranzato et al., 2007; Jain and Seung, 2008; Ranzato
et al., 2008; Vincent et al., 2008; Kavukcuoglu et al., 2009; Rifai et al., 2011b,a; Gregor et al.,
2011). An auto-encoder reconstructs the input through two stages, an encoder function f
1
arXiv:1211.4246v5  [cs.LG]  19 Aug 2014
(which outputs a learned representation h = f(x) of an example x) and a decoder function
g, such that g(f(x)) ≈x for most x sampled from the data-generating distribution. These
feature learning algorithms can be stacked to form deeper and more abstract representations.
Deep learning algorithms learn multiple levels of representation, where the number of levels is
data-dependent. There are theoretical arguments and much empirical evidence to suggest that
when they are well-trained, deep learning algorithms (Hinton et al., 2006; Bengio, 2009; Lee
et al., 2009; Salakhutdinov and Hinton, 2009; Bengio and Delalleau, 2011; Bengio et al., 2013b)
can perform better than their shallow counterparts, both in terms of learning features for the
purpose of classiﬁcation tasks and for generating higher-quality samples.
Here we restrict ourselves to the case of continuous inputs x∈Rd with the data-generating
distribution being associated with an unknown target density function, denoted p. Manifold
learning algorithms assume that pis concentrated in regions of lower dimension (Cayton, 2005;
Narayanan and Mitter, 2010), i.e., the training examples are by deﬁnition located very close to
these high-density manifolds. In that context, the core objective of manifold learning algorithms
is to identify where the density concentrates.
Some important questions remain concerning many of feature learning algorithms based on
reconstruction error. Most importantly, what is their training criterion learning about the input
density? Do these algorithms implicitly learn about the whole density or only some aspect? If
they capture the essence of the target density, then can we formalize that link and in particular
exploit it to sample from the model ? The answers may help to establish that these algorithms
actually learn implicit density models, which only deﬁne a density indirectly, e.g., through the
estimation of statistics or through a generative procedure. These are the questions to which
this paper contributes.
The paper is divided in two main sections, along with detailed appendices with proofs of the
theorems. Section 2 makes a direct link between denoising auto-encoders (Vincent et al., 2008)
and contractive auto-encoders (Rifai et al., 2011b), justifying the interest in the contractive
training criterion studied in the rest of the paper. Section 3 is the main contribution and
regards the following question: when minimizing that criterion, what does an auto-encoder
learn about the data generating density ? The main answer is that it estimates the score (ﬁrst
derivative of the log-density), i.e., the direction in which density is increasing the most, which
also corresponds to the local mean, which is the expected value in a small ball around the
current location. It also estimates the Hessian (second derivative of the log-density).
Finally, Section 4 shows how having access to an estimator of the score can be exploited to
estimate energy diﬀerences, and thus perform approximate MCMC sampling. This is achieved
using a Metropolis-Hastings MCMC in which the energy diﬀerences between the proposal and
the current state are approximated using the denoising auto-encoder. Experiments on artiﬁcial
datasets show that a denoising auto-encoder can recover a good estimator of the data-generating
distribution, when we compare the samples generated by the model with the training samples,
projected into various 2-D views for visualization.
2. Contractive and Denoising Auto-Encoders
Regularized auto-encoders (see Bengio et al. (2012b) for a review and a longer exposition)
capture the structure of the training distribution thanks to the productive opposition between
reconstruction error and a regularizer. An auto-encoder maps inputs xto an internal represen-
tation (or code) f(x) through the encoder function f, and then maps back f(x) to the input
2
Figure 1: Regularization forces the auto-encoder to become less sensitive to the input, but
minimizing reconstruction error forces it to remain sensitive to variations along the
manifold of high density. Hence the representation and reconstruction end up cap-
turing well variations on the manifold while mostly ignoring variations orthogonal to
it.
space through a decoding function g. The composition of f and g is called the reconstruction
function r, with r(x) = g(f(x)), and a reconstruction loss function ℓ penalizes the error made,
with r(x) viewed as a prediction of x. When the auto-encoder is regularized, e.g., via a sparsity
regularizer, a contractive regularizer (detailed below), or a denoising form of regularization
(that we ﬁnd below to be very similar to a contractive regularizer), the regularizer basically
attempts to make r (or f) as simple as possible, i.e., as constant as possible, as unresponsive
to xas possible. It means that f has to throw away some information present in x, or at least
represent it with less precision. On the other hand, to make reconstruction error small on the
training set, examples that are neighbors on a high-density manifold must be represented with
suﬃciently diﬀerent values of f(x) or r(x). Otherwise, it would not be possible to distinguish
and hence correctly reconstruct these examples. It means that the derivatives of f(x) or r(x) in
the x-directions along the manifold must remain large, while the derivatives (of f or r) in the
x-directions orthogonal to the manifold can be made very small. This is illustrated in Figure 1.
In the case of principal components analysis, one constrains the derivative to be exactly 0 in the
directions orthogonal to the chosen projection directions, and around 1 in the chosen projection
directions. In regularized auto-encoders, f is non-linear, meaning that it is allowed to choose
diﬀerent principal directions (those that are well represented, i.e., ideally the manifold tangent
directions) at diﬀerent x’s, and this allows a regularized auto-encoder with non-linear encoder
to capture non-linear manifolds. Figure 2 illustrates the extreme case when the regularization
is very strong (r(·) wants to be nearly constant where density is high) in the special case where
the distribution is highly concentrated at three points (three training examples). It shows the
compromise between obtaining the identity function at the training examples and having a
ﬂat r near the training examples, yielding a vector ﬁeld r(x) −x that points towards the high
density points.
Here we show that the denoising auto-encoder (Vincent et al., 2008) with very small
Gaussian corruption and squared error loss is actually a particular kind of contractive auto-
encoder (Rifai et al., 2011b), contracting the whole auto-encoder reconstruction function rather
than just the encoder, whose contraction penalty coeﬃcient is the magnitude of the perturba-
tion. This was ﬁrst suggested in (Rifai et al., 2011c).
The contractive auto-encoder, or CAE (Rifai et al., 2011b), is a particular form of regularized
auto-encoder which is trained to minimize the following regularized reconstruction error:
LCAE = E
[
ℓ(x,r(x)) + λ

∂f(x)
∂x

2
F
]
(1)
where r(x) = g(f(x)) and ||A||2
F is the sum of the squares of the elements of A. Both the
squared loss ℓ(x,r) = ||x−r||2 and the cross-entropy loss ℓ(x,r) = −xlog r−(1 −x) log(1−r)
3
x"
r(x)"
x1" x2" x3"
Figure 2: The reconstruction function r(x) (in turquoise) which would be learned by a high-capacity
auto-encoder on a 1-dimensional input, i.e., minimizing reconstruction error at the training
examples xi (with r(xi) in red) while trying to be as constant as possible otherwise. The
ﬁgure is used to exagerate and illustrate the eﬀect of the regularizer (corresponding to a large
σ2 in the loss function Llater described by (6)). The dotted line is the identity reconstruction
(which might be obtained without the regularizer). The blue arrows shows the vector ﬁeld
of r(x) −x pointing towards high density peaks as estimated by the model, and estimating
the score (log-density derivative), as shown in this paper.
have been used, but here we focus our analysis on the squared loss because of the easier
mathematical treatment it allows. Note that success in minimizing the CAE criterion strongly
depends on the parametrization of f and g and in particular on the tied weights constraint
used, with f(x) = sigmoid( Wx + b) and g(h) = sigmoid( WTh+ c). The above regularizing
term forces f (as well as g, because of the tied weights) to be contractive, i.e., to have singular
values less than 1 1. Larger values of λ yield more contraction (smaller singular values) where
it hurts reconstruction error the least, i.e., in the local directions where there are only little or
no variations in the data. These typically are the directions orthogonal to the manifold of high
density concentration, as illustrated in Fig. 2.
The denoising auto-encoder, or DAE (Vincent et al., 2008), is trained to minimize the
following denoising criterion:
LDAE = E[ℓ(x,r(N(x)))] (2)
where N(x) is a stochastic corruption of xand the expectation is over the training distribution
and the corruption noise source. Here we consider mostly the squared loss and Gaussian noise
corruption, again because it is easier to handle them mathematically. In many cases, the exact
same proofs can be applied to any kind of additive noise, but Gaussian noise serves as a good
frame of reference.
1. Note that an auto-encoder without any regularization would tend to ﬁnd many leading singular values near
1 in order to minimize reconstruction error, i.e., preserve input norm in all the directions of variation present
in the data.
4
Theorem 1 Let p be the probability density function of the data. If we train a DAE using the
expected quadratic loss and corruption noise N(x) = x+ ϵ with
ϵ∼N
(
0,σ2I
)
,
then the optimal reconstruction function r∗(x) will be given by
r∗(x) = Eϵ[p(x−ϵ)(x−ϵ)]
Eϵ[p(x−ϵ)] (3)
for values of x where p(x) ̸= 0.
Moreover, if we consider how the optimal reconstruction function r∗
σ(x) behaves asymptoti-
cally as σ→0, we get that
r∗
σ(x) = x+ σ2 ∂log p(x)
∂x + o(σ2) as σ→0. (4)
The proof of this result is found in the Appendix. We make use of the small o notation
throughout this paper and assume that the reader is familiar with asymptotic notation. In the
context of Theorem 1, it has to be understood that all the other quantities except for σ are
ﬁxed when we study the eﬀect of σ→0. Also, note that the σ in the index of r∗
σ is to indicate
that r∗
σ was chosen based on the value of σ. That σ should not be mistaken for a parameter to
be learned.
Equation (3) reveals that the optimal DAE reconstruction function at every point xis given
by a kind of convolution involving the density function p, or weighted average from the points
in the neighbourhood of x, depending on how we would like to view it. A higher noise level σ
means that a larger neighbourhood of x is taken into account. Note that the total quantity of
“mass” being included in the weighted average of the numerator of (3) is found again at the
denominator.
Gaussian noise is a simple case in the sense that it is additive and symmetrical, so it avoids
the complications that would occur when trying to integrate over the density of pre-images x′
such that N(x′) = x for a given x. The ratio of those quantities that we have in equation (3),
however, depends strongly on the decision that we made to minimize the expected square error.
When we look at the asymptotic behavior with equation (4), the ﬁrst thing to observe is
that the leading term in the expansion of r∗
σ(x) is x, and then the remainder goes to 0 as σ→0.
When there is no noise left at all, it should be clear that the best reconstruction target for any
value x would be that x itself.
We get something even more interesting if we look at the second term of equation (4)
because it gives us an estimator of the score from
∂log p(x)
∂x = (r∗
σ(x) −x) /σ2 + o(1) as σ→0. (5)
This result is at the core of our paper. It is what allows us to start from a trained DAE, and
then recover properties of the training density p(x) that can be used to sample from p(x).
Most of the asymptotic properties that we get by considering the limit as the Gaussian
noise level σ goes to 0 could be derived from a family of noise distribution that approaches a
point mass distribution in a relatively “nice” way.
5
An interesting connection with contractive auto-encoders can also be observed by using a
Taylor expansion of the denoising auto-encoder loss and assuming only that rσ(x) = x+ o(1)
as σ→0. In that case we get the following proposition.
Proposition 1 Let p be the probability density function of the data. Consider a DAE using
the expected quadratic loss and corruption noise N(x) = x+ ϵ, with ϵ ∼N
(
0,σ2I
)
. If we
assume that the non-parametric solutions rσ(x) satistﬁes
rσ(x) = x+ o(1) as σ→0,
then we can rewrite the loss as
LDAE (rσ) = E
[
∥rσ(x) −x∥2
2 + σ2

∂rσ(x)
∂x

2
F
]
+ o(σ2) as σ→0
where the expectation is taken with respect to X, whose distribution is given by p.
The proof is in Appendix and uses a simple Taylor expansion around x.
Proposition 1 shows that the DAE with small corruption of variance σ2 is similar to a
contractive auto-encoder with penalty coeﬃcient λ= σ2 but where the contraction is imposed
explicitly on the whole reconstruction function r(·) = g(f(·)) rather than on f(·) alone2.
This analysis motivates the deﬁnition of thereconstruction contractive auto-encoder(RCAE),
a variation of the CAE where loss function is instead the squared reconstruction loss plus con-
tractive penalty on the reconstruction:
LRCAE = E
[
∥r(x) −x∥2
2 + σ2

∂r(x)
∂x

2
F
]
. (6)
This is an analytic version of the denoising criterion with small noise σ2, and also corresponds
to a contractive auto-encoder with contraction on both f and g, i.e., on r.
Because of the similarity between DAE and RCAE when taking λ = σ2 and because the
semantics of σ2 is clearer (as a squared distance in input space), we will denoteσ2 for the penalty
term coeﬃcient in situations involving RCAE. For example, in the statement of Theorem 2,
this σ2 is just a positive constant; there is no notion of additive Gaussian noise, i.e., σ2 does
not explicitly refer to a variance, but using the notation σ2 makes it easier to intuitively see
the connection to the DAE setting.
The connection between DAE and RCAE established in Proposition 1 motivates the follow-
ing Theorem 2 as an alternative way to achieve a result similar to that of Theorem 1. In this
theorem we study the asymptotic behavior of the RCAE solution.
Theorem 2 Let pbe a probability density function that is continuously diﬀerentiable once and
with support Rd (i.e. ∀x∈Rd we have p(x) ̸= 0). Let Lσ be the loss function deﬁned by
Lσ(r) =
∫
Rd
p(x)
[
∥r(x) −x∥2
2 + σ2

∂r(x)
∂x

2
F
]
dx (7)
2. In the CAE there is a also a contractive eﬀect on g(·) as a side eﬀect of the parametrization with weights
tied between f(·) and g(·).
6
for r: Rd →Rd assumed to be diﬀerentiable twice, and 0 ≤σ2 ∈R used as factor to the penalty
term.
Let r∗
σ(x) denote the optimal function that minimizes Lσ. Then we have that
r∗
σ(x) = x+ σ2 ∂log p(x)
∂x + o(σ2) as σ→0.
Moreover, we also have the following expression for the derivative
∂r∗
σ(x)
∂x = I+ σ2 ∂2 log p(x)
∂x2 + o(σ2) as σ→0. (8)
Both these asymptotic expansions are to be understood in a context where we consider
{r∗
σ(x)}σ≥0 to be a family of optimal functions minimizing Lσ for their corresponding value
of σ. The asymptotic expansions are applicable point-wise in x, that is, with any ﬁxed x we
look at the behavior as σ→0.
The proof is given in the appendix and uses the Euler-Lagrange equations from the calculus
of variations.
3. Minimizing the Loss to Recover Local Features ofp(·)
One of the central ideas of this paper is that in a non-parametric setting (without parametric
constraints on r), we have an asymptotic formula (as the noise level σ →0) for the optimal
reconstruction function for the DAE and RCAE that allows us to recover the score ∂log p(x)
∂x .
A DAE is trained with a method that knows nothing about p, except through the use of
training samples to minimize a loss function, so it comes as a surprise that we can compute the
score of p at any point x.
In the following subsections we explore the consequences and the practical aspect of this.
3.1 Empirical Loss
In an experimental setting, the expected loss (7) is replaced by the empirical loss
ˆL= 1
N
N∑
n=1
(r(x(n)) −x(n)

2
2
+ σ2

∂r(x)
∂x
⏐⏐⏐⏐
x=x(n)

2
F
)
based on a sample
{
x(n)}N
n=1 drawn from p(x).
Alternatively, the auto-encoder is trained online (by stochastic gradient updates) with a
stream of examples x(n), which corresponds to performing stochastic gradient descent on the
expected loss (7). In both cases we obtain an auto-encoder that approximately minimizes the
expected loss.
An interesting question is the following: what can we infer from the data generating density
when given an auto-encoder reconstruction function r(x)?
The premise is that this auto-encoder r(x) was trained to approximately minimize a loss
function that has exactly the form of (7) for some σ2 > 0. This is assumed to have been
done through minimizing the empirical loss and the distribution pwas only available indirectly
through the samples
{
x(n)}N
n=1. We do not have access to p or to the samples. We have only
r(x) and maybe σ2.
7
We will now discuss the usefulness of r(x) based on diﬀerent conditions such as the model
capacity and the value of σ2.
3.2 Perfect World Scenario
As a starting point, we will assume that we are in a perfect situation, i.e., with no constraint
on r (non-parametric setting), an inﬁnite amount of training data, and a perfect minimization.
We will see what can be done to recover information about pin that ideal case. Afterwards, we
will drop certain assumptions one by one and discuss the possible paths to getting back some
information about p.
We use notation rσ(x) when we want to emphasize the fact that the value of r(x) came
from minimizing the loss with a certain ﬁxed σ.
Suppose that rσ(x) was trained with an inﬁnite sample drawn from p. Suppose also that it
had inﬁnite (or suﬃcient) model capacity and that it is capable of achieving the minimum of
the loss function (7) while satisfying the requirement that r(x) be twice diﬀerentiable. Suppose
that we know the value of σ and that we are working in a computing environment of arbitrary
precision (i.e. no rounding errors).
As shown by Theorem 1 and Theorem 2, we would be able to get numerically the values of
∂log p(x)
∂x at any point x∈Rd by simply evaluating
rσ(x) −x
σ2 →∂log p(x)
∂x as σ→0. (9)
In the setup described, we do not get to pick values ofσso as to take the limit σ→0. However,
it is assumed that σ is already suﬃciently small that the above quantity is close to ∂log p(x)
∂x for
all intents and purposes.
3.3 Simple Numerical Example
To give an example of this in one dimension, we will show what happens when we train a
non-parametric model ˆr(x) to minimize numerically the loss relative to p(x). We train both a
DAE and an RCAE in this fashion by minimizing a discretized version of their losses deﬁned
by equations (2) and (6). The goal here is to show that, for either a DAE or RCAE, the
approximation of the score that we get through equation (5) gets arbitrarily close to the actual
score ∂
∂x log p(x) as σ→0.
The distribution p(x) studied is shown in Figure 3 (left) and it was created to be simple
enough to illustrate the mechanics. We plot p(x) in Figure 3 (left) along with the score of p(x)
(right).
The model ˆr(x) is ﬁtted by dividing the interval [ −1.5,1.5] into M = 1000 partition points
x1,...,x M evenly separated by a distance ∆. The discretized version of the RCAE loss function
is
M∑
i=1
p(xi)∆ (ˆr(xi) −xi)2 + σ2
M−1∑
i=1
p(xi)∆
(ˆr(xi+1) −ˆr(xi)
∆
)2
. (10)
Every value ˆr(xi) for i = 1,...,M is treated as a free parameter. Setting to 0 the derivative
with respect to the ˆr(xi) yields a system of linear equations in M unknowns that we can solve
exactly. From that RCAE solution ˆr we get an approximation of the score of p at each point
xi. A similar thing can be done for the DAE by using a discrete version of the exact solution
(3) from Theorem 1. We now have two ways of approximating the score of p.
8
(a) p(x) = 1
Z exp(−E(x))
 (b) ∂
∂x log p(x) = − ∂
∂xE(x)
Figure 3: the density p(x) and its score for a simple one-dimensional example.
In Figure 4 we compare the approximations to the actual score of p for decreasingly small
values of σ∈{1.00,0.31,0.16,0.06}.
3.4 Vector Field Around a Manifold
We extend the experimentation of section 3.3 to a 1-dimensional manifold in 2-D space, in which
one can visualize r(x)−xas a vector ﬁeld, and we go from the non-parametric estimator of the
previous section to an actual auto-encoder trained by numerically minimizing the regularized
reconstruction error.
Two-dimensional data points (x,y) were generated along a spiral according to the following
equations:
x= 0.04 sin(t), y = 0.04 cos(t), t ∼Uniform (3,12) .
A denoising auto-encoder was trained with Gaussian corruption noiseσ= 0.01. The encoder
is f(x) = tanh( b+ Wx) and the decoder is g(h) = c+ Vh. The parameters ( b,c,V,W ) are
optimized by BFGS to minimize the average squared error, using a ﬁxed training set of 10 000
samples (i.e. the same corruption noises were sampled once and for all). We found better
results with untied weights, and BFGS gave more accurate models than stochastic gradient
descent. We used 1000 hiddens units and ran BFGS for 1000 iterations.
The non-convexity of the problem makes it such that the solution found depends on the
initialization parameters. The random corruption noise used can also inﬂuence the ﬁnal out-
come. Moreover, the fact that we are using a ﬁnite training sample size with reasonably small
noise may allow for undesirable behavior of r in regions far away from the training samples.
For those reasons, we trained the model multiple times and selected two of the most visually
appealing outcomes. These are found in Figure 5 which features a more global perspective
along with a close-up view.
Figure 5 shows the data along with the learned score function (shown as a vector ﬁeld). We
see that that the vector ﬁeld points towards the nearest high-density point on the data manifold.
The vector ﬁeld is close to zero near the manifold (i.e. the reconstruction error is close to zero),
also corresponding to peaks of the implicitly estimated density. The points on the manifolds
9
Figure 4: Comparing the approximation of the score of pgiven by discrete versions of optimally
trained auto-encoders with inﬁnite capacity. The approximations given by the RCAE
are in orange while the approximations given by the DAE are in purple. The results
are shown for decreasing values of σ∈{1.00,0.31,0.16,0.06}that have been selected
for their visual appeal.
As expected, we see in that the RCAE (orange) and DAE (purple) approximations
of the score are close to each other as predicted by Proposition 1. Moreover, they are
also converging to the true score (green) as predicted by Theorem 1 and Theorem 2.
play the role of sinks for the vector ﬁeld. Other places where reconstruction error may be low,
but where the implicit density is not high, are sources of the vector ﬁeld. In Figure 5(b) we
can see that we have that kind of behavior halfway between two sections of the manifold. This
shows that reconstruction error plays a very diﬀerent role as what was previously hypothesized:
whereas Ranzato et al. (2008) viewed reconstruction error as an energy function, our analysis
suggests that in regularized auto-encoders, it is the norm of an approximate score, i.e., the
derivative of the energy w.r.t. input. Note that the norm of the score should be small near
training examples (corresponding to local maxima of density) but it could also be small at
other places corresponding to local minima of density . This is indeed what happens in the
spiral example shown. It may happen whenever there are high-density regions separated by
a low-density region: tracing paths from one high-density region to another should cross a
“median” lower-dimensional region (a manifold) where the density has a local maximum along
10
(a) r(x) − x vector ﬁeld, acting as sink, zoomed out
 (b) r(x) − x vector ﬁeld, close-up
Figure 5: The original 2-D data from the data generating density p(x) is plotted along with the
vector ﬁeld deﬁned by the values ofr(x)−xfor trained auto-encoders (corresponding
to the estimation of the score ∂log p(x)
∂x ).
the path direction. The reason such a median region is needed is because at these points the
vectors r(x) −x must change sign: on one side of the median they point to one of the high-
density regions while on the other side they point to the other, as clearly visible in Figure 5(b)
between the arms of the spiral.
We believe that this analysis is valid not just for contractive and denoising auto-encoders,
but for regularized auto-encoders in general. The intuition behind this statement can be ﬁrmed
up by analyzing Figure 2: the score-like behavior of r(x) −xarises simply out of the opposing
forces of (a) trying to make r(x) = x at the training examples and (b) trying to make r(x) as
regularized as possible (as close to a constant as possible).
Note that previous work (Rifai et al., 2012; Bengio et al., 2013b) has already shown that
contractive auto-encoders (especially when they are stacked in a way similar to RBMs in a
deep belief net) learn good models of high-dimensional data (such as images), and that these
models can be used not just to obtain good representations for classiﬁcation tasks but that
good quality samples can be obtained from the model, by a random walk near the manifold of
high-density. This was achieved by essentially following the vector ﬁeld and adding noise along
the way.
3.5 Missing σ2
When we are in the same setting as in section 3.2 but the value of σ2 is unknown, we can
modify (9) a bit and avoid dividing by σ2. That is, for a trained reconstruction function r(x)
given to us we just take the quantity r(x) −xand it should be approximatively the score up to
a multiplicative constant.
r(x) −x∝∂log p(x)
∂x
Equivalently, if one estimates the density via an energy function (minus the unnormalized log
density), then x−r(x) estimates the derivative of the energy function.
11
We still have to assume that σ2 is small. Otherwise, if the unknown σ2 is too large we
might get a poor estimation of the score.
3.6 Limited Parameterization
We should also be concerned about the fact that r(x) −x is trying to approximate −∂E(x)
∂x as
σ→0 but we have not made any assumptions about the space of functions thatrcan represent
when we are dealing with a speciﬁc implementation.
When using a certain parameterization of r such as the one from section 3.3, there is
no guarantee that the family of functions in which we select r each represent a conservative
vector ﬁeld (i.e. the gradient of a potential function). Even if we start from a density p(x) ∝
exp(−E(x)) and we have that r(x) −xis very close to −∂
∂xE(x) in terms of some given norm,
there is not guarantee that there exists an associated function E0(x) for which r(x) −x ∝
−∂
∂xE0(x) and E0(x) ≈E(x).
In fact, in many cases we can trivially show the non-existence of such a E0(x) by computing
the curl of r(x). The curl has to be equal to 0 everywhere if r(x) is indeed the derivative of a
potential function. We can omit the xterms from the computations because we can easily ﬁnd
its antiderivative by looking at x= ∂
∂x ∥x∥2
2.
Conceptually, another way to see this is to argue that if such a function E0(x) existed, its
second-order mixed derivatives should be equal. That is, we should have that
∂2E0(x)
∂xi∂xj
= ∂2E0(x)
∂xj∂xi
∀i,j,
which is equivalent to
∂ri(x)
∂xj
= ∂rj(x)
∂xi
∀i,j.
Again in the context of section 3.3, with the parameterization used for that particual kind
of denoising auto-encoder, this would yield the constraint that VT = W. That is, unless we
are using tied weights, we know that no such potential E0(x) exists, and yet when running the
experiments from section 3.3 we obtained much better results with untied weights. To make
things worse, it can also be demonstrated that the energy function that we get from tied weights
leads to a distribution that is not normalizable (it has a divergent integral over Rd). In that
sense, this suggests that we should not worry too much about the exact parameterization of
the denoising auto-encoder as long as it has the required ﬂexibility to approximate the optimal
reconstruction function suﬃciently well.
3.7 Relation to Denoising Score Matching
There is a connection between our results and previous research involving score matching for
denoising auto-encoders. We will summarize here the existing results from Vincent (2011)
and show that, whereas they have shown that denoising auto-encoders with a particular form
estimated the score, our results extend this to a very large family of estimators (including the
non-parametric case). This will provide some reassurance given some of the potential issues
mentioned in section 3.6.
Motivated by the analysis of denoising auto-encoders, Vincent (2011) are concerned with
the case where we explicitly parametrize an energy function E(x), yielding an associated score
function ψ(x) = −∂E(x)
∂x and we stochastically corrupt the original samples x ∼p to obtain
12
noisy samples ˜x∼qσ(˜x|x). In particular, the article analyzes the case where qσ adds Gaussian
noise of variance σ2 to x. The main result is that minimizing the expected square diﬀerence
between ψ(˜x) and the score of qσ(˜x|x),
Ex,˜x[||ψ(˜x) −∂log qσ(˜x|x)
∂˜x ||2],
is equivalent to performing score matching (Hyv¨ arinen, 2005) with estimatorψ(˜x) and target
density qσ(˜x) =
∫
qσ(˜x|x)p(x)dx, where p(x) generates the training samples x. Note that when
a ﬁnite training set is used, qσ(˜x) is simply a smooth of the empirical distribution (e.g. the
Parzen density with Gaussian kernel of width σ). When the corruption noise is Gaussian,
qσ(˜x|x)
∂˜x = x−˜x
σ2 , from which we can deduce that if we deﬁne a reconstruction function
r(˜x) = ˜x+ σ2ψ(˜x), (11)
then the above expectation is equivalent to
Ex,˜x[||r(˜x) −˜x
σ2 −x−˜x
σ2 ||2] = 1
σ2 Ex,˜x[||r(˜x) −x||2]
which is the denoising criterion. This says that when the reconstruction functionris parametrized
so as to correspond to the score ψof a model density (as per eq. 11, and where ψis a derivative
of some log-density), the denoising criterion onrwith Gaussian corruption noise is equivalent to
score matching with respect to a smooth of the data generating density, i.e., a regularized form
of score matching. Note that this regularization appears desirable, because matching the score
of the empirical distribution (or an insuﬃciently smoothed version of it) could yield undesirable
results when the training set is ﬁnite. Since score matching has been shown to be a consistent
induction principle (Hyv¨ arinen, 2005), it means that this denoising score matching (Vincent,
2011; Kingma and LeCun, 2010; Swersky et al., 2011) criterion recovers the underlying density,
up to the smoothing induced by the noise of variance σ2. By making σ2 small, we can make
the estimator arbitrarily good (and we would expect to want to do that as the amount of
training data increases). Note the correspondance of this conclusion with the results presented
here, which show (1) the equivalence between the RCAE’s regularization coeﬃcient and the
DAE’s noise variance σ2, and (2) that minimizing the equivalent analytic criterion (based on
a contraction penalty) estimates the score when σ2 is small. The diﬀerence is that our result
holds even when r is not parametrized as per eq. 11, i.e., is not forced to correspond with the
score function of a density.
3.8 Estimating the Hessian
Since we have r(x)−x
σ2 as an estimator of the score, we readily obtain that the Hessian of the
log-density, can be estimated by the Jacobian of the reconstruction function minus the identity
matrix:
∂2 log p(x)
∂x2 ≈(∂r(x)
∂x −I)/σ2
as shown by equation (8) of Theorem 2.
In spite of its simplicity, this result is interesting because it relates the derivative of the
reconstruction function, i.e., a Jacobian matrix, with the second derivative of the log-density
(or of the energy). This provides insights into the geometric interpretation of the reconstruction
13
function when the density is concentrated near a manifold. In that case, near the manifold the
score is nearly 0 because we are near a ridge of density, and the density’s second derivative
matrix tells us in which directions the ﬁrst density remains close to zero or increases. The
ridge directions correspond to staying on the manifold and along these directions we expect the
second derivative to be close to 0. In the orthogonal directions, the log-density should decrease
sharply while its ﬁrst and second derivatives would be large in magnitude and negative in
directions away from the manifold.
Returning to the above equation, keep in mind that in these derivationsσ2 is near 0 and r(x)
is near x, so that ∂r(x)
∂x is close to the identity. In particular, in the ridge (manifold) directions,
we should expect ∂r(x)
∂x to be closer to the identity, which means that the reconstruction remains
faithful (r(x) = x) when we move on the manifold, and this corresponds to the eigenvalues of
∂r(x)
∂x that are near 1, making the corresponding eigenvalues of ∂2 log p(x)
∂x2 near 0. On the other
hand, in the directions orthogonal to the manifold, ∂r(x)
∂x should be smaller than 1, making the
corresponding eigenvalues of ∂2 log p(x)
∂x2 negative.
Besides ﬁrst and second derivatives of the density, other local properties of the density are
its local mean and local covariance, discussed in the Appendix, section 6.4.
4. Sampling with Metropolis-Hastings
4.1 Estimating Energy Diﬀerences
One of the immediate consequences of equation (5) is that, while we cannot easily recover the
energy E(x) itself, it is possible to approximate the energy diﬀerence E(x∗) −E(x) between
two states x and x∗. This can be done by using a ﬁrst-order Taylor approximation
E(x∗) −E(x) = ∂E(x)
∂x
T
(x∗−x) + o(∥x∗−x∥).
To get a more accurate approximation, we can also use a path integral from x to x∗ that we
can discretize in suﬃciently many steps. With a smooth path γ(t) : [0,1] →Rd, assuming that
γ stays in a region where our DAE/RCAE can be used to approximate ∂E
∂x well enough, we
have that
E(x∗) −E(x) =
∫ 1
0
[∂E
∂x (γ(t))
]T
γ′(t)dt. (12)
The simplest way to discretize this path integral is to pick points {xi}n
i=1 spread at even
distances on a straight line from x1 = x to xn = x∗. We approximate (12) by
E(x∗) −E(x) ≈1
n
n∑
i=1
[∂E
∂x (xi)
]T
(x∗−x) (13)
4.2 Sampling
With equation (12) from section 4.1 we can perform approximate sampling from the estimated
distribution, using the score estimator to approximate energy diﬀerences which are needed
in the Metropolis-Hastings accept/reject decision. Using a symmetric proposal q(x∗|x), the
acceptance ratio is
α= p(x∗)
p(x) = exp(−E(x∗) + E(x))
14
which can be computed with (12) or approximated with (13) as long as we trust that our
DAE/RCAE was trained properly and has enough capacity to be a suﬃciently good estimator
of ∂E
∂x. An example of this process is shown in Figure 6 in which we sample from a density
concentrated around a 1-d manifold embedded in a space of dimension 10. For this particular
task, we have trained only DAEs and we are leaving RCAEs out of this exercise. Given that the
data is roughly contained in the range [ −1.5,1.5] along all dimensions, we selected a training
noise level σtrain = 0 .1 so that the noise would have an appreciable eﬀect while still being
relatively small. As required by Theorem 1, we have used isotropic Gaussian noise of variance
σ2
train.
The Metropolis-Hastings proposal q(x∗|x) = N(0,σ2
MHI) has a noise parameter σMH that
needs to be set. In the situation shown in Figure 6, we used σMH = 0.1. After some hyperpa-
rameter tweaking and exploring various scales for σtrain,σMH, we found that setting both to be
0.1 worked well.
When σtrain is too large, the DAE trained learns a “blurry” version of the density that
fails to represent the details that we are interested in. The samples shown in Figure 6 are very
convincing in terms of being drawn from a distribution that models well the original density. We
have to keep in mind that Theorem 1 describes the behavior as σtrain →0 so we would expect
that the estimator becomes worse when σtrain is taking on larger values. In this particular case
with σtrain = 0.1, it seems that we are instead modeling something like the original density to
which isotropic Gaussian noise of variance σ2
train has been added.
In the other extreme, when σtrain is too small, the DAE is not exposed to any training
example farther away from the density manifold. This can lead to various kinds of strange
behaviors when the sampling algorithm falls into those regions and then has no idea what to do
there and how to get back to the high-density regions. We come back to that topic in section
4.3.
It would certainly be possible to pick both a very small value for σtrain = σMH = 0.01 to
avoid the spurius maxima problem illustrated in section 4.3. However, this leads to the same
kind of mixing problems that any kind of MCMC algorithm has. Smaller values of σMH lead
to higher acceptance ratios but worse mixing properties.
4.3 Spurious Maxima
There are two very real concerns with the sampling method discussed in section 4.2. The
ﬁrst problem is with the mixing properties of MCMC and it is discussed in that section. The
second issue is with spurious probability maxima resulting from inadequate training of the
DAE. It happens when an auto-encoder lacks the capacity to model the density with enough
precision, or when the training procedure ends up in a bad local minimum (in terms of the
DAE parameters).
This is illustrated in Figure 7 where we show an example of a vector ﬁeld r(x) −x for a
DAE that failed to properly learn the desired behavior in regions away from the spiral-shaped
density.
5. Conclusion
Whereas auto-encoders have long been suspected of capturing information about the data
generating density, this work has clariﬁed what some of them are actually doing, showing that
they can actually implicitly recover the data generating density altogether. We have shown that
15
original sampled original sampled
Figure 6: Samples drawn from the estimate of ∂E
∂x given by a DAE by the Metropolis-
Hastings method presented in section 4. By design, the data density distribu-
tion is concentrated along a 1-d manifold embedded in a space of dimension 10.
This data can be visualized in the plots above by plotting pairs of dimensions
(x0,x1),..., (x8,x9),(x9,x0), going in reading order from left to right and then line
by line. For each pair of dimensions, we show side by side the original data (left)
with the samples drawn (right).
regularized auto-encoders such as the denoising auto-encoder and a form of contractive auto-
encoder are closely related to each other and estimate local properties of the data generating
density: the ﬁrst derivative (score) and second derivative of the log-density, as well as the local
mean. This contradicts the previous interpretation of reconstruction error as being an energy
function (Ranzato et al., 2008) but is consistent with our experimental ﬁndings. Our results do
not require the reconstruction function to correspond to the derivative of an energy function as
in Vincent (2011), but hold simply by virtue of minimizing the regularized reconstruction error
training criterion. This suggests that minimizing a regularized reconstruction error may be an
alternative to maximum likelihood for unsupervised learning, avoiding the need for MCMC in
the inner loop of training, as in RBMs and deep Boltzmann machines, analogously to score
16
(a) DAE misbehaving when away from manifold
 (b) sampling getting trapped into bad attractor
Figure 7: (a) On the left we show a r(x) −xvector ﬁeld similar to that of the earlier Figure 5.
The density is concentrated along a spiral manifold and we should have the recon-
struction function r bringing us back towards the density. In this case, it works well
in the region close to the spiral (the magnitude of the vectors is so small that the
arrows are shown as dots). However, things are out of control in the regions outside.
This is because the level of noise used during training was so small that not enough
of the training examples were found in those regions.
(b) On the right we sketch what may happen when we follow a sampling procedure
as described in section 4.2. We start in a region of high density (in purple) and we
illustrate in red the trajectory that our samples may take. In that situation, the
DAE/RCAE was not trained properly. The resulting vector ﬁeld does not reﬂect the
density accurately because it should not have this attractor (i.e. stable ﬁxed point)
outside of the manifold on which the density is concentrated. Conceptually, the sam-
pling procedure visits that spurious attractor because it assumes that it corresponds
to a region of high probability. In some cases, this eﬀect is regrettable but not catas-
trophic, but in other situations we may end up with completely unusable samples.
In the experiments, training with enough of the examples involving suﬃciently large
corruption noise typically eliminates that problem.
matching (Hyv¨ arinen, 2005; Vincent, 2011). Toy experiments have conﬁrmed that a good
estimator of the density can be obtained when this criterion is non-parametrically minimized.
The experiments have also conﬁrmed that an MCMC could be setup that approximately samples
from the estimated model, by estimating energy diﬀerences to ﬁrst order (which only requires
the score) to perform approximate Metropolis-Hastings MCMC.
Many questions remain open and deserve futher study. A big question is how to general-
ize these ideas to discrete data, since we have heavily relied on the notions of scores, i.e., of
derivatives with respect to x. A natural extension of the notion of score that could be ap-
plied to discrete data is the notion of relative energy, or energy diﬀerence between a point x
and a perturbation ˜x of x. This notion has already been successfully applied to obtain the
equivalent of score matching for discrete models, namely ratio matching (Hyv¨ arinen, 2007).
More generally, we would like to generalize to any form of reconstruction error (for example
17
many implementations of auto-encoders use a Bernouilli cross-entropy as reconstruction loss
function) and any (reasonable) form of corruption noise (many implementations use masking
or salt-and-pepper noise, not just Gaussian noise). More fundamentally, the need to rely on
σ →0 is troubling, and getting rid of this limitation would also be very useful. A possible
solution to this limitation, as well as adding the ability to handle both discrete and continuous
variables, has recently been proposed while this article was under review (Bengio et al., 2013a).
It would also be interesting to generalize the results presented here to other regularized
auto-encoders besides the denoising and contractive types. In particular, the commonly used
sparse auto-encoders seem to ﬁt the qualitative pattern illustrated in Figure 2 where a score-like
vector ﬁeld arises out of the opposing forces of minimizing reconstruction error and regularizing
the auto-encoder.
We have mostly considered the harder case where the auto-encoder parametrization does
not guarantee the existence of an analytic formulation of an energy function. It would be
interesting to compare experimentally and study mathematically these two formulations to
assess how much is lost (because the score function may be somehow inconsistent) or gained
(because of the less constrained parametrization).
Acknowledgements
The authors thank Salah Rifai, Max Welling, Yutian Chen and Pascal Vincent for fruitful
discussions, and acknowledge the funding support from NSERC, Canada Research Chairs and
CIFAR.
References
Yoshua Bengio. Learning deep architectures for AI. Foundations and Trends in Machine Learning , 2
(1):1–127, 2009. Also published as a book. Now Publishers, 2009.
Yoshua Bengio and Olivier Delalleau. On the expressive power of deep architectures. In ALT’2011,
2011.
Yoshua Bengio, Pascal Lamblin, Dan Popovici, and Hugo Larochelle. Greedy layer-wise training of deep
networks. In NIPS 19, pages 153–160. MIT Press, 2007.
Yoshua Bengio, Guillaume Alain, and Salah Rifai. Implicit density estimation by local moment matching
to sample from auto-encoders. Technical report, arXiv:1207.0057, 2012a.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new
perspectives. Technical report, arXiv:1206.5538, 2012b.
Yoshua Bengio, Yao Li, Guillaume Alain, and Pascal Vincent. Generalized denoising auto-encoders as
generative models. Technical Report arXiv:1305.6663, Universite de Montreal, 2013a.
Yoshua Bengio, Gr´ egoire Mesnil, Yann Dauphin, and Salah Rifai. Better mixing via deep representations.
In ICML’13, 2013b.
Lawrence Cayton. Algorithms for manifold learning. Technical Report CS2008-0923, UCSD, 2005.
B. Dacorogna. Introduction to the Calculus of Variations . World Scientiﬁc Publishing Company, 2004.
Karol Gregor, Arthur Szlam, and Yann LeCun. Structured sparse coding via lateral inhibition. In
Advances in Neural Information Processing Systems (NIPS 2011) , volume 24, 2011.
18
Geoﬀrey E. Hinton, Simon Osindero, and Yee Whye Teh. A fast learning algorithm for deep belief nets.
Neural Computation, 18:1527–1554, 2006.
Aapo Hyv¨ arinen. Estimation of non-normalized statistical models using score matching. Journal of
Machine Learning Research, 6:695–709, 2005.
Aapo Hyv¨ arinen. Some extensions of score matching. Computational Statistics and Data Analysis , 51:
2499–2512, 2007.
Viren Jain and Sebastian H. Seung. Natural image denoising with convolutional networks. In NIPS’08,
pages 769–776, 2008.
Koray Kavukcuoglu, Marc’Aurelio Ranzato, Rob Fergus, and Yann LeCun. Learning invariant features
through topographic ﬁlter maps. In Proceedings of the Computer Vision and Pattern Recognition
Conference (CVPR’09), pages 1605–1612. IEEE, 2009.
Diederik Kingma and Yann LeCun. Regularized estimation of image statistics by score matching. In
J. Laﬀerty, C. K. I. Williams, J. Shawe-Taylor, R.S. Zemel, and A. Culotta, editors, Advances in
Neural Information Processing Systems 23 , pages 1126–1134, 2010.
Honglak Lee, Roger Grosse, Rajesh Ranganath, and Andrew Y. Ng. Convolutional deep belief net-
works for scalable unsupervised learning of hierarchical representations. In L´ eon Bottou and Michael
Littman, editors, ICML 2009. ACM, Montreal (Qc), Canada, 2009.
Hariharan Narayanan and Sanjoy Mitter. Sample complexity of testing the manifold hypothesis. In
NIPS’2010. 2010.
B. A. Olshausen and D. J. Field. Sparse coding with an overcomplete basis set: a strategy employed
by V1? Vision Research, 37:3311–3325, December 1997. URL http://view.ncbi.nlm.nih.gov/
pubmed/9425546.
Marc’Aurelio Ranzato, Christopher Poultney, Sumit Chopra, and Yann LeCun. Eﬃcient learning of
sparse representations with an energy-based model. In NIPS’06, pages 1137–1144. MIT Press, 2007.
Marc’Aurelio Ranzato, Y-Lan Boureau, and Yann LeCun. Sparse feature learning for deep belief net-
works. In NIPS’07, pages 1185–1192, Cambridge, MA, 2008. MIT Press.
Salah Rifai, Yann Dauphin, Pascal Vincent, Yoshua Bengio, and Xavier Muller. The manifold tangent
classiﬁer. In NIPS’2011, 2011a. Student paper award.
Salah Rifai, Pascal Vincent, Xavier Muller, Xavier Glorot, and Yoshua Bengio. Contractive auto-
encoders: Explicit invariance during feature extraction. In ICML’11, 2011b.
Salah Rifai, Pascal Vincent, Xavier Muller, Xavier Glorot, and Yoshua Bengio. Contractive auto-
encoders: Explicit invariance during feature extraction. In ICML’2011, 2011c.
Salah Rifai, Yoshua Bengio, Yann Dauphin, and Pascal Vincent. A generative process for sampling
contractive auto-encoders. In ICML’12, 2012.
R. Salakhutdinov and G.E. Hinton. Deep Boltzmann machines. In Proceedings of the Twelfth Interna-
tional Conference on Artiﬁcial Intelligence and Statistics (AISTATS 2009) , volume 8, 2009.
Kevin Swersky, Marc’Aurelio Ranzato, David Buchman, Benjamin Marlin, and Nando de Freitas. On
autoencoders and score matching for energy based models. In ICML’2011. ACM, 2011.
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Computation,
23(7):1661–1674, 2011.
19
Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and Pierre-Antoine Manzagol. Extracting and com-
posing robust features with denoising autoencoders. In ICML’08, pages 1096–1103. ACM, 2008.
6. Appendix
6.1 Optimal DAE solution
Theorem 1 Let p be the probability density function of the data. If we train a DAE using the
expected quadratic loss and corruption noise N(x) = x+ ϵ with
ϵ∼N
(
0,σ2I
)
,
then the optimal reconstruction function r∗(x) will be given by
r∗(x) = Eϵ[p(x−ϵ)(x−ϵ)]
Eϵ[p(x−ϵ)] (14)
for values of x where p(x) ̸= 0.
Moreover, if we consider how the optimal reconstruction function r∗
σ(x) behaves asymptoti-
cally as σ→0, we get that
r∗
σ(x) = x+ σ2 ∂log p(x)
∂x + o(σ2) as σ→0.
Proof The ﬁrst part of this proof is to get to equation (14) without assuming that σ→0.
By using an auxiliary variable ˜x = x+ ϵ, we can rewrite this loss in a way that puts the
quantity r(˜x) in focus and allows us to perform the minimization with respect to each choice
of r(˜x) independantly. That is, we have that
LDAE =
∫
Rd
Eϵ∼N(0,σ2I)
[
p(˜x−ϵ) ∥r(˜x) −˜x+ ϵ∥2
2
]
d˜x
which can be diﬀerentiated with respect to the quantity r(˜x) and set to be equal to 0. Denoting
the optimum by r∗(˜x), we get
0 = Eϵ∼N(0,σ2I) [p(˜x−ϵ) (r∗(˜x) −˜x+ ϵ)]
Eϵ∼N(0,σ2I) [p(˜x−ϵ)r∗(˜x)] = Eϵ∼N(0,σ2I) [p(˜x−ϵ)(˜x−ϵ)]
r∗(˜x) = Eϵ∼N(0,σ2I) [p(˜x−ϵ)(˜x−ϵ)]
Eϵ∼N(0,σ2I) [p(˜x−ϵ)] . (15)
We used ˜x out of convenience in equation (15). Theorem 1 is just nicer to read when using
r(x), and we switch back to using r(x) for the rest of this proof.
Now, for the second part of this proof, we study the behavior of the solution r∗
σ(x) as σ
approaches 0. We start with equation (15) that we rewrite in a way to be able to pull out the
leading term x. Given the symmetry of the distribution of ϵ∼N
(
0,σ2I
)
, we can also simplify
equation (15) by converting all the −ϵ into +ϵ.
20
Eϵ[p(x+ ϵ)(x+ ϵ)]
Eϵ[p(x+ ϵ)] = Eϵ[p(x+ ϵ)x]
Eϵ[p(x+ ϵ)] + Eϵ[p(x+ ϵ)ϵ]
Eϵ[p(x+ ϵ)]
= x+ Eϵ[p(x+ ϵ)ϵ]
Eϵ[p(x+ ϵ)]
We now look at the Taylor expansion of p(x+ ϵ) around x. We perform this expansion
inside of the expectation, where ϵ just represents as small quantity of scale σ.
p(x+ ϵ) = p(x) + ∇p(x)Tϵ+ 1
2ϵT ∂2p(x)
∂x2 ϵ+ o(ϵ2) as ϵ→0 (16)
When taking the expectation of p(x+ ϵ) with respect to ϵ, we get a zero for all the terms
containing an odd power of ϵ, for reasons of symmetry. When we take the expectation of
p(x+ ϵ)ϵ instead, all the terms with an even power of ϵ in the above expectation will vanish.
Thus, we get that
Eϵ[p(x+ ϵ)ϵ] = Eϵ[ϵϵT]∇p(x) + Eϵ[o(ϵ2)] = σ2∇p(x) + o(σ2)
Eϵ[p(x+ ϵ)] = p(x) + O(σ2)
Provided that p(x) ̸= 0, our quotient can be rewritten as
σ2∇p(x) + o(σ2)
p(x) + O(σ2) = σ2∇log p(x) + o(σ2)
1 + O(σ2) . (17)
We can use the basic geometric series expansion to write
1
1 + O(σ2) = 1 −O(σ2) = 1 + O(σ2)
so that equation (17) is now
σ2∇log p(x) + o(σ2)
1 + O(σ2) = σ2∇log p(x) + o(σ2).
This is the result that we wanted.
Note that p(x) was treated as a constant when we studied the asymptotic behavior as
σ →0. In the Taylor expansion around some given x, we want to stuﬀ all the higher-order
derivatives of p(x) into the asymptotic remainder term. It is quite possible that the size of the
σ required to do so would depend on the particular x, that there would not be a uniform σ >0
suitable for all the x. Therefore we can only say that we are dealing with pointwise convergence
(not uniform convergence) in our formula for the asymptotic expansion of r∗
σ(x). For practical
applications, we do not need more than that.
The assumption that p(x) ̸= 0 could also be viewed as problematic, but if we go back to
the deﬁnition of the DAE loss, then any point where p(x) = 0 would never produce a training
example and would not even contribute in the deﬁnition of the expectation of DAE loss. The
correct value to assign to r at such a point is not well-deﬁned.
21
6.2 Relationship between Contractive Penalty and Denoising Criterion
Proposition 1 Let p be the probability density function of the data. Consider a DAE using
the expected quadratic loss and corruption noise N(x) = x+ ϵ, with ϵ ∼N
(
0,σ2I
)
. If we
assume that the non-parametric solutions rσ(x) satistﬁes
rσ(x) = x+ o(1) as σ→0,
then we can rewrite the loss as
LDAE (rσ) = E
[
∥rσ(x) −x∥2
2 + σ2

∂rσ(x)
∂x

2
F
]
+ o(σ2) as σ→0
where the expectation is taken with respect to X, whose distribution is given by p.
Proof We can drop the σ index from rσ to lighten the notation if we just keep in mind that
we are considering the particular family of solutions such that r(x) −x= o(1). With a Taylor
expansion around x we have that
r(x+ ϵ) = r(x) + ∂r(x)
∂x ϵ+ O(ϵ2).
The DAE loss involves taking the expectation with respect to X and with respect to the noise
ϵ. We substitute the Taylor expansion into LDAE, we express the norm as a dot product, and
we show how the expectation with respect to ϵ cancels out certain terms.
EXEϵ
[x−
(
r(x) + ∂r(x)
∂x ϵ+ O(ϵ2)
)
2
2
]
= EX
[
Eϵ
[
∥x−r(x)∥2
2
]
−2 (x−r(x))T Eϵ
[∂r(x)
∂x ϵ+ O(ϵ2)
]
+ Eϵ
[
ϵT ∂r(x)
∂x
T ∂r(x)
∂x ϵ+ O(ϵ3)
]]
= EX
[
[
∥r(x) −x∥2
2
]
+ o(1)O(σ2) + σ2Tr
(
∂r(x)
∂x
T ∂r(x)
∂x
)]
+ O(σ3)
= EX
[
∥r(x) −x∥2
2
]
+ o(σ2) + σ2EX
[
∂r(x)
∂x

2
F
]
The fact that we get a trace (rewritten as Frobenius norm) comes from the assumption
that the noise is Gaussian and is independant from X. This explains why E
[
ϵϵT]
= σ2I and
E[ϵ] = 0. It also motivates many of our substitutions such as Eϵ
[
O(ϵk)
]
= O(σk).
6.3 Calculus of Variations
Theorem 2 Let pbe a probability density function that is continuously diﬀerentiable once and
with support Rd (i.e. ∀x∈Rd we have p(x) ̸= 0). Let Lσ be the loss function deﬁned by
Lσ(r) =
∫
Rd
p(x)
[
∥r(x) −x∥2
2 + σ2

∂r(x)
∂x

2
F
]
dx
22
for r: Rd →Rd assumed to be diﬀerentiable twice, and 0 ≤σ∈R used as factor to the penalty
term.
Let r∗
σ(x) denote the optimal function that minimizes Lσ. Then we have that
r∗
σ(x) = x+ σ2 ∂log p(x)
∂x + o(σ2) as σ→0.
Moreover, we also have the following expression for the derivative
∂r∗
σ(x)
∂x = I+ σ2 ∂2 log p(x)
∂x2 + o(σ2) as σ→0.
Both these asymptotic expansions are to be understood in a context where we consider{r∗
σ(x)}σ≥0
to be a family of optimal functions minimizing Lσ for their corresponding value of σ. The
asymptotic expansions are applicable point-wise in x, that is, with any ﬁxed x we look at the
behavior as σ→0.
Proof
This proof is done in two parts. In the ﬁrst part, the objective is to get to equation (20)
that has to be satisﬁed for the optimum solution.
We will leave out the σ indices from the expressions involving r(x) to make the notation
lighter. We have a more important need for indices k in rk(x) that denote the d components
of r(x) ∈Rd.
We treat σ as given and constant for the ﬁrst part of this proof.
In the second part we work out the asymptotic expansion in terms of σ. We again work
with the implicit dependence of r(x) on σ.
(part 1 of the proof)
We make use of the Euler-Lagrange equation from the Calculus of Variations. We would
refer the reader to either (Dacorogna, 2004) or Wikipedia for more on the topic. Let
f(x1,...,x n,r,rx1 ,...,r xn) = p(x)
[
∥r(x) −x∥2
2 + σ2

∂r(x)
∂x

2
F
]
where x= (x1,...,x d) , r(x) = (r1(x),...,r d(x)) and rxi = ∂f
∂xi
.
We can rewrite the loss L(r) more explicitly as
L(r) =
∫
Rd
p(x)


d∑
i=1
(ri(x) −xi)2 + σ2
d∑
i=1
d∑
j=1
∂ri(x)
∂xj
2

dx
=
d∑
i=1
∫
Rd
p(x)

(ri(x) −xi)2 + σ2
d∑
j=1
∂ri(x)
∂xj
2

dx (18)
to observe that the components r1(x),...,r d(x) can each be optimized separately.
The Euler-Lagrange equation to be satisﬁed at the optimal r: Rd →Rd is
∂f
∂r =
d∑
i=1
∂
∂xi
∂f
∂rxi
.
23
In our situation, the expressions from that equation are given by
∂f
∂r = 2(r(x) −x)p(x)
∂f
∂rxi
= 2σ2p(x)
[
∂r1
∂xi
∂r2
∂xi
··· ∂rd
∂xi
]T
∂
∂xi
( ∂f
∂rxi
)
= 2 σ2 ∂p(x)
∂xi
[
∂r1
∂xi
∂r2
∂xi
··· ∂rd
∂xi
]T
+2σ2p(x)
[
∂2r1
∂x2
i
∂2r2
∂x2
i
··· ∂2rd
∂x2
i
]T
and the equality to be satisﬁed at the optimum becomes
(r(x) −x)p(x) = σ2
d∑
i=1


∂p(x)
∂xi
∂r1
∂xi
+ p(x)∂2r1
∂x2
i
...
∂p(x)
∂xi
∂rd
∂xi
+ p(x)∂2rd
∂x2
i

. (19)
As equation (18) hinted, the expression (19) can be decomposed into the diﬀerent compo-
nents rk(x) : Rd →R that make r. For k= 1,...,d we get
(rk(x) −xk)p(x) = σ2
d∑
i=1
(∂p(x)
∂xi
∂rk(x)
∂xi
+ p(x)∂2rk(x)
∂x2
i
)
.
As p(x) ̸= 0 by hypothesis, we can divide all the terms by p(x) and note that ∂p(x)
∂xi
/p(x) =
∂log p(x)
∂xi
.
We get
rk(x) −xk = σ2
d∑
i=1
(∂log p(x)
∂xi
∂rk(x)
∂xi
+ ∂2rk(x)
∂x2
i
)
. (20)
This ﬁrst thing to observe is that when σ2 = 0 the solution is just rk(x) = xk, which
translates into r(x) = x. This is not a surprise because it represents the perfect reconstruction
value that we get when we the penalty term vanishes in the loss function.
(part 2 of the proof)
This linear partial diﬀerential equation (20) can be used as a recursive relation for rk(x) to
obtain a Taylor series in σ2. The goal is to obtain an expression of the form
r(x) = x+ σ2h(x) + o(σ2) as σ2 →0 (21)
where we can solve for h(x) and for which we also have that
∂r(x)
∂x = I+ σ2 ∂h(x)
∂x + o(σ2) as σ2 →0.
24
We can substitute in the right-hand side of equation (21) the value for rk(x) that we get from
equation (21) itself. This substitution would be pointless in any other situation where we are
not trying to get a power series in terms of σ2 around 0.
rk(x) = xk +σ2
d∑
i=1
(∂log p(x)
∂xi
∂rk(x)
∂xi
+ ∂2rk(x)
∂x2
i
)
(22)
= xk +σ2
d∑
i=1

∂log p(x)
∂xi
∂
∂xi

xk + σ2
d∑
j=1
(
∂log p(x)
∂xj
∂rk(x)
∂xj
+ ∂2rk(x)
∂x2
j
)


(23)
+σ2
d∑
i=1
∂2rk(x)
∂x2
i
(24)
= xk +σ2
d∑
i=1
∂log p(x)
∂xi
I(i= k) + σ2
d∑
i=1
∂2rk(x)
∂x2
i
(25)
+(σ2)2
d∑
i=1
d∑
j=1
(
∂log p(x)
∂xi
∂
∂xi
(
∂log p(x)
∂xj
∂rk(x)
∂xj
+ ∂2rk(x)
∂x2
j
))
(26)
rk(x) = xk +σ2 ∂log p(x)
∂xk
+ σ2
d∑
i=1
∂2rk(x)
∂x2
i
+ σ4ρ(σ2,x) (27)
Now we would like to get rid of that σ2 ∑d
i=1
∂2rk(x)
∂x2
i
term by showing that it is a term that
involves only powers of σ4 or higher. We get this by showing what we get by diﬀerentiating
the expression for rk(x) in line (27) twice with respect to some l.
∂rk(x)
∂xl
= I(i= l) + σ2 ∂2 log p(x)
∂xl∂xk
+ σ2 ∂
∂xl
( d∑
i=1
∂2rk(x)
∂x2
i
+ σ2ρ(σ2,x)
)
∂2rk(x)
∂x2
l
= σ2 ∂3 log p(x)
∂x2
l∂xk
+ σ2 ∂
∂x2
l
( d∑
i=1
∂2rk(x)
∂x2
i
+ σ2ρ(σ2,x)
)
Since σ2 is a common factor in all the terms of the expression of ∂2rk(x)
∂x2
l
we get what we
needed. That is,
rk(x) = xk + σ2 ∂log p(x)
∂xk
+ σ4η(σ2,x).
This shows that
r(x) = x+ σ2 ∂log p(x)
∂x + o(σ2) as σ2 →0
and
∂r(x)
∂x = I+ σ2 ∂2 log p(x)
∂x2 + o(σ2) as σ2 →0
which completes the proof.
25
6.4 Local Mean
In preliminary work (Bengio et al., 2012a), we studied how the optimal reconstruction could
possibly estimate so-called local moments. We revisit this question here, with more appealing
and precise results.
What previous work on denoising and contractive auto-encoders suggest is that regularized
auto-encoders can capture the local structure of the density through the value of the encoding
(or reconstruction) function and its derivative. In particular, Rifai et al. (2012); Bengio et al.
(2012a) argue that the ﬁrst and second derivatives tell us in which directions it makes sense
to randomly move while preserving or increasing the density, which may be used to justify
sampling procedures. This motivates us here to study so-called local moments as captured
by the auto-encoder, and in particular the local mean, following the deﬁnitions introduced
in Bengio et al. (2012a).
6.4.1 Definitions for Local Distributions
Let p be a continuous probability density function with support Rd. That is, ∀x∈Rd we have
that p(x) ̸= 0. We deﬁne below the notion of a local ball Bδ(x0), along with an associated local
density, which is the normalized product of p with the indicator for the ball:
Bδ(x0) = {x s.t. ∥x−x0∥2 <δ}
Zδ(x0) =
∫
Bδ(x0)
p(x)dx
pδ(x|x0) = 1
Zδ(x0)p(x)I(x∈Bδ(x0))
where Zδ(x0) is the normalizing constant required to make pδ(x|x0) a valid pdf for a dis-
tribution centered on x0. The support of pδ(x|x0) is the ball of radius δ around x0 denoted
by Bδ(x0). We stick to the 2-norm in terms of deﬁning the balls Bδ(x0) used, but everything
could be rewritten in terms of another p-norm to have slightly diﬀerent formulas.
We use the following notation for what will be referred to as the ﬁrst two local moments
(i.e. local mean and local covariance) of the random variable described by pδ(x|x0).
mδ(x0)
def
=
∫
Rd
xpδ(x|x0)dx
Cδ(x0)
def
=
∫
Rd
(x−mδ(x0))(x−mδ(x0))Tpδ(x|x0)dx
Based on these deﬁnitions, one can prove the following theorem.
Theorem 3 Let p be of class C3 and represent a probability density function. Let x0 ∈Rd
with p(x0) ̸= 0. Then we have that
mδ(x0) = x0 + δ2 1
d+ 2
∂log p(x)
∂x
⏐⏐⏐⏐
x0
+ o
(
δ3)
.
This links the local mean of a density with the score associated with that density. Combining
this theorem with Theorem 2, we obtain that the optimal reconstruction function r∗(·) also
estimates the local mean:
mδ(x) −x= δ2
σ2(d+ 2) (r∗(x) −x) + A(δ) + δ2B(σ2) (28)
26
for error terms A(δ),B(σ2) such that
A(δ) ∈o(δ3) as δ→0,
B(σ2) ∈o(1) as σ2 →0.
This means that we can loosely estimate the direction to the local mean by the direction of
the reconstruction:
mδ(x) −x ∝ r∗(x) −x. (29)
6.5 Asymptotic formulas for localised moments
Proposition 4 Let p be of class C2 and let x0 ∈Rd. Then we have that
Zδ(x0) = δd πd/2
Γ (1 +d/2)
[
p(x0) + δ2 Tr(H(x0))
2(d+ 2) + o(δ3)
]
where H(x0) = ∂2p(x)
∂x2
⏐⏐⏐
x=x0
. Moreover, we have that
1
Zδ(x0) = δ−dΓ (1 +d/2)
πd/2
[ 1
p(x0) −δ2 1
p(x0)2
Tr(H(x0))
2(d+ 2) + o(δ3)
]
.
Proof
Zδ(x0) =
∫
Bδ(x0)
[
p(x0) + ∂p(x)
∂x
⏐⏐⏐⏐
x0
(x−x0) + 1
2!(x−x0)TH(x0)(x−x0)
+ 1
3!D(3)p(x0)(x−x0) + o(δ3)
]
dx
= p(x0)
∫
Bδ(x0)
dx + 0 + 1
2
∫
Bδ(x0)
(x−x0)TH(x0)(x−x0)dx + 0 + o(δd+3)
= p(x0)δd πd/2
Γ (1 +d/2) + δd+2 πd/2
4Γ (2 +d/2)Tr (H(x0)) + o(δd+3)
= δd πd/2
Γ (1 +d/2)
[
p(x0) + δ2 Tr(H(x0))
2(d+ 2) + o(δ3)
]
We use Proposition 10 to get that trace come up from the integral involving H(x0). The
expression for 1/Zδ(x0) comes from the fact that, for any a,b> 0 we have that
1
a+ bδ2 + o(δ3) = a−1
1 + b
aδ2 + o(δ3) = 1
a
(
1 −( b
aδ2 + o(δ3)) + o(δ4)
)
= 1
a − b
a2 δ2 + o(δ3) as δ→0.
by using the classic result from geometric series where 1
1+r = 1 −r+ r2 −... for |r|<1.
Now we just apply this to
1
Zδ(x0) = δ−dΓ (1 +d/2)
πd/2
1[
p(x0) + δ2 Tr(H(x0))
2(d+2) + o(δ3)
]
27
and get the expected result.
Theorem 5 Let p be of class C3 and represent a probability density function. Let x0 ∈Rd
with p(x0) ̸= 0. Then we have that
mδ(x0) = x0 + δ2 1
d+ 2
∂log p(x)
∂x
⏐⏐⏐⏐
x0
+ o
(
δ3)
.
Proof
The leading term in the expression for mδ(x0) is obtained by transforming the x in the
integral into a x−x0 to make the integral easier to integrate.
mδ(x0) = 1
Zδ(x0)
∫
Bδ(x0)
xp(x)dx= x0 + 1
zδ(x0)
∫
Bδ(x0)
(x−x0)p(x)dx.
Now using the Taylor expansion around x0
mδ(x0) = x0 + 1
Zδ(x0)
∫
Bδ(x0)
(x−x0)
[
p(x0) + ∂p(x)
∂x
⏐⏐⏐⏐
x0
(x−x0)
+1
2(x−x0)T ∂2p(x)
∂x2
⏐⏐⏐⏐
x0
(x−x0) + o(∥x−x0∥2)
]
dx.
Remember that
∫
Bδ(x0) f(x)dx= 0 whenever we have a function f is anti-symmetrical (or
“odd”) relative to the point x0 (i.e. f(x−x0) = f(−x−x0)). This applies to the terms
(x−x0)p(x0) and (x−x0)(x−x0) ∂2p(x)
∂x2
⏐⏐⏐
x=x0
(x−x0)T. Hence we use Proposition 9 to get
mδ(x0) = x0 + 1
Zδ(x0)
∫
Bδ(x0)
[
(x−x0)T ∂p(x)
∂x
⏐⏐⏐⏐
x0
(x−x0) + o(∥x−x0∥3)
]
dx
= x0 + 1
Zδ(x0)
(
δd+2 π
d
2
2Γ
(
2 + d
2
)
)
∂p(x)
∂x
⏐⏐⏐⏐
x0
+ o(δ3).
Now, looking at the coeﬃcient in front of ∂p(x)
∂x
⏐⏐⏐
x0
in the ﬁrst term, we can use Proposition
4 to rewrite it as
1
Zδ(x0)
(
δd+2 π
d
2
2Γ
(
2 + d
2
)
)
= δ−dΓ (1 +d/2)
πd/2
[ 1
p(x0) −δ2 1
p(x0)2
Tr(H(x0))
2(d+ 2) + o(δ3)
]
δd+2 π
d
2
2Γ
(
2 + d
2
)
= δ2 Γ
(
1 + d
2
)
2Γ
(
2 + d
2
)
[ 1
p(x0) −δ2 1
p(x0)2
Tr(H(x0))
2(d+ 2) + o(δ3)
]
= δ2 1
p(x0)
1
d+ 2 + o(δ3).
28
There is no reason the keep the −δ4 Γ(1+ d
2 )
2Γ(2+ d
2 )
1
p(x0)2
Tr(H(x0))
2(d+2) in the above expression because
the asymptotic error from the remainder term in the main expression is o(δ3). That would
swallow our exact expression for δ4 and make it useless.
We end up with
mδ(x0) = x0 + δ2 1
d+ 2
∂log p(x)
∂x
⏐⏐⏐⏐
x0
+ o(δ3).
6.6 Integration on balls and spheres
This result comes from Multi-dimensional Integration : Scary Calculus Problems from Tim
Reluga (who got the results from How to integrate a polynomial over a sphere by Gerald B.
Folland).
Theorem 6 Let B =
{
x∈Rd
⏐⏐⏐
∑d
j=1 x2
j ≤1
}
be the ball of radius 1 around the origin. Then
∫
B
d∏
j=1
|xj|aj dx=
∏Γ
(aj+1
2
)
Γ
(
1 + d
2 + 1
2
∑aj
)
for any real numbers aj ≥0.
Corollary 7 Let B be the ball of radius 1 around the origin. Then
∫
B
d∏
j=1
xaj
j dx=



∏Γ
(aj+1
2
)
Γ(1+ d
2 + 1
2
∑aj) if all the aj are even integers
0 otherwise
for any non-negative integers aj ≥0. Note the absence of the absolute values put on the xaj
j
terms.
Corollary 8 Let Bδ(0) ⊂Rd be the ball of radius δ around the origin. Then
∫
Bδ(0)
d∏
j=1
xaj
j dx=



δd+∑aj
∏Γ
(aj+1
2
)
Γ(1+ d
2 + 1
2
∑aj) if all the ajare even integers
0 otherwise
for any non-negative integers aj ≥0. Note the absence of the absolute values on the xaj
j
terms.
Proof
We take the theorem as given and concentrate here on justifying the two corollaries.
Note how in Corollary 7 we dropped the absolute values that were in the original Theorem 6.
In situations where at least one aj is odd, we have that the function f(x) = ∏d
j=1 xaj
j becomes
odd in the sense that f(−x) = −f(x). Because of the symmetrical nature of the integration on
the unit ball, we get that the integral is 0 as a result of cancellations.
29
For Corollary 8, we can rewrite the integral by changing the domain with yj = xj/δ so that
δ−∑aj
∫
Bδ(0)
d∏
j=1
xaj
j dx=
∫
Bδ(0)
d∏
j=1
(xj/δ)aj dx=
∫
B1(0)
d∏
j=1
yajδddy.
We pull out the δd that we got from the determinant of the Jacobian when changing from
dx to dy and Corollary 8 follows.
Proposition 9 Let v∈Rd and let Bδ(0) ⊂Rd be the ball of radius δ around the origin. Then
∫
Bδ(0)
y <v,y >dy =
(
δd+2 π
d
2
2Γ
(
2 + d
2
)
)
v
where <v,y > is the usual dot product.
Proof
We have that
y <v,y > =


v1y2
1
...
vdy2
d


which is decomposable into d component-wise applications of Corollary 8. This yields the
expected result with the constant obtained from Γ
(3
2
)
= 1
2 Γ
(1
2
)
= 1
2
√π.
Proposition 10 Let H ∈Rd×d and let Bδ(x0) ⊂Rd be the ball of radius δ around x0 ∈Rd.
Then
∫
Bδ(x0)
(x−x0)TH(x−x0)dx= δd+2 πd/2
2Γ (2 +d/2) trace(H) .
Proof
First, by substituting y= (x−x0) /δ we have that this is equivalent to showing that
∫
B1(0)
yTHydy = πd/2
2Γ (2 +d/2) trace (H) .
This integral yields a real number which can be written as
∫
B1(0)
yTHydy =
∫
B1(0)
∑
i,j
yiHi,jyjdy=
∑
i,j
∫
B1(0)
yiyjHi,jdy.
Now we know from Corollary 8 that this integral is zero when i̸= j. This gives
30
∑
i,j
Hi,j
∫
B1(0)
yiyjdy=
∑
i
Hi,i
∫
B1(0)
y2
idy= trace (H) πd/2
2Γ (2 +d/2).
31


Peer review:
[Review]:
This paper shows that we can relate the solution of specific autoencoder to the data generating distribution. Specifically solving for general reconstruction function with regularizer that is the L2 penalty of reconstruction contraction relates the reconstruction function derivative of the data probability log likelihood. This is in the limit of small regularization. The paper also shows that in the limit of small penalty this autoencoder is equivalent to denoising autoencoder with small noise.

Section 3.2.3: You get similar attractive behavior using almost any autoencoder with limited capacity. The point of your work is that with the specific form of regularization - square norm of contraction of r - the r(x)-x relates to derivative of log probability (proof seem to require it - it would be interesting to know what can be said about other regularizers). It would be good to compare these plots with other regularizers and show that getting log(p) for contractive one is somehow advantageous. Otherwise this section doesn't support this paper in any way.

As authors point out, it would be good to know something not in the limit of penalty going to zero. At least have some numerical experiments, for example in 1d or 2d. 

Figure 4. - 'Top plots are for one model and bottom plots for another' - what are the two models? It would be good to specify this in the figure, e.g. denosing autoencoders with different initial conditions and parameter settings.

Section 3.2.5 is important and should be written a little more clearly. 

I would suggest deriving (13) in the appendix directly from (11) without having the reader recall or read about Euler-Lagrange equations, and it might actually turn out to be simpler. Differentiating the first term with r(x) gives r(x)-x. For the second term one moves the derivative to the other size using integration by parts (and droping the boundary term) and then just applying it to the product p(x)dr/dx resulting in (13).

Minor - twice you say in the appending that the proof is in the appendinx (e.g. after statement of theorem 1)

The second last sentence in the abstract is uncomfortable to read.

This is probably not important, but can we assume that r given by (11) actually has a taylor expansion in lambda? (probably, but in the spirit of prooving things).

You don't actually derive formulas the second moments in the appendix like you do for the first moment, you mean they can similarly be derived?