Introduction to Probabilistic Graphical Models and Deep Generative Models

Contents

• Introduction and overview
• 0. Review: probability tools for probabilistic models
• 1. Maximum likelihood, linear models and K-means
• 2. Latent variables: EM, Gaussian mixtures and probabilistic PCA
• 3. Bayesian linear models, Gaussian processes and model selection
• 4. Directed and undirected graphical models
• 5. Exact inference: the sum-product algorithm and hidden Markov models
• 6. Approximate inference I: Monte Carlo and MCMC
• 7. Approximate inference II: variational methods, stochastic block models and expectation propagation
• 8. Deep latent variable models, variational auto-encoders and amortised inference
• 9. Deep generative models beyond VAEs: autoregressive models, normalizing flows and GANs
• 10. Energy-based models and score-based diffusion models
• 11. Causal inference: potential outcomes, treatment effects and policy learning
• 12. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, statistics, computer science or physics who knows probability and linear algebra, and wants one language for latent variables, inference and generative modelling.

Empirical first: in the observational study simulated in ch. 11, the naive difference of means estimates a treatment effect of 1.95 when the true effect is 1.00. On the two-moons data of scikit-learn, the best Gaussian reaches an average log-likelihood of −2.73 per point, a six-layer normalizing flow −1.21 (ch. 9). Both gaps come from a model that ignores structure: hidden confounders in one case, a curved low-dimensional support in the other.

The course builds probabilistic modelling from one idea: write the joint distribution of everything, observed or hidden, as a product of simple pieces, then compute or approximate what you need from it. Graphs encode the factorisation. Inference is exact on trees (sum-product), and approximate elsewhere (MCMC, variational methods, expectation propagation). Deep generative models replace the simple pieces by neural networks and keep the same inference questions.

It shares tools with Computational Statistics (Monte Carlo), Stochastic Calculus in Machine Learning (diffusions) and Deep Learning. In semester 2 it leads to Bayesian Machine Learning, Graphical Models: Discrete Inference and Learning and Generative Models for Images.

What you will be able to do

• Derive maximum-likelihood estimators for linear, logistic and Gaussian models, and compute them by closed form or Newton’s method.
• Derive and run the EM algorithm for a latent-variable model, and prove that it never decreases the likelihood.
• Compute posterior and predictive distributions in Bayesian linear regression and Gaussian-process regression, and select hyperparameters by the evidence.
• Read conditional independences off a directed or undirected graph (d-separation, Hammersley–Clifford).
• Run sum-product, forward–backward and Viterbi, and learn a hidden Markov model with Baum–Welch.
• Design and check MCMC samplers (Metropolis–Hastings, Gibbs, Langevin) and variational approximations (mean field, variational EM, expectation propagation).
• Fit a stochastic block model to a network by variational EM.
• Train a variational auto-encoder, explain amortisation, and tighten its bound by importance weighting.
• Compare autoregressive models, normalizing flows, GANs, energy-based models and score-based diffusion models by what they can compute.
• Estimate an average treatment effect under ignorability with inverse-propensity weighting and a doubly robust estimator, and state what policy learning optimises.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com presents the course as a unifying introduction to probabilistic modelling through graphical models, with their learning and inference algorithms. It calls the course historical and at the core of the MVA programme, and states that recent developments are largely covered: deep latent variable models, variational auto-encoders, energy-based models, score-based diffusion models, amortised inference and causal inference.

The master-mva.com page lists the topics below (“Thèmes abordés”); the Université Paris-Saclay course detail lists no programme beyond the assessment type. Each row points to the section where the topic is taught.

Beyond the syllabus, the page adds the review of ch. 0, Langevin dynamics as the bridge between MCMC and diffusion models (ch. 6.4), and a self-study chapter on research directions (ch. 12).

Where this course fits in the programme

The official prerequisite is a course on probability. The probability and statistics refresher covers it, including conditional expectation and the multivariate Gaussian. The linear algebra and optimisation refresher covers eigendecompositions, Schur complements and Newton’s method. The scientific Python refresher prepares the NumPy and PyTorch code.

• Related courses in semester 1 — Introduction to Statistical Learning (risk, regression, model selection), Computational Statistics (Monte Carlo, MCMC), Convex Optimization (Newton, duality), Stochastic Calculus in Machine Learning (Langevin, diffusion), Deep Learning, Foundations of Large Language Models (autoregressive models), Optimal Transport for Machine Learning.
• Courses that build on it in semester 2 — Bayesian Machine Learning, Graphical Models: Discrete Inference and Learning, Generative Models for Images, Diffusion and Multiscale Models in AI, Graphs in Machine Learning, Kernel Methods (Gaussian processes), Algorithms for Speech and NLP (HMMs), Sequential Learning and Biostatistics (treatment effects).

The instructors listed on master-mva.com are Pierre Latouche (UCA), Pierre-Alexandre Mattei (Inria) and Mathieu Even (Inria).

Course outline

• Review: probability tools — conditional independence, Gaussian conditioning, Kullback–Leibler divergence, exponential families.
• Maximum likelihood, linear models and K-means — the estimation principle and its first four instances.
• Latent variables — Gaussian mixtures, the EM algorithm and its monotonicity, probabilistic PCA.
• Bayesian linear models, Gaussian processes and model selection — posteriors in closed form, EM revisited, evidence, BIC.
• Directed and undirected graphical models — factorisations, d-separation, Markov properties, Hammersley–Clifford.
• Exact inference — variable elimination, sum-product on trees, forward–backward, Viterbi, Baum–Welch.
• Approximate inference I — importance sampling, Metropolis–Hastings, Gibbs and Langevin samplers.
• Approximate inference II — mean field, variational EM for stochastic block models, expectation propagation.
• Deep latent variable models — VAEs, reparameterisation, amortisation, importance-weighted bounds.
• Beyond VAEs — autoregressive models, normalizing flows, GANs.
• Energy-based and score-based models — contrastive divergence, score matching, diffusion models.
• Causal inference — potential outcomes, identification, IPW and doubly robust estimators, policy learning.
• Research directions — tighter bounds, missing data, diffusion theory, network models and causal machine learning (self-study).

Session plan

30 h in class = 10 sessions of 3 h.

The volume is the Université Paris-Saclay programme table’s (30 h of lectures, 5 ECTS). The master-mva.com page announces “9 lectures of 3 hours each”, on Zoom and recorded, with sessions at ENS Paris-Saclay for questions. The split into ten sessions below is a proposal.

Level 1 = foundations, 2 = core models and algorithms, 3 = deep generative models, diffusion and causal inference.

Self-study, outside the 30 h (about 5 h in total): ch. 12 Research directions and open problems (≈ 3 h) · Ex. 16–17, the two open-ended computational exercises (≈ 2 h).

How to study this module

• Read ch. 0 before the first session, even if the tools look familiar: Gaussian conditioning and the KL divergence return in every later chapter.
• After sessions 2–4, run the code of ch. 1–3 and change one thing (initialisation, noise level, kernel length-scale) before doing the exercises.
• After session 6, redo the HMM example of ch. 5 by hand (Exercise 5) before the approximate-inference sessions (7–8): every later algorithm generalises one of its two passes.
• For the deep-generative sessions (8–10), the scientific Python refresher is enough; all models train on a CPU in under a minute.
• Keep the ELBO decomposition (Key result, ch. 2) in view: EM, variational EM, VAEs and IWAE are four ways of closing the same gap.

0. Review: probability tools for probabilistic models

Every later chapter manipulates joint densities, conditions them on observations and compares two distributions. This review fixes the four tools used throughout: conditional independence, the Gaussian, the Kullback–Leibler divergence and exponential families. The probability refresher gives the background.

0.1 Conditional independence

Graphical models (ch. 4) are nothing but a bookkeeping of such factorisations.

0.2 The multivariate Gaussian

A useful corollary, the linear-Gaussian model: if z \sim \mathcal N(m, P) and x \mid z \sim \mathcal N(Az + b, R), then x \sim \mathcal N(Am + b,\; R + APA^\top).

The posterior of z given x is also Gaussian, by Proposition 0.3 applied to the pair (z, x). Probabilistic PCA, Bayesian regression and Gaussian processes (ch. 2–3) are three instances.

0.3 The Kullback–Leibler divergence

The KL divergence is not symmetric. Minimising it in q (variational inference, ch. 7) and in p (maximum likelihood, ch. 1; expectation propagation, ch. 7.4) give different answers.

0.4 Exponential families

1. Maximum likelihood, linear models and K-means

Given data and a parametric family, which parameter should we pick? Maximum likelihood answers in one line and covers regression, classification and, after a limit, clustering. This chapter states the principle, then applies it to linear regression, logistic regression and K-means.

1.1 Maximum likelihood

Maximising \ell_n / n is minimising an empirical version of \mathrm{KL}(p^\star \| p_\theta), where p^\star is the data distribution. Indeed, by the law of large numbers: \frac1n \ell_n(\theta) \;\to\; \mathbb E_{p^\star}[\log p_\theta(X)] = -\mathrm{KL}(p^\star \| p_\theta) - H(p^\star).

1.2 Linear regression

The model is y_i = w^\top x_i + \varepsilon_i with \varepsilon_i \sim \mathcal N(0, \sigma^2) i.i.d. Stack the inputs in X \in \mathbb R^{n \times d} and the outputs in y \in \mathbb R^n.

Geometrically, X\hat w is the orthogonal projection of y onto the column space of X. Laplace noise would give least absolute deviations instead: the loss is the negative log-likelihood of the noise model.

1.3 Logistic regression

If the two classes are linearly separable, \ell_n increases along the separating direction without bound, and no MLE exists (Exercise 7). A Gaussian prior on w, that is a ridge penalty, restores existence (ch. 3).

import numpy as np
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(0)
n, d = 500, 3
X = np.c_[np.ones(n), rng.normal(size=(n, d - 1))]   # intercept + 2 features
w_true = np.array([-0.5, 2.0, -1.0])
y = rng.random(n) < 1 / (1 + np.exp(-X @ w_true))

w = np.zeros(d)
for it in range(7):                       # Newton = IRLS
    p = 1 / (1 + np.exp(-X @ w))
    grad = X.T @ (y - p)                  # gradient of the log-likelihood
    H = X.T @ (X * (p * (1 - p))[:, None])
    w = w + np.linalg.solve(H, grad)
    print(it, np.round(w, 4), f"|grad| = {np.linalg.norm(grad):.2e}")

sk = LogisticRegression(C=np.inf, fit_intercept=False, tol=1e-10, max_iter=1000).fit(X, y)
print("scikit-learn:", np.round(sk.coef_.ravel(), 4))

The gradient norm falls from 150 to 10^{-14} in six steps, the quadratic convergence of Newton’s method. The estimate (-0.449, 1.638, -1.103) agrees with scikit-learn’s unpenalised fit.

1.4 K-means

The fixed point is only a local minimum: minimising J exactly is NP-hard in general. In practice one uses several restarts and k-means++ seeding (Arthur and Vassilvitskii, 2007). Ch. 2.3 shows K-means is the zero-variance limit of EM for a Gaussian mixture.

2. Latent variables: EM, Gaussian mixtures and probabilistic PCA

Many models are simple only once an unobserved variable is added: the cluster of a point, the low-dimensional cause of a high-dimensional vector. The likelihood of the observed data then contains a sum or an integral over the hidden variable, and has no closed-form maximiser. The EM algorithm turns this into a sequence of complete-data problems.

2.1 Gaussian mixtures

The logarithm of a sum has no closed-form maximiser. Two pathologies deserve a warning.

• Unbounded likelihood: put \mu_k = x_1 and let \Sigma_k \to 0; then \ell_n \to +\infty. The MLE is a local maximum away from these singularities, or one adds a small floor or a prior on \Sigma_k.
• Label switching: permuting the components leaves p(x; \theta) unchanged, so the parameter is identifiable only up to a permutation.

2.2 The EM algorithm

For any distribution q on the latent variables, define the evidence lower bound (ELBO): \mathcal L(q, \theta) = \mathbb E_{q(z)}\big[\log p(x, z; \theta)\big] - \mathbb E_{q(z)}\big[\log q(z)\big].

EM (Dempster, Laird and Rubin, 1977) alternates two steps from \theta_0:

1. E-step: q_{t+1}(z) = p(z \mid x; \theta_t), which closes the gap;
2. M-step: maximise the expected complete log-likelihood, \theta_{t+1} \in \arg\max_\theta Q(\theta, \theta_t), \qquad Q(\theta, \theta_t) = \mathbb E_{p(z \mid x; \theta_t)}\big[\log p(x, z; \theta)\big].

Monotonicity does not give convergence to a global maximum. Under smoothness conditions, limit points are stationary points of the likelihood (Wu, 1983). In practice EM is run from several initialisations, and the best final likelihood is kept.

2.3 EM for Gaussian mixtures

import numpy as np
from scipy.stats import multivariate_normal as mvn

rng = np.random.default_rng(1)
X = np.r_[rng.normal([0, 0], 1.0, (200, 2)), rng.normal([4, 3], 0.7, (150, 2))]
n, K = len(X), 2
pi, mu = np.full(K, 1 / K), X[rng.choice(n, K, replace=False)]
Sigma = np.array([np.eye(2)] * K)

trace = []
for it in range(30):
    dens = np.column_stack([pi[k] * mvn.pdf(X, mu[k], Sigma[k]) for k in range(K)])
    trace.append(np.log(dens.sum(1)).sum())          # log-likelihood at current theta
    tau = dens / dens.sum(1, keepdims=True)          # E-step: responsibilities
    Nk = tau.sum(0)                                  # M-step
    pi, mu = Nk / n, (tau.T @ X) / Nk[:, None]
    for k in range(K):
        R = X - mu[k]
        Sigma[k] = (tau[:, k, None] * R).T @ R / Nk[k]

print("log-likelihood:", np.round(trace[::6], 2))
print("never decreases:", np.all(np.diff(trace) >= -1e-9))
print("weights", np.round(pi, 3), " means", np.round(mu, 2).tolist())

The log-likelihood rises from −3456 to −1094.66 and never decreases, as Theorem 2.3 predicts. The estimated weights are (0.57, 0.43) for true proportions (0.57, 0.43), and the means are close to (0, 0) and (4, 3).

2.4 Probabilistic PCA

As \sigma^2 \to 0, the posterior mean projects x onto the principal subspace: classical PCA is the noiseless limit. EM for PPCA costs O(ndq) per iteration instead of a full eigendecomposition, and handles missing coordinates naturally.

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(2)
n, d, q = 2000, 5, 2
W_true = rng.normal(size=(d, q))
X = rng.normal(size=(n, q)) @ W_true.T + 0.3 * rng.normal(size=(n, d)) + 1.0

S = np.cov(X, rowvar=False, bias=True)          # ML covariance (divides by n)
lam, U = np.linalg.eigh(S)
lam, U = lam[::-1], U[:, ::-1]                   # decreasing eigenvalues
sigma2 = lam[q:].mean()                          # Theorem 2.7: discarded variance
W = U[:, :q] * np.sqrt(lam[:q] - sigma2)         # rotation R = I

C = W @ W.T + sigma2 * np.eye(d)
ll = -0.5 * n * (d * np.log(2 * np.pi) + np.linalg.slogdet(C)[1]
                 + np.trace(np.linalg.solve(C, S)))
pca = PCA(n_components=q).fit(X)
print(f"sigma2 = {sigma2:.4f} (true 0.09); scikit-learn: {pca.noise_variance_:.4f}")
print(f"average log-likelihood {ll / n:.4f}; scikit-learn PCA.score: {pca.score(X):.4f}")

The estimate \hat\sigma^2 = 0.0899 recovers the true noise variance 0.09. The average log-likelihood, -5.1007, equals scikit-learn’s PCA.score, which implements the PPCA likelihood.

3. Bayesian linear models, Gaussian processes and model selection

Maximum likelihood returns one parameter and no uncertainty, and it always prefers the most complex model. A prior on the parameters gives a posterior, a predictive distribution with error bars and a marginal likelihood to compare models. Gaussian priors keep every computation in closed form.

3.1 Bayesian linear regression

The model takes a prior w \sim \mathcal N(0, \alpha^{-1} I_d) and a likelihood y \mid w \sim \mathcal N(Xw, \beta^{-1} I_n), where \alpha and \beta are precisions.

The posterior mean is the ridge estimator with penalty \lambda = \alpha / \beta, which is also the MAP estimate. The predictive variance has two terms: the noise \beta^{-1}, and the parameter uncertainty x_\ast^\top S_N x_\ast, which shrinks as data accumulate.

3.2 Gaussian processes

Replace x by features \phi(x). Under the prior of ch. 3.1, the function f(x) = w^\top \phi(x) has Gaussian finite-dimensional marginals with covariance \alpha^{-1}\phi(x)^\top \phi(x'). A Gaussian process specifies this covariance directly, through a kernel.

The evidence balances a data-fit term against a complexity term, the log-determinant. Hyperparameters such as \ell are chosen by maximising it (type-II maximum likelihood). Exact inference costs O(n^3) through a Cholesky factorisation.

import numpy as np

rng = np.random.default_rng(3)
x = np.sort(rng.uniform(-3, 3, 20))
y = np.sin(2 * x) + 0.1 * rng.normal(size=20)
xs = np.linspace(-3, 3, 7)

def k(a, b, ell):                       # squared-exponential kernel
    return np.exp(-0.5 * (a[:, None] - b[None, :]) ** 2 / ell ** 2)

def gp(ell, s2=0.01):
    Ky = k(x, x, ell) + s2 * np.eye(len(x))
    L = np.linalg.cholesky(Ky)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    mean = k(xs, x, ell) @ alpha
    V = np.linalg.solve(L, k(x, xs, ell))
    var = 1.0 - (V ** 2).sum(0)          # k(x*, x*) = 1
    logev = -0.5 * y @ alpha - np.log(np.diag(L)).sum() - 0.5 * len(x) * np.log(2 * np.pi)
    return mean, var, logev

for ell in [0.1, 0.5, 2.0]:
    m, v, le = gp(ell)
    print(f"ell = {ell:3.1f}   log evidence = {le:8.2f}   max posterior sd = {np.sqrt(v.max()):.2f}")
m, v, _ = gp(0.5)
print("x      ", xs.round(1))
print("mean   ", m.round(2))
print("sd     ", np.sqrt(v).round(2))
print("sin(2x)", np.sin(2 * xs).round(2))

The log evidence is −17.7 for a length-scale of 0.1 (too flexible), −6.0 for 0.5 and −97.6 for 2 (too rigid). With length-scale 0.5, the posterior standard deviation is between 0.08 and 0.22 at the test inputs from −2 to 2, where the data lie. It reaches 0.48 at x = 3 and 0.66 at x = -3, beyond the extreme training inputs −2.49 and 2.74.

3.3 EM revisited

The proof of Theorem 2.3 used only the function \mathcal L(q, \theta). Neal and Hinton (1998) call \mathcal L the (negative) free energy and read EM as coordinate ascent on it.

This view justifies several variants, since any step that increases \mathcal L keeps the monotonicity.

• Generalised EM: the M-step only increases \mathcal L, for example by one gradient step.
• Incremental EM: the E-step updates the posterior of one data point at a time.
• MAP-EM: adding \log p(\theta) to \mathcal L gives the posterior mode instead of the MLE.
• Variational EM: q is restricted to a tractable family, and the bound is no longer tight (ch. 7.2).

EM for the evidence of Bayesian regression. In ch. 3.1, treat w as the latent variable and \theta = (\alpha, \beta) as the parameter. The E-step is Proposition 3.1, q(w) = \mathcal N(m_N, S_N). With d weights and n observations, maximising the expected complete log-likelihood gives \alpha \leftarrow \frac{d}{\|m_N\|^2 + \mathrm{tr}\, S_N}, \qquad \beta \leftarrow \frac{n}{\|y - Xm_N\|^2 + \mathrm{tr}(X S_N X^\top)} .

import numpy as np
from scipy.stats import multivariate_normal as mvn
from sklearn.mixture import GaussianMixture

# EM for the hyperparameters (alpha, beta) of Bayesian linear regression
rng = np.random.default_rng(4)
N, M = 30, 10
x = rng.uniform(0, 1, N)
y = np.sin(2 * np.pi * x) + 0.2 * rng.normal(size=N)        # true beta = 25
Phi = np.vander(x, M, increasing=True)                       # degree-9 polynomial
alpha, beta = 1.0, 1.0
for it in range(200):
    S = np.linalg.inv(alpha * np.eye(M) + beta * Phi.T @ Phi)  # E-step: q(w) = N(m, S)
    m = beta * S @ Phi.T @ y
    alpha = M / (m @ m + np.trace(S))                          # M-step
    beta = N / (np.sum((y - Phi @ m) ** 2) + np.trace(Phi @ S @ Phi.T))
    if it in (0, 9, 199):
        C = Phi @ Phi.T / alpha + np.eye(N) / beta
        print(f"iter {it + 1:3d}: alpha {alpha:7.3f}  beta {beta:6.2f}"
              f"  log evidence {mvn.logpdf(y, np.zeros(N), C):7.3f}")

# Model selection: BIC for the number of mixture components
X = np.r_[rng.normal(0, 1, (150, 2)), rng.normal([5, 0], 1, (150, 2)), rng.normal([2, 4], 1, (150, 2))]
bic = [GaussianMixture(K, n_init=3, random_state=0).fit(X).bic(X) for K in range(1, 7)]
print("BIC for K = 1..6:", np.round(bic, 1), "-> K =", 1 + int(np.argmin(bic)))

The log evidence climbs from -19.7 to -8.7. The noise precision reaches \beta \approx 35, of the order of the true value 25 (the flexible fit absorbs part of the noise), and \alpha falls to 0.01. A degree-9 polynomial fitted to 30 noisy points is thus regularised without any validation set.

3.4 Model selection

Comparing models \mathcal M_1, \mathcal M_2, \dots by their maximised likelihood always favours the largest. Three corrections are standard.

• Evidence: the Bayesian answer compares the marginal likelihoods p(\mathcal D \mid \mathcal M), which penalise unused flexibility automatically (ch. 3.2–3.3).

  p(\mathcal D \mid \mathcal M) = \int p(\mathcal D \mid \theta, \mathcal M)\, p(\theta \mid \mathcal M)\, d\theta .

• BIC: a large-sample approximation of the log evidence (Schwarz, 1978).

• AIC and cross-validation: estimates of the expected predictive log-likelihood (Akaike, 1974), and its direct empirical estimate.

AIC replaces k \log n by 2k. It targets prediction, while BIC is consistent for selecting the true model in a finite list. In the code above, BIC over Gaussian mixtures with K = 1, \dots, 6 components is minimal at K = 3, the true number. A full-covariance mixture in dimension 2 has 6K - 1 free parameters.

4. Directed and undirected graphical models

A joint distribution of n binary variables has 2^n - 1 free parameters, far too many to store or learn. Graphical models restrict the joint to a product of local factors, and the graph tells which conditional independences follow. Directed graphs encode generative stories; undirected graphs encode symmetric interactions.

4.1 Directed graphical models

For binary variables the factorisation needs \sum_i 2^{|\pi_i|} parameters instead of 2^n - 1. The models of ch. 1–3 are all directed graphical models.

• Naive Bayes: a class y with arrows to each feature x_j, so the features are independent given the class.
• Gaussian mixture: z_i \to x_i, repeated over i; a plate drawn around the pair stands for the n copies.
• Bayesian regression: w \to y_i \leftarrow x_i, with the plate over i.
• Hidden Markov model (ch. 5.3): a chain z_1 \to z_2 \to \dots \to z_T, with an arrow z_t \to x_t at each step.

4.2 d-separation and explaining away

Three structures on three nodes a, c, b contain all the cases.

• Chain a \to c \to b: the joint is p(a)\,p(c \mid a)\,p(b \mid c), of the form f(a, c)\, g(b, c). So a \perp b \mid c by Proposition 0.2.
• Fork a \leftarrow c \to b: the joint p(c)\,p(a \mid c)\,p(b \mid c) has the same form, so again a \perp b \mid c.
• Collider a \to c \leftarrow b: summing over c gives p(a, b) = p(a)\, p(b), so a \perp b. Conditioning on c (or on a descendant of c) creates a dependence.

4.3 Undirected graphical models

The normalising constant Z, the partition function, is a sum over exponentially many states. It is the main computational obstacle of undirected models, and it returns with energy-based models in ch. 10.

Three Markov properties can hold with respect to G, from weakest to strongest:

• pairwise: two non-adjacent variables are independent given all the others;
• local: a variable is independent of the rest given its neighbours;
• global: x_A \perp x_B \mid x_S whenever S separates A from B in G.

4.4 Hammersley–Clifford, Ising and Gaussian graphical models

Positivity is needed: there are distributions with zeros that satisfy the Markov properties without factorising.

The Ising model. For spins x \in \{-1, +1\}^n on a graph: p(x) = \frac1Z \exp\Big(\sum_i h_i x_i + \sum_{(i,j) \in E} J_{ij}\, x_i x_j\Big).

It is the prototype of a pairwise Markov random field. Computing its Z on a general graph is #P-hard, which motivates the samplers of ch. 6 (Exercise 10).

So the graph of a Gaussian model is read off the zeros of its precision matrix. Estimating a sparse precision matrix with an \ell_1 penalty is the graphical lasso (Friedman, Hastie and Tibshirani, 2008).

4.5 From directed to undirected: moralisation

Each factor p(x_i \mid x_{\pi_i}) is a potential on the set \{i\} \cup \pi_i, which is a clique of the moral graph. So a directed factorisation is a Gibbs factorisation on the moral graph, with Z = 1. The converse fails: moralisation hides the marginal independence of the two parents of a collider.

5. Exact inference: the sum-product algorithm and hidden Markov models

Once a model is written, every question about it is a sum: a marginal, a conditional given evidence, the partition function. Done naively, these sums have exponentially many terms. On trees, the distributive law reduces them to messages passed along the edges, at linear cost.

5.1 Variable elimination and factor graphs

Take a chain x_1 - x_2 - \dots - x_n with K states each and pairwise potentials. The marginal of x_n is a sum over K^{n-1} configurations. Pushing each sum inside the product removes one variable at a time: \sum_{x_1, \dots, x_{n-1}} \prod_{t=1}^{n-1} \psi_t(x_t, x_{t+1}) = \sum_{x_{n-1}} \psi_{n-1}(x_{n-1}, x_n) \cdots \sum_{x_1} \psi_1(x_1, x_2), which costs O(nK^2) operations instead of O(K^n).

Directed and undirected models both convert to factor graphs. On a general graph, the cost of elimination is exponential in the size of the largest clique created during elimination, which depends on the elimination order. The best order makes this size the treewidth plus one. Exact inference is NP-hard in general (Cooper, 1990).

5.2 Sum-product and max-product on trees

The sum-product algorithm sends messages along the edges of a factor graph. A variable sends to a factor the product of its other incoming messages. A factor sends to a variable the sum, over its other variables, of its potential times their incoming messages: m_{i \to f}(x_i) = \prod_{g \ni i,\, g \ne f} m_{g \to i}(x_i), \qquad m_{f \to i}(x_i) = \sum_{x_{f \setminus i}} \psi_f(x_f) \prod_{j \in f,\, j \ne i} m_{j \to f}(x_j).

Replacing sums by maxima gives max-product. It computes max-marginals, and backtracking the argmaxes returns a most probable configuration (MAP). On graphs with cycles, the same updates define loopy belief propagation, which often works well but has no general guarantee (Murphy, Weiss and Jordan, 1999).

5.3 Hidden Markov models: forward–backward and Viterbi

The HMM factor graph is a chain, so sum-product applies. Its two passes are the forward and backward recursions, for \alpha_t(k) = p(x_{1:t}, z_t = k) and \beta_t(k) = p(x_{t+1:T} \mid z_t = k): \alpha_{t+1}(k) = b_k(x_{t+1}) \sum_j \alpha_t(j) A_{jk}, \qquad \beta_t(j) = \sum_k A_{jk}\, b_k(x_{t+1})\, \beta_{t+1}(k).

The cost is O(TK^2). The \alpha_t underflow for long sequences, so one normalises them at each step and keeps the normalisers c_t; then \log p(x_{1:T}) = \sum_t \log c_t. The Viterbi algorithm is max-product in the log domain. When the states and observations are Gaussian with linear dynamics, the same recursions become the Kalman filter and smoother.

import itertools
import numpy as np

A = np.array([[0.9, 0.1], [0.2, 0.8]])        # transitions  A[i, j] = p(z_t+1 = j | z_t = i)
B = np.array([[0.7, 0.2, 0.1], [0.1, 0.3, 0.6]])  # emissions  B[i, o] = p(x_t = o | z_t = i)
pi0 = np.array([0.5, 0.5])
x = [0, 0, 2, 2, 1, 2]
T, K = len(x), 2

alpha, c = np.zeros((T, K)), np.zeros(T)       # scaled forward pass
for t in range(T):
    a = pi0 * B[:, x[t]] if t == 0 else (alpha[t - 1] @ A) * B[:, x[t]]
    c[t] = a.sum(); alpha[t] = a / c[t]
beta = np.ones((T, K))                         # scaled backward pass
for t in range(T - 2, -1, -1):
    beta[t] = A @ (B[:, x[t + 1]] * beta[t + 1]) / c[t + 1]
gamma = alpha * beta                           # p(z_t | x_1:T)

# brute force over the K**T hidden paths
joint = {z: pi0[z[0]] * np.prod([A[z[t - 1], z[t]] for t in range(1, T)])
            * np.prod([B[z[t], x[t]] for t in range(T)])
         for z in itertools.product(range(K), repeat=T)}
px = sum(joint.values())
g_bf = np.array([[sum(p for z, p in joint.items() if z[t] == k) / px for k in range(K)] for t in range(T)])
print(f"log p(x): forward {np.log(c).sum():.6f}   brute force {np.log(px):.6f}")
print("max |gamma - brute force| =", np.abs(gamma - g_bf).max())

delta, back = np.log(pi0 * B[:, x[0]]), []    # Viterbi (max-product in logs)
for t in range(1, T):
    s = delta[:, None] + np.log(A)
    back.append(s.argmax(0)); delta = s.max(0) + np.log(B[:, x[t]])
path = [int(delta.argmax())]
for bp in reversed(back):
    path.append(int(bp[path[-1]]))
print("Viterbi path:", path[::-1], "  brute force:", list(max(joint, key=joint.get)))

The forward pass gives \log p(x_{1:6}) = -6.651829, equal to the sum over the 2^6 = 64 hidden paths. The posteriors agree to 10^{-15}, and Viterbi returns the path (0, 0, 1, 1, 1, 1) found by enumeration.

5.4 Learning an HMM: Baum–Welch

Baum–Welch is EM (ch. 2.2) for the HMM. The E-step computes \gamma and \xi by forward–backward. As for mixtures (Proposition 2.4), the M-step replaces counts of states and transitions by their expectations: \pi_k = \gamma_1(k), \qquad A_{jk} = \frac{\sum_{t=1}^{T-1} \xi_t(j, k)}{\sum_{t=1}^{T-1} \gamma_t(j)}, \qquad b_k(o) = \frac{\sum_t \gamma_t(k)\, \mathbb 1[x_t = o]}{\sum_t \gamma_t(k)} .

With Gaussian emissions, the last update becomes a \gamma-weighted mean and covariance. By Theorem 2.3, each iteration does not decrease the likelihood. As for mixtures, the result depends on the initialisation, and several restarts are standard.

6. Approximate inference I: Monte Carlo and MCMC

When the graph has cycles or the variables are continuous and non-Gaussian, exact sums are out of reach. Monte Carlo methods replace them by averages over samples. Markov chain Monte Carlo (MCMC) produces such samples from a density known only up to its normalising constant, which is exactly what posteriors and Gibbs distributions provide.

6.1 Monte Carlo and importance sampling

For x_1, \dots, x_n i.i.d. from p, the average of f(x_i) is unbiased for \mathbb E_p[f], with variance \mathrm{Var}_p(f)/n, whatever the dimension.

When the target is hard to sample, draw from a proposal q and reweight.

The quality of the weights is summarised by the effective sample size (\sum_i w_i)^2 / \sum_i w_i^2. In high dimension, the weights of a fixed proposal degenerate: a few samples carry all the mass. Importance sampling returns in ch. 8.3, inside the IWAE bound.

6.2 Metropolis–Hastings

Invariance is not convergence. If the chain is also irreducible and aperiodic, the law of x_t converges to \pi, and ergodic averages satisfy a law of large numbers (Robert and Casella, 2004). For a random-walk proposal in high dimension, a classical result recommends tuning the step to an acceptance rate near 0.234 (Roberts, Gelman and Gilks, 1997).

6.3 Gibbs sampling

The Gibbs sampler updates one variable at a time from its full conditional p(x_i \mid x_{-i}), cycling over i.

In a graphical model, the full conditional depends only on the Markov blanket: the neighbours in an undirected graph, or the parents, children and co-parents in a DAG. For the Ising model of ch. 4.4, with \sigma the logistic function: p(x_i = +1 \mid x_{-i}) = \sigma\Big(2\big(h_i + \textstyle\sum_{j} J_{ij} x_j\big)\Big).

Gibbs mixes slowly when variables are strongly coupled: in a highly correlated Gaussian, coordinate moves are tiny. Blocking correlated variables or reparameterising helps.

6.4 Langevin dynamics

Write the target as \pi \propto e^{-U}. The overdamped Langevin diffusion follows the gradient of \log \pi plus noise: dX_t = -\nabla U(X_t)\, dt + \sqrt 2\, dW_t .

The Euler discretisation with step h is the unadjusted Langevin algorithm (ULA): x_{k+1} = x_k - h\, \nabla U(x_k) + \sqrt{2h}\; \xi_k, \qquad \xi_k \sim \mathcal N(0, I).

ULA is biased: its stationary law differs from \pi by an amount controlled by h (Durmus and Moulines, 2017). Using the ULA move as the proposal of a Metropolis–Hastings step removes the bias: this is MALA.

import numpy as np

rng = np.random.default_rng(5)
U = lambda x: 0.5 * x ** 2            # target pi = N(0, 1), potential U = -log pi
gradU = lambda x: x

def run(h, n=200_000, adjust=False):
    x, xs, acc = 0.0, np.empty(n), 0
    for i in range(n):
        y = x - h * gradU(x) + np.sqrt(2 * h) * rng.normal()      # Langevin proposal
        if adjust:                                                 # MALA: MH correction
            logq = lambda a, b: -(a - b + h * gradU(b)) ** 2 / (4 * h)
            if np.log(rng.random()) < U(x) - U(y) + logq(x, y) - logq(y, x):
                x, acc = y, acc + 1
        else:
            x = y                                                  # ULA: always accept
        xs[i] = x
    return xs.var(), acc / n

for h in [0.1, 0.5]:
    v_ula, _ = run(h)
    v_mala, a = run(h, adjust=True)
    print(f"h = {h}: ULA variance {v_ula:.3f} (theory {1 / (1 - h / 2):.3f})"
          f" | MALA variance {v_mala:.3f}, acceptance {a:.2f}")

With h = 0.5, ULA’s empirical variance is 1.338, as Example 6.6 predicts, while MALA gives 0.995 with 92 % acceptance. Langevin dynamics need only the gradient of \log \pi, the score. Score-based diffusion models (ch. 10.3) run it with a learned score.

7. Approximate inference II: variational methods, stochastic block models and expectation propagation

Variational methods turn inference into optimisation: choose, in a tractable family, the distribution closest to the posterior. They are deterministic and fast, and they plug directly into EM. This chapter covers mean field, variational EM, its application to stochastic block models for networks, and expectation propagation, which minimises the KL divergence in the other direction.

7.1 Variational inference and mean field

By Proposition 2.2, maximising the ELBO \mathcal L(q) over a family \mathcal Q is minimising \mathrm{KL}(q \,\|\, p(z \mid x)) over \mathcal Q. The mean-field family makes the latent coordinates independent: q(z) = \prod_{j=1}^m q_j(z_j).

Cycling over j never decreases the ELBO. It converges to a stationary point, in practice a local optimum, so the result depends on the initialisation.

Because \mathrm{KL}(q \| p) penalises mass of q where p is small, the solution is mode-seeking: mean field underestimates posterior variances. For a Gaussian target with precision \Lambda, the optimal factor q_j has variance 1/\Lambda_{jj}, never larger than the true marginal variance (Exercise 11).

7.2 Variational EM

When the E-step of EM is intractable, restrict q to \mathcal Q and alternate:

1. VE-step: maximise \mathcal L(\cdot, \theta_t) over q \in \mathcal Q, for example by CAVI, giving q_{t+1};
2. M-step: maximise \mathcal L(q_{t+1}, \cdot) over \theta, giving \theta_{t+1}.

By Proposition 3.5, this is coordinate ascent on the free energy over a smaller set. The ELBO increases at each step. The likelihood itself need not increase, because the gap \mathrm{KL}(q \| p(z \mid x; \theta)) no longer vanishes.

7.3 Stochastic block models

Diagonal-dominant \pi gives communities, but the model also captures hubs, bipartite and core–periphery structures. Exact EM fails: every edge X_{ij} is a common child of z_i and z_j, so moralisation couples all the labels. The posterior p(z \mid X) does not factorise, and its normaliser sums over K^n labellings.

The number of classes is chosen by the integrated classification likelihood (ICL). For an undirected graph, Daudin et al. (2008) penalise the complete-data log-likelihood, evaluated at the estimated labels, by \frac{K - 1}{2} \log n + \frac{K(K+1)}{4} \log \frac{n(n-1)}{2} .

The first term counts the K - 1 free proportions; the second counts the K(K+1)/2 connection probabilities, estimated from n(n-1)/2 pairs.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score as ari

rng = np.random.default_rng(6)
n, K = 300, 3
z = rng.choice(K, n, p=[0.5, 0.3, 0.2])
P = np.full((K, K), 0.04) + 0.10 * np.eye(K)                 # assortative SBM
A = np.triu(rng.random((n, n)) < P[z][:, z], 1).astype(float)
A = A + A.T                                                   # undirected, no self-loop

U = np.linalg.eigh(A)[1][:, -K:]                              # spectral initialisation
z0 = KMeans(K, n_init=10, random_state=0).fit_predict(U)
tau = np.eye(K)[z0] * 0.9 + 0.1 / K
off = 1 - A - np.eye(n)                                       # non-edges, j != i
for it in range(50):
    alpha = tau.mean(0)                                       # M-step
    Pi = (tau.T @ A @ tau) / (tau.T @ (1 - np.eye(n)) @ tau)
    Pi = Pi.clip(1e-6, 1 - 1e-6)
    for _ in range(5):                                        # E-step: fixed point on tau
        S = np.log(alpha) + A @ tau @ np.log(Pi).T + off @ tau @ np.log(1 - Pi).T
        tau = np.exp(S - S.max(1, keepdims=True))
        tau /= tau.sum(1, keepdims=True)

print("ARI spectral init:", round(ari(z, z0), 3))
print("ARI after VEM:    ", round(ari(z, tau.argmax(1)), 3))
print("estimated alpha", np.round(np.sort(alpha), 2), "\nestimated Pi\n", np.round(Pi, 3))

Spectral clustering of the adjacency matrix gives an adjusted Rand index of 0.86. Variational EM, started from it, reaches 0.975, and recovers the proportions (0.23, 0.29, 0.48) and connection probabilities within 0.005 of the truth (0.14 inside, 0.04 between).

7.4 Expectation propagation

Expectation propagation (EP; Minka, 2001) targets posteriors of the form p(\theta \mid \mathcal D) \;\propto\; p_0(\theta) \prod_{i=1}^n t_i(\theta), and approximates each site t_i by a member \tilde t_i of an exponential family. Its building block is a KL projection in the direction opposite to variational inference.

For a Gaussian family, q^\star matches the mean and the covariance of p. EP iterates over the sites:

1. cavity: remove site i, giving q^{\setminus i} \propto q / \tilde t_i;
2. tilted distribution: put back the exact site, giving \hat p_i \propto q^{\setminus i}\, t_i;
3. projection: set q_{\text{new}} to the moment-matched projection of \hat p_i onto \mathcal F (Proposition 7.4);
4. site update: \tilde t_i \propto q_{\text{new}} / q^{\setminus i}.

Each step needs only one-dimensional or low-dimensional integrals against one exact site, which are often available in closed form, as for probit likelihoods. A single pass is assumed density filtering. EP has no general convergence guarantee. Where it converges, it tends to cover the posterior mass better than mean field, which is why it is a standard choice for Gaussian-process classification (Rasmussen and Williams, 2006, chapter 3).

8. Deep latent variable models, variational auto-encoders and amortised inference

Replace the linear map of probabilistic PCA (ch. 2.4) by a neural network, and the model can generate images, sounds or molecules. The price is that neither the likelihood nor the posterior has a closed form any more. This chapter combines the ELBO of ch. 2, an inference network shared by all data points, and a gradient estimator that goes through sampling.

8.1 Deep latent variable models

The integral has no closed form, and neither does the posterior p_\theta(z \mid x), so the E-step of EM is impossible. A further subtlety: with a Gaussian decoder whose variance is also learned, the likelihood of a DLVM can be unbounded, like that of a Gaussian mixture (ch. 2.1). Mattei and Frellsen (2018) analyse this degeneracy.

8.2 Amortised variational inference and the VAE

Variational EM would optimise one distribution q_i(z) per data point. Amortised inference replaces them by a single network, the encoder, that maps x to the parameters of q_\phi(z \mid x), typically a Gaussian with mean \mu_\phi(x) and diagonal variances \sigma^2_\phi(x). The variational auto-encoder (VAE; Kingma and Welling, 2014; Rezende, Mohamed and Wierstra, 2014) maximises over (\theta, \phi) the sum over data points of

The first term is a reconstruction error; the second keeps the codes close to the prior. For a diagonal Gaussian q and a standard normal prior, the KL term is explicit: \mathrm{KL}\big(\mathcal N(\mu, \mathrm{diag}\,\sigma^2) \,\|\, \mathcal N(0, I)\big) = \frac12 \sum_{j=1}^q \big(\mu_j^2 + \sigma_j^2 - 1 - \log \sigma_j^2\big).

A single draw of \varepsilon gives an unbiased gradient of the ELBO, usually with much lower variance than the score-function estimator \mathbb E_q[h(z)\, \nabla_\phi \log q_\phi(z \mid x)], which also works for discrete z.

The gap between \log p_\theta(x) and the ELBO splits in two parts (Cremer, Li and Duvenaud, 2018):

• the approximation gap: even the best Gaussian q for this x differs from the posterior;
• the amortisation gap: the encoder’s output is not that best Gaussian.

A known failure is posterior collapse: with a very flexible decoder, the optimum can ignore z and set q_\phi(z \mid x) = p(z), making the KL term zero.

8.3 Importance-weighted bounds

The ELBO is \mathbb E_q[\log w] with w = p_\theta(x, z) / q_\phi(z \mid x). Averaging several weights inside the logarithm gives a tighter bound.

Training on \mathcal L_K gives the IWAE; even a VAE trained on \mathcal L_1 is best evaluated with a large K. The code trains a small VAE on the 8×8 digits bundled with scikit-learn, binarised, then evaluates \mathcal L_K on held-out images.

import math
import torch
from torch import nn
from sklearn.datasets import load_digits

torch.manual_seed(0); torch.set_num_threads(1)
X = torch.tensor(load_digits().data / 16.0, dtype=torch.float32)   # 1797 images 8x8 in [0, 1]
X = (X > 0.5).float()                                              # binarise: Bernoulli pixels
Xtr, Xte = X[:1500], X[1500:]
d, q, h = 64, 4, 128
enc = nn.Sequential(nn.Linear(d, h), nn.Tanh(), nn.Linear(h, 2 * q))   # amortised q(z|x)
dec = nn.Sequential(nn.Linear(q, h), nn.Tanh(), nn.Linear(h, d))       # logits of p(x|z)
opt = torch.optim.Adam([*enc.parameters(), *dec.parameters()], lr=1e-3)

def log_weights(x, K):
    mu, logvar = enc(x).chunk(2, -1)
    z = mu + torch.exp(0.5 * logvar) * torch.randn(K, *mu.shape)      # reparameterisation
    log_px_z = -nn.functional.binary_cross_entropy_with_logits(
        dec(z), x.expand(K, *x.shape), reduction="none").sum(-1)
    log_pz = -0.5 * (z ** 2 + math.log(2 * math.pi)).sum(-1)
    log_qz = -0.5 * (((z - mu) ** 2) / logvar.exp() + logvar + math.log(2 * math.pi)).sum(-1)
    return log_px_z + log_pz - log_qz                                  # shape (K, batch)

for epoch in range(60):
    for xb in Xtr[torch.randperm(1500)].split(100):
        loss = -log_weights(xb, 1).mean()                              # negative ELBO
        opt.zero_grad(); loss.backward(); opt.step()

with torch.no_grad():
    for K in [1, 10, 100]:
        lw = log_weights(Xte, K)
        bound = (torch.logsumexp(lw, 0) - math.log(K)).mean()          # IWAE bound
        print(f"K = {K:3d}: test bound on log p(x) = {bound:.2f} nats")

On the 297 test images, the bound rises from −19.84 nats (K = 1) to −19.66 (K = 10) and −19.62 (K = 100): the ELBO underestimates the log-likelihood by at least 0.2 nats per image. Training takes a few seconds on a CPU.

9. Deep generative models beyond VAEs: autoregressive models, normalizing flows and GANs

A VAE gives only a bound on the likelihood. Three other families make different trade-offs. Autoregressive models and normalizing flows give the exact likelihood, at a price on sampling or on architecture. GANs give up the likelihood altogether and learn to sample by competing with a classifier.

9.1 Autoregressive models

By the chain rule, any distribution on x = (x_1, \dots, x_D) factorises without assumption: p(x) = \prod_{d=1}^D p(x_d \mid x_{<d}).

An autoregressive model parameterises every conditional with one network, whose weights are shared across positions. As a directed graphical model, it is the complete DAG in a fixed order.

• Likelihood: exact, a sum of D log-conditionals; training is maximum likelihood (teacher forcing).
• Architectures: masks enforce the ordering, in MADE (Germain et al., 2015) and PixelCNN (van den Oord, Kalchbrenner and Kavukcuoglu, 2016); causal transformers do the same for text.
• Sampling: sequential, D network evaluations per sample.

The large language models of Foundations of Large Language Models are autoregressive models of token sequences.

9.2 Normalizing flows

The determinant of a general D \times D Jacobian costs O(D^3). Flows use layers whose Jacobian is triangular.

Alternating which half is transformed lets every coordinate depend on every other after a few layers. Flows keep the dimension and the topology of the base distribution. To send a Gaussian onto two separate moons, a flow must stretch space strongly between them.

import math
import torch
from torch import nn
from sklearn.datasets import make_moons

torch.manual_seed(0); torch.set_num_threads(1)
X = torch.tensor(make_moons(2000, noise=0.05, random_state=0)[0], dtype=torch.float32)
X = (X - X.mean(0)) / X.std(0)

class Coupling(nn.Module):                     # affine coupling: x_b -> x_b * exp(s(x_a)) + t(x_a)
    def __init__(self, swap):
        super().__init__()
        self.swap = swap
        self.net = nn.Sequential(nn.Linear(1, 64), nn.ReLU(), nn.Linear(64, 2))
    def forward(self, x):                      # data -> latent, returns log|det J|
        a, b = (x[:, 1:], x[:, :1]) if self.swap else (x[:, :1], x[:, 1:])
        s, t = self.net(a).chunk(2, -1)
        s = torch.tanh(s)
        b = b * torch.exp(s) + t
        out = torch.cat([b, a], 1) if self.swap else torch.cat([a, b], 1)
        return out, s.sum(-1)

layers = nn.ModuleList([Coupling(i % 2 == 1) for i in range(6)])
opt = torch.optim.Adam(layers.parameters(), lr=3e-3)

def log_prob(x):
    logdet = torch.zeros(len(x))
    for f in layers:
        x, ld = f(x)
        logdet = logdet + ld
    return -0.5 * (x ** 2).sum(-1) - math.log(2 * math.pi) + logdet   # change of variables

for step in range(1500):
    loss = -log_prob(X[torch.randint(0, 2000, (256,))]).mean()
    opt.zero_grad(); loss.backward(); opt.step()
    if step % 500 == 0 or step == 1499:
        with torch.no_grad():                   # evaluate on all 2000 points, not one batch
            print(f"step {step:4d}: average log-likelihood {log_prob(X).mean().item():.3f}")
C = torch.cov(X.T)                              # best Gaussian (maximum likelihood) for comparison
print(f"best Gaussian: average log-likelihood {-0.5 * (2 * (1 + math.log(2 * math.pi)) + torch.logdet(C).item()):.3f}")

On all 2000 points, six coupling layers raise the average log-likelihood from −2.88 after the first step to −1.21 after 1500 steps. The best Gaussian on the same standardised data reaches only −2.73. Each layer here maps data to the latent space, which is the direction the likelihood needs; sampling uses the inverse of Proposition 9.2.

9.3 Generative adversarial networks

A generator G maps noise z \sim p(z) to samples with law p_G. A discriminator D : x \mapsto [0, 1] tries to tell data from samples. The GAN of Goodfellow et al. (2014) solves \min_G \max_D \; V(D, G) = \mathbb E_{x \sim p_{\text{data}}}\big[\log D(x)\big] + \mathbb E_{z \sim p(z)}\big[\log\big(1 - D(G(z))\big)\big].

In practice D is never optimal, and the generator usually maximises \log D(G(z)) instead, which gives stronger gradients early in training. GANs need no likelihood and produce sharp samples. They can also miss modes (mode collapse) and are hard to evaluate. The Wasserstein GAN (Arjovsky, Chintala and Bottou, 2017) replaces the JS divergence by a transport distance.

10. Energy-based models and score-based diffusion models

An energy-based model specifies a density up to its normalising constant, like the undirected models of ch. 4.3, with a neural network as the energy. Its likelihood gradient needs samples from the model, which MCMC provides only slowly. Score matching avoids the constant altogether. Learning the score of noisy versions of the data, then running a reverse diffusion, gives today’s diffusion models.

10.1 Energy-based models

The Ising model (ch. 4.4) is an EBM with an energy linear in its parameters. A neural energy can express any shape, but Z_\theta is intractable.

Averaged over the data, the first term lowers the energy of observed points (positive phase). The second raises the energy of model samples (negative phase), and needs MCMC from p_\theta, typically Langevin (ch. 6.4). Contrastive divergence (Hinton, 2002) starts the chains at the data and runs only a few steps. This is cheap but biased.

10.2 Score matching

The score \nabla_x \log p(x) does not involve the normalising constant. Score matching fits a model score s_\theta(x) to the data score by minimising the Fisher divergence J(\theta) = \frac12\, \mathbb E_{p}\Big[\big\|s_\theta(x) - \nabla_x \log p(x)\big\|^2\Big].

The trace of a Jacobian is expensive for networks in high dimension. Denoising avoids it. Perturb the data with Gaussian noise, \tilde x = x + \sigma \varepsilon, and let p_\sigma be the law of \tilde x.

So the network learns to predict the noise: -\sigma s_\theta(\tilde x) regresses \varepsilon. The minimiser is linked to the best denoiser by Tweedie’s formula (Exercise 14): \mathbb E[x \mid \tilde x] = \tilde x + \sigma^2\, \nabla \log p_\sigma(\tilde x).

10.3 Score-based diffusion models

Diffusion models (Sohl-Dickstein et al., 2015; Ho, Jain and Abbeel, 2020; Song et al., 2021) learn scores at all noise levels at once. A forward SDE carries the data distribution p_0 to a simple law p_T: dx = f(x, t)\, dt + g(t)\, dW_t .

In the variance-exploding case, f = 0 and the noise scale \sigma(t) increases, with g(t)^2 = \frac{d}{dt}\,\sigma(t)^2 .

Then x_t = x_0 + \sigma(t)\,\varepsilon when \sigma(0) = 0, and p_T is nearly \mathcal N(0, \sigma(T)^2 I). On a grid of noise scales, g^2\,dt becomes \sigma_i^2 - \sigma_{i+1}^2, as in the code below.

In the variance-preserving case, f(x, t) = -\frac12 \beta(t)\, x, and p_T is nearly \mathcal N(0, I); its discretisation is DDPM.

The unknown score is replaced by a network s_\theta(x, t), trained by denoising score matching (Theorem 10.4) over noise levels. Sampling integrates the reverse SDE or the ODE. Because the ODE is a continuous normalizing flow (ch. 9.2), it also gives a likelihood.

import math
import torch
from torch import nn

torch.manual_seed(0); torch.set_num_threads(1)
# data: mixture of 8 Gaussians on a circle of radius 2
ang = 2 * math.pi * torch.randint(0, 8, (4000,)) / 8
X = 2 * torch.stack([ang.cos(), ang.sin()], 1) + 0.1 * torch.randn(4000, 2)

sig_min, sig_max = 0.01, 3.0                       # variance-exploding noise scales
net = nn.Sequential(nn.Linear(3, 128), nn.SiLU(), nn.Linear(128, 128), nn.SiLU(), nn.Linear(128, 2))
opt = torch.optim.Adam(net.parameters(), lr=2e-3)
score = lambda x, s: net(torch.cat([x, s.log()], 1)) / s   # s_theta(x, sigma)

for step in range(3000):                           # denoising score matching
    x0 = X[torch.randint(0, 4000, (512,))]
    s = sig_min * (sig_max / sig_min) ** torch.rand(512, 1)
    eps = torch.randn_like(x0)
    loss = ((s * score(x0 + s * eps, s) + eps) ** 2).sum(1).mean()   # sigma^2-weighted DSM
    opt.zero_grad(); loss.backward(); opt.step()

with torch.no_grad():                              # reverse-time SDE, Euler-Maruyama
    N = 500
    sig = sig_min * (sig_max / sig_min) ** torch.linspace(1, 0, N + 1)
    x = sig_max * torch.randn(2000, 2)
    for i in range(N):
        g2 = sig[i] ** 2 - sig[i + 1] ** 2          # g(t)^2 dt for this step
        s = torch.full((2000, 1), sig[i].item())
        x = x + g2 * score(x, s) + g2.sqrt() * torch.randn_like(x)
r = x.norm(dim=1)
print(f"radius of samples: mean {r.mean():.3f}, sd {r.std():.3f} (data: 2.0, ~0.1)")
k = torch.remainder(torch.round(torch.atan2(x[:, 1], x[:, 0]) / (2 * math.pi / 8)), 8)
print("samples per mode:", torch.bincount(k.long(), minlength=8).tolist())

The 2000 samples have mean radius 1.986 and standard deviation 0.115, against 2 and about 0.1 for the data. All eight modes are populated, with between 220 and 347 samples each (250 expected): no mode collapse, but the weights are not perfectly balanced.

11. Causal inference: potential outcomes, treatment effects and policy learning

A predictive model learns p(y \mid x, w); a decision needs to know what would happen if we set the treatment w. The two coincide in a randomised experiment, and differ when treatment depends on variables that also affect the outcome. This chapter states when causal effects are identified from observational data, how to estimate them, and how to learn whom to treat.

11.1 Potential outcomes

The fundamental problem of causal inference is that Y_i(0) and Y_i(1) are never observed together. Causal effects must therefore be identified from assumptions on how treatment was assigned.

11.2 Identification: randomisation and ignorability

In a randomised experiment, W is independent of (Y(0), Y(1)). Then \mathbb E[Y \mid W = 1] = \mathbb E[Y(1)], and the difference of the two group means is unbiased for \tau.

In observational data, write the propensity score e(x) = P(W = 1 \mid X = x) and the outcome regressions \mu_w(x) = \mathbb E[Y \mid X = x, W = w].

In the language of ch. 4, ignorability holds when X blocks every back-door path from W to Y in a causal DAG, as in W \leftarrow X \to Y; the first formula is then Pearl’s back-door adjustment (Pearl, 2009). Rosenbaum and Rubin (1983) show that conditioning on the scalar e(X) suffices.

11.3 Inverse-propensity weighting and doubly robust estimation

Plugging estimates into Theorem 11.2 gives two estimators.

• Inverse-propensity weighting (IPW) averages W_i Y_i / \hat e(X_i) - (1 - W_i) Y_i / (1 - \hat e(X_i)). Its variance explodes when \hat e approaches 0 or 1.
• Regression adjustment averages \hat\mu_1(X_i) - \hat\mu_0(X_i). It is biased when \hat\mu_w is misspecified.

The augmented IPW (AIPW) estimator combines them (Robins, Rotnitzky and Zhao, 1994). It averages the scores \Gamma_i = \hat\mu_1(X_i) - \hat\mu_0(X_i) + \frac{W_i\,(Y_i - \hat\mu_1(X_i))}{\hat e(X_i)} - \frac{(1 - W_i)\,(Y_i - \hat\mu_0(X_i))}{1 - \hat e(X_i)} .

With machine-learning estimates of the nuisances, fitted on separate folds (cross-fitting), AIPW is \sqrt n-consistent and asymptotically normal as soon as both nuisance errors decrease faster than n^{-1/4}: this is double/debiased machine learning (Chernozhukov et al., 2018). The code simulates confounded data with a true ATE of 1.

import numpy as np
from sklearn.linear_model import LogisticRegression, LinearRegression

rng = np.random.default_rng(8)
n = 20_000
X = rng.normal(size=(n, 2))
e = 1 / (1 + np.exp(-(1.5 * X[:, 0] - X[:, 1])))          # propensity: confounded assignment
W = rng.random(n) < e
tau = 1.0 + X[:, 1]                                        # heterogeneous effect, ATE = 1
Y = 2 * X[:, 0] + X[:, 1] + W * tau + rng.normal(size=n)

naive = Y[W].mean() - Y[~W].mean()
e_hat = LogisticRegression().fit(X, W).predict_proba(X)[:, 1]
ipw = np.mean(W * Y / e_hat - (~W) * Y / (1 - e_hat))
m1 = LinearRegression().fit(X[W], Y[W]).predict(X)        # outcome models mu_1, mu_0
m0 = LinearRegression().fit(X[~W], Y[~W]).predict(X)
psi = m1 - m0 + W * (Y - m1) / e_hat - (~W) * (Y - m0) / (1 - e_hat)
aipw, se = psi.mean(), psi.std() / np.sqrt(n)               # AIPW and its standard error
print(f"true ATE 1.000 | naive {naive:.3f} | IPW {ipw:.3f} | AIPW {aipw:.3f} +/- {1.96 * se:.3f}")

The naive difference of means gives 1.946, because treated units have larger X_1, which also raises Y. IPW gives 1.109, unbiased but noisy. AIPW gives 1.002, with a 95 % confidence half-width of 0.054.

11.4 Policy learning

A policy \pi : x \mapsto \{0, 1\} decides whom to treat. Its value is V(\pi) = \mathbb E[Y(\pi(X))]. The best unrestricted policy treats exactly when \tau(x) > 0. Under the assumptions of Theorem 11.2, the advantage of \pi over treating nobody is identified and estimated with the AIPW scores: V(\pi) - V(0) = \mathbb E\big[\pi(X)\, \tau(X)\big], \qquad \hat A(\pi) = \frac1n \sum_{i=1}^n \pi(X_i)\, \Gamma_i .

Policy learning maximises \hat A over an interpretable class \Pi, such as shallow decision trees. Athey and Wager (2021) show that, with cross-fitted doubly robust scores, the regret of the learned policy against the best policy in \Pi is of order \sqrt{\mathrm{VC}(\Pi)/n}. In the simulation above, \tau(x) = 1 + x_2, so the best policy treats when x_2 > -1 (Exercise 15).

12. Research directions and open problems

This self-study chapter follows five threads of the course into current research. Each paragraph names the question, a few papers that shaped it, and what remains open. It assumes ch. 0–11.

12.1 Tighter bounds and better amortisation

The importance-weighted bound of Theorem 8.3 raised a question about training. Rainforth and coauthors (2018) show that, as the number of samples K grows, the signal-to-noise ratio of the encoder’s gradient estimate decreases. A tighter bound for the decoder can thus mean a worse-trained encoder.

Cremer, Li and Duvenaud (2018) measure the amortisation gap separately from the approximation gap. On their benchmarks, the amortisation gap is often the larger of the two.

• Open: how to train inference networks whose gap is controlled, and how much the encoder family matters once the decoder is expressive.

12.2 Deep latent variable models with missing data

Work by Mattei and Frellsen (2018) shows that the likelihood of a DLVM with Gaussian outputs and learned variances can be unbounded, as for Gaussian mixtures (ch. 2.1). They also give conditions under which the maximum-likelihood problem is well posed.

The same authors’ MIWAE (Mattei and Frellsen, 2019) trains a DLVM on data with entries missing at random, by an importance-weighted bound over the observed coordinates only. It then imputes the missing entries by self-normalised importance sampling (Proposition 6.1).

When missingness depends on the missing values themselves (missing not at random), not-MIWAE (Ipsen, Mattei and Frellsen, 2021) adds a model of the missingness mechanism.

• Open: identifiability of the data model when the missingness mechanism is also learned, and guarantees for the imputations.

12.3 Network models

For stochastic block models, Celisse, Daudin and Pierre (2012) prove that the maximum-likelihood and variational estimators of the parameters are consistent in dense graphs. Latouche, Birmelé and Ambroise (2012) develop variational Bayes inference for the SBM, with a model selection criterion derived from the variational bound.

In sparse graphs, recovery has sharp thresholds. Decelle, Krzakala, Moore and Zdeborová (2011) predicted, by statistical-physics arguments, the threshold for two equal communities with edge probabilities a/n inside and b/n across, so that the average degree is (a + b)/2. Detection is possible only when (a - b)^2 > 2(a + b). Abbe (2018) surveys the proofs that followed.

• Open: guarantees for selecting the number of classes, and models for networks with node attributes, weights or time stamps.

12.4 Theory of diffusion models

Convergence guarantees for diffusion models now hold under weak assumptions. Chen, Chewi, Li, Li, Salim and Zhang (2023) show that, if the learned score is accurate in L^2, the sampled distribution is close in total variation to the data distribution. Their bound grows polynomially in the dimension, without log-concavity. De Bortoli (2022) gives bounds when the data lie on a low-dimensional manifold.

A second question is why trained models generalise instead of reproducing their training set, since the exact score of the empirical distribution would only return training points. Kadkhodaie, Guth, Simoncelli and Mallat (2024) observe that two networks trained on disjoint halves of a data set generate nearly identical images, and relate this to adaptive harmonic bases.

Training without simulating SDEs is a third thread. Flow matching (Lipman et al., 2023), rectified flows (Liu, Gong and Liu, 2023) and stochastic interpolants (Albergo and Vanden-Eijnden, 2023) regress a velocity field directly. Consistency models (Song, Dhariwal, Chen and Sutskever, 2023) aim at one-step sampling.

• Open: a theory of generalisation for score networks, and sharp trade-offs between the number of sampling steps and accuracy.

12.5 Causal inference with machine learning

Double/debiased machine learning (Chernozhukov et al., 2018) made the AIPW construction of ch. 11.3 the default for using flexible regressions as nuisances. For heterogeneous effects, causal forests (Wager and Athey, 2018), the R-learner (Nie and Wager, 2021) and the DR-learner (Kennedy, 2023) estimate \tau(x) with guarantees. Policy learning with doubly robust scores is analysed by Athey and Wager (2021).

All of these rely on ignorability, which data cannot test. Sensitivity analysis quantifies how much hidden confounding would overturn a conclusion. Causal representation learning (Schölkopf et al., 2021) asks how to find causal variables when only high-dimensional observations, such as images, are available.

• Open: inference with weak overlap, combining randomised and observational data, and identifiability results for learned causal representations.

Exercises

Levels: ★ apply the course, ★★ combine or prove a short result, ★★★ open-ended. Exercises 10, 14, 16 and 17 are computational; the code given for 10 and 14 runs in a few seconds.

1. ★ Conditioning a bivariate Gaussian. Let (x_a, x_b) be Gaussian with mean (1, 2), variances 4 and 3, and covariance 2. Find the law of x_a given x_b = 5, and compare its variance with that of x_a.

   Solution: by Proposition 0.3, the conditional mean is 1 + \frac23 (5 - 2) = 3, and the conditional variance is 4 - \frac{2^2}{3} = \frac83.

   x_a \mid x_b = 5 \;\sim\; \mathcal N(3,\; 8/3 \approx 2.67).

   Observing x_b removes a third of the variance of x_a, whatever the observed value.

2. ★ Heteroscedastic regression. Observe y_i = w x_i + \varepsilon_i with independent \varepsilon_i \sim \mathcal N(0, \sigma_i^2), the variances known. Derive the MLE of w and its variance.

   Evaluate for inputs (1, 2, 3), outputs (2, 3, 7) and noise variances (1, 1, 4).

   Solution: the log-likelihood is -\frac12 \sum_i (y_i - w x_i)^2 / \sigma_i^2 plus a constant. It is maximal at the weighted least-squares estimate, with Fisher information \sum_i x_i^2/\sigma_i^2: \hat w = \frac{\sum_i x_i y_i / \sigma_i^2}{\sum_i x_i^2 / \sigma_i^2} = \frac{13.25}{7.25} \approx 1.828, \qquad \mathrm{Var}(\hat w) = \frac{1}{7.25} \approx 0.138 .

   Ordinary least squares would give 29/14 \approx 2.071, because it trusts the noisy third point as much as the others.

3. ★ One EM step by hand. A one-dimensional mixture of two Gaussians with unit variances, weights (\frac12, \frac12) and means (0, 2). The data are x = (0, 1, 3). Compute the responsibilities, then the updated weights and means (variances kept at 1).

   Solution: the odds of component 2 for a point x are \exp(2x - 2). So the responsibilities of component 1 are 0.881, 0.5 and 0.018, and N_1 = 1.399, N_2 = 1.601. The M-step of Proposition 2.4 gives \pi = (0.466,\ 0.534), \qquad \mu_1 = \frac{0.5 + 3 \times 0.018}{1.399} \approx 0.396, \qquad \mu_2 \approx 2.152 .

4. ★ d-separation. Take the DAG with edges a \to c, b \to c, c \to d and b \to e. Decide whether each statement below holds for every distribution that factorises according to it: \text{(i) } a \perp b; \quad \text{(ii) } a \perp b \mid d; \quad \text{(iii) } a \perp e; \quad \text{(iv) } a \perp e \mid c; \quad \text{(v) } a \perp e \mid \{b, c\}.

   Solution: (i) yes, the collider c is unobserved and so is its descendant. (ii) no, d is a descendant of the collider, which opens the path.

   (iii) yes, the only path a \to c \leftarrow b \to e is blocked at c. (iv) no, observing c opens the collider, and b is unobserved. (v) yes, b is observed and tail-to-tail on that path.

5. ★ Forward–backward by hand. Use the HMM of the code in ch. 5.3 (same initial law, transitions and emission table) with the observations x = (0, 2). Compute the forward variables, p(x), the backward variable \beta_1 and the posterior of z_1.

   Solution: \alpha_1 = (0.35, 0.05). Then \alpha_2(0) = 0.1 \times (0.315 + 0.01) = 0.0325 and \alpha_2(1) = 0.6 \times (0.035 + 0.04) = 0.045, so p(x) = 0.0775.

   Backward: \beta_1 = (0.9 \times 0.1 + 0.1 \times 0.6,\; 0.2 \times 0.1 + 0.8 \times 0.6) = (0.15,\ 0.50).

   Hence, consistently with p(x) = 0.0525 + 0.025: p(z_1 \mid x) = \frac{(0.35 \times 0.15,\; 0.05 \times 0.50)}{0.0775} \approx (0.677,\ 0.323).

6. ★ Gaussian KL and the reparameterisation gradient. (a) Compute the KL divergence from \mathcal N(1, 0.5^2) to \mathcal N(0, 1). (b) For z = \mu + \sigma\varepsilon, give the one-sample reparameterised gradient of \mathbb E[z^2] in (\mu, \sigma) and check that it is unbiased.

   Solution: (a) \frac12 (1 + 0.25 - 1 - \log 0.25) \approx 0.818.

   (b) The estimator is 2(\mu + \sigma\varepsilon)\,(1, \varepsilon). Its expectation is (2\mu, 2\sigma), which is the gradient of \mathbb E[z^2] = \mu^2 + \sigma^2.

7. ★★ No logistic MLE on separable data. Suppose some w_0 satisfies w_0^\top x_i > 0 when y_i = 1 and w_0^\top x_i < 0 when y_i = 0.

   Show that the logistic log-likelihood has no maximiser, and that adding -\frac\lambda2 \|w\|^2 with \lambda > 0 restores a unique one.

   Solution: every term of \ell_n is negative, so \ell_n < 0. Along t w_0, each term tends to 0 as t \to \infty, so \sup \ell_n = 0 is not attained.

   With the penalty, the objective is \lambda-strongly concave (Proposition 1.5) and tends to -\infty at infinity. It has a maximiser, unique by strict concavity.

8. ★★ Ridge, MAP and shrinking uncertainty. In ch. 3.1, show that m_N is the ridge estimator with \lambda = \alpha/\beta, and that adding an observation (x, y) never increases the predictive variance at any x_\ast. Quantify the decrease.

   Solution: m_N = (X^\top X + \frac\alpha\beta I)^{-1} X^\top y, the ridge solution. Adding (x, y) replaces S_N^{-1} by S_N^{-1} + \beta x x^\top. By the Sherman–Morrison formula, the new S is S' = S_N - \frac{\beta\, S_N x x^\top S_N}{1 + \beta\, x^\top S_N x}, so the variance term drops by \beta (x_\ast^\top S_N x)^2 / (1 + \beta\, x^\top S_N x) \ge 0. It does not depend on y.

9. ★★ Counting parameters for BIC. Count the free parameters of a K-component Gaussian mixture in dimension d with full, diagonal and spherical covariances. For which sample sizes does BIC penalise each parameter more than AIC?

   Solution: there are K - 1 weights and Kd means. Full covariances add Kd(d+1)/2, diagonal ones Kd, spherical ones K.

   For d = 2 and full covariances, this gives 6K - 1, as used in ch. 3.4.

   BIC charges \log n per parameter and AIC charges 2, so BIC is stricter as soon as n > e^2 \approx 7.4, that is for n \ge 8.

10. ★★ Gibbs sampling for the Ising model. On a 16 \times 16 periodic lattice with h = 0 and couplings J_{ij} = \beta, run a Gibbs sampler and estimate the mean absolute magnetisation for \beta \in \{0.2, 0.44, 0.6\}. Explain why all sites of one colour of a checkerboard can be updated at once.

    Solution: each site’s Markov blanket is its four neighbours, all of the other colour. Given the other colour, the sites of one colour are therefore independent (Proposition 4.7), and updating them together is an exact block Gibbs step.

    import numpy as np
    
    rng = np.random.default_rng(10)
    L = 16
    ii, jj = np.indices((L, L))
    colour = (ii + jj) % 2                                  # checkerboard: two conditionally independent halves
    
    def sweep(x, beta):
        for c in (0, 1):
            nb = np.roll(x, 1, 0) + np.roll(x, -1, 0) + np.roll(x, 1, 1) + np.roll(x, -1, 1)
            p_up = 1 / (1 + np.exp(-2 * beta * nb))         # p(x_i = +1 | neighbours), h = 0
            new = np.where(rng.random((L, L)) < p_up, 1, -1)
            x = np.where(colour == c, new, x)
        return x
    
    for beta in [0.2, 0.44, 0.6]:
        x = rng.choice([-1, 1], (L, L))
        mags = []
        for t in range(2000):
            x = sweep(x, beta)
            if t >= 500:
                mags.append(abs(x.mean()))
        print(f"beta = {beta:4.2f}: mean |magnetisation| = {np.mean(mags):.3f}")

    The output is about 0.085, 0.73 and 0.975: disordered at high temperature, ordered at low temperature. The change happens near Onsager’s critical value \beta_c = \frac12 \log(1 + \sqrt 2) \approx 0.441.

11. ★★ Mean field underestimates variances. For a Gaussian target \mathcal N(\mu, \Lambda^{-1}) in dimension 2, show that the CAVI factors are q_j = \mathcal N(m_j, 1/\Lambda_{jj}), that the fixed point has m = \mu, and that 1/\Lambda_{11} \le \Sigma_{11}. Evaluate for unit variances and correlation 0.9.

    Solution: \log p is quadratic. Taking its expectation in z_2 (Proposition 7.1) leaves a quadratic in z_1 with precision \Lambda_{11}, and mean \mu_1 - \Lambda_{12}(m_2 - \mu_2)/\Lambda_{11}. The fixed point solves a linear system whose solution is m = \mu.

    By the block-inverse formula, \Sigma_{11} = \frac{1}{\Lambda_{11} - \Lambda_{12}^2/\Lambda_{22}} \;\ge\; \frac{1}{\Lambda_{11}} .

    With correlation 0.9, \Lambda_{11} = 1/0.19 \approx 5.26: \mathrm{Var}_{q}(z_1) = 0.19 \quad \text{versus} \quad \mathrm{Var}_p(z_1) = 1 .

12. ★★ The two KL projections of a bimodal law. Let p = \frac12 \mathcal N(-3, 1) + \frac12 \mathcal N(3, 1).

    (a) Find the Gaussian q minimising \mathrm{KL}(p \| q). (b) Show that q = \mathcal N(3, 1) satisfies \mathrm{KL}(q \| p) \le \log 2. (c) Compare the two candidates numerically, in both directions.

    Solution: (a) by Proposition 7.4, the projection matches the mean 0 and the variance 1 + 9 = 10. So q = \mathcal N(0, 10): it covers both modes and puts its peak where p is small.

    (b) Pointwise, p \ge \frac12 \mathcal N(3, 1) = \frac12 q. So \log(q/p) \le \log 2 everywhere; integrate against q.

    (c) Numerical integration gives:

    • q = \mathcal N(0, 10): \mathrm{KL}(p \| q) \approx 0.462 and \mathrm{KL}(q \| p) \approx 0.938;
    • q = \mathcal N(3, 1): \mathrm{KL}(p \| q) \approx 8.31 and \mathrm{KL}(q \| p) \approx 0.689.

    Each direction prefers a different candidate. Over all Gaussians, the reverse KL is smallest near \mathcal N(\pm 2.98, 1.05), one mode each: it is mode-seeking. A wider symmetric local minimum, \mathcal N(0, 7.5) with value 0.84, also exists.

    The modes must be well separated for this. With \pm 2 instead of \pm 3, the reverse KL has a single minimiser, \mathcal N(0, 4.19). This is the contrast between expectation propagation and variational inference.

13. ★★ An optimal discriminator in closed form. Take p_{\text{data}} = \mathcal N(0, 1) and p_G = \mathcal N(1, 1). Compute D^\star, and the value of the GAN game at D^\star, given that \mathrm{JS}(p_{\text{data}} \| p_G) \approx 0.111 nats.

    Solution: the density ratio is p_G/p_{\text{data}} = \exp(x - \frac12), so by Theorem 9.3 D^\star(x) = \frac{1}{1 + e^{x - 1/2}} = \sigma\big(\tfrac12 - x\big).

    The value is -\log 4 + 2 \times 0.111 \approx -1.163, above the global minimum -\log 4 \approx -1.386 reached when p_G = p_{\text{data}}.

14. ★★ Tweedie’s formula. Let y = x + \sigma\varepsilon with \varepsilon \sim \mathcal N(0, 1) independent of x, and let p_\sigma be the density of y. Prove the formula below, then check it numerically for a two-component Gaussian mixture. \mathbb E[x \mid y] = y + \sigma^2 \frac{d}{dy}\log p_\sigma(y).

    Solution: p_\sigma(y) = \int \varphi_\sigma(y - x)\, p(x)\, dx, and \varphi_\sigma'(u) = -u\, \varphi_\sigma(u)/\sigma^2. Differentiating under the integral: p_\sigma'(y) = \int \frac{x - y}{\sigma^2}\, \varphi_\sigma(y - x)\, p(x)\, dx = \frac{\mathbb E[x \mid y] - y}{\sigma^2}\; p_\sigma(y).

    import numpy as np
    from scipy.stats import norm
    
    rng = np.random.default_rng(14)
    w, m, s, sigma = np.array([0.3, 0.7]), np.array([-2.0, 1.0]), np.array([0.5, 0.8]), 0.7
    n = 400_000
    k = rng.choice(2, n, p=w)
    x = rng.normal(m[k], s[k])                       # clean data: two-component mixture
    y = x + sigma * rng.normal(size=n)               # noisy observation
    
    def score_noisy(y):                              # exact score of p_sigma (mixture of N(m, s^2 + sigma^2))
        v = s ** 2 + sigma ** 2
        dens = w * norm.pdf(y[:, None], m, np.sqrt(v))
        return (dens * (m - y[:, None]) / v).sum(1) / dens.sum(1)
    
    for y0 in [-2.0, 0.0, 1.5]:
        near = np.abs(y - y0) < 0.02                 # Monte Carlo estimate of E[x | y = y0]
        tweedie = y0 + sigma ** 2 * score_noisy(np.array([y0]))[0]
        print(f"y = {y0:4.1f}: E[x|y] by Monte Carlo {x[near].mean():6.3f}, by Tweedie {tweedie:6.3f}")

    The two estimates agree to within 0.012 at the three points, for example 1.279 and 1.283 at y = 1.5.

15. ★★ Double robustness and the value of a policy. (a) Show that when the outcome models are set to zero and \hat e = e, the AIPW score reduces to the IPW score. (b) In the simulation of ch. 11.3, where \tau(x) = 1 + x_2 and x_2 \sim \mathcal N(0, 1), compute the gain of the policy “treat when x_2 > -1” over “treat everyone”.

    Solution: (a) with \hat\mu_w = 0, the score is WY/e(X) - (1-W)Y/(1 - e(X)). It is unbiased by Theorem 11.2, consistent with Theorem 11.3.

    (b) The gain is \mathbb E[(-\tau(X))_+] = \mathbb E[(Z - 1)_+] with Z standard normal, by symmetry. With \varphi and \Phi the normal density and distribution function: \mathbb E[(Z - 1)_+] = \varphi(1) - \big(1 - \Phi(1)\big) \approx 0.242 - 0.159 = 0.083 .

16. ★★★ Block models near the detection threshold. Simulate n = 2000 nodes in two equal communities, with edge probabilities a/n inside and b/n across, and average degree (a + b)/2 = 5. Vary a - b across the threshold (a - b)^2 = 2(a + b) of ch. 12.3, and compare spectral clustering of the adjacency matrix with variational EM started from it.

    Hints and expected result: reuse the code of ch. 7.3 with sparse matrices (scipy.sparse). Below the threshold, no method should beat chance; report the adjusted Rand index with error bars over 10 graphs. Above it, look at how the eigenvectors of the adjacency matrix localise on high-degree nodes in sparse graphs, and whether VEM corrects this. Conclude on which of the two is limited by the information in the graph and which by the algorithm.

17. ★★★ Reverse SDE or probability-flow ODE? Modify the sampler of ch. 10.3 to integrate the probability-flow ODE of Theorem 10.5 (replace the update by x = x + 0.5 * g2 * score(x, s)). Compare both samplers for 10, 50 and 500 steps, using the radius statistics and the mode counts.

    Hints and expected result: in one run with the same trained network, the standard deviation of the sample radius (data: about 0.1) was 0.185 for the SDE and 0.257 for the ODE with 10 steps, 0.124 and 0.141 with 50, and 0.115 and 0.122 with 500. Plot the samples; explain what the stochastic term of the SDE corrects and what it costs. Relate the ODE to a continuous normalizing flow (ch. 9.2), and say how you would compute its likelihood.

Projects and assessment

The master-mva.com page states that “the students will be evaluated by a written exam”, and its organisation section mentions in-person sessions to answer questions on the lectures and “the project”. The Université Paris-Saclay course detail lists “Evaluation Continue non Intégrale”. No weights are displayed. The projects below train the skills of the project and of the exam; the exam proposal shows the expected level of a written test on this page’s content.

Solo projects

1. [Solo · ≈ 20 h] Paper study: importance-weighted auto-encoders. Reproduce the central comparison of Burda, Grosse and Salakhutdinov (2016) at small scale: does training on \mathcal L_K improve the test log-likelihood over a VAE?
   • Work: extend the code of ch. 8.3; train with K \in \{1, 5, 50\} on the binarised scikit-learn digits; evaluate every model with \mathcal L_{5000} on held-out images; count active latent units.
   • Extension: measure the signal-to-noise ratio of the encoder gradient as K grows (Rainforth et al., 2018).
   • Deliverables: code, 4–6 page report, 10-min talk. Chapters: 2, 6.1, 8.
   • Success: a table of test bounds with standard errors over 3 seeds, and a clear verdict.
2. [Solo · ≈ 20 h] Stochastic block models: variational EM, spectral clustering and ICL. When does variational EM improve on spectral clustering, and does ICL find the number of classes?
   • Work: implement Proposition 7.3 and the ICL criterion; test on synthetic assortative, disassortative and core–periphery SBMs; apply to Zachary’s karate club (bundled with NetworkX).
   • Deliverables: code, 4–6 page report, 10-min talk. Chapters: 4, 7 (and 3.4 for model selection).
   • Success: adjusted Rand index curves against the signal strength for both methods; ICL selects the true number of classes on at least 8 of 10 graphs per family at n = 300, with the block probabilities 0.10 apart as in ch. 7.3.
3. [Solo · ≈ 15 h] Hidden Markov models for regime switching. How much does modelling time help clustering?
   • Work: implement Baum–Welch with Gaussian emissions and scaled forward–backward (ch. 5.3–5.4); simulate a 3-regime series; compare Viterbi decoding with a Gaussian mixture that ignores time (ch. 2.3); choose the number of states by BIC.
   • Deliverables: code, 4-page report, 10-min talk. Chapters: 2, 3.4, 5.
   • Success: decoding accuracy above the mixture’s, and BIC selecting 3 states on most of 20 simulated series.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Paper study: score-based generative modelling at small scale. Reproduce the main mechanisms of Song et al. (2021), with Ho, Jain and Abbeel (2020) as the discrete-time reference, on 2D data and the 8×8 digits of scikit-learn.
   • Split: member A trains variance-exploding and variance-preserving models by denoising score matching; member B implements the samplers (reverse SDE, probability-flow ODE, Langevin corrector from ch. 6.4); member C evaluates (likelihood through the ODE with a trace estimator, a digit classifier for sample quality).
   • Deliverables: shared code, 8–10 page report, 20-min talk. Chapters: 6.4, 9.2, 10.
   • Success: on the eight-Gaussian data of ch. 10.3, 500 reverse-SDE steps give all eight modes and a mean sample radius within 0.02 of 2. Curves of the radius statistics and of the classifier’s class balance against the number of steps (10 to 1000), for each sampler. The ODE log-likelihood of held-out 2D points, compared with a coupling flow (ch. 9.2) trained on the same data.
2. [Group of 2–3 · ≈ 30 h each] Treatment effects and policies from observational data. How reliable are IPW, AIPW and policy learning when overlap weakens or models are misspecified?
   • Split: member A builds a semi-synthetic generator with known effects and controllable overlap; member B implements regression adjustment, IPW and cross-fitted AIPW with scikit-learn nuisances and confidence intervals; member C learns depth-2 tree policies from the AIPW scores (Athey and Wager, 2021).
   • Deliverables: shared code, 8–10 page report, 20-min talk. Chapters: 4.2, 11.
   • Success: empirical coverage of the 95 % intervals over 200 replications, and policy regret against the known optimum.

Written exam (proposal)

Format (proposal, not an official rule): 3 hours, one handwritten A4 sheet allowed, no electronic device, four exercises.

1. Question (5 points). State the ELBO decomposition of \log p(x; \theta) for a latent-variable model, describe the two steps of EM, and prove that EM never decreases the likelihood. When is the bound tight?

   Expected: Proposition 2.2 with its two-line proof; the E-step sets q to the posterior and the M-step maximises the expected complete log-likelihood; the chain of inequalities of Theorem 2.3; tightness exactly when q is the posterior.

2. Question (6 points). For the HMM of ch. 5.3 and the observations x = (0, 2, 2): compute the forward variables and p(x); compute p(z_1 \mid x); find the Viterbi path. Comment on the two answers.

   Expected: the forward variables are (0.35, 0.05), (0.0325, 0.045) and (0.003825, 0.02355), so p(x) = 0.027375. Then \beta_1 = (0.0435, 0.243) and p(z_1 \mid x) \approx (0.556, 0.444).

   The Viterbi path is (1, 1, 1), with joint probability 0.01152. The most probable path does not start with the most probable first state: marginal and joint maximisation differ.

3. Question (4 points). Prove that the AIPW estimator of the average treatment effect is unbiased when the propensity score is correct, whatever the outcome models.

   Expected: condition on X; use ignorability to write the expectation of the treated correction as e/\hat e times \mu_1 - \hat\mu_1; with \hat e = e the outcome model cancels; the control part is symmetric; conclude with Theorem 11.2.

4. Question (5 points). Draw the stochastic block model as a directed graphical model, moralise it, and explain why the E-step of EM is intractable. Derive the mean-field update of \tau_i.

   Expected: nodes z_i with arrows into each X_{ij}; moralisation links every pair (z_i, z_j), so the posterior couples all labels and its normaliser sums over K^n labellings. The update follows from Proposition 7.1, giving the fixed point of Proposition 7.3 with a short justification.

Grading rubric (proposal)

References

Official references (master-mva.com):

• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — free PDF.
• C. M. Bishop and H. Bishop, Deep Learning: Foundations and Concepts, Springer Nature, 2023 — book site.
• K. P. Murphy, Probabilistic Machine Learning: An Introduction, MIT Press, 2022 — book site.
• K. P. Murphy, Probabilistic Machine Learning: Advanced Topics, MIT Press, 2023 — book site.
• S. Wager, Causal Inference: A Statistical Learning Approach, 2024 — PDF.

Textbooks:

• D. Koller and N. Friedman, Probabilistic Graphical Models: Principles and Techniques, MIT Press, 2009 — ch. 4–5.
• S. L. Lauritzen, Graphical Models, Oxford University Press, 1996 — ch. 4.
• M. J. Wainwright and M. I. Jordan, “Graphical models, exponential families, and variational inference”, Foundations and Trends in Machine Learning 1(1–2), 2008 — ch. 0.4, 5, 7.
• C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, MIT Press, 2006 — ch. 3, 7.4.
• C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004 — ch. 6.
• A. W. van der Vaart, Asymptotic Statistics, Cambridge University Press, 1998 — ch. 1.
• G. W. Imbens and D. B. Rubin, Causal Inference for Statistics, Social, and Biomedical Sciences, Cambridge University Press, 2015 — ch. 11.
• J. Pearl, Causality, 2nd ed., Cambridge University Press, 2009 — ch. 4, 11.

Papers, by chapter:

• H. Akaike, “A new look at the statistical model identification”, IEEE Transactions on Automatic Control 19(6), 1974; G. Schwarz, “Estimating the dimension of a model”, Annals of Statistics 6(2), 1978 — ch. 3.4.
• D. Arthur and S. Vassilvitskii, “k-means++: the advantages of careful seeding”, SODA, 2007 — ch. 1.4.
• A. P. Dempster, N. M. Laird and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm”, Journal of the Royal Statistical Society B 39(1), 1977; C. F. J. Wu, “On the convergence properties of the EM algorithm”, Annals of Statistics 11(1), 1983 — ch. 2.
• R. M. Neal and G. E. Hinton, “A view of the EM algorithm that justifies incremental, sparse, and other variants”, in M. I. Jordan (ed.), Learning in Graphical Models, Kluwer, 1998 — ch. 3.3.
• M. E. Tipping and C. M. Bishop, “Probabilistic principal component analysis”, Journal of the Royal Statistical Society B 61(3), 1999 — ch. 2.4.
• J. Besag, “Spatial interaction and the statistical analysis of lattice systems”, Journal of the Royal Statistical Society B 36(2), 1974; D. Geiger, T. Verma and J. Pearl, “Identifying independence in Bayesian networks”, Networks 20(5), 1990 — ch. 4.
• J. Friedman, T. Hastie and R. Tibshirani, “Sparse inverse covariance estimation with the graphical lasso”, Biostatistics 9(3), 2008 — ch. 4.4.
• G. F. Cooper, “The computational complexity of probabilistic inference using Bayesian belief networks”, Artificial Intelligence 42(2–3), 1990; F. R. Kschischang, B. J. Frey and H.-A. Loeliger, “Factor graphs and the sum-product algorithm”, IEEE Transactions on Information Theory 47(2), 2001; K. P. Murphy, Y. Weiss and M. I. Jordan, “Loopy belief propagation for approximate inference: an empirical study”, UAI, 1999 — ch. 5.
• L. R. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition”, Proceedings of the IEEE 77(2), 1989 — ch. 5.3–5.4.
• G. O. Roberts, A. Gelman and W. R. Gilks, “Weak convergence and optimal scaling of random walk Metropolis algorithms”, Annals of Applied Probability 7(1), 1997; G. O. Roberts and R. L. Tweedie, “Exponential convergence of Langevin distributions and their discrete approximations”, Bernoulli 2(4), 1996; A. Durmus and É. Moulines, “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm”, Annals of Applied Probability 27(3), 2017 — ch. 6.
• D. M. Blei, A. Kucukelbir and J. D. McAuliffe, “Variational inference: a review for statisticians”, Journal of the American Statistical Association 112(518), 2017 — ch. 7.1.
• P. W. Holland, K. B. Laskey and S. Leinhardt, “Stochastic blockmodels: first steps”, Social Networks 5(2), 1983; K. Nowicki and T. A. B. Snijders, “Estimation and prediction for stochastic blockstructures”, Journal of the American Statistical Association 96(455), 2001; J.-J. Daudin, F. Picard and S. Robin, “A mixture model for random graphs”, Statistics and Computing 18(2), 2008 — ch. 7.3.
• T. P. Minka, “Expectation propagation for approximate Bayesian inference”, UAI, 2001 — ch. 7.4.
• D. P. Kingma and M. Welling, “Auto-encoding variational Bayes”, ICLR, 2014; D. J. Rezende, S. Mohamed and D. Wierstra, “Stochastic backpropagation and approximate inference in deep generative models”, ICML, 2014 — ch. 8.2.
• Y. Burda, R. Grosse and R. Salakhutdinov, “Importance weighted autoencoders”, ICLR, 2016; C. Cremer, X. Li and D. Duvenaud, “Inference suboptimality in variational autoencoders”, ICML, 2018; T. Rainforth, A. R. Kosiorek, T. A. Le, C. J. Maddison, M. Igl, F. Wood and Y. W. Teh, “Tighter variational bounds are not necessarily better”, ICML, 2018 — ch. 8, 12.1.
• M. Germain, K. Gregor, I. Murray and H. Larochelle, “MADE: masked autoencoder for distribution estimation”, ICML, 2015; A. van den Oord, N. Kalchbrenner and K. Kavukcuoglu, “Pixel recurrent neural networks”, ICML, 2016 — ch. 9.1.
• L. Dinh, J. Sohl-Dickstein and S. Bengio, “Density estimation using Real NVP”, ICLR, 2017; G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed and B. Lakshminarayanan, “Normalizing flows for probabilistic modeling and inference”, Journal of Machine Learning Research 22(57), 2021 — ch. 9.2.
• I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville and Y. Bengio, “Generative adversarial nets”, NeurIPS, 2014; M. Arjovsky, S. Chintala and L. Bottou, “Wasserstein generative adversarial networks”, ICML, 2017 — ch. 9.3.
• G. E. Hinton, “Training products of experts by minimizing contrastive divergence”, Neural Computation 14(8), 2002; A. Hyvärinen, “Estimation of non-normalized statistical models by score matching”, Journal of Machine Learning Research 6, 2005; P. Vincent, “A connection between score matching and denoising autoencoders”, Neural Computation 23(7), 2011 — ch. 10.1–10.2.
• B. D. O. Anderson, “Reverse-time diffusion equation models”, Stochastic Processes and their Applications 12(3), 1982; J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan and S. Ganguli, “Deep unsupervised learning using nonequilibrium thermodynamics”, ICML, 2015; J. Ho, A. Jain and P. Abbeel, “Denoising diffusion probabilistic models”, NeurIPS, 2020; Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon and B. Poole, “Score-based generative modeling through stochastic differential equations”, ICLR, 2021 — ch. 10.3.
• P. R. Rosenbaum and D. B. Rubin, “The central role of the propensity score in observational studies for causal effects”, Biometrika 70(1), 1983; J. M. Robins, A. Rotnitzky and L. P. Zhao, “Estimation of regression coefficients when some regressors are not always observed”, Journal of the American Statistical Association 89(427), 1994 — ch. 11.2–11.3.
• V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey and J. Robins, “Double/debiased machine learning for treatment and structural parameters”, The Econometrics Journal 21(1), 2018; S. Athey and S. Wager, “Policy learning with observational data”, Econometrica 89(1), 2021 — ch. 11.3–11.4.

Research directions (ch. 12):

• P.-A. Mattei and J. Frellsen, “Leveraging the exact likelihood of deep latent variable models”, NeurIPS, 2018; P.-A. Mattei and J. Frellsen, “MIWAE: deep generative modelling and imputation of incomplete data sets”, ICML, 2019; N. B. Ipsen, P.-A. Mattei and J. Frellsen, “not-MIWAE: deep generative modelling with missing not at random data”, ICLR, 2021.
• A. Celisse, J.-J. Daudin and L. Pierre, “Consistency of maximum-likelihood and variational estimators in the stochastic block model”, Electronic Journal of Statistics 6, 2012; P. Latouche, E. Birmelé and C. Ambroise, “Variational Bayesian inference and complexity control for stochastic block models”, Statistical Modelling 12(1), 2012.
• A. Decelle, F. Krzakala, C. Moore and L. Zdeborová, “Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications”, Physical Review E 84, 066106, 2011; E. Abbe, “Community detection and stochastic block models: recent developments”, Journal of Machine Learning Research 18(177), 2018.
• S. Chen, S. Chewi, J. Li, Y. Li, A. Salim and A. R. Zhang, “Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions”, ICLR, 2023; V. De Bortoli, “Convergence of denoising diffusion models under the manifold hypothesis”, Transactions on Machine Learning Research, 2022; Z. Kadkhodaie, F. Guth, E. P. Simoncelli and S. Mallat, “Generalization in diffusion models arises from geometry-adaptive harmonic representations”, ICLR, 2024.
• Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel and M. Le, “Flow matching for generative modeling”, ICLR, 2023; X. Liu, C. Gong and Q. Liu, “Flow straight and fast: learning to generate and transfer data with rectified flow”, ICLR, 2023; M. S. Albergo and E. Vanden-Eijnden, “Building normalizing flows with stochastic interpolants”, ICLR, 2023; Y. Song, P. Dhariwal, M. Chen and I. Sutskever, “Consistency models”, ICML, 2023.
• S. Wager and S. Athey, “Estimation and inference of heterogeneous treatment effects using random forests”, Journal of the American Statistical Association 113(523), 2018; X. Nie and S. Wager, “Quasi-oracle estimation of heterogeneous treatment effects”, Biometrika 108(2), 2021; E. H. Kennedy, “Towards optimal doubly robust estimation of heterogeneous causal effects”, Electronic Journal of Statistics 17(2), 2023; B. Schölkopf, F. Locatello, S. Bauer, N. R. Ke, N. Kalchbrenner, A. Goyal and Y. Bengio, “Toward causal representation learning”, Proceedings of the IEEE 109(5), 2021.