We sincerely appreciate your careful review of our work and your positive feedback. Below please find our point-by-point responses.

Weakness 1. We appreciate the perspective that the depth and technical nature of our results might be well-suited to a journal format. However, we expect that our contributions will be of strong interest to some conference audience as well, since quite a few similar works [1, 2] have been presented at conferences like NeurIPS and ICLR. We hope our paper will likewise contribute meaningfully to ongoing discussions in this area.

Weakness 2. We thank the reviewer for the thoughtful question. Theorem 5 establishes the convergence guarantees for STMH when sampling from a mixture of Gaussians with shared covariance matrix. This result can be naturally extended to target distributions that are sufficiently close to such Gaussian mixtures, as shown in [1]. In such cases, the time complexity would additionally depend on closeness between the actual distribution and the Gaussian mixture approximation. We also anticipate that similar techniques can be adapted to handle mixtures of log-concave distributions that are translates to each other, or more broadly, distributions that are well-approximated by such mixtures.

Additionally, there are realistic applied settings where our assumptions may approximately hold, making the proposed framework potentially applicable. For instance, [3] showed that the posterior distribution for Bayesian variable selection (which is one of the most common model selection problems in statistics) can be approximated by a mixture of Gaussians under certain conditions.

Another relevant application arises in the context of (Bayesian) deep neural networks. Consider the case when the target is $e^{-f}$. In the context of DNN, $f$ is the loss composed with a neural network function. In this setting, the complexity of sampling will depend on the structure of $f$. If the widths of the neural networks tend to infinity, several works have shown that $f$ satisfied PL-type condition and hence by recent works like [4, 5, 6] $e^{-f}$ satisfies Poincare or Log-Sobolev type isoperimetric inequalities. Therefore, algorithms like Langevin Monte Carlo, proximal samplers or MALA could be used to sample from the posterior. However, in the practical regime where the width is finite, the corresponding $f$ does not satisfy any isoperimetric inequality and so $e^{-f}$ is highly multimodal, and the proposed algorithm is likely to be useful.

Question 1. This has been addressed in the response to Weakness 2.

Limitation. We thank the reviewer for pointing out the more relevant references related to sequential Monte Carlo methods. We will incorporate them into the revised version of the manuscript.

[1] H. Lee, A. Risteski, and R. Ge. "Beyond log-concavity: Provable guarantees for sampling multi-modal distributions using simulated tempering Langevin Monte Carlo." Advances in neural information processing systems 31 (2018).

[2] O. Chehab, A. Korba, A. J. Stromme, and A. Vacher. "Provable convergence and limitations of geometric tempering for Langevin dynamics." In The Thirteenth International Conference on Learning Representations (2025).

[3] I. Castillo, J. Schmidt-Hieber and A. Van der Vaart. "Bayesian linear regression with sparse priors." The Annals of Statistics 43.5 (2015): 1986-2018.

[4] Y. Gong, N. He and Z. Shen. "Poincare inequality for local log-Polyak-Łojasiewicz measures: Non-asymptotic analysis in low-temperature Regime." arXiv preprint arXiv:2501.00429 (2024).

[5] A. Y. Chen and K. Sridharan. "Optimization, isoperimetric inequalities, and sampling via Lyapunov potentials." arXiv preprint arXiv:2410.02979 (2024).

[6] S. Chewi and A. J. Stromme. ``The ballistic limit of the log-Sobolev constant equals the Polyak-Łojasiewicz constant." arXiv preprint arXiv:2411.11415 (2024).

First of all, we sincerely appreciate your careful review of our manuscript and your thoughtful feedback. Below please find our point-by-point responses.

Weakness 1. We agree that Assumption 2 can be difficult to verify in general and typically needs to be assessed on a case-by-case basis. However, since our analysis is based on the restricted spectral gap, it is sufficient to verify this assumption only within a subset of the state space—often a compact set—rather than globally. This makes the assumption more tractable in practice. Specifically, for other Metropolis-Hastings samplers such as MALA (Metropolis-adjusted Langevin algorithm), inequality (2) in Assumption 2 essentially reduces to bounding the acceptance and proposal probabilities. Since we are focusing on compact subsets of the state space, obtaining such bounds is usually feasible. Inequality (3) in Assumption 2 essentially assumes a bound on the restricted spectral gap of each component chain. Again, on compact spaces, there are standard techniques for obtaining such bounds, such as the path methods of [1]. Whether the resulting mixing time bound is sharp has to be checked on a case-by-case basis.

Regarding the dependence on dimension $d$, we acknowledge that the constants in our complexity bound scale exponentially with $d$ when using the random walk Metropolis–Hastings (RWMH) kernel. This is a known limitation of RWMH in high-dimensional settings and has already been discussed in Remark 5. We will make it more explicit in the revised version.

[1] Yuen, Wai Kong. "Applications of geometric bounds to the convergence rate of Markov chains on Rn." Stochastic Processes and their Applications 87.1 (2000): 1-23.

Weakness 2. Thank you for the feedback. The goal of our experiment was to illustrate the theoretical result, and we found that the empirical behavior aligned well with the theory. However, we agree that the example used is a toy setting, and we have not included comparison with baseline methods or error bars. In the revised version, we will follow your suggestion to improve the presentation of this toy 2D Gaussian mixture example and strengthen the empirical evaluation by adding a more realistic example.

Question 1. This has been addressed in the response to weakness 1.

We sincerely appreciate your careful review of our work and your insightful and constructive feedback. Below please find our point-by-point responses.

A bit more discussion on whether the analysis has methodological recommendations associated with it would be good, if something can be said.

We thank the reviewer for the thoughtful question. First, the main methodological contribution is that our theory provides practical guidance on how to tune algorithm parameters, such as the temperature ladder, as described in Appendix C.4. The conditions we have on the temperature ladder are not artifacts of our proof techniques (though a proof is not provided): intuitively, we need (1) the two adjacent temperatures to be close enough so that temperature swaps are likely to be accepted, and (2) the highest temperature to be large enough so that the tempered target can be explored efficiently by an algorithm like random walk Metropolis-Hastings or MALA. The conditions we have on $\beta_i$ make these intuitions precise. Second, our theory characterizes the dependence of the complexity of STMH on key parameters such as mode separation, component weights, and condition number. In practice, one can use it to approximately estimate how the runtime would change if the problem becomes more difficult.

The paper briefly mentions in Remark 5 that the new bound has an exponential dependence in the dimensionality. I feel this should be flagged in the introduction, as well as by adding a row in Table 1 to include the proposed bound.

We thank the reviewer for the suggestion. We agree that the exponential dependence on dimension should be made more transparent, and we will highlight this point in the revised version.

The paper falls in the usual fallacy of using upper bound comparisons to make statements about comparing algorithms, for example while 254 STMH achieves a more favorable logarithmic dependence on $\varepsilon - what is being compared there are upper bounds, so no statement should be made about STMH achieving a more favourable dependence.

We agree with the reviewer that the original phrasing may be misleading. Our intention was simply to highlight that the bound we obtain for STMH represents a state-of-the-art theoretical result among existing known upper bounds. We will revise the wording to make this distinction clear in the revision.

Question 1. We thank the reviewer for raising this truly insightful question. We did consider this very subtle issue when developing our proof, and the current version has already addressed this, as we now explain. We recognize that the writing needs to be improved at a few places to further clarify the related argument.

First, to clarify, Algorithm 1 is always run for a fixed number of steps $N$ and returns the sample one collects at the last step (we will clarify this in the revision.) If the resulting sample does not come from the desired temperature level, then Algorithm 1 is simply re-run for another $N$ steps. We do not run Algorithm 1 for a random number of steps until the desired temperature is reached (so no random stopping time is involved). Next, in Lemma 4, we analyze the total variation distance between the conditional distributions given the temperature index, which can be applied to the output of Algorithm 1 conditioning on the temperature level. Finally, to properly estimate the overall complexity of the whole algorithm including estimation of normalizing constants, in Lemma 27 we analyze how many times we need to re-run Algorithm 1 with $N$ steps (for fixed $N$) in order to obtain a sample from the desired temperature level. Note that we essentially followed the simulated tempering algorithm construction given in [1] but this final step (i.e., Lemma 27) apparently was missing in the proof of [1]. We will make this reasoning more explicit in the revised version.

Using a random stopping time (e.g. the first hitting time after N steps of a given temperature level) could be potentially problematic, since due to the Markovian nature of the algorithm, the distribution of such a sample collected can be severely biased. Of course, we are aware that our construction of Algorithm 1 is not practical; it is designed to facilitate the proof by producing independent samples for estimating normalizing constants.

[1] Ge, Rong, Holden Lee, and Andrej Risteski. ``Simulated tempering Langevin Monte Carlo ii: An improved proof using soft Markov chain decomposition." arXiv preprint arXiv:1812.00793 (2018).

First of all, we sincerely appreciate your careful reading of our work and your detailed feedback. Below please find our point-by-point responses.

Weakness 1. Thanks for the question. Simulated tempering (ST) is a widely used MCMC sampling technique (see, e.g., the textbook [1]). The use of STMH dates back at least to the seminal work of [2] on the Ising model, and it remains as one of the most popular methods for tackling metastability in the chemical/statistical physics literature [3]. In Bayesian statistics, the theory and methodology of tempering-based methods have also been studied extensively since the work of [4].

ST can be combined with any local MCMC kernel. In this work, we choose to focus on the random walk MH algorithm, which is a simple but illustrative and analytically tractable example. We analyze the complexity of STMH using a new decomposition theorem, which can also be used to analyze the convergence of ST combined with other local schemes such as MALA. Further, we expect that this technique can be adapted for analyzing other tempering-based algorithms such as parallel tempering (which has been utilized in many deep learning tasks [5, 6]), and those more advanced variants of ST such as non-reversible ST [7]. Although random walk MH could be a naive choice in some applied settings, the broader contribution of this work lies in the general theoretical framework we develop, which can inform the analysis and design of general tempering-based samplers.

[1] S. Brooks, A. Gelman, G. Jones, and X. Meng, eds. "Handbook of Markov Chain Monte Carlo.'' CRC press, 2011.

[2] E. Marinari, and P. Giorgio. "Simulated tempering: a new Monte Carlo scheme.'' Europhysics Letters 19.6 (1992): 451.

[3] Z. You, L. Li, J. Lu, and H. Ge. ``Integrated tempering enhanced sampling method as the infinite switching limit of simulated tempering." The Journal of Chemical Physics 149.8 (2018).

[4] C. J. Geyer, and E. A. Thompson. "Annealing Markov chain Monte Carlo with applications to ancestral inference." Journal of the American Statistical Association 90.431 (1995): 909-920.

[5] C. Smith, Q. T. Campbell, and T. Albash. ``Optimizing temperature distributions for training neural quantum states using parallel tempering." Physical Review E 111.5 (2025): 055306.

[6] W. Deng, Q. Feng, L. Gao, F. Liang, and G. Lin. "Non-convex learning via replica exchange stochastic gradient mcmc." In International Conference on Machine Learning, pp. 2474-2483. PMLR, 2020.

[7] M. Biron-Lattes, T. Campbell, and A. Bouchard-Côté. "Automatic regenerative simulation via non-reversible simulated tempering." Journal of the American Statistical Association 120.549 (2025): 318-330.

Weakness 2. We thank the reviewer for the suggestion. Roughly speaking, the spectral gap of a Markov chain quantifies the average-case convergence rate towards the stationary distribution; a larger gap implies faster mixing and thus more efficient sampling. However, in many scenarios, the spectral gap can be difficult to bound or vanish quickly, and the restricted spectral gap serves as a more robust and informative technique.

To explain simulated tempering more intuitively: the key idea is to augment the state space with a temperature index that allows the chain to move between a family of distributions, ranging from the original target distribution (low temperature) to a flattened version (high temperature). At higher temperatures, the probability landscape becomes smoother, allowing the chain to move more easily between modes. Transitions between temperature levels are governed by a Metropolis–Hastings rule, and samples from the target distribution are obtained by collecting states only at the lowest temperature (target distribution). This helps overcome the energy barriers between modes and promotes global exploration.

Weakness 3. Thanks for the suggestion. Assumptions 1 and 2 are motivated by the structure of mixture models, also used in prior works such as [8]. Intuitively, the sampler's convergence can be arbitrarily slow if (i) the mixture components are very far away from each other, or (ii) the mixing of the local chain targeting a single component is very slow (e.g. due to non-log-concavity). So some conditions are always needed to achieve a meaningful bound. We also want to clarify that these assumptions are only made for the decomposition theorem, and when applying the decomposition theorem to Gaussian mixtures, we rigorously prove that these assumptions hold. We agree that a more intuitive explanation would improve clarity, and in the revised version, we will include more discussion on the intuition behind the assumptions and interpretation of key constants and bounds.

[8] R. Ge, H. Lee, and A. Risteski. "Simulated tempering Langevin Monte Carlo ii: An improved proof using soft Markov chain decomposition." arXiv preprint arXiv:1812.00793 (2018).

Weakness 4. Thanks for the question. First, our result characterizes the complexity of the algorithm (e.g. polynomial dependence on D), which tells practitioners how the runtime of the algorithm scales as the problem becomes more difficult. Note that the dependence on other key parameters such as component weights and condition number is given in our proof of Theorem 5 in Appendix C.4.3. Second, our result also provides practical guidance for how the tuning parameters of ST should be selected such as the temperature ladder; see Appendix C.4.

Weakness 5. We thank the reviewer for the thoughtful question. In Figure 1, we empirically study the convergence behavior of STMH as a function of mode separation. The observed relationship is approximately linear rather than exact, likely because (i) we use the rate at which the empirical sample mean approaches the target mean as a proxy for convergence, rather than measuring the total variation distance directly, and (ii) our theoretical result provides only an upper bound on the convergence time, rather than an exact characterization.

Minor 1. We will correct it in the revision. Thanks for pointing it out.

Minor 2. Our result is similar in spirit to the decomposition theorem in Ge et al. (2018) but not directly comparable, as our results applies in a discrete-time setting under different assumptions. We will clarify this in the revision.

Minor 3. Thanks for the suggestion. We will update the section title accordingly.

Question 1. Thank you for the question. To clarify, $r_i$ is not needed to run the algorithm, and we define $r_i \propto Z_i / \hat{Z}_i$ exactly because $Z_i$ is typically unknown in practice. When implementing the STMH algorithm, one needs to evaluate the acceptance probability given in Equation (1), and this definition of $r_i$ ensures that the acceptance probability is independent of $Z_i$ and relies solely on the estimate $\hat{Z}_i$. In Appendix C.4, we show how sufficiently good estimates for $\hat{Z}_i$ can be determined, which makes the whole STMH algorithm fully implementable. We understand that the writing might have caused confusion and will clarify this in the revision.

Question 2. Thanks for the thoughtful question. Thm 5 establishes the convergence guarantees for STMH when sampling from a mixture of Gaussians with shared covariance matrix. This result can be naturally extended to target distributions that are sufficiently close to such Gaussian mixtures, as shown in Ge et al. (2018). In such cases, the time complexity would additionally depend on closeness between the actual distribution and the Gaussian mixture approximation. We also anticipate that similar techniques can be adapted to handle mixtures of log-concave distributions that are translates to each other, or more broadly, distributions that are well-approximated by such mixtures. However, as demonstrated by Ge et al. (2018), some seemingly mild violations of the assumptions, such as component covariance matrices differing by a constant factor, can lead to exponential time complexity for any reasonable algorithm with similar guarantees.

Question 3. Following up on our answer to Q2, there are realistic applied settings where our assumptions may approximately hold, making the proposed framework useful. For instance, [9] showed that the posterior distribution for Bayesian variable selection (which is one of the most common model selection problems in statistics) can be approximated by a mixture of Gaussians under certain conditions. Another relevant application arises in the context of (Bayesian) DNNs. Let the target be $e^{-f}$. In the context of DNN, $f$ is the loss composed with a neural network function, and the complexity of sampling will depend on the structure of $f$. If the widths of the neural networks tend to infinity, several works have shown that $f$ satisfied PL-type conditions and hence recent works like [10, 11, 12] could be used to show $e^{-f}$ satisfies Poincare or Log-Sobolev type inequalities. Therefore, algorithms like LMC, proximal sampler or MALA can be used to sample from the posterior efficiently. However, in the practical regime where the width is finite, the $e^{-f}$ does not satisfy any isoperimetric inequalities and is highly multimodal, where the proposed algorithm is likely to be useful.

[9] I. Castillo, J. Schmidt-Hieber and A. Van der Vaart. "Bayesian linear regression with sparse priors." The Annals of Statistics 43.5 (2015): 1986-2018.

[10] Y. Gong, N. He and Z. Shen. "Poincare inequality for local log-Polyak-Łojasiewicz measures: Non-asymptotic analysis in low-temperature Regime." arXiv:2501.00429 (2024).

[11] A. Y. Chen and K. Sridharan. "Optimization, isoperimetric inequalities, and sampling via Lyapunov potentials." arXiv:2410.02979 (2024).

[12] S. Chewi and A. J. Stromme. ``The ballistic limit of the log-Sobolev constant equals the Polyak-Łojasiewicz constant." arXiv:2411.11415 (2024).