Improved Autoregressive Modeling
with Distribution Smoothing

An autoregressive model (Larochelle & Murray, 2011; Salimans et al., 2017) decomposes a joint distribution $p_{\theta}(\mathbf{x})$ into the product of univariate conditionals:

where $x_{i}$ stands for the $i$-th component of $\mathbf{x}$, and $\mathbf{x}_{<i}$ refers to the components with indices smaller than $i$. In general, an autoregressive model parameterizes each conditional $p_{\theta}(x_{i}|\mathbf{x}_{<i})$ using a pre-specified density function (e.g. mixture of logistics). This bounds the capacity of the model by limiting the number of modes for each conditional.

Although autoregressive models have achieved top likelihoods amongst all types of density based models, their sample quality is still lacking compared to energy-based models (Du & Mordatch, 2019) and score-based models (Song & Ermon, 2019). We believe this can be caused by the following two reasons.

Several existing methods (Roweis & Saul, 2000; Tenenbaum et al., 2000) rely on the manifold hypothesis, i.e. that real-world high-dimensional data tends to lie on a low-dimensional manifold (Narayanan & Mitter, 2010). If the manifold hypothesis is true, then the density of the data distribution is not well defined in the ambient space; if the manifold hypothesis holds only approximately and the data lies in the vicinity of a manifold, then only points that are very close to the manifold would have high density, while all other points would have close to zero density. Thus we may expect the data density around the manifold to have large first-order derivatives, i.e. the density function has a high Lipschitz constant (if not infinity).

To see this, let us consider a 2-d example where the data distribution is a thin ring distribution (almost a unit circle) formed by rotating the 1-d Gaussian distribution $\mathcal{N}(1,0.01^{2})$ around the origin. The density function of the ring has a high Lipschitz constant near the “boundary”. Let us focus on a data point travelling along the diagonal as shown in the leftmost panel in figure 2. We plot the first-order directional derivatives of the density for the point as it approaches the boundary from the inside, then lands on the ring, and finally moves outside the ring (see figure 2). As we can see, when the point is far from the boundary, the derivative has a small magnitude. When the point moves closer to the boundary, the magnitude increases and changes significantly near the boundary even with small displacements in the trajectory. However, once the point has landed on the ring, the magnitude starts to decrease. As it gradually moves off the ring, the magnitude first increases and then decreases just like when the point approached the boundary from the inside. It has been observed that certain likelihood models, such as normalizing flows, exhibit pathological behaviors on data distributions whose densities have high Lipschitz constants (Cornish et al., 2019). Since each conditional of an autoregressive model is a 1-d normalizing flow given a fixed context, a high Lipschitz constant on data density could also hinder learning of autoregressive models.

Autoregressive models can also be susceptible to compounding errors from the conditional distributions (Lamb et al., 2016) during sampling time. We notice that an autoregressive model $p_{\theta}(\mathbf{x})$ learns the joint density $p_{\text{data}}(\mathbf{x})$ by matching each of the conditional $p_{\theta}(x_{i}|\mathbf{x}_{<i})$ with $p_{\text{data}}(x_{i}|\mathbf{x}_{<i})$. In practice, we typically have access to a limited amount of training data, which makes it hard for an autoregressive model to capture all the conditional distributions correctly due to the curse of dimensionality. During sampling, since a prediction is made based on the previously generated context, once a mistake is made at a previous step, the model is likely to make more mistakes in the later steps, eventually generating a sample $\hat{\mathbf{x}}$ that is far from being an actual image, but is mistakenly assigned a high-likelihood by the model.

The generated image $\hat{\mathbf{x}}$, being unrealistic but assigned a high-likelihood, resembles an adversarial example, i.e., an input that causes the model to make mistakes. Recent works (Cohen et al., 2019) in adversarial defense have shown that random noise can be used to improve the model’s robustness to adversarial perturbations — a process during which adversarial examples that are close to actual data are generated to fool the model. We hypothesize that such approach can also be applied to improve an autoregressive modeling process by making the model less vulnerable to compounding errors occurred during density estimation. Inspired by the success of randomized smoothing in adversarial defense (Cohen et al., 2019), we propose to apply smoothing to autoregressive modeling to address the problems mentioned above.

Unlike Cohen et al. (2019) where randomized smoothing is applied to a model, we apply smoothing directly to the data distribution $p_{\text{data}}(\mathbf{x})$. To do this, we introduce a smoothing distribution $q(\tilde{\mathbf{x}}|\mathbf{x})$ — a distribution that is symmetric and stationary (e.g. a Gaussian or Laplacian kernel) — and convolve it with $p_{\text{data}}(\mathbf{x})$ to obtain a new distribution $q(\tilde{\mathbf{x}})\triangleq\int q(\tilde{\mathbf{x}}|\mathbf{x})p_{\text{data}}(\mathbf{x})d\mathbf{x}$. When $q(\tilde{\mathbf{x}}|\mathbf{x})$ is a normal distribution, this convolution process is equivalent to perturbing the data distribution with Gaussian noise, which, intuitively, will make the data distribution smoother. In the following, we formally prove that convolving a 1-d distribution $p_{\text{data}}(x)$ with a suitable noise can indeed “smooth” $p_{\text{data}}(x)$.

Theorem 1 shows that convolving a 1-d data distribution $p_{\text{data}}(x)$ with a suitable noise distribution $q(\tilde{x}|x)$ (e.g. $\mathcal{N}(\tilde{x}|x,\sigma^{2})$) can reduce the Lipschitzness (i.e. increase the smoothness) of $p_{\text{data}}(x)$. We provide the proof of Theorem 1 in Appendix A.

Given $p_{\text{data}}(\mathbf{x})$ with a high Lipschitz constant, we empirically verify that density estimation becomes an easier task on the smoothed distribution $q(\tilde{\mathbf{x}})$ than directly on $p_{\text{data}}(\mathbf{x})$. To see this, we visualize a 1-d example in figure 3(a), where we want to model a ten-mode data distribution with a mixture of logistics model. If our model has three logistic components, there is almost no way for the model, which only has three modes, to perfectly fit this data distribution, which has ten separate modes with sharp transitions. The model, after training (see figure 3(a)), mistakenly assigns a much higher density to the low density regions between nearby modes. If we convolve the data distribution with $q(\tilde{x}|x)=\mathcal{N}(\tilde{x}|x,0.5^{2})$, the new distribution becomes smoother (see figure 3(b)) and can be captured reasonably well by the same mixture of logistics model with only three modes (see figure 3(b)). Comparing the same model’s performance on the two density estimation tasks, we can see that the model is doing a better job at modeling the smoothed version of the data distribution than the original data distribution, which has a high Lipschitz constant.

Proposition 1 shows that our smoothing process provides a regularization effect on the original objective $\mathbb{E}_{p_{\text{data}}(\mathbf{x})}[\log p_{\theta}(\mathbf{x})]$ when no noise is added, where the regularization aims to maximize $\frac{\eta}{2}\sum_{i}\frac{\partial^{2}\log p_{\theta}}{\partial x_{i}^{2}}$. Since samples from $p_{\rm{data}}$ should be close to a local maximum of the model, this encourages the second order gradients computed at a data point $\mathbf{x}$ to become closer to zero (if it were positive then $\mathbf{x}$ will not be a local maximum), creating a smoothing effect. This extra term is also the trace of the score function (up to a multiplicative constant) that can be found in the score matching objective (Hyvärinen, 2005), which is closely related to many denoising methods (Vincent, 2011; Hyvärinen, 2008). This regularization effect can, intuitively, increase the generalization capability of the model. In fact, it has been demonstrated empirically that training with noise can lead to improvements in network generalization (Sietsma & Dow, 1991; Bishop, 1995). Our argument is also similar to that used in (Bishop, 1995) except that we consider a more general generative modeling case as opposed to supervised learning with squared error. We provide the formal statement and proof of Proposition 1 in Appendix A.

Motivated by the previous 1-d example, instead of directly modeling $p_{\text{data}}(\mathbf{x})$, which can have a high Lipschitz constant, we propose to first train an autoregressive model on the smoothed version of the data distribution $q(\tilde{\mathbf{x}})$. Although the smoothing process makes the distribution easier to learn, it also introduces bias. Thus, we need an extra step to debias the learned distribution by reverting the smoothing process.

If our goal is to generate approximate samples for $p_{\rm{data}}(\mathbf{x})$, when $q(\tilde{\mathbf{x}}|\mathbf{x})=\mathcal{N}(\tilde{\mathbf{x}}|\mathbf{x},\sigma^{2}I)$ and $\sigma$ is small, we can use the gradient of $p_{\theta}(\tilde{\mathbf{x}})$ for denoising (Alain & Bengio, 2014). More specifically, given smoothed samples $\tilde{\mathbf{x}}$ from $p_{\theta}(\tilde{\mathbf{x}})$, we can “denoise” samples via:

which only requires the knowledge of $p_{\theta}(\tilde{\mathbf{x}})$ and the ability to sample from it. However, this approach does not provide a likelihood estimate and Eq. (2) only works when $q(\tilde{\mathbf{x}}|\mathbf{x})$ is Gaussian (though alternative denoising updates for other smoothing processes could be derived under the Empirical Bayes framework (Raphan & Simoncelli, 2011)). Although Eq. (2) could provide reasonable denoising results when the smoothing distribution has a small variance, $\bar{\mathbf{x}}$ obtained in this way is only a point estimation of $\bar{\mathbf{x}}=\mathbb{E}[\mathbf{x}|\tilde{\mathbf{x}}]$ and does not capture the uncertainty of $p(\mathbf{x}|\tilde{\mathbf{x}})$.

To invert more general smoothing distributions (beyond Gaussians) and to obtain likelihood estimations, we introduce a second autoregressive model $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$. The parameterized joint density $p_{\theta}(\mathbf{x},\tilde{\mathbf{x}})$ can then be computed as $p_{\theta}(\mathbf{x},\tilde{\mathbf{x}})=p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})p_{\theta}(\tilde{\mathbf{x}})$. To obtain our approximation of $p_{\text{data}}(\mathbf{x})$, we need to integrate over $\tilde{\mathbf{x}}$ on the joint distribution $p_{\theta}(\mathbf{x},\tilde{\mathbf{x}})$ to obtain $p_{\theta}(\mathbf{x})=\int p_{\theta}(\mathbf{x},\tilde{\mathbf{x}})d\tilde{\mathbf{x}}$, which is in general intractable. However, we can easily obtain an evidence lower bound (ELBO):

Note that when $q(\tilde{\mathbf{x}}|\mathbf{x})$ is fixed, the entropy term $\mathbb{E}_{q(\tilde{\mathbf{x}}|\mathbf{x})}[\log q(\tilde{\mathbf{x}}|\mathbf{x})]$ is a constant with respect to the optimization parameters. Maximizing ELBO on $p_{\text{data}}(\mathbf{x})$ is then equivalent to maximizing:

From equation 4, we can see that optimizing the two models $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ separately via maximum likelihood estimation is equivalent to optimizing $J(\theta)$.

Thus, the key here is to appropriately choose a smoothing level so that both $q(\tilde{\mathbf{x}})$ and $p(\mathbf{x}|\tilde{\mathbf{x}})$ can be approximated relatively well by existing autoregressive models. In general, the optimal variance might be hard to find. Although one can train $q(\tilde{\mathbf{x}}|\mathbf{x})$ by jointly optimizing ELBO, in practice, we find this approach often assigns a very large variance to $q(\tilde{\mathbf{x}}|\mathbf{x})$, which can trade-off sample quality for better likelihoods on high dimensional image datasets. We find empirically that a pre-specified $q(\tilde{\mathbf{x}}|\mathbf{x})$ chosen by heuristics (Saremi & Hyvarinen, 2019; Garreau et al., 2017) is able to generate much better samples than training $q(\tilde{\mathbf{x}}|\mathbf{x})$ via ELBO. In this paper, we will focus on the sample quality and leave the training of $q(\tilde{\mathbf{x}}|\mathbf{x})$ for future work.

To help us build insights into the selection of the smoothing distribution $q(\tilde{\mathbf{x}}|\mathbf{x})$, we first focus on a 1-d multi-modal distribution (see figure 4 leftmost panel). We use model-based methods to invert the smoothed distribution and provide analysis on “single-step denoising” in Appendix B.1. We start with the exploration of three different types of smoothing distributions – Gaussian distribution, Laplace distribution, and uniform distribution. For each type of distribution, we perform a grid search to find the optimal variance. Since our approach requires the modeling of both $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$, we stack $\tilde{\mathbf{x}}$ and $\mathbf{x}$ together, and use a MADE model (Germain et al., 2015) with a mixture of two logistic components to parameterize $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ at the same time. For the baseline model, we train a mixture of logistics model directly on $p_{\text{data}}(\mathbf{x})$. We compare the results in the middle two panels in figure 4.

We find that although the baseline with eight logistic components has the capacity to perfectly model the multi-modal data distribution, which has six modes, the baseline model still fails to do so. We believe this can be caused by optimization or initialization issues for modeling a distribution with a high Lipschitz constant. Our method, on the other hand, demonstrates more robustness by successfully modeling the different modes in the data distribution even when using only two mixture components for both $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$.

For all the three types of smoothing distributions, we observe a reverse U-shape correlation between the variance of $q(\tilde{\mathbf{x}}|\mathbf{x})$ and ELBO values — with ELBO first increasing as the variance increases and then decreasing as the variance grows beyond a certain point. The results match our discussion on the trade-off between modeling $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ in Section 3.3. We notice from the empirical results that Gaussian smoothing is able to obtain better ELBO than the other two distributions. Thus, we will use $q(\tilde{\mathbf{x}}|\mathbf{x})=\mathcal{N}(\tilde{\mathbf{x}}|\mathbf{x},\sigma^{2}I)$ for the later experiments.

In this section, we consider two challenging 2-d multi-modal synthetic datasets (see figure 5). We focus on model-based denoising methods and present discussion on “single-step denoising” in Appendix B.2. We use a MADE model with comparable number of total parameters for both the baseline and our approach. For the baseline, we train the MADE model directly on the data. For our randomized smoothing model, we choose $q(\tilde{\mathbf{x}}|\mathbf{x})=\mathcal{N}(\tilde{\mathbf{x}}|\mathbf{x},0.3^{2}I)$ to be the smoothing distribution. We observe that with this randomized smoothing approach, our model is able to generate better samples than the baseline (according to a human observer) even when using less logistic components (see figure 5). We provide more analysis on the model’s performance in Appendix B.2. We also provide the negative log-likelihoods in Tab. 1.

In this section, we focus on three common image datasets, namely MNIST, CIFAR-10 (Krizhevsky et al., 2009) and CelebA (Liu et al., 2015). We select $q(\tilde{\mathbf{x}}|\mathbf{x})=\mathcal{N}(\tilde{\mathbf{x}}|\mathbf{x},\sigma^{2}I)$ to be the smoothing distribution. We use PixelCNN++ (Salimans et al., 2017) as the model architecture for both $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$. We provide more details about settings in Appendix C.

Image generation. For image datasets, we select the $\sigma$ of $q(\tilde{\mathbf{x}}|\mathbf{x})=\mathcal{N}(\tilde{\mathbf{x}}|\mathbf{x},\sigma^{2}I)$ according to analysis in (Saremi & Hyvarinen, 2019) (see Appendix C for more details). Since $q(\tilde{\mathbf{x}}|\mathbf{x})$ is a Gaussian distribution, we can apply “single-step denoising” to reverse the smoothing process for samples drawn from $p_{\theta}(\tilde{\mathbf{x}})$. In this case, the model $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ is not required for sampling since the gradient of $p_{\theta}(\tilde{\mathbf{x}})$ can be used to denoise samples (also from $p_{\theta}(\tilde{\mathbf{x}})$) (see equation 2). We present smoothed samples from $p_{\theta}(\tilde{\mathbf{x}})$, reversed smoothing samples processed by “single-step denoising” and processed by $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ in figure 6. For comparison, we also present samples from a PixelCNN++ with parameters comparable to the sum of total parameters of $p_{\theta}(\tilde{\mathbf{x}})$ and $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$. We find that by using this randomized smoothing approach, we are able to drastically improve the sample quality of PixelCNN++ (see the rightmost panel in figure 6). We note that with only $p_{\theta}(\tilde{\mathbf{x}})$, a PixelCNN++ optimized on the smoothed data, we already obtain more realistic samples compared to the original PixelCNN++ method. However, $p_{\theta}(\mathbf{x}|\tilde{\mathbf{x}})$ is needed to compute the likelihood lower bounds. We report the sample quality evaluated by Fenchel Inception Distance (FID (Heusel et al., 2017)), Kernel Inception Distance (KID  (Bińkowski et al., 2018)), and Inception scores (Salimans et al., 2016) in Tab. 2. Although our method obtains better samples compared to the original PixelCNN++, our model has worse likelihoods as evaluated in BPDs. We believe this is because likelihood and sample quality are not always directly correlated as discussed in Theis et al. (2015). We also tried training the variance for $q(\tilde{\mathbf{x}}|\mathbf{x})$ by jointly optimizing ELBO. Although training the variance can produce better likelihoods, it does not generate samples with comparable quality as our method (i.e. choosing variance by heuristics). Thus, it is hard to conclusively determine what is the best way of choosing $q(\tilde{\mathbf{x}}|\mathbf{x})$. We provide more image samples in Appendix C.4 and nearest neighbors analysis in Appendix C.5.