Telescoping Density-Ratio Estimation

Introduction

Unsupervised learning via density-ratio estimation is a powerful paradigm in machine learning that continues to be a source of major progress in the field. It consists of estimating the ratio $p/q$ from their samples without separately estimating the numerator and denominator. A common way to achieve this is to train a neural network classifier to distinguish between the two sets of samples, since for many loss functions the ratio $p/q$ can be extracted from the optimal classifier . This discriminative approach has been leveraged in diverse areas such as covariate shift adaptation , energy-based modelling , generative adversarial networks , bias correction for generative models , likelihood-free inference , mutual-information estimation , representation learning , Bayesian experimental design and off-policy reward estimation in reinforcement learning . Across this diverse set of applications, density-ratio based methods have consistently yielded state-of-the-art results.

Despite the successes of discriminative density-ratio estimation, many existing loss functions share a severe limitation. Whenever the ‘gap’ between $p$ and $q$ is large, the classifier can obtain almost perfect accuracy with a relatively poor estimate of the density ratio. We refer to this failure mode as the density-chasm problem—see Figure 1(a) for an illustration. We observe empirically that the density-chasm problem manifests whenever the KL-divergence $D_{KL}(p\;\|\;q)$ exceeds $\sim 20$ nats‘nat’ being a unit of information measured using the natural logarithm (base $e$). This observation accords with recent findings in the mutual information literature regarding the limitations of density-ratio based estimators of the KL . In high dimensions, it can easily occur that two densities $p$ and $q$ will have a KL-divergence measuring in the hundreds of nats, and so the ratio may be virtually intractable to estimate with existing techniques.

In this paper, we propose a new framework for estimating density-ratios that can overcome the density-chasm problem. Our solution uses a ‘divide-and-conquer’ strategy composed of two steps. The first step is to gradually transport samples from $p$ to samples from $q$, creating a chain of intermediate datasets. We then estimate the density-ratio between consecutive datasets along this chain, as illustrated in the top row of Figure 1(b). Unlike the original ratio $p/q$, these ‘chained ratios’ can be accurately estimated via classification (see bottom row). Finally, we combine the chained ratios via a telescoping product to obtain an estimate of the original density-ratio $p/q$. Thus, we refer to the method as Telescoping density-Ratio Estimation (TRE).

We empirically demonstrate that TRE can accurately estimate density-ratios using deep neural networks on high-dimensional problems, significantly outperforming existing single-ratio methods. We show this for two important applications: representation learning via mutual information (MI) estimation and the learning of energy-based models (EBMs).

In the context of mutual information estimation, we show that TRE can accurately estimate large MI values of 30+ nats, which is recognised to be an outstanding problem in the literature . However, obtaining accurate MI estimates is often not our sole objective; we also care about learning representations from e.g. audio or image data that are useful for downstream tasks such as classification or clustering. To this end, our experimental results for representation learning confirm that TRE offers substantial gains over a range of existing single-ratio baselines.

In the context of energy-based modelling, we show that TRE can be viewed as an extension of noise-contrastive estimation that more efficiently scales to high-dimensional data. Whilst energy-based modelling has been a topic of interest in the machine learning community for some time , there has been a recent surge of interest, with a wave of new methods for learning deep EBMs in high dimensions . These methods have shown promising results for image and 3D shape synthesis , hybrid modelling , and modelling of exchangeable data .

However, many of these methods result in expensive/challenging optimisation problems, since they rely on approximate Markov chain Monte Carlo (MCMC) sampling during learning , or on adversarial optimisation . In contrast, TRE requires no MCMC during learning and uses a well-defined, non-adversarial, objective function. Moreover, as we show in our mutual information experiments, TRE is applicable to discrete data, whereas all other recent EBM methods only work for continuous random variables. Applicability to discrete data makes TRE especially promising for domains such as natural language processing, where noise-contrastive estimation has been widely used .

Discriminative ratio estimation and the density-chasm problem

Suppose $p$ and $q$ are two densities for which we have samples, and that $q({\mathbf{x}})>0$ whenever $p({\mathbf{x}})>0$. We can estimate the density-ratio $r({\mathbf{x}})=p({\mathbf{x}})/q({\mathbf{x}})$ by training a classifier to distinguish samples from $p$ and $q$ . There are many choices for the loss function of the classifier , but in this paper we concentrate on the widely used logistic loss

where $r({\mathbf{x}};{\boldsymbol{\theta}})$ is a non-negative ratio estimating model. To enforce non-negativity, $r$ is typically expressed as the exponential of an unconstrained function such as a neural network. For a correctly specified model, the minimiser of this loss, ${\boldsymbol{\theta}}^{*}$, satisfies $r({\mathbf{x}};{\boldsymbol{\theta}}^{*})=p({\mathbf{x}})/q({\mathbf{x}})$, without needing any normalisation constraints . Other classification losses do not always have this self-normalising property, but only yield an estimate proportional to the true ratio—see e.g. .

We experimentally find that density-ratio estimation via classification typically works well when $p$ and $q$ are ‘close’ e.g. the KL divergence between them is less than $\sim 20$ nats. However, for sufficiently large gaps, which we refer to as density-chasms, the ratio estimator is often severely inaccurate. This raises the obvious question: what is the cause of such inaccuracy?

There are many possible sources of error: the use of misspecified models, imperfect optimisation algorithms, and inaccuracy stemming from Monte Carlo approximations of the expectations in (1). We argue that this mundane final point—Monte Carlo error due to finite sample size—is actually sufficient for inducing the density-chasm problem. Figure 1(a) depicts a toy problem for which the model is well-specified, and because it is 1-dimensional (w.r.t. $\theta$), optimisation is straightforward using grid-search. And yet, if we use a sample size of $n=10,000$ and minimise the finite-sample loss

we obtain an estimate $\hat{\theta}$ that is far from the asymptotic minimiser $\theta^{*}={\operatorname{arg}\,\operatorname{min}}\;\mathcal{L}(\theta)$. Repeating this same experiment for different sample sizes, we can empirically measure the method’s sample efficiency, which is plotted as the blue curve in Figure 2. For the regime plotted, we see that an exponential increase in sample size only yields a linear decrease in estimation error. This empirical result is concordant with theoretical findings that density-ratio based lower bounds on KL divergences are only tight for sample sizes exponential in the the number of nats .

Whilst we focus on the logistic loss, we believe the density chasm problem is a broader phenomenon. As shown in the appendix, the issues identified in Figure 1 and the sample inefficiency seen in Figure 2 also occur for other commonly used discriminative loss functions.

Thus, when faced with the density-chasm problem, simply increasing the sample size is a highly inefficient solution and not always possible in practice. This begs the question: is there a more intelligent way of using a fixed set of samples from $p$ and $q$ to estimate the ratio?

Telescoping density-ratio estimation

We introduce a new framework for estimating density-ratios $p/q$ that can overcome the density-chasm problem in a sample-efficient manner. Intuitively, the density-chasm problem arises whenever classifying between $p$ and $q$ is ‘too easy’. This suggests that it may be fruitful to decompose the task into a collection of harder sub-tasks.

For convenience, we make the notational switch $p\equiv p_{0},\ \ q\equiv p_{m}$ (which we will keep going forward), and expand the ratio via a telescoping product

where, ideally, each $p_{k}$ is chosen such that a classifier cannot easily distinguish it from its two neighbouring densities. Instead of attempting to build one large ‘bridge’ (i.e. density-ratio) across the density-chasm, we propose to build many small bridges between intermediate ‘waymark’ distributions. The two key components of the method are therefore:

Waymark creation. We require a method for gradually transporting samples $\{{\mathbf{x}}^{1}_{0},\ldots,{\mathbf{x}}^{n}_{0}\}$ from $p_{0}$ to samples $\{{\mathbf{x}}^{1}_{m},\ldots,{\mathbf{x}}^{n}_{m}\}$ from $p_{m}$. At each step in the transportation, we obtain a new dataset $\{{\mathbf{x}}^{1}_{k},\ldots,{\mathbf{x}}^{n}_{k}\}$ where $k\in\{0,\ldots m\}$. Each intermediate dataset can be thought of as samples from an implicit distribution $p_{k}$, which we refer to as a waymark distribution.

Bridge-building: A method for learning a set of parametrised density-ratios between consecutive pairs of waymarks $r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})\approx p_{k}({\mathbf{x}})/p_{k+1}({\mathbf{x}})$ for $k=0,\ldots,m-1$, where each bridge $r_{k}$ is a non-negative function. We refer to these ratio estimating models as bridges. Note that the parameters of the bridges, $\{{\boldsymbol{\theta}}_{k}\}_{k=0}^{m-1}$, can be totally independent or they can be partially shared.

An estimate of the original ratio is then given by the product of the bridges

where ${\boldsymbol{\theta}}$ is the concatenation of all ${\boldsymbol{\theta}}_{k}$ vectors. Because of the telescoping product in (4), we refer to the method as Telescoping density-Ratio Estimation (TRE).

TRE has conceptual ties with a range of methods in optimisation, statistical physics and machine learning that leverage sequences of intermediate distributions, typically between a complex density $p$ and a simple tractable density $q$. Of particular note are the methods of Simulated Annealing , Bridge Sampling & Path Sampling and Annealed Importance Sampling (AIS) . Whilst none of these methods estimate density ratios, and thus serve fundamentally different purposes, they leverage similar ideas. In particular, AIS also computes a chain of density-ratios between artificially constructed intermediate distributions. It typically does this by first defining explicit expressions for the intermediate densities, and then trying to obtain samples via MCMC. In contrast, TRE implicitly defines the intermediate distributions via samples and then tries to learn the ratios. Additionally, in TRE we would like to evaluate the learned ratios in (4) at the same input ${\mathbf{x}}$ while AIS should only evaluate a ratio $r_{k}$ at ‘local’ samples from e.g. $p_{k}$.

In this paper, we consider two simple, deterministic waymark creation mechanisms: linear combinations and dimension-wise mixing. We find these mechanisms yield good performance and are computationally cheap. However, we note that other mechanisms are possible, and are a promising topic for future work.

Linear combinations. Given a random pair ${\mathbf{x}}_{0}\sim p_{0}$ and ${\mathbf{x}}_{m}\sim p_{m}$, define the $k$th waymark via

where the $\alpha_{k}$ form an increasing sequence from to $1$, which control the distance of ${\mathbf{x}}_{k}$ from ${\mathbf{x}}_{0}$. For all of our experiments (except, for illustration purposes, those depicted in Figure 1), each dimension of $p_{0}$ and $p_{m}$ has the same varianceFor MI estimation this always holds, for energy-based modelling this is enforceable via the choice of $p_{m}$. and the coefficients in (5) are chosen to preserve this variance, with the goal being to match basic properties of the waymarks and thereby make consecutive classification problems harder.

Dimension-wise mixing. An alternative way to ‘mix’ two vectors is to concatenate different subsets of their dimensions. Given a $d$-length vector ${\mathbf{x}}$, we can partition it into $m$ sub-vectors of length $d/m$, assuming $d$ is divisible by $m$. We denote this as ${\mathbf{x}}=({\mathbf{x}},\ldots,{\mathbf{x}}[m])$, where each ${\mathbf{x}}[i]$ has length $d/m$. Using this notation, define the $k$th waymark via

where, again, ${\mathbf{x}}_{0}\sim p_{0}$ and ${\mathbf{x}}_{m}\sim p_{m}$ are randomly paired.

Number and spacing. Given these two waymark generation mechanisms, we still need to decide the number of waymarks, $m$, and, in the case of linear combinations, how the $\alpha_{k}$ are spaced in the unit interval. We treat these quantities as hyperparameters, and demonstrate in the experiments (Section 4) that tuning them is feasible with a limited search budget.

2 Bridge-building

Each bridge $r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})$ in (4) can be learned via binary classification using a logistic loss function as described in Section 2. Solving this collection of classification tasks is therefore a multi-task learning (MTL) problem—see for a review. Two key questions in MTL are how to share parameters and how to define a joint objective function.

Parameter sharing. We break the construction of the bridges $r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})$ into two stages: a (mostly) shared body computing hidden vectors $f_{k}({\mathbf{x}})$For simplicity, we suppress the parameters of $f_{k}$, and will do the same for $r_{k}$ in the experiments section., followed by bridge-specific heads. The body $f_{k}$ is a deep neural network with shared parameters and pre-activation per-hidden-unit scales and biases for each bridge (see appendix for details). Similar parameter sharing schemes have been successfully used in the multi-task learning literature . The heads map the hidden vectors $f_{k}({\mathbf{x}})$ to the scalar $\log r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})$. We use either linear or quadratic mappings depending on the application; the precise parameterisation is stated in each experiment section.

TRE loss function. The TRE loss function is given by the average of the $m$ logistic losses

This simple unweighted average works well empirically. More sophisticated multi-task weighting schemes exist , but preliminary experiments suggested they were not worth the extra complexity.

An important aspect of this loss function is that each ratio estimator $r_{k}$ sees different samples during training. In particular, $r_{0}$ sees samples close to the real data i.e. from $p_{0}$ and $p_{1}$, while the final ratio $r_{m-1}$ sees data from $p_{m-1}$ and $p_{m}$. This creates a potential mismatch between training and deployment, since after learning, we would like to evaluate all ratios at the same input ${\mathbf{x}}$. In our experiments, we do not find this mismatch to be a problem, suggesting that each ratio, despite seeing different inputs during training, is able to generalise to new test points. We speculate that this generalisation is encouraged by parameter sharing, which allows each ratio-estimator to be indirectly influenced by samples from all waymark distributions. Nevertheless, we think a deeper analysis of this issue of generalisation deserves further work.

3 TRE applied to mutual information estimation

The mutual information (MI) between two random variables ${\mathbf{u}}$ and ${\mathbf{v}}$ can be written as

Given samples from the joint density $p({\mathbf{u}},{\mathbf{v}})$, one obtains samples from the product-of-marginals $p({\mathbf{u}})p({\mathbf{v}})$ by shuffling the ${\mathbf{v}}$ vectors across the dataset. This then enables standard density-ratio estimation to be performed.

For TRE, we require waymark samples. To generate these, we take a sample from the joint, ${\mathbf{x}}_{0}=({\mathbf{u}},{\mathbf{v}}_{0})$, and a sample from the product-of-marginals, ${\mathbf{x}}_{m}=({\mathbf{u}},{\mathbf{v}}_{m})$, where ${\mathbf{u}}$ is held fixed and only ${\mathbf{v}}$ is altered. We then apply a waymark construction mechanism from Section 3.1 to generate ${\mathbf{x}}_{k}=({\mathbf{u}},{\mathbf{v}}_{k})$, for $k=0,\ldots,m$.

4 TRE applied to energy-based modelling

An energy-based model (EBM) is a flexible parametric family $\{\phi({\mathbf{x}};{\boldsymbol{\theta}})\}$ of non-negative functions, where each function is proportional to a probability-density. Given samples from a data distribution with density $p({\mathbf{x}})$, the goal of energy-based modelling is to find a parameter ${\boldsymbol{\theta}}^{*}$ such that $\phi({\mathbf{x}};{\boldsymbol{\theta}}^{*})$ is ‘close’ to $cp({\mathbf{x}})$, for some positive constant $c$.

In this paper, we consider EBMs of the form $\phi({\mathbf{x}};{\boldsymbol{\theta}})=r({\mathbf{x}};{\boldsymbol{\theta}})q({\mathbf{x}})$, where $q$ is a known density (e.g. a Gaussian or normalising flow) that we can sample from, and $r$ is an unconstrained positive function. Given this parameterisation, the optimal $r$ simply equals the density-ratio $p({\mathbf{x}})/q({\mathbf{x}})$, and hence the problem of learning an EBM becomes the problem of estimating a density-ratio, which can be solved via TRE. We note that, since TRE actually estimates a product of ratios as stated in Equation 4, the final EBM will be a product-of-experts model of the form $\phi({\mathbf{x}};{\boldsymbol{\theta}})=\prod_{k=0}^{m-1}r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})q({\mathbf{x}})$.

The estimation of EBMs via density-ratio estimation has been studied in multiple prior works, including noise-contrastive estimation (NCE) , which has many appealing theoretical properties . Following NCE, we will refer to the known density $q$ as the ‘noise distribution’.

Experiments

We include two toy examples illustrating both the correctness of TRE and the fact that it can solve problems which verge on the intractable for standard density ratio estimation. We then demonstrate the utility of TRE on two high-dimensional complex tasks, providing clear evidence that it substantially improves on standard single-ratio baselines.

For experiments with continuous random variables, we use the linear combination waymark mechanisms in (5); otherwise, for discrete variables, we use dimension-wise mixing (6). For the linear combination mechanism, we collapse the $\alpha_{k}$ into a single spacing hyperparameter, and grid-search over this value, along with the number of waymarks. Details are in the appendix.

Figure 2 shows the full results. These sample efficiency curves clearly demonstrate that, across all sample sizes, TRE is significantly more accurate than single ratio estimation. In fact, TRE obtains a better solution with 100 samples than single-ratio estimation does with 100,000 samples: a three orders of magnitude improvement.

2 High-dimensional ratio with large MI

We apply TRE using quadratic bridges of the form: $\log r_{k}({\mathbf{x}})={\mathbf{x}}^{T}{\mathbf{W}}_{k}{\mathbf{x}}+b_{k}$. The results in Figure 3 show that single ratio estimation becomes severely inaccurate for MI values greater than 20 nats. In contrast, TRE can accurately estimate MI values as large as 80 nats for 320 dimensional variables. To our knowledge, TRE is the first discriminative MI estimation method that can scale this gracefully.

3 MI estimation & representation learning on SpatialMultiOmniglot

We applied TRE to the SpatialMultiOmniglot problem taken from We mirror their experimental setup as accurately as possible, however we were unable to obtain their code. where characters from Omniglot are spatially stacked in an $n\times n$ grid, where each grid position contains characters from a fixed alphabet. Following , the individual pixel values of the characters are not considered random variables; rather, we treat the grid as a collection of $n^{2}$ categorical random variables whose realisations are the characters from the respective alphabet. Pairs of grids, $({\mathbf{u}},{\mathbf{v}})$, are then formed such that corresponding grid-positions contain alphabetically consecutive characters. Given this setup, the ground truth MI can be calculated (see appendix).

Each bridge in TRE uses a separable architecture given by $\log r_{k}({\mathbf{u}},{\mathbf{v}})=g({\mathbf{u}})^{T}{\mathbf{W}}_{k}f_{k}({\mathbf{v}})$, where $g$ and $f_{k}$ are 14-layer convolutional ResNets and $f_{k}$ uses the parameter-sharing scheme described in Section 3.2. We note that separable architectures are standard in the MI-based representation learning literature . We construct waymarks using the dimension-wise mixing mechanism (6) with $m=n^{2}$ (i.e. one dimension is mixed at a time).

After learning, we adopt a standard linear evaluation protocol (see e.g. ), where we train supervised linear classifiers on top of the output layer $g({\mathbf{u}})$ to predict the alphabetic position of each character in ${\mathbf{u}}$. We compare our results to those reported in . Specifically, we report their baseline method—contrastive predictive coding (CPC) , a state-of-the-art representation learning method based on single density-ratio estimation—along with their variant, Wasserstein predictive coding (WPC).

Figure 4 shows the results. The left plot shows that only TRE can accurately estimate high MI values of $\sim 35$ nats do not provide MI estimates for CPC & WPC, but shows that they are bounded by log batch-size.. The representation learning results (right) show that all single density-ratio baselines degrade significantly in performance as we increase the number of characters in a grid (and hence increase the MI). In contrast, TRE always obtains greater than $97\%$ accuracy.

4 Energy-based modelling on MNIST

As explained in Section 3.4, TRE can be used estimate an energy-based model of the form $\phi({\mathbf{x}};{\boldsymbol{\theta}})=\prod_{k=0}^{m-1}r_{k}({\mathbf{x}};{\boldsymbol{\theta}}_{k})q({\mathbf{x}})$, where $q$ is a pre-specified ‘noise’ distribution from which we can sample, and the product of ratios is given by TRE. In this section, we demonstrate that such an approach can scale to high-dimensional data, by learning energy-based models of the MNIST handwritten digit dataset . We consider three choices of the noise distribution: a multivariate Gaussian, a Gaussian copula and a rational-quadratic neural spline flow (RQ-NSF) with coupling layers . Each distribution is first fitted to the data via maximum likelihood estimation—see appendix for details.

For a Gaussian, $F$ is linear, and hence (10) is identical to the original waymark mechanism in (5).

We use the parameter sharing scheme from Section 3.2 together with quadratic heads. This gives $\log r_{k}({\mathbf{x}})=-f_{k}({\mathbf{x}})^{T}{\mathbf{W}}_{k}f_{k}({\mathbf{x}})-f_{k}({\mathbf{x}})^{T}{\mathbf{b}}_{k}-c_{k}$, where we set $f_{k}$ to be an 18-layer convolutional Resnet and constrain ${\mathbf{W}}_{k}$ to be positive definite. This constraint enforces an upper limit on the log-density of the EBM, which has been useful in other work , and improves results here.

We evaluate the learned EBMs quantitatively via estimated log-likelihood in Table 1 and qualitatively via random samples from the model in Figure 5. For both of these evaluations, we employ NUTS to perform annealed MCMC sampling as explained in the appendix. This annealing procedure provides two estimators of the log-likelihood: the Annealed Importance Sampling (AIS) estimator and the more conservative Reverse Annealed Importance Sampling Estimator (RAISE) .

The results in Table 1 and Figure 5 show that single ratio estimation performs poorly in high-dimensions for simple choices of the noise distribution, and only works well if we use a complex neural density-estimator (RQ-NSF). This illustrates the density-chasm problem explained in Section 2. In contrast, TRE yields improvements for all choices of the noise, as measured by the approximate log-likelihood and the visual fidelity of the samples. TRE’s improvement over the Gaussian noise distribution is particularly large: the bits per dimension (bpd) is around 0.66 lower, corresponding to an improvement of roughly $360$ nats. Moreover, the samples are significantly more coherent, and appear to be of higher fidelity than the RQ-NSF samplesWe emphasise here that the quality of the RQ-NSF model depends on the exact architecture. A larger model may yield better samples. Thus, we do not claim that TRE generally yields superior results in any sense., despite the fact that TRE (with Gaussian noise) has a worse log-likelihood. This final point is not contradictory since log-likelihood and sample quality are known to be only loosely connected .

Finally, we analysed the sensitivity of our results to the construction of the waymarks and include the results in the appendix. Using TRE with a copula noise distribution as an illustrative case, we found that varying the number of waymarks between 5-30 caused only minor changes in the approximate log-likelihoods, no greater than $0.03$ bpd. We also found that if we omit the ${\mathbf{z}}$-space waymark mechanism in (10), and work in ${\mathbf{x}}$-space, then TRE’s negative log-likelihood increases to $1.33$ bpd, as measured by RAISE. This is still significantly better than single-ratio estimation, but does show that the quality of the results depends on the exact waymark mechanism.

Conclusion

We introduced a new framework—Telescoping density-Ratio Estimation (TRE)—for learning density-ratios that, unlike existing discriminative methods, can accurately estimate ratios between extremely different densities in high-dimensions.

TRE admits many exciting directions for future work. Firstly, we would like a deeper theoretical understanding of why it is so much more sample-efficient than standard density-ratio estimation. The relationship between TRE and standard methods is structurally similar to the relationship between annealed importance sampling and standard importance sampling. Thus, exploring this connection further may be fruitful. Relatedly, we believe that TRE would benefit from further research on waymark mechanisms. We presented simple mechanisms that have clear utility for both discrete and continuous-valued data. However, we suspect more sophisticated choices may yield improvements, especially if one can leverage domain or task-specific assumptions to intelligently decompose the density-ratio problem. Lastly, whilst this paper has focused on the logistic loss, it would be interesting to more deeply investigate TRE with other discriminative loss functions.

Broader Impact

As outlined in the introduction, density-ratio estimation is a foundational tool in machine learning with diverse applications. Our work, which improves density-ratio estimation, may therefore increase the scope and power of a wide spectrum of techniques used both in research and real-world settings. The broad utility of our contribution makes it challenging to concretely assess the societal impact of the work. However, we do discuss here two applications of density-ratio estimation with obvious potential for positive & negative impacts on society.

Generative Adversarial Networks are a popular class of models which are often trained via density-ratio estimation and are able to generate photo-realistic image/video content. To the extent that TRE can enhance GAN training (a topic we do not treat in this paper), our work could conceivably lead to enhanced ‘deepfakes’, which can be maliciously used in fake-news or identity fraud.

More positively, density-ratio estimation is being used to correct for dataset bias, including the presence of skewed demographic factors like race and gender . While we are excited about such applications, we emphasise that density-ratio based methods are not a panacea; it is entirely possible for the technique to introduce new biases when correcting for existing ones. Future work should continue to be mindful of such a possibility, and look for ways to address the issue if it arises.

Acknowledgments and Disclosure of Funding

Benjamin Rhodes was supported in part by the EPSRC Centre for Doctoral Training in Data Science, funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016427/1) and the University of Edinburgh. Kai was supported by Edinburgh Huawei Research Lab in the University of Edinburgh, funded by Huawei Technologies Co. Ltd.

References

Appendix A ResNet architectures with parameter sharing

In Figure 6, we give the exact architectures for the $f_{k}$ used in the two high-dimensional experiments on SpatialMultiOmniglot and MNIST. These $f_{k}$ output a hidden vector for the $k$th bridge, which is then mapped to the scalar value of the log-ratio, as stated in each experiment section. All convolution operations share their parameters across the bridges, and are thus independent of $k$.

The only difference between our conditional residual blocks (i.e. ‘CondResBlocks’) and a standard residual block is the use of ‘ConditionalScaleShift’ layers. These layers map a hidden vector ${\mathbf{z}}_{k}$ to a hidden vector of the same size, ${\mathbf{z}}^{\prime}_{k}$, via

where ${\mathbf{s}}_{k}$ and ${\mathbf{b}}_{k}$ are bridge-specific parameters and $\odot$ denotes element-wise multiplication. This operation could be thought of as class-conditional Batch Normalisation (BN) without the normalisation. We did not investigate the use of BN, since many energy-based modelling papers (e.g. ) found it to harm performance. We did perform preliminary experiments with Instance Normalisation in the context of energy-based modelling, finding it to be harmful to performance.

For the MNIST energy-based modelling experiments, we use average pooling operations since other work has found this to produce higher quality samples than max pooling. For the SpatialMultiOmniglot experiments, we grid-search over average pooling and max pooling. For both sets of experiments, we use LeakyRelu activations with a slope of $0.3$.

The MNIST architecture includes an attention block which has been used in GANs to model long-range dependencies in the input image. We found that this attention layer did not yield improvements in estimated log-likelihood, but we think it may yield slightly more globally coherent samples. We note that that another commonly used feature in recent GAN and EBM architectures is Spectral Normalisation (SN) . Our preliminary experiments suggested that SN was not beneficial for performance. That said, all of our negative results should be taken with a grain of salt, given the preliminary nature of the experiments.

Appendix B Waymark number and spacing

As stated in the main text, the number and (in the case of linear combinations) the spacing of the waymarks are treated as hyperparameters. Finding good values of these hyperparameters is made simpler by the following observations.

If any of the TRE logistic losses saturate close to 0 during learning, then this indicates that the density-chasm problem has occured for that bridge, and we can terminate the run.

As illustrated by our sensitivity analysis for MNIST (see Figure 10) it seems that, past a certain point, performance plateaus with the addition of extra waymarks. The fact that it plateaus, and does not decrease, is good news since it indicates that there is little risk of ‘overshooting’, and obtaining a bad model by having too many waymarks.

We now recall the linear combinations waymark mechanism, given by

where $m$ is the number of waymarks. We consider two ways of reducing the coefficients $\alpha_{k}$ to a function of a single spacing hyperparameter $p$ via

Both mechanisms yield linearly spaced $\alpha_{k}$ when $p=1$. For the first mechanism in (13), setting $p>1$ means the gaps between waymarks increase with $k$ (and conversely decrease if $p<1$). The spacing mechanism in (16) is a kind of symmetrised version of (13).

Table 2 shows the grid-searches we performed for all experiments. We note that these weren’t always all performed in parallel. When using linear combinations, we typically set $p=1$ initially and searched over values of $m$. If, for all values of $m$ tested, one of the TRE logistic losses saturated close to 0, then we would expand our search space and test different values of $p$.

Appendix C Minibatching

Recall that the TRE loss is a sum of logistic losses:

When generating minibatch estimates of this loss, we can either sample from each $p_{k}$ independently, or we can couple the samples. By ‘couple’, we mean first drawing $B$ samples each from $p_{0}$ and $p_{m}$, randomly pairing members from each set, and then, for each pair, constructing all possible intermediate waymark samples to obtain a final minibatch of size $B\times M$. Coupling in this way means that the gradient of (17) w.r.t. ${\boldsymbol{\theta}}$ is estimated using shared sources of randomness, which can act as a form of variance reduction .

In all of our experiments, we use coupling when forming minibatches, since we found it to be useful in some preliminary investigations. However, coupling does have memory costs: the number of independent samples drawn from the data distribution, $B$, may need to be very small for the full minibatch, $B\times M$, to fit into memory. We speculate that as $B$ becomes sufficiently small, coupled minibatches will produce inferior results to non-coupled minibatches (which can use a greater number of independent real data samples). Empirical investigation of this claim is left to future work.

Appendix D 1d peaked ratio toy experiment

In this experiment we estimate the ratio $p_{0}/p_{m}$, where both densities are Gaussian, $p_{0}=\mathcal{N}(0,\sigma_{0}^{2})$ and $p_{m}=\mathcal{N}(0,\sigma_{m}^{2})$, where $\sigma_{0}=10^{-6}$ and $\sigma_{m}=1$. We generate waymarks using the linear combinations mechanism (12), which implies that each waymark distribution is Gaussian, since linear combinations of Gaussian random variables are also Gaussian. Specifically, the waymark distributions have the form

where the $\sigma_{k}$ form an increasing sequence between $\sigma_{0}$ and $\sigma_{m}$. The log-ratio between two waymark distributions is therefore given by

where the quadratic coefficient $-\exp(\theta_{k})$ is always negative. We note that this model is well-specified since it contains the ground-truth solution in (20).

The bridges can then be combined via summation to provide an estimate of the original log-ratio

Where $\theta_{TRE}=\log(\sum_{k=0}^{m-1}\exp(\theta_{k}))$. We observe that (24) has the same form as (21) if we were to set $m=1$ in (21) (i.e. if we use a single bridge). Hence $\theta_{TRE}$ can be directly compared to the parameter value we would obtain if we used single density-ratio estimation. This is precisely the comparison we make in Figure 1a and Figure 2 of the main text.

In the main paper, we illustrated the density-chasm problem for the logistic loss using the 1d peaked ratio experiment. Here, we illustrate precisely the same phenomenon for the NWJ/MINE-f loss and a Least Squares (LSQ) loss used by . The loss functions are given by

where the $\sigma$ in (26) denotes the sigmoid function.

In Figures 8 & 8, we can see how single-density ratio estimation performs when using the NWJ and LSQ loss functions for 10,000 samples. the loss curves display the same ‘saturation’ effect seen for the logistic loss, where many settings of the parameter yield an almost identical value of the loss. Moreover, the minimiser of these saturated objectives is far from the ‘true’ minimiser (black dotted lines).

Figures 8 & 8 also show the performance of TRE when each bridge is estimated using the NWJ/LSQ losses. Each TRE loss has a quadratic bowl shape, where the finite-sample minimisers almost perfectly overlap with the true minimisers.

Finally, we plot sample efficiency curves for both the NWJ and LSQ losses, showing the results in Figure 9. We see that single density-ratio estimation with NWJ or LSQ performs poorly, with at best linear gains for exponential increases in sample size. In contrast, if we perform TRE using NWJ or LSQ losses, then we obtain significantly better performance with orders of magnitude fewer samples. These findings are essentially the same as those presented in the main paper for the logistic loss.

Appendix E High-dimensional ratio with large MI toy experiment

We generate $100,000$ samples for each of the train/validation/test splits. We use a total batch size of $1024$, which includes all samples from the waymark trajectories. The bridges in TRE have the form $\log r_{k}({\mathbf{x}})={\mathbf{x}}^{T}{\mathbf{W}}_{k}{\mathbf{x}}+b_{k}$, where we enforce that the diagonal entries of ${\mathbf{W}}_{k}$ are positive and that the matrix is symmetric. We use the Adam optimiser with an initial learning rate of $0.0001$ for TRE, and $0.0005$ for single ratio estimation. We use the default Tensorflow settings for $\beta_{1},\beta_{2}$ and $\epsilon$. We train the models for $40,000$ iterations, which takes at most 1 hour.

Appendix F MI estimation & representation learning on SpatialMultiOmniglot

We here describe how we created the SpatialMultiOmniglot dataset and give the derivation for the ground truth mutual information values presented in the main paperThe original work from which we borrow this experiment did not not provide a detailed explanation or code.. We will share the dataset, along with code for the paper, upon publication. We also state the hyperparameter settings used in our experiments.

We take the Tensorflow version of the Omniglot dataset (https://www.tensorflow.org/datasets/catalog/omniglot) and resize it to $28\times 28$ using the tf.image.resize function. We arrange the data into alphabets $\{A_{i}\}_{i=1}^{l}$, where each alphabet contains $n_{i}$ characters. The alphabets are sorted by size, so that $n_{1}>n_{2}>\ldots>n_{l}$. Each character in a alphabet has $20$ different versions (e.g. there are 20 different images depicting the letter ‘w’). Hence, we can express each alphabet as a set $A_{i}=\{\{a_{j,k}^{i}\}_{k=1}^{20}\}_{j=1}^{n_{i}}$, where $a_{j,k}^{i}$ refers to the $k$th version of the $j$th character of the $i$th alphabet.

In order to construct the $d$-dimensional version of the SpatialMultiOmniGlot dataset, we restrict ourselves to the $d$ largest alphabets $\{A_{i}\}_{i=1}^{d}$. We then sample a vector of categorical random variables

where the $i$th categorical distribution is uniform over the set $\{1,\ldots,n_{i}\}$ and is independent from the other categorical distributions. The vector ${\mathbf{j}}$ should be thought of as an index vector that specifies a particular character from each of the $d$ alphabets.

We then sample two i.i.d random variables ${\mathbf{k}}$ and ${\mathbf{k}}^{\prime}$, via

where, again, each Categorical distribution is independent from the rest. These vectors should be thought of as index vectors that specify a particular version of a character.

Now, we define a datapoint as a tuple ${\mathbf{x}}=({\mathbf{u}},{\mathbf{v}})$, where

In words, we construct ${\mathbf{u}}$ and ${\mathbf{v}}$ such that $u_{i}$ and $v_{i}$ are consecutive characters within their alphabet (whilst the precise versions of the characters are randomised). Finally, we arrange ${\mathbf{u}}$ and ${\mathbf{v}}$ into a grid using raster ordering. This is possible since we assume $d$ to be a square number.

Importantly, we emphasise that ${\mathbf{u}},{\mathbf{v}}\in\prod_{i=1}^{d}A_{i}$ are discrete random variables defined over a set of template images. They are not defined over a space of pixel values, as is usually the case in image-modelling.

F.2 Derivation of ground truth MI

By construction, we have that ${\mathbf{u}}$ and ${\mathbf{v}}$ are conditionally independent given ${\mathbf{j}}$. This means

Furthermore, will assume that, for all ${\mathbf{u}}$ there exists a unique ${\mathbf{j}}_{{\mathbf{u}}}$ such that

Similarly, for any ${\mathbf{v}}$, there exists a unique ${\mathbf{j}}_{{\mathbf{v}}}$ satisfying the same condition. In words, this simply means that, given a grid of Omniglot images, we assume there is no ambiguity about which characters are present. Using Bayes’ rule, and the fact that for a given ${\mathbf{j}}$, ${\mathbf{u}}$ is uniquely determined by ${\mathbf{k}}$, one can then deduce that

We now proceed to derive an analytical formula for the ground truth mutual information between ${\mathbf{u}}$ and ${\mathbf{v}}$. We show that the mutual information is equal to the sum of the log alphabet sizes $\mathcal{I}({\mathbf{u}},{\mathbf{v}})=\sum_{i=1}^{d}\log n_{i}$.

F.3 Experimental settings

We generate 3 versions of the SpatialMultiOmniglot dataset for $d=1,4,9$. For each version, we sample $50,000$ training points and $10,000$ validation and test points. As stated in the main text, we use a separable architecture given by

where $f_{k}$ is a convolutional ResNet whose architecture is given in Figure 6(a). The function $g$ is also a convolutional ResNet with almost the same architecture, except that none of its parameters are bridge-specific, and hence the ‘ConditionalScaleShift’ layers simply become ‘ScaleShift’ layers, with no dependence on $k$.

To construct a mini-batch, we first sample a batch from the joint distribution $p({\mathbf{u}},{\mathbf{v}})$. We then obtain samples from $p({\mathbf{u}})p({\mathbf{v}})$ by sampling a second batch from the joint distribution (which could overlap with the first batch), and shuffling the ${\mathbf{v}}$ vectors across this second batch. Finally, we construct waymark trajectories as described in the main text. For all experiments, the ‘total’ batch size is $\sim 512$, which includes all samples from the waymark trajectories. Thus, as the number of waymarks increases, the number of trajectories in a batch decreases.

We use the Adam optimiser with an initial learning rate of $10^{-4}$ with default Tensorflow settings for $\beta_{1},\beta_{2}$ and $\epsilon$. We gradually decrease the learning rate over the course of training with cosine annealing . All models are trained using a single NVIDIA Tesla P100 GPU card for $200,000$ iterations, which takes at most a day.

We grid-searched over the type of pooling (max vs. average) and the size of the final dense layer ($150n,300n$ and $450n$, where $d=n^{2}$). Interestingly, average pooling was less prone to overfitting and often yielded better final performance, however it was often ‘slow to get started’, with the TRE losses hardly making any progress during the first quarter of training.

For the representation learning evaluations, we first obtained the hidden representations $g({\mathbf{u}})$ for the entire dataset. We then trained a collection of independent supervised linear classifiers on top of these representations, in order to predict the alphabetic position of each character in ${\mathbf{u}}$. We used the L-BFGS optimiser to fit these classifiers via the tfp.optimizer.lbfgs_minimize function, setting the maximum iteration number to $10,000$.

Appendix G Energy-based modelling on MNIST

We here discuss the parameterisation of the noise distributions used in the experiments, the exact method for sampling from the learned EBMs, and the experimental settings used for TRE.

For all noise distributions and TRE models, we use the Adam optimiser with an initial learning rate of $10^{-4}$ with default Tensorflow settings for $\beta_{1},\beta_{2}$ and $\epsilon$. We gradually decrease the learning rate over the course of training with cosine annealing . All models are trained using a single NVIDIA Tesla P100 GPU card.

As stated in the main text, we consider three noise distributions: a multivariate Gaussian, a Gaussian copula and a rational-quadratic neural spline flow (RQ-NSF), all of which are pre-trained via maximum likelihood estimation.

The full-covariance multivariate Gaussian is by far the simplest, and can be fitted in around a minute via np.cov. The Gaussian copula is slightly more complicated. Its density can be written as $p({\mathbf{x}})=\mathcal{N}([s_{1}(x_{1}),\ldots,s_{d}(x_{d})];\mu,\Sigma)\prod_{i=1}^{d}|s^{\prime}_{i}(x_{i})|$. The $s_{i}$ are given by the composition of the inverse CDF of a standard normal and the CDF of the univariate $x_{i}$. It is possible to exploit this to learn the $s_{i}$—as well as $\mu$ and $\Sigma$—however, we found it slightly simpler to directly parametrise the $s_{i}$ via flexible rational-quadratic spline functions of which there are official implementations in Tensorflow and Pytorch and to jointly learn all parameters via maximum likelihood. We follow the basic hyperparameter recommendations in . The hyperparameters that required tuning were the number of bins (we use 128) and the interval widths (which we set to 3 times the standard deviation of the data). For optimisation, we used a batch size of 512 and trained for $40,000$ iterations.

Finally, we turn to the RQ-NSF model . We largely adopt the architectural choices of , and so for a more detailed explanation, we refer the reader to their work. We use a multi-scale convolutional architecture comprised of 2 levels, where each level contains 8 ‘steps’. A step consists of an actnorm layer, an invertible $1\times 1$ convolution, and a rational-quadratic coupling transform. The coupling transforms are parameterised by a block of convolution operations following , which use 64 feature maps. The spline functions use 8 bins and the interval width is set to . We do not ‘factor out’ half of the variables at the end of each level, but do perform ‘squeeze’ operation and an additional $1\times 1$ convolution. For optimisation, we set the batch size to 256, the dropout rate to 0.1, and train for $200,000$ iterations, which takes under a day.

G.2 Annealed MCMC Sampling

We here describe how we leverage the specific products-of-experts structure of the TRE model to perform annealed MCMC sampling. Firstly, we initialise a set of MCMC chains with i.i.d samples from the noise distribution $p_{m}$. We could then run an MCMC sampler with the full TRE model as the target distribution. However, we instead use an annealing procedure, whereby we iteratively sample from a sequence of distributions that interpolate between $p_{m}$ and $p_{0}$. Such distributions can be obtained by multiplying $p_{m}$ with an increasing number of bridges

To obtain an even smoother interpolation, we further define exponentially-averaged intermediate distributions $p_{k,t}({\mathbf{x}})=p_{k}({\mathbf{x}})^{\beta_{t}}p_{k+1}({\mathbf{x}})^{1-\beta_{t}}$, where $\{\beta_{t}\}$ is a decreasing sequence of numbers ranging from 1 to 0.

In addition to obtaining samples, we can simultaneously use this annealing procedure for estimating the log-likelihood of the model via annealed importance sampling (AIS) . We may also run the annealing procedure ‘in reverse’, initialising a chain at a datapoint and iteratively removing bridges until the target distribution of the MCMC sampler is the noise distribution. Using this reverse sampling procedure, we can obtain a second, more conservative, estimate of the log-likelihood via the reverse annealed importance sampling estimator (RAISE) .

Whilst in principle any MCMC sampler could be used, the efficiency of different samplers can vary greatly. We choose to use the gradient-based No-U-turn sampler (NUTS) , which is a highly efficient method for many applications. We use the official Tensorflow implementation along with most of the default hyperparameter settings. We set the target acceptance rate to 0.6, and use a max tree depth of 6 during the annealed sampling. We also continue to run the sampler after the annealing phase is finished, using a max tree depth of 10. We use a total of $1000$ intermediate distributions with $100$ parallel chains.

Finally, recall from the main text that each noise distribution in our experiments can be expressed as invertible transformation $F$ of a standard normal distribution. We use this $F$ to further enhance the efficiency of the NUTS sampler, by performing the sampling in the ${\mathbf{z}}$-space, and then mapping the final results back to ${\mathbf{x}}$-space. Working in ${\mathbf{z}}$-space, by the rules of transformations of random variables, the intermediate distributions of (45) become

AIS and RAISE can still be applied, just as before, to obtain an estimate of the log-likelihood in ${\mathbf{z}}$-space. The change of variables formula for probability density functions can then be applied to obtain estimated log-likelihoods for the original TRE model in ${\mathbf{x}}$-space. We note that when the noise distribution is a normalising flow, prior work has demonstrated that ${\mathbf{z}}$-space MCMC sampling can be significantly more effective than working in the original data space .

G.3 Experimental settings

We use the standard version of the MNIST dataset , with $50,000$ training points, and $10,000$ validation and test points. We follow the same preprocessing steps as , ‘dequantizing’ the dataset with uniform noise, re-scaling to the unit interval, and then mapping to the real line via a logit transformation.

The architecture for the TRE bridges is given in Figure 6(b). The waymark mechanism and associated grid-search is given in Table 2. A consistent observation across all our MNIST experiments was that the first ratio-estimator between the data distribution $p_{0}$ and a slightly perturbed data distribution $p_{1}$ was extremely prone to overfitting. We found that the only way to mitigate this problem was to simply drop the ratio by setting the $\alpha_{0}$ in (12) to a very small value (0.01) rather than exactly 0. Equivalently, this can be viewed as applying standard TRE to a very slightly perturbed data distribution. We note that this perturbation is small enough that is barely visible to the human eye when comparing samples. We conjecture that this problem may stem from the fact that the original MNIST dataset is actually discrete not continuous and the ‘dequantizing’ perturbation used to make the data continuous is perhaps not sufficient.

To form mini-batches, we sample $25$ datapoints each from $p_{0}$ and $p_{m}$, and then generate waymark trajectories as described in the main text. Thus, the total batch size is $25\times(m+1)$. We use the optimisation settings described at the beginning of this section, training for $200,000$ iterations, which takes about a day.

G.4 Additional results

In Figure 10, we present a sensitivity analysis showing how the quality of the learned EBM varies as we alter the number of waymarks, as well as the space in which the waymarks are generated. We found that working in ${\mathbf{x}}$-space yielded lower performance compared to working in ${\mathbf{z}}$-space, as measured by the most conservative estimator, RAISE. In particular, we found that the ${\mathbf{x}}$-space mechanism required more waymarks (around 15) to avoid any of the logistic losses saturating close to 0, and it was significantly harder to tune the spacing of the waymarks as indicated by Table 2.

Finally, for the models whose results were given in the main paper, we display extended image samples in Figure 11. Note that these samples are ordered by log-density (lowest density in top left corner, highest in bottom right).