Probability assigns events a number in \([0,1]\):
Example: dice roll, \(P(\{4\})=1/6\).
For continuous variables, we require an integral for it to make sense:
\[ P(3 \le x \le 10) = \int_3^{10}p(x)\,dx =0.95. \]
Probability = area under the probability density
\[ (\Omega,\mathcal F,P) \]
\(\Omega\): outcome / parameter space
\(\mathcal F\): measurable events \(A\subseteq\Omega\)
\(P\): probability measure
\[ P(A)\ge0,\qquad P(\Omega)=1. \]
For disjoint \(A_i\):
\[ P\!\left(\bigcup_i A_i\right)=\sum_iP(A_i). \]
\[ \Omega=\mathbb R^n. \]
We integrate with respect to Lebesgue measure on \(\mathbb R^n\), written \(dx\).
From now on, densities are with respect to this measure:
\[ p(x)=\frac{dP}{dx}(x), \qquad P(A)=\int_A p(x)\,dx. \]
Define a Bayesian model: prior + likelihood \(\rightarrow\) joint density.
Prior choice: \(p(x\mid M)\) \(\propto\) \(\sigma^{-1}\)
Likelihood for \(D\): \(p(D\mid x,M)\) \(=\) \((2\pi\sigma^2)^{-3/2}\prod_{i=1}^3\exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]\)
Joint density:
| \(p(D,x\mid M)\) | \(=\) | \(p(D\mid x,M)\) \(p(x\mid M)\) |
| \(p(D,x\mid M)\) | \(\propto\) | \(\sigma^{-3}\prod_{i=1}^3\exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]\) \(\sigma^{-1}\) |
Parameter space:
\[ x=(\mu,\sigma)\in\mathbb R\times\mathbb R_+. \]
Factorization:
\[ p(D,x\mid M) = p(D\mid x,M)p(x\mid M). \]
Likelihood contribution:
\[ p(D\mid x,M)\propto_x \sigma^{-3} \prod_{i=1}^3 \exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]. \]
Prior contribution:
\[ p(x\mid M)\propto\sigma^{-1}. \]
Therefore:
\[ p(D,x\mid M) \propto_x \sigma^{-4} \prod_{i=1}^3 \exp\!\left[-\frac{(y_i-\mu)^2}{2\sigma^2}\right]. \]
Bayes, “An essay towards solving a problem in the doctrine of chances”, Philos. Trans. R. Soc. Lond., 1763.
Same joint density:
\[ p(D,x\mid M). \]
Condition on \(x\): likelihood.
Condition on \(D\): posterior.
\(p(D\mid x,M)\) \(=\) \(p(D,x\mid M)\) \(p(x\mid M)\)
But this is the term we defined before.
\(p(x\mid D,M)\) \(=\) \(p(D\mid x,M)\) \(p(x\mid M)\) \(p(D\mid M)\)
Bayes
\(p(D\mid M)\) \(=\) \(\int\) \(p(D\mid x,M)\) \(p(x\mid M)\) \(dx\)
Joint:
\[ p(D,x\mid M) = p(D\mid x,M)p(x\mid M). \]
Likelihood:
\[ p(D\mid x,M) = \frac{p(D,x\mid M)} {\int p(D,x\mid M)\,dD}. \]
Parameter conditional:
\[ p(x\mid D) = \frac{p(D,x\mid M)} {\int p(D,x\mid M)\,dx}. \]
For this model:
\[ dx=d\mu\,d\sigma. \]
Bayes, “An essay towards solving a problem in the doctrine of chances”, Philos. Trans. R. Soc. Lond., 1763.
Example: is variable \(x_3\) larger than \(0\)?
\[ P(x_3>0\mid D) = \int_{\{x:x_3>0\}} p(x\mid D)\,dx. \]
Many posterior questions become calculations of the form:
\[ \mathbb E_p[f] = \int_\Omega f(x)\,p(x\mid D)\,dx. \]
\(f(x) = x\) calculate the mean of the variable \(x\).
Variable/parameter:
\[ x\in\Omega. \]
Observable:
\[ f:\Omega\to\mathbb R. \]
Second moment:
\[ \mathbb E_p[x_j^2] = \int x_j^2\,p(x\mid D)\,dx. \]
Interval probability:
\[ f(x)=\mathbf 1_A(x), \qquad P(x\in A\mid D)=\mathbb E_p[f]. \]
Credible interval \(C_{1-\alpha}\):
\[ F_j(t) = P(x_j\le t\mid D). \]
\[ \begin{aligned} F_j(a)&=\alpha/2,\\ F_j(b)&=1-\alpha/2 \end{aligned} \quad\text{solve for }a,b. \]
\[ C_{1-\alpha}=[a,b]. \]
Easiest approach: Riemann sum.
Replace area by small rectangles:
\[
\int f(x)p(x)\,dx
\approx
\sum_i f(x_i)p(x_i)\Delta x.
\] 
More intelligent nodes and weights: Gaussian quadrature.
Excellent in low dimension when the integrand is polynomial-like.
Integrand:
\[ g(x)=f(x)p(x\mid D). \]
Riemann:
\[ \int g(x)\,dx \approx \sum_i g(x_i)\Delta x. \]
Gaussian quadrature:
\[ \int_a^b g(x)\rho(x)\,dx \approx Q_n(g) := \sum_{i=1}^n \lambda_i g(x_i). \]
Choose \(x_i\) as roots of the degree-\(n\) polynomial orthogonal w.r.t. \(\rho\).
Optimal in polynomial class: orthogonality cancels the interpolation-error moments, so \(Q_n\) matches \(2n\) moments:
\[ Q_n(g)=\int_a^b g(x)\rho(x)\,dx \quad\text{for }\operatorname{degree}(g)\le 2n-1. \]
Riemann, “Ueber die Darstellbarkeit einer Function durch eine trigonometrische Reihe”, Abh. K. Ges. Wiss. Gottingen, 1868. Gauss, “Methodus nova integralium valores per approximationem inveniendi”, Comment. Soc. Regiae Sci. Gotting. Recent., 1814.
Better deterministic idea: refine intervals/cells only when needed.
Use a local error estimate:
On a cell \(I\):
\[ Q_I(g)\approx\int_I g(x)\,dx. \]
Embedded estimate:
\[ \widehat e_I = |Q_I^{\mathrm{high}}(g)-Q_I^{\mathrm{low}}(g)|. \]
Refine if:
\[ \widehat e_I>\tau_I. \]
Examples:
Still partitions ambient space.
Kuncir, “Algorithm 103: Simpson’s rule integrator”, Commun. ACM, 1962. Gander et al., “Adaptive Quadrature - Revisited”, BIT Numer. Math., 2000.
Mesh depends on several choices/parameters.
Simpler/more elegant: random delta peaks.
\[ p(x)\,dx \approx \sum_i w_i\delta_{x_i}(dx). \]
Draw from an easier proposal:
\[ x_i\sim q. \]
Correct by density ratios:
\[ w_i = \frac{p(x_i)/q(x_i)} {\sum_j p(x_j)/q(x_j)}. \]
Failure: proposal misses support or tails.
For \(N\) proposals:
\[ x_i\sim q, \qquad r_i=\frac{p(x_i)}{q(x_i)}. \]
Normalized weights:
\[ \bar w_i = \frac{r_i}{\sum_j r_j}. \]
Useful fraction:
\[ \eta_{\mathrm{IS},N} = \frac{N_{\mathrm{eff}}}{N} \approx \frac{(\sum_i r_i)^2} {N\sum_i r_i^2}. \]
Large \(N\):
\[ \eta_{\mathrm{IS}} \to \frac{Z^2}{\mathbb E_q r^2}, \qquad Z=\mathbb E_q r. \]
If \(p\) normalized: \(\mathbb E_q r=1\).
Kahn, “Random sampling (Monte Carlo) techniques in neutron attenuation problems–I”, Nucleonics, 1950. Kish, “Survey Sampling”, Wiley, 1965.
Weighted samples: correct, but not a real sampler from \(p\).
Unweighted samples by accept/reject:
\[ x\sim q, \qquad u\sim \mathrm{Unif}(0,1), \qquad u\le \frac{p(x)}{c\,q(x)}. \]
Need valid probabilities:
\[ c\,q(x)\ge p(x) \quad\text{for all }x. \]
Retained fraction \(\approx Z/c\); cost \(\approx c/Z\) proposals per sample.
For \(N\) proposals:
\[ x_i\sim q, \qquad r_i=\frac{p(x_i)}{q(x_i)}, \qquad Z=\mathbb E_q r. \]
Envelope:
\[ r_i\le c \quad\text{for all proposals.} \]
Binary thinning:
\[ B_i\mid x_i \sim \mathrm{Bernoulli}(r_i/c), \qquad A=\sum_i B_i. \]
Useful fraction:
\[ \eta_{\mathrm{rej},N} = \frac{A}{N}, \qquad \mathbb E[\eta_{\mathrm{rej},N}\mid x_{1:N}] = \frac{1}{Nc}\sum_i r_i. \]
Large \(N\):
\[ \eta_{\mathrm{rej}} \to \frac{Z}{c}, \qquad Z=\mathbb E_q r. \]
Compare to IS:
\[ \eta_{\mathrm{IS}} = \frac{Z^2}{\mathbb E_q r^2}, \qquad \eta_{\mathrm{rej}} = \frac{Z}{c}. \]
von Neumann, “Various techniques used in connection with random digits”, NBS Appl. Math. Ser., 1951.
Function surrogate:
\[ q_\phi(x)\approx p(x). \]
Classical: Laplace/Gaussian, polynomial or spline quadrature.
Gaussian quadrature: exact for polynomial test functions.
Normalizing flow: invertible map of simple \(q_0\).
\[ z\sim q_0,\qquad x=T_\phi(z). \]
If approximate: reweight samples, \(\bar w_i\propto p(x_i)/q_\phi(x_i)\).
Problem: missed support stays invisible.
Quadrature:
\[ \int g(x)\rho(x)\,dx \approx \sum_i \lambda_i g(x_i). \]
Exact on a chosen function class.
Family:
\[ q_\phi\in\mathcal Q. \]
Flow:
\[ q_\phi=T_{\phi\#}q_0. \]
Fit:
\[ q_\phi = \arg\min_{q\in\mathcal Q} \mathrm{KL}(q\|p/Z). \]
\[ \mathrm{KL}(q\|p/Z) = \int q(x)\log\frac{q(x)}{p(x)/Z} \,dx. \]
Zero support:
\[ q_\phi(x)=0 \Rightarrow \text{no samples there}. \]
Laplace, “Memoire sur les approximations des formules qui sont fonctions de tres grands nombres”, Mem. Acad. R. Sci. Paris, 1785. Jordan et al., “An Introduction to Variational Methods for Graphical Models”, Mach. Learn., 1999. Rezende et al., “Variational Inference with Normalizing Flows”, ICML, 2015.
Global proposals might miss probability mass. Solution: move locally.
Define only a proposal kernel:
\[ x_{t+1}^\star\sim q(\cdot\mid x_t), \qquad q(x_{t+1}^\star\mid x_t). \]
Add Metropolis-Hastings (MH) correction (like rejection sampling):
\[ \alpha(x_t,x_{t+1}^\star) = \min\!\left( 1, \frac{p(x_{t+1}^\star)q(x_t\mid x_{t+1}^\star)} {p(x_t)q(x_{t+1}^\star\mid x_t)} \right). \]
\(\alpha\): acceptance probability.
Current state:
\[ x_t=x. \]
Proposal:
\[ y\sim q(\cdot\mid x). \]
Proposal ratio:
\[ r(x,y) = \frac{p(y)q(x\mid y)} {p(x)q(y\mid x)}. \]
Accept:
\[ \alpha(x,y)=\min(1,r(x,y)). \]
If not accepted, add the current value again:
\[ x_{t+1} = \begin{cases} y, & \text{with probability }\alpha(x,y),\\ x, & \text{otherwise}. \end{cases} \]
Transition:
\[ K(x,dy) = q(y\mid x)\alpha(x,y)dy +\rho(x)\delta_x(dy). \]
Sample measure:
\[ \widehat\Pi_N = \frac1N\sum_{t=1}^N\delta_{x_t}. \]
Metropolis et al., “Equation of State Calculations by Fast Computing Machines”, J. Chem. Phys., 1953. Hastings, “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”, Biometrika, 1970.
Pairwise detailed balance: no net local flow.
\[ p(x)K(x,dy) = p(y)K(y,dx). \]
Integrate all balanced pairs: target density invariant.
\[ pK=p. \]
If the chain is ergodic, it converges to \(p\).
Then one long path gives posterior averages.
Target measure:
\[ \Pi(dx)=p(x)\,dx. \]
For any event \(A\):
\[ \begin{aligned} (\Pi K)(A) &= \int_\Omega p(x)K(x,A)\,dx\\ &= \int_\Omega\int_A p(x)K(x,dy)\,dx\\ &\overset{\mathrm{balance}}{=} \int_A\int_\Omega p(y)K(y,dx)\,dy\\ &= \int_A p(y)\,dy = \Pi(A). \end{aligned} \]
For the Markov chain:
\[ \begin{gathered} X_0=x_0,\qquad X_t\sim K^t(x_0,\cdot),\\ K^t(x_0,A) = \mathbb P(X_t\in A\mid X_0=x_0),\\ \text{balance} \Rightarrow \Pi K=\Pi \quad\text{(invariant target)},\\ \text{ergodicity conditions:}\\ \Pi\text{-irreducible},\quad \text{aperiodic},\quad \text{positive Harris}\\ \Rightarrow \lim_{t\to\infty}K^t(x_0,A)=\Pi(A) \quad\text{(independent of start)},\\ \Rightarrow \lim_{N\to\infty} \frac1N\sum_{t=1}^N f(X_t) = \int f\,d\Pi . \end{gathered} \]
Metropolis et al., “Equation of State Calculations by Fast Computing Machines”, J. Chem. Phys., 1953. Birkhoff, “Proof of the Ergodic Theorem”, Proc. Natl. Acad. Sci. U.S.A., 1931.
We already know: the transition kernel preserves the target.
\[ p K=p. \]
But local random-walk moves are not efficient:
\[ Y=X+\sigma\xi, \qquad \alpha(x,y)= \min\!\left\{1,\frac{p(y)}{p(x)}\right\}. \]
→ We want the proposal distribution to follow nearby probability mass.
Balance implies the invariant target:
\[ p(x)K(x,dy) = p(y)K(y,dx) \Rightarrow pK=p. \]
\[ Y=X+\sigma_d\xi, \qquad \sigma_d=\ell/\sqrt d. \]
For product targets:
\[ p_d(x)=\prod_{i=1}^d p_1(x_i), \qquad I=\mathbb E[(\partial_x\log p_1(X))^2]. \]
Then
\[ \bar\alpha_d(\ell) \to 2\Phi\!\left(-\frac{\ell\sqrt I}{2}\right), \qquad \ell_{\mathrm{opt}}: \bar\alpha\approx0.234. \]
So useful local moves require:
\[ \text{follow mass} \quad+\quad \text{small correction}. \]
Roberts et al., “Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms”, Ann. Appl. Probab., 1997.
Solution: gradient drift.
Drift the proposal toward higher density:
\[ \begin{aligned} Y &= X +\frac{h}{2}\nabla\log p(X) +\sqrt{h}\,\xi,\\ \xi&\sim N(0,I). \end{aligned} \]
After MH correction: MALA = Metropolis-adjusted Langevin algorithm.
Forward proposal:
\[ q_h(y\mid x) = N\!\left(y;\, x+\frac h2\nabla\log p(x),\, hI\right). \]
Balance compares flows:
\[ p(x)q_h(y\mid x) \quad\text{vs.}\quad p(y)q_h(x\mid y). \]
Therefore:
\[ \alpha(x,y) = \min\left( 1, \frac{p(y)q_h(x\mid y)} {p(x)q_h(y\mid x)} \right). \]
Not symmetric: contains drift \(\nabla\log p(x)\).
\[ q_h(y\mid x)\ne q_h(x\mid y) \]
Drift and diffusion flows cancel:
\[ \begin{gathered} dX_t = \frac12\nabla\log p(X_t)\,dt+dW_t,\\ \rho_t=\text{density}(X_t), \qquad W_t=\text{Brownian motion}. \end{gathered} \]
Probability current \(J_t\): net mass flux.
\[ J_t := \frac12\rho_t\nabla\log p - \frac12\nabla\rho_t, \qquad \partial_t\rho_t=-\nabla\cdot J_t. \]
At \(\rho_t=p\): \(J_t=\frac12p\nabla\log p-\frac12\nabla p=0\).
Rossky et al., “Brownian dynamics as smart Monte Carlo simulation”, J. Chem. Phys., 1978.
Remaining limitation: how far can one gradient be extrapolated?
Useful step length depends on how that gradient changes with \(x\).
Use \(U(x)=-\log p(x)\). Around the MAP point:
\[ x_{\mathrm{MAP}} := \underset{x}{\operatorname{arg\,max}}\,p(x) = \underset{x}{\operatorname{arg\,min}}\,U(x), \qquad \nabla U(x_{\mathrm{MAP}})=0. \]
The simplest local model: gradient is linear in \(\delta x\) with slope 1 (unit local scale).
\[ \begin{gathered} \delta x=x-x_{\mathrm{MAP}},\\ \nabla U(x_{\mathrm{MAP}}+\delta x) \approx \delta x. \end{gathered} \]
This is the Taylor / Gaussian local model.
Potential:
\[ U(x)=-\log(p(x)). \]
MAP point / minimum:
\[ x_{\mathrm{MAP}} = \underset{x}{\operatorname{arg\,max}}\,p(x) = \underset{x}{\operatorname{arg\,min}}\,U(x). \]
\[ \nabla U(x_{\mathrm{MAP}})=0, \qquad \delta x=x-x_{\mathrm{MAP}}. \]
Unit local scale:
\[ \nabla U(x_{\mathrm{MAP}}+\delta x) \approx \delta x. \]
Taylor at a minimum:
\[ U(x_{\mathrm{MAP}}+\delta x) \approx U(x_{\mathrm{MAP}}) + \frac12\|\delta x\|^2. \]
Local Gaussian:
\[ p_{\mathrm{loc}}(\delta x) \propto \exp\!\left(-\frac12\|\delta x\|^2\right). \]
If \(\nabla U\) changes like \(\delta x\), the measured gradient stays predictive over a local step.
Easy case for:
If the local model is anisotropic, “slopes” \(H\) differ in different directions:
\[ \nabla U(x_{\mathrm{MAP}}+\delta x)\approx H\delta x, \qquad H:=\nabla^2 U(x_{\mathrm{MAP}}). \]
with \(H\) the Hessian of \(U\), only diagonal if slopes are aligned with the axes.
Preconditioning problem with \(G = H^{-1}\) makes it look isotropic:
\[ G\approx H^{-1} \qquad\Rightarrow\qquad GH\approx I. \]
Equivalent view: change the distance scale.
Raw gradient change:
\[ \delta x\mapsto H\delta x. \]
Preconditioned change:
\[ \delta x\mapsto GH\delta x. \]
Nice units:
\[ GH\approx I. \]
Distance view:
same points, different unit lengths.
Choose a distance rule at each position:
\[ \|\delta x\|_{g(x)}^2 = (\delta x)^\top g(x)\,\delta x, \qquad g(x)\ \text{is symmetric positive definite (SPD)}. \]
The inverse \(G(x)=g(x)^{-1}\) scales gradients and noise.
If the distance rule changes with \(x\): Riemannian metric.
Infinitesimal distance:
\[ ds^2=dx^\top g(x)\,dx. \]
Inverse distance rule:
\[ G(x)=g(x)^{-1}. \]
Gradient in localized units:
\[ G(x)\nabla U(x). \]
Local quadratic prototype:
\[ g(x)\approx \nabla^2 U(x) \quad \text{if SPD.} \]
Otherwise use SPD proxies:
Fisher, Gauss–Newton, covariance, regularized Hessian.
Riemann, “Ueber die Hypothesen, welche der Geometrie zu Grunde liegen”, Abh. K. Ges. Wiss. Gottingen, 1868.
Insert the preconditioner into Langevin:
\[ Y = X+\frac h2G\nabla\log p(X)+\sqrt h\,G^{1/2}\xi, \qquad G\succ0. \]
\(G\) rescales drift and noise; the gradient is still measured at \(X\).
The preconditioner matches proposal length scales to local target length scales.
Drift: preconditioned by \(G\).
Noise: covariance \(hG\).
Target gradient: still \(\nabla\log p(X)\).
For the local quadratic, the next slide studies:
\[ GH. \]
Still one-step Langevin.
Matched local geometry would give \(G H(x)\approx I\).
Assume the local scales are unequal:
\[ H(x)=-\nabla^2\log\!p(x), \quad \theta_i\in\operatorname{eig}(G H(x)). \]
Fast direction fixes the safe step size:
\[ h=O(\theta_{\mathrm{fast}}^{-1}). \]
Slow gradient drift covers only
\[ h\theta_{\mathrm{slow}} = O\!\left( \frac{\theta_{\mathrm{slow}}}{\theta_{\mathrm{fast}}} \right) =O(1/\kappa), \qquad \kappa=\frac{\theta_{\mathrm{fast}}}{\theta_{\mathrm{slow}}} \quad\text{(condition number)}. \]
\(\rightarrow\) Efficiency decreases with condition number \(\kappa\).
Other problem: the noise term is sampled again every step.
Local one-step update:
\[ \delta x' \approx \delta x-\frac h2G H(x)\,\delta x +\sqrt h\,G^{1/2}\xi. \]
Slow deterministic fraction:
\[ \frac{|\Delta_{\mathrm{drift,slow}}|} {|\delta x_{\mathrm{slow}}|} \sim h\theta_{\mathrm{slow}} =O(1/\kappa). \]
In standardized slow units:
\[ z_{k+1} \approx \left(1-\frac h2\theta_{\mathrm{slow}}\right)z_k +\sqrt{h\theta_{\mathrm{slow}}}\,\eta_k. \]
\[ \eta_k\sim N(0,1), \qquad \eta_k \perp \eta_{k+1}. \]
Roberts et al., “Optimal scaling of discrete approximations to Langevin diffusions”, J. R. Stat. Soc. B, 1998.
Measure exploration by root mean square distance:
\[ \operatorname{RMS}_i(N) := \sqrt{\mathbb E[(z_{i,N}-z_{i,0})^2]}. \]
Independent noise kicks add like a random walk: \(N\) steps give a factor \(\sqrt N\).
With \(h\theta_{\mathrm{fast}}=O(1)\):
\[ \operatorname{RMS}_{\mathrm{fast}}(N) \sim\sqrt N, \qquad \operatorname{RMS}_{\mathrm{slow}}(N) \sim \sqrt{N/\kappa}. \]
Thus a target-scale slow refresh needs \(N=O(\kappa)\) local kicks.
To improve this, stop redrawing direction.
Roberts et al., “Optimal scaling of discrete approximations to Langevin diffusions”, J. R. Stat. Soc. B, 1998.
Fresh noise directions do not line up.
Keeping direction makes slow increments reinforce:
\[ (\text{slow coordinate }s) \qquad \Delta s_k := s_{k+1}-s_k . \]
\[ \left|\sum_{k=1}^L \Delta s_k\right| =O(L) \quad\text{instead of}\quad O(\sqrt L). \]
The gradient still bends the path back toward the density.
Slow increment:
\[ \Delta s_k=s_{k+1}-s_k. \]
Fresh directions:
\[ \mathbb E[\Delta s_k\Delta s_{k+1}]\approx0 \]
\[ \operatorname{RMS}\!\left(\sum_{k=1}^L\Delta s_k\right) \approx \sqrt L\,\sigma_s. \]
Keep direction:
\[ \Delta s_k\approx \Delta s \quad\Longrightarrow\quad \left|\sum_{k=1}^L \Delta s_k\right| \approx L|\Delta s|. \]
The coherent part lasts until the path is bent back.
Normal deviation \(n\perp s\):
\[ \begin{aligned} \lambda_n&:=\partial_n^2U(s,0)>0,\\ U(s,n)&\approx U_0+\frac12\lambda_n n^2. \end{aligned} \]
Normal steering:
\[ -\partial_n U(s,n)\approx-\lambda_n n. \]
Thus:
\[ s\text{-increments add,} \qquad n\text{ is corrected.} \]
Duane et al., “Hybrid Monte Carlo”, Phys. Lett. B, 1987.
A persistent path needs one more variable:
\[ x_t \quad\leadsto\quad (x_t,m_t). \]
We call \(m_t\) the momentum.
Momentum determines velocity through an SPD (symmetric positive definite) mass matrix \(M\):
\[ v_t=\dot x_t=M^{-1}m_t. \]
The target acts on momentum:
\[ \begin{aligned} \frac{d}{dt} \begin{bmatrix} x\\ m \end{bmatrix} &= \begin{bmatrix} M^{-1}m\\ -\nabla U(x) \end{bmatrix}, \qquad U(x)=-\log p_{\mathrm{tar}}(x). \end{aligned} \]
Duane et al., “Hybrid Monte Carlo”, Phys. Lett. B, 1987.
For each position \(x\), we add the same momentum distribution \(\varphi_M(m)\), because the Markov kernel depends on \(x\) and \(m\) now:
That means it can be factorized:
\[ \pi_{\mathrm{aug}}(x,m)=\pi(x)\,\varphi_M(m). \]
The \(x\)-marginal is unchanged:
\[ \int \pi_{\mathrm{aug}}(x,m)\,dm=\pi(x). \]
The chain may move in \((x,m)\); summaries/sample results still use only \(x\).
Normalize the added variable:
\[ \int \varphi_M(m)\,dm=1. \]
Then:
\[ \int \pi_{\mathrm{aug}}(x,m)\,dm = \pi(x)\int\varphi_M(m)\,dm = \pi(x). \]
Therefore, for any function of \(x\):
\[ \mathbb E_\pi[f(x)] = \mathbb E_{\pi_{\mathrm{aug}}}[f(x)]. \]
Momentum can be refreshed independently because it is not part of the requested marginal.
\[ m\sim \varphi_M \quad\text{given }x. \]
Duane et al., “Hybrid Monte Carlo”, Phys. Lett. B, 1987.
The previous slide chose the augmented density. The flow must now preserve it.
Write the momentum law as a potential:
\[ \varphi(m)\propto e^{-K(m)} \quad\Longrightarrow\quad H(x,m)=U(x)+K(m). \]
The compatible local coupling is:
\[ \dot x=\nabla_mK(m), \qquad \dot m=-\nabla_xU(x). \]
For Gaussian momentum,
\[ K(m)=\frac12m^\top M^{-1}m, \quad\Longrightarrow\quad \dot x=\nabla_mK(m)=M^{-1}m. \]
Energy cancellation:
\[ \frac{dH}{dt} = \nabla_xU^\top\nabla_mK - \nabla_mK^\top\nabla_xU =0. \]
Volume preservation:
\[ \nabla\cdot(\dot x,\dot m) = \nabla_x\cdot\nabla_mK - \nabla_m\cdot\nabla_xU =0. \]
Thus the exact continuous flow preserves:
\[ e^{-H(x,m)}\,dx\,dm \propto \pi(x)\varphi(m)\,dx\,dm. \]
This is the minimal compatible deterministic coupling; choosing \(M\) tunes its geometry.
\[ \text{momentum law}\Rightarrow K, \qquad K\Rightarrow \dot x. \]
Duane et al., “Hybrid Monte Carlo”, Phys. Lett. B, 1987.
Calculate slow-direction exploration speed:
\[ x_0\ \text{arbitrary}, \qquad \delta x=x-x_0, \qquad H_0=\nabla^2U(x_0). \]
Linearize the coupled dynamics around \(x_0\):
\[ \frac{d}{dt} \begin{bmatrix} \delta x\\ \delta m \end{bmatrix} = B_0 \begin{bmatrix} \delta x\\ \delta m \end{bmatrix}, \qquad B_0= \begin{bmatrix} 0 & M^{-1}\\ -H_0 & 0 \end{bmatrix}. \]
Eigenvalues of \(B_0\) identify fastest and slowest bends:
\[ M^{-1}H_0 u_i=\lambda_i u_i \quad\Longrightarrow\quad \rho_i^\pm\in\operatorname{eig}(B_0) =\pm i\sqrt{\lambda_i} =\pm i\omega_i. \]
Complex eigenvalues mean bending. Only the square-root ratio matters:
\[ \frac{\omega_{\max}}{\omega_{\min}} = \sqrt{\frac{\lambda_{\max}}{\lambda_{\min}}} =O(\sqrt{\kappa}), \qquad \text{vs. Langevin }O(\kappa). \]
Start from the left matrix:
\[ B_0 \begin{bmatrix}u\\w\end{bmatrix} = \rho \begin{bmatrix}u\\w\end{bmatrix}. \]
This means:
\[ \rho u=M^{-1}w, \qquad \rho w=-H_0u. \]
Eliminate \(w\):
\[ \rho^2u = -M^{-1}H_0u. \]
So for \(M^{-1}H_0u_i=\lambda_i u_i\):
\[ \rho_i^2=-\lambda_i, \qquad \rho_i^\pm=\pm i\sqrt{\lambda_i}. \]
Thus curvature rates become angular frequencies:
\[ \omega_i=|\operatorname{Im}\rho_i|=\sqrt{\lambda_i}. \]
Resolve the fastest bend:
\[ \epsilon\omega_{\max}=O(1). \]
Move through the slow bend:
\[ T\omega_{\min}=O(1). \]
Leapfrog steps:
\[ L=\frac{T}{\epsilon} = O\!\left( \frac{\omega_{\max}}{\omega_{\min}} \right). \]
Apers et al., “Hamiltonian Monte Carlo for efficient Gaussian sampling: long and random steps”, JMLR, 2024.
The square-root gain used one frozen local system:
\[ B_t\approx B_0, \qquad H_t\approx H_0. \]
While that is true, slow motion stays coherent:
\[ \left|\sum_{k=1}^L\Delta s_k\right|=O(L). \]
When curvature changes or modes mix, the increments stop reinforcing:
\[ \left|\sum_{k=1}^L\Delta s_k\right|=O(\sqrt L). \]
Time-varying linearization:
\[ B_t= \begin{bmatrix} 0 & M^{-1}\\ -H_t & 0 \end{bmatrix}, \qquad H_t=\nabla^2U(x_t). \]
Square-root scaling assumes:
\[ \int_0^T\!\|B_t-B_0\|\,dt \ll 1. \]
Equivalently, the force is still locally linear:
\[ r_t := \nabla U(x_t)-\nabla U(x_0)-H_0(x_t-x_0) \quad\text{small}. \]
Mode mixing starts when
\[ C_t=M^{-1/2}H_tM^{-1/2} \]
changes eigenvalues or rotates eigenvectors enough to create off-diagonal coupling.
NUTS detects one symptom:
\[ (x_t-x_0)^\top M^{-1}m_t\le0. \]
It detects backtracking, not linearization error. Either one can happen first.
Hoffman et al., “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo”, JMLR, 2014.
