Faster Diffusion Sampling with Randomized Midpoints: Sequential and Parallel

## Abstract

Abstract Sampling algorithms play an important role in controlling the quality and runtime of diffusion model inference. In recent years, a number of works [ CCL + 23b , CCL + 23a , BDD23 , LLT22 ] have proposed schemes for diffusion sampling with provable guarantees; these works show that for essentially any data distribution, one can approximately sample in polynomial time given a sufficiently accurate estimate of its score functions at different noise levels.
In this work, we propose a new scheme inspired by Shen and Lee’s randomized midpoint method for log-concave sampling [ SL19 ] . We prove that this approach achieves the best known dimension dependence for sampling from arbitrary smooth distributions in total variation distance ( O ~ ⁢ ( d 5 / 12 ) ~ 𝑂 superscript 𝑑 5 12 \widetilde{O}(d^{5/12}) over~ start_ARG italic_O end_ARG ( italic_d start_POSTSUPERSCRIPT 5 / 12 end_POSTSUPERSCRIPT ) compared to O ~ ⁢ ( d ) ~ 𝑂 𝑑 \widetilde{O}(\sqrt{d}) over~ start_ARG italic_O end_ARG ( square-root start_ARG italic_d end_ARG ) from prior work). We also show that our algorithm can be parallelized to run in only O ~ ⁢ ( log 2 ⁡ d ) ~ 𝑂 superscript 2 𝑑 \widetilde{O}(\log^{2}d) over~ start_ARG italic_O end_ARG ( roman_log start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d ) parallel rounds, constituting the first provable guarantees for parallel sampling with diffusion models. As a byproduct of our methods, for the well-studied problem of log-concave sampling in total variation distance, we give an algorithm and simple analysis achieving dimension dependence O ~ ⁢ ( d 5 / 12 ) ~ 𝑂 superscript 𝑑 5 12 \widetilde{O}(d^{5/12}) over~ start_ARG italic_O end_ARG ( italic_d start_POSTSUPERSCRIPT 5 / 12 end_POSTSUPERSCRIPT ) compared to O ~ ⁢ ( d ) ~ 𝑂 𝑑 \widetilde{O}(\sqrt{d}) over~ start_ARG italic_O end_ARG ( square-root start_ARG italic_d end_ARG ) from prior work.

## 1 Introduction

Diffusion models [, , , , , , ] have emerged as the *de facto* approach to generative modeling across a range of data modalities like images [, ], audio [], video [], and molecules []. In recent years a slew of theoretical works have established surprisingly general convergence guarantees for this method [, , , , , , , , ]. They show that for essentially any data distribution, assuming one has a sufficiently accurate estimate for its *score function*, one can approximately sample from it in polynomial time.

While these results offer some theoretical justification for the empirical successes of diffusion models, the upper bounds they furnish for the number of iterations needed to generate a single sample are quite loose relative to what is done in practice. The best known provable bounds scale as $O(\sqrt{d}/\varepsilon)$, where $d$ is the dimension of the space in which the diffusion is taking place (e.g. $d=16384$ for Stable Diffusion) [], and $\varepsilon$ is the target error. Even ignoring the dependence on $\varepsilon$ and the hidden constant factor, this is at least $2-3\times$ larger than the default value of $50$ inference steps in Stable Diffusion.

In this work we consider a new approach for driving down the amount of compute that is provably needed to sample with diffusion models. Our approach is rooted in the *randomized midpoint method*, originally introduced by Shen and Lee [] in the context of Langevin Monte Carlo for log-concave sampling. At a high level, this is a method for numerically solving differential equations where within every discrete window of time, one forms an unbiased estimate for the drift by evaluating it at a random “midpoint” (see Section 2.2 for a formal treatment). For sampling from log-concave densities, the number of iterations needed by their method scales with $d^{1/3}$, and this remains the best known bound in the “low-accuracy” regime.

While this method is well-studied in the log-concave setting [, , , ], its applicability to diffusion models has been unexplored both theoretically and empirically. Our first result uses the randomized midpoint method to obtain an improvement over the prior best known bound of $O(\sqrt{d}/\varepsilon)$ for sampling arbitrary smooth distributions with diffusion models:

###### Theorem 1.1 (Informal, see Theorem A.10 ) .

Suppose that the data distribution $q$ has bounded second moment, its score functions $\nabla\ln q_{t}$ along the forward process are $L$-Lipschitz, and we are given score estimates which are $L$-Lipschitz and $\widetilde{O}(\frac{\varepsilon}{d^{1/12}\sqrt{L}})$(^1^11$\widetilde{O}(\cdot)$ hides polylogarithmic factors in $d,L,\varepsilon$ and $\operatorname*{\mathbb{E}}_{x\sim q}[\|x\|^{2}]$)-close to $\nabla\ln q_{t}$ for all $t$. Then there is a diffusion-based sampler using these score estimates (see Algorithm 1) which outputs a sample whose law is $\varepsilon$-close in total variation distance to $q$ using $\widetilde{O}(L^{5/3}d^{5/12}/\varepsilon)$ iterations.

Our algorithm is based on the ODE-based predictor-corrector algorithm introduced in [], but in place of the standard exponential integrator discretization in the predictor step, we employ randomized midpoint discretization. We note that in the domain of log-concave sampling, the result of Shen and Lee only achieves recovery in *Wasserstein distance*. Prior to our work, it was actually open whether one can achieve the same dimension dependence in total variation or KL divergence, for which the best known bound was $\widetilde{O}(\sqrt{d})$ [, , ]. In contrast, our result circumvents this barrier by carefully trading off time spent in the corrector phase of the algorithm for time spent in the predictor phase. We defer the details of this, as well as other important technical hurdles, to Section 3.1.

Next, we turn to a different computational model: instead of quantifying the cost of an algorithm in terms of the total number of iterations, we consider the *parallel* setting where one has access to multiple processors and wishes to minimize the total number of parallel rounds needed to generate a single sample. This perspective has been explored in a recent empirical work [], but to our knowledge, no provable guarantees were known for parallel sampling with diffusion models (see Section 1.1 for discussion of concurrent and independent work). Our second result provides the first such guarantee:

###### Theorem 1.2 (Informal, see Theorem B.13 ) .

Under the same assumptions on $q$ as in Theorem 1.1, and assuming that we are given score estimates which are $\widetilde{O}(\frac{\varepsilon}{\sqrt{L}})$-close to $\nabla\ln q_{t}$ for all $t$, there is a diffusion-based sampler using these score estimates (see Algorithm 9) which outputs a sample whose law is $\varepsilon$-close in total variation distance to $q$ using $\widetilde{O}(L\cdot\mathrm{polylog}(Ld/\varepsilon))$ parallel rounds.

This result follows in the wake of several recent theoretical works on parallel sampling of log-concave densities using Langevin Monte Carlo [, , ]. A common thread among these works is the observation that differential equations can be numerically solved via fixed point iteration (see Section 2.3 for details), and we adopt a similar perspective in the context of diffusions. To our knowledge this is the first provable guarantee for parallel sampling beyond the log-concave setting.

Finally, we show that, as a byproduct of our methods, we can actually obtain a similar dimension dependence of $\widetilde{O}(d^{5/12})$ as in Theorem 1.1 for log-concave sampling in TV, superseding the previously best known bound of $\widetilde{O}(\sqrt{d})$ mentioned above.

###### Theorem 1.3 (Informal, see Theorem C.2 ) .

Suppose distribution $q$ is $m$-strongly-log-concave, and its score function $\nabla\ln q$ is $L$-Lipschitz. Then, there is a underdamped-Langevin-based sampler that uses this score (Algorithm 11) and outputs a sample whose law is $\varepsilon$-close in total variation to $q$ using $\widetilde{O}\left(d^{5/12}\left(\frac{L^{4/3}}{\varepsilon^{2/3}m^{4/3}}+
\frac{1}{\varepsilon}\right)\right)$ iterations.

###### Theorem 1.1 (Informal, see Theorem A.10 ) .

###### Theorem 1.2 (Informal, see Theorem B.13 ) .

###### Theorem 1.3 (Informal, see Theorem C.2 ) .

### 1.1 Related work

Our discretization scheme is based on the randomized midpoint method of [], which has been studied at length in the domain of log-concave sampling [, , ].

The proof of our parallel sampling result builds on the ideas of [, , ] on parallelizing the collocation method. These prior results were focused on Langevin Monte Carlo, rather than diffusion-based sampling. We review these ideas in Section 2.3.

In [], the authors proposed the predictor-corrector framework that we also use for analysing convergence guarantee of the probability flow ODE and which achieved iteration complexity scaling with $\widetilde{O}(\sqrt{d})$. In addition to this, there have been many works in recent years giving general convergence guarantees for diffusion models [, , , , , , , , , , , , , , , ]. Of these, one line of work [, , , ] analyzed DDPM, the stochastic analogue of the probability flow ODE, and showed $\tilde{O}(d)$ iteration complexity bounds. Another set of works [, , , ] studied the probability flow ODE, for which our work provides a new discretization scheme for the probability flow ODE, that achieves a state-of-the-art $\widetilde{O}(d^{5/12})$ dimension dependence for sampling from a diffusion model.

#### Concurrent work.

Here we discuss the independent works of [] and []. [] gave an analysis for parallel sampling with diffusion models that also achieves a $\mathrm{polylog}(d)$ number of parallel rounds like in the present work. [] showed an improved dimension dependence of $\widetilde{O}(d^{5/12})$ for log-concave sampling in total variation, similar to our analogous result, but via a different proof technique. In addition to this, they show a similar result when the distribution only satisfies a log-Sobolev inequality. They also show empirical results for diffusion models, showing that an algorithm inspired by the randomized midpoint method outperforms ODE based methods with similar compute. While their work builds on the randomized midpoint method, they do not theoretically analyze the diffusion setting and do not study parallel sampling.

## 2 Preliminaries

### 2.1 Probability flow ODE

In this section we review basics about deterministic diffusion-based samplers; we refer the reader to [] for a more thorough exposition.

Let $q^{*}$ denote the data distribution over $\mathbb{R}^{d}$. We consider the standard Ornstein-Uhlenbeck (OU) *forward process*, i.e. the “VP SDE,” given by

$$ $\mathrm{d}x^{\rightarrow}_{t}=-x^{\rightarrow}_{t}\,\mathrm{d}t+\sqrt{2}\, \mathrm{d}B_{t}\,\qquad x_{0}^{\rightarrow}\sim q^{*}\,,$ (1) $$

where $(B_{t})_{t\geq 0}$ denotes a standard Brownian motion in $\mathbb{R}^{d}$. This process converges exponentially quickly to its stationary distribution, the Gaussian distribution $\mathcal{N}(0,\mathrm{Id})$.

Suppose the OU process is run until terminal time $T>0$, and for any $t\in[0,T]$, let $q^{*}_{t}\triangleq\text{law}(x^{\rightarrow}_{t})$, i.e. the law of the forward process at time $t$. We will consider the *reverse process* given by the *probability flow ODE*

$$ $\mathrm{d}x_{t}=(x_{t}+\nabla\ln q_{T-t}(x_{t}))\,\mathrm{d}t\,.$ (2) $$

This is a time-reversal of the forward process, so that if $x_{0}\sim q_{T}$, then $\text{law}(x_{t})=q^{*}_{T-t}$. In practice, one initializes at $x_{0}\sim\mathcal{N}(0,\mathrm{Id})$, and instead of using the exact *score function* $\nabla\ln q_{T-t}$, one uses estimates $\widehat{s}_{T-t}\approx\nabla\ln q_{T-t}$ which are learned from data. Additionally, the ODE is solved numerically using any of a number of discretization schemes. The theoretical literature on diffusion models has focused primarily on *exponential integration*, which we review next before turning to the discretization scheme, the *randomized midpoint method* used in the present work.

### 2.2 Discretization schemes

Suppose we wish to discretize the following semilinear ODE:

$$ $\mathrm{d}x_{t}=(x_{t}+f_{t}(x_{t}))\,\mathrm{d}t\,.$ (3) $$

For our application we will eventually take $f_{t}\triangleq\widehat{s}_{T-t}$, but we use $f_{t}$ in this section to condense notation.

Suppose we want to discretize Equation 3 over a time window $[t_{0},t_{0}+h]$. The starting point is the integral formulation for this ODE:

$$ $x_{t_{0}+h}=e^{h}x_{t_{0}}+\int^{t_{0}+h}_{t_{0}}e^{t_{0}+h-t}f_{t}(x_{t})\, \mathrm{d}t\,.$ (4) $$

Under the standard *exponential integrator* discretization, one would approximate the integrand by $e^{t_{0}+h-t}f_{t_{0}}(x_{t_{0}})$ and obtain the approximation

$$ $x_{t_{0}+h}\approx e^{h}x_{t_{0}}+(e^{h}-1)f_{t_{0}}(x_{t_{0}})\,.$ (5) $$

The drawback of this discretization is that it uses an inherently *biased* estimate for the integral in Eq. (4). The key insight of [] was to replace this with the following unbiased estimate

$$ $\int^{t_{0}+h}_{t_{0}}e^{t_{0}+h-t}f_{t}(x_{t})\,\mathrm{d}t\approx he^{(1- \alpha)h}f_{t_{0}+\alpha h}(x_{t_{0}+\alpha h})\,,$ (6) $$

where $\alpha$ is a uniformly random sample from $[0,1]$. While this alone does not suffice as the estimate depends on $x_{t_{0}+\alpha h}$, naturally we could iterate the above procedure again to obtain an approximation to $x_{t_{0}+\alpha h}$. It turns out though that even if we simply approximate $x_{t_{0}+\alpha h}$ using exponential integrator discretization,
we can obtain nontrivial improvements in discretization error (e.g. our Theorem 1.1). In this case, the above sequence of approximations takes the following form:

$$ $\displaystyle x_{t_{0}+\alpha h}$ $\displaystyle\approx e^{\alpha h}x_{t_{0}}+(e^{\alpha h}-1)f_{t_{0}}(x_{t_{0}})$ (7) $\displaystyle x_{t_{0}+h}$ $\displaystyle\approx e^{h}x_{t_{0}}+he^{(1-\alpha)h}f_{t_{0}+\alpha h}(x_{t_{0 }+\alpha h})\,.$ (8) $$

Note that a similar idea can be used to discretize stochastic differential equations, but in this work we only use it to discretize the probability flow ODE.

#### Predictor-Corrector.

For important technical reasons, in our analysis we actually consider a slightly different algorithm than simply running the probability flow ODE with approximate score, Gaussian initialization, and randomized midpoint discretization. Specifically, we interleave the ODE with *corrector* steps that periodically inject noise into the sampling trajectory. We refer to the phases in which we are running the probability flow ODE as *predictor* steps.

The corrector step will be given by running *underdamped Langevin dynamics*. As our analysis of this will borrow black-box from bounds proven in [], we refer to Section A.2 for details.

### 2.3 Parallel sampling

The scheme outlined in the previous section is a simple special case of the *collocation method*. In the context of the semilinar ODE from Eq. (3), the idea behind the collocation method is to solve the integral formulation of the ODE in Eq. (4) via fixed point iteration. For our parallel sampling guarantees, instead of choosing a single randomized midpoint $\alpha$, we break up the window $[t_{0},t_{0}+h]$ into $R$ sub-windows, select randomized midpoints $\alpha_{1},\ldots,\alpha_{R}$ for these sub-windows, and approximate the trajectory of the ODE at any time $t_{0}+i\delta$, where $\delta\triangleq h/R$, by

$$ $x_{t_{0}+\alpha_{i}h}\approx e^{\alpha_{i}h}x_{t_{0}}+\sum^{i}_{j=1}\Bigl{(}e^ {\alpha_{i}h-(j-1)\delta}-\max(e^{\alpha_{i}h-j\delta},1)\Bigr{)}\cdot f_{t_{0 }+\alpha_{j}h}(x_{t_{0}+\alpha_{j}h})\,.$ (9) $$

One can show that as $R\to\infty$, this approximation tends to an equality. For sufficiently large $R$, Eq. (9) naturally suggests a fixed point iteration that can be used to approximate each $x_{t_{0}+\alpha_{i}h}$, i.e. we can maintain a sequence of estimates $\widehat{x}^{(k)}_{t_{0}+\alpha_{i}h}$ defined by the iteration

$$ $\widehat{x}^{(k)}_{t_{0}+\alpha_{i}h}\leftarrow e^{\alpha_{i}h}\widehat{x}^{(k -1)}_{t_{0}}+\sum^{i}_{j=1}\Bigl{(}e^{\alpha_{i}h-(j-1)\delta}-\max(e^{\alpha_ {i}h-j\delta},1)\Bigr{)}\cdot f_{t_{0}+\alpha_{j}h}(\widehat{x}^{(k-1)}_{t_{0} +\alpha_{j}h})\,,$ (10) $$

for $k$ ranging from $1$ up to some sufficiently large $K$. Finally, analogously to Eq. (8), we can estimate $x_{t_{0}+h}$ via

$$ $x_{t_{0}+h}\approx e^{h}\widehat{x}^{(K)}_{t_{0}}+\delta\sum^{R}_{i=1}e^{(1- \alpha_{i})h}f_{t_{0}+\alpha_{i}h}(\widehat{x}^{(K)}_{t_{0}+\alpha_{i}h})\,.$ (11) $$

The key observation, made in [] and also in related works of [, , ], is that for any fixed round $k$, all of the iterations Eq. (10) for different choices of $i=1,\ldots,R$ can be computed in parallel. With $R$ parallel processors, one can thus compute the estimate for $x_{t_{0}+h}$ in $K$ parallel rounds, with $O(KR)$ total work.

### 2.4 Assumptions

Throughout the paper, for our diffusion results, we will make the following standard assumptions on the data distribution and score estimates.

###### Assumption 2.1 (Bounded Second Moment) .

$$ $\displaystyle\mathfrak{m}_{2}^{2}:=\operatorname*{\mathbb{E}}_{x\sim q_{0}} \left[\|x\|^{2}\right]<\infty.$ $$

###### Assumption 2.2 (Lipschitz Score) .

For all $t$, the score $\nabla\ln q_{t}$ is $L$-Lipschitz.

###### Assumption 2.3 (Lipschitz Score estimates) .

For all $t$ for which we need to estimate the score function in our algorithms, the score estimate $\widehat{s}_{t}$ is $L$-lipschitz.

###### Assumption 2.4 (Score Estimation Error) .

For all $t$ for which we need to estimate the score function in our algorithms,

$$ $\displaystyle\operatorname*{\mathbb{E}}_{x_{t}\sim q_{t}}\left[\|\widehat{s}_{ t}(x_{t})-\nabla\ln q_{t}(x_{t})\|^{2}\right]\leq\varepsilon_{\mathrm{sc}}^{2}.$ $$

###### Assumption 2.1 (Bounded Second Moment) .

###### Assumption 2.2 (Lipschitz Score) .

###### Assumption 2.3 (Lipschitz Score estimates) .

###### Assumption 2.4 (Score Estimation Error) .

## 3 Technical overview

Here we provide an overview of our sequential and parallel algorithms, along with the analysis of our iteration complexity bounds. We begin with a description of the sequential algorithm.

### 3.1 Sequential algorithm

Following the framework of [], our algorithm consists of “predictor” steps interspersed with “corrector” steps, with the time spent on each carefully tuned to obtain our final $\widetilde{O}(d^{5/12})$ dimension dependence. We first describe our predictor step – this is the piece of our algorithm that makes use of the Shen and Lee’s randomized midpoint method [].

Figure: Algorithm 1 PredictorStep (Sequential)

The main difference between the above and the predictor step of [] are steps $1(b)$ – $1(d)$. $1(b)$ and $1(c)$ together compute a randomized midpoint, and $1(d)$ uses this midpoint to obtain an approximate solution to the integral of the ODE. We describe these steps in more detail in Section 3.3.

Next, we describe the “corrector” step, introduced in []. First, recall the underdamped Langevin ODE:

$$ $\displaystyle\begin{split}\mathrm{d}\widehat{x}_{t}&=\widehat{v}_{t}\,\mathrm{ d}t\\ \mathrm{d}\widehat{v}_{t}&=(\widehat{s}(\widehat{x}_{\lfloor t/h\rfloor h})- \gamma\widehat{v}_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}B_{t}\end{split}$ (12) $$

Here $\widehat{s}$ is our $L^{2}$ accurate score estimate for a fixed time (say $t$). Then, the corrector step is described below.

Figure: Algorithm 2 CorrectorStep (Sequential)

Finally, Algorithm 3 below puts the predictor and corrector steps together to give our final sequential algorithm.

Figure: Algorithm 3 SequentialAlgorithm

For the final setting of parameters in Algorithm 3, see Theorem A.10. Now, we describe the analysis of the above algorithm in detail.

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Starting time $t_{0}$, Number of steps $N$, Step sizes $h_{n\in[0,\dots,N-1]}$, Score estimates $\widehat{s}_{t}$

- 1.
For $n=0,\dots,N-1$:
(a)
Let $t_{n}=t_{0}-\sum_{i=0}^{n-1}h_{i}$
(b)
Randomly sample $\alpha$ uniformly from $[0,1]$.
(c)
Let $\widehat{x}_{n+\frac{1}{2}}=e^{\alpha h_{n}}\widehat{x}_{n}+\left(e^{\alpha h_
{n}}-1\right)\widehat{s}_{t_{n}}(\widehat{x}_{n})ds$
(d)
Let $\widehat{x}_{n+1}=e^{h_{n}}\widehat{x}_{n}+h_{n}\cdot e^{(1-\alpha){h_{n}}}
\widehat{s}_{t_{n}-\alpha h_{n}}(\widehat{x}_{n+\frac{1}{2}})$
- 2.
Let $t_{N}=t_{0}-\sum_{i=0}^{N-1}h_{i}$
- 3.
Return $\widehat{x}_{N},t_{N}$.

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Total time $T_{\mathrm{corr}}$, Step size $h_{\mathrm{corr}}$, Score estimate $\widehat{s}$

- 1.
Run underdamped Langevin Monte Carlo in (12) for total time $T_{\mathrm{corr}}$ using step size $h_{\mathrm{corr}}$, and let the result be $\widehat{x}_{N}$.
- 2.
Return $\widehat{x}_{N}$.

#### Input parameters:

- •
Start time $T$, End time $\delta$, Corrector steps time $T_{\mathrm{corr}}\lesssim 1/\sqrt{L}$, Number of predictor-corrector steps $N_{0}$, Predictor step size $h_{\mathrm{pred}}$, Corrector step size $h_{\mathrm{corr}}$, Score estimates $\widehat{s}_{t}$

- 1.
Draw $\widehat{x}_{0}\sim\mathcal{N}(0,I_{d})$.
- 2.
For $n=0,\dots,N_{0}-1$:
(a)
Starting from $\widehat{x}_{n}$, run Algorithm 4 with starting time $T-n/L$ using step sizes $h_{\mathrm{pred}}$ for all $N$ steps, with $N=\frac{1}{Lh_{\mathrm{pred}}}$, so that the total time is $1/L$. Let the result be $\widehat{x}_{n+1}^{\prime}$.
(b)
Starting from $\widehat{x}_{n+1}^{\prime}$, run Algorithm 2 for total time $T_{\mathrm{corr}}$ with step size $h_{\mathrm{corr}}$ and score estimate $\widehat{s}_{T-(n+1)/L}$ to obtain $\widehat{x}_{n+1}$.
- 3.
Starting from $\widehat{x}_{N_{0}}$, run Algorithm 4 with starting time $T-N_{0}/L$ using step sizes $h_{\mathrm{pred}}/2,h_{\mathrm{pred}}/4,h_{\mathrm{pred}}/8,\dots,\delta$ to obtain $\widehat{x}_{N_{0}+1}^{\prime}$.
- 4.
Starting from $\widehat{x}^{\prime}_{N_{0}+1}$, run Algorithm 2 for total time $T_{\mathrm{corr}}$ with step size $h_{\mathrm{corr}}$ and score estimate $\widehat{s}_{\delta}$ to obtain $\widehat{x}_{N_{0}+1}$.
- 5.
Return $\widehat{x}_{N_{0}+1}$.

### 3.2 Predictor-corrector framework

The general framework of our algorithm closely follows that of [], which proposed to run the (discretized) reverse ODE but interspersed with “corrector” steps given by running underdamped Langevin dynamics. The idea is that the “predictor” steps where the discretized reverse ODE is being run keep the sampler close to the true reverse process in Wasserstein distance, but they cannot be run for too long before potentially incurring exponential blowups. The main purpose of the corrector steps is then to inject stochasticity into the trajectory of the sampler in order to convert closeness in Wasserstein to closeness in KL divergence. This effectively allows one to “restart the coupling” used to control the predictor steps. For technical reasons that are inherited from [], for most of the reverse process the predictor steps (Step 2(a)) are run with a fixed step size, but at the end of the reverse process (Step 3), they are run with exponentially decaying step sizes.

We follow the same framework, and the core of our result lies in refining the algorithm and analysis for the predictor steps by using the randomized midpoint method. Below, we highlight our key technical steps.

### 3.3 Predictor step – improved discretization error with randomized midpoints

Here, we explain the main idea behind why randomized midpoint allows us to achieve improved dimension dependence. We first focus on the analysis of the predictor (Algorithm 1) and restrict our attention to running the reverse process for a small amount of time $h\ll 1/L$.

We begin by recalling the dimension dependence achieved by the standard exponential integrator scheme. One can show (see e.g. Lemma 4 in []) that if the true reverse process and the discretized reverse process are both run for small time $h$ starting from the same initialization, the two processes drift by a distance of $O(d^{1/2}h^{2})$. By iterating this coupling $O(1/h)$ times, we conclude that in an $O(1)$ window of time, the processes drift by a distance of $O(d^{1/2}h)$. To ensure this is not too large, one would take the step size $h$ to be $O(1/\sqrt{d})$, thus obtaining an iteration complexity of $O(\sqrt{d})$ as in [].

The starting point in the analysis of randomized midpoint is to instead track the *squared* displacement between the two processes instead. Given two neighboring time steps $t-h$ and $t$ in the algorithm, let $x_{t}$ denote the true reverse process at time $t$, and let $\widehat{x}_{t}$ denote the algorithm at time $t$ (in the notation of Algorithm 2, this is $\widehat{x}_{n}$ for some $n$, but we use $t$ in the discussion here to make the comparison to the true reverse process clearer). Note that $\widehat{x}_{t}$ depends on the choice of randomized midpoint $\alpha$ (see Step 1(b)). One can bound the squared displacement $\operatorname*{\mathbb{E}}\,\left\|x_{t}-\widehat{x}_{t}\right\|^{2}$ as follows. Let $y_{t}$ be the result of running the reverse process for time $h$ starting from $\widehat{x}_{t-h}$. Then by writing $x_{t}-\widehat{x}_{t}$ as $(x_{t}-y_{t})-(\widehat{x}_{t}-y_{t})$ and applying Young’s inequality, we obtain

$$ $\operatorname*{\mathbb{E}}_{\widehat{x}_{t-h},\alpha}\|x_{t}-\widehat{x}_{t}\| ^{2}\leq\Bigl{(}1+\frac{Lh}{2}\Bigr{)}\operatorname*{\mathbb{E}}_{\widehat{x}_ {t-h}}\|x_{t}-y_{t}\|^{2}+\frac{2}{Lh}\operatorname*{\mathbb{E}}_{\widehat{x}_ {t-h}}\|\operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{t}-y_{t}\|^{2}+ \operatorname*{\mathbb{E}}_{\widehat{x}_{t-h}}\operatorname*{\mathbb{E}}_{ \alpha}\|\widehat{x}_{t}-y_{t}\|^{2}\,.$ (13) $$

For the first term, because $x_{t}$ and $y_{t}$ are the result of running the same ODE on initializations $x_{t-h}$ and $\widehat{x}_{t-h}$ , the first term is close to $\operatorname*{\mathbb{E}}\|x_{t-h}-\widehat{x}_{t-h}\|^{2}$ provided $h\ll 1/L$. The upshot is that the squared displacement at time $t$ is at most the squared displacement at time $t-h$ plus the remaining two terms on the right of Equation 13.

The main part of the proof lies in bounding these two terms, which can be thought of as “bias” and “variance” terms respectively. The variance term can be shown to scale with the square of the aforementioned $O(d^{1/2}h^{2})$ displacement bound that arises in the exponential integrator analysis, giving $O(dh^{4})$:

###### Lemma 3.1 (Informal, see Lemma A.4 for formal statement) .

If $h\lesssim\frac{1}{L}$ and $T-t\geq(T-t-h)/2$, then

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\widehat{x}_{t-h}}\operatorname*{ \mathbb{E}}_{\alpha}\left\|\widehat{x}_{t}-y_{t}\right\|^{2}\lesssim L^{2}dh^{ 4}\left(L\lor\frac{1}{T-(t-h)}\right)+h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2} h^{2}\operatorname*{\mathbb{E}}_{\widehat{x}_{t-h}}\|x_{t-h}-\widehat{x}_{t-h} \|^{2}\,.$ $$

Note that in this bound, in addition to the $O(dh^{4})$ term and a term for the score estimation error, there is an additional term which depends on the squared displacement from the previous time step. Because the prefactor $L^{2}h^{2}$ is sufficiently small, this will ultimately be negligible.

The upshot of the above Lemma is that if the bias term is of lower order, then this means that the squared displacement essentially increases by $O(dh^{4})$ with every time step of length $h$. Over $O(1/h)$ such steps, the total squared displacement is $O(dh^{3})$, so if we take the step size $h$ to be $O(1/d^{1/3})$, this suggests an improved iteration complexity of $O(d^{1/3})$.

Arguing that the bias term $\frac{2}{Lh}\operatorname*{\mathbb{E}}_{\widehat{x}_{t-h}}\left\|\operatorname
*{\mathbb{E}}_{\alpha}\widehat{x}_{t}-y_{t}\right\|^{2}$ is dominated by the variance term is where it is crucial that we use randomized midpoint instead of exponential integrator. But recall that the randomized midpoint method was engineered so that it would give an unbiased estimate for the true solution to the reverse ODE if the estimate of the trajectory at the randomized midpoint were exact. In reality we only have an approximation to the latter, but as we show, the error incurred by this is indeed of lower order (see Lemma A.3). One technical complication that arises here is that the relevant quantity to bound is the distance between the true process at the randomized midpoint versus the algorithm, when both are initialized at an intermediate point in the algorithm’s trajectory. Bounding such quantities in expectation over the randomness of the algorithm’s trajectory can be difficult, but our proof identifies a way of “offloading” some of this difficulty by absorbing some excess terms into a term of the form $\|x_{t-h}-\widehat{x}_{t-h}\|^{2}$, i.e. the squared displacement from the previous time step. Concretely, we obtain the following bound on the bias term:

###### Lemma 3.2 (Informal, see Lemma A.2 for formal statement) .

$$ $\operatorname*{\mathbb{E}}_{\widehat{x}_{t-h}}\|\operatorname*{\mathbb{E}}_{ \alpha}\widehat{x}_{t}-y_{t}\|^{2}\lesssim L^{4}dh^{6}\left(L\lor\frac{1}{T-t+ h}\right)+h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{4}h^{4}\operatorname*{\mathbb{ E}}_{\widehat{x}_{t-h}}\|x_{t-h}-\widehat{x}_{t-h}\|^{2}$ $$

###### Lemma 3.1 (Informal, see Lemma A.4 for formal statement) .

###### Lemma 3.2 (Informal, see Lemma A.2 for formal statement) .

### 3.4 Shortening the corrector steps

While we have outlined how to improve the predictor step in the framewok of [], it is quite unclear whether the same can be achieved for the corrector step. Whereas the the former is geared towards closeness in Wasserstein distance, the latter is geared towards closeness in KL divergence, and it is a well-known open question in the log-concave sampling literature to obtain analogous discretization bounds in KL for the randomized midpoint method [].

We will sidestep this issue and argue that even using exponential integrator discretization of the underdamped Langevin dynamics will suffice for our purposes, by simply shortening the amount of time for which each corrector step is run.

First, let us briefly recall what was shown in [] for the corrector step. If one runs underdamped Langevin dynamics with stationary distribution $q$ for time $T$ and exponential integrator discretization with step size $h$ starting from two distributions $p$ and $q$, then the resulting distributions $p^{\prime}$ and $q$ satisfy

$$ $\mathsf{TV}(p^{\prime},q)\lesssim\frac{W_{2}(p,q)}{L^{1/4}T^{3/2}}+L^{3/4}T^{1 /2}d^{1/2}h\,,$ (14) $$

where $L$ is the Lipschitzness of $\nabla\ln q$ (see Theorem A.6). At first glance this appears insufficient for our purposes: because of the $d^{1/2}h$ term coming from the discretization error, we would need to take step size $h=1/\sqrt{d}$, which would suggest that the number of iterations must scale with $\sqrt{d}$.

To improve the dimension dependence for our overall predictor-corrector algorithm, we observe that if we take $T$ itself to be smaller, then we can take $h$ to be larger while keeping the discretization error in Equation 14 sufficiently small. Of course, this comes at a cost, as $T$ also appears in the term $\frac{W_{2}(p,q)}{L^{1/4}T^{3/2}}$ in Equation 14. But in our overall proof, the $W_{2}(p,q)$ term is bounded by the predictor analysis. There, we had quite a bit of slack: even with step size as large as $1/d^{1/3}$, we could achieve small Wasserstein error. By balancing appropriately, we get our improved dimension dependence.

### 3.5 Parallel algorithm

Now, we summarize the main proof ideas for our parallel sampling result. In Section 2.3, we described how to approximately solve the reverse ODE over time $h$ by running $K$ rounds of the iteration in Equation 10. In our final algorithm, we will take $h$ to be *dimension-independent*, namely $h=\Theta(1/\sqrt{L})$, so that the main part of the proof is to bound the discretization error incurred over each of these time windows of length $h$. As in the sequential analysis, we will interleave these “predictor” steps with corrector steps given by (parallelized) underdamped Langevin dynamics.

We begin by describing the parallel predictor step. Suppose we have produced an estimate for the reverse process at $t_{0}$ and now wish to solve the ODE from time $t_{0}$ to $t_{0}+h$. We initialize at $\{\widehat{x}^{(0)}_{t_{0}+\alpha_{i}h}\}_{i\in[R]}$ via exponential integrator steps starting from the beginning of the window – see Line 1(c) in Algorithm 9 (this can be thought of as the analogue of Equation 7 used in the sequential algorithm). The key difference relative to the sequential algorithm is that here, because the length of the window is dimension-free, the discretization error incurred by this initialization is too large and must be refined using the fixed point iteration in Equation 10. The main step is then to show that with each iteration of Equation 10, the distance to the true reverse process contracts:

###### Lemma 3.3 (Informal, see Lemma B.2 for formal statement) .

Suppose $h\lesssim 1/L$. If $y_{t}$ denotes the solution of the true ODE starting at $\widehat{x}_{t_{0}}$ and running until time $t_{0}+\alpha_{i}h$, then for all $k\in\{1,\cdots K\}$ and $i\in\{1,\cdots,R\}$,

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\hat{x}_{t_{0}},\alpha_{1},\cdots \alpha_{R}}\left\|\widehat{x}^{(k)}_{t_{0}+\alpha_{i}h}-y_{t_{0}+\alpha_{i}h} \right\|^{2}$ $\displaystyle\lesssim\left(8h^{2}L^{2}\right)^{k}\cdot\left(\frac{1}{R}\sum_{j =1}^{R}\operatorname*{\mathbb{E}}_{\hat{x}_{t_{0}},\alpha_{j}}\left\|\widehat{ x}^{(0)}_{t_{0}+\alpha_{j}h}-y_{t_{0}+\alpha_{j}h}\right\|^{2}\right)$ $\displaystyle+h^{2}\left(\varepsilon^{2}_{\mathrm{sc}}+\frac{L^{2}dh^{2}}{R^{2 }}(L\lor\frac{1}{T-t_{0}+h})+L^{2}\cdot\operatorname*{\mathbb{E}}_{\hat{x}_{t_ {0}}}\left\|\widehat{x}_{t_{0}}-x_{t_{0}}\right\|^{2}\right)\,,$ (15) $$

where $\widehat{x}_{t_{0}}$ is the iterate of the algorithm from the previous time window, and $x_{t_{0}}$ is the corresponding iterate in the true ODE.

In particular, because $h$ is at most a small multiple of $1/L$, the prefactor $(8h^{2}L^{2})^{k}$ is exponentially decaying in $k$, so that the error incurred by the estimate $\widehat{x}^{(k)}_{t_{0}+\alpha_{i}h}$ is contracting with each fixed point iteration. Because the initialization is at distance $\mathrm{poly}(d)$ from the true process, $O(\log d)$ rounds of contraction thus suffice, which translates to $O(\log d)$ parallel rounds for the sampler. The rest of the analysis of the predictor step is quite similar to the analogous proofs for the sequential algorithm (i.e. Lemma B.4 and Lemma B.5 give the corresponding bias and variance bounds).

One shortcoming of the predictor analysis is that the contraction achieved by fixed point iteration ultimately bottoms out at error which scales with $d/R^{2}$ (see the second term in Equation 15). In order for the discretization error to be sufficiently small, we thus have to take $R$, and thus the total work of the algorithm, to scale with $O(\sqrt{d})$. So in this case we do not improve over the dimension dependence of [], and instead the improvement is in obtaining a parallel algorithm.

For the corrector analysis, we mostly draw upon the recent work of [] which analyzed a parallel implementation of the underdamped Langevin dynamics. While their guarantee focuses on sampling from log-concave distributions, implicit in their analysis is a bound for general smooth distributions on how much the law of the algorithm and the law of the true process drift apart in a bounded time window (see Lemma B.8). This bound suffices for our analysis of the corrector step, and we can conclude the following:

###### Theorem 3.4 (Informal, see Theorem B.12 ) .

Let $\beta\geq 1$ be an adjustable parameter. Let $p^{\prime}$ denote the law of the output of running the parallel corrector (see Algorithm 8) for total time $1/\sqrt{L}$ and step size $h$, using an $\varepsilon_{\mathrm{sc}}$-approximate estimate for $\nabla\ln q$ and starting from a sample from another distribution $p$.

$$ $\textup{{TV}}(p^{\prime},q)\lesssim\sqrt{\textup{{KL}}(p^{\prime},q)}\lesssim \frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}}+\frac{\varepsilon}{\beta}+\frac{ \varepsilon}{\beta\sqrt{d}}\cdot W_{2}(p,q)\,.$ $$

Furthermore, this algorithm uses $\widetilde{\Theta}(\beta\sqrt{d}/\varepsilon)$ score evaluations over $\Theta(\log(\beta^{2}d/\varepsilon^{2}))$ parallel rounds.

Overall, parallel algorithm is somewhat different from the parallel sampler developed in the empirical work of [], even apart from the fact that we use randomized midpoint discretization and corrector steps. The reason is that our algorithm applies collocation to fixed windows of time, whereas the algorithm of [] utilizes a sliding window approach that proactively shifts the window forward as soon as the iterates at the start of the previous window begin to converge. We leave rigorously analyzing the benefits of this approach as an interesting future direction.

###### Lemma 3.3 (Informal, see Lemma B.2 for formal statement) .

###### Theorem 3.4 (Informal, see Theorem B.12 ) .

### 3.6 Log-concave sampling in total variation

Finally, we briefly summarize the simple proof for our result on log-concave sampling in TV, which achieves the best known dimension dependence of $\widetilde{O}(d^{5/12})$. Our main observation is that Shen and Lee’s randomized midpoint method [] applied to the underdamped Langevin process gives a *Wasserstein* guarantee for log-concave sampling, while the corrector step of [] can convert a Wasserstein guarantee to closeness in TV. Thus, we can simply run the randomized midpoint method, followed by the corrector step to achieve closeness in TV. Carefully tuning the amount of time spend and step sizes for each phase of this algorithm yields our improved dimension dependence – see Appendix C for the full proof.

## 4 Discussion and Future Work

In this work, we showed that it is possible to leverage Shen and Lee’s randomized midpoint method [] to achieve the best known dimension dependence for sampling from arbitrary smooth distributions in TV using diffusion. We also showed how to parallelize our algorithm, and showed that $\widetilde{O}(\log^{2}d)$ parallel rounds suffice for sampling. These constitute the first provable guarantees for parallel sampling with diffusion models. Finally, we showed that our techniques can be used to obtain an improved dimension dependence for log-concave sampling in TV.

We note that relative to [], our result requires a slightly stronger guarantee on the score estimation error, by a $d^{1/12}$ factor; we believe this is an artifact of our analysis, and it would be interesting to remove this dependence in future work. Importantly however, it was not known how to achieve an improvement over the $O(\sqrt{d})$ dependence shown in that paper even in case that the scores are known *exactly* prior to the present work. Moreover, another line of work [, , ] analyzing diffusion sampling makes the *stronger* assumption that the score estimation error is $\widetilde{O}\left(\frac{\varepsilon}{\sqrt{d}}\right)$, an assumption stronger than ours by a $d^{5/12}$ factor; this does not detract from the importance of these works, and we feel the same is true in our case.

We also note that our diffusion results require smoothness assumptions – we assume that the true score, as well as our score estimates are $L$-Lipschitz. Although this assumption is standard in the literature, recent work [, ] has analyzed DDPM in the absence of these assumptions, culminating in a $\widetilde{O}(d)$ dependence for sampling using a discretization of the reverse SDE. However, unlike in the smooth case, it is not known whether even a *sublinear* in $d$ dependence is possible without smoothness assumptions via any algorithm. We leave this as an interesting open question for future work.

## Acknowledgements

S.C. thanks Sinho Chewi for a helpful discussion on log-concave sampling. S.G. was funded by NSF Award CCF-1751040 (CAREER) and the NSF AI Institute for Foundations of Machine Learning (IFML).

## Appendix A Sequential algorithm

In this section, we describe our sequential randomized-midpoint-based algorithm in detail. Following the framework of [], we begin by describing the *predictor Step* and show in Lemma A.5 that in $\widetilde{O}(d^{1/3})$ steps (ignoring other dependencies), when run for time $t$ at most $O(\frac{1}{L})$ starting from $t_{n}$, it produces a sample that is close to the true distribution at time $t_{n}-t$. Then, we show that the *corrector step* can be used to convert our $W_{2}$ error to error in TV distance by running the underdamped Langevin Monte Carlo algorithm, as described in []. We show in A.8 that if we run our predictor and corrector steps in succession for a careful choice of times, we obtain a sample that is close to the true distribution in TV using just $\widetilde{O}(d^{5/12})$ steps, but covering a time $O\left(\frac{1}{L}\right)$. Finally, in Theorem A.10, we iterate this bound $\widetilde{O}(\log^{2}d)$ times to obtain our final iteration complexity of $\widetilde{O}(d^{5/12})$.

### A.1 Predictor step

To show the $\widetilde{O}(d^{1/3})$ dependence on dimension for the predictor step, we will, roughly speaking, show that its bias after *one step* is bounded by $\approx O_{L}\left(dh^{6}\right)$ in Lemma A.3, and that the variance is bounded by $O_{L}\left(dh^{4}\right)$ in Lemma A.4. Then, iterating these bounds $\approx\frac{1}{h}$ times as shown in Lemma A.5 will give error $O_{L}\left(dh^{4}+dh^{3}\right)$ in squared Wasserstein Distance.

Figure: Algorithm 4 PredictorStep (Sequential)

###### Lemma A.1 (Naive ODE Coupling) .

Consider two variables $x_{0},x_{0}^{\prime}$ starting at time $t_{0}$, and consider the result of running the true ODE for time $h$, and let the results be $x_{1},x_{1}^{\prime}$. For $L\geq 1$, $h\leq 1/L$, we have

$$ $\displaystyle\|x_{1}-x_{1}^{\prime}\|^{2}\leq\exp(O(Lh))\|x_{0}-x_{0}^{\prime} \|^{2}$ $$

###### Proof.

Recall that the true ODE is given by

$$ $\displaystyle dx_{t}=(x_{t}+\nabla\ln q_{T-t}(x_{t}))dt$ $$

So,

$$ $\displaystyle\partial_{t}\|x_{t}-x_{t}^{\prime}\|^{2}$ $\displaystyle=2\langle x_{t}-x_{t}^{\prime},\partial_{t}x_{t}-\partial_{t}x_{t }^{\prime}\rangle$ $\displaystyle=2\langle x_{t}-x_{t}^{\prime},x_{t}-x_{t}^{\prime}+\nabla\ln q_{ T-t}(x_{t})-\nabla\ln q_{T-t}(x_{t}^{\prime})\rangle$ $\displaystyle\lesssim L\|x_{t}-x_{t}^{\prime}\|^{2}$ $$

So,

$$ $\displaystyle\|x_{1}-x_{1}^{\prime}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh)\right)\|x_{0}-x_{0}^{\prime}\|^{2}\,.\qed$ $$

###### Lemma A.2 .

Suppose $L\geq 1$. In Algorithm 4, for all $n\in\{0,\dots,N-1\}$, let $x^{*}_{n}(t)$ be the solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$, running until time $t_{n}-t$. If $h_{n}\lesssim\frac{1}{L}$ and $t_{n}-h_{n}\geq t_{n}/2$, we have

$$ $\operatorname*{\mathbb{E}}\|h_{n}e^{(1-\alpha)h_{n}}\widehat{s}_{t_{n}-\alpha h _{n}}(\widehat{x}_{n+\frac{1}{2}})-h_{n}e^{(1-\alpha)h_{n}}\nabla\ln q_{t_{n}- \alpha h_{n}}(x_{n}^{*}(\alpha h_{n}))\|^{2}\\ \lesssim h_{n}^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{4}dh_{n}^{6}\left(L\lor \frac{1}{t_{n}}\right)+L^{4}h_{n}^{4}\operatorname*{\mathbb{E}}\|\widehat{x}_{ n}-x_{t_{n}}\|^{2}\,,$ $$

where $\operatorname*{\mathbb{E}}$ refers to the expectation over the initial choice $\widehat{x}_{0}\sim q_{t_{0}}$.

###### Proof.

For the proof, we will let $h:=h_{n}$. It suffices to show that

$$ $\|\widehat{s}_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-\nabla\ln q_{t_{n} -\alpha h}(x_{n}^{*}(\alpha h))\|^{2}\lesssim\varepsilon_{\mathrm{sc}}^{2}+L^{ 4}dh^{4}\left(L\lor\frac{1}{t_{n}}\right)+L^{4}h^{2}\operatorname*{\mathbb{E}} \|\widehat{x}_{n}-x_{t_{n}}\|^{2}\,.$ (16) $$

Now,

$$ $\displaystyle\|\widehat{s}_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})- \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))\|^{2}$ $\displaystyle\lesssim\|\widehat{s}_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2} })-\nabla\ln q_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})\|^{2}+\|\nabla\ln q _{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-\nabla\ln q_{t_{n}-\alpha h}(x_ {n}^{*}(\alpha h))\|^{2}$ $\displaystyle\lesssim\varepsilon_{\mathrm{sc}}^{2}+L^{2}\|\widehat{x}_{n+\frac {1}{2}}-x_{n}^{*}(\alpha h)\|^{2}\,.$ (17) $$

Now, note $x_{n}^{*}(\alpha h)$ is the solution to the following ODE run for time $\alpha h$, starting at $\widehat{x}_{n}$ at time $t_{n}$:

$$ $dx_{t}=\left(x_{t}+\nabla\ln q_{t}(x_{t})\right)dt$ $$

Similarly, $\widehat{x}_{n+\frac{1}{2}}$ is the solution to the following ODE run for time $\alpha h$, starting at $\widehat{x}_{n}$ at time $t_{n}$:

$$ $d\widehat{x}_{t}=\left(\widehat{x}_{t}+\widehat{s}_{t_{n}}(\widehat{x}_{n}) \right)dt$ $$

So, we have

$$ $\displaystyle\partial_{t}\|x_{t}-\widehat{x}_{t}\|^{2}$ $\displaystyle=2\langle x_{t}-\widehat{x}_{t},\partial_{t}x_{t}-\partial_{t}{ \widehat{x}}_{t}\rangle$ $\displaystyle=2\left(\|x_{t}-\widehat{x}_{t}\|^{2}+\langle x_{t}-\widehat{x}_{ t},\nabla\ln q_{t}(x_{t})-\widehat{s}_{t_{n}}(\widehat{x}_{n})\rangle\right)$ $\displaystyle\leq\left(2+\frac{1}{h}\right)\|x_{t}-\widehat{x}_{t}\|^{2}+h\| \nabla\ln q_{t}(x_{t})-\widehat{s}_{t_{n}}(\widehat{x}_{n})\|^{2}$ $$

where the last line is by Young’s inequality.
So, by Grönwall’s inequality,

$$ $\displaystyle\|x_{n}^{*}(\alpha h)-\widehat{x}_{n+\frac{1}{2}}\|^{2}$ $\displaystyle\leq\exp\left(\left(2+\frac{1}{h}\right)\cdot\alpha h\right)\int_ {0}^{h}h\|\nabla\ln q_{t_{n}-s}(x_{n}^{*}(s))-\widehat{s}_{t_{n}}(\widehat{x}_ {n})\|^{2}\ \mathrm{d}s$ $\displaystyle\lesssim h\int_{0}^{h}\|\nabla\ln q_{t_{n}-s}\left(x_{n}^{*}(s) \right)-\widehat{s}_{t_{n}}(\widehat{x}_{n})\|^{2}\ \mathrm{d}s$ $\displaystyle\lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+h\int_{0}^{h}\|\nabla \ln q_{t_{n}-s}(x_{n}^{*}(s))-\nabla\ln q_{t_{n}}(\widehat{x}_{n})\|^{2}\ \mathrm{d}s$ $$

Now, we have

$$ $\displaystyle\|\nabla\ln q_{t_{n}-s}(x_{n}^{*}(s))-\nabla\ln q_{t_{n}}( \widehat{x}_{n})\|^{2}$ $\displaystyle\lesssim\|\nabla\ln q_{t_{n}-s}(x_{t_{n}-s})-\nabla\ln q_{t_{n}}( x_{t_{n}})\|^{2}$ $\displaystyle\qquad+\|\nabla\ln q_{t_{n}-s}(x_{n}^{*}(s))-\nabla\ln q_{t_{n}-s }(x_{t_{n}-s})\|^{2}+\|\nabla\ln q_{t_{n}}(x_{t_{n}})-\nabla\ln q_{t_{n}}( \widehat{x}_{n})\|^{2}\,.$ $$

By Corollary D.1,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-s}(x_{t_{n}-s})- \nabla\ln q_{t_{n}}(x_{t_{n}})\|^{2}\lesssim L^{2}dh^{2}\left(L\lor\frac{1}{t_ {n}}\right)$ $$

By Lipschitzness of $\nabla\ln q_{t}$ and Lemma A.1, for $s\leq h_{n}\leq 1/L$,

$$ $\displaystyle\|\nabla\ln q_{t_{n}-s}(x_{n}^{*}(s))-\nabla\ln q_{t_{n}-s}(x_{t_ {n}-s})\|^{2}$ $\displaystyle\leq L^{2}\|x_{n}^{*}(s)-x_{t_{n}-s}\|^{2}$ $\displaystyle\lesssim L^{2}\exp\left(O(Lh_{n})\right)\|\widehat{x}_{n}-x_{t_{n }}\|^{2}$ $\displaystyle\lesssim L^{2}\|\widehat{x}_{n}-x_{t_{n}}\|^{2}$ $$

and similarly,

$$ $\displaystyle\|\nabla\ln q_{t_{n}}(x_{t_{n}})-\nabla\ln q_{t_{n}}(\widehat{x}_ {n})\|^{2}\leq L^{2}\|\widehat{x}_{n}-x_{t_{n}}\|^{2}$ $$

So, we have shown that

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-s}(x_{n}^{*}(s))- \nabla\ln q_{t_{n}}(\widehat{x}_{n})\|^{2}\lesssim L^{2}dh^{2}\left(L\lor\frac {1}{t_{n}}\right)+L^{2}\operatorname*{\mathbb{E}}\|\widehat{x}_{n}-x_{t_{n}}\| ^{2}$ $$

so that

$$ $\displaystyle\operatorname*{\mathbb{E}}\|x_{n}^{*}(\alpha h)-\widehat{x}_{n+ \frac{1}{2}}\|^{2}\lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2}dh^{4}\left (L\lor\frac{1}{t_{n}}\right)+L^{2}h^{2}\operatorname*{\mathbb{E}}\|\widehat{x} _{n}-x_{t_{n}}\|^{2}\,.$ $$

Combining this with the bound in (17) and recalling that $h\lesssim 1/L$ yields the desired inequality in (16).
∎

###### Lemma A.3 (Sequential Predictor Bias) .

Suppose $L\geq 1$. In Algorithm 4, for all $n\in\{0,\dots,N-1\}$, let $x_{n}^{*}(t)$ be the solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$, and let $x_{t}\sim q_{t}$ be the solution of the true ODE, starting at $\widehat{x}_{0}\sim q_{t_{0}}$. If $h_{n}\lesssim\frac{1}{L}$, and $t_{n}-h_{n}\geq t_{n}/2$, we have

$$ $\operatorname*{\mathbb{E}}\|\operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{n+ 1}-x_{n}^{*}(h_{n})\|^{2}\lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{4}dh^{ 6}\left(L\lor\frac{1}{t_{n}}\right)+L^{4}h^{4}\operatorname*{\mathbb{E}}\| \widehat{x}_{n}-x_{t_{n}}\|^{2}$ $$

where $\operatorname*{\mathbb{E}}_{\alpha}$ is the expectation with respect to the $\alpha$ chosen in the $n^{th}$ step, and $\operatorname*{\mathbb{E}}$ is the expectation with respect to the choice of the initial $\widehat{x}_{0}\sim q_{t_{0}}$.

###### Proof.

For the proof, we wil fix $n$, and let $h:=h_{n}$.
By the integral formulation of the true ODE,

$$ $x_{n}^{*}(h)=e^{h}\widehat{x}_{n}+\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n}-h)}\nabla \ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\,.$ $$

Thus, we have

$$ $\displaystyle\|\operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{n+1}-x_{n}^{*}( h)\|^{2}$ $\displaystyle=\|h\operatorname*{\mathbb{E}}_{\alpha}e^{(1-\alpha)h}\widehat{s} _{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-\int_{t_{n}-h}^{t_{n}}e^{s-(t_{ n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}_{\alpha}\|he^{(1-\alpha)h} \widehat{s}_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-he^{(1-\alpha)h} \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))\|^{2}$ $\displaystyle\qquad+\|h\cdot\operatorname*{\mathbb{E}}_{\alpha}e^{(1-\alpha)h} \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))-\int_{t_{n}-h}^{t_{n}}e^{s-( t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}\,.$ $$

The second term is 0 0 since

$$ $\displaystyle h\operatorname*{\mathbb{E}}_{\alpha}e^{(1-\alpha)h}\nabla\ln q_{ t_{n}-\alpha h}(x_{n}^{*}(\alpha h))$ $\displaystyle=h\int_{0}^{1}e^{(1-\alpha)h}\nabla\ln q_{t_{n}-\alpha h}(x_{n}^{ *}(\alpha h))\,\mathrm{d}\alpha$ $\displaystyle=\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t _{n}-s))\,\mathrm{d}s\,.$ $$

For the first term, we have, by Lemma A.2

$$ $\operatorname*{\mathbb{E}}\|he^{(1-\alpha)h}\widehat{s}_{t_{n}-\alpha h}( \widehat{x}_{n+\frac{1}{2}})-he^{(1-\alpha)h}\nabla\ln q_{t_{n}-\alpha h}(x_{n }^{*}(\alpha h))\|^{2}\\ \lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{4}dh^{6}\left(L\lor\frac{1}{t_{ n}}\right)+L^{4}h^{4}\operatorname*{\mathbb{E}}\|\widehat{x}_{n}-x_{t_{n}}\|^{ 2}\,.$ $$

The claimed bound follows.
∎

###### Lemma A.4 (Sequential Predictor Variance) .

Suppose $L\geq 1$. In Algorithm 4, for all $n\in\{0,\dots,N-1\}$, let $x_{n}^{*}(t)$ be the solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$, and let $x_{t}\sim q_{t}$ be the solution of the true ODE starting at $\widehat{x}_{0}\sim q_{t_{0}}$. If $h_{n}\lesssim\frac{1}{L}$ and $t_{n}-h_{n}\geq t_{n}/2$, we have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\widehat{x}_{n+1}-x_{n}^{*}(h_{n})\|^ {2}\lesssim h_{n}^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2}dh_{n}^{4}\left(L\lor \frac{1}{t_{n}}\right)+L^{2}h_{n}^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}}- \widehat{x}_{n}\|^{2}$ $$

where $\operatorname*{\mathbb{E}}$ refers to the expectation wrt the random $\alpha$ in the $n^{th}$ step, along with the initial choice $\widehat{x}_{0}\sim q_{t_{0}}$.

###### Proof.

Fix $n$ and let $h:=h_{n}$. We have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\widehat{x}_{n+1}-x_{n}^{*}(h)\|^{2}$ $\displaystyle=\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h}\widehat{s}_{t _{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n}- h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h} \widehat{s}_{t_{n}-\alpha h}(\widehat{x}_{n+\frac{1}{2}})-he^{(1-\alpha)h} \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))\|^{2}$ $\displaystyle\qquad\qquad+\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h} \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))-\int_{t_{n}-h}^{t_{n}}e^{(1- \alpha)h}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $\displaystyle\qquad\qquad+\operatorname*{\mathbb{E}}\|\int_{t_{n}-h}^{t_{n}}e^ {(1-\alpha)h}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s-\int_{t_{n}-h}^{ t_{n}}e^{s-(t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $$

The first term was bounded in Lemma A.2:

$$ $\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h}\widehat{s}_{t_{n}-\alpha h} (\widehat{x}_{n+\frac{1}{2}})-he^{(1-\alpha)h}\nabla\ln q_{t_{n}-\alpha h}(x_{ t_{n}-\alpha h})\|^{2}\\ \lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{4}dh^{6}\left(L\lor\frac{1}{t_{ n}}\right)+L^{4}h^{4}\operatorname*{\mathbb{E}}\|\widehat{x}_{n}-x_{t_{n}}\|^{ 2}\,.$ $$

For the second term,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h}\nabla\ln q_{t_ {n}-\alpha h}(x_{n}^{*}(\alpha h))-\int_{t_{n}-h}^{t_{n}}e^{(1-\alpha)h}\nabla \ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $\displaystyle=\operatorname*{\mathbb{E}}\|\int_{t_{n}-h}^{t_{n}}e^{(1-\alpha)h }\cdot\left(\nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))-\nabla\ln q_{s}( x_{n}^{*}(t_{n}-s))\right)\,\mathrm{d}s\|^{2}$ $\displaystyle\lesssim h\int_{t_{n}-h}^{t_{n}}\operatorname*{\mathbb{E}}\| \nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))-\nabla\ln q_{s}(x_{n}^{*}(t_ {n}-s)\|^{2}\,\mathrm{d}s\,.$ $$

Now,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-\alpha h}(x_{n}^{* }(\alpha h))-\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\|^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-\alpha h}( x_{t_{n}-\alpha h})-\nabla\ln q_{s}(x_{s})\|^{2}$ $\displaystyle\qquad+\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-\alpha h}(x _{t_{n}-\alpha h})-\nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))\|^{2}$ $\displaystyle\qquad+\operatorname*{\mathbb{E}}\|\nabla\ln q_{s}(x_{s})-\nabla \ln q_{s}(x_{n}^{*}(t_{n}-s))\|^{2}$ $$

The first of these terms is bounded in Corollary D.1:

$$ $\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-\alpha h}(x_{t_{n}-\alpha h})- \nabla\ln q_{s}(x_{s})\|^{2}\lesssim L^{2}dh^{2}\left(L\lor\frac{1}{t_{n}}\right)$ $$

For the remaining two terms, note that by the Lipschitzness of $\nabla\ln q_{t}$ and Lemma A.1,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\nabla\ln q_{t_{n}-\alpha h}(x_{t_{n} -\alpha h})-\nabla\ln q_{t_{n}-\alpha h}(x_{n}^{*}(\alpha h))\|^{2}$ $\displaystyle\leq L^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}-\alpha h}-x_{n}^{* }(\alpha h)\|^{2}$ $\displaystyle\lesssim L^{2}\exp(O(Lh))\operatorname*{\mathbb{E}}\|x_{t_{n}}- \widehat{x}_{n}\|^{2}$ $\displaystyle\lesssim L^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}}-\widehat{x}_{ n}\|^{2}$ $$

and similarly, for $t_{n}-h\leq s\leq t_{n}$,

$$ $\operatorname*{\mathbb{E}}\|\nabla\ln q_{s}(x_{s})-\nabla\ln q_{s}(x_{n}^{*}(t _{n}-s))\|^{2}\lesssim L^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}}-\widehat{x}_ {n}\|^{2}\,.$ (18) $$

Thus, we have shown that the second term in our bound on $\operatorname*{\mathbb{E}}\|\widehat{x}_{n+1}-x^{*}_{n}(h)\|^{2}$ is bounded as follows:

$$ $\operatorname*{\mathbb{E}}\|h\cdot e^{(1-\alpha)h}\nabla\ln q_{t_{n}-\alpha h} (x_{n}^{*}(\alpha h))-\int_{t_{n}-h}^{t_{n}}e^{(1-\alpha)h}\nabla\ln q_{s}(x_{ n}^{*}(t_{n}-s))\,,\mathrm{d}s\|^{2}\\ \lesssim L^{2}dh^{4}\left(L\lor\frac{1}{t_{n}}\right)+L^{2}h^{2}\operatorname* {\mathbb{E}}\|x_{t_{n}}-\widehat{x}_{n}\|^{2}\,.$ $$

For the third term,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\int_{t_{n}-h}^{t_{n}}e^{(1-\alpha)h} \nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,,\mathrm{d}s-\int_{t_{n}-h}^{t_{n}}e^{s-( t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,,\mathrm{d}s\|^{2}$ $\displaystyle=\operatorname*{\mathbb{E}}\|\int_{t_{n}-h}^{t_{n}}\left(e^{(1- \alpha)h}-e^{s-(t_{n}-h)}\right)\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,,\mathrm{ d}s\|^{2}$ $\displaystyle\lesssim h\int_{t_{n}-h}^{t_{n}}\operatorname*{\mathbb{E}}_{ \alpha}\left(e^{(1-\alpha)h}-e^{s-(t_{n}-h)}\right)^{2}\operatorname*{\mathbb{ E}}_{\widehat{x}_{0}\sim q_{t_{0}}}\|\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\|^{2} \,\mathrm{d}s\,.$ $$

Now, we have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\| ^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}\|\nabla\ln q_{s}(x_{s})\|^{2}+ \operatorname*{\mathbb{E}}\|\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))-\nabla\ln q_{s }(x_{s})\|^{2}$ $\displaystyle\lesssim\frac{d}{s}+L^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}}- \widehat{x}_{n}\|^{2}\,.$ $$

where the last step follows by Lemma D.4 and (18).
So,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\int_{t_{n}-h}^{t_{n}}e^{(1-\alpha)h} \nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s-\int_{t_{n}-h}^{t_{n}}e^{s-(t _{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\|^{2}$ $\displaystyle\lesssim h\int_{t_{n}-h}^{t_{n}}\operatorname*{\mathbb{E}}_{ \alpha}\left(e^{(1-\alpha)h}-e^{s-(t_{n}-h)}\right)^{2}\cdot\left(\frac{d}{s}+ L^{2}\operatorname*{\mathbb{E}}\|x_{t_{n}}-\widehat{x}_{n}\|^{2}\right)\, \mathrm{d}s$ $\displaystyle\lesssim h^{4}\cdot\left(\frac{d}{t_{n}}+L^{2}\operatorname*{ \mathbb{E}}\|x_{t_{n}}-\widehat{x}_{n}\|^{2}\right)$ $$

Thus, noting that $h\leq\frac{1}{L}$, we obtain the claimed bound on $\operatorname*{\mathbb{E}}\|\widehat{x}_{n+1}-x^{*}_{n}(h_{n})\|^{2}$.
∎

Finally, we put together the bias and variance bounds above to obtain a bound on the Wasserstein error at the end of the Predictor Step.

###### Lemma A.5 (Sequential Predictor Wasserstein Guarantee) .

Suppose that $L\geq 1$, and that for our sequence of step sizes $h_{0},\dots,h_{N-1}$, $\sum_{i}h_{i}\leq 1/L$. Let $h_{\text{max}}=\max_{i}h_{i}$. Then, at the end of Algorithm 4,

- 1.
If $t_{N}\gtrsim 1/L$,
$\displaystyle W_{2}^{2}(\widehat{x}_{N},x_{t_{N}})\lesssim\frac{\varepsilon_{
\mathrm{sc}}^{2}}{L^{2}}+L^{3}dh_{\max}^{4}+L^{2}dh_{\max}^{3}$
- 2.
If $t_{N}\lesssim 1/L$ and $h_{n}\lesssim\frac{t_{n}}{2}$ for each $n$,
$\displaystyle W_{2}^{2}(\widehat{x}_{N},x_{t_{N}})\lesssim\frac{\varepsilon_{
\mathrm{sc}}^{2}}{L^{2}}+\left(L^{3}dh_{\max}^{4}+L^{2}dh_{\max}^{3}\right)\cdot
N$

Here, $x_{t_{N}}\sim q_{t_{N}}$ is the solution of the true ODE beginning at $x_{t_{0}}=\widehat{x}_{0}\sim q_{t_{0}}$.

###### Proof.

For all $n\in[1,\dots,N]$, let $y_{n}$ be the solution of the exact one step ODE starting from $\widehat{x}_{n-1}$. Let the operator $\operatorname*{\mathbb{E}}_{\alpha}$ be the expectation over the random choice of $\alpha$ in the $n^{th}$ iteration. Note that only $\widehat{x}_{N}$ depends on $\alpha$. We have

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\left[\|x_{t_{N}}-\widehat{x}_ {N}\|^{2}\right]$ $\displaystyle=\operatorname*{\mathbb{E}}_{\alpha}\left[\|(x_{t_{N}}-y_{N})-( \widehat{x}_{N}-y_{N})\|^{2}\right]$ $\displaystyle=\|x_{t_{N}}-y_{N}\|^{2}-2\langle x_{t_{N}}-y_{N},\operatorname*{ \mathbb{E}}_{\alpha}\widehat{x}_{N}-y_{N}\rangle+\operatorname*{\mathbb{E}}_{ \alpha}\|\widehat{x}_{N}-y_{N}\|^{2}$ $\displaystyle\leq\left(1+\frac{Lh_{N-1}}{2}\right)\|x_{t_{N}}-y_{N}\|^{2}+ \frac{2}{Lh_{N-1}}\|\operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{N}-y_{N}\| ^{2}+\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{N}-y_{N}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\|x_{t_{N-1}}-\widehat{x}_{N-1}\| ^{2}+\frac{2}{Lh_{N-1}}\|\operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{N}-y_ {N}\|^{2}+\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{N}-y_{N}\|^{2}\,,$ $$

where the third line is by Young’s inequality, and the fourth line is by Lemma A.1.
Taking the expectation wrt $\widehat{x}_{0}\sim q_{t_{0}}$, by Lemmas A.3 and A.4,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}+\frac{2}{Lh_{N-1}}\operatorname*{\mathbb{E}}\| \operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{N}-y_{N}\|^{2}+\operatorname*{ \mathbb{E}}\|\widehat{x}_{N}-y_{N}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}$ $\displaystyle\qquad+O\left(\frac{1}{Lh_{N-1}}\left(h_{N-1}^{2}\varepsilon_{ \mathrm{sc}}^{2}+L^{4}dh_{N-1}^{6}\left(L\lor\frac{1}{t_{N-1}}\right)+L^{4}h_{ N-1}^{4}\operatorname*{\mathbb{E}}\|x_{t_{N-1}}-\widehat{x}_{N-1}\|^{2}\right)\right)$ $\displaystyle\qquad+O\left(h_{N-1}^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2}dh_{N -1}^{4}\left(L\lor\frac{1}{t_{N-1}}\right)+L^{2}h_{N-1}^{2}\operatorname*{ \mathbb{E}}\|x_{t_{N-1}}-\widehat{x}_{N-1}\|^{2}\right)$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}$ $\displaystyle\qquad+O\left(\frac{h_{N-1}\varepsilon_{\mathrm{sc}}^{2}}{L}+h_{N -1}^{2}\varepsilon_{\mathrm{sc}}^{2}+\left(L^{3}dh_{N-1}^{5}+L^{2}dh_{N-1}^{4} \right)\left(L\lor\frac{1}{t_{N-1}}\right)\right)$ $$

By induction, noting that $x_{t_{0}}=\widehat{x}_{0}$, we have

$$ $\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2}\lesssim\sum_{n=0}^ {N-1}\left(\frac{h_{n}\varepsilon_{\mathrm{sc}}^{2}}{L}+h_{n}^{2}\varepsilon_{ \mathrm{sc}}^{2}+\left(L^{3}dh_{n}^{5}+L^{2}dh_{n}^{4}\right)\cdot\left(L\lor \frac{1}{t_{n}}\right)\right)\cdot\exp\left(O\left(L\sum_{i=n+1}^{N-1}h_{i} \right)\right)\,.$ $$

By assumption, $\sum_{i}h_{i}\leq\frac{1}{L}$. In the first case, $L\vee\frac{1}{t_{n}}\lesssim L$ for all $n$, so

$$ $W_{2}^{2}(\widehat{x}_{N},x_{t_{N}})\lesssim\frac{\varepsilon_{\mathrm{sc}}^{2 }}{L^{2}}+L^{3}dh_{\max}^{4}+L^{2}dh_{\max}^{3}\,.$ $$

In the second case,

$$ $\displaystyle W_{2}^{2}(\widehat{x}_{N},x_{t_{N}})$ $\displaystyle\lesssim\frac{\varepsilon_{\mathrm{sc}}^{2}}{L^{2}}+\left(L^{3}dh _{\max}^{4}+L^{2}dh_{\max}^{3}\right)\cdot\sum_{n=0}^{N-1}\frac{h_{n}}{t_{n}}$ $\displaystyle\lesssim\frac{\varepsilon_{\mathrm{sc}}^{2}}{L^{2}}+\left(L^{3}dh _{\max}^{4}+L^{2}dh_{\max}^{3}\right)\cdot N\,.\qed$ $$

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Starting time $t_{0}$, Number of steps $N$, Step sizes $h_{n\in[0,\dots,N-1]}$, Score estimates $\widehat{s}_{t}$

- 1.
For $n=0,\dots,N-1$:
(a)
Let $t_{n}=t_{0}-\sum_{i=0}^{n-1}h_{i}$
(b)
Randomly sample $\alpha$ uniformly from $[0,1]$.
(c)
Let $\widehat{x}_{n+\frac{1}{2}}=e^{\alpha h_{n}}\widehat{x}_{n}+\left(e^{\alpha h_
{n}}-1\right)\widehat{s}_{t_{n}}(\widehat{x}_{n})ds$
(d)
Let $\widehat{x}_{n+1}=e^{h_{n}}\widehat{x}_{n}+h_{n}\cdot e^{(1-\alpha){h_{n}}}
\widehat{s}_{t_{n}-\alpha h_{n}}(\widehat{x}_{n+\frac{1}{2}})$
- 2.
Let $t_{N}=t_{0}-\sum_{i=0}^{N-1}h_{i}$
- 3.
Return $\widehat{x}_{N},t_{N}$.

###### Lemma A.1 (Naive ODE Coupling) .

###### Proof.

###### Lemma A.2 .

###### Proof.

###### Lemma A.3 (Sequential Predictor Bias) .

###### Proof.

###### Lemma A.4 (Sequential Predictor Variance) .

###### Proof.

###### Lemma A.5 (Sequential Predictor Wasserstein Guarantee) .

###### Proof.

### A.2 Corrector step

For the sequential algorithm, we make use of the underdamped Langevin corrector step and analysis from []. We reproduce the same here for convenience.

The underdamped Langevin Monte Carlo process with step size $h$ is given by:

$$ $\displaystyle\begin{split}\mathrm{d}\widehat{x}_{t}&=\widehat{v}_{t}\,\mathrm{ d}t\\ \mathrm{d}\widehat{v}_{t}&=(\widehat{s}(\widehat{x}_{\lfloor t/h\rfloor h})- \gamma\widehat{v}_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}B_{t}\end{split}$ (19) $$

where $B_{t}$ is Brownian motion, and $\widehat{s}$ satisfies

$$ $\operatorname*{\mathbb{E}}_{x\sim q}\left\|\widehat{s}(x)-\nabla\ln q(x)\right \|^{2}\leq\varepsilon_{\mathrm{sc}}^{2}\,.$ (20) $$

for some target measure $q$.
Here, we set the friction parameter $\gamma=\Theta(\sqrt{L})$.

Then, our corrector step is as follows.

Figure: Algorithm 5 CorrectorStep (Sequential)

###### Theorem A.6 (Theorem 5 of [ CCL + 23a ] , restated) .

Suppose Eq. (20) holds. For any distribution $p$ over $\mathbb{R}^{d}$, and total time $T_{\mathrm{corr}}\lesssim 1/\sqrt{L}$, if we let $p_{N}$ be the distribution of $\widehat{x}_{N}$ resulting from running Algorithm 5 initialized at $\widehat{x}_{0}\sim p$, then we have

$$ $\displaystyle\textup{{TV}}(p_{N},q)\lesssim\frac{W_{2}(p,q)}{L^{1/4}T_{\mathrm {corr}}^{3/2}}+\frac{\varepsilon_{\mathrm{sc}}T_{\mathrm{corr}}^{1/2}}{L^{1/4} }+L^{3/4}T_{\mathrm{corr}}^{1/2}d^{1/2}h_{\mathrm{corr}}\,.$ $$

###### Corollary A.7 (Underdamped Corrector) .

For $T_{\mathrm{corr}}=\Theta\left(\frac{1}{\sqrt{L}d^{1/18}}\right)$

$$ $\displaystyle\textup{{TV}}(p_{N},q)\lesssim W_{2}(p,q)\cdot d^{1/12}\cdot\sqrt {L}+\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}d^{1/36}}+\sqrt{L}d^{17/36}h_{ \mathrm{corr}}\,.$ $$

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Total time $T_{\mathrm{corr}}$, Step size $h_{\mathrm{corr}}$, Score estimate $\widehat{s}$

- 1.
Run underdamped Langevin Monte Carlo in (19) for total time $T_{\mathrm{corr}}$ using step size $h_{\mathrm{corr}}$, and let the result be $\widehat{x}_{N}$.
- 2.
Return $\widehat{x}_{N}$.

###### Theorem A.6 (Theorem 5 of [ CCL + 23a ] , restated) .

###### Corollary A.7 (Underdamped Corrector) .

### A.3 End-to-end analysis

Finally, we put together the analysis of the predictor and corrector step to obtain our final $\widetilde{O}(d^{5/12})$ dependence on sampling time. We first show that carefully choosing the amount of time to run the corrector results in small TV error after successive rounds of the predictor and corrector steps in Lemma A.8. Finally, we iterate this bound to obtain our final guarantee, given by Theorem A.10.

Figure: Algorithm 6 SequentialAlgorithm

###### Lemma A.8 (TV error after one round of predictor and corrector) .

Let $x_{t}\sim q_{t}$ be a sample from the true distribution at time $t$. Let $t_{n}=T-n/L$ for $n\in[0,\dots,N_{0}]$. If we set $T_{\mathrm{corr}}=\Theta\left(\frac{1}{\sqrt{L}d^{1/18}}\right)$, we have,

- 1.
For $n\in[0,\dots,N_{0}-1]$, if $t_{n}\gtrsim 1/L$,
$\textup{{TV}}(\widehat{x}_{n+1},x_{t_{n+1}})\leq\textup{{TV}}(\widehat{x}_{n},
x_{t_{n}})+O\left(L^{2}d^{7/12}h_{\mathrm{pred}}^{2}+L^{3/2}d^{7/12}h_{\mathrm
{pred}}^{3/2}+\sqrt{L}d^{17/36}h_{\mathrm{corr}}+\frac{\varepsilon_{\mathrm{sc
}}d^{1/12}}{\sqrt{L}}\right)$
- 2.
If $t_{N_{0}}\lesssim 1/L$,
$\displaystyle\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{\delta})$
$\displaystyle\leq\textup{{TV}}(\widehat{x}_{N_{0}},x_{t_{N_{0}}})$
$\displaystyle\quad+O\left(\left(L^{2}d^{7/12}h_{\mathrm{pred}}^{2}+L^{3/2}d^{7
/12}h_{\mathrm{pred}}^{3/2}\right)\cdot\sqrt{\log\frac{h_{\mathrm{pred}}}{
\delta}}+\sqrt{L}d^{17/36}h_{\mathrm{corr}}+\frac{\varepsilon_{\mathrm{sc}}d^{
1/12}}{\sqrt{L}}\right)$

###### Proof.

For $n\in[0,\ldots,N_{0}]$, let $\widehat{y}_{n+1}$ be the result of a single predictor-corrector sequence as described in step 2 of Algorithm 6, but starting from $x_{t_{n}}\sim q_{t_{n}}$ instead of $\widehat{x}_{n}$. Additionally, let $\widehat{y}_{N_{0}+1}$ be the result of running steps $3$ and $4$ starting from $x_{t_{N_{0}}}\sim q_{t_{N_{0}}}$ instead of $\widehat{x}_{N_{0}}$. Similarly, let $\widehat{y}_{n+1}^{\prime}$ be the result of only applying the predictor step starting from $x_{t_{n}}\sim q_{t_{n}}$, analogous to $\widehat{x}_{n+1}^{\prime}$ defined in step 2a.

We have, by the triangle inequality and the data-processing inequality, for $n\in[0,\dots,N_{0}-1]$,

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{n+1},x_{t_{n+1}})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{n+1},\widehat{y}_{n+1})+\textup{{ TV}}(\widehat{y}_{n+1},x_{t_{n+1}})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{n},x_{t_{n}})+\textup{{TV}}( \widehat{y}_{n+1},x_{t_{n+1}})$ $$

By Corollary A.7,

$$ $\textup{{TV}}(\widehat{y}_{n+1},x_{t_{n+1}})\lesssim W_{2}(\widehat{y}_{n+1}^{ \prime},x_{t_{n+1}})\cdot d^{1/12}\cdot\sqrt{L}+\frac{\varepsilon_{\mathrm{sc} }}{\sqrt{L}d^{1/36}}+\sqrt{L}d^{17/36}h_{\mathrm{corr}}$ $$

Now, for $t_{n}\gtrsim 1/L$, by Lemma A.5,

$$ $W_{2}(\widehat{y}_{n+1}^{\prime},x_{t_{n+1}})\lesssim\frac{\varepsilon_{ \mathrm{sc}}}{L}+L^{3/2}\sqrt{d}h_{\mathrm{pred}}^{2}+L\sqrt{d}h_{\mathrm{pred }}^{3/2}\,.$ $$

Combining the above gives the first claim.

For the second claim, similar to above, we have

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{\delta})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{N_{0}+1},\widehat{y}_{N_{0}+1})+ \textup{{TV}}(\widehat{y}_{N_{0}+1},x_{\delta})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{N_{0}},x_{t_{N_{0}}})+\textup{{TV} }(\widehat{y}_{N_{0}+1},x_{\delta})$ $$

By Corollary A.7,

$$ $\textup{{TV}}(\widehat{y}_{N_{0}+1},x_{\delta})\leq W_{2}(\widehat{y}^{\prime} _{N_{0}+1},x_{\delta})\cdot d^{1/12}\cdot\sqrt{L}+\frac{\varepsilon_{\mathrm{ sc}}}{\sqrt{L}d^{1/36}}+\sqrt{L}d^{17/36}h_{\mathrm{corr}}$ $$

For $t_{N_{0}}\lesssim 1/L$, by Lemma A.5, noting that the number of predictor steps in this case is $O\left(\log\frac{h_{\mathrm{pred}}}{\delta}\right)$,

$$ $W_{2}(\widehat{y}^{\prime}_{N_{0}+1},x_{\delta})\lesssim\frac{\varepsilon_{ \mathrm{sc}}}{L}+\left(L^{3/2}\sqrt{d}h_{\mathrm{pred}}^{2}+L\sqrt{d}h_{ \mathrm{pred}}^{3/2}\right)\cdot\sqrt{\log\frac{h_{\mathrm{pred}}}{\delta}}$ $$

The second claim follows by combining the above.
∎

We recall the following lemma on the convergence of the OU process from []

###### Lemma A.9 (Lemma 13 of [ CCL + 23a ] ) .

Let $q_{t}$ denote the marginal law of the OU process started at $q_{0}=q$. Then, for all $T\gtrsim 1$,

$$ $\displaystyle\textup{{TV}}(q_{T},\mathcal{N}(0,I_{d}))\lesssim(\sqrt{d}+ \mathfrak{m}_{2})\exp(-T)$ $$

Finally, we prove our main theorem on the convergence of our sequential algorithm.

###### Theorem A.10 (Convergence bound for sequential algorithm) .

Suppose Assumptions 2.1-2.4 hold. If $\widehat{x}$ denotes the output of Algorithm 6, for $T=\Theta\left(\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{\varepsilon^{2}}
\right)\right),T_{\mathrm{corr}}=\Theta\left(\frac{1}{\sqrt{L}d^{1/18}}\right)$ and $\delta=\Theta\left(\frac{\varepsilon^{2}}{L^{2}(d\lor\mathfrak{m}_{2}^{2})}\right)$, if we set $h_{\mathrm{pred}}=\widetilde{\Theta}\left(\min\left(\frac{\varepsilon^{1/2}}{d
^{1/3}L^{3/2}},\frac{\varepsilon^{2/3}}{d^{5/12}L^{5/3}}\right)\cdot\frac{1}{
\log(\mathfrak{m}_{2})}\right)$, $h_{\mathrm{corr}}=\widetilde{\Theta}\left(\frac{\varepsilon}{d^{17/36}L^{3/2}
\log\mathfrak(m_{2})}\right)$, and if the score estimation satisfies $\varepsilon_{\mathrm{sc}}\leq\widetilde{O}\left(\frac{\varepsilon}{\sqrt{L}d^{
1/12}\log\mathfrak{m_{2}}}\right)$, we have that

$$ $\displaystyle\textup{{TV}}(\widehat{x},x_{0})\lesssim\varepsilon$ $$

with iteration complexity $\widetilde{\Theta}\left(\frac{L^{5/3}d^{5/12}}{\varepsilon}\cdot\log^{2}(
\mathfrak{m}_{2})\right)$

###### Proof.

We will let $t_{n}=T-n/L$.
First, note that by Lemma A.9

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{0},x_{t_{0}})\lesssim(\sqrt{d}+ \mathfrak{m}_{2})\exp(-T)$ $$

We divide our analysis into two steps. For the first $N_{0}=O(LT)$ steps, we iterate the first part of Lemma A.8 to obtain

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N_{0}},x_{t_{N_{0}}})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{0},x_{t_{0}})+O\left(L^{2}d^{7/12} h_{\mathrm{pred}}^{2}+L^{3/2}d^{7/12}h_{\mathrm{pred}}^{3/2}+\sqrt{L}d^{17/36} h_{\mathrm{corr}}+\frac{\varepsilon_{\mathrm{sc}}d^{1/12}}{\sqrt{L}}\right) \cdot N_{0}$ $\displaystyle\lesssim\left(\sqrt{d}+\mathfrak{m}_{2}\right)\exp(-T)+L^{3}d^{7/ 12}h_{\mathrm{pred}}^{2}T+L^{5/2}d^{7/12}h_{\mathrm{pred}}^{3/2}T+L^{3/2}d^{17 /36}h_{\mathrm{corr}}T+\varepsilon_{\mathrm{sc}}d^{1/12}T\sqrt{L}$ $$

Applying the second part of Lemma A.8 for the second stage of the algorithm, we have

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{\delta})$ $\displaystyle\lesssim\left(\sqrt{d}+\mathfrak{m}_{2}\right)\exp(-T)+\left(L^{3 }d^{7/12}h_{\mathrm{pred}}^{2}+L^{5/2}d^{7/12}h_{\mathrm{pred}}^{3/2}\right) \left(T+\sqrt{\log\frac{h_{\mathrm{pred}}}{\delta}}\right)$ $\displaystyle+L^{3/2}d^{17/36}h_{\mathrm{corr}}T+\varepsilon_{\mathrm{sc}}d^{1 /12}T\sqrt{L}$ $$

Setting $T=\Theta\left(\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{\varepsilon^{2}}
\right)\right)$, $h_{\mathrm{pred}}=\widetilde{\Theta}\left(\min\left(\frac{\varepsilon^{1/2}}{d
^{1/3}L^{3/2}},\frac{\varepsilon^{2/3}}{d^{5/12}L^{5/3}}\right)\frac{1}{\log(
\mathfrak{m}_{2})}\right)$, and $h_{\mathrm{corr}}=\widetilde{\Theta}\left(\frac{\varepsilon}{d^{17/36}L^{3/2}}
\cdot\frac{1}{\log(\mathfrak{m}_{2})}\right)$, if the score estimation error satisfies $\varepsilon_{\mathrm{sc}}\leq\widetilde{O}\left(\frac{\varepsilon}{\sqrt{L}d^{
1/12}\log(\mathfrak{m}_{2})}\right)$, with iteration complexity $\widetilde{\Theta}\left(\frac{L^{5/3}d^{5/12}\log^{2}\mathfrak{m}_{2}}{
\varepsilon}\right)$, we obtain $\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{\delta})\leq\varepsilon$.
∎

#### Input parameters:

- •
Start time $T$, End time $\delta$, Corrector steps time $T_{\mathrm{corr}}\lesssim 1/\sqrt{L}$, Number of predictor-corrector steps $N_{0}$, Predictor step size $h_{\mathrm{pred}}$, Corrector step size $h_{\mathrm{corr}}$, Score estimates $\widehat{s}_{t}$

- 1.
Draw $\widehat{x}_{0}\sim\mathcal{N}(0,I_{d})$.
- 2.
For $n=0,\dots,N_{0}-1$:
(a)
Starting from $\widehat{x}_{n}$, run Algorithm 4 with starting time $T-n/L$ using step sizes $h_{\mathrm{pred}}$ for all $N$ steps, with $N=\frac{1}{Lh_{\mathrm{pred}}}$, so that the total time is $1/L$. Let the result be $\widehat{x}_{n+1}^{\prime}$.
(b)
Starting from $\widehat{x}_{n+1}^{\prime}$, run Algorithm 5 for total time $T_{\mathrm{corr}}$ with step size $h_{\mathrm{corr}}$ and score estimate $\widehat{s}_{T-(n+1)/L}$ to obtain $\widehat{x}_{n+1}$.
- 3.
Starting from $\widehat{x}_{N_{0}}$, run Algorithm 4 with starting time $T-N_{0}/L$ using step sizes $h_{\mathrm{pred}}/2,h_{\mathrm{pred}}/4,h_{\mathrm{pred}}/8,\dots,\delta$ to obtain $\widehat{x}_{N_{0}+1}^{\prime}$.
- 4.
Starting from $\widehat{x}^{\prime}_{N_{0}+1}$, run Algorithm 5 for total time $T_{\mathrm{corr}}$ with step size $h_{\mathrm{corr}}$ and score estimate $\widehat{s}_{\delta}$ to obtain $\widehat{x}_{N_{0}+1}$.
- 5.
Return $\widehat{x}_{N_{0}+1}$.

###### Lemma A.8 (TV error after one round of predictor and corrector) .

###### Proof.

###### Lemma A.9 (Lemma 13 of [ CCL + 23a ] ) .

###### Theorem A.10 (Convergence bound for sequential algorithm) .

###### Proof.

## Appendix B Parallel algorithm

### B.1 Predictor step

In this section, we will apply a parallel version of randomized midpoint for the predictor step, where only $\widetilde{\Theta}(\log^{2}(\frac{Ld}{\varepsilon}))$ iteration complexity will be required to attain our desired error bound for one predictor step.

In each iteration $n$, we will first sample $R$ randomized midpoints that are in expectation evenly spaced with $\delta_{n}=\frac{h_{n}}{R_{n}}$ time intervals between consecutive midpoints. Next, in our step (c), we provide an initial estimate on the $x$ value of midpoints using our estimate of position $\widehat{x}_{n}$ at time $t_{n}$ provided by iteration $n-1$. This step is analogous to step (c) in Algorithm 4 for the sequential predictor step. Then, in step (d) we refine our initial estimates by using a discrete version of Picard iteration, where for round $k$, we compute a new estimate of $x_{t_{n}-\alpha_{i}h_{n}}$ based on the estimates of $x_{t_{n}-\alpha_{j}h_{n}}$ for $j\leq i$ in round $k-1$. Note that a trajectory $x(t)$ that starts from time $t_{0}$ and follows the true ODE is a fix point of operator $\tau$ that maps continuous function to continuous function, where

$$ $\tau(x)(t)=e^{t-t_{0}}x(t_{0})+\int_{t_{0}}^{t}e^{s-t_{0}}\cdot\nabla q_{s}(x( s))\mathrm{d}s.$ $$

By smoothness of the true ODE, the continuous Picard iteration converges exponentially to the true trajectory, and we will show that discretization error for Picard iteration is controlled.
After the refinements have sufficiently reduced the estimation error for our randomized midpoints, we make a final calculation, estimating the value of $x_{t_{n}+h_{n}}$ based on the estimated value at all the randomized midpoints.

Figure: Algorithm 7 PredictorStep (Parallel)

In our analysis, we follow the same notation as in Section A.1. We first establish a $\operatorname{poly}(L,d)$ bound on the initial estimation error incurred in step (c) of each iteration in Algorithm 7.

###### Claim B.1 .

Suppose $L\geq 1$. Assume $h_{n}\lesssim 1/L$. For any $n=0,\cdots,N-1$, suppose we draw $\widehat{x}_{n}$ from an arbitrary distribution $p_{n}$, then run step (a) - (e) in Algorithm 7. Then for any $i=1,\cdots,R_{n}$,

$$ $\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(0)}_{n,i}-x_{n}^{*}(\alpha_{i}h _{n})\right\|^{2}\lesssim h_{n}^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2}dh_{n}^{ 4}(L\lor\frac{1}{t_{n}-h_{n}})+L^{2}h^{2}\left\|\widehat{x}_{n}-x_{t_{n}} \right\|^{2},$ $$

where $x_{n}^{*}(t)$ is solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$.

###### Proof.

Notice that in step (c) of Algorithm 7, the initial estimate of the randomized midpoint is done with the exact same formula as in step (c) of Algorithm 6, except we calculate this initial estimate for $R_{n}$ different randomized midpoints. Notice also that the bound for discretization error in Lemma A.2 is not dependent on specific value of $\alpha$, as long as the randomed value is at most $1$. Hence we can use identical calculation to yield the claim.
∎

Next, we show how to drive the initialization error from Lemma B.1 down (exponentially) using the Picard iterations in step (d) of Algorithm 7.

###### Lemma B.2 .

Suppose $L\geq 1$. Assume $h_{n}\lesssim 1/L$. For all iterations $n\in\{0,\cdots,N-1\}$, suppose we draw $\widehat{x}_{n}$ from an arbitrary distribution $p_{n}$, then run step (a) - (e) in Algorithm 7. Then for all $k\in\{1,\cdots K_{n}\}$ and $i\in\{1,\cdots,R_{n}\}$,

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(k)}_{n,i}-x_{n}^{* }(\alpha_{i}h_{n})\right\|^{2}$ $\displaystyle\lesssim\left(8h_{n}^{2}L^{2}\right)^{k}\cdot\left(\frac{1}{R} \sum_{j=1}^{R}\left\|\widehat{x}^{(0)}_{n,j}-x_{n}^{*}(\alpha_{j}h_{n})\right \|^{2}\right)$ $\displaystyle\quad+h_{n}^{2}\left(\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh_ {n}^{2}}{R_{n}^{2}}(L\lor\frac{1}{t_{n}-h_{n}})+L^{2}\cdot\operatorname*{ \mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right),$ (21) $$

where $x_{n}^{*}(t)$ is solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$.

###### Proof.

Fixing iteration $n$, we will let $h:=h_{n}$, $R:=R_{n}$ and $\delta:=\delta_{n}$. The formula of $\widehat{x}^{(k)}_{n,i}$ and $x_{n}^{*}(\alpha_{i}h)$ has the same coefficient for $\hat{x}_{n}$, thus we can bound the difference as follows:

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(k)}_{n,i}-x_{n}^{* }(\alpha_{i}h)\right\|^{2}$ $\displaystyle\leq$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\sum_{j=1}^{i}\left(\int_{t_{n}- \min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\, \mathrm{d}s\cdot\widehat{s}_{t_{n}-\alpha_{j}h}(\widehat{x}^{(k-1)}_{n,j})- \int_{t_{n}-\alpha_{i}h}^{t_{n}}e^{s-(t_{n}-\alpha_{i}h)}\nabla\ln q_{s}(x_{n} ^{*}(t_{n}-s))\,\mathrm{d}s\right)\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2\operatorname*{\mathbb{E}}\left\|\sum_{j=1}^{i}\int_{t_{n}-\min (j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\,\mathrm{d }s\cdot\left(\widehat{s}_{t_{n}-\alpha_{j}h}(\widehat{x}^{(k-1)}_{n,j})-\nabla \ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))\right)\right\|^{2}$ (22) $\displaystyle\quad+2\operatorname*{\mathbb{E}}\left\|\sum_{j=1}^{i}\int_{t_{n} -\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\left( \nabla q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n}^{*} (t_{n}-s))\right)\,\mathrm{d}s\right\|^{2}.$ (23) $$

The first to second line is by definition and the second to third line is by Young’s inequality. Now, we will bound Equation 22 and Equation 23 separately. By 2.2 and 2.4,

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{s}_{t_{n}-\alpha_{j}h}( \widehat{x}^{(k-1)}_{n,j})-\nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j }h))\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2\operatorname*{\mathbb{E}}\left\|\widehat{s}_{t_{n}-\alpha_{j}h }(\widehat{x}^{(k-1)}_{n,j})-\nabla\ln q_{t_{n}-\alpha_{j}h}(\widehat{x}^{(k-1 )}_{n,j})\right\|^{2}+2\operatorname*{\mathbb{E}}\left\|\nabla\ln q_{t_{n}- \alpha_{j}h}(\widehat{x}^{(k-1)}_{n,j})-\nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^ {*}(\alpha_{j}h))\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2\varepsilon_{\mathrm{sc}}^{2}+2L^{2}\cdot\left\|\widehat{x}^{(k -1)}_{n,j}-x_{n}^{*}(\alpha_{j}h)\right\|^{2}.$ $$

The term in Equation 22 can now be bounded as follows

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\sum_{j=1}^{i}\int_{t_{n}-\min(j \delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\,\mathrm{d}s \cdot\left(\widehat{s}_{t_{n}-\alpha_{j}h}(\widehat{x}^{(k-1)}_{n,j})-\nabla \ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))\right)\right\|^{2}$ $\displaystyle\leq$ $\displaystyle R\cdot\sum_{j=1}^{i}\operatorname*{\mathbb{E}}\left\|\int_{t_{n} -\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\, \mathrm{d}s\cdot\left(\widehat{s}_{t_{n}-\alpha_{j}h}(\widehat{x}^{(k-1)}_{n,j })-\nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))\right)\right\|^{2}$ $\displaystyle\leq$ $\displaystyle R\cdot\delta^{2}\cdot e^{2\alpha_{i}h}\sum_{j=1}^{i} \operatorname*{\mathbb{E}}\left\|\widehat{s}_{t_{n}-\alpha_{j}h}(\widehat{x}^{ (k-1)}_{n,j})-\nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))\right\|^ {2}$ $\displaystyle\leq$ $\displaystyle 2R\cdot\delta^{2}\cdot e^{2\alpha_{i}h}\sum_{j=1}^{i}\left( \varepsilon_{\mathrm{sc}}^{2}+L^{2}\cdot\left\|\widehat{x}^{(k-1)}_{n,j}-x_{n} ^{*}(\alpha_{j}h)\right\|^{2}\right)$ $\displaystyle\leq$ $\displaystyle 2R^{2}\cdot\delta^{2}\cdot e^{2\alpha_{i}h}\frac{1}{R}\sum_{j=1} ^{R}\left(\varepsilon_{\mathrm{sc}}^{2}+L^{2}\cdot\left\|\widehat{x}^{(k-1)}_{ n,j}-x_{n}^{*}(\alpha_{j}h)\right\|^{2}\right)$ $\displaystyle\leq$ $\displaystyle 4h^{2}\varepsilon_{\mathrm{sc}}^{2}+4h^{2}L^{2}\cdot\frac{1}{R} \sum_{j=1}^{R}\cdot\left\|\widehat{x}^{(k-1)}_{n,j}-x_{n}^{*}(\alpha_{j}h) \right\|^{2}.$ $$

The first to second line is by inequality $(\sum_{i=1}^{n}a_{i})^{2}\leq n\sum_{i=1}^{n}a_{i}^{2}$, the second to third line is by the fact that $\int_{t_{n}-\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{
i}h)}\,\mathrm{d}s\leq\delta\cdot e^{\alpha_{i}h}$, the fifth to sixth line is by $R\delta=h$ and that $e^{2\alpha_{i}h}$ is at most $2$ when $h<1/4$.

Similarly, the term in Equation 23 can be bounded as follows

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\sum_{j=1}^{i}\int_{t_{n}-\min(j \delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\left(\nabla \ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n}^{*}(t_ {n}-s))\right)\,\mathrm{d}s\right\|^{2}$ $\displaystyle\leq$ $\displaystyle R\cdot\sum_{j=1}^{i}\operatorname*{\mathbb{E}}\left\|\int_{t_{n} -\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}e^{s-(t_{n}-\alpha_{i}h)}\left( \nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n}^ {*}(t_{n}-s))\right)\,\mathrm{d}s\right\|^{2}$ $\displaystyle\leq$ $\displaystyle R\cdot\delta\cdot e^{2\alpha_{i}h}\cdot\sum_{j=1}^{i}\int_{t_{n} -\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}\operatorname*{\mathbb{E}}\left \|\nabla\ln q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n }^{*}(t_{n}-s))\right\|^{2}\,\mathrm{d}s$ $$

Since $|x_{n}^{*}(\alpha_{j}h)-s|\leq\delta$,

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j}h}(x_{ n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 3\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j}h}(x _{n}^{*}(\alpha_{j}h))-\nabla q_{t_{n}-\alpha_{j}h}(x_{t_{n}-\alpha_{j}h}) \right\|^{2}+3\operatorname*{\mathbb{E}}\left\|\nabla q_{s}(x_{s})-\nabla q_{s }(x_{n}^{*}(t_{n}-s))\right\|^{2}$ $\displaystyle\quad+3\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j }h}(x_{t_{n}-\alpha_{j}h})-\nabla q_{s}(x_{s})\right\|^{2}$ $\displaystyle=$ $\displaystyle 3\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j}h}(x _{n}^{*}(\alpha_{j}h))-\nabla q_{t_{n}-\alpha_{j}h}(x_{t_{n}-\alpha_{j}h}) \right\|^{2}+3\operatorname*{\mathbb{E}}\left\|\nabla q_{s}(x_{s})-\nabla q_{s }(x_{n}^{*}(t_{n}-s))\right\|^{2}$ $\displaystyle\quad+3\operatorname*{\mathbb{E}}\left\|\int_{s}^{t_{n}-\alpha_{j }h}\partial_{u}\nabla\ln q_{u}(x_{u})du\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 3L^{2}\cdot\operatorname*{\mathbb{E}}\left\|x_{n}^{*}(\alpha_{j} h)-x_{t_{n}-\alpha_{j}h}\right\|^{2}+3L^{2}\cdot\operatorname*{\mathbb{E}} \left\|x_{n}^{*}(t_{n}-s)-x_{s}\right\|^{2}+3\operatorname*{\mathbb{E}}\left\| \int_{s}^{t_{n}-\alpha_{j}h}\partial_{u}\nabla\ln q_{u}(x_{u})du\right\|^{2}.$ $$

The second to third line is by Young’s inequality, and the fourth to fifth line is by 2.2.
By Lemma 3 in [],

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\int_{s}^{t_{n}-\alpha_{j}h} \partial_{u}\nabla\ln q_{u}(x_{u})du\right\|^{2}$ $\displaystyle\leq\delta\cdot\int_{s}^{t_{n}-\alpha_{j}h}\operatorname*{\mathbb {E}}\left\|\partial_{u}\nabla\ln q_{u}(x_{u})\right\|^{2}du$ $\displaystyle\leq\delta\cdot\int_{s}^{t_{n}-\alpha_{j}h}L^{2}d\max(L,\frac{1}{ u})du\leq L^{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h}).$ $$

Now, by Lemma A.1 and the fact that $h\lesssim\frac{1}{L}$,

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j}h}(x_{ n}^{*}(\alpha_{j}h))-\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\right\|^{2}$ $\displaystyle\leq 12L^{2}\exp(Lh)\operatorname*{\mathbb{E}}\left\|\widehat{x}_ {n}-x_{t_{n}}\right\|^{2}+3L^{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h})$ $\displaystyle\leq 36L^{2}\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t _{n}}\right\|^{2}+3L^{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h}).$ $$

We conclude that the term in Equation 23 can be bounded by

$$ $\displaystyle R\cdot\delta\cdot e^{2\alpha_{i}h}\cdot\sum_{j=1}^{i}\int_{t_{n} -\min(j\delta,\alpha_{i}h)}^{t_{n}-(j-1)\delta}36L^{2}\operatorname*{\mathbb{E }}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+3L^{2}d\delta^{2}(L\lor\frac{1} {t_{n}-h})\,\mathrm{d}s$ $\displaystyle\leq$ $\displaystyle 2R^{2}\delta^{2}\cdot\left(36L^{2}\operatorname*{\mathbb{E}} \left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+3L^{2}d\delta^{2}(L\lor\frac{1}{t _{n}-h})\right)$ $\displaystyle=$ $\displaystyle 4h^{2}L^{2}\left(18\operatorname*{\mathbb{E}}\left\|\widehat{x}_ {n}-x_{t_{n}}\right\|^{2}+\frac{3}{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h})\right).$ $$

Combining the bounds for Equation 22 and Equation 23, we get

$$ $\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(k)}_{n,i}-x_{n}(\alpha_{i}h) \right\|^{2}\\ \leq 8h^{2}L^{2}\cdot\left(\frac{1}{R}\sum_{j=1}^{R}\left\|\widehat{x}^{(k-1)} _{n,j}-x_{n}^{*}(\alpha_{j}h)\right\|^{2}+\frac{\varepsilon_{\mathrm{sc}}^{2}} {L^{2}}+\frac{3}{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h})+18\operatorname*{ \mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right).$ $$

Given sufficiently small constant $h$, $e^{\alpha_{i}h}\leq 2$. Moreover, by the definition of $\delta$, $R\delta=h$.
By unrolling the recursion, we get

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(k)}_{n,i}-x_{t_{n} -\alpha_{i}h}\right\|^{2}$ $\displaystyle\lesssim\left(8h^{2}L^{2}\right)^{k}\cdot\left(\frac{1}{R}\sum_{j =1}^{R}\left\|x^{(0)}_{n,j}-x_{n}^{*}(\alpha_{j}h)\right\|^{2}\right)$ $\displaystyle\quad+\frac{8h^{2}L^{2}}{1-8h^{2}L^{2}}\left(\frac{\varepsilon_{ \mathrm{sc}}^{2}}{L^{2}}+d\delta^{2}(L\lor\frac{1}{t_{n}-h})+\operatorname*{ \mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right)$ $\displaystyle\lesssim\left(8h^{2}L^{2}\right)^{k}\cdot\left(\frac{1}{R}\sum_{j =1}^{R}\left\|\widehat{x}^{(0)}_{n,j}-x_{n}^{*}(\alpha_{j}h)\right\|^{2}\right)$ $\displaystyle\quad+h^{2}\left(\varepsilon_{\mathrm{sc}}^{2}+L^{2}d\delta^{2}(L \lor\frac{1}{t_{n}-h})+L^{2}\cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_ {n}-x_{t_{n}}\right\|^{2}\right)\,.\qed$ $$

As a consequence of Lemma B.2, if we take the number $K_{n}$ of Picard iterations sufficiently large, the error incurred in our estimate for $x_{t_{n}-\alpha_{i}h_{n}}$ is dominated by the terms in the second line of Eq. (21).

###### Corollary B.3 .

Assume $L\geq 1$. For all $n\in\{0,\cdots,N-1\}$, suppose we draw $\widehat{x}_{n}$ from an arbitrary distribution $p_{n}$, then run step (a) - (e) in Algorithm 7.
In addition, suppose $h_{n}<\frac{1}{3L}$ and $K_{n}\gtrsim\log(R_{n})$. Then for any $i\in\{1,\cdots,R_{n}\}$,

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(K_{n})}_{n,i}-x_{n }^{*}(\alpha_{i}h_{n})\right\|^{2}$ $\displaystyle\lesssim h_{n}^{2}\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh_{n} ^{4}}{R_{n}^{2}}\left(L\lor\frac{1}{t_{n}-h_{n}}\right)+L^{2}h_{n}^{2}\cdot \operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2},$ $$

where $x_{n}^{*}(t)$ is solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$.

###### Proof.

Fixing iteration $n$, we will let $h:=h_{n}$, $R:=R_{n}$, $\delta:=\delta_{n}$ and $K:=K_{n}$.
Notice that when $K\geq\frac{2}{\log\frac{1}{8h^{2}L^{2}}}\cdot\log(R)$, $\left(8h^{2}L^{2}\right)^{K}$ is at most $\frac{1}{R^{2}}$. Now by plugging B.1 into Lemma B.2, we get

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\widehat{x}^{(K)}_{n,i}-x_{n}^{* }(\alpha_{i}h)\right\|^{2}$ $\displaystyle\lesssim\left(8h^{2}L^{2}\right)^{K}\cdot h^{2}\cdot\left( \varepsilon_{\mathrm{sc}}^{2}+L^{2}dh^{2}(L\lor\frac{1}{t_{n}-h})+L^{2}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right)$ $\displaystyle\quad+h^{2}\left(\varepsilon_{\mathrm{sc}}^{2}+L^{2}d\delta^{2}(L \lor\frac{1}{t_{n}-h})+L^{2}\cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_ {n}-x_{t_{n}}\right\|^{2}\right)$ $\displaystyle\lesssim\left((8h^{2}L^{2})^{K}+1\right)h^{2}\cdot\left( \varepsilon_{\mathrm{sc}}^{2}+L^{2}\cdot\operatorname*{\mathbb{E}}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right)$ $\displaystyle\qquad+((8h^{2}L^{2})^{K}+\frac{1}{R^{2}})\cdot L^{2}dh^{4}\left( L\lor\frac{1}{t_{n}-h}\right)$ $\displaystyle\lesssim h^{2}\varepsilon_{\mathrm{sc}}^{2}+L^{2}h^{2}\cdot \operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+\frac{L ^{2}dh^{4}}{R^{2}}\left(L\lor\frac{1}{t_{n}-h}\right)$ $$

The first to second inequality by rearrangement of terms, while the second to third inequality is by the fact that the terms $\frac{1}{R^{2}}$ dominates $\left(8h^{2}L^{2}\right)^{K}$.
∎

We can now prove the parallel analogue of Lemma A.3 and Lemma A.4. Note that the bounds in Lemma B.4 and Lemma B.5 are identical to the bounds in Lemma A.3 and Lemma A.4, except from an additional $\frac{1}{R_{n}^{2}}$ factor for the middle term. This additional factor stems from using $R_{n}$ midpoints in each iteration $n$ (compared to using one midpoint each iteration in Algorithm 6).

###### Lemma B.4 (Parallel Predictor Bias) .

Assume $L\geq 1$. For all $n\in\{0,\cdots,N-1\}$, suppose we draw $\widehat{x}_{n}$ from an arbitrary distribution $p_{n}$, then run step (a) - (e) in Algorithm 7.
In addition, suppose $h_{n}<\frac{1}{3L}$ and $K_{n}\gtrsim\log(R_{n})$. Then we have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\operatorname*{\mathbb{E}}_{\alpha} \widehat{x}_{n+1}-x_{n}^{*}(h_{n})\|^{2}\lesssim h_{n}^{2}\cdot\varepsilon_{ \mathrm{sc}}^{2}+\frac{L^{4}h_{n}^{6}d}{R_{n}^{2}}\cdot(L\lor\frac{1}{t_{n}-h_ {n}})))+L^{4}h_{n}^{4}\cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_ {t_{n}}\right\|^{2},$ $$

where $x_{n}^{*}(t)$ is solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$.

###### Proof.

Fixing iteration $n$, we will let $h:=h_{n}$, $R:=R_{n}$, $\delta:=\delta_{n}$ and $K:=K_{n}$. We have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|\operatorname*{\mathbb{E}}_{\alpha} \widehat{x}_{n+1}-x_{n}^{*}(h)\|^{2}$ $\displaystyle\leq$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\operatorname*{\mathbb{E}}_{ \alpha}\left[\delta\cdot\sum_{i=1}^{R}e^{h-\alpha_{i}h}\cdot\widehat{s}_{t_{n} -\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})\right]-\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n }-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2\operatorname*{\mathbb{E}}\left\|\operatorname*{\mathbb{E}}_{ \alpha}\left[\sum_{i=1}^{R}\delta e^{h-\alpha_{i}h}\cdot\left(\widehat{s}_{t_{ n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})-\nabla\ln q_{t_{n}-\alpha_{i}h}(x_{t_ {n}-\alpha_{i}h})\right)\right]\right\|^{2}$ (24) $\displaystyle+2\operatorname*{\mathbb{E}}\left\|\operatorname*{\mathbb{E}}_{ \alpha}\sum_{i=1}^{R}\delta e^{h-\alpha_{i}h}\cdot\nabla\ln q_{t_{n}-\alpha_{i }h}(x_{t_{n}-\alpha_{i}h})-\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n}-h)}\nabla\ln q_{s }(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\right\|^{2}\,.$ (25) $$

Since $\alpha_{i}$ is drawn uniformly from $[(i-1)\delta,i\delta]$,

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\left[\delta e^{h-\alpha_{i}h} \cdot\nabla\ln q_{t_{n}-\alpha_{i}h}(x_{t_{n}-\alpha_{i}h})\right]=\int_{t_{n} -i\delta}^{t_{n}-(i-1)\delta}e^{s-(t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s) )\,\mathrm{d}s,$ $$

and the term in Equation 25 is equal to 0 0. Now we need to bound Equation 24.

$$ $\displaystyle\left\|\widehat{s}_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})- \nabla\ln q_{t_{n}-\alpha_{i}h}(x_{t_{n}-\alpha_{i}h})\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2\left\|\widehat{s}_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i}) -\nabla\ln q_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})\right\|^{2}+2\left\| \nabla\ln q_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})-\nabla\ln q_{t_{n}- \alpha_{i}h}(x_{t_{n}-\alpha_{i}h})\right\|^{2}$ $\displaystyle\lesssim$ $\displaystyle\varepsilon_{\mathrm{sc}}^{2}+L^{2}\cdot\left\|x^{(K)}_{n,j}-x_{n }^{*}(\alpha_{j}h)\right\|^{2}$ $\displaystyle\lesssim$ $\displaystyle\varepsilon_{\mathrm{sc}}^{2}+L^{2}\cdot\left(h^{2}\varepsilon_{ \mathrm{sc}}^{2}+\frac{L^{2}dh^{4}}{R^{2}}\left(L\lor\frac{1}{t_{n}-h}\right)+ L^{2}h^{2}\cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}} \right\|^{2}\right).$ $$

The first to second line by inequality $(\sum_{i=1}^{n}a_{i})^{2}\leq n\sum_{i=1}^{n}a_{i}^{2}$, the second to third line is by 2.4 and 2.2, and the third to fourth line is by Corollary B.3 ($K\gtrsim\log(R)$, which satisfies the condition in Corollary B.3). Hence

$$ $\displaystyle 2\operatorname*{\mathbb{E}}\left\|\operatorname*{\mathbb{E}}_{ \alpha}\left[\sum_{i=1}^{R}\delta e^{h-\alpha_{i}h}\cdot\left(\widehat{s}_{t_{ n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})-\nabla\ln q_{t_{n}-\alpha_{i}h}(x_{t_ {n}-\alpha_{i}h})\right)\right]\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 2R(2\delta)^{2}\sum_{i=1}^{R}\operatorname*{\mathbb{E}}_{\alpha} \left\|\widehat{s}_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i})-\nabla\ln q_{t _{n}-\alpha_{i}h}(x_{t_{n}-\alpha_{i}h})\right\|^{2}$ $\displaystyle\lesssim$ $\displaystyle 8\delta^{2}R^{2}\cdot\left(\varepsilon_{\mathrm{sc}}^{2}+L^{2} \cdot\left(h^{2}\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh^{4}}{R^{2}}\left(L \lor\frac{1}{t_{n}-h}\right)+L^{2}h^{2}\cdot\operatorname*{\mathbb{E}}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right)\right)$ $\displaystyle\lesssim$ $\displaystyle h^{2}\cdot\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{4}dh^{6}}{R^{2} }(L\lor\frac{1}{t_{n}-h})+L^{4}h^{4}\cdot\operatorname*{\mathbb{E}}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}.$ $$

The first to second step is by inequality $(\sum_{i=1}^{n}a_{i})^{2}\leq n\sum_{i=1}^{n}a_{i}^{2}$ and Young’s inequality, the second to third line is by plugging in our previous calculation, and the third to forth line is by $h=\delta R$ and that $h\lesssim 1/L$.
∎

###### Lemma B.5 (Parallel Predictor Variance) .

Assume $L\geq 1$. For all $n\in\{0,\cdots,N-1\}$, suppose we draw $\widehat{x}_{n}$ from an arbitrary distribution $p_{n}$, then run step (a) - (e) in Algorithm 7.
In addition, suppose $h_{n}<\frac{1}{3L}$ and $K_{n}\gtrsim\log(R_{n})$. Then we have

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{n+1}-x_{n}^{*}( h_{n})\|^{2}\lesssim h_{n}^{2}\cdot\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh _{n}^{4}}{R_{n}^{2}}\left(L\lor\frac{1}{t_{n}-h_{n}}\right)+L^{2}h_{n}^{2} \cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2},$ $$

where $x_{n}^{*}(t)$ is solution of the true ODE starting at $\widehat{x}_{n}$ at time $t_{n}$ and running until time $t_{n}-t$.

###### Proof.

Fixing iteration $n$, we will let $h:=h_{n}$, $R:=R_{n}$, $\delta:=\delta_{n}$ and $K:=K_{n}$. We will separate $\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{n+1}-x_{n}^{*}(h)\|^{2}$ into several terms and bound each term separately.

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{n+1}-x_{n}^{*}( h)\|^{2}$ $\displaystyle\leq$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\left\|\delta\cdot\sum_{i=1}^{ R}e^{h-\alpha_{i}h}\cdot\widehat{s}_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{n,i })-\int_{t_{n}-h}^{t_{n}}e^{s-(t_{n}-h)}\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\, \mathrm{d}s\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 3\operatorname*{\mathbb{E}}_{\alpha}\left\|\sum_{i=1}^{R}\delta e ^{h-\alpha_{i}h}\cdot\left(\widehat{s}_{t_{n}-\alpha_{i}h}(\widehat{x}^{(K)}_{ n,i})-\nabla\ln q_{t_{n}-\alpha_{i}h}(x_{t_{n}-\alpha_{i}h})\right)\right\|^{2}$ (26) $\displaystyle+3\operatorname*{\mathbb{E}}\left\|\sum_{i=1}^{R}\int_{t_{n}-i \delta}^{t_{n}-(i-1)\delta}e^{h-\alpha_{i}h}\cdot\left(\nabla\ln q_{t_{n}- \alpha_{i}h}(x_{t_{n}-\alpha_{i}h})-\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\right) \,\mathrm{d}s\right\|^{2}$ (27) $\displaystyle+3\operatorname*{\mathbb{E}}\left\|\sum_{i=1}^{R}\int_{t_{n}-i \delta}^{t_{n}-(i-1)\delta}\left(e^{h-\alpha_{i}h}-e^{s-(t_{n}-h)}\right)\cdot \nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\right\|^{2}.$ (28) $$

Equation 26 is identical to Equation 24 in Lemma B.4, and can be bounded by

$$ $\displaystyle\lx@cref{creftypecap~refnum}{eq:parallel_predictor_var_term_1} \lesssim h^{2}\cdot\varepsilon_{\mathrm{sc}}^{2}+L^{4}h^{4}\cdot\operatorname* {\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+\frac{L^{4}dh^{6}}{R^ {2}}(L\lor\frac{1}{t_{n}-h}))).$ $$

Next we will bound Equation 27. By Lemma D.1 and Lemma A.1,

$$ $\operatorname*{\mathbb{E}}\left\|\nabla q_{t_{n}-\alpha_{j}h}(x_{n}^{*}(\alpha _{j}h))-\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\right\|^{2}\lesssim L^{2}d\delta^{ 2}(L\lor\frac{1}{t_{n}-h})+L^{2}\cdot\operatorname*{\mathbb{E}}\left\|\widehat {x}_{n}-x_{t_{n}}\right\|^{2},$ $$

hence Equation 27 can be bounded with similar calculations as for Equation 23 in Lemma B.2, by the following term:

$$ $12R^{2}\delta^{2}\cdot\left(L^{2}d\delta^{2}(L\lor\frac{1}{t_{n}-h})+O(L^{2}) \cdot\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2} \right)\lesssim\frac{L^{2}dh^{4}}{R^{2}}(L\lor\frac{1}{t_{n}-h})+L^{2}h^{2} \operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}.$ $$

Finally we will bound Equation 28
Since both $\alpha_{i}$ and $s$ belong to the range $[(i-1)\delta,i\delta]$,

$$ $e^{h-\alpha_{i}h}-e^{s-(t_{n}-h)}\leq e^{h}(e^{-(i-1)\delta}-e^{-i\delta})\leq e ^{h}\cdot\delta\leq 2\delta.$ $$

Moreover, by the fact that $\operatorname*{\mathbb{E}}\left\|\nabla\ln q_{s}(x_{t_{n}})\right\|^{2}\leq Ld$ (by integration by parts), Lemma D.1, Lemma A.1 and the fact that $L\delta=o(1)$, we have

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\nabla\ln q_{s}(x_{n}^{*}(t_{n}- s))\right\|^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}\left\|\nabla\ln q_{s}(x_{t_{n} })\right\|^{2}+\operatorname*{\mathbb{E}}\left\|\nabla\ln q_{s}(x_{n}^{*}(t_{n }-s))-\nabla\ln q_{s}(x_{t_{n}})\right\|^{2}$ $\displaystyle\lesssim Ld+L^{2}\exp(L\delta)\left\|\widehat{x}_{n}-x_{t_{n}} \right\|^{2}\lesssim Ld+L^{2}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}.$ $$

Hence Equation 28 can be bounded by

$$ $\displaystyle 3\operatorname*{\mathbb{E}}\left\|\sum_{i=1}^{R}\int_{t_{n}-i \delta}^{t_{n}-(i-1)\delta}\left(e^{h-\alpha_{i}h}-e^{s-(t_{n}-h)}\right)\cdot \nabla\ln q_{s}(x_{n}^{*}(t_{n}-s))\,\mathrm{d}s\right\|^{2}$ $\displaystyle\leq$ $\displaystyle 3R\delta\sum_{i=1}^{R}\int_{t_{n}-i\delta}^{t_{n}-(i-1)\delta} \operatorname*{\mathbb{E}}\left\|2\delta\nabla\ln q_{s}(x_{n}^{*}(t_{n}-s)) \right\|^{2}$ $\displaystyle\lesssim$ $\displaystyle 3R\delta\cdot 4\delta^{2}\sum_{i=1}^{R}\int_{t_{n}-i\delta}^{t_{ n}-(i-1)\delta}\left(Ld+L^{2}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}\right)$ $\displaystyle\lesssim$ $\displaystyle R^{2}\delta^{4}\left(Ld+L^{2}\left\|\widehat{x}_{n}-x_{t_{n}} \right\|^{2}\right)=\frac{Ldh^{4}}{R^{2}}+\frac{L^{2}h^{4}}{R^{2}}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}.$ $$

By adding together Equation 26, Equation 27 and Equation 28, and combining terms, we conclude that

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\|\widehat{x}_{n+1}-x_{n}^{*}( h)\|^{2}$ $\displaystyle\lesssim h^{2}\cdot\varepsilon_{\mathrm{sc}}^{2}+L^{4}h^{4}\cdot \operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+\frac{L ^{4}dh^{6}}{R^{2}}(L\lor\frac{1}{t_{n}-h})))$ $\displaystyle\quad+\frac{L^{2}dh^{4}}{R^{2}}(L\lor\frac{1}{t_{n}-h})+L^{2}h^{2 }\operatorname*{\mathbb{E}}\left\|\widehat{x}_{n}-x_{t_{n}}\right\|^{2}+\frac{ Ldh^{4}}{R^{2}}+\frac{L^{2}h^{4}}{R^{2}}\left\|\widehat{x}_{n}-x_{t_{n}}\right \|^{2}$ $\displaystyle\lesssim h^{2}\cdot\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh^{4 }}{R^{2}}(L\lor\frac{1}{t_{n}-h})+L^{2}h^{2}\operatorname*{\mathbb{E}}\left\| \widehat{x}_{n}-x_{t_{n}}\right\|^{2}\,.\qed$ $$

We can now prove our main guarantee for the parallel predictor step, which states that with logarithmically many parallel rounds and $\widetilde{O}(\sqrt{d})$ score estimate queries over a short time interval (of length $O(1/L)$) of the reverse process, the algorithm does not drift too far from the true ODE. Our proof follows a similar flow as Lemma A.5.

Note that due to the existence of $R_{n}$ midpoints in each step, we can set $h_{n}=\Theta(1)$. We will set $h_{n}$ to be $\Theta(\frac{1}{L})$, unless $t_{n}$ is close to the end time $\delta$ (see Algorithm 9 for the global algorithm and timeline). If $t_{n}$ is close to $\delta$, we will repeated half $h_{n}$ as $n$ increases, until we reach the end time.

###### Theorem B.6 .

Assume $L\geq 1$. Let $\beta\geq 1$ be an adjustable parameter. When we set $h_{n}=\min\{\frac{1}{4L},t_{n}/2,t_{n}-\delta\}$, $K\gtrsim\log(\frac{\beta\sqrt{d}}{\varepsilon})$ and $R_{n}=h_{n}\cdot\beta\cdot\frac{L\sqrt{d}}{\varepsilon}=O\left(\frac{\beta
\sqrt{d}}{\varepsilon}\right)$, the Wasserstein distance between the true ODE process and the process in Algorithm 7, both starting from $\hat{x}_{0}\sim q_{t_{0}}$ and run for total time $T_{N}\lesssim\frac{1}{L}$ is bounded by

$$ $W_{2}(\widehat{x}_{N},x_{T_{N}})\lesssim\frac{\varepsilon_{\mathrm{sc}}}{L}+ \frac{\varepsilon}{\beta\sqrt{L}}.$ $$

###### Proof.

To avoid confusion, we will reserve $\widehat{x}_{n}$ as the result of running Algorithm 7 for $n$ steps, starting at $\widehat{x}_{0}\sim q_{0}$. We will use $\widehat{y}_{n}$ to denote the result from running the true ODE process, starting at $\widehat{x}_{n-1}$. Then by an identical calculation as in Lemma A.5,

$$ $\displaystyle\operatorname*{\mathbb{E}}_{\alpha}\left\|x_{t_{N}}-\widehat{x}_{ N}\right\|^{2}\leq\exp\left(O(Lh_{N-1})\right)\left\|x_{t_{N-1}}-\widehat{x}_{ N-1}\right\|^{2}+\frac{2}{Lh_{N-1}}\left\|\operatorname*{\mathbb{E}}_{\alpha} \widehat{x}_{N}-y_{N}\right\|^{2}+\operatorname*{\mathbb{E}}_{\alpha}\left\| \widehat{x}_{N}-y_{N}\right\|^{2}.$ $$

Now we do a similar calculation as in Lemma A.5, but utilizes the bias and variance of one step in the parallel algorithm instead of sequential. Taking the expectation wrt $\widehat{x}_{0}\sim q_{t_{0}}$, by Lemmas B.4 and B.5,

$$ $\displaystyle\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}+\frac{2}{Lh_{N-1}}\operatorname*{\mathbb{E}}\| \operatorname*{\mathbb{E}}_{\alpha}\widehat{x}_{N}-y_{N}\|^{2}+\operatorname*{ \mathbb{E}}\|\widehat{x}_{N}-y_{N}\|^{2}$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}$ $\displaystyle\qquad+O\left(\frac{1}{Lh_{N-1}}\left(h_{N-1}^{2}\varepsilon_{ \mathrm{sc}}^{2}+\frac{L^{4}dh_{N-1}^{6}}{R_{n}^{2}}\left(L\lor\frac{1}{t_{N-1 }}\right)+L^{4}h_{N-1}^{4}\operatorname*{\mathbb{E}}\|x_{t_{N-1}}-\widehat{x}_ {N-1}\|^{2}\right)\right)$ $\displaystyle\qquad+O\left(h_{N-1}^{2}\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2 }dh_{N-1}^{4}}{R_{n}^{2}}\left(L\lor\frac{1}{t_{N-1}}\right)+L^{2}h_{N-1}^{2} \operatorname*{\mathbb{E}}\|x_{t_{N-1}}-\widehat{x}_{N-1}\|^{2}\right)$ $\displaystyle\leq\exp\left(O(Lh_{N-1})\right)\operatorname*{\mathbb{E}}\|x_{t_ {N-1}}-\widehat{x}_{N-1}\|^{2}$ $\displaystyle\qquad+O\left(\frac{h_{N-1}\varepsilon_{\mathrm{sc}}^{2}}{L}+h_{N -1}^{2}\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{3}dh_{N-1}^{5}+L^{2}dh_{N-1}^{4} }{R_{n}^{2}}\left(L\lor\frac{1}{t_{N-1}}\right)\right)$ $$

Since $h_{N-1}<\frac{1}{L}$, the term $L^{3}dh_{N-1}^{5}$ is dominated by the term $L^{2}dh_{N-1}^{4}$ and the term $h_{N-1}^{2}\varepsilon_{\mathrm{sc}}^{2}$ os dominated by $\frac{h_{N-1}\varepsilon_{\mathrm{sc}}^{2}}{L}$. By induction, noting that $x_{t_{0}}=\widehat{x}_{0}$, we have

$$ $\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2}\lesssim\sum_{n=0}^ {N-1}\left(\frac{h_{n}}{L}\varepsilon_{\mathrm{sc}}^{2}+\frac{L^{2}dh_{n}^{4}} {R_{n}^{2}}\cdot\left(L\lor\frac{1}{t_{n}}\right)\right)\cdot\exp\left(O\left( L\sum_{i=n+1}^{N-1}h_{i}\right)\right)\,.$ $$

Since $\sum_{i=n+1}^{N-1}h_{i}\leq\frac{1}{L}$, $\exp\left(O\left(L\sum_{i=n+1}^{N-1}h_{i}\right)\right)$ is a constant. Moreover, by our choice of $h_{n}$ , it is always the case that $h_{n}\leq\frac{t_{n}}{2}\leq t_{n}-h_{n}$ and that $h_{n}\leq\frac{1}{L}$. Therefore

$$ $\displaystyle\frac{L^{2}dh_{n}^{4}}{R_{n}^{2}}\left(L\lor\frac{1}{t_{n}-h_{n}} \right)\leq\frac{L^{2}dh_{n}^{3}}{R_{n}^{2}}\leq\frac{L^{2}dh_{n}}{\frac{\beta ^{2}L^{2}d}{\varepsilon^{2}}}=\frac{h_{n}\varepsilon^{2}}{\beta^{2}},$ $$

and thus

$$ $\displaystyle\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2} \lesssim\sum_{n=0}^{N-1}\frac{h_{n}}{L}\varepsilon_{\mathrm{sc}}^{2}+\frac{h_{ n}\varepsilon^{2}}{\beta^{2}}\lesssim\frac{\varepsilon_{\mathrm{sc}}^{2}}{L^{2 }}+\frac{\varepsilon^{2}}{L\beta^{2}}.$ $$

We conclude that

$$ $\displaystyle W_{2}(x_{t_{N}},\widehat{x}_{N})$ $\displaystyle=\sqrt{\operatorname*{\mathbb{E}}\|x_{t_{N}}-\widehat{x}_{N}\|^{2 }}\lesssim\frac{\varepsilon_{\mathrm{sc}}}{L}+\frac{\varepsilon}{\beta\sqrt{L} }.\qed$ $$

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Starting time $t_{0}$, Number of steps $N$, Step size $\{h_{n}\}_{n=0}^{N-1}$, Number of midpoint estimates $\{R_{n}\}_{n=0}^{N-1}$, Number of parallel iteration $\{K_{n}\}_{n=0}^{N-1}$, Score estimates $\widehat{s}_{t}$
- •
For all $n=0,\cdots,N-1$:
let $\delta_{n}=\frac{h_{n}}{R_{n}}$

- 1.
For $n=0,\dots,N-1$:
(a)
Let $t_{n}=t_{0}-\sum_{w=0}^{n-1}h_{w}$
(b)
Randomly sample $\alpha_{i}$ uniformly from $[(i-1)/R_{n},i/R_{n}]$ for all $i\in\{1,\cdots,R_{n}\}$
(c)
For $i=1,\cdots,R_{n}$ in parallel: Let $\widehat{x}^{(0)}_{n,i}=e^{\alpha_{i}h_{n}}\widehat{x}_{n}+\left(e^{\alpha_{i}
h_{n}}-1\right)\cdot\widehat{s}_{t_{n}}(\widehat{x}_{n})$
(d)
For $k=1,\cdots,K_{n}$:
For $i=1,\cdots,R_{n}$ in parallel:
$\widehat{x}^{(k)}_{n,i}:=e^{\alpha_{i}h_{n}}\widehat{x}_{n}+\sum_{j=1}^{i}
\left(e^{\alpha_{i}h_{n}-(j-1)\delta_{n}}-\max(e^{\alpha_{i}h_{n}-j\delta_{n}}
,1)\right)\cdot\widehat{s}_{t_{n}-\alpha_{j}h_{n}}(\widehat{x}^{(k-1)}_{n,j})$
(e)
$\widehat{x}_{n+1}=e^{h_{n}}\widehat{x}_{n}+\delta_{n}\cdot\sum_{i=1}^{R}e^{h_{
n}-\alpha_{i}h_{n}}\widehat{s}_{t_{n}-\alpha_{i}h_{n}}(\widehat{x}^{(K_{n})}_{
n,i})$
- 2.
Let $t_{N}=t_{0}-\sum_{n=0}^{N-1}h_{n}$
- 3.
Return $\widehat{x}_{N},t_{N}$.

###### Claim B.1 .

###### Proof.

###### Lemma B.2 .

###### Proof.

###### Corollary B.3 .

###### Proof.

###### Lemma B.4 (Parallel Predictor Bias) .

###### Proof.

###### Lemma B.5 (Parallel Predictor Variance) .

###### Proof.

###### Theorem B.6 .

###### Proof.

### B.2 Corrector step

In this step we will be using the parallel algorithm in [] to estimate the underdamped Langevin diffusion process. Since we will be fixing the score function in time, we will use $\nabla\ln q$ to denote the true score function for the diffusion process, and $\widehat{s}$ to denote the estimated score function. We will choose the friction parameter to be $\gamma\asymp\sqrt{L}$.

Figure: Algorithm 8 Corrector Step (Parallel) []

Let $T_{N}$ denote the total time the parallel corrector step is run (namely, $T_{N}=nh$). Consider two continuous underdamped Langevin diffusion processes $u^{*}(t)=(x^{*}(t),v^{*}(t))$ and $u_{t_{0}+t}=(x_{t_{0}+t},v_{t_{0}+t})$ with coupled brownian motions. The first one start from position $x^{*}(0)=\widehat{x}_{0}\sim p$ and the second one start from position $x_{t_{0}}\sim q$. Both processes start with velocity $v^{*}(0)=v_{t_{0}}\sim\mathcal{N}(0,I_{d})$.
We will bound both the distance measure between $x^{*}(t)$ and the true sample $x_{t_{0}+t}$, and the distance measure between $x^{*}(t_{N})$ and outputs of Algorithm 8.
First, [] gives the following bound on the total variation error between $x^{*}(T_{N})$ and $x_{t_{N}}$.

###### Lemma B.7 ( [ CCL + 23a ] , Lemma 9) .

If $h\lesssim\frac{1}{\sqrt{L}}$, then

$$ $\textup{{TV}}(x^{*}(T_{N}),x_{t_{N}})\lesssim\frac{W_{2}(p,q)}{L^{1/4}T_{N}^{3 /2}}\,.$ $$

Next, [] bounds the discretization error in Algorithm 8 in terms of quantities that relates to the supremum of $\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x^{*}(t))\right\|^{2}$ and $\operatorname*{\mathbb{E}}\left\|v^{*}(t)\right\|^{2}$ where $t\in[0,T_{N}]$.

###### Lemma B.8 ( [ ACV24 ] , Theorem 20, Implicit) .

Assume $L\geq 1$. In Algorithm 8, assume $K\gtrsim\log(d)$ (for sufficiently large constant), $K\lesssim 4\log R$ and $h\lesssim\frac{1}{\sqrt{L}}$. Then

$$ $\displaystyle\textup{{KL}}(\widehat{x}_{N},x^{*}(T_{N}))\lesssim\frac{T_{N}}{ \sqrt{L}}\cdot\left(\varepsilon_{\mathrm{sc}}^{2}+L^{2}(\frac{\gamma dh^{3}}{R ^{4}}+\frac{h^{2}}{R^{2}}{\cal P}+\frac{h^{4}}{R^{4}}{\cal Q})\right),$ $$

where ${\cal P}=\sup_{t\in[0,T_{N}]}\operatorname*{\mathbb{E}}[\left\|v^{*}(t)\right
\|^{2}]$ and ${\cal Q}=\sup_{t\in[0,T_{N}]}\operatorname*{\mathbb{E}}[\left\|\nabla\ln q(x^{
*}(t))\right\|^{2}]$.

To reason about the value of ${\cal P}$ and ${\cal Q}$, we will use the following lemma in [].

###### Lemma B.9 ( [ CCL + 23a ] , Lemma 10) .

For any $t\lesssim\frac{1}{\sqrt{L}}$,

$$ $\operatorname*{\mathbb{E}}\left\|u^{*}(t)-u_{t_{0}+t}\right\|^{2}\lesssim W_{2 }^{2}(p,q).$ $$

###### Lemma B.10 .

Assume $L\geq 1$. For any $T_{N}\lesssim\frac{1}{\sqrt{L}}$,

$$ $\sup_{t\in[0,T_{N}]}\operatorname*{\mathbb{E}}[\left\|\nabla\ln q(x^{*}(t)) \right\|^{2}]\lesssim L^{2}W_{2}^{2}(p,q)+Ld$ $$

and

$$ $\sup_{t\in[0,T_{N}]}\operatorname*{\mathbb{E}}\left\|v^{*}(t)\right\|^{2} \lesssim W_{2}^{2}(p,q)+d.$ $$

###### Proof.

Note that $(q,\mathcal{N}(0,I_{d}))$ is a stationary distribution of the underdamped Langevin diffusion process, hence $x_{t}\sim q$ and $v_{t}\sim\mathcal{N}(0,I_{d})$. Hence $\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x_{t})\right\|^{2}\leq Ld$ by integration by parts. Similarly, $\operatorname*{\mathbb{E}}\left\|v_{t}\right\|^{2}=\operatorname*{\mathbb{E}}[
\left\|\mathcal{N}(0,I_{d})\right\|^{2}]\lesssim d$. Since $T_{N}\lesssim\frac{1}{\sqrt{L}}$, for any $t\in[0,T_{N}]$, we can now bound $\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x^{*}(t))\right\|^{2}$ and $\operatorname*{\mathbb{E}}\left\|v^{*}(t)\right\|^{2}$ by Lemma B.9 as follows:

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x^{*}(t))\right\|^{2}$ $\displaystyle\leq 2\cdot\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x_{t}) \right\|^{2}+2\cdot\operatorname*{\mathbb{E}}\left\|\nabla\ln q(x^{*}(t))- \nabla\ln q(x_{t})\right\|^{2}$ $\displaystyle\lesssim Ld+L^{2}\operatorname*{\mathbb{E}}\left\|x^{*}(t)-x_{t} \right\|^{2}\lesssim Ld+L^{2}\operatorname*{\mathbb{E}}\left\|u^{*}(t)-u_{t} \right\|^{2}\lesssim Ld+L^{2}W_{2}^{2}(p,q),$ $$

and

$$ $\displaystyle\operatorname*{\mathbb{E}}\left\|v^{*}(t)\right\|^{2}$ $\displaystyle\leq 2\cdot\operatorname*{\mathbb{E}}\left\|v_{t}\right\|^{2}+2 \cdot\operatorname*{\mathbb{E}}\left\|v^{*}(t)-v_{t}\right\|^{2}$ $\displaystyle\lesssim\operatorname*{\mathbb{E}}\left\|v_{t}\right\|^{2}+ \operatorname*{\mathbb{E}}\left\|u^{*}(t)-u_{t}\right\|^{2}\lesssim d+W_{2}^{2 }(p,q).\qed$ $$

###### Theorem B.11 .

Let $\beta\geq 1$ be an adjustable parameter. Algorithm 8 with parameter $h=\frac{1}{\sqrt{8L}}$,$R=\beta\cdot\Theta(\frac{\sqrt{d}}{\varepsilon})$, $K=4\cdot\log(R)$ and $T_{N}\lesssim\frac{1}{\sqrt{L}}$ has discretization error

$$ $\textup{{TV}}(\widehat{x}_{N},x^{*}(T_{N}))\lesssim\sqrt{\textup{{KL}}( \widehat{x}_{N},x^{*}(T_{N}))}\lesssim\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L }}+\frac{\varepsilon}{\beta}+\frac{\varepsilon}{\beta\sqrt{d}}\cdot W_{2}(p,q).$ $$

###### Proof.

Since $T_{N}\lesssim\frac{1}{\sqrt{L}}$ and $h=\Theta(\frac{1}{\sqrt{L}})$, $N=O(1)$. Plugging Lemma B.10 into Lemma B.8, we get that

$$ $\displaystyle\textup{{KL}}(\widehat{x}_{N},x^{*}(T_{N}))$ $\displaystyle\lesssim\frac{T_{N}}{\sqrt{L}}\cdot\left(\varepsilon_{\mathrm{sc} }^{2}+L^{2}\left(\frac{\gamma dh^{3}}{R^{4}}+\frac{h^{2}}{R^{2}}\cdot(d+W_{2}^ {2}(p,q))+\frac{h^{4}}{R^{4}}\cdot(L^{2}W_{2}^{2}(p,q)+Ld)\right)\right)$ $\displaystyle\lesssim\frac{1}{L}\cdot\left(\varepsilon_{\mathrm{sc}}^{2}+\frac {d}{R^{4}}+\frac{Ld}{R^{2}}+\frac{Ld}{R^{4}}+\left(\frac{L}{R^{2}}+\frac{L^{2} }{R^{4}}\right)\cdot W_{2}^{2}(p,q)\right)$ $\displaystyle=\frac{\varepsilon_{\mathrm{sc}}^{2}}{L}+\frac{\varepsilon^{2}}{ \beta^{2}}+\frac{\varepsilon^{2}}{\beta^{2}d}W_{2}^{2}(p,q).$ $$

The first to second line is by combining terms and setting $h=\Theta(\frac{1}{\sqrt{L}})$, $\gamma=\Theta(\sqrt{L})$, and the second to third line is by setting $R=\beta\cdot\Theta(\frac{\sqrt{d}}{\varepsilon})$. Taking the square root of $\textup{{KL}}(\widehat{x}_{N},x^{*}(T_{N}))$ yields the claim.
∎

###### Theorem B.12 .

Let $\beta\geq 1$ be an adjustable parameter.
When Algorithm 8 is initialized at $(\widehat{x}_{0},\widehat{v}_{0})\sim p\otimes\mathcal{N}(0,I_{d})$, there exists parameters $h=\frac{1}{\sqrt{8L}}$,$R=\beta\cdot\Theta(\frac{\sqrt{d}}{\varepsilon})$, $K=\Theta(\log(\frac{\beta^{2}d}{\varepsilon^{2}}))$ and $T_{N}\lesssim\frac{1}{\sqrt{L}}$ such that the total variation distance between the final output of Algorithm 8 and the true distribution can be bounded as

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N},x_{t_{N}})\lesssim\frac{\varepsilon _{\mathrm{sc}}}{\sqrt{L}}+\frac{\varepsilon}{\beta}+\sqrt{L}\cdot W_{2}(p,q).$ $$

###### Proof.

By triangle inequality, $\textup{{TV}}(\widehat{x}_{N},x_{t_{N}})\leq\textup{{TV}}(\widehat{x}_{N},x^{*
}(T_{N}))+\textup{{TV}}(x^{*}(T_{N}),x_{t_{N}})$. Combining Lemma B.7 and Theorem B.11 yields the claim.
∎

#### Input parameters:

- •
Starting sample $(\widehat{x}_{0},\widehat{v}_{0})\sim p\otimes\mathcal{N}(0,I_{d})$, Number of steps $N$, Step size $h$, Score estimates $\widehat{s}\approx\nabla\ln q$, Number of midpoint estimates $R$, $\delta:=\frac{h}{R}$

- 1.
For $n=0,\dots,N-1$:
(a)
Let $t_{n}=nh$
(b)
Let $(\widehat{x}^{(k)}_{n,i},\widehat{v}^{(k)}_{n,i})$ represent the algorithmic estimate of $(x_{t_{n}+ih/R},v_{t_{n}+ih/R})$ at iteration $k$.
(c)
Let $(\zeta^{x},\zeta^{v})$ be a correlated gaussian vector corresponding to change caused by the Brownian motion term in $h/R$ time (see more detail in [])
(d)
For $i=0,\cdots,R$ in parallel: Let $(\widehat{x}^{(0)}_{n,i},\widehat{v}^{(0)}_{n,i})=(\widehat{x}_{n},\widehat{v}
_{n})$
(e)
For $k=1,\cdots,K$:
For $i=1,\cdots,R$ in parallel:
$\widehat{x}^{(k)}_{n,i}:=\widehat{x}^{(k-1)}_{n,i-1}+\frac{1-\exp(-\gamma h/R)
}{\gamma}\cdot\widehat{v}^{(k-1)}_{n,i-1}-\frac{h/R-(1-\exp(-\gamma h/R))/
\gamma}{\gamma}\cdot\widehat{s}(x^{k-1}_{n,j})+\zeta^{x}$
$\widehat{v}^{(k)}_{n,i}=\exp(-\gamma h/R)\cdot\widehat{v}^{(k-1)}_{n,{i-1}}-
\frac{1-\exp(-\gamma h/R)}{\gamma}\cdot\widehat{s}(x^{(k-1)}_{n,i-1})+\zeta^{v}$
(f)
$(\widehat{x}_{n+1},\widehat{v}_{n+1})=(\widehat{x}^{K}_{n,R},\widehat{v}^{K}_{
n,R})$
- 2.
Let $t_{N}=Nh$
- 3.
Return $\widehat{x}_{N},t_{N}$.

###### Lemma B.7 ( [ CCL + 23a ] , Lemma 9) .

###### Lemma B.8 ( [ ACV24 ] , Theorem 20, Implicit) .

###### Lemma B.9 ( [ CCL + 23a ] , Lemma 10) .

###### Lemma B.10 .

###### Proof.

###### Theorem B.11 .

###### Proof.

###### Theorem B.12 .

###### Proof.

### B.3 End-to-end analysis

Figure: Algorithm 9 ParallelAlgorithm

###### Theorem B.13 (Parallel End to End Error) .

By setting $T=\Theta\left(\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{\varepsilon^{2}}
\right)\right)$, $T_{corr}=\frac{1}{\sqrt{L}}$, $\delta=\Theta\left(\frac{\varepsilon^{2}}{L^{2}(d\lor\mathfrak{m}_{2}^{2})}\right)$, and $\beta_{1}=\beta_{2}=\Theta\left(L\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{
\varepsilon^{2}}\right)\right)$ in Algorithm 7 and Algorithm 8, when $\varepsilon_{\mathrm{sc}}\lesssim\widetilde{\Theta}(\frac{\varepsilon}{\sqrt{L
}})$, the total variation distance between the output of Algorithm 9 and the target distribution $x_{0}\sim q^{*}$ is

$$ $\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{0})\lesssim\varepsilon,$ $$

with iteration complexity $\widetilde{\Theta}(L\cdot\log^{2}\left(\frac{Ld\lor\mathfrak{m}_{2}^{2}}{
\varepsilon}\right))$.

###### Proof.

Let $x_{t_{n}}$ be the result of running the true ODE for time $T-t_{n}$, starting from $x_{T}\sim q^{*}$. Let $y^{\prime}_{n}$ be the result of running the predictor step in step $n-1$ of Algorithm 9, starting from $x_{t_{n-1}}\sim q_{t_{n-1}}$ and start time $t_{n-1}$. In addition, let $\widehat{y}_{n}$ be the result of the corrector step in step $n-1$ of Algorithm 9, starting from $y^{\prime}_{n}$.

We will first bound the error in one predictor + corrector step that starts at $t_{n-1}=T-(n-1)/L$. By triangle inequality of TV distance and data processing inequality (applied to $\widehat{x}_{n}$ and $\widehat{y}_{n}$),

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{n},x_{t_{n}})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{n},\widehat{y}_{n})+\textup{{TV}}( \widehat{y}_{n},x_{t_{n}})$ $\displaystyle\leq\textup{{TV}}(\widehat{x}_{n-1},x_{t_{n-1}})+\textup{{TV}}( \widehat{y}_{n},x_{t_{n}})$ (29) $$

By Theorem B.6 parametrized by $\beta_{1}$ and Theorem B.12 parametrized by $\beta_{2}$,

$$ $\displaystyle\textup{{TV}}(\widehat{y}_{n},x_{t_{n}})$ $\displaystyle\lesssim\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}}+\frac{ \varepsilon}{\beta_{2}}+\sqrt{L}\cdot W_{2}(y^{\prime}_{n},x_{t_{n}})$ $\displaystyle\lesssim\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}}+\frac{ \varepsilon}{\beta_{2}}+\sqrt{L}\left(\frac{\varepsilon_{\mathrm{sc}}}{L}+ \frac{\varepsilon}{\beta_{1}\sqrt{L}}\right)$ $\displaystyle\lesssim\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}}+\frac{ \varepsilon}{\min(\beta_{1},\beta_{2})}.$ $$

The first line is by Theorem B.12, and the first to second line is by Theorem B.12. Next, note that at the beginning of the process, $t_{0}=T$, and at the end of the process, $t_{N_{0}+1}=\delta$. By induction on Equation 29,

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{0})$ $\displaystyle\leq\textup{{TV}}(x_{0},x_{t_{N_{0}+1}})+\textup{{TV}}(\widehat{x }_{N_{0}+1},x_{t_{N_{0}+1}})$ $\displaystyle\leq\textup{{TV}}(x_{0},x_{\delta})+\textup{{TV}}(x_{T},\mathcal{ N}(0,I_{d}))+\sum_{n=1}^{N_{0}+1}\textup{{TV}}(\widehat{y}_{n},x_{t_{n}})$ $\displaystyle\leq\textup{{TV}}(x_{0},x_{\delta})+\textup{{TV}}(x_{T},\mathcal{ N}(0,I_{d}))+N_{0}\cdot\left(\frac{\varepsilon_{\mathrm{sc}}}{\sqrt{L}}+\frac{ \varepsilon}{\min(\beta_{1},\beta_{2})}\right).$ $$

By Lemma A.9, $\textup{{TV}}(x_{T},\mathcal{N}(0,I_{d}))\lesssim(\sqrt{d}+\mathfrak{m}_{2})
\exp(-T)$. By [, Lemma 6.4], $\textup{{TV}}(x_{0},x_{\delta})\leq\varepsilon$. Therefore by setting $T=\Theta\left(\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{\varepsilon^{2}}
\right)\right)$, $N_{0}=\Theta\left(L\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{\varepsilon^{2}}
\right)\right)$ and $\beta_{1}=\beta_{2}=\Theta\left(L\log\left(\frac{d\lor\mathfrak{m}_{2}^{2}}{
\varepsilon^{2}}\right)\right)$ in Algorithm 7 and Algorithm 8, when $\varepsilon_{\mathrm{sc}}\lesssim\widetilde{\Theta}(\frac{\varepsilon}{\sqrt{L
}})$, we obtain $\textup{{TV}}(\widehat{x}_{N_{0}+1},x_{0})\lesssim\varepsilon$.

The iteration complexity of Algorithm 9 given above parameters is roughly number of predictor-corrector steps times the iteration complexity in one predictor-corrector step. Note that in any corrector step and any predictor step except the last one, only $N=O(1)$ number of sub-steps are taken, therefore the iteration complexity of one predictor step (except the last step) is
$\Theta(\log(\frac{\beta_{1}\sqrt{d}}{\varepsilon}))$ and iteration complexity of one corrector step is $\Theta(\log(\frac{\beta_{2}^{2}d}{\varepsilon^{2}}))$. In the last predictor step, the number of steps taken is $O\left(\log\left(\frac{1}{\delta L}\right)\right)=O(\log(L)+T)$, and thus the iteration complexity is $\Theta((\log(L)+T)\cdot\log(\frac{\beta_{1}\sqrt{d}}{\varepsilon}))$. We conclude that the total iteration complexity of Algorithm 9 is

$$ $LT\cdot\left(\Theta(\log(\frac{\beta_{1}\sqrt{d}}{\varepsilon}))+\Theta(\log( \frac{\beta_{2}^{2}d}{\varepsilon^{2}}))\right)+\Theta\left((\log(L)+T)\cdot \log\left(\frac{\beta_{1}\sqrt{d}}{\varepsilon}\right)\right)=\widetilde{ \Theta}\left(L\log^{2}\left(\frac{Ld\lor\mathfrak{m}_{2}^{2}}{\varepsilon} \right)\right)\,.$ $$

∎

#### Input parameters:

- •
Start time $T$, End time $\delta$, Corrector steps time $T_{\mathrm{corr}}\lesssim 1/\sqrt{L}$, Number of predictor-corrector steps $N_{0}$, Score estimates $\widehat{s}_{t}$

- 1.
Draw $\widehat{x}_{0}\sim\mathcal{N}(0,I_{d})$.
- 2.
For $n=0,\dots,N_{0}$:
(a)
Starting from $\widehat{x}_{n}$, run Algorithm 7 with starting time $T-n/L$ with total time $\min(1/L,T-n/L-\delta)$. Let the result be $\widehat{x}_{n+1}^{\prime}$.
(b)
Starting from $\widehat{x}_{n+1}^{\prime}$, run Algorithm 8 for total time $T_{\mathrm{corr}}$ and score estimate $\widehat{s}_{T-(n+1)/L}$ to obtain $\widehat{x}_{n+1}$.
- 3.
Return $\widehat{x}_{N_{0}+1}$.

###### Theorem B.13 (Parallel End to End Error) .

###### Proof.

## Appendix C Log-concave sampling in total variation

In this section, we give a simple proof, using our observation about trading off the time spent on the predictor and corrector steps, of an improved bound for sampling from a log-concave distribution in total variation. Note that for this section, we assume that $\widehat{s}$ is the *true score* of the distribution and is known, as is standard in the log-concave sampling literature.

We begin by recalling Shen and Lee’s randomized midpoint method applied to approximate the underdamped Langevin process, for log-concave sampling in the Wasserstein metric [] in Algorithm 10.

Figure: Algorithm 10 RandomizedMidpointMethod []

###### Theorem C.1 (Theorem 3 of [ SL19 ] , restated) .

Let $\widehat{s}=\nabla\ln p$, the score function of a log-concave distribution $p$ be such that $0\preccurlyeq m\cdot I_{d}\preccurlyeq J_{\widehat{s}}(x)\preccurlyeq L\cdot I
_{d}$, for the Jacobian $J_{\widehat{s}}$ of $\widehat{s}$. Let $\widehat{x}_{0}$ be the root of $\widehat{s}$, and $v_{0}=0$. Let $\kappa=\frac{L}{m}$ be the condition number. For any $0<\varepsilon<1$, if we set the step size of Algorithm 10 as $h=C\min\left(\frac{\varepsilon^{1/3}m^{1/6}}{d^{1/6}\kappa^{1/6}}\log^{-1/6}
\left(\frac{d}{\varepsilon m}\right),\frac{\varepsilon^{2/3}m^{1/3}}{d^{1/3}}
\log^{-1/3}\left(\frac{d}{\varepsilon m}\right)\right)$ for some small constant $C$ and run the algorithm for $N=\frac{4\kappa}{h}\log\frac{20d}{\varepsilon^{2}m}\leq\widetilde{O}\left(
\frac{\kappa^{7/6}d^{1/6}}{\varepsilon^{1/3}m^{1/6}}+\frac{\kappa d^{1/3}}{
\varepsilon^{2/3}m^{1/3}}\right)$ iterations, then Algorithm 10 after $N$ iterations can generate $\widehat{x}_{N}$ such that

$$ $\displaystyle W_{2}(\widehat{x}_{N},x)\leq\varepsilon$ $$

where $x\sim p$.

Now, we make the following simple observation – if we run the corrector step from Section A.2 for a short time, we can convert the above Wasserstein guarantee to a TV guarantee. We carefully trade off the time spent on the Randomized Midpoint step above and the corrector step to obtain the improved dimension dependence. Our final algorithm is given in Algorithm 11.

Figure: Algorithm 11 LogConcaveSampling []

We obtain the following guarantee with our improved dimension dependence of $\widetilde{O}(d^{5/12})$.

###### Theorem C.2 (Log-Concave Sampling in Total Variation) .

Let $\widehat{s}=\nabla\ln p$ be the score function of a log-concave distribution $p$ such that $0\preccurlyeq m\cdot I_{d}\preccurlyeq J_{\widehat{s}}(x)\preccurlyeq L\cdot I
_{d}$ for the Jacobian $J_{\widehat{s}}$ of $\widehat{s}$. Let $\kappa=\frac{L}{m}$ be the condition number. For any $\varepsilon<1$, if we set $h_{\text{rand}}=C\left(\frac{\varepsilon^{2/3}}{d^{5/12}\kappa^{1/3}}\log^{-1/
3}\left(\frac{d\kappa}{\varepsilon}\right)\right)$ for a small constant $C$, $N_{\text{rand}}=\frac{4\kappa}{h}\log\frac{20d\kappa}{\varepsilon^{2}}\leq
\widetilde{O}\left(\frac{\kappa^{4/3}d^{5/12}}{\varepsilon^{2/3}}\right)$, $h_{\mathrm{corr}}=\widetilde{O}\left(\frac{\varepsilon}{d^{17/36}\sqrt{L}}\right)$ and $T_{\mathrm{corr}}=O\left(\frac{1}{\sqrt{L}d^{1/18}}\right)$, we have that Algorithm 11 returns $\widehat{x}_{N_{\text{rand}}}$ with

$$ $\displaystyle\textup{{TV}}(\widehat{x}_{N_{\text{rand}}},x)\lesssim\varepsilon$ $$

for $x\sim p$. Furthemore, the total iteration complexity is $\widetilde{O}\left(d^{5/12}\left(\frac{\kappa^{4/3}}{\varepsilon^{2/3}}+\frac{
1}{\varepsilon}\right)\right)$.

###### Proof.

By Theorem C.1, we have, for our setting of $N_{\text{rand}}$ and $h_{\text{rand}}$ that, at the end of step $2$ of Algorithm 11,

$$ $W_{2}(\widehat{x}_{N_{\text{rand}}}^{\prime},x)\leq\frac{\varepsilon}{d^{1/12} \sqrt{L}}\,.$ $$

Then, by the first part of Corollary A.7,

$$ $\textup{{TV}}(\widehat{x}_{N_{\text{rand}}},x)\lesssim\varepsilon+\sqrt{L}d^{1 7/36}\cdot\left(\frac{\varepsilon}{d^{17/36}\sqrt{L}}\right)\lesssim \varepsilon\,.$ $$

Our iteration complexity is bounded by $N_{\text{rand}}+\frac{T_{\mathrm{corr}}}{h_{\mathrm{corr}}}=\widetilde{O}\left
(\frac{\kappa^{4/3}d^{5/12}}{\varepsilon^{2/3}}+\frac{d^{5/12}}{\varepsilon}\right)$ as claimed.
∎

#### Input parameters:

- •
Starting sample $\widehat{x}_{0}$, Starting $v_{0}$, Number of steps $N$, Step size $h$, Score function $\widehat{s}$, $u=\frac{1}{L}$.

- 1.
For $n=0,\dots,N-1$:
(a)
Randomly sample $\alpha$ uniformly from $[0,1]$.
(b)
Generate Gaussian random variable $(W_{1}^{(n)},W_{2}^{(n)},W_{3}^{(n)})\in\mathbb{R}^{3d}$ as in Appendix A of [].
(c)
Let $\widehat{x}_{n+\frac{1}{2}}=\widehat{x}_{n}+\frac{1}{2}\left(1-e^{-2\alpha h}
\right)v_{n}-\frac{1}{2}u\left(\alpha h-\frac{1}{2}\left(1-e^{-2(h-\alpha h)}
\right)\right)\widehat{s}(x_{n})+\sqrt{u}W_{1}^{(n)}$.
(d)
Let $\widehat{x}_{n+1}=\widehat{x}_{n}+\frac{1}{2}\left(1-e^{-2h}\right)v_{n}-\frac
{1}{2}uh\left(1-e^{-2(h-\alpha h)}\right)\widehat{s}(x_{n+\frac{1}{2}})+\sqrt{
u}W_{2}^{(n)}$.
(e)
Let $v_{n+1}=v_{n}e^{-2h}-uhe^{-2(h-\alpha h)}\widehat{s}(x_{n+\frac{1}{2}})+2\sqrt
{u}W_{3}^{(n)}$.
- 2.
Return $\widehat{x}_{N}$.

###### Theorem C.1 (Theorem 3 of [ SL19 ] , restated) .

#### Input parameters:

- •
Number of Randomized Midpoint steps $N_{\text{rand}}$, Corrector steps Time $T_{\mathrm{corr}}\lesssim\frac{1}{\sqrt{L}}$, Randomized Midpoint Step size $h_{\text{rand}}$, Corrector step size $h_{\mathrm{corr}}$, Score function $\widehat{s}$.

- 1.
Let $\widehat{x}_{0}$ be the root of $\widehat{s}$, and let $v_{0}=0$.
- 2.
Run Algorithm 10 with $N_{\text{rand}}$ steps and step size $h_{\text{rand}}$, using $\widehat{x}_{0},v_{0}$, and let the result be $\widehat{x}_{N_{\text{rand}}}^{\prime}$.
- 3.
Run Algorithm 5 starting from $\widehat{x}_{N_{\text{rand}}}^{\prime}$ for time $T_{\mathrm{corr}}$, using step size $h_{\mathrm{corr}}$. Let the result be $\widehat{x}_{N_{\text{rand}}}$.
- 4.
Return $\widehat{x}_{N_{\text{rand}}}$.

###### Theorem C.2 (Log-Concave Sampling in Total Variation) .

###### Proof.

## Appendix D Helper lemmas

###### Lemma D.1 (Corollary $1$ of [ CCL + 23a ] ) .

For the ODE

$$ $\mathrm{d}x_{t}=\left(x_{t}+\nabla\ln q_{t}(x_{t})\right)\mathrm{d}t\,,$ $$

if $L\geq 1$ and $\operatorname*{\mathbb{E}}\left[\|\nabla^{2}\log q_{t}(x)\|^{2}\right]\leq L$, we have, for $0<s<t$ and $h=t-s$,

$$ $\operatorname*{\mathbb{E}}\left[\|\nabla\ln q_{t}(x_{t})-\nabla\ln q_{s}(x_{s} )\|^{2}\right]\lesssim L^{2}dh^{2}\left(L\lor\frac{1}{t}\right)\,.$ $$

###### Lemma D.2 (Implicit in Lemma $4$ of [ CCL + 23a ] ) .

Suppose $L\geq 1$, $h\lesssim\frac{1}{L}$ and $t_{0}-h\geq t_{0}/2$. For ODEs starting at $x_{t_{0}}=\widehat{x}_{t_{0}}$, where

$$ $\displaystyle\mathrm{d}x_{t}$ $\displaystyle=\left(x_{t}+\nabla\ln q_{t}(x_{t})\right)\mathrm{d}t$ $\displaystyle\mathrm{d}\widehat{x}_{t}$ $\displaystyle=\left(x_{t}+\widehat{s}_{t_{0}}(\widehat{x}_{t_{0}})\right) \mathrm{d}t,$ $$

we have

$$ $\displaystyle\operatorname*{\mathbb{E}}\|x_{t_{0}-h}-\widehat{x}_{t_{0}-h}\|^{ 2}\lesssim h^{2}\left(L^{2}dh^{2}\left(L\lor\frac{1}{t_{0}}\right)+\varepsilon _{\mathrm{sc}}^{2}\right).$ $$

###### Lemma D.3 (Lemma $B.1.$ of [ GLP23 ] , restated) .

Let $p_{0}$ be a distribution over $\mathbb R^{d}$. For $x_{0}\sim p_{0}$, let $x_{t}=x_{0}+z_{t}\sim p_{t}$ for $z_{t}\sim\mathcal{N}(0,tI_{d})$ independent of $x_{0}$. Then,

$$ $\displaystyle\frac{p_{t}(x_{t}+\varepsilon)}{p_{t}(x_{t})}=\operatorname*{ \mathbb{E}}_{z_{t}|x_{t}}[e^{\frac{\varepsilon^{T}z_{t}}{t}-\frac{\| \varepsilon\|^{2}}{2t}}]$ $$

and

$$ $\displaystyle\nabla\ln p_{t}(x_{t})=\operatorname*{\mathbb{E}}_{z_{t}|x_{t}} \left[-\frac{z_{t}}{t}\right]$ $$

###### Lemma D.4 .

For $q_{t}(y_{t})\propto p_{e^{2t}-1}(e^{t}y_{t})$, for $z_{t}\sim\mathcal{N}(0,(e^{2t}-1)I_{d})$, we have

$$ $\displaystyle\nabla\ln q_{t}(y_{t})=e^{t}\nabla\ln p_{e^{2t}-1}(e^{t}y)=e^{t} \operatorname*{\mathbb{E}}_{z_{t}|e^{t}y_{t}}\left[\frac{-z_{t}}{e^{2t}-1}\right]$ $$

Furthermore,

$$ $\displaystyle\operatorname*{\mathbb{E}}_{y_{t}\sim q_{t}}\left[\|\nabla\ln q_{ t}(y_{t})\|^{2}\right]\lesssim\frac{d}{t}$ $$

###### Proof.

The first claim is an immediate consequence of the definition of $q_{t}$ and Lemma D.3. For the second claim, note that

$$ $\displaystyle\operatorname*{\mathbb{E}}_{y_{t}\sim q_{t}}\left[\|\nabla\ln q_{ t}(y_{t})\|^{2}\right]$ $\displaystyle=e^{2t}\operatorname*{\mathbb{E}}_{y_{t}\sim q_{t}}\left[\left\| \operatorname*{\mathbb{E}}_{z_{t}|e^{t}y_{t}}\left[\frac{-z_{t}}{e^{2t}-1} \right]\right\|^{2}\right]$ $\displaystyle\leq e^{2t}\operatorname*{\mathbb{E}}_{y_{t}\sim q_{t}}\left[ \operatorname*{\mathbb{E}}_{z_{t}|e^{t}y_{t}}\left[\frac{\|z_{t}\|^{2}}{(e^{2t }-1)^{2}}\right]\right]$ $\displaystyle=e^{2t}\operatorname*{\mathbb{E}}_{z_{t}}\left[\frac{\|z_{t}\|^{2 }}{(e^{2t}-1)^{2}}\right]$ $\displaystyle=\frac{e^{2t}\cdot d}{e^{2t}-1}\quad\text{since $z_{t}\sim \mathcal{N}(0,(e^{2t}-1)I_{d})$}$ $\displaystyle\lesssim\frac{d}{t}\,.\qed$ $$

###### Lemma D.1 (Corollary 1 1 1 1 of [ CCL + 23a ] ) .

###### Lemma D.2 (Implicit in Lemma 4 4 4 4 of [ CCL + 23a ] ) .

###### Lemma D.3 (Lemma B ⁢ .1 . 𝐵 .1 B.1. italic_B .1 . of [ GLP23 ] , restated) .

###### Lemma D.4 .

###### Proof.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: - is it possible to improve the total query complexity of paralell version to $d^{4/9}$?

- What are the primary challenges associated with removing the smoothness assumption?

- In the parallel version of sampling for diffusion models and the log-concave sampling algorithm, the predictor-corrector framework runs only once, whereas in the sequential version, the predictor-corrector loop is repeated $N_0$ times. Why does this difference occur?
soundness: 3
strengths: - Their method achieves significantly better sequential complexity and parallel complexity.

- Overall I found the paper well-written and clear.
summary: The authors explore sampling for diffusion with mild assumptions and introduce both sequential and parallel algorithms, which extend the ODE-based predictor-corrector framework initially presented in [CCL+23a]. In their sequential approach, they adeptly utilize the randomized midpoint method in[SL19] to discretize the probability flow ODE during the predictor phase, and use the corrector method outlined in [CCL+23a] with a short time $T$. Impressively, this method achieves a complexity of $O(L^{5/3}d^{4/9}/\varepsilon)$, which marks a substantial improvement over the previous best complexity of $O(L^2d^{1/2}/\varepsilon)$. For the parallel algorithm, the authors also apply a similar parallel randomized midpoint method from [SL19] for the predictor phase and use the parallel underdamped Langevin method from [ACV24] as the corrector. This innovative combination gives the first provable guarantees for parallel sampling within diffusion models. Furthemore, this paper gives an improved bound for log-concave sampling.
weaknesses: - Their sequential result requires a slightly stronger assumption on the score estimation error, but I think it is fine.

- The technical contributions are limited to combining existing ODE-based predictor-corrector framework and discrete schemes.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: **Space Complexity and Scalability.** Can the authors elaborate on how the proposed method would scale with large data distributions in practical applications? How feasible is this method for scenarios requiring substantial memory or computation resources?

**Empirical Benchmarks.** Could the authors provide insights on the types of empirical tests they recommend or any preliminary results? Testing on standard datasets or benchmarks in score-based generative models would be particularly valuable for evaluating this method’s efficiency.
soundness: 3
strengths: **Theoretical Innovation.** The manuscript extends a randomized midpoint technique, traditionally used in other domains, to diffusion sampling, presenting significant advancements in sample complexity for smooth distributions.

**Parallel Sampling Framework.** Including a parallel sampling approach with provable efficiency bounds is particularly notable, potentially paving the way for scalable implementations in large-scale generative modeling tasks.

**Solid Analysis.** The authors offer a rigorous mathematical foundation for their method, with thorough proofs and well-defined assumptions on score functions and dimensionality dependence.
summary: The manuscript introduces a novel approach for enhancing sampling efficiency in diffusion models by incorporating a randomized midpoint method. This method, inspired by Shen and Lee's work on log-concave sampling, is adapted to diffusion processes to reduce sample complexity and improve iteration bounds. The authors present theoretical guarantees for both sequential and parallel sampling, addressing the limitations of previous diffusion model sampling techniques.
weaknesses: **Assumptions and Practical Applicability.** While the theoretical advancements are clear, the reliance on strong Lipschitz assumptions and high accuracy in score estimates might limit practical applicability. The authors could consider discussing the feasibility of these assumptions in real-world diffusion models.

**Lack of Experimental Validation.** Despite the paper's focus on theoretical results, experimental validation would have been beneficial. Empirical benchmarks could highlight practical gains and offer insights into scenarios where the proposed method outperforms existing approaches.

**Parallelization Overhead.** While the approach achieves polylogarithmic round complexity, it requires significant parallel resources to fully realize these benefits. The paper could improve by addressing potential bottlenecks and clarifying hardware implications.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Can the authors provide simulations of their methods? It would be interesting to see if their method works in practice.  This concerns both the runtime and the quality of the generated samples.
soundness: 4
strengths: The paper obtains guarantees, in the TV norm, for solving the reverse SDE in diffusion models, with improved dimension dependence on the runtime.

The paper obtains improved guarantees, which are logarithmic in dimension, for parallelized solving of reverse SDE in diffusion models.

As a by-product of their analysis the authors also obtain improved bounds for sampling from log-concave distributions using the Langevin MCMC algorithm.

The paper introduces a new technique (the randomized midpoint method, originally developed for log-concave sampling via MCMC algorithms) to the area of generative diffusion models.
summary: This paper uses the randomized midpoint method to solve the reverse-diffusion SDE in generative diffusion models, to generate samples.  The randomized midpoint method is a more efficient method of solving SDEs than methods which incur a first-order error, such as the Euler-Maruyama integrator, as long as the coefficients are smooth functions.   The randomized midpoint method was originally developed for implementing the Hamiltonian Monte Carlo Markov Chain algorithm, but the here the authors apply the same technique in the setting of generative diffusion models.  This allows the authors to obtain faster runtime bounds for sampling (in the total variation norm) from the reverse-diffusion SDE in diffusion models, provided the certain smoothness assumptions are satisfied.  As a by-product of their analysis the authors also obtain improved bounds for sampling from log-concave distributions using the Langevin MCMC algorithm.
weaknesses: The paper does not include any numerical simulations.  Numerical simulations would help evaluate whether their algorithm works well in practice (both in terms of runtime and in terms of sample generation quality).

Also, while the authors adapt a new technique to the generative diffusion models, this is still an existing technique, but was used in a different area (MCMC for log-concave sampling).  That being said, I still feel the technical contribution is good since the authors introduce an existing technique to a new area.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: 1. As claimed in [1], the linear dependency on $d$ is optimal for KL divergence, thus $\sqrt{d}$ for TV. Can you remark on this compared with the better than $\sqrt{d}$ result in this paper?

2. In the proof of Theorem A.10, given the value of the hyperparameters in Line 1109, can you explain more about how you can get the iteration complexity? Is it max $\\{1/h_{pred}, 1/h_{corr}\\}$? Then why the dependency on $d$ is not $d^{17/36}$, which is larger than $d^{5/12}$?
soundness: 3
strengths: 1. The organization is clear, highlighting the main contribution.

2. The analysis method is very novel, and it’s an interesting result to see that the dimensional dependency for sampling can be improved to be better than $\tilde O(\sqrt{d})$.

3. The paper shows that the analysis can also be extended to parallel algorithms.
summary: This paper studies diffusion models and introduces a faster sampling algorithm utilizing the randomized midpoint method within a predictor-corrector framework. In the predictor step, the paper introduces the randomized midpoint method, which allows for bounding the squared displacement into “bias” and“variance”terms. The corrector step follows the analysis in Chen et al.,2023b. This approach reduces discretization error, leading to an algorithm with improved theoretical dependency on the dimension $d$. Similar analysis can also be used in the setting of log-concave sampling in total variation distance, achieving a better dependency in $d$.
weaknesses: 1.The assumptions are a little strong.  It has been shown in [1][2] that the assumption of Lipschitz score can be relaxed. Moreover, it is often enough to assume that the true score is Lipschitz (Assumption 2.2). However, in this paper, the estimated score is also required to be Lipschitz. Although I don’t think it’s a big problem, can the authors remark about where this additional assumption is used?

2.There are issues with the proof of Lemma A.2. The argument in (Line 736 equation 17) holds only under pointwise estimation (the estimated score is accurate for any $x$). However, in Assumption 2.4, it is defined as the expectation. Similar conditions occur in Line 755. As all the lemmas after this depend on Lemma A.2, the main result is correct only with a stronger assumption, which is not reflected in the paper. 

3.It misses some background knowledge and discussion:

(1) The paper discusses log-concave sampling in section 3.6. However, there is no background knowledge in this paper. It is unclear what is the exact setting. Thus, the analysis in Appendix C is confusing. 

(2) In the main content, all the theorems are informal. Therefore, in the proof technique sections, like Sections 3.3, and 3.4, the target of the argument is ambiguous. I think at least a formal version of the theorem should be provided, and more details of the proof should be provided to help readers understand.

4.Some typos: I think Algorithms 1, 2, and 3 are exactly the same as Algorithms 4,5,6 in the appendix. Thus, the reference in Algorithm 3 Line 3 should be Algorithm 1 instead of Algorithm 4, right?

I'm willing to increase my score if the authors can solve my concerns, especially weakness 2.

----
[1] Benton et al. 2024, Nearly $d$-Linear Convergence Bounds for Diffusion Models via Stochastic Localization, ICLR2024

[2] Chen et al. 2023, Improved Analysis of Score-based Generative Modeling: User-Friendly Bounds under Minimal Smoothness Assumptions, ICML2023

Meta-review: This paper considers diffusion sampling via the randomized midpoint method for log-concave sampling. Authors improve the existing dimension dependence for sampling from arbitrary smooth distributions in total variation distance. Authors also consider a parallel version of this algorithm and provide respective guarantees. 


This paper was reviewed by three expert reviewers and received the following Scores/Confidence: 6/2, 8/3, 6/4, 6/4. Even though certain results are not surprising, I still think the paper is studying an interesting topic and the results are relevant to ICLR community. The following concerns were brought up by the reviewers:

- The idea is straightforward, and the technical contribution seems limited. Authors should explain the novelty in their work more clearly.
  
- The feasibility of the conditions should be discussed in detail.
  
- Authors discuss a 'new' intergrator and provide improved guarantees for it. I expected to see extensive numerical results that validates the improved performance. I suggest authors to extend their experiments in camera ready.


Authors should carefully go over reviewers' suggestions and address any remaining concerns in their final revision. Based on the reviewers' suggestion, as well as my own assessment of the paper, I recommend including this paper to the ICLR 2025 program.

Decision: Accept (Poster)