Comment on "Solving Statistical Mechanics Using VANs": Introducing saVANt - VANs Enhanced by Importance and MCMC Sampling
Kim Nicoli Pan Kessel Nils Strodthoff Wojciech Samek Klaus-Robert Müller Shinichi Nakajima





Machine Learning Group, Technische Universität Berlin, 10587 Berlin, Germany


Machine Learning Group, Technische Universität Berlin, 10587 Berlin, Germany
Berlin Big Data Center, 10587 Berlin, Germany


Fraunhofer Heinrich Hertz Institute, 10587 Berlin, Germany


wojciech.samek@hhi.fraunhofer.de
Fraunhofer Heinrich Hertz Institute, 10587 Berlin, Germany
Berlin Big Data Center, 10587 Berlin, Germany
Berliner Zentrum für Maschinelles Lernen, 10587 Berlin, Germany


klaus-robert.mueller@tu-berlin.de
Machine Learning Group, Technische Universität Berlin, 10587 Berlin, Germany
Berlin Big Data Center, 10587 Berlin, Germany
Department of Brain and Cognitive Engineering, Korea University, Anam-dong, Seongbuk-gu, Seoul 136-713, South Korea
Max-Planck-Institut für Informatik, Saarbrücken, Germany
Berliner Zentrum für Maschinelles Lernen, 10587 Berlin, Germany


Machine Learning Group, Technische Universität Berlin, 10587 Berlin, Germany
Berlin Big Data Center, 10587 Berlin, Germany
RIKEN Center for AIP, 103-0027, Tokyo, Japan







05.10.Ln, 
05.20.−y, 
02.50.Ng, 
11.15.Ha 


Comment on “Solving Statistical Mechanics Using VANs”: 
 Introducing saVANt—VANs Enhanced by Importance and MCMC Sampling
    Shinichi Nakajima
    Received ; accepted 
=========================================================================================================================




 
In statistical physics and strongly-interacting field theories, physical observables are often computed using Monte-Carlo methods; most interestingly near criticality. 
In a recent exciting publication <cit.>, 
the authors have demonstrated the promising potential of generative neural samplers (GNSs)
in this context. 
In this comment, we propose a subtle yet crucial modification which enhances their approach drastically by giving theoretical guarantees and fast computation. We envision our method to be helpful for broad applications in high-energy and condensed matter physics.

The Variational Autoregressive Network (VAN)  <cit.>
is a GNS, successfully trained so as to provide independent samples from an approximate distribution q(s) of the Boltzmann distribution p(s)= 1/Z  exp(-β H(s)).
As the authors showed, thermodynamic observables (expectation values of an operator 𝒪) can be estimated by the sample mean ⟨𝒪(s) ⟩_p≈1/N∑_i𝒪(s_i) of the samples  drawn from a VAN s_i ∼ q.




Our proposal, the sampled VAN training (saVANt) is motivated by an exceptional feature of VANs, namely that the exact sampling probability q(s) is accessible unlike in most of the popular GNSs <cit.>. 
Taking advantage of this fact, we posit that the sampling error can be corrected by using neural network-based MCMC or importance sampling which leads to asymptotically unbiased estimators for physical quantities.

Practically, we train saVANt samplers using the same approach as Wu et al. but with output probabilities bounded within [ϵ, 1 - ϵ] for small ϵ >0. Then, Neural Importance Sampling (saVANt-NIS) defined below gives an unbiased estimator:

    ⟨𝒪(s) ⟩_p≈∑_i w_i  𝒪(s_i)       with        s_i ∼ q   ,

where w_i = ŵ_i∑_i ŵ_i for ŵ_i = e^-β H(s_i)q(s_i)
is the importance weight.
Alternatively, Neural MCMC Sampling (saVANt-NMCMC) 
with the acceptance probability

    min(1, p_0(s'| s)   p(s')p_0(s| s')   p(s))
        = min(1, q(s')  exp(-β H(s'))q(s)  exp(-β H(s)))

also provides an unbiased estimator.
Note that, in saVANt-NMCMC, we can draw samples from the sampler independently, and apply the Metropolis rejection step to the samples arranged in an arbitrary order.
This is
because the trial distribution p_0(s'| s) = q(s') is independent from the previous sample s.  

The asymptotic unbiasedness is guaranteed even if q does not accurately approximate the target distribution p (see e.g. <cit.>),
which allows transfer across parameter space. For example,
a sampler trained at a certain temperature can be used to estimate physical observables at other temperatures.

We demonstrate the effectiveness of our proposal for the Ising model with H(s)=- ∑_<i, j> s_i s_j. We trained all samplers on a Tesla P100 GPU by using the implementation provided by Wu et al <cit.>. The error bars were determined using the method described in <cit.>, and all runs used 500k samples. 
Figure <ref> clearly shows that saVANt with a single trained sampler accurately predicts the energies for different temperatures and reduces the runtime considerably.  We also observe that the estimator by original VAN is not fully compatible with the reference value provided by the Wolff algorithm <cit.>, whereas saVANt is.  




For analyzing complex physical systems,

the unbiasedness of the estimators is of particular importance 
since any GNS can hardly be expected to capture the underlying physics perfectly at reasonable training costs.  Practically, the saVANt framework could allow to predict observables close to the critical point based on samplers trained further away from it. More importantly, we consider saVANt a viable strategy to explore expensive parts in parameter space based on models trained in cheaper regions. Furthermore, our saVANt-NMCMC can be easily combined with well-established Monte-Carlo methods as a generalized overrelaxation step. We believe that this may alleviate the local minima problem and large autocorrelation times—key issues in
studying critical phenomena and the continuum limit of lattice field theories. 






Acknowledgment–
We thank Wu et al. for providing their implementation. This work was supported by the German Ministry for Education and Research as Berlin Big Data Center (01IS18025A) and Berlin Center for Machine Learning (01IS18037I). This work was also supported by the Information & Communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (No. 2017-0-00451).  

















§ ADDITIONAL DETAILS ON ALGORITHM



 §.§ Lightning Review of VAN

Wu et al approximate Boltzmann distributions p(s)=1/Z exp(-β H(s) ) with an auto-regressive generative model q by minimizing the KL divergence

    KL( q | p)    = ∑_s q(s) lnq(s)/p(s) = ∑_s q(s) (  lnq(s) + β H(s) + β Z )  .

A PixelCNN is used which allows exact evaluation of the probability q(s) and relatively efficient sampling. Also note that the partition function Z is a constant and  therefore the last summand leads to no contribution to the gradient. As a result, VAN can be trained by sampling from the model q, evaluating the probability q(s) and the Hamiltonian H(s) for these samples and using gradient descent.  



 §.§ Bounding the output probabilities of VAN


We can simply interpret the original network output q' ∈ [0, 1] as the probability q ∈ [ϵ, 1-ϵ] by the following mapping:

    q = (q' - 1/2) (1 - 2 ϵ) + 1/2.





 §.§ Proof: Estimators are Asymptotically Unbiased

Assume that the support of the sampling distribution q contains the support of the target distribution p. This property is ensured by ensuring that the probability takes values in q ∈ [ϵ, 1-ϵ].


  §.§.§ Neural Importance Sampling

Then, importance sampling with respect to q, i.e. 

    ⟨𝒪(s) ⟩≈∑_i=1^N w_i  𝒪(s_i)  ,         s_i ∼ q(s)  ,        w_i = ŵ_i/∑_i ŵ_i ,        ŵ_i = e^-β H(s_i)/q(s_i) .

is an asymptotically unbiased estimator of the expectation value ⟨𝒪(s) ⟩ because

    ⟨𝒪(s) ⟩_p = ∑_s p(s) 𝒪(s) = ∑_s q(s)  p(s)/q(s) 𝒪(s) = 1/Z∑_s q(s) exp(-β H(s))/q(s)_=ŵ(s) 𝒪(s) ≈1/Z N∑_i=1^N ŵ(s_i) 𝒪(s_i)  ,

where s_i ∼ q. The partition function Z can be similarly determined

    Z = ∑_s exp(-β H(s)) = ∑_s q(s) exp(-β H(s))/q(s)≈1/N∑_i=1^N  ŵ(s_i)  .

Combining the previous equations, we obtain 

    ⟨𝒪(s) ⟩_p≈∑_i w_i  𝒪(s_i)       with       w_i = ŵ_i/∑_i ŵ_i .





  §.§.§ Neural MCMC Sampling

The sampler q can be used as a trial distribution p_0(s' | s) = q(s') for a Markov-Chain which uses the following acceptance probability in its Metropolis step

    p_a(s'| s) = min(1, p_0(s| s') p(s')/p_0(s'| s) p(s))  = min(1, q(s)  exp(-β H(s'))/q(s')  exp(-β H(s)))   .

This fulfills the detailed balance condition

    p_t(s'|s) exp(-β H(s)) = p_t(s|s') exp(-β H(s))

because the total transition probability is given by p_t(s'|s)=q(s')   p_a(s'|s) and therefore

    p_t(s'|s) exp(-β H(s))    = q(s')  min(1, q(s)  exp(-β H(s'))/q(s')  exp(-β H(s)))  exp(-β H(s)) 
       = min{q(s') exp(-β H(s)),   q(s)exp(-β H(s'))}
       = p_t(s|s') exp(-β H(s'))  ,

where we have used the fact that the min operator is symmetric and that all factors are strictly positive. The latter property is ensured by the fact that q(s) ∈ [ϵ, 1-ϵ]. 





§ ADDITIONAL DETAILS ON EXPERIMENTS



 §.§ Setup

We use a Tesla P100 GPU with 16GB of memory both for training and sampling. A 16×16 lattice is considered. The reference values are generated using the Wolff algorithm using 2M steps with 100k warm-up steps. We use the ResNet version of VAN with the following hyperparameter choices (chosen to match with the ones used by Wu et al.):


The reference implementation of Wu et al. is used to train the VANs. For estimating the results of VAN and saVANt-NIS, we use 1000 iterations sampling 500 configurations each. For saVANt-NMCMC, we use 100k steps sampling 500 candidate configurations in a batch. No warm-up steps are required because candidates are sampled from a pre-trained VAN. As is demonstrated in Table <ref>, all algorithms have roughly the same runtime for sampling but saVANt leads to significant reduction in training time as explained in the main text.




 §.§ Additional Results



  §.§.§ Properties of saVANt-NMCMC

saVANt-NMCMC is less likely to get stuck in local minima and has reduced autocorrelation times as compared to Metropolis. This is demonstrated in Table <ref> and <ref>. 





  §.§.§ Using Models Trained at Different Temperatures


As already shown in the main text, Fig. <ref> and Fig. <ref> demonstrate that the bias for VANs holds for different betas. In particular, we looked at β=0.48 and β=0.57. As for the main study shown in the manuscript for β=0.36, here we trained a VAN with the aforementioned setup at the reference β, and we subsequently used this model to predict energies at different β. We note that VAN does not reproduce the reference values for almost all values of β. This discrepancy is particularly pronounced as one approaches the critical temperature. Fig. <ref> shows the aforementioned transfer property from a different perspective.



 

 



