In-and-Out: Algorithmic Diffusion for Sampling Convex Bodies

$**What practical advantage of the method?**$

The main advantage of our method is the simplicity of the algorithm and the analysis with provable guarantees in commonly used probabilistic metrics, more general than previously known. Not only is our algorithm much simpler to analyze, but it also provides some tangible improvements (such as improving the guarantee from $\mathsf{TV}$ to $\mathsf{KL}$ and $\mathcal{R}_q$); it also shows a clear and direct connection between isoperimetry and convergence. In practice, we might expect comparable performance between all these approaches, although this would depend on the details of the implementation.

$**Does the result follow from the proximal sampler by taking some suitable limit**$

As demonstrated in Lemma 12 and Theorem 3, the mixing guarantee of our sampler (i.e., how many outer loops are needed) does follow from a limiting argument. However, what matters in the end is a bound on query complexity (not just mixing rate), which essentially requires us to bound the number of rejections throughout backward steps.

The approach required for the constrained case is not comparable to the unconstrained case. Previous work on proximal sampling works under a well-conditioned setting without hard constraints, where the number of trials for rejection sampling (for the backward step) is always $\mathcal{O}(1)$ regardless of the forward step. However, in the presence of constraints, this type of analysis for the backward step is no longer possible. Analyzing proximal-type sampling for this setting has been a well-known open problem. Instead, we carried out the analysis for the backward step in a more careful way. As a result, there are a number of new ideas in the rejection analysis.

$**Line 146-149: How would INO compare to projection-based methods**$

In general, a projection oracle is stronger than a membership oracle, and its implementation using a membership oracle requires $O(d^2)$ membership calls per projection in the worst case. One can ask instead what happens if we assume a projection oracle which is roughly the same cost as a membership oracle. In this case, leveraging projection might be faster, but this is a non-trivial problem. We leave this for future work.

$**Failure = Lazy? How should one take ergodic averages along the chain for computing expectation of some observable?**$

We thank the reviewer for this question. In our view, failure of our chain is not the same as taking lazy steps in e.g. Ball walk. Lazy chains have a 1/2 chance to remain stationary at each instant, potentially doubling the iterations needed for convergence. By contrast, our algorithm has an arbitrarily small failure probability (moreover, the mixing guarantee has only a poly-log dependence $\text{polylog}(1/\eta)$). In the event of failure, we restart the algorithm until we succeed, which does not introduce any bias.

Regarding the second question, it is still possible to take ergodic averages along the chain (assuming you pick the failure probability at each iteration to be sufficiently small so that your run rarely fails during the entire horizon of your run). This is because our Markov chain still has a spectral gap, and so the same techniques used to obtain guarantees for ergodic averages continue to work in our case.

$**How to obtain a warm start in practice**$

Obtaining warm-starts in practice is non-trivial, and requires some type of annealing algorithm in general. Indeed, warm-start generation has been studied for more than two decades in theoretical computer science. See e.g., [1, 2] for how this can be done. For simplicity, we assume that a warm start is provided in our algorithm, although one can use the guarantees in [1, 2] to provide a rigorous justification.

[1] Gaussian Cooling and $O^*(n^3)$ Algorithms for Volume and Gaussian Volume, Ben Cousins and Santosh S. Vempala.

[2] Reducing Isotropy and Volume to KLS: An $O(n^3\psi^2)$ Volume Algorithm, Jia et al.

We are unclear on the reviewer's "reservation". The problem of sampling constrained convex bodies has been widely studied for decades and so far did not have guarantees beyond TV distance; moreover, in recent years, it has been a well-known open problem to see if diffusion-based methods could be used to obtain a polytime algorithm. Our main results show stronger guarantees for the classical problem using a simple and clean diffusion approach. The resulting analysis brings together several well-known analysis components, along with some extensions (to the constrained setting) and some new ideas (rejection analysis). We feel that the fact the solution is relatively simple and does not need substantial technical sophistication is an attractive feature.

$**Initialization process and parameter selection**$

• Initialization process: We assume that the initial start is $M$-warm in this work. Obtaining this warm-start is non-trivial, and requires some type of annealing algorithm in general. Warm-start generation has been studied for more than two decades in theoretical computer science. See e.g. [1] for how this can be done.
• Parameter selection: The details of $N, h, \eta$ are given in the proof (see Lemma 3 and 14 for example). We will make this more clear in the initial Theorem statement as well. Note that we do not track the absolute constants hidden in big-O notation, but this can be obtained from a careful reading of the proof; our analysis was not optimal with respect to numerical constants, which can likely be sharpened with a more precise argument.

$**How FI parameters affect the convergence rate for specific geometry**$

The Poincare constant in general is on the same order as the maximum eigenvalue of the covariance, while the log-Sobolev constant is bounded by the squared diameter of the convex body (please refer to Appendix D). However, for a more structured body, it is not clear what the sharpest possible constant is for these functional inequalities. In general, we expect that it would be very difficult to estimate unless the body is something simple like an $\ell_p$ ball. We can also obtain the constants under a linear map $T: x \mapsto Ax$ by multiplying the constants by $||A||^2$.

$**Extension to a general setting**$

We thank the reviewer for their insightful suggestion. While our framework could potentially accommodate it, the rates after incorporating a first-order oracle for sampling $e^{-f}1_K$ are not immediately clear from our analysis and would take some more effort. We leave this for future work.

[1] Gaussian Cooling and $O^*(n^3)$ Algorithms for Volume and Gaussian Volume, Ben Cousins and Santosh S. Vempala.

$**Techniques are already existing**$

As noted by Reviewer W5VJ, our paper is not “just” putting together known components. Previous work on proximal sampling works under a well-conditioned setting without hard constraints, where the number of trials for rejection sampling (for the backward step) is always $\mathcal{O}(1)$ regardless of the forward step. However, in the presence of constraints, this type of analysis for the backward step is no longer possible. Analyzing proximal-type sampling for this setting has been a well-known open problem. Instead, we carried out the analysis for the backward step in a more careful way. As a result, there are a number of new ideas in the rejection analysis.

$**Extension to a general setting**$

We thank the reviewer for their insightful suggestion. While our framework could potentially accommodate it, the rates after incorporating a first-order oracle for sampling $e^{-f}1_K$ are not immediately clear from our analysis and would take some more effort. We leave this for future work.

$**Polytope constraint?**$

There are several ways of leveraging the explicit structure of the log-concave sampling problem. For instance, one can easily implement a projection oracle for polytope constraints (which is stronger than a membership one), so one may consider a projection-based algorithm. Another approach for exploiting the structure is to use barrier-related information. In general, a convex constraint admits a self-concordant barrier $\phi$, and samplers with the local metric given by the Hessian $\nabla^2 \phi$ are provably fast and practical (due to condition-number independence, for instance).

$**Application**$

Our work mainly aims to provide a new algorithmic/analytic framework for the uniform sampling problem under the membership oracle model. The applications of this problem are already widespread and well-known, therefore we do not feel the need to propose any new applications. Instead, we note that this problem is used in a number of fundamental settings: volume computation for convex bodies; as a core subroutine in general convex body sampling; metabolic flux sampling in systems biology; Bayesian inference in statistics; and other domains of scientific computing. Thus, new results for the core algorithmic problem would immediately imply theoretical improvements for all these settings.

$**Extension to potentials**$

We thank the reviewer for their insightful suggestion. While our framework could potentially accommodate it, the rates after incorporating a first-order oracle for sampling $e^{-f}1_K$ are not immediately clear from our analysis and would take some more effort. We leave this for future work.

$**Confusion around Line 175**$

We apologize for the confusion. Corollary 1 holds for the near-isotropic case as well, since the operator norm of the covariance is $\mathcal{O}(1)$.

$**How to make it isotropic**$

The question of obtaining an isotropic position has been studied for more than two decades in theoretical computer science. For references, please consult [1, 2] for an idea of how this can be done. In summary, their approach is to draw a few samples from the body and then obtain a rough estimate of the covariance $T$. Then, applying a proper linear map $T^{-1/2}$ to the body will reduce the skewness. Repeating these steps $\mathcal{O}(\log d)$ many times (along with some type of annealing algorithm) ensures that a transformed body is nearly isotropic with high probability. Although we had mentioned this procedure in our related works, we will make a comment earlier in our paper for clarity. We thank the reviewer for their question.

$**Regarding Line 139-141 (time-reversal of SDE)**$

The final paragraph of Appendix E discusses why the pointwise reversal property is needed. In particular, it allows us to claim that given a starting point $z$, the reverse SDE from such a starting point has distribution $\pi^{X|Y=y}$, a property used in the initial paragraphs of Appendix B. This will be clarified in our revision.

[1] Gaussian Cooling and $O^*(n^3)$ Algorithms for Volume and Gaussian Volume, Ben Cousins and Santosh S. Vempala.

[2] Reducing Isotropy and Volume to KLS: An $O(n^3\psi^2)$ Volume Algorithm, Jia et al.