Ensemble Monte Carlo Calculations with Five Novel Moves

To move walker $i$ from location $\vec{r}_{i}$ to $\vec{r}_{i}^{\prime}$, Goodman & Weare (2010) employed the affine invariant moves that were originally introduced by Christen (2007). From the ensemble, one selects one additonal walker, $j$, to set the stretch direction,

To make such moves reversible, the stretch factor, $\lambda$, must be sampled from the interval $\left[\frac{1}{a},a\right]$ where $a>1$ is constant parameter that we controls the step size. To sample the stretch factor, $\lambda$, one has a bit of a choice. Goodman & Weare (2010) followed Christen (2007) when they adopted a function, $P_{S}(\lambda)$, that satisfies,

The acceptance probability is given by

where the factor $\lambda^{\alpha}$ with $\alpha=N-1$ emerges because one performs a one-dimensional move to sample points in a $N$-dimensions space (see derivations in Christen (2007) and Militzer (2023a)).

In Fig. 1, we illustrate the affine as well as the quadratic moves. For a quadratic move, Militzer (2023a) selects two walkers $j$ and $k$ from the ensemble to move walker $i$ from $\vec{r}_{i}$ to $\vec{r}_{i}^{\prime}$.

The interpolation weights $w$ are chosen from,

where $L_{2}$ is the typical Lagrange weighting function that guarantees a proper quadratic interpolation so that $\vec{r}_{i}^{\prime}=\vec{r}_{i}$ if $t_{i}^{\prime}=t_{i}$; $\vec{r}_{i}^{\prime}=\vec{r}_{j}$ if $t_{i}^{\prime}=t_{j}$; and $\vec{r}_{i}^{\prime}=\vec{r}_{k}$ if $t_{i}^{\prime}=t_{k}$. One sets $t_{j}=-1$ and $t_{k}=+1$ to introduce a scale into the parameter space, $t$.

To satisfy the detailed balance condition, $T(\vec{r}_{i}\to\vec{r}_{i}^{\prime})=T(\vec{r}_{i}^{\prime}\to\vec{r}_{i})$, it is key that we sample the parameters $t_{i}$ and $t_{i}^{\prime}$ from the same distribution $\mathcal{P_{S}}(t)$. The acceptance probability then becomes,

Again a factor of $\left|w_{i}\right|^{N}$ is needed because we sample the one-dimensional $t$ space but then switch to the $N$-dimensional parameter space, $\vec{r}$. In appendix of Militzer (2023a), this factor is derived rigorously from the generalized detailed balance equation by Green & Mira (2001).

Goodman & Weare (2010) also introduced walk moves. To move walker $k$ from $\vec{r}_{k}$ to $\vec{r}_{k}^{\prime}=\vec{r}_{k}+W$, one chooses at random a subset of walkers, $S$. The subset size, $N_{S}$, is a free parameter that one needs to choose within $2\leq N_{S}<N_{W}$ where $N_{W}$ is the total number of walkers. Militzer (2023a) followed Goodman & Weare (2010) in computing the average location all walkers in the subset,

but then modifoed their formula for computing the step size, $W$, by introduding a scaling factor $a$:

$Z_{j}$ are univariate standard normal random numbers. By setting $a=1$, one obtains the original walk moves, for which the covariance of the step size, $W$, is the same as the covariance of subset $S$. The introdcution of the scaling parameter, $a$, enabled Militzer (2023a) to make smaller (or larger) steps in situations where the covariance of the instantaneous walker distribution is a not an optimal representation of local structure of the sampling function. It was demonstrated that the scaling factor $a$ significantly improves the sampling efficiency of the Rosenbrock function and for the ring potential in high dimensions.

For a quadratic MC move one selects two guide points to move walker $i$ from $\vec{r}_{i}$ to the new location $\vec{r}^{\prime}_{i}$. The new location is generated by quadratic interpolation in parameter space $t$. Here we now generalize this concept to arbitrary order, $N_{O}$. At random, we first select $N_{O}$ guide points from the ensemble of walkers. Together with the moving walker, we now have $N_{O}+1$ points to perform a Lagrange interpolation of order $N_{O}$ to derive a new location for the moving walker, $\vec{r}_{i}^{\prime}$,

The interpolation weights, $w$, are given by

where $\vec{t}$ is an $N_{O}+1$ dimensional vector $\vec{t}=\{t_{0},t_{1}\ldots t_{N}\}$ that begins with index 0. $L_{N_{O}}$ are the Lagrange polynomials,

where the modular division in $t_{(i+s)\%(N_{O}+1)}$ guarantees that index $i+s$ referes to vector element $t_{i+s}$ if $i+s<N_{O}+1$ and to element $t_{i+s-(N_{O}+1)}$ if $i+s$ is larger. This leads to the usual Lagrange interpolation, $\vec{r}(t^{\prime})$, that goes through all $N_{O}+1$ points so that $\vec{r}_{i}^{\prime}=\vec{r}_{i}$ if $t^{\prime}=t_{0}$ and $\vec{r}_{i}^{\prime}=\vec{r}_{j}$ if $t^{\prime}=t_{j}$. For the quadratic moves in Militzer (2023a), we set $t_{j=1}=-1$ and $t_{j=2}=+1$ to introduce a scale into the parameter space $t$. The $t$ arguments for the original and new locations, $t_{0}$ and $t^{\prime}$, are both sampled at random.

Militzer (2023a) proposed two alternative methods and sampled them either a uniform distribution from $[-a,+a]$ or from a Gaussian function with standard deviation $\sigma=a$ that is centered at $t=0$. For the acceptance probabilities, we use again Eq. 12. For both sampling functions, the parameter $a$ controls the average size of the MC steps in $t$ and $\vec{r}$ spaces. While the standard deviation of resulting $t$ values is equal to $a$ for the Gaussian $t$ sampling method, it is equal to only $a/\sqrt{3}$ for the linear $t$ sampling. So the favorable ranges of $a$ tend to be a little larger for linear $t$ sampling than for the Gaussian method when the best settings of $a$ for both methods are compared for the same application, as we will later see.

When we extend the quadratic moves to higher orders here, we still sample $t_{j=0}$ and $t^{\prime}$ as before but we need to specify the arguments, $t_{j>0}$, for the remaining interpolation points. For a third order interpolation, a natural choice would be $t_{1}=-1$, $t_{2}=0$, and $t_{3}=+1$. Following this concept, we distribute the $t$ arguments uniformly,

so that $t_{j=1}=-1$ and $t_{j=N_{O}}=+1$. We study how this MC method performs when we choose the order $N_{O}=\{3,4,6,10\}$ and combined it with linear and Gaussian $t$ sampling and different values of $a=\{0.1,0.3,0.5,1.0,1.2,1.5,2.0,3.0\}$. While we see a dependence on $a$ when we compare the performance of this method, the performance of the linear $t$ sampling method with the best choice $a$ was always comparable to that of the Gaussian $t$ sampling method with its best $a$. So for simplicity we combine the results with linear and Gaussian $t$ sampling into one dataset when we later compare the performance this method for different orders, $N_{O}$.

None of the moves in Goodman & Weare (2010) or Militzer (2023a) nor any other moves in this article make use of the probability values of the other walkers, $\pi(\vec{r}_{j})$, when a new location, $\vec{r}_{i}^{\prime}$, for walker $i$ is proposed. Currently only their locations, $\vec{r}_{j}$, are employed. Some valuable information might be harnessed by utilizing that a particular walker, $k$, resides in an unlikely location with a probability value, $\pi(\vec{r}_{k})$, that much lower than that of the others. (One may recall that the regula falsi method is more efficient in finding the roots of a function than the bisection method because it makes use of the function value while the bisection method only relies on the sign of the values.)

Here we go back the quadratic moves that relies on the walker, $\vec{r}(t_{0})$, and two guide points $\vec{r}(t_{1}=-1)$ and $\vec{r}(t_{2}=+1)$. We assume their probabilities, $\pi(\vec{r}(t^{\prime}))$, came from an energy function, $E(\vec{r}(t^{\prime}))=-k_{B}T*\log(\pi(\vec{r}(t^{\prime})))$. The factor $k_{B}T$ may be set to 1. We further assume that we can approximate the energy function by the quadratic function, $E(t^{\prime})=At^{\prime 2}+Bt^{\prime}+C$, whose coefficients are chosen so that it interpolates the three known points, $E(t_{0})$, $E(t_{1})$, and $E(t_{2})$. If $A$ is positive, it has a minimum at $t_{\rm min}=-B/2A$ and we can sample $t^{\prime}$ from a Gaussian function,

with the standard deviation, $\sigma=\left(2A/k_{B}T\right)^{-1/2}$.

For the reverse move, one first needs to consider the probability of sampling $t^{\prime}$ from the distribution, $P_{S}(t^{\prime})$. To sample $t$, one constructs a different quadratic function, $P^{\prime}_{G}(t)$, that interpolates $E(t^{\prime})$, $E(t_{1})$, and $E(t_{2})$. This introduces an additional factor into the acceptance ratio,

The approach has one caveat because the interpolation coefficient, $A$, may occasionally become negative during the forward or reverse move, in which case the distribution $P_{G}(t^{\prime})$ cannot be properly normalized. When this happens if first try to reorder the elements of the vector $\vec{t}$. If this is not successful we overwrite $A$ with $|A|$, which leads to a stable MC algorithm. We tested it with with linear and Gaussian sampling, $P_{S}(t)$ and a series of $a$ values, $a=\{0.1,0.3,0.5,1.0,1.2,1.5,2.0,3.0\}$.

Here we introduce three extensions of the affine and quadratic moves. First for our modified affine moves, we linearly interpolate between walkers $j$ and $k$ to pick a direction for the stretch move of walker, $i$,

where $c$ is a random number chosen uniformly between 0 and 1. In all other respects, the move proceeds like the original affine move. Illustrations of this and the two following moves are given in Fig. 1.

For our affine simplex moves, we set the direction of the stretch move in yet a different way. As guide points, we select $N_{G}$ walkers from the ensemble, compute their center of mass, and insert it as $\vec{r}_{*}$ into Eq. 21. Finally for our quadratic simplex move, we select two separate sets of $N_{G}$ walkers at random and compute their respective centers of mass before inserting them as of $\vec{r}_{j}$ and $\vec{r}_{j}$ into Eq. 7. The move then proceeds like our quadratic moves.

We test all our MC methods by applying them to two test problems: the ring potential and Rosenbrock density. The ring potential is defined as,

where $\vec{r}=\{r_{1},\ldots,r_{N}\}$ is a vector in the $N\geq 2$ dimensions. $\rho=\sqrt{r_{1}^{2}+r_{2}^{2}}$ is the distance from the origin in the $r_{1}$-$r_{2}$ plane. The first term ensures that the potential is small along a ring of constant radius, $R=1$. The second term keeps the magnitude of all remaining parameters, $r_{3\ldots N}$, small. Increasing the positive integer, $m$, makes the walls of the potential around the ring steeper. For this article, we performed calculations with $N=12$ and 24 dimensions while keeping $m=6$ fixed. The last term breaks axial symmetry. We set $C$ to small value of 0.01 so that the potential minimum is approximately located at point $\vec{A}=(+R,0,\ldots)$ while the potential is raised at opposing point $\vec{B}=(-R,0,\ldots)$. Figure 3 of Militzer (2023a) shows an illustration of this effect. The prefactor of the first term in Eq. 22 is introduced so that the location of potential minimum does not shift much with increasing $m$.

We employ the Boltzmann factor, Eq. 1, to convert the energy of the ring potential into a probability function that we can sample. We set $k_{B}T=10^{-4}$ so that the equilibrium distribution is reasonably well confined around point $\vec{A}$. So when we compute the autocorrelation time and the blocked error bar (Allen & Tildesley, 1987) of the potential energy, we initialized the ensemble of walkers around this point. When we want to determine how long it takes the ensemble of walkers to travel around the ring, we initialize the ensemble near the high-energy point, $\vec{B}$. We define the travel time to be the number of MC moves that are required for ensemble average of $\left<r_{1}\right>$ to change sign from negative to positive. Based on such a travel time analysis, Militzer (2023a) concluded that employing between $2N+1$ and $3N+1$ walkers was a reasonable ensemble size while it took larger ensembles a disproportionally long time to travel around the ring.

The ring potential also allows to define and study the cohesion among an ensemble of walkers, which is relevant because all MC methods in this article perform local moves that rely on information of walkers in the ensemble. This means that if the walkers have split up in a highly dimensional and fragmented fitness landscape, or it were initialized that way, it will be challenging to reunited the walkers later. This will likely to lead to a degradation of performance because a distant walker cannot provide useful guidance for a local move. Therefore, we consider cohesion among the walkers to be a favorable trait of an ensemble MC algorithm.

For the ring potential, we define cohesion in the following, simple way. As for the travel time calculation, we initialize the ensemble near the high-energy state of point $\vec{B}$ and then monitor whether all walkers find to their way to the low-energy state near point $\vec{A}$. Cohesion is defined by the fracton of walkers with $r_{1}\geq 0$ at the end of the MC calculation. This average could in principle depend on the duration of the MC calculations but Fig. 2 provides an illustration why this dependence is weak. Conservatively, we end our MC calculations after they completed enough moves equal to twice the ring travel time.

To demonstrate why it is informative to compare the walker cohesion between different MC algorithms, we plot the $r_{1}$ coordinates of all walkers at the end of every block in Fig. 2. For the ring potential in $N=12$ dimensions, we initialized two ensembles of with $N_{W}=24$ walkers near $r_{1}=-1$ and propagated them with our QMC and with the affine methods. The QMC ensemble remained together and traveled more efficiently towards the low-energy near $r_{1}=+1$. The affine ensemble traveled more slowly and split up so that only 17 out of 24 walkers reached the low-energy state. A careful inspection of the first 2000 blocks reveals that even some walkers in the QMC ensemble fell behind but eventually they all caught up. Still one may not expect the QMC ensemble to show perfect cohesion in all circumstances especially if sampling parameters are chose poorly. While Fig. 2 shows just one calculation, we always average the cohesion and the travel time over 1000 independent MC calculations when we compare prediction from different MC moves and parameter settings in this article.

Goodman & Weare (2010) tested their methods by sampling the 2d Rosenbrock density,

which carves a narrow curved channel into the $(r_{1},r_{2})$ landscape. $B$ effectively plays the role of temperature while $A$ controls the width of the channel. Here we set $A=100$ and $B=5$ to be consistent with Goodman & Weare (2010). There are three different ways one can generalize the Rosenbrock density to higher dimensions. First one can simply treat it as $N/2$ independent 2d Rosenbrock problems by setting,

Alternatively one can “connect” the individual coordinates by defining,

Finally one can make the problem periodic by also connecting the coordinates $r_{1}$ and $r_{N}$,

In Fig. 3, we compare error bars and autocorrelation times that we computed with all the different MC moves and parameters that we described in Sec. 2 for these three types of Rosenbrock densities in $N=10$ dimensions. We find that the simple Rosenbrock density to be the most challenging one to sample because it leads to the largest error bars and the longest autocorrelation times. Since we are interesting here in designing algorithms that can solve challenging sampling problems, we focus all following investigations on the simple Rosenbrock density and exclude the two others from further consideration. The simple Rosenbrock density was also studied Huijser et al. (2022) who demonstrated that the affine invariant method exhibits undesirable properties in more than 50 dimensions.

Our code is open source. Examples and installation instruction are available here: Militzer (2024). A simpler source code for our quadratic Monte Carlo method is available here: Militzer (2023b).

In Figs. 4 and 5, we respectively compare the results from different MC sampling methods for the Rosenbrock function in $N=2$ and 20 dimensions. In Fig. 6, we show similar results for the ring potential in $N=24$. We also computed results for $N=12$ dimensions but we do not provide a plot here because they are very similar, with one exception: The order 3 results did not lag behind the QMC results in $N=12$ case as they do for $N=24$ in Fig. 6. One sees the same trend if one compares the $N=2$ and 20 results for the Rosenbrock function.

For any application of MC methods, one wants the autocorrelation time and computed averages as small as possible. We compute both for the energy because this is the central quantity of many physical system. The results from different sampling methods are compared in panels (a)-(d) in Figs. 4, 5, and Fig. 6. For the ring potential, an efficient algorithm should fulfill two additional criteria: It should travel efficiently from unfavorable to favorable regions of the parameter space and exhibit a high level of cohesion. Both are plotted in panels (e)-(h) of Fig. 6.

The QMC method with linear and Gaussian $t$ sampling yield very good results across all panels in Figs. 4, 5, and 6 but there is also quite a bit of variance among the predictions. So with the following discussion, we identify favorable choice for the different sampling parameters and compare for the two applications with goal of identifying some reasonable default settings for future applications.

For Rosenbrock function in $N=2$ dimensions, the QMC method with linear $t$ sampling works well if the scaling parameter $a$ is set between 1.0 and 3.0 for any number of walkers between 3 and 7. If $a$ is chosen poorly between 0.1 and 0.5, the energy error bar and autocorrelation time increase drastically, which is the trend we see in Fig. 4a. Similarly, the QMC method with Gaussian $t$ sampling works well for $a$ value between 0.5 and 3.0 and not so well for $a=0.1$ and 0.3.

If we increase the number of dimensions from 2 to 20, the favorable parameter ranges are slightly modified. The QMC method with linear $t$ sampling works well if $a$ is set only between 0.3 and 1.0 but not so well for larger $a$ values, nor for $a=0.1$. The QMC method with Gaussian $t$ sampling works well for $a$ values between 0.1 and 0.5 but not so well for larger $a$ values. As one increases the dimensionality of the Rosenbrock density, sampling the search space becomes more challenging (Huijser et al., 2022) and one needs to make smaller steps to still move efficiently.

When we compare the performance of the QMC method with linear $t$ sampling for the ring potential in $N=24$ dimension, we find the method to yield small error bars, short autocorrelation and travel times as well as a high level of cohesion of over 90% if $a$ is set between 0.3 and 0.5 for between 26 and 73 walkers.

For $a=1.0$, the method works well for between 27 and 49 walkers. If we switch to Gaussian $t$ sampling, $a$ values between 0.1 and 0.5 yield favorable results. The QMC simplex method yields results that are very similar to that of the original QMC method as Figs. 4, 5, and 6 show consistently. This also means that these test cases did not reveal any particular advantage of this method.

We find that the affine invariant method does not perform as well as the QMC method for all three test cases. It lags behind by factors between 2 and 4 behind the best QMC results as panels (b) of Figs. 4, 5, and 6 show. We obtained the best results for the ring potential by setting $a=1.2$ and using between 28 and 73 walkers. For these simulations, the cohesion level varied between 80 and 97% (see Fig. 6f). For the Rosenbrock density in 2 and 20 dimensions, we needed to set $a$ to 2.5 and 1.2 respectively to obtain the best results. This confirms the trend that sampling the $N=20$ space efficiently requires smaller steps.

Our modification to the affine method did not yield any improvements. In Figs. 5 and 6, it slightly lags behind the original affine method while Figs. 4 shows the autocorrelation time for the 2d Rosenbrock density is about twice as long. The affine simplex performed yet a bit worth. It lagged behind the other affine methods in term of the travel time and cohesion as Fig. 6f shows.

We obtained the best results for the ring potential in $N=24$ dimensions in Fig. 6 when we employed the QMC method with linear $t$ sampling with $a=0.3\ldots 0.5$ and Gaussian $t$ sampling with $a=0.1\ldots 0.3$. Using between 27 and 73 walkers gave consistently good results. In comparison, the quadratic simplex method did not quite perform as well for any sampling parameters.

Affine and modified affine methods did slightly worse than the quadratic simplex method. The autocorrelation time, error bars, and travel times were not as good and the cohesion varied substantially between 50% and 95% (see Fig. 6b and f). The best choice for the $a$ value was 1.2. The affine simplex method performed poorly on all counts.

In Fig. 6c and g, we compare results from interpolations with different orders. Some simulations of fourth and sixth order performed well if we set $a$ to 0.1 or 0.3. In comparison, calculations of third and tenth order were not competitive.

In Fig. 6d and h, we compare results obtained with the directed quadratic method and with walk moves. With the directed quadratic method, we obtained the best results by setting $a$ to 0.3 and 0.5 but they were not as good as those from original quadratic method. On the other hand, the results from the walk method were consistent better as long as we set the scaling parameter $a$ to 0.1 or 0.5, which underlines the importance this parameter that was introduced only recently by Militzer (2023a).

The preceeding analysis relied on a detailed comparision of results spread across three figures. So here we introduce the relative inverse efficiency to greatly simplify the comparison of different methods across various applications. Such an efficiency measure should be able to handle the noise that we see in Figs. 4, 5, and 6. Furthermore, one cannot expect that adjustable parameters like $a$ and $N_{W}$ have been optimized to perfection. Instead it should be sufficient to have results for good number of reasonable $a$ and $N_{W}$ values in the dataset. Also a method should not be penalized very much if in addition to results for reasonable parameters there are also some findings for unreasonable values like $a=3$ in the dataset. Finally if a method like the modified walk method has more adjustable parameters and more calculations have been performed, this should not lead to an advantage over other methods that have fewer adjustable parameters. With these tree goals in mind, we define the relative inverse efficiency to be the ratio, $A/B$. For given application (e.g. the Rosenbrock density in 2 dimensions), we sort all autocorrelation times obtained with all methods and all parameters by magnitude and then compute the average time of the best $f=20\%$ results to obtain $B$. (This value would not change if a completely useless method were to be added to the dataset.) We derive $A$ by averaging the best 20% result derived with one particular method only. Employing the fraction $f$ of the results rather than single data points helps reduce the effect of the noise and it favors methods to yield good results over a wider range of parameters, which is helpful for real challenging application for which one may not be able to perform many test calculations. The ratio, $A/B$, is thus a measure how a particular method compares to the best of its peers. In Fig. 7, we plot a relative inverse efficiency for twelve MC methods and four applications. We not only computed it for the autocorrelation time but also for the squared error bar of the energy and for the ring potential, also for the ring travel time and inverse of the cohesion so that small values of all four measures imply a high level of efficiency. (Setting $f=10\%$ did not yield any meaningful change.)

Fig. 7 illustrates that the QMC method with linear and Gaussian $t$ sampling has the highest efficiency overall, which is consistently above average. The performance of the QMC simplex method is also very good but it does not do quite as well for the ring potential. Fig. 7 also illustrates that the original affine method nor its modified or simplex forms can compete with the QMC method. The efficieny is lower by factors between 2 and 10.

Furthermore Fig. 7 shows that increasing the interpolation order from 2 (QMC) to 10 leads to a gradual decrease in efficiency but this trend has one exceptions. The fourth-order method does notable better than the third-order method, which only does well for the Rosenbrock density in 2 dimensions.

The directed QMC method does very well but, for the four applications presented here, it does not offer any improvement over the original QMC method. Finally the modified walk method does markedly better than the QMC moves for the ring potential according to all four efficiency criteria but it does not sample the Rosenbrock density very well.