Unadjusted Langevin algorithm, irreversible perturbations, expected squared jump distance, Markov chain Monte Carlo

1 Introduction↩︎

In modern Bayesian inference and scientific computing, a central challenge is to compute expectations under a target distribution \(\pi\) known only up to normalization [1]–[4]. Markov chain Monte Carlo (MCMC) methods address this problem by constructing a Markov process with stationary distribution \(\pi\), enabling estimation via ergodic averages without explicit integration (e.g., quadrature) in high dimensions [5]–[9]. This paper focuses on the unadjusted Langevin algorithm (ULA), a gradient-based MCMC method derived from discretizing overdamped Langevin dynamics. ULA is widely used for its simplicity and scalability, and its mixing rate and discretization bias have been studied extensively [10]–[13].

Despite its simplicity, ULA faces two key limitations: slow mixing on challenging targets and discretization bias that grows with step size [11], [14]. To address these issues, two complementary strategies have been developed to accelerate convergence by enlarging spectral gaps and reducing asymptotic variance. Reversible perturbations reshape the dynamics using local or global geometric information [15]–[24], while irreversible perturbations introduce additional drift to enhance exploration of the state space [25]–[35]. Both strategies are well studied in continuous time: reversible perturbations, also understood as preconditioning, are supported by rigorous optimality principles, and irreversible perturbations have been shown to accelerate a wide range of diffusion-based samplers. For irreversible perturbations, however, most analyses remain qualitative and/or confined to the continuous-time limit. [36] provide a principled characterization of optimal perturbations in the continuous-time Gaussian setting. [32] propose a way of constructing geometry-informed irreversible perturbations for more general non-Gaussian targets, but do not invoke a particular notion of optimality. It is also important to note that irreversible perturbations can introduce stiff, multiscale dynamics, which in turn degrade the accuracy of standard time discretization schemes, as observed in numerical studies of irreversible Langevin samplers [37]. A principled framework for designing irreversible perturbations that simultaneously accounts for discretization error and extends beyond Gaussian target distributions \(\pi\) remains largely unrealized.

This paper develops a principled framework for designing state-independent irreversible perturbations of ULA, explicitly balancing two competing goals: accelerating mixing and controlling discretization bias. Our framework first imposes a constraint that enforces a mixing-improvement criterion, motivated by spectral-gap optimality in the Gaussian setting but extended to the non-Gaussian case via the Fisher information matrix. Then we optimize within this admissible class of perturbations, minimizing a discretization error proxy based on a weighted expected squared jump distance (ESJD\(^W\)). The solution set of this optimization problem can be characterized explicitly, and yields an irreversible perturbation that achieves the best trade-off between these two objectives.

We demonstrate our method on several representative target distributions: anisotropic Gaussians, Gaussian mixtures, and the posterior distributions of Bayesian logistic regression and Bayesian independent component analysis problems. In these experiments, ULA equipped with our optimal irreversible perturbation achieves a favorable trade-off between mixing efficiency (estimator variance) and estimator bias, and lower mean-squared estimation errors overall. Additionally, it remains stable over a broad range of step sizes, in regimes where alternative constructions of irreversible perturbations can diverge or stagnate.

The rest of the paper is organized as follows: 2 reviews irreversible perturbations of the overdamped Langevin dynamics. 3 presents our method for constructing optimal irreversible perturbations, informed by the Fisher information of the target distribution, together with a practical approximation algorithm based on online Fisher information matrix estimation. 4 reports numerical experiments comparing our method to alternative designs of irreversible perturbation. 5 provides a discussion of future research directions.

The appendix contains derivations of formulas, proofs of the theoretical results of the paper, and additional numerical examples. In particular, Appendix 6 discusses reversible perturbations in the context of choosing a preconditioner that minimizes the expected square jump distance (ESJD), in the spirit of [23], which has partially motivated the approach in this paper. Appendix 7 contains explicit derivations of various formulas used throughout the paper, including the mean and covariance of Gaussian ULA iterates and the discretization error originating from the Itô–Taylor expansion. Appendix 8 contains rigorous and complete proofs of the theoretical results of the paper. Appendix 9 discusses aspects of Fisher information-based preconditioning. Appendix 10 contains additional experiments and numerical details. Appendix 11 offers an empirical and target-dependent comparison of irreversible and reversible perturbations. Code for our algorithms and experiments is available at https://github.com/qianyu-zhu/IrrevULA.

Acknowledgments↩︎

JZ and YM acknowledge support from the US Department of Energy, Office of Advanced Scientific Computing Research, under grants DE-SC0023188 and DE-SC0026245. KS was partially supported by the National Science Foundation (NSF) under grant DMS-2311500. BZ was supported by AFOSR grant FA9550-21-1-0354.

2 Background on sampling from equilibrium↩︎

Let \(\pi\propto e^{-V}\) denote a target density on \(\Omega \subset \mathbb{R}^d\). The overdamped Langevin diffusion targeting \(\pi\) is \[\label{eq:overdamped95LD} \mathrm{d} X_t = -\nabla V(X_t) \, \mathrm{d}t + \sqrt{2} \, \mathrm{d}W_t,\tag{1}\] where \(W_t\) is a standard Brownian motion and the initial condition of the process is \(X_0\). Let \(\pi_t\) denote the law of \(X_t\). (When \(\pi_t\) has a density with respect to the Lebesgue measure on \(\mathbb{R}^d\), we use the same notation for both measure and density.) Under mild regularity conditions on the potential \(V\) (e.g., it is smooth and confining), the process admits a unique invariant measure \(\pi\) and \(\pi_t\) converges to \(\pi\) weakly [38]. The ergodic theorem then ensures \[\mathop{\mathrm{\mathbb{E}}}_\pi[f(X)] = \int_{\mathbb{R}^d}f(x)\pi(x)\text{d}x = \lim_{T\rightarrow \infty}\frac{1}{T}\int_0^T f(X(t))\text{d}t.\] For a detailed treatment, see [38].

As discussed in the introduction, we focus on irreversible perturbations, which introduce an additional drift term \(\phi: \mathbb{R}^d \to \mathbb{R}^d\): \[\label{eq:irr-LD} \mathrm{d}X_t \;=\; \big[\,\nabla\log\pi(X_t)\;+\;\phi(X_t)\,\big]\mathrm{d}t \;+\; \sqrt{2}\,\mathrm{d}W_t.\tag{2}\] The process leaves \(\pi\) invariant whenever \(\phi\) is \(\pi\)–divergence free, i.e.\(\nabla\cdot\big(\phi(x)\,\pi(x)\big)=0\). A natural choice satisfying this constraint is \(\phi(x)=J\,\nabla\log\pi(x)\) with a constant skew-symmetric matrix \(J^\top=-J\), since \[\nabla\!\cdot\!\big(\phi\,\pi\big)=\nabla\!\cdot\!\big(J\nabla\pi\big)=\mathrm{tr}\big(J\,\nabla^2\pi\big)=0.\] We call this the state-independent case; see Remark [rmk:state-dependent-perturbation] for the general case.

Such perturbations break detailed balance by inducing a rotational vector field which enhances state-space exploration [26], [31]. Prior work has shown that irreversible drift accelerates convergence to equilibrium, with optimal rates and corresponding optimality criteria characterized for Gaussian targets [25], [29]. Building on this line of work, [36] provide an explicit construction of optimal irreversible perturbations achieving these rates. Behavior of the generator under strongly scaled nonreversible drift has also been analyzed [39]–[41], and subsequent work established that irreversible perturbations for overdamped Langevin dynamics reduce the asymptotic variance of the ergodic estimators [33]–[35], [42].

In this work, we focus on state-independent case and identify a class of perturbations \(J\) that jointly reduces mixing time and discretization error. For brevity, we refer to \(J\) itself as the irreversible perturbation.

2.1 Optimal irreversible perturbation for Gaussian targets in continuous time↩︎

For Gaussian targets, the impact of irreversible perturbations can be quantified exactly through the spectral gap, enabling explicit characterization of optimal perturbations. Consider \(\pi \propto \exp(-x^\top S x / 2)\) with irreversible perturbations \(\phi(x) = JSx\), \(J = -J^\top\) [25], [36]. The infinitesimal generator of the corresponding Langevin dynamics is \[\mathcal{L} = -(I + J) S x \cdot \nabla + \Delta \, .\] The spectrum of \(\mathcal{L}\) in \(L^p(\pi)\), \(p>1\), consists of negative integer linear combinations of the eigenvalues of \((I + J) S\) [44]: \[\sigma(\mathcal{L}) = \left\{ -\sum_{j=1}^r n_j \sigma_j \; : \; \forall\, n_j \in \mathbb{N}_0 \right\},\] where \(\{\sigma_j\}_{j=1}^r\) are the (generally complex) eigenvalues of \((I + J) S\), all satisfying \(\Re(\sigma_j)>0\) (where \(\Re\) denotes the real part). The spectral gap is thus determined by the smallest among \(\{\Re(\sigma_j)\}_{j=1}^r\). Since \(J\) is skew-symmetric, \(\mathrm{Tr}((I+J)S) = \mathrm{Tr}(S)\), so the sum of the real parts of the eigenvalues \(\{\sigma_j\}_{j=1}^r\) is fixed. Thus, maximizing the spectral gap corresponds to maximizing the smallest real part under a fixed sum constraint, which is achieved when all \(\Re(\sigma_j)\) are equal [36]: \[\label{eq:J95spec} \mathcal{J}_{\mathrm{s}}(S) := \left\{J\in\mathbb{R}^{d\times d}: J^\top=-J,\; \min_i \Re \big( \sigma_i( (I + J) S ) \big) = \frac{\mathrm{Tr}(S)}{d}\right\}.\tag{3}\] where the subscript “s” denotes spectral. When \(S\) is clear, we write \(\mathcal{J}_{\mathrm{s}}\). The maximizer of the spectral gap in continuous time is generally not unique; [36] further proposed a constructive algorithm that yields admissible maximizers within this class. However, this analysis is confined to continuous time and does not account for discretization error, a critical consideration in practice. We address this gap in 3.

3 Design of optimal position-independent irreversible perturbations↩︎

In practice, the continuous-time Langevin dynamics are approximated by a discrete-time scheme, typically Euler–Maruyama [45]: \[\label{eq:ULA} X_{t} = X_{t-1} + h(I+J)\nabla\log\pi(X_{t-1})+ \sqrt{2h} \, Z_{t-1},\qquad t = 1,2,\dots\tag{4}\] where \(h>0\) is the step size and \(\{Z_t\}_{t\geq 0}\) are i.i.d.standard Gaussian. This first-order discretization, also referred to as the unadjusted Langevin algorithm (ULA), admits an invariant distribution \(\pi_h\) that generally differs from \(\pi\). While a Metropolis–Hastings correction can recover \(\pi\), the asymmetry induced by irreversible perturbations lowers acceptance rates and adds computational cost, especially in high dimensions [46], [47]. We therefore focus on unadjusted dynamics and design \(J\) to balance mixing acceleration against discretization bias.

Our approach is one of constrained optimization. Among choices of \(J\) that satisfy the Gaussian spectral optimality defined by [36], we will minimize a discretization-error proxy based on the weighted expected squared jump distance. We begin with the Gaussian setting, where spectral gap can be calculated exactly, before extending these principles to general targets.

For a Gaussian target with mean \(\mu\) and covariance matrix \(S^{-1}\), the ULA update in 4 yields \[\label{eq:ULA-gaussian} X_{t} = X_{t-1} - h (I+J)S(X_{t-1}-\mu) + \sqrt{2h}\,Z_{t-1}.\tag{5}\] With \(X_0 \sim \mathcal{N}(\mu_0, \tau^2 I)\), the iterates remain Gaussian with mean and covariance \[\begin{align} \mu_t(h) = \mu+B^t(\mu_0-\mu),\;\; \Sigma_t(h) = \tau^2B^t(B^\top)^t + 2h\sum_{i=0}^{t-1}B^i(B^\top)^i, \;\; B:=I-h(I+J)S. \end{align}\] See 7.1 for a detailed derivation. As \(t\rightarrow\infty\), \(\pi_t\) converges to the Gaussian distribution \(\pi_h\), whose moments as functions of \(h\) are \[\begin{align} \label{eq:momentswithbias} \mu(h) = \mu,\qquad \Sigma(h) = S^{-1}+h\Sigma^{(1)}+ \mathcal{O}(h^2), \end{align}\tag{6}\] where \(\Sigma^{(1)}\) satisfies \[\begin{align} A\Sigma^{(1)}+\Sigma^{(1)}A^\top=(I+J) S (I+J)^\top, \;\text{with} \;A := (I+J)S. \end{align}\] While the mean is unbiased, the covariance incurs bias driven by \((I+J)S(I+J)^\top\), reflecting the discrepancy between the invariant distribution \(\pi_h\) of the discretized dynamics and the target distribution \(\pi\).

Moment estimation error provides useful intuition for how discretization affects downstream statistics, but it is not a sufficiently general design objective. Directly optimizing the bias of a particular moment estimator is both observable-dependent and less principled as a general design criterion. Since our goal is to choose perturbations that improve the overall performance of the discretized sampler, we instead introduce tractable proxies for discretization error. These proxies capture the leading-order dependence on \(J\) while remaining general enough to support analysis and optimization.

3.1 Discretization error proxies↩︎

We introduce three proxies for the discretization error, evaluated in closed form for the ULA 5 with \(X_t\sim\pi=\mathcal{N}(\mu, S^{-1})\).

3.1.0.1 (i) ESJD.

The expected squared jump distance, used as a discretization error criterion in [23], measures the average squared displacement per iteration: \[\label{eq:ESJD-def} \mathrm{ESJD}(X_t, h)\;:=\;\mathbb{E}\big[\|X_{t+h}-X_t\|^2\big].\tag{7}\] Substituting the Gaussian ULA step 5 at stationarity, \[\begin{align} \mathrm{ESJD}(X_t,h) &= \mathrm{Tr}\!\left(h^2 A\, S^{-1}A^\top\right)+2h d = -h^2\mathrm{Tr}\big[ J S J \big] +g_{1}(S,h),\qquad A:=(I+J)S, \end{align}\] where \(g_{1}(S,h):=h^2\mathrm{Tr}[S]+2h d\) collects \(J\)-independent terms and we used \(\text{Tr}[JS+SJ^\top]=0\) by skew-symmetry. Minimizing ESJD over \(J\) thus reduces to minimizing \(\mathcal{E}_1(J) :=-\mathrm{Tr}\big[ J S J \big]\).

This proxy has a direct connection to the moment estimation error in 6 : the \(J\)-dependent part of ESJD exactly equals the \(J\)-dependent part of the forcing term driving the leading covariance bias. In particular, \(\mathrm{Tr}\bigl( (I+J) S (I+J)^\top \bigr ) = \mathrm{Tr}(S) + \mathrm{Tr}(JSJ^\top)\). Larger irreversible perturbations produce larger jumps and, simultaneously, larger bias.

3.1.0.2 (ii) ESJD\(^W\).

For anisotropic targets, a large step along a high-curvature direction inflates bias more than the same step along a flat direction. A natural remedy is to measure jumps in coordinates that isotropize the target. This generally motivates a weighted ESJD with weight \(W \succ 0\): \[\label{eq:w-ESJD-def} \mathrm{ESJD}^W(X_t, h) \;:=\;\mathbb{E}\big[\|X_{t+1}-X_t\|_{W}^2\big].\tag{8}\] where we use the notation \(\|a \|_W^2 := a^\top W a\). At stationarity, \[\begin{align} \mathrm{ESJD}^W(X_t,h) &= h^2\mathrm{Tr}\!\left(W(I+J)S(I+J)^\top\right)+2h\mathrm{Tr}[W]. \end{align}\] Taking \(W=S\) makes the cross terms vanish. In particular, we get that \[\begin{align} \mathrm{ESJD}^S(X_t,h) &= - h^2\mathrm{Tr}(JSJS)+g_{2}(S,h), \end{align}\] with \(g_{2}(S,h)=-h^2\mathrm{Tr}(SS)+2h\mathrm{Tr}[S]\). This calculation gives the proxy \(\mathcal{E}_2(J) := -\mathrm{Tr}(JSJS)\). Compared with \(\mathcal{E}_1\), this proxy favors larger effective steps along directions where the target geometry is flatter.

3.1.0.3 (iii) Itô–Taylor remainder.

A complementary viewpoint comes from the numerical analysis of SDEs: the local truncation error of the Euler–Maruyama scheme directly quantifies the per-step approximation quality. The Itô–Taylor expansion at \(X_t\) (the stochastic analogue of a Taylor series) gives the leading error term as \[X_{t+h} - \tilde{X}_{t+h} = -\sqrt{2}(I+J)S\,Z_{t,h}+ \mathcal{O}(h^2),\] where \(\tilde{X}_{t+h}\) denotes the EM approximation of \(X_{t+h}\) and \(Z_{t,h}\sim\mathcal{N}(0,\tfrac{h^3}{3}I)\). Its expected squared norm \(\mathbb{E}\big[\|X_{t+h}-\tilde{X}_{t+h}\|^2\big]\) is proportional to \(\mathrm{Tr}[(I+J)SS(I+J)^\top]\), yielding \(\mathcal{E}_3(J) := -\mathrm{Tr}(J S S J)\) after isolating \(J\)-dependent terms; see 7.2 for the full derivation.

To summarize, we proposed the following three discretization error proxies: \[\mathcal{E}_1(J) = -\mathrm{Tr}(J S J), \qquad \mathcal{E}_2(J) =-\mathrm{Tr}(J S J S),\qquad \mathcal{E}_3(J) = -\mathrm{Tr}(J S S J).\] The next section compares these functions and examines how well they track the mean-square estimation error that ultimately determines sampling performance.

3.2 Comparison of proxies and relationship with MSE↩︎

As a numerical illustration, we fix a Gaussian target with precision \(S=\operatorname{diag}(1,1/2,1/5)\), and randomly sample \(J\in\mathcal{J}_\mathrm{s}\) using the constructive algorithm of [36]. 1 and 2 reveal the following key findings.

• The histograms in the top row of 1 show that none of \(\mathcal{E}_{\{1,2,3\}}(J)\) is bounded over \(\mathcal{J}_\mathrm{s}\). Thus, spectral optimality alone does not control discretization error: arbitrary choices of spectrally optimal perturbation can yield arbitrarily large discretization error. This motivates a second-stage optimization of \(\mathcal{E}_i(J)\) within \(\mathcal{J}_\mathrm{s}\).

• The pairwise scatter plots in the bottom row of 1 show strong positive correlation among all three criteria. This suggests that optimizing a single representative proxy is sufficient in practice. Proposition [prop:proxy-equivalence] below formalizes this equivalence.

• 2 connects the proxies to estimator accuracy. In this example, the proxies show little correlation with the MSE of the estimator of the first moment \(\mathbb{E}[\mathrm{sum}(X^{(i)})]\), but strong positive correlation with the MSE of the estimator of the second moment \(\mathbb{E}[\|X\|_2^2]\). This matches the discretization-bias analysis in the opening of 3: mean estimation is asymptotically unbiased even after discretization, whereas covariance-related quantities are affected by discretization.

The next proposition establishes equivalence of \(\mathcal{E}_{ \{1, 2, 3\} }\) up to constants depending only on the spectrum of \(S\). See 8.1 for the proof.

Together, these findings justify optimizing a single discretization-error proxy within \(\mathcal{J}_\mathrm{s}\). We will see in 3.4 that only \(\mathcal{E}_2\) admits a tractable closed-form optimization under the spectral constraint.

3.3 Illustration: two-dimensional Gaussian↩︎

Before presenting our general method for selecting \(J\in\mathcal{J}_\mathrm{s}\) while controlling discretization error, we first illustrate the key trade-off in a two-dimensional setting with a Gaussian target, where every skew-symmetric \(J\in\mathbb{R}^{2\times 2}\) is of the form \[J = \begin{pmatrix} 0 & a \\ -a & 0 \end{pmatrix} \qquad\text{ with }\quad a \in \mathbb{R}.\] For the precision matrix \(S = \text{diag}(1,\varepsilon)\), the spectral optimality condition in [36] is attained by every \(a\) satisfying \(|a|\geq (1-\varepsilon)/2\sqrt{\varepsilon}\). Thus, continuous-time spectral optimality alone does not identify a unique perturbation. Once discretization error is taken into account, however, the ambiguity is removed: minimizing \(\mathcal{E}_2\) over this spectrally optimal class selects the two boundary values \(a^* =\pm (1-\varepsilon)/2\sqrt{\varepsilon}\).

3 visualizes the first part of this trade-off. The blue shaded region denotes the spectrally optimal class: once \(|a|\) exceeds the threshold, the spectral gap has reached its maximal value, so all larger perturbations are equivalent from the continuous-time spectral perspective. The discretization proxy, however, continues to increase with \(|a|\). Thus, minimizing \(\mathrm{ESJD}^{S}\), equivalently \(\mathcal{E}_2\), selects the smallest perturbation that is still spectrally optimal, namely the bifurcation point where the spectral constraint first becomes active. In this two-dimensional Gaussian example, this same boundary value also minimizes \(\mathcal{E}_1\) and \(\mathcal{E}_3\), although this coincidence is special to this setting.

4 then shows why this distinction matters for finite-step estimation. For linear statistics, such as \(\mathbb{E}[\operatorname{sum}(X)]\), larger \(|a|\) primarily improves mixing and therefore reduces variance and hence MSE. For nonlinear statistics, however, the behavior is no longer monotone: increasing \(|a|\) reduces variance at first, but excessive perturbation strength introduces discretization bias. Thus, perturbations that are equivalent in continuous time can lead to different finite-sample accuracy after discretization.

This example identifies the basic selection problem. Spectral optimality determines a class of accelerating perturbations \(\mathcal{J}_s\), but discretization error determines which member of that class is preferable in practice. In two dimensions this reduces to choosing the boundary point \(a^\ast\); in higher dimensions, the spectrally optimal set has many more degrees of freedom because the space of skew-symmetric matrices has dimension \(d(d-1)/2\). This motivates the general optimization framework developed next.

3.4 Optimal irreversible perturbations as solution to constrained optimization↩︎

We now formalize the design of optimal irreversible perturbations as a constrained optimization problem that minimizes discretization error while maximizing spectral gap. Although all three proxies correlate positively with MSE and are equivalent up to spectral constants, as shown in the preceding subsections, only \(\mathcal{E}_2\) leads to an explicit characterization of the optimal solutions.

We write \(\gamma_\star=\mathrm{Tr}(S)/d\) and minimize \(\mathcal{E}_2(J)\) subject to the spectral gap constraint of [36]: \[\label{eq:P1} \begin{align} \min_{J:J = -J^\top} \quad &\mathcal{E}_2(J) := \mathrm{Tr}(-JSJS) \\ \text{s.t.} \qquad &\Re\!\left[\sigma_i\big((I+J)S\big)\right] = \gamma_\star, \;\; \forall i. \end{align}\qquad{(1)}\] The constraint enforces spectral optimality: all eigenvalues of \((I+J)S\) share the same real part, maximizing the spectral gap for the continuous-time process. The objective penalizes large jumps in high-curvature directions, controlling discretization bias. Together, they identify perturbations that accelerate mixing without sacrificing numerical accuracy.

3.4.0.1 Equivalent spectral coordinate formulation.

Now define \[\tilde{J}=S^{1/2}JS^{1/2}, \qquad \tilde{B}_J=S+\tilde{J} =S^{1/2}(I+J)S\,S^{-1/2}.\] Then \(\tilde{B}_J\) is similar to \((I+J)S\), and \[\mathrm{Tr}(-JSJS)=\|\tilde{J}\|_F^2.\] For any orthonormal matrix \(\Psi=[\psi_1,\ldots,\psi_d]\), write \[B_\Psi=\Psi^\top S\Psi, \qquad M=\Psi^\top \tilde{J}\Psi .\] Orthogonal invariance gives \(\|\tilde{J}\|_F=\|M\|_F\), and \(\Psi^\top \tilde{B}_J\Psi=B_\Psi+M\). Hence ?? is equivalent to the following problem in different coordinates: \[\label{eq:P2} \begin{align} \min_{\Psi,M}\qquad &\|M\|_F^2\\ \text{s.t.}\qquad &\Psi^\top\Psi=I,\quad M^\top=-M,\\ &\Re\!\left[\sigma_i(B_\Psi+M)\right]=\gamma_\star,\quad i=1,\ldots,d . \end{align}\qquad{(2)}\] Given a feasible pair \((\Psi,M)\), the corresponding perturbation is \[J=S^{-1/2}\Psi M\Psi^\top S^{-1/2}.\] Conversely, every feasible \(J\) in ?? gives feasible coordinates \((\Psi,M)\) in ?? by setting \(M=\Psi^\top S^{1/2}JS^{1/2}\Psi\).

3.4.0.2 Optimal value and solution set.

The form of ?? makes the attained optimum explicit.

The proof is in 8.2.

We define the solution set to ?? as \[\mathcal{J}_\mathrm{se} := \operatorname*{argmin}_{J\in\mathcal{J}_\mathrm{s}}\;\mathrm{Tr}(-JSJS).\] [algo:algorithm] provides a constructive procedure to obtain an element of \(\mathcal{J}_{\mathrm{se}}\). It is inspired by the constructive approach of [36], whose goal is to build an orthonormal basis \(\Psi=(\psi_1, \dots, \psi_d)\) satisfying \(\psi_k^\top S\psi_k = \gamma_\star\). This condition ensures that the resulting perturbation can satisfy the optimal spectral constraint. However, this construction is aimed at spectral optimality alone, whereas our goal is to select, among spectrally optimal perturbations, one that also optimizes the ESJD\(^W\) discretization-error proxy. Our algorithm thus modifies the previous constructive procedure in two respects.

First, instead of repeatedly selecting directions with maximal and minimal values of \((\psi_i,S\psi_i)\), we use a threshold-based selection rule. This avoids discontinuous changes in the selected directions when \(S\) is slightly perturbed, which is useful in adaptive settings where the precision matrix is estimated and updated online; see Section 3.5.1.2.

Second, after constructing such a basis \(\Psi\), spectral optimality still leaves degrees of freedom in the skew-symmetric perturbation. These degrees of freedom do not affect the spectral constraint, but they do affect the discretization error proxy. We resolve this ambiguity by choosing the triangular perturbation in [thm:P1-P2-relation], namely \[M_{jk}=\mathrm{sign}(k-j)\psi_j^\top S\psi_k,\quad j\neq k, \qquad M_{kk}=0.\] This choice preserves spectral optimality while attaining the minimum of \(\mathcal{E}_2(J)\).

Returning to 1, the orange points show perturbations \(J\in\mathcal{J}_{\mathrm{se}}\) generated by [algo:algorithm]. These perturbations, which minimize \(\mathcal{E}_2\) by construction, simultaneously achieve low values of \(\mathcal{E}_1\) and \(\mathcal{E}_3\), confirming the equivalence of [prop:proxy-equivalence]: for well-conditioned \(S\), optimizing any single proxy effectively controls the others.

3.5 Extension to non-Gaussian targets↩︎

The constrained optimization problem ?? is formulated for Gaussian targets, where the precision matrix \(S\) fully determines the geometry of the distribution. For a general target \(\pi\), we can still ask what the best position-independent irreversible perturbation should be. A natural way to characterize the target distribution for this purpose is through its Fisher information matrix: \[\mathcal{F} \;=\; \mathbb{E}_{\pi}\!\left[\nabla\log\pi(X)\,\nabla\log\pi(X)^\top\right].\label{Eq:FisherGeneral}\tag{9}\] This choice is principled in two ways: first, \(\mathcal{F}\) captures the average local curvature of \(\log\pi\), and for a Gaussian with precision \(S\) it recovers \(\mathcal{F}=S\) exactly. Second, the weighted ESJD for general \(\pi\) is a function of \(\mathcal{F}\): \[\begin{align} \mathrm{ESJD}^W(X_t,h) :&= \mathrm{Tr} \, \mathbb{E}\big[ W^{1/2}(X_{t+h} - X_t)(X_{t+h} - X_t)^\top W^{1/2}\big] \\ &=h^2\mathrm{Tr}\!\left(W(I+J)\mathcal{F}(I+J)^\top\right)+2h\mathrm{Tr}[W]. \end{align}\] Setting \(W = \mathcal{F}\) makes the cross terms vanish. Writing \(\gamma_{\mathcal{F}}=\mathrm{Tr}(\mathcal{F})/d\), minimizing this \(\mathcal{F}\)-weighted ESJD over \(J \in \mathcal{J}_s(\mathcal{F})\) yields the following optimization problem: \[\label{eq:P3} \begin{align} &\min_{J = -J^\top} \qquad-\mathrm{Tr}\;\left( J \mathcal{F} J \mathcal{F}\right)\\ &\text{s.t.}\qquad \qquad \Re[\sigma_i\left((I+J)\mathcal{F}\right)] = \gamma_{\mathcal{F}},\qquad\forall i. \end{align}\qquad{(3)}\] This is precisely the Gaussian problem ?? with \(S\) replaced by \(\mathcal{F}\), so the same closed-form solution and [algo:algorithm] apply directly. We note that the Fisher information matrix 9 is also used by [23] to construct optimal position-independent reversible perturbations; here, we use the same target summary for the irreversible setting.

3.5.1 When does the Fisher-based approach work?↩︎

Since \(\mathcal{F}\) averages local geometry over the entire state space, it is most informative when this geometry is sufficiently consistent. We now analyze two regimes where design based on the global Fisher information matrix \(\mathcal{F}\) is most effective and describe how this matrix can be estimated adaptively in practice.

3.5.1.1 Near-Gaussian targets.

We first consider target distributions \(\pi \propto \exp(-V)\) that are close to Gaussian, in the sense that \(V\) is a small perturbation of a quadratic form. Write \(\mathcal{F}_\pi\) for the Fisher information matrix 9 of such a \(\pi\). The natural question is whether an irreversible perturbation induced by \(J^* \in \mathcal{J}_\mathrm{s}(\mathcal{F}_\pi)\), i.e., designed from \(\mathcal{F}_\pi\), will still accelerate the convergence of Langevin dynamics for \(\pi\). The next result shows that, under a diagonalizability assumption, this spectral acceleration is indeed stable to perturbations.

The proof is provided in 8.5.

3.5.1.2 Adaptive estimation of the Fisher information matrix.

In practice, the Fisher information matrix \(\mathcal{F}\) is typically unknown and must be estimated during sampling. We employ a streaming estimator that updates the matrix online, using score observations: \[\label{eq:adapted-fisher} \begin{align} \text{Initialization:} \quad &\widehat{\mathcal{F}}_0 = I_d, \\ \text{Iteration:} \quad &\widehat{\mathcal{F}}_n = \widehat{\mathcal{F}}_{n-1} + \frac{1}{K+n}\left( s_n s_n^\top - \widehat{\mathcal{F}}_{n-1} \right), & n \geq 1, \end{align}\tag{10}\] where \(s_n = \nabla\log\pi(X_n)\) is the score at iteration \(n\) and \(K>0\) controls the transition from initialization to online updates. This estimator is consistent and numerically stable, as illustrated in the numerical experiments of 4.

3.5.1.3 Aligned local geometries.

Beyond the near-Gaussian regime, the Fisher information matrix can still be informative in cases where multiple modes share similar local geometry, even though obtaining formal guarantees in this setting is challenging. When the target distribution has multiple modes with similar shapes and orientations, the global Fisher information matrix aggregates their local curvature information in a coherent manner. 6 illustrates this phenomenon in two dimensions: when modes are misaligned (left), their local contributions cancel and \(\mathcal{F}\) loses directional information; when modes are geometrically aligned (middle and right), \(\mathcal{F}\) reflects their shared local structure.

Together, these results suggest that the \(\mathcal{F}\)-based optimization is most effective for targets that are either (i) approximately Gaussian, or (ii) multi-modal with geometrically aligned modes. We will verify these guidelines empirically in 4.

4 Numerical experiments↩︎

In this section, we evaluate our method against several alternative irreversible perturbations on four benchmark problems: a non-isotropic Gaussian (4.1), a Gaussian mixture (4.2), Bayesian logistic regression (4.3), and a Bayesian formulation of independent component analysis (4.4). For each problem, we compare MCMC estimates of representative observables against ground truth values obtained analytically or from high-accuracy MALA simulations.

4.0.0.1 Methods.

We run ULA with the six types of irreversible perturbations summarized in 1. These comprise the unperturbed case (\(J=0\)) and five other choices of \(J\).

In this table, \(C\) is the Frobenius norm of \(J\) from spec-E, averaged over randomness in Algorithm [algo:algorithm]. \(\mathcal{F}\) is given by 9 . Note that neither spec nor spec-E is unique: both depend on the random initial basis supplied to the constructive algorithm, so each run produces a different \(J\) from the solution sets \(\mathcal{J}_{\mathrm{s}}\) or \(\mathcal{J}_{\mathrm{se}}\).

4.0.0.2 Evaluation metrics.

For an observable \(\varphi:\mathbb{R}^d \to \mathbb{R}\), the ULA estimator of \(\mathbb{E}_\pi [\varphi(X)]\) based on \(K\) steps is \(\bar{\varphi}_K \;=\; \frac{1}{K} \sum_{k=0}^{K-1} \varphi(X_k)\). We evaluate the performance of the estimator using bias and variance estimated from \(M\) independent chains: \[\begin{align} \widehat{\mathrm{Bias}}(\bar{\varphi}_K) &= \frac{1}{M}\sum_{m=1}^{M}\bar{\varphi}_K^{(m)} - \mathbb{E}_\pi[\varphi(X)], \quad\quad \widehat{\mathrm{Var}}(\bar{\varphi}_K) = \frac{1}{M}\sum_{m=1}^{M}\left(\bar{\varphi}_K^{(m)}\right)^2 - \left(\frac{1}{M}\sum_{m=1}^{M}\bar{\varphi}_K^{(m)}\right)^2, \end{align}\] where \(\bar{\varphi}_K^{(m)}\) denotes the estimate from the \(m^\mathrm{th}\) chain. The MSE of the estimator is then simply \(\widehat{\mathrm{Bias}}(\bar{\varphi}_K)^2 + \widehat{\mathrm{Var}}(\bar{\varphi}_K)\). The values of \(K\) and of the ULA step size \(h\) will be varied systematically in the experiments below.

The chains for the first two experiments are initialized randomly from the target distribution; those for the last two experiments are initialized at \(\mathbf{0}\) and at a random diagonal sign matrix \(\text{diag}(\pm1)\), with burn-in intervals of \(0\) and \(5000\) steps, respectively, discarded before computing any \(\bar{\varphi}_K\).

In addition, we evaluate numerical stability by examining the probability that an estimator fails to converge, in the sense that it exceeds a prescribed threshold chosen to be significantly larger than the true value of \(\mathbb{E}_\pi[\varphi(X)]\). In plots of MSE versus step size \(h\), we exclude points for which the averaged estimates exceed this threshold, as this provides a sufficient (though not necessary) indication of divergence. Across the various experiments, we find that spec, rand-M, and rand-L are substantially more numerically unstable, a phenomenon that is particularly evident in Experiment 2. Accordingly, we classify these methods as “discretization-sensitive,” in contrast to the remaining “discretization-robust” methods.

4.1 Non-isotropic Gaussian↩︎

We begin with a zero-mean Gaussian target with covariance \(S^{-1} \;=\;\text{diag}(1,4,16,64)\), providing a controlled setting where \(\mathcal{F} = S\) is known exactly. Let the components of the state vector \(X\) be denoted as \(X = (X^{(1)}, X^{(2)}, X^{(3)}, X^{(4)})\).

4.1.0.1 Experiment 1: MSE versus step size.

7 and 8 compare the MSE, bias, and variance of estimators of the expected \(\ell_1\) norm \(\|X\|_1=\sum_{i=1}^4 |X^{(i)}|\) and the tail probability \(P(X^{(4)}>16)\) , across step sizes \(h \in [0.02, 0.4]\), with either fixed simulation time (\(T = hK =10\,000\)) or fixed computational budget (\(K = 10^5\)). We also test our method spec-E with adaptive estimation of the Fisher information matrix as described in 3.5.1.2, which is labeled as adaptive. Results for additional observables \(\phi\) are reported in 10.1.

The methods split into two groups. The discretization-unaware methods spec, rand-M, and rand-L, which do not involve \(\text{ESJD}^W\) minimization, diverge for \(h>0.05\) (red/orange curves in both figures). Other methods, namely rand-S, spec-E, and unperturbed, remain well-behaved across all \(h\) (blue curves) with spec-E converging the most quickly, achieving MSE roughly an order of magnitude lower than the other stable methods. Note that the \(J\) of rand-M has exactly the Frobenius norm as that of spec-E, showing that stability is not only a matter of the scale of the perturbation. Also notably, adaptive nearly matches the performance of the latter despite having no access to the true Fisher information matrix, estimating it on-the-fly from samples alone. The fact that the MSE of adaptive remains only marginally larger than that of spec-E demonstrates that our method is robust to Fisher information matrix misspecification—an issue that arises in practice when the target distribution is not analytically tractable.

The fixed-\(T\) and fixed-\(K\) experiments highlight complementary trade-offs. Under fixed simulation time (7), spec-E achieves lower variance for all \(h\) (right column) without sacrificing bias (middle column), reflecting its maximized spectral gap. Under fixed computational budget (8), increasing \(h\) reduces variance through longer simulation time but increases bias through discretization error. The resulting U-shaped MSE curves (left column) show that spec-E achieves the lowest MSE across a wide range of \(h\), demonstrating robustness to step size selection.

4.1.0.2 Experiment 2: Proxy comparison

9 compares the performance of spec-E with that of two alternative discretization-aware pertubations, obtained by solving ?? but with its objective replaced by \(\mathcal{E}_1\) or \(\mathcal{E}_3\). Since the latter optimization problems are not analytically tractable, we approximate their minimizers by Monte Carlo search over \(\mathcal{J}_s\). Using \(\mathcal{E}_2\) yields the most accurate estimation, for most values of \(h\). Improved performance over \(\mathcal{E}_{ \{ 1,3 \} }\) may be due to two complementary factors: first, minimizing \(\mathcal{E}_2\) admits an explicit and numerically stable construction; second, \(\mathcal{E}_2\) itself seems to better correlate with discretization error and thus provides a more effective criterion for selecting \(J\) in this example.

4.2 Mixture of Gaussians↩︎

We next consider a random vector \(X=(X^{(1)},X^{(2)},X^{(3)})\in\mathbb{R}^3\) whose target distribution is a mixture of five Gaussian components (21). The modes are separated along the first coordinate, so efficient exploration requires transitions between distant regions. At the same time, the local geometry within each mode is similar, and the global Fisher information matrix remains aligned with this shared local structure. This example therefore isolates a favorable multi-modal setting in which the Fisher information matrix provides an effective global summary of the target geometry, consistent with the discussion in 3.5.1.2. Further details on the problem setup are provided in 10.2.

4.2.0.1 Experiment 1: Stability and trajectory behavior.

10 reports the number of divergent chains among 128 independent runs, where divergence is declared if the empirical estimator of \(\mathop{\mathrm{\mathbb{E}}}_\pi [\|X\|_1]\) or \(\mathop{\mathrm{\mathbb{E}}}_\pi [\max_i(|X^{(i)}|)]\) exceeds \(50\) after \(T=10{,}000\) iterations, a threshold chosen to be more than twice the corresponding true value. The discretization-sensitive methods spec, rand-M, and rand-L become increasingly unstable as the step size grows: by \(h=0.3\), nearly all chains diverge for spec and rand-L. In contrast, the discretization-robust methods unperturbed, rand-S, and spec-E rarely diverge.

The trajectory visualizations in 12 further explain these stability trends. The unperturbed and rand-S/M chains often fail to visit all modes, while rand-L/spec can over-explore and produce unstable trajectories. The properly regularized spec-E achieves the desired balance: it improves inter-modal exploration while controlling discretization-induced instability.

4.2.0.2 Experiment 2: MSE versus step size.

11 reports MSE, bias, and variance for estimators of \(\mathbb{E}_\pi[ \max_i |X^{(i)}| ]\) and \(\Pr(X^{(1)}>20)\) under fixed simulation time \(T = 4000\). The pattern mirrors the Gaussian results: spec-E achieves the lowest MSE across a broad range of \(h\), while unperturbed/rand-S mix poorly at small \(h\) and rand-M/rand-L/spec accumulate bias at large \(h\). Additional statistics and fixed-budget experiments are in 10.2.

4.3 Bayesian logistic regression↩︎

Next we apply our method to Bayesian logistic regression on the german dataset from the benchmark collection used by [48] and compiled by [49] (\(d = 20\) features, \(N = 400\) observations). Given data \(\{ (Z_i, t_i) \}_{i=1}^N\) with \(Z = [z_1,\dots, z_N]\in\mathbb{R}^{d\times N}\), \(t_i \in \{0, 1\}\), the log posterior for weights \(X \in \mathbb{R}^d\) under a Gaussian prior \(X\sim\mathcal{N}(0,\alpha^{-1} I)\) is \[\log\pi(X\mid Z) = -\frac{\alpha}{2}\|X\|^2 + \sum_{i=1}^N \left[ t_i z_i^\top X - \log(1 + \exp(z_i^\top X)) \right] + \mathrm{const}.\] Define the logistic function \(\varphi(y) = (1+\exp(-y))^{-1}\). The gradient of the log-posterior is \[\nabla_{X} \log \pi(X \mid Z) = -\alpha X + \sum_{i=1}^N t_i z_i - \sum_{i=1}^N \varphi(z_i^\top X) \, z_i.\] The posterior of \(X\) is unimodal but anisotropic (see 24 in 10.3).

For the numerical experiments, we use stochastic gradient Langevin dynamics [50], [51] and subsample at a rate of \(n=10\) data points per likelihood computation, leading to the following drift \[\nabla_{X} \log \tilde{\pi}(X \mid Z) = -\alpha X + \frac{N}{n}\sum_{i=1}^n t_{\tau_i} z_{\tau_i} - \frac{N}{n}\sum_{i=1}^n \varphi(z_{\tau_i}^\top X) \, z_{\tau_i}\] where \(\{\tau_i\}_{i=1}^n\) are the indices of the subsampled data.

4.3.0.1 Experiment: MSE versus step size.

[fig:blr95statistics] reports MSE, bias, and variance of estimators of \(\mathop{\mathrm{\mathbb{E}}}_\pi[ |X^{(5)}|]\) and \(P(X^{(5)}>2)\) under fixed computational budget \(K=100,000\); estimator is estimated using \(M=128\) independent chains. Unlike the Gaussian targets considered above, this posterior has asymmetric tail behavior: along some directions, the density has flatter tails on one side and steeper decay on the other. Thus, while \(\mathcal{F}\) no longer exactly characterizes the target geometry, it still captures directions in which the log-density changes rapidly.

All methods remain numerically stable in this example, unlike in the Gaussian experiments where several unregularized methods diverged. Nevertheless, the unregularized spectral perturbation spec performs poorly. The random perturbations and unperturbed exhibit moderate MSE, limited by higher variance. spec-E retains much of the variance reduction from other irreversible perturbations while controlling bias in steep directions, and therefore achieves the lowest MSE across step sizes. Results for additional statistics and under fixed simulation time \(T\) are in 10.3.

4.4 Independent component analysis↩︎

The last application is a blind signal separation problem studied in [50], [52]. Given \(N\) observations of \(m\) mixed signals \(x_t = Ms_t\in\mathbb{R}^m\), where \(M\in\mathbb{R}^{m\times m}\) is an unknown mixing matrix and \(s_t\in\mathbb{R}^m\) contains the original source signals, the goal is to infer a de-mixing matrix \(W\in\mathbb{R}^{m\times m}\) such that \(Wx_t\) recovers the sources up to scaling and permutation.

We adopt a Bayesian formulation with Gaussian prior \(W_{ij} \stackrel{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0, \lambda^{-1})\). Following standard ICA practice, each recovered component \(y = [Wx]_i\) is assumed to follow a hyperbolic secant distribution \(p(y) = \frac{1}{4} \,\mathrm{sech}^2\!\left(\frac{y}{2}\right)\), yielding the dataset likelihood \[p(X| W) = \prod_{t=1}^N|\det W| \prod_{i=1}^m p([W x_t]_i) .\] The log-posterior gradient is \[\nabla_{W} \log p(W|X) = N (W^\top)^{-1} - \sum_{t=1}^N \tanh\!\left(\frac{1}{2} Wx_t\right) x_t^\top - \lambda W ,\] where the \(\tanh\) is applied element-wise, and \(X = [x_1,\dots, x_N]\in\mathbb{R}^{m\times N}\). We use \(m=3\) independent source signals \(s_t^1, s_t^2, s_t^3\) and \(N=400\) observations, yielding a \(d = m^2 = 9\) dimensional posterior. Stochastic gradients use minibatch size \(n = 40\). The initial condition is a diagonal matrix with entries \(\pm 1\) chosen uniformly at random. The posterior is strongly multimodal with symmetric modes (14), making this the most challenging test case. Reference samples and expectations are obtained from long MALA chains with \(K=10^7\) steps and step size \(h=5\times 10^{-5}\), using Riemannian manifold preconditioning together with a fixed irreversible perturbation. Our test runs use \(K=10^{7}\) steps with step size \(h=2\times10^{-4}\).

The Fisher information matrix \(\mathcal{F}\), estimated offline for this example, exhibits a \(9\times 9\) block-diagonal structure with strong intra-block correlations and weak inter-block coupling, reflecting the three independent source signals in the ICA model. See 10.4 for more details.

4.4.0.1 Experiment 1: Trajectory visualization.

15 shows estimated marginal distributions and chain trajectories along the first coordinate \(W_{11}\). unperturbed stays trapped in a subset of the modes; rand-S/M/L switch frequently between nearby modes but still exhibit visible bias; and spec leads to blowing-up chains. spec-E better captures the global symmetry and promotes more balanced transitions among the modes, yielding a closer approximation to the target distribution. These qualitative results show, again, that irreversibility alone does not guarantee global exploration; the perturbation must also be chosen to avoid discretization-induced bias.

4.4.0.2 Experiment 2: MSE versus iteration.

16 shows the evolution of MSE, squared bias, and variance of estimators of various posterior expectations, for increasing MCMC simulation time \(T\) (and fixed step size \(h\)). Across the statistics considered, spec-E is consistently among the best-performing methods, achieving fast convergence and low MSE. This quantitative comparison reinforces the qualitative findings above.

5 Conclusion↩︎

We have presented a systematic method for constructing irreversible perturbations of time-discretized Langevin dynamics, designed to accelerate mixing while controlling numerical error. The perturbations we propose are solutions to a constrained optimization problem, wherein we minimize a weighted expected squared jump distance (with weight equal to the global Fisher information matrix) subject to a spectral-gap constraint motivated by the Gaussian setting. We show how to construct explicit solutions to this problem. A key insight is that spectral optimality alone is insufficient for good performance of ULA: spectrally optimal perturbations that do not account for discretization can amplify bias and often lead to divergence at even moderate step sizes. By jointly optimizing for mixing speed and numerical stability, our method achieves substantially lower mean squared error than alternative irreversible perturbations, across diverse target distributions.

Several directions merit further investigation. First, the irreversible perturbations \(J\) developed here are state-independent. Extending to state-dependent perturbations \(J(x)\) opens up numerous possibilities: adapting to local geometry, controlling higher-order discretization errors, and leveraging structure-preserving integrators. Second, though the present paper is focused solely on the design of irreversible perturbations, it is natural to consider combining these perturbations with reversible preconditioning. These different classes of perturbations may have complementary strengths; indeed, our preliminary experiments in 11.0.0.4 suggest that irreversible perturbations can be more robust than reversible preconditioning when local and global geometry diverge. A complete theoretical analysis is beyond the scope of this work, but would be interesting in the future. Third, our analysis has focused on overdamped Langevin dynamics; lifting to the underdamped setting could exploit hypocoercivity through skew-symmetric position-momentum couplings for further mixing improvements. Finally, interacting particle systems offer another opportunity for the design of more effective algorithms. Ensemble-informed irreversible perturbations, e.g., using ensemble estimates of the Fisher information, could enhance exploration of complex multi-modal landscapes.

6 Reversible perturbations↩︎

Reversible perturbations precondition the Langevin dynamics by introducing a symmetric positive-definite matrix \(A(x) \succ 0\) that rescales and reshapes the state space. The preconditioned overdamped Langevin diffusion is \[\label{eq:preconditioned-LD} \mathrm{d}X_t = \big[-A(X_t)\nabla V(X_t) + \nabla \cdot A(X_t)\big]\mathrm{d}t + \sqrt{2A(X_t)}\,\mathrm{d}W_t,\tag{11}\] where the divergence correction \(\nabla \cdot A(X_t)\) ensures that \(\pi \propto e^{-V}\) remains the invariant measure. The diffusion coefficient \(\sqrt{2A(X_t)}\) denotes any matrix square root satisfying \(\sqrt{2A}\,\sqrt{2A}^\top = 2A\). When \(A(x) = I\), equation 11 reduces to the standard overdamped Langevin dynamics 1 . This framework is also known as Riemannian manifold Langevin dynamics (RMLD) [17].

For a constant preconditioner \(A(x) \equiv A\), the divergence term vanishes and the dynamics simplify to \[\label{eq:const-precond-LD} \mathrm{d}X_t = -A\nabla V(X_t)\,\mathrm{d}t + \sqrt{2A}\,\mathrm{d}W_t.\tag{12}\] The corresponding ULA discretization is \[X_{t+1} = X_t + hA\nabla\log\pi(X_t) + \sqrt{2hA}\,Z_t, \qquad Z_t \sim \mathcal{N}(0,I).\] The preconditioner \(A\) reshapes the effective geometry of the target: choosing \(A\) to align with the local curvature of \(V\) can significantly accelerate convergence. [23] showed that the ESJD-optimal constant preconditioner is \(A^* \propto \mathcal{F}^{-1}\), where \[\mathcal{F} = \mathbb{E}_\pi\!\left[\nabla\log\pi(X)\,\nabla\log\pi(X)^\top\right]\] is the Fisher information matrix, which captures the average local curvature of \(\log\pi\). For a Gaussian target with precision \(S\), the Fisher information matrix recovers \(\mathcal{F} = S\) exactly.

Our approach extends this ESJD-based framework from reversible to irreversible perturbations. While reversible preconditioning rescales the geometry, irreversible perturbations introduce rotational dynamics that break detailed balance and induce probability currents. The two strategies are complementary, and both benefit from Fisher-informed design.

7 Derivations↩︎

7.1 Mean and covariance of Gaussian ULA↩︎

Consider the ULA iteration for a Gaussian target with mean \(\mu\) and covariance \(S^{-1}\): \[X_t = X_{t-1} - h(I+J)S(X_{t-1} - \mu) + \sqrt{2h}\,Z_{t-1}, \qquad Z_{t-1} \sim \mathcal{N}(0, I).\] Define \(A :=(I+J)S\) and \(B := I - h A\). Subtracting \(\mu\) and unrolling the recursion gives \[X_t - \mu = \big( I - h (I+J) S \big)^t (X_0 - \mu) + \sqrt{2h} \sum_{i=0}^{t-1} \big( I - h (I+J) S \big)^{t-1-i} \, Z_i.\]

7.1.0.1 Mean.

For Gaussian initialization \(X_0 =\mathcal{N}(\mu_0,\tau^2 I)\), taking expectations yields: \[\mathbb{E}[X_t] - \mu = B^t (\mu_0 - \mu).\] The mean converges geometrically to \(\mu\) at a rate determined by the spectral radius \(\rho(B)\), whenever the initial mean \(\mu_0\) differs from \(\mu\). For small \(h\), \(\rho(B)\) is governed to first order by the real parts of the eigenvalues of \(A=(I+J)S\), matching the continuous-time spectral-gap criterion.

7.1.0.2 Covariance.

\[\begin{align} \Sigma_t &= \mathop{\mathrm{\mathbb{E}}}[(X_t-\mu_t)(X_t-\mu_t)^\top]\\ &=\mathop{\mathrm{\mathbb{E}}}\left[\Bigl(B^t(X_0-\mu)+\sqrt{2h}\sum_{i=0}^{t-1}B^{t-1-i}\xi_i-(\mu_t-\mu)\Bigr)\Bigl(\cdots\Bigr)^\top\right]\\ &=\mathop{\mathrm{\mathbb{E}}}\left[\Bigl(B^t(X_0-\mu_0)+\sqrt{2h}\sum_{i=0}^{t-1}B^{t-1-i}\xi_i\Bigr) \Bigl(\cdots \Bigr)^\top\right]\\ &=\tau^2B^t(B^\top)^t + 2h\sum_{i=0}^{t-1}B^i(B^\top)^i. \end{align}\] The covariance satisfies the recursion \(\Sigma_{t+1} = B \Sigma_t B^\top + 2h I\). As \(t\rightarrow\infty\), if \(\rho(B)< 1\), then \(\lim_{t\rightarrow\infty}\Sigma_t = \Sigma_h\), which satisfies the discrete Lyapunov equation \[\label{eq:discrete-lyapunov} \Sigma_h = B \Sigma_h B^\top + 2h I.\tag{13}\]

7.1.0.3 Expansion in \(h\).

Substituting \(B = I-h A\) into 13 yields \[\Sigma_h = \Sigma_h - h(A\Sigma_h + \Sigma_h A^\top) + h^2 A\Sigma_h A^\top + 2h I.\] Rearranging and dividing by \(h\): \[A\Sigma_h + \Sigma_h A^\top = 2I + h A\Sigma_h A^\top.\] Expanding \(\Sigma_h = \Sigma^{(0)}+h\Sigma^{(1)}+\mathcal{O}(h^2)\) and matching powers of \(h\), \[\begin{align} \mathcal{O}(1): &\quad A\Sigma^{(0)} + \Sigma^{(0)}A^\top = 2I, \\ \mathcal{O}(h): &\quad A\Sigma^{(1)} + \Sigma^{(1)}A^\top = A\Sigma^{(0)}A^\top. \end{align}\] For the \(\mathcal{O}(1)\) equation, one can verify \(\Sigma^{(0)} = S^{-1}\) is a valid (and unique) solution. For the \(\mathcal{O}(h)\) equation, substituting \(\Sigma^{(0)} = S^{-1}\), \[\begin{align} A\Sigma^{(1)} + \Sigma^{(1)}A^\top &= AS^{-1}A^\top = (I+J)S(I+J)^\top =: \Delta,\\ \Sigma^{(1)} &= \int_0^\infty e^{-tA} \Delta \, e^{-tA^\top} \, dt. \end{align}\] In summary, \(\Sigma_h = S^{-1} + h\Sigma^{(1)} + O(h^2)\), where the leading discretization error \(\Sigma^{(1)}\) is driven by \(\Delta = (I+J)S(I+J)^\top\). Taking the trace: \[\mathrm{Tr}(\Delta) = \mathrm{Tr}(S) + \mathrm{Tr}(JSJ^\top),\] where \(\mathrm{Tr}(JSJ^\top)\) is the \(J\)-dependent term. This motivates minimizing \(\mathcal{E}_1(J)=\mathrm{Tr}(JSJ^\top)\) or related quantities to control discretization bias.

7.2 Discretization error in Itô-Taylor expansion↩︎

In this section we derive the local truncation error of Euler–Maruyama scheme using the Itô–Taylor expansion, providing an alternative perspective on discretization bias.

7.2.0.1 General Itô–Taylor expansion

Consider the SDE \[\mathrm{d}X_t = b(X_t) \, \mathrm{d}t + \sigma(X_t) \, \mathrm{d}W_t, \qquad X_0 \sim \pi_0.\] For a smooth test function \(f\), define \[\mathcal{L}_0 f(x) = b(x) \cdot \nabla f(x) + \frac{1}{2} \sigma(x)\sigma(x)^\top : \Delta f(x), \qquad \mathcal{L}_1 f(x)= \sigma(x) \cdot \nabla f(x).\] For a fixed time \(t\) and step size \(h>0\), let \(\Delta W_t:=W_{t+h}-W_t\). Applying Itô’s lemma to the identity function \(f(x) = x\) over the interval \([t,t+h]\) gives \[\begin{align} \label{eq:ito95formula} X_{t+h} &= X_t + b(X_t)h+\sigma(X_t)\Delta W_t+R_1,\\ R_1 &= \int_0^h \int_0^s \mathcal{L}_0 b(X_{t+z})\,\mathrm{d}z\,\mathrm{d}s + \int_0^h \int_0^s \mathcal{L}_1 b(X_{t+z})\,\mathrm{d}W_{t+z}\,\mathrm{d}s \notag\\ &\quad + \int_0^h \int_0^s \mathcal{L}_0 \sigma(X_{t+z})\,\mathrm{d}z\,\mathrm{d}W_{t+s} + \int_0^h \int_0^s \mathcal{L}_1 \sigma(X_{t+z})\,\mathrm{d}W_{t+z}\,\mathrm{d}W_{t+s}. \end{align}\tag{14}\] The four terms have \(L^2\) orders \(\mathcal{O}(h^2)\), \(\mathcal{O}(h^{3/2})\), \(\mathcal{O}(h^{3/2})\), and \(\mathcal{O}(h)\), respectively.

7.2.0.2 Application to Langevin dynamics

For the Langevin dynamics considered in this paper, \(\sigma(x) = \sqrt{2}I\) is constant, so \(\mathcal{L}_0\sigma = 0\) and \(\mathcal{L}_1\sigma = 0\). Consequently, the last two terms vanish, and we obtain \[X_{t+h} = X_t + b(X_t)\,h + \sqrt{2}\Delta W_t + \int_0^h \int_0^s \mathcal{L}_1 b(X_{t+z})\,\mathrm{d}W_{t+z}\,\mathrm{d}s + O(h^2).\] The dominant remainder term is \(\mathcal{O}(h^{3/2})\).

7.2.0.3 Linear Langevin dynamics

For the irreversibly perturbed linear Langevin SDE, we have \(b(X)=-(I+J) SX\), \(\sigma = \sqrt{2}I\). Therefore, \(\nabla b = -(I+J)S\) and \(\mathcal{L}_1 b(X_0) = -\sqrt{2}(I+J)S\). Since \(L_1 b\) is constant, the \(\mathcal{O}(h^{3/2})\) term simplifies to \[\int_0^h \int_0^s \mathcal{L}_1 b(X_{t+z})\,\mathrm{d}W_{t+z}\,\mathrm{d}s = -\sqrt{2}(I+J)S \int_0^h (W_{t+s}-W_t)\,\mathrm{d}s =: -\sqrt{2}(I+J)S\,Z_{t,h},\] where \(Z_{t,h} := \int_0^h (W_{t+s}-W_t)\,\mathrm{d}s \sim \mathcal{N}\!\left(0,\frac{h^3}{3}I\right).\)

Therefore, the local expansion from time \(t\) is \[X_{t+h} = X_t - (I+J)SX_t\,h + \sqrt{2}\Delta W_t - \sqrt{2}(I+J)S Z_{t,h} + \mathcal{O}(h^2).\]

7.2.0.4 Local truncation error

The Euler–Maruyama scheme over one step is \(\widetilde{X}_{t+h} = X_t + b(X_t)h + \sqrt{2}\Delta W_t\). For the linear Langevin dynamics, the leading local truncation error is therefore \[X_{t+h}-\widetilde{X}_{t+h} = -\sqrt{2}(I+J)S Z_{t,h} + \mathcal{O}(h^2).\] Consequently, \[\begin{align} \mathbb{E}\!\left[ \left\| \sqrt{2}(I+J)S Z_{t,h} \right\|^2 \right] &= 2\,\mathrm{Tr}\!\left[ (I+J)S \mathbb{E}[Z_{t,h}Z_{t,h}^\top] S(I+J)^\top \right] \\ &= \frac{2h^3}{3} \mathrm{Tr}\!\left[(I+J)S^2(I+J)^\top\right] \\ &= \frac{2h^3}{3}\mathrm{Tr}(S^2) + \frac{2h^3}{3}\mathrm{Tr}(JS^2J^\top), \end{align}\] where the cross terms vanish because \(S^2\) is symmetric and \(J\) is skew-symmetric. The first term is independent of \(J\), so the \(J\)-dependent contribution is \(\mathcal{E}_3(J) := \mathrm{Tr}(JS^2J^\top) =-\mathrm{Tr}\!\left[ J S S J \right]\). This motivates minimizing \(\mathcal{E}_3(J)\) as a proxy for controlling the leading local discretization error.

8 Proofs↩︎

8.1 Proof of Proposition [prop:proxy-equivalence]↩︎

8.2 Proof of Theorem [thm:P1-P2-relation]↩︎

8.3 Spectral stability under small perturbations↩︎

The following proposition shows that the spectral gap is stable under target perturbations with small Hessian. Although it is not used directly in the main paper, it provides a useful construction and key intuition for the proof of [prop:stability:spectrum] in the next subsection.

8.4 Preparatory results for Theorem [thm:spectral-gap-stability]↩︎

We begin by introducing some useful notation and preliminary results. Under the hypotheses of Theorem [thm:spectral-gap-stability], \[V(x)=\frac{1}{2}x^\top S x+r_S(x), \qquad \nabla r_S(0)=0, \qquad \delta_S=\sup_x\|\nabla^2 r_S(x)\|_2 ,\] the Fisher information matrix satisfies \[\mathcal{F}_\pi =\mathbb{E}_\pi[\nabla^2V(X)] =S+\mathbb{E}_\pi[\nabla^2r_S(X)], \qquad \|\mathcal{F}_\pi-S\|_2\leq \delta_S .\] Rewrite the same potential around the \(\mathcal{F}_\pi\)-induced quadratic reference: \[V(x)=\frac{1}{2}x^\top\mathcal{F}_\pi x+\widetilde{r}(x), \qquad \widetilde{r}(x):=r_S(x)+\frac{1}{2}x^\top(S-\mathcal{F}_\pi)x .\] Then \[\widetilde{\delta}:=\sup_x\|\nabla^2\widetilde{r}(x)\|_2 \leq \delta_S+\|S-\mathcal{F}_\pi\|_2 \leq 2\delta_S .\] The centering is preserved: \(\nabla\widetilde{r}(0)=\nabla r_S(0)=0\). For the rest of the proof, write \(\mathcal{F}=\mathcal{F}_\pi\), \(r=\widetilde{r}\), and \(\delta=\widetilde{\delta}\). This reduces the theorem to the perturbative argument around the Gaussian reference \(\mathcal{N}(0,\mathcal{F}_\pi^{-1})\). For any antisymmetric matrix \(J\), let \[\mathcal{L}_J = \Delta - (I+J)\nabla V\cdot\nabla = \Delta - \nabla V \cdot \nabla - (J \nabla V) \cdot \nabla\] be the nonreversible Langevin generator on \(L^2(\pi)\), where \(\pi \propto e^{-V}\). The theorem concerns \(J=J^*\). The Courant-Fischer variational principle used in the proof of Proposition [prop:stability:spectrum] applies only to self-adjoint (reversible) generators, so for \(\mathcal{L}_{J^*}\) we use resolvent perturbation techniques for sectorial operators [55].

As in the proof of Proposition [prop:stability:spectrum], set \(R = \nabla r\). Since \(\nabla r(0)=0\), we have \(\nabla V = \mathcal{F} x + R\) and \(\|R(x)\| \leq \delta\|x\|\). Let \(U:L^2(\pi)\to L^2(dx)\) denote the ground-state transform \((Uf)(x)=e^{-V(x)/2}f(x)\), and define \[\tilde{L}_J := -U\mathcal{L}_J U^{-1}.\] Then \[\label{eq:conjugated-op} \tilde{L}_J = -\Delta + \frac{1}{4}|\nabla V|^2 - \frac{1}{2}\Delta V + (J \nabla V) \cdot \nabla.\tag{18}\] Here the zeroth-order term from conjugating \((J\nabla V)\cdot\nabla\) vanishes because \((J\nabla V)\cdot\nabla V=0\). The spectra of \(-\mathcal{L}_J\) and \(\tilde{L}_J\) are identical, with \(0\) as the eigenvalue of smallest real part. Moreover, the irreversible part is skew-adjoint in \(L^2(\pi)\), since \(\nabla\cdot(J\nabla V\,e^{-V})=0\); hence \(-\mathcal{L}_J\) is accretive and every eigenvalue of \(\tilde{L}_J\) has nonnegative real part. Finally, decompose \(\tilde{L}_J = \tilde{L}_J^{(0)} + \mathcal{W}_J\), where \[\begin{align} \tilde{L}_J^{(0)} &= -\Delta + \frac{1}{4}|\mathcal{F} x|^2 - \frac{1}{2}\mathrm{Tr}(\mathcal{F}) + (J\mathcal{F} x) \cdot \nabla,\\ \mathcal{W}_J &= \frac{1}{2}(\mathcal{F} x) \cdot R + \frac{1}{4}|R|^2 - \frac{1}{2}\Delta r +(JR) \cdot \nabla . \end{align}\] Here \(\tilde{L}_J^{(0)}\) is the Gaussian reference operator associated with \(V_0(x)=\frac{1}{2} x^\top\mathcal{F} x\), i.e., the Gaussian distribution \(\mathcal{N}(0,\mathcal{F}^{-1})\), and \(\mathcal{W}_J\) collects the non-Gaussian perturbation terms.

We recall the key spectral properties of this Gaussian reference operator from [36].

We also need the relative boundedness of the non-Gaussian part.

8.5 Proof of Theorem [thm:spectral-gap-stability]↩︎

With these ingredients in hand, we can now write the proof of [thm:spectral-gap-stability].

9 Relating the Fisher information matrix to local or global geometry↩︎

For a mixture \(\pi = \sum_{k} w_k \pi_k\) with well-separated components, the density at most points is dominated by a single component \(\pi_i\), so \(\nabla \log \pi(x) \approx \nabla \log \pi_i(x)\) and the Fisher information matrix decomposes as \(\mathcal{F} \approx \sum_k w_k \mathcal{F}_k\), where \(\mathcal{F}_k := \mathbb{E}_{\pi_k}[\nabla \log \pi_k \, \nabla \log \pi_k^\top]\). As the modes overlap, this decomposition breaks down and \(\mathcal{F}\) transitions toward the geometry of the merged distribution.

17 visualizes this transition using an ellipse associated with \(\mathcal{F}\). In this example, the off-diagonal entries vanish, so \(\mathcal{F}\) is diagonal and the ellipse is aligned with the coordinate axes. The lengths of its two semi-axes are determined by the two diagonal entries of \(\mathcal{F}\), thereby illustrating how the local geometry encoded by the Fisher matrix changes as the modes merge.

This also highlights the limitation of Fisher information-based preconditioning. Although \(\mathcal{F}^{-1}\) adapts to local anisotropy, \(\mathcal{F}\) primarily reflects within-mode geometry and does not directly encode the inter-modal directions needed for global exploration.

10 Additional experiments and numerical details↩︎

10.1 Non-isotropic Gaussian↩︎

Additional results for MSE, bias, and variance across different sampling methods under fixed budget and fixed simulation time are provided in [fig:app:exp1:1] and [fig:app:exp1:2]. [fig:app:exp1:3] further reports estimation error for other statistics corresponding to different proxies.

10.2 Mixture of Gaussians↩︎

The target density is \[\pi(x) \;=\; \sum_{k=1}^{5} w_k \, \mathcal{N}(x \mid \mu_k, \Sigma_k),\] where the mixture parameters are: \[\begin{align} &\text{Means:} \quad \mu_1 = (-30, 0, 0)^\top, \quad \mu_2 = (-15, 0, 0)^\top, \quad \mu_3 = (0, 0, 0)^\top, \\ & \mu_4 = (15, 0, 0)^\top, \quad \mu_5 = (30, 0, 0)^\top, \\[0.5em] &\text{Weights:} \quad w = (0.3,\, 0.1,\, 0.2,\, 0.1,\, 0.3), \\[0.5em] &\text{Covariances:} \quad \Sigma_1 = \begin{bmatrix} 5 & 0 & 0 \\ 0 & 1 & \frac{1}{5} \\ 0 & \frac{1}{5} & \frac{1}{5} \end{bmatrix}, \quad \Sigma_2 = \begin{bmatrix} 10 & 0 & 0 \\ 0 & 1 & \frac{1}{5} \\ 0 & \frac{1}{5} & \frac{1}{5} \end{bmatrix}, \\ & \Sigma_3 = \begin{bmatrix} 5 & 0 & 0 \\ 0 & 1 & \frac{1}{5} \\ 0 & \frac{1}{5} & \frac{1}{5} \end{bmatrix}, \quad \Sigma_4 = \begin{bmatrix} 10 & 0 & 0 \\ 0 & 1 & \frac{1}{5} \\ 0 & \frac{1}{5} & \frac{1}{5} \end{bmatrix}, \Sigma_5 = \begin{bmatrix} 5 & 0 & 0 \\ 0 & 1 & \frac{1}{5} \\ 0 & \frac{1}{5} & \frac{1}{5} \end{bmatrix}. \end{align}\] See 21 for demonstration.

Additional results for MSE, bias, and variance across different sampling methods under fixed budget and fixed simulation time are provided in [fig:app:exp2:1] and [fig:app:exp2:2].

10.3 Bayesian logistic regression↩︎

The posterior of Bayesian logistic regression task is illustrated in 24, which is uni-modal but non-symmetric. Additional results for MSE, bias, and variance across different sampling methods under fixed budget and fixed simulation time are provided in [fig:app:exp3:1] and [fig:app:exp3:2].

10.4 Independent component analysis↩︎

The estimated Fisher information matrix \(\hat{\mathcal{F}} \in \mathbb{R}^{9\times 9}\) (with \(D = d^2 = 9\)) used for constructing the spectral perturbation is \[\hat{\mathcal{F}} = \begin{pmatrix} 1376.6 & 821.98 & 601.77 & 18.56 & 26.61 & 1.08 & 0.57 & -5.64 & 7.88 \\ 821.98 & 1048.80 & 387.15 & -0.57 & 5.98 & -14.91 & -7.88 & -8.41 & -16.20 \\ 601.77 & 387.15 & 544.71 & 12.76 & 17.55 & 3.13 & 24.72 & 17.63 & 16.26 \\ 18.56 & -0.57 & 12.76 & 1379.31 & 817.66 & 596.03 & 20.69 & 26.94 & -43.39 \\ 26.61 & 5.98 & 17.55 & 817.66 & 1036.79 & 380.94 & 8.62 & 13.46 & -28.55 \\ 1.08 & -14.91 & 3.13 & 596.03 & 380.94 & 539.21 & 34.54 & 32.86 & -8.04 \\ 0.57 & -7.88 & 24.72 & 20.69 & 8.62 & 34.54 & 1397.90 & 819.09 & 602.28 \\ -5.64 & -8.41 & 17.63 & 26.94 & 13.46 & 32.86 & 819.09 & 1030.88 & 380.12 \\ 7.88 & -16.20 & 16.26 & -43.39 & -28.55 & -8.04 & 602.28 & 380.12 & 536.94 \end{pmatrix}.\]

11 Irreversible versus reversible perturbations↩︎

A natural question is whether the Fisher information matrix should instead be used to construct a reversible preconditioner. The answer depends on how well the global Fisher information matrix \(\mathcal{F}\) represents the local geometry of the target. For the Gaussian and mixture-of-Gaussian targets in [sec:Gaussian,sec:mix:of:Gaussian], local and global curvature coincide, or nearly so. In these settings, reversible preconditioning with \(\mathcal{F}^{-1}\) increases the spectral gap while introducing relatively little discretization error, and is therefore expected to perform well.

The BLR posterior is different: its curvature varies across the state space, so the local Hessian \(-\nabla^2\log\pi(x)\) can differ substantially from \(\mathcal{F}\). A constant reversible preconditioner \(\mathcal{F}^{-1}\) modifies both the drift and the diffusion, so this mismatch can introduce persistent discretization bias. An irreversible perturbation, by contrast, uses \(\mathcal{F}\) only to construct a skew-symmetric drift rotation and leaves the diffusion unchanged, making it more robust to global–local geometry mismatch. The experiments below illustrate this distinction empirically on the BLR target: even with the true Fisher information matrix, our irreversible perturbation outperforms reversible preconditioning, and the gap widens in the adaptive setting, where inaccurate Fisher information matrix estimates can distort the reversible noise covariance.

All experiments use Euler–Maruyama discretization with step size \(\Delta t=2\times 10^{-3}\), \(K=20{,}000\) steps, \(M=32\) chains, and stochastic gradients with mini-batch size \(n=10\). We compare standard ULA with oracle and adaptive versions of the optimal position-independent reversible perturbations proposed by [23] and with oracle and adaptive versions of our optimal position-independent irreversible perturbations spec-E. For adaptive methods, \(\mathcal{F}\) is estimated online using the streaming estimator from 3.5.1.2, and the preconditioner or perturbation is updated every 100 steps. For a paired comparison, all methods use the same Gaussian noise and mini-batch indices within each chain. We report the MSE of running estimates of \(\mathbb{E}[\mathrm{sum}(X^{(i)})]\) and \(\mathbb{E}_\pi[\|X\|^2]\).

We emphasize that this comparison is empirical and target-dependent. A more complete theoretical characterization of when irreversible perturbations are preferable to reversible preconditioning is beyond the scope of this paper and remains an interesting direction for future work.

11.0.0.1 Robustness to initialization

27 compares standard ULA, oracle/adaptive reversible preconditioning, and oracle/adaptive irreversible perturbations, initialized at \(X_0=s\mathbf{1}_d\) with \(s\in\{0,1,5,10\}\). Adaptive reversible preconditioning stalls across initialization distances, while adaptive irreversible perturbation converges reliably and tracks its oracle version. Even with the true Fisher information matrix, the irreversible perturbation outperforms reversible preconditioning, indicating that \(\mathcal{F}^{-1}\) is a poor constant reversible preconditioner for this spatially varying posterior.

11.0.0.2 Trajectory behavior

[fig:traj-oracle,fig:traj-adaptive] visualize representative trajectories projected onto the top two principal components of the posterior. The oracle irreversible perturbation produces a spiraling trajectory that explores posterior level sets more broadly than reversible preconditioning. In the adaptive case, reversible preconditioning collapses into a trapped cluster for distant initializations, whereas the adaptive irreversibly perturbed chain continues to explore and moves toward the posterior mean.

11.0.0.3 Robustness to inaccurate Fisher estimates

To isolate the effect of FIM quality, we corrupt the Fisher information matrix as \[\mathcal{F}_\alpha = \alpha \, \mathcal{F}_{\mathrm{true}} + (1 - \alpha) \, I_d, \qquad \alpha \in \{0, 0.25, 0.5, 0.75, 1\},\] and initialize all chains at \(X_0 = \mathbf{0}\). [fig:robustness-fim-curves-blr,fig:robustness-fim-summary-blr] show that irreversible perturbations consistently outperform reversible preconditioning, both with oracle and corrupted Fisher matrices. The adaptive irreversible perturbation achieves the lowest MSE across corruption levels, while adaptive reversible preconditioning again stalls.

11.0.0.4 Analysis of adaptive failure modes

Failure of the adaptive reversible method (and success of the adaptive irreversible method) can be ascribed to how the two methods use an inaccurate Fisher information matrix estimate. Reversible preconditioning inverts \(\hat{\mathcal{F}}\) and uses \(\hat{\mathcal{F}}^{-1}\) for both drift and noise; when \(\hat{\mathcal{F}}\) is ill-conditioned, this suppresses exploration in important directions and prevents recovery. An irreversible perturbation based on \(\hat{\mathcal{F}}\) may also produce a poorly calibrated rotational term in the drift, but the diffusion remains isotropic, so the chain continues to explore and the Fisher estimate can improve. This drift–diffusion decoupling can explain the greater robustness of irreversible perturbations in adaptive schemes.

References↩︎