{'fraction_non_alphanumeric': 0.04168156931796297, 'fraction_numerical': 0.01372878091726965, 'mean_word_length': 4.438874756845728, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic gradient descent and fast relaxation to thermodynamic equilibrium: a stochastic control approach

We study the convergence to equilibrium of an underdamped Langevin equation that is controlled by a linear feedback force. Specifically, we are interested in sampling the possibly multimodal invariant probability distribution of a Langevin system at small noise (or low temperature), for which the dynamics can easily get trapped inside metastable subsets of the phase space. We follow [Chen et al., J. Math. Phys. 56, 113302, 2015] and consider a Langevin equation that is simulated at a high temperature, with the control playing the role of a friction that balances the additional noise so as to restore the original invariant measure at a lower temperature. We discuss different limits as the temperature ratio goes to infinity and prove convergence to a limit dynamics. It turns out that, depending on whether the lower ("target") or the higher ("simulation") temperature is fixed, the controlled dynamics converges either to the overdamped Langevin equation or to a deterministic gradient flow. This implies that (a) the ergodic limit and the large temperature separation limit do not commute in general, and that (b) it is not possible to accelerate the speed of convergence to the ergodic limit by making the temperature separation larger and larger. We discuss the implications of these observation from the perspective of stochastic optimisation algorithms and enhanced sampling schemes in molecular dynamics.

1. Introduction. Methods to accelerate the convergence of Monte Carlo algorithms that sample the equilibrium distribution of a possibly high-dimensional system have been around since at least the invention of the Metropolis algorithm [40]. Yet, they have recently enjoyed an enormous boost in attention, partly driven by the latest developments in machine learning [30,47] and theoretical advances in hypocoercive systems [23,45]. Depending on the scientific community, accelerated sampling algorithms come under various different names, e.g. enhanced sampling [3], smart sampling [29], guided sampling [43], to mention just a few. They come into play when the target probability distribution is not just high-dimensional, but also multimodal, which is a challenge for any Monte Carlo algorithm [31]. Similar issues arise in machine learning and the (stochastic) optimisation of high-dimensional, non-convex functions, and it has been argued that Monte Carlo sampling can improve the training of large-scale networks by reducing overfitting and accelerating convergence [25,30].
A model for Monte Carlo simulations that is popular in statistical mechanics and that shares features of Brownian dynamics and ballistic motion dynamics is the underdamped Langevin model [4]. The underdamped Langevin dynamics is a second-order stochastic differential equation (SDE) and is based on Hamilton's equations perturbed by friction and noise. Besides the physical appeal of the model, it has been conjectured that the underdamped Langevin dynamics shows superior convergence over its first-order counterpart, the so-called overdamped Langevin equation, and this conjecture has been recently proved using the 2-Wasserstein distance [9]. The fact that the (underdamped) Langevin equation is a second-order SDE makes it harder to analyse from a theoretical point of view, because the underlying infinitesimal generator is degenerate and non-symmetric. However, it is precisely the degeneracy and the lack of symmetry that offers the possibility to tweak the dynamics by changing the drift in a way that the convergence to equilibrium is accelerated without altering the invariant measure [12,14]; see also [22,26,39] for related ideas.
In view of the aforementioned, two recent developments are relevant: sampling methods based on either optimal control or mean field approaches [20,21], and the application of Monte Carlo methods from statistical mechanics-especially molecular dynamics-to problems in machine learning or Bayesian inference [16,33]. Some of the modern stochastic optimisation methods from machine learning, like ADAM, AdaGrad or RMSProp adaptively control the learning rate so as to improve the convergence to a local minimum, but they also share many features with adaptive versions of the Langevin equation [17,25]. It therefore comes as no surprise that the Langevin model is the basis for the stochastic modified equation approach that can be used to analyse state of the art momentum-based stochastic optimisation algorithms like ADAM [1,27]. The common feature of all these algorithms is that they are designed to adaptively sample only relevant parts of a complicated energy landscape or loss function in order to accelerate convergence to equilibrium (i.e. to an equilibrium distribution or an approximation of the global optimum). Control theory plays a key role here in that many of the aforementioned adaptive algorithms can be interpreted and analysed as solutions to stochastic optimal control problems [15,18,35]; cf. also [32].
In this paper, we take up an idea that was put forward in [6] and that is reminiscent of the Schrödinger bridge approach (e.g. [10,38]). Inspired by recent ideas [7,8] for solving the Schrödinger bridge problem for a linear SDE, the idea here is to augment the (nonlinear) Langevin equation by a linear feedback control that plays the role of bridging the equilibrium distribution of an uncontrolled Langevin equation at a high ("simulation") temperature to a non-equilibrium target probability measure at a lower ("target") temperature on an infinite time horizon. The control target (i.e. the non-equilibrium invariant measure) is reached asymptotically in the long time limit by letting the control act as a friction force that balances the additional noise in the high-temperature Langevin equation so as to restore the invariant measure at the target temperature. The approach rests on standard arguments from equilibrium statistical mechanics (specifically: fluctuation-dissipation relations), combined with functional inequalities (specifically: FIR inequalities) from which an ergodic optimal control problem is derived. The latter leads to the minimum dissipation approach of [6].
Adopting a slightly different perspective, the control approach can be thought of as a nonasymptotic version of simulated annealing, in which the control target is to efficiently sample the modes of a probability density without getting trapped in local modes too often. The question of finding good temperature protocols in simulated annealing is intimately related to the question of accelerating the speed of convergence towards the stationary distribution by optimising the spectral gap of the underlying generator [19,28]. Having these two sides of a coin, it is thus natural to look at the limit as the ratio between the simulation temperature and the target temperature goes to infinity. It turns out that, depending on whether the target or the simulation temperature is fixed, the controlled dynamics converges either to the overdamped Langevin equation or to a deterministic gradient flow. This implies that the ergodic limit and the large temperature separation limit do not commute in general. As a consequence, the controlled Langevin dynamics interpolates between sampling and stochastic gradient descent (under an appropriate temperature scaling), and this fact can be used to devise algorithms for non-convex optimisation problems. On the other hand, it implies that it is not possible to accelerate the system indefinitely by increasing the simulation temperature under the constraint of asymptotically reaching the stationary distribution, which is in agreement with related results in [12,14].
The article is outlined as follows. In Section 2, we introduce the Langevin model with linear feedback controls and derive optimality conditions for the target invariant measure. The large temperature limits of the model are studied in Section 3. All theoretical findings are discussed and illustrated with numerical examples in Section 4. We conclude with a short discussion in Section 5. The article contains an appendix that records various technical identities, derivations and proofs.
2. High-temperature Langevin dynamics. Let (X, Y ) = (X t , Y t ) t≥0 be the solution of the controlled Langevin equation Here V : R n → R is a smooth potential energy function that is bounded from below and coercive, W t is a standard Brownian motion in R n , γ ∈ R n×n is a symmetric positive definite matrix, and the noise coefficient σ ∈ R n×n satisfies the fluctuationdissipation relation for someβ > 0. The control u assumes values in R n and will be specified below; for u ≡ 0 the dynamics admits the unique invariant measure (omitting normalisation constants here and in what follows) with the Hamiltonian We call 1/β the "simulation temperature". Our aim is to simulate (1) at inverse temperatureβ, and choose a control so that the dynamics converges to a probability distribution µ ∞ at a lower temperature 1/β that we call the "target temperature". We follow the approach proposed in [6], and consider controls of the form u t = B Y t for some suitable matrix B ∈ R n×n . The following Lemma characterises the admissible matrices B that achieve the aforementioned aim.
In what follows the operators ∇ x , ∇ y are the gradients and ∇ 2 y is the Hessian with respect to the variable in the subscript, where x, y ∈ R n . For two matrices D, C, we use D : C = tr(D C).  Bσ + σB = (β − β)σσ .
Proof. The infinitesimal generator of (1) on C 2 has the form L = L 1 + L 2 , with where the Hamiltonian part, L 1 , is skew-symmetric with respect to the standard L 2 scalar product. Hence it suffices to show that, for any smooth and integrable function f , if and only if (5) holds. Using integration by parts, we obtain Since ρ ∞ > 0, it follows by (2) that g = 0 if and only if the symmetric part of the matrix σB T is equal to 1 2 (β − β)σσ or, in other words, Hence the assertion is proved.
Equation (5) is a T -Sylvester equation for which existence and uniqueness results can be found in [46]. Unlike the standard Lyapunov equation for a symmetric unknown, its solution need not be unique. In fact, note that B = 1 2 (β − β)σ + σM is a solution to (5) for any skew-symmetric matrix M = −M . However, since σσ is positive definite and, hence, σ is invertible there is also a unique solution such that σB is symmetric. This particular solution corresponds to M = 0 and will be of interest in the next section.
Roughly speaking, the matrix σB plays the role of a friction coefficient that can cause positive or negative dissipation, depending on whether β >β or β ≤β. Typically, we want to speed up the convergence to the target distribution µ ∞ and so the system is simulated at a temperature 1/β that is higher than the target temperature 1/β, i.e. we have β >β. Then, the eigenvalues of Bσ + σB have strictly negative real part and therefore the control term σu dissipates energy. If β =β, equation (5) implies that σB must be skew-symmetric.
We will henceforth refer to u as either "control", "friction" or "dissipation".

A stochastic control problem.
There are many matrices B satisfying the relation (5). For example, if β =β, any matrix B with the property Bσ = −σB preserves the invariant measure µ ∞ =μ ∞ , especially the zero matrix B = 0. In fact, the symmetric part of Bσ is unique, whereas the anti-symmetric or skew-symmetric part is not, which is in accordance with previous works on non-reversible perturbations of reversible diffusions, e.g. [22] We now identify a unique control that has a property that it minimises the relative entropy between the controlled dynamics and an uncontrolled Langevin system that is initialised at the correct equilibrium distribution, with the aim of accelerating the convergence to equilibrium. The stochastic control approach leads to a Langevin dynamics that has a minimum-dissipation property, in that it is as close as possible to a reversible system in a wide sense (i.e. including momentum flips [12]), since the friction term becomes symmetric again.
Letting ρ t denote the probability density function of the marginal law of the controlled process (1) at time t ≥ 0. Then, by Lemma 5 in Appendix A, it follows that where KL(·, ·) as defined in (31) denotes the Kullback-Leibler divergence (or relative entropy) between two probability densities, and is the asymptotic entropy production (AEP) rate when the uncontrolled process is started in equilibrium; see Appendix A.2 for details [5,42]. Here, |z| g = z gz for a symmetric and positive definite matrix g, and E[·] denotes the expectation with respect to the target density ρ ∞ . Our aim is to minimise the AEP rate R over the set of admissible controls which boils down to a minimisation over matrices B ∈ R n×n that are of the form for some skew-symmetric matrix M = −M . The AEP functionalR(B) = R(B T y) is strictly convex with the unique minimumR(B * ) = 0 that is attained at The optimal matrix B * has been computed in [6] by minimising a slightly different functional, however. Combining equation (11) with (2), we see that the controlled system that is simulated at temperature 1/β satisfies a fluctuation-dissipation relation with parameter β: 3. Large temperature separation limit. We are interested in the limit β/β → ∞, which can either correspond to the high simulation temperature limitβ → 0 when the target temperature 1/β is held fixed, or to the zero target temperature limit β → ∞ when the simulation temperature 1/β is fixed.
Using the relation (2), it is convenient to write (1) under the constraint (11) as Note that the scaled friction matrix D satisfies the fluctuation-dissipation relation at the correct target temperature β, Fixed simulation temperature. We introduce the dimensionless parameter ε =β/β ∈ (0, 1]. Upon rescaling time according to t → t/ε, we can recast our Langevin equation (13) as We are interested in the limit ε → 0 when both σ and γ are held fixed. It turns out that, in this limit, the x-component of (14) mimics a stochastic gradient descent algorithm for small ε in the sense that, as ε → 0, it converges to the solution of the gradient system (15) d dt This can be seen as follows. Integrating the second equation of (14) from zero to t yields which, together with the first equation leads to and therefore formally passing ε → 0 we arrive at the limit equation. We now give the precise statement.
Theorem 1. Let ∇V be Lipschitz continuous on R n with constant L V . Furthermore let X = X ε and x be the solutions of (14) and (15) with X 0 = x(0). Then, for any T > 0, Moreover we have the pre-asymoptotic bound The proof is based on [41, Prop. 2.14] and we outline the main steps in Appendix B.
Fixed target temperature. By construction, under the limit β/β → ∞ for fixedβ the target measure is changed. If, however β/β → ∞ while β is held fixed, then the target measure is obviously unaltered. We define ς := σ β /β, with βςς =βσσ = 2γ such that ς is independent ofβ. After rescaling time again according to t → t/ε, we obtain Note the 1/ε scaling of the noise term that is different from the 1/ √ ε scaling in (14). We are interested in the limit ε → 0 when both ς and γ are held fixed. (Note that the noise coefficient ς is independent ofβ.) It is well-known that, under this scaling, the x-component converges pathwise to the solution of the overdamped Langevin equation (19) γdX The following result is due to Nelson [34, Thm. 10.1].
Theorem 2 (High-friction limit). Let ∇V be Lipschitz continuous on R n . Further let X = X ε andX be the solutions of (18) and (19) respectively with X 0 =X 0 . Then for any

Numerical examples.
We consider three different examples, with increasing degree of complexity: a linear system where we can explicitly control the speed of convergence, a highdimensional bistable system, for which we want to compute a slowly converging marginal, and a Lennard-Jones cluster for which we seek a minimum energy configuration.
4.1. Optimal temperature ratio. We consider a linear Langevin system of the form (20) dZ t = AZ t dt + CdW t with the matrices where γ, K ∈ R n×n are symmetric positive definite matrices. Using (2) together with (11), the coefficient matrices can be recast as Note that (20) is basically the controlled Langevin equation (13) with a quadratic Hamiltonian. Hence, by construction, it has a unique Gaussian invariant measure We are interested in the optimal temperature ratio β/β that maximizes the speed of convergence to µ ∞ . Using [2, Theorem 4.9] and the Cziszar-Kullback-Pinsker inequality, the convergence is exponential with a rate given by the real part of the principal eigenvalue of A (spectral abscissa). Specifically, letting −r = max{Re(λ) : Av = λv} and calling η t the probability density associated with the law of Z t at time t > 0, we have, for some C > 0, where ρ ∞ = η ∞ denotes the density associated with the invariant measure µ ∞ and KL(η|ρ) denotes the Kullback-Leibler (KL) divergence or relative entropy between two probability measures with densities η and ρ; see Definition 4 below. We want to determine the inverse temperature ratio β/β that brings the principal eigenvalue of A as far away as possible from the imaginary axis, i.e., we want to minimize the largest real part of the eigenvalues. To begin with, we consider the case n = 1. The two eigenvalues A are given by and they are either negative real or form a conjugate complex pair, since tr(A) < 0 and The optimal rate 2r * = tr(A) is therefore achieved when the two eigenvalues are equal and real, which is exactly the case when (tr(A)) 2 = 4 det(A) or, equivalently, This is to say that, for given K, the optimal ratio between simulation temperature and target temperature scales as 2/γ with the friction coefficient γ; for example, when K = γ = 1, the optimal simulation temperature is twice the target temperature. Note, however, that the optimal temparature ratio can be smaller than 1 if the friction coefficient is large. Illustration. We illustrate the previous considerations with (20) where we choose the coefficients for β = 1 and random initial conditions Z 0 , with The KL divergence between the two centered Gaussians η t and ρ ∞ is readily seen to be [36] KL(η t |ρ t ) = Figure 1 shows the results of a comparison for different values ofβ. It is clearly visible that the optimal valueβ = β/2 beats the equilibrium simulationβ = β by several orders of magnitude (purple curve vs. blue curve) whereas the deviation between the distributions is clearly larger for other suboptimal temperature ratios (light blue and orange curves). Furthermore, a larger simulation temperature need not necessarily increase the speed of convergence as shown by the orange curve.
Multidimensional linear system. We briefly discuss a generalisation of the previous considerations to the case n > 1. Generally, since the eigenvalues of A are either real or occur in complex conjugate pairs, the optimal eigenvalue bound and thus the maximally achievable rate of convergence in (23) is when all eigenvalues have the same real part, i.e. when To see this we decompose the matrix A ∈ R 2n×2n in (22) into symmetric and antisymmetric parts according to The fact that tr(A) is a real number now implies that where m j denotes the algebraic multiplicity of the eigenvalue λ j . By construction, the sum of m j equals 2n, and so Note, however, that this bound does not imply that the temperature ratio can be made arbitrarily large so as to improve the rate of convergence indefinitely. By Theorem 2, the convergence rate for infinitely large temperature ratios equals the convergence rate of the overdamped Langevin equation at the target temperature. Moreover, the optimal convergence rate cannot always be reached. For example, let A ∈ R 4×4 be defined as in (22), with Introducing the shorthand α = β/β, the eigenvalues of the matrix are The optimal temperature ratio α that minimises the real part of the maximum eigenvalue is α * = 4/3 ≈ 1.1547, and it leads to the eigenvalues (see Figure 2) Clearly, for any given simulation temperature (i.e. when α is constant), the spectral abscissa cannot be smaller than r * = tr(A)/4 ≈ −0.8660, which determines the theoretically optimal speed of convergence (see page 8). We observe that the optimal rate of convergence is not reached in this case. It is unclear, under which conditions this rate can be generally reached, and we refer to [44] for numerical algorithms for spectral abscissa minimization.
The following can be shown if the potential energy matrix K and friction matrix γ have a special form. Lemma 2. Let K and γ commute. Then there exists a matrix S ∈ R n×n such that SKS −1 = diag(k 1 , . . . , k n ) and SγS −1 = diag(g 1 , . . . , g n ), where without loss of generality g 1 ≤ g 2 ≤ . . . ≤ g n . The optimal convergence rate is then attained for Proof. Following the above considerations (see also [2, Theorem 4.9]) we are looking for the maximum real part of the eigenvalues of the drift matrix A that will be minimised over α = β/β. To this end we consider the eigenvalue problem for the matrix A Rewriting the last equation in the basis given by the simultaneous diagonalization S yields the following system of equations Hence, for each i we find a pair of eigenvalues λ 1,2 Since we are looking for the maximum real part, we only consider  Note that in case r = 4k 1 /g 2 1 there is no intersection at α * = r = argmin α Re(λ 1 (α)), but Re(λ 1 (α)) attains its global minimum at α * .
Remark 3. The optimal temperature ratio (24) is a special case of Lemma 2 for n = 1.

4.2.
Fast sampling of a bistable system. We now illustrate the previous considerations for linear systems with a nonlinear toy system at low friction. To this end, we consider the Langevin equation of dimension n = d + 1 (d = 100), with the potential for a k > 0. The friction coefficient γ = γ T > 0 is given by a tridiagonal matrix having 0.04 on the main diagonal and 0.02 above and below. In what follows we abbreviate the configuration components as x = (q, ξ) and the velocities or momenta (strictly speaking) as y = (u, ζ).
We run the controlled Langevin dynamics forβ = 1.0 and compare with an uncontrolled simulation at target temperature β = 5.0. Note that in this situation-other than in the linear case-no analytical result regarding the optimal temperature is available, so we run several numerical tests first, from which we determine a temperature ratio β/β ≈ 5 that shows good convergence towards the equilibrium distribution in the eyeball norm.
For all numerical simulations, we use a variant of the BAOAB scheme with step size ∆t = 5 · 10 −3 (see [24]), and we record the empirical q-marginals as well as the velocity autocorrelation function for the first coordinate after 4 000 000 steps using a subsample of 20 000 points. Figure 3 shows a typical trajectory for the first and the third components q and ξ 1 . While ξ 1 is essentially Gaussian, we observe that the mixing rate between the metastable sets, i.e. the 2 wells of the double-well potential is much higher for the controlled simulation atβ < β than for the uncontrolled simulation withβ = β; this is in accordance with the well converged empirical distribution shown in the left panel of Figure 4 that should be contrasted with the result of the uncontrolled simulation that has not yet converged.
Decay of correlations. For further comparison, we also compute the velocity autocorrelation function for the q-component, i.e. U = dQ/dt where we use the shorthands Q = X 1 and U = Y 1 . Since the Langevin equation is degenerate, with the noise driving only the ycomponents, the velocity autocorrelation exists. For all practical purposes the corresponding time series is weakly stationary with mean zero and variance ϑ 2 > 0, and thus the velocity autocorrelation is unambigously defined by where the expectation is understood over the (stationary) marginal distribution of the component U t . Using ergodicity, we can approximate the autocorrelation function by for some sufficiently large T ∈ (0, ∞), with the integrals being replaced by sums when discrete time series data are used (as is the case here). Interestingly enough, despite the much higher mixing rate in the controlled simulation due to the higher noise level, even the velocity autocorrelation functions of the two simulations forβ < β andβ = β are in good qualitative agreement (see right panel of Figure 4). The fact that the autocorrelation function of the controlled simulation decays faster than for the uncontrolled simulation is in line with the observation in Section 3 that the temperature ratio plays the role of a time scale separation parameter; cf. equation (18).  where ε, σ > 0 are parameters and r ij = |x (i) − x (j) |, i = j denotes the distance between the i-th and the j-th particle with coordinates x (i) , x (j) ∈ R d , respectively. The function V is a sum of all pairwise interaction terms.
Suppose we want to find the unique minimum energy configuration x * = (x (1) , . . . , x (N ) ) ∈ R n by gradient descent in V where n = dN with N denoting the number of particles and d being the spatial dimension-typically 1, 2 or 3. (Note that x * is unique up to translations and rotations, i.e. up to rigid body symmetries.) Since the potential has many irrelevant local minima a gradient descent algorithm will get stuck in the local minimum that is closest to its starting value. The Langevin dynamics (14) approximates a gradient descent as ε → 0 that gets stuck in local minima, but for every finite ε > 0, the dynamics has has a unique ergodic invariant measure µ ε ∞ with density whereH =βH is independent of ε. (The representation (26) readily follows from Lemma 1 by rescaling the Hamiltonian.) By Laplace's principle, we have (e.g., see [11,Theorem 4.3.1]) In other words the invariant measure µ ε ∞ concentrates at the unique global minimum x * = argmin x V (x) as ε → 0. Since, by Theorem 1, the limit ε of the Langevin equation on any finite time interval [0, T ] converges to the deterministic gradient system that has the property to converge to the nearest local minimum as t → ∞, it is a straight consequence that the limits t → ∞ and ε → 0 in the two-temperature Langevin dynamics do not commute.
We will exploit this property to simulate the two-temperature Langevin dynamics (14) at large, but finite temperature separation limit, so as to mimic a stochastic gradient descent that is likely to converge to a minimum that is close to the global one in the sense that V (X T ) ≈ V (x * ) whenever T is sufficiently large. Cooling a Lennard-Jones system. We seek the unique minimum energy configuration of a Lennard-Jones cluster of N = 7 particles-or a close approximation thereof. The idea now is to exploit the fact that, for finite ε, the Langevin dynamics (14) is close to a stochastic gradient descent with momentum (like ADAM or its variants) that has the property of being absorbed by a very low minimum on O(1) time scales and to remain there over O(exp(1/ε)) time scales. We expect that for sufficiently small ε > 0, the dynamics (14) will essentially fluctuate around the global minimum x * .
As before we employ the BAOAB scheme with step size ∆t = 5 · 10 −3 and a total number of 2 000 000 time steps, where we record only every 200th step. The simulation temperature is fixed atβ = 1.0, whereas the target temperature is varied between β =β and β =β/100. Figure 5 shows the results of two simulations for β = 1.0 (left panel) and β = 0.01 (right panel), both starting from the same high energy initial configuration. It can be seen that the simulation in the stochastic gradient descent regime for β = 0.01 quickly goes down the energy ladder by roughly 11 to 12 orders of magnitude, similar to the uncontrolled equilibrium simulation, but then stabilises in a low-lying local minimum (see energy plot in the right panel of Figure 5). Since the target temperature is small compared to the typical energy barrier, the probability to escape from the attractor of the local minimum during simulation time is negligible. The uncontrolled simulation even reaches a yet lower energy state, but it oscillates wildly between local energy minima that span three orders of magnitude in energy, with a high escape probability (see energy plot in the left panel).

Discussion.
We have studied the asymptotic properties of a feedback-controlled underdamped Langevin equation at high (simulation) temperature that is controlled by a linear feedback force. The feedback force acts as a friction and balances the noise so that the dynamics weakly converges to a thermal equilibrium at a lower (target) temperature. The rationale of choosing a simulation temperature that is higher than the target temperature is that the dynamics is less likely to get trapped in local minima of the Hamiltonian compared to simulating the system at the target temperature. The approach is reminiscent of the Schrödinger bridge problem on an infinite time horizon for diffusions with different noise coefficients. By using a functional inequality approach (FIR inequality) we have identified a control that enjoys a minimum dissipation property and that has been derived in [6] using different arguments.
The fact that the simulation temperature can be different from the target temperature opens up various possibilities for designing sampling or global optimisation algorithms. For example, it is possible to ask what is the optimal temperature ratio that maximises the speed of the convergence to the desired target distribution. We have given only a partial answer to this question for Gaussian processes and collected further empirical evidence based on numerical simulations; cf. also [14] for a theoretical analysis. It should be mentioned that by choosing a too high simulation temperature, the speed of convergence can even be reduced as compared to the case when the simulation and the target temperatures are equal. On the other hand, if one is interested in exploring the minima of the underlying Hamiltonian, one could lower the target temperature, while holding the simulation temperature fixed, so that the invariant probability measure concentrates in the low-lying minima. We have numerically demonstrated with a benchmark example from molecular dynamics that the controlled Langevin dynamics can be used to efficiently solve high-dimensional and non-convex minimisation problems.
From an algorithmic point of view, a relevant question then is whether it is possible to reach a (theoretically) arbitrary speed of convergence for sampling or optimisation by sending the temperature ratio to infinity. The answer is clearly negative as has been already pointed out in [12,14]. More specifically, it turns out that the controlled dynamics converges either to the overdamped Langevin equation or to a deterministic gradient flow, depending on whether the target or the simulation temperature is held fixed. As a consequence, the ergodic limit and the large temperature separation limit do not commute in general, and the feedback-controlled Langevin dynamics interpolates between sampling and stochastic gradient descent.
Future work shall address the question of optimising temperature ratios beyond the Gaussian setting or choosing more general feedback laws. The latter may be done by rephrasing the problem as an optimal control problem for the Fokker-Planck equation which relates this line of research to recent work on mean-field approaches for optimisation (e.g. [37]).
Acknowledgement. This work has been partially supported by the Collaborative Research Center Scaling Cascades in Complex Systems (DFG-SFB 1114) through project A05 and by the MATH+ Cluster of Excellence (DFG-EXC 2046) through the projects EP4-4 and EF4-6.
Appendix A. FIR inequality and minimum dissipation functional. We want to study the density evolution of the controlled Langevin dynamics (1) and compare it to an uncontrolled simulation at the target temperature 1/β: A.1. FIR inequality. Obviously, the path measures of (1) and (28) are mutually singular because the equations are not driven by the same Brownian motion, and therefore we cannot apply Girsanov's theorem to do a pathwise comparison. (The noise coefficients are different.) Nevertheless it makes sense to compare the time t marginals of the two equations. To this end, we write the Fokker-Planck equation associated with (1) as a continuity equation (29) ∂ρ ∂t + ∇ · (vρ) = 0 , We use the shorthand notation ρ t = ρ(·, t) for the smooth probability density of (X t , Y t ) ∈ R n × R n at time t > 0 and U = (0, σu) ∈ R n × R n to denote the control that is acting on the momentum equation only, and we have introduced the (2n × 2n)-matrices J = 0 n I n −I n 0 n , Γ = 0 n 0 n 0 n γ , Σ = 0 n σ .
We further let denote the Fokker-Planck equation associated with the reference Langevin dynamics (28) at the target (inverse) temperature β withΣ = 0 n σ .
Definition 4 (Relative entropy). We define the relative entropy between ρ t and η t as and otherwise we set KL(η t |ρ t ) = ∞.
Assuming that η t is sufficiently decaying at infinity 1 and that KL(η t |ρ t ) < ∞ for all t > 0, it is straightforward to show that since (using a dot to denote the partial derivative with respect to t) Here we have substituted (29) and (30) in the second equality and used integration by parts in the third equality, assuming fast decay of the density η t which allowed us to ignore boundary terms.
Theorem 5 (FIR inequality, cf. [42]). It holds, for all τ > 0, is the relative Fisher information, with |z| G = z Gz denoting the G-weighted Euclidean vector norm with a semidefinite weight matrix G = σσ > 0, and is a variant of the Donsker-Varadhan large deviations rate functional for the empirical measure generated by an SDE (e.g. see [13,Theorem 2.5]).
Proof. We start from the representation (32). Substituting v andv we obtain by adding a zero: where in the last line we have used that only the lower right block of the matrices Σ,Σ is nonzero. We set δ = σσ −σσ and consider the second integral. Using Cauchy-Schwarz and Young's inequality, we find for all τ > 0. Taking a closer look at the leftmost integral, we see that As a consequence, A.2. Entropy production rate. In (33) we may tighten the upper bound by minimising over τ > 0. From a practical point of view, however, it is more convenient to set τ = 1, as a consequence of which the Fisher information term vanishes.
We now want to compute the entropy production of the controlled system. To this end, we set η 0 = ρ ∞ which implies that η t = ρ ∞ for all t > 0. Using that δ = (1 −β/β)σσ and integrating (33) we obtain the inequality (36) KL(ρ ∞ |ρ T ) − KL(ρ ∞ |ρ 0 ) ≤ 1 2 T 0 E σu − 1 2 δ∇ y log ρ ∞ 2 (σσ ) −1 dt , where the expectation is taken over the invariant density ρ ∞ ∝ exp(−βH). The term inside the integral is the (asymptotic) entropy production rate, and we define Note that when u is admissible, then R(u) shows up as asymptotic entropy production rate as T → ∞, even without the assumption η 0 = ρ ∞ , because the law η t of the uncontrolled dynamics (28) asymptotically approaches the target density ρ ∞ as t → ∞, independently of the initial distribution. Hence the name asymptotic entropy production rate.
For the third term we use, denoting the eigenvalues of γ by λ i , that |∇V (x r )|, where |γ −1/2 | 2 F = i≤n 1/λ i . Note that max r∈[0,t] |∇V (x r )| takes a deterministic finite value since ∇V (·) and x · are both continuous and hence attain their minimum and maximum on compact intervals. The last term can be bounded using integration by parts as Since W r is almost surely continuous, we have that |W r | ≤ K < ∞ almost surely so that evaluating the integral t 0 |e − t−r ε 2 γ | F dr yields the almost sure bound Putting everything together gives |∇V (x r )| which, using Gronwall's inequality, yields the assertion. The second assertion follows by the modified estimate for the last term, for which we use Doob's inequality and Itô's isometry. This, together with the previously derived bounds, Young's inequality, and finally applying Gronwall's inequality gives |∇V (x t )| 2 +εn|y| 2 and C 2 = 4n|γ −1 | 2 F L 2 V T .