Chapter 9
Variational inference

Lectures 17–18 · inference as optimisation, the ELBO, mean-field and full-rank families, and knowing when to trust the answer

fast answers that are a little too sure of themselves

Every posterior in this course so far has been computed by sampling. Hamiltonian Monte Carlo draws a cloud of points from the posterior and lets the cloud stand in for the distribution; given enough draws, it is exact. Its weakness is cost. Every iteration evaluates the gradient of the log posterior, which means touching every data point, and a few thousand iterations across four chains means tens of thousands of passes through the data. With a few thousand rows this is invisible. With hundreds of thousands of rows, or a model that must be refitted every night, or dozens of candidate models to compare, it becomes the bottleneck of the whole analysis.

Variational inference (VI) replaces sampling with optimisation. Instead of drawing from the posterior, it chooses a family of simple distributions (independent normals, say) and searches that family for the member closest to the posterior. The search is a gradient-based optimisation, the same kind that fits a neural network, and it is usually one or two orders of magnitude faster than MCMC. The price is an approximation whose errors are not random but systematic and predictable: variational posteriors are typically centred in the right place and too narrow. This chapter is about the payoff and the price, and about how to tell, for a given model, whether the approximation can be trusted.

Lecture 17 is about the idea. We will turn inference into optimisation via the Kullback–Leibler divergence and the evidence lower bound (ELBO), see why the direction of the divergence makes VI overconfident, why the mean-field family cannot represent correlations, and watch the stochastic-gradient algorithm that Stan uses (ADVI) find a posterior step by step.
Lecture 18 is about practice. We will predict cancellations among 119,000 hotel bookings with HMC and with VI, time both, compare the answers, watch VI fail on a multilevel model, diagnose the failure with Pareto- $\hat{k}$ , and survey the modern alternatives (Laplace, Pathfinder, stochastic VI) and the rules for reporting approximate inference.

Think About It: Before You Start

Suppose the true posterior for a parameter is bimodal: two separate humps of equal height. You are only allowed to report a single normal distribution. Would you place it over one of the humps, narrow, or between them, wide enough to cover both? Which choice gives better point estimates? Which gives more honest uncertainty? Variational inference makes this choice for you, and Section 9.3 shows which one.

Lecture 17 — Inference as optimisation

9.1 Will this booking be cancelled?

Hotels overbook because guests cancel, and a good estimate of each booking's probability of cancellation is worth a great deal of money. Nuno Antonio, Ana de Almeida and Luis Nunes (2019) published the booking records of two Portuguese hotels, a city hotel in Lisbon and a resort in the Algarve, for July 2015 to August 2017: 119,390 bookings, of which 37% were cancelled. The data were distributed through TidyTuesday. Each record has dozens of variables, and the first modelling decision is a causal one.

library(tidyverse)
library(brms)
library(posterior)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

Warning: Leakage Is a Causal Error

Several variables in the hotel file are recorded after the outcome: the reservation status ("Canceled", "Check-Out"), the room type actually assigned, and the number of changes made to the booking (which keeps counting until the booking ends). In a causal graph they are descendants of cancellation, not causes of it. A model that uses them predicts cancellations superbly in the data set and is useless for a new booking, where they are not yet known. This is called leakage, and it is one of the commonest errors in published prediction models. The remedy is the causal graph of Chapter 1, drawn with attention to time: include when each variable is measured, and use only what is known at the moment the prediction is needed.

We keep variables known when the booking is made: how far ahead it was made (lead time), the deposit type, the price per night, the number of special requests, whether a parking space was requested, whether the guest has cancelled before or stayed before, which hotel, and the guest's country.

url <- "https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-02-11/hotels.csv"
hotels <- read_csv(url, show_col_types = FALSE) |>
  filter(adr >= 0, adr < 1000, adults + children + babies > 0) |>        # drop data errors
  transmute(cancel   = is_canceled,
            hotel    = if_else(hotel == "City Hotel", "city", "resort"),
            lead     = as.numeric(scale(log1p(lead_time))),              # standardised log lead time
            deposit  = factor(deposit_type, c("No Deposit", "Non Refund", "Refundable")),
            prev     = as.integer(previous_cancellations > 0),           # has cancelled before
            repeated = is_repeated_guest,
            requests = pmin(total_of_special_requests, 3),
            parking  = as.integer(required_car_parking_spaces > 0),
            price    = as.numeric(scale(log1p(adr))),                    # standardised log price per night
            country)
c(bookings = nrow(hotels), cancelled = mean(hotels$cancel), countries = n_distinct(hotels$country))
hotels |> group_by(deposit) |> summarise(bookings = n(), cancelled = mean(cancel))

output

 bookings cancelled countries 
1.192e+05 3.707e-01 1.780e+02 
# A tibble: 3 x 3
  deposit    bookings cancelled
  <fct>         <int>     <dbl>
1 No Deposit   104456     0.284
2 Non Refund    14586     0.994
3 Refundable      162     0.222

After removing a few hundred records with impossible values (negative prices, no guests), 119,204 bookings remain, from 178 countries, and 37% were cancelled. One pattern stands out at once: bookings with a non-refundable deposit, 14,586 of them, were cancelled 99.4% of the time, against 28% for bookings with no deposit. That is not how non-refundable deposits are supposed to work, and Exercise 1 investigates it; for now, notice that a model will find this variable enormously predictive, and that a prediction learned from it may not transfer to other hotels.

9.2 The cost of sampling

A logistic regression with eight predictors is the simplest model in this course, and HMC handles it without any difficulty, except time. Each leapfrog step evaluates the log likelihood and its gradient over all 119,000 bookings; each iteration takes several leapfrog steps; and we ask for 2,000 iterations in each of four chains. Here is the fit, with its timing.

f <- cancel ~ lead + deposit + prev + repeated + requests + parking + price + hotel
pr <- c(prior(normal(0, 1.5), class = Intercept), prior(normal(0, 1), class = b))
m0 <- brm(f, data = hotels, family = bernoulli(), prior = pr, chains = 0)   # compile only
time_hmc <- system.time(
  m_hmc <- update(m0, chains = 4, seed = 9, refresh = 0)                    # 4 chains x 2,000 iterations
)["elapsed"]
fixef(m_hmc) |> round(3)
c(seconds = time_hmc)

output

                  Estimate Est.Error   Q2.5  Q97.5
Intercept           -0.686     0.011 -0.708 -0.665
lead                 0.475     0.009  0.458  0.493
depositNonRefund     5.095     0.103  4.896  5.300
depositRefundable   -0.511     0.199 -0.896 -0.121
prev                 3.895     0.091  3.719  4.078
repeated            -2.127     0.101 -2.328 -1.934
requests            -0.450     0.010 -0.470 -0.430
parking             -6.210     0.389 -7.029 -5.503
price                0.436     0.012  0.413  0.459
hotelresort          0.048     0.017  0.016  0.079
seconds.elapsed 
          237.3

Every coefficient is estimated precisely, with standard deviations of about 0.01 for the continuous predictors, because there are so many bookings. Longer lead times and higher prices go with more cancellations; special requests, parking requests and being a repeat guest go with fewer; a previous cancellation raises the odds enormously. On this machine, with two cores, sampling took about four minutes. That is tolerable once, but not if the model is to be refitted every night with new bookings, or if twenty candidate models are to be compared, or if the data had ten million rows instead of a hundred thousand. Each of those multiplies the cost.

9.3 Inference as optimisation

Variational inference starts by choosing a family $𝒬$ of distributions that are easy to work with: for example, all normal distributions with independent components. It then looks for the member $q^{⋆} \in 𝒬$ that is closest to the posterior $p (θ ∣ y)$ . "Closest" is measured by the Kullback–Leibler divergence,

KL (q ‖ p) = \int q (θ) \log \frac{q (θ)}{p (θ ∣ y)} d θ = 𝔼_{q} [\log q (θ)] - 𝔼_{q} [\log p (θ ∣ y)],

which is never negative and is zero only when $q = p$ . There is an obvious obstacle: the divergence involves the posterior, $p (θ ∣ y) = p (y, θ) / p (y)$ , and the normalising constant $p (y)$ is exactly the integral we could not compute. The way around it is a rearrangement. Substituting and noting that $\log p (y)$ does not depend on $θ$ ,

\log p (y) = \underset{ELBO (q)}{\underset{⏟}{𝔼_{q} [\log p (y, θ)] - 𝔼_{q} [\log q (θ)]}} + KL (q ‖ p (θ ∣ y)) .

The left side is a constant. So making the ELBO larger makes the divergence smaller by exactly the same amount: minimising the KL divergence is the same as maximising the ELBO. And the ELBO involves only the joint density $p (y, θ)$ , which is likelihood times prior, and $q$ , which we chose. Because the divergence is never negative, the ELBO is a lower bound on the log evidence $\log p (y)$ : the evidence lower bound. The two terms of the ELBO have a natural reading. The first rewards $q$ for putting its mass where the joint density is high; the second, the entropy of $q$ , rewards it for being spread out. The best $q$ balances the two.

Definition: The Evidence Lower Bound

For a family $𝒬$ and a model with joint density $p (y, θ)$ , the ELBO of $q \in 𝒬$ is $ELBO (q) = 𝔼_{q} [\log p (y, θ)] + ℍ (q)$ , where $ℍ (q) = - 𝔼_{q} [\log q (θ)]$ is the entropy of $q$ . It satisfies $ELBO (q) = \log p (y) - KL (q ‖ p (θ ∣ y)) \leq \log p (y)$ . Variational inference reports $q^{⋆} = arg {max}_{q \in 𝒬} ELBO (q)$ as the approximate posterior.

The direction of the divergence matters. $KL (q ‖ p)$ averages $\log (q / p)$ over the places where q puts its mass. Wherever $q$ has mass and $p$ has almost none, $\log (q / p)$ is huge, so $q$ is severely punished for straying into regions the posterior considers implausible. Wherever $p$ has mass and $q$ has none, the integrand is zero and $q$ pays nothing. The safe strategy for $q$ is therefore to stay inside the posterior's bulk and leave the tails (or a whole second mode) uncovered. This is called mode-seeking or zero-forcing behaviour. The reverse divergence, $KL (p ‖ q)$ , behaves the opposite way: it punishes $q$ for missing any region where $p$ has mass, so its minimiser spreads out to cover everything (it matches the mean and variance of $p$ ). We cannot use the reverse divergence directly, because it requires expectations under the posterior we do not have. The widget shows the difference.

When the two modes overlap, both approximations are similar: one wide normal over both. As they separate, there is an abrupt switch, at a separation of about five standard deviations with the default weights (a little more when the weights are equal), after which the variational solution (blue) locks onto the heavier mode and becomes narrow: its standard deviation is close to that of a single mode, and it reports nothing about the other. The moment-matching solution sits between the modes, in a region where the target has little mass, and is wide enough to cover both. Neither is a good description of a bimodal posterior; the lesson is about what VI does when its family is too simple. It prefers being confidently partly right to being vaguely right everywhere, and its uncertainty is therefore too small.

Check Yourself

Why can variational inference maximise the ELBO even though it cannot compute the posterior's normalising constant

p (y)

?Because it estimates

p (y)

by Monte Carlo first.Because

p (y) = 1

for any proper prior.Because

\log p (y)

does not depend on

q

: the ELBO uses only the joint density

p (y, θ)

, and differs from

\log p (y)

by exactly the KL divergence, so maximising one minimises the other.

9.4 Mean-field and full-rank families

The family $𝒬$ decides how good the approximation can possibly be. The most common choice is the mean-field family: $q$ is a product of independent one-dimensional normals, one per parameter (on the unconstrained scale, so a standard deviation becomes a normal on its logarithm). It is fast, scales to millions of parameters, and has one structural blind spot: it cannot represent correlation. The full-rank family uses a single multivariate normal with an arbitrary covariance matrix, which can tilt to follow correlations, at a cost that grows with the square of the number of parameters.

For a correlated normal posterior, the mean-field solution can be found exactly, and it is instructive. If the posterior is normal with precision matrix $Λ = Σ^{- 1}$ , the best mean-field $q$ has the right means, but its variances are $1 / Λ_{i i}$ , the conditional variances of each parameter given the others, rather than the marginal variances $Σ_{i i}$ . For two parameters with correlation $ρ$ , that shrinks each standard deviation by a factor $\sqrt{1 - ρ^{2}}$ .

The mean-field ellipse is axis-aligned and centred correctly; it reaches only a fraction of the way along the posterior's long direction (and slightly past it across the short one), and its marginal is narrower than the true marginal by $\sqrt{1 - ρ^{2}}$ : about 60% of the true width at $ρ = 0.8$ , and 30% at $ρ = 0.95$ . Correlated posteriors are the rule, not the exception: the coefficients of correlated predictors, intercepts and slopes when predictors are not centred, group effects and their standard deviation in multilevel models. Notice that the full-rank family fixes the correlation problem but not the mode-seeking problem of Section 9.3, which comes from the direction of the divergence, not the family.

Check Yourself

Two regression coefficients have a normal posterior with correlation $- 0.9$ . By what factor does mean-field VI shrink the standard deviation of each?

ratio =

9.5 ADVI, step by step

How is the ELBO maximised in practice? It is an expectation under $q$ , which usually cannot be computed exactly, but it can be estimated by drawing from $q$ . The key device is the reparameterisation trick: write a draw from $q$ as a fixed transformation of standard normal noise, $θ = m + L ε$ with $ε \sim Normal (0, I)$ , where $m$ is the mean and $L$ is diagonal (mean-field) or lower-triangular (full-rank). Then

\nabla_{m, L} ELBO = 𝔼_{ε} [\nabla_{m, L} \log p (y, m + L ε)] + \nabla_{m, L} ℍ (q),

and the expectation can be estimated with a handful of draws of $ε$ , using the same gradient of the log posterior that HMC uses. The entropy of a normal is $\sum_{i} \log L_{i i}$ plus a constant, whose gradient is exact. Stochastic gradient ascent on $(m, L)$ then climbs the ELBO. Kucukelbir et al. (2017) combined this with automatic differentiation and automatic transformation of constrained parameters to the real line, and called it automatic differentiation variational inference (ADVI). It is what Stan and brms run when you ask for algorithm = "meanfield" or "fullrank".

The widget runs ADVI live, in your browser, on a two-parameter version of the hotel model fitted to a random sample of 2,000 bookings: the probability of cancellation as a function of lead time alone, with an intercept and a slope.

Run 500 steps with each family. The ELBO climbs steeply from the poor starting point and then fluctuates around a plateau; the fluctuations are the Monte Carlo noise in the gradient, and they are why ADVI's answer changes from run to run and why the step size must shrink near the end. Both families find the right centre. The full-rank ellipse tilts to match the posterior; the mean-field ellipse cannot, and it is slightly narrower than the posterior in both directions (the correlation here is only about −0.3, so the squeeze is small). Try a large learning rate: the optimiser overshoots and the ELBO trace becomes erratic. Diagnosing convergence of a stochastic optimiser is harder than diagnosing MCMC, because there is no equivalent of $\hat{R}$ ; Stan stops when the relative change in the ELBO falls below a tolerance, which can happen too early.

History: From Physics to Big Data

The mean-field idea comes from statistical physics: to study a magnet with billions of interacting atoms, Pierre Weiss (1907) replaced each atom's interactions with its neighbours by their average, a "mean field", turning an impossible coupled problem into many independent ones. Richard Feynman's variational principle in quantum mechanics bounds an intractable free energy by an expectation under a simpler distribution, which is the ELBO in another notation. In the 1990s, Michael Jordan, Zoubin Ghahramani, Tommi Jaakkola and Lawrence Saul (1999) brought these ideas into machine learning as a way to fit large graphical models. David Blei, Andrew Ng and Jordan's latent Dirichlet allocation (2003), the topic model of Chapter 12, was fitted by variational inference, and Hoffman et al. (2013) made it scale to millions of documents with stochastic VI, which estimates the gradient from a random minibatch of the data at each step. The reparameterisation trick (Kingma and Welling, 2014) made VI the engine of deep generative models such as the variational autoencoder, and ADVI (Kucukelbir et al., 2017) brought it to general Bayesian models.

Key Takeaways: Lecture 17

Variational inference finds the member of a simple family closest to the posterior in KL divergence, by maximising the ELBO, which needs only the joint density.
The direction $KL (q ‖ p)$ makes VI mode-seeking: it stays inside the posterior and underestimates its spread.
The mean-field family cannot represent correlation and shrinks each standard deviation by up to $\sqrt{1 - ρ^{2}}$ ; full-rank can, at quadratic cost.
ADVI maximises the ELBO by stochastic gradient ascent with the reparameterisation trick. Its answer is noisy and its convergence hard to diagnose.

Lecture 18 — Variational inference in practice

9.6 VI and HMC on 119,000 bookings

In brms, the only change needed is the algorithm argument. We fit the hotel model with both variational families, reusing the compiled model so that the timings compare the algorithms alone.

time_mf <- system.time(
  m_mf <- update(m0, algorithm = "meanfield", seed = 9, refresh = 0)
)["elapsed"]
time_fr <- system.time(
  m_fr <- update(m0, algorithm = "fullrank", seed = 9, refresh = 0)
)["elapsed"]
c(hmc = time_hmc, meanfield = time_mf, fullrank = time_fr)

output

      hmc.elapsed meanfield.elapsed  fullrank.elapsed 
          237.325             6.038             7.507

summ <- function(m, method) as_draws_df(m) |> select(starts_with("b_")) |>
  pivot_longer(everything(), names_to = "coef") |>
  group_by(coef) |> summarise(mean = mean(value), sd = sd(value),
                              lo = quantile(value, .025), hi = quantile(value, .975)) |>
  mutate(method = method)
cmp <- bind_rows(summ(m_hmc, "HMC"), summ(m_mf, "mean-field"), summ(m_fr, "full-rank")) |>
  mutate(coef = str_remove(coef, "^b_"), method = factor(method, c("HMC", "mean-field", "full-rank")))
ratio <- cmp |> group_by(coef) |> mutate(sd_ratio = sd / sd[method == "HMC"]) |> filter(method != "HMC")
cols3 <- c(HMC = harbour[["ink"]], `mean-field` = harbour[["sky"]], `full-rank` = harbour[["amber"]])
p1 <- ggplot(cmp, aes(y = coef, colour = method)) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, position = position_dodge(width = 0.6)) +
  geom_point(aes(x = mean), position = position_dodge(width = 0.6), size = 1.3) +
  scale_colour_manual(values = cols3) + labs(x = "coefficient (log odds)", y = NULL, colour = NULL)
p2 <- ggplot(ratio, aes(sd_ratio, coef, colour = method)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = harbour[["muted"]]) +
  geom_point(size = 2) + scale_colour_manual(values = cols3, guide = "none") +
  scale_x_log10() + labs(x = "sd relative to HMC (log)", y = NULL)
patchwork::wrap_plots(p1, p2, widths = c(1.5, 1))
cmp |> select(coef, method, mean, sd) |> pivot_wider(names_from = method, values_from = c(mean, sd)) |>
  mutate(across(where(is.numeric), ~ round(.x, 3)))

output

# A tibble: 10 x 7
   coef         mean_HMC `mean_mean-field` `mean_full-rank` sd_HMC `sd_mean-field`
   <chr>           <dbl>             <dbl>            <dbl>  <dbl>           <dbl>
 1 Intercept      -0.686            -0.646           -0.822  0.011           0.014
 2 depositNonR~    5.10              5.08             5.20   0.103           0.044
 3 depositRefu~   -0.511            -0.542           -0.475  0.199           0.243
 4 hotelresort     0.048             0.014            0.077  0.017           0.018
 5 lead            0.475             0.478            0.469  0.009           0.01 
 6 parking        -6.21             -5.94            -6.92   0.389           0.103
 7 prev            3.90              3.90             4.03   0.091           0.051
 8 price           0.436             0.467            0.384  0.012           0.011
 9 repeated       -2.13             -2.08            -2.09   0.101           0.069
10 requests       -0.45             -0.472           -0.428  0.01            0.011
# i 1 more variable: `sd_full-rank` <dbl>

Posterior means and 95% intervals for the ten coefficients of the hotel model (including the intercept), from HMC (reference), mean-field VI and full-rank VI. The centres agree; the interval widths, shown in the right panel as a ratio to HMC, do not.

Mean-field ADVI took about 6 seconds and full-rank about 7.5, against 237 for HMC: roughly 30 to 40 times faster. (Both were fitted with the same compiled model, so the comparison is fair.) Mean-field ADVI found the centres well: every posterior mean is within about 0.05 of HMC's except the extreme parking coefficient (−5.9 against −6.2). Its standard deviations are a mixed picture. For the precisely estimated continuous predictors they are about right, but for the coefficients of rare categories they are much too small: 0.044 against 0.103 for the non-refundable deposit, 0.051 against 0.091 for previous cancellations, and 0.10 against 0.39 for parking, where only a handful of bookings with a parking request were cancelled and the posterior is skewed, which no normal approximation can capture. Full-rank ADVI did worse, not better: its standard deviations are too large for most coefficients (0.19 against 0.011 for the intercept) and its means drift. With ten correlated coefficients and the default convergence tolerance, the stochastic optimiser stopped before it had converged. A richer family is only better if the optimiser can actually find its best member.

Check Yourself

A variational fit and an HMC fit give the same posterior means but the VI intervals are half as wide. For deciding which predictors matter most, VI is

9.7 When VI fails, and how to tell

The fixed-effects model is easy for VI: lots of data, few parameters, a nearly normal posterior. Multilevel models are harder. Add a varying intercept for the guest's country (there are 178 countries, from Portugal with tens of thousands of bookings to dozens with one or two), and the posterior acquires the funnel of Chapter 6: when the between-country standard deviation $τ$ is small, all the country effects must be small too. That is a strong correlation between $τ$ and every country effect, which the mean-field family cannot represent.

pr_ml <- c(pr, prior(exponential(1), class = sd))
f_ml <- update(f, ~ . + (1 | country))
m_ml_hmc <- brm(f_ml, data = hotels, family = bernoulli(), prior = pr_ml, seed = 9, refresh = 0,
                file = "fits/ch9_ml_hmc")                  # about 35 minutes on a two-core machine
m0_ml <- brm(f_ml, data = hotels, family = bernoulli(), prior = pr_ml, chains = 0)
time_ml_mf <- system.time(m_ml_mf <- update(m0_ml, algorithm = "meanfield", seed = 9, refresh = 0))["elapsed"]
c(meanfield_seconds = time_ml_mf)
rbind(HMC = posterior_summary(m_ml_hmc, variable = "sd_country__Intercept"),
      meanfield = posterior_summary(m_ml_mf, variable = "sd_country__Intercept")) |> round(3)

output

meanfield_seconds.elapsed 
                    8.659 
                      Estimate Est.Error  Q2.5 Q97.5
sd_country__Intercept    1.096     0.099 0.924 1.315
sd_country__Intercept    0.620     0.008 0.604 0.635

The difference is dramatic. HMC needed about 35 minutes for this model, and even then mixed slowly for the intercept (because Portugal, with tens of thousands of bookings, is tightly linked to the overall intercept). Mean-field ADVI took 9 seconds. But it estimates the between-country standard deviation $τ$ at 0.62 (95% interval 0.60 to 0.64), against HMC's 1.10 (0.92 to 1.31): too small by almost half, and with an interval so narrow that it excludes the right answer by a wide margin. Countries with few bookings would be shrunk far too strongly toward the average, and their predicted cancellation probabilities would be wrong and overconfident (Exercise 5). This is the typical failure of mean-field VI on hierarchical models, and it is exactly the kind of quantity, a variance component, that a multilevel analysis exists to estimate.

How could we know that the approximation is poor, without running HMC to compare? Yao, Vehtari, Simpson and Gelman (2018) proposed a diagnostic that uses only the variational fit. Draw $θ^{(s)}$ from $q$ , and compute the importance ratios $r_{s} = p (y, θ^{(s)}) / q (θ^{(s)})$ . If $q$ were exactly the posterior, the ratios would all be equal. If $q$ is too narrow, a few draws in the tails of $q$ will have enormous ratios, because the posterior has much more mass there than $q$ expects. Fit a generalised Pareto distribution to the largest ratios, the same device that Pareto-smoothed importance sampling (PSIS) uses to check approximate leave-one-out cross-validation in loo(); its shape parameter $\hat{k}$ measures how heavy the tail of the ratios is. Values below 0.7 mean that the approximation is good enough for importance sampling to correct it; values above 0.7 mean that $q$ misses important parts of the posterior. Stan computes $\hat{k}$ at the end of every ADVI run, and brms keeps it in the fitted object.

pareto_k <- function(m) m$fit@sim$diagnostics$psis$pareto_k   # computed by Stan's ADVI, stored by brms
c(fixed_meanfield = pareto_k(m_mf), fixed_fullrank = pareto_k(m_fr),
  multilevel_meanfield = pareto_k(m_ml_mf)) |> round(2)

output

     fixed_meanfield       fixed_fullrank multilevel_meanfield 
                2.06                17.10                 2.87

All three variational fits fail the diagnostic, two of them badly. The mean-field fixed-effects model has $\hat{k} \approx 2.1$ , which is consistent with what we saw: good means, some standard deviations far too small. The full-rank fit has $\hat{k} \approx 17$ , signalling that its optimiser did not converge. The mean-field multilevel model has $\hat{k} \approx 2.9$ . In every case the diagnostic, which uses only the variational fit, correctly warns that the uncertainty cannot be trusted. Stan prints the same value as a warning, which is easy to ignore; do not. Tightening the convergence tolerance (tol_rel_obj = 0.001) and running more iterations helps a little (in our trials, mean-field $\hat{k}$ fell to 1.4), but here no setting brought it below 0.7, and the honest conclusion is that this posterior needs MCMC, or at least validation against it.

9.8 Beyond ADVI

ADVI is one of several ways to approximate a posterior quickly, and in recent years it has lost ground to methods that are more robust or more accurate.

The Laplace approximation finds the posterior mode by optimisation and uses the curvature (Hessian) there as the covariance of a normal approximation. It is the oldest method of all (Laplace used it in 1774), deterministic, and often better than ADVI for models with many observations per parameter. INLA (Rue, Martino and Chopin, 2009) builds an extremely fast and accurate method for latent Gaussian models (including many spatial and time-series models) by nesting Laplace approximations.
Pathfinder (Zhang et al., 2022) runs a quasi-Newton optimiser toward the mode, fits a normal approximation at each point along the path, and picks the one with the highest ELBO. It is usually faster and more reliable than ADVI and is increasingly used in Stan to initialise MCMC chains; it is available in cmdstanr and, through it, in recent versions of brms.
Stochastic VI estimates the gradient from a random minibatch of the data, so that each step costs a fraction of a pass through the data. It is what makes VI feasible for millions of observations, and it is the default in probabilistic programming languages built on deep-learning frameworks (Pyro and NumPyro on PyTorch and JAX, for example).
Richer families, from mixtures of normals to normalising flows (invertible neural networks that transform a simple distribution into a flexible one), reduce the approximation error of the family, at the cost of much harder optimisation.

Research Practice: Reporting Approximate Inference

State the algorithm (MCMC, ADVI mean-field or full-rank, Laplace, Pathfinder, INLA) and the software, with its convergence settings.
For VI, report the Pareto- $\hat{k}$ diagnostic. If it exceeds 0.7, treat the uncertainty as unreliable, and if possible improve the approximation (full-rank, importance sampling) or switch method.
Validate against MCMC on a subsample or a simplified model, and say so: "on a 10% subsample, mean-field VI recovered the HMC posterior means to within 0.01 and the standard deviations to within 15%".
Use VI's point estimates freely for exploration; do not base claims about uncertainty (intervals, tail probabilities, predictive intervals, multilevel variance components) on VI without validation.
Final results in a paper should come from MCMC whenever it is computationally possible. "It would have taken a week" is a reason to use VI; "it would have taken ten minutes" is not.

9.9 Approximate inference in practice

A sensible workflow uses both kinds of inference, each for what it is good at. During model building, when a model is fitted, criticised and changed dozens of times, fast approximate fits (VI, Pathfinder or Laplace) let you iterate in seconds; their point estimates are usually reliable enough to steer the next step, and their failures are usually visible in the diagnostics. For the final model, run MCMC if it is at all feasible, and use its results in the paper. When MCMC is not feasible, because the data are enormous or the model must be refitted continuously in production, VI is what makes Bayesian inference possible at all, and a fast approximate posterior, honestly labelled and validated on a subsample, is much better than an exact posterior that never finishes.

Project Connections: Logs, Text and Fast Iteration

The hostile-traffic detection project is the natural home for this chapter: authentication and web-server logs arrive at millions of lines a day, attack patterns shift, and a detector must be retrained often. A Bayesian logistic or Poisson model of request patterns with stochastic VI can be retrained continuously, and still attaches a probability to each alert; validate it periodically against MCMC on a sample of a day's traffic, and report the Pareto- $\hat{k}$ . In the AI-generated text detection project, models with thousands of word or embedding features and many candidate specifications benefit from fast VI fits during feature selection, with MCMC reserved for the final model on the chosen features. In both, the multilevel lesson of Section 9.7 applies: if the model has group-level effects (by user, by IP block, by course), check the variance components against MCMC before trusting them.

Key Takeaways: Lecture 18

In brms, algorithm = "meanfield" or "fullrank" switches to ADVI. On a large, simple model it is many times faster than HMC, with nearly identical posterior means.
VI's standard deviations are typically too small, and for multilevel models its variance components can be badly wrong.
The Pareto- $\hat{k}$ of the importance ratios $p / q$ diagnoses a poor approximation without MCMC; values above 0.7 are a warning.
Laplace, INLA and Pathfinder are fast alternatives; stochastic VI scales to millions of observations. Use fast methods to explore and MCMC to report, whenever possible.

9.10 Exercises

Each exercise has a hint and an answer button.

Exercises

The non-refundable puzzle. Bookings with a non-refundable deposit are cancelled almost 100% of the time, which seems absurd: why would someone pay in full and then cancel? Look at the lead times, market segments and countries of these bookings. What is a plausible explanation, and what does it imply for using the hotel data to predict cancellations elsewhere?

Compare the distribution of lead times and the share from travel agencies (market_segment) for non-refundable and other bookings in the raw data.

The 14,586 non-refundable bookings are nothing like the rest. Their median lead time is 183 days against 57; 97% came through groups or offline travel agents and tour operators, against 29% of other bookings; 97% are from Portugal, against a third; and 88% are at the city hotel. The most plausible explanation is that these are blocks of rooms reserved in advance by tour operators and groups, recorded with a non-refundable deposit type, most of which were later released. Whatever the mechanism, the "deposit" variable here is a proxy for a booking channel peculiar to these two hotels, not a measure of a guest's financial commitment. A model trained on it would predict nonsense for a hotel where non-refundable deposits mean what they say. This is a question about the causal and measurement process behind a variable, and no amount of fitting can answer it.
Mean-field variance by formula. From the HMC draws of m_hmc, compute the posterior covariance matrix of the coefficients, invert it to get the precision matrix, and compute the mean-field standard deviations $1 / \sqrt{Λ_{i i}}$ predicted by Section 9.4. Compare them with the mean-field VI standard deviations from m_mf.

S <- cov(as_draws_df(m_hmc) |> select(starts_with("b_"))), then 1 / sqrt(diag(solve(S))).

The formula $1 / \sqrt{Λ_{i i}}$ applied to the HMC posterior predicts some of mean-field VI's standard deviations well: for the correlated pair of previous cancellation and repeat guest (correlation −0.71), it gives 0.063 and 0.069 against VI's 0.051 and 0.069, and HMC's 0.091 and 0.101. For the rare categories it does not: for the non-refundable deposit the formula gives 0.101 and VI 0.044; for parking 0.387 and VI 0.103. Those posteriors are far from normal (skewed by the near-separation of rare categories), and in addition the optimiser had not fully converged (VI's intercept standard deviation, 0.014, is larger than either). The correlation formula explains part of mean-field's narrowness; non-normality and imperfect optimisation explain the rest.
Run it again. Refit the mean-field model five times with different seeds. How much do the posterior means and standard deviations change between runs? Compare this with the Monte Carlo standard errors of the HMC fit.

Loop over seed = 1:5 with update(m_mf, seed = s).

Across five seeds, the mean-field posterior means vary by up to 0.07 for the precisely estimated coefficients (lead time, price, requests) and by 0.17 for parking; the posterior standard deviations vary by factors of up to two; and Pareto- $\hat{k}$ ranges from 1.6 to 5.5. HMC's Monte Carlo standard errors for the same means are about 0.0001 for lead time and 0.009 for parking. The run-to-run variation of ADVI is hundreds of times larger than HMC's Monte Carlo error, and comparable to the posterior uncertainty itself. If you use ADVI, run it several times and report the spread; a single run's numbers are not reproducible to the digits printed.
Subsample validation. Fit the model with HMC to a random 10% of the bookings and with mean-field VI to the full data. Using the 10% fit, predict how wide the HMC posterior on the full data would be (standard deviations scale as $1 / \sqrt{n}$ ), and compare with VI's standard deviations. Is this a useful validation strategy?

Multiply the 10%-sample HMC standard deviations by $\sqrt{0.1}$ .

Scaling the 10%-sample HMC standard deviations by $\sqrt{0.1}$ predicts the full-data HMC values well for most coefficients (lead time 0.009 against 0.009; price 0.012 against 0.012; previous cancellations 0.086 against 0.091), and so exposes the mean-field values that are too small (previous cancellations 0.051, non-refundable deposit 0.044). It fails for the rare categories: it predicts 0.139 for parking, where the truth is 0.389, because with few events the posterior does not shrink like $1 / \sqrt{n}$ . Subsample validation is useful and cheap, and it catches most problems, but it is least reliable exactly where VI is also least reliable: sparse categories and skewed posteriors.
Predictive probabilities. For the multilevel model, compare the predicted cancellation probabilities of a booking from a country with only two bookings in the data under HMC and under mean-field VI (posterior mean and 90% interval). Which is more honest?

Pick a country with n == 2 and use posterior_epred(..., newdata = ...) with typical values of the other predictors.

For a booking from Aruba (two bookings in the data), with typical values of the other predictors, the multilevel HMC fit gives a predicted cancellation probability of 0.26 with a 90% interval from 0.06 to 0.56; mean-field VI gives 0.30 with a 90% interval from 0.15 to 0.49. HMC's interval is about 40% wider, because it knows that the country effect of a country with two bookings is very uncertain (and that $τ$ is large). VI, with $τ$ underestimated and every uncertainty squeezed, is more confident than two bookings can justify. The HMC answer is the honest one.
The ELBO as a model comparison. The ELBO is a lower bound on $\log p (y)$ , the quantity in a Bayes factor. Why is it a poor tool for comparing two models, even if both are fitted by VI?

Think about what the gap between the ELBO and $\log p (y)$ is, and whether it is the same for two models.

The gap between the ELBO and $\log p (y)$ is the KL divergence between $q$ and the posterior, which is unknown and differs between models: a model whose posterior is close to normal has a small gap, and one with a skewed, correlated or multimodal posterior has a large one. Comparing ELBOs therefore compares $\log p (y)$ minus two different unknown amounts, and can prefer the model whose posterior happens to be easier to approximate. In addition, the ELBO from ADVI is a noisy stochastic estimate, and $\log p (y)$ itself is sensitive to the prior in ways that make Bayes factors fragile. For predictive model comparison, use leave-one-out or leave-future-out cross-validation instead.
Laplace by hand. For the two-parameter model of the widget in Section 9.5 (intercept and lead-time slope, on 2,000 bookings), find the posterior mode with optim(), compute the Hessian there, and use its inverse as the covariance of a normal approximation. Compare it with HMC on the same 2,000 bookings.

optim(c(0, 0), negative_log_posterior, hessian = TRUE), then solve(fit$hessian).

The Laplace approximation (mode and inverse Hessian) gives intercept −0.580 (sd 0.051) and slope 0.811 (sd 0.061), with correlation −0.30. HMC on the same 2,000 bookings gives −0.582 (0.050), 0.814 (0.061) and −0.30. The two agree to two or three decimal places, because with 2,000 observations and two parameters the posterior is very nearly normal. The Laplace approximation is deterministic, cost a fraction of a second, and has no tuning; for models with many observations per parameter and a nearly normal posterior, it is often the best fast approximation, better than ADVI.
A methods paragraph. You fitted a model to ten million rows with mean-field VI because MCMC would take weeks. Write the methods paragraph (four or five sentences) reporting the inference and its validation, following the Research Practice box in Section 9.8.

Algorithm, software and settings; the Pareto- $\hat{k}$ ; a subsample comparison with MCMC; what the reader should and should not rely on.

Example. Because the data set contained ten million rows and full MCMC was estimated to require several weeks, we fitted the model by mean-field automatic differentiation variational inference (ADVI) in Stan 2.32 via brms 2.20, with a relative convergence tolerance of 0.001 and 50,000 maximum iterations. The Pareto- $\hat{k}$ diagnostic of Yao et al. (2018) was 0.45, below the threshold of 0.7, indicating that the variational approximation was adequate for importance-sampling correction. To validate it, we fitted the same model with the No-U-Turn sampler (4 chains of 2,000 iterations, all $\hat{R} \leq 1.01$ ) to a random 1% subsample, and to the same subsample with ADVI: posterior means agreed to within 0.02 and posterior standard deviations to within 15% for all coefficients. Repeating ADVI with five random seeds changed posterior means by at most 0.01. We report posterior means and 95% intervals from ADVI; because variational approximations can understate uncertainty, intervals and tail probabilities should be regarded as approximate.

Chapter 9
Variational inference

Think About It: Before You Start

Lecture 17 — Inference as optimisation

9.1 Will this booking be cancelled?

Warning: Leakage Is a Causal Error

9.2 The cost of sampling

9.3 Inference as optimisation

Definition: The Evidence Lower Bound

Check Yourself

9.4 Mean-field and full-rank families

Check Yourself

9.5 ADVI, step by step

History: From Physics to Big Data

Key Takeaways: Lecture 17

Lecture 18 — Variational inference in practice

9.6 VI and HMC on 119,000 bookings

Check Yourself

9.7 When VI fails, and how to tell

9.8 Beyond ADVI

Research Practice: Reporting Approximate Inference

9.9 Approximate inference in practice

Project Connections: Logs, Text and Fast Iteration

Key Takeaways: Lecture 18

9.10 Exercises

Exercises

Further reading and references

Further Reading

Chapter 9Variational inference

Think About It: Before You Start

Lecture 17 — Inference as optimisation

9.1 Will this booking be cancelled?

Warning: Leakage Is a Causal Error

9.2 The cost of sampling

9.3 Inference as optimisation

Definition: The Evidence Lower Bound

Check Yourself

9.4 Mean-field and full-rank families

Check Yourself

9.5 ADVI, step by step

History: From Physics to Big Data

Key Takeaways: Lecture 17

Lecture 18 — Variational inference in practice

9.6 VI and HMC on 119,000 bookings

Check Yourself

9.7 When VI fails, and how to tell

9.8 Beyond ADVI

Research Practice: Reporting Approximate Inference

9.9 Approximate inference in practice

Project Connections: Logs, Text and Fast Iteration

Key Takeaways: Lecture 18

9.10 Exercises

Exercises

Further reading and references

Further Reading

Chapter 9
Variational inference