An Easy to Interpret Diagnostic for Approximate Inference: Symmetric Divergence Over Simulations

1 Introduction

This paper considers the probabilistic inference problem. Given a known distribution $p (z, x)$ and observing some specific value $x$ , one wishes to infer $z$

. (E.g. predict its mean or variance.) Unless

p

is simple, there is no simple form for the posterior

p (z | x)

. Many approximate methods exist, including variants of MCMC, message-passing, Laplace’s method, and variational inference (VI). All these methods produce large errors on some problems. Thus, diagnostic techniques are of high interest to understand when a given inference method will perform well on a given problem.

For MCMC, there are several widely-used diagnostics. The potential scale reduction factor $^R$ diagnostic [Gelman_1992_InferenceIterativeSimulation] runs multiple chains, and then compares within-chain and between-chain variances. The expected sample size diagnostic considers correlations in a single chain. Diagnostics of this type are an active research area [Vehtari_2020_Ranknormalizationfoldinglocalization].

While successful, these diagnostics are grounded in the fact that MCMC is asymptotically exact. That is, under mild conditions MCMC will converge to the stationary distribution if run long enough. Informally, diagnostics for MCMC only need to diagnose “has the chain converged?”, rather than “has it converged to the correct distribution?”.

For inference methods that are asymptotically approximate, different diagnostics are needed. This paper is in the line of simulation-based diagnostics. These are a fairly radical departure. Rather than measuring how well inference performs on the given $x$ , these estimate how well inference performs on average over data generated from the model. These diagnostics repeatedly sample $(z, x) \sim p (z, x)$ and then do inference on the simulated $x$ . The power of this approach is that the true latent $z$ corresponding to the observed $x$ is known.

To the best of our knowledge, this simulation-based approach was first pursued by Cook_2006_ValidationSoftwareBayesian, who sample $x \sim p (x)$ and then perform inference to approximately sample $z \sim p (z | x) .$

The quantiles of each component

z_{i}

generated this way are compared to those generated directly from the prior

p (z) .

This can be done visually (looking at histograms), or by using a Kolmogorov-Smirnov test. More recently, Yao_2018_YesDidIt suggest testing for symmetry. These error measures may not be appropriate for all situations. First, these measures can be challenging to automate or interpret, since they do not provide a scalar quantity but rather a procedure to perform in each dimension. Second, there could be inference errors not detected by looking at univariate distributions.

In this paper, we observe that some inference methods, such as Laplace’s method and variational inference (VI), do not simply give a set of samples, but an approximate distribution $q (z | x)$ . This turns out to enable diagnostics that would be impossible with MCMC.

Our central idea is simple. Suppose that on input $x$ , inference returns a distribution $q (z | x) .$

Define the joint distribution

q (z, x) = p (x) q (z | x) .

Then, our diagnostic is an estimate of

SKL (p (z, x) ∥ q (z, x)),

the symmetric KL-divergence between

p (z, x)

and

q (z, x)

. This essentially measures how far

q (z | x)

is from

p (z | x)

on average, over simulated datasets.

The key observation is that the symmetric divergence induces cancellations induces cancellations between among the unknown normalization terms. Specifically, if $(z, x) \sim p (z, x)$ and $~ z \sim q (z | x)$ then we can simulate

d = log \frac{p (z, x)}{q (z | x)} - log \frac{p (~ z, x)}{q (~ z | x)},

and the expected value of $d$ is the symmetric divergence. To compute this diagnostic, one must: (1) simulate $(z, x) \sim p (z, x)$ and $~ z \sim q (z | x)$ and (2) compute $p (z, x)$ and $q (z | x)$ . We do not need to be able to evaluate $p (x),$ even though it is part of the definitions of $p (z, x)$ and $q (z, x)$ .

We also show that this idea can be extended to situations with conditional or hidden variables. As an example of the latter, we show that it can be used with importance-weighted inference methods that generate many samples $z \sim q (z | x)$ and select then one according to the importance weights $p (z, x) / q (z | x)$ [Burda_2015_ImportanceWeightedAutoencoders]. Experiments show that the diagnostic gives practical measures of performance, for regular VI, for Laplace’s method, and for importance-weighted variants of both.

1.1 Notation

The KL-divergence is $KL(q(z)∥p(z))=\operatornamewithlimitsEq(z)log(q(z)/p(z))$

. Sans-serif font marks random variables. This disambiguates conflicting conventions in machine learning and information theory.

KL(q(z|x)∥p(z|x))=\operatornamewithlimitsEq(z|x)log(q(z|x)/p(z|x))

is a divergence over

z

for a fixed

x

. Meanwhile,

KL(q(z|x)∥p(z|x))=\operatornamewithlimitsEq(z,x)log(q(z|x)/p(z|x))

is the conditional divergence [Cover_2006_Elementsinformationtheory], with an expectation over both

z

and

x

. In all cases, symmetric divergences are defined as

SKL (q ∥ p) = KL (q ∥ p) + KL (p ∥ q)

2 A New Simulation-Based Diagnostic

This section gives a novel simulation-based diagnostic based on the symmetric KL-divergence. The key idea is that some inference methods (e.g. VI or Laplace’s method) do not just give approximate samples, but an approximate distribution that can be evaluated at any point. This enables certain diagnostics that would be impossible with just a set of samples. Again, let $p (z, x)$ be the target. We consider approximate inference methods that input some $x$ and produce a distribution over $z$ that approximates $p (z | x) .$ We denote that approximation as $q (z | x) .$

One might hope to use the KL divergence $KL (q (z | x) ∥ p (z | x))$ as a diagnostic. This is almost never tractable since $p (x)$ is unknown. One can instead compute the evidence lower bound (ELBO) $\operatornamewithlimitsEq(z|x)log(q(z|x)/p(z,x))$ which is equal to the KL-divergence plus $log p (x) .$ The ELBO precisely measures the relative error for different algorithms, but gives little information about the absolute error, since $p (x)$ is unknown.

Instead, our diagnostic is based on the symmetric KL-divergence. The basic idea of the diagnostic is to define a joint distribution $q (z, x) = p (x) q (z | x)$ with the same distribution over $x$ as $p (z | x) .$ Then, cancellations make it possible to estimate the joint symmetric divergence between $p (z, x)$ and $q (z, x)$ . This is formalized in the following result.

Theorem 1.

Given $p (z, x)$ and $q (z | x),$ define $q (z, x) = p (x) q (z | x)$ . Then,

SKL (q (z, x) ∥ p (z, x))

=\operatornamewithlimitsE[logp(z,x)q(z|x)−logp(~z,x)q(~z|x)],

where $(z, x) \sim p (z, x)$ is sampled from the model distribution and $~ z \sim q (z | x)$ is sampled from the approximating distribution.

For

k = 1, 2, \dots, K

Simulate

(z, x) \sim p (z, x)

. Infer

p (z, x)

to get

q (z | x)

. (fix

x

) Simulate

~ z \sim q (z | x)

. (fix

x

)

d_{k} \leftarrow log \frac{p (z, x)}{q (z | x)} - log \frac{p (~ z, x)}{q (~ z | x)}

Use

ave (d_{1} \dots d_{K}) \approx SKL (q (z, x) ∥ p (z, x)) .

Algorithm 1 Computing the proposed diagnostic. Input

x

. For

k = 1, 2, \dots, K

Simulate

(z, y) \sim p (z, y | x)

. Infer

p (z, y | x)

to get

q (z | y, x)

. (fix

x, y

) Simulate

~ z \sim q (z | y, x)

. (

x, y

fixed)

d_{k} \leftarrow log \frac{p (z, y | x)}{q (z | y, x)} - log \frac{p (~ z, y | x)}{q (~ z | y, x)}

Use

ave (d_{1} \dots d_{K}) \approx SKL (q (z, y | x) ∥ p (z, y | x)) .

Algorithm 2 Diagnostic for conditional models.

Pseudocode for how one would use this result is given in Alg. 1. One attractive aspect is that the output is the mean of a set of $K$

independent quantities. This makes it is easy to produce uncertainty measures such as confidence intervals. These bound how far the estimated diagnostic may be from the true symmetric divergence.

2.1 Inference in Conditional Models

Many inference problems are conditional, meaning one is given a model $p (z, y | x)$ , with no distribution specified over $x$ . After observing $x$ and $y$ , the goal is to predict $z$ . For example, in regression or classification problems, $x$ would represent the input features, $y$ the output values/labels, and $z$ the latent parameters.

In these cases, inference takes as input a pair $(y, x)$ and produces a distribution $q (z | x, y)$ approximating $p (z | x, y)$ . It’s easy to see that the following generalization of Thm. 1 holds. This is given by taking Thm. 1, substituting $y$ for $x$ and then conditioning all distributions on the fixed value $x$ .

Corollary 2.

SKL (q (z, y | x) ∥ p (z, y | x))

=\operatornamewithlimitsElogp(z,y|x)q(z|y,x)−logp(~z,y|x)q(~z|y,x)

where $(z, y) \sim p (z, y | x)$ is sampled from the model distribution and $~ z \sim q (z | y, x)$ is sampled from the approximating distribution.

Pseudocode for how the diagnostic would be used with conditional models is given as Alg. 2. It’s critical that $x$ is not a random variable– it is the actual observed input data. The simulated latent variables $z$ and datasets $y$ are conditioned on $x$ .

In order to run use this algorithm, one must be able to perform the following operations: (1) Simulate $(z, y) \sim p (z, y | x)$ and $~ z \sim q (z | y, x)$ . (2) Compute $p (z, y | x)$ for a given $x$ , $y$ , and $z$ . (3) Compute $q (z | x)$ for a given $x$ , $y$ , and $z$ . It is not necessary to be able to evaluate $p (y | x),$ despite the fact that it is part of the definition of $q (z, x)$ .

Example Results. Before moving on to more complex cases, we give some examples of the use of this diagnostic. Fig. 1 shows an example of running the diagnostic on five example models using variational inference (VI) and Laplace’s method that maximizes $log p (z, x)$ to get $^z$ and uses a Gaussian centered at $^z$ with a covariance matching the Hessian $H$ of $log p$ . We also compare to an “adjusted” Laplace’s method that better matches the curvature if $^z$ is only an approximate maxima. This instead uses a mean of $H^{- 1} g$ where $g$ is that gradient of $log p$ at $^z .$ A more full description of the models and inference algorithms is given in Sec. 5.

Figure 1: The diagnostic gives plausible measures of the value of different inference algorithms and of optimizing for different numbers of iterations. The diagnostic computed with $K = 100$ repetitions. Lines show the mean while the colored areas show 95% confidence intervals. Confidence intervals are computed before the log-transform and therefore appear large for the lower-bounds seen on concrete. Laplace’s method fails due to numerical problems with few iterations on hospitals. Adjusted Laplace’s method is exact for concrete.

3 Inference with Augmented Variables

Many approximate inference methods used the idea of augmentation. The idea is to create an extra variable $h$ and then approximate $p (z, h | x)$ with $q (z, h | x)$ . Why would this be useful? The basic reason is that many powerful approximating families are obtained by integrating out other random variables. Such families often do not have tractable densities $q (z | x)$ , but can be represented as the marginal of some density $q (z, h | x) .$ If we choose $p (h | z, x)$ in a way that is “easy” for $q$ to match, then augmented inference might be nearly as accurate as directly approximating $p (z | x)$ with $q (z | x) .$ Agakov_2004_AuxiliaryVariationalMethod introduced the idea of auxiliary variational inference, which fits this form.

To apply the diagnostic to inference with hidden variables, we need another version of Thm. 1. This can be proven by taking Thm. 1 and substituting $(z, h)$ for $z$ .

Corollary 3.

Given $p (z, h, x)$ and $q (z, h | x),$ define $q (z, h, x) = p (x) q (z, h | x) .$ Then

SKL (q (z, h, x) ∥ p (z, h, x))

=\operatornamewithlimitsElogp(z,h,x)q(z,h|x)−logp(~z,~h,x)q(~z,~h|x),

where $(z, x) \sim p (z, x)$ is sampled from the model distribution, $h \sim p (h | z, x)$ is sampled from the augmenting distribution, and $(~ z, ~ h) \sim q (z, h | x)$ is sampled from the approximating distribution.

The diagnostic is useful because (by the chain rule of KL-divergence),

SKL (q (z, h, x) ∥ p (z, h, x)) = SKL (q (z, x) ∥ p (z, x)) + SKL (q (h | z, x) ∥ p (h | z, x)),

i.e. the diagnostic yields an upper bound on the error of $q (z, x) .$

The corresponding algorithm is given as Alg. 3. This shows an important computational constraint. We assume that $p (z, x)$ is given as a “black box”. Thus, the augmented $p (z, h, x)$ must be chosen so that it is tractable to simulate $p (h | z, x) .$ The next section will give a special case of this result for a certain class of approximate augmented inference methods.

For $k = 1, 2, \dots, K$ :

Simulate $(z, x) \sim p (z, x)$ .
Infer $p (z, x) p (h | z, x)$ to get $q (z, h | x)$ . (fix $x$ )
Simulate $h \sim p (h | z, x) .$
Simulate $(~ z, ~ h) \sim q (z, h | x)$ . (fix $x$ )
$d_{k} \leftarrow log \frac{p (z, h, x)}{q (z, h | x)} - log \frac{p (~ z, ~ h, x)}{q (~ z, ~ h | x)}$

Use $ave (d_{1} \dots d_{K}) \approx SKL (q (z, h, x) ∥ p (z, h, x)) .$

Algorithm 3 Diagnostic with augmentation.

For $k = 1, 2, \dots, K$ :

Simulate $(z_{1}, x) \sim p (z, x)$ .
Run inference to find a base distribution $q (z | x)$ .
Simulate $z_{2}, \dots, z_{M} \sim q (z | x) .$
Simulate ${~ z}_{1}, \dots, {~ z}_{M} \sim q (z | x)$ . (fix $x$ )
$d_{k} \leftarrow log M \sum m = 1 \frac{p (z_{m}, x)}{q (z_{m} | x)}$ $- log M \sum m = 1 \frac{p ({~ z}_{m}, x)}{q ({~ z}_{m} | x)}$

Use $ave (d_{1} \dots d_{K}) \approx$

$SKL (q_{I W} (z_{1} \dots z_{M}, x) ∥ p_{I W} (z_{1} \dots z_{M}, x)) .$

Algorithm 4 With importance-weighting.

4 Importance Sampling

Self-normalized importance sampling is a classic Monte-Carlo method [Owen_2013_MonteCarlotheory]. Given any distribution $q (z | x)$ , one can approximately sample from the posterior $p (z | x)$ by drawing a set of $M$ samples ${^z}_{1}, \dots, {^z}_{M} \sim q (z | x)$

, selecting one with probability proportional to the importance weights

p ({^z}_{m}, x) / q ({^z}_{m})

and then returning the final sample

z_{1} = {^z}_{m} .

Call the resulting density

q_{I W} (z_{1} | x) .

One might hope to directly apply the diagnostic to self-normalizing importance sampling. However, this cannot be done because it is intractable to to evaluate $q_{I W} (z_{1} | x) .$ However, we can identify augmented distributions, and thereby upper-bound the symmetric divergence. Define the distribution

p_{I W} (z_{1}, \dots, z_{M}, x) = p (z_{1}, x) M \prod m = 2 q (z_{m} | x) .

(1)

This can be seen as an augmented distribution with $z_{2}, \dots, z_{M}$ the hidden variables augmenting the original $z_{1} .$ Define also

q_{I W} (z_{1}, \dots, z_{M} | x) = \frac{p_{I W} (z_{1}, \dots, z_{M}, x)}{\frac{1}{M} \sum_{m = 1}^{M} \frac{p (z_{m}, x)}{q (z_{m} | x)}} .

(2)

It is not immediately obvious that this augments the self-normalized importance sampling density $q_{I W} (z_{1} | x)$ introduced at the beginning of this section (or indeed that this is a valid density at all). However, it was recently shown [Domke_2018_ImportanceWeightingVariational] that the following algorithm samples from $q_{I W} .$

Claim 4.

The following process yields a sample from $q_{I W} (z_{1}, \dots, z_{M} | x)$ as defined in Eq. 2.

Draw ${^z}_{1}, \dots {^z}_{M} \sim q (z | x)$ .
Choose $m \in {1, \dots, M}$ with $\operatornamewithlimitsP[m]∝p(^zm,x)q(^zm)$
Set $(z_{1}, \dots, z_{M}) = ({^z}_{m}, {^z}_{1}, \dots, {^z}_{m - 1}, {^z}_{m + 1}, \dots, {^z}_{M})$

Informally, this algorithm draws $m$ samples from $q (z | x)$ then swaps one to be in the first position, chosen according to importance weights. Thus, this is a valid augmentation of the self-normalized importance sampling distribution.

This shows that Eq. 1 and Eq. 2 augment the target and self-normalized importance-sampling density, and thus that $SKL (q_{I W} (z_{1}, \dots, z_{M}, x) ∥ p_{I W} (z_{1}, \dots, z_{M}, x))$ upper-bounds $SKL (q_{I W} (z_{1}, x) ∥ p (z_{1}, x)) .$ It can be shown that the (non-symmetric) divergence from $q_{I W}$ to $p_{I W}$ asymptotically decreases at a $1 / M$ rate [Maddison_2017_FilteringVariationalObjectives, Domke_2018_ImportanceWeightingVariational].

Figure 2: The diagnostic applied to self-normalized importance-weighting with $M$ samples. The proposal distribution is estimated by Laplace’s method with adjustment. The $M = 1$ curve is the same as the “laplace (adjusted)” curve in Fig. 1.

4.1 Diagnosing Importance Weighted Inference

The above discussion shows that adding importance-sampling can improve any distribution $q .$ Even better results can be obtained by explicitly optimizing $q$ to work well after augmentation. Importance-weighted variational inference directly performs an optimization to maximize the ELBO from $q_{I W}$ to $p_{I W},$ equivalent to minimizing the KL-divergence $KL (q_{I W} ∥ p_{I W})$ between the augmented distributions. This ELBO can be simplified using the relationship that

\operatornamewithlimitsEqIW(z1,⋯,zM|x)logpIW(z1,⋯,zM,x)qIW(z1,⋯,zM|x)=\operatornamewithlimitsEq(z1|x)⋯\operatornamewithlimitsEq(zM|x)log1MM∑m=1logp(zm,x)q(zm|x),

(3)

which follows from cancellations between $p_{I W}$ and $q_{I W}$ , followed from the observation that the argument of the expectation is constant with respect to permutations of $z_{1}, \dots, z_{m} .$ This objective was originally introduced in the context of importance-weighted auto-encoders [Burda_2015_ImportanceWeightedAutoencoders, Cremer_2017_ReinterpretingImportanceWeightedAutoencoders, Domke_2018_ImportanceWeightingVariational, Naesseth_2018_VariationalSequentialMonte, Le_2018_AutoEncodingSequentialMonte] (without explicitly identifying $q_{I W}$ and $p_{I W}$ ) and subsequently studied by various others [Cremer_2017_ReinterpretingImportanceWeightedAutoencoders, Domke_2018_ImportanceWeightingVariational, Bachman_2015_TrainingDeepGenerative, Naesseth_2018_VariationalSequentialMonte, Le_2018_AutoEncodingSequentialMonte].

The following result specializes Cor. 3 to the case of importance-weighted inference.

Corollary 5.

Given $p (z, x)$ and $q (z | x),$ define $p_{I W}$ and $q_{I W}$ as in Eq. 1 and Eq. 2. Further, set $q_{I W} (z_{1}, \dots, z_{M}, x) = p (x) q_{I W} (z_{1}, \dots, z_{M} | x) .$ Then

SKL (q_{I W} (z_{1}, \dots, z_{M}, x) ∥ p_{I W} (z_{1}, \dots, x))

=\operatornamewithlimitsElog1MM∑m=1p(zm,x)q(zm|x)−log1MM∑m=1p(~zm,x)q(~zm|x),

where $(z_{1}, x) \sim p (z, x)$ is sampled from the model distribution, $z_{2}, \dots, z_{M} \sim q (z | x)$ is sampled from the approximating distribution and ${~ z}_{1}, \dots, {~ z}_{M} \sim q (z | x)$ are also sampled from the approximating distribution.

A proof is in the supplement. Unlike Cor. 2 and Cor. 3, the result is not trivial. The main idea is to substitute $p_{I W} (z_{1}, \dots, z_{M}, x)$ for $p (z, h, x)$ and $q_{I W} (z_{1}, \dots, z_{M} | x)$ for $q (z, h | x) .$ Then, many expressions can be simplified based on the particular forms of $p_{I W}$ and $q_{I W}$ . Finally, we can observe that the argument of the expectation is independent of permutations of ${~ z}_{1}, \dots, {~ z}_{M}$ . This allows a final simplification.

5 Experiments

Figure 3: The diagnostic as computed using Alg. 4 where parameters are determined by optimizing Eq. 3 with the same value $M$ . The $M = 1$ curve is the same as the “vi” curve in Fig. 1.

Models. We use five models, described in detail in Sec. Sec. 7 (Supplement). glm_binomial is a hierarchical model of the number of a bird population over time. heart_transplants models the survival times of patients after surgery. hospitals measures the number of deaths in different hospitals. ionosphere

is a Bayesian logistic regression model on a classic dataset.

concrete

is a Bayesian linear regression model – included as a baseline because the exact posterior is Gaussian.

Optimization. The first inference algorithm we consider is Laplace’s method which produces a multivariate Gaussian approximation $q (z | x)$ . For this method, we run Adam with a step-size of 0.01 for the first half of iterations, and 0.001 for the second half. The resulting point $^z$ is the mean. Finally, the Hessian of $log p$ at $^z$ is used to estimate the covariance of $q$ as $Σ = - H^{- 1} .$ This method fails to match the local curvature of $log p$ when $^z$ is far from an optima. We also consider an “adjusted” Laplace’s method that instead uses a mean of $H^{- 1} g$ where $g$ is that gradient of $log p$ at $^z$ . This guarantees that $log q$ has the same gradient at $^z$ as $log p$ .

We also consider variational inference. We initialize to a standard Gaussian and optimize with Adam with a step size of 0.001 for the first half of iterations, and 0.0001 for the second half. We estimate the gradient of the ELBO using the reparameterization trick, using the “sticking the landing” estimator [Roeder_2017_StickingLandingSimple].

Constraints and Transformations. Often, random variables have constraints – they are not supported over the reals. As is common [Kucukelbir_2017_AutomaticDifferentiationVariational] we deal with this through a process of transformations. For our models, it is sufficient to consider two cases:

Random variables $x$ that are either defined over the non-negative reals $[0, \infty)$ . In this case, we replace $x$ with a new random variable $x^{'} = log x$ that is unconstrained.
Random variables $x$ defined on a closed interval $[a, b] .$ Here, we transform to $x^{'} = L o g i t (\frac{x - a}{a - b}) .$

Results. Results comparing VI and the two variants of Laplace’s method are shown in Fig. 1, averaging over $K = 100$ simulated datasets. Laplace’s method is reasonably accurate in many cases, but usually has a “floor” of accuracy it does not exceed. The adjustment to Laplace’s method is often helpful and never harmful. VI performs better with many iterations. For these models, the diagnostic shows that inference error is reasonably low with many iterations, but not quite “exact”.

Fig. 2 shows the results of importance sampling with a proposal computed using Laplace’s method with adjustment. Using more samples yields a clear improvement.

Finally, results with full importance weighted variational inference (optimizing the Eq. 3 rather than the standard ELBO) is shown in Fig. 3. The same value $M$ is used during optimization and at test time.

6 Discussion

This paper proposed a new diagnostic for approximate inference methods. This is a simulation-based diagnostic, meaning it is computed by repeatedly simulating latent variables along with datasets and running inference on each dataset. The central idea is that cancellations in unknown constants make it possible to estimate a symmetric divergence. This is notable in being a simple, scalar quantity with a clear information-theoretic interpretation. It can also be computed in a fully automated way along with error measures like confidence intervals. We showed that the diagnostic can be extended to augmented inference methods, in particular importance-weighted inference. Empirically, the method gives reasonable diagnostic information on several test models.

While Thm. 1 is quite simple, there are numerous points worth clarifying in its use as a diagnostic:

What the diagnostic measures. One possibly counter-intuitive aspect of this diagnostic (like all simulation-based diagnostics) is that it does not use the actual observed data $x$ . Rather, it measures the typical error, averaged over $x$ simulated from the model $p (x)$ . It is this essential that the prior $p (z)$ and likelihood $p (x | z)$ be selected so that $p (x) = \int p (z) p (x | z) d z$ yields realistic simulated datasets. In particular, very broad priors are might lead to “nonsense” observations $x$ that are unrepresentative of the data that would be seen in practice.

Computational considerations. In order to compute this diagnostic, one must be able to perform several operations: (1) Simulate $(z, x) \sim p (z, x)$ and $~ z \sim q (z | x)$ . (2) Compute $p (z, x)$ for a given $z$ and $x$ . (3) Compute $q (z | x)$ for a given $z$ and $x$ . Crucially, it is not necessary to be able to evaluate $p (x)$ . This is true despite the fact that $p (x)$ is part of the definition of $q (z, x)$ .

Other representations of the diagnostic. The quantity that the diagnostic is representing can be written in a different form that emphasizes that it measures errors over $z$ . This uses the notion of a conditional divergence $KL (q (z | x) ∥ p (z | x)) .$ It is not hard to show that in the setting of Thm. 1 that

SKL (q (z, x) ∥ p (z, x)) = SKL (q (z | x) ∥ p (z | x)) .

This is true because of cancellations between $q$ and $p$ due to the inclusion of the $p (x)$ term in $q (x, z) .$

Use with randomized inference methods. In practice, approximate inference algorithms are often non-deterministic. This is not reflected by the notation $q (z | x)$ . With non-determinism, no modification to the diagnostic technique in Alg. 1 is needed. Only the interpretation is slightly different. To formalize what the diagnostic measures in this case, define $q (z | x, ω)$ to be the approximate posterior produced where $ω$ are the random numbers underlying the algorithm. Then, a diagnostic can be defined as the expected divergence between $q$ and $p$ , i.e., $\operatornamewithlimitsEωSKL(q(z,x|ω)∥p(z,x)),$ where $q (z, x | ω) = p (x) q (z | x, ω) .$ This is still a very reasonable measure of the accuracy of inference. For simplicity, our presentation mostly neglects the issue of non-determinism in approximating distributions.

6.1 Related Work

There are several lines of related work not mentioned so far. One recent line of work explores inference diagnostics based on Stein’s method [Gorham_2017_MeasuringSampleQuality, Gorham_2015_MeasuringSampleQuality]. The idea is to create sets of functions whose true expectations must be zero. Any deviation from zero in those functions indicates inference failure. This diagnostic can be used with methods that are asymptotically approximate. However, it is intended for cases when error decreases to zero. There is no claim that the magnitude of the diagnostic is a good measure of the usefulness of an approximate posterior.

The test proposed by Geweke_2004_GettingItRight is an interesting early diagnostic that repeatedly simulates datasets in a non-independent manner. The idea is to iteratively sample $x \sim p (x | z)$ , the run inference to produce $q (z | x)$ and then sample $z \sim q (z | x) .$ Then, one compares the expectation of some function $g (z, x)$ to those on exact samples $(z, x) \sim p (z, x) .$ If $q$ is exact, these expectations should match. We prefer an approach where each simulation is independent since this is easier to parallelize, avoids correlations between simulations, and makes it easier to compute error measures like confidence intervals.

Bidirectional MCMC [Grosse_2015_Sandwichingmarginallikelihood] runs MCMC on repeated simulated datasets to get upper and lower bounds on the marginal likelihood $log p (x) .$ This is intended as a technique to evaluate the quality of a model, not as a diagnostic for inference. Still, in principle one could use these to transform an ELBO into bounds on the KL-divergence. One drawback is the expense of repeatedly running MCMC. Typically, variational inference is used in settings where MCMC would be too expensive.

6.2 Limitations and Future Work

This work has several limitations shared with all simulation-based diagnostics: First, computing them requires repeating inference numerous times. This comes with an associated cost. Second, these methods can be overly pessimistic when used with extremely broad or uninformative priors. It is important that the model is chosen so that simulated data are representative of the datasets one cares about. Third, the diagnostic measures average accuracy over data simulated the prior, as opposed to the expected accuracy for a particular dataset. (Put another way, the diagnostic is arguably frequentist rather than Bayesian.)

One might be concerned about the success of this diagnostic when used with variational inference methods. Namely, VI typically minimizes $KL (q ∥ p)$ while the diagnostic is based on the symmetric divergence. Informally, VI cares about finding a distribution that is close in a “mode finding” divergence, while the diagnostic measures both “mode finding” and “mode spanning”. It is possible that a distribution could be close in VI’s objective, yet yield a high diagnostic value. This is arguably a flaw not of the diagnostic, but of variational inference. One interesting future direction would be to investigate recent VI variants that try to minimize other divergences [Li__RenyiDivergenceVariational, Dieng_2017_VariationalInferencechi].

In future work, it would be interesting to address MCMC methods. Of course, most MCMC methods are not suitable for this framework. However, some methods like annealed importance sampling [Neal_1998_AnnealedImportanceSampling] formally create augmented target and proposal densities at a variety of “temperatures”. It may be possible to use the diagnostic proposed here to measure the symmetric divergence between these augmented distributions. This could potentially offer a diagnostic for MCMC with the unusual property that the diagnostic going to zero is both necessary and sufficient to guarantee convergence to the stationary distribution.

References

7 Models

GLM Binomial. This is a model of the number of peregrine pairs $c_{i}$ in the French Jura in year $x_{i} .$ The data is from between $1964$ and $2003,$ but $x_{i}$ is scaled to between $- 1$ and $+ 1$ . The model and data are from Kery_2012_Bayesianpopulationanalysis.

	$\upalpha$	$\sim N (0, (10)^{2})$
	$\upbeta1$	$\sim N (0, (10)^{2})$
	$\upbeta2$	$\sim N (0, (10)^{2})$
	$c_{i}$	$\sim B i n o m i a l (n_{i}, α + β_{1} x_{i} + β_{2} x_{i}^{2})$

Heart Transplants. This is a model of a hypothetical population of patients who underwent a surgery, of whom $y_{T} = 8$ survived [Lunn_2013_BUGSbookpractical, Ex. 3.5.1]. These were tracked to see the number of years $s_{i}$ ( $i \in {1, \dots, 8}$ ) that the $i$

-th patient who survived surgery lived post-surgery. This is assumed to be determined by an exponential distribution with parameter

θ

. Thus, the model is:

	$p_{T}$	$\sim U n i f o r m (0, 1)$
	$y_{T}$	$\sim B i n o m i a l (N, p_{T})$
	$\uptheta$	$\sim G a m m a (1 / 3, 1 / 3)$
	$s_{i}$	$\sim E x p o n e n t i a l (θ)$

Note that, when generating synthetic datasets for this model, we always use the same set of variables $s_{1}, \dots, s_{8}$ , independent of the value of $y_{T}$ . This is done because of the difficulties posed by having different dimensionality in different realizations of the posterior. While not fully in keeping with the spirit of the original model, this still defines a perfectly valid probabilistic model and test of the diagnostic.

Hospitals. This is a hierarchical model of the mortality rate of English hospitals performing heart surgery [Lunn_2013_BUGSbookpractical, Ex. 10.1.1]. The data is ${(n_{i}, y_{i})}$ where $n_{i}$ is the number of operations in hospital $i$ and $y_{i}$

is the corresponding number of deaths. The logit of the true mortality rate

θ_{i}

of hospital

i

is a Gaussian with unknown mean

μ

and standard deviation

ω

. The latent variables are

ω

μ

, and

{θ_{i}}

	$ω$	$\sim U n i f o r m (.25, 1)$
	$μ$	$\sim U n i f o r m (- 3, 3)$
	$L o g i t (θ_{i})$	$\sim N (μ, ω^{2})$
	$y_{i}$	$\sim B i n o m i a l (n_{i}, θ_{i})$

The original model has a very wide prior on $ω$ and $μ$ , which leads to the problems discussed in LABEL:subsec:new-diagnostic-discussion. We use the above model with more modest priors.

Ionosphere. This is a classic dataset for binary classification. We model it as a Bayesian logistic regression problem with a standard Gaussian prior over the weights $w .$

Concrete. This is a well-known dataset for linear regression. We model it as a Bayesian linear regression problem with a standard Gaussian prior over the weights $w .$ This model is particularly notable because the true posterior is exactly Gaussian. Since both Laplace’s method and variational inference can exactly represent such a posterior, this provides an important test if the diagnostic can correctly recognize inference success when it occurs.

8 Theory

See Thm. 1

Proof.

The divergence $SKL (q (z, x) ∥ p (z, x))$ is equal to

	$\operatornamewithlimitsEp(z,x)logp(z,x)q(z,x)+\operatornamewithlimitsEq(z,x)logq(z,x)p(z,x)$	$=\operatornamewithlimitsEp(x)[\operatornamewithlimitsEp(z\|x)logp(z,x)q(z,x)+\operatornamewithlimitsEq(z\|x)logq(z,x)p(z,x)]$
		$=\operatornamewithlimitsEp(x)[\operatornamewithlimitsEp(z\|x)logp(z,x)q(z\|x)+\operatornamewithlimitsEq(z\|x)logq(z\|x)p(z,x)]$

In the first line we use the fact that $q (x) = p (x),$ while in the second line we pull out a factor of $log q (x)$ from each term, which cancel. The claimed result is the same as the last line with a sign change. ∎

See Cor. 5

Proof.

Start with the result of Cor. 3.

	$SKL (q (z, h, x) ∥ p (z, h, x))$	$=\operatornamewithlimitsElogp(z,h,x)q(z,h\|x)−logp(~z,~h,x)q(~z,~h\|x)$
	$(z, x)$	$\sim p (z, x)$
	$h$	$\sim p (h \| z, x)$
	$(~ z, ~ h)$	$\sim q (z, h \| x) .$

Now, make the following transformations

$q$	$\Rightarrow$	$q_{I W}$
$p$	$\Rightarrow$	$p_{I W}$
$z$	$\Rightarrow$	$z_{1}$
$h$	$\Rightarrow$	$(z_{2}, \dots, z_{M}) .$

Then, we get

	$SKL (q_{I W} (z_{1} \dots, z_{M}, x) ∥ p (z_{1}, \dots, z_{M}, x))$	$=\operatornamewithlimitsElogpIW(z1⋯,zM,x)qIW(z1⋯,zM\|x)−logpIW(~z1⋯,~zM,x)qIW(~z1⋯,~zM\|x)$
	$(z_{1}, x)$	$\sim p_{I W} (z_{1}, x)$
	$z_{2} \dots, z_{M}$	$\sim p_{I W} (z_{2} \dots, z_{M} \| z, x)$
	$({~ z}_{1}, \dots, {~ z}_{M})$	$\sim q_{I W} (z_{1} \dots, z_{M} \| x) .$

Now, note that

$\frac{p_{I W} (z_{1} \dots, z_{M}, x)}{q_{I W} (z_{1} \dots, z_{M} \| x)}$	$=$	$\frac{1}{M} M \sum m = 1 \frac{p (z_{m}, x)}{q (z_{m} \| x)}$
$(z_{1}, x)$	$\sim$	$p_{I W} (z_{1}, x)$
	$=$	$p (z_{1}, x)$
$z_{2} \dots, z_{M}$	$\sim$	$p_{I W} (z_{2} \dots, z_{M} \| z, x)$
	$=$	$M \prod m = 1 q (z_{m} \| x)$

This leaves us with the result of

	$SKL (q_{I W} (z_{1}, \dots, z_{M}, x) ∥ p_{I W} (z_{1}, \dots, x))$
	$=\operatornamewithlimitsElog1MM∑m=1p(zm,x)q(zm\|x)−log1MM∑m=1p(~zm,x)q(~zm\|x),$
	$=\operatornamewithlimitsElogM∑m=1p(zm,x)q(zm\|x)−logM∑m=1p(~zm,x)q(~zm\|x),$

	$(z_{1}, x)$	$\sim p (z, x)$
	$z_{2}, \dots, z_{M}$	$\sim q (z \| x)$
	$({~ z}_{1}, \dots, {~ z}_{M})$	$\sim q_{I W} (z_{1} \dots, z_{M} \| x) .$

Now, finally, note that the expectation is unchanged under permutations of the order of ${~ z}_{1}, \dots, {~ z}_{M}$ . Thus, the expectation is unchanged if we replace the distribution with

{~ z}_{1}, \dots, {~ z}_{M} \sim q (z | x) .

∎

An Easy to Interpret Diagnostic for Approximate Inference: Symmetric Divergence Over Simulations

Authors

AIDE: An algorithm for measuring the accuracy of probabilistic inference algorithms

A Contrastive Divergence for Combining Variational Inference and MCMC

Divergence Triangle for Joint Training of Generator Model, Energy-based Model, and Inference Model

Diagnostics for Monte Carlo Algorithms for Models with Intractable Normalizing Functions

Alpha-Beta Divergence For Variational Inference

Recursive Monte Carlo and Variational Inference with Auxiliary Variables

Quantifying the probable approximation error of probabilistic inference programs

1 Introduction

1.1 Notation

2 A New Simulation-Based Diagnostic

Theorem 1.

2.1 Inference in Conditional Models

Corollary 2.

3 Inference with Augmented Variables

Corollary 3.

4 Importance Sampling

Claim 4.

4.1 Diagnosing Importance Weighted Inference

Corollary 5.

5 Experiments

6 Discussion

6.1 Related Work

6.2 Limitations and Future Work

References

7 Models

8 Theory

Proof.

Proof.

An Easy to Interpret Diagnostic for Approximate Inference: Symmetric Divergence Over Simulations

Authors

Related Research

AIDE: An algorithm for measuring the accuracy of probabilistic inference algorithms

A Contrastive Divergence for Combining Variational Inference and MCMC

Divergence Triangle for Joint Training of Generator Model, Energy-based Model, and Inference Model

Diagnostics for Monte Carlo Algorithms for Models with Intractable Normalizing Functions

Alpha-Beta Divergence For Variational Inference

Recursive Monte Carlo and Variational Inference with Auxiliary Variables

Quantifying the probable approximation error of probabilistic inference programs

This week in AI

1 Introduction

1.1 Notation

2 A New Simulation-Based Diagnostic

Theorem 1.

2.1 Inference in Conditional Models

Corollary 2.

3 Inference with Augmented Variables

Corollary 3.

4 Importance Sampling

Claim 4.

4.1 Diagnosing Importance Weighted Inference

Corollary 5.

5 Experiments

6 Discussion

6.1 Related Work

6.2 Limitations and Future Work

References

7 Models

8 Theory

Proof.

Proof.