Subpolynomial query complexity for well-conditioned log-concave sampling

Proof. Strong convexity implies Gaussian tails and existence of a unique minimizer. Integration by parts gives \E_\mu (X-z_0)\cdot\nabla H(X)=n; strong monotonicity of \nabla H gives the moment bound. The globally Lipschitz drift gives a nonexplosive diffusion. For synchronous solutions, differentiating their squared distance and using strong monotonicity gives the contraction.

First suppose H is smooth and has bounded derivatives of orders at least two. For a smooth compactly supported h, let P_th(x)=\E[h(X_t)\mid X_0=x]. Differentiating the synchronous flow in the initial state gives bounded first and second derivatives on bounded time intervals. Itô’s formula and the semigroup identity yield

The equality with \mathcal AP_th also follows by differentiating P_{t+s}h=P_s(P_th) at s=0. Its at-most-linear growth is integrable against \mu. Integration by parts gives \int\mathcal AP_th\,d\mu=0, hence invariance. Coupling to an invariant starting point and using contraction shows P_th\to\mu h in L^2(\mu). Therefore

The flow derivative has operator norm at most e^{-mt}. Jensen’s inequality and invariance thus bound the inner integral by e^{-2mt}\int|\nabla h|^2d\mu, proving the Poincaré inequality.

For C^2 potentials, convolve with smooth compactly supported mollifiers of radii tending to zero. Both Hessian bounds are preserved; all higher derivatives are bounded for each mollification. The normalized laws converge in total variation by the uniform quadratic tail bounds (subtract the potential’s value at zero). The gradients converge uniformly since the original gradient is globally Lipschitz. Synchronous solutions of the mollified and original diffusions consequently converge uniformly in their common starting point on each bounded time interval, by Gronwall’s inequality. Testing invariance against bounded Lipschitz functions passes it to the limit. Poincaré first passes on smooth compact tests, then by cutoff and approximation to the stated class.

For (2.5), start the two diffusions with their respective invariant marginals, with any initial coupling, and use common Brownian motion. Writing a(t)=(\E|X_t-Y_t|^2)^{1/2} and \beta=(\E_\eta|\nabla H-\nabla\bar H|^2)^{1/2}, strong monotonicity and Cauchy–Schwarz give, in the usual upper-derivative sense,

The marginals remain fixed, so \Wtwo(\mu,\eta)\le a(t). Letting t tend to infinity proves the claim.□

Proof. Use conditional near-optimal couplings for the first inequality and their independent product for the second, followed by the triangle inequality in L^2. Measurability of such conditional comparisons can also be seen without an exact optimal-kernel selection: quantize each output law on a finite grid inside a box, sending the complement to zero. A minimizing plan on the finite sets can be chosen measurably by a fixed ordering of the finitely many basic feasible solutions. Lift it using the conditional measures in the cells, and refine the grids and boxes. Finite unconditional second moments and the triangle inequality give the asserted limit.

For deterministic centers x,y, the relative entropy of N(x,s^2I) from N(y,s^2I) is |x-y|^2/(2s^2). Pinsker’s inequality gives TV at most |x-y|/(2s). Integrate this estimate over any coupling and apply Cauchy–Schwarz to obtain (2.7), with room in its constant. The total-variation rules follow directly by testing measurable functions taking values in [0,1].□

Proof. We first prove a dimension-independent covariance bound for centered polynomial tensors. We then recover cumulants from those covariances. If \beta=0, applying the hypothesis to coordinate functions shows that S is almost surely constant, so every cumulant in the statement vanishes. We may therefore suppose \beta>0.

Define the symmetric Appell tensor polynomials by

For a coefficient tensor T of order j, let Q_T(s)=T:P_j(s), where the colon denotes contraction on all slots. Differentiating the generating function in s multiplies it by copies of \theta. Consequently

Here \operatorname{Sym}(T) is the average over all permutations of the slots; hence |\operatorname{Sym}(T)|_F\le|T|_F. Apply Poincaré to each component of \nabla^h Q_T and sum. Iterating from h=0 to h=j-1 yields

Only first derivatives in the displayed Poincaré inequality are used at each application.

In particular, if T and V have orders j,k\ge1, then

This is the direct covariance estimate. It concerns arbitrary coefficient tensors, rather than products of one-dimensional test directions, which is why it controls a matrix split.

To pass from covariance to cumulants while preserving a specified split, consider the two-variable generating function

Its derivative of orders j\ge1 in \theta and k\ge1 in \omega at (0,0) is the cross-moment tensor \E[P_j(S)\otimes P_k(S)]. Across the split between the \theta and \omega slots its operator norm is at most \beta^{(j+k)/2}j!k! by (3.3).

Both \mathcal G(\theta,0) and \mathcal G(0,\omega) are identically one. Thus every nonzero term of \mathcal G-1 contains variables from both groups. In the formal identity

only m\le n can contribute at total order n. Every differentiated product is a tensor product of the preceding cross-moment operators, up to permutations within the two sides. Its norm is therefore at most the product of their norms. To bound the total coefficient, distribute the n labeled derivative slots among at most n ordered nonempty factors. There are at most n^{n+1} such assignments; the product of the factorials within any assignment is at most n!. Enlarging an absolute constant bounds their total by C^n n^{Cn}.

The mixed derivatives of \log\mathcal G are precisely the corresponding derivatives of \log\E\exp((\theta+\omega)\cdot S), since the other two logarithms have only one variable group. They are the cumulants, and the two groups can represent any prescribed nontrivial split of the n slots. This proves the claim.□

Proof. The conditional potential in (3.1) has Hessian between (a^{-1}-l)I and (a^{-1}+l)I. By Lemma 2.2, its Poincaré constant is at most Ca. Under this conditional law the map

has derivative norm at most \sqrt2, since f is L-Lipschitz. Pulling back an ambient test function shows that the possibly singular law of S satisfies Poincaré with constant Ca. It has exponential moments of every linear functional: f has at most linear growth and the conditional law has Gaussian tails. Lemma 3.1 therefore applies. This use of repeated Poincaré does not differentiate f repeatedly. Each application pulls back a new ambient polynomial and uses only Df.

In the conditional density, the linear tilt parameter for z is \rho y/a. Repeated differentiation of an expectation with respect to a linear tilt gives joint cumulants. In coordinates this says

It can also be checked directly by differentiating the logarithm of the conditional moment-generating function with one auxiliary parameter coupled to f_i. Restrict the cumulant of S to one f/L coordinate and q ordinary coordinates. Lemma 3.1 gives a factor L a^{(q+1)/2}, and the prefactor a^{-q} gives (3.4).

All the conditional integrals used here are smooth in (\rho,y) for \rho<1. This follows by differentiating their densities: on compact parameter sets, every derivative inserts a polynomial with an integrable Gaussian majorant. Formula (3.7) gives the smooth extension at zero without division by \rho. For later integrations by parts we record a slightly stronger regularity fact. On each compact interval in \rho<1, conditional moments of any fixed order grow at most polynomially in y. Indeed, the conditional mode satisfies

so |z_0|\le C(|y|+r|g(x)|). The upper and lower conditional Hessian bounds compare the density and its normalizing integral around this mode to Gaussian integrals. This proves the claimed moment growth. The same argument applies to time derivatives of the conditional integrals. Dimension-dependent constants are harmless for this regularity justification; the quantitative estimates above remain dimension-uniform.

Conditional integration by parts gives \E[Z\mid Y=y]=\rho y-arM_\rho(y). Differentiating the marginal density then yields the first identity in (3.5). For \rho>0 it implies symmetry of DM_\rho, hence of U_1; continuity handles zero. Its derivative, together with (3.4), gives the bounded Hessian and (3.6). In particular, all these conclusions use only the original C^2 regularity of F.□

Proof. For the random interpolation used to define \xi_\rho,

The conditional mean of this derivative is -rM_\rho(Y) by the conditional integration-by-parts identity above. Thus, for compactly supported smooth \phi,

On the other hand, (3.4) with q=1 gives r\|DM_\rho\|_{\op}\le Cl. The velocity has at most linear growth, so its smooth ordinary differential equation has a global-in-space flow and inverse on every compact interval in \rho<1. To identify its pushforward, fix a smooth compactly supported terminal test and transport it backward along this flow. The transported test has zero material derivative. Its supports remain in a bounded set on the interval, by the linear-growth bound. The last display then shows its expectation under \xi_\rho is constant, which proves the pushforward assertion.□

Proof. Differentiate \E[f_i(Z)\phi(Y)]=\int (M_\rho)_i\phi\,d\xi_\rho. The right side differentiates using the transporting velocity from Proposition 3.3. On the left the random velocity is (Z-\rho Y)/a. Subtract its conditional mean from the flux multiplying \nabla\phi. The remaining coefficient is \Cov(f_i(Z),Z_k\mid Y)/a=(U_1)_{i,k} by (3.7). Integration by parts proves (4.1).

For the higher identity, coordinate differentiation gives

Insert (4.1) and multiply by \rho^{-q}. When all derivatives pass inside the adjoint divergence, the result is \nabla^*U_{q+1}. In the other terms they hit the score y+r\rho M_\rho. A derivative hitting its linear part contributes qU_q/\rho, canceled by differentiating the prefactor \rho^{-q}. The nonlinear score terms and the material/spatial commutator terms are identical after scaling: each is

This proves (4.2). Each contraction retains at least one free slot in each tensor.

Finally differentiate the adjoint identity against a transported scalar test, componentwise in the other slots. The terms involving \mathcal D\phi cancel. The gradient commutator is

Moving it to the other side by the adjoint relation gives (4.3). All coefficients in the final formulas are nonsingular at zero, so the extensions follow from Proposition 3.2.□

Proof. Call the two sides of the required output split the input and output sides. Orient the contracted index from one factor to the other. If all free slots of one factor lie on the same side, its contracted slot must be placed on the opposite side of that factor’s unfolding. These requirements are compatible: incompatibility would put every free slot of the result on one side, contrary to the hypothesis. If a factor already has free slots on both sides, it imposes no requirement on the orientation. With the resulting orientation, the contraction is composition of the two matrix actions, with identities on the unused input and output coordinates and with coordinate permutations. Identities and these permutations have operator norm one. Submultiplicativity proves the assertion.□

Proof. We induct on the finite construction of E, proving all derivative and moment orders simultaneously. Permutations and single contractions follow directly from the preceding lemma. For a divergence, we first integrate by parts in a high trace moment. The resulting coordinate sums can be arranged as a composition of matrix unfoldings with only one closing trace. We then choose the moment large enough to absorb that trace’s dimension factor.

For a normalized jet, \nabla^jU_q=\rho^jU_{q+j}, so Proposition 3.2 proves the assertion after enlarging a constant depending on q. Permutations cause no change. For a contraction, distribute the derivatives by Leibniz, use Lemma 4.2 on each summand, and use the inductive estimates at exponent 2p on both factors. The number of Leibniz summands is at most 2^j. Their weights and derivative orders add, giving (4.4) with an enlarged fixed constant.

It remains to consider E=\nabla^*V, where V has weight k_V and E has at least two slots. Commute the j derivatives past the divergence. Apart from

the terms are single contractions of \nabla^{j-h}V with an order-h derivative of the score, 1\le h\le j. Both factors retain slots. The score derivative is bounded pointwise by a constant of the form C^h h^{Ch}(\sqrt a)^{-(h-1)}, by (3.5)–(3.6). Lemma 4.2, the inductive hypothesis, and Leibniz multiplicities therefore give the required bound for all these terms, even with exponent k_V+j-1 in place of the larger allowed exponent k_V+j+1. It remains to bound T.

Removing the adjoints. Fix an unfolding of T, with both sides nonempty, and an integer s\ge1. Pointwise,

Expand the trace by coordinate indices. It contains 2s labeled primary factors, each an adjoint divergence of \nabla^jV. Their free slots are joined in the alternating trace cycle; the two bundles incident to each primary factor are nonempty.

Remove adjoint divergences one at a time by integration by parts. For a chosen primary, its adjoint differentiates the product of the other factors, never its own inside field. If a derivative hits a remaining adjoint, use

The first term leaves that adjoint to be processed later. The second consumes it and creates a Hessian factor joining the two divergence indices. A later derivative may also hit a previously created Hessian factor. After all adjoints have been removed, a typical term has primary factors \nabla^{j+e_i}V, 1\le i\le2s, and n_H factors \nabla^{h_v}\nabla^2H_\rho. The bookkeeping identity is

Indeed, a removed adjoint either adds one derivative or pairs with another adjoint to create a Hessian. Each step has at most Cs possible target factors and at most two commutator terms, so there are at most (Cs)^{Cs} resulting terms.

We use a graph only to describe these explicit coordinate sums. Its primary vertices are the 2s labeled tensor factors; a Hessian vertex represents one factor \nabla^{h_v}\nabla^2H_\rho. An edge represents a pair of tensor slots summed over a common coordinate, and a bundle represents several such edges. Three properties will control its contraction. All original trace bundles remain. Every new direct edge connects distinct primary factors, because an adjoint never differentiates its own inside field. Finally, each Hessian factor has at least two distinct primary neighbors, the pair that created it, and there are no Hessian-to-Hessian edges. Later derivatives only add further primary-to-Hessian edges. These statements remain true if some derivatives had already been commuted inside an unprocessed adjoint.

Opening the trace. Choose one bundle in the original trace cycle and cut it, leaving its b\ge1 indices unsummed at the two ends. All remaining original bundles form a chain through the 2s primary factors. Order these factors along that chain. For each Hessian factor, its two distinct creating primaries have different positions, so we can place it strictly between its earliest and latest primary neighbors. No edge joins two Hessians, so these placements impose no competing condition. Order Hessians occupying the same interval arbitrarily, and orient every uncut edge from its earlier factor to its later factor. The b cut indices enter the first primary and leave the last one.

Every tensor has both an input and an output in this ordering. An interior primary receives one original chain bundle and sends the next; at the two endpoints the cut bundle supplies the missing side. A Hessian receives an edge from an earlier neighbor and sends an edge to a later one. Thus all the matrix unfoldings used below are nontrivial splits, exactly as required by the induction hypothesis.

Here is an explicit matrix-composition interpretation. Immediately before processing a factor, retain one copy of \R^d for each index whose earlier endpoint has been processed but whose later endpoint has not. The initial space is (\R^d)^{\otimes b}, belonging to the cut inputs. Permute these tensor coordinates so that the incoming indices of the next factor are adjacent. Apply that factor’s matrix unfolding to them and the identity to all other retained coordinates; then permute the outputs into the new retained list. Matrix multiplication performs precisely the sums over the incoming indices. Induction over the factors therefore accounts for every uncut index once and leaves only the b cut outputs at the end. If Q denotes the resulting map from (\R^d)^{\otimes b} to itself, then

Indeed coordinate permutations and the intervening identity maps have operator norm one. Closing the original cut is now exactly a single matrix trace, so the absolute value of this completely contracted term is at most

This argument explains why other undirected cycles introduce no further dimension factor: their coordinate sums have already been performed by the matrix compositions. Figure 1 shows a possible ordering.

Uniform high-moment bookkeeping. The coordinate expansion and its contraction estimate are now complete. It remains to verify that their constants allow s to grow with d and j. Apply the induction hypothesis to every primary factor with moment exponent 2s, and use the pointwise Hessian estimates on all other factors. There are exactly 2s random primary factors, so Hölder’s inequality uses no larger exponent. The power of 1/\sqrt a before taking the 2s-th root is at most

Since e_i,h_v\le2s, all polynomial bases in the inductive estimates are at most C(s+j+\log(d+1)). The total exponent of those bases is at most

using (4.5). The term count (Cs)^{Cs} is absorbed in the same bound. The cut bundle uses only the original free slots of T=\nabla^*(\nabla^jV), not the new indices introduced by integration by parts. Consequently b\le j+C_V, where C_V depends only on the fixed number of slots of V, and not on s. We have proved

Choose s to be the smallest positive integer satisfying 2s\ge p and s\ge(j+C_V)\log(d+1), increasing C_V if necessary. Then d^{(j+C_V)/(2s)}\le e^{1/2} and

The expression on the right is bounded by a fixed multiple of (p+j+\log(d+1))^2. Taking the 2s-th root, using monotonicity of probability-space L^p norms, and increasing C_E now gives (4.4) with k=k_V+1.

All integration-by-parts identities above are legitimate for each finite j,s,d and each \rho<1. Indeed, the conditional-integral regularity in Proposition 3.2, together with the score formula, gives polynomial growth for every field and derivative appearing in the finite expansion. The laws \xi_\rho have Gaussian tails. Cutoffs therefore justify the identities with vanishing boundary terms. Uniformity in growing s is supplied by the explicit estimates after the identities, not by an unproved uniform regularity assertion. This completes the induction.□

Proof. First take a compactly supported smooth vector field, so that the output is scalar. Integrating by parts twice gives

The first term is bounded in absolute value by |\nabla V|_F^2, and the second by \|\nabla^2H\|_{\op}|V|^2. For tensor output, apply the scalar identity to each fixed collection of remaining indices and sum. Taking a square root proves (4.7) for compactly supported fields. For general V, multiply by smooth cutoffs \chi_R equal to one on the ball of radius R, with |\nabla\chi_R|\le C/R. Then \chi_RV\to V and \nabla(\chi_RV)\to\nabla V in L^2(\mu), while their adjoint divergences converge pointwise. Fatou’s lemma proves the stated estimate. All fields to which we apply it have polynomially growing derivatives and hence satisfy the required integrability.□

Proof. Iterate Lemma 4.1, using the product rule on contractions. Replacing U_q by \nabla^*U_{q+1} raises weight by two. Replacing it by one of the contracted products in (4.2) does not raise weight: the two jet weights sum to q-1. The commutator (4.3) inserts DM_\rho=\rho U_1 of weight zero on the divergence slot. Every additional jet factor L is accompanied by a factor r, so, after retaining one overall L, the extra factors are powers of l=rL\le1. Scalar coefficients are fixed polynomials in \rho and remain bounded under the fixed number of differentiations. Starting from a jet U_q, all divergences retain at least two free slots, and all contractions retain a slot of each factor. The resulting terms are therefore expressions of Theorem 4.4, of weight at most q-1+2m. Its case j=0 proves (4.9).

Starting instead from M, the first identity is \mathcal DM=\nabla^*U_1. Its outermost divergence remains a divergence to a vector under further material differentiations. Inside it, the same construction gives finitely many normalized expressions of weight at most 2(m-1), each with one overall factor L. Apply Lemma 4.5 at this outermost operation. For a tensor with at least two slots, its Frobenius norm is at most \sqrt d times the operator norm of an unfolding that puts one coordinate slot on one side: that matrix has at most d singular values. Thus Theorem 4.4, for j=0,1 and p=2, bounds the right side of (4.7) by

This is (4.8).□

Proof. Differentiate the transport equation with respect to its starting position. Gronwall’s inequality gives a bounded first derivative and an O(l) bound for its difference from I_d. For order k>1, the chain rule has terms consisting of an order-j derivative of the velocity and derivatives of the state of orders i_1,\ldots,i_j\ge1, where i_1+\cdots+i_j=k. By (3.4), the outer tensor has every split norm at most C_j lR^{-(j-1)}. Attach the inner tensors one at a time. Each attachment is a single contraction retaining the outer output and the inner derivative slots, so Lemma 4.2 applies without a dimension factor. Bound a first state derivative by a constant and use the induction hypothesis for all higher state derivatives. Since l\le1, the total power of R^{-1} is at most

and each forcing term retains the outer factor l. The term involving the unknown order-k state derivative is linear, with coefficient bounded by Cl, and this derivative starts at zero. Gronwall proves the assertion. The same argument with absolute interval lengths applies to transport in the reverse direction.

For the second assertion, write C_\rho=rDM_\rho(Y_\rho)=r\rho U_1(\rho,Y_\rho), so that B_\rho'=B_\rho C_\rho. The product rule expresses B_\rho^{(m)} as a finite sum of ordered products of B_\rho and k\ge1 factors C_\rho^{(j_1)},\ldots,C_\rho^{(j_k)}, with j_1+\cdots+j_k=m-k. The zero-order jet bound and (4.9) imply, for every fixed moment order,

Here the path estimates apply because inverse transport sends y\sim\xi_* to Y_\rho(y)\sim\xi_\rho. Apply Hölder with sufficiently high fixed moments to each product. Using l^k\le l and (1-\rho^2)^{-(m-k)}\le(1-\rho^2)^{-m} gives the claimed conservative exponent. Moments below one, if desired, follow from the L^1 bound.□

Proof. Fix an output coordinate i and a smooth compactly supported scalar function \phi. The transported test \phi\circ T_\rho has zero material derivative: it is constant along each probability-transport path. Differentiate its product with M_{\rho,i} under the transported laws. Lemma 4.1 gives \mathcal DM_\rho=\nabla_{\xi_\rho}^*U_1, so integration by parts and the chain rule yield

Smooth transport on the compact time interval and Gaussian tails justify differentiation; alternatively the test can first be cut off and the resulting identities passed to the limit. Push the expectation forward to y\sim\xi_*. Since U_1 is symmetric, its coefficient of \partial_j\phi(y) is

The terminal expectation is \E_{\xi_*}[M_{*,i}\phi]. At zero, M_0=u and (T_0)_\#\xi_0=\xi_*, so the initial expectation is u_i\E_{\xi_*}\phi. Integrating the displayed derivative in \rho proves

This proves the adjoint-divergence identity weakly, and smoothness makes it a pointwise identity.

For the bounds, \|B_\rho\|_{\op}\le C and \|U_1\|_{\op}\le CL immediately give \|K\|_{\op}\le CL. For a positive-order spatial derivative, use Lemma 5.1 on both DT_\rho and Y_\rho. An order-j derivative of DT_\rho has all split norms at most C_jlR^{-j}, and an order-j derivative of U_1 has all split norms at most C_jLR^{-j} by (3.4). When these fields are composed with Y_\rho, attach each chain-rule tensor by a single edge as in the preceding lemma. First derivatives of Y_\rho are bounded, and a derivative of order i>1 costs C_i lR^{-(i-1)}. Every order-j term in the product has total cost at most C_jLR^{-j}, since l\le1. Each factor retains free slots, so the single-contraction bound preserves every nontrivial split. Integrating over an interval of length at most one proves the assertion. All orders are smooth by Gaussian conditional integration and the smooth transport on \rho\le b_*<1; this uses no higher derivatives of the original unsmoothed gradient.□

Proof. The vector field v in (5.3) is smooth and hence locally Lipschitz. The bound on K and the global Lipschitz bound on M_* imply |v(y,g)|\le C_{x,r,s}(1+|y|+|g|). Gronwall’s inequality therefore prevents a solution from escaping to infinity in finite positive or negative time. Local uniqueness gives a complete smooth flow \Phi_t, with inverse \Phi_{-t}. We only need local Lipschitz continuity here; the factor G does not in general give a globally Lipschitz joint vector field.

Write the product density as

Differentiating in the position and momentum variables separately, with the row convention of the kernel lemma, gives

Thus \operatorname{div}(pv)=0. The Jacobian formula for the complete flow now yields, pointwise,

The determinant starts at one and stays positive. Change of variables therefore proves that \Phi_t preserves the probability density p.

For later differentiation it is useful to record the same cancellation as an operator identity. Relative to the product law, the Lie derivative is

Indeed the second-order part vanishes by skew symmetry. The first-order coefficients are -K^Tg/s in position and \nabla_{\xi_*}^*K/s=(M_*-u)/s in momentum, exactly as in (5.3).

Finally, integrating the momentum equation over [0,1] gives the pathwise identity

Stationarity makes G_1 a standard Gaussian, proving (5.4).□

Proof. Put \phi_k=\mathcal L^k\phi_0. We prove, for j\ge0 and k\ge1, the stronger estimate

For k=0 the same bound is assumed only for j\ge1, which suffices: the iteration differentiates \phi_{k-1} before multiplying by A. All L^2 norms below are for the invariant density.

Set V=A\nabla\phi_{k-1}. Leibniz’s rule and the induction hypothesis give, for each h\ge0,

To see why no dimension factor appears, in each coefficient tensor \nabla^vA isolate the slot contracted against the gradient of \phi_{k-1}. Its one-slot operator norm multiplies the Frobenius norm of the other factor. Summing the finitely many Leibniz terms preserves the displayed estimate.

Now commute j derivatives past the adjoint divergence in \phi_k=\nabla^*V. The main term is \nabla^*(\nabla^jV). The Frobenius estimate (4.7) bounds its L^2 norm by a constant times \|\nabla^jV\|_{F,2}+\|\nabla^{j+1}V\|_{F,2}. These two terms have powers R^{-(j+2k-1)} and R^{-(j+2k)}, respectively. Each remaining commutator contracts an order-v derivative of the score, 1\le v\le j, with \nabla^{j-v}V. The coefficient’s one-slot norm again suffices; the total exponent of R^{-1} is

All these terms are bounded by (5.7) because 0<R\le1. This closes the induction, including j=0 for k\ge1. Gaussian tails and polynomial growth justify the integrations by parts through cutoffs.

Along an autonomous trajectory, the kth time derivative of \phi_0 is \phi_k evaluated at the current state. Stationarity identifies its RMS norm at every time with the fixed-measure norm in (5.7) for j=0.□

Proof. Write Y_\rho for the exact path from the given random input and Q_\rho=M_\rho(Y_\rho). The marginal identity in Proposition 3.3 and (4.8) give

On a cell I=[a_c,b_c] of length \ell_c, let P be the Taylor polynomial of Q of degree m-1 at a_c. Minkowski’s inequality applied to the integral remainder gives, at each v\in I,

The last step uses \ell_c\le Ch_0(1-z) on that cell and 1-z^2\ge1-z. The same estimate holds at the interpolation nodes. Because interpolation reproduces P and its basis has bounded absolute sum, \|Q_v-\mathcal I_cQ(v)\|_2 obeys the same bound. Integrate this pointwise estimate over the cells, whose total length is at most one. This proves the required quadrature error, for partial integrals and for either orientation. No supremum of the random path has been taken.

Let e_k=\max_i\|y_i^{[k]}-Y_{\rho_i}\|_2. The global bound \Lip(M_\rho)\le CL, together with the absolute row sums, gives

On each exact transported marginal, conditional expectation and the primitive moment bound give \|M_\rho\|_2\le CD. Integrating the velocity shows e_0\le CrD. Since C_ml<1/2 for large d, summing the geometric recurrence proves (6.4). The norm here is the maximum of the nodewise RMS errors; it is not the RMS of their maximum.

For deterministic inputs the same recurrence, now in maximum node norm, gives a bounded Lipschitz constant for every iterate. In the correction y_i^{[k+1]}-y_{\mathrm{in}}, the factor r and the CL-Lipschitz mean give the stronger constant C_ml.□

Proof. For exact inputs (Y_\rho(y),G), harmonic motion \dot W=P, \dot P=-\nabla H_\rho(W) preserves \xi_\rho\otimes N(0,I_d). The vector field is globally Lipschitz, and its weighted divergence vanishes. Its Lie derivative is \nabla^*A\nabla with constant skew coefficient whose (W,P) block is -I. Lemma 5.4 therefore applies with \alpha=1, with a_0=l for T_\rho-\mathrm{id}, and with a_0=L for M_\rho. Indeed Lemma 5.1 and (3.4) give the required spatial bounds; changing a one-slot split bound to a Frobenius bound costs only \sqrt d. To apply the angle stencil, subtract the Taylor polynomial of z\mapsto(T_\rho-\mathrm{id})(W_z) of degree m-1 at zero. Its derivative at zero is exact, and its remainder at each stencil node has RMS norm at most C_ml\sqrt d\,\psi^mR^{-2m}. Multiplication by the absolute stencil sum C_m/\psi proves (6.9). This uses pointwise RMS remainders on the stationary harmonic path.

For the harmonic state equation, interpolate M_\rho(W_v) in (6.6). Since z\le(m-1)\psi, its Taylor remainder, multiplied by the prefactor r, is bounded by C_ml\sqrt d\,\psi^mR^{-2m} =C_ml\sqrt d\,\delta^{2tm}; the additional integration factors are bounded and can be discarded. The initial free-path error is at most CrD, and Picard contracts by C_ml. The computed harmonic state is globally Lipschitz with bounded constant in its two inputs, by the same row-sum argument as before.

Next compare the forward probability correction started at the exact harmonic state. That state has marginal law \xi_\rho, so Proposition 6.1 applies. Errors in its starting state then propagate through the approximate correction with factor C_ml, not a factor of order one. This also controls the backward-position error transmitted through the harmonic computation. Combining these estimates proves (6.10) and the asserted global Lipschitz bound.□

Proof. By Lemma 5.1 and independence of G, the RMS norm of the order-m derivative of C_\rho(y)G is at most \polylog(d)l\sqrt d(1-\rho^2)^{-m}. Here multiplication by the independent Gaussian converts the matrix norm to a Frobenius norm and costs at most \sqrt d.

On an ordinary cell, Taylor expansion and scaling the interpolation basis bound the derivative error by this derivative bound times C_m\ell_c^{m-1}. It is therefore at most C_m\polylog(d)l\sqrt d h_0^{m-1}/(1-\rho). After integration against 1/\rho, all ordinary cells contribute at most this prefactor times C(1+|\log R|+|\log h_0|), which fits the claimed more conservative h_0^{m-2} bound.

For the first cell, subtract the Taylor polynomial at zero of degree m-1. Its derivative at zero agrees with the exact zero derivative. The second derivative of the remainder, and that of the Hermite interpolant of the remainder, are bounded in RMS by C_m\polylog(d)l\sqrt d\,\ell_c^{m-2}. Both first derivatives vanish at zero. If E is the difference between the exact function and its Hermite interpolant on this cell, then

Consequently the integral of E'(\rho)/\rho over the first cell has RMS norm at most C_mB_dl\sqrt d\,\ell_c^{m-1}. This is within the stated, more conservative h_0^{m-2} bound.

On an ordinary cell the absolute coefficient sum of the interpolant derivative is at most C_m/\ell_c. Integrating and using \ell_c\ge c'h_0(1-\rho) gives a total bounded by

This gives the logarithms in (6.11). On the first cell, scale to the unit interval. The polynomial derivative vanishes at zero and is bounded by distance to zero times the second derivative, so the integral of its absolute value divided by that distance is bounded. Scaling back gives a weight sum C_m/\ell_c\le C_mh_0^{-1}. This proves the claim.□

Proof. Combine Lemmas 6.2 and 6.3, and divide by rs. The value errors are amplified only by C_m\psi^{-1} and the weight sum (6.11). The exact stencil truncation has factor \psi^{m-1}R^{-2m}=\delta^{2tm-4t} and the same probability-quadrature weight sum. Crucially,

Thus neither division introduces a negative power of the potentially small scale r. The resulting estimate is (6.14).

Each forward correction has Lipschitz constant C_ml in its inputs, and the backward and harmonic maps have bounded constants. Applying the same absolute weight sums and using l/(r\sigma)=L/\sigma gives (6.15); the momentum component has the smaller Lipschitz bound CL/\sigma.

For the last assertion, increase interpolation orders until the first line of (6.14) is of the desired order. Under (6.2), sufficiently many O(J+1) Picard layers make the second line of that order as well. The explicit common choice given below verifies independence from interpolation orders and from w.□

Proof. First compare the exact flow to polynomial integration of its exact velocity. Lemma 5.4, applied to the state with one additional derivative, bounds this quadrature error by

The approximation error of \widehat v at every exact trajectory node is bounded by Proposition 6.4, by stationarity. Previous state errors propagate with its global Lipschitz constant and the absolute row sums. The true centering velocity need not be globally Lipschitz in the joint variables; it is the approximate velocity’s bound that is used in this recurrence.

By (6.1), b\ge1/2 and w<t give b-4t-w>1/3. For sufficiently large d, fixed constants and polylogarithms in (6.15) can therefore be absorbed so that the recurrence contracts by at most \delta^{1/3}. Its initial error is at most C(L/s)\sqrt d by (5.8). Choose the quadrature and velocity error orders sufficiently high, and then O(J+1) Picard layers yield the state error claimed in the statement.

For the output, apply Lemma 5.4 to M_* with a_0=L. Its path quadrature error is at most C_mL\sqrt d\,((L/s)R^{-2})^m, which can be made at most \polylog(d)\sigma\delta^JD. The additional error from evaluating on the approximate positions is bounded by their RMS error times C_mL, using the CL-Lipschitz bound for M_*. Since L/\sigma is small, this is within the same budget.□

Proof. For large d, probability and harmonic Picard contractions C_ml\le\polylog(d)\delta^b can be bounded by \delta^{1/3}. The only additional inverse powers in kernel evaluation are \delta^{-4t-w}, besides the divisions already canceled in Proposition 6.4. The centering iteration also contracts by \delta^{1/3}. The probability iteration remainder is at most C\delta^{N/3}D; after the kernel weights the iteration remainder is at most C\delta^{N/3-4t-w}D. Both exponents exceed J with a fixed positive margin. For the truncation terms it suffices to check the exponents

The first four come respectively from probability quadrature, the zero-endpoint rule, its amplified probability-value error, and the harmonic and angle errors. For the last one use (L/s)R^{-2}\le\delta^{1/3} for large d in both outer quadratures, dividing the mean-output error by \sigma and using L/\sigma\le\delta^{1/3}. Every displayed exponent exceeds J by (6.17). The outer Picard remainder contracts by the same \delta^{1/3} factor. Thus fixed constants and polylogarithms, including those depending on m, fit within the strict exponent margins as d\to\infty.

Within a probability or harmonic Picard layer, all nodes depend on only the preceding layer, regardless of the number of nodes. A path through one kernel evaluation passes through a backward probability computation, a harmonic computation, and a forward probability computation, each having at most N layers. Such a path traverses at most N layers of the outer centering iteration. This gives depth O((N+1)^2); final quadratures add a bounded number of layers. All weights integrate known polynomials, or polynomials against sine or 1/\rho with the stated zero-endpoint condition. They require no oracle information and are allowable real constants in the model.□

Proof. We expand only the additive state formulas, stopping whenever an earlier mean value is reached. This convention is what makes the coefficient bounds local: the center expression inside that earlier mean is kept as its own node.

The three transport formulas. A probability Picard state has its input position as direct term and mean coefficients with absolute sum at most C_mr. A harmonic state has direct term \cos z times its input position plus \sin z times its input momentum, with mean coefficients of absolute sum at most C_mr. The sine-weighted row sums are bounded because their integration intervals have bounded length. These assertions hold at every iteration: its additive direct term is the original input, not the preceding Picard iterate.

For a kernel evaluation, the scalar expansion (6.13) has absolute coefficient sum at most B_d\sigma^{-1}\delta^{-4t-w}, by the derivative-quadrature, angle-stencil, and forward-integral weight bounds. The momentum formula (M_*-M_0)/s has coefficient sum at most C/\sigma. A centering Picard update has the original anchor as direct term and bounded outer integral rows. Its position and momentum corrections thus have the same respective coefficient bounds.

Centers and their coefficients. Consider a mean evaluation inside any of these formulas. Its physical center is x+r\rho v when its normalized position is v, by (3.1). Substitute its direct starting states: a backward probability state starts at the current centering position; a harmonic state starts at that backward position and the current momentum; and a forward probability state starts at the harmonic position. Stop each substitution at the mean terms just described. Only the free terms pass through this chain. Thus their anchor coefficients are exactly (\rho,0) or (\rho\cos z,\rho\sin z). For sampling there is only the first probability computation, so the coefficient is (\rho,0). Multiplication by r\rho, with \rho\le1, bounds the center coefficient sum in the mean case by

The term C_mr^2 accounts for the backward, harmonic, and forward integral corrections. In the sampling case the bound is simply C_mr^2. Enlarging these positive upper bounds to (6.21) is convenient for the implementation. The mean output is the single quadrature in (6.18), so its coefficient sum is at most C_m and its means all have probability time b_*. For sampling, dividing a probability correction by \rho_e\ge2/\sqrt5 gives output coefficient sum at most C_mr. These observations also give the stated final direct terms.

The scale at probability time \rho is r\sqrt{1-\rho^2}. In the mean computation \rho\le b_*, giving \tau\ge R, including equality for the output means. In the sampling computation,

This proves every assertion about the anchors and scales.

Expression size and depth. Let P\le B_d\delta^{-w} bound the number of nodes in any one probability mesh. Angle stencils, harmonic quadratures, and outer centering quadratures have fixed numbers of nodes once m is fixed. Formula (6.13) has at most C_mP^2 mean terms: one factor P comes from probability-time derivative quadrature, and the other from the forward transport integral. A direct centering position has the same bound, and its momentum has only a fixed number of terms. Substituting their free states into a backward, harmonic, or forward center adds at most one backward list, one harmonic list, and one forward list. No substitution continues through a mean value. Hence B_d\delta^{-2w} already bounds the direct expression size; the exponent c_1w=10w is a uniform conservative choice.

All Picard updates at a fixed layer depend only on earlier layers. Order an internal kernel computation as backward transport, harmonic motion, and forward transport, after its input states. A mean path within this computation adds at most 3N nodes above a path in those inputs. There are at most N outer Picard layers, and one final mean layer in (6.18). Thus 3N^2+1 bounds the depth. The sampling path has depth at most N. This gives an acyclic order for (6.20). Meshes, weights, angles and partial integral rows depend only on the prescribed scalar parameters; no exact mean is needed to form any of them.

Ideal-law errors. Couple the numerical mean formula and the exact centering path integral by the same ideal anchors. Corollary 6.6 bounds their RMS difference by B_d\sigma\delta^JD(x), while Proposition 5.3 identifies the exact output as N(u,\sigma^2I_d/4). This coupling proves (6.22). For sampling, the exact probability transport from a standard Gaussian has law \xi_{\rho_e}. After division by \rho_e this is the law of

The same-anchor RMS estimate and \rho_e\ge2/\sqrt5 prove (6.23).□

Proof. The eigenvalues of P lie in [\sqrt3/2,1]. The translation identity rewrites the first output as

We suppress the tensor product with I_d on slot-space matrices. The seed array in this expression is centered Gaussian and has slot covariance

Its slots are therefore independent standard Gaussian vectors. They may remain correlated with Z; that is harmless because the rewritten call is at the deterministic center x.□

Proof. Before 4K_0+1 decreasing transitions, the total loss in b is less than 10t(4K_0+1)<1/4, so b>3/4. At that many transitions, p\le K_0-(4K_0+1)/4<0, forcing termination. Therefore no branch can have more such transitions. Positive transitions increase b by t, while a nonterminal satisfies b<p\le K_0; their total number is consequently bounded in terms of K_0,t and the already bounded total decrease in b. A neutral transition is an \mathsf S-to-\mathsf M transition. Another neutral transition requires first reaching another \mathsf S, which can occur only through a decreasing transition. This bounds the neutral transitions as well and proves both claims.□

Proof. Let D_0 be the dependency-depth bound in Theorem 6.7, and let F_0 be its direct fan-out bound. Expanding a circuit into an expression tree costs at most (D_0+1)(1+F_0)^{D_0} nodes. Its size is thus \polylog(d)\delta^{-C_2w}, where C_2 depends on the Picard counts but not on w or the interpolation orders. This bounds its direct child invocations and its own anchor and residual slots. Expanding through the height H from Lemma 7.2 gives \polylog(d)\delta^{-C_3w} slots and terminal queries. Interpolation orders change only fixed multiplicative constants in these bounds, not C_2 or C_3.

The terminal shift has norm 1. Assuming the shift bound for the children, the bounded ratios and inverse norms above show that a parent’s shift norm is at most a fixed constant times one plus the sum of its direct children’s shift norms, including the anchor call. The extra reserved-slot shift has bounded coefficient. The triangle inequality is sufficient even if ordinary and extra anchor shifts affect the same slot. Iterating through the fixed routine height gives the same bound \polylog(d)\delta^{-C_3w}, after increasing C_3 if necessary.

Choose w>0 with

and only then choose the fixed interpolation orders needed in Theorem 6.7. For sufficiently large d, the last shift bound is smaller than D_s=\delta^{-t}. This is a valid bottom-up induction: it proves each child’s bound before using its inverse transformation in the parent. The count bound is also at most \polylog(d)\delta^{-t}.

To make the parameter order explicit, first fix the Picard counts, hence D_0,H,C_3, then w, then the interpolation orders. The universal scalar D_s fixes every center-noise budget. The circuit shapes, labels, requested noises, and slot counts are determined by these scalar parameters before any shifts are computed. Compute child shifts before their matrices and the parent shift. No center value, oracle response, numerical-error estimate, or assumption (7.2) enters this finite construction or the count. Only terminal routines query g, once each; hence the bound is a cap on every execution, not an expected query count.□

Proof. For the standard Gaussian seed Z, the primitive moment bounds give

Since L/\sigma\le\polylog(d)\delta^b and p\le b, this is within the mean error budget; couple the added Gaussian identically.

Let z_0 be the mode of \nu_{x,r}, so z_0=-rg(x+rz_0). Then

Also W_2(N(z_0,I_d),\nu_{x,r})\le Cl\sqrt d. Apply the synchronous comparison in Lemma 2.2 with Gaussian potential H_0(z)=|z-z_0|^2/2 (contraction constant 1) and tilted potential H(z)=|z|^2/2+F(x+rz). Their gradient difference is r[g(x+rz)-g(x+rz_0)], whose L^2(\nu_{x,r}) norm is at most l\|Z-z_0\|_{L^2(\nu_{x,r})}\le Cl\sqrt d. This proves the displayed Wasserstein bound under exactly the C^2 and bounded-Hessian hypotheses of that lemma. Adding the same independent comparison noises proves the sampling error ClD(x)\le\polylog(d)\eta\delta^bD(x), which is at most the requested bound because \eta\le1 and p\le b. Omitting G_2 proves the pre-reserved assertion. The Lipschitz bounds follow directly from the two displayed terminal formulas.□

Proof. For a final term of an \mathsf M parent, the output uses only means of scale rR, so r_i=rR, while \sigma_i is a fixed positive multiple of \sigma. Thus L_i/\sigma_i gains a factor R and l_i gains R^2, allowing labels (b+t,p). For a final term of an \mathsf S parent, \sigma_i is a fixed positive multiple of \eta/r. Therefore L_i/\sigma_i\le Cl/\eta and l_i\le l, allowing the unchanged labels.

For a nonfinal term, let r\tau_{\rm prev} be the scale of the node whose center is currently being supplied. Its expression is evaluated with budget r\tau_{\rm prev}/(2D_s), irrespective of the tier at which it occurs. Consequently

The coefficient estimates give \lambda A\le\polylog(d)\delta^{b-4t-c_1w} for an \mathsf M parent. For an \mathsf S parent, \lambda A=C_1l\le\polylog(d)\delta^b, since \eta\le1. Using D_s/R=\delta^{-2t} and c_1w<t costs less than the allowed 10t decrease in b. Also l_i\le l. For the position-anchor \mathsf S call of \mathsf M, l/\eta\le l/R, which has the same allowance. This proves the size conditions by the bounded routine height. All accumulated fixed and polylogarithmic factors are still one polylogarithmic factor.

As b>1/2 and 4t+c_1w<5t<1/6, the coefficient bounds imply (7.9). Finally, the positive declared bound for A includes a fixed multiple of r/\sigma in an \mathsf M parent, so \sigma_i\le r/(4D_sA)\le C\sigma. For \mathsf S, A=C_1r^2 gives r\sigma_i\le1/(4C_1D_s)\le C.□

Proof. Induct on routine height, with the terminal case supplied by Lemma 7.4. First hold the anchors fixed. A direct center term has Lipschitz constant 1 in x and at most Cr in a position anchor. The center-supplying expression combines child outputs with total absolute coefficient at most A. Using the child \mathsf M bound, the contribution of earlier supplied-center sensitivities is at most \polylog(d)\lambda A times their maximum. By (7.9), this is small. Induction through the bounded expression depth therefore bounds center sensitivities by a constant times their direct sensitivities, with an allowed polylogarithmic factor.

At the final expression the coefficient sum is A_f. Its direct term has no center dependence when anchors are held fixed; for \mathsf M it also has no position-anchor dependence. Hence its center sensitivity is at most \polylog(d)A_f\lambda, and its \mathsf M position-anchor sensitivity at most \polylog(d)A_f\lambda r. These are the asserted bounds with anchors treated as separate inputs. The \mathsf S Gaussian anchor is indeed center independent. The \mathsf M position anchor is supplied by \mathsf S at the same center, whose center Lipschitz constant by induction is \polylog(d)r\lambda. Its additional contribution is at most \polylog(d)Lr\lambda=\polylog(d)L^2, and L^2\le\lambda for r\le1 and \lambda\le1. All seed transformations depend only on scalar parameters, not on x. The proof uses only fixed-seed identities and Lipschitz bounds, and therefore requires no seed independence.□

Proof. The terminal assertion is Lemma 7.4. All random variables below have finite second moments: the actual algorithms are finite compositions of affine maps and globally Lipschitz oracle evaluations, while ideal tilted variables have the primitive moment bounds. Thus all W_2 comparisons are defined.

For a nonterminal \mathsf M, first couple its approximate \mathsf S position anchor to an ideal anchor, independently of all remaining seed blocks. The child labels are (b-10t,p-1/4). By induction and Lemma 7.6, replacing this anchor costs at most

After division by the mean error unit \sigma\delta^pD(x), the cost is at most \polylog(d)\delta^{b-1/4}, hence is bounded for large d. We may now assume that the position anchor has law b_*Z+RN, where Z\sim\nu_{x,r} and N is standard Gaussian, and that the momentum anchor is an independent standard Gaussian. These are exactly the input laws in Theorem 6.7.

Let x_i be the centers of the deterministic exact-mean circuit under these anchors, and put D_i=\sqrt d+|g(x_i)|. Their RMS values over anchors obey

Indeed, both position and momentum anchors have RMS at most CD(x). Use |g(x_i)|\le|g(x)|+\lambda|x_i-x| and the circuit’s direct terms. An exact primitive mean at a preceding center has magnitude at most CD_j. Their contribution to this maximum-of-RMS recurrence is at most C\lambda A times its previous maximum. This is small by (7.9); the direct terms cost at most CD(x). Induction through the node ordering proves (7.10) without a maximum inside the expectation.

Condition on the anchors. For an expression E, write \epsilon_E for the W_2 error of \operatorname{Eval}(E,n) relative to its deterministic exact value plus N(0,n^2I_d). Couple the supplied center of its ith term to x_i+n_iZ with its conditional W_2 error \epsilon_{E_i}. Use the same transformed child seeds on these two centers, drawn independently of this coupling. The fixed-seed Lipschitz bound costs at most \polylog(d)\lambda\epsilon_{E_i}. At the ideal noisy center, Lemma 7.1 identifies the child’s law with its standalone law at x_i. The induction hypothesis therefore supplies the additional error \polylog(d)\sigma_i\delta^{p_i}D_i.

Different term occurrences have disjoint blocks, so their outputs are independent conditional on the anchors. Product couplings compare them to independent noised exact means. Together with the independent residual in (7.7), the ideal sum has exactly covariance n^2I_d. We conclude that

Notice the order of this argument: inaccurate supplied centers use only pathwise stability, and standalone accuracy is invoked only after Gaussian absorption. No assertion is made that the actual algorithm’s intermediate states have the ideal input laws.

Take RMS in the anchors and expand (7.11) through the finite expression tree. At its root, direct-child labels have p_i=p and

Thus their intrinsic errors fit the requested budget. At tier j\ge2, absolute coefficient sums and Lipschitz factors contribute at most

All such children have the single accuracy-label loss p_i=p-1/4, not one loss per expression tier. Their noise units are at most C\sigma in an \mathsf M parent and C/r in an \mathsf S parent by Lemma 7.5; these units combine with A_f=C_1 or C_1r, respectively. Already the first factor \polylog(d)\lambda A\le\polylog(d)\delta^{1/3} pays for the loss \delta^{-1/4}, leaving a positive power \delta^{1/12}. The bounded expression depth allows all further tiers and polylogarithmic factors to be summed. Leaf affine expressions have zero error. This proves the required final evaluation error.

Conditionally on the anchors, the ideal final evaluation is its deterministic exact-mean circuit value plus a Gaussian of the specified fixed covariance. Its noise may therefore be taken independent of the anchors. By Lemma 7.5, Theorem 6.7 applies. It compares the ideal mean circuit to N(u,\sigma^2I_d/4) and the ideal sampling circuit to \nu_{x,r}*N(0,\eta^2I_d/4), with the respective errors \polylog(d)\sigma\delta^JD(x) and \polylog(d)\delta^JD(x). Since J=K_0+2 and p\le K_0, both errors fit the requested bounds. Couple the independent evaluation noise identically in each comparison. For \mathsf M, its variance 3\sigma^2/4 raises the Gaussian target variance from \sigma^2/4 to \sigma^2. For \mathsf S, its variance \eta^2/2 raises the Gaussian convolution variance from \eta^2/4 to 3\eta^2/4. This closes the standalone induction. Adding the common reserved noise preserves the W_2 bound; Gaussian smoothing by its standard deviation \eta/2 gives the stated total variation bound by (2.7).□

Proof of Theorem 2.1. Use the top \mathsf S call with (b,p,\eta)=(1,K_0,1). Because l=\delta r^2\le\delta, its size condition holds with constant 1. Proposition 7.7 gives both W_2 and total variation error at most \polylog(d)\delta^{K_0}D(x) from \nu_{x,r}*N(0,I_d), uniformly in x and 0<r\le1. Lemma 7.3 gives at most \polylog(d)\delta^{-t} evaluations of g on every execution. The output has the form X_{\rm pre}+G/2, where the reserved standard Gaussian G is independent of X_{\rm pre} and all queries, and the pre-reserved conclusion of Proposition 7.7 gives W_2(\mathcal L(X_{\rm pre}),\nu_{x,r}*N(0,3I_d/4)) \le\polylog(d)\delta^{K_0}D(x). The exact means and numerical comparison laws serve only to analyze this finite algorithm; its only oracle operations are the terminal evaluations.□

Proof. Under the joint law \pi_V(dx)N(x,hI_d)(dy), the conditional law of x given y is \kappa_y. Drawing this conditional and then redrawing y given x preserves the joint law, whose y-marginal is \widetilde\pi.

The potential in (8.1) has curvature at least 1+1/h. Replacing y_1 by y_2 changes its gradient by the constant vector (y_1-y_2)/h. The synchronous comparison (2.5) therefore bounds the conditional \Wtwo distance by |y_1-y_2|/(1+h). Couple the subsequent Gaussian increments identically. Integrating these conditional couplings against an input coupling, as in Lemma 2.3, proves the kernel contraction.□

Proof. Lemma 2.2, applied at the mode zero, gives \E_{\pi_V}|X|^2\le d. Hence the root mean square under \widetilde\pi is at most C\sqrt d. Let \mu_k be the actual law after k steps, and set a_k=\Wtwo(\mu_k,\widetilde\pi). Then

All these moments are finite. Indeed they are finite at k=0, and the local \Wtwo bound in (8.3), together with the conditional moment bounds for P_h, propagates finite second moments to the next step.

Take the root mean square of the first local error in (8.3). Conditional coupling and Lemma 8.1 give

For all large d, the coefficient of a_k is at most 1-ch for an absolute c>0. Since a_0\le C\sqrt d, summing this geometric recurrence yields

The second term is o(\sqrt h): after division by \sqrt h, it is bounded by a constant times B_d\delta^{K_0-1}\sqrt d, which tends to zero because e(K_0-1)>1/2. At k=\lceil Ch^{-1}\log d\rceil, a sufficiently large fixed C makes the first term o(\sqrt h) as well.

Use one more transition to obtain TV control. The mixture bound in Lemma 2.3 gives

For the remaining comparison, the conditional x'-laws from \mu_k and \widetilde\pi can be coupled with root mean square at most a_k/(1+h). Their common independent \sqrt h-Gaussian then gives, by (2.7),

Invariance of \widetilde\pi completes the proof.□

Proof. The initial estimate follows from |g_y(0)|\le2\sqrt h\,|y| and the moment bound for \widetilde\pi. Conditional on an ideal history, the next observation has the law of Z+N at center x_j and scale r_j=2^{-j/2}. Equation (2.6) therefore gives

The center update and \Lip(g_y)\le\delta imply, by the conditional and then unconditional L^2 triangle inequalities,

Both \sum_{j\ge0}r_j and \sum_{j\ge0}r_j^2 are bounded by absolute constants. In particular the products of 1+C\delta r_j^2 are bounded by \exp(C\delta\sum_jr_j^2). Expanding the recurrence proves the stated uniform bound.□

Proof. Translate by z_0. Up to an additive constant, the potential of w=z-z_0 is |w|^2/2+J(w), where

The mode equation gives \nabla J(0)=0, and \Lip(\nabla J)\le l. In particular J(0)=0 and |J(w)|\le l|w|^2/2. For 0\le s\le1, the potential |w|^2/2+sJ(w) has mode zero and Hessian between (1-l)I_d and (1+l)I_d. Its law \mu_s satisfies, by Lemma 2.2,

The same quadratic bounds justify differentiating

twice under the integral, and give \Phi''(s)=\Var_{\mu_s}J. Since \mu_0=N(0,I_d),

Pinsker’s inequality, followed by translation back by z_0, proves the assertion.□

Proof. Given \varepsilon>0, first choose e\in(0,1) with 2e<\varepsilon, then an integer K_0>2+1/e, and then t>0 satisfying (2.3). Fix the auxiliary construction’s orders in its prescribed order, and finally the constants in the chain, descent depth, and base iteration count. These choices may depend on \varepsilon, but they are fixed before observing the oracle or the random seeds.

Every invocation of Theorem 2.1 uses at most \polylog(d)\delta^{-t} evaluations of the gradient oracle g_y. By (8.2), each evaluation is one original query, including calls at random descent centers. The cap is uniform in those centers. There are O(\delta^{-1}\log d) proximal calls, O(\log d) descent calls, and O(\log d) additional gradient queries in the base rule. Summing gives the deterministic cap

The query counts and recursion shapes of this construction are fixed by the scalar parameters; they do not depend on oracle magnitudes or approximation events. All maps are measurable by the primitive theorem and the explicit finite operations above. Since t<1, we have e(1+t)<2e<\varepsilon. The fixed polylogarithmic factors fit C_\varepsilon d^\varepsilon for all sufficiently large d. Equation (8.10) gives TV error at most 1/10 above a uniform dimension threshold.

For each of the finitely many smaller dimensions, use independent standard Gaussian proposals. Accept a proposal X with probability

Integrating the Hessian bounds from the known minimizer gives |x|^2/2\le V(x)\le|x|^2. The acceptance probability is therefore at most one, and its average is at least \E e^{-|X|^2/2}=2^{-d/2}. The proposal density times its acceptance probability is proportional to e^{-V(x)}, so an accepted proposal has law \pi_V.

Cap the number of trials at \lceil4\cdot2^{d/2}\rceil. Return the first accepted proposal, or zero if none is accepted. The failure probability is at most e^{-4}<1/10; conditional on success, the returned proposal has law \pi_V. The fallback thus has TV error at most e^{-4} and makes at most the prescribed number of original first-order queries on every execution, one per trial. Increasing C_\varepsilon covers these finitely many caps. This proves the assertion for all d\ge2.□

Proof. Include all the algorithm’s randomness, including any randomness used at output, in a private seed independent of O. Fix this seed for the moment. Pad executions with zero queries so that there are exactly q query slots.

We use the following conditional property of Haar measure. For a fixed restriction T:E\to W of an orthogonal map, the conditional law of its remaining restriction is uniform among isometries from E^\perp to W^\perp. For a specified unit vector in either complement, its partner is uniform on the unit sphere of the other complement. After revealing that pair, the map between the smaller complements is again uniform. This follows from invariance under the orthogonal transformations fixing E and W pointwise. Applying this rule successively also permits each specified direction to be a measurable function of the already revealed data.

All conditional laws below are laws of finite-dimensional data after the seed has been fixed. To represent the subspaces and partial isometries measurably, apply Gram–Schmidt to the revealed vectors in their order of appearance, discard zero residuals, and complete a basis by scanning the standard coordinate vectors in order. These operations are Borel on each rank stratum. Thus regular conditional laws and the successive Gaussian draws can be chosen as measurable kernels; no continuity of the algorithm’s query maps is required.

Initially E=W=\{0\}. After each reply, reveal the restriction T=O|_E, where E is the span of the eigenframe coordinates of all queries and replies so far, and W=OE is their span in the original coordinates. Our conditional invariant is that the unrevealed map is Haar on the two complements, given the full latent history: the seed, the transcript, this restriction, and all auxiliary variables constructed so far.

Let x be the next query, which is fixed given that history. Write x=x_W+x_\perp according to W\oplus W^\perp. The preimage of x_W is known. If x_\perp\ne0, put

Conditionally on the history, v is uniform on the complementary unit sphere. Let h be a standard Gaussian vector in E, and let a have the distribution of the Euclidean norm of a standard Gaussian in \R^{d-\dim E}. Draw h and a independently of one another and of the unrevealed map, conditional on the history, and define

The Gaussian polar decomposition shows that g_i is a standard Gaussian in \R^d conditional on the entire previous history. Moreover, a>0 almost surely, so v=(g_i-h)/a\in E+\operatorname{span}\{g_i\}. Reveal Ov=w and extend T to E+\operatorname{span}\{v\}. Conditioning on the full vector g_i does not reveal any further part of O: its parallel component h and complementary radius a are independent auxiliary variables, and its normalized complementary component is precisely v. Thus the conditional Haar invariant persists. If x_\perp=0, instead draw a fresh independent standard Gaussian g_i and make no extension at this step.

In either case the query preimage u=O^Tx is now known and lies in the old space E+\operatorname{span}\{g_i\}. Its reply has preimage \Lambda u. Project \Lambda u onto the enlarged eigenframe space. The image of this projection is known. If the residual is nonzero, reveal the image under O of its unit direction and extend the isometry once more. This is the other form of the same Haar-revelation rule: the direction in the eigenframe is fixed by the current history, and its image is uniform in the unrevealed complementary sphere. The resulting reply is exactly Ax, and the conditional invariant is preserved.

Every g_i has conditional law N(0,I_d) given all earlier Gaussian vectors and the earlier history. Consequently these vectors are independent. To check the polynomial degrees, let E_i be the revealed eigenframe space after reply i. Induction gives

Indeed, the new query preimage is in E_{i-1}+\operatorname{span}\{g_i\}, and adjoining its reply applies \Lambda once. All old directions remain in the indicated larger span.

Finally treat the algorithm’s output as one more query direction, with no reply: repeat the preimage construction using g_{q+1}. This gives O^TX\in E_q+\operatorname{span}\{g_{q+1}\}, proving (9.1). At most two dimensions were exposed per oracle query and one for the output, so the assumed dimension inequality ensures that every required complement is available. The construction started with the genuine Haar input and never altered it. Its added variables are obtained by measurable kernels from the original variables and independent auxiliary draws. The construction therefore applies jointly measurably for the fixed-seed algorithms; integrating over the original independent seed proves the assertion for randomized algorithms.□

Proof. Let m=(q+1)^2 and form a d\times m matrix B, whose columns are indexed by (i,j) with 1\le i\le q+1 and 0\le j\le q, by

The functions \cos(j\theta) are degree-j polynomials in \cos\theta with nonzero leading coefficients. Since \cos\theta_k=2\lambda_k-3, the range of B is precisely \mathcal K_q.

Discrete cosine orthogonality gives

For completeness, the equally spaced half-grid satisfies

Here the finite cosine sum follows by summing a geometric progression, and its denominator is nonzero in the indicated range. Expand products of two cosines for the first identity in (9.4). For the second, expand the product of three cosines with frequencies n,j,j'. Its frequencies are nonzero because n>j+j', and have absolute value at most 4q+1<2d. Different Gaussian blocks have zero cross means in both cases.

An entry of either Gram matrix is a sum of independent row terms. The deterministic coefficient in each term has absolute value at most 2/d, and the variance of a product of two standard Gaussian entries is at most 2. Thus every entry has variance at most 8/d, and

By Markov’s inequality,

Since \|C\|_{\op}\le\|C\|_{\HS}, outside this event both operator norms are at most 40m/\sqrt d=\tau_d. Whenever these bounds hold, for every a\in\R^m,

These inequalities imply (9.3) on the entire range of B.□

Proof. Independence of the coordinates gives \mu_{d,q}=\sum_k p_k/\lambda_k. Substitute (9.6). The constant term vanishes by (9.5). For 1\le j<2d-n, use

If j\ne n, both frequencies have nonzero absolute value below 2d, so their sums vanish. If j=n, the constant term contributes d/2 and the frequency 2n<2d contributes zero. Thus only j=n contributes in this range. Each remaining cosine product has absolute sum at most d, so summing the geometric tail proves (9.7). Finally, a centered normal variable with variance v satisfies \Var(U_k^2)=2v^2. Since |p_k|\le1 and \lambda_k\ge1, summing these coordinate variances proves (9.8).□

Proof of Theorem 9.1. Take A=O\Lambda O^T with O Haar and \Lambda as in (9.2). The reduction at the start of the section identifies an exact value-and-gradient algorithm on V_A(x)=x^TAx/2 with an algorithm querying x\mapsto Ax, with the same query cap and output law.

Fix any algorithm using at most q queries on every execution, where

All dimension restrictions in Lemmas 9.2, 9.3, and 9.4 hold for sufficiently large d, uniformly over these q. Let \P_{\mathrm{alg}} be the joint law of the Haar rotation and the algorithm’s output on this input. Define the event on pairs (O,x)

This event is measurable. It is also determined by (A,x): if T_n is the polynomial defined by T_n(\cos\theta)=\cos(n\theta), then OD_pO^T=T_n(2A-3I). We retain O in the notation to use eigenframe coordinates throughout the comparison.

Apply Lemma 9.2 to the algorithm’s output and then Lemma 9.3 to its auxiliary Gaussian vectors. On the latter lemma’s event, the norm condition in (9.9) forces the quadratic form to be at most 2\tau_d(2d)=160(q+1)^2\sqrt d. Hence

For the target joint law \P_{\mathrm{target}}, sample O from the same prior and then X from N(0,O\Lambda^{-1}O^T). Conditional on O, the vector U=O^TX has law N(0,\Lambda^{-1}). Our choice of c_0 gives

The remainder in (9.7) satisfies

uniformly, since n\le1+2c_0\log d. Thus it is negligible relative to the leading term. In particular, with T_{d,q}=160(q+1)^2\sqrt d, we have |\mu_{d,q}|>T_{d,q} for large d, and

uniformly, by Lemma 9.4 and Chebyshev’s inequality. Thus the probabilities of \mathcal E_{d,q} under the two joint experiments differ by at least 0.99-o(1).

To compare these joint laws with the required guarantee on each fixed input, let K_O denote the algorithm’s output law on A=O\Lambda O^T and let G_O=N(0,A^{-1}). For a measurable joint event \mathcal E, write \mathcal E_O=\{x:(O,x)\in\mathcal E\}. Both laws use the same Haar prior, so

The pointwise sampling guarantee would bound the last integral by 1/10. The event used to disprove that guarantee is a test of the joint laws; the algorithm need not know or compute it. This is impossible for (9.9). Therefore no such algorithm with q\le c_0\log d can meet the required accuracy for all inputs. Decreasing the constant if necessary proves the theorem.□