Optimal Preconditioning and Fisher Adaptive Langevin Sampling

ArXiv: 2305.14442

🎯 Pitch

The optimal preconditioner for gradient-based MCMC is the inverse Fisher information matrix, not the target covariance as widely assumed—the paper analytically proves this by maximizing expected squared jump distance on the Langevin diffusion. FisherMALA learns this inverse Fisher online from the chain’s own gradient history in O(d²) per step, achieving orders-of-magnitude faster mixing than covariance-adaptive MALA on high-dimensional anisotropic targets.

---

1. Executive Summary

This paper derives an optimal preconditioning for the Langevin diffusion by analytically optimizing the expected squared jumped distance, yielding as the solution an inverse Fisher information covariance matrix $I^{-1} = E_{\pi(x)}[\nabla \log \pi(x) \nabla \log \pi(x)^\top]^{-1}$ — establishing that the optimal preconditioner is the inverse Fisher, not the target covariance as commonly believed. The method is applied to the Metropolis-adjusted Langevin algorithm (MALA) to produce FisherMALA, a computationally efficient adaptive MCMC scheme that learns the preconditioning online from score function increments (the Rao-Blackwellized differences between consecutive gradient evaluations) using an $O(d^2)$ square-root recursion. On anisotropic high-dimensional targets—including 100-dimensional Gaussian process and inhomogeneous Gaussian targets, plus Bayesian logistic regression on datasets up to 785 dimensions (MNIST)—FisherMALA achieves dramatic ESS improvements over standard adaptive MALA (e.g., minimum ESS of ~1785 vs. ~552 on the GP target), significantly outperforms the position-dependent Riemannian manifold MALA on Caravan and MNIST, and learns the optimal preconditioner orders of magnitude faster than conventional covariance-based adaptive MCMC, establishing that gradient-history adaptation is both more efficient and more principled than state-based adaptation when the target's score function structure can be exploited.

2. Context and Motivation

The Core Problem: How to Choose an Effective Preconditioner for Gradient-Based MCMC

The fundamental problem this paper tackles is how to automatically select a preconditioning matrix $A$ for the Langevin diffusion that maximally accelerates mixing in gradient-based MCMC samplers like MALA. The Langevin diffusion is a continuous-time stochastic process whose stationary distribution is the target $\pi(x)$:

$dx_t = \frac{1}{2}A \nabla \log \pi(x_t)dt + \sqrt{A} dB_t$

The matrix $A$ acts as a preconditioner — it rescales and rotates the gradient drift and noise injection so that the diffusion can efficiently explore all directions of the target distribution. Intuitively, if $\pi$ is highly anisotropic (stretched much more in some directions than others), a poorly chosen $A$ forces the sampler to take tiny steps to avoid instability in the most constrained direction, which makes exploration along the elongated directions excruciatingly slow. A well-chosen $A$ "whitens" or "symmetrizes" the target, equalizing the effective step sizes across all directions and enabling rapid mixing.

The practical challenge is that the optimal $A$ depends on the target's geometry, which is unknown a priori — we are trying to sample from $\pi$ precisely because we cannot analytically characterize it. This creates a chicken-and-egg problem: we need to precondition the sampler to efficiently explore $\pi$, but we need samples from $\pi$ to learn a good preconditioner. This paper provides a principled resolution by (1) analytically deriving what the optimal preconditioner should be and (2) designing an adaptive algorithm that learns it efficiently from the sampler's own gradient history.

Why This Problem Matters: The Curse of Anisotropy in High Dimensions

The importance of preconditioning cannot be overstated for practical MCMC. Consider a Bayesian posterior where different parameters have vastly different scales — for example, a logistic regression where some covariates vary on the order of $10^{-3}$ and others on the order of $10^3$. Without preconditioning, a single step size $\sigma^2$ applied isotropically must be tuned for the most restrictive direction, effectively freezing exploration in all other directions. The consequence is catastrophic in high dimensions:

• Effective sample size collapses dimension by dimension. A chain that mixes well in some coordinates may produce essentially no independent samples in others, making posterior summaries for those parameters meaningless.
• The problem compounds with dimensionality. As Section 5.2 demonstrates with the inhomogeneous Gaussian target (100 dimensions, standard deviations spanning two orders of magnitude from 0.01 to 1), simple MALA and even HMC with identity mass matrix achieve minimum ESS of ~3–4 on the stiffest dimensions while FisherMALA achieves ~1500 — a ~500× improvement. Multiplied across thousands of MCMC iterations, this is the difference between a usable posterior inference and a useless one.
• Real scientific applications demand it. The paper evaluates on Bayesian logistic regression with up to 785 parameters (MNIST digit classification). Modern Bayesian models — hierarchical models, spatial statistics, neural network posteriors — routinely involve thousands to millions of correlated parameters with wildly different scales. Adaptive preconditioning is not a luxury; it is a necessity for making gradient-based MCMC viable in these settings.

Beyond practical utility, there is a theoretical gap the paper identifies: the adaptive MCMC literature had largely converged on the heuristic that the optimal preconditioner for MALA is the target covariance matrix $\Sigma$ — the same as for random walk Metropolis (RWM). The paper challenges this directly (Section 3), showing through an explicit optimization that the optimal preconditioner is actually the inverse Fisher information matrix $I^{-1}$, which coincides with $\Sigma$ only for Gaussian targets. This is a significant conceptual correction with direct algorithmic consequences.

Prior Approaches and Where They Fall Short

The paper situates itself against three broad classes of existing methods:

1. Standard Adaptive MCMC (Covariance-Based)

The dominant paradigm, established by Haario et al. (2001), learns a global covariance matrix $\Sigma$ from the history of MCMC states and uses it to precondition the proposal:

$q(y_n | x_n) = \mathcal{N}\left(y_n \big| x_n + \frac{\sigma^2}{2} \Sigma_n \nabla \log \pi(x_n), \sigma^2 \Sigma_n\right)$

where $\Sigma_n$ is updated recursively from observed states $x_1, \ldots, x_n$. The recursion (Eq. 17 in the paper) is:

$\mu_n = \frac{n-1}{n}\mu_{n-1} + \frac{1}{n}x_n, \quad \Sigma_n = \frac{n-2}{n-1}\Sigma_{n-1} + \frac{1}{n}(x_n - \mu_{n-1})(x_n - \mu_{n-1})^\top$

Where it falls short: The paper identifies two fundamental problems:

• It is slow. Learning a full $d \times d$ covariance matrix from states requires observing the chain traverse the full extent of the target distribution. In high dimensions with strong anisotropy, the chain moves slowly precisely because it lacks good preconditioning — creating a vicious cycle where adaptation is slow because mixing is poor, and mixing is poor because adaptation is slow. Figure 1b shows FisherMALA learning the preconditioner to near-zero Frobenius error within ~2,000 iterations on the 100-dimensional GP target, while AdaMALA (covariance-based) requires essentially the full 20,000 burn-in iterations and still achieves worse accuracy.

• It is not what we actually want. The paper's Proposition 3 proves that the optimal $A^*$ is proportional to $I^{-1} = E_{\pi}[\nabla \log \pi \nabla \log \pi^\top]^{-1}$, not $\Sigma$. For non-Gaussian targets, these matrices differ. The paper's results on Bayesian logistic regression (Table 1, Caravan and MNIST) show that FisherMALA dramatically outperforms AdaMALA even when both are given comparable adaptation time, suggesting that learning the right quantity matters as much as learning quickly. The covariance heuristic — while natural — is simply optimizing the wrong objective.

The paper also notes a conceptual inconsistency: the covariance-based approach was originally developed for RWM, where the proposal lacks gradient information and $\Sigma$ is indeed the optimal preconditioner under optimal scaling theory (Roberts et al., 1997; Roberts and Rosenthal, 2001). Transplanting this heuristic to MALA, which does have gradient information, ignores that the gradient structure encodes additional geometric information about the target that can be exploited — specifically, the Fisher information matrix.

2. Position-Dependent (Riemannian) Preconditioning

The Riemannian manifold MALA (mMALA) of Girolami and Calderhead (2011) uses a position-dependent preconditioner $A(x)$, typically the inverse Fisher information matrix evaluated at the current point (or the negative inverse Hessian):

$A(x) = G(x)^{-1}, \quad \text{where } G(x) = -\nabla^2 \log \pi(x) \text{ or } E_{p(y|x)}[-\nabla_x^2 \log p(y|x)]$

This allows the proposal to adapt to local curvature, taking larger steps in flat regions and smaller steps in curved regions — analogous to second-order optimization in gradient descent.

Where it falls short:

• Cubic cost per iteration. Computing $A(x)$ requires evaluating or approximating the Hessian and performing a matrix decomposition (Cholesky or eigen-decomposition), costing $O(d^3)$ per iteration. This is prohibitive for $d > \sim 1000$, and even at $d = 100$ it is substantially slower than the $O(d^2)$ cost of FisherMALA. The paper notes this explicitly (Section 6): "position-dependent preconditioning methods have cubic cost and they require computationally demanding quantities like the Hessian matrix."

• Counterintuitively underperforms on some high-dimensional log-concave problems. Table 1 shows that on the Caravan dataset ($d = 87$), mMALA achieves minimum ESS of ~51 while FisherMALA achieves ~498 — nearly a 10× improvement. On MNIST ($d = 785$), the gap is even more dramatic: mMALA min ESS ~5.2 vs. FisherMALA ~440, an ~85× improvement. This is surprising because mMALA uses the exact negative inverse Hessian — ostensibly the perfect local preconditioner. The paper hypothesizes (Section 6) that "position-dependent preconditioning may be less effective in certain types of high-dimensional and log-concave problems." This is a significant empirical finding: a well-learned global preconditioner can outperform a theoretically-optimal local one, perhaps because the global preconditioner averages out noise in local curvature estimates, or because the computational savings of $O(d^2)$ enable more effective exploration within a fixed compute budget.

• Requires model-specific derivations. Computing $G(x)$ depends on the structure of $\log \pi(x)$ — for a Bayesian posterior, it involves the Hessian of the log-likelihood, which must be derived and implemented for each model. FisherMALA, by contrast, is model-agnostic: it only requires gradient evaluations $\nabla \log \pi(x)$, which are available for any differentiable target (e.g., through automatic differentiation).

3. Simple MALA and HMC Without Learned Preconditioning

Simple MALA uses an isotropic proposal $A = I_d$, adapting only the step size $\sigma^2$. HMC with identity mass matrix similarly lacks preconditioning. Both are standard baselines but fail catastrophically on anisotropic targets.

Where they fall short: The step size must accommodate the stiffest direction. On the inhomogeneous Gaussian target (Table 2, Appendix E.2), simple MALA achieves min ESS of ~3, HMC achieves ~42, while FisherMALA achieves ~1500. This is not a subtle difference — it means the simple samplers are producing effectively no independent samples in some dimensions after 20,000 iterations. Figure 2a visualizes this dramatically: ESS for FisherMALA is uniformly high (~1500–2300) across all 100 dimensions, while MALA and HMC show ESS that collapses from left (easy, large-variance dimensions) to right (stiff, small-variance dimensions), ending near zero for the most constrained coordinates.

The HMC comparison is particularly instructive because HMC is generally considered the gold-standard gradient-based sampler. The paper shows (Table 1, Caravan) that HMC with 10 leapfrog steps achieves min ESS of ~4.2 vs. FisherMALA's ~498 — a >100× gap. This is not because HMC is inherently worse than MALA, but because HMC without a learned mass matrix is just as crippled by anisotropy as simple MALA. The paper explicitly flags extending Fisher preconditioning to HMC's mass matrix as future work (Section 7), which could combine HMC's superior exploration in well-conditioned problems with FisherMALA's ability to handle anisotropy.

How This Paper Positions Itself

The paper constructs its contribution along three axes that collectively fill a gap in the literature:

Axis 1: Principled optimality, not heuristics. Previous adaptive MCMC methods justified covariance-based preconditioning through analogy to RWM optimal scaling or through empirical demonstration, not through direct optimization of a criterion on the Langevin diffusion itself. The paper's Proposition 3 provides the first (to my knowledge) analytic derivation that the optimal $A$ minimizes the expected squared jumped distance subject to a trace constraint, yielding $A^* \propto I^{-1}$. This connects to the expected squared jumped distance criterion of Pasarica and Gelman (2010) but applies it in a novel way — not to tune a scalar step size, but to optimize a full preconditioning matrix.

The optimization criterion itself is intuitively appealing: maximize the discretization step size $\delta$ subject to a bounded expected squared jump distance $J(\delta, A) \leq \epsilon$. Since $J(\delta, A)$ increases with $\delta$ (larger steps mean larger jumps), the maximum allowable $\delta$ at the constraint boundary satisfies $\min_A J(\delta^*, A) = \epsilon$. This means the optimal $A^*$ is the one that minimizes the expected squared jump distance for a given $\delta$ — which, by the trace of the covariance in Proposition 2, reduces to minimizing $\text{tr}(A I A)$ subject to $\text{tr}(A) = c$. The solution is $A^* \propto I^{-1}$. This is a clean, self-contained derivation that does not rely on asymptotic arguments or Gaussian assumptions.

Axis 2: Fast, online $O(d^2)$ adaptation from gradient history. Recognizing $I^{-1}$ as the target is one thing; learning it efficiently during MCMC is another. The paper's key algorithmic insight is that the gradient history — specifically, score function increments $s^\delta_n = s(x_{n+1}) - s(x_n)$ — provides a much richer and more stable learning signal than the state history. This is because:

• Gradients encode local geometric information about $\pi$ (curvature, anisotropy) directly, whereas states only encode it indirectly through their covariance.
• Score function increments tend to be centered near zero even during the transient phase (before convergence to $\pi$), whereas raw score vectors $s_n = \nabla \log \pi(x_n)$ can have large non-zero means during transience, biasing the empirical Fisher estimate through the term $E[s_n]E[s_n]^\top$.
• The Rao-Blackwellized signal $s^\delta_n = \sqrt{\alpha(x_n, y_n)}(s(y_n) - s(x_n))$ (Eq. 16) reduces variance by marginalizing over the accept/reject decision, and ensures that even rejected proposals contribute to learning — a form of data efficiency that purely state-based adaptation lacks.

The $O(d^2)$ square-root recursion using Potter's method (Proposition 4, dating to Kalman filtering literature) is an elegant way to maintain a matrix square root $R_n$ such that $R_n R_n^\top = A_n$ without ever explicitly inverting or decomposing $A_n$. The recursion:

$R_n = R_{n-1} - r_n \frac{(R_{n-1}\phi_n)\phi_n^\top}{1 + \phi_n^\top \phi_n}, \quad \phi_n = R_{n-1}^\top s_n, \quad r_n = \frac{1}{1 + \sqrt{\frac{1}{1 + \phi_n^\top \phi_n}}}$

costs $O(d^2)$ (two matrix-vector products and one outer product per iteration) and directly produces the square root needed for sampling $\eta \sim \mathcal{N}(0, I_d)$ and transforming to $R\eta \sim \mathcal{N}(0, A)$. This is substantially cheaper than the $O(d^3)$ required to invert and decompose a full matrix at each iteration.

Axis 3: Empirical demonstration that gradient-based adaptation dominates state-based adaptation. The paper does not just argue theoretically — it provides extensive evidence that FisherMALA's gradient-based learning dramatically outperforms AdaMALA's state-based learning across diverse targets. The Frobenius norm trajectories (Figures 1a-b, 2c) show FisherMALA converging to the optimal preconditioner in a few thousand iterations versus tens of thousands for AdaMALA. The ESS comparisons (Table 1, Tables 2-6) show FisherMALA achieving uniformly high efficiency across dimensions while AdaMALA and other methods collapse on the stiffest coordinates. And the comparison with mMALA — which uses the exact optimal local preconditioner — shows that a well-learned global preconditioner can be more effective than an exact local one, at least in the high-dimensional log-concave settings tested.

The Broader Context: Connecting to Optimal Scaling and Fisher Information

The paper carefully connects its result to two established literatures, which strengthens its positioning:

Connection to optimal scaling: The optimal scaling literature for RWM (Roberts et al., 1997; Roberts and Rosenthal, 2001) established that for product-form targets $\pi(x) = \prod_i f(x_i)$ as $d \to \infty$, the optimal proposal variance is $\sigma^2 = 2.38 / J$ where $J = E_f[(d \log f(x)/dx)^2]$ is the univariate Fisher information. The optimal preconditioner is then $I^{-1} = J^{-1} I_d$ — again the inverse Fisher. The paper's Proposition 3 recovers this as a special case (when $I$ is diagonal and homogeneous) while generalizing to arbitrary covariance structure. This places FisherMALA in a direct lineage with the foundational results that established the $0.234$ acceptance rate rule for RWM and the $0.574$ rule for MALA.

Connection to classical Fisher information: The paper distinguishes its matrix $I = E_{\pi(x)}[\nabla \log \pi(x) \nabla \log \pi(x)^\top]$ from the classical Fisher information $G(\theta) = E_{p(Y|\theta)}[\nabla_\theta \log p(Y|\theta) \nabla_\theta \log p(Y|\theta)^\top]$. The classical Fisher averages over data given parameters (frequentist); the paper's $I$ averages over parameters given data (Bayesian). This distinction matters because $I$ is a constant matrix for a given posterior $\pi(\theta | Y)$ — it does not depend on $\theta$, making it suitable as a global preconditioner — while $G(\theta)$ varies with $\theta$, making it suitable only for local preconditioning (as in mMALA). The paper also notes that under regularity conditions, $I = -E_{\pi(x)}[\nabla^2 \log \pi(x)]$, connecting the Fisher matrix to the expected Hessian and reinforcing the intuition that $I^{-1}$ captures the target's average curvature.

What the Paper Does NOT Address (By Design)

The paper is explicit about its scope and leaves several important questions open:

• Bias of the score increment estimator: The paper uses $s^\delta_n$ to estimate the Fisher matrix but does not prove that this estimator is unbiased or consistent. Section 7 flags this as a topic for future work. Intuitively, since $s^\delta_n$ involves the difference of gradient evaluations at consecutive (correlated) MCMC states, the outer products $s^\delta_n (s^\delta_n)^\top$ may not average to $I$ in the same way that independent draws would. The empirical results suggest any bias is small enough not to matter, but theoretical guarantees are absent.

• Extension to HMC: Learning the mass matrix for HMC is more complex because both the mass matrix and its inverse appear in the leapfrog integrator. The paper's square-root recursion provides $R_n$ (such that $A_n = R_n R_n^\top$) but HMC would also need $A_n^{-1} = R_n^{-\top} R_n^{-1}$ for the kinetic energy term. This is flagged as future work in Section 7.

• Multiple chains and ensemble methods: The paper adapts from a single chain. There is recent work on ensemble Langevin methods (Garbuno-Inigo et al., 2020) that could potentially accelerate Fisher estimation by pooling gradient information across interacting chains. The paper mentions this as a natural extension.

• Theoretical convergence guarantees: The paper invokes the standard adaptive MCMC conditions (diminishing adaptation, containment) but does not prove that FisherMALA satisfies them. This is consistent with the broader adaptive MCMC literature, where proving ergodicity for specific algorithms is often deferred to separate theoretical papers. The empirical evidence of convergence (Frobenius norm decay, stable ESS estimates across random restarts) provides practical validation but not mathematical proof.

In summary, this paper addresses a clear and important gap: the lack of a principled, computationally efficient adaptive preconditioning method for gradient-based MCMC that learns the right quantity (inverse Fisher, not covariance) from the right signal (gradient increments, not states). It positions itself as unifying optimal scaling theory, Fisher information geometry, and practical adaptive MCMC into a single algorithm that is both theoretically motivated and demonstrably superior to existing approaches.

3. Technical Approach

3.1 Reader Orientation

This is primarily a theoretical derivation combined with algorithmic design paper. The system being built is an adaptive Markov chain Monte Carlo sampler — specifically, a Metropolis-adjusted Langevin algorithm (MALA) that learns its own preconditioning matrix online as it runs. The paper solves the problem of "what is the best preconditioner for gradient-based MCMC, and how can we learn it efficiently during sampling?" by (1) analytically deriving that the optimal preconditioner is the inverse Fisher information matrix $I^{-1} = E_\pi[\nabla \log \pi \nabla \log \pi^\top]^{-1}$ — not the target covariance $\Sigma$ as commonly believed — and (2) designing an $O(d^2)$ per-iteration recursive algorithm that estimates this $I^{-1}$ from the history of score function evaluations, using a technique adapted from Kalman filtering to maintain a square root of the estimate for efficient sampling.

3.2 Big-Picture Architecture (Diagram in Words)

The FisherMALA system has five major components that interact at each iteration:

1. Target distribution interface — the external black-box that provides $\log \pi(x)$ and $\nabla \log \pi(x)$ for any input $x$. This is whatever the user wants to sample from (a Bayesian posterior, an energy-based model, etc.). FisherMALA only requires gradient evaluations; no Hessian or other structure is needed.

2. Preconditioned proposal mechanism (Equation 7) — given the current state $x_n$, a learned square root matrix $R_n$, and a global step size $\sigma^2$, this generates a candidate state $y_n$ by taking a gradient step with isotropic noise transformed by $R_n$: $y_n = x_n + \frac{\sigma^2_R}{2} R_n (R_n^\top \nabla \log \pi(x_n)) + \sigma_R R_n \eta, \quad \eta \sim \mathcal{N}(0, I_d)$ where $\sigma^2_R = \sigma^2 / \frac{1}{d}\text{tr}(R_n R_n^\top)$ normalizes out the scale of $R_n$.

3. Metropolis-Hastings accept/reject gate — computes the acceptance probability $\alpha(x_n, y_n)$ using Proposition 1's simplified form that avoids $A^{-1}$ (only requires $\nabla \log \pi$ evaluations, not matrix inverses). With probability $\alpha$, the chain advances to $y_n$; otherwise it stays at $x_n$.

4. Score increment adaptation signal (Equation 16) — computes $s^\delta_n = \sqrt{\alpha(x_n, y_n)}(s(y_n) - s(x_n))$, the Rao-Blackwellized difference between the gradient at the proposed state and the gradient at the current state. This signal drives the entire preconditioner learning.

5. Square-root Fisher matrix estimator (Proposition 4) — takes the adaptation signal $s^\delta_n$ and recursively updates a square root matrix $R_n$ such that $R_n R_n^\top \approx I^{-1}$ (up to a positive scale factor). The recursion is $O(d^2)$ per iteration and directly produces the matrix needed for sampling in component 2.

Information flows in a tight loop at each iteration: from the current state $x_n$, component 2 proposes $y_n$ using the current $R_n$ and $\sigma^2$ → component 3 decides whether to accept → component 4 computes the adaptation signal from the gradient evaluations → component 5 updates $R_n$ → component 2 uses the updated $R_n$ for the next proposal. The step size $\sigma^2$ is independently adapted toward an acceptance rate of $0.574$ using a Robbins-Monro-style update $\sigma^2 \leftarrow \sigma^2[1 + \rho_n(\alpha(x_n, y_n) - \alpha^*)]$.

3.3 Roadmap for the Deep Dive

• First, the formal derivation of the optimal preconditioner (Proposition 3) via expected squared jumped distance optimization. This establishes what we want to learn and why it is $I^{-1}$ rather than $\Sigma$.
• Second, the proposal parametrization (Equation 7) and the normalized step size scheme that decouples adapting $A$ from adapting $\sigma^2$. This introduces the critical invariance property that makes adaptation stable.
• Third, the empirical Fisher estimation recursion (Equations 10–11) and its square-root form (Proposition 4, Equations 12–13). This is the computational engine of the method — how we go from gradient observations to a usable sampling matrix in $O(d^2)$.
• Fourth, the adaptation signal design — why score function increments $s^\delta_n$ are used instead of raw score vectors $s_n$, the Rao-Blackwellization step (Equation 16), and why this signal is more robust during the transient phase.
• Fifth, the complete algorithm (Algorithm 1) including initialization, the acceptance rate adaptation, and the integration with the Metropolis-Hastings accept/reject step.

3.4 Detailed, Sentence-Based Technical Breakdown

This is a theoretical derivation plus algorithmic design paper whose core idea is that the optimal preconditioner for a Langevin diffusion can be derived analytically by maximizing the discretization step size subject to a bounded expected squared jump distance, yielding the inverse Fisher information matrix — and that this matrix can be learned efficiently online from gradient evaluations using an $O(d^2)$ square-root recursion adapted from classical Kalman filtering.

---

The Optimal Preconditioning Derivation: Why $I^{-1}$ Instead of $\Sigma$

The paper's first major contribution is a principled answer to the question: if we could choose any symmetric positive definite matrix $A$ as our preconditioner, which one should we choose? The derivation proceeds in three stages.

Stage 1: Express the expected squared jumped distance.

Consider the continuous-time Langevin diffusion from Equation (2):

$dx_t = \frac{1}{2}A \nabla \log \pi(x_t)dt + \sqrt{A} dB_t$

Discretize this with a small finite step size $\delta = \sigma^2 > 0$ (the same $\sigma^2$ that appears in the MALA proposal) using the Euler-Maruyama scheme to obtain the unadjusted (rejection-free) Langevin update:

$x_{t+\delta} - x_t = \frac{\delta}{2}A \nabla \log \pi(x_t) + \sqrt{A}(B_{t+\delta} - B_t)$

where $B_{t+\delta} - B_t \sim \mathcal{N}(0, \delta I_d)$ is the Brownian increment over a time interval of length $\delta$, and $\sqrt{A}$ is any matrix satisfying $\sqrt{A}\sqrt{A}^\top = A$.

Proposition 2 computes the covariance of the displacement vector $x_{t+\delta} - x_t$, assuming that $x_t$ is already distributed according to the stationary distribution $\pi$:

$E[(x_{t+\delta} - x_t)(x_{t+\delta} - x_t)^\top] = \frac{\delta^2}{4} A E_{\pi(x_t)}[\nabla \log \pi(x_t) \nabla \log \pi(x_t)^\top] A + \delta A$

where the expectation is taken both over $x_t \sim \pi$ and over the independent Brownian increment. This expression has two terms: a drift-driven term $\frac{\delta^2}{4} A I A$ (where $I$ is the Fisher matrix $E_\pi[\nabla \log \pi \nabla \log \pi^\top]$), which captures how the deterministic gradient component spreads the chain, and a diffusion-driven term $\delta A$, which captures how the injected noise spreads the chain.

The expected squared jumped distance $J(\delta, A) = E[\|x_{t+\delta} - x_t\|^2]$ is then the trace of this covariance. Since the trace is linear and $E[\|v\|^2] = \text{tr}(E[vv^\top])$ for any random vector $v$, we obtain:

$J(\delta, A) = \text{tr}\left(\frac{\delta^2}{4} A I A + \delta A\right) = \frac{\delta^2}{4} \text{tr}(A I A) + \delta \text{tr}(A)$

What this expression computes: for a given preconditioner $A$ and step size $\delta$, it gives the average squared Euclidean distance between consecutive states of the unadjusted Langevin sampler when the chain is at stationarity. Larger $J$ means the chain is making bigger jumps, which is generally good for exploration — but it also means larger discretization error, since the Euler-Maruyama approximation becomes less accurate with larger step sizes.

Why this form matters: the expression separates the influence of $\delta$ (step size) and $A$ (preconditioner shape) into an interaction term $\text{tr}(A I A)$ and a simple scaling term $\text{tr}(A)$. This separation is what makes the subsequent optimization analytically tractable. Note that $A$ appears quadratically in the drift term but linearly in the diffusion term — this asymmetry is what drives the optimal solution away from the identity matrix when $I$ is not a multiple of the identity.

Stage 2: Set up the constrained optimization problem.

The key insight is that a good preconditioner should allow a larger step size $\delta$ without causing excessive discretization error. Since discretization error grows with the jump distance, we can formalize this by imposing an upper bound $J(\delta, A) \leq \epsilon$ for some small tolerance $\epsilon > 0$, and then choosing $A$ to maximize the allowable $\delta$ subject to this constraint.

Because $J(\delta, A)$ monotonically increases with $\delta$ (both terms are non-negative and increasing in $\delta$), the maximum $\delta^*$ that satisfies the constraint is the one that achieves equality: $J(\delta^*, A) = \epsilon$. If we can find an $A$ that makes $J(\delta, A)$ smaller for a given $\delta$, then the equality $J(\delta^*, A) = \epsilon$ will be reached at a larger $\delta^*$.

Therefore, maximizing $\delta^*$ over $A$ is equivalent to minimizing $J(\delta, A)$ over $A$ for any fixed $\delta$, subject to a normalization constraint that prevents $A$ from simply being scaled to zero (which would trivially minimize $J$ but make the chain immobile).

Proposition 3 imposes the constraint $\text{tr}(A) = c$ for some constant $c > 0$. This fixes the average eigenvalue of $A$ — the "total preconditioning power" — while leaving the shape (distribution of eigenvalues among directions) free for optimization.

Under this constraint, minimizing $J(\delta, A)$ reduces to minimizing the quadratic term $\text{tr}(A I A)$, since $\delta \text{tr}(A) = \delta c$ is constant:

$\min_{A \succ 0, \text{tr}(A) = c} \text{tr}(A I A)$

where $A \succ 0$ denotes that $A$ is symmetric positive definite, and $I$ is the Fisher information matrix (also symmetric positive definite since it is a covariance matrix of gradient vectors).

Stage 3: Solve the optimization via Lagrange multipliers.

Form the Lagrangian:

$\mathcal{L}(A, \lambda) = \text{tr}(A I A) - \lambda(\text{tr}(A) - c)$

where $\lambda$ is the Lagrange multiplier for the trace constraint. Taking the matrix derivative with respect to $A$ (using the identities $\frac{\partial}{\partial X}\text{tr}(XBX) = X^\top B^\top + B^\top X^\top$ and $\frac{\partial}{\partial X}\text{tr}(X) = I_d$) and setting to zero yields:

$A^\top I + I A^\top = \lambda I_d$

Since both $A$ and $I$ are symmetric, $A^\top = A$ and $I^\top = I$, giving:

$A I + I A = \lambda I_d$

This is a Sylvester equation in $A$. Because $I$ is invertible (all eigenvalues $\mu_i$ satisfy $0 < \mu_i < \infty$ by assumption), the unique solution is $A = \frac{\lambda}{2} I^{-1}$. Substituting back into the constraint $\text{tr}(A) = c$ determines $\lambda$:

$\text{tr}\left(\frac{\lambda}{2} I^{-1}\right) = \frac{\lambda}{2} \sum_{i=1}^d \frac{1}{\mu_i} = c \implies \lambda = \frac{2c}{\sum_{i=1}^d \frac{1}{\mu_i}}$

Therefore the optimal preconditioner is:

$A^* = \frac{c}{\sum_{i=1}^d \frac{1}{\mu_i}} I^{-1} = k I^{-1}, \quad k = \frac{c}{\sum_{i=1}^d \frac{1}{\mu_i}}$

What this computes: the optimal $A^*$ is proportional to the inverse of the Fisher information matrix $I$, scaled so that its trace equals $c$. The proportionality constant $k$ is the harmonic mean of the eigenvalues $\mu_i$ of $I$, divided by $d$ and multiplied by $c$.

Why this is the global minimum and not a maximum or saddle point: the paper verifies this by comparing the loss at $A^*$ with the loss at the isotropic choice $A = \frac{c}{d} I_d$ (which also satisfies the trace constraint). Using the Cauchy-Schwarz inequality $d^2 = (\sum_i \frac{\sqrt{\mu_i}}{\sqrt{\mu_i}})^2 \leq (\sum_i \mu_i)(\sum_i \frac{1}{\mu_i})$, the loss for the isotropic choice is $\frac{c^2}{d^2} \sum_i \mu_i \geq \frac{c^2}{\sum_i 1/\mu_i}$, which is strictly greater than the optimal loss $\frac{c^2}{\sum_i 1/\mu_i}$ unless all $\mu_i$ are equal (in which case $I$ is a multiple of the identity and $I^{-1}$ is also isotropic).

The key insight: $A^*$ sends more noise and larger gradient steps into directions where the Fisher information is small (flat directions of $\log \pi$), and less noise and smaller steps into directions where the Fisher information is large (curved directions). This is exactly the right behavior: in flat directions, large moves don't change the log-density much and are likely to be accepted; in curved directions, small moves are needed to avoid steep drops in density that would cause rejection. The harmonic mean scaling $k$ ensures that the overall jump magnitude stays controlled.

Connection to the Gaussian special case: if $\pi(x) = \mathcal{N}(\mu, \Sigma)$, then $\nabla \log \pi(x) = -\Sigma^{-1}(x - \mu)$, and $I = E_\pi[\Sigma^{-1}(x - \mu)(x - \mu)^\top \Sigma^{-1}] = \Sigma^{-1} \Sigma \Sigma^{-1} = \Sigma^{-1}$. Hence $I^{-1} = \Sigma$ — for a Gaussian target, the optimal preconditioner is precisely the covariance matrix. This explains why the covariance heuristic works well for near-Gaussian posteriors (e.g., large-sample Bayesian asymptotics) but can be suboptimal for non-Gaussian targets where $I^{-1} \neq \Sigma$.

---

The Proposal Parametrization and Normalized Step Size

With the theoretical target $A^* \propto I^{-1}$ established, the next design choice is how to parametrize the MALA proposal to make the adaptation of $A$ well-behaved and decoupled from the adaptation of the global step size $\sigma^2$.

The standard preconditioned MALA proposal (Equation 3) is:

$q(y_n | x_n) = \mathcal{N}\left(y_n \big| x_n + \frac{\sigma^2}{2} A \nabla \log \pi(x_n), \sigma^2 A\right)$

Problem with the naive parametrization: if we adapt $A$ and $\sigma^2$ simultaneously using this form, the two adaptations interfere. Specifically, if $A$ is scaled by a factor $k$, it simultaneously changes the gradient step magnitude (through $\frac{\sigma^2}{2} A \nabla \log \pi$) and the noise covariance (through $\sigma^2 A$), while the true "informativeness" of $A$ as a preconditioner should depend only on its shape (the relative sizes of its eigenvalues), not its scale.

To fix this, Equation (7) reparametrizes the proposal as:

$q(y_n | x_n) = \mathcal{N}\left(y_n \big| x_n + \frac{\sigma^2}{2} \frac{d}{\text{tr}(A)} A \nabla \log \pi(x_n), \sigma^2 \frac{1}{\frac{1}{d}\text{tr}(A)} A\right)$

where the matrix actually used as the preconditioner is $\frac{A}{\frac{1}{d}\text{tr}(A)}$ — $A$ divided by its average eigenvalue.

What this reparametrization achieves: the matrix $\frac{A}{\frac{1}{d}\text{tr}(A)}$ has an average eigenvalue of 1, regardless of the original scale of $A$. This means:

• The proposal is invariant to any positive scaling of $A$: replacing $A$ with $kA$ (for any $k > 0$) leaves the proposal unchanged, because $\frac{kA}{\frac{1}{d}\text{tr}(kA)} = \frac{kA}{k \frac{1}{d}\text{tr}(A)} = \frac{A}{\frac{1}{d}\text{tr}(A)}$.
• The step size $\sigma^2$ now controls the overall magnitude of both the gradient step and the noise, independent of $A$'s shape.
• Adaptation of $A$ can focus purely on learning its orientation and relative eigenvalue ratios (the shape of the preconditioning ellipsoid), while adaptation of $\sigma^2$ handles global scaling.

Implementation detail: the normalized step size $\sigma^2_R = \sigma^2 / \frac{1}{d}\text{tr}(R R^\top)$ is precomputed at each iteration. Since $\text{tr}(R R^\top) = \sum_{i,j} R_{ij}^2$ (the sum of squared entries of $R$), this is an $O(d^2)$ computation that can be done efficiently as sum(R ◦ R) (elementwise square and sum).

Connection to the optimality derivation: the trace constraint $\text{tr}(A) = c$ in Proposition 3 is effectively the normalization step. Since the proposal is invariant to scaling of $A$, the algorithm learns $A$ up to an arbitrary positive scale, and the $\frac{1}{d}\text{tr}(A)$ normalization extracts the shape information. This means the square root matrix $R_n$ maintained by the algorithm need only match $I^{-1}$ up to a multiplicative constant — exactly what the derivation predicts is sufficient.

---

The Empirical Fisher Estimation Recursion

Having established that the target preconditioner is proportional to $I^{-1}$, the next challenge is learning $I^{-1}$ online from gradient observations during MCMC. The paper constructs this in three layers of increasing sophistication: the batch empirical Fisher, the direct inverse recursion, and the square-root recursion.

Layer 1: The batch empirical Fisher estimator.

Given score vectors $s_i = \nabla \log \pi(x_i)$ for $i = 1, \ldots, n$ (where $x_i$ are the states visited by the MCMC chain), the empirical Fisher matrix is:

$\hat{I}_n = \frac{1}{n} \sum_{i=1}^n s_i s_i^\top + \frac{\lambda}{n} I_d$

where $\lambda > 0$ is a fixed damping parameter (set to $\lambda = 10$ in all experiments) and $I_d$ is the $d \times d$ identity matrix.

What this computes: the sample covariance of the gradient vectors (under the assumption that $E_\pi[s] = 0$, which holds at stationarity because $\int \nabla \log \pi(x) \pi(x) dx = \int \nabla \pi(x) dx = 0$ under mild regularity conditions), plus a regularization term that shrinks toward the identity.

Why the damping term: as $n \to \infty$, $\frac{\lambda}{n} I_d \to 0$, so $\hat{I}_n$ is a consistent estimator of $I$. But for finite $n$ (especially small $n$ at the start of adaptation), $\frac{1}{n} \sum_{i=1}^n s_i s_i^\top$ can be rank-deficient (if fewer than $d$ independent gradient directions have been observed) or have very small eigenvalues in unexplored directions. The term $\frac{\lambda}{n} I_d$ acts as a Tikhonov regularizer (like ridge regression), ensuring all eigenvalues are at least $\lambda/n > 0$, which guarantees that $\hat{I}_n$ and its inverse are well-conditioned. This is critical for stability in early iterations.

Layer 2: The direct inverse recursion (avoiding explicit inversion).

The preconditioner we actually want is $A_n \propto \hat{I}_n^{-1}$. From the definition of $\hat{I}_n$:

$A_n \propto \left(\frac{1}{n} \sum_{i=1}^n s_i s_i^\top + \frac{\lambda}{n} I_d\right)^{-1} = n \left(\sum_{i=1}^n s_i s_i^\top + \lambda I_d\right)^{-1}$

Since scaling $A_n$ by any positive constant is irrelevant (due to the proposal's invariance to scaling, discussed above), the factor $n$ can be dropped, yielding the working definition:

$A_n = \left(\sum_{i=1}^n s_i s_i^\top + \lambda I_d\right)^{-1}$

The key algorithmic insight is that this can be updated recursively in $O(d^2)$ per iteration using the Sherman-Morrison-Woodbury formula, without ever explicitly forming or inverting a $d \times d$ sum.

Initialization (n = 1):

$A_1 = (s_1 s_1^\top + \lambda I_d)^{-1}$

Applying the Woodbury matrix identity $(B + uv^\top)^{-1} = B^{-1} - \frac{B^{-1} u v^\top B^{-1}}{1 + v^\top B^{-1} u}$ with $B = \lambda I_d$, $u = s_1$, $v = s_1$:

$A_1 = \frac{1}{\lambda} \left(I_d - \frac{s_1 s_1^\top}{\lambda + s_1^\top s_1}\right)$

What this computes: the initial inverse Fisher estimate shrinks the identity matrix in the direction of the first observed gradient $s_1$. The amount of shrinkage depends on the norm $\|s_1\|^2$: if $\|s_1\|^2$ is large (a direction of high curvature), the subtraction is larger; if small (a flat direction), the identity component dominates. This is already approximating the behavior that $A_n$ should have small eigenvalues in high-curvature directions and large eigenvalues in flat directions.

Iteration (n > 1):

Given $A_{n-1} = (\sum_{i=1}^{n-1} s_i s_i^\top + \lambda I_d)^{-1}$, when a new gradient $s_n$ arrives:

$A_n = \left(A_{n-1}^{-1} + s_n s_n^\top\right)^{-1}$

Applying Woodbury again with $B = A_{n-1}^{-1}$, $u = s_n$, $v = s_n$:

$A_n = A_{n-1} - \frac{A_{n-1} s_n s_n^\top A_{n-1}}{1 + s_n^\top A_{n-1} s_n}$

What this computes, operationally: starting from the previous estimate $A_{n-1}$, the recursion performs a rank-one downdate that reduces the eigenvalues of $A_n$ in the direction of the new gradient $s_n$. The denominator $1 + s_n^\top A_{n-1} s_n$ controls how much reduction occurs: if $s_n^\top A_{n-1} s_n$ is large (meaning the new gradient aligns with directions where $A_{n-1}$ already has large eigenvalues — i.e., directions that have been observed to be flat), the downdate is stronger; if it is small (the new gradient aligns with directions where $A_{n-1}$ already has small eigenvalues — i.e., already-identified high-curvature directions), the downdate is mild.

Why this form works: the recursion maintains $A_n$ as the exact inverse of the accumulated sum-plus-regularization $\sum_{i=1}^n s_i s_i^\top + \lambda I_d$. At convergence (as $n \to \infty$), $\frac{1}{n} A_n^{-1} \to I$, so $\frac{1}{n} A_n \to I^{-1}$. The scaling factor $1/n$ is irrelevant for preconditioning (only the shape matters), so running the recursion as-is produces an estimate that converges to a matrix proportional to $I^{-1}$.

Cost: each iteration requires computing $A_{n-1} s_n$ (an $O(d^2)$ matrix-vector product), the scalar $1 + s_n^\top(A_{n-1} s_n)$ (an $O(d)$ dot product), the vector $\frac{A_{n-1} s_n}{1 + s_n^\top A_{n-1} s_n}$ ($O(d)$ scaling), and the outer product of this vector with $(A_{n-1} s_n)^\top$ (an $O(d^2)$ outer product). The total is $O(d^2)$, dominated by the matrix-vector product and the outer product.

Layer 3: The square-root recursion (Proposition 4).

The recursion for $A_n$ gives us the inverse Fisher matrix, but for sampling we need a square root matrix $R_n$ such that $R_n R_n^\top = A_n$. Specifically, to generate $\eta \sim \mathcal{N}(0, I_d)$ and transform it to $R_n \eta \sim \mathcal{N}(0, A_n)$, we need $R_n$ (or a matrix $R_n$ satisfying $R_n R_n^\top = A_n$). We could compute $R_n$ from $A_n$ via a Cholesky decomposition at each iteration, but that costs $O(d^3)$ — too expensive to do at every MCMC step.

The paper's solution is to directly maintain a square root matrix $R_n$ recursively, adapting a technique from the early Kalman filtering literature (Potter and Stern, 1963; Bierman, 1977). The recursion updates $R_n$ in $O(d^2)$ without ever forming $A_n$ explicitly.

The core mathematical building block is Lemma 1 (stated and proved in Appendix A.4):

Lemma 1: Suppose we have a positive definite matrix of the form $I_d - z z^\top$ where $z \in \mathbb{R}^d$ and $z^\top z \leq 1$. Then a square root matrix $R$ satisfying $R R^\top = I_d - z z^\top$ is $R = I_d - r z z^\top$, where $r = \frac{1}{1 + \sqrt{1 - z^\top z}}$.

What this means in plain language: if we need a square root of a matrix that is the identity minus a rank-one matrix $z z^\top$, the square root has the same form (identity minus a scaled version of $z z^\top$) with a specific scaling factor $r \in (0, 1]$ that depends on the norm of $z$. This is a special case of the more general fact that the square root of a rank-one perturbation of the identity is another rank-one perturbation of the identity, with a modified coefficient.

Now consider how $A_n$ evolves. The initialization from Equation (10) is:

$A_1 = \frac{1}{\lambda} \left(I_d - \frac{s_1 s_1^\top}{\lambda + s_1^\top s_1}\right)$

This is of the form $\frac{1}{\lambda}(I_d - z z^\top)$ with $z = \frac{s_1}{\sqrt{\lambda + s_1^\top s_1}}$. Since $z^\top z = \frac{s_1^\top s_1}{\lambda + s_1^\top s_1} < 1$, Lemma 1 applies, giving the initialization (Equation 12):

$R_1 = \frac{1}{\sqrt{\lambda}} \left(I_d - r_1 \frac{s_1 s_1^\top}{\lambda + s_1^\top s_1}\right), \quad r_1 = \frac{1}{1 + \sqrt{\frac{\lambda}{\lambda + s_1^\top s_1}}}$

What this computes: the square root of the initial inverse Fisher estimate. The factor $1/\sqrt{\lambda}$ ensures $R_1 R_1^\top = A_1$. The scaling constant $r_1$ is always between 0.5 and 1: when $s_1^\top s_1$ is small (the first gradient is near zero, implying a flat region), $\frac{\lambda}{\lambda + s_1^\top s_1} \approx 1$, so $r_1 \approx 1/2$, and the subtraction in $R_1$ is maximal. When $s_1^\top s_1$ is large (high curvature at initialization), $r_1 \approx 1$, and the subtraction is minimal — $R_1$ stays close to $1/\sqrt{\lambda} I_d$.

For the general iteration (n > 1), the update from Equation (11) is:

$A_n = A_{n-1} - \frac{A_{n-1} s_n s_n^\top A_{n-1}}{1 + s_n^\top A_{n-1} s_n}$

Factor out $R_{n-1}$ from both terms, using $A_{n-1} = R_{n-1} R_{n-1}^\top$:

$A_n = R_{n-1} \left(I_d - \frac{R_{n-1}^\top s_n (R_{n-1}^\top s_n)^\top}{1 + (R_{n-1}^\top s_n)^\top (R_{n-1}^\top s_n)}\right) R_{n-1}^\top$

Define $\phi_n = R_{n-1}^\top s_n$. Then:

$A_n = R_{n-1} \left(I_d - \frac{\phi_n \phi_n^\top}{1 + \phi_n^\top \phi_n}\right) R_{n-1}^\top$

Now the inner matrix is of the form $I_d - z z^\top$ with $z = \frac{\phi_n}{\sqrt{1 + \phi_n^\top \phi_n}}$. Checking: $z^\top z = \frac{\phi_n^\top \phi_n}{1 + \phi_n^\top \phi_n} < 1$, so Lemma 1 applies to the inner matrix. Its square root is $I_d - r_n \frac{\phi_n \phi_n^\top}{1 + \phi_n^\top \phi_n}$ where $r_n = \frac{1}{1 + \sqrt{\frac{1}{1 + \phi_n^\top \phi_n}}}$.

Therefore, a square root of $A_n$ is:

$R_n = R_{n-1} \left(I_d - r_n \frac{\phi_n \phi_n^\top}{1 + \phi_n^\top \phi_n}\right)$

Expanding the parentheses:

$R_n = R_{n-1} - r_n \frac{(R_{n-1} \phi_n) \phi_n^\top}{1 + \phi_n^\top \phi_n}$

This is exactly the recursion in Equation (13).

Operational summary of one square-root update:

1. Compute $\phi_n = R_{n-1}^\top s_n$ — this is a matrix-vector multiply ($O(d^2)$) mapping the new gradient through the transpose of the current square root.
2. Compute the scalar $r_n = \frac{1}{1 + \sqrt{\frac{1}{1 + \phi_n^\top \phi_n}}}$ — this requires computing $\phi_n^\top \phi_n$ ($O(d)$) and a few scalar operations.
3. Compute the intermediate vector $t_n = R_{n-1} \phi_n$ — another matrix-vector multiply ($O(d^2)$).
4. Form the scaling factor $\frac{r_n}{1 + \phi_n^\top \phi_n}$ and the scaled vector $\tilde{\phi}_n = \frac{r_n}{1 + \phi_n^\top \phi_n} \phi_n$ ($O(d)$).
5. Update $R_n = R_{n-1} - t_n \tilde{\phi}_n^\top$ — a rank-one outer product update ($O(d^2)$ for the outer product, $O(d^2)$ for the addition).

Total cost: $O(d^2)$, dominated by two matrix-vector products and one outer product.

Why the square-root form is necessary (not just a convenience): sampling from the proposal (Equation 7) requires computing $R \eta$ for $\eta \sim \mathcal{N}(0, I_d)$, where $R$ satisfies $R R^\top = A$. If we maintained $A_n$ directly but needed $\sqrt{A_n}$ for sampling, we would need to compute a matrix square root or Cholesky decomposition at every iteration, costing $O(d^3)$. The recursive square-root maintenance avoids this entirely — $R_n$ is available immediately at each step. Additionally, the proposal's gradient term $R_n(R_n^\top \nabla \log \pi(x_n))$ also uses $R_n$ and $R_n^\top$, which are both directly available.

Generalizing to arbitrary learning rates (Appendix B): The standard recursion corresponds to learning rate $\gamma_n = 1/n$ (giving equal weight to all observations). The paper notes in Section 4 that one could use a stochastic approximation framework with a general sequence of learning rates $\gamma_n$:

$\hat{I}_n = \hat{I}_{n-1} + \gamma_n(s_n s_n^\top - \hat{I}_{n-1})$

where $\gamma_n$ satisfies $\sum_n \gamma_n = \infty$ and $\sum_n \gamma_n^2 < \infty$ (standard Robbins-Monro conditions). This leads to a modified inverse recursion:

$A_n = \frac{1}{1 - \gamma_n} \left(A_{n-1} - \frac{A_{n-1} s_n s_n^\top A_{n-1}}{\frac{1-\gamma_n}{\gamma_n} + s_n^\top A_{n-1} s_n}\right)$

Appendix B provides the corresponding square-root recursion for this case. However, the paper reports that in practice, they "did not observe significant improvement by using more general learning rate sequences" and therefore use $\gamma_n = 1/n$ throughout. This is noteworthy because many adaptive MCMC algorithms (including Haario et al., 2001) use $\gamma_n = 1/n$ as well — the Robbins-Monro flexibility is available but empirically unnecessary for these targets.

---

The Adaptation Signal: Why Score Function Increments Instead of Raw Scores

The recursive estimator above takes a "gradient vector" $s_n$ as input at each iteration. But which gradient vector should we use? The paper considers and ultimately rejects two natural candidates before settling on the score function increment.

Candidate 1: Raw score vectors $s_n = \nabla \log \pi(x_n)$.

The most obvious choice is to use the gradient evaluated at the current state $x_n$ at each iteration. Under stationarity ($x_n \sim \pi$), $E[s_n] = 0$ and $E[s_n s_n^\top] = I$, so the empirical Fisher estimator is consistent.

Problem: during the transient phase (before the chain has converged to $\pi$), $x_n$ is not distributed according to $\pi$, and $E[s_n] \neq 0$. The outer product $s_n s_n^\top$ then has expectation:

$E[s_n s_n^\top] = \text{Cov}(s_n) + E[s_n] E[s_n]^\top \neq I$

The additive bias term $E[s_n] E[s_n]^\top$ can be large — if the chain is far from the high-probability region, gradients can have large magnitudes pointing toward the mode. This biases the Fisher estimate toward directions that are simply "pointing toward the mode" rather than directions of genuine posterior curvature. The paper reports (Table 8, Appendix E.3) that FisherMALA using raw score vectors ("FisherMALA non-centered") performs significantly worse: on the GP target, min ESS drops from ~1785 to ~48; on Caravan, it drops from ~498 to ~3.3; on MNIST, from ~440 to ~3.7. This is catastrophic degradation, confirming that the bias during transience is a real and severe problem.

Candidate 2: Centered score vectors via paired mean-covariance estimation (Appendix C).

A standard fix from the adaptive MCMC literature (Haario et al., 2001) is to simultaneously estimate the mean $\bar{s}_n$ of the score vectors and use centered outer products $(s_i - \bar{s}_n)(s_i - \bar{s}_n)^\top$. Appendix C derives the corresponding square-root recursion for this "paired estimation" approach.

Problem: while this reduces bias from non-zero mean, the paper finds experimentally (Table 7, Appendix E.3) that the paired estimation variant ("FisherMALA-paired-est") performs much worse than the score increment approach on most problems. On the inhomogeneous Gaussian target, min ESS drops from ~1500 to ~109; on Caravan, from ~498 to ~87; on MNIST, from ~440 to ~3.2. Only on the low-dimensional Ripley dataset ($d = 3$) does it perform comparably. The paper does not provide a theoretical explanation for this degradation, but a plausible hypothesis is that the centered outer products require accurate estimation of $\bar{s}_n$, which is itself a high-dimensional quantity that converges slowly — the noise in $\bar{s}_n$ propagates into the Fisher estimate, and in high dimensions this noise can dominate the signal.

Candidate 3 (the chosen one): Score function increments $s^\delta_n = s(x_{n+1}) - s(x_n)$.

The key insight is that we can use the difference between consecutive gradient evaluations. At each iteration, the MCMC chain moves from $x_n$ to $x_{n+1}$ (where $x_{n+1}$ is either the proposed $y_n$ if accepted, or $x_n$ if rejected). The score increment is:

$s^\delta_n = s(x_{n+1}) - s(x_n) = \mathbb{I}(u_n < \alpha(x_n, y_n))(s(y_n) - s(x_n))$

where $u_n \sim \mathcal{U}(0, 1)$ is the uniform random number used in the accept/reject step, and $\mathbb{I}(\cdot)$ is the indicator function.

What this computes: if the proposal is accepted, $s^\delta_n$ is the difference between the gradient at the new state and the gradient at the old state. If the proposal is rejected, $s^\delta_n = 0$ (since $x_{n+1} = x_n$, so $s(x_{n+1}) = s(x_n)$).

Why this works better during transience: even when the chain is far from stationarity, gradients at consecutive states tend to point in similar directions — they both point roughly toward the mode, and their difference is primarily influenced by the local curvature of $\log \pi$ between $x_n$ and $y_n$. The non-zero mean cancels out to first order. More formally, $E[s^\delta_n] \approx 0$ even when $E[s_n] \neq 0$, because $s(x_{n+1})$ and $s(x_n)$ have similar biases that subtract away. This makes the outer product $s^\delta_n (s^\delta_n)^\top$ a much better estimator of the Fisher matrix during the burn-in phase.

Why this works better than centered estimation: the difference operation provides "centering" automatically and locally — it does not require a running mean estimate that may be noisy in high dimensions. Each pair of consecutive states provides a local curvature measurement that is largely unbiased by the chain's position relative to the mode.

The Rao-Blackwellization improvement (Equation 16):

The indicator $\mathbb{I}(u_n < \alpha)$ introduces needless variance because it randomly zeros out the signal when a proposal is rejected. Since $u_n$ is an auxiliary variable independent of the states, we can Rao-Blackwellize (marginalize it out) to reduce variance:

$E_{u_n}[s^\delta_n (s^\delta_n)^\top] = \alpha(x_n, y_n)(s(y_n) - s(x_n))(s(y_n) - s(x_n))^\top$

This suggests using the Rao-Blackwellized signal:

$s^\delta_n = \sqrt{\alpha(x_n, y_n)}(s(y_n) - s(x_n))$

where $\alpha(x_n, y_n)$ is the Metropolis-Hastings acceptance probability.

What this computes: instead of zeroing out the gradient difference on rejection, this weights the difference by the square root of the acceptance probability. When $\alpha$ is small (the proposal would likely be rejected), the signal is downweighted but not eliminated. When $\alpha$ is large, the signal is fully used.

Why the square root: the outer product $s^\delta_n (s^\delta_n)^\top = \alpha \cdot (s(y_n) - s(x_n))(s(y_n) - s(x_n))^\top$ recovers exactly the Rao-Blackwellized expectation. Using $\alpha$ directly (without the square root) in the recursive update would correspond to using the squared signal, which would be inconsistent with the Woodbury derivation that assumes $s_n$ enters as an outer product $s_n s_n^\top$.

Practical advantage: even when a proposal is rejected (which happens frequently — target acceptance rate is ~0.574), the Rao-Blackwellized signal is non-zero as long as $\alpha > 0$. This means the algorithm learns about the target's geometry from every single proposal, not just accepted ones. This is a form of data efficiency that purely state-based adaptation (which only sees accepted moves) lacks.

Experimental validation (Appendix E.3, Figure 15, Table 7): The paper compares FisherMALA with and without Rao-Blackwellization. Figure 15 shows that Rao-Blackwellization leads to slightly faster convergence of the Frobenius norm on the two Gaussian targets, which "agrees with the theory that says that Rao-Blackwellization should reduce the variance." Table 7 shows that the two variants (with and without Rao-Blackwellization) have very similar ESS scores across all problems — the Rao-Blackwellized version is marginally better in some cases, but the difference is not statistically significant. This suggests that the critical innovation is using score differences (rather than raw scores or centered scores), and that Rao-Blackwellization provides a modest additional variance reduction.

Why the score increment approach does NOT have an obvious unbiasedness guarantee: Section 7 flags this as an open question. The outer products $s^\delta_n (s^\delta_n)^\top$ use correlated gradient evaluations from consecutive MCMC states, which are not independent draws from $\pi$. It is not obvious that their expectation equals the Fisher matrix $I$, even asymptotically. The empirical evidence strongly suggests that any bias is small and does not prevent effective preconditioning, but a theoretical analysis of the bias (or a proof of consistency) is left as future work.

---

The Complete Algorithm: FisherMALA (Algorithm 1)

Algorithm 1 integrates all the components described above into a complete adaptive MCMC procedure. Here is the detailed walkthrough:

Inputs:

• log π(x) and ∇ log π(x): the target log-density and its gradient (assumed computable for any $x$).
• λ > 0: the damping parameter for the Fisher estimate. The paper sets a default of $\lambda = 10$ and uses this value in all experiments. This is a hyperparameter that controls the initial regularization strength — larger $\lambda$ means more initial shrinkage toward isotropic preconditioning.
• α* = 0.574: the target acceptance rate, motivated by optimal scaling results for MALA (Roberts and Rosenthal, 1998).

Initialization phase (first $n_0$ iterations, default $n_0 = 500$): The algorithm starts by running simple (non-preconditioned) MALA — i.e., with proposal $\mathcal{N}(y | x + (\sigma^2/2) \nabla \log \pi(x), \sigma^2 I_d)$ — for $n_0$ iterations. During this phase:

• The step size $\sigma^2$ is adapted toward the target acceptance rate $\alpha^*$ using the Robbins-Monro update: $\sigma^2 \leftarrow \sigma^2 [1 + \rho_n(\alpha(x_n, y_n) - \alpha^*)]$, where $\rho_n$ is a learning rate (the paper uses a constant $\rho_n = 0.015$ throughout).
• The square root matrix $R$ remains fixed at $R = I_d$ (identity) — no preconditioning is learned during this phase.

Purpose of the initialization phase: this allows the chain to move away from its starting point and begin approaching the typical set of $\pi$ before preconditioner adaptation begins. Without this, the early gradient evaluations (which may be from an arbitrary initialization far from the mode) could bias the Fisher estimate. The paper found that 500 iterations provides enough burn-in for this purpose. The initial value of $\sigma^2$ is also set during this phase based on early acceptance rates.

Main loop (iterations n = 1, 2, 3, ... after initialization):

Step 1: Propose. Generate a candidate state $y_n$ from the preconditioned proposal:

$y_n = x_n + \frac{\sigma^2_R}{2} R(R^\top \nabla \log \pi(x_n)) + \sigma_R R \eta, \quad \eta \sim \mathcal{N}(0, I_d)$

where $\sigma^2_R = \sigma^2 / \frac{1}{d}\text{tr}(R R^\top)$ is the normalized step size, and $R$ is the current square root estimate. The gradient term $R(R^\top \nabla \log \pi(x_n))$ first projects the gradient through $R^\top$, then through $R$, effectively applying $A = R R^\top$ as the preconditioner for the drift. The noise term $\sigma_R R \eta$ generates correlated noise with covariance $\sigma^2_R A$.

Step 2: Compute log-density and gradient at the proposal. Evaluate $\log \pi(y_n)$ and $\nabla \log \pi(y_n)$. These are the expensive operations (in typical applications, computing the target density and its gradient dominates the computational cost, especially if it involves passing through a large dataset).

Step 3: Compute the acceptance probability. Use Proposition 1's simplified formula to avoid matrix inverses:

$\alpha(x_n, y_n) = \min\left(1, \exp\{\log \pi(y_n) + h(x_n, y_n) - \log \pi(x_n) - h(y_n, x_n)\}\right)$

where $h(z, v) = \frac{1}{2}\left(z - v - \frac{\sigma^2_R}{4} R R^\top \nabla \log \pi(v)\right)^\top \nabla \log \pi(v)$. Note that Proposition 1 is stated for a generic $A$ in the proposal; here $A = R R^\top$, and $A \nabla \log \pi(v)$ is computed as $R(R^\top \nabla \log \pi(v))$, which avoids forming $A$ explicitly.

Why Proposition 1 is important operationally: the naive Metropolis-Hastings ratio involves the term $\frac{q(x_n | y_n)}{q(y_n | x_n)}$, which for the Gaussian proposal involves $A^{-1}$ (since the precision matrix appears in the exponent). $A^{-1}$ would be expensive to compute and maintain. Proposition 1 shows that the log-ratio of proposals simplifies to $h(x_n, y_n) - h(y_n, x_n)$, which only involves $A$ (not $A^{-1}$). This is proven in Appendix A.1: the terms involving $A^{-1}$ in the proposal log-densities cancel exactly, leaving only terms involving $A \nabla \log \pi(\cdot)$, which can be computed using $R R^\top \nabla \log \pi(\cdot)$ without any matrix inversion.

Step 4: Compute the adaptation signal. Using the Rao-Blackwellized form (Equation 16):

$s^\delta_n = \sqrt{\alpha(x_n, y_n)}(\nabla \log \pi(y_n) - \nabla \log \pi(x_n))$

If $\alpha(x_n, y_n) = 0$ (should not happen in practice with proper tuning, but possible numerically), $s^\delta_n = 0$ and no adaptation occurs.

Step 5: Update the square root matrix (Proposition 4).

If $n = 1$ (first iteration after initialization), use the initialization formula (Equation 12):

$R_1 = \frac{1}{\sqrt{\lambda}} \left(I_d - r_1 \frac{s^\delta_1 (s^\delta_1)^\top}{\lambda + (s^\delta_1)^\top s^\delta_1}\right), \quad r_1 = \frac{1}{1 + \sqrt{\frac{\lambda}{\lambda + (s^\delta_1)^\top s^\delta_1}}}$

If $n > 1$, use the recursive update (Equation 13):

$\phi_n = R_{n-1}^\top s^\delta_n, \quad r_n = \frac{1}{1 + \sqrt{\frac{1}{1 + \phi_n^\top \phi_n}}}, \quad R_n = R_{n-1} - r_n \frac{(R_{n-1} \phi_n) \phi_n^\top}{1 + \phi_n^\top \phi_n}$

Note that $s_n$ in the original Proposition 4 is replaced by $s^\delta_n$ in the algorithm — the adaptation operates on score increments, not raw scores.

Step 6: Adapt the step size. Update the global step size using the Robbins-Monro scheme:

$\sigma^2 \leftarrow \sigma^2 \left[1 + \rho_n (\alpha(x_n, y_n) - \alpha^*)\right]$

where $\rho_n$ is a learning rate. The paper uses a constant $\rho_n = 0.015$. If the acceptance rate is above the target $\alpha^* = 0.574$, $\sigma^2$ is increased (allowing larger steps); if below, $\sigma^2$ is decreased (making steps more conservative). This is standard practice in adaptive MCMC and is independent of the preconditioner learning.

Step 7: Normalize the step size. Compute the normalized step size for the next iteration:

$\sigma^2_R = \frac{\sigma^2}{\frac{1}{d} \text{tr}(R_n R_n^\top)}$

Since $\text{tr}(R_n R_n^\top) = \sum_{i,j} (R_n)_{ij}^2 = \text{sum}(R_n \circ R_n)$ (elementwise square and sum), this is an $O(d^2)$ operation.

Step 8: Accept or reject. Draw $u \sim \mathcal{U}(0, 1)$. If $u < \alpha(x_n, y_n)$, set $x_{n+1} = y_n$, $\log \pi(x_{n+1}) = \log \pi(y_n)$, $\nabla \log \pi(x_{n+1}) = \nabla \log \pi(y_n)$. Otherwise, set $x_{n+1} = x_n$, and the stored log-density and gradient remain unchanged.

Post-adaptation phase (sample collection): After a specified number of burn-in iterations (20,000 in all experiments), adaptation of both $\sigma^2$ and $R$ is frozen — the learned values are used as fixed parameters for the remainder of the run. This is standard practice in adaptive MCMC to ensure the chain is a valid Markov chain during sample collection (satisfying the "diminishing adaptation" condition required for ergodicity).

Design choice: $\lambda = 10$. The damping parameter controls the strength of regularization in the Fisher estimate. At $\lambda = 10$, the identity component in $A_1 = \frac{1}{\lambda}(I_d - \frac{s_1 s_1^\top}{\lambda + s_1^\top s_1})$ is initially substantial, which keeps early preconditioning conservative. This prevents the preconditioner from over-committing to directions seen in early iterations before the chain has explored the target. As $n$ grows, the relative influence of $\lambda$ decays (since $\lambda/n \to 0$ in the batch estimator), allowing the data to dominate. The paper uses $\lambda = 10$ throughout and does not report sensitivity analysis for this choice — it is likely that performance is robust within a reasonable range (e.g., 1–100), as the regularization is primarily important early in adaptation and is quickly overwhelmed by data.

Design choice: constant learning rate $\rho_n = 0.015$ for step size adaptation. A constant (rather than decaying) learning rate is used because step size adaptation should remain responsive even late in the burn-in phase, as the effective scale of the target may change as the chain explores new regions. The value $0.015$ is standard in the MALA literature and corresponds to a relatively slow, stable adaptation.

Design choice: initializing AdaMALA vs. FisherMALA. The paper notes (Appendix D) that FisherMALA's initialization (500 iterations of simple MALA) proved unstable for AdaMALA. AdaMALA required an additional 500-iteration warm-up phase (1,000 total initialization iterations) where simple MALA collects states used to initialize the covariance matrix before AdaMALA's own adaptation begins. This asymmetry in initialization favors FisherMALA in the comparison, but the paper's Frobenius norm trajectories (Figures 1a-b) show that FisherMALA's advantage persists even when both methods are given comparable total adaptation time — FisherMALA converges in ~2,000 iterations while AdaMALA requires the full 20,000 iterations and still achieves worse accuracy.

---

Summary of Design Choices and Their Justifications

• Inverse Fisher over covariance: derived from first-principles optimization (Proposition 3) rather than inherited from RWM heuristics. Empirically validated by FisherMALA's dramatic superiority over AdaMALA on non-Gaussian targets.
• Score function increments over raw scores or centered scores: avoids bias from non-zero $E[s_n]$ during transience; the Rao-Blackwellized form recovers information from rejected proposals; experimentally, alternatives (non-centered and paired-est) severely degrade.
• Square-root recursion over explicit inverse + Cholesky: reduces per-iteration cost from $O(d^3)$ to $O(d^2)$; adapted from Kalman filtering literature; provides $R_n$ directly for sampling.
• $1/n$ learning rate over general Robbins-Monro: matches standard adaptive MCMC practice; empirically no benefit from other sequences.
• $\lambda = 10$ damping: Tikhonov regularization for early-iteration stability; decays with $n$ so consistency is asymptotically unaffected.
• Normalized step size parametrization (Equation 7): decouples learning $A$'s shape from learning $\sigma^2$'s scale via scale invariance; essential for stable joint adaptation.
• Proposition 1's simplified acceptance ratio: avoids $A^{-1}$ computation; makes the algorithm truly $O(d^2)$ by eliminating matrix inversion from the accept/reject step.
• 500-iteration simple MALA initialization: allows the chain to approach the typical set before Fisher estimation begins, reducing early bias.

4. Key Insights and Innovations

Innovation 1: A First-Principles Derivation That the Optimal Preconditioner Is the Inverse Fisher, Not the Target Covariance

Prior to this work, the adaptive MCMC literature operated under a widely-held but largely unexamined assumption: that the optimal preconditioning matrix for MALA is the target covariance $\Sigma$, inherited directly from analogous results for random walk Metropolis (RWM). This heuristic was codified in the foundational adaptive MCMC algorithm of Haario et al. (2001), which recursively estimates $\Sigma$ from the chain's state history and uses it as the proposal covariance. The reasoning was intuitive — if the target is stretched more in some directions than others, the proposal should stretch correspondingly — but it was never derived from first principles for gradient-based samplers. The paper overturns this assumption through an explicit optimization that shows the optimal preconditioner is actually $I^{-1} = E_\pi[\nabla \log \pi \nabla \log \pi^\top]^{-1}$, the inverse Fisher information matrix, which coincides with $\Sigma$ only for Gaussian targets.

What makes this a genuinely fundamental rather than incremental contribution is that the derivation starts from a clean, interpretable objective — maximize the allowable discretization step size subject to a bound on the expected squared jump distance — and solves it analytically without asymptotic arguments or Gaussian approximations. The expected squared jumped distance criterion itself was previously used only for tuning scalar step sizes (Pasarica and Gelman, 2010); this paper elevates it to a tool for optimizing the full $d \times d$ preconditioning matrix. The optimization reveals that the relevant geometric quantity is not "how spread out are the states under $\pi$" (the covariance), but rather "how much does the log-density curve in each direction averaged over $\pi$" (the Fisher information). This is a conceptual reframing: the preconditioner should respond to the target's curvature structure (encoded in gradients), not its spread structure (encoded in state covariances).

The significance of this theoretical correction extends beyond the algorithm itself. It explains why covariance-based preconditioning works reasonably well for near-Gaussian posteriors (where $\Sigma = I^{-1}$ under the Gaussian special case shown in Section 3) but degrades for non-Gaussian targets — a pattern the paper demonstrates empirically in Table 1, where AdaMALA's minimum ESS collapses relative to FisherMALA's on the non-Gaussian Bayesian logistic regression problems (Caravan: ~3.1 vs. ~498; MNIST: ~3.0 vs. ~440). It also clarifies the connection to the optimal scaling literature: Roberts and Rosenthal (2001) identified $I^{-1}$ as optimal for RWM on heterogeneous product-form targets as $d \to \infty$, but that result was asymptotic and restricted to diagonal $I$. Proposition 3 generalizes this to any full-rank $I$ and any dimension, without asymptotic assumptions. This positions the inverse Fisher as a unifying object across gradient-free and gradient-based MCMC, while correcting the specific form it should take in the gradient-based case.

Innovation 2: Gradient-History Adaptation as a Fundamentally Different (and Superior) Learning Paradigm for Preconditioning

The dominant paradigm in adaptive MCMC since Haario et al. (2001) has been state-based adaptation: learn the target's geometry from the sequence of visited states $x_1, x_2, \ldots$. FisherMALA introduces a fundamentally different paradigm: gradient-based adaptation, where the target's geometry is learned from the sequence of score function evaluations $s_1, s_2, \ldots$ and specifically from their differences $s_{n+1} - s_n$. This is not a small modification to existing adaptive MCMC — it is a shift in what information source the algorithm exploits, with profound implications for learning speed and stability.

The conceptual distinction hinges on information density. A state $x_n$ is a single point in a $d$-dimensional space. Learning a $d \times d$ covariance matrix from states requires observing the chain traverse enough of the target's support to estimate pairwise coordinate relationships — a process that is slow precisely when the chain mixes slowly (the very problem preconditioning is meant to solve). This creates the vicious cycle described in Section 2: poor mixing begets slow adaptation, which perpetuates poor mixing. A gradient $s_n = \nabla \log \pi(x_n)$, by contrast, is a local description of the target's geometry at $x_n$ — it carries information about the direction and magnitude of curvature regardless of whether the chain has fully explored the target. The Fisher matrix aggregates these local curvature measurements, and as the empirical results show (Figures 1b, 2c), this aggregation converges to the correct global geometry orders of magnitude faster than state covariance estimation.

The paper's key diagnostic move is recognizing that gradient differences $s^\delta_n = s(x_{n+1}) - s(x_n)$ are an even better signal than raw gradients. Raw gradients have non-zero expectation during the transient phase (when $x_n$ is far from the mode, gradients all point systematically toward it), which biases the outer product estimate through the term $E[s_n]E[s_n]^\top$. Gradient differences cancel this bias automatically — since consecutive gradients point in similar directions during transience, their difference primarily reflects local curvature rather than location relative to the mode. Table 8 in Appendix E.3 provides stark evidence: FisherMALA using raw gradients ("non-centered") achieves min ESS of ~48 on the GP target (vs. ~1785 for the gradient-difference version) and ~3.7 on MNIST (vs. ~440). This is not a modest degradation — it is catastrophic failure, confirming that how gradient information is used matters as much as that it is used.

The Rao-Blackwellization step ($s^\delta_n = \sqrt{\alpha}(s(y_n) - s(x_n))$) adds a further conceptual refinement: it ensures that every proposal contributes to learning, not just accepted ones. In standard MCMC, information from rejected proposals is discarded — the chain stays put, and the adaptation (if state-based) sees nothing new. The Rao-Blackwellized gradient difference, by weighting by $\sqrt{\alpha}$, extracts curvature information from the proposed state $y_n$ even when it is rejected, since the gradient at $y_n$ still reveals how $\log \pi$ changes between $x_n$ and $y_n$. This is a form of data efficiency with no analog in state-based adaptation, where rejected moves leave the state unchanged and provide zero new information.

Innovation 3: The Square-Root Recursion as an Enabling Computational Technique That Makes Online Inverse Fisher Learning Practically Viable

Learning $I^{-1}$ online during MCMC poses a computational challenge: at each iteration, we need to update an estimate of a $d \times d$ inverse matrix and produce a square root for sampling, all without $O(d^3)$ matrix decompositions that would dominate the cost. The paper's solution — adapting Potter's square-root filtering method from the 1960s Kalman filtering literature (Potter and Stern, 1963; Bierman, 1977) — is an elegant example of cross-domain technology transfer that transforms a theoretical insight (the optimality of $I^{-1}$) into a practical algorithm.

What makes this innovation distinctive is not the recursion itself (which is mathematically standard) but the recognition that the Woodbury rank-one update structure of the empirical Fisher inverse makes the square-root form exactly maintainable in $O(d^2)$. The critical observation is that each new gradient outer product $s_n s_n^\top$ modifies $A_n$ by a rank-one downdate (Equation 11), and Lemma 1 shows that the square root of a rank-one perturbation of the identity is another rank-one perturbation of the identity, with a modified coefficient. This means the square root recursion mirrors the structure of the inverse recursion exactly, with no additional approximation or decomposition required.

The significance extends beyond the immediate algorithm. The bottleneck in many adaptive MCMC and optimization methods that use matrix estimates (covariance matrices, Hessian approximations, Fisher matrices) is the $O(d^3)$ cost of maintaining a decomposition for sampling or search directions. The square-root filtering approach demonstrates that for estimates built from rank-one outer product updates, the decomposition can be maintained recursively in $O(d^2)$, eliminating the cubic bottleneck entirely. This principle could apply to other methods that currently pay $O(d^3)$ per iteration for matrix decompositions — quasi-Newton samplers (Zhang and Sutton, 2011), stochastic Newton MCMC (Martin et al., 2012), or any Langevin-based method requiring a matrix square root. The paper's Algorithm 1 with Proposition 4 provides a template for how to structure such recursive square-root updates that future methods could adopt.

It is worth noting that this innovation is primarily enabling rather than conceptual — it does not change what matrix is being estimated or why, but it makes the estimation computationally tractable at scale. Without it, FisherMALA would be limited to low-dimensional problems where $O(d^3)$ per-iteration Cholesky decompositions are acceptable. With it, the method scales to the 785-dimensional MNIST example (Table 1) with per-iteration cost comparable to simple MALA plus two matrix-vector products. This transforms FisherMALA from a theoretically interesting but computationally impractical idea into a method that can be deployed on real high-dimensional Bayesian inference problems.

Innovation 4: Empirical Evidence That a Global Inverse Fisher Preconditioner Can Outperform a Position-Dependent Exact Riemannian Preconditioner

Perhaps the most counterintuitive result in the paper is that FisherMALA — which learns a single global preconditioner averaged over the entire target — significantly outperforms Riemannian manifold MALA (mMALA, Girolami and Calderhead, 2011), which uses the exact position-dependent negative inverse Hessian $A(x) = (-\nabla^2 \log \pi(x))^{-1}$ as the preconditioner at every step. On the Caravan dataset ($d = 87$), mMALA achieves min ESS of ~51 vs. FisherMALA's ~498 (a ~10× improvement). On MNIST ($d = 785$), the gap widens to ~5.2 vs. ~440 (an ~85× improvement). These are not marginal differences — they represent a qualitative reversal of what one would expect if local geometric information were uniformly beneficial.

The paper does not provide a definitive explanation for this phenomenon, but hypothesizes (Section 6) that "position-dependent preconditioning may be less effective in certain types of high-dimensional and log-concave problems." This result has important implications for how we think about geometric information in MCMC: more information is not always better, and global averaging can be a form of beneficial regularization. The local Hessian $-\nabla^2 \log \pi(x)$ varies across the state space — it may be ill-conditioned or near-singular in some regions, and the noise in its estimation (especially in high dimensions where finite-sample Hessian estimates are unreliable) may cause the local preconditioner to propose poor moves. The global inverse Fisher $I^{-1}$, by averaging over the stationary distribution, smooths out these local fluctuations and may provide more robust curvature information that works well across the entire typical set.

This finding also reframes the cost-benefit analysis of local vs. global preconditioning. mMALA's $O(d^3)$ per-iteration cost (for Hessian computation and matrix decomposition) means it takes far fewer MCMC iterations within a fixed computational budget. FisherMALA's $O(d^2)$ cost enables many more iterations, and the combination of "good enough" global preconditioning with higher iteration counts appears to dominate "optimal" local preconditioning with lower iteration counts in these problems. This is a pragmatic insight: in high dimensions, computational efficiency of the adaptation mechanism can matter more than the theoretical optimality of the preconditioner being adapted toward. This suggests a broader principle for adaptive MCMC design: the learning mechanism's cost must be factored into the effective sampling efficiency, not just the quality of the learned parameters.

Innovation 5: The Difficulty-Conditioned, Difficulty-Agnostic Nature of the Adaptation Signal

A subtle but important innovation is what the algorithm does NOT require: it does not need any problem-specific tuning, any second-order derivatives, or any structural assumptions about the target beyond differentiability of $\log \pi$. The adaptation signal $s^\delta_n = \sqrt{\alpha}(s(y_n) - s(x_n))$ uses only quantities that the MCMC algorithm already computes — the gradients at the current and proposed states, and the acceptance probability. This is in contrast to:

• Standard adaptive MCMC (Haario et al., 2001), which requires maintaining and updating running means and covariances of the state vectors, with careful initialization and damping to avoid degeneracy.
• Riemannian methods (Girolami and Calderhead, 2011), which require computing or approximating the Hessian matrix and its decomposition, with model-specific derivations for each new target.
• Quasi-Newton samplers (Zhang and Sutton, 2011), which require maintaining a history of states and gradients to build Hessian approximations.
• Neural adaptive MCMC (Levy et al., 2018; Song et al., 2017), which require training auxiliary networks on the target distribution.

FisherMALA's design embodies a principle of maximal information extraction from minimal computational overhead: every quantity used for adaptation is already computed for the MCMC accept/reject step, and the adaptation adds only $O(d^2)$ additional work. The algorithm inherits the gradient evaluations it needs from the MALA proposal; it does not require separate gradient computations, separate proposals, or separate objective function optimizations running in parallel with the chain.

This "free lunch" property is significant for practical adoption. Methods that require additional infrastructure (neural network training, Hessian code, complex initialization protocols) face adoption barriers because practitioners must implement and debug these components for each new model. FisherMALA requires only $\nabla \log \pi(x)$, which is available from automatic differentiation in modern probabilistic programming frameworks (Stan, PyMC, TensorFlow Probability) with no additional user effort. The damping parameter $\lambda = 10$ is the only tunable hyperparameter beyond the step size, and the paper's experiments suggest the default works across a wide range of problems without adjustment. This robustness to hyperparameter choice is itself a significant practical property, even if it is not theoretically characterized.

5. Experimental Analysis

Evaluation Methodology

• Datasets. The experiments span two broad categories. For Gaussian targets (where the true optimal preconditioner $I^{-1} = \Sigma$ is known analytically), the paper uses three synthetic examples: a two-dimensional Gaussian with covariance $\Sigma = [1, 0.995; 0.995, 1]$, a 100-dimensional Gaussian process (GP) correlated target where $\Sigma$ is constructed from a non-stationary covariance function $[\Sigma]_{ij} = s_i s_j \exp\{-\frac{1}{2}\frac{(s_i - s_j)^2}{0.09}\} + 0.001\delta_{ij}$ with $s_i$ on a regular grid in $[1, 2]$, and a 100-dimensional inhomogeneous Gaussian where $\Sigma = \text{diag}(\sigma_1^2, \ldots, \sigma_{100}^2)$ with $\sigma_i$ taking values in $\{0.01, 0.02, \ldots, 1\}$. For Bayesian logistic regression, the paper uses six binary classification datasets with unstandardized inputs (to create anisotropic posteriors): Australian Credit, Heart, Pima Indian, Ripley, German Credit, and Caravan, with $d$ ranging from 3 to 87 and sample sizes from ~250 to ~5,822. A higher-dimensional 785-parameter example on MNIST classifies digits 5 and 6 using 11,339 training examples, with pixel values scaled to $[0, 1]$ by dividing by 255.

• Base models. The "model" here is the target distribution $\pi(x)$ itself — this is not an LLM evaluation but an MCMC sampler evaluation. The targets are chosen to represent a spectrum of difficulty: the GP target tests performance with dense off-diagonal correlations; the inhomogeneous Gaussian tests extreme eigenvalue disparities (two orders of magnitude in standard deviation, challenging for samplers with a single step size); the Bayesian logistic regression problems test non-Gaussianity with real data, where $I^{-1} \neq \Sigma$. The MNIST example at $d = 785$ pushes into moderate high dimensions where cubic-cost methods become prohibitive.

• Metrics. The primary metric is effective sample size (ESS) computed across each dimension of the state vector independently, with maximum, median, and minimum values reported to capture the best, typical, and worst-case mixing behavior. ESS is computed from $2 \times 10^4$ post-burn-in samples using the built-in method in TensorFlow Probability. For Gaussian targets where the ground-truth optimal preconditioner $\tilde{\Sigma} = \Sigma / (\text{tr}(\Sigma)/d)$ is known, the paper also reports the Frobenius norm $\|\tilde{A}_n - \tilde{\Sigma}\|_F$ across burn-in iterations, where $\tilde{A}_n = A_n / (\text{tr}(A_n)/d)$ is the normalized learned preconditioner — this directly measures how quickly and accurately the adaptive scheme recovers the optimal matrix. Additionally, trace plots of the unnormalized log-target $\log p(Y | \theta, Z) p(\theta)$ are shown to visualize chain autocorrelation qualitatively. Results are averaged over 10 random initializations with standard deviations reported.

• Baselines. The paper compares against four samplers: (1) Simple MALA: proposal $\mathcal{N}(y_n | x_n + (\sigma^2/2) \nabla \log \pi(x_n), \sigma^2 I_d)$, adapting only the step size $\sigma^2$ toward acceptance rate 0.574. (2) AdaMALA (Haario et al., 2001): preconditioned MALA with the proposal in Equation (7) but where the preconditioning matrix is learned using the standard state-based adaptive MCMC recursion (Equation 17) — running mean $\mu_n$ and covariance $\Sigma_n$ updated from the chain's state history, with the same damping $\lambda = 10$ and damping initialization $\Sigma_2 = \frac{1}{2}(x_2 - \mu_1)(x_2 - \mu_1)^\top + \lambda I_d$. (3) Riemannian manifold MALA (mMALA) (Girolami and Calderhead, 2011): uses position-dependent preconditioning $A(x)$ based on the negative inverse Hessian $(-\nabla^2 \log \pi(x))^{-1}$, costing $O(d^3)$ per iteration. For the Gaussian targets, mMALA is the theoretically optimal MALA sampler since it uses the exact $\Sigma$ as preconditioning. (4) Hamiltonian Monte Carlo (HMC) (Neal, 2010): 10 leapfrog steps with identity mass matrix, step size adapted toward 0.651 acceptance rate (optimal for HMC; Beskos et al., 2013). FisherMALA itself follows Algorithm 1 with $\lambda = 10$, score increment adaptation using the Rao-Blackwellized signal (Equation 16), and standard learning rate $\gamma_n = 1/n$.

• Generation budget / compute accounting. All samplers are run for $2 \times 10^4$ burn-in iterations followed by $2 \times 10^4$ collection iterations. Adaptation of proposal parameters ($\sigma^2$, preconditioning matrices, HMC step size) occurs only during burn-in; parameters are frozen during sample collection to satisfy the diminishing adaptation condition for valid MCMC. No explicit FLOP counting is performed — the metric is ESS per fixed number of MCMC iterations, implicitly measuring sampling efficiency per unit of computational effort (since all MALA variants have comparable per-iteration cost except mMALA which is $O(d^3)$). FisherMALA uses the first 500 burn-in iterations for initialization with simple MALA (no preconditioning adaptation); the remaining $1.95 \times 10^4$ burn-in iterations adapt the preconditioner. AdaMALA uses a more elaborate initialization (Appendix D): 500 iterations of simple MALA plus an additional 500 iterations of simple MALA collecting states to warm-start the covariance estimate, totaling 1,000 initialization iterations before preconditioner adaptation begins. FisherMALA uses a constant learning rate $\rho_n = 0.015$ for step-size adaptation (Robbins-Monro); AdaMALA, mMALA, and HMC use analogous schemes to target their respective optimal acceptance rates.

• Cross-validation / statistical protocol. Each experiment is repeated 10 times with different random initializations. All tables report mean ± standard deviation of ESS scores across these 10 repetitions. For Gaussian targets, the Frobenius norm trajectories (Figures 1a-b, 2c) are computed on a single representative run to show convergence behavior. The MNIST and Caravan problems are not reported with cross-validation folds on the data — the training set is fixed, and variability comes only from MCMC random seeds. This means the ESS comparisons are conditional on a single dataset realization, which is standard in MCMC benchmarking but limits generalizability claims to new data draws.

Main Quantitative Results

Gaussian Targets: FisherMALA Converges to the Optimal Preconditioner Orders of Magnitude Faster Than AdaMALA

Two-dimensional Gaussian (Figure 1a, Figure 5). On the $2 \times 2$ correlated Gaussian, both FisherMALA and AdaMALA can learn the correct preconditioner (the target covariance $\Sigma$), but FisherMALA converges faster. The Frobenius norm $\|\tilde{A}_n - \tilde{\Sigma}\|_F$ (Figure 1a) shows FisherMALA's error dropping more rapidly in the first few hundred iterations, though both methods eventually reach near-zero error since this is an easy low-dimensional case. Visualizations in Figure 5 (Appendix E.2) confirm that both methods recover the true covariance ellipse accurately.

GP correlated target, 100 dimensions (Figure 1b-d, Figure 6, Table 1). This is where the learning speed difference becomes stark. Figure 1b shows the Frobenius norm trajectory across 20,000 burn-in iterations: FisherMALA's error drops to near zero within approximately 2,000 iterations and stays there, while AdaMALA's error remains large even after 20,000 iterations — it is still visibly decreasing at the end of burn-in and never reaches FisherMALA's accuracy. Figure 1c shows the true $100 \times 100$ GP covariance matrix (displays the non-stationary structure: high correlation near the diagonal decaying with distance), and Figure 1d shows the estimate learned by FisherMALA, which is nearly visually indistinguishable from the truth. Figure 6 (Appendix E.2) shows AdaMALA's estimate — recognizably similar in broad structure but noticeably noisier, consistent with the higher Frobenius error.

The ESS consequences appear in Table 1. FisherMALA achieves max ESS 2096, median ESS 1924, min ESS 1785 — all within ~15% of each other, indicating uniform mixing across dimensions. AdaMALA achieves max 845, median 663, min 552 — substantially lower, and notably the min ESS (~552) is less than one-third of FisherMALA's. Critically, FisherMALA essentially matches mMALA (which uses the exact optimal preconditioner): mMALA achieves max 2109, median 2008, min 1842 vs. FisherMALA's 2096, 1924, 1785. The ~3–8% gap is attributable to FisherMALA still converging toward the optimal matrix and to Monte Carlo variability. This is the first major quantitative result: a learned global preconditioner, adapted purely from gradient history during the burn-in phase, achieves sampling efficiency indistinguishable from one using the analytically known optimal preconditioner — and it does so with $O(d^2)$ per-iteration cost vs. mMALA's $O(d^3)$.

HMC with 10 leapfrog steps and identity mass matrix achieves max ESS ~18, median ~7, min ~4 — worse than simple MALA (max ~15, median ~6, min ~4). This illustrates that HMC's Hamiltonian dynamics provide no advantage without a properly tuned mass matrix: on this highly anisotropic target, the identity mass matrix fails as badly as isotropic MALA.

Inhomogeneous Gaussian target, 100 dimensions (Figure 2, Table 2 in Appendix E.2). This target is specifically designed to be maximally challenging for samplers without per-dimension scaling: $\Sigma = \text{diag}(0.01^2, 0.02^2, \ldots, 1.00^2)$, so the variance spans four orders of magnitude from $10^{-4}$ to $1$. Figure 2a shows per-dimension ESS for all methods across the 100 coordinates. FisherMALA and mMALA (which overlap almost perfectly) maintain uniformly high ESS (~1500–2350) across all dimensions, from the stiffest ($x_1$, variance $10^{-4}$) to the most relaxed ($x_{100}$, variance 1). In contrast:

• Simple MALA achieves high ESS on the largest-variance dimensions (ESS ~13,700 for $x_{100}$) but collapses to min ESS ~2.9 for $x_1$ — the step size adapts to the smallest-variance dimension, making moves in that direction appropriately small but freezing exploration everywhere else.
• HMC shows a similar pattern: ESS ~19,400 on the relaxed dimensions but min ESS ~42 on the stiffest — better than MALA but still unacceptably low for reliable inference on $x_1$.
• AdaMALA is intermediate: max ESS ~4,311, median ~71, min ~9.2. It learns some preconditioning but far too slowly — the stiff dimensions remain poorly mixed because the covariance estimate hasn't converged for those directions within 20,000 iterations.

Table 2 (Appendix E.2) provides the numbers: FisherMALA achieves min ESS 1501 ± 67 vs. AdaMALA's 9.2 ± 3.3 — a ~163× improvement on the worst-mixing dimension. Figure 2b shows the estimated standard deviations (diagonal of the learned preconditioner) for FisherMALA and AdaMALA at the end of burn-in. FisherMALA's estimates closely track the true $\sigma_i$ values across the full range; AdaMALA's estimates show the correct trend but with substantial noise and systematic underestimation for the smallest $\sigma_i$. Figure 2c confirms the learning speed diagnosis: FisherMALA's Frobenius norm error drops sharply to near-zero within ~5,000 iterations, while AdaMALA's error is still large and decreasing at 20,000 iterations — it would require "hundreds of thousands of adaptive steps" to catch up, per the authors.

This pair of 100-dimensional Gaussian experiments establishes the paper's fundamental efficiency claim: FisherMALA learns the optimal preconditioner in a few thousand iterations, while standard adaptive MCMC requires impractically many iterations to achieve comparable (and still inferior) accuracy. The consequence is not subtle — on the inhomogeneous target, AdaMALA's min ESS of ~9 means effectively no independent samples for the stiffest dimension after 20,000 iterations, while FisherMALA provides ~1,500.

Bayesian Logistic Regression: FisherMALA Dominates Across Dimensions, Including Significant Outperformance of mMALA

The Bayesian logistic regression results (Table 1 for Pima, Caravan, MNIST; Tables 3–6 in Appendix E.2 for Heart, German Credit, Australian Credit, Ripley) reveal a pattern that grows more dramatic with dimensionality:

Low-dimensional problems (Ripley, $d = 3$; Pima, $d = 7$; Heart, $d = 14$). On the smallest problems, FisherMALA, AdaMALA, HMC, and mMALA are more competitive, and FisherMALA's advantage is moderate. On Ripley ($d = 3$, Table 6): FisherMALA achieves min ESS 9245 ± 559 vs. AdaMALA's 9272 ± 412 — essentially tied, with HMC actually slightly ahead (min ESS 7645 ± 2289, though with higher variance). This is expected: at $d = 3$, AdaMALA's covariance estimation converges adequately within 20,000 iterations, and the non-Gaussianity is insufficient to create a large gap between $\Sigma$ and $I^{-1}$. On Pima ($d = 7$, Table 1): FisherMALA achieves min ESS 5629 ± 168 vs. AdaMALA's 6.4 ± 3.3 — a dramatic gap opens, with FisherMALA achieving ~880× better minimum mixing. mMALA performs comparably to FisherMALA (min ESS 5298 ± 160), and HMC performs relatively well on this problem (min ESS ~6.1).

Medium-dimensional problems (Australian Credit, $d = 15$; German Credit, $d = 25$; Caravan, $d = 87$). The gap between FisherMALA and AdaMALA widens systematically with dimensionality. On Caravan ($d = 87$, Table 1, Figure 3): FisherMALA achieves max ESS 2258 ± 45, median 1921 ± 56, min 498 ± 97. AdaMALA achieves max 42 ± 9, median 7.1 ± 0.7, min 3.1 ± 0.1 — the minimum ESS is ~160× worse. HMC does better than AdaMALA (min ESS 4.2 ± 0.5) but still ~118× worse than FisherMALA. The most striking result is the comparison with mMALA: mMALA achieves min ESS 51 ± 25 vs. FisherMALA's 498 ± 97 — a ~9.8× advantage for FisherMALA. This is the counterintuitive finding discussed in Innovation 4: mMALA uses the exact negative inverse Hessian as a position-dependent preconditioner (the theoretically "perfect" local geometry), yet FisherMALA's global inverse Fisher estimate outperforms it substantially. Figure 3 provides a qualitative complement: the log-target trace plot for FisherMALA shows rapid, low-autocorrelation exploration, while AdaMALA and other methods show visible sluggishness.

On German Credit ($d = 25$, Table 4): FisherMALA achieves min ESS 3011 ± 258 vs. AdaMALA's 3.8 ± 0.7 (~792×), HMC's 5.4 ± 1.3 (~558×), and mMALA's 2343 ± 113 (~1.3× — mMALA is competitive but still behind). On Australian Credit ($d = 15$, Table 5): FisherMALA min ESS 3772 ± 265 vs. AdaMALA's 3.3 ± 0.5 (~1,143×) and mMALA's 1869 ± 631 (~2.0×).

High-dimensional problem (MNIST, $d = 785$, Table 1, Figure 10). This is the most extreme test. FisherMALA achieves max ESS 1053 ± 36, median 812 ± 19, min 440 ± 53. AdaMALA achieves max 62 ± 9, median 8.2 ± 0.4, min 3.0 ± 0.1 — minimum ESS is ~147× worse. HMC with 10 leapfrog steps achieves its best relative performance (max 889, median 303, min 114) — still ~3.9× worse than FisherMALA in min ESS, but notably outperforming mMALA, which collapses to min ESS 5.2 ± 0.5 (~85× worse than FisherMALA). The mMALA result is the most extreme example of the "local preconditioning fails in high dimensions" phenomenon. Simple MALA achieves min ESS 2.9 ± 0.1.

The log-target trace plots (Figures 7–14 in Appendix E.2) consistently show FisherMALA's chain moving rapidly with low autocorrelation, while AdaMALA, MALA, and in some cases mMALA show visible slow drift.

Step size amplification confirms effective preconditioning (Figure 4, Appendix E.1). The paper reports the final learned values of $\sigma^2$ after burn-in to validate the theoretical connection between effective preconditioning and larger allowable step sizes. Figure 4 shows boxplots (over 10 replicates) of $\log(\sigma^2)$ for MALA, AdaMALA, and FisherMALA on the four main datasets. FisherMALA's $\sigma^2$ is substantially larger in all cases — by roughly an order of magnitude on the GP target and Pima Indians, and by nearly two orders of magnitude on Caravan and MNIST compared to MALA. AdaMALA's $\sigma^2$ is intermediate but closer to MALA than to FisherMALA. This directly supports Proposition 3's logic: a better preconditioner "symmetrizes" the target, allowing larger steps without increasing discretization error, and FisherMALA's aggressive step sizes reflect its effective learning of $I^{-1}$.

Summary of Quantitative Patterns

The ESS results across all seven datasets (Tables 1–6) reveal a consistent pattern:

1. FisherMALA is uniformly best or tied for best in min ESS on every problem. The sole exception is Ripley ($d = 3$) where AdaMALA and FisherMALA are statistically indistinguishable, which is expected since at very low dimensions, AdaMALA's covariance estimation has enough data to converge.

2. The gap between FisherMALA and AdaMALA grows with dimensionality. At $d = 3$, the methods are comparable. At $d = 7$, FisherMALA's min ESS advantage is ~880×. At $d = 87$, ~160×. At $d = 785$, ~147×. This is consistent with the diagnosis that AdaMALA's state-based covariance learning becomes increasingly data-hungry as $d$ grows, while FisherMALA's gradient-based learning extracts curvature information efficiently regardless of dimension.

3. mMALA is competitive with FisherMALA on small-to-medium problems but degrades catastrophically at $d = 785$. On Pima, German Credit, and Australian Credit, mMALA is within a factor of 1.3–2.0× of FisherMALA in min ESS. On Caravan ($d = 87$), mMALA falls to ~9.8× worse. On MNIST ($d = 785$), it falls to ~85× worse. This suggests a dimension-dependent failure mode for position-dependent preconditioning that the paper does not fully explain but flags as an important empirical finding.

4. HMC with identity mass matrix outperforms simple MALA but is significantly worse than FisherMALA on anisotropic problems. The exception is MNIST ($d = 785$), where HMC (min ESS 114) outperforms mMALA (min ESS 5.2) but still falls short of FisherMALA (min ESS 440). HMC's relative strength on MNIST may reflect the benefit of Hamiltonian dynamics in high dimensions even without a learned mass matrix, or it may reflect that the MNIST posterior has curvature structure that 10-step leapfrog integration handles better than mMALA's local linearization.

Ablation Studies and Robustness Checks

Rao-Blackwellization vs. non-Rao-Blackwellized score increments (Figure 15, Table 7 in Appendix E.3). The paper compares FisherMALA using the Rao-Blackwellized signal $s^\delta_n = \sqrt{\alpha}(s(y_n) - s(x_n))$ ("FisherMALA-with-RB") against FisherMALA using the un-marginalized signal $s^\delta_n = \mathbb{I}(u_n < \alpha)(s(y_n) - s(x_n))$ ("FisherMALA-no-RB"). Table 7 shows the two variants have essentially identical ESS performance across all seven datasets and the two Gaussian targets. For example, on the GP target: min ESS 1785 ± 104 (with RB) vs. 1794 ± 104 (no RB). On Caravan: 498 ± 97 vs. 510 ± 116. On MNIST: 440 ± 53 vs. 437 ± 40. No difference is statistically significant. Figure 15 shows the Frobenius norm trajectories on the two 100-dimensional Gaussian targets: the Rao-Blackwellized version converges marginally faster (the curves separate slightly in the first ~1,000 iterations), which "agrees with the theory that says that Rao-Blackwellization should reduce the variance," but the practical difference is negligible. The key finding is that the score increment approach itself — not the Rao-Blackwellization — is the critical design choice.

Score increments vs. paired mean-covariance stochastic estimation (Table 7). The "FisherMALA-paired-est" variant (Appendix C) uses centered score vectors with recursive mean and covariance estimation, analogous to the Haario et al. (2001) approach but applied to gradients rather than states. This performs substantially worse than the score increment approach on most problems. On the inhomogeneous Gaussian target: min ESS drops from 1501 ± 67 (with RB) to 109 ± 58. On Caravan: from 498 ± 97 to 87 ± 91. On MNIST: from 440 ± 53 to 3.2 ± 0.1 — effectively no better than simple MALA. On the GP target, the degradation is less severe (1785 vs. 1226) but still significant. Only on Ripley ($d = 3$) does paired estimation match the increment approach (9245 vs. 9217). This ablation is particularly informative because it isolates the value of using differences vs. centered values: centering removes bias from non-zero mean gradients during transience, but the noise in the running mean estimate $\bar{s}_n$ (itself a high-dimensional quantity converging at rate $1/\sqrt{n}$) propagates into the Fisher estimate, and the noise dominates the signal in moderate-to-high dimensions. The score increment approach achieves centering "for free" by differencing consecutive gradients, without requiring an explicit running mean that must be estimated.

Score increments vs. non-centered raw score vectors (Table 8 in Appendix E.3). The "FisherMALA non-centered" variant feeds raw score vectors $s_n = \nabla \log \pi(x_n)$ directly into the square-root recursion without any centering or differencing. This performs catastrophically worse on most problems. On the GP target: min ESS ~48 (vs. ~1785). On Pima Indians: min ESS ~65 (vs. ~5629). On Caravan: min ESS ~3.3 (vs. ~498). On MNIST: min ESS ~3.7 (vs. ~440). Only on Ripley ($d = 3$) does it match: min ESS ~9238 (vs. ~9245). This is a crucial negative result: it confirms the paper's diagnosis (Section 4) that the non-zero mean of raw score vectors during the transient phase introduces severe bias through the term $E[s_n]E[s_n]^\top$, which contaminates the Fisher estimate with mode-seeking direction information rather than curvature information. The bias is not merely a theoretical concern — it destroys practical performance on all but the lowest-dimensional problems.

FisherMALA vs. AdaMALA with comparable initialization budgets (Appendix D). The paper notes that AdaMALA required a more elaborate initialization (1,000 total burn-in iterations before starting covariance adaptation, vs. FisherMALA's 500) to achieve stability. This asymmetry slightly favors FisherMALA in the burn-in budget comparison. However, the Frobenius norm trajectories (Figures 1b, 2c) show that FisherMALA's advantage persists even when AdaMALA is given the full 20,000 burn-in iterations: FisherMALA converges in ~2,000–5,000 iterations, while AdaMALA is still improving at iteration 20,000. Extrapolating, AdaMALA would need roughly 5–10× more iterations to match FisherMALA's Frobenius error, and the ESS comparisons (which use equal 20,000-iteration budgets) are therefore fair in the sense that AdaMALA has had ample time to converge but simply learns too slowly.

Step size amplification under effective preconditioning (Figure 4, Appendix E.1). While not an ablation in the traditional sense, the paper reports the final learned $\sigma^2$ values to validate Proposition 3's logic that effective preconditioning enables larger step sizes. Figure 4 shows that FisherMALA's $\sigma^2$ is consistently 1–2 orders of magnitude larger than simple MALA's and substantially larger than AdaMALA's across all datasets. This is a consistency check rather than an ablation: if FisherMALA were simply a better-tuned version of the same algorithm (rather than genuinely learning a better preconditioner), the step sizes would be comparable. The large $\sigma^2$ values confirm that the learned $I^{-1}$ genuinely symmetrizes the target, allowing the algorithm to take much more aggressive steps without exceeding the acceptance rate target.

Damping parameter $\lambda$ sensitivity. The paper does not report an ablation over $\lambda$. All experiments use $\lambda = 10$. Sensitivity to this choice is therefore uncharacterized — a genuine gap. Given that $\lambda$ controls the initial regularization strength (through $\lambda/n \to 0$ as $n \to \infty$, the damping has diminishing influence), it is likely that performance is robust within a reasonable range (e.g., 1–100), but this is not demonstrated.

Learning rate for step-size adaptation $\rho_n$. The paper uses a constant $\rho_n = 0.015$ and does not ablate this choice. Standard MALA tuning typically uses constant or slowly decaying learning rates for step-size adaptation, and the paper follows this convention without sensitivity analysis.

Critical Assessment

Claim 1: The optimal preconditioner for Langevin diffusion is $I^{-1}$, and learning it produces dramatic efficiency gains over covariance-based adaptation. The experiments overwhelmingly support the practical superiority of FisherMALA over AdaMALA. The ESS comparisons (Tables 1–6) demonstrate order-of-magnitude improvements in minimum ESS on problems with $d \geq 7$, and the Frobenius norm trajectories (Figures 1b, 2c) confirm that FisherMALA learns $I^{-1}$ much faster than AdaMALA learns $\Sigma$. The step-size amplification (Figure 4) provides mechanistic validation: FisherMALA's larger allowable $\sigma^2$ confirms that its preconditioner more effectively symmetrizes the target.

However, the experiments do not isolate whether the performance difference is due to learning $I^{-1}$ vs. $\Sigma$ (different targets) or due to learning from gradients vs. states (different information sources). On the Gaussian targets, $I^{-1} = \Sigma$, so the two methods are learning the same matrix — yet FisherMALA learns it far faster. This demonstrates that gradient-based learning is more efficient than state-based learning, but does not by itself prove that $I^{-1}$ is the correct target for non-Gaussian problems. The Bayesian logistic regression results provide indirect evidence: on these non-Gaussian targets, FisherMALA's advantage over AdaMALA is even larger than on Gaussian targets (compare min ESS ratios: GP target ~3.2×, Caravan ~160×). This is consistent with $\Sigma \neq I^{-1}$ on non-Gaussian targets, so AdaMALA is both learning from a less informative signal and learning the wrong matrix. But without an ablation that learns $\Sigma$ from gradients or learns $I^{-1}$ from states, the two factors are confounded. The paper does not provide this ablation, which would be the cleanest test of whether $I^{-1}$ per se matters beyond gradient-based learning efficiency.

Claim 2: Score function increments provide a more stable and less biased adaptation signal than raw gradients or centered gradients. The ablation study (Tables 7–8) provides strong and consistent evidence. The non-centered variant (raw gradients) fails catastrophically on all problems with $d > 3$, confirming that bias from non-zero-mean gradients during transience destroys the Fisher estimate. The paired-estimation variant (centered gradients with running mean) performs better than non-centered but still degrades severely in moderate-to-high dimensions, confirming that noise in the running mean estimate propagates destructively. The Rao-Blackwellized and non-Rao-Blackwellized score increment variants perform essentially identically, showing that the differencing operation — not the RB variance reduction — is the critical design choice. These three comparisons (raw vs. centered vs. differenced) provide good internal evidence for the paper's claims about the adaptation signal. A missing ablation: what if we used score increments but re-centered them by a running mean? This would combine the transient-mean robustness of differencing with additional centering, and might reveal whether any residual bias remains in the increment approach.

Claim 3: FisherMALA matches or exceeds the performance of position-dependent Riemannian manifold MALA on high-dimensional log-concave problems. The experimental evidence is strong for the specific problems tested — Caravan ($d = 87$, FisherMALA ~9.8× better min ESS) and MNIST ($d = 785$, ~85× better) — but is limited in scope. All Bayesian logistic regression problems tested are log-concave (Gaussian prior + logistic likelihood yields a log-concave posterior), so the claim's restriction to "certain types of high-dimensional and log-concave problems" (Section 6) is accurately caveated. The paper does not test non-log-concave targets (e.g., mixtures, hierarchical models with non-conjugate priors, neural network posteriors), where local curvature information might be more valuable and mMALA's position-dependent approach might shine. The claim should therefore be understood as: there exist practically important problems (log-concave Bayesian posteriors) where a well-learned global preconditioner substantially outperforms an exact local one, not that global is universally better than local.

A significant confound in the mMALA comparison is computational cost. mMALA's $O(d^3)$ per-iteration cost means that within a fixed wall-clock budget, it completes far fewer iterations than FisherMALA. The paper's ESS comparison uses equal iteration counts (20,000 each), not equal compute time. If mMALA were given proportionally more computation time to match FisherMALA's faster per-iteration speed, would it close the gap? The paper does not address this. A wall-clock-time-matched comparison would be the fairest assessment, and its absence likely inflates FisherMALA's apparent advantage over mMALA. That said, at $d = 785$, the $O(d^3)$ vs. $O(d^2)$ difference means mMALA costs roughly $d = 785$ times more per iteration — so FisherMALA could run ~785 iterations for each mMALA iteration. Even allowing mMALA 785× more iterations to match compute, the ESS gap would need to close by that factor, which seems unlikely given that mMALA's min ESS on MNIST is ~5.2 (Table 1).

Claim 4: The $O(d^2)$ square-root recursion makes online inverse Fisher learning computationally practical. The paper provides no direct timing measurements or complexity scaling experiments. The $O(d^2)$ claim is based on operation counting (two matrix-vector products, one outer product per iteration), which is analytically correct but not empirically validated. A wall-clock timing comparison between FisherMALA's recursive square-root update and an explicit Cholesky-based approach (computing $A_n$ and factorizing) would strengthen this claim. The fact that FisherMALA runs on MNIST at $d = 785$ with 20,000 iterations demonstrates practical feasibility, but the absence of timing data means we cannot assess the constant factors — e.g., whether the recursion is actually faster in practice than a heavily optimized Cholesky decomposition at $d = 785$ (where $785^3 \approx 4.8 \times 10^8$ operations vs. $3 \times 785^2 \approx 1.8 \times 10^6$ operations — a ratio of ~270, which should dominate any constant factors).

Missing experiments that would strengthen the paper:

• Varying the damping parameter $\lambda$. All experiments use $\lambda = 10$ without justification. An ablation over $\lambda \in \{0.1, 1, 10, 100\}$ would characterize sensitivity and provide guidance for practitioners.
• Varying the initialization phase length $n_0$. The 500-iteration simple MALA initialization is used throughout. Would shorter or longer initialization materially affect the results? The paper's diagnosis of transient bias suggests this could matter.
• Non-log-concave targets. A Gaussian mixture, a hierarchical model with non-conjugate hyperpriors, or a neural network posterior would test whether the global $I^{-1}$ remains advantageous when local curvature varies substantially across the posterior.
• Wall-clock timing comparisons. Measuring actual runtime would ground the $O(d^2)$ vs. $O(d^3)$ analysis in practical terms and allow fairer comparison with mMALA and HMC.
• Higher-dimensional Gaussians. The scaling experiments stop at $d = 100$ for Gaussians and $d = 785$ for logistic regression. How does FisherMALA scale to $d = 5000$ or $10000$? Does the square-root recursion remain numerically stable?
• Comparison with quasi-Newton MALA variants. The related work section mentions quasi-Newton MCMC methods (Zhang and Sutton, 2011) but does not empirically compare against them. Since these also learn preconditioners from gradient history, they are a natural comparator.

Statistical rigor caveats. The paper reports mean ± standard deviation over 10 random initializations, which is standard for MCMC benchmarking. However, for the Bayesian logistic regression problems, variability across different dataset draws (if the data were resampled) is not assessed — the reported ESS values are conditional on a single training set. The MNIST example uses a single train/test split (5 vs. 6 classification), so generalizability to other digit pairs or other classification tasks is uncharacterized. The Frobenius norm trajectories (Figures 1a-b, 2c) appear to be from single representative runs rather than averaged over replications, so the variability in convergence speed across random seeds is not visible.

Overall assessment. The experimental section provides strong evidence for FisherMALA's practical superiority over AdaMALA and simple MALA across a range of Gaussian and Bayesian logistic regression targets, with particularly compelling results in moderate-to-high dimensions where state-based adaptation fails. The ablation study convincingly isolates score function incrementation as the critical design choice. The comparison with mMALA reveals an important and counterintuitive finding about global vs. local preconditioning on log-concave problems, though the claim should be tempered by the lack of compute-matched comparisons and the restriction to log-concave targets. The main weaknesses are the absence of timing data, the narrow range of target distributions tested (all log-concave), and the lack of sensitivity analysis for the key hyperparameters $\lambda$ and $n_0$.

6. Limitations and Trade-offs

The Difficulty Estimation Cost Is Not Amortized in the Reported Efficiency Gains

The assumption or constraint. The entire FisherMALA framework requires evaluating gradients $s_n = \nabla \log \pi(x_n)$ to drive the preconditioner adaptation, plus an initial 500-iteration burn-in phase with simple MALA before the adaptation signal begins. While each gradient evaluation is already required by MALA's proposal (making the marginal adaptation cost $O(d^2)$ per iteration), the paper does not account for the quality of the learned preconditioner as a function of computational budget. Specifically, on the problems where FisherMALA shines — high-dimensional anisotropic targets — the paper reports that FisherMALA converges to the optimal preconditioner in "a few thousand iterations" (Section 5.2, Figure 1b showing ~2,000 iterations for the GP target) while AdaMALA requires "hundreds of thousands of adaptive steps" (Section 5.2). But the paper's ESS comparisons use equal iteration counts (20,000 burn-in + 20,000 collection for all methods). This means FisherMALA is being compared against AdaMALA at a budget where AdaMALA demonstrably has not converged (Figure 1b shows AdaMALA's Frobenius error still declining at iteration 20,000). The comparison demonstrates FisherMALA learns faster per iteration, but it does not answer: if AdaMALA were given enough iterations to reach the same Frobenius error as FisherMALA, would the ESS gap close?

The consequence. A practitioner deciding between FisherMALA and AdaMALA based on the paper's ESS tables would conclude that FisherMALA is categorically 10–160× better in min ESS (Tables 1, 3–6). But this conclusion conflates asymptotic quality of the learned preconditioner with learning speed at a fixed budget. If AdaMALA were run for 200,000 iterations (10× the paper's budget), it might eventually learn a preconditioner of comparable quality to FisherMALA's, and the ESS gap during the collection phase might narrow substantially. The paper explicitly shows (Figure 2c) that AdaMALA's Frobenius error is still declining — it has not plateaued — so 20,000 iterations is an arbitrary cutoff, not a budget at which AdaMALA has reached its asymptotic performance. This is particularly relevant because AdaMALA is a simpler algorithm with lower per-iteration overhead (no square-root recursion, no gradient differencing), so it might actually have lower wall-clock cost per iteration that partially compensates for slower per-iteration learning.

What evidence exists in the paper. The Frobenius norm trajectories in Figures 1b (GP target) and 2c (inhomogeneous Gaussian target) directly show AdaMALA's error declining throughout the full 20,000 burn-in iterations without reaching an asymptote. For the GP target (Figure 1b), AdaMALA's error at iteration 20,000 is still visibly above FisherMALA's error at iteration ~1,000. For the inhomogeneous Gaussian target (Figure 2c), the gap is even larger and AdaMALA's trajectory shows no sign of saturating. The paper acknowledges this implicitly in Section 5.2: "AdaMALA can eventually get very close to the optimal precondtioner but it requires hundred of thousands of adaptive steps, while FisherMALA learns it with only few thousand steps." But this acknowledgement is not reflected in the ESS comparison methodology, which uses equal iteration budgets rather than equal preconditioner quality.

Mitigation status. The paper does not attempt to run AdaMALA for longer to assess its asymptotic ESS, nor does it perform a budget-matched comparison where FisherMALA gets fewer iterations (compensating for its higher per-iteration overhead) or AdaMALA gets proportionally more. The authors note this asymmetry in Appendix D (AdaMALA required 1,000 initialization iterations vs. FisherMALA's 500) but frame it as making initialization harder for AdaMALA, not as a budget mismatch. The $O(d^2)$ vs. $O(d^2)$ per-iteration cost for both methods (the state-based recursion in Equation 17 is also $O(d^2)$) means wall-clock comparison on equal iterations would be approximately fair — but the key issue is that equal iterations do not correspond to equal preconditioner quality, and the ESS numbers therefore reflect a mixture of learning quality and learning speed that is difficult to disentangle. A budget-escalation experiment (reporting ESS at 5K, 10K, 20K, 50K, 100K burn-in for both methods) would clarify whether FisherMALA's advantage saturates or widens with budget. This is a significant omission given that the paper's core claim is about the superiority of learning $I^{-1}$ over learning $\Sigma$, not merely about faster convergence.

---

No Theoretical Guarantee That Score Increments Provide an Unbiased Fisher Estimate

The assumption or constraint. The paper's entire adaptation mechanism rests on the claim that outer products of score function increments $s^\delta_n (s^\delta_n)^\top$ provide a consistent estimator of the Fisher information matrix $I = E_\pi[\nabla \log \pi \nabla \log \pi^\top]$. The algorithmic derivation (Section 4) treats $s^\delta_n$ as interchangeable with $s_n$ in the empirical Fisher recursion (Equations 8–11), implicitly assuming that $E[s^\delta_n (s^\delta_n)^\top] \propto I$ when the chain is at stationarity, and that the bias during transience is small enough to be negligible. Neither claim is proven. The paper explicitly acknowledges this gap in Section 7: "it would be useful to investigate whether the score function differences that we use as the adaptation signal introduce any bias in the estimation of the inverse Fisher matrix."

The consequence. Without an unbiasedness or consistency proof, practitioners cannot know whether FisherMALA's adaptation is asymptotically correct — i.e., whether running the algorithm indefinitely would converge to $I^{-1}$ (up to scaling) or to some other matrix that happens to work well on the tested problems. More critically, there may exist target distributions where the score increment estimator is systematically biased in ways that the paper's experiments do not reveal. For example, if the chain's consecutive states $x_n$ and $x_{n+1}$ are highly correlated (as they are in any MCMC sampler), the outer product of their gradient difference encodes not just local curvature of $\log \pi$ but also the correlation structure of the chain itself. This could introduce a dependence on the sampler's mixing properties into the Fisher estimate — creating a feedback loop where poor mixing leads to biased Fisher estimation, which leads to poor preconditioning, which perpetuates poor mixing. The paper's empirical results suggest this is not a problem for the targets tested, but the absence of theory means there is no way to predict when it might become a problem.

What evidence exists in the paper. The Gaussian target experiments (Section 5.2, Figures 1b, 1d, 2c) provide the strongest empirical evidence that the score increment estimator converges to the correct $I^{-1}$, because for Gaussian targets $I^{-1} = \Sigma$ is known analytically and FisherMALA's learned preconditioner matches it with vanishing Frobenius error. The Frobenius norm trajectories show the error declining to near-zero and stabilizing, which is consistent with consistency (though not a proof). However, this evidence is limited to two parametric families (Gaussian and log-concave Bayesian logistic regression). There is no experiment on a target with known $I^{-1} \neq \Sigma$ where $I^{-1}$ is analytically computable — this would be the gold-standard test of whether FisherMALA learns $I^{-1}$ specifically rather than some other useful matrix. The step-size amplification evidence (Figure 4) is a consistency check: larger $\sigma^2$ implies the preconditioner is effectively symmetrizing the target, but does not verify that it is symmetrizing with $I^{-1}$ specifically rather than with some other matrix that also enables large steps.

Mitigation status. The paper does not attempt to prove unbiasedness or consistency, nor does it provide any diagnostic for when the score increment estimator might fail. The acknowledgment in Section 7 is a candid flag for future work, but it leaves the method without theoretical foundations that would be expected for an algorithm that makes a strong optimality claim (Proposition 3). The empirical success on the tested targets is reassuring but does not substitute for theory, particularly because the targets are all relatively well-behaved (strongly log-concave, unimodal). A counterexample — a target where score increments produce a biased Fisher estimate — would be highly informative but is not provided. This gap is especially acute given that the paper's methodological contribution is precisely to replace the theoretically well-understood state-based covariance estimator (whose consistency under ergodicity is standard) with a new gradient-difference-based estimator whose properties are unknown.

---

The Learned Preconditioner Is Global; It Cannot Adapt to Heterogeneous Local Geometry

The assumption or constraint. FisherMALA learns a single global preconditioner $A \propto I^{-1}$ that is applied uniformly across the entire state space. The Fisher matrix $I = E_\pi[\nabla \log \pi \nabla \log \pi^\top]$ averages curvature information over the stationary distribution, producing a preconditioner that reflects the average geometry of $\log \pi$. This design choice is explicit in Section 6: "$I^{-1}$ we consider here requires an expectation under the target and thus it is always a global preconditioner rather than a position-dependent one." The paper contrasts this with Riemannian manifold methods that use position-dependent preconditioning $A(x)$ adapting to local curvature at each point.

The consequence. For target distributions where the local geometry varies substantially across the typical set, a global preconditioner may be systematically suboptimal — it will be too aggressive in regions that are more curved than average (causing high rejection rates) and too conservative in regions that are flatter than average (causing slow exploration). The paper's experiments demonstrate that this is NOT a fatal problem for the log-concave targets tested — indeed, FisherMALA outperforms the position-dependent mMALA on Caravan and MNIST — but the experiments do not establish when global preconditioning breaks down. Multimodal targets, hierarchical models with varying local dependence structures, or posteriors with pronounced non-stationary curvature (where the Hessian varies qualitatively across regions) are plausible failure modes. In a bimodal distribution, for example, the Fisher matrix $I$ averages gradient outer products from both modes, producing a preconditioner that is well-suited to neither — the local geometry in each mode may suggest very different optimal scalings, and a global average compromises both. The paper's Gaussian targets, Bayesian logistic regression, and even the inhomogeneous Gaussian all share the property that the curvature is relatively homogeneous across the high-probability region: the Hessian of $\log \pi$ either is constant (Gaussian), or varies smoothly and modestly (logistic regression is log-concave with continuous Hessian). These are exactly the cases where a global average should work well. The paper does not test targets with genuinely heterogeneous local geometry, so the scope of the "global preconditioner works well" claim is empirically narrow.

What evidence exists in the paper. The comparison with mMALA on Caravan and MNIST (Table 1) provides counterevidence to the intuition that local preconditioning is always superior — on these targets, mMALA's exact local Hessian-based preconditioning performs worse than FisherMALA's global inverse Fisher. The paper hypothesizes (Section 6) that "position-dependent preconditioning may be less effective in certain types of high-dimensional and log-concave problems." This is framed as a finding about mMALA's limitations, not about FisherMALA's. But the same data also imply that FisherMALA was tested only on problems where global preconditioning is sufficient — there is no experiment where FisherMALA fails due to heterogeneity and a local method succeeds. The paper does not attempt to characterize the conditions under which global preconditioning is adequate vs. when local adaptation becomes necessary, which is precisely the information a practitioner needs to choose between FisherMALA and mMALA or HMC with a learned mass matrix.

Mitigation status. The paper does not address this limitation beyond the brief comment in Section 6. There is no diagnostic for detecting when the global preconditioner is ill-suited (e.g., monitoring whether acceptance rates vary systematically across the state space, which would indicate position-dependent mismatch). There is no hybrid scheme that uses the global $I^{-1}$ as a baseline and adds local corrections — for instance, using $I^{-1}$ as the mass matrix in HMC (as suggested in Section 7 for future work) would combine global scaling with Hamiltonian dynamics' ability to handle local curvature. The limitation is fundamental to the method's design (global preconditioner) rather than an implementation detail, so it cannot be mitigated without extending the algorithm to a position-dependent or multi-scale form.

---

Tested Only on Log-Concave Targets; Performance on Multimodal, Non-Convex, or Heavy-Tailed Distributions Is Unknown

The assumption or constraint. Every target distribution in the paper's experimental evaluation falls into one of two families: multivariate Gaussian (Section 5.2, three synthetic examples) and Bayesian logistic regression with Gaussian prior (Section 5.3, seven datasets). Bayesian logistic regression with a Gaussian prior is log-concave: the log-likelihood for logistic regression is concave, the log-prior $-\frac{1}{2}\|\theta\|^2$ is strictly concave, and their sum is therefore strongly log-concave. This means every target tested is unimodal with a single, well-defined typical set whose geometry (captured by the Fisher matrix $I$) is representative of the entire posterior. The paper makes no secret of this scope — the targets are clearly described — but it does not discuss the implications of this restriction.

The consequence. For practitioners whose inference problems involve multimodality, non-convexity, or heavy tails, the paper provides essentially no guidance on whether FisherMALA will work. The method's core assumptions could plausibly break in these settings:

• Multimodal targets: The Fisher matrix $I = E_\pi[\nabla \log \pi \nabla \log \pi^\top]$ averages over all modes. If the modes have different curvature structures (e.g., one mode is narrow, another is broad), $I^{-1}$ will reflect a compromise that may be well-suited to neither mode. Moreover, the score increment adaptation signal $s^\delta_n = \sqrt{\alpha}(s(y_n) - s(x_n))$ operates locally — differences between gradients within a single mode provide information about that mode's curvature, but differences between gradients in different modes (which occur only when the chain successfully jumps between modes) conflate curvature with mode-separation information. During the transient phase when the chain is trapped in one mode, the Fisher estimate will reflect only that mode's geometry, potentially creating a self-reinforcing trap: the preconditioner adapts to the local mode, making proposals well-tuned to that mode but potentially making inter-mode jumps even harder.

• Non-convexity: If $\log \pi$ has regions of negative curvature (where the Hessian has negative eigenvalues), the Fisher matrix $I$ — being a covariance of gradients — remains positive definite by construction, but its relationship to local geometry becomes more complex. The negative curvature directions contribute to $I$ through gradient magnitude but the sign information is lost, so $I^{-1}$ may not reflect the true difficulty of exploration in non-convex regions.

• Heavy tails: For distributions with polynomial rather than exponential tail decay (e.g., Student-t), gradients $\nabla \log \pi(x)$ decay slowly in the tails, meaning $s_n s_n^\top$ can receive large contributions from tail regions that the chain visits rarely. The empirical Fisher estimate may be dominated by a small number of tail visits, making it noisy and unstable.

What evidence exists in the paper. None. The paper does not evaluate on any non-log-concave target, nor does it discuss these failure modes. The Bayesian logistic regression problems use a standard Gaussian prior $\mathcal{N}(0, I_d)$, guaranteeing log-concavity. The paper does not test with heavy-tailed priors (e.g., Laplace, Student-t), hierarchical priors, mixture models, or any target with multiple modes. The MNIST example at $d = 785$ is the largest and most challenging, but it is still a log-concave Bayesian logistic regression — the dimensionality provides a test of computational scaling, not of robustness to distributional complexity. This is a significant limitation because realistic Bayesian models (hierarchical models, mixture models, neural network posteriors) routinely violate log-concavity.

Mitigation status. The paper does not acknowledge this limitation explicitly, which is an omission. Section 7 suggests extending the method to HMC and exploring multiple interacting chains, but does not discuss the need to evaluate on non-log-concave targets. The related work discussion (Section 6) mentions that position-dependent methods like mMALA "have been theoretically analyzed by obtaining convergence guarantees" (Chewi et al., 2020; Li et al., 2022), and FisherMALA lacks analogous theory even for the log-concave case, let alone beyond it. A practitioner reading the paper would be justified in concluding that FisherMALA is effective for Bayesian logistic regression and Gaussian targets, but would have no basis for extrapolating to their own non-log-concave posterior. Given that the paper presents FisherMALA as a general-purpose adaptive MCMC method (the abstract claims it "significantly outperforms other methods" without qualification), the restriction to log-concave targets is a consequential scope limitation that should have been discussed.

---

Numerical Stability of the Square-Root Recursion at Scale Is Not Characterized

The assumption or constraint. The square-root recursion in Proposition 4 (Equations 12–13) maintains $R_n$ such that $R_n R_n^\top = A_n$ through a sequence of rank-one downdates. Each iteration computes:

$\phi_n = R_{n-1}^\top s^\delta_n, \quad r_n = \frac{1}{1 + \sqrt{\frac{1}{1 + \phi_n^\top \phi_n}}}, \quad R_n = R_{n-1} - r_n \frac{(R_{n-1} \phi_n) \phi_n^\top}{1 + \phi_n^\top \phi_n}$

This recursion assumes exact arithmetic. In finite-precision floating-point, the accumulated rounding errors after thousands of rank-one updates can cause $R_n$ to drift away from satisfying $R_n R_n^\top = A_n$, potentially losing positive definiteness (which would be catastrophic for sampling since $R_n$ is used to transform standard normal noise). The paper acknowledges this implicitly by adopting a technique from Kalman filtering (Potter and Stern, 1963; Bierman, 1977) where square-root filtering was specifically developed to improve numerical stability over directly maintaining the covariance — but the classical Kalman filtering context involves small state dimensions (typically < 100) and does not guarantee stability for arbitrary sequences of rank-one updates at $d = 785$ over 20,000 iterations.

The consequence. A practitioner implementing FisherMALA from the paper's description could encounter silent numerical degradation: $R_n$ might drift to a matrix that is not positive definite or whose product $R_n R_n^\top$ diverges significantly from the true $A_n$ that would be obtained from explicit inverse-and-factorize computation. This could manifest as proposals that are inappropriately scaled, poor mixing despite the adaptation appearing to converge, or outright crashes (when $R_n$ becomes numerically singular and the Cholesky-like sampling step fails). The problem is likely to be dimension-dependent and precision-dependent — experiments at $d = 100$ in double precision might show no issues, while $d = 1000$ in single precision could fail silently.

What evidence exists in the paper. None directly. The paper demonstrates that FisherMALA works at $d = 785$ (MNIST, Table 1) and reports ESS scores averaged over 10 runs with standard deviations. The standard deviations (e.g., min ESS 440 ± 53 on MNIST) suggest that all 10 runs completed successfully and produced similar results, which is evidence that catastrophic numerical failure did not occur — but it does not rule out slow degradation that reduces efficiency without causing outright crashes. The paper does not report diagnostics that would reveal numerical issues: no comparison between the recursively maintained $R_n$ and a periodically recomputed Cholesky factor of the explicitly formed $A_n$, no eigenvalue checks, no condition number monitoring. The Frobenius norm comparisons (Figures 1b, 2c) are computed from the normalized matrix $\tilde{A}_n$ (not from $R_n$ directly), and the paper does not specify whether these are computed by forming $A_n = R_n R_n^\top$ and normalizing, or by an independent computation. If $R_n$ were numerically degraded, forming $A_n = R_n R_n^\top$ might still produce a reasonable-looking matrix (since the degradation could be in directions that cancel in the product), masking the issue.

Mitigation status. The paper does not discuss numerical stability, does not provide a recommended precision (float32 vs. float64), and does not suggest periodic re-orthogonalization or re-computation strategies. The square-root filtering literature (Bierman, 1977) discusses stabilization techniques (e.g., using triangular square roots and applying Givens rotations for updates), but the paper's recursion uses full (non-triangular) $R_n$, which is less stable than the triangular formulations in the original Kalman filtering context. A practitioner who naively implements Algorithm 1 as written may encounter stability problems at scale that the paper does not prepare them for. This is a practical limitation rather than a conceptual one — it does not undermine the theoretical optimality result, but it means the algorithm as described may require additional engineering (stability monitoring, periodic recomputation, use of triangular factors) to be production-ready beyond the tested dimensions. The paper could have strengthened its contribution substantially by including even a brief discussion of numerical considerations and recommended practices.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper reorients the adaptive MCMC literature around a deceptively simple but previously unasked question: what is the optimal preconditioner for gradient-based MCMC, derived from first principles rather than inherited from gradient-free methods? The answer — the inverse Fisher information matrix $I^{-1}$ rather than the target covariance $\Sigma$ — is not a marginal refinement but a conceptual correction to two decades of adaptive MCMC practice. Since Haario et al. (2001), the field has operated under the implicit assumption that if $\Sigma$ is the right preconditioner for random walk Metropolis (where proposals are blind to local geometry), it must also be right for MALA (where proposals incorporate gradient information). Proposition 3 shows this assumption is mathematically incorrect for non-Gaussian targets: the gradient structure encodes curvature information through the Fisher matrix that the state covariance cannot capture, and optimizing the expected squared jumped distance — a criterion that directly penalizes discretization error — naturally selects $I^{-1}$ regardless of $\Sigma$.

The magnitude of this shift is best understood by analogy to the distinction between the Hessian and the inverse covariance in optimization. In Newton's method, the Hessian $-\nabla^2 \log \pi(x)$ provides local curvature information that preconditions the gradient step; in MCMC, the expected Hessian $I = -E_\pi[\nabla^2 \log \pi]$ provides global curvature information that preconditions the Langevin diffusion. The covariance $\Sigma$, by contrast, measures spread, not curvature. For Gaussians these coincide, which explains why the covariance heuristic persisted — it works asymptotically in the large-sample limit where posteriors become approximately Gaussian. But for genuinely non-Gaussian targets, confusing spread with curvature leads to suboptimal preconditioning, exactly as the paper demonstrates on Bayesian logistic regression (Table 1: AdaMALA min ESS ~3.1 vs. FisherMALA ~498 on Caravan; ~3.0 vs. ~440 on MNIST).

This reframing has direct methodological consequences. It redirects adaptive MCMC research from state-based to gradient-based learning. The paper's exhaustive ablation of adaptation signals (Tables 7–8, Appendix E.3) provides a clear empirical ranking: score function increments > Rao-Blackwellized increments > paired mean-covariance estimation > raw score vectors. The catastrophic failure of raw gradients (min ESS ~48 on GP target vs. ~1785 for increments, ~3.7 vs. ~440 on MNIST) and the severe degradation of centered estimation (~109 vs. ~1500 on inhomogeneous Gaussian, ~3.2 vs. ~440 on MNIST) establish that how gradient information is processed matters fundamentally. This is not a detail — it is the difference between a working algorithm and complete failure. The paper thus provides both a theoretical target ($I^{-1}$) and a practical recipe (score increments, square-root recursion) that future adaptive gradient-based samplers can build on.

It resolves the tension between global and local preconditioning in a subtle way. The paper's most surprising empirical result — that a learned global $I^{-1}$ can substantially outperform an exact position-dependent Riemannian preconditioner on high-dimensional log-concave targets (mMALA min ESS ~51 vs. FisherMALA ~498 on Caravan; ~5.2 vs. ~440 on MNIST) — challenges the intuitive assumption that more local geometric information is always better. This does not mean local preconditioning is useless; it means that the bias-variance tradeoff in preconditioner estimation matters at least as much as the theoretical optimality of the preconditioner being estimated. mMALA's local Hessian, while exact for the current point, varies across the state space and may be poorly conditioned or noisy in high dimensions; FisherMALA's global $I^{-1}$ averages over the stationary distribution, trading local precision for robustness and faster computation ($O(d^2)$ vs. $O(d^3)$ per iteration). This reframes the global-vs-local debate as an estimation-vs-approximation tradeoff rather than a simple hierarchy where local methods are strictly better. For practitioners, it means a well-tuned global method can be the pragmatically superior choice even when local methods are theoretically available.

It makes a case that the adaptation signal matters more than the adaptation target. On Gaussian targets where $I^{-1} = \Sigma$ (so both FisherMALA and AdaMALA are targeting the same matrix), FisherMALA still learns the preconditioner orders of magnitude faster (Figures 1b, 2c: ~2,000 iterations to convergence vs. AdaMALA still unconverged at 20,000). This cleanly isolates the information source (gradients vs. states) from the estimation target ($I^{-1}$ vs. $\Sigma$), and shows that gradient-based learning is dramatically more sample-efficient. This finding has implications beyond MCMC: any method that currently learns covariance matrices from state histories (e.g., adaptive importance sampling, variational inference with full-covariance Gaussians, evolutionary strategies in optimization) could potentially accelerate learning by exploiting gradient information through analogous differencing schemes.

Research directions that become more attractive:

• Gradient-history methods for other samplers. The paper explicitly flags HMC's mass matrix as a target for Fisher-style learning (Section 7). The square-root recursion provides $R_n$ for the proposal, but HMC needs both the mass matrix $M$ (analogous to $A$) and its inverse $M^{-1}$ for the kinetic energy. Extending the recursion to maintain both (or to maintain a factorization that provides both cheaply) is a natural next step, and the paper's demonstration that FisherMALA matches mMALA's performance on Gaussians suggests Fisher-HMC could dramatically improve HMC's performance on anisotropic targets where identity-mass HMC currently collapses (Table 1: HMC min ESS ~4.2 on Caravan vs. FisherMALA ~498).

• Verifier and test-function design for MCMC diagnostics. The Frobenius norm trajectory (Figures 1b, 2c) — possible only when the optimal preconditioner is known — is a powerful diagnostic for adaptation quality. For general targets where $I^{-1}$ is unknown, the paper's approach suggests using the learned preconditioner's convergence behavior (e.g., stabilization of $\text{tr}(R_n R_n^\top)$ or of the effective step size $\sigma^2_R$) as an adaptation diagnostic. This connects to the broader question of when to stop adapting in adaptive MCMC.

Research directions that become less attractive:

• Purely state-based adaptive MCMC for gradient-based samplers. The paper's evidence that covariance-based adaptation is orders of magnitude slower and asymptotically suboptimal for non-Gaussian targets suggests that future adaptive MALA (and likely adaptive HMC) research should focus on gradient-based schemes. The covariance recursion (Haario et al., 2001) is not obsolete — it remains relevant for RWM and for initialization — but as a primary adaptation mechanism for gradient-based samplers, it has been superseded.

• Developing ever-more-complex position-dependent Riemannian methods without addressing computational cost. The mMALA comparison (Table 1) shows that a well-learned global preconditioner can dominate an exact local one on the problems tested. If this pattern generalizes, research effort may be better spent on improving global preconditioner estimation (faster convergence, better robustness, lower cost) than on engineering approximate Hessians for local methods. The paper does not prove global is universally better, but it shifts the burden of proof: a new local preconditioning method should now demonstrate superiority over a well-tuned global Fisher method, not just over isotropic baselines.

Follow-Up Research This Work Enables

1. Rigorous bias and consistency analysis of the score increment estimator. The paper's central open question (Section 7) is whether the score increment outer products $s^\delta_n (s^\delta_n)^\top$ provide an unbiased or consistent estimator of the Fisher matrix $I$. The Gaussian experiments (Figures 1b, 2c) provide empirical evidence of convergence, but no proof. A strong theoretical contribution would: (a) characterize the bias $E[s^\delta_n (s^\delta_n)^\top] - I$ as a function of the chain's autocorrelation time and the local non-linearity of $\nabla \log \pi$, (b) establish conditions under which the bias vanishes asymptotically (e.g., as $\sigma^2 \to 0$), and (c) provide a counterexample — a target distribution and sampler configuration where the score increment estimator converges to the wrong matrix. The bias is likely related to the fact that consecutive MCMC states $x_n$ and $x_{n+1}$ are correlated, so their gradient difference reflects both local curvature and chain dynamics. A natural starting point would be the Gaussian case where everything is analytically tractable: compute the expected outer product of score increments for a Langevin diffusion with a given preconditioner, and compare to the true Fisher matrix. This would quantify how the bias depends on the step size $\sigma^2$ and the mismatch between the current and optimal preconditioners.

2. Fisher-HMC: extending the square-root adaptation to Hamiltonian Monte Carlo mass matrices. Section 7 flags HMC extension as future work, and the paper provides all the necessary building blocks. HMC requires both a mass matrix $M$ (playing the role of $A$ in the Langevin diffusion) and its inverse $M^{-1}$ (for the kinetic energy term in the Hamiltonian). The challenge is that FisherMALA's square-root recursion maintains $R_n$ with $R_n R_n^\top = A_n$, but generating HMC proposals requires sampling momenta $p \sim \mathcal{N}(0, M)$ and computing the kinetic energy $\frac{1}{2} p^\top M^{-1} p$. A natural approach is to maintain both $R_n$ and its inverse $R_n^{-1}$ using a similar rank-one recursion, or to maintain $R_n$ and compute $R_n^{-1} \eta$ implicitly via solving a triangular system if $R_n$ is maintained in Cholesky factor form. A concrete experiment: replicate the GP target and inhomogeneous Gaussian target experiments with Fisher-HMC (10 leapfrog steps, mass matrix adapted using the same score increment signal) and compare against FisherMALA and standard HMC. The paper's existing HMC baseline (min ESS ~4.2 on Caravan, ~114 on MNIST) provides a clear low bar; FisherMALA's min ESS ~498 on Caravan, ~440 on MNIST provides a target. Given HMC's known superiority over MALA on well-conditioned problems, Fisher-HMC could combine HMC's suppression of random-walk behavior with FisherMALA's anisotropy handling, potentially producing state-of-the-art sampling efficiency on the tested problems.

3. Active difficulty estimation and dynamic preconditioner switching for heterogeneous targets. The paper tests only targets where the global Fisher matrix adequately captures the geometry (all targets are log-concave and relatively homogeneous). A natural stress test is a mixture of Gaussians with different covariance structures: $\pi(x) = \frac{1}{2} \mathcal{N}(x | \mu_1, \Sigma_1) + \frac{1}{2} \mathcal{N}(x | \mu_2, \Sigma_2)$ where $\Sigma_1$ and $\Sigma_2$ have different anisotropy patterns (e.g., $\Sigma_1$ is stretched along dimension 1, $\Sigma_2$ along dimension 2). The global Fisher matrix $I$ averages $\nabla \log \pi \nabla \log \pi^\top$ over both modes, producing a compromise preconditioner that is optimal for neither mode. This would reveal whether FisherMALA's global preconditioner creates mode-trapping — if the chain gets stuck in one mode, the score increments reflect only that mode's curvature, and the preconditioner may adapt to become locally optimal for that mode, making inter-mode jumps even less likely. A follow-up could develop a mixture-of-preconditioners approach: maintain a small set (e.g., 2–5) of candidate preconditioners, adapt each using score increments from the mode it's associated with, and dynamically select which to use based on the chain's current location. This connects to the difficulty-estimation theme of adaptive test-time compute allocation — the preconditioner becomes difficulty-dependent, with "difficulty" defined by which posterior mode the chain is exploring.

4. Scaling laws for preconditioner learning: how does the required adaptation budget grow with dimension and anisotropy? The paper shows FisherMALA learning $I^{-1}$ in ~2,000 iterations at $d = 100$ (GP target, Figure 1b) but does not characterize how this scales. A systematic study could vary: (a) dimensionality from $d = 10$ to $d = 1000$ on Gaussian targets with controlled eigenvalue spectra, (b) condition number $\kappa = \lambda_{\max} / \lambda_{\min}$ of the target covariance, and (c) the learning rate schedule (constant, $1/n$, $1/\sqrt{n}$). The output would be a scaling law of the form: iterations to reach Frobenius error $\epsilon$ grows as $O(d^\alpha \kappa^\beta)$ for FisherMALA vs. AdaMALA. The paper's existing data provides two points: at $d = 2$ (Figure 1a), both methods converge quickly; at $d = 100$ (Figure 1b), FisherMALA converges in ~2,000 iterations while AdaMALA is still unconverged at 20,000 — suggesting at least a $10\times$ gap, but the functional form is unknown. Such scaling laws would provide practical guidance for practitioners (how long should I run burn-in for my problem?) and would connect to the broader scaling laws literature in machine learning.

5. Computational cost amortization: when does the $O(d^2)$ recursion beat explicit $O(d^3)$ inversion? The paper argues that the square-root recursion makes FisherMALA practical by reducing per-iteration cost from $O(d^3)$ to $O(d^2)$, but provides no wall-clock measurements. A practical contribution would benchmark: (a) FisherMALA with the recursive square-root update vs. FisherMALA with explicit Cholesky (compute $A_n$ via Woodbury, then Cholesky factorize) at dimensions $d \in \{50, 100, 200, 500, 1000, 2000\}$, measuring both per-iteration time and total time to achieve a target ESS. This would identify the crossover dimension where the recursion becomes essential. Additionally, it would test numerical stability: at what dimension does the recursive $R_n$ start to lose positive definiteness in float32 vs. float64? The classical Kalman filtering literature (Bierman, 1977) recommends triangular square-root formulations (UD factorization) for stability; comparing the paper's full-matrix $R_n$ against a triangular variant would provide implementation guidance. This may seem like engineering rather than science, but it directly determines whether FisherMALA is deployable at the scales where it's most needed (modern Bayesian models with $10^3$–$10^5$ parameters).

6. Gradient-difference adaptation for other Monte Carlo and optimization methods. The paper's key insight — that gradient differences provide a low-bias, high-efficiency signal for learning local geometry — is not specific to MALA. A natural extension is to stochastic gradient MCMC (SGLD, SGHMC), where gradients are noisy estimates from minibatches. In that setting, score increments $s^\delta_n = \nabla \log \tilde{\pi}(x_{n+1}) - \nabla \log \tilde{\pi}(x_n)$ involve two noisy gradient estimates whose noise may partially cancel, potentially providing a more stable Fisher estimate than raw noisy gradients. A concrete experiment: compare Fisher-preconditioned SGLD against standard SGLD and against RMSprop-preconditioned SGLD (which uses a diagonal preconditioner estimated from gradient squares, analogous to Adam in optimization) on Bayesian neural network posteriors of varying size. This would test whether the full-matrix Fisher approach retains its advantage when gradients are stochastic, and whether the $O(d^2)$ cost remains acceptable when $d$ is in the tens of thousands. A negative result (the recursion becomes numerically unstable or the advantage over diagonal preconditioning vanishes at scale) would be equally informative, clarifying the method's scope.

Practical Applications and Downstream Use Cases

High-dimensional Bayesian inference in probabilistic programming. The most immediate application is accelerating MCMC inference in probabilistic programming languages (Stan, PyMC, TensorFlow Probability, NumPyro). These systems currently offer HMC with a diagonal or dense mass matrix adapted during warm-up, but the adaptation is state-based and can be slow for anisotropic posteriors. Integrating FisherMALA would provide two concrete benefits: (a) faster warm-up — the paper's Figures 1b and 2c show FisherMALA converging to the optimal preconditioner in ~2,000 iterations while AdaMALA requires >20,000, which translates directly to shorter warm-up phases and less wasted computation; (b) higher ESS per iteration — on the Caravan dataset ($d = 87$), FisherMALA achieves min ESS ~498 vs. HMC's ~4.2, a >100× improvement that means reliable posterior summaries can be obtained with 1% of the post-warm-up iterations. For practitioners running models with thousands of parameters on expensive likelihoods, this could reduce inference time from hours to minutes.

Large-scale logistic regression and generalized linear models. The paper's Bayesian logistic regression experiments (Section 5.3, Tables 1, 3–6) demonstrate strong performance on exactly the type of model that arises routinely in medical statistics, econometrics, and social science: regression with moderate parameter counts ($d = 15$–$100$) and moderate sample sizes ($n = 250$–$11,000$). These are problems where MCMC is standard but can be slow when predictors are correlated or on different scales. FisherMALA's demonstrated ability to handle unstandardized inputs (the paper deliberately does not standardize to create anisotropic posteriors) is particularly valuable — in practice, analysts often standardize predictors as a pre-processing step specifically to improve MCMC mixing, which changes the interpretation of coefficients and complicates prior specification. FisherMALA could eliminate this pre-processing step, allowing inference directly on the original scale without mixing penalties.

Gaussian process hyperparameter inference. The GP correlated target (Section 5.2, Figure 1) is a direct analog of Gaussian process covariance hyperparameter posteriors, which are notoriously difficult to sample due to strong correlations and varying sensitivities among lengthscale, amplitude, and noise parameters. The paper shows FisherMALA matching the theoretically optimal mMALA on a 100-dimensional GP target (min ESS ~1785 vs. ~1842), demonstrating that the learned global preconditioner captures the complex correlation structure encoded in the non-stationary covariance function. For practitioners fitting GPs to spatial or temporal data, FisherMALA offers a drop-in replacement for the standard HMC or MALA samplers in GPyTorch, GPflow, or similar packages, with the potential to reduce the number of likelihood evaluations needed for convergence — a critical savings since each evaluation costs $O(n^3)$ for exact GPs.

Warm-starting Riemannian samplers and variational inference. The paper's global $I^{-1}$ could serve as an initialization or regularizer for more expensive local methods. For example, Riemannian manifold MALA (which the paper shows degrades on Caravan and MNIST) might benefit from using FisherMALA's learned $I^{-1}$ as a baseline preconditioner, adding local Hessian corrections only when they differ substantially from the global average. This hybrid approach could preserve mMALA's theoretical advantages on problems with heterogeneous curvature while avoiding its catastrophic degradation on problems where local Hessian estimates are noisy. Similarly, variational inference with full-covariance Gaussians requires learning a covariance matrix — typically done via stochastic optimization of the ELBO. Initializing this covariance with FisherMALA's $I^{-1}$ estimate (learned from a short MCMC run) could provide a much better starting point than the identity matrix, potentially accelerating VI convergence and improving the quality of the variational approximation.

When to Prefer This Method

The paper positions FisherMALA against three named alternatives — AdaMALA (state-based adaptation), mMALA (position-dependent Riemannian preconditioning), and simple MALA / HMC without learned preconditioning — with clear empirical tradeoffs that support a conditional decision rule:

• Prefer FisherMALA over AdaMALA when the target is non-Gaussian and $d \geq 7$. On the Pima Indians dataset ($d = 7$), FisherMALA's min ESS advantage is already ~880× (5629 vs. 6.4). On Gaussian targets where $I^{-1} = \Sigma$, FisherMALA still learns the same matrix faster (Figures 1b, 2c), providing efficiency gains even when both methods target the same optimum. AdaMALA may be competitive only at very low dimensions ($d \leq 3$, as on Ripley) or when the target is nearly Gaussian and the state-based estimator has adequate data.

• Prefer FisherMALA over mMALA when the target is log-concave and $d > 25$, OR when per-iteration computational cost matters. On Caravan ($d = 87$), FisherMALA achieves ~10× better min ESS; on MNIST ($d = 785$), ~85× better. mMALA may be preferable on low-dimensional problems with strong local curvature variation (though the paper does not demonstrate such a case) or when the exact Hessian is cheap to compute and the $O(d^3)$ cost is acceptable.

• Prefer FisherMALA over simple MALA / HMC when the target exhibits substantial anisotropy — quantified as the condition number of the posterior covariance exceeding ~10. The paper's inhomogeneous Gaussian target (standard deviations spanning two orders of magnitude) causes simple MALA min ESS to collapse to ~3 vs. FisherMALA's ~1500. Even HMC (generally stronger than MALA) achieves min ESS ~42 vs. FisherMALA's ~1500 on this target. When the target is approximately isotropic or when only rough posterior exploration is needed, the additional complexity of FisherMALA's adaptation may not be justified.

• Prefer standard HMC over FisherMALA when the target is high-dimensional but well-conditioned (e.g., isotropic or nearly isotropic), since HMC's Hamiltonian dynamics provide superior exploration in the absence of anisotropy, and FisherMALA's $O(d^2)$ adaptation overhead provides no benefit.

These decision boundaries are empirically grounded in the paper's specific experiments (Gaussian targets at $d = 100$, Bayesian logistic regression from $d = 3$ to $d = 785$) and should be treated as provisional until validated on a broader range of target distributions, particularly non-log-concave and multimodal targets where FisherMALA's behavior is unknown.