The Wang-Landau Algorithm as Stochastic Optimization and Its Acceleration
Chenguang Dai Jun S. Liu
stat.CO cs.LG math.OC







chenguangdai@g.harvard.edu

jliu@stat.harvard.edu
Department of Statistics, Harvard University




We show that the Wang-Landau algorithm can be formulated as a stochastic gradient descent algorithm minimizing a smooth and convex objective function, of which the gradient is estimated using Markov chain Monte Carlo iterations. The optimization formulation provides us a new way to establish the convergence rate of the Wang-Landau algorithm, by exploiting the fact that almost surely, the density estimates (on the logarithmic scale) remain in a compact set, upon which the objective function is strongly convex.
The optimization viewpoint motivates us to improve the efficiency of the  Wang-Landau algorithm using popular tools including the momentum method and the adaptive learning rate method. We demonstrate the accelerated Wang-Landau algorithm on a two-dimensional Ising model and a two-dimensional ten-state Potts model. 


Valid PACS appear here
                             

                              
The Wang-Landau Algorithm as Stochastic Optimization and Its Acceleration
    Jun S. Liu
    December 30, 2023
=========================================================================






§ INTRODUCTION


The Wang-Landau (WL) algorithm <cit.> has been proven useful in solving a wide range of computational problems in statistical physics, including spin-glass models <cit.>, fluid phase equilibria <cit.>, polymers <cit.>, lattice gauge theory <cit.>, protein folding <cit.>, free energy profile <cit.>, and numerical integration <cit.>. Its successful applications in statistics have also been documented <cit.>. The WL algorithm directly targets the density of states (the number of all possible configurations for an energy level of a system), thus allowing us to calculate thermodynamic quantities over an arbitrary range of temperature within a single run of the algorithm.

Much effort has been made to understand the dynamics of the WL algorithm, along with numerous proposed improvements, of which we highlight three here. 
(i) Optimizing the modification factor (flatness criterion) <cit.>. Belardinelli and Pereyra <cit.> proposed that instead of reducing the modification factor exponentially, the log modification factor should be scaled down at the rate of 1/t in order to avoid the saturation in the error. (ii) Employing a Parallelization scheme. Wang and Landau <cit.> suggested that multiple random walkers working simultaneously on the same density of states can accelerate the convergence of the WL algorithm. The efficiency of the parallelization scheme can be further enhanced using the replica-exchange framework <cit.>. 
(iii) Incorporating efficient Monte Carlo trial moves <cit.>. 

In this paper, we consider the WL algorithm from an optimization perspective and formulate it as a first-order method. We derive the corresponding smooth and convex objective function, of which the gradient involves the unknown density of states. 
Wang and Landau <cit.> used a random-walk based Metropolis algorithm <cit.>
to estimate the gradient. In general, any suitable Markov chain Monte Carlo (MCMC) strategies <cit.> can be employed for this purpose. Therefore, the WL algorithm is essentially a stochastic gradient descent algorithm. 

The optimization viewpoint enables us to establish the convergence rate of the WL algorithm. Following <cit.> and using the standard stochastic approximation theory <cit.>, we first show that the density estimates (on the logarithmic scale) almost surely stay in a compact set. Based on this, we exploit the strong convexity of the objective function, restricted on this compact set, to prove the convergence rate. We note that the gradient estimator output from the MCMC iterations is generally biased, thus a critical step is to show that the bias vanishes properly as t→∞.

The optimization framework also provides us with a new direction for improving the WL algorithm. We explore one possible improvement, by combining the momentum method <cit.> and the adaptive learning rate method <cit.>. The general goal is to accelerate the transient phase <cit.> of the WL algorithm before it enters the fine local convergence regime. The effectiveness of the acceleration method is demonstrated on a two-dimensional Ising model and a two-dimensional ten-state Potts model, in which the learning in the transient phase is considerably demanding.

The rest of the paper is organized as follows. Section <ref> discusses the optimization formulation of the WL algorithm, and establishes the convergence rate from an optimization perspective.  Section <ref>  introduces possible strategies to accelerate the WL algorithm using optimization tools. Section <ref> demonstrates the accelerated WL algorithm on two benchmark examples. Finally, Section <ref>  concludes with a few remarks. 






§ AN OPTIMIZATION FORMULATION


Let the space of all microscopic configurations be 𝖷. Suppose there are totally N energy levels, E_1 < ⋯ < E_N, for the underlying physical model. For a microscopic configuration x∈𝖷, we use E(x) to denote its energy. Let {g(E_n)}_n = 1^N be the normalized density of states, i.e.,

    g(E_n) ∝#{x∈𝖷, E(x) = E_n},   ∑_n = 1^N g(E_n) = 1.

After initializing g_0(E_n) as 1/N, the WL algorithm iterates between the following two steps: (i) Propose a transition configuration and accept it with probability min{1, g_t(E_i)/g_t(E_j)}, where E_i and E_j refer to the energy levels before and after this transition, respectively. This is essentially a step of the Metropolis algorithm <cit.> with the corresponding stationary distribution: 

    π_t(x) ∝∑_n = 1^N 1/g_t(E_n)1(E(x) = E_n).

(ii) Update the density of states. If E(x_t + 1) = E_n, multiply g_t(E_n) by a modification factor f_t + 1 > 1. That is, g_t + 1(E_n) ← g_t(E_n) × f_t + 1. 

The modification factor f_t should be properly scaled down in order to guarantee the convergence of the algorithm. There is a rich literature on how to adapt  f_t online, including the flat/minimum histogram criterion,
and the 1/t rule <cit.> with its various extensions <cit.>. Under a proper scaling rule, the magnitude of the modification factor f_t is informative of the estimation error <cit.>. Thus, a commonly used stopping criteria for the WL algorithm is that f_t is small enough (say, below exp(10^-8)).  

In the following, we will work on the logarithmic scale of the density of states. Denote u_n^(t) = log(g_t(E_n)) for n ∈ [N], and let u = (u_1,⋯, u_N). 
The density update in the WL algorithm can be rewritten as

    u_n^(t + 1)← u_n^(t) + η_t + 11(E(x_t + 1) = E_n),

where η_t + 1 = log f_t + 1, which will be referred to as the learning rate henceforth. The intermediate target distribution π_t(x) defined in Equation (<ref>) can also be formulated in terms of u^(t). We define

    π_u(x) ∝∑_n = 1^Nexp(-u_n)1(E(x) = E_n),

and denote P_u as a general transition kernel invariant to π_u(x).
For notational convenience, we use π_t(x) to refer to  π_u^(t)(x), and use P_t to refer to the transition kernel invariant to π_t(x). After each density update, we normalize u^(t) to sum to 0, i.e.,
u_n^(t)← u_n^(t) - ∑_i = 1^N u_i^(t)/N,
so that u^(t) stays in a compact set (see Proposition <ref>).
The WL algorithm can be slightly rephrased as in Algorithm <ref>.


Let us consider the following optimization problem:

    min_u∈R^N h(u) = log(∑_n = 1^N  exp(u^⋆_n - u_n)),
       subject to  ∑_n = 1^N u_n = 0,

in which
u_n^⋆ = log(g(E_n)) - 1/N∑_i = 1^Nlog(g(E_i)).
We write u^⋆ = (u^⋆_1, ⋯, u^⋆_N).
It is not difficult to see that this is a convex optimization problem because the objective function h(u) is a log-sum-exp function and the constraint is linear.
It has a unique solution at u_n = u^⋆_n for n∈[N], in which exp(u^⋆_n) equals to the density of states g(E_n) up to an multiplicative constant.

The projected gradient descent algorithm is a standard approach to solve the constrained optimization problem (<ref>). The gradient of the objective function h(u) is

    ∂ h(u)/∂ u_n 
    = -exp(u_n^⋆ - u_n)/∑_i = 1^N exp(u_i^⋆ - u_i),   n∈[N],

which is not directly available because it involves the unknown density of states. However, one can think of approximating the gradient function defined in Equation (<ref>) by one-step or multiple-step Monte Carlo simulations, leading to a stochastic version of the projected gradient descent algorithm.

More precisely, a gradient descent step for minimizing h(u) takes the following form:

    u_n^(t + 1)← u_n^(t) +  η_t + 1exp(u^⋆_n - u_n^(t))/∑_i = 1^N exp(u^⋆_i - u_i^(t)).

Denote the probability of the set {x∈𝖷: E(x) = E_n} with respect to π_t(x) as π_t(E_n). Since the probability π_t(E_n) is proportional to exp(u^⋆_n - u_n^(t)), the density update in Equation (<ref>) is essentially

    u^(t + 1)_n ← u_n^(t) + η_t + 1π_t(E_n).

A crude approximation to π_t(E_n) is the indicator function 1(E(x_t + 1) = E_i), given that after several steps of Monte Carlo simulations according to the transition kernel P_t invariant to π_t(x), x_t + 1 is approximately a sample from π_t(x). This corresponds to the density update in Equation (<ref>).

We note that the projection step to the set Π = {u∈R^N, ∑_n = 1^N u_n = 0} is equivalent to the normalization step (see Algorithm <ref> step 2(c)).
Thus, we have shown that the stochastic projected gradient descent algorithm solving the constrained optimization problem (<ref>), which estimates the probability π_t(E_n) by 1(E(x_t + 1) = E_n) using the output from Monte Carlo simulations, is equivalent to the WL algorithm.

The above optimization formulation has the following immediate implications. 
First, the parallel WL algorithm estimates the negative gradient π_t(E_n) by 1/m∑_k = 1^m[1(E(x^(k)_t) = E_n)], in which m denotes the total number of random walkers, and x^(k)_t denotes the kth random walker. Therefore, it reduces the variance of the gradient estimate by a factor m. Second, instead of implementing a single transition step, the separation strategy mentioned in <cit.> implements multiple transition steps within each iteration, so that the law of the random walker gets closer to the intermediate target distribution π_t(x) defined in Equation (<ref>). Therefore, it reduces the bias of the gradient estimate. 

The optimization formulation also points out a new approach to establish the convergence rate of the WL algorithm. We first state a required assumption, which assumes that the transition kernels are (uniformly) geometrically ergodic over the space Π.


There exists a constant ρ∈ (0, 1) such that for all u∈Π, x∈𝖷, k∈N, we have

    sup_u∈Πsup_x∈𝖷||P^k_u(x, ·) - π_u||_TV≤ 2(1 - ρ)^k,

in which for a signed measure μ, the total variation norm is defined as

    ||μ||_TV = sup_|q|≤ 1|∫_𝖷q(x)μ(dx)|.


We note that sufficient conditions for Assumption <ref> exist in the literature (e.g., condition A2 in <cit.>), and relaxation of Assumption <ref> is also possible <cit.>. 
We have the following result.


Under Assumption <ref>, if we scale down the learning rate η_t in the order of O(1/t), the following two statements hold.


  * Almost surely convergence.


  * There exists a compact set 𝒦⊆Π such that for any t ≥ 0, u^(t)∈𝒦 almost surely.

  * P(lim_t→∞u^(t) = u^⋆) = 1.


  * Convergence rate.
There exists a constant C > 0 such that

    E||u^(t) - u^⋆||^2 ≤ C/t.



The proof of Proposition <ref> is given in the Supplemental Material.

The first part of Proposition <ref> follows similarly as <cit.>. The main idea is to rewrite the WL update, including the density update and the normalization step, as 

    u^(t + 1)←u^(t) + η_t + 1r(u^(t)) + η_t + 1(R(x_t + 1) - r(u^(t))),

in which R_n(x)  = 1(E(x) = E_n) - 1/N, and r(u) is the mean-field function defined as

    r(u) = ∫_𝖷R(x)π_u(x)dx = exp(u^⋆ - u)/∑_n = 1^Nexp(u_n^⋆ - u_n) - 1/N.

The proof of the almost-sure convergence concludes by applying the standard stochastic approximation theory (Theorem 2.2 and Theorem 2.3 in <cit.>) after we establish the following two facts. (1) The remainder term η_t + 1(R(x_t + 1) - r(u^(t))) vanishes properly as t→∞. (2) There exists a Lyapunov function V(u) specified below, 

    V(u) = 1/N∑_n = 1^N exp(u^⋆_n - u_n) - 1,

with respect to the mean-field function r(u), such that ⟨∇ V(u), r(u)⟩ < 0, ∀ u≠u^⋆, and ⟨∇ V(u^⋆), r(u^⋆)⟩ = 0.

The second part of Proposition <ref> is our main theoretical contribution. There are two essential ingredients in establishing the convergence rate. (i) Strong convexity. The objective function h(u) is only convex but not strongly convex on R^N. However, because u^(t) stays in a compact set 𝒦⊆Π almost surely (see Proposition <ref>, part 1(a)), we are able to establish the strong convexity of h(u) restricted on this compact set 𝒦.


Under Assumption <ref>, there exists a constant ℓ > 0 such that for any t ≥ 0, almost surely, it holds

    ⟨∇ h(u^(t)), u^(t) - u^⋆⟩≥ℓ ||u^(t) - u^⋆||^2.


(ii) Vanishing bias. Because x_t + 1 is only an approximate sample from the intermediate target distribution π_t(x), the indicator 1(E(x_t + 1) = E_n) is not an unbiased estimator to the negative gradient π_t(E_n). The following Lemma <ref> shows that the bias of the gradient estimator vanishes properly, as fast as the learning rate, when t→∞.


Under Assumption <ref>, there exists a constant C > 0 such that 

    E||π_t - P_t(x_t, ·)||_TV≤ Cη_t + 1.



The convergence rate of the WL algorithm has been established in different forms in the literature. Zhou and Bhatt <cit.> show that the discrete probability distribution {π_t(E_n)}_n = 1^N will be attracted, in terms of the KL-divergence, to the vicinity of the uniform distribution (π_∞(E_n) = 1/N) as t→∞. In addition, they show that the standard deviation of exp(u_n^⋆ - u_n^(t)) roughly scales like √(log f_t) when the modification factor f_t is close to 1. Although we are looking at the L^2 error of u^(t), which is slightly different from the aforementioned standard deviation, their convergence rate is consistent with our result because √(log f_t)  = √(η_t) is in the order of O(1/√(t)) if we scale down the learning rate η_t in the order of O(1/t).
It is also worthwhile to mention that a corresponding central limit theorem in the original density space is provided in <cit.>. 







§ ACCELERATING WANG-LANDAU ALGORITHM


The optimization formulation motivates us to further improve the WL algorithm using optimization tools <cit.>. Our goal in this paper is to accelerate the convergence in the transient phase. The transient phase <cit.> generally refers to the initial stage of running a stochastic gradient descent algorithm. For instance, if we scale down the learning rate according to the flat/minimum histogram criterion, we can refer to the transient phase as the running period from the beginning up to the time when the flat/minimum histogram criterion is first satisfied.

When the transient phase appears noticeable, the acceleration tools can be very effective in practice, and have been widely used in large-scale systems such as deep neural networks <cit.>. In this paper, we restrict ourselves on the first-order acceleration methods, and leave other possibilities for future explorations. In particular, we find that both the momentum method and the adaptive learning rate method are effective in accelerating the WL algorithm. Before we go into details, we note that improvement in the asymptotic convergence rate of the stochastic gradient descent algorithm is hard to achieve (or even impossible) <cit.> except for some well-structured objective functions such as finite sums. 

The momentum method exponentially accumulates a momentum vector, denoted as m_t in the following, to amplify the persistent gradient across iterations. 
The basic momentum update operates as follows:

    m^(t)←βm^(t - 1) + η_t + 1∇ h(u^(t)), 
       u^(t + 1)←u^(t) - m^(t),

where we initialize the momentum vector to be m^(0) = 0.
We note that the momentum update essentially adds a fraction β of the previously accumulated gradients m^(t - 1) into the current update vector m^(t). The weighting factor β is a tuning parameter, and is commonly set to be 0.9 or higher. 

In the setting of the WL algorithm, the momentum update in Equation (<ref>) becomes

    m_n^(t)←β m_n^(t - 1) - η_t + 11(E(x_t + 1) = E_n), 
        u_n^(t + 1)← u_n^(t) - m_n^(t),     ∀ n∈[N].

The intuition behind the momentum acceleration for the WL algorithm can be heuristically described as follows. The event E(x_t + 1) = E_n suggests that π_t(E_n) is likely larger than 1/N, thus the Markov kernel P_t has a better chance to transit the microscopic configuration x_t into the energy level E_n. Therefore, in order to push π_t(E_n) towards 1/N, that is, downweight the probability mass in the energy level E_n, we increase u_n^(t) by η_t + 1, which corresponds to the density update in Equation (<ref>). 
In contrast to the WL algorithm, which only increases u^(t)_n by η_t + 1 at the current iteration t, we keep increasing u^(t)_n for a few more iterations by an exponentially decay momentum m^(t)_n to achieve a faster convergence. 

The adaptive learning rate method helps standardize the gradient across different coordinates of the parameter u, so that they scale in a similar magnitude. Otherwise, it can be challenging to find a suitable global learning rate η_t over different coordinates.
Popular algorithms along this research direction include AdaGrad <cit.>, AdaDelta <cit.>, and RMSprop (an unpublished method proposed by Geoffrey Hinton). The RMSprop update operates as follows:

    G^(t)←γG^(t - 1) + (1 - γ)∇ h(u^(t))^2, 
       u^(t + 1)←u^(t) - η_t + 1[G^(t)]^-1/2∇ h(u^(t)),

in which both the square and the square root are taken elementwise. G^(t) represents the moving average of the squared gradients, so that the current gradient ∇ h(u^(t)), standardized by [G^(t)]^1/2, is in a similar magnitude across different coordinates. The weighting factor γ is a tuning parameter, which is commonly set to be 0.9 in order to prevent the updates from diminishing too fast. 
In the setting of the WL algorithm, the RMSprop update in Equation (<ref>) becomes

    G^(t)_n ←γ G_n^(t - 1) + (1 - γ)1(E(x_t + 1) = E_n), 
        u_n^(t + 1)← u_n^(t) - η_t + 1[G^(t)_n]^-1/21(E(x_t + 1) = E_n).


The combination of the momentum method and the adaptive learning rate method leads to the Adaptive Moment Estimation (Adam) method <cit.>. The Adam update operates as follows:

    m^(t)←βm^(t - 1) + (1 - β) ∇ h(u^(t)), 
       G^(t)←γG^(t - 1) + (1 - γ)∇ h(u^(t))^2, 
       u^(t + 1)←u^(t) - η_t + 1[G^(t)]^-1/2m^(t).

In the setting of the WL algorithm, we note that,  although β and γ can be potentially two tuning parameters, if we set β = γ and initialize m^(0) and G^(0) to be 0, we have G^(t) = - m^(t), since -∇ h(u^(t)) is approximated by a one-hot vector, which contains only a single “1" with the remaining elements being 0. This simplification leads to Algorithm <ref>, which we refer to as the AWL algorithm henceforth.


We remark that for large-scale systems, a naive implementation of Equation (<ref>) can be very inefficient, as we have to loop over every coordinate of m^(t) and u^(t) in each iteration. A simple solution is to introduce a vector s = (s_1, ⋯, s_N), in which s_n records the last time when m_n and u_n are updated. 
With the help of s_n, instead of updating m_n and u_n in each iteration, we shall update them only when the energy level E_n is involved in the Monte Carlo simulations.




§ ILLUSTRATIONS


We compare the AWL algorithm with the original WL algorithm on two benchmark examples: (a) a nearest-neighbour Ising model; (b) a nearest-neighbour ten-state Potts model. Both models are defined on a two-dimensional L× L square lattice equipped with the periodic boundary condition. 

For the Ising model, the energy E(x) is given by the Hamiltonian:

    E(x) = -∑_<i,j>J_ijx_ix_j - ψ∑_jb_jx_j,

where x_i ∈{±1}. The subscripts i,j denote the lattice sites, and the notation <i,j> implies that the site i and the site j are nearest neighbors. For the ten-state Potts model, the energy E(x) is given by:

    E(x) = -∑_<i,j>J_ij1(x_i = x_j) - ψ∑_jb_jx_j,

where x_i∈{1,⋯, 10}. For both models, we assume that J_ij≡ 1 and b_j≡0 (no external magnetic field). If b_j≡0, the two-dimensional Ising model exhibits a second-order phase transition. Otherwise, in the presence of an external magnetic field, the two-dimensional Ising model exhibits a first-order phase transition. When b_j≡0, the two-dimensional Potts model exhibits a first-order phase transition when the number of states is larger than 4.

Let {H_t(E_n)}_n = 1^N be the histogram of all energy levels at iteration t. We initialize H_0(E_n) = 0 for n∈[N]. At each iteration t, the AWL algorithm and the WL algorithm update u^(t) according to Algorithm <ref> and Algorithm <ref>, respectively. In addition, we update the energy histogram as H_t(E_n) = H_t - 1(E_n) + 1(E(x_t + 1) = E_n). 

The adaptation of the learning rate η_t follows <cit.>, which is detailed in the following.


  * After every 1,000 MC sweeps, we check {H_t(E_n)}. If min_n H_t(E_n) > 0, we set η_t + 1 = η_t/2, and reset H_t(E_n) = 0 for each energy level E_n. Otherwise if min_n H_t(E_n) = 0, we keep η_t + 1 = η_t.

  * If η_t + 1≤ N/t, then η_t = N/t for all the subsequent iterations. H_t(E_n) is discarded and the above step is not executed any more.

We note that each MC sweep contains L^2 iterations, in which each iteration refers to a single round of parameter update. That is, step 2(a)–2(c) in Algorithm <ref> and Algorithm <ref>.
The energy histogram {H_t(E_n)} essentially represents the number of visits to each energy level up to iteration t, since the last update of the learning rate. 

We implement one step of the Metropolis algorithm to estimate the gradient, i.e., step 2(a) in Algorithm <ref> and Algorithm <ref>. The proposal schemes for the Ising model and the Potts model are described as follows. Given the current configuration x_t, we randomly pick up a site and change its value. For the Ising model, we filp its sign. For the ten-state Potts model, we set it to be a number uniformly sampled from {1, ⋯, 10}.

To illustrate the efficiency of the AWL algorithm, we investigate the following four perspectives. (i) The scaling of the first equilibration time, in terms of the number of MC sweeps, with respect to the dimension L. The first equilibration time, which corresponds to the transient phase as we discussed in Section <ref>, is defined to be min{t: min_n H_t(E_n)} > 0. That is, the first time when the energy histogram becomes nonzero everywhere. According to the adaptation rule of the learning rate η_t, the equilibration time is also the first time we decrease the learning rate. (ii) The scaling of the first equilibration time, in terms of the CPU time, with respect to the dimension L. Because the AWL algorithm requires additional computations in updating the momentum vector, the comparison between the two algorithms on the actual CPU time is necessary to see whether the implementation of the acceleration method is indeed worthwhile. (iii) The dynamics of the estimation error ϵ(t) defined as below following <cit.> for L = 80,

    ϵ(t) = 1/N - 1∑_n = 1^N |1 - log(g_t(E_n))/log(g(E_n))|.

For the Ising model, the exact density of states g(E_n) is available, and can be calculated using a publicly available Mathematica program <cit.>. For the Potts model, no exact solution of g(E_n) is available, thus we pre-run a 1/t WL simulation for 5×10^7 MC sweeps, in which the final learning rate is 2×10^-8. We then treat the density estimates as an approximation to the exact density of states. 
(iv) The accuracy in the task of estimating the specific heat for the Ising model with L = 80.

We compare the AWL algorithm and the WL algorithm with different initializations of the learning rate, η_0 = 0.05, 0.10 and 1.00. We test out the two algorithms for different sizes of the two-dimensional square lattice, L = 50, 60, 70, 80, 90, 100. 
The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

Figure <ref> summarizes the computational overheads of the two algorithms for the Ising model. The reported results are based on 50 independent runs of both algorithms, in which the dot represents the empirical mean and the error bar represents the empirical standard deviation. We see that the AWL algorithm takes significantly fewer MC sweeps as well as less CPU time to reach the first equilibration among all settings with different lattice sizes and different initializations of the learning rate. 

Figure <ref> summarizes the computational overheads of the two algorithms on the Potts model. Similar to the case of Ising model, the AWL algorithm is more efficient than the WL algorithm in terms of the first equilibration time measured by the number of MC sweeps and the CPU time. 


Figure <ref> shows the empirical dynamics of ϵ(t), averaged over 50 independent runs of both algorithms. The first 100×10^3 MC sweeps for the Ising model and the first 1500× 10^3 MC sweeps for the Potts model are representative for the transient phase. We see that in the transient phase, the convergence speed of the AWL algorithm, in terms of the number of MC sweeps,  is significantly faster than the convergence speed of the WL algorithm with different initializations of the learning rate. 


For the Ising model with L = 80, Table <ref> compares the accuracy of the two algorithms in the calculation of the specific heat defined as:

    C(T) = ⟨ E^2⟩_T - ⟨ E⟩_T^2/T^2,

in which T denotes the temperature. We test out temperatures ranging from 0.4 to 8 incremented by 0.1. The internal energy ⟨ E⟩_T is defined as 

    ⟨ E⟩_T     = ∑_nE_ng(E_n)exp(-E_n/T)/∑_ng(E_n)exp(-E_n/T).

The fluctuation expression ⟨ E^2⟩_T is defined similarly. We note that the theoretical value of the specific heat at a given temperature T can be evaluated exactly when the exact density of states is available, which is the case for the two-dimensional Ising model. We independently run each algorithm 50 times to obtain 50 independent estimates of the specific heat at each temperature. The relative error at each temperature is calculated based on the mean of the 50 independent estimates. Table <ref> summarizes the quantiles of the relative errors for T∈[0.4, 8], by running each algorithm for 100×10^3, 150×10^3, and 200×10^3 MC sweeps, respectively. Compared to the WL algorithm, the AWL algorithm yields significantly more accurate estimates of the specific heat especially in the transient phase.
 


More details of this numerical study can be found in the Supplemental Material. First, within the first 2×10^5 MC sweeps and 2×10^6 MC sweeps for the Ising model and the Potts model, respectively, we report the number of equilibrations that the AWL algorithm and the WL algorithm have reached (equivalently, the number of changes of the learning rate η_t), for different lattice sizes L and different initializations of the learning rate η_0. We also report the corresponding first 8 equilibration time in terms of the number of MC sweeps. Second, for the Ising model with L = 80, we provide a graphical comparison of the estimated specific heat obtained by the AWL algorithm and the WL algorithm, over the temperature region T∈[0.4, 8].







§ CONCLUSION


To summarize, in this paper we present a new interpretation of the WL algorithm from the optimization perspective. We show that the WL algorithm is essentially a stochastic (projected) gradient descent algorithm minimizing a smooth and convex function, in which MCMC steps are used to estimate the unknown gradient. 
The optimization formulation intuitively explains that because of using more accurate gradient estimates, some notable modifications of the algorithm, such as utilizing multiple random walkers, can improve the WL algorithm. In addition, using the (strong) convexity of the objective function, we provide a new approach to establish the convergence rate of the WL algorithm, which is more explicit compared to the existing results <cit.>. We expect that our contributions are useful for further theoretical investigations of the WL algorithm.

The optimization interpretation also opens a new way to improve the efficiency of the WL algorithm. There are rich tools in the optimization literature to accelerate the stochastic gradient descent algorithm, including but not restricted to the methods we mentioned in Section <ref>. Different methods can be favorable for different applications. In the presence of noisy gradients, it usually requires some careful tuning to successfully apply the acceleration tools. We demonstrate  one possible acceleration approach, using the momentum method and the adaptive learning rate strategy, on a two-dimensional Ising model and a two-dimensional ten-state Potts model.

Acknowledgments. This work was partially supported by NSF DMS-1613035
and DMS-1712714.





unsrt

