Reweighted Interacting Langevin Diffusions: an Accelerated Sampling Method for Optimization

In this section, we introduce the classical Langevin Dynamics which is related to GLD and SGLD algorithms.

Suppose $V:\mathbb{R}^{d}\to\mathbb{R}$ that is twice differentiable and has a simple global minimum ${\bm{x}}^{\ast}=\arg\min_{x\in\mathbb{R}^{d}}V({\bm{x}})$ for simplicity.

The overdamped Langevin Dynamics is defined as in 1. Such a dynamic has plenty of relationships with the Fokker-Planck equation. Denote

the infinitesimal generator (Øksendal, 1987) of the Markov process $({\bm{x}}_{t})_{t\geq 0}$, and its $L^{2}$ adjoint operator:

Let us assume the law of ${\bm{x}}_{t}$ at time $t$ has a density $p_{t}(x)$ for Lebesgure measure. Then $p_{t}({\bm{x}})$ satisfies the Fokker-Planck equation:

We may denote $p_{t}({\bm{x}})=e^{\mathcal{L}^{\dagger}t}p_{0}({\bm{x}})$ where we denote $p_{0}({\bm{x}})$ the initial distribution of overdamped Langevin dynamics, and it is well-known that the Markov operator $e^{\mathcal{L}^{\dagger}t}$ admits a unique invariant probability measure $\nu({\bm{x}})=Z_{\nu}^{-1}e^{-2\sigma^{-2}V({\bm{x}})},Z_{\nu}=\int_{\mathbb{R}^{d}}e^{-2\sigma^{-2}V({\bm{x}})}dx$. The rate converging to $\nu({\bm{x}})$ has been wildly studied,

Note that $\lambda_{1}$ is real because $\mathcal{L}$ is self-adjoint²²2We say a linear operator $A$ is self-adjoint on a Hilbert space $\{\mathcal{H},\langle\cdot\rangle_{\mathcal{H}}\}$, if for any $f,g\in\mathcal{H},\langle f,Ag\rangle_{\mathcal{H}}=\langle Af,g\rangle_{\mathcal{H}}$ on $L^{2}(\nu)$.

When for sampling purposes, individual samples are sampled from independent paths by Monte Carlo Markov Chain (MCMC, (Berg, 2004)) corresponding to the discretized Langevin dynamics. This method is somehow slow for large dimensional, non-convex problems. Such a weakness make it less attractive for optimization purposes. Fortunately, we can explore the connection between individual samples, as they did in Garbuno-Iñigo et al. (2019). We will introduce their work in the next subsection.

In Garbuno-Iñigo et al. (2019), the authors introduced a modified Langevin dynamics and analyzed its property, proving its superiority in designing sampling algorithms. We now introduce the main result of their work.

The convergence rate of classical Langevin dynamic based algorithm can be really slow when $V$ varies highly anisotropic, . A common approach for canceling out this effect is to introduce a $d\times d$ preconditioning positive semi-definite matrix $C$ in the corresponding gradient scheme,

Here $\sqrt{C}$ can be any $d\times r$ real matrix $U$ such that $UU^{T}=C$ (Note that in this case, $W_{t}$ can be reduced to $r$-dimensional standard Brownian motion as the essential rank of $C$ is no larger than r).

The corresponding Fokker-Planck equation now becomes

and the infinitesimal generator

One can easily verify that $e^{-2\sigma^{-2}V({\bm{x}})}$ is the invariant measure to the above system for all semi-definite $C$, and the only one positive invariant measure if $C$ is strictly positive definite.

To find a suitable $C$ is of general interest. One of the best choices is to let $C=\text{Hess}V$, as a counterpart of Newton’s scheme in optimization, which is unfriendly for computation. One intrinsic benefit of choosing $\text{Hess}V$ is its affine-invariant property when designing numerical schemes.

Such a property can also be preserved if we take $C=\mathcal{C}(p)$, the covariance matrix under the probability measure $p({\bm{x}})$,

The dynamic now becomes a nonlinear flow

To simulate from such a mean-field model, the finite interacting particle system $\mathcal{X}_{t}=\{{\bm{x}}^{i}_{t}\}_{i=1}^{N}$ is introduced:

where $\mathcal{B}_{t}^{i}$ are i.i.d. standard Brownian motions, and

Particles in this system are then no longer independent to each other, in contrast to independent simulations in SGLD algorithm.

Now suppose $V$ is a least squares functional³³3For any positive-definite matrix $A$, we define $\langle a,a^{\prime}\rangle_{A}=\langle a,A^{-1}a^{\prime}\rangle=a^{T}A^{-1}a^{\prime}$, and $\|a\|_{A}=\|A^{-\frac{1}{2}}a\|$. with Tikhonov-Phillips regularization (Engl et al., 1996):

where $y\in\mathbb{R}^{k},\mathcal{G}:\mathbb{R}^{d}\to\mathbb{R}^{k}$ is the forward map, $\Gamma$ and $\Gamma_{0}$ are two positive definite matrixs. In this situation, the system (13) can be re-written as

Using the $1^{\text{st}}$ order Taylor approximation

one may approximate Eq. 15 in a derivative-free manner as

Keeping the invariant measure to be exactly $e^{-2\sigma^{-2}V({\bm{x}})}$ is not necessary for optimization. It would be preferable if a new process can converge faster to its invariant measure with its invariant measure still concentrating on the global optimum as $\sigma\to 0$.

Now we introduce an additional source term to modify Eq. (6) into

Here we mainly consider the situation when $C=I$ corresponding to Eq. 4, or $C=\mathcal{C}(p_{t})$ corresponding to Eq. 12, but the analysis remains valid for all positive-definite $C$. We assume in this paper that the function $W$ is smooth and upper bounded.

Let us look inside what the role $W$ plays in the evolution of this process. If we take $W$ as a function to $f$ such that $W({\bm{x}})$ becomes larger when $f({\bm{x}})$ becomes smaller, intuitively the ratio of the mass in the better-fitness region (closer to the global minimum) becomes larger, thus we expect the invariant measure concentrates more on the global optimum.

Note that such a process (17), (18) is no longer a gradient flow structure, which does not preserve total mass, thus a normalization is added to keep the total mass equal to 1. Such a normalizing process has been well-studied as the so-called Feynman-Kac Semigroup when considering the linear case: $C$ is a constant matrix. As the spectral analysis is wildly studied for linear operators (Kato, 1966), and for numerical discretization the $C$ is fixed at each time step, we thus conduct analysis for any fixed matrix $C$.

We introduce the solution operator corresponding to Eq. 17 and Eq. 18: recall the infinite generator $\mathcal{L}^{\dagger}_{C}$, the solution in Eq. 17 can be represented by $\tilde{p}_{t}=e^{t(\mathcal{L}^{\dagger}_{C}+W)}p_{0}$. We denote the corresponding reweighted operator, which is also well-known as Feynman-Kac Semigroup, as $\Phi^{t}_{\mathcal{L}_{C}+W}:$

then $p_{t}=\Phi^{t}_{\mathcal{L}_{C}+W}(p_{0})$ is exactly the solution to Eq. 18.

The operator $\Phi^{t}_{\mathcal{L}_{C}+W}$ is nonlinear in general, but it has many similarity with the linear operator $e^{t(\mathcal{L}^{\dagger}_{C}+W)}$. In specific, the unique positive fixed-point $p_{\infty}$ of $\Phi^{t}_{\mathcal{L}_{C}+W}$ is proportional to the principle eigenfunction of $(\mathcal{L}^{\dagger}_{C}+W)$, and the convergence rate of $\Phi^{t}_{\mathcal{L}_{C}+W}(p_{0})$ to $p_{\infty}$ is controlled by the spectral gap of $(\mathcal{L}_{C}+W)$. For the Feynman-Kac semigroup, one can refer to Ferré & Stoltz (2017) for time-invariant case, Lyu et al. (2021) for time-periodic case, and Moral (2004) for systematical details.

We are interested in what benefits the source term $W$ can contribute. We show that, when taking $W=\varepsilon m(V)$ for some small $\varepsilon>0$ where $m:\mathbb{R}\to\mathbb{R}$ is a monotonic decreasing function of $V$, this brings mainly two benefits:

• •

  the spectral gap to the operator $\mathcal{L}_{C}+W$ is larger than $\mathcal{L}_{C}$ considered in same functional space;

• •

  the invariant measure of $\Phi^{t}_{\mathcal{L}_{C}+W}$ concentrates more on the global minimum compared to the invariant measure of $e^{t(\mathcal{L}^{\dagger}_{C})}$.

These two benefits show that, in the same noise level, process (18) can converge faster and have more concentrated invariant measure than process (4) or (27). We will further explain the benefits, giving proof in Section 4.1.

Now let us use interacting particle methods (Moral & Miclo, 2000; Moral, 2013) and mean field theory (Kac et al., 1960) to design algorithms. We introduce the so-called Reweighted Interacting Langevin Diffusion algorithm, by simply approximating $p_{t}$ with a population of particles, and discretize the evolution of time by operator-splitting technique and forward-Euler scheme.

First, we convert Eq. 17 and 18 into a discrete-time version: let $\tau>0$ be a fixed timestep:

then we use the operator splitting technique⁴⁴4Note that higher order splitting technique can be applied, such as Strange splitting (Strang, 1968) $e^{\frac{\tau}{2}W}\circ e^{\tau\mathcal{L}^{\dagger}_{C}}\circ e^{\frac{\tau}{2}W}$ with splitting error $O(\tau^{3})$. However, the approximation error can still be induced in later steps, thus we only choose a simplest one here.: approximate $e^{\tau(\mathcal{L}^{\dagger}_{C}+W)}$ by $e^{\tau W}\circ e^{\tau\mathcal{L}^{\dagger}_{C}}$ with splitting error $O(\tau^{2})$ .

Let’s now approximate $e^{\tau\mathcal{L}^{\dagger}_{C}}p_{n\tau}$. Suppose $p_{n\tau}$ is approximated by a weighted empirical measure that can be expressed as $\hat{p}_{n\tau}({\bm{x}})=\sum_{i=1}^{N}w^{i}_{n}\delta_{{\bm{x}}_{n}^{i}}({\bm{x}})$ generated from a sort of sample-weight pair $\{{\bm{x}}_{n}^{i},w_{n}^{i}\}_{i=1}^{N}$. We use the simplest Euler-Maruyama (Faniran, 2015) approximation,

Here $C$ can be $I$ or the covariance matrix of current step: if we take the matrix $C$ to be the covariance matrix, it is computed as (denote $\Lambda=\text{diag}(w_{n}^{i}),\quad\bar{\bm{x}}_{n}=\frac{1}{N}\sum_{i=1}^{N}w_{n}^{i}{\bm{x}}_{n}^{i}$)

Then, after normalizing, we get the natural approximation

Note that elements in $\{w^{i}_{n}\}_{i=1}^{N}$ may be polarized, we use a resampling technique for better approximation: If $\frac{\max_{i}w^{i}_{n}}{\min_{i}w^{i}_{n}}$ reaches to a threshold, we resample the replicas $\{\bm{x}_{n}^{i}\}_{i=1}^{N}$ according to the multinomial distribution associated with $\{w_{n}^{i}\}_{i=1}^{N}$, which defines a new set of replicas $\{\bm{\tilde{x}}_{n}^{i}\}_{i=1}^{N}$ and the empirical distribution

Then we replace $\hat{p}_{n\tau}$ by $\tilde{p}_{n\tau}$ in Eq. 21 and conduct further computation.

When $V$ is a least square functional with Tikhonov-Phillips regularization like in Eq. 14, we can follow the same analysis, approximating $\mathcal{C}(\hat{p}_{n\tau})\nabla V({\bm{x}}_{n}^{i})$ in a similar derivative-free manner as in Eq. 16: ( $\mathcal{\bar{G}}_{n}:=\frac{1}{N}\sum_{k=1}^{N}w_{n}^{i}\mathcal{G}(\bm{x}^{k}_{n})$ )

We now conclude in algorithm 1. Note that we fix the population $N$, stepsize $\tau$, and noise level $\sigma$ for simple, but these can be adjusted dynamically for faster convergence.

Note that our methods can also be interpreted as a mutation-selection genetic particle algorithm with MCMC mutations: the update rule for $\bm{x}_{n}^{i}\to\bm{x}_{n+1}^{i}$ can be regarded as mutation, and the resampling step can be regarded as selection. Such a type of algorithm acts like a ridge connecting Genetic Algorithm and PDEs, which can conduct better convergence analysis.

In practice, many problems are hard to get the exact gradients for optimizing: gradient information is infeasible, or computationally expensive. We thus suggest a gradient-free variant, which corresponds to a process only with diffusion and source term. We will prove that such a process can exponentially converge to its invariant measure, and it also has an invariant measure that generally concentrates on the global minimum as $\sigma\to 0$.

We introduce the formula of the modified process as

Note that we only delete the term related to $\nabla V$, the term $W$ that can still be chosen to relate to $V$. We will state in Thm. 4.6 that, the system will finally converge to a distribution concentrating on the maximum of $W$.

The discrete algorithm is designed similarly as in Algorithm 1, but just ignores the gradient term $C\nabla V$. This can be seen by taking $V(\bm{x})\equiv\text{const}$.

Now we analyze how $W$ helps in improving the convergence rate, as well as the sharpness of the invariant measure.

First let us restrict the comparison in the $L^{2}(\nu)$, where $\nu(x)=e^{-\frac{2V(x)}{\sigma^{2}}}$. The reason to choose this space is mainly because the following property,

To continue our analysis, we need to make following assumption for $V$ and $W$, which is necessary for $L_{C}+W$ and $L_{C}$ have a discrete and up-bounded spectrum (see Pankov (2001)).

Next, we need to prove that the Feynman-Kac semigroup $\Phi^{t}_{\mathcal{L}_{C}+W}$ does map any initial density $p_{0}\in L^{2}(\nu)$ to a limit density. Such a result is concluded as follows:

Next, let us analyze how the convergence speed of $\Phi^{t}_{\mathcal{L}_{C}+W}p_{0}$ is improved compared to the original process $e^{t\mathcal{L}_{C}}(p_{0})$, which is heavily related to the spectral gap. To exactly analyze the spectral gap to a differential operator is hard in general. Many existing analysis of spectral gap only consider the simplest case: $\nabla V$ is a constant matrix. In contrast to these existing techniques, we analyze the contribution of $W$ to the spectral gap by perturbation theory: we will prove that the new process has a better convergence rate compared to the old one, if $W=\varepsilon m(V)$ for a small $\varepsilon>0$, where $m:\mathbb{R}\to\mathbb{R}$ is a monotonic decreasing function to $V$.

In the gradient-free situation, the original operator $\mathcal{D}_{C}:=\text{div}\big{(}\frac{\sigma^{2}}{2}C\nabla\cdot\big{)}$ has a trivial invariant measure: uniform distribution, thus no optimization property can be expected when $W$ is not included. We thus turn our results to another way, showing that $\Phi_{\mathcal{D}_{C}+W}^{t}p_{0}$ exponentially converges to a distribution that concentrates on the neighborhood of the global minimum, and as $\sigma\to 0$, the invariant distribution gets more and more concentrate on the global minimum. We now state the result as follows:

Let us consider a 1 $d$ periodic problem. We first verify our results in Thm. 4.5. Let $x\in[0,1)$⁵⁵5Although our analysis considers the space $\mathbb{R}$, the results can be easily transferred to the periodic case with Asmp. 4.2 be removed., we use Fourier Spectral method (Shen et al., 2011) to discretize the differential operator $\mathcal{L}$ and $\mathcal{L}-\varepsilon V$,where

for any periodic $f\in C^{\infty}([0,1))$, and

with boundary smoothly modified.

In Fig. 1, We plot the graphs of the principal eigenfunction $\nu$ to the operator $\mathcal{L}$, the principal eigenfunction $\mu$ to the operator $\mathcal{L}-0.1V$, function $V$ and $\frac{\mu}{\nu}$ respectively. As we expected, $\mu$ concentrates more on the global minimum of $V$ than $\nu$. In Fig. 1, we can see $\lambda_{0}(\varepsilon)-\lambda_{1}(\varepsilon)$ increasing as $\varepsilon$ increasing, showing the enhancement of spectral gap. This experiment gives a visual explanation to Thm. 4.5: specifically, the eigenfunction to $\mathcal{L}$ is the Gibbs measure $\nu(\bm{x})=\exp(-2V(\bm{x})/\sigma^{2})$, getting larger when $V$ gets smaller; the eigenfunction $\mu$ to $\mathcal{L}-0.1V$ has no doubt more mass than $\nu$ that concentrating nearby the global minimum of $V$, and the staircase-like graph of the quotient of $\mu/\nu$ gives a direct explanation to (42).

Next, let us verify our result in Thm. 4.6. We use the same space $x\in[0,1)$ and the source term $W(\bm{x})=-V(\bm{x})$ where $V(x)$ is the same as in (30), and test the analytical property of the following operator

for any periodic $f\in C^{\infty}([0,1))$. We plot in Fig. 2 the principal eigenfunction, or invariant distribution density $\mu_{\sigma}(x)$ of $\mathcal{D}_{\sigma}$ with different $\sigma$. We also calculate the mass nearby the global minimum: we take the interval $I=[0.44,0.68]$, and calculate $\int_{I}\mu_{\sigma}(x)dx$. The results is plotted in Fig. 2, showing the concentrating tendency when $\sigma\to 0$. This coincides with the result in Thm. 4.6: the invariant measure concentrates more and more on the global maximum of $W$ as $\sigma\to 0$.

In this section, we design some numerical tests in inverse problem solving fields. We compare our RILD algorithm with EKS (Garbuno-Iñigo et al., 2019) and EKI (Kovachki & Stuart, 2018). These problems can be seen as an optimization problem of a least square functional with Tickhonov-Phillips regularization (14), thus a derivative-free approximation to the updating rule ((16) for EKS and EKI, or (26) for RILD) can be used to design derivative-free schemes.

We first try to solve a low-dimensional inverse problem. The numerical experiment considered here is the example originally presented by Ernst et al. (2015), and also used in Herty & Visconti (2018). We compare with the result from Garbuno-Iñigo et al. (2019), and the experimental settings are exactly the same. The forward map is given by the solution of a one-dimensional elliptic boundary value problem as defined in Garbuno-Iñigo et al. (2019),

with $f(0)=0,f(1)=x_{2}$. The explicit solution is given by

Thus we define the forward map

Here $\bm{x}=(x_{1},x_{2})^{T}$ is a constant vector we want to solve, and we have noisy measurements $\bm{y}=(27.5,79.7)^{T}$ of $f(\cdot)$ at locations $u_{1}=0.25,u_{2}=0.75$. This can be solved by minimizing the least-square functional defined as in Eq. 14. We assume observation noise $\Gamma=0.1^{2}I_{2}$, prior matrix $\Gamma_{0}=10^{2}I_{2}$, and initial ensemble drawn from $N(0,1)\times U(90,110).$ The ensemble size is $N=10^{3}$. We fix $\sigma=\sqrt{2}$. The stepsize $\tau$ is updated adaptively as in Garbuno-Iñigo et al. (2019). We take $W(\bm{x})=-\|\mathcal{G}(\bm{x})-\bm{y}\|^{2}_{\Gamma}$.

We compare our RILD algorithm 1 with EKS and EKI algorithms. The key difference between our RILD and EKS in this situation is the use of reweighting/resampling. The results are plotted in Fig. 3, 3, and 3. From the figure, we can see that our RILD algorithm converges much faster than EKI or EKS algorithms using the same super-parameters, and as our problem is settled as a posterior sampling problem, the ensemble of RILD stops shrinking to the minimum point after some iterations, unless one decrease the diffusion parameter $\sigma$. One could expect, using a decay schedule to $\sigma$, RILD algorithm will perform much better than EKS or EKI algorithms for optimizations.

We also tested in a high-dimensional case. Specifically, we define the map $\mathcal{G}:\mathbb{R}^{d}\to\mathbb{R}^{2d-2}$

$\bm{y}=(0,\cdots,0,1,\cdots,1)^{T}$ with $0$ repeats $d-1$ times and $1$ repeats $d-1$ times. One can verify that

is exactly a $Rosenbrock$ function. We choose $d=100$, observation noise matrix $\Gamma=0.1^{2}I_{198}$, prior distribution matrix $\Gamma_{0}=10^{2}I_{100}$, and initial ensemble drawn from $N(2,0.3^{2}I_{100})$. The global solution of $\mathcal{G}(\bm{x})=\bm{y}$ is $\bm{x}=(1,\cdots,1)^{T}.$ The ensemble size is fixed to $N=10^{3}$. We take $\sigma=\sqrt{2}$ and the stepsize $\tau$ is updated adaptively as before. For RILD algorithm, we choose $W(\bm{x})=-5*10^{-3}\|\mathcal{G}(\bm{x})-\bm{y}\|^{2}$. Similar to the test before, one can find RILD converges much faster than EKI or EKS in these high-dimensional sampling tasks, and thus perform better in optimization situations.

We now test our RILD algorithm in a highly nonconvex high-dimensional situation. we test with 100 dimensional $Ackley$ function, which is defined as follows:

where $\bm{x}=(x_{1},\cdots,x_{d})^{T},d=100,a=20,b=0.2,c=2\pi$. As the difficulties mainly raise from the numerious local minimum, covariance modification that is designed for ill-posed problems is not suitable here, we just take $C=I$.

We compare our algorithm RILD in Alg. 1 with the classical Gradient Langevin Dynamics (GLD) algorithm. GLD algorithm discretizes a single path of the Langevin Dynamics:

The difference between RILD and GLD is: RILD maintains an ensemble of size $N$ while GLD only maintains 1 individual, and at each step, RILD calculates a weight associated with each individual, then resamples the ensemble according to the weight, see Alg. 1. Now we test if RILD has better ability getting out of local minimums, compared to GLD.

For RILD, we take ensemble size $N=50$, and randomly pick the initial ensemble⁶⁶6Such an initial setting creates many difficulties to find the global minimum, as the first term in the $Ackley$ function becomes quickly dominated when $\bm{x}$ gets away from the origin. from $N(0,30^{2}I_{100})$. For GLD, the initial point is randomly chosen from the initial ensemble of RILD. We test a wide range of the stepsize $\tau\in[2,32]$, $\sigma\in[1,16]$, and for each fixed $\tau$ and $\sigma$, we repeat 10 trials to calculate the pass rate: we say one trial is passed, if the RILD or GLD algorithm can find a point $\bm{x}$ that $V(\bm{x})<17$ in $5*10^{4}$ evaluations⁷⁷7If one finds a point smaller than 17, the remaining task will be trivial as the first term in the $Ackley$ function gets dominant.. All trials in all super-parameter settings begin with the same initial ensemble. In 5 One can see that RILD has a wide range of super-parameter settings to find the true decay directions, while GLD algorithm cannot make it in any tested settings. We also tested GA and PSO under the same initial condition with different super-parameters, and reported the best searching result in Fig. 5 together with RILD and GLD.