Learning deep kernels for exponential family densities

Introduction

Density estimation is a foundational problem in statistics and machine learning (Devroye & Györfi, 1985; Wasserman, 2006), lying at the core of both supervised and unsupervised machine learning problems. Classical techniques such as kernel density estimation, however, struggle to exploit the structure inherent to complex problems, and thus can require unreasonably large sample sizes for adequate fits. For instance, assuming only twice-differentiable densities, the $L_{2}$ risk of density estimation with $N$ samples in $D$ dimensions scales as $\mathcal{O}(N^{-4/(4+D)})$ (Wasserman, 2006, Section 6.5).

One promising approach for incorporating this necessary structure is the kernel exponential family (Canu & Smola, 2006; Fukumizu, 2009; Sriperumbudur et al., 2017). This model allows for any log-density which is suitably smooth under a given kernel, i.e. any function in the corresponding reproducing kernel Hilbert space. Choosing a finite-dimensional kernel recovers any classical exponential family, but when the kernel is sufficiently powerful the class becomes very rich: dense in the family of continuous probability densities on compact domains in KL, TV, Hellinger, and $L^{r}$ distances (Sriperumbudur et al., 2017, Corollary 2). The normalization constant is not available in closed form, making fitting by maximum likelihood difficult, but the alternative technique of score matching (Hyvärinen, 2005) allows for practical usage with theoretical convergence guarantees (Sriperumbudur et al., 2017).

The choice of kernel directly corresponds to a smoothness assumption on the model, allowing one to design a kernel corresponding to prior knowledge about the target density. Yet explicitly deciding upon a kernel to incorporate that knowledge can be complicated. Indeed, previous applications of the kernel exponential family model have exclusively employed simple kernels, such as the Gaussian, with a small number of parameters (e.g. the length scale) chosen by heuristics or cross-validation (e.g. Sasaki et al., 2014; Strathmann et al., 2015). These kernels are typically spatially invariant, corresponding to a uniform smoothness assumption across the domain. Although such kernels are sufficient for consistency in the infinite-sample limit, the induced models can fail in practice on finite datasets, especially if the data takes differently-scaled shapes in different parts of the space. Figure 1 (left) illustrates this problem when fitting a simple mixture of Gaussians. Here there are two “correct” bandwidths, one for the broad mode and one for the narrow mode. A translation-invariant kernel must pick a single one, e.g. an average between the two, and any choice will yield a poor fit on at least part of the density.

In this work, we propose to learn the kernel of an exponential family directly from data. We can then achieve far more than simply tuning a length scale, instead learning location-dependent kernels that adapt to the underlying shape and smoothness of the target density. We use kernels of the form

where the deep network $\bm{\phi}$ extracts features of the input and $\kappa$ is a simple kernel (e.g. a Gaussian) on those features. These types of kernels have seen success in supervised learning (Wilson et al., 2016; Jean et al., 2018) and critic functions for implicit generative models (Li et al., 2017; Bellemare et al., 2017; Bińkowski et al., 2018; Arbel et al., 2018), among other settings. We call the resulting model a deep kernel exponential family (DKEF).

We can train both kernel parameters (including all the weights of the deep network) and, unusually, even regularization parameters directly on the data, in a form of meta-learning. Normally, directly optimizing regularization parameters would always yield 0, since their beneficial effect in preventing overfitting is by definition not seen on the training set. Here, though, we can exploit the closed-form fit of the kernel exponential family to optimize a “held-out” score (Section 3). Figure 1 (right) demonstrates the success of this model on the same mixture of Gaussians; here the learned, location-dependent kernel gives a much better fit.

We compare the results of our new model to recent general-purpose deep density estimators, primarily autoregressive models (Uria et al., 2013; Germain et al., 2015; van den Oord et al., 2016) and normalizing flows (Jimenez Rezende & Mohamed, 2015; Dinh et al., 2017; Papamakarios et al., 2017). These models learn deep networks with structures designed to compute normalized densities, and are fit via maximum likelihood. We explore the strengths and limitations of both deep likelihood models and deep kernel exponential families on a variety of datasets, including artificial data designed to illustrate scenarios where certain surprising problems arise, as well as benchmark datasets used previously in the literature. The models fit by maximum likelihood typically give somewhat higher likelihoods, whereas the deep kernel exponential family generally better fits the shape of the distribution.

Background

Under mild regularity conditions, this is equal to

up to an additive constant depending only on $p_{0}$, which can be ignored during training. Here $\partial_{d}^{n}f({\bm{x}})$ denotes $\frac{\partial^{n}}{(\partial{y}_{d})^{n}}f({\bm{y}})\rvert_{{\bm{y}}={\bm{x}}}$. We can estimate (3) with $\hat{J}(p_{\bm{\theta}},\mathcal{D})$:

Kernel exponential families

Using a simple finite-dimensional $\mathcal{H}$, we can recover any standard exponential family, e.g. normal, gamma, or Poisson; if $\mathcal{H}$ is sufficiently rich, the family can approximate any continuous distribution with tails like $q_{0}$ arbitrarily well (Sriperumbudur et al., 2017, Example 1 and Corollary 2).

Sutherland et al. (2018) used a small $\lambda_{\alpha}$ for numerical stability but primarily regularized with $\lambda_{H}\lVert f_{{\bm{\alpha}},{\bm{z}}}^{k}\rVert_{\mathcal{H}}^{2}$. As we change $k$, however, $\lVert f\rVert_{\mathcal{H}}$ changes meaning, and we found empirically that this regularizer tends to harm the fit. The $\lambda_{C}$ term was recommended by Kingma & LeCun (2010), encouraging the learned log-density to be smooth without much extra computation; it provides some empirical benefit in our context. Given $k$, ${\bm{z}}$, and ${\bm{\lambda}}$, 3 (Appendix B) shows we can find the optimal ${\bm{\alpha}}$ by solving an $M\times M$ linear system in $\mathcal{O}(M^{2}ND+M^{3})$ time: the ${\bm{\alpha}}$ which minimizes $\hat{J}(f_{{\bm{\alpha}},{\bm{z}}}^{k},{\bm{\lambda}},\mathcal{D})$ is

Fitting Deep Kernels

All previous applications of score matching in the kernel exponential family of which we are aware (e.g. Sriperumbudur et al., 2017; Strathmann et al., 2015; Sun et al., 2015; Sutherland et al., 2018) have used kernels of the form $k({\bm{x}},{\bm{y}})=\exp\left(-\frac{1}{2\sigma^{2}}\lVert{\bm{x}}-{\bm{y}}\rVert^{2}\right)+r\left({\bm{x}}^{\mathsf{T}}{\bm{y}}+c\right)^{2}$, with kernel parameters and regularization weights either fixed a priori or selected via cross-validation. This simple form allows the various kernel derivatives required in (6) to be easily computed by hand, and the small number of parameters makes grid search adequate for model selection. But, as discussed in Section 1, these simple kernels are insufficient for complex datasets. Thus we wish to use a richer class of kernels $\{k_{\bm{w}}\}$, with a large number of parameters ${\bm{w}}$ – in particular, kernels defined by a neural network. This prohibits model selection via simple grid search.

One could attempt to directly minimize $\hat{J}(f_{{\bm{\alpha}},{\bm{z}}}^{k_{\bm{w}}},{\bm{\lambda}},\mathcal{D})$ jointly in the kernel parameters ${\bm{w}}$, the model parameters ${\bm{\alpha}}$, and perhaps the inducing points ${\bm{z}}$. Consider, however, the case where we simply use a Gaussian kernel and $\{{\bm{z}}_{m}\}=\mathcal{D}$. Then we can achieve arbitrarily good values of (3) by taking $\sigma\to 0$, drastically overfitting to the training set $\mathcal{D}$.

Solving for ${\bm{\alpha}}$ and computing the loss (4) require matrices of kernel second derivatives, but current deep learning-oriented automatic differentiation systems are not optimized for evaluating tensor-valued higher-order derivatives at once. We therefore implement backpropagation to compute ${\bm{G}}$, ${\bm{U}}$, and ${\bm{b}}$ of (6) as TensorFlow operations (Abadi et al., 2015) to obtain the scalar loss $\hat{J}$, and used TensorFlow’s automatic differentiation only to optimize ${\bm{w}}$, ${\bm{z}}$, ${\bm{\lambda}}$, and $q_{0}$ parameters.

Backpropagation to find these second derivatives requires explicitly computing the Hessians of intermediate layers of the network, which becomes quite expensive as the model grows; this limits the size of kernels that our model can use. A more efficient implementation based on Hessian-vector products might be possible in an automatic differentiation system with better support for matrix derivatives.

Kernel architecture

We will choose our kernel $k_{\bm{w}}({\bm{x}},{\bm{y}})$ as a mixture of $R$ Gaussian kernels with length scales $\sigma_{r}$, taking in features of the data extracted by a network $\bm{\phi}_{{\bm{w}}_{r}}(\cdot)$:

Combining $R$ components makes it easier to account for both short-range and long-range dependencies. We constrain $\rho_{r}\geq 0$ to ensure a valid kernel, and $\sum_{r=1}^{R}\rho_{r}=1$ for simplicity. The networks $\bm{\phi}_{{\bm{w}}}$ are made of $L$ fully connected layers of width $W$. For $L>1$, we found that adding a skip connection from data directly to the top layer speeds up learning. A softplus nonlinearity, $\log(1+\exp(x))$, ensures that the model is twice-differentiable so (3) is well-defined.

1 Behavior on Mixtures

If all modes are connected by regions of positive density, then the log density gradient in between components will determine their relative weight, and indeed score matching is then consistent. But when $p_{0}$ is nearly zero between two dense components, so that there are no or few samples in between, score matching will generally have insufficient evidence to weight nearly-separate components.

4 (Appendix C) studies the the kernel exponential family in this case. For two components that are completely separated according to $k$, (6) fits each as it would if given only that component, except that the effective $\lambda_{\alpha}$ is scaled: smaller components are regularized more.

Section C.1 studies a simplified case where the kernel length scale $\sigma$ is far wider than the component; then the density ratio between components, which should be $\frac{\pi}{1-\pi}$, is approximately $\exp\left(\frac{D}{2\sigma^{2}\lambda_{\alpha}}\left(\pi-\frac{1}{2}\right)\right)$. Depending on the value of $\frac{D}{2\sigma^{2}\lambda_{\alpha}}$, this ratio will often either be quite extreme, or nearly $1$. It is unclear, however, how well this result generalizes to other settings.

A heuristic workaround when disjoint components are suspected is as follows: run a clustering algorithm to identify disjoint components, separately fit a model to each cluster, then weight each model according to its sample count. When the components are well-separated, this clustering is straightforward, but it may be difficult in high-dimensional cases when samples are sparse but not fully separated.

2 Model Evaluation

As all our models are twice-differentiable, we also compare the score-matching loss (4) on held-out test data. A lower score-matching loss implies a smaller Fisher divergence between model and data distributions.

Finally, we compare test log-likelihoods, using importance sampling estimates of the normalizing constant $Z_{\bm{\theta}}$:

(Proof in Appendix D.) We can find $a$, $s$ because we propose from $q_{0}$, and thus we can effectively estimate the bound (Section D.1). This estimate of the upper bound is itself biased upwards (6), so it is likely, though not guaranteed, that the estimate overstates the amount of bias.

3 Previous Attempts at Deep Score Matching

The scalar in brackets is fit analytically, so $\log p$ is determined almost entirely by the network $\phi_{\bm{w}}$ plus $\log q_{0}({\bm{x}})$.

Saremi et al. (2018) recently also attempted parameterizing the unnormalizing log-density as a deep network, using an approximation called Parzen score matching. This approximation requires a global constant bandwidth to define the Parzen window size for the loss function, fit to the dataset before learning the model. This is likely only sensible on datasets for which simple fixed-bandwidth kernel density estimation is appropriate; on more complex datasets, the loss may be poorly motivated. It also leads to substantial oversmoothing visible in their results. As they did not provide code for their method, we do not compare to it empirically.

Experimental Results

In our experiments, we compare to several alternative methods. The first group are all fit by maximum likelihood, and broadly fall into (at least) one of two categories: autoregressive models decompose $p({\textnormal{x}}_{1},\dots,{\textnormal{x}}_{D})=\prod_{d=1}^{D}p({\textnormal{x}}_{d}|{\mathbf{x}}_{\leq d})$ and learn a parametric density model for each of these conditionals. Normalizing flows instead apply a series of invertible transformations to some simple initial density, say standard normal, and then compute the density of the overall model via the Jacobian of the transformation. We use implementationsgithub.com/gpapamak/maf of the following several models in these categories by Papamakarios et al. (2017):

MADE (Germain et al., 2015) masks the connections of an autoencoder so it is autoregressive. We use two hidden layers, and each conditional a Gaussian. MADE-MOG is the same but with each conditional a mixture of 10 Gaussians.

Real NVP (Dinh et al., 2017) is a normalizing flow; we use a general-purpose form for non-image datasets.

MAF (Papamakarios et al., 2017). A combination of a normalizing flow and MADE, where the base density is modeled by MADE, with 5 autoregressive layers. MAF-MOG instead models the base density by MADE-MOG.

For the models above, we use layers of width 30 for experiments on synthetic data, and 100 for benchmark datasets. Larger values did not improve performance.

KCEF (Arbel & Gretton, 2018). Inspired by autoregressive models, the density is modeled by a cascade of kernel conditional exponential family distributions, fit by score matching with Gaussian kernels.github.com/MichaelArbel/KCEF

DKEF. On synthetic datasets, we consider four variants of our model with one kernel component, $R=1$. KEF-G refers to the model using a Gaussian kernel with a learned bandwidth. DKEF-G-15 has the kernel (7), with $L=3$ layers of width $W=15$. DKEF-G-50 is the same with $W=50$. To investigate whether the top Gaussian kernel helps performance, we also train DKEF-L-50, whose kernel is $k_{\theta}({\bm{x}},{\bm{y}})=\bm{\phi}_{\bm{w}}({\bm{x}})\cdot\bm{\phi}_{\bm{w}}({\bm{y}})$, where $\bm{\phi}_{\bm{w}}$ has $W=50$. To compare with the architecture of Kingma & LeCun (2010), DKEF-L-50-1 has the same architecture as DKEF-L-50 except that we add an extra layer with a single neuron, and use $M=1$. In all experiments, $q_{0}({\bm{x}})=\prod_{d=1}^{D}\exp\left(-\lvert x_{d}-\mu_{d}\rvert^{\beta_{d}}/(2\sigma_{d}^{2})\right)$, with $\beta_{d}>1$. On benchmark datasets, we use DKEF-G-50 and KEF-G with three kernel components, $R=3$. Code for DKEF is at github.com/kevin-w-li/deep-kexpfam.

We first demonstrate the behavior of the models on several two-dimensional synthetic datasets Funnel, Banana, Ring, Square, Cosine, Mixture of Gaussians (MoG) and Mixture of Rings (MoR). Together, they cover a range of geometric complexities and multimodality.

We visualize the fit of various methods by showing the log density function in Figure 2. For each model fit on each distribution, we report the normalized log likelihood (left) and Fisher divergence (right). In general, the kernel score matching methods find cleaner boundaries of the distributions, and our main model KDEF-G-15 produces the lowest Fisher divergence on many of the synthetic datasets while maintaining high likelihoods.

Among the kernel exponential families, DKEF-G outperformed previous versions where ordinary Gaussian kernels are used for either joint (KEF-G) or autoregressive (KCEF) modeling. DKEF-G-50 does not substantially improve over DKEF-G-15; we omit it for space. We can gain additional insights into the model by looking at the shape of the learned kernel, shown by the colored lines; the kernels do indeed adapt to the local geometry.

DKEF-L-50 and DKEF-L-50-1 show good performance when the target density has simple geometries, but had trouble in more complex cases, even with much larger networks than used by DKEF-G-15. It seems that a Gaussian kernel with inducing points provides much stronger representational features than a using linear kernel and/or a single inducing point. A large enough network $\bm{\phi}_{\bm{w}}$ would likely be able to perform the task well, but, using currently available software, the second derivatives in the score matching loss limit our ability to use very large networks. A similar phenomenon was observed by Bińkowski et al. (2018) in the context of GAN critics, where combining some analytical RKHS optimization with deep networks allowed much smaller networks to work well.

As expected, models fit by DKEFs generally have smaller Fisher divergences than likelihood-based methods. For Funnel and Banana, the true densities are simple transformations of Gaussians, and the normalizing flow models perform relatively well. But on Ring, Square, and Cosine, the shape of the learned distribution by likelihood-based methods exhibits noticeable artifacts, especially at the boundaries. These artifacts, particularly the “breaks” in Ring, may be caused by a normalizing flow’s need to be invertible and smooth. The shape learned by DKEF-G-15 is much cleaner.

On multimodal distributions with disjoint components, likelihood-based and score matching-based methods show interesting failure modes. The likelihood-based models often connect separated modes with a “bridge”, even for MADE-MOG and MAF-MOG which use explicit mixtures. On the other hand, DKEF is able to find the shapes of the components, but the weighting between modes is unstable. As suggested in Section 3.1, we also fit mixtures of all models (except KCEF) on a partition of MoR found by spectral clustering (Shi & Malik, 2000); DKEF-G-15 produced an excellent fit.

Another challenge we observed in our experiments is that the estimator of the objective function, $\hat{J}$, tends to be more noisy as the model fit improves. This happens particularly on datasets where there are “sharp” features, such as Square (see Figure 4 in Section E.1), where the model’s curvature becomes extreme at some points. This can cause higher variance in the gradients of the parameters, and more difficulty in optimization.

2 Results on Benchmark Datasets

Following recent work on density estimation (Uria et al., 2013; Arbel & Gretton, 2018; Papamakarios et al., 2017; Germain et al., 2015), we trained DKEF and the likelihood-based models on five UCI datasets (Dheeru & Karra Taniskidou, 2017); in particular, we used RedWine, WhiteWine, Parkinson, HepMass, and MiniBoone. All performances were measured on held-out test sets. We did not run KCEF due to its computational expense. Section E.2 gives further details.

Results for models with a fixed Gaussian $q_{0}$ were similar (Figure 6, in Section E.2).

Discussion

Learning deep kernels helps make the kernel exponential family practical for large, complex datasets of moderate dimension. We can exploit the closed-form fit of the ${\bm{\alpha}}$ vector to optimize kernel and even regularization parameters using a “held-out” loss, in a particularly convenient instance of meta-learning. We are thus able to find smoothness assumptions that fit our particular data, rather than arbitrarily choosing them a priori.

Computational expense makes score matching with deep kernels difficult to scale to models with large kernel networks, limiting the dimensionality of possible targets. Combining with the kernel conditional exponential family might help alleviate that problem by splitting the model up into several separate but marginally complex components. The kernel exponential family, and score matching in general, also struggles to correctly balance datasets with nearly-disjoint components, but it seems to generally learn density shapes better than maximum likelihood-based deep approaches.

References

Appendix A DKEFs can be normalized

Consider the kernel $k({\bm{x}},{\bm{y}})=\kappa(\bm{\phi}({\bm{x}}),\bm{\phi}({\bm{y}}))$, where $\kappa$ is a kernel such that $\kappa({\bm{a}},{\bm{a}})\leq L_{\kappa}\lVert{\bm{a}}\rVert^{2}+C_{\kappa}$ and $\bm{\phi}$ a function such that $\lVert\bm{\phi}({\bm{x}})\rVert\leq L_{\bm{\phi}}\lVert{\bm{x}}\rVert+C_{\bm{\phi}}$. Let $q_{0}({\bm{x}})=Q\,r_{0}({\bm{V}}^{-1}({\bm{x}}-{\bm{\mu}}))$, where $Q>0$ is any scalar and $r_{0}$ is a product of independent generalized Gaussian densities, with each $\beta_{d}>1$:

(For example, $\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$ for strictly positive definite ${\bm{\Sigma}}$ could be achieved with $\beta_{d}=2$ and ${\bm{V}}$ the Cholesky factorization of ${\bm{\Sigma}}$.) Then, for any function $f$ in the RKHS $\mathcal{H}$ corresponding to $k$,

First, we have that $f({\bm{x}})=\langle f,k({\bm{x}},\cdot)\rangle_{\mathcal{H}}\leq\lVert f\rVert_{\mathcal{H}}\sqrt{k({\bm{x}},{\bm{x}})},$ and

defining $C_{1}:=\lVert f\rVert_{\mathcal{H}}\sqrt{L_{\kappa}L_{\bm{\phi}}}$, $C_{0}:=\lVert f\rVert_{\mathcal{H}}\sqrt{L_{\kappa}C_{\bm{\phi}}+C_{\kappa}}$.

Let ${\bm{z}}={\bm{V}}^{-1}({\bm{x}}-{\bm{\mu}})$, and let $C_{r}$ be the normalizing constant of $r_{0}$, $C_{q}:=\prod_{d=1}^{D}\frac{{\beta}_{d}}{2\alpha_{d}\operatorname{\Gamma}\left(\frac{1}{\beta_{d}}\right)}$. Then

We can now show that each of these expectations is finite: letting $C=C_{1}\lVert{\bm{V}}\rVert$,

for any $s\in(0,\infty)$. The first integral is clearly finite. Picking $s=\left(2\lvert C\rvert\right)^{\frac{1}{\beta-1}}$, so that $\lvert Cz\rvert<\frac{1}{2}z^{\beta}$ for $z>s$, gives that

The condition on $\bm{\phi}$ holds for any $\bm{\phi}$ given by a deep network with Lipschitz activation functions, such as the softplus function we use in this work. The condition on $\kappa$ also holds for a linear kernel (where $L_{\kappa}=1$, $C_{\kappa}=0$), any translation-invariant kernel ($L_{\kappa}=0$, $C_{\kappa}=\kappa(0,0)$), or mixtures thereof. If $\kappa$ is bounded, the integral is finite for any function $\bm{\phi}$.

The given proof would not hold for a quadratic $\kappa$, which has been used previously in the literature; indeed, it is clearly possible for such an $f$ to be unnormalizable.

Appendix B Finding the optimal 𝜶𝜶{\bm{\alpha}}

We will show a slightly more general result than we need, also allowing for an $\lVert f\rVert_{\mathcal{H}}^{2}$ penalty. This result is related to Lemma 4 of Sutherland et al. (2018), but is more elementary and specialized to our particular needs while also allowing for more types of regularizers.

For fixed $k$, ${\bm{z}}$, and ${\bm{\lambda}}$, as long as $\lambda_{\bm{\alpha}}>0$ then the optimal ${\bm{\alpha}}$ is

We will show that the loss is quadratic in ${\bm{\alpha}}$. Note that

The $\lambda_{C}$ term is of the same form, but with second derivatives:

Because $\lambda_{\alpha}>0$ and ${\bm{G}}$, ${\bm{K}}$, ${\bm{U}}$ are all positive semidefinite, the matrix in parentheses is strictly positive definite, and the claimed result follows directly from standard vector calculus. ∎

Appendix C Behavior on mixtures

Let $\mathcal{D}=\bigcup_{i=1}^{I}\mathcal{D}_{i}$, where $\mathcal{D}_{i}\subset\mathcal{X}_{i}$, $\lvert\mathcal{D}_{i}\rvert=\pi_{i}N$, $\sum_{i=1}^{I}\pi_{i}=1$. Also suppose that the inducing points are partitioned as ${\bm{Z}}=\left[{\bm{Z}}_{1};\dots;{\bm{Z}}_{I}\right]$, with ${\bm{Z}}_{i}\subset\mathcal{X}_{i}$. Further let the kernel $k$ be such that $k({\bm{x}}_{1},{\bm{x}}_{2})=0$ when ${\bm{x}}_{1}\in\mathcal{X}_{i}$, ${\bm{x}}_{2}\in\mathcal{X}_{j}$ for $i\neq j$, with its first and second derivatives also zero. Then the kernel exponential family solution of 3 is

Let ${\bm{G}}_{i}$, ${\bm{b}}_{i}$ be the ${\bm{G}}$, ${\bm{b}}$ of 3 when using only ${\bm{Z}}_{i}$ and $\mathcal{D}_{i}$. Then, because the kernel values and derivatives are zero across components, if $m$ and $m^{\prime}$ are from separate components then

as at least one of the kernel derivatives will be zero for each term of the sum. When $m$ and $m^{\prime}$ are from the same component, the total will be the same except that $N$ is bigger, giving

${\bm{U}}$ is of the same form and factorizes in the same way. ${\bm{K}}$ does not scale:

Each term in the sum for which ${\bm{x}}_{n}$ is in a different component than ${\bm{z}}_{m}$ will be zero, giving ${\bm{b}}=\left(\pi_{1}{\bm{b}}_{1},\cdots,\pi_{I}{\bm{b}}_{I}\right)$. Thus ${\bm{\alpha}}({\bm{\lambda}},k,{\bm{z}},\mathcal{D})$ becomes

Thus the fits for the components are essentially added together, except that each component uses a different $\lambda_{\alpha}$ and $\lambda_{\mathcal{H}}$; smaller components are regularized more. $\lambda_{C}$, interestingly, is unscaled.

It is difficult in general to tell how two components will be weighted relative to one another; the problem is essentially equivalent to computing the overall normalizing constant of a fit. However, we can gain some insight by analyzing a greatly simplified case, in Section C.1.

Consider, for the sake of our study of mixture fits, one of the simplest possible situations for a kernel exponential family: $p_{0}=\operatorname{\mathcal{N}}({\bm{0}},{\bm{I}})$, with a kernel $k({\bm{x}},{\bm{y}})=\exp\left(-\frac{1}{\sigma^{2}}\lVert{\bm{x}}-{\bm{y}}\rVert^{2}\right)$ for $\sigma\gg\sqrt{D}$, so that $k({\bm{x}},{\bm{y}})\approx 1$ for all ${\bm{x}},{\bm{y}}$ sampled from $p_{0}$. Let $q_{0}$ be approximately uniform, $q_{0}=\operatorname{\mathcal{N}}({\bm{0}},q{\bm{I}})$ for $q\gg\sigma^{2}$, so that $\nabla\log q_{0}({\bm{x}})=\frac{-1}{q^{2}}{\bm{x}}\approx 0$. Also assume that $N\to\infty$, but $M$ is fixed. Assume that $\lambda_{\mathcal{H}}=\lambda_{C}=0$, and refer to $\lambda_{\alpha}$ as simply $\lambda$. Then we have

Because $k({\bm{z}},{\bm{x}})\approx k({\bm{z}}^{\prime},{\bm{x}})$ for any ${\bm{z}},{\bm{z}}^{\prime}$ near the data in this setup, it’s sufficient to just consider a single ${\bm{z}}={\bm{0}}$. In that case,

Thus, if we attempt to fit the mixture $\pi\operatorname{\mathcal{N}}({\bm{0}},{\bm{I}})+(1-\pi)\operatorname{\mathcal{N}}({\bm{r}},{\bm{I}})$ with $q^{2}\gg\lVert{\bm{r}}\rVert^{2}\gg\sigma^{2}\gg D$, we are approximately in the regime of 4 and so the ratio between the two components in the fit is

If $\pi=\frac{1}{2}$, the density ratio is correctly 1. If we further assume that $\lambda\gg D/\sigma^{4}$, so that the denominator is dominated by the last term, then the ratio becomes approximately

Thus, depending on the size of $D/(2\sigma^{2}\lambda)\ll\sigma^{2}/2$, the ratio will usually either remain too close to $\frac{1}{2}$ or become too extreme as $\pi$ changes; only in a very narrow parameter range is it approximately correct.

Appendix D Upper bound on normalizer bias

Recall the importance sampling setup of Section 3.2:

First note that the following form of a Taylor expansion holds identically:

We can thus write the bias as the following, where $\varphi(x)=-\log(x)$, $P$ is the distribution of $\hat{Z}$, and $t\geq a$:

is convex in $x$ and thus its supremum is $\max\left(\log\frac{Z}{a}+\frac{a}{Z}-1,\log\frac{Z}{t}+\frac{t}{Z}-1\right)$, with the term at $a$ being necessarily larger as long as $t<Z$.

Picking $t=(s+a)/2$ gives the desired bound on $\Pr(\hat{Z}\leq t)$ via Lemma 5.

Lemma 1 of Liao & Berg (2018) shows that since $\varphi^{\prime}(x)=-1/x$ is concave, $h(x,Z)$ is decreasing in $x$. Thus $\sup_{x\geq t}h(x,Z)=h(t,Z)$, giving the claim. ∎

Let $a$ and $s$ be such that $\Pr({\textnormal{r}}_{u}\geq a)=1$ and $\Pr({\textnormal{r}}_{u}\leq s)\leq\rho<\frac{1}{2}$, with $a<s$. Then

Let $K$ denote the number of samples of ${\textnormal{r}}_{u}$ which are smaller than $s$, so that $U-K$ samples are at least $s$. Then we have

$K$ is distributed binomially with probability of success at most $\rho<\frac{1}{2}$, so applying Theorem 1 of Arratia & Gordon (1989) yields

Let $r:=t/x$, so $x\in[t,\infty)$ corresponds to $r\in(0,1]$, and $x\in(0,t]$ corresponds to $r\in[1,\infty)$. Then

We can evaluate $\lim_{r\to 1}\chi^{\prime\prime}(t/r)=\frac{3}{2t^{4}}>0$. For $r\neq 1$, $\chi_{t}^{\prime\prime}>0$ if and only if $f(r)>0$, where

Clearly $\lim_{r\to 0}f(r)=\infty$ and $f(1)=0$. But notice that

so that $f(r)$ is strictly decreasing on $(0,1)$, and strictly increasing on $(1,\infty)$. Thus $f(r)>0$ for all $r\in(0,1)\cup(1,\infty)$, and $\chi_{t}^{\prime\prime}(x)>0$ for all $x>0$. The claim follows by Jensen’s inequality. ∎

For a kernel such as (7) bounded in $$,$a:=\exp\left(\sum_{m=1}^{M}\min(\alpha_{m},0)\right)$is a uniform lower bound on${\textnormal{r}}_{u}$.

For large $U$, essentially any $\rho<\frac{1}{2}$ will make the second term practically zero, so we select $s$ as slightly less than the 40th percentile of an initial sample of ${\textnormal{r}}_{u}$, and confirm a high-probability ($0.999$) Hoeffding upper bound $\rho$ on $\Pr({\textnormal{r}}_{u}\leq s)$ with another sample. (We use $s$ as $\exp(-0.001)\approx 0.999$ times the estimate of the 40th percentile, to avoid ties.) We use $10^{7}$ samples for each of these.

To finally estimate the bound, we estimate $Z_{\bm{\theta}}$ on yet another independent sample, again usually of size $10^{9}$ but $10^{10}$ for MiniBoone.

Crucially, the function $\psi(t,x)/(x-t)^{2}$ is convex (6); because the variance is unbiased, our estimate of the bias bound is itself biased upwards. As 1’s bound is also not tight, our estimate thus likely overstates the actual amount of bias.

Appendix E Additional experimental details

For each synthetic distribution, we sample $10\,000$ random points from the distribution, $1\,000$ of which are used for testing; of the rest, 90% ($8\,100$) are used for training, and 10% (900) are used for validation. Training was early stopped when validation cost does not improve for 200 minibatches. The current implementation of KCEF does not include a Nyström approximation, and trains via full-batch L-BFGS-B, so we down-sampled the training data to 1000 points. We used the Adam optimizer (Kingma & Ba, 2015) for all other models. For MADE, RealNVP, and MAF, we used minibatches of size 200 and the learning rate was $10^{-3}$ For KEF-G and DKEF, we used 200 inducing points, used $|\mathcal{D}_{t}|=|\mathcal{D}_{v}|=100$, and learning rate $10^{-3}$. The same parameters are used for each component for mixture models trained on MoR.

To show that learning is stable, we ran the experiments on 5 random draws of training, validation and test sets from the synthetic distributions, trained KDEF initialized using 5 random seeds and calculated validation score at each iteration until convergence in the first phase of training (before optimizing for $\lambda$’s). The traces are shown in Figure 4.

The same data for benchmark datasets are shown in Figure 5. There is no overfitting except for the small Redwine dataset. Runs on Hepmass and Miniboone do not seem to fully converge, despite having met the early stopping criterion.

E.2 Benchmark datasets

RedWine and WhiteWine are quantized, and thus problematic for modeling with continuous densities; we added to each dimension uniform noise with support equal to the median distances between two adjacent values. For HepMass and MiniBoone, we removed ill-conditioned dimensions as did Papamakarios et al. (2017). For all datasets except HepMass, 10% of the entire data was used as testing, and 10% of the remaining was used for validation with an upper limit of $1\,000$ due to time cost of validation at each iteration. For HepMass, we used the same splitting as done in Papamakarios et al. (2017) and with the same upper limit on validation set. The data is then whitened before fitting and the whitening matrix was computed on at most $10\,000$ data points.

Likelihood-based models

We set MADE, MADE-MOG, each autoregressive layer of MAF and each scaling and shifting layers of real NVP to have two hidden layers of 100 neurons. For real NVP, MAF and MAF-MOG, five autoregressive layers were used; MAF-MOG and MADE-MOG has a mixture of 10 Gaussians for each conditional distribution. Learning rate was $10^{-3}$ The size of a minibatch is 200.

Deep kernel exponential family

We set the DKEF model to have three kernels ($R=3$), each a Gaussian on features of a 3-layer network with 30 neurons in each layer. There was also a skip-layer connection from data directly to the last layer which accelerated learning. Length scales $\sigma_{r}$ were initialized to 1.0, 3.3 and 10.0. Each $\lambda$ was initialized to 0.001. The weights of the network were initialized from a Gaussian distribution with standard deviation equal to $1/\sqrt{30}$. We also optimized the inducing points ${\bm{z}}_{m}$ which were initialized with random draws from training data. The number of inducing points $M=300$, and $|\mathcal{D}_{t}|=|\mathcal{D}_{v}|=100$. The learning rate was $10^{-2}$. We found that our initialization on the weight std and $\sigma_{r}$’s are importance for fast and stable learning; other parameters did not significantly change the results under similar computational budget (time and memory).

FSSD tests were conducted using 100 points ${\bm{v}}_{b}$ selected at random from the test set, with added normal noise of standard deviation $0.2$, using code provided by the authors.

We estimated $\log Z_{\bm{\theta}}$ with $10^{10}$ samples proposed from $q_{0}$, as in Section 3.2, and estimated the bias as in Section D.1.

We added independent $\operatorname{\mathcal{N}}(0,0.05^{2})$ noise to the data in training. This is similar to the regularization applied by (Kingma & LeCun, 2010; Saremi et al., 2018), except that the noise is added directly to the data instead of the model.

For all models, we stopped training when the objective ((4) or log likelihood) did not improve for 200 minibatches. We also set a time budget of 3 hours on each model; this was fully spent by MAF, MOG-MAF and Real NVP on HepMass. We found that MOG-MADE had unstable runs on some datasets; out of 15 runs on each dataset, 7 on WhiteWine, 4 on Parkinsons and 9 on MiniBoone produced invalid log likelihoods. These results were discarded in Figure 3 log likelihood panels.

The DKEF in our main results (Figure 3) has an adaptive $q_{0}$ which is a generalized normal distribution. We also trained DKEF with $q_{0}$ being an isotropic multivariate normal of standard deviation 2.0. These results Figure 6 are similar to Figure 3 but exhibit much smaller bias estimates in the log normalizer for RedWine and Parkinsons.