Thank you for your positive and constructive feedback. We appreciate your recognition of the optimal bounds in specific regimes, the improved algorithm enabled by parallelism, and the implications regarding the separation between sampling and optimization, as well as the impossibility of general high-accuracy samplers. Below, we address your specific suggestions and questions:

"Notation is at times confusing...Clarify the notation."

Thank you for pointing this out. We acknowledge that the notation in some parts of the paper may be confusing, and we will revise the manuscript to ensure consistency and clarity. Specifically, we will:

• Use $V$ to represent the potential function consistently throughout the paper.
• Explicitly state that $L$ is a positive scalar.
• Use $\mathsf{Or}$ to denote the oracle and reserve $\mathcal{O}$ exclusively for complexity notation.

"The initialisation assumption on line 129 ... it might be nice to explicitly mention this to give further intuition about the assumption."

Thank you for pointing out that it would be helpful to make the intuition behind the initialization assumption more explicit. In fact, by Lemma 32 in [CEL+24], it suffices to initialize at a point x such that the optimality gap $V(x) - \min V = O(d)$. This condition ensures that the relative Fisher information is finite and meaningful for convergence analysis. Various methods from the optimization literature, such as gradient descent or gradient flow trapping [BM24, ZHTS25], can be used to efficiently find such a point. We will clarify this intuition in the revised manuscript.

[CEL+24] Analysis of Langevin Monte Carlo from Poincaré to Log-Sobolev, Sinho Chewi, Murat A. Erdogdu, Mufan Bill Li, Ruoqi Shen, Matthew Zhang, COLT 2022

[BM24] How to Trap a Gradient Flow, Sébastien Bubeck and Dan Mikulincer, COLT 2020

[ZHTS25] The adaptive complexity of finding a stationary point, Huanjian Zhou, Andi Han, Akiko Takeda, and Masashi Sugiyama, COLT 2025

"why the last line above 425 needs to be used....Why the bound on the exponential function needs to be applied."

Thank you for the suggestion. We agree with your observation and will revise it to keep the exponential form and fold it.

"The paper does not contain any experiments whatsoever. Is it possible to somehow plot empirical bounds for non log concave densities?"

Thank you for the suggestion. We have completed a standard experiment on a 2-dimensional Gaussian mixture, which demonstrates that the sequential method from [BCE+22] requires approximately 2000 iterations to converge, while our parallel method achieves comparable accuracy in just 300 iterations. We will include these results in the revised version and also plan to conduct additional experiments on more realistic and representative non-log-concave targets.

Regarding typos and confusing notations: We will carefully proofread the manuscript and correct all grammatical and typographical errors to improve clarity and overall presentation, and ensure that the notation is consistent throughout.

Thank you for your detailed and positive feedback. We greatly appreciate your recognition of our very nice results—with matching upper and lower bounds—as well as your appreciation of the insightful discussion and the paper’s clear structure, strong motivation, and well-organized presentation of the main results and their proofs. Below, we address your specific suggestions and questions:

"...issue on readability...For readers that are unfamiliar with parallel sampling...pose a challenge..."

Thank you for your thoughtful feedback and for highlighting the importance of improving the readability of the algorithm section. We agree that the parallel Picard method, as introduced in [ZS24], may not yet be familiar to all readers in the sampling community.

To improve clarity, we will revise Section 3.1 to include a more detailed explanation of the algorithm's intuition, motivation, and structure. Specifically, we will enhance the explanation of Algorithm 1 by providing more intuition behind its design, clarifying the purpose of each line in the pseudocode, and including an illustrative example of the update step to address potential confusion. Additionally we will add a section in Appendix to review existing parallel methods for sampling and diffusion models.

For your specific questions, we address as below:

• Q1: What's the purpose of $J$ and $N$? A1: The variable $j$ indexes the number of Picard updates applied to each time slice, and $J$ is the total number of such updates required to ensure convergence along the Picard direction. The parameter $N$ denotes the number of time slices used to discretize the time direction.

• Q2: Why we need to split the same update rule (lines 11 & 15 in Alg. 1) into two parts with different indices? A2: Our algorithm follows a diagonal-style update pattern. The first loop corresponds to updating the upper-left rectangles in Figure 1, while the second loop handles the remaining rectangles. While it is possible to merge both loops into a single unified loop, we chose to separate them for clarity and to better illustrate the structure of the computation graph.

• Q3: How does it improve upon [ACV24]? A more interesting question to think about is that, if one applies the averaged LMC into the algorithm in [ACV24], how would the complexity bound be compared with Thm. 3.1? A3: The key distinction between [ZS24] and [ACV24] is that [ACV24] updates time slices sequentially, whereas [ZS24] performs parallel updates across time slices. As a result, [ZS24] achieves a better complexity of $O(\log d)$ compared to the $O(\log^2 d)$ bound in [ACV24] for the strongly log-concave case. Our paper builds on the parallel framework of [ZS24], leading to an adaptive complexity of

N \cdot \log\left(\frac{L d}{\varepsilon^2}\right) = O\left(\frac{L K_0}{\varepsilon^2} \cdot \log\left(\frac{L d}{\varepsilon^2}\right)\right),

Thank you for your detailed and positive feedback. We appreciate your recognition of our improved upper and matching lower bounds for the important problem of parallel non-log-concave sampling. We are glad you highlighted the broader implications, such as the impossibility of a general high-accuracy sampler via parallelism, and the separation between optimization and sampling when using parallelism in strongly log-concave versus non-log-concave settings. Below, we address your specific suggestions and questions:

"...The algorithm is not clearly presented...Currently, Figure 1 + Algorithm 1 do not have enough explanations...arrows in Figure 1 are quite abstract and confusing."

Thank you for your thoughtful feedback. We acknowledge that the current presentation of Algorithm 1 and Figure 1 may lack clarity. To improve this, we will revise Figure 1 by replacing rectangles with lines and points to represent trajectories and grid points, respectively, and use color coding to distinguish computation graphs across different iterations. In Section 3.1, we will enhance the explanation of Algorithm 1 by providing more intuition behind its design, clarifying the purpose of each line in the pseudocode, and including an illustrative example of the update step to address potential confusion. Furthermore, we will add a discussion in Appendix that reviews existing parallel methods for sampling and diffusion models.

"..how practical... poly(d) queries can be made per iteration... especially in high dimensions."

Thank you for raising this important point. For the upper bound, the total number of oracle queries per iteration is $O(d^2)$, assuming $K_0 = O(d)$, which remains practical in moderate to high-dimensional settings. To enhance scalability in practice, we consider a sliding window strategy with window size $W$, which limits the number of points updated per iteration and thereby reduces memory usage while preserving parallelism.

Under this framework, when $W = 1$, the algorithm reduces to the sequential setting; when $W = \frac{d^2}{\varepsilon^4}$ (assuming $K_0 = O(d)$), it recovers our fully parallel method. The resulting space complexity is $Wd$, and the adaptive complexity becomes

where the space complexity is measured in terms of word usage rather than bit count.

As for the lower bound, assuming algorithms that make a polynomial number of queries in $d$ is standard in the literature (e.g., [ZWS24,ZHTS25]), and this aligns with prior work on adaptive complexity lower bounds in high-dimensional settings. We will incorporate the above discussion into the revised version of the manuscript.

"The current exposition of the lower bound is not very friendly... page 8 line 231-233, why is the probability at least $1/2$?...Lemma 4.5 is not proved/referenced... outline of the proof is not well-stated... Donsker-Varadhan variational principle..."

Thank you for the thoughtful feedback. We acknowledge that the exposition of the lower bound may not be sufficiently accessible to non-experts. In the revision, we will improve the structure of the lower bound section by adding a high-level overview at the beginning to guide the reader through the main proof strategy and clarify the logical flow.

We will expand the explanation of Theorem 4.4 and add the proper citation and context for Lemma 4.5 (which corresponds to Lemma 17 in [BCE+22]).

Regarding the probability lower bound of at least $1/2$ in Lines 231–233, this follows directly from Markov's inequality. This bound is used in the reduction from minimizing the relative Fisher information to the task of locating a stationary point.

Regarding the use of the Donsker-Varadhan variational principle, when $\mathbb{E}[\exp{g(X)}]<+\infty$ we can still apply the Donsker-Varadhan variational principle (Theorem 4.6 in [PW22]).

[PW22] Polyanskiy, Yuri, and Wu, Yihong. Information Theory: From Coding to Learning. Cambridge University Press.

"...averaging is not mentioned in the statement of theorem 3.1."

Thank you for your comment. We would like to clarify that the output in Line 17 of Algorithm 1 corresponds to the averaged distribution.

"Algorithm 1: does line 8 initializes all $x_{n,m}^0$ identically? Also, I have difficulty parsing the setup of $n$ and $j$ in line 10 and line 14 of the algorithm. Maybe some explanations will help readers."

Thank you for the question. In Line 8 of Algorithm 1, all $x_{n,m}^0$ are initialized identically for all $n$ and $m$. Regarding the indices, $n$ corresponds to the index of the time slice, while $j$ denotes the index of Picard iterations. We will clarify this in the revised version and improve the explanation of Algorithm 1 in Section 3.1 to aid readability.

"...In the proof of Lemma B.5, the equation after line 424, I am not sure how those two binomials pop out in the second inequality..."

Thank you for catching that. You're right and the first displayed inequality is incorrect. The second inequality is the correct one, and it follows directly from the decomposition established in Lemma B.4. We will remove the first inequality and explicitly refer to Lemma B.4 in the revised proof to make this clear.

Regarding typos: We will carefully proofread the manuscript and correct all grammatical and typographical errors to improve the clarity and overall presentation in the revised version.

Thank you for your clear response! Two minor points:

1. For the Donsker-Varadhan principle, it might be better to refer to the version you used in the rebuttal, which aligns with the setting in the paper.
2. I don't see why there is an averaging in line 17 in Algorithm 1. Is there a typo? I don't see dependence on $n$ and $m$ in the expression of $x_t$, but the for loop is wrt $n$ and $m$.

We sincerely thank the reviewer for their thoughtful evaluation and insightful questions. We are encouraged by your recognition of our novel application of the parallel Picard method to achieve efficient sampling. Below, we address your questions and comments:

"...focus on the specific regime of large FI...Is it possible to extend the results of Theorem 4.1 to smaller values $\varepsilon$?"

We focus on the large Fisher Information (FI) regime because the reduction from minimizing FI to minimizing the gradient norm, formalized in Lemma 4.3, is valid only when FI is sufficiently large. This mirrors the sequential setting, where a similar reduction underpins Theorem 9 in [CGLL22]. Our Theorem 4.1 strengthens that result by extending it to the parallel setting.

As for the possibility of extending the result to the small FI regime, as in Theorem 10 of [CGLL22], we believe this would require fundamentally different techniques. Specifically, if we adopt the bump function construction from [CGLL22] and place it in an unreachable region as in [ZWS25], one can show that with $d$ iterations and $\mathrm{poly}(d)$ queries per iteration, it is not possible to produce a sample with FI smaller than $\exp(-cd)$, for some constant $c > 0$. However, this is strictly weaker than Theorem 10 in [CGLL22], which establishes a lower bound of $1/\varepsilon^2$ for the total query complexity. In particular, Theorem 10 in [CGLL22] implies that achieving such accuracy would require exponentially many queries in $d$, and since the number of queries per iteration is polynomial in $d$, the total number of iterations must also be exponential in $d$. Such result is much stronger hardness guarantee than the $d$ iteration regime we can obtain.

[CGLL22] Fisher information lower bounds for sampling, Sinho Chewi, Patrik Gerber, Holden Lee, Chen Lu, ALT 2023,

[ZWS25] The adaptive complexity of parallelized log-concave sampling, Huanjian Zhou, Baoxiang Wang, Masashi Sugiyama, ICLR 2025

"...a discussion of the space complexity of the parallel algorithm is necessary..."

Thank you for highlighting the importance of the space complexity. In practice, fully parallel updates may require excessive memory. To address this, we consider a slicing window strategy, as proposed in [SBE+24], where in each iteration we update at most $W > 0$ points in parallel. When $W = 1$, this recovers the sequential method; when $W = \frac{d^2}{\varepsilon^4}$ (assuming $K_0 = O(d)$), it recovers our fully parallel method.

Under this framework, we obtain the following tradeoff: the space complexity is $Wd$, and the adaptive complexity becomes$O\left(\frac{d}{\varepsilon^2} + \frac{d^2}{\varepsilon^4 W} \log \frac{d}{\varepsilon^2}\right),$ where we define the space complexity in terms of word usage rather than bit count.

As noted by Reviewer 2U73, an interesting intermediate case arises when applying averaged Langevin Monte Carlo to the algorithm in [ACV24], which updates time slices sequentially. This corresponds to $W = \frac{d}{\varepsilon^2}$, yielding an adaptive complexity of $O\left(\frac{d}{\varepsilon^2} \log \frac{d}{\varepsilon^2} \right),$ which fits naturally within the above tradeoff curve. We will incorporate the above discussion into the revised version of the manuscript.

[SBE+24] Parallel Sampling of Diffusion Models, Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, Nima Anari, NeurIPS 2024

"The absence of experimental results limits the validation and practical support of the theoretical findings."

Thank you for the suggestion. We have completed a standard experiment on a 2-dimensional Gaussian mixture, which demonstrates that the sequential method from [BCE+22] requires approximately 2000 iterations to converge, while our parallel method achieves comparable accuracy in just 300 iterations. We will include these results in the revised version and also plan to conduct additional experiments on more realistic and representative non-log-concave targets.

"Theorem 4.1...does not explicitly show the dependence on $\varepsilon$, Would it be clearer to retain th dependence in the expression on line 248..."

Thank you for the suggestion. We agree it would improve clarity to include the $\varepsilon$-dependence, and we will revise the text accordingly.