What Regularized Auto-Encoders Learn from the Data Generating Distribution

Introduction

Machine learning is about capturing aspects of the unknown distribution from which the observed 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 configurations of the observed variables are plausible.

Unsupervised representation-learning algorithms attempt to characterize the data-generating distribution through the discovery of a set of features or latent variables whose variations capture 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; Rifai et al. 2011a; Gregor et al. 2011). An auto-encoder reconstructs the input through two stages, an encoder function $f$ (which outputs a learned representation $h=f(x)$ of an example $x$) and a decoder function $g$, such that $g(f(x))\approx 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 classification tasks and for generating higher-quality samples.

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 define 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 (first 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 differences, and thus perform approximate MCMC sampling. This is achieved using a Metropolis-Hastings MCMC in which the energy differences between the proposal and the current state are approximated using the denoising auto-encoder. Experiments on artificial 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.

Contractive and Denoising Auto-Encoders

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 coefficient is the magnitude of the perturbation. This was first 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:

The denoising auto-encoder, or DAE (Vincent et al. 2008), is trained to minimize the following denoising criterion:

where $N(x)$ is a stochastic corruption of $x$ and 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.

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+\epsilon$ with

then the optimal reconstruction function $r^{*}(x)$ will be given by

Moreover, if we consider how the optimal reconstruction function $r^{*}_{\sigma}(x)$ behaves asymptotically as $\sigma\rightarrow 0$, we get that

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 $\sigma$ are fixed when we study the effect of $\sigma\rightarrow 0$. Also, note that the $\sigma$ in the index of $r_{\sigma}^{*}$ is to indicate that $r_{\sigma}^{*}$ was chosen based on the value of $\sigma$. That $\sigma$ should not be mistaken for a parameter to be learned.

Equation (3) reveals that the optimal DAE reconstruction function at every point $x$ is 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 $\sigma$ 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^{\prime}$ such that $N(x^{\prime})=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 first thing to observe is that the leading term in the expansion of $r^{*}_{\sigma}(x)$ is $x$, and then the remainder goes to $0$ as $\sigma\rightarrow 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

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 $\sigma$ goes to $0$ could be derived from a family of noise distribution that approaches a point mass distribution in a relatively “nice” way.

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_{\sigma}(x)=x+o(1)$ as $\sigma\rightarrow 0$. In that case we get the following proposition.

Let $p$ be the probability density function of the data. Consider a $DAE$ using the expected quadratic loss and corruption noise $N(x)=x+\epsilon$, with $\epsilon\sim\mathcal{N}\left(0,\sigma^{2}I\right)$. If we assume that the non-parametric solutions $r_{\sigma}(x)$ satistfies

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 $\sigma^{2}$ is similar to a contractive auto-encoder with penalty coefficient $\lambda=\sigma^{2}$ but where the contraction is imposed explicitly on the whole reconstruction function $r(\cdot)=g(f(\cdot))$ rather than on $f(\cdot)$ alone In the CAE there is a also a contractive effect on $g(\cdot)$ as a side effect of the parametrization with weights tied between $f(\cdot)$ and $g(\cdot)$..

This analysis motivates the definition of the reconstruction contractive auto-encoder (RCAE), a variation of the CAE where loss function is instead the squared reconstruction loss plus contractive penalty on the reconstruction:

This is an analytic version of the denoising criterion with small noise $\sigma^{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 $\lambda=\sigma^{2}$ and because the semantics of $\sigma^{2}$ is clearer (as a squared distance in input space), we will denote $\sigma^{2}$ for the penalty term coefficient in situations involving RCAE. For example, in the statement of Theorem 2, this $\sigma^{2}$ is just a positive constant; there is no notion of additive Gaussian noise, i.e., $\sigma^{2}$ does not explicitly refer to a variance, but using the notation $\sigma^{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 following 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.

Let $r_{\sigma}^{*}(x)$ denote the optimal function that minimizes $\mathcal{L}_{\sigma}$. Then we have that

Moreover, we also have the following expression for the derivative

Both these asymptotic expansions are to be understood in a context where we consider $\left\{r_{\sigma}^{*}(x)\right\}_{{\sigma}\geq 0}$ to be a family of optimal functions minimizing $\mathcal{L}_{\sigma}$ for their corresponding value of ${\sigma}$. The asymptotic expansions are applicable point-wise in $x$, that is, with any fixed $x$ we look at the behavior as ${\sigma}\rightarrow 0$.

The proof is given in the appendix and uses the Euler-Lagrange equations from the calculus of variations.

Minimizing the Loss to Recover Local Features of p⁡(⋅)p(\cdot)

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 ${\sigma}\rightarrow 0$) for the optimal reconstruction function for the DAE and RCAE that allows us to recover the score $\frac{\partial\log p(x)}{\partial 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.

In an experimental setting, the expected loss (7) is replaced by the empirical loss

based on a sample $\left\{x^{(n)}\right\}_{n=1}^{N}$ 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 ${\sigma^{2}}>0$. This is assumed to have been done through minimizing the empirical loss and the distribution $p$ was only available indirectly through the samples $\left\{x^{(n)}\right\}_{n=1}^{N}$. We do not have access to $p$ or to the samples. We have only $r(x)$ and maybe ${\sigma^{2}}$.

We will now discuss the usefulness of $r(x)$ based on different conditions such as the model capacity and the value of ${\sigma^{2}}$.

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 infinite amount of training data, and a perfect minimization. We will see what can be done to recover information about $p$ in 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_{\sigma}(x)$ when we want to emphasize the fact that the value of $r(x)$ came from minimizing the loss with a certain fixed ${\sigma}$.

Suppose that $r_{\sigma}(x)$ was trained with an infinite sample drawn from $p$. Suppose also that it had infinite (or sufficient) 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 differentiable. Suppose that we know the value of ${\sigma}$ and that we are working in a computing environment of arbitrary precision (i.e. no rounding errors).

In the setup described, we do not get to pick values of ${\sigma}$ so as to take the limit ${\sigma}\rightarrow 0$. However, it is assumed that ${\sigma}$ is already sufficiently small that the above quantity is close to $\frac{\partial\log p(x)}{\partial x}$ for all intents and purposes.

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 $\hat{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 defined 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 $\frac{\partial}{\partial x}\log p(x)$ as $\sigma\rightarrow 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 $\hat{r}(x)$ is fitted by dividing the interval $\left[-1.5,1.5\right]$ into $M=1000$ partition points $x_{1},\ldots,x_{M}$ evenly separated by a distance $\Delta$. The discretized version of the RCAE loss function is

Every value $\hat{r}(x_{i})$ for $i=1,\ldots,M$ is treated as a free parameter. Setting to $0$ the derivative with respect to the $\hat{r}(x_{i})$ yields a system of linear equations in $M$ unknowns that we can solve exactly. From that RCAE solution $\hat{r}$ we get an approximation of the score of $p$ at each point $x_{i}$. 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$.

In Figure 4 we compare the approximations to the actual score of $p$ for decreasingly small values of $\sigma\in\{1.00,0.31,0.16,0.06\}$.

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)-x$ as a vector field, 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:

A denoising auto-encoder was trained with Gaussian corruption noise $\sigma=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 fixed 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 influence the final outcome. Moreover, the fact that we are using a finite 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 field). We see that that the vector field points towards the nearest high-density point on the data manifold. The vector field 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 play the role of sinks for the vector field. Other places where reconstruction error may be low, but where the implicit density is not high, are sources of the vector field. 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 different 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 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 firmed up by analyzing Figure 2: the score-like behavior of $r(x)-x$ arises 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 classification 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 field and adding noise along the way.

5 Missing σ2{\sigma^{2}}

When we are in the same setting as in section 3.2 but the value of ${\sigma^{2}}$ is unknown, we can modify (9) a bit and avoid dividing by ${\sigma^{2}}$. That is, for a trained reconstruction function $r(x)$ given to us we just take the quantity $r(x)-x$ and it should be approximatively the score up to a multiplicative constant.

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.

We still have to assume that ${\sigma^{2}}$ is small. Otherwise, if the unknown ${\sigma^{2}}$ is too large we might get a poor estimation of the score.

6 Limited Parameterization

We should also be concerned about the fact that $r(x)-x$ is trying to approximate $-\frac{\partial E(x)}{\partial x}$ as ${\sigma}\rightarrow 0$ but we have not made any assumptions about the space of functions that $r$ can represent when we are dealing with a specific 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 field (i.e. the gradient of a potential function). Even if we start from a density $p(x)\propto\exp(-E(x))$ and we have that $r(x)-x$ is very close to $-\frac{\partial}{\partial x}E(x)$ in terms of some given norm, there is not guarantee that there exists an associated function $E_{0}(x)$ for which $r(x)-x\propto-\frac{\partial}{\partial x}E_{0}(x)$ and $E_{0}(x)\approx E(x)$.

In fact, in many cases we can trivially show the non-existence of such a $E_{0}(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 $x$ terms from the computations because we can easily find its antiderivative by looking at $x=\frac{\partial}{\partial x}\left\|x\right\|^{2}_{2}$.

Conceptually, another way to see this is to argue that if such a function $E_{0}(x)$ existed, its second-order mixed derivatives should be equal. That is, we should have that

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.

then the above expectation is equivalent to

which is the denoising criterion. This says that when the reconstruction function $r$ is parametrized so as to correspond to the score $\psi$ of a model density (as per eq. 11, and where $\psi$ is a derivative of some log-density), the denoising criterion on $r$ with 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 insufficiently smoothed version of it) could yield undesirable results when the training set is finite. Since score matching has been shown to be a consistent induction principle (Hyvärinen 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 $\sigma^{2}$. By making $\sigma^{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 coefficient and the DAE’s noise variance $\sigma^{2}$, and (2) that minimizing the equivalent analytic criterion (based on a contraction penalty) estimates the score when ${\sigma^{2}}$ is small. The difference 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.

8 Estimating the Hessian

Since we have $\frac{r(x)-x}{{\sigma^{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:

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 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 first 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 first 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 $\sigma^{2}$ is near 0 and $r(x)$ is near $x$, so that $\frac{\partial r(x)}{\partial x}$ is close to the identity. In particular, in the ridge (manifold) directions, we should expect $\frac{\partial r(x)}{\partial 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 $\frac{\partial r(x)}{\partial x}$ that are near 1, making the corresponding eigenvalues of $\frac{\partial^{2}\log p(x)}{\partial x^{2}}$ near 0. On the other hand, in the directions orthogonal to the manifold, $\frac{\partial r(x)}{\partial x}$ should be smaller than 1, making the corresponding eigenvalues of $\frac{\partial^{2}\log p(x)}{\partial x^{2}}$ negative.

Besides first 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.

Sampling with Metropolis-Hastings

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 difference $E(x^{*})-E(x)$ between two states $x$ and $x^{*}$. This can be done by using a first-order Taylor approximation

The simplest way to discretize this path integral is to pick points $\left\{x_{i}\right\}_{i=1}^{n}$ spread at even distances on a straight line from $x_{1}=x$ to $x_{n}=x^{*}$. We approximate (12) by

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 differences which are needed in the Metropolis-Hastings accept/reject decision. Using a symmetric proposal $q(x^{*}|x)$, the acceptance ratio is

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 sufficiently good estimator of $\frac{\partial E}{\partial 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 $\sigma_{\textrm{train}}=0.1$ so that the noise would have an appreciable effect while still being relatively small. As required by Theorem 1, we have used isotropic Gaussian noise of variance $\sigma_{\textrm{train}}^{2}$.

The Metropolis-Hastings proposal $q(x^{*}|x)=\mathcal{N}(0,\sigma_{\textrm{MH}}^{2}I)$ has a noise parameter $\sigma_{\textrm{MH}}$ that needs to be set. In the situation shown in Figure 6, we used $\sigma_{\textrm{MH}}=0.1$. After some hyperparameter tweaking and exploring various scales for $\sigma_{\textrm{train}},\sigma_{\textrm{MH}}$, we found that setting both to be $0.1$ worked well.

When $\sigma_{\textrm{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 $\sigma_{\textrm{train}}\rightarrow 0$ so we would expect that the estimator becomes worse when $\sigma_{\textrm{train}}$ is taking on larger values. In this particular case with $\sigma_{\textrm{train}}=0.1$, it seems that we are instead modeling something like the original density to which isotropic Gaussian noise of variance $\sigma_{\textrm{train}}^{2}$ has been added.

In the other extreme, when $\sigma_{\textrm{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 $\sigma_{\textrm{train}}=\sigma_{\textrm{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 $\sigma_{\textrm{MH}}$ lead to higher acceptance ratios but worse mixing properties.

3 Spurious Maxima

There are two very real concerns with the sampling method discussed in section 4.2. The first 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 field $r(x)-x$ for a DAE that failed to properly learn the desired behavior in regions away from the spiral-shaped density.

Conclusion

Whereas auto-encoders have long been suspected of capturing information about the data generating density, this work has clarified what some of them are actually doing, showing that they can actually implicitly recover the data generating density altogether. We have shown that 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 first 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 findings. 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 matching (Hyvärinen 2005; Vincent 2011). Toy experiments have confirmed that a good estimator of the density can be obtained when this criterion is non-parametrically minimized. The experiments have also confirmed that an MCMC could be setup that approximately samples from the estimated model, by estimating energy differences to first order (which only requires the score) to perform approximate Metropolis-Hastings MCMC.

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 fit the qualitative pattern illustrated in Figure 2 where a score-like vector field 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).

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

Appendix

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+\epsilon$ with

then the optimal reconstruction function $r^{*}(x)$ will be given by

Moreover, if we consider how the optimal reconstruction function $r^{*}_{\sigma}(x)$ behaves asymptotically as $\sigma\rightarrow 0$, we get that

The first part of this proof is to get to equation (14) without assuming that $\sigma\rightarrow 0$.

Now, for the second part of this proof, we study the behavior of the solution $r_{\sigma}^{*}(x)$ as $\sigma$ 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 $\epsilon\sim\mathcal{N}\left(0,\sigma^{2}I\right)$, we can also simplify equation (15) by converting all the $-\epsilon$ into $+\epsilon$.

We now look at the Taylor expansion of $p(x+\epsilon)$ around $x$. We perform this expansion inside of the expectation, where $\epsilon$ just represents as small quantity of scale $\sigma$.

When taking the expectation of $p(x+\epsilon)$ with respect to $\epsilon$, we get a zero for all the terms containing an odd power of $\epsilon$, for reasons of symmetry. When we take the expectation of $p(x+\epsilon)\epsilon$ instead, all the terms with an even power of $\epsilon$ in the above expectation will vanish. Thus, we get that

Provided that $p(x)\neq 0$, our quotient can be rewritten as

We can use the basic geometric series expansion to write

Note that $p(x)$ was treated as a constant when we studied the asymptotic behavior as $\sigma\rightarrow 0$. In the Taylor expansion around some given $x$, we want to stuff all the higher-order derivatives of $p(x)$ into the asymptotic remainder term. It is quite possible that the size of the $\sigma$ required to do so would depend on the particular $x$, that there would not be a uniform $\sigma>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_{\sigma}^{*}(x)$. For practical applications, we do not need more than that.

The assumption that $p(x)\neq 0$ could also be viewed as problematic, but if we go back to the definition of the DAE loss, then any point where $p(x)=0$ would never produce a training example and would not even contribute in the definition of the expectation of DAE loss. The correct value to assign to $r$ at such a point is not well-defined.

2 Relationship between Contractive Penalty and Denoising Criterion

Let $p$ be the probability density function of the data. Consider a $DAE$ using the expected quadratic loss and corruption noise $N(x)=x+\epsilon$, with $\epsilon\sim\mathcal{N}\left(0,\sigma^{2}I\right)$. If we assume that the non-parametric solutions $r_{\sigma}(x)$ satistfies

where the expectation is taken with respect to $X$, whose distribution is given by $p$.

We can drop the $\sigma$ index from $r_{\sigma}$ 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

The DAE loss involves taking the expectation with respect to $X$ and with respect to the noise $\epsilon$. We substitute the Taylor expansion into ${\cal L}_{DAE}$, we express the norm as a dot product, and we show how the expectation with respect to $\epsilon$ cancels out certain terms.

3 Calculus of Variations

Let $r_{\sigma}^{*}(x)$ denote the optimal function that minimizes $\mathcal{L}_{\sigma}$. Then we have that

Moreover, we also have the following expression for the derivative

Both these asymptotic expansions are to be understood in a context where we consider $\left\{r_{\sigma}^{*}(x)\right\}_{{\sigma}\geq 0}$ to be a family of optimal functions minimizing $\mathcal{L}_{\sigma}$ for their corresponding value of ${\sigma}$. The asymptotic expansions are applicable point-wise in $x$, that is, with any fixed $x$ we look at the behavior as ${\sigma}\rightarrow 0$.

This proof is done in two parts. In the first part, the objective is to get to equation (22) that has to be satisfied for the optimum solution.

We treat ${\sigma}$ as given and constant for the first part of this proof.

In the second part we work out the asymptotic expansion in terms of ${\sigma}$. We again work with the implicit dependence of $r(x)$ on ${\sigma}$.

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

where $x=(x_{1},\ldots,x_{d})$ , $r(x)=(r_{1}(x),\ldots,r_{d}(x))$ and $r_{x_{i}}=\frac{\partial f}{\partial x_{i}}$.

We can rewrite the loss $\mathcal{L}(r)$ more explicitly as

to observe that the components $r_{1}(x),\ldots,r_{d}(x)$ can each be optimized separately.

In our situation, the expressions from that equation are given by

and the equality to be satisfied at the optimum becomes

As $p(x)\neq 0$ by hypothesis, we can divide all the terms by $p(x)$ and note that $\frac{\partial p(x)}{\partial x_{i}}/p(x)=\frac{\partial\log p(x)}{\partial x_{i}}$.

This first thing to observe is that when ${\sigma^{2}}=0$ the solution is just $r_{k}(x)=x_{k}$, 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.

This linear partial differential equation (22) can be used as a recursive relation for $r_{k}(x)$ to obtain a Taylor series in ${\sigma^{2}}$. The goal is to obtain an expression of the form

where we can solve for $h(x)$ and for which we also have that

We can substitute in the right-hand side of equation (23) the value for $r_{k}(x)$ that we get from equation (23) itself. This substitution would be pointless in any other situation where we are not trying to get a power series in terms of ${\sigma^{2}}$ around $0$.

Now we would like to get rid of that ${\sigma^{2}}\sum_{i=1}^{d}\frac{\partial^{2}r_{k}(x)}{\partial x_{i}^{2}}$ term by showing that it is a term that involves only powers of ${\sigma^{4}}$ or higher. We get this by showing what we get by differentiating the expression for $r_{k}(x)$ in line (29) twice with respect to some $l$.

Since ${\sigma^{2}}$ is a common factor in all the terms of the expression of $\frac{\partial^{2}r_{k}(x)}{\partial x_{l}^{2}}$ we get what we needed. That is,

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 first 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 definitions introduced in Bengio et al. 2012a.

where $Z_{\delta}(x_{0})$ is the normalizing constant required to make $p_{\delta}(x|x_{0})$ a valid pdf for a distribution centered on $x_{0}$. The support of $p_{\delta}(x|x_{0})$ is the ball of radius $\delta$ around $x_{0}$ denoted by $B_{\delta}(x_{0})$. We stick to the $2$-norm in terms of defining the balls $B_{\delta}(x_{0})$ used, but everything could be rewritten in terms of another $p$-norm to have slightly different formulas.

We use the following notation for what will be referred to as the first two local moments (i.e. local mean and local covariance) of the random variable described by $p_{\delta}(x|x_{0})$.

Based on these definitions, one can prove the following theorem.

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^{*}(\cdot)$ also estimates the local mean:

for error terms $A(\delta),B({\sigma^{2}})$ such that

This means that we can loosely estimate the direction to the local mean by the direction of the reconstruction:

5 Asymptotic formulas for localised moments

where $H(x_{0})=\left.\frac{\partial^{2}p(x)}{\partial x^{2}}\right|_{x=x_{0}}$. Moreover, we have that

We use Proposition 10 to get that trace come up from the integral involving $H(x_{0})$. The expression for $1/Z_{\delta}(x_{0})$ comes from the fact that, for any $a,b>0$ we have that

by using the classic result from geometric series where $\frac{1}{1+r}=1-r+r^{2}-\ldots$ for $|r|<1$.

The leading term in the expression for $m_{\delta}(x_{0})$ is obtained by transforming the $x$ in the integral into a $x-x_{0}$ to make the integral easier to integrate.

Now using the Taylor expansion around $x_{0}$

Remember that $\int_{B_{\delta}(x_{0})}f(x)dx=0$ whenever we have a function $f$ is anti-symmetrical (or “odd”) relative to the point $x_{0}$ (i.e. $f(x-x_{0})=f(-x-x_{0})$). This applies to the terms $(x-x_{0})p(x_{0})$ and $(x-x_{0})(x-x_{0})\left.\frac{\partial^{2}p(x)}{\partial x^{2}}\right|_{x=x_{0}}(x-x_{0})^{T}$. Hence we use Proposition 9 to get

Now, looking at the coefficient in front of $\left.\frac{\partial p(x)}{\partial x}\right|_{x_{0}}$ in the first term, we can use Proposition 4 to rewrite it as

There is no reason the keep the $-\delta^{4}\frac{\Gamma\left(1+\frac{d}{2}\right)}{2\Gamma\left(2+\frac{d}{2}\right)}\frac{1}{p(x_{0})^{2}}\frac{\textrm{Tr}(H(x_{0}))}{2(d+2)}$ in the above expression because the asymptotic error from the remainder term in the main expression is $o(\delta^{3})$. That would swallow our exact expression for $\delta^{4}$ and make it useless.

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).

Let $B$ be the ball of radius $1$ around the origin. Then

for any non-negative integers $a_{j}\geq 0$. Note the absence of the absolute values put on the $x_{j}^{a_{j}}$ terms.

for any non-negative integers $a_{j}\geq 0$. Note the absence of the absolute values on the $x_{j}^{a_{j}}$ terms.

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 $a_{j}$ is odd, we have that the function $f(x)=\prod_{j=1}^{d}x_{j}^{a_{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.

For Corollary 8, we can rewrite the integral by changing the domain with $y_{j}=x_{j}/\delta$ so that

We pull out the $\delta^{d}$ that we got from the determinant of the Jacobian when changing from $dx$ to $dy$ and Corollary 8 follows.

which is decomposable into $d$ component-wise applications of Corollary 8. This yields the expected result with the constant obtained from $\Gamma\left(\frac{3}{2}\right)=\frac{1}{2}\Gamma\left(\frac{1}{2}\right)=\frac{1}{2}\sqrt{\pi}.$

First, by substituting $y=\left(x-x_{0}\right)/\delta$ we have that this is equivalent to showing that

This integral yields a real number which can be written as

Now we know from Corollary 8 that this integral is zero when $i\neq j$. This gives