A Variational Approach to Monte Carlo Renormalization Group
Yantao Wu Roberto Car
cond-mat.stat-mech











^1The Department of Physics, Princeton University

^2The Department of Chemistry, Princeton University 





             


We present a Monte Carlo method for computing the renormalized coupling constants and the critical exponents within renormalization theory. 
The scheme, which derives from a variational principle, overcomes critical slowing down, by means of a bias potential that renders the coarse grained variables uncorrelated. The 2D Ising model is used to illustrate the method. 













Valid PACS appear here
                             

                              
A Variational Approach to Monte Carlo Renormalization Group
    Yantao Wu^1 and Roberto Car^1, 2
    December 30, 2023
===========================================================






Since the introduction of renormalization group (RG) theory <cit.>, there has been strong interest in methods to compute the renormalized coupling constants and the critical exponents in a non-perturbative fashion. 
This goal has been achieved with the Monte Carlo (MC) RG approach of Swendsen. 
In 1979, he introduced a method to compute the critical exponents, which did not require explicit knowledge of the renormalized Hamiltonian <cit.>. 
A few years later, he solved the problem of calculating the renormalized coupling constants, using an equality due to Callen <cit.> to write the correlation functions in a form explicitly depending on the couplings. 
By imposing that the standard MC expression of a correlation function and its corresponding Callen form be equal, he derived equations whose iterative solution led to the coupling constants <cit.>. 
Finding the renormalized Hamiltonian is an example of inverse statistical mechanical problem <cit.>. 
MCRG has been used successfully in many applications but difficulties related to sampling efficiency may be severe. 
Typically, the evaluation of the correlation functions near a critical point suffers from critical slowing down and is affected by large sampling errors in large systems. 
This difficulty can be alleviated with ingenious cluster algorithms <cit.>, which, however, are limited to specific models.   

Here we present an MCRG framework based on a variational principle for a biasing potential acting on the coarse grained degrees of freedom of a RG transformation. 
In our approach, the coupling constants and the critical exponents derive from the same unifying principle. 
Swendsen's formulae emerge as a special case, but our scheme also leads to formulations exempt from critical slowing down. 
In addition, it permits to estimate variationally the effect of truncating the Hamiltonian.  

Although the approach is rather general, here we limit ourselves, for concreteness, to lattice models with discrete spin degrees of freedom, {σ}. 
A generic Hamiltonian has the form  

    H(σ) = ∑_α K_α S_α(σ),

where the K_α are coupling constants and the S_α are operators acting on the spins σ, such as sums or products of spins or combinations thereof.  

RG considers a flux in the space of Hamiltonians (<ref>) under scale transformations that reduce the linear size of the original lattice by a factor b. 
The rescaled degrees of freedom take the same discrete values of the original spins, to which they are related by a coarse graining transformation, σ' = τ(σ). 
For example, τ can be the block spin transformation of Kadanoff <cit.>.   

The distribution of the σ' is obtained from the distribution of the σ by tracing out the original degrees of freedom while keeping the σ' fixed:  

    p(σ') = ∑_σδ_τ(σ), σ' e^-H(σ)/Z =  e^-H'(σ')/Z' .

Here δ is the discrete Kroneker-delta function, 
Z and Z' are partition functions that ensure the normalization of the corresponding distributions. 
While the partition function Z' is invariant under RG transformations, the renormalized Hamiltonian H' is not, except at fixed points of the RG flow:    

    Z = ∑_σ e^-H(σ) = ∑_σ' e^-H(σ') = Z'

and 

    H'(σ') = -log∑_σδ_τ(σ), σ' e^-H(σ)

Repeated at infinitum, the RG transformations generate a flux in the space of Hamiltonians, in which all possible coupling terms appear, unless forbidden by symmetry. 
For example, in an Ising model with no magnetic field, only even spin products appear. 
The space of the coupling terms is, in general, infinite. 
However, perturbative and non-perturbative calculations suggest that only a finite number of couplings should be sufficient for a given degree of accuracy. 

In the proximity of a critical point, the distribution (<ref>) of the block spins σ' displays a divergent correlation length, originating critical slowing down of local MC updates. 
This can be avoided by modifying the distribution of the σ' by adding to the Hamiltonian H'(σ') a biasing potential V(σ') to force the biased distribution of the block spins, p_V(σ'), to be equal to a chosen target distribution, p_t(σ'). For instance, p_t can be the constant probability distribution.  
Then the σ' have the same probability at each lattice site and act as uncorrelated spins, even in the vicinity of a critical point.           

It turns out that V(σ') obeys a powerful variational principle that facilitates the sampling of the Landau free energy <cit.>. 

In the present context, we define the functional Ω[V] of the biasing potential V(σ') by: 

    Ω [V] = log∑_σ' e^-[H'(σ') + V(σ')]/∑_σ' e^-H'(σ')+ ∑_σ' p_t (σ') V(σ'),

where p_t(σ') is a normalized known target probability distribution. As demonstrated in <cit.>, the following properties hold:  


  * Ω [V] is a convex functional with a lower bound. 

  * The minimizer, V_min(σ'), of Ω is unique up to a constant and is such that: 

    H'(σ') = - V_min(σ') - log p_t (σ') + constant


  * The probability distribution of the σ' under the action of V_min is:   

    p_V_min(σ') = e^-(H'(σ') + V_min(σ'))/∑_σ' e^-(H'(σ') + V_min (σ')) =  p_t(σ')


The above three properties lead to the following MCRG scheme.  

First, we approximate V(σ') with V_J⃗(σ'), a linear combination of a finite number of terms S_α(σ') with unknown coefficients J_α, forming a vector J⃗ = {J_1, ..., J_α, ..., J_n}.    

    V_J⃗(σ') = ∑_α J_α S_α(σ')

Then the functional Ω[V] becomes a convex function of J⃗, due to the linearity of the expansion, and the minimizing vector, J⃗_min, and the corresponding V_min(σ') can be found with a local minimization algorithm using the gradient and the Hessian of Ω: 

    ∂Ω(J⃗)/∂ J_α = - ⟨S_α(σ')|_⟩V_J⃗ + ⟨S_α(σ')|_⟩p_t


    ∂^2 Ω(J⃗)/∂ J_α∂ J_β = ⟨S_α(σ') S_β(σ')|_⟩V_J⃗ - ⟨S_α(σ')|_⟩V_J⃗⟨S_β(σ')|_⟩V_J⃗

Here ⟨·|_⟩V_J⃗ is the biased ensemble average under V_J⃗ and ⟨·|_⟩p_t is the ensemble average under the target probability distribution p_t. 
The first average is associated to the Boltzmann factor exp{-(H'(σ') + V(σ'))} = ∑_σδ_τ(σ), σ'exp(-H(σ)) exp(-V(τ(σ))) and can be computed with MC sampling. 
The second average can be computed analytically if p_t is simple enough.      
⟨·|_⟩V_J⃗ always has inherent random noise, or even inaccuracy, and some sophistication is required in the optimization problem. 
Following <cit.>, we adopt the stochastic optimization procedure of <cit.>, and improve the statistics by running independent MC simulations, called multiple walkers, in parallel. 
For further details, consult <cit.> and the  Supplementary Material (SM) <cit.>. 

The renormalized Hamiltonian H'(σ') is given by Eq. <ref> in terms of V_min(σ'). Taking a constant p_t, we have modulo a constant: 

    H'(σ') = -V_min(σ') = ∑_α (-J_min, α) S_α(σ')

In this finite approximation the renormalized Hamiltonian has exactly the same terms of V_min(σ') with renormalized coupling constants 

    K'_α = -J_min, α.

The relative importance of an operator S_α in the renormalized Hamiltonian can be estimated variationally in terms of the relative magnitude of the coefficient J_min, α. 
When J_min, α is much smaller than the other components of J⃗_min, the corresponding S_α(σ') is comparably unimportant and can be ignored. The accuracy of this approximation could be quantified by measuring the deviation of p_V_min(σ') from p_t(σ').  

To illustrate the method, we present a study of the Ising model on a 2D square lattice in the absence of a magnetic field. 
We adopt 3 × 3 block spins with the majority rule. 
26 coupling terms were chosen initially, including 13 two-spin and 13 four-spin products. 
One preliminary iteration of variational RG (VRG) was performed on a 45× 45 lattice starting from the nearest-neighbor Hamiltonian.
The coupling terms with renormalized coupling constants smaller than 0.001 in absolute value were deemed unimportant and dropped from further calculations. 
13 coupling terms, including 7 two-spin and 6 four-spin products, survived this criterion and were kept in all subsequent calculations <cit.>. 
Each calculation consisted of 5 VRG iterations starting with nearest-neighbor coupling, K_nn, only. 
All the subsequent iterations used the same lattice of the initial iteration. 
Standard Metropolis MC sampling <cit.> was adopted, and the calculations were done at least twice to ensure that statistical noise did not alter the results significantly. 

In Fig. <ref>, results are shown for a 300 × 300 lattice with two initial K_nn, equal to 0.4355 and to 0.4365, respectively. When K_nn = 0.4365, the renormalized coupling constants increase over the five iterations shown, and would increase more dramatically with further iterations. 
Similarly, they decrease when K_nn = 0.4355. 
Thus, the critical coupling K_c should belong to the window 0.4355-0.4365. 
The same critical window is found for the 45×45, 90× 90, 150× 150, and 210× 210 lattices <cit.>.  
Because each iteration is affected by truncation and finite size errors, less iterations for the same rescaling factor would reduce the error.
For example, 4 VRG iterations with a 2×2 block have the rescaling factor of a 16 × 16 block. 
The latter is computationally more costly than a calculation with 2× 2 blocks, but can still be performed with modest computational resources. Indeed, with a 16 × 16 block, RG iterations on a 128 × 128 lattice gave a critical window 0.4394-0.4398 <cit.>, very close to the exact value, K_c ∼0.4407, due to Onsager <cit.>.

The statistical uncertainty of the renormalized couplings from the variational method is small. 
Using the standard approach, Ref. <cit.> found a renormalized nearest neighbor coupling equal to 0.408 ± 0.002 after the first RG iteration on a 36× 36 lattice using a 3 × 3 block spin, starting with K_nn = 0.4407. 
This result required 5.76 × 10^5 MC sweeps. 
With our method, applied to a 300 × 300 lattice, starting with K_nn = 0.4365, we found a renormalized nearest-neighbor coupling equal to 0.38031 ± 0.00002 after 3.398 × 10^5 MC sweeps.
The standard error in our case was computed with the block averaging method <cit.>. Because <cit.> used only seven coupling terms and a different initial K_nn, the renormalized couplings should not be expected to be the same in the two calculations, but a comparison of the corresponding statistical uncertainties should be meaningful.  


According to theory <cit.>, the critical exponents are obtained from the leading eigenvalues of ∂ K'_α/∂ K_β, the Jacobian matrix of the RG transformation, at a critical fixed point. 
In order to find ∂ K'_α/∂ K_β near a fixed point, we need to know how the renormalized coupling constants K'_α from a RG iteration on the Hamiltonian H = ∑_β K_β S_β, change when K_β is perturbed to K_β + δ K_β, for fixed target probability p_t and operators S_α. 

The minimum condition, Eq. <ref>, 
implies dΩ/d J_α = 0, i.e. for all γ:  

    ∑_σ S_γ(σ') e^- ∑_β (K_β S_β(σ) - K'_β S_β(σ'))/∑_σ e^- ∑_β (K_β S_β(σ) - K'_β S_β(σ'))  = ⟨S_γ(σ')|_⟩p_t,

and  

    ∑_σ S_γ(σ') e^- ∑_β ((K_β + δ K_β) S_β(σ) - (K'_β + δ K'_β) S_β(σ'))/∑_σ e^- ∑_β ((K_β + δ K_β) S_β(σ) - (K'_β + δ K'_β)S_β(σ'))
    
    = ⟨S_γ(σ')|_⟩p_t.

Expanding Eq. <ref> to linear order in δ K'_α and δ K_β, we obtain (<cit.>) 

    A_βγ = ∑_α∂ K'_α/∂ K_β·  B_αγ,

where 

    A_βγ = ⟨S_β(σ) S_γ(σ')|_⟩V - ⟨S_β(σ)|_⟩V⟨S_γ(σ')|_⟩V,

and 

    B_αγ = ⟨S_α(σ') S_γ(σ')|_⟩V - ⟨S_α(σ')|_⟩V⟨S_γ(σ')|_⟩V.

Here ⟨·|_⟩V denotes average under the biased Hamiltonian, H = ∑_β K_β S_β(σ) - K'_β S_β(σ').

If we require the target average of S_γ(σ') to coincide with the unbiased average under H = ∑_β K_β S_β, K' would necessarily vanish and Eqs. <ref>-<ref> would coincide with Swendsen's formulae <cit.>. 

If we use a uniform target probability, the σ' at different sites would be uncorrelated, and critical slowing down would be absent.  

In practice, in order to compute the critical exponents, we first need to locate K_c. 
From the above calculations on the 45 × 45, 90 × 90, and 300 × 300 lattices with a 3 × 3 block spin, we expect that K_c = 0.436 should approximate the critical nearest-neighbor coupling in our model. 
Indeed an RG iteration starting from this value gives couplings that remain essentially constant, as illustrated in Figs. S11-S13 of the SM <cit.>. 

Then, we use Eqs. <ref>-<ref> to compute the Jacobian of the RG transformation by setting K_c = 0.436. 
The renormalized coupling constants after the first RG iteration represent K_α, and those after the second RG iteration represent K'_α. 
The results for biased and unbiased ensembles are shown in Table <ref>, which reports the leading even (e) and odd (o) eigenvalues of ∂ K'_α/∂ K_β when including 13 coupling terms for the three L× L lattices with L = 45, 90, and 300.  
As seen from the table, biased and unbiased calculations give slightly different eigenvalues, as one should expect, given that the respective calculations are different embodiments of the truncated Hamiltonian approximation. 
For L = 300 the results are well converged in the biased ensemble.  
By contrast, we were not able to obtain converged results for this lattice in the unbiased ensemble on the time scale of our simulation.    
The absence of critical slowing down in the biased simulation is demonstrated in Fig. <ref>, which displays  time decay of a correlation function in the biased and unbiased ensembles. See also Figs. S14-S15 of the SM <cit.>. 



The fixed point used for Table <ref> is approximate, and we did not make any effort to fine tune the approximation. 
Refinements could be done iteratively using Eqs. <ref>-<ref>, as we will discuss in a future paper. 
There is an important benefit in knowing accurately the location of the fixed point, because then a single RG iteration, instead of multiple implicit iterations would suffice to compute the Jacobian.  
Moreover, one could use small block spins, having a smaller statistical uncertainty than larger block spins. 

In summary, we have unified the calculation of critical exponents and renormalized couplings within the same framework.
A key feature of our approach is that we adopt a biased ensemble, ⟨·|_⟩V, for the averages. 
This not only simplifies the algorithm, but also enhances the sampling. 
In fact, the original motivation for the variational principle <cit.> was to overcome the long correlation time in first-order phase transitions. 
The bias potential constructed by optimizing the functional acquires a history-dependence that discourages the sampling of previously visited configurations <cit.>, thereby breaking the long correlation time of the unbiased simulation. 
In the RG context, enhanced sampling eliminates critical slowing down.  
We expect that it should be also helpful in systems with deep local free energy minima, as the variational method was originally designed to deal precisely with such systems. 

The finite size of the numerical samples is a source of error.
If the RG iterations are carried out on a single L× L lattice, the coarse grained lattice will have size L/b×L/b. 
Then, as noted in <cit.>, the calculated renormalized couplings will have different size errors on the L × L and L/b×L/b lattices. 
A better way, as suggested in <cit.>, would be to perform calculations on two lattices, L × L and L/b×L/b, so that the coarse grained lattice rescaled by b^n, at the nth iteration starting from L × L, would coincide with the lattice rescaled by b^n-1, at the (n-1)th iteration starting from L/b×L/b. 
In this way, two successive RG iterations have the same lattice size, with a significant cancellation of finite size errors. 
We plan to discuss in a future paper how this idea could be implemented within VRG.  

In the present paper we have used a constant probability distribution p_t, but there is no reason to always do so. 
For example, in systems with continuous and unbounded degrees of freedom, like molecular systems or lattice field theory, it may be convenient to use a Gaussian distribution for p_t. 

Finally, we note that a regular term g(K) always appears as the inhomogeneous part of a RG transformation <cit.>:  

    exp[H'(K'; σ') + Ng(K)] = ∑_σδ_τ(σ), σ'exp[H(K; σ)]

The g(K) in this equation is precisely the thermodynamic free energy per site in the biased ensemble ⟨·|_⟩V, as shown in the SM <cit.>.
It is then interesting, and somewhat surprising, that the information on the critical behavior is fully contained in the statistical behavior of ⟨·|_⟩V, even though g(K) is a regular function and ⟨·|_⟩V does not show singular behavior.           

All the codes used in this project were written in C, and would be available upon request. The authors would like to thank C. Castellani and L. Pietronero for discussions. Partial support for this work was provided by the Department of Energy under Grant no. DE-FG02-05ER46201. 
apsrev

17
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Wilson(1971)]wilson_rg
authorK. G. Wilson,
  journalPhys. Rev. B volume4,
  pages3174 (year1971).

[Swendsen(1979)]mcrg
authorR. H. Swendsen,
  journalPhys. Rev. Lett. volume42,
  pages859 (year1979).

[Callen(1963)]callen
authorH. Callen,
  journalPhysics Letters volume4,
  pages161 (year1963), ISSN issn0031-9163.

[Swendsen(1984a)]mcrg_rc
authorR. H. Swendsen,
  journalPhys. Rev. Lett. volume52,
  pages1165 (year1984a).

[Albert and Swendsen(2014)]invising
authorJ. Albert and
  authorR. H. Swendsen,
  journalPhysics Procedia volume57,
  pages99  (year2014), ISSN issn1875-3892.

[Swendsen and Wang(1987)]clustermc
authorR. H. Swendsen and
  authorJ.-S. Wang,
  journalPhys. Rev. Lett. volume58,
  pages86 (year1987).

[Kadanoff(1966)]block2
authorL. P. Kadanoff,
  journalPhysics volume2, pages263
  (year1966).

[Valsson and Parrinello(2014)]varyfes
authorO. Valsson and
  authorM. Parrinello,
  journalPhys. Rev. Lett. volume113,
  pages090601 (year2014).

[Bach and Moulines(2013)]bach
authorF. Bach and
  authorE. Moulines, in
  booktitleAdvances in Neural Information Processing Systems
  26, edited by editorC. J. C.
  Burges,
  editorL. Bottou,
  editorM. Welling,
  editorZ. Ghahramani,
  and editorK. Q.
  Weinberger (publisherCurran Associates, Inc.,
  year2013), pp. pages773–781.

[sm()]sm
noteSupplementary Material.

[Metropolis et al.(1953)Metropolis,
  Rosenbluth, Rosenbluth, Teller, and Teller]metropolis
authorN. Metropolis,
  authorA. W. Rosenbluth,
  authorM. N. Rosenbluth,
  authorA. H. Teller,
  and authorE. Teller,
  journalThe Journal of Chemical Physics
  volume21, pages1087 (year1953).

[Onsager(1944)]onsager
authorL. Onsager,
  journalPhys. Rev. volume65,
  pages117 (year1944).

[Swendsen(1984b)]mcrg2dising
authorR. H. Swendsen,
  journalPhys. Rev. B volume30,
  pages3866 (year1984b).

[Flyvbjerg and Petersen(1989)]block_method
authorH. Flyvbjerg and
  authorH. G. Petersen,
  journalThe Journal of Chemical Physics
  volume91, pages461 (year1989).

[Wilson(1975)]wilsonkondo
authorK. G. Wilson,
  journalRev. Mod. Phys. volume47,
  pages773 (year1975).

[Landau and Swendsen(1981)]two_lattice
authorD. P. Landau and
  authorR. H. Swendsen,
  journalPhys. Rev. Lett. volume46,
  pages1437 (year1981).

[Nauenberg and Nienhuis(1974)]nauenberg
authorM. Nauenberg and
  authorB. Nienhuis,
  journalPhys. Rev. Lett. volume33,
  pages1598 (year1974).


