Multilevel Generative Samplers for Investigating Critical Phenomena

W4: The appendix introduces unrelated terminologies (e.g., Swendsen-Wang algorithm, Wolff algorithm) without meaningful connections to the core work or further experimental analysis. Moreover, the lack of theoretical support or detailed experimental setups in the appendix limits its usefulness.

• Appendix on Cluster Algorithms: Appendix B.2 is meant to be a self-contained introduction to two popular Cluster algorithms, one of which is used in our experiment as the “cluster algorithm”. As the referee pointed out, the connection is unclear because we did not mention that the Wolff algorithm, effectively a cluster algorithm, is chosen as the cluster method in the main text. We will clarify these connections in the revised manuscript.

• Lack of theoretical support: We do not see what kind of theoretical support is necessary to support our method. As mentioned above, the training of RiGCS is guaranteed to converge because the training of ARNN is known to converge. The generative MC simulation with importance re-weighting is guaranteed to be unbiased [1] unless mode-dropping occurs. Our experiments empirically verify that RiGCS does not show mode-dropping. We request the reviewer to please concretely specify what kind of “lack of theoretical support” limits the usefulness of our method.

• Detailed experimental setups: We detailed the experimental setup for our numerical analysis in Appendix F. In the revised version of the manuscript we will also add pseudocodes for RiGCS. We would appreciate it if the reviewer could specify what relevant details are missing and we will be happy to add them in the revised manuscript.

W5: The manuscript cites more than 100 references, many of which seem only loosely related to the topic.

We agree that we cited many papers, including weakly related works. However, this is because we decided to have an “Extended related work” section in the Appendix, where the general influence of renormalization group theory (RGT) on machine learning is discussed. We argue that this part is useful for readers who are not necessarily familiar with the topic but may be interested in RGT and would consider incorporating it for their future machine learning research. By quoting the words of reviewer cr5x, we stress that our goal was to provide readers with “​​[..] an impressive bibliography on the subject which is indeed quite vast and old. A synthetic and pedagogical section in the supplementary material, which should be useful to non-physicists, is devoted to the RG theory and should make the whole paper readable for non-physicists”. However, we found that we did not clearly describe some of the references, and therefore, we revised the appendix so that the (either strong or weak) relations between RiGCS and other works become clear.

We appreciate the reviewer’s comments, which highlighted readability issues in the submitted version. We are committed to addressing these and will make every effort to improve the manuscript. We first refer the reviewer to the general reply, which includes an important correction to the experimental results in Fig.2 (left). Below we address the specific questions and concerns. Please note that we have not uploaded a revision yet, but we plan on doing it before the discussion phase ends.

Weakness 1: The core contributions (lines 90–99) are not sufficiently clear. See Questions 1-3.

Q 1: Why is addressing SIC important in this context? While the authors claim that SIC reduces sampling efficiency, they do not clearly explain how RiGCS specifically addresses this issue. A more detailed justification is required.

• Why is addressing SIC important: The reason is nicely described by Reviewer cr5x in the summary part. Taking their words, we will reinforce our explanation in Section 1 by adding the following: “At the criticality, physical systems typically exhibit self-similarity with respect to the change of scale, i.e., the physics of the coarse-grained system is similar to that of the fine-grained one. This requires us to deal with arbitrarily long-range correlations for which standard MCMC samplers based on local moves undergo critical slowing down, meaning that the auto-correlation time becomes arbitrarily large with system size.”

• How RiGCS specifically addresses this issue: In general, any distribution can be decomposed as Eq.(7). Although the left-hand side $p(s)$ exhibits long-range correlations, each of the factors \{p(s^l|s^{\leq l-1}\)} and $p(s^0)$ in the right-hand side exhibit only short-range correlations (with respect to the corresponding resolution scale), according to the renormalization group theory. RiGCS approximates not directly $p(s)$ with long-range correlations but each factor in the right-hand side with short-range correlations, using generative models with small receptive fields. We can intuitively say that RiGCS addresses the issue, i.e., long-range correlations induced by SIC, by capturing the long-range correlations as short-range correlations in the coarse level marginal distributions.

Q2: Although the warm start is listed as a key contribution, the implementation details are unclear. How is the warm start applied across different resolution levels, and why is it considered a significant innovation?

• Clarity: We agree with the reviewer that the implementation details of the training procedure were unclear. We will improve the explanation and add an illustration of the workflow as well as a detailed pseudocode.
• Innovation: We argue that our warm start (sequential) training with model transfer is a significant innovation as a novel physics-informed training procedure. First of all, model initialization is crucial when training a model by minimizing the reverse KL divergence because samples from badly initialized models do not provide informative gradient signals for training. The situation gets exponentially worse for larger lattice systems. Our sequential training with model transfer leverages SIC, i.e., the self-similarity across scales, and starts from training models for small lattice sizes, e.g., $2\times 2$ and $4\times 4$. Then, to train a larger model, e.g., $8 \times 8$, each conditional model in Eq.(7) has a similar conditional model already trained for smaller lattice sizes, which can be used for initialization. This idea is novel to the best of our knowledge and effective according to our experiment shown in Fig. 2 right. We will improve the description of this training procedure (see below).

Q3: If the authors only employ standard metrics (e.g., effective sample size, free energy bias), it is unclear why performance evaluation is listed as a major contribution. Were any new metrics introduced? If not, the framing of performance evaluation as a contribution should be reconsidered.

We moved the corresponding sentence outside the list of contributions in the revised manuscript.

I appreciate the authors' detailed response, which clarified some aspects of the submission. I have a few follow-up questions:

• The authors mentioned a novel "physics-informed training procedure" in their response. Could the authors further clarify what "physics-informed" entails? Does it refer to incorporating physical constraints, governing equations, or other principles? I could not find related details in the manuscript or the rebuttal.

• Could the authors explain how a trained $2 \times 2$ model is used to initialize a $4 \times 4$ model?

• I have not yet seen the updated revision of the manuscript. I will reconsider my score once the updated version is available.

We thank the referee for their comments, which informed us that we did not explicitly outline the crucial requirements for generative models in the context of our target problem—Monte Carlo simulations for estimating physical observables. First of all, we have an important correction to the experimental results in Fig.2 (left). Please see our general reply and also note that we have not uploaded a revision yet, but we plan on doing it before the discussion phase ends.

My first criticism is that the authors do not mention these works or position their contribution in the backdrop of this existing work. Can the authors comment on how their problem setting is substantively different from what is tackled in these existing works, and/or why the approaches from this line of work cannot be applied "zero-shot" to this application domain?

As the referee suspected, there is a fundamental difference in the problem setting from the typical applications in computer vision and natural language processing, which prevents us from using the state-of-the-art architectures (e.g., transformers, GANs, diffusion models) for general machine learning applications. The reasons are twofold.

• 1. The exact sampling probability needs to be efficiently accessible: When estimating physical observables from Monte Carlo (MC) simulations, one estimates both the mean value and the uncertainty, which is crucial to provide confidence intervals, e.g., for validating theory (e.g., the standard model) with experiments (e.g, high energy physics experiments in a collider). The confidence interval is estimated by making the estimators unbiased, and evaluating their variances (the variance is easily evaluated, while the bias is not unless the true value is known). Unlike MCMC methods, which asymptotically provide samples exactly from the target distribution, generative neural samplers (GNSs) can not directly sample from the exact target distribution, which makes the MC estimator of GNSs biased. As mentioned in the first paragraph of Section 2.2, one can make the estimator unbiased by re-weighting (Neural Importance sampling or NeuralMCMC [1]) if and only if the sampling probability is accessible. For this reason, modern research in the field focuses on using generative models with accessible sampling probabilities, such as autoregressive models, normalizing flows, and variants thereof, which are standard choices for sampling discrete and continuous variables. Diffusion models are actually potential candidates because the stochastic differential equation (SDE) for the reverse diffusion has an equivalent ordinary differential equation (ODE) (i.e., probability flow, continuous flow, or neural ODE). However, the SDE is solved approximately, and an efficient method to compute the exact sampling probability is yet to be established, to the best of our knowledge.
• 2. Training needs to be done without a large number of samples from the target distribution: Another crucial difference from the more typical ML problems is that we cannot assume the existence of the training data drawn from the target distribution. If a large amount of data from the target distribution, on which the state-of-the-art architectures like transformers with large receptive fields can be trained, are available, one can simply compute the MC estimator for the observables without training a generative model. Although one could assume some situations where sufficiently large training data for a few physical systems, e.g., Ising models for a few coupling parameter values, are available, and a conditional generative model is trained only on those data to generalize it for many other parameter settings, this is not a common scenario and indeed is not our setting. We refer to a recent review paper, and references therein, where the problem of using generative models for Lattice QCD is tackled [1].

Because of these two requirements, state-of-the-art ML models like transformers with large receptive fields cannot be used for our target application. Following previous state-of-the-art generative neural samplers for the Ising model, VAN [2] and HAN [3], we adopted the PixelCNN architecture as the basic component of our generative model. Consequently, we need to limit the receptive field size of the PixelCNN for tractability. As shown in Fig. 2 (left), the plain VAN with receptive fields containing all lattice sites can be trained only up to systems of size $32 \times 32$ (training time explodes to several weeks/months for larger systems).

We will make these requirements explicit in the revision.

We have now uploaded a revision. We appreciate the reviewer’s follow-up questions, giving us a further chance to convince them. We agree that the model and the architecture can be distinguished, and we understand the reviewer’s request to evaluate the VAN model with a transformer architecture.

We explored three possibilities to conduct the requested experiment, which could be finished by the deadline.

• 1. Looking for existing methods with their code available and testing one of them.
• 2. Implement a VAN with a transformer by ourselves.
• 3. Test our RiGCS with a transformer as the basic building block.

For Option 1, we found Transformer Neural Autoregressive Flows [1], which we could try to train with self-generated samples (via the reverse KL divergence) or with limited numbers of training data drawn by MCMC methods from the target Boltzmann distribution (via the forward KL divergence or the maximum likelihood training). However, the paper has not been published in conference nor journal (although presented in a workshop), and therefore they have not pushed their code into their GitHub repository. We further explored the literature, but so far we could not find a suitable existing method. We are happy to try this option by the (extended!) deadline if the reviewer would suggest an existing method with the code publicly available, such that we can easily test its performance.

For Option 2, we considered if there would be an easy way to implement VAN with a transformer, e.g., by replacing the PixelCNN with a transformer in the VAN implementation [2]. However, we found that it is not straightforward as we would need to carefully design a transformer-based layer that can deal with masked weights. A priori, it is hard to determine whether it is more efficient to use the standard encoder-decoder transformer architecture or to use transformer layers with masked weights like in the MADE model. Implementing this approach would require substantial time and rigorous testing, making it unrealistic to complete within the limited timeframe of the rebuttal period.

Option 3 is viable, as in the conditional models for RiGCS one masks the lattices rather than the weights, as usually done with standard autoregressive models like MADE or PixelCNN. We thus replaced CNNs with transformer layers and tested them by keeping the PixelCNN at the coarsest level $l=0$ (because it is the plain VAN). The result for $16 \times 16$ lattice can be found at the ANONYMIZED LINK, where RiGCS with Transformer (blue) and RiGCS with CNN (orange) are compared in terms of ESS during training. Note that the plot also shows the pre-training phase, i.e., the first phase (1-1001 epochs) corresponds to the pretraining of unconditional VAN at $l=0$ (where RiGCS with transformer also uses CNN and thus the performance is comparable). The second (1002-11002 epochs) and the third (11003-21003 epochs) phases involve training of the conditional VANs at $l=1, 2$, and the last phase (21004- epochs) corresponds to the training of the full RiGCS after pre-training for $N=16$. We clearly see that RiGCS with CNN outperforms RiGCS with Transformer. Note that the ESS shown in the figure is for $2 \times 2$, $4 \times 4$, $8 \times 8$, and $16 \times 16$, respectively, in the first, second, third pretraining phases, and the last full training phase.

We understand that Option 3 is not exactly what the reviewer requested. However, we found it to be the only feasible implementation given the time-constraints of the rebuttal period. More importantly, our result is consistent with the observed preference for CNNs over transformers in the field of generative neural samplers (GNSs). For example, in [3], where GNSs are used to simulate lattice QCD, the authors concluded in the last paragraph of Section V-A that “[...] While the attention mechanism is at the heart of some of the most powerful artificial intelligence applications to date, we find that in our experiments, convolution architectures generally train faster and result in higher performances. [...]”. Considering that negative results are less likely to be published, the current preference for CNNs over transformers in lattice sampling is also evident from the fact that public implementations of GNSs using transformers are rare (as seen in our challenges with Option 1), despite that CNN-based normalizing flows and ARNNs are widely available.

While we may not fully address the reviewer’s request, we kindly ask for consideration of the practical constraints (Options 1 and 2 are challenging within a short timeframe) and the prevailing consensus in the GNS field that CNNs are currently the more practical and effective choice of architecture.

We thank the referee for raising important points and for their thoughtful questions. We have an important correction to the experimental results in Fig.2 (left). Please see our general reply and also note that we have not uploaded a revision yet, but we plan on doing it before the discussion phase ends. Below we reply to the reviewer’s enquiries.

Weakness: I found some aspects of the experimental not completely satisfactory. The claim that RiGCS performs better than Wolff is a strong claim that should be better sustained in my opinion. Authors concentrate in the experiments on the computation of the free energy, but characterization of physical properties involve many other possible indicators, which the Wolff algorithm is known to directly compute very efficiently, like the magnetization, susceptibility [...].

In the original submission, we claimed that RiGCS outperforms the state-of-the-art generative model baselines (in the abstract), and that it outperforms the cluster method in free energy estimation (Line 402). Note that this claim is not as strong as claiming that RiGCS performs better than the Wolff algorithm in general. However, we need to modify our claims because of the error reported in the general reply. In the revision, we do not claim that RiGCS outperforms the cluster method in any of our experiments. As the referee mentioned, the cluster method is a very strong sampler for the Ising model, and no generative neural sampler has been able to outperform it. In the revised manuscript, we will make it clearer that our RiGCS outperforms only the state-of-the-art generative samplers.

The combination of Wolff with AIS even though reasonable seems a bit like of an ad-hoc procedure to yield a baseline comparison for the computation of the free energy, I would expect that integrating the mean energy with respect to beta or simply using directly the empirical distribution of energy delivered by the Wolff algorithm would yield more robust estimation but I maybe wrong?

We chose Cluster-AIS as the state-of-the-art method for free energy computation with MCMC samplers. Recent works [1-3] showed that AIS combined with Cluster/MCMC significantly outperforms the integral method, where fine discretization of the temperature parameter is required to ensure overlapping distribution supports between the subsequent steps. We do not know how to compute the free energy directly from the “empirical distribution of energy” using the Wolff algorithm. We would appreciate it if the referee could provide more explanations or references about it.

The size itself (128x128) seems rather limited compared to what can commonly be achieved with the Wolff algorithm. It seems to be a limitation in term of computation time but this is not discussed. Actually the scaling of the the computational time of RiGCS w.r.t to L is actually not given. Question: How does the algorithm scale system size (even empirically) ?

We agree that $L=128$ is limited compared to the capability of the Wolff algorithm. In the revision, we made explicit that we claim state-of-the-art performance only with respect to the generative neural sampler (GNS) baselines, i.e., VAN and HAN. In the context of GNSs, sizes of $128\times128$ already represent a large and challenging problem in comparison to other recent works [4-6].

Leveraging physical information, i.e., scale invariance at criticality (SIC), we designed RiGCS with multilevel conditional samplers, where the conditional samplers are convolutions with a fixed filter size (receptive field) for all levels. This allows linear model complexity (the number of network parameters) with respect to the number of levels, e.g., complexity scales linearly with the lattice sites $(V = L^2)$. We will add this argument and empirical computation time, which verifies the linear scaling, in the revised manuscript.

We thank the reviewer for their encouraging comments and their thoughtful questions. We have an important correction to the experimental results in Fig.2 (left), please see our general reply. Please also note that we have not uploaded a revision yet, but we plan on doing it before the discussion phase ends.

First, we thank the reviewer for the positive assessment of our work.

Strengths: The proposed method utilizes SIC effectively, enabling ancestral sampling from low-resolution to high-resolution lattice configurations with site-wise-independent conditional HB sampling. This innovation enhances both scalability and efficiency, making it a significant contribution to the field of generative models.

Below, we reply to the reviewer’s inquiries.

Weaknesses (1): It appears that the training of the algorithm at each resolution follows a simultaneous sampling and training paradigm. However, this is not clearly articulated in the manuscript. I recommend providing complete pseudocode for the algorithm, which would greatly aid readers in understanding the implementation and workflow.

Question (1): Provide a detailed algorithm or flowchart that outlines the key steps of RiGCS, including how training proceeds at each resolution level. This would help clarify the simultaneous sampling and training processes.

In the revised manuscript, which we will upload by the end of the discussion phase, we will include a figure explaining the training workflow in the main text, as well as a pseudocode in the appendix F.6.

Weaknesses (2): The loss function presented in equation (11) is known to have a mode-seeking tendency when training generative models. This could lead to the neglect of certain peaks in multimodal distributions. It would be beneficial for the authors to address whether RiGCS is susceptible to similar issues and, if so, how they might mitigate them. Question (2) Address whether RiGCS is susceptible to the mode-seeking, if so, how they might mitigate them.

We agree that the mode-seeking tendency is a concern. In the manuscript we evaluated the effective mode-dropping measure (EDMD), a metric with which the bias caused by mode-dropping can be bounded [1]. We observe in Fig. 3 left that RiGCS is hardly prone to mode-dropping, unlike the baseline GNSs (VAN, HAN). This result is consistent with the free energy estimation (Fig.2 left). When RiGCS would suffer from the mode-dropping, standard remedies from the literature, such as annealing or forward-KL training [1], should be applied to mitigate these effects.

Weaknesses (3): The paper does not mention an important class of methods for the Ising model that rely on relaxation to continuous variable representations. Approaches such as the Neural Network Renormalization Group by Li and Wang use this relaxation technique to enable generative sampling in a continuous variable space. Including a discussion of these relaxation-based methods would provide a more comprehensive comparison and allow readers to better appreciate where RiGCS stands in relation to existing methods.

Question (3): Compare the proposed approach to methods like the Neural Network Renormalization Group, highlighting potential advantages or trade-offs of RiGCS in relation to these continuous relaxation based techniques.

NeuralRG [2] aims to learn the renormalization operations with hierarchical bijective mapping, where the depth of the bijective layers depends on the level of the latent variables in the hierarchy. Although this approach reduces the degrees of freedom of the model (i.e., the number of the network parameters) for efficiency, the reduction is not as drastic as for our RiGCS, because long-range correlations are still needed to be captured by the high-resolution layers. In contrast, our RiGCS significantly limits the receptive field size to be constant for all layers, based on the function-forms of the renormalized Hamiltonians, leading to significant reduction of the model complexity. As a quick evaluation, we tested the code for Neural RG [2] on a local GPU (using device=mps on M1 MacBookPro with 32G RAM). We observed that training Neural RG takes 1h 30min for L=16, which is much slower than training RiGCS (which takes less than 15 min), supposedly due to its higher model complexity. We include NeuralRG [2] and discuss the comparison with our method in the related work section of the revision.

References

• [1] Nicoli, Kim A., et al. "Detecting and mitigating mode-collapse for flow-based sampling of lattice field theories." Physical Review D 108.11 (2023): 114501.
• [2] Li, Shuo-Hui, and Lei Wang. "Neural network renormalization group." Physical review letters 121.26 (2018): 260601.

• [1] Cranmer, Kyle, et al. "Advances in machine-learning-based sampling motivated by lattice quantum chromodynamics." Nature Reviews Physics 5.9 (2023): 526-535.
• [2] Nicoli, Kim A., et al. "Estimation of thermodynamic observables in lattice field theories with deep generative models." Physical review letters 126.3 (2021): 032001.
• [3] Kanwar, Gurtej, et al. "Equivariant flow-based sampling for lattice gauge theory." Physical Review Letters 125.12 (2020): 121601.
• [4] Wu, Dian, Lei Wang, and Pan Zhang. "Solving statistical mechanics using variational autoregressive networks." Physical review letters 122.8 (2019): 080602.
• [5] Li, Shuo-Hui, and Lei Wang. "Neural network renormalization group." Physical review letters 121.26 (2018): 260601.

We sincerely thank all reviewers for the active discussion in the rebuttal phase. We would like to inform the reviewers that we have uploaded a revision, where most of the reviewers’ suggestions have been incorporated, except some parts related to additional experiments. Revised parts except minor corrections, e.g., typos, are highlighted in blue. A major revision was done in Section 3.3, where the sequential training with model transfer is explained. We added Fig. 2 that intuitively explains the procedure, including how the trained model parameters are used to initialize the parameters in the subsequent steps (levels). Due to the space limit, we moved a large part of the mathematical descriptions to Appendix E.1, where the explanations on how the parameters are initialized and on how SIC justifies the sequential training process are further elaborated. Pseudocode for training and sampling is also added there.

Some experimental parts have not been fully updated yet as they will depend on the outcome of the suggested additional experiments, for which we could only conduct preliminary experiments so far. Specifically, we implemented a continuous version of RiGCS (as well as VAN) with continuous autoregressive neural network (ARNN) using a mixture of Gaussians [1], and applied it to the two-dimensional $\phi^4$ lattice field theory up to $N=32$, see the ANONYMIZED LINK. Similarly to the Ising model case, RiGCS significantly outperforms the plain VAN, which could be trained only up to $N=16$ within reasonable computation time. Note that the ESS of both models are worse in the continuous $\phi^4$ theory compared to the discrete Ising model (Fig. 4 right), not only because the theory is harder, but also because the continuous ARNN might not be the optimal choice as the basic building block. We conjecture that RiGCS based on other models, e.g., affine-coupling normalizing flows and continuous normalizing flows, could significantly improve the ESS. We would like to leave this exploration to our future work.

Refs.

• [1] Uria, Benigno, et al. "Neural autoregressive distribution estimation." Journal of Machine Learning Research 17.205 (2016): 1-37.