Gradient-based Discrete Sampling with Automatic Cyclical Scheduling

Authors: patrick  pynadath
riddhiman  bhattacharya
arun  hariharan
ruqi  zhang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-10 21:55:31

Abstract: Discrete distributions, particularly in high-dimensional deep models, are often highly multimodal due to inherent discontinuities. While gradient-based discrete sampling has proven effective, it is susceptible to becoming trapped in local modes due to the gradient information. To tackle this challenge, we propose an automatic cyclical scheduling, designed for efficient and accurate sampling in multimodal discrete distributions. Our method contains three key components: (1) a cyclical step size schedule where large steps discover new modes and small steps exploit each mode; (2) a cyclical balancing schedule, ensuring "balanced" proposals for given step sizes and high efficiency of the Markov chain; and (3) an automatic tuning scheme for adjusting the hyperparameters in the cyclical schedules, allowing adaptability across diverse datasets with minimal tuning.  We prove the non-asymptotic convergence and inference guarantee for our method in general discrete distributions. Extensive experiments demonstrate the superiority of our method in sampling complex multimodal discrete distributions.

Main paper: 
 introduction: 
 Discrete variables are common in many machine learning problems, highlighting the crucial need for efficient discrete samplers. Recent advances [Grathwohl et al., 2021 , Zhang et al., 2022b , Sun et al., 2021 , 2023b ,a, Xiang et al., 2023] have leveraged gradient information in discrete distributions to improve proposal distributions, significantly boosting their efficiency. These advancements have set new benchmarks in discrete sampling tasks across graphical models, energy-based models, and combinatorial optimization [Goshvadi et al., 2023 ].However, one major limitation of gradient-based methods is their susceptibility to becoming trapped in local modes [Ruder, 2016 , Ziyin et al., 2021] , which significantly reduces the accuracy and efficiency of sampling results. In continuous spaces, several strategies such as cyclical step sizes [Zhang et al., 2020] , parallel tempering [Swendsen and Wang, 1986, Deng et al., 2020a] , and flat histograms [Berg and Neuhaus, 1991, Deng et al., 2020b] , have been proposed to address this issue. When it comes to discrete distributions, which are inherently more multimodal due to their discontinuous nature, the problem becomes even more severe. Despite the pressing need, there is a lack of methodology for gradient-based discrete samplers to effectively explore multimodal distributions. Current methods 38th Conference on Neural Information Processing Systems (NeurIPS 2024). [Zhang et al., 2022b] with a manually tuned step size. (d): results from AB [Sun et al., 2023a] . (e): results from our method ACS. While the random walk sampler can find all modes, its characterization is noisy and lacks details for each mode.Gradient-based samplers (b) and (c) effectively characterize a specific mode but are easily trapped in some local modes. Our method (d) can find all modes efficiently and characterize each mode accurately.often fall far short in traversing the complex landscapes of multimodal distributions, as illustrated in Figure 1 .In this paper, we propose automatic cyclical scheduling for gradient-based discrete sampling to efficiently and accurately sample from multimodal distributions. To balance between uncovering new modes and characterizing the current mode, we parameterize a family of gradient-based proposals that span a spectrum from local to global proposals. The parameterized proposal dynamically adjusts according to cyclical schedules of both step size and the balancing parameter, smoothly transitioning from global exploratory moves to more localized moves within each cycle. These cyclical schedules are automatically tuned by a specially designed algorithm, which identifies optimal step sizes and balancing parameters for discrete distributions. Our contributions are summarized as follows:• We present the first gradient-based discrete sampling method that targets multimodal distributions. Our method incorporates cyclical schedules for both step size and balancing parameter to facilitate the exploration and exploitation in discrete distributions.• We propose an automatic tuning algorithm to configure the cyclical schedule, enabling effortless and customized adjustments across various datasets without much manual intervention.• We offer non-asymptotic convergence and inference guarantees for our method in general discrete distributions. To our knowledge, this is the first non-asymptotic convergence bound of gradient-based discrete sampling to the target distribution with inference guarantees, which could be of independent interest.• We demonstrate the superiority of our method for both sampling and learning tasks including restricted Boltzmann machines, deep energy-based models, and large language models. 
related work: 
 Gradient-based Discrete Sampling Zanella [2017] introduced a family of locally informed proposals, laying the foundation for recent developments in efficient discrete sampling. Building upon this, Grathwohl et al. [2021] further incorporates gradient approximation, significantly reducing computational costs. Following these pioneering efforts, numerous studies have proposed various gradient-based discrete sampling techniques [Rhodes and Gutmann, 2022 , Sun et al., 2021 , 2022 , 2023b , Xiang et al., 2023] . Zhang et al. [2022b] develops a discrete Langevin proposal, translating the powerful Langevin algorithm to discrete spaces. Sansone [2022] introduces a self-balancing method to optimize the balancing functions in locally balanced proposals. While our work also utilizes an adaptive phase, it differs in that our parameterization extends beyond the local regime, and our proposal parameterization is considerably simpler.Perhaps the most closely related study is the any-scale balanced sampler [Sun et al., 2023a] . This method uses a non-local balancing proposal and adaptively tunes it. Our work, however, differs in several key aspects: (1) We focus on combining both local and non-local proposals to effectively characterize multimodal discrete distributions, as opposed to focusing on a single optimal proposal.(2) Our automatic tuning algorithm adjusts the step size and balancing parameter by considering the special discrete structures and targets a specific Metropolis-Hastings acceptance rate, rather than maximizing the average coordinates changed per step.(3) Our method can be applied to learning energy-based models (EBM) and sampling large language models, whereas their approach cannot.Sampling on Multimodal Distributions There exist several sampling methods targeting discrete multimodal distributions, such as simulated tempering [Marinari and Parisi, 1992] , the Swendsen-Wang algorithm [Swendsen and Wang, 1987] , and the Wolff algorithm [Wolff, 1989] . However, these methods usually use random walk or Gibbs sampling as their proposals. It is unclear how these methods can be adapted for gradient-based discrete sampling.In continuous spaces, various gradient-based methods have been developed specifically for multimodal distributions [Zhang et al., 2020 , Deng et al., 2020a,b] . Our method distinguishes from the cyclical step size in Zhang et al. [2020] by incorporating an additional cyclical balancing parameter schedule and an automatic tuning scheme, which are crucial for efficient exploration in discrete distributions. Furthermore, our theoretical analysis of convergence is different from that in Zhang et al. [2020] which relies on continuous stochastic processes.3 Preliminaries 
problem definition: 
 We consider the task of sampling from some target distribution defined over a discrete spaceπ(θ) = 1 Z exp(U (θ)), θ ∈ Θ.Here, θ is a d dimensional discrete variable in domain Θ, U is the energy function, and Z is the normalizing constant. We make the following assumptions of the domain and the energy function, following the literature of gradient-based discrete sampling [Grathwohl et al., 2021 , Sun et al., 2021 , Zhang et al., 2022b] : (1) The domain is coordinatewisely factorized, Θ = Π d i=1 Θ i .(2) The energy function U can be extended to a differentiable function in R d .Zanella [2017] introduces a family of informed proposals, which is defined below:EQUATIONHere, K α is a kernel that determines the scale of the proposal where α plays a similar role as the step size. g(t) is a balancing function that determines how to incorporate the information about π. If g(t) = tg( 1 t ), the proposal becomes a locally balanced proposal, which is asymptotically optimal in the local regime, that is, when the step size α → 0. 
automatic cyclical sampler: 
 We aim to develop a sampler capable of escaping local modes in general multimodal discrete distributions, including those that appear in deep energy-based models and large language models. First, we motivate using the cyclical schedule by demonstrating the issue of gradient-based samplers getting stuck in local modes on a toy dataset. We then present our sampler's parameterization of the step size and balancing function. Next, we introduce a cyclical schedule for the proposal distribution that enables effective exploration and characterization of discrete multimodal distributions. Finally, we develop an automatic tuning method that simplifies the process of identifying hyperparameters in cyclical schedules.To demonstrate the crucial issue of local modes trapping gradient-based samplers, we construct a 2dimensional dataset consisting of integers. We define Θ = {0, 1, • • • N } 2 , where N is the maximum value for each coordinate. Given a set of modes {µ 1 , µ 2 , . . . µ l }, we define the energy as follows:EQUATIONThis distribution enables easy comparison between different methods in terms of their ability to both explore and exploit the target distribution. We demonstrate the results of various samplers in Figure 1 .More experimental details can be found in Appendix D.1.A visual comparison reveals that while gradient-based samplers (DMALA [Zhang et al., 2022b] and AB [Sun et al., 2023a] ) are very effective at characterizing a given mode, they tend to get trapped in some small neighborhood, preventing a proper characterization of the distribution as a whole.We can understand this behavior of gradient-based samplers by comparing them to a random walk sampler (RW), which is able to explore all the modes but unable to fully characterize the detail of each one. While the RW sampler proposes movements uniformly over the sample space, gradient-based samplers propose movement based on the geometry of the distribution as captured by the gradient.Because the proposed movements are in the direction of increasing density, these proposals are able to characterize a given mode in detail. At the same time, these proposals hinder escape to more distant modes as the gradient points away from their direction. For this reason, we observe that local modes are able to "trap" gradient-based samplers.To derive an automatic schedule for the proposal, we need to parameterize the proposal first. We define K α and g(t) in the informed proposal [Zanella, 2017] as follows:EQUATIONwhere β is called a balancing parameter. α → 0, β = 0.5 correspond to a locally-balanced proposal and α → ∞, β = 1 correspond to a globally-balanced proposal. Values in between result in interpolations between locally-balanced and globally-balanced proposals. Note that β ∈ (0, 1) in Sun et al. [2023a] while our range is narrower.We substitute these definitions into Equation (1) and apply the first-order Taylor expansion:EQUATIONAs in Zhang et al. [2022b] , we use the assumption of coordinate-wise factorization to obtain the following coordinate-wise proposal function:EQUATIONIn order to make the resulting Markov chain reversible, we apply the Metropolis-Hastings correction, where we accept the proposed step with the following probability:EQUATIONIn summary, we parameterize our proposal as in Equation ( 5) which includes a spectrum of local and global proposals. Our proposal is determined by two hyperparameters, the step size α and the balancing parameter β.Step Size Schedule In order to effectively explore the whole target distribution while retaining the ability to exploit local modes, we adopt the cyclical step size schedule from Zhang et al. [2020] . The definition of step size α for iteration k is as follows: EQUATIONθ ← samples[-1] 5: for coordinate i in range d do 6: construct Q i α k ,β k (•|θ) as in (5) 7: sample θ ′ i ∼ Q i α k ,β k (•|θ) 8:end for 9:samples ← θ ′ with probability (6) 10:end for 11: end for 12: return samples where α max is the initial step size, α min is the minimum step size, and s is the number of sampling steps per cycle. Differing from the cyclical schedule in Zhang et al. [2020] , we additionally add α min to make sure that even the smallest step size remains effective in discrete spaces.Cyclical Balancing Schedule Using large step sizes in (7) can easily result in very low acceptance rates, removing any benefit of exploration. To address this issue, we introduce a balancing parameter schedule, which enables reasonable acceptance rates for large step sizes. As discussed in Zanella [2017] , Sun et al. [2023a] , the balancing parameter should vary with different step sizes to achieve a "balanced" proposal. A balanced proposal ensures that the Markov chain is reversible with respect to a certain distribution, which will converge weakly to the target distribution. For example, when the step size α → 0, the optimal balancing parameter is β = 0.5, whereas for α → ∞, the ideal balancing parameter becomes β = 1.Thus for a schedule of step sizes, each α i requires a different β i ∈ [.5, 1), with larger step sizes having β i closer to 1 and smaller step sizes having β i closer to 0.5. Using the Metropolis-Hastings acceptance rate to characterize the quality of a given α, β pair, we define the value of β i as follows:EQUATIONIntuitively, this definition means that the best β i for a given step size α i maximizes the average acceptance rate for the proposal function Q α,β . It also conveys that larger step sizes will have larger balancing parameters.We include a visualization of the resulting schedules in Figure 2a and outline our algorithm using the α, β schedules in Algorithm 1. Note that it incurs no extra overhead compared to previous gradient-based discrete sampling methods as it only adjusts hyperparameters α and β. By using a combination of large and small α and β, we enable the sampler to explore the distribution fully without sacrificing the ability to characterize each mode. This is demonstrated in Figure 1e .Require: β min = .5, β max , target acceptance rate ρ * , initial state θ init , steps per cycle s, initial largest step size α ceil = 60, initial smallest step size α floor = .05 1: θ ← InitBurnin(α ceil , β max , θ init ) 2: α min ← EstAlpha(α floor , β min , θ, ρ * , MAX=False) 3: α max ← EstAlpha(α ceil , β max , θ, ρ * , MAX=True) 4: Construct α-sched of length s using (7) 5: β-sched ← EstBalSched(α-sched, β max , β min , θ) 6: return α-sched, β-schedFor schedules in Equations ( 7) and ( 8), we have parameters α max , α max , and {β 1 , β 2 . . . β s } to be decided. In this section, we will introduce an automatic tuning algorithm to easily find suitable values.Main Idea Our automatic tuning algorithm depends on the initial balancing parameter β max , the final balancing parameter β min , a target acceptance rate ρ * , and the number of steps per cycle s. These values are relatively easy to select, as detailed in Appendix A. Below, we assume they are already determined. The tuning algorithm first estimates the optimal choices for α max and α min based on ρ * , which can then be used to construct the full step-size schedule using (7). We then construct the balancing parameter schedule using (8). The method is summarized in Algorithm 2 with details regarding subroutines in Appendix A. Our automatic tuning introduces minimal overhead relative to the more expensive sampling process. For example, in Section 6, we use 500 steps as the budget for Algorithm 2 where the total number of sampling steps is at least 5000. We further demonstrate that our algorithm is relatively robust to hyperparameters in Appendix A.1.In short, our tuning algorithm adopts an alternative optimization strategy, leveraging existing knowledge about hyperparameter values (e.g. β min and β max should be around 0.5 and 1 respectively). While estimating the best pair α, β is challenging due to their interdependence, it is much easier to fix one and optimize the other [Sun et al., 2023a] .Estimating α max , α min For a given β max , β min , our goal is to find step sizes α max , α min that enable an acceptance rate close to ρ * . We can formally state this goal as follows:EQUATIONGiven β max , β min , we construct the following objectives to pick the corresponding α max , α min :EQUATION)By defining the initial and final step sizes in this manner, we ensure that our cyclical schedule includes a wide range of hyperparameter pairs with different trade-offs in exploration and exploitation.To solve (10), we estimate α max by starting with a large step size and gradually decreasing it to find the step size that yields ρ * . Unlike existing works that start with small step sizes, we observed that multiple α values can yield the same acceptance rate for a given β, as shown in Figure 2b . Therefore, we start with an upper limit α ceil and reduce the step size to avoid missing any larger α values that meet our criteria. Detailed implementation is provided in Algorithm 4 in the Appendix. α min can be obtained similarly.Estimating Balancing Schedule After setting the start and end pairs for the α and β schedules, we now define intermediate β values. As the entire step size schedule is fixed by ( 7), the problem is to determine the best balancing parameter for each step size. A simple strategy is to test different β spaced out evenly throughout the interval [.5, β i-1 ] and select the best value in terms of acceptance rate. This approach leverages the observation that smaller step sizes tend to have smaller optimal balancing parameters. Detailed implementation is given in Algorithm 5 in Appendix. 
theoretical analysis: 
 In this section, we present a convergence rate analysis for Algorithm 1. For general step size and balancing parameter schedules, i.e., at each cycle, the algorithm will go through s steps in which it will use step sizes α 1 , α 2 , • • • , α s and balancing parameters β 1 , β 2 , • • • , β s . Note that for each pair (α i , β i ), we have a Markov transition operator which we label P i for i = 1, 2, • • • , s. The Markov operator for a single cycle is given by P = P 1 P 2 • • • P s . We have the following two assumptions: Assumption 5.1. The functionU (•) ∈ C 2 (R d ) has M -Lipschitz gradient. That is ∥∇U (θ) -∇U (θ ′ )∥ ≤ M ∥θ -θ ′ ∥ .Note that it implicitly assumes that the set in domain Θ is finite. We define conv(Θ) as the convex hull of the set Θ. Assumption 5.2. For each θ ∈ R d , there exists an open ball containing θ of some radius r θ , denoted by B(θ, r θ ), such that the function U (•) is m θ -strongly concave in B(θ, r θ ) for some m θ > 0.Assumptions 5.1 and 5.2 are standard in optimization and sampling literature [Bottou et al., 2018 , Dalalyan, 2017] . Under Assumption 5.2, U (•) is m-strongly concave on conv(Θ), following Lemma C.3 in Appendix.We define diam(Θ) = sup θ,θ ′ ∈Θ ∥θ -θ ′ ∥ and ϵ αi,βi to be exp -1 2α i + β i M - β i m 2 diam(Θ) 2 -∥∇U (a)∥ diam(Θ) .The Markov kernel corresponding to each P i in each step of the cycle in Algorithm 1 isEQUATIONwhereL(θ) = θ ′ ∈Θ π(θ ′ )Q αi,βi (θ|θ ′ ) π(θ)Q αi,βi (θ ′ |θ) ∧ 1 Q αi,βi (θ ′ |θ)is the total rejection probability from θ. Finally, recall that the total variation distance between two probability measures µ and ν, defined on some spaceΘ ⊂ R d is ∥µ -ν∥ T V = sup A∈B(Θ) |µ(A) -ν(A)|where B(Θ) is the set of all measurable sets in Θ.To analyze Algorithm 1 with step size and balancing parameter schedules, we first solve a simpler problem where the step size and balancing parameter are fixed and then extend the analysis to the setting of Algorithm 1.Our main method of proof is to establish uniform ergodicity of the Markov chain P , for a single α, β, by establishing a uniform minorization for P . We denote the transition kernel for this Markov chain P as p(• | •), which is given in (11) with α i , β i replaced by a fixed α, β. Lemma 5.3. Let Assumptions 5.1-5.2 with α < 1 βM hold. Then for the Markov chain P we have, for any θ, θ ′ ∈ Θ,p(θ | θ ′ ) ≥ ϵ β,α exp {βU (θ ′ )} θ ′ ∈Θ exp {βU (θ ′ )}, whereϵ β,α = exp - 1 2α + βM - β m 2 diam(Θ) 2 -∥∇U (a)∥diam(Θ)}with a ∈ arg min θ∈Θ ∥∇U (θ)∥ .Proof. The proof is provided in Appendix C.1.Theorem 5.4. Let Assumptions 5.1-5.2 hold with α < 1/βM . Then for the Markov chain P , the following hold: i. P is uniformly ergodic with∥P n -π∥ T V ≤ (1 -ϵ β,α ) n .ii. For any real-valued function f and samples X 1 , X 2 , X 3 , • • • , X n from P , one has√ n 1 n n i=1 f (X i ) - θ∈Θ f (θ)π(θ) d → N (0, σ2 * )for some σ * > 0 as n → ∞.Proof. The proof directly follows from our Lemma 5.3 and Jones [2004][Corollary 5 ].Note that as α → 0, we have ϵ β,α → 1 which implies that small step sizes result in low convergence rates. This is intuitive as the algorithm could not explore much in this case. Furthermore, our results suggest that large β restricts α to small values. Given that large β generally requires large α, our findings imply an upper bound for the step size.Step Size and Balancing Parameter Now we tackle the case of varying step sizes and balancing parameters. Each cycle has s steps with step sizes α 1 , α 2 , • • • , α s and balancing parametersβ 1 , β 2 , • • • , β s .Note that this case is more challenging as at each step the transition operator changes and the Markov chain is no longer homogeneous. However, the marginal chain for each cycle is indeed homogeneous and can be analyzed. We present our results in this setting as follows:Theorem 5.5. Let Assumptions 5.1 and 5.2 hold withα i < 1/β i M , i = 1, 2, • • • s.Then for the Markov chain P , the following hold i. P is uniformly ergodic withP n -π T V ≤ (1 -ϵ βs,αs ) n .ii. For any real-valued function f and samplesX 1 , X 2 , X 3 , • • • , X n from P , one has √ n 1 n n i=1 f (X i ) - θ∈Θ f (θ)π(θ) d → N (0, σ2 * )for some σ * > 0 as n → ∞, where,ϵ βs,αs = exp - 1 2α s + β s M - β s m 2 diam(Θ) 2 • exp {-∥∇U (a)∥diam(Θ)} with a ∈ arg min θ∈Θ ∥∇U (θ)∥.Proof. The proof follows from our Lemma 5. Both Theorems 5.4 and 5.5 hold uniformly over all functions in the class of functions with at least a local minima in Θ. The Central Limit Theorem results in Theorems 5.4 and 5.5 imply that we may perform inference on the target distribution π(•) even though the asymptotic variances are unknown, as we may perform batch-means to estimate these variances Vats et al. [2019] .In summary, we have established a geometric convergence rate to the target distribution for our sampler. Previous research has only established asymptotic convergence [Zhang et al., 2022b] or relative convergence rate bounds [Grathwohl et al., 2021] for gradient-based discrete samplers. To the best of our knowledge, our results present the first non-asymptotic convergence bounds that explicitly quantify the distance between the estimated and target distributions. Further, our convergence bound also shows that discrete spaces play a fundamental part in the ergodic nature of these algorithms.  
experiments: 
 We call our method that combines Algorithm 1 and 2 Automatic Cyclical Sampler (ACS). For RBM and EBM sampling tasks, we compare our method to Gibbs-with-Gradient (GWG) [Grathwohl et al., 2021] , Any-scale sampler (AB) [Sun et al., 2023a] , and Discrete Metropolis Adjusted Langevin Algorithm (DMALA) [Zhang et al., 2022b] , which are popular and recent gradient-based discrete samplers. For learning tasks, we omit AB sampler as it is not originally applied to the model learning tasks. More experimental details are in Appendix D. We released our code at the following link: https://github.com/patrickpynadath1/automatic_cyclical_sampling.We evaluate our sampling method on both Restricted Boltzmann Machines (RBMs) and deep convolutional Energy-Based Models (EBMs). For RBMs, we measure accuracy by comparing the Maximum Mean Divergence (MMD) between samples generated by our method and Block Gibbs, which can be considered the ground truth. We sample on EBMs to demonstrate our method's scalability to more complex distributions. Experimental details are provided in Appendices D.2 and D.3 for RBM and EBM sampling, respectively.In Figure 3 , our proposed ACS method performs competitively for both RBMs and EBMs across all datasets. For RBM sampling, ACS is able to converge to the ground truth quicker than other methods due to the ability to capture the multi-modal nature of the target distribution. We see that this performance generalizes to more complex distributions as represented by deep EBMs.Table 1 : Deep Convolution EBM Log likelihood scores on test data as estimated by AIS. GWG results are taken from [Grathwohl et al., 2021] . ACS is able to achieve better results than the baselines. Results Table 1 demonstrates that ACS is capable of learning better quality EBMs given the same computational budget as DMALA. Furthermore, ACS learns better quality models with less computational budget than GWG.One challenging application of discrete MCMC methods is text-infilling, where the goal is to complete a sentence with some missing words. Given a dataset of sentences, we randomly mask our 50% of the words and fill them in using the distribution given by a pretrained RoBERTa model. We include experiment details in Appendix D.6.Results Table 2 demonstrates that ACS is capable of generating more diverse sentences, as ACS has a lower self-BLEU and higher percentage of unique n-grams. While the perplexity results seem to imply that ACS generates lower quality than DMALA, we note that the ACS generations are more likely to be predicted as linguistically acceptable as shown by the CoLA scores. We discuss the results more extensively in Appendix D.6.  
conclusion and limitations: 
 In this work, we propose Automatic Cyclical Sampler (ACS) to more effectively characterize multimodal distributions in discrete spaces. First, we demonstrate that gradient-based samplers are prone to getting trapped in local modes, preventing a full characterization of target distributions. To address this issue, we combine a cyclical step size schedule with a cyclical balancing parameter schedule along with an automatic tuning algorithm to configure these schedules. We also theoretically establish the non-asymptotic convergence bound of our method to the target distribution in addition to providing extensive experimental results.While our proposed ACS method generates impressive results on a wide range of experiments, there are some limitations to our work that should be mentioned. Specifically, though we have proven a geometric convergence rate and the relationship between α and β in our theoretical analysis, we require U (•) to be twice differentiable as well as locally strongly concave and the proof is not based on the specific tuning algorithm implemented. This is why we provide extensive experimentation to demonstrate that our algorithm is capable of picking good α, β schedules. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We have read the Ethics Guidelines, and our submission aligns with all the points listed.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: This is a MCMC sampling technique which does not have a direct societal impact. The research presented is either theoretical, or based on common benchmarks such as MNIST. The most relevant impact would be due to text infilling, which can be used for LLMs text generation.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: This foundational research that does not directly have a societal impact, as it is primarily an MCMC algorithm for discrete spaces. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [No] Justification: We were unable to find the license for the datasets we used. However, they are fairly popular and well known datasets, and we cite the relevant paper where necessary. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: We do not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This research does not involve crowdsourcing or human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? 
 
 Appendix: Here we include more details regarding the Automatic Cyclical Sampler algorithm. We discuss all the individual sub-routines that compose the algorithm shown in Algorithm 2. We also include an ablation study to demonstrate the robustness of our algorithm to various hyper-parameter configurations.InitialBurnin We find that in order to produce meaningful estimates for the objective in (9), it is necessary to burn in the MCMC sampling chain. This is due to the dependence of the acceptance rate on current sample θ. If we use θ very low in density with respect to the target distribution, the acceptance rates estimated by the tuning algorithm will lose accuracy as the sampler converges to the target distribution. In order to avoid this issue, we run a quick burn-in stage with two distinct stages.The first stage uses the gradient information to move the sampler away from the initialized point as quickly as possible. We use the parameterized proposal from Equation ( 4) with stepsize α ceil , β max without any Metropolis-Hastings correction as this enables very large movements from the initial sample.For some datasets, this enables a very quick burn-in. This can be noticed in Figure 3 for Static/Dynamic MNIST and Omniglot. We hypothesize that this is due to the distribution having a relatively simple structure that enables the gradient to provide meaningful information for very large sampling steps. It is impossible to determine a priori whether a given distribution will have this property, so we include a following stage that uses a Metropolis-Hastings correction to increase the chance of arriving at a reasonable sample θ.For this stage, we construct a naive step size schedule and balancing constant schedule using the values of α ceil , α floor , β max , β min . We then run the parameterized sampler from Equation ( 4) with the Metropolis-Hastings correction. Our goal is to move the sampler to samples θ that are more likely in the target distribution. This will enable the acceptance rates computed during the tuning algorithm to be closer to the acceptance rates for the steady-state chain.For all the sampling experiments, these two stages combined use 100 sampling steps.EstimateAlpha Here we discuss the algorithm used to calculate both α max , α min as defined in Equation ( 10). When calculating α max , the goal is to pick the largest stepsize α max that acheives the acceptance rate ρ * for a given β max . When calculating α min , the goal is to determine the smallest step-size capable of acheiving the target acceptance rate. We put the full pseudo-code in Algorithm 4.For calculating α max and α min , the algorithm follows the general pattern of automatically shifting the range of potential α based on the best values calculated from the previous iteration. When calculating α max , the algorithm starts with an upper-bound initialized to α bound = α ceil and iteratively decreases the range of proposed α. For α min , the algorithm starts with a lower bound α bound = α floor and iteratively increases the range. For both, the other bound is calculated by the following learning rule:α prop = α bound ± ζ|ρ -ρ * |.Here, ζ is the learning rate that determines how much we can adjust the step size in one tuning step. We found ζ insensitive and set ζ = .5 in all tasks. Additionally, ρ is the best acceptance rate computed from the previous iteration of the algorithm. For the first step of the algorithm, we set ρ = 0.The algorithm uses α prop , α bound to determine the range of α to test. For calculating α max , the algorithm searches in the range of[α prop , α bound ]. For calculating α min , the range is [α bound , α prop ].Given the appropriate range of α and an initial θ, we test t potential α and calculate their respective acceptance rates using Equation ( 6). Once we have computed all the acceptance rates, we set α bound to the value that resulted in the most optimal acceptance rate as determined by Equation ( 9), θ to the corresponding θ ′ , and ρ to the corresponding acceptance rate.Choice of β max , β min , ρ * , s The automatic tuning algorithm depends on an initial choice of β max , , β min , ρ * , s that enable it to automatically configure an effective hyper-parameter schedule.Here we describe the general approach to picking these values.For some target distributions, it is possible that the best possible acceptance rate with a very high β max , such as β max = .95, will not be close to the target acceptance rate ρ * . In this case, the EstimateAlphaMax algorithm will keep on decreasing the proposed α max , which will result in a very small α max . In order to avoid this behavior, we recommend starting with β max = .95, and decreasing it by .05 if the resulting α max is reasonable.We always set β min = 0.5 which is the smallest value β can take.We determine the target ρ * by starting with a value of .5 and increasing it by .05 until desirable performance metrics are obtained. While this process is essentially the same as a grid search, we note that we only needed to apply this process in the specific case of training a deep EBM on the Caltech Silhouettes dataset. For all other tasks and datasets, the target acceptance rate of ρ * = .5 was effective. We discuss the unique difficulty presented within the Caltech Silhouettes dataset in D.5.To determine the steps per cycle s, we required a similar approach to determine the optimal value.In our experiments, we only look at two different values: either s = 8, or s = 20. Having a longer cycle length tends to enable more exploitation of the target distribution, whereas having a shorter cycle enables more exploration. While we do not have an algorithm for automatically configuring this value, we were able to achieve good results across all tasks and datasets by choosing either of these two values. For more details on the resulting hyper-parameters used for each experiment, see Appendix D.Our method introduces the following hyperparameters: β max , β min , α ceil , α floor , learning rate for tuning γ, steps per cycle s, target acceptance rate ρ * , and budget B. This may seem like many additional hyperparameters, but the majority of these are introduced due to the automatic tuning mechanism and are not changed across all tasks and datasets in the paper: γ = .5, α floor = .05, α ceil = 5, β min = .5, B = 200. Thus the only hyperparameters requiring tuning in practice are β max , ρ * , and s. Note that the existing adaptive discrete sampler, any-scale sampler introduced in [Sun et al., 2023a] , has a similar number of hyper-parameters: initial step size σ, initial balancing parameter α, update rate γ, decay rate β, buffer size N , initial Hessian matrix W , and initial diagonal matrix D. Like our method, most of these hyperparameters are fixed across experiments.We conduct an ablation study to evaluate the sensitivity of our tuning algorithm to these hyperparameters choices. We choose one hyperparameter at a time to ablate and keep the rest at default values of the hyperparameters at their default setting. We run the RBM sampling experiment over multiple datasets, each for 10 random seeds, and report the average results in Figure 4 . We omit the standard error as that would harm the interpretability of the graph as many of the plots are quite close together.We can summarize the main takeaways as follows:1. The sensitivity of our algorithm to the hyperparameters depends on the dataset. For example, the sensitivity of our algorithm is low on MNIST, kMNIST, eMNIST, and Omniglot while the sensitivity is relatively high on Caltech.2. The optimal hyperparameter values depend on the dataset. For example, high values of s generally yield superior results, except for Caltech, where lower values excel. Similarly, low β max values are usually less effective, though Caltech is an exception, showcasing decent outcomes. In general, the hyperparameter values we selected to generate the final results in the experiment section were the ones that generalized across the datasets.3. For each ablation, the values tested demonstrate reasonable results when compared with the baselines. While not all hyperparameter values result in equally competitive performance, all of them outperform the Gibbs-With-Gradient sampler Grathwohl et al. [2021] . This demonstrates that our method performs well with a wide range of hyperparameters and can achieve even better performance with careful hyperparameter tuning.In conclusion, we believe these results demonstrate that our algorithm is relatively robust to choice in hyperparameters.Require: α ceil , α floor , β max , β min , steps per cycle s, steps to take without MH correction l = 50, steps to take with MH correction l MH = 50, initial state θ 1: for step i in range(l) do We note that all configurations to exhibit convergence to the ground truth as indicated by the maximum mean discrepancy (log MMD), albeit with varying convergence speeds. In some cases, specific hyper-parameter configurations are able to achieve better performance than what we report in the RBM sampling experiment. Overall, we can observe that our algorithm is reasonably robust to various hyper-parameter configurations as it will still demonstrate convergent behavior towards the ground truth.Energy Based Models (EBMs) are a class of generative models that learn some unnormalized distribution over a sample space. As discussed in Hinton [2002] , these models can be trained via the following Maximum Likelihood objective:EQUATIONThe gradient updates for this loss function are known to be as follows:EQUATIONWhile the expectation on the left is straight forward to calculate given a dataset, calculating the right expectation is not as clear. Here we will mention the two methods that are relevant towards our experiments with EBMs.Contrastive Divergence (CD) In order to estimate the second term, we initialize some sampler using the x in the first term and run it for a set number of sampling steps. For a more detailed description, refer to Hinton [2002] .The expectation on the right can be calculated using samples from a persistent Markov Chain that approximates the true distribution Tieleman [2008] .Require:Step size schedule {α max , α 1 , . . . α min }, β max , β min , number of proposals per step t = 10, initial state θ, target acceptance rate ρ * 1:β floor = β min , β ceil = β max 2: β-sched ← {β max } 3: for i in {1, 2, . . . s -1} do 4:proposed-params ← LinSpace(β floor , β max , t) ▷ Create t potential balance parameters for index i in the schedule β ceil ← β i ▷ Shrink the range of potential balancing parameters by using assumption β i > β i+113: θ = θ ′ correspending to β i 14: end for 15: β-sched.append(β min ) 16: return β-sched Instead of resetting the chain each training iteration, we maintain a buffer of the generated samples that we use to calculate the second expectation. This method relies on the intuition that the model distribution does not vary too widely within one iteration. Using the intuition provided by [Du and Mordatch, 2019] , we can view this process as updating the model parameters ϕ to put more weight on true samples and less weight on fake samples. By doing so, the model will in turn generate samples that closer to those from the true distribution.Main Idea We can apply the ACS algorithm combining the automatic tuning of the cyclical schedule with the original PCD learning algorithm. Our goal in doing so is to improve PCD through better characterization of the entire model distribution. During training, we can view PCD as adjusting the model parameters to "push down" the probability of samples from the model distribution while "pushing up" samples from the true data distribution. Because our sampling method is able to explore the model's distribution more effectively than other samplers, we can adjust more regions of the model distribution at a quicker rate than previous sampling methods, which should improve the quality of gradient updates and thus lead to better model parameters. We adapt ACS to work within PCD by having the step size depend on the training iteration as opposed to the sampling iteration, with the corresponding α, β pair being used for all the sampling steps within the iteration. We include the complete learning algorithm in Algorithm 6.Cyclical Scheduling We find that the learning task requires a different approach to the cyclical scheduling than the sampling task. Rather than having a relative equal amount of exploration and exploitation, we find that it is more effective to use a cyclical schedule biased towards exploitation. However, exploration is still important as it enables the model to better represent the distribution as a whole rather than a few local modes. Given this, we construct a cyclical schedule consisting of one iteration that uses α max , β max with the rest using α min , β min .Tuning One of the advantages of using the simplified cyclical schedule is that it only requires two pairs of hyper-parameters to be optimized. Thus we can leverage the EstimateAlphaMax and EstimateAlphaMin algorithm to both tune the respective α, β pair while also updating the persistent buffer. Not only does this reduce the additional overhead of the tuning component, but it allows the hyper-parameters to adapt to the changing EBM distribution. end if19: Construct Q = Q α,β (•|X f ) using (4) 20:for sampling step in range(S big ) do 21:X ∼ Q(•|X f ) 22: if i mod s = 0 then 23: X f ← X 24:continue ▷ If i is the first step of the cycle, omit the MH correction 25:end if 26:X f ← X with acceptance probability as calculated in (6) 27:end for 28:end if 29: Calculate E x∼p ϕ [∇ ϕ E ϕ (x)] using X f 30: Calculate E x∼pdata [∇ ϕ E ϕ (x)] using X t 31: ∇L(ϕ) = E x∼p ϕ [∇ ϕ E ϕ (x)] -E x∼pdata [∇ ϕ E ϕ (x)] ▷ϕ = ϕ -γ ϕ ∇L(ϕ) 33: i+ = 1 34:end for 35: end whileWe define the problem setting in more detail. We have a target that is of the formπ(θ) = 1 Z exp(U (θ)).We consider the proposal kernel asQ α,β (θ ′ |θ) ∝ exp β∇U (θ) T (θ ′ -θ) - 1 2α ∥θ ′ -θ∥ 2and consider the transition kernel asp(θ ′ | θ) = π(θ ′ )Q α,β (θ | θ ′ ) π(θ)Q α,β (θ ′ | θ) ∧ 1 Q α,β (θ ′ | θ) + (1 -L(θ)) δ θ (θ ′ )where δ θ (θ ′ ) is the Kronecker delta function and L(θ) is the total acceptance probability from the point θ withL(θ) = θ ′ ∈Θ π(θ ′ )Q α,β (θ|θ ′ ) π(θ)Q α,β (θ ′ |θ) ∧ 1 Q α,β (θ ′ |θ).We also defineZ α,β (θ) = x∈Θ exp β ∇U (θ) T (x -θ) - 1 2α ∥x -θ∥ 2which is the normalizing constant for the proposal kernel.C.1 Proof of Lemma 5.3Proof. By including the balancing parameter, we start by noting thatEQUATIONConsider the term,EQUATIONSubstituting ( 15) in ( 14), the numerator ofQ α,β (θ, θ ′ ) β∇U (θ) T (θ ′ -θ) - 1 2α ∥θ ′ -θ∥ 2 =β (-U (θ) + U (θ ′ )) - β 2 (θ -θ ′ ) T 1 0 ∇ 2 U ((1 -s)θ + sθ ′ ) ds (θ -θ ′ ) - 1 2α (θ -θ ′ ) T I(θ -θ ′ ) =β (-U (θ) + U (θ ′ )) - 1 2 (θ -θ ′ ) T β 1 0 ∇ 2 U ((1 -s)θ + sθ ′ ) ds + 1 α I (θ -θ ′ ).From Assumption 5.1 (U is M -gradient Lipschitz), we haveβ 1 0 ∇ 2 U ((1 -s)θ + sθ ′ ) ds)(θ -θ ′ ) + 1 α I ≥ 1 α -βM ISince α < 1/βM , the matrix 1 2α -βM I is positive definite. We note thatEQUATION) Z α,β (θ) = x∈Θ exp β ∇U (θ) T (x -θ) - 1 2α ∥x -θ∥ 2 = x∈Θ exp -β (U (θ) -U (x)) - 1 2 (θ -x) T (β 1 0 ∇ 2 U ((1 -s)θ + sx) ds)(θ -x) + 1 α I)(θ -x) .This can be seen asπ(θ)Q α,β (θ ′ |θ) = 1 Z Z α,β (θ) exp β (U (θ) + U (θ ′ )) -(θ ′ -θ) T 1 2α I + β 2 1 0 ∇ 2 U ((1 -s)θ + sθ ′ )ds (θ ′ -θ) .Since Assumption 5.2 holds true in this setting, we have an m > 0 such that for any θ ∈ conv(Θ)-∇ 2 U (θ) ≥ m I. From this, one notes that exp -βU (θ) - 1 2 1 α -β m diam(Θ) 2 x∈Θ exp (βU (x)) ≤ Z α,β (θ) ≤ exp (-βU (θ)) x∈Θ exp (βU (x))where the right-hand side follows from the fact that α < 1/(βM ). Therefore,Z α,β (θ) Z α,β (θ ′ ) ≥ exp {β (-U (θ) + U (θ ′ ))} exp 1 2 1 α -βm diam(Θ) 2 Also note that Q α,β (θ ′ |θ) = exp β (-U (θ) + U (θ ′ )) -(θ -θ ′ ) T 1 2α I + β 2 1 0 ∇ 2 U ((1 -s)θ + sθ ′ ) (θ -θ ′ ) θ ′ ∈Θ exp β (-U (θ) + U (θ ′ )) -(θ -θ ′ ) T 1 2α I + β 2 1 0 ∇ 2 U ((1 -s)θ + sθ ′ ) (θ -θ ′ ) ≥ exp β ⟨∇U (θ), θ ′ -θ⟩ -1 2α ∥θ -θ ′ ∥ 2 θ ′ Θ exp {β (-U (θ) + U (θ ′ ))} .We also note that-β ⟨∇U (θ), θ ′ -θ⟩ + 1 2α ∥θ -θ ′ ∥ 2 = β ⟨-∇U (θ) + ∇U (a), θ ′ -θ⟩ + β ⟨-∇U (a), θ ′ -θ⟩ + 1 2α ∥θ -θ ′ ∥ 2 ≤ β ⟨-∇U (θ) + ∇U (a), θ ′ -θ⟩ + β ⟨-∇U (a), θ ′ -θ⟩ + 1 2α diam(Θ) 2 ≤ β ∥-∇U (θ) + ∇U (a)∥∥θ ′ -θ∥ + β ∥∇U (a)∥∥θ ′ -θ∥ + 1 2α diam(Θ) 2 ≤ β∥ -∇U (θ) + ∇U (a)∥diam(Θ) + β∥∇U (a)∥diam(Θ) + 1 2α diam(Θ) 2 ≤ βM + 1 2α diam(Θ) 2 + β∥∇U (a)∥ diam(Θ).Combining, we getp(θ ′ |θ) ≥ ϵ β,α exp {βU (θ ′ )} θ ′ Θ exp {βU (θ ′ )} where ϵ β,α = exp - 1 α + βM - β m 2 diam(Θ) 2 -∥∇U (a)∥ diam(Θ) .Synthetic Distribution In order to construct a distribution that is easy to visualize, we first must define a few experiment parameters. We must define the space between the modes, the total number of modes, and the variance of each mode σ. For convenience, we have the number of modes as 25, which is a perfect square. We define the space between modes as 75, and the variance for each mode σ 2 as .15. Given this, we can calculate the maximum value for each coordinate as follows:MaxVal = ( √ NumModes + 1) * SpaceBetweenModesWe can calculate the coordinate value for each mode µ i,j as follows:µ i,j [0] = MaxVal √ NumModes + 2 (i + 1) µ i,j [1] = MaxVal √ NumModes + 2 (j + 1)Sampler Configuration Our goal in this experiment is to demonstrate how gradient-based samplers typically behave when faced with a distribution with modes that are far apart. In order for this experiment to be meaningful, it is important that the representation of each sample respect the notion of distance between the integer values. For this reason, we cannot use a categorical distribution or represent each coordinate with a one-hot encoding, as every sample in this representation would be within a 2-hamming ball of every other point. In order to determine the step sizes for the baselines, we tune each until we reach an acceptance rate around .574. For DMALA, this ends up being around α = 53. For the any-scale sampler, we set the initial step size to be the same and use their implemented adaptive algorithm.For the cyclical sampler, we set α max = 1575, α min = 3, and steps per cycle s = 20. Because the goal of the experiment is to demonstrate the need for larger step sizes along with smaller step sizes, we do not use the automatic tuning algorithm on this example as restricting the space to be ordinal changes the optimal setting for α ceil . In most practical cases, the samples would be represented by a categorical or binary form, which the proposed tuning algorithm is able to handle as demonstrated by the performance on real data distributions.Uneven Multi-modal Distributions Not only does a cyclical step-size enable more accurate sampling in highly multi-modal distributions, but it is also able to handle distributions where the modes are weighted unevenly. This problem is more difficult since this requires not only exploring all the modes of a distribution, but ensuring that the less likely modes are not over represented in the generated samples. We provide a visual comparison between the target distribution, the estimated distribution from DMALA, and the estimated distribution from ACS in Figure 5 .Since the modes may not be clear due to the nature of this specific problem, we also include a quantitative comparison between DMALA and ACS in Table 3 by computing the KL divergence between the estimated distribution and the target distribution in addition to the average energy of the generated samples. Through both Table 3 and Figure 5 , we observe that a cyclical step-size enables accurate sampling from uneven multi-modal distributions.Table 3 : Quantitative comparison of sampler performance on uneven multi-modal distribution. ACS retains the ability to accurately capture all the modes within the distribution despite the uneven weighting.DMALA ACS KL Divergence 0.70 0.13 Average Energy -2.66 ± 1.68 -3.39 ± 1.63 Furthermore, it is interesting to observe that generating more highprobability samples does not necessarily correspond to accurate sampling, highlighting the difference in goals between generating very likely samples and accurately sampling the target distribution.RBM Overview We will give a brief overview of the Block-Gibbs sampler used to represent the ground truth of the RBM distribution. For a more in-depth explanation, see Grathwohl et al. [2021] . Given the hidden units h and the sample x, we define the RBM distribution as follows:log p(x, h) = h T W x + b T x + c T -log Z (19)As before, Z is the normalizing constant for the distribution. The sample x is represented by the visible layer with units corresponding to the sample space dimension and h represents the model capacity. It can be shown that the marginal distributions are as follows:p(x|h) = Bernoulli(W x + c) p(h|x) = Bernoulli(W t h + b)The Block-Gibbs sampler updates x and h alternatively, allowing for many of the coordinates to get changed at the same time, due to utilizing the specific structure of the RBM model.Experiment Setup Similar to the experimental setup of Zhang et al. [2022a] , we use RBM models with 500 hidden units and 784 visible units. We adopt the same training protocol, except we train the RBM with 100 steps of Contrastive Divergence as opposed to 10. We also train the models for 1000 iterations as opposed to a single pass through the dataset. We find that this enables the RBMs to generate more realistic samples. We include the generated images in Figure 6 to demonstrate that these models have learned the dataset reasonably well. Sampler Configuration For GWG, we use the same settings as Grathwohl et al. [2021] , for DMALA, we set step size to .2, and for AB we use the default hyper-parameters for the first order sampler.For ACS, we use ρ * = .5, β max = .95, ζ = .5, cycle length s = 20 for all the datasets. We also fix the total overhead of the tuning algorithm to 10% of the total sampling steps.Escape from Local Modes In addition to using the same initialization as Zhang et al. [2022a] , Grathwohl et al. [2021] , we extend the experiment to measure the ability of a sampler to escape from local modes. We initialize the sampler within the most likely mode, as measured by unnormalized energy of the RBM. Samplers that are less prone to getting trapped in local modes will be able to converge quicker to the ground truth, as measured by log MMD. We include the performance of the various samplers across 11 random seeds in 7. ACS demonstrates superior robustness to mode-specific initialization due to its capability to escape from local modes. ACS demonstrates more robust sampling behavior across the datasets than other methods, as evidenced by superior convergence on all datasets except KMNIST. We do note that ACS performance is still competitive on KMNIST with the added benefit of a smaller standard error.We found that a visual inspection of the generated images demonstrates the ability of ACS to escape local modes. We include the generated images in Figure 8 . Figure 8 : Images sampled from RBM trained on MNIST when the sampler is initialized to most likely mode. ACS is able to generate a diverse range of digits, demonstrating its ability to escape from modes. It should also noted that while AB is able to generate a diverse range of digits as well, the images are slightly less clear than those generated by ACS.We can make two primary inferences from the generated images: the first being that ACS is able to escape from local modes and explore the distribution as a whole, as demonstrated by the wide range of generated images; and that ACS does not compromise on the ability to characterize each mode as evidenced by the quality of generated samples.Sampling Speed While the run time can vary depending on the specific implementation of a given sampling algorithm, we illustrate the efficiency of ACS in Figure 9 . ACS is able to outperform DMALA in terms of convergence with respect to time, even including the overhead of the tuning algorithm. Figure 9 : Log MMD for the RBM sampling task across time for both DMALA and ACS for the kMNIST dataset. We observe that while ACS is offset slightly in the beginning due to the initial tuning algorithm, it is quickly able to demonstrate superior convergence.In order to train the EBMs, we use Gibbs-with-Gradient to sample the EBM distribution during PCD, following the same training protocol as Grathwohl et al. [2021] . We train these models for 50,000 iterations total with 40 sampling steps per iteration and use the parameters corresponding to the best log likelihood scores on the validation dataset.Experimental Design For each of the trained models, we evaluate the samplers based on how quickly the average energy of the generated samples rises. This gives an estimate of the speed at which a sampler is able to reach a stationary distribution.Sampler Configuration For GWG, we use the same settings as Grathwohl et al. [2021] , for DMALA, we set step size to .2, and for AB we use the default hyper-parameters for the diagonal variant of the AB sampler. We choose this variant as this is what they evaluate for their experiments when measuring mixing speed of samplers on EBMs.For ACS, we use ρ * = .5, β max = .8, ζ = .5, cycle length s = 20 for all the datasets. As in RBM Sampling, we fix the total overhead of the tuning algorithm to 10% of the total sampling steps.Sampler Performance It is worth commenting on the similarity in performance between ACS and DMALA when sampling from Caltech. We find that when sampling from the EBM trained on the Caltech dataset, ACS finds a α max similar to α min , thus making the ACS sampler similar to DMALA for this specific case. We hypothesize that small step sizes are most effective for this dataset. The results in Figure 3 demonstrate that ACS can handle such cases automatically: while the step size for DMALA must be hand-tuned, the ACS method can automatically adapt to a suitable step size schedule.We include the images generated by ACS when sampling from deep EBMs in Figure 10 . Experiment Design We use the same RBM structure as the sampling task, with 500 hidden units and 784 visible units. However, we apply the samplers of interest to the PCD algorithm introduced by Tieleman [2008] . The model parameters are tuned via the Adam optimizer with a learning rate of .001.In order to evaluate the learned RBMs, we run AIS with Block-Gibbs as the sampler to calculate the log likelihood values for the models Neal [2001] . We run AIS for 100,000 steps, which is adequate given the efficiency of Block Gibbs for this specific model.For DMALA, we use a step size of .2. For the ACS algorithm, we set β max = .9, ρ * = .5 for all the data-sets. We do modify the number of cycles for each data-set as different distributions require different amounts of exploration and exploitation. We use cycle length of 8 for MNIST, eMNIST, and kMNIST; we use 20 for Omniglot and Caltech silhouettes. This difference reflects the specific needs for each dataset in terms of exploration and exploitation -more complex datasets tend to need longer cycles in order to better exploit each region, while simpler datasets tend to need shorter cycles in order to capture all the modes of the learned distribution. In Figure 11 , we show the samples generated from AIS for 100,000 steps as opposed to the persistent buffer as this forms a longer MCMC chain, thus giving a better visual of what the learned distribution represents.In order to ensure that the overhead for the tuning algorithm does not add to the overall computational cost, we spread out the computations of the EstimateAlphaMin algorithm throughout the training process. We keep a running list of α min and set α floor to be one standard deviation below the mean of this list. By doing this, we start closer to what the ideal α min . For EstimateAlphaMax, we simply call the tuning function every 50 cycles containing 50 training iterations, with α ceil = 5. As the initial step does not use the Metropolis-Hastings correction and has half the sampling steps, the budget for each call of EstimateAlphaMax can be seen as coming in part by the computation saved.We include the AIS results for RBMs trained with different sampling methods in Table 4 . We see that ACS achieves superior log likelihood results when compared to other sampling methods across all datasets. Furthermore, the AIS results are consistently close to those achieved by Block-Gibbs, which can be considered close to ideal for RBMs since it leverages the known structure of the model.We include the generated images from the RBMs trained using different samplers in Figure 11 . In general, the images generated from the ACS-trained RBM capture more modes than other methods, except for the Caltech Silhouettes dataset. In particular, all the methods struggle to generate reasonable images for this dataset. We hypothesize that this is due to the increased complexity of the distribution relative to the other datasets -Caltech Silhouettes is composed of the silhouettes from real objects, whereas the other datasets are hand-written symbols. This hypothesis is supported by the generated images in Figure 6 , where the images generated when using Block-Gibbs on Caltech Silhouettes also seem less reasonable than the samples obtained from different datasets. Since Block-Gibbs is the best sampler for this specific model as it leverages the known structure of the RBM, this appears to be unavoidable as a result of limited model capacity. This motivates our experiments with deep convolutional EBMs, where we can understand how our method does when using a model architecture with sufficient capacity.Experiment Design We use the same EBM model architecture as Zhang et al. [2022a] , Grathwohl et al. [2021] and follow the same experimental design, with the only change being to the number of sampling steps alotted for each sampler.In order to determine the number of sampling steps that we could use for ACS-PCD, we tested different sampling steps. For Static/Dynamic MNIST and Omniglot, we found that we only needed to use 10 sampling steps to achieve good models. However, we observed divergence when training models on Caltech. In order to determine what number of sampling steps to use for ACS, we do a grid search over the number of sampling steps and ρ * with all other values remaining the same. We test sampling steps of 10, 20, and 30; and we use ρ * .5, .6, .7, .8. We decide which hyper-parameters to use based on when training diverged the latest; and we use the best model parameters as indicated by validation log likelihood. We use 10 sampling steps for Static, Dynamic MNIST, and Omniglot, while we found 30 was the minimum we could use for Caltech Silhouettes and obtain reasonable results. We apply this number of sampling steps for both DMALA and ACS to demonstrate how the methods compare when facing a similar budget constraint.In order to evaluate these learned models, we use the same evaluation protocol as Zhang et al.[2022a], Grathwohl et al. [2021] . We run AIS for 300,000 iterations using Gibbs-With-Gradient as the evaluation sampler. By following the same experimental design as previous works, we can draw meaningful comparisons from previous results in Grathwohl et al. [2021] .Sampler Configuration For DMALA, we use a step size of .15 as used in Zhang et al. [2022b] . For ACS, we use 200 sampling steps for EstimateAlphaMax and EstimateAlphaMin. For Static MNIST, Dynamic MNIST, and Omniglot, we set the algorithm to tune α max and α min every 25 cycles, where each cycle has 50 training iterations. The additional overhead of this is 16,000 extra sampling steps, which is a 3.2% of the total budget of 500,000 sampling steps. For Caltech Silhouettes, we have to adapt every 10 cycles with the same number of training iterations. This results in 40,000 additional sampling steps due to the tuning algorithm. For this specific dataset, because we use 30 sampling steps, the additional cost is 2.6% of the total sampling steps 1,500,000.In terms of the final parameters for cycle length and sampling steps, we find that we can use the same ρ * across all datasets, with the exception of Caltech Silhouettes. For Static/Dynamic MNIST and Omniglot, we were able to use ρ * = .5 and For this dataset, we found good results by setting ρ * = .7.We hypothesize that the need for a higher acceptance rate is due to the fundamental difference between Caltech Silhouettes and the other datasets, as previously mentioned. Because Caltech Silhouettes contain samples are derived from real objects, they are more complex than the hand-written figures.In addition to the empirical results in Table 1 , we provide some qualitative data in the form of the generated images from the PCD buffer when using ACS. We choose to include the buffer images for this experiment as the chain from the persistent buffer is much longer than the chain from AIS due to the increased training duration: the chain from AIS is obtained using 300,000 sampling steps whereas the persistent buffer is obtained from 500,000 sampling steps on Static/Dynamic MNIST and Omniglot, 1,500,000 sampling steps for Caltech Silhouettes. By visualizing the generated images from the longer chain, we get a better understanding of the quality of the trained distribution. We put the images in Figure 12 .We also observe that this behavior is not unique to ACS and does occur when Gibbs-With-Gradient and DMALA are used with 40 sampling steps as indicated. Instability is common when training deep EBMs, and this is most likely why the original experimental design included check-pointing throughout the training process as well as comparisons based on the validation set. We also note that despite this behavior, the trained models are able to generate fairly realistic images. We present the images from the PCD buffer for ACS below in figure.When the images in Figure 12 are taken in context of the improvements in log likelihoods as presented in Table 1 , the results indicate the benefits of using ACS when learning multi-modal discrete distributions. Experimental Design For both datasets, we sample 100 sentences randomly and mask 50% of the tokens. We use a pretrained RoBERTa model available through the Hugging Face API Liu et al. [2019] . We run 25 separate chains for each example and take the final state of each chain to be a sample. We then take the top-5 most likely samples and use these for empirical comparisons.We define the energy function the same as in Zhang et al. [2022a] . Let us define a sentence of length d θ = {θ 1 , θ 2 , . . . θ d }, where θ i is a one hot vector over vocabulary V . Let M ⊂ {1, 2, . . . d} the set of indices we wish to sample. We define the function f (θ i ¬i ) to be the log probability distribution over V for the i position conditioned on all other positions. Given this, we define the energy function for the sentence θ to be as follows:EQUATIONSampler Configuration For DMALA, we tune the step-size to achieve an acceptance rate of 50%, which ends up being α = .5. For ACS, we use a cycle length of 20. We use the same hyper-parameters for the tuning algorithms as in previous tasks, demonstrating that our algorithm can be applied across domains and tasks with little modification. We include example generations for both ACS and DMALA in Figure 13 .In Table 2 , we observe that the ACS generations have higher perplexity than the DMALA generations. While perplexity is a popular means of evaluating language generations, it is important to recognize that perplexity is based on the likelihood of the generation under the language model. This metric is biased towards frequent patterns and does not account for diverse modes. Therefore, it does not completely align with the goal of MCMC, which is to accurately characterize the target distribution.Minimizing the average perplexity of the sample corresponds to maximizing the average likelihood of the generations, which can be at odds with the goal of accurately capturing the target distribution. We illustrate this in Appendix D.1, where we compare the performance of DMALA and ACS when sampling from a uneven multi-modal distribution.• The comedy that follows feels hackneyed or just plain crude, calculated to provoke shocked stares, without opening up to a deeper truth.• This comedy could either be hacky, or just plain crude, calculated to provoke shocked curiosity, without opening up to a deeper insight.• A comedy that feels slightly hacky, or just plain crude, calculated to achieve shocked results, without coming up with a deeper message.• The comedy was either unnecessarily hacky, or just plain crude, calculated to create shocked humor, without leading up to a deeper plot.• Most comedy is flat or hack-ish or just plain crude, calculated to evoke shocked laughs, without opening up to a deeper audience.(a) ACS• This comedy is either plain hacky, or just plain crude, calculated to provoke shocked discussion, without linking up to a deeper message.• And comedy that can be hacky, or just plain crude, calculated to provoke shocked discussion, without opening up to a deeper topic.• Modern comedy can be deliberately hacklish, or just plain crude, calculated to provoke shocked disbelief, without opening up to a deeper meaning.• Simple comedy has all things hacky, or just plain crude, calculated to be shocked away, without opening up to a deeper meaning.• A comedy usually ranges from hacky, or just plain crude, calculated to provoke shocked surprise, without opening up to a deeper reality.(b) DMALA Justification: We include the specific GPU at the beginning of Appendix D.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).Justification: There are no study participants.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 


Decision: This paper proposes an automatic cyclical scheduling, designed for efficient and accurate sampling in multimodal discrete distributions. The authors prove a non-asymptotic convergence guarantee in general discrete distributions. The method is also verified by experiments.

All reviewers agree to accept this paper because (1) using a cyclic step size schedule is well motivated, (2) an automated tuning scheme is used, making it easy to configure the two cyclic schedules, (3) strong empirical results, (4) non-asymptotic convergence guarantees, etc. I agree with the opinions and recommend an accept.

Reviews: review: paper_summary: paper_summary: This paper proposes a new discrete sampling method called ACS that addresses a common problem for existing gradient-based  approach where they are susceptible to becoming trapped in local modes. ACS combines local-balancing proposals with a cyclic step size to balance local exploitation and global exploration; it is in essence an extension of cyclic stochastic-gradient MCMC to discrete distributions.  To ensure proposals are still balanced with a varying step size, ACS uses a cyclic balancing
schedule along with an automatic tuning scheme to easily adapt the schedules.  Non-asymptotic convergence guarantees are provided.  Results demonstrate ACS to outperform prior approaches for sampling from energy based models, training RBMs, and text-infilling.
strength_weakness: strength_weakness: - Experimental results for RBM have discrepancies wrt results reported in previous papers.  In particular the ranges differ for average energy from other papers and the curves for log MMD for ACS show unexpected curvature.  Please respond to the questions in the section below to clarify.
- Inadequate discussion of text-infilling results.  It is not clear why higher perplexity and diversity is good for ACS since the goal is to be able to efficiently sample from a target discrete distribution.
- Error bars missing for results in Figure 3 and Table 1.

Typos & formatting:
- Inline citation format should show author name
- Figure 3 caption: cpnvergence -> convergence
questions: questions: - Why is average energy in Figure 3 negative?  Also, the results for ACS on dynamic_mnist and omniglot are unexpected with better performance on fewer iterations before converging.
- Why is the scale for average energy in Figure 3 different from that reported in Figure 4 of the [AB paper](https://openreview.net/pdf?id=lEkl0jdSb7B)?
limitations: limitations: Limitations are adequately discussed.
paper_summary: paper_summary: The paper introduces a novel method for sampling from multimodal discrete distributions, which presents an innovative approach to address the challenge of local modes trapping in gradient-based discrete sampling, together with non-asymptotic convergence guarantee and empirical validation of the proposed method.
strength_weakness: strength_weakness: 1. In the experiments, there seems to be no error bars in Figure 1 and Table 1. 
2. The quality of the samples seems worse than DMALA in Table 2 and Figure 12. Is it possible that the proposed method might sacrifice sample quality to achieve higher diversity?
questions: questions: 1. The paper mentions the proof is not consistent with the specific tuning algorithm used in the experiments. Could you elaborate on the reasons?
limitations: limitations: N/A
paper_summary: paper_summary: The paper proposes a solution to the challenge of sampling from high-dimensional discrete spaces, where conventional discrete samplers often get trapped in local modes. To address this, the authors introduce a discrete Langevin sampler with automatic cyclical scheduling. This method comprises three components: a cyclical step size schedule, a cyclical balancing schedule, and an automatic hyperparameter tuning scheme. The authors provide theoretical guarantees for non-asymptotic convergence and inference, and extensive experiments demonstrate the method's superiority in sampling complex multimodal discrete distributions.
strength_weakness: strength_weakness: - My primary concern lies in the complexity of the proposed methods. The automatic schedule tuning scheme appears to be quite time-consuming, particularly the grid search required for the balancing parameters \beta_i, demanding significant computational resources.
- Another concern pertains to the theoretical assumptions underlying the analysis. The non-asymptotic convergence of the proposed samplers relies on the strong convexity of the negative energy function, an assumption that may not hold in practical deep EBM scenarios. Despite this potentially restrictive condition, the analysis provides valuable insights, and empirically, the proposed method demonstrates effective performance, as supported by extensive studies.
questions: questions: - In equations 8 and 9, how do you estimate the acceptance rate A? Is it estimated by averaging across the training batch? If so, that implies the complexity would increase to n*s times compared to the original DMALA samplers at each step, where n is the number of grids for parameters \beta_i and s is the number of sampling steps per cycle. This would be highly time-consuming.
- Could you provide a comparison of the running times between the proposed ACS sampler and the DMALA samplers?
- Could you elaborate on the rationale behind setting the target acceptance rate $\rho^*$ to 0.5 in your experiments? What are the implications of setting it to 1 or 0.234 instead?
limitations: limitations: Despite the strong assumptions in the theoretical analysis and the perceived complexity of the algorithm, this paper presents a substantial contribution to the field by offering a well-explained, theoretically sound, and empirically validated approach to enhancing discrete sampling in high-dimensional spaces.
paper_summary: paper_summary: The paper presents a novel gradient-based algorithm to sample from complex multimodal discrete distributions based on differentiable energy functions. Overall, the method is based on locally balanced proposals, previously introduced, and instantiates it with parametrized functions and a cyclical schedule for the "learning rate" that promotes the alternation of modes discovery and modes refinement. 
The authors also provide an algorithm to automatically tune the introduced parameters based on an input acceptance rate and initial and final balancing parameters.
The paper features a theoretical analysis that includes concrete convergence rates under (somewhat restrictive) assumptions and finishes with some experiments in learning and sampling from RBMs and EBMs and finally in text infilling with masked language models.
strength_weakness: strength_weakness: - I think too much of the paper is in the appendix, and the information presented in the paper is not sufficient to fully follow the logical flow (see also points below). In my opinion, some details regarding the development of the method and the theory could be moved to the appendix to make more space for both preliminaries (like how are these methods used in practice) and developing more intuition.  
- the algorithm is fairly complex and the paper fails in provide enough intuition for some parts of its functioning. Like the authors claim that "it is fairly easy to choose initial and final balancing factors", but do not elaborate why (in the main paper)
- the only "real-world" experiment is only sketched in the main paper, not providing enough information to appreciate the task. What is the precise difficulty here? Since roberta is a masked language model I believe it is possible to derive "pseudo-distributions" like in [1]  
- the theoretical analysis seems to require strong assumptions that may be violated in compelling real-world use cases (like LLMs)


[1] Hennigen, Lucas Torroba, and Yoon Kim. "Deriving language models from masked language models." arXiv preprint arXiv:2305.15501 (2023).
questions: questions: - Can you please discuss in which cases the locally concave hypothesis holds for realistic models such as LLMs? 
- Even if this is probably more of a sanity check, I'd like to see how the method behave on unimodal, or mildly multimodal distributions 
- Can you please describe some other applications of the method to real-world problems (like text infilling)
limitations: limitations: Partially discussed, see questions

Comments: title: Post author response
Thanks
comment: Thank you for responding to my questions.  I have changed my score from 4 to 5.
Thanks for the authors' response. It has solved my concern about the complexity. I keep my opinion that this work is good to be in, and I highly recommend including the discussion of complexity in the camera ready.
Thank you for the response and clarifications. I will keep my score.
Dear authors. Thank you very much for your rebuttal. I appreciate the additional experiments clarifications on the infilling task. I keep my opinion that this work is a valid addition to the conference.