\documentclass[11pt]{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb,amsthm,mathtools}
\usepackage{microtype}
\usepackage{enumitem}
\usepackage{graphicx}
\usepackage[numbers,sort&compress]{natbib}
\usepackage{xurl}
\usepackage{xcolor}
\usepackage[colorlinks=true,linkcolor=blue!50!black,citecolor=blue!50!black,urlcolor=blue!50!black]{hyperref}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newcommand{\R}{\mathbb R}
\newcommand{\E}{\mathbb E}
\renewcommand{\P}{\mathbb P}
\newcommand{\TV}{\mathrm{TV}}
\newcommand{\Wtwo}{W_2}
\newcommand{\Lip}{\operatorname{Lip}}
\newcommand{\polylog}{\operatorname{polylog}}
\newcommand{\op}{\mathrm{op}}
\newcommand{\HS}{\mathrm F}
\newcommand{\tr}{\operatorname{tr}}
\newcommand{\Cov}{\operatorname{Cov}}
\newcommand{\Var}{\operatorname{Var}}
\newcommand{\KL}{\operatorname{KL}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\norm}[1]{\lVert#1\rVert}
\newcommand{\abs}[1]{\lvert#1\rvert}
\setlist[enumerate]{itemsep=3pt,topsep=5pt}
\setlength{\emergencystretch}{2em}
\begin{document}
The problem is to sample from a probability density proportional to
$e^{-V}$ on $\R^d$, under a uniform second-order bound and a known
minimizer, when the only access to $V$ is exact first-order queries.
For $d\ge2$ let $\mathcal V_d$ be the set of $C^2$ potentials
$V:\R^d\to\R$ with
\begin{equation}\label{eq:original-class}
V(0)=0,\qquad \nabla V(0)=0,\qquad
I_d\preceq\nabla^2V(x)\preceq2I_d\quad\text{for every }x\in\R^d .
\end{equation}
The origin is therefore the unique minimizer and is known to the
algorithm. The target law is
\[
\pi_V(dx)=Z_V^{-1}e^{-V(x)}\,dx,\qquad
Z_V=\int_{\R^d}e^{-V(x)}\,dx .
\]
\emph{Oracle model.} An exact first-order query at a point $x$
returns the pair $(V(x),\nabla V(x))$. Queries may be adaptive.
Random real numbers may be generated. Computation between queries is
unrestricted: only the number of queries is measured. An algorithm
may depend on $d$ but learns about $V$ only through oracle answers.
All algorithms are measurable.
Write $\norm{\mu-\eta}_{\TV}=\sup_B|\mu(B)-\eta(B)|$ for total
variation. Let $Q(d)$ be the smallest integer $q$ such that some
algorithm makes at most $q$ queries on \emph{every} execution and
returns $X$ with $\norm{\mathcal L(X)-\pi_V}_{\TV}\le1/10$ for every
$V\in\mathcal V_d$. The budget is a uniform cap, not an expectation.
Define the optimal dimension exponent
\[
\gamma_*=\inf\{\gamma\ge0:Q(d)\le d^{\gamma+o(1)}\}
=\limsup_{d\to\infty}\frac{\log Q(d)}{\log d}.
\]
\begin{theorem}[Main theorem]\label{thm:main}
In the oracle model above:
(a) for every $\varepsilon>0$ there is $C_\varepsilon<\infty$ with
$Q(d)\le C_\varepsilon d^\varepsilon$ for all $d\ge2$;
(b) there are absolute constants $c>0$, $d_0$ with $Q(d)\ge c\log d$
for every $d\ge d_0$, even when inputs are restricted to quadratic
potentials $V(x)=x^TAx/2$ with $I\preceq A\preceq2I$.
The lower bound holds for arbitrary randomized adaptive algorithms.
In particular, $\gamma_*=0$.
\end{theorem}
The upper bound (a) is Theorem~\ref{thm:upper}, the lower bound (b) is
Theorem~\ref{thm:lower}. Together they fix the dimension exponent but
not the finer growth of $Q(d)$: the upper bound does not assert a
polylogarithmic query count. The constants and the algorithmic
approximation orders may depend on $\varepsilon$. No arithmetic or
bit-complexity bound is imposed or deduced.
Familiar samplers -- Langevin, Hamiltonian Monte Carlo, and their
discretizations -- are governed by a step-size analysis whose error,
at fixed conditioning and accuracy, typically grows like a positive
power of $d$. Two recent lines of work break this pattern at fixed
accuracy: shifted composition (Altschuler, Chewi, and Zhang) and exact
diffusion simulation (Chen, Chewi, Rakhlin, and Zhang) achieve
$\widetilde O(d^{1/3})$ queries; smoothed Picard methods (Chen, Chewi,
Lu, and Zhang) achieve $\widetilde O(d^{1/6})$ expected gradient
queries and a Gaussian-cloud variant achieves $\widetilde O(d^{1/5})$
at high accuracy. All of these are positive powers. The question is whether \emph{any} positive power is necessary.
\url{https://arxiv.org/pdf/2506.23062}
\url{https://arxiv.org/pdf/2608.05022}
\url{https://arxiv.org/pdf/2609.06906}
\url{https://arxiv.org/pdf/2609.38710}
The answer is no: the exponent can be driven to zero. The cost of
this improvement is that the approximation order must grow with the
requested exponent $\varepsilon$; the algorithm is aimed at oracle
complexity, not practical running time. The lower bound says that,
nevertheless, some logarithmic dependence on $d$ is unavoidable.
Historical context for the individual tools is given where they are
used. The Gaussian conditional-sampling framework is due to Lee, Shen,
and Tian and to Chen, Chewi, Salim, and Wibisono; the smoothed Picard
precedent to Chen, Chewi, Lu, and Zhang; the all-split tensor techniques
build on the classical Appell polynomials of Anshelevich, the Poincar\'e
concentration bounds of Adamczak and Wolff, and the relaxed tensor norms
of Kothari and Steinhardt. The lower-bound strategy adapts the
block-Krylov method of Chewi, de Dios Pont, Li, Lu, and Narayanan.
\url{https://arxiv.org/pdf/2010.03106}
\url{https://proceedings.mlr.press/v178/chen22c/chen22c.pdf}
\url{https://arxiv.org/pdf/2609.06906}
\url{https://arxiv.org/pdf/1711.07465}
\url{https://arxiv.org/pdf/2304.02599}
The proof has five movements. Each is a construction plus an
estimate, and each estimates only one kind of error.
\paragraph{Movement 1: reduce to a near-quadratic problem.}
Condition a target sample $X\sim\pi_V$ on $X+\sqrt h\,G=y$, where $G$
is a standard Gaussian and $h>0$ is small. The conditional density is
proportional to $\exp(-V(x')-|x'-y|^2/(2h))$. Writing
$x'=y+\sqrt h\,z$, the normalized potential is $|z|^2/2+F_y(z)$ with
$F_y(z)=V(y+\sqrt h\,z)$, and $\Lip(\nabla F_y)\le2h$. Thus the
hard part is the following smaller, uniform problem. For a $C^2$
perturbation $F$ with $\Lip(\nabla F)\le\delta$ small, a center $x$,
and $0<r\le1$, sample or compute means of
\begin{equation}\label{eq:primitive-law-intro}
\nu_{x,r}(dz)\propto
e^{-|z|^2/2-F(x+rz)}\,dz .
\end{equation}
Both the sample and the conditional mean
$u_{x,r}=\E_{\nu_{x,r}}\nabla F(x+rZ)$ are constructed with an added
Gaussian, whose noise is later removed by conditioning.
\paragraph{Movement 2: Gaussian smoothing and dimension control.}
A Gaussian convolution makes the conditional mean smooth although
$F$ is only $C^2$. Its derivatives are cumulant tensors of the
conditional law. Two facts keep them in hand: a Poincar\'e
inequality controls every tensor split, and keeping at least two free
tensor slots until the very end ensures each additional derivative
costs at most a polylogarithmic factor in $d$, never a fresh power.
Only the final passage from a tensor to a vector costs one $\sqrt d$.
\paragraph{Movement 3: two flows.}
The probability transport
$Y_\rho=\rho Z+\sqrt{1-\rho^2}\,N$ moves a Gaussian to $\xi_\rho$.
Its velocity is a conditional mean, hence another instance of the
auxiliary problem. A complementary \emph{stationary} flow in position
and momentum, built from a skew-matrix drift, turns a path integral of
gradient means into an exact Gaussian centered at the desired mean.
High-order quadrature of both flows uses the tensor bounds of Movement
2 and a singularity-cancellation at the endpoint.
\paragraph{Movement 4: finite evaluation with Gaussian noise.}
Both flows become finite expressions in exact conditional means.
The difficulty: a mean is needed at a center that is itself built from
means. A prescribed Gaussian is added to every computed center and
absorbed, via an exact covariance identity on the random seeds, into
the lower-level routine. This restores the law of a call at the
intended center without changing any physical query. Every routine
carries an error label and a size label; the labels make the recursion
finite and cap the number of queries on \emph{every} execution.
\paragraph{Movement 5: return to the target and the lower bound.}
A contracting Gaussian proximal chain, run for a prescribed number of
steps, approaches $\pi_V*N(0,hI_d)$. The noised conditional sample is
then unnoised by variance-halving descent, with a mode-centered
Gaussian at the bottom. The query count of the descent and of the
chain, multiplied by the per-call cap, gives $d^{e(1+t)+o(1)}$, and
$t$ is chosen so $e(1+t)<\varepsilon$. The lower bound is independent:
a random rotation hides the input, and a block-Krylov statistic forces
$\log d$ queries for quadratic potentials.
The conditional density is a Gaussian times a factor with small
Hessian. A uniform guarantee is needed for that larger class; it is
stated below, together with the comparison and smoothing
facts used everywhere after. The center is arbitrary: later routines
are called at centers produced by earlier random computations, so no
gradient bound at the center may be assumed.
Vectors use Euclidean norm; matrices the induced operator norm unless
stated otherwise. For a random vector or tensor $X$,
$\norm{X}_2=(\E|X|^2)^{1/2}$ with the chosen pointwise norm, and
$\norm{X}_{\HS,2}$ specifies the Frobenius norm. $\Wtwo$ is quadratic
Wasserstein distance, and
$\norm{\mu-\eta}_{\TV}=\sup_B|\mu(B)-\eta(B)|$. Laws and oracle
procedures are measurable.
Fix $0<e<1$ and put $\delta=d^{-e}$. Let $F:\R^d\to\R$ be $C^2$ with
$g=\nabla F$ Lipschitz with constant at most $\lambda=\delta$; values
of $F$ need not be available. For $x\in\R^d$ and $0<r\le1$ define
\begin{equation}\label{eq:primitive-law}
\nu_{x,r}(dz)=
\frac{\exp(-|z|^2/2-F(x+rz))\,dz}
{\int_{\R^d}\exp(-|v|^2/2-F(x+rv))\,dv},\quad
f(z)=g(x+rz),\quad u=\E_{\nu_{x,r}}f .
\end{equation}
Often $\nu=\nu_{x,r}$, and write
\begin{equation}\label{eq:primitive-scales}
L=\lambda r,\qquad l=\lambda r^2=rL,\qquad
D(x)=\sqrt d+|g(x)| .
\end{equation}
The center $x$ need not be a minimizer. Since $\lambda$ is only a
declared bound, affine and constant $F$ are included. The full
potential in \eqref{eq:primitive-law} has Hessian between $(1-l)I$ and
$(1+l)I$; as $l\le\delta<1$ for $d\ge2$ the normalizing integral is
finite even when $F$ is not convex. For all quantitative work it suffices to take
$d$ large enough that $l\le1/2$.
Throughout, $\polylog(d)$ is a factor bounded by $C(1+\log d)^C$ whose
constants may depend on the fixed parameters ($e$, accuracy orders,
and fixed orders introduced below) but never on $F,x,r$ or the
answers. Different occurrences may use different constants. Every
threshold ``for large $d$'' is uniform.
\begin{theorem}[Noisy near-quadratic sampling]\label{thm:primitive}
Fix $0<e<1$, an integer $K_0\ge2$, and $t>0$ with
\begin{equation}\label{eq:t-choice}
40(K_0+2)t<\tfrac14 .
\end{equation}
There are constants and a threshold depending only on $e,K_0,t$ such
that for every integer $d$ above the threshold, with $\delta=d^{-e}$,
one measurable randomized oracle procedure works for every $C^2$
function $F$ with $\Lip(\nabla F)\le\delta$, every center $x$, and
every $0<r\le1$. Given access to $g=\nabla F$ it makes at most
$\polylog(d)\delta^{-t}$ evaluations of $g$ on every execution and
returns
\[
X=X_{\mathrm{pre}}+\tfrac12G,
\]
where $G\sim N(0,I_d)$ is independent of $X_{\mathrm{pre}}$ and of all
queries, and
\begin{equation}\label{eq:primitive-pre-noise}
\Wtwo\!\left(\mathcal L(X_{\mathrm{pre}}),
\nu_{x,r}*N(0,\tfrac34I_d)\right)
\le \polylog(d)\delta^{K_0}D(x).
\end{equation}
Consequently both
\[
\Wtwo\!\left(\mathcal L(X),\nu_{x,r}*N(0,I_d)\right)
\quad\text{and}\quad
\norm{\mathcal L(X)-\nu_{x,r}*N(0,I_d)}_{\TV}
\]
are at most $\polylog(d)\delta^{K_0}D(x)$. All constants are uniform
in $F,x,r$. The scalar parameters, computation shapes, and query cap
of the procedure do not depend on $x$ or on $F$, and it uses no
evaluations of $F$ itself.
\end{theorem}
The proof is carried out in the derivations below. An
independent Gaussian is reserved until all queries are made: adding the
same reserve to a coupling keeps its transportation error, and smoothing
by the reserve gives the total-variation conclusion through
Lemma~\ref{lem:coupling-smoothing}. The target carries a full
Gaussian convolution; that convolution is removed later.
The elementary diffusion facts used throughout are collected next, in the
form actually needed, with regularity requirements explicit. These are
familiar forms of the curvature method for functional inequalities
(Bakry, Gentil, and Ledoux); the argument is included for completeness.
\url{https://doi.org/10.1007/978-3-319-00227-9}
\begin{lemma}[Uniformly convex laws]\label{lem:convex-preliminaries}
Let $H\in C^2(\R^n)$ satisfy $mI\preceq\nabla^2H\preceq MI$ with
$0<m\le M<\infty$. Write $\mu(dx)\propto e^{-H(x)}dx$ and let $z_0$
be the unique mode. Then
\[
\E_\mu|X-z_0|^2\le n/m,\qquad
\Var_\mu h\le m^{-1}\E_\mu|\nabla h|^2
\]
for $C^1$ $h$ with polynomial growth of value and gradient. The
diffusion $dX_t=-\nabla H(X_t)\,dt+\sqrt2\,dB_t$ has invariant law
$\mu$, and synchronous solutions contract in distance by $e^{-mt}$.
If $\bar H$ is another strongly convex $C^2$ potential with bounded
Hessian and law $\eta$, then
\begin{equation}\label{eq:synchronous-comparison}
\Wtwo(\mu,\eta)
\le \frac1m\left(\E_\eta|\nabla H-\nabla\bar H|^2\right)^{1/2}.
\end{equation}
\end{lemma}
\begin{proof}
Strong convexity yields Gaussian tails and 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. A globally Lipschitz drift makes the diffusion
nonexplosive, and differentiating the squared distance of two
synchronous solutions against strong monotonicity gives the $e^{-mt}$
contraction.
\emph{Poincar\'e.} First take $H$ smooth with bounded derivatives of
order at least two. For smooth compactly supported $h$ let
$P_th(x)=\E[h(X_t)\mid X_0=x]$. It\^o's formula and the semigroup
identity give $\partial_tP_th=P_t\mathcal Ah=\mathcal AP_th$, where
$\mathcal A=\Delta-\nabla H\cdot\nabla$, and the second equality also
follows by differentiating $P_{t+s}h=P_s(P_th)$ at $s=0$. Then
$\int\mathcal AP_th\,d\mu=0$ by parts, so $\mu$ is invariant; coupling
to an invariant start and using contraction gives $P_th\to\mu h$ in
$L^2(\mu)$. Hence
$\Var_\mu h=2\int_0^\infty\!\int|\nabla P_th|^2\,d\mu\,dt$, the flow
derivative has norm at most $e^{-mt}$, and Jensen bounds the inner
integral, proving Poincar\'e.
For $C^2$ potentials, convolve with smooth compactly supported
mollifiers of radii tending to zero: both Hessian bounds survive and all
higher derivatives are bounded for each mollification. Normalized laws
converge in total variation by the quadratic tail bounds (subtract the
value at zero), and gradients converge uniformly by the global
Lipschitz bound. Synchronous solutions of mollified and original
diffusions converge uniformly on bounded time intervals by Gronwall, so
invariance passes to the limit; Poincar\'e first passes on smooth
compact tests, then by cutoff.
For \eqref{eq:synchronous-comparison}, start the two diffusions from
their invariant marginals with common Brownian motion and let
$a(t)=(\E|X_t-Y_t|^2)^{1/2}$, $\beta=(\E_\eta|\nabla H-\nabla\bar H|^2)^{1/2}$.
Strong monotonicity and Cauchy--Schwarz give, in the upper-derivative
sense, $a'(t)\le-ma(t)+\beta$. The marginals are fixed, so
$\Wtwo(\mu,\eta)\le a(t)$; letting $t\to\infty$ proves it.
\end{proof}
These apply to the auxiliary laws without convexity of $F$. The mode
of $\nu_{x,r}$ satisfies $z_0=-r g(x+rz_0)$; for $l\le1/2$,
\begin{equation}\label{eq:primitive-moments}
|z_0|\le\frac{r|g(x)|}{1-l},\qquad
\bigl(\E_\nu|Z|^2\bigr)^{1/2}\le C(\sqrt d+r|g(x)|),\qquad
|u|\le CD(x).
\end{equation}
Indeed $|g(x+rz_0)-g(x)|\le L|z_0|$ gives the first, the lemma gives
$\norm{Z-z_0}_2\le\sqrt{d/(1-l)}$, and $|u|\le|g(x)|+L\norm{Z}_2$
gives the third. Thus $D(x)$ controls all moments uniformly in the
center, even where the gradient is huge.
Errors of a deterministic map on a fixed input law are measured in root
mean square; errors of randomized outputs in $\Wtwo$. The next rules
keep these two notions apart.
\begin{lemma}[Conditional comparisons and smoothing]\label{lem:coupling-smoothing}
Let $(\mu_z),(\eta_z)$ be probability kernels on Euclidean space and $Z$
have law $\zeta$. Whenever second moments are finite,
\[
\Wtwo^2\!\left(\int\mu_z\,d\zeta(z),\int\eta_z\,d\zeta(z)\right)
\le\int\Wtwo^2(\mu_z,\eta_z)\,d\zeta(z).
\]
For finite families, independent sums satisfy
\[
\Wtwo\!\left(\mathcal L\Bigl(\sum_i a_iX_i\Bigr),
\mathcal L\Bigl(\sum_i a_iY_i\Bigr)\right)
\le\sum_i|a_i|\Wtwo(\mathcal L(X_i),\mathcal L(Y_i)),
\]
independence being required within each family. For $s>0$,
\begin{equation}\label{eq:gaussian-smoothing}
\norm{\mu*N(0,s^2I)-\eta*N(0,s^2I)}_{\TV}\le s^{-1}\Wtwo(\mu,\eta).
\end{equation}
A common probability kernel cannot increase total variation, and
$\bigl\|\int\mu_z\,d\zeta-\int\eta_z\,d\zeta\bigr\|_{\TV}
\le\int\norm{\mu_z-\eta_z}_{\TV}\,d\zeta$.
\end{lemma}
\begin{proof}
Use conditional near-optimal couplings for the first inequality and
their independent product for the second, then the $L^2$ triangle
inequality. Measurability of the kernels is obtained by quantizing each
output law on a finite grid inside a box, sending the complement to
zero, choosing a minimizing plan on the finite sets by a fixed ordering
of basic feasible solutions, lifting via the conditional cell measures,
and refining. Finite second moments give the limit.
For centers $x,y$ the relative entropy of $N(x,s^2I)$ from
$N(y,s^2I)$ is $|x-y|^2/(2s^2)$, so Pinsker gives TV at most
$|x-y|/(2s)$; integrate over a coupling and apply Cauchy--Schwarz for
\eqref{eq:gaussian-smoothing}. The total-variation rules follow by
testing functions valued in $[0,1]$.
\end{proof}
Deterministic root-mean-square errors are also combined using
absolute coefficient sums; this uses only the ordinary $L^2$ triangle
inequality and needs no independence.
A standard Gaussian is transported to the auxiliary law $\xi$ by
Gaussian interpolation, whose velocity is a conditional gradient mean.
The point of the Gaussian is that it smooths: although $F$ is only
$C^2$, the conditional mean is smooth in the conditioning variable, and
every derivative tensor can be bounded.
The bounds are stated as matrix norms under every split of the tensor
slots, so later contractions never pay a dimension factor. The
interpolation itself is the stochastic interpolant of
Albergo, Boffi, and Vanden-Eijnden (Section 5.1, equations (5.5)--(5.6)).
\url{https://jmlr.org/papers/volume26/23-1605/23-1605.pdf}
Retain the notation fixed earlier for the near-quadratic law and assume
$l\le1/2$. For independent $Z\sim\nu$ and $N\sim N(0,I_d)$ set
\[
a=1-\rho^2,\qquad Y=\rho Z+\sqrt a\,N,\qquad
\xi_\rho=\mathcal L(Y),\qquad 0\le\rho<1 ,
\]
write $\xi_\rho(dy)=e^{-H_\rho(y)}dy$ (absorbing the constant in
$H_\rho$), and define $M_\rho(y)=\E[f(Z)\mid Y=y]$, using the
conditional density at every $y$. Completing the square,
\begin{equation}\label{eq:conditional-tilt}
\mathcal L(Z\mid Y=y)(dz)\propto
\exp\!\left(-\frac{|z-\rho y|^2}{2a}-F(x+rz)\right)dz .
\end{equation}
Under $z=\rho y+\sqrt a\,w$, $M_\rho(y)$ is exactly the primitive mean
at center $x+r\rho y$ and scale $r\sqrt a$, so $M_0=u$. Thus each
conditional mean is again an instance of \eqref{eq:primitive-law}, at a
prescribed center and smaller positive scale. This is a density
identity; evaluating it with finitely many queries is taken up later.
A tensor $T\in(\R^D)^{\otimes n}$ becomes a matrix by partitioning its
$n$ slots into two nonempty sets $A,B$ and contracting the $A$ slots:
a linear map $(\R^D)^{\otimes|A|}\to(\R^D)^{\otimes|B|}$ between
Frobenius-normed spaces. Its operator norm is the \emph{split norm}
$|T|_{\mathrm{split}}$ for that partition; $|T|_{\mathrm{split}}$ means
the bound holds for every partition. Permuting slots permutes these
norms, and $|T|_F$ is the plain Frobenius norm. The relevance: if
every factor in a contraction is bounded in split norm, the product is
bounded by matrix submultiplicativity with no dimension factor.
Multivariate Appell polynomials, introduced by Appell (1880), are used
(Anshelevich, Section 2.3); their vanishing expected lower derivatives
link Poincar\'e to tensor
covariance (Adamczak and Wolff, Proposition 3.2 and Theorem 3.3;
Kothari and Steinhardt, Section 4.2).
\url{https://arxiv.org/pdf/math/0311043v3}
\url{https://arxiv.org/pdf/1304.1826v1}
\url{https://arxiv.org/pdf/1711.07465v1}
\url{https://www.numdam.org/item/10.24033/asens.186.pdf}
\begin{lemma}[All-split cumulant bound]\label{lem:appell-cumulants}
Let $S$ be a random vector in $\R^D$ with moment-generating function
finite near zero. Suppose that for some $\beta\ge0$,
\[
\Var(P(S))\le\beta\E|\nabla P(S)|^2
\]
for every polynomial $P$ on $\R^D$. For $n\ge2$ let
$\kappa_n(S)=\nabla_\theta^n\log\E e^{\theta\cdot S}|_{\theta=0}$ be
the cumulant tensor. Then for an absolute $C$,
\[
|\kappa_n(S)|_{\mathrm{split}}\le C^n n^{Cn}\beta^{n/2}
\]
for every nontrivial split, independently of $D$ and of the mean.
\end{lemma}
\begin{proof}
Covariances of centered polynomial tensors are first bounded without a
dimension factor, and cumulants are then recovered. If $\beta=0$, testing
coordinate functions shows $S$ is almost surely constant and every
cumulant vanishes; assume $\beta>0$.
Define symmetric Appell tensor polynomials
\[
P_j(s)=\nabla_\theta^j
\frac{\exp(\theta\cdot s)}{\E\exp(\theta\cdot S)}\Big|_{\theta=0}.
\]
For a coefficient tensor $T$ of order $j$ put $Q_T(s)=T:P_j(s)$, colon
contracting all slots. Differentiating the generating function in $s$
multiplies it by copies of $\theta$, so
\[
\E\nabla^h Q_T(S)=0\ (0\le h<j),\qquad
\nabla^j Q_T=j!\operatorname{Sym}(T),
\]
with $|\operatorname{Sym}(T)|_F\le|T|_F$. Apply Poincar\'e to each
component of $\nabla^hQ_T$ and sum, iterating $h=0,\ldots,j-1$:
\begin{equation}\label{eq:appell-l2}
\E|T:P_j(S)|^2\le\beta^j(j!)^2|T|_F^2 .
\end{equation}
Hence for $T,V$ of orders $j,k\ge1$,
\begin{equation}\label{eq:appell-covariance}
\left|\E\bigl[(T:P_j(S))(V:P_k(S))\bigr]\right|
\le\beta^{(j+k)/2}j!k!\,|T|_F|V|_F .
\end{equation}
This controls arbitrary coefficient tensors, not just products of
one-dimensional test directions -- which is precisely what a matrix
split needs.
To recover cumulants while keeping a prescribed split, use the
two-variable generating function
\[
\mathcal G(\theta,\omega)=
\frac{\E\exp((\theta+\omega)\cdot S)}
{\E\exp(\theta\cdot S)\,\E\exp(\omega\cdot S)}.
\]
Its derivative of orders $j,k\ge1$ 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 \eqref{eq:appell-covariance}. Since
$\mathcal G(\theta,0)=\mathcal G(0,\omega)=1$, every nonzero term of
$\mathcal G-1$ contains both groups, and in
$\log\mathcal G=\sum_{m\ge1}\frac{(-1)^{m+1}}m(\mathcal G-1)^m$ only
$m\le n$ contributes at total order $n$. Each differentiated product is
a tensor product of cross-moment operators up to permutations within a
side, so its norm is at most the product of theirs. Distributing $n$
labeled slots among at most $n$ ordered nonempty factors gives at most
$n^{n+1}$ assignments, and the factorials inside any one multiply to at
most $n!$; enlarging the constant bounds the total by $C^nn^{Cn}$.
The mixed derivatives of $\log\mathcal G$ equal those of
$\log\E\exp((\theta+\omega)\cdot S)$, since the other two logarithms
involve a single group. These are the cumulants, and the two groups
represent any prescribed nontrivial split.
\end{proof}
In the conditional density $y$ enters through the linear tilt
$\rho y/a$, and each $y$-derivative of $M_\rho$ supplies a factor
$\rho$; dividing it out keeps the bound meaningful at $\rho=0$.
\begin{proposition}[Gaussian jet bounds]\label{prop:gaussian-jets}
Under the above, in particular $l\le1/2$, define for $q\ge1$
\[
U_q(\rho,y)=\rho^{-q}\nabla^qM_\rho(y)\quad(\rho>0).
\]
These extend smoothly to $\rho=0$. With the vector-output slot counted
among the tensor slots, there is an absolute $C$ such that for every
$0\le\rho<1$, $y\in\R^d$, and every split,
\begin{equation}\label{eq:jets}
|U_q(\rho,y)|_{\mathrm{split}}
\le C^q q^{Cq}L(\sqrt a)^{-(q-1)},
\end{equation}
uniformly in dimension, $F$, center, scale, and order $q$. Moreover
$DM_\rho$ and $U_1$ are symmetric, and
\begin{equation}\label{eq:interpolation-score}
\nabla H_\rho(y)=y+r\rho M_\rho(y),\qquad
\|\nabla^2H_\rho\|_{\op}\le C .
\end{equation}
For $j\ge1$ every split of $\nabla^j\nabla^2H_\rho$ is bounded by
\begin{equation}\label{eq:score-derivatives}
C^{j+1}(j+1)^{C(j+1)}(\sqrt a)^{-j}.
\end{equation}
\end{proposition}
\begin{proof}
The conditional potential in \eqref{eq:conditional-tilt} has Hessian
between $(a^{-1}-l)I$ and $(a^{-1}+l)I$, so by
Lemma~\ref{lem:convex-preliminaries} its Poincar\'e constant is at most
$Ca$. Under this conditional law the map
$z\mapsto S(z)=(z,f(z)/L)$ has derivative norm at most $\sqrt2$ because
$f$ is $L$-Lipschitz; pulling back an ambient test function shows the
possibly singular law of $S$ satisfies Poincar\'e with constant $Ca$.
It has exponential moments of every linear functional ($f$ grows at
most linearly and the conditional law has Gaussian tails), so
Lemma~\ref{lem:appell-cumulants} applies. This repeated use of
Poincar\'e never differentiates $f$ repeatedly: each application pulls
back a fresh ambient polynomial and uses only $Df$.
Repeated differentiation of the expectation in the linear tilt
$\rho y/a$ gives joint cumulants, i.e.
\begin{equation}\label{eq:jet-cumulant}
(U_q)_{i,k_1,\ldots,k_q}
=a^{-q}\operatorname{cum}
(f_i(Z),Z_{k_1},\ldots,Z_{k_q}\mid Y=y).
\end{equation}
Restricting the cumulant of $S$ to one $f/L$ coordinate and $q$ ordinary
coordinates, Lemma~\ref{lem:appell-cumulants} gives a factor
$L a^{(q+1)/2}$, and the prefactor $a^{-q}$ yields \eqref{eq:jets}.
All conditional integrals are smooth in $(\rho,y)$ for $\rho<1$:
differentiating densities inserts polynomials with integrable Gaussian
majorants on compact parameter sets, and \eqref{eq:jet-cumulant} gives
the extension at zero without dividing by $\rho$. For the integrations
by parts below, polynomial growth is also recorded: the conditional mode
satisfies $z_0=\rho y-ar g(x+rz_0)$, so $|z_0|\le C(|y|+r|g(x)|)$, and
comparing the density and its normalizing integral to Gaussian integrals
around the mode bounds moments. The same holds for time derivatives;
dimension-dependent constants are harmless.
Conditional integration by parts gives
$\E[Z\mid Y=y]=\rho y-arM_\rho(y)$; differentiating the marginal density
gives the first identity in \eqref{eq:interpolation-score}. For
$\rho>0$ it implies symmetry of $DM_\rho$, hence of $U_1$, and continuity
handles $\rho=0$. Its derivative with \eqref{eq:jets} gives the bounded
Hessian and \eqref{eq:score-derivatives}. All of this uses only the
$C^2$ regularity of $F$.
\end{proof}
\begin{proposition}[Probability transport]\label{prop:transport}
The ODE
\begin{equation}\label{eq:transport}
\frac{dY_\rho}{d\rho}=-rM_\rho(Y_\rho)
\end{equation}
transports $\xi_0=N(0,I_d)$ to $\xi_\rho$ for every $\rho<1$. Its velocity
is globally $Cl$-Lipschitz in space, and on every compact probability-time
interval it has a smooth flow in both directions.
\end{proposition}
\begin{proof}
Along the random interpolation,
$\frac{d}{d\rho}(\rho Z+\sqrt a\,N)=(Z-\rho Y)/a$, whose conditional
mean is $-rM_\rho(Y)$ by the identity above. Thus for smooth compactly
supported $\phi$,
\[
\frac{d}{d\rho}\int\phi\,d\xi_\rho=-r\int M_\rho\cdot\nabla\phi\,d\xi_\rho .
\]
By \eqref{eq:jets} with $q=1$, $r\|DM_\rho\|_{\op}\le Cl$, so the
velocity has at most linear growth and its smooth ODE has a
global-in-space flow and inverse on every compact interval in $\rho<1$.
To identify the pushforward, fix a smooth compactly supported terminal
test and transport it backward; it has zero material derivative and its
support stays bounded by linear growth, and the display shows its
$\xi_\rho$-expectation is constant.
\end{proof}
The remaining difficulty is not spatial smoothness (supplied by the
Gaussian) but control of derivatives \emph{along} the transport while
paying only one $\sqrt d$ for a vector output. That control is the next task.
Gaussian interpolation gave spatial derivatives of all orders from a
$C^2$ potential. Differentiation now proceeds \emph{along} the
transporting trajectories, since that controls the error when a path
integral is replaced by quadrature. A naive Euclidean norm after each
differentiation loses a power of $d$ each time. The fix: keep at least
two tensor slots free throughout, convert to a vector norm only at the
last step, and bound every matrix split.
Keep $L=\lambda r$, $l=rL\le1/2$, and the interpolation
$\xi_\rho=\mathcal L(\rho Z+\sqrt a\,N)$, $a=1-\rho^2$, $0\le\rho<1$.
Expectations and $L^p$ norms use $\xi_\rho$ at the displayed time.
Along $Y'_\rho=-rM_\rho(Y_\rho)$ the chain rule gives
\[
\frac d{d\rho}Q(\rho,Y_\rho)=(\mathcal DQ)(\rho,Y_\rho),\qquad
\mathcal D=\partial_\rho-rM_\rho\cdot\nabla ,
\]
so path derivatives are iterates of $\mathcal D$. If a tensor field
$V$ has a distinguished coordinate slot, its \emph{adjoint divergence}
on that slot is
\[
(\nabla^*V)_I=\sum_{k=1}^d\{-\partial_kV_{I,k}+(\partial_kH_\rho)V_{I,k}\},
\]
characterized by
$\int(\nabla^*V)_I\phi\,d\xi_\rho
=\int\sum_kV_{I,k}\partial_k\phi\,d\xi_\rho$
for smooth compactly supported $\phi$. The identities below express
material derivatives through this operation and contractions of the
normalized jets $U_q=\rho^{-q}\nabla^qM_\rho$.
\begin{lemma}\label{lem:material-identities}
The mean and its jets satisfy
\begin{align}
\mathcal D M&=\nabla^*U_1, \label{eq:material-mean-identity}\\
(\mathcal D U_q)_{i,I}
&=(\nabla^*U_{q+1})_{i,I}
+2r\rho\!\!\sum_{\varnothing\ne S\subseteq I}\!\!
(U_{|S|})_{k,S}(U_{q-|S|+1})_{i,k,I\setminus S},
\quad q\ge1, \label{eq:material-jet-identity}
\end{align}
$I$ a labeled list of $q$ slots, subsets referring to positions and $k$
summed. For any smooth tensor field $V$,
\begin{equation}\label{eq:material-adjoint-identity}
\mathcal D(\nabla^*V)=\nabla^*(\mathcal DV)+r\nabla^*((DM_\rho)V),
\end{equation}
the matrix acting on the slot removed by the divergence. All hold at
$\rho=0$ by continuous extension.
\end{lemma}
\begin{proof}
Differentiate $\E[f_i(Z)\phi(Y)]=\int(M_\rho)_i\phi\,d\xi_\rho$; the
right side differentiates using the transport velocity, the left has
random velocity $(Z-\rho Y)/a$. Subtracting its conditional mean from
the flux gives coefficient
$\Cov(f_i(Z),Z_k\mid Y)/a=(U_1)_{i,k}$ by \eqref{eq:jet-cumulant}, and
integration by parts proves \eqref{eq:material-mean-identity}.
Coordinate differentiation gives
\[
\mathcal D\partial_IM_i=\partial_I\mathcal DM_i
+r\!\!\sum_{\varnothing\ne S\subseteq I}\!\!
(\partial_SM_k)\partial_k\partial_{I\setminus S}M_i .
\]
Insert \eqref{eq:material-mean-identity} and multiply by $\rho^{-q}$:
passing all derivatives inside the divergence gives
$\nabla^*U_{q+1}$; where they hit the score $y+r\rho M_\rho$, a
derivative of the linear part contributes $qU_q/\rho$, cancelled by
the $\rho^{-q}$ prefactor, and the remaining nonlinear and commutator
terms are both
$r\rho\sum_{\varnothing\ne S\subseteq I}(U_{|S|})_{k,S}
(U_{q-|S|+1})_{i,k,I\setminus S}$. That proves
\eqref{eq:material-jet-identity}; each contraction retains a free slot
in each factor.
Finally differentiate the adjoint identity against a transported
scalar test, componentwise in the other slots: the $\mathcal D\phi$
terms cancel, and the gradient commutator
$\mathcal D\partial_k\phi-\partial_k\mathcal D\phi
=r(\partial_kM_j)\partial_j\phi$ moved across by the adjoint relation
gives \eqref{eq:material-adjoint-identity}. All coefficients are
nonsingular at zero, so extensions follow from
Proposition~\ref{prop:gaussian-jets}.
\end{proof}
Before iterating, the estimate making product terms
harmless is isolated. Contracting one slot of each of two tensors means summing
over one common coordinate, keeping all other slots free; restricting
to one slot is what makes it matrix composition.
\begin{lemma}[A single contraction]\label{lem:split-contraction}
Let $A,B$ be tensors, contract one slot of each leaving at least one
free slot in each. For any nontrivial split of the resulting free
slots there are nontrivial splits of $A$ and of $B$ with
$|A*B|_{\mathrm{split}}\le|A|_{\mathrm{split}}|B|_{\mathrm{split}}$.
\end{lemma}
\begin{proof}
Call the two output sides input and output. 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 sit on the opposite side
of that factor's unfolding; these requirements are compatible, since
incompatibility would put every output free slot on one side. If a
factor already has free slots on both sides it imposes none. With this
orientation the contraction is composition of the two matrix actions,
plus identities on unused coordinates and coordinate permutations, all
of operator norm one; submultiplicativity proves it.
\end{proof}
The following records exactly which tensor formulas are estimated and
which contractions are allowed (contracting two slots of one tensor is
not).
\begin{definition}[Tensor expressions and weights]\label{def:divergence-expressions}
A tensor expression is a finite formula built by:
(1) for each $q\ge1$, $U_q/L$ is an expression with $q+1$ free slots and
weight $q-1$; (2) permuting free slots preserves weight; (3) contracting
one chosen slot of two expressions, each keeping a free slot, adds the
weights; (4) an adjoint divergence on a chosen free slot is allowed
while at least two slots remain and raises the weight by one. The
formula fixes its leaves, orders, and slots, a finite syntactic datum
independent of $d,\rho,x,r$ and $F$. Repeated occurrences are separate
factors even when equal in value.
\end{definition}
For instance $\nabla^*(U_2/L)$ has two slots and weight two, while
$\nabla^*(U_1/L)$ is not an expression (vector output) -- the final
$\sqrt d$. If $E$ is an expression, $\nabla^jE$ has $j$ more slots, a
split placing any on either side. The next theorem bounds all splits at
once, with explicit dependence on $j$ and $p$ because a later proof
takes a moment order growing with $\log d$.
\begin{theorem}[Tensor divergence bound]\label{thm:tensor-divergence}
Let $E$ be a tensor expression of weight $k$. There is $C_E\ge1$,
depending only on the finite formula, such that for every integer
$j\ge0$, every real $2\le p<\infty$, every $0\le\rho<1$, and every
nontrivial split,
\begin{equation}\label{eq:network}
\bigl\|\,|\nabla^jE|_{\mathrm{split}}\,\bigr\|_{L^p(\xi_\rho)}
\le C_E^{j+1}
\bigl(p+j+\log(d+1)\bigr)^{C_E(j+1)}(\sqrt a)^{-(k+j)} .
\end{equation}
The same $C_E$ works for all $j,p,d,\rho,x,r$.
\end{theorem}
\begin{proof}
Induct on the finite construction of $E$, proving all derivative and
moment orders at once. Permutations and single contractions follow from
Lemma~\ref{lem:split-contraction}. A normalized jet has
$\nabla^jU_q=\rho^jU_{q+j}$, so Proposition~\ref{prop:gaussian-jets}
handles leaves after enlarging a constant in $q$. For a contraction,
distribute derivatives by Leibniz (at most $2^j$ summands), apply
Lemma~\ref{lem:split-contraction} to each, and use the inductive
estimates at exponent $2p$, giving \eqref{eq:network} with an enlarged
constant.
It remains to consider $E=\nabla^*V$, $V$ of weight $k_V$, $E$ with at
least two slots. Commuting $j$ derivatives past the divergence leaves a
main term $T=\nabla^*(\nabla^jV)$; the others are single contractions of
$\nabla^{j-h}V$ with an order-$h$ score derivative, $1\le h\le j$, each
bounded pointwise by $C^hh^{Ch}(\sqrt a)^{-(h-1)}$ via
\eqref{eq:interpolation-score}--\eqref{eq:score-derivatives}, so
Lemma~\ref{lem:split-contraction}, the inductive hypothesis, and
Leibniz multiplicities bound them by the required power. Only $T$
remains.
\emph{Removing the adjoints.} Fix an unfolding of $T$ (both sides
nonempty) and $s\ge1$. Pointwise $\|T\|_{\op}^{2s}\le
\tr((TT^T)^s)$; expand the trace by indices, giving $2s$ labeled primary
factors, each an adjoint divergence of $\nabla^jV$, joined in the
alternating trace cycle with nonempty incident bundles. Remove adjoints
one at a time by parts; a chosen primary differentiates the product of
the \emph{others}. If a derivative hits a remaining adjoint, use
$\partial_i(\nabla^*W)=\nabla^*(\partial_iW)
+(\partial_i\partial_kH_\rho)W_k$: the first term leaves that adjoint,
the second consumes it and creates a Hessian factor joining the two
divergence indices. After all removals 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$, with
\begin{equation}\label{eq:network-degree-count}
\sum_{i=1}^{2s}e_i+\sum_{v=1}^{n_H}h_v+2n_H=2s,
\end{equation}
since a removed adjoint either adds one derivative or pairs with another
to create a Hessian. Each step has at most $Cs$ targets and two
commutator terms, so at most $(Cs)^{Cs}$ terms result.
Use a graph only to describe these sums: primary vertices are the $2s$
labeled tensor factors, a Hessian vertex is one factor
$\nabla^{h_v}\nabla^2H_\rho$, an edge is a summed slot pair, a bundle
several such edges. Three properties control contraction: all original
trace bundles remain; every new direct edge joins distinct primaries (an
adjoint never differentiates its own field); every Hessian has at least
two distinct primary neighbors (its creators) and no Hessian-Hessian
edges.
\emph{Opening the trace.} Cut one bundle of the original trace cycle,
leaving $b\ge1$ indices unsummed at its ends; the rest form a chain
through the primaries. Order primaries along it. Each Hessian has two
distinct creating primaries at different positions; place it strictly
between its earliest and latest neighbors (possible because no
Hessian-Hessian edge). Order Hessians in the same gap arbitrarily and
orient every uncut edge from earlier to later. The $b$ cut indices
enter the first primary and leave the last.
Every tensor now has an input and an output. An interior primary
receives one original chain bundle and sends the next, with the cut
bundle supplying the missing side at the endpoints; a Hessian receives
from an earlier and sends to a later neighbor. Thus all matrix
unfoldings are nontrivial splits, as required. Concretely, retain one
copy of $\R^d$ per index whose earlier endpoint is processed but later is
not, start from $(\R^d)^{\otimes b}$ on the cut inputs, permute the
incoming indices of the next factor adjacent, apply that factor's
unfolding and the identity on the rest, then permute the outputs; matrix
multiplication performs the sums. If $Q$ maps
$(\R^d)^{\otimes b}\to(\R^d)^{\otimes b}$ is the induced map,
\[
\|Q\|_{\op}\le\prod_{i=1}^{2s}|\nabla^{j+e_i}V|_{\mathrm{split}}
\prod_{v=1}^{n_H}|\nabla^{h_v}\nabla^2H_\rho|_{\mathrm{split}},
\]
permutations and identities having norm one, and closing the cut is one
matrix trace:
\begin{equation}\label{eq:divergence-single-trace}
|\tr Q|\le d^b\|Q\|_{\op}.
\end{equation}
That is why other cycles add no dimension factor -- their sums are
already done by the compositions.
\emph{Uniform high-moment bookkeeping.} Apply the inductive hypothesis to
the $2s$ primaries at exponent $2s$ and the pointwise Hessian estimates
elsewhere (H\"older uses no larger exponent). The power of
$1/\sqrt a$ before the $2s$-th root is at most
$2s(k_V+j)+\sum_ie_i+\sum_vh_v\le2s(k_V+j+1)$; polynomial bases are at
most $C(s+j+\log(d+1))$, with total exponent at most
\[
C\Bigl(2s(j+1)+\sum_ie_i+\sum_v(h_v+1)\Bigr)\le C's(j+1)
\]
by \eqref{eq:network-degree-count}, absorbing $(Cs)^{Cs}$. The cut
bundle uses only the original free slots of $T$, so $b\le j+C_V$, $C_V$
depending only on $V$'s fixed slot count. Thus
\[
\E\|T\|_{\op}^{2s}\le d^{j+C_V}
[C(s+j+\log(d+1))]^{Cs(j+1)}(\sqrt a)^{-2s(k_V+j+1)}.
\]
Take $s$ smallest with $2s\ge p$ and $s\ge(j+C_V)\log(d+1)$; then
$d^{(j+C_V)/(2s)}\le e^{1/2}$ and
$s\le C_V(p+(j+1)\log(d+1))\le C(p+j+\log(d+1))^2$. Taking the
$2s$-th root and using $L^p$ monotonicity gives \eqref{eq:network} with
$k=k_V+1$. All parts are legitimate for each finite $j,s,d$ and each
$\rho<1$ by the polynomial growth of Proposition~\ref{prop:gaussian-jets}
together with Gaussian tails; cutoffs remove boundary terms. Uniformity
in $s$ comes from the explicit bounds after the identities, not from an
unproved regularity claim.
\end{proof}
The theorem controls every tensor expression while at least two slots
remain; only the last divergence to a vector costs a $\sqrt d$.
\begin{lemma}[Frobenius divergence estimate]\label{lem:divergence-l2}
Let $\mu(dy)\propto e^{-H(y)}dy$, $H\in C^2(\R^d)$ with bounded Hessian
and Gaussian tails, and $V$ smooth with a distinguished divergence slot
and $V,\nabla V\in L^2(\mu)$ in Frobenius norm. Then
\begin{equation}\label{eq:divergence-l2}
\|\nabla^*V\|_{F,2}\le C\bigl(\|V\|_{F,2}+\|\nabla V\|_{F,2}\bigr),
\end{equation}
$C$ depending only on $\|\nabla^2H\|_{\op,\infty}$, independent of $d$
and of the number of remaining slots. For $H=H_\rho$ it is absolute.
\end{lemma}
\begin{proof}
For compactly supported smooth $V$ outputting a scalar, parts twice:
\[
\int(\nabla^*V)^2d\mu
=\int\sum_{i,j}\bigl[(\partial_iV_j)(\partial_jV_i)
+V_i(\partial_{ij}H)V_j\bigr]d\mu ,
\]
the first term bounded by $|\nabla V|_F^2$ and the second by
$\|\nabla^2H\|_{\op}|V|^2$. For tensor output apply to each fixed
collection of remaining indices and sum; a square root gives
\eqref{eq:divergence-l2}. For general $V$ multiply by cutoffs $\chi_R$
equal to one on the ball of radius $R$ with $|\nabla\chi_R|\le C/R$;
then $\chi_RV\to V$, $\nabla(\chi_RV)\to\nabla V$ in $L^2(\mu)$, and
adjoint divergences converge pointwise, so Fatou proves the estimate.
Every field used has polynomially growing derivatives, hence the required
integrability.
\end{proof}
\begin{corollary}[Derivatives along transported paths]\label{cor:path-derivatives}
For each fixed integer $m\ge1$,
\begin{equation}\label{eq:path-mean}
\|\mathcal D^mM_\rho\|_{L^2(\xi_\rho)}
\le\polylog(d)L\sqrt d\,a^{-m}.
\end{equation}
For fixed integers $q\ge1$, $m\ge0$, fixed real $2\le p<\infty$, every
nontrivial split satisfies
\begin{equation}\label{eq:path-jet}
\bigl\|\,|\mathcal D^mU_q|_{\mathrm{split}}\,\bigr\|_{L^p(\xi_\rho)}
\le\polylog(d)L(\sqrt a)^{-(q-1+2m)}.
\end{equation}
The $\polylog$ constants may depend on fixed $m,q,p$ but not on
$\rho$, center, scale, or oracle.
\end{corollary}
\begin{proof}
Iterate Lemma~\ref{lem:material-identities} with the product rule on
contractions: replacing $U_q$ by $\nabla^*U_{q+1}$ raises weight by two;
a contracted product in \eqref{eq:material-jet-identity} does not (jet
weights sum to $q-1$); the commutator
\eqref{eq:material-adjoint-identity} inserts $DM_\rho=\rho U_1$ of weight
zero. Each extra jet factor $L$ carries an $r$, so after one overall $L$
the extras are powers of $l=rL\le1$; scalar coefficients are bounded
polynomials in $\rho$. Starting from $U_q$ all divergences keep at least
two slots and all contractions keep a slot of each factor, so the terms
are expressions of Theorem~\ref{thm:tensor-divergence} of weight at most
$q-1+2m$, and $j=0$ gives \eqref{eq:path-jet}.
Starting from $M$, the first identity is $\mathcal DM=\nabla^*U_1$, and
further material derivatives keep the outermost divergence a divergence
to a vector; inside it the construction gives finitely many normalized
expressions of weight at most $2(m-1)$, each with one overall $L$. Apply
Lemma~\ref{lem:divergence-l2} there. A tensor with at least two slots
has Frobenius norm at most $\sqrt d$ times the operator norm of an
unfolding putting one coordinate slot on one side (that matrix has at
most $d$ singular values), so Theorem~\ref{thm:tensor-divergence} at
$j=0,1$, $p=2$ bounds the right side of \eqref{eq:divergence-l2} by
\[
\polylog(d)L\sqrt d\,(\sqrt a)^{-(2m-1)}
\le\polylog(d)L\sqrt d\,a^{-m}.
\]
\end{proof}
For quadrature, start the transport at its correct marginal and set
$Q_\rho=M_\rho(Y_\rho)$, which has law $\xi_\rho$ by
Proposition~\ref{prop:transport}; \eqref{eq:path-mean} bounds
$\|Q_\rho^{(m)}\|_2$ at every time. On a cell
$[\rho_0,\rho_1]\subset[0,1)$ of length $h_c$, Taylor's integral
remainder with the $L^2$ triangle inequality gives, for
$\rho\in[\rho_0,\rho_1]$,
\[
\Bigl\|Q_\rho-\sum_{h<m}\frac{(\rho-\rho_0)^h}{h!}Q_{\rho_0}^{(h)}\Bigr\|_2
\le\frac{h_c^m}{m!}\sup_{t\in[\rho_0,\rho_1]}\|Q_t^{(m)}\|_2
\le C_m\polylog(d)L\sqrt d\,h_c^m(1-\rho_1^2)^{-m}.
\]
So mesh ratio $h_c\le\theta(1-\rho_1)$ makes the remainder at most
$C_m\polylog(d)L\sqrt d\,\theta^m$: a supremum of deterministic RMS
bounds, no expected supremum over a random path. This cancels cell
length against endpoint singularity in the quadrature below.
The probability transport turns conditional means into samples. The
complementary operation is a path integral whose law is
$N(u,s^2I_d)$ for prescribed $s>0$. Represent the centered observable
$M_*-u$ as a matrix adjoint divergence; that matrix couples a position
with law $\xi_*$ to a Gaussian momentum while preserving their product
law, and integrating the momentum equation yields the noisy mean. This
is still exact mathematics -- both its velocity and its ideal initial
law are implemented later. Keep the earlier interpolation notation:
$L=\lambda r$, $l=\lambda r^2\le1/2$, $M_0=u$,
$DM_\rho=\rho U_1$. Fix $0<R<1/2$, $b_*=(1-R^2)^{1/2}$, and use a star
for evaluation at $\rho=b_*$. The positive radius $R$ keeps all
coefficients smooth although $F$ is only $C^2$. The symbol $\rho$ is the
interpolation parameter; the flow uses a separate time $t$.
Let $T_\rho=T_{\rho\to b_*}$ be the probability transport
\eqref{eq:transport} from $\rho$ to $b_*$. Given terminal $y$ define
\[
Y_\rho(y)=T_\rho^{-1}(y),\qquad B_\rho(y)=DT_\rho(Y_\rho(y))=(DY_\rho(y))^{-1},
\]
so $Y_{b_*}(y)=y$. The backward path satisfies
$Y_\rho'=-rM_\rho(Y_\rho)$ with this terminal condition, and
$J_\rho=DY_\rho$ solves $J_\rho'=-rDM_\rho(Y_\rho)J_\rho$, $J_{b_*}=I_d$.
Inverting gives
\begin{equation}\label{eq:backward-jacobian}
B_\rho'=rB_\rho DM_\rho(Y_\rho)=r\rho B_\rho U_1(\rho,Y_\rho),\qquad B_0'=0.
\end{equation}
The last is an exact endpoint identity, not an approximation. Also
$\|B_\rho\|_{\op}\le C$ since $r\|DM_\rho\|_{\op}\le Cl$ and the time
interval has length at most one.
\begin{lemma}[Transport derivatives]\label{lem:transport-derivatives}
For each fixed $k\ge1$, order-$k$ spatial derivatives of any transport
between times in $[0,b_*]$, after subtracting the identity, have all
nontrivial split norms at most $C_klR^{-(k-1)}$. For each fixed $m\ge1$
and finite $p\ge1$, uniformly in $0\le\rho\le b_*$,
\[
\bigl\|\|B_\rho^{(m)}(y)\|_{\op}\bigr\|_{L^p(y\sim\xi_*)}
\le\polylog(d)\,l(1-\rho^2)^{-m}.
\]
Constants and the $\polylog$ may depend on fixed orders and $p$, not on
center or scale.
\end{lemma}
\begin{proof}
Differentiate the transport equation in the starting position.
Gronwall 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 with an
order-$j$ velocity derivative and state derivatives of orders
$i_1,\ldots,i_j\ge1$, $\sum i_\nu=k$. By \eqref{eq:jets} the outer
tensor has split norms at most $C_jlR^{-(j-1)}$; attach inner tensors one
at a time, each a single contraction retaining the outer output and inner
slots, so Lemma~\ref{lem:split-contraction} applies with no dimension
factor. Bound first state derivatives by a constant and use induction
for higher ones. Since $l\le1$, the total power of $R^{-1}$ is at most
$(j-1)+\sum_\nu(i_\nu-1)=k-1$ with each forcing bounded by $Cl$; the
unknown order-$k$ derivative enters linearly with coefficient $\le Cl$
and starts at zero, so Gronwall proves it. The same holds with absolute
interval lengths for reverse transport.
For the second assertion write $C_\rho=rDM_\rho(Y_\rho)=r\rho U_1(\rho,Y_\rho)$,
so $B_\rho'=B_\rho C_\rho$. The product rule writes $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 $\sum j_\ell=m-k$. The
zero-order jet bound and \eqref{eq:path-jet} give
$\|\|C_\rho^{(j)}\|_{\op}\|_{L^p}\le\polylog(d)l(1-\rho^2)^{-j}$
(inverse transport sends $y\sim\xi_*$ to $Y_\rho(y)\sim\xi_\rho$). H\"older
with high fixed moments, $l^k\le l$, and
$(1-\rho^2)^{-(m-k)}\le(1-\rho^2)^{-m}$ finish; moments below one follow
from $L^1$.
\end{proof}
For a matrix field $K$ define its row-wise adjoint divergence by
\[
(\nabla_{\xi_*}^*K)_i=\sum_j\{-\partial_jK_{ij}+(\partial_jH_*)K_{ij}\},
\]
i.e. $\nabla_{\xi_*}^*K=h$ means
$\E_{\xi_*}[h_i\phi]=\E_{\xi_*}\sum_jK_{ij}\partial_j\phi$ for smooth
compactly supported scalar $\phi$. The goal is this identity with $h=M_*-u$.
This is a Stein-type integration-by-parts identity; classical coordinate
Stein kernels represent $y-\E_\xi y$ instead, and Courtade, Fathi, and
Pananjady build such kernels from a spectral gap (arXiv v2, Definition 2.1
and Theorem 2.4).
\url{https://arxiv.org/pdf/1703.07707v2}
The observable
$M_*-u$ has mean zero by conditional expectation; the transport formula
below supplies its adjoint-divergence representation \emph{with pointwise
derivative bounds} needed for quadrature. Define
\begin{equation}\label{eq:stein-kernel}
K(y)^T=\int_0^{b_*}B_\rho(y)U_1(\rho,Y_\rho(y))\,d\rho
=\frac1r\int_0^{b_*}B_\rho'(y)\,\frac{d\rho}{\rho};
\end{equation}
the first integrand extends smoothly to zero in the second form
($B_\rho'/(r\rho)=B_\rho U_1$), and the transpose matters (the adjoint
divergence acts on rows of $K$, the velocity uses $K^T$).
\begin{lemma}[The centering kernel]\label{lem:centering-kernel}
$K$ is smooth, $\nabla_{\xi_*}^*K=M_*-u$, $\|K(y)\|_{\op}\le CL$, and
every nontrivial split norm of $\nabla^jK$ ($j\ge1$) is bounded pointwise
by $C_jLR^{-j}$.
\end{lemma}
\begin{proof}
Fix output coordinate $i$ and smooth compactly supported $\phi$. The
transported test $\phi\circ T_\rho$ is constant along transport paths.
Differentiate its product with $M_{\rho,i}$ under $\xi_\rho$ and use
$\mathcal DM_\rho=\nabla_{\xi_\rho}^*U_1$
(Lemma~\ref{lem:material-identities}); integration by parts and the chain
rule give
\[
\frac{d}{d\rho}\E_{\xi_\rho}[M_{\rho,i}\phi\circ T_\rho]
=\E_{\xi_\rho}\sum_{j,k}(U_1)_{ik}(DT_\rho)_{jk}(\partial_j\phi)\circ T_\rho .
\]
Push to $y\sim\xi_*$; symmetry of $U_1$ makes the coefficient of
$\partial_j\phi(y)$ equal $(B_\rho(y)U_1(\rho,Y_\rho(y)))_{ji}$. The
terminal expectation is $\E_{\xi_*}[M_{*,i}\phi]$ and the initial one
$u_i\E_{\xi_*}\phi$ since $M_0=u$. Integrating the derivative proves
$\E_{\xi_*}[(M_{*,i}-u_i)\phi]=\E_{\xi_*}\sum_jK_{ij}\partial_j\phi$; the
identity is weak and, by smoothness, pointwise.
For bounds, $\|B_\rho\|_{\op}\le C$, $\|U_1\|_{\op}\le CL$ give
$\|K\|_{\op}\le CL$. For positive-order derivatives use
Lemma~\ref{lem:transport-derivatives} on $DT_\rho$ and $Y_\rho$: an
order-$j$ derivative of $DT_\rho$ has split norms $\le C_jlR^{-j}$, of
$U_1$ at most $C_jLR^{-j}$ by \eqref{eq:jets}. Composing with $Y_\rho$,
attach each chain-rule tensor by a single edge as before; first
derivatives of $Y_\rho$ are bounded and order $i>1$ costs
$C_ilR^{-(i-1)}$, so every order-$j$ term costs at most $C_jLR^{-j}$
since $l\le1$, and each factor keeps free slots so the contraction
preserves every split. Integration over an interval of length at most
one finishes; all orders are smooth by Gaussian conditional integration
on $\rho\le b_*<1$, with no higher derivative of the unsmoothed gradient.
\end{proof}
The next dynamics specialize the skew-matrix construction of
Ma, Chen, and Fox (Section 2.1, equation (3), and Theorem 1).
\url{https://proceedings.neurips.cc/paper/2015/file/9a4400501febb2a95e79248486a5f6d3-Paper.pdf}
The transport-derived
$K$ makes the momentum displacement encode $M_*-u$. Completeness and
invariance are proved directly, since a formal stationary-density
calculation alone would not ensure the solution exists on the whole
interval.
\begin{proposition}[A noisy mean from a path integral]\label{prop:centering-law}
For $s>0$, let $(Y_0,G_0)$ have law $\xi_*\otimes N(0,I_d)$ and solve
\begin{equation}\label{eq:centering-flow}
\dot Y_t=-K(Y_t)^TG_t/s,\qquad \dot G_t=(M_*(Y_t)-u)/s.
\end{equation}
The solution exists for every finite time and its law stays
$\xi_*\otimes N(0,I_d)$. In particular,
\begin{equation}\label{eq:centering-output}
sG_0+\int_0^1M_*(Y_t)\,dt\ \stackrel{\mathrm{law}}=
u+sG,\qquad G\sim N(0,I_d).
\end{equation}
\end{proposition}
\begin{proof}
The vector field $v$ in \eqref{eq:centering-flow} is smooth and locally
Lipschitz; the bound on $K$ and global Lipschitz bound on $M_*$ give
$|v(y,g)|\le C_{x,r,s}(1+|y|+|g|)$, so Gronwall prevents escape in finite
time in either direction. Local uniqueness gives a complete smooth flow
$\Phi_t$ with inverse $\Phi_{-t}$; only local Lipschitz is needed (the
factor $G$ does not give a globally Lipschitz joint field).
Write the product density $p(y,g)=(2\pi)^{-d/2}\exp(-H_*(y)-|g|^2/2)$.
Differentiating in position and momentum separately (row convention),
\[
\operatorname{div}_y(pv_Y)=\frac p s\,g\cdot(M_*-u),\qquad
\operatorname{div}_g(pv_G)=-\frac p s\,g\cdot(M_*-u),
\]
so $\operatorname{div}(pv)=0$, and the Jacobian formula for the complete
flow gives
$\frac{d}{dt}\{p(\Phi_t(z))\det D\Phi_t(z)\}=0$ pointwise. The
determinant starts at one and stays positive; change of variables proves
$\Phi_t$ preserves $p$.
The same cancellation is the operator identity
\begin{equation}\label{eq:centering-skew-operator}
\mathcal L=\nabla^*A\nabla,\qquad
A=\begin{pmatrix}0&K^T/s\\-K/s&0\end{pmatrix},
\end{equation}
the second-order part vanishing by skew symmetry and the first-order
coefficients being $-K^Tg/s$ and $\nabla_{\xi_*}^*K/s=(M_*-u)/s$.
Integrating the momentum equation over $[0,1]$ gives the pathwise
identity $sG_0+\int_0^1M_*(Y_t)\,dt=u+sG_1$, and stationarity makes $G_1$
standard Gaussian.
\end{proof}
Stationarity also makes $G_1$ independent of $Y_1$, but not of the
initial anchors; the exact mean $u=M_0$ stays in the velocity and the
proposition assumes an ideal sample of $\xi_*$. So it is a comparison
construction, not yet a routine. Finitely many exact means approximate
it, and every such mean, including $M_0$, is subsequently replaced.
High-order quadrature needs time-derivative bounds at each fixed time;
stationarity allows them to be proved in the fixed product measure, without a
bound on the supremum of a random path.
\begin{lemma}[Stationary time derivatives]\label{lem:stationary-derivatives}
Fix $0<R\le1$. Let a smooth probability density on $\R^n$, $n\le Cd$,
have Gaussian tails and negative log $H$ with bounded Hessian and every
one-slot split norm of $\nabla^j\nabla H$ at most $C_jR^{-(j-1)}$
($j\ge1$; a one-slot split isolates one slot as matrix input). Suppose a
stationary autonomous flow has Lie derivative
$\mathcal L=\nabla^*A\nabla$ for this density, $A$ skew,
$|\nabla^jA|_{\text{one-slot}}\le C_j\alpha R^{-j}$ ($j\ge0$). Let
$\phi_0$ be smooth, scalar or vector, with
$|\nabla^j\phi_0|_F\le C_ja_0\sqrt d\,R^{-j}$ ($j\ge1$), and assume
polynomial growth. Then for each fixed $k\ge1$ the order-$k$ time
derivative of $\phi_0$ on the flow has RMS norm at most
\begin{equation}\label{eq:stationary-derivatives}
C_ka_0\sqrt d\,(\alpha R^{-2})^k.
\end{equation}
\end{lemma}
\begin{proof}
Put $\phi_k=\mathcal L^k\phi_0$. For $j\ge0$, $k\ge1$ the
following stronger estimate is proved
\begin{equation}\label{eq:stationary-spatial-induction}
\|\nabla^j\phi_k\|_{F,2}\le C_{j,k}a_0\sqrt d\,\alpha^kR^{-(j+2k)}.
\end{equation}
For $k=0$ it is assumed only for $j\ge1$, which suffices (the iteration
differentiates $\phi_{k-1}$ before multiplying by $A$). All $L^2$ norms
are in the invariant measure.
Set $V=A\nabla\phi_{k-1}$. Leibniz and the hypothesis give, for each
$h\ge0$,
$\|\nabla^hV\|_{F,2}\le C_{h,k}a_0\sqrt d\,\alpha^kR^{-(h+1+2(k-1))}$.
No dimension factor appears because in each coefficient tensor
$\nabla^vA$ one isolates the slot contracted with $\nabla\phi_{k-1}$; its
one-slot operator norm multiplies the Frobenius norm of the rest, and
finitely many Leibniz terms preserve the bound.
Commute $j$ derivatives past the divergence in $\phi_k=\nabla^*V$, the
main term $\nabla^*(\nabla^jV)$ bounded in $L^2$ by a constant times
$\|\nabla^jV\|_{F,2}+\|\nabla^{j+1}V\|_{F,2}$ via
\eqref{eq:divergence-l2} (powers $R^{-(j+2k-1)}$, $R^{-(j+2k)}$). Each
remaining commutator contracts an order-$v$ score derivative,
$1\le v\le j$, with $\nabla^{j-v}V$; again the coefficient's one-slot
norm suffices, and the total power is
$(v-1)+(j-v+1+2(k-1))=j+2k-2$, bounded by
\eqref{eq:stationary-spatial-induction} since $0<R\le1$. The induction
closes, including $j=0$ for $k\ge1$; Gaussian tails and polynomial growth
justify the integrations by parts through cutoffs. On an autonomous
trajectory the $k$th time derivative is $\phi_k$ at the current state, so
stationarity identifies its RMS norm at every time with the fixed-measure
norm at $j=0$.
\end{proof}
For \eqref{eq:centering-flow} the ambient dimension is $2d$ and
\eqref{eq:centering-skew-operator} has coefficient bounds with
$\alpha=L/s$. The product potential $H_*(y)+|g|^2/2$ has score-derivative
bounds from \eqref{eq:jets} and the Gaussian factor. For the joint state
take $a_0=1$ (first derivative Frobenius norm $\sqrt{2d}$, higher
derivatives vanish); for $M_*$ take $a_0=L$ (isolating any slot in
\eqref{eq:jets} gives a matrix of input dimension $d$, so operator to
Frobenius costs at most $\sqrt d$). Hence at every stationary time
\[
\Bigl\|\frac{d^k}{dt^k}(Y_t,G_t)\Bigr\|_2\le C_k\sqrt d\,((L/s)R^{-2})^k,\qquad
\Bigl\|\frac{d^k}{dt^k}M_*(Y_t)\Bigr\|_2\le C_kL\sqrt d\,((L/s)R^{-2})^k,
\]
the strict radius $0<R<1/2$ satisfying the lemma's $0<R\le1$.
An initial displacement estimate is also needed. Conditional expectation
and the primitive Poincar\'e inequality give
$\E_{\xi_*}|M_*-u|^2\le\E_\nu|f-u|^2\le CL^2d$; at every stationary time
$G_t$ is independent of $Y_t$, so
$\E|K(Y_t)^TG_t|^2=\E\|K(Y_t)\|_F^2\le CL^2d$. Integrating the velocity
and using Minkowski,
\begin{equation}\label{eq:centering-initial-displacement}
\|(Y_t,G_t)-(Y_0,G_0)\|_2\le Ct(L/s)\sqrt d\le C(L/s)\sqrt d,
\qquad 0\le t\le1.
\end{equation}
These fixed-time RMS estimates feed the discretization below.
The probability transport and the stationary centering flow are now turned
into finite formulas in conditional means, meeting two requirements at
once: small error on the prescribed input law, and each new mean
evaluated at a center that is an affine expression in earlier means with
scalar coefficients and controlled size. The latter is what permits the
recursive implementation described later.
The main conclusion is Theorem~\ref{thm:numerical-interface}. Three
numerical steps: global Picard iteration discretizes probability
transport without dependency depth proportional to the number of time
cells; short harmonic trajectories give a formula for $K(y)^TG$; Picard
iteration discretizes the stationary flow. Two cancellations matter: a
zero derivative removes the singularity at the lower endpoint of the
kernel integral, and keeping only a forward transport \emph{correction}
cancels an apparent inverse power of $r$. Polynomial collocation of
Picard iteration is due to Lee, Song, and Vempala (Algorithm 2 and
Theorem 2.3); Gaussian smoothing and high-order Picard to Chen, Chewi,
Lu, and Zhang (Appendix E, Algorithm E.1, and Theorem E.2).
\url{https://arxiv.org/pdf/1812.06243v1}
\url{https://arxiv.org/pdf/2609.06906v1}
Scalar
affine formulas are needed whose dependency depth can be fixed before their
interpolation order. Exact conditional means are analytical nodes
throughout; finite queries replace them later.
Fix an integer $K_0\ge2$ and $t>0$ with
\begin{equation}\label{eq:numerical-parameters}
40(K_0+2)t<\tfrac14,\qquad R=\delta^t,\qquad \psi=\delta^{4t},\qquad
h_0=\delta^w,\quad 0<w<t.
\end{equation}
The iteration count first fixes a depth bound; then $w$ is chosen; then
the interpolation orders. For the centering flow set $s=\sigma/2$,
$\sigma>0$ the parent mean-noise parameter. Throughout its analysis
assume
\begin{equation}\label{eq:numerical-size}
\max\{l,L/\sigma\}\le\polylog(d)\delta^b,\qquad b\ge\tfrac12 ;
\end{equation}
for probability transport alone only $l\le\polylog(d)\delta^b$ is
needed. All probability-time intervals lie in
$[0,(1-cR^2)^{1/2}]$ for fixed $c>0$. Constants may depend on fixed
orders but not on center, dimension, scale, noise, or oracle. A symbol
$B_d$ denotes $C(1+\log d)^C$ with this uniformity. All assertions hold
for sufficiently large $d$.
Partition a probability-time interval into steps of size at most $h_0$ in
$-\log(1-\rho)$; for a nonzero interval take the ceiling of the
transformed length over $h_0$ as the number of cells. Then
$\ell_c\le Ch_0(1-\rho)$ on each cell, and the count is
$O(1+h_0^{-1}(1+|\log R|))$; a zero-length interval does nothing.
On each cell interpolate at $m\ge2$ equally spaced nodes including the
ends and integrate the degree-$(m-1)$ polynomial, with partial integrals
to every state node. The scaled fixed Lagrange basis on $[0,1]$ has
absolute value sum at most $C_m$ and order-$k$ derivative absolute sum at
most $C_m\ell_c^{-k}$, so every partial-integral row has absolute weight
sum at most $C_m$ uniformly in the number of cells.
Let $w_{ij}$ be the partial-integral weights from starting time to node
$\rho_i$ (signed backward), duplicating shared endpoints as labeled
nodes if convenient. From starting position $y_{\mathrm{in}}$ use
\begin{equation}\label{eq:probability-update}
y_i^{[0]}=y_{\mathrm{in}},\qquad
y_i^{[k+1]}=y_{\mathrm{in}}-r\sum_jw_{ij}M_{\rho_j}(y_j^{[k]}),
\end{equation}
a global iteration that integrates the whole interval at every layer.
\begin{proposition}[Probability-flow circuit]\label{prop:probability-circuit}
If the input position has marginal law $\xi_{\rho_{\mathrm{in}}}$, then
after $N_p$ iterations the maximum over all state nodes of the RMS error
relative to the exact transport is at most
\begin{equation}\label{eq:probability-picard}
C_m\polylog(d)\,l\sqrt d\,h_0^m+(C_ml)^{N_p}CrD .
\end{equation}
Each computed state is globally Lipschitz in its initial position with
constant $C_m$, and after subtracting the initial position with constant
at most $C_ml$.
\end{proposition}
\begin{proof}
Write $Y_\rho$ for the exact path and $Q_\rho=M_\rho(Y_\rho)$; the
marginal identity in Proposition~\ref{prop:transport} and
\eqref{eq:path-mean} give
$\|Q_\rho^{(m)}\|_2\le B_dL\sqrt d(1-\rho^2)^{-m}$. On a cell
$I=[a_c,b_c]$ of length $\ell_c$ let $P$ be the Taylor polynomial of
degree $m-1$ at $a_c$; Minkowski with the integral remainder gives, at
$v\in I$,
\begin{equation}\label{eq:numerical-taylor-remainder}
\|Q_v-P(v)\|_2\le\frac1{(m-1)!}\int_{a_c}^v(v-z)^{m-1}\|Q_z^{(m)}\|_2dz
\le C_mB_dL\sqrt d\,h_0^m,
\end{equation}
using $\ell_c\le Ch_0(1-z)$ and $1-z^2\ge1-z$. The same holds at nodes,
and since interpolation reproduces $P$ with bounded basis sum,
$\|Q_v-\mathcal I_cQ(v)\|_2$ obeys the same bound. Integrating over the
cells (total length at most one) gives the quadrature error for partial
integrals and either orientation, with no supremum of the random path.
Let $e_k=\max_i\|y_i^{[k]}-Y_{\rho_i}\|_2$. The bound
$\Lip(M_\rho)\le CL$ and absolute row sums give
$e_{k+1}\le C_ml\,e_k+C_mB_dl\sqrt d\,h_0^m$. On each exact transported
marginal, conditional expectation and the primitive moment bound give
$\|M_\rho\|_2\le CD$, and integrating the velocity shows $e_0\le CrD$.
Since $C_ml<1/2$ for large $d$, summing the geometric recurrence proves
\eqref{eq:probability-picard}. The norm is the maximum of the nodewise
RMS errors, not the RMS of their maximum.
For deterministic inputs the same recurrence 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 $CL$-Lipschitz mean give
$C_ml$.
\end{proof}
The position velocity in \eqref{eq:centering-flow} contains $K(y)^TG$.
Approximating $K$ then multiplying the error by an arbitrary momentum
would not give the global Lipschitz property Picard needs. Instead the action is
approximated directly, using short harmonic trajectories
and forward transport corrections.
For the kernel integral use the mesh on $[0,b_*]$ above, whose cell
lengths also satisfy $\ell_c\ge c'h_0(1-\rho)$ for large $d$ (from the
ceiling construction and the transformed length tending to infinity).
Set $C_\rho(y)=B_\rho(y)-I$. A circuit is described for $C_\rho(y)G$ at
one probability node $\rho$.
\paragraph{Backward position.} Apply Proposition~\ref{prop:probability-circuit}
backward from input $y$ at time $b_*$ to approximate $Y_\rho(y)$ by
$\widehat Y_\rho$.
\paragraph{Short harmonic trajectories.} At fixed time $\rho$, Hamiltonian
equations with potential $H_\rho$ have position integral equation
\begin{equation}\label{eq:harmonic-integral}
W_z=\cos(z)Y_\rho+\sin(z)G-r\rho\int_0^z\sin(z-v)M_\rho(W_v)\,dv .
\end{equation}
For each requested angle $z$ interpolate on the single cell $[0,z]$.
With nodes $z_i$, Lagrange basis $\ell_j$, and
$w_{ij}^{\sin}=\int_0^{z_i}\sin(z_i-v)\ell_j(v)\,dv$, use
\begin{equation}\label{eq:harmonic-picard}
W_i^{[0]}=\cos(z_i)\widehat Y_\rho+\sin(z_i)G,\quad
W_i^{[k+1]}=\cos(z_i)\widehat Y_\rho+\sin(z_i)G
-r\rho\sum_jw_{ij}^{\sin}M_\rho(W_j^{[k]}).
\end{equation}
At angle zero there is no integral; the free sine and cosine terms are
kept exactly.
\paragraph{Differentiating the transport correction.} Let $d_j$ be the
weights differentiating at zero the polynomial interpolant on angle nodes
$0,\psi,\ldots,(m-1)\psi$, with absolute sum $C_m\psi^{-1}$. At each
angle node compute $T_\rho(W_z)-W_z$ by the forward probability-flow
circuit, keeping only its integral correction at the output. For exact
paths,
\begin{equation}\label{eq:kernel-action-derivative}
\Bigl.\frac d{dz}(T_\rho-\mathrm{id})(W_z)\Bigr|_{z=0}=(B_\rho-I)G=C_\rho G,
\end{equation}
since $W_0=Y_\rho$, $W'_0=G$; the weighted sum with $d_j$ approximates
$C_\rho G$ using exact means only. Use a common interpolation order and
the same number $N_p\ge1$ of Picard layers in the backward, harmonic, and
forward computations.
\begin{lemma}[Errors in the angle computation]\label{lem:kernel-angle-errors}
Let $y\sim\xi_*$ and $G\sim N(0,I_d)$ independently. The exact-input
angle differentiation in \eqref{eq:kernel-action-derivative} has RMS
truncation error at most
\begin{equation}\label{eq:angle-stencil-error}
C_ml\sqrt d\,\psi^{m-1}R^{-2m}=C_ml\sqrt d\,\delta^{2tm-4t}.
\end{equation}
Before applying the angle weights, the computed forward-correction values
have RMS error at most
\begin{equation}\label{eq:kernel-value-error}
C_m\polylog(d)\,l\sqrt d\,(h_0^m+\delta^{2tm})+(C_ml)^{N_p}CrD .
\end{equation}
Each approximate correction is globally Lipschitz in $(y,G)$ with
constant at most $C_ml$.
\end{lemma}
\begin{proof}
For exact inputs $(Y_\rho(y),G)$, harmonic motion preserves
$\xi_\rho\otimes N(0,I_d)$; the field is globally Lipschitz with vanishing
weighted divergence, and its Lie derivative is $\nabla^*A\nabla$ with
constant skew $(W,P)$ block $-I$. So
Lemma~\ref{lem:stationary-derivatives} applies with $\alpha=1$, $a_0=l$
for $T_\rho-\mathrm{id}$ and $a_0=L$ for $M_\rho$, using
Lemma~\ref{lem:transport-derivatives} and \eqref{eq:jets} (operator to
Frobenius costs only $\sqrt d$). Subtracting 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}$; the absolute stencil sum
$C_m/\psi$ gives \eqref{eq:angle-stencil-error}.
Interpolate $M_\rho(W_v)$ in \eqref{eq:harmonic-integral}; since
$z\le(m-1)\psi$, the Taylor remainder times $r$ is at most
$C_ml\sqrt d\,\psi^mR^{-2m}=C_ml\sqrt d\,\delta^{2tm}$, integration
factors bounded. The free-path error is at most $CrD$, Picard contracts
by $C_ml$, and the harmonic state is globally Lipschitz in its inputs by
the row-sum argument.
Compare the forward probability correction started at the exact harmonic
state: that state has marginal law $\xi_\rho$, so
Proposition~\ref{prop:probability-circuit} applies, and errors in it
propagate with $C_ml$, not order one, also controlling the backward-position
error through the harmonic computation. Combining gives
\eqref{eq:kernel-value-error} and the Lipschitz bound.
\end{proof}
It remains to integrate the derivative of $C_\rho G$ against $1/\rho$.
The exact identity $B'_0=0$ removes the singularity, and it is imposed on
the interpolant too. On all but the first probability cell use ordinary
degree-$(m-1)$ interpolation; on the first cell use $m-1$ equally spaced
value nodes including both endpoints plus the Hermite condition that the
derivative at zero is zero, for $m\ge3$ conditions. This determines a
unique polynomial of degree at most $m-1$: a null-space polynomial would
have a double zero at zero and the remaining $m-2$ distinct zeros. Let
$q_i$ be the resulting value weights when the derivative of this piecewise
interpolant is integrated against $1/\rho$.
\begin{lemma}[Derivative quadrature at zero]\label{lem:kernel-endpoint-quadrature}
Let $y\sim\xi_*$, $G\sim N(0,I_d)$ independently. For exact values of
$C_\rho(y)G$, the RMS error of this quadrature for
$\int_0^{b_*}B'_\rho(y)G\,d\rho/\rho$ is at most
$C_m\polylog(d)\,l\sqrt d\,h_0^{m-2}$, and
\begin{equation}\label{eq:kernel-weight-sum}
\sum_i|q_i|\le C_mh_0^{-1}(1+|\log R|+|\log h_0|).
\end{equation}
These weights remain well defined for approximate values: the derivative
condition at zero stays exact.
\end{lemma}
\begin{proof}
By Lemma~\ref{lem:transport-derivatives} and independence of $G$, the
order-$m$ derivative of $C_\rho(y)G$ has RMS norm at most
$\polylog(d)l\sqrt d(1-\rho^2)^{-m}$ (the Gaussian converts matrix to
Frobenius norm, cost $\sqrt d$). On an ordinary cell, Taylor and basis
scaling bound the derivative error by this times $C_m\ell_c^{m-1}$, i.e.
at most $C_m\polylog(d)l\sqrt dh_0^{m-1}/(1-\rho)$; after
$1/\rho$-integration all ordinary cells contribute this prefactor times
$C(1+|\log R|+|\log h_0|)$, within $h_0^{m-2}$.
For the first cell subtract the degree-$(m-1)$ Taylor polynomial at zero,
whose derivative at zero matches the exact zero derivative. The second
derivative of the remainder and of its Hermite interpolant are bounded in
RMS by $C_m\polylog(d)l\sqrt d\,\ell_c^{m-2}$, and both first derivatives
vanish at zero. If $E$ is the difference from the Hermite interpolant,
$\|E'(\rho)\|_2\le\int_0^\rho\|E''(v)\|_2dv\le C_mB_dl\sqrt d\,\rho\ell_c^{m-2}$,
so the $1/\rho$-integral over the first cell is at most
$C_mB_dl\sqrt d\,\ell_c^{m-1}$, within $h_0^{m-2}$.
On an ordinary cell the interpolant-derivative coefficient sum is at most
$C_m/\ell_c$; integrating with $\ell_c\ge c'h_0(1-\rho)$ gives at most
$C_mh_0^{-1}\int_{\rho_{\mathrm{first}}}^{b_*}\frac{d\rho}{\rho(1-\rho)}$,
the logarithms of \eqref{eq:kernel-weight-sum}. On the first cell, scale
to $[0,1]$: the polynomial derivative vanishes at zero and is bounded by
distance-to-zero times the second derivative, so the $1/\rho$-integral of
its absolute value is bounded; scaling back gives $C_m/\ell_c\le
C_mh_0^{-1}$.
\end{proof}
Let $\mathrm{corr}_{ij}(y,G)$ be the forward correction at probability
node $i$ and angle node $j$. The approximate position velocity is
\begin{equation}\label{eq:kernel-approximate-velocity}
\widehat v_Y(y,G)=-\frac1{rs}\sum_iq_i\sum_jd_j\,\mathrm{corr}_{ij}(y,G).
\end{equation}
Each correction keeps only the integral part of the forward Picard formula;
write its weights $w_{ijk}$ and means $M_{ijk}$, so
$\mathrm{corr}_{ij}=-r\sum_kw_{ijk}M_{ijk}$ and
\begin{equation}\label{eq:kernel-scalar-expansion}
\widehat v_Y=\frac1s\sum_{i,j,k}q_id_jw_{ijk}M_{ijk}.
\end{equation}
The apparent $1/r$ cancels algebraically; there is no additive forward
state left, and all coefficients multiplying means are deterministic
scalars. Use the exact expression
$\widehat v_G(y,G)=(M_*(y)-M_0)/s$ for the momentum velocity; the momentum
occurs inside angle states and is never multiplied by a separately
approximated matrix at the output.
\begin{proposition}[Kernel-action circuit]\label{prop:kernel-circuit}
On the stationary product input law, the RMS error of
$\widehat v=(\widehat v_Y,\widehat v_G)$ relative to
\eqref{eq:centering-flow} is bounded by $\polylog(d)$ times
\begin{equation}\label{eq:kernel-velocity-error}
\frac L\sigma\sqrt d\bigl[h_0^{m-2}+h_0^{-1}\delta^{-4t}(h_0^m+\delta^{2tm})\bigr]
+\delta^{-4t-w}(C_ml)^{N_p-1}r\frac L\sigma D .
\end{equation}
It is globally Lipschitz with constant at most
\begin{equation}\label{eq:velocity-lipschitz}
C_m\polylog(d)\,(L/\sigma)\psi^{-1}h_0^{-1}.
\end{equation}
For every fixed $J\ge1$ its RMS error can be made at most
$\polylog(d)\delta^JD$ with $N_p=O(J+1)$, independently of $w$ and the
interpolation orders.
\end{proposition}
\begin{proof}
Combine Lemmas~\ref{lem:kernel-angle-errors} and~\ref{lem:kernel-endpoint-quadrature}
and divide by $rs$. Value errors are amplified only by $C_m\psi^{-1}$ and
the weight sum \eqref{eq:kernel-weight-sum}; the exact stencil truncation
has factor $\psi^{m-1}R^{-2m}=\delta^{2tm-4t}$ and the same probability
weight sum. Crucially $l/(r\sigma)=L/\sigma$ and
$l^{N_p}rD/(r\sigma)=l^{N_p-1}r(L/\sigma)D$, so neither division
introduces a negative power of $r$; the estimate is
\eqref{eq:kernel-velocity-error}.
Each forward correction is $C_ml$-Lipschitz in its inputs and the backward
and harmonic maps bounded; the same absolute weight sums with
$l/(r\sigma)=L/\sigma$ give \eqref{eq:velocity-lipschitz} (the momentum
component smaller, $CL/\sigma$). For the last assertion, increase the
interpolation order until the first line of \eqref{eq:kernel-velocity-error}
has the desired order; under \eqref{eq:numerical-size} $O(J+1)$ Picard
layers make the second line of that order too (the common choice below
verifies independence from $w$ and orders).
\end{proof}
Discretize the centering flow on the single cell $[0,1]$ with ordinary
degree-$(m-1)$ interpolation and partial-integral rows $w_{ij}$. From
anchor $Z_0=(Y_0,G_0)$ set
\[
Z_i^{[0]}=Z_0,\qquad Z_i^{[k+1]}=Z_0+\sum_jw_{ij}\widehat v(Z_j^{[k]}),
\]
and at the end keep the position components $y_j$ and output
\begin{equation}\label{eq:mean-output}
sG_0+\sum_j\omega_jM_*(y_j),\qquad \sum_j|\omega_j|\le C_m,
\end{equation}
$\omega_j$ the full-interval integration weights.
\begin{proposition}[Mean-flow circuit]\label{prop:mean-circuit}
For ideal independent anchors $Y_0\sim\xi_*$, $G_0\sim N(0,I_d)$, and
every fixed $J\ge1$, the circuit \eqref{eq:mean-output} can be chosen to
differ in RMS by at most $\polylog(d)\sigma\delta^JD$ from the exact path
integral \eqref{eq:centering-output}; the approximate joint states have
RMS error at most $\polylog(d)\delta^JD$. Every integration uses $O(J+1)$
Picard layers, a bound independent of $w$ and all interpolation orders.
\end{proposition}
\begin{proof}
Compare the exact flow to polynomial integration of its exact velocity:
Lemma~\ref{lem:stationary-derivatives} (with one more derivative) bounds
this quadrature error by $C_m\sqrt d\,((L/s)R^{-2})^{m+1}$. The
approximation error of $\widehat v$ at each exact trajectory node is
bounded by Proposition~\ref{prop:kernel-circuit} by stationarity; previous
state errors propagate with its global Lipschitz constant and absolute row
sums. (The true centering velocity need not be globally Lipschitz; the
approximate one is.) By \eqref{eq:numerical-parameters}, $b\ge1/2$ and
$w<t$ give $b-4t-w>1/3$, so for large $d$ the constants and polylogarithms
of \eqref{eq:velocity-lipschitz} let the recurrence contract by at most
$\delta^{1/3}$; its initial error is at most $C(L/s)\sqrt d$ by
\eqref{eq:centering-initial-displacement}. Choose the quadrature and
velocity error orders high, then $O(J+1)$ Picard layers give the stated
state error.
For the output apply Lemma~\ref{lem:stationary-derivatives} to $M_*$ with
$a_0=L$: the path quadrature error is at most
$C_mL\sqrt d\,((L/s)R^{-2})^m\le\polylog(d)\sigma\delta^JD$. The extra
error from evaluating at approximate positions is at most their RMS error
times $C_mL$ (using $M_*$ is $CL$-Lipschitz), within the same budget since
$L/\sigma$ is small.
\end{proof}
\begin{corollary}[Accuracy without increasing dependency depth]\label{prop:exact-mean-circuits}
Fix $J\ge1$. In all the probability, harmonic, and centering integrations
one may take at most $N=\lceil10(J+2)\rceil$ Picard layers. After any
fixed $w\in(0,t)$, a common fixed interpolation order satisfying
\begin{equation}\label{eq:interpolation-order-choice}
w(m-3)>J+10,\quad t(m-3)>J+10,\quad m>4(J+10)
\end{equation}
suffices, for all large $d$: probability transport from its exact starting
marginal has endpoint RMS error at most $\polylog(d)\delta^JD$ (also after
division by an endpoint parameter bounded away from zero), and the mean
circuit with ideal anchors has RMS error at most
$\polylog(d)\sigma\delta^JD$ relative to \eqref{eq:centering-output}. The
number of layers of exact-mean dependencies within either circuit is
$O((N+1)^2)$, independent of $m$, $w$, and the number of mesh cells. All
nodes, weights, scalar parameters, and counts are deterministic functions
of the parameters and dimension.
\end{corollary}
\begin{proof}
For large $d$, probability and harmonic Picard contractions
$C_ml\le\polylog(d)\delta^b$ are at most $\delta^{1/3}$. The only extra
inverse powers in kernel evaluation are $\delta^{-4t-w}$; the centering
iteration also contracts by $\delta^{1/3}$. The probability iteration
remainder is at most $C\delta^{N/3}D$, and after kernel weights at most
$C\delta^{N/3-4t-w}D$, both exponents exceeding $J$ with a fixed margin.
For truncation it suffices that the exponents
\[
wm,\quad w(m-2),\quad w(m-1)-4t,\quad 2tm-4t-w,\quad (m+1)/3
\]
exceed $J$: they come respectively from probability quadrature, the
zero-endpoint rule, its amplified probability-value error, the harmonic and
angle errors, and the outer quadratures, where $(L/s)R^{-2}\le\delta^{1/3}$
for large $d$ and $L/\sigma\le\delta^{1/3}$. This is exactly
\eqref{eq:interpolation-order-choice}. Fixed constants and
$m$-dependent polylogarithms fit within the strict margins as $d\to\infty$.
Within a probability or harmonic Picard layer all nodes depend only on the
preceding layer, regardless of node count. A path through one kernel
evaluation passes a backward probability, a harmonic, and a forward
probability computation, each with at most $N$ layers, traversing at most
$N$ outer centering layers: depth $O((N+1)^2)$, final quadratures adding a
bounded number. All weights integrate known polynomials (or polynomials
against sine or $1/\rho$ with the zero-endpoint condition), require no
oracle information, and are allowable real constants.
\end{proof}
The order of choices is essential: first the target error and Picard
depth, then the mesh exponent $w$, finally the interpolation order. The
count of evaluations may depend on that order; the number of nested
dependencies does not.
To display a primitive's parameters write, for $x'\in\R^d$, $0<r'\le1$,
$u_{x',r'}=\int g(x'+r'z)\,\nu_{x',r'}(dz)$, so $u_{x,r}=u$. Next,
$\mathcal C_M,\mathcal C_S$ are deterministic maps of their
anchors with exact primitive means as analytical nodes; the parameter
$\sigma$ of $\mathcal C_M$ is the parent noise, so it sets $s=\sigma/2$
once. Define
\begin{align}
\mathcal C_M(x,r,\sigma;y,G)&=\frac\sigma2G+\sum_j\omega_jM_*(y_j),
\label{eq:exact-mean-circuit}\\
\mathcal C_S(x,r,\eta;y)&=\rho_e^{-1}\widehat T_{0\to\rho_e}(y),\quad
\rho_e=(1+\eta^2/4)^{-1/2},\ R\le\eta\le1,
\label{eq:exact-sample-circuit}
\end{align}
where in the first line $y_j$ are the last centering Picard iterates with
anchors $(y,G)$ using \eqref{eq:kernel-approximate-velocity}, and the
second uses \eqref{eq:probability-update}. Neither includes an extra
Gaussian; the noise and the position anchor are allocated later.
A \emph{mean-dependency path} follows a mean node to a mean value used in
its center expression, repeating. Direct expression size counts
occurrences in one expression, before expanding center expressions.
\begin{theorem}[Numerical circuit interface]\label{thm:numerical-interface}\label{lem:compiler-circuit}
Fix $J\ge1$, the parameters \eqref{eq:numerical-parameters}, and
$N=\lceil10(J+2)\rceil$, and choose the interpolation order by
\eqref{eq:interpolation-order-choice}. For all sufficiently large $d$:
\begin{enumerate}[label=(\roman*),leftmargin=*]
\item \emph{Deterministic scalar formulas.} The exact mean nodes can be
ordered $M_1,\ldots,M_Q$ with
\begin{equation}\label{eq:exact-affine-centers}
M_i=u_{E_i,r_i},\quad E_i=b_i+\sum_{j<i}a_{ij}M_j,\quad
E_f=b_f+\sum_ja_{fj}M_j,
\end{equation}
all $a_{ij},a_{fj}$ deterministic scalars formable without any mean or
query.
\item \emph{Direct terms and scales.} $r_i=r\tau_i$,
$\tau_i=\sqrt{1-\rho_i^2}\ge R/\sqrt5$, $b_i=x+r(\ell_iy+m_iG)$. For
$\mathcal C_S$ the $G$ term is absent, $\rho_i\le\rho_e$,
$(\ell_i,m_i)=(\rho_i,0)$. For $\mathcal C_M$, $\rho_i\le b_*$ and the
pair is $(\rho_i,0)$ or $(\rho_i\cos z,\rho_i\sin z)$, $|z|\le C_m\psi$;
both coefficients have absolute value at most one. A node for $M_0=u$
has center $x$, scale $r$. Every mean occurring directly in the final
$\mathcal C_M$ expression is $M_*$ at a computed position, scale exactly
$rR$. Final direct terms are $b_f=\sigma G/2$ and $b_f=y/\rho_e$.
\item \emph{Coefficient sums and direct expression size.} There are
$A,A_f>0$ with $\sum_{j<i}|a_{ij}|\le A$, $\sum_j|a_{fj}|\le A_f$:
\begin{equation}\label{eq:coefficient-bounds}
\begin{gathered}
A_M=C_1(1+\log d)^{C_1}\bigl[(r/\sigma)\delta^{-4t-c_1w}+r^2\bigr],\quad
A_S=C_1r^2,\\
A_{f,M}=C_1,\quad A_{f,S}=C_1r,\quad c_1=10,
\end{gathered}
\end{equation}
$C_1\ge1$ depending only on fixed parameters and orders. Each direct
expression has at most $B_d\delta^{-c_1w}$ terms, and every mean-dependency
path length is at most $D_0=3N^2+1$, independent of $w,m$, cell count.
\item \emph{Errors on the ideal anchor laws.} Under
\eqref{eq:numerical-size}, for independent $Y\sim\xi_*$, $G\sim N(0,I_d)$,
\begin{equation}\label{eq:ideal-mean-circuit-error}
\Wtwo\bigl(\mathcal L(\mathcal C_M(x,r,\sigma;Y,G)),N(u,\sigma^2I_d/4)\bigr)
\le B_d\sigma\delta^JD(x),
\end{equation}
and (only $l\le B_d\delta^b$, $b\ge1/2$, needed) with
$Y\sim N(0,I_d)$,
\begin{equation}\label{eq:ideal-sample-circuit-error}
\Wtwo\bigl(\mathcal L(\mathcal C_S(x,r,\eta;Y)),\nu_{x,r}*N(0,\eta^2I_d/4)\bigr)
\le B_d\delta^JD(x).
\end{equation}
The coefficient and dependency assertions are deterministic for every
anchor choice; the distributional estimates need the ideal anchor laws.
\end{enumerate}
\end{theorem}
\begin{proof}
Expand only the additive state formulas, stopping at any earlier mean
value. This keeps coefficient bounds local.
\emph{Transport formulas.} A probability Picard state has its input
position as direct term and mean coefficients of absolute sum at most
$C_mr$; a harmonic state has direct term $\cos z$ times its input plus
$\sin z$ times its momentum, mean coefficients of absolute sum at most
$C_mr$ (sine-weighted rows bounded because integration intervals are
bounded). Both hold at every iteration, the additive direct term always
being the original input. A kernel evaluation, \eqref{eq:kernel-scalar-expansion},
has absolute coefficient sum at most $B_d\sigma^{-1}\delta^{-4t-w}$; the
momentum formula $(M_*-M_0)/s$ has sum at most $C/\sigma$; a centering
Picard update has the anchor as direct term and bounded outer rows.
\emph{Centers and coefficients.} A mean evaluation's physical center is
$x+r\rho v$ for normalized position $v$ by \eqref{eq:conditional-tilt}.
Substituting the direct starting states (backward probability at the
centering position, harmonic from it and the momentum, forward probability
from the harmonic position) and stopping at the means just described, only
free terms pass through, so anchor coefficients are $(\rho,0)$ or
$(\rho\cos z,\rho\sin z)$; for sampling only the first probability stage
occurs, so $(\rho,0)$. Multiplying by $r\rho$ ($\rho\le1$) bounds the
center coefficient sum by
$B_d(r/\sigma)\delta^{-4t-w}+C_mr^2$ in the mean case (the $C_mr^2$ from
the backward, harmonic, forward corrections) and by $C_mr^2$ in the
sampling case; enlarging to \eqref{eq:coefficient-bounds} is convenient.
The mean output is the single quadrature of \eqref{eq:exact-mean-circuit},
coefficient sum at most $C_m$, means all at time $b_*$; dividing a
probability correction by $\rho_e\ge2/\sqrt5$ gives output sum at most
$C_mr$.
The scale at time $\rho$ is $r\sqrt{1-\rho^2}$; in the mean computation
$\rho\le b_*$ gives $\tau\ge R$ (equality for output means), and in
sampling $\sqrt{1-\rho_e^2}=\eta/\sqrt{4+\eta^2}\ge R/\sqrt5$.
\emph{Expression size and depth.} Let $P\le B_d\delta^{-w}$ bound nodes
in one probability mesh; angle, harmonic, outer centering quadratures have
fixed node counts once $m$ is fixed. Formula
\eqref{eq:kernel-scalar-expansion} has at most $C_mP^2$ mean terms (one
factor $P$ from time-derivative quadrature, one from the forward integral).
A direct centering position has the same bound, its momentum a fixed number.
Substituting free states into a backward, harmonic, or forward center adds
at most one such list of each; no substitution continues through a mean
value. So $B_d\delta^{-2w}$ already bounds direct expression size and
$c_1w=10w$ is a uniform conservative choice.
Picard updates at a fixed layer depend only on earlier layers. Order an
internal kernel computation as backward transport, harmonic motion,
forward transport, after its input states: a mean path adds at most $3N$
nodes over a path in those inputs; at most $N$ outer Picard layers and one
final mean layer give depth $3N^2+1$ (sampling: at most $N$). This gives
an acyclic order for \eqref{eq:exact-affine-centers}; meshes, weights,
angles, and rows depend only on the prescribed scalars.
\emph{Ideal-law errors.} Couple the numerical mean formula and the exact
centering path integral by the same ideal anchors; Corollary~\ref{prop:exact-mean-circuits}
bounds their RMS difference by $B_d\sigma\delta^JD(x)$ while
Proposition~\ref{prop:centering-law} identifies the exact output as
$N(u,\sigma^2I_d/4)$, proving \eqref{eq:ideal-mean-circuit-error}. For
sampling the exact transport from a standard Gaussian has law
$\xi_{\rho_e}$; after division by $\rho_e$ it is
$Z+(\sqrt{1-\rho_e^2}/\rho_e)G=Z+(\eta/2)G$ with $Z\sim\nu_{x,r}$,
$G\sim N(0,I_d)$ independent, and $\rho_e\ge2/\sqrt5$ plus the same-anchor
RMS estimate prove \eqref{eq:ideal-sample-circuit-error}.
\end{proof}
For the implementation use $J=K_0+2$. The theorem fixes $N$ and $D_0$
before $w$; the implementation then chooses $w$ from its recursion count
and finally $m$ by \eqref{eq:interpolation-order-choice}. No random
failure event is introduced: errors are RMS estimates under ideal laws,
and the algebraic structure holds for every input.
A deterministic formula that calls for a conditional mean at a center
built from other conditional means is awkward: its error changes the next
center, and the next call is analyzed at a random center. A prescribed
Gaussian is instead added to every computed center. A translation
identity of the lower-level routine absorbs that Gaussian into its random
seeds, restoring the law of a call at the intended center, which is fixed
after conditioning on the parent's random anchors.
The identity is proved first, then the full recursive construction
(termination, per-execution query cap, distributional error). The
analytical input is the family of affine formulas in exact means supplied
by Theorem~\ref{thm:numerical-interface}. These are implemented below,
proving Theorem~\ref{thm:primitive}.
A seeded mean routine at scale $r>0$ is a measurable deterministic map
$\mathsf Q(x;z_1,\ldots,z_m)\in\R^d$ of a center $x\in\R^d$ and $m$ seed
slots $z_j\in\R^d$, possibly querying a fixed oracle. A vector
$v=(v_1,\ldots,v_m)\in\R^m$ gives a \emph{query-preserving translation} if
for every center, displacement, and fixed seeds,
\begin{equation}\label{eq:mean-shift-warmup}
\mathsf Q(x+\Delta;z+v\otimes\Delta/r)=\mathsf Q(x;z),
\end{equation}
with $(v\otimes a)_j=v_ja$, and the two computations ask exactly the same
physical queries. This is deterministic; no independence or Gaussianity
of the seeds is used.
\begin{lemma}[Absorbing a Gaussian center displacement]\label{lem:compiler-absorption}
Let $D_s>0$ and $\mathsf Q$ satisfy \eqref{eq:mean-shift-warmup} with
$\|v\|_2\le D_s$. Put $n=r/(2D_s)$ and
$P=(I_m-vv^T/(4D_s^2))^{1/2}$. For deterministic $x$, with independent
standard Gaussian $Z,G_1,\ldots,G_m$, the output of
$\mathsf Q(x+nZ;(P\otimes I_d)G)$ has exactly the law of $\mathsf Q(x;G)$.
\end{lemma}
\begin{proof}
The eigenvalues of $P$ lie in $[\sqrt3/2,1]$. The translation identity
rewrites the first output as
$\mathsf Q(x;PG-v\otimes Z/(2D_s))$ (suppressing $\otimes I_d$). This seed
array is centered Gaussian with slot covariance
$PP^T+vv^T/(4D_s^2)=I_m$, so its slots are independent standard Gaussians;
they may correlate with $Z$, harmless because the rewritten call is at the
deterministic center $x$.
\end{proof}
The covariance reduction $P$ is essential: with unmodified fresh seeds,
translating the noisy center back would add variance, whereas with $P$ the
center noise restores exactly the variance removed from the seeds. The
construction maintains the deterministic translation identity at every
level, so this applies recursively.
Retain the earlier notation: $g=\nabla F$,
$\Lip(g)\le\lambda=\delta=d^{-e}$, $0<r\le1$, $\nu_{x,r}$ as in
\eqref{eq:primitive-law}, $u_{x,r}=\E_{\nu_{x,r}}g(x+rZ)$, $L=\lambda r$,
$l=\lambda r^2$, $D(x)=\sqrt d+|g(x)|$, with the moment bounds
\eqref{eq:primitive-moments} and the coupling rules
Lemma~\ref{lem:coupling-smoothing}. Only evaluations of $g$ are used.
Fix $K_0,t$ as in Theorem~\ref{thm:primitive} and the numerical parameters
$R=\delta^t$, $b_*=(1-R^2)^{1/2}$, $\psi=\delta^{4t}$, $J=K_0+2$.
Theorem~\ref{thm:numerical-interface} supplies finite acyclic circuits
$\mathcal C_M,\mathcal C_S$ with exact conditional means as nodes; use its
$A$ (center coefficient bound), $A_f$ (output bound), $D_0$ (dependency
depth), and $c_1$ (fan-out exponent). A center has the form
$x+r(\ell_iy+m_iG)+\sum_{j<i}a_{ij}u_{x_j,r_j}$ with $r_i=r\tau_i$; all
scalar data are deterministic functions of the prescribed parameters,
independent of centers, anchors, and oracle answers.
For the mean circuit with parent noise $\sigma$, the underlying centering
flow uses $s=\sigma/2$ once, so its direct output is $(\sigma/2)G$, with
ideal-anchor error $\Wtwo(\mathcal L(\mathcal C_M),N(u_{x,r},\sigma^2I_d/4))
\le\polylog(d)\sigma\delta^JD(x)$. For the sampling circuit the direct
output is $y/\rho_e$, $\rho_e=(1+\eta^2/4)^{-1/2}$, ideal target
$\nu_{x,r}*N(0,\eta^2I_d/4)$. These need the anchor laws and size
assumptions of Theorem~\ref{thm:numerical-interface}, verified at their
uses.
The remaining task: replace every exact mean by finitely many queries while
preserving enough independent Gaussian noise for absorption. The
construction first establishes a query cap on every execution, then uses
pathwise stability at inaccurate centers, then absorption at ideal noisy
centers, then a standalone law for a child call.
Construct a mean routine $\mathsf M(x,r,\sigma;b,p)$ targeting
$u_{x,r}+\sigma G$ and a sampling routine $\mathsf S(x,r,\eta;b,p)$
targeting $Z+\eta G$ ($Z\sim\nu_{x,r}$, $G$ standard Gaussian). The top
call is $\mathsf S(x,r,1;1,K_0)$. A \emph{reachable call} is any descendant
of this call including itself. The first label $b$ records size
conditions
\begin{equation}\label{eq:routine-sizes}
\max\{l,L/\sigma\}\le\polylog(d)\delta^b\ (\mathsf M),\qquad
l/\eta\le\polylog(d)\delta^b\ (\mathsf S),
\end{equation}
$R\le\eta\le1$. The second label $p$ records the requested error power,
intended errors $\polylog(d)\sigma\delta^pD(x)$ and
$\polylog(d)\delta^pD(x)$. A call is terminal iff $p\le b$, nonterminal
otherwise. Every sampling routine reserves a final noise $\eta G/2$ unused
in its own queries; before adding it, its target is
$\nu_{x,r}*N(0,3\eta^2I_d/4)$.
For a nonterminal sampling call, use $\mathcal C_S$ with a standard
Gaussian position anchor, supply additional final noise $n=\eta/\sqrt2$,
then add the reserve. For a nonterminal mean call, pass the parent
$\sigma$ to $\mathcal C_M$ and supply final noise $n=\sqrt3\sigma/2$;
obtain its position anchor from an $\mathsf S$ call at the same center and
scale with $\eta=R/b_*$, multiplying that output by $b_*$ (an ideal such
anchor has law $b_*Z+RN$); its momentum anchor is an independent standard
Gaussian. For large $d$, $R/b_*\in[R,1]$, so the anchor call has an
admissible noise parameter.
An expression dependency and a routine call play different roles: if a
parent circuit contains $E_f=a\,u_{x_2,r_2}$ and $x_2=x+c\,u_{x_1,r_1}$,
evaluating $E_f$ invokes two new mean routines, both children of the same
parent routine, whose own internal circuits begin the next routine level.
The two nested affine expressions do not create two successive losses in
the parent's accuracy label.
To evaluate an expression, expand its mean dependencies as an uncached
tree, duplicating every mean occurrence with its entire center-supplying
subtree. All mean calls at any expression tier are direct children of the
routine whose circuit is being evaluated. Children of a parent $(b,p)$:
\begin{equation}
\begin{array}{c|c}
\text{child's role}&\text{child labels}\\ \hline
\text{term of }E_f\text{ in an }\mathsf M\text{ parent}&(b+t,p)\\
\text{term of }E_f\text{ in an }\mathsf S\text{ parent}&(b,p)\\
\text{any other term, or the }\mathsf M\text{ position anchor}
&(b-10t,p-1/4).
\end{array}
\label{eq:compiler-labels}
\end{equation}
In the last row labels are those of the parent routine, not successively
decreased while descending its expression tree.
\begin{lemma}[Bounded routine height]\label{lem:compiler-height}
Starting at $(b,p)=(1,K_0)$, \eqref{eq:compiler-labels} gives a finite
routine tree of height at most $H=H(K_0,t)$, and every reachable label has
$b>1/2$.
\end{lemma}
\begin{proof}
Before $4K_0+1$ decreasing transitions the 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$ forces termination, so no branch has more.
Positive transitions raise $b$ by $t$, while a nonterminal has
$b<p\le K_0$, bounding their total by $K_0,t$ and the bounded decrease of
$b$. A neutral transition is $\mathsf S\to\mathsf M$, and another needs a
fresh $\mathsf S$, reached only through a decreasing one; this bounds
neutral transitions too.
\end{proof}
A routine is a deterministic function of its center, scalar parameters,
oracle, and a finite ordered list of $d$-dimensional seed slots; standalone
use fills slots with independent standard Gaussians. Each routine has a
deterministic vector $v$ (one scalar per slot) with
\begin{equation}\label{eq:shift-budget}
\|v\|_2\le D_s:=\delta^{-t}.
\end{equation}
For every fixed seed list, changing $x$ to $x+\Delta$ and adding
$v\otimes\Delta/r$ to the seeds preserves every physical query point: the
output of $\mathsf M$ is unchanged, the normalized output of $\mathsf S$
decreases by $\Delta/r$. These are pathwise properties.
Within each call write $u=u_{x,r}$, $\nu=\nu_{x,r}$. At a terminal call
\begin{align*}
\mathsf M&=g(x+rZ)+\sigma G,\\
\mathsf S&=Z-rg(x+rZ)+(\sqrt3\eta/2)G_1+(\eta/2)G_2,
\end{align*}
separate seed slots, the last term of $\mathsf S$ its reserved noise.
Taking $v=-1$ on the $Z$ slot proves the shift property with $\|v\|_2=1$.
Fix a parent circuit and its anchors, write $M_i=u_{x_i,r_i}$ for its
exact mean nodes. For each expression $E=b_0+\sum_ia_iM_i$ the operation
$\operatorname{Eval}(E,n)$ approximates its exact value plus independent
$N(0,n^2I_d)$ noise. Assign disjoint seed blocks to every term
occurrence, its center-supplying computation, and each residual noise;
these are independent of the anchors, so term blocks are independent
conditional on the anchors in a standalone call.
Recursive rule (children and their shift vectors already built). Take
$A'=A$ for a center expression and $A'=A_f$ at the output, put
$\sigma_i=n/(2A')$. For the term of scale $r_i=r\tau_i$, first evaluate
its center expression with noise
\begin{equation}\label{eq:compiler-center-noise}
n_i=r_i/(2D_s).
\end{equation}
At this supplied center run the child $\mathsf M$ with scale $r_i$, noise
$\sigma_i$, labels from \eqref{eq:compiler-labels}. If its shift vector is
$v_i$, transform its fresh seed block $G_i$ by
\begin{equation}\label{eq:compiler-seed-transform}
P_i=(I-v_iv_i^T/(4D_s^2))^{1/2},\qquad\text{seeds }(P_i\otimes I_d)G_i,
\end{equation}
and denote its output $Q_i$. Return
\begin{equation}\label{eq:compiler-eval}
b_0+\sum_ia_iQ_i+\Bigl(n^2-\sum_ia_i^2\sigma_i^2\Bigr)^{1/2}G_E,
\end{equation}
$G_E$ a separate slot; the square root is defined since
$\sum_ia_i^2\sigma_i^2\le(\sum_i|a_i|)^2n^2/(4(A')^2)\le n^2/4$. With no
terms return $b_0+nG_E$. The mean-node ordering makes the recursion
finite; at the parent's output use the specified $n$, adding the reserved
noise for $\mathsf S$. Correlated seeds perform the same deterministic
operations.
At an ideal supplied center $x_i+n_iZ$, Lemma~\ref{lem:compiler-absorption}
identifies this child with its standalone law at $x_i$; the shift vectors
and the accuracy of the noisy centers are built next.
The shift vectors of Lemma~\ref{lem:compiler-absorption} are shown to exist
with the universal bound \eqref{eq:shift-budget}, which also makes the
definition noncircular and yields a deterministic query cap.
For a nonterminal $\mathsf S$ shift its Gaussian position anchor by
$-\rho_e\Delta/r$. For a nonterminal $\mathsf M$ its position anchor
$y_*$ must shift by $-\Delta/(b_*r)$ with momentum anchor fixed: apply
the shift of its anchor $\mathsf S$ call and additionally shift that
call's reserved seed by $-2R/b_*$ times $\Delta/r$. The ordinary shift
changes its output by $-\Delta/r$; since the reserved standard deviation
is $R/(2b_*)$, the extra change is $-R^2\Delta/(b_*^2r)$, and multiplying
by $b_*$ gives
$-(b_*+R^2/b_*)\Delta/r=-\Delta/(b_*r)$. The extra slot is used in no
query, so query invariance is preserved.
For a term of scale $r\tau_i$, the direct part of its supplied center then
changes by $k_i\Delta$ with
$k_i=1-\ell_i\rho_e$ ($\mathsf S$) or $k_i=1-\ell_i/b_*$ ($\mathsf M$),
satisfying $|k_i|/\tau_i\le C$: for $\mathsf S$, $\rho_i\le\rho_e$ gives
$1-\rho_i\rho_e\le2(1-\rho_i)\le2\tau_i$; for $\mathsf M$,
$0\le1-\rho_i\cos z/b_*\le C(1-\rho_i)+Cz^2$ with $z=0$ when there is no
cosine, using $1-\rho_i\le\tau_i$, $\tau_i\ge cR$, and
$z^2/R\le C\delta^{7t}\to0$. Shift that child's untransformed seed block
by
\begin{equation}\label{eq:compiler-child-shift}
\frac{k_i}{\tau_i}P_i^{-1}v_i\otimes\frac{\Delta}{r}.
\end{equation}
If \eqref{eq:shift-budget} holds for the child, the least eigenvalue of
$P_i$ is at least $\sqrt3/2$, so $\|P_i^{-1}\|_{\op}\le2$; the formula
induces exactly the child's shift for center displacement $k_i\Delta$.
Leave residual noise slots fixed. Induction from leaf expressions upward
shows every mean output is unchanged, so the actual supplied-center change
is its direct-term change; all queries are preserved. The parent's final
direct term is unchanged for $\mathsf M$, changes by $-\Delta/r$ for
$\mathsf S$. The center-supplying subtree and the invoked child have
disjoint slots, so $k_i/\tau_i$ multiplies only that child's shift.
\begin{lemma}[Well-defined construction and query cap]\label{lem:compiler-count}
One can choose $w>0$, depending only on $K_0,t$ and the fixed Picard
counts, and then the interpolation orders, so the whole construction is
well defined for all large $d$; every routine reached from a top call
satisfies \eqref{eq:shift-budget}; and its total number of seed slots and
oracle queries is at most $\polylog(d)\delta^{-t}$ on every execution.
\end{lemma}
\begin{proof}
Let $D_0$ be the dependency-depth bound and $F_0$ the direct fan-out bound
of Theorem~\ref{thm:numerical-interface}. Expanding a circuit into an
expression tree costs at most $(D_0+1)(1+F_0)^{D_0}$ nodes, hence
$\polylog(d)\delta^{-C_2w}$ slots with $C_2$ independent of $w$ and the
orders; this bounds its direct child invocations and its own anchor and
residual slots. Expanding through the height $H$ of
Lemma~\ref{lem:compiler-height} gives $\polylog(d)\delta^{-C_3w}$ slots
and terminal queries. Interpolation orders change only fixed constants.
The terminal shift has norm $1$. Given the bound for the children, the
bounded ratios and inverse norms show a parent's shift norm is at most a
fixed constant times one plus the sum of its direct children's, including
the anchor call (the extra reserved-slot shift has bounded coefficient;
the triangle inequality covers coincident ordinary and extra shifts).
Iterating through the fixed height gives the same
$\polylog(d)\delta^{-C_3w}$.
Choose $w>0$ with $C_3w<t/2$, $w<t$, $c_1w<t$, then the fixed
interpolation orders of Theorem~\ref{thm:numerical-interface}. For large
$d$ the last shift bound is below $D_s=\delta^{-t}$; this is a valid
bottom-up induction, proving each child's bound before using its inverse
in the parent, and the count is at most $\polylog(d)\delta^{-t}$.
The parameter order: fix Picard counts (hence $D_0,H,C_3$), then $w$, then
orders. $D_s$ fixes every center-noise budget; circuit shapes, labels,
noises, and slot counts are determined by these scalars before any shift.
Compute child shifts before their matrices and the parent shift. No center
value, oracle answer, error estimate, or \eqref{eq:routine-sizes} enters
this construction or the count; only terminal routines query $g$, once
each, so it is a cap on every execution.
\end{proof}
The finite construction and its cap are now independent of accuracy; it remains
to verify terminal errors and the size conditions.
\begin{lemma}[Terminal accuracy]\label{lem:compiler-terminal}
In standalone use a terminal call satisfying \eqref{eq:routine-sizes}
obeys
\begin{align*}
\Wtwo(\mathcal L(\mathsf M),N(u,\sigma^2I_d))
&\le\polylog(d)\sigma\delta^pD(x),\\
\Wtwo(\mathcal L(\mathsf S_{\rm pre}),\nu_{x,r}*N(0,3\eta^2I_d/4))
&\le\polylog(d)\delta^pD(x),
\end{align*}
$\mathsf S_{\rm pre}$ the output before adding $\eta G_2/2$; the sampling
bound also holds after adding that reserve to both laws. The fixed-seed
center Lipschitz constant is at most $\lambda$ for $\mathsf M$ and
$r\lambda$ for $\mathsf S$.
\end{lemma}
\begin{proof}
For standard Gaussian $Z$, the moment bounds give
$\|g(x+rZ)-u\|_{L^2}\le CL(\sqrt d+r|g(x)|)\le CLD(x)$; since
$L/\sigma\le\polylog(d)\delta^b$ and $p\le b$ this is within budget, and the
added Gaussian is coupled identically.
Let $z_0$ be the mode, $z_0=-rg(x+rz_0)$; then
$\|Z-rg(x+rZ)-(Z+z_0)\|_{L^2}\le l\|Z-z_0\|_{L^2}\le ClD(x)$ and
$\Wtwo(N(z_0,I_d),\nu_{x,r})\le Cl\sqrt d$. Apply the synchronous
comparison \eqref{eq:synchronous-comparison} with Gaussian potential
$H_0=|z-z_0|^2/2$ (contraction constant $1$) and $H=|z|^2/2+F(x+rz)$,
whose gradient difference is $r[g(x+rz)-g(x+rz_0)]$, of
$L^2(\nu_{x,r})$ norm at most $l\|Z-z_0\|\le Cl\sqrt d$. Adding the same
comparison noises proves $ClD(x)\le\polylog(d)\eta\delta^bD(x)\le$ the
requested bound ($\eta\le1$, $p\le b$); omitting $G_2$ gives the
pre-reserved assertion. Lipschitz follows from the two terminal formulas.
\end{proof}
\begin{lemma}[Size conditions descend]\label{lem:compiler-sizes}
Starting from a top call satisfying \eqref{eq:routine-sizes}, every
descendant satisfies it with its labels up to a uniform polylogarithmic
factor; at each nonterminal circuit
\begin{equation}\label{eq:small-coefficients}
\lambda A\le\polylog(d)\delta^{1/3}.
\end{equation}
For nonfinal expression terms, $\sigma_i\le C\sigma$ in an $\mathsf M$
parent, $r\sigma_i\le C$ in an $\mathsf S$ parent.
\end{lemma}
\begin{proof}
A final term of an $\mathsf M$ parent uses means of scale $rR$, so
$r_i=rR$ while $\sigma_i$ is a fixed positive multiple of $\sigma$;
$L_i/\sigma_i$ gains a factor $R$ and $l_i$ gains $R^2$, allowing
$(b+t,p)$. A final term of an $\mathsf S$ parent has $\sigma_i$ a fixed
positive multiple of $\eta/r$, so $L_i/\sigma_i\le Cl/\eta$, $l_i\le l$,
allowing $(b,p)$.
For a nonfinal term, let $r\tau_{\rm prev}$ be the scale of the node whose
center is being supplied; its expression is evaluated with budget
$r\tau_{\rm prev}/(2D_s)$ regardless of tier, so
$\sigma_i=r\tau_{\rm prev}/(4D_sA)$ and
$L_i/\sigma_i\le CD_s\lambda A/R$. The coefficient estimates give
$\lambda A\le\polylog(d)\delta^{b-4t-c_1w}$ for an $\mathsf M$ parent and
$\lambda A=C_1l\le\polylog(d)\delta^b$ for $\mathsf S$ ($\eta\le1$); with
$D_s/R=\delta^{-2t}$ and $c_1w<t$ this costs under the allowed $10t$
decrease in $b$, and $l_i\le l$. For the position-anchor $\mathsf S$ call
of $\mathsf M$, $l/\eta\le l/R$ has the same allowance. This proves the
size conditions by the bounded height, all factors still one polylog.
Since $b>1/2$ and $4t+c_1w<5t<1/6$, the coefficient bounds give
\eqref{eq:small-coefficients}. The positive $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$.
\end{proof}
\begin{lemma}[Fixed-seed stability]\label{lem:compiler-lipschitz}
For every seed array, including correlated or transformed arrays, the
compiled $\mathsf M$ output has $x$-Lipschitz constant at most
$\polylog(d)\lambda$; the normalized $\mathsf S$ output at most
$\polylog(d)r\lambda$. Treating the anchors of a nonterminal $\mathsf M$
as separate inputs, the constant in the position anchor is at most
$\polylog(d)L$.
\end{lemma}
\begin{proof}
Induct on routine height; the terminal case is
Lemma~\ref{lem:compiler-terminal}. Hold anchors fixed. A direct center
term has 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 bound, earlier supplied-center
sensitivities contribute at most $\polylog(d)\lambda A$ times their
maximum, small by \eqref{eq:small-coefficients}. Induction through the
bounded expression depth bounds center sensitivities by a constant times
their direct ones with a polylog factor.
At the final expression the coefficient sum is $A_f$; the direct term has
no center dependence with anchors fixed and (for $\mathsf M$) no
position-anchor dependence, so center sensitivity is at most
$\polylog(d)A_f\lambda$ and $\mathsf M$ position-anchor sensitivity at most
$\polylog(d)A_f\lambda r$. The $\mathsf S$ Gaussian anchor is center
independent. The $\mathsf M$ position anchor, supplied by $\mathsf S$ at
the same center, contributes at most $\polylog(d)Lr\lambda=\polylog(d)L^2
\le\polylog(d)\lambda$. All seed transformations depend only on scalars,
not $x$; no seed independence is needed.
\end{proof}
The law guarantees are proved by induction on height, the only stage
assuming independent standard Gaussian seeds, used for a child only after
absorption restores its standalone seed law.
\begin{proposition}[Accuracy of the compiled routines]\label{prop:compiler-accuracy}
For every call reached from the top, the standalone routines satisfy
\begin{align*}
\Wtwo(\mathcal L(\mathsf M),N(u,\sigma^2I_d))
&\le\polylog(d)\sigma\delta^pD(x),\\
\Wtwo(\mathcal L(\mathsf S_{\rm pre}),\nu_{x,r}*N(0,3\eta^2I_d/4))
&\le\polylog(d)\delta^pD(x),
\end{align*}
$\mathsf S_{\rm pre}$ before its reserved noise; the latter after adding
the reserve too, and
$\|\mathcal L(\mathsf S)-\nu_{x,r}*N(0,\eta^2I_d)\|_{\TV}
\le\polylog(d)\delta^pD(x)/\eta$.
\end{proposition}
\begin{proof}
Terminal case: Lemma~\ref{lem:compiler-terminal}. All variables below
have finite second moments (actual algorithms are finite compositions of
affine maps and globally Lipschitz evaluations; ideal variables obey the
primitive moment bounds), so the $\Wtwo$ comparisons are defined.
For a nonterminal $\mathsf M$, first couple its approximate $\mathsf S$
position anchor to an ideal one, independently of the remaining blocks.
Child labels are $(b-10t,p-1/4)$; by induction and
Lemma~\ref{lem:compiler-lipschitz} this costs at most
$\polylog(d)L\delta^{p-1/4}D(x)$, i.e. after dividing by the mean unit
$\sigma\delta^pD(x)$, at most $\polylog(d)\delta^{b-1/4}$, bounded for
large $d$. So the position anchor may be assumed to have law $b_*Z+RN$
($Z\sim\nu_{x,r}$) and the momentum anchor is independent standard
Gaussian, the inputs in Theorem~\ref{thm:numerical-interface}.
Let $x_i$ be the centers of the deterministic exact-mean circuit under
these anchors and $D_i=\sqrt d+|g(x_i)|$. Their RMS values over anchors
obey
\begin{equation}\label{eq:compiler-center-moments}
\max_i\|D_i\|_{L^2}\le\polylog(d)D(x).
\end{equation}
Both anchors have RMS at most $CD(x)$; $|g(x_i)|\le|g(x)|+\lambda|x_i-x|$
with the circuit's direct terms, a preceding exact mean of magnitude at
most $CD_j$ contributing at most $C\lambda A$ times the previous maximum,
small by \eqref{eq:small-coefficients}; indirect terms cost at most
$CD(x)$. Induction proves \eqref{eq:compiler-center-moments} without a
maximum inside the expectation.
Condition on the anchors. For expression $E$ let $\epsilon_E$ be the
$\Wtwo$ error of $\operatorname{Eval}(E,n)$ relative to its exact value plus
$N(0,n^2I_d)$. Couple the supplied center of its $i$th term to $x_i+n_iZ$
with error $\epsilon_{E_i}$, using the same transformed child seeds drawn
independently; the fixed-seed Lipschitz bound costs at most
$\polylog(d)\lambda\epsilon_{E_i}$. At the ideal noisy center,
Lemma~\ref{lem:compiler-absorption} identifies the child's law with its
standalone law at $x_i$, and induction supplies
$\polylog(d)\sigma_i\delta^{p_i}D_i$.
Distinct occurrences have disjoint blocks, so their outputs are
independent conditional on the anchors; product couplings compare to
independent noised exact means, and with the independent residual in
\eqref{eq:compiler-eval} the ideal sum has covariance $n^2I_d$. Hence
\begin{equation}\label{eq:compiler-error-recursion}
\epsilon_E\le\sum_i|a_i|\bigl[\polylog(d)\lambda\epsilon_{E_i}
+\polylog(d)\sigma_i\delta^{p_i}D_i\bigr].
\end{equation}
Note the order: inaccurate centers use only pathwise stability; standalone
accuracy is invoked only after absorption; no claim is made that the
algorithm's intermediate states have ideal laws.
Take RMS in the anchors and expand \eqref{eq:compiler-error-recursion}
through the finite expression tree. At the root direct children have
$p_i=p$ and $\sum_i|a_i|\sigma_i\le A_f n/(2A_f)=n/2$, fitting the budget.
At tier $j\ge2$ the coefficient sums and Lipschitz factors contribute at
most $A_f(\polylog(d)\lambda A)^{j-1}$, all such children having the one
loss $p_i=p-1/4$; their noise units are at most $C\sigma$ ($\mathsf M$)
or $C/r$ ($\mathsf S$) by Lemma~\ref{lem:compiler-sizes}, combining with
$A_f=C_1$ or $C_1r$. Already
$\polylog(d)\lambda A\le\polylog(d)\delta^{1/3}$ pays for the loss
$\delta^{-1/4}$, leaving $\delta^{1/12}$; the bounded depth sums the rest.
Leaf affine expressions have zero error.
Conditionally on the anchors, the ideal final evaluation is the exact
circuit value plus a Gaussian of the specified covariance, taken
independent of the anchors. By Lemma~\ref{lem:compiler-sizes},
Theorem~\ref{thm:numerical-interface} 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 errors $\polylog(d)\sigma\delta^JD(x)$
and $\polylog(d)\delta^JD(x)$; since $J=K_0+2\ge p$, both fit. Couple the
evaluation noise identically: 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 convolution
from $\eta^2/4$ to $3\eta^2/4$. Adding the common reserve preserves the
$\Wtwo$ bound; smoothing by standard deviation $\eta/2$ gives the total
variation bound by \eqref{eq:gaussian-smoothing}.
\end{proof}
\begin{proof}[Proof of Theorem~\ref{thm:primitive}]
Use the top $\mathsf S$ call with $(b,p,\eta)=(1,K_0,1)$. Since
$l=\delta r^2\le\delta$, its size condition holds with constant $1$.
Proposition~\ref{prop:compiler-accuracy} gives $\Wtwo$ 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~\ref{lem:compiler-count} gives at
most $\polylog(d)\delta^{-t}$ evaluations of $g$ on every execution. The
output $X_{\rm pre}+G/2$ has $G$ a reserved standard Gaussian independent
of $X_{\rm pre}$ and all queries, and the pre-reserved conclusion gives
$\Wtwo(\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 comparison laws serve
only to analyze this finite algorithm, whose only oracle operations are the
terminal evaluations.
\end{proof}
The auxiliary algorithm returns a conditional sample with an added
Gaussian. It is used to sample the original unnoised law $\pi_V$. A
contracting proximal chain first approaches
$\widetilde\pi=\pi_V*N(0,hI_d)$; conditional on its output, variance-halving
steps reduce to a law close to a mode-centered Gaussian, which completes
the descent. Call numbers are prescribed in advance.
Throughout, $V\in C^2(\R^d)$ satisfies $V(0)=0$, $\nabla V(0)=0$, and
$I_d\preceq\nabla^2V\preceq2I_d$; the origin is the known minimizer. Each
query returns $(V,\nabla V)$; the construction uses only the gradient, the
fallback also the value. Fix $0<e<1$, $\delta=d^{-e}$, $h=\delta/2$, and
an integer $K_0>2+1/e$. Choose $t>0$ as in \eqref{eq:t-choice} and fix the
orders required by Theorem~\ref{thm:primitive}. All thresholds are uniform
in $V$.
The Gaussian proximal sampler alternates a Gaussian step with a
restricted Gaussian conditional draw
(Lee, Shen, and Tian, Definition 1 and Lemmas 2 and 4); see also Chen,
Chewi, Salim, and Wibisono (Section 3 and Theorem 1).
\url{https://arxiv.org/pdf/2010.03106v4}
\url{https://proceedings.mlr.press/v178/chen22c/chen22c.pdf}
For $y\in\R^d$ let $\kappa_y$ be
proportional to
\begin{equation}\label{eq:proximal-conditional}
\exp\!\left(-V(x')-\frac{|x'-y|^2}{2h}\right).
\end{equation}
The ideal transition $P_h$ draws $x'\sim\kappa_y$ and returns
$x'+\sqrt h\,N$, so $P_h(y,\cdot)=\kappa_y*N(0,hI_d)$.
\begin{lemma}[Proximal contraction]\label{lem:proximal-contraction}
$\widetilde\pi$ is invariant for $P_h$, and
$\Wtwo(\kappa_{y_1},\kappa_{y_2})\le|y_1-y_2|/(1+h)$. Consequently for
input laws $\mu,\eta$ of finite second moment,
$\Wtwo(\mu P_h,\eta P_h)\le\Wtwo(\mu,\eta)/(1+h)$.
\end{lemma}
\begin{proof}
Under $\pi_V(dx)N(x,hI_d)(dy)$ the conditional law of $x$ given $y$ is
$\kappa_y$; drawing it and redrawing $y$ given $x$ preserves the joint law,
whose $y$-marginal is $\widetilde\pi$.
The potential in \eqref{eq:proximal-conditional} has curvature at least
$1+1/h$; replacing $y_1$ by $y_2$ changes its gradient by the constant
$(y_1-y_2)/h$. The synchronous comparison \eqref{eq:synchronous-comparison}
bounds the conditional $\Wtwo$ by $|y_1-y_2|/(1+h)$; coupling the Gaussian
increments identically and integrating against an input coupling
(Lemma~\ref{lem:coupling-smoothing}) proves the contraction.
\end{proof}
To implement the transition substitute $x'=y+\sqrt h\,z$: the normalized
conditional potential becomes $|z|^2/2+F_y(z)$ with
\begin{equation}\label{eq:proximal-primitive}
F_y(z)=V(y+\sqrt h\,z),\qquad g_y(z)=\sqrt h\,\nabla V(y+\sqrt h\,z),
\end{equation}
and the Hessian bound with $\nabla V(0)=0$ give
$\Lip(g_y)\le2h=\delta$, $|g_y(0)|\le2\sqrt h\,|y|$. An evaluation of
$g_y$ is one original query. Apply Theorem~\ref{thm:primitive} to $F_y$ at
center zero, scale one, with fresh seeds, and map the output $S$ to
$y+\sqrt h\,S$; this defines the actual transition $\widehat P_h$, whose
ideal comparison law is exactly $P_h(y,\cdot)$, retaining the primitive's
reserved $G/2$. Its pre-reserve $\Wtwo$ guarantee and Gaussian smoothing
give
\begin{equation}\label{eq:proximal-local-errors}
\Wtwo(\widehat P_h(y,\cdot),P_h(y,\cdot))\le\sqrt h\,B_d\delta^{K_0}
(\sqrt d+2\sqrt h\,|y|),\quad
\norm{\widehat P_h-P_h}_{\TV}\le B_d\delta^{K_0}(\sqrt d+2\sqrt h|y|).
\end{equation}
The constant $B_d=C(1+\log d)^C$ may be increased by a fixed factor or power, its
constants depending only on the fixed parameters.
\begin{proposition}[Mixing of the approximate noised chain]\label{prop:noised-mixing}
With $K_0>2+1/e$ and $t$ as above, start the $\widehat P_h$-chain at zero
and run $\lceil Ch^{-1}\log d\rceil+1$ steps. For a sufficiently large
fixed $C$, its output has total-variation distance $o(1)$ from
$\widetilde\pi$, uniformly over admissible $V$.
\end{proposition}
\begin{proof}
Lemma~\ref{lem:convex-preliminaries} at the mode zero gives
$\E_{\pi_V}|X|^2\le d$, so the RMS under $\widetilde\pi$ is at most
$C\sqrt d$. Let $\mu_k$ be the actual law after $k$ steps and
$a_k=\Wtwo(\mu_k,\widetilde\pi)$; then $({\textstyle\int}|y|^2\mu_k(dy))^{1/2}
\le a_k+C\sqrt d$. Moments are finite: finite at $k=0$, and the local
$\Wtwo$ bound with $P_h$'s conditional moment bounds propagates them.
Taking RMS of the first local error, conditional coupling and the
contraction give
$a_{k+1}\le(1/(1+h)+2hB_d\delta^{K_0})a_k+C\sqrt h\,B_d\delta^{K_0}\sqrt d$.
For large $d$ the $a_k$ coefficient is at most $1-ch$, and $a_0\le C\sqrt d$,
so summing the geometric recurrence yields
\begin{equation}\label{eq:proximal-global-error}
a_k\le C\sqrt d\,e^{-chk}+Ch^{-1/2}B_d\delta^{K_0}\sqrt d.
\end{equation}
The second term is $o(\sqrt h)$ since after division by $\sqrt h$ it is
bounded by a constant times $B_d\delta^{K_0-1}\sqrt d\to0$ (because
$e(K_0-1)>1/2$); at $k=\lceil Ch^{-1}\log d\rceil$ a large fixed $C$ makes
the first term $o(\sqrt h)$ too.
One more transition for TV control. The mixture bound gives
$\norm{\mu_k\widehat P_h-\mu_kP_h}_{\TV}\le
B_d\delta^{K_0}[\sqrt d+2\sqrt h(a_k+C\sqrt d)]=o(1)$. The conditional
$x'$-laws from $\mu_k$ and $\widetilde\pi$ couple with RMS at most
$a_k/(1+h)$, and their common $\sqrt h$-Gaussian gives, by
\eqref{eq:gaussian-smoothing},
$\norm{\mu_kP_h-\widetilde\pi P_h}_{\TV}\le a_k/((1+h)\sqrt h)=o(1)$.
Invariance of $\widetilde\pi$ completes it.
\end{proof}
An exact draw from $\kappa_y$, $y\sim\widetilde\pi$, has marginal
$\pi_V$. This conditional draw is implemented; for analysis first replace
the incoming approximate $y$ by $y\sim\widetilde\pi$, costing at most the
TV error of Proposition~\ref{prop:noised-mixing} by data processing. For a
related restricted-Gaussian-oracle recursion see Chen, Chewi, Lu, and
Zhang (Algorithm 3.3 and Section 6.3);
\url{https://arxiv.org/pdf/2609.06906v1}
the scheme uses fixed variance
halving, mode-centered Gaussian termination, and an every-execution cap.
Fix $y$; keep $F_y,g_y$ throughout, changing only center and scale. At
center $x$, scale $r$, the desired law is
$\nu_{x,r}(dz)\propto\exp(-|z|^2/2-F_y(x+rz))dz$. If $Z\sim\nu_{x,r}$ and
$N\sim N(0,I_d)$, then given $Y=Z+N$ the conditional density of $Z$ is
$\propto\exp(-|z-Y/2|^2-F_y(x+rz))$, and $z=Y/2+z_{\rm new}/\sqrt2$ gives
\begin{equation}\label{eq:conditional-descent}
Z=\frac Y2+\frac{Z_{\rm new}}{\sqrt2},\qquad
Z_{\rm new}\mid Y\sim\nu_{x+rY/2,\,r/\sqrt2}.
\end{equation}
The algorithm starts at $x_0=0,r_0=1$ with $J=\lceil C_1\log d\rceil$
levels; at each $0\le j<J$ it calls Theorem~\ref{thm:primitive} with fresh
seeds at $(x_j,r_j)$ to get a noised observation $Y_j$, setting
$x_{j+1}=x_j+r_jY_j/2$, $r_{j+1}=r_j/\sqrt2$; at level $J$ it produces a
base sample $Z_J$ and returns $y+\sqrt h\,(x_J+r_JZ_J)$. With backward
variables $Z_j=Y_j/2+Z_{j+1}/\sqrt2$, the center update gives
$x_j+r_jZ_j=x_{j+1}+r_{j+1}Z_{j+1}$, so since $x_0=0,r_0=1$ the output is
$y+\sqrt h\,Z_0$; by \eqref{eq:conditional-descent}, ideal observations and
an exact base make it an exact draw from $\kappa_y$.
Errors at random centers are controlled using the ideal-history moments.
\begin{lemma}[Moments along the ideal descent]\label{lem:descent-moments}
If $y\sim\widetilde\pi$ and every conditioning observation preceding level
$j$ is ideal, then $b_j:=(\E|g_y(x_j)|^2)^{1/2}\le C\sqrt d$ for
$0\le j\le J$, with absolute $C$ for large $d$.
\end{lemma}
\begin{proof}
The initial bound follows from $|g_y(0)|\le2\sqrt h|y|$ and the moment of
$\widetilde\pi$. Conditional on an ideal history the next observation has
law $Z+N$ at center $x_j$, scale $r_j=2^{-j/2}$, so
\eqref{eq:primitive-moments} gives
$(\E[|Y_j|^2\mid x_j,y])^{1/2}\le C(\sqrt d+r_j|g_y(x_j)|)$. The center
update with $\Lip(g_y)\le\delta$ and the $L^2$ triangles give
$b_{j+1}\le(1+C\delta r_j^2)b_j+C\delta r_j\sqrt d$; $\sum_jr_j$ and
$\sum_jr_j^2$ are bounded, so $\prod_j(1+C\delta r_j^2)\le
\exp(C\delta\sum_jr_j^2)$, and expanding the recurrence proves the bound.
\end{proof}
For the TV sum, let $H_j$ ($0\le j\le J$) be the output law with the first
$j$ observations ideal, all later ones actual, base rule actual. For
$0\le j<J$ the laws $H_j,H_{j+1}$ share an ideal prefix of length $j$;
conditional on it only the next observation kernel changes, subsequent
computation being the same kernel, so data processing and the mixture TV
bound give
$\norm{H_j-H_{j+1}}_{\TV}\le B_d\delta^{K_0}\E_{\rm ideal}[\sqrt d+|g_y(x_j)|]$.
Lemma~\ref{lem:descent-moments} then yields
\begin{equation}\label{eq:descent-tv}
\norm{H_0-H_J}_{\TV}\le CJB_d\delta^{K_0}\sqrt d=o(1),
\end{equation}
summing errors using only ideal-history moments. It remains to compare the
actual and exact base laws under that history.
At depth $J$ write $l=\delta r_J^2$ and start $z^{(0)}=0$; make
$m_0=\lceil C_2\log d\rceil$ gradient iterations
$z^{(k+1)}=-r_Jg_y(x_J+r_Jz^{(k)})$, $0\le k<m_0$, and output
$Z_J=z^{(m_0)}+G$ for a fresh standard Gaussian $G$. The iteration has
Lipschitz constant at most $l$ and unique fixed point the mode $z_0$;
\eqref{eq:primitive-moments} gives
\begin{equation}\label{eq:mode-accuracy}
|z^{(m_0)}-z_0|\le l^{m_0}|z_0|\le l^{m_0}\frac{r_J|g_y(x_J)|}{1-l}.
\end{equation}
\begin{lemma}[Mode-centered Gaussian approximation]\label{lem:mode-gaussian}
For any $C^2$ $F$ with $\Lip(\nabla F)\le\lambda$, center $x$, scale
$0<r\le1$, $l=\lambda r^2\le1/2$, mode $z_0$ of $\nu_{x,r}$,
$\norm{\nu_{x,r}-N(z_0,I_d)}_{\TV}\le Cl\sqrt d$.
\end{lemma}
\begin{proof}
Translate by $z_0$; up to a constant the potential of $w=z-z_0$ is
$|w|^2/2+J(w)$ with $J(w)=z_0\cdot w+F(x+r(z_0+w))-F(x+rz_0)$. The mode
gives $\nabla J(0)=0$, $\Lip(\nabla J)\le l$, so $J(0)=0$,
$|J(w)|\le l|w|^2/2$. For $0\le s\le1$, $|w|^2/2+sJ(w)$ has mode zero and
Hessian between $(1-l)I$ and $(1+l)I$; its law $\mu_s$ satisfies by
Lemma~\ref{lem:convex-preliminaries}
$\Var_{\mu_s}J\le\frac1{1-l}\E_{\mu_s}|\nabla J|^2
\le\frac{l^2}{1-l}\E_{\mu_s}|w|^2\le\frac{l^2d}{(1-l)^2}\le4l^2d$. The
quadratic bounds justify differentiating
$\Phi(s)=\log\int e^{-|w|^2/2-sJ(w)}dw$ twice, giving
$\Phi''(s)=\Var_{\mu_s}J$; since $\mu_0=N(0,I_d)$,
$\KL(\mu_0\Vert\mu_1)=\int_0^1(1-s)\Phi''(s)ds\le2l^2d$. Pinsker, then
translating back by $z_0$, proves it.
\end{proof}
Under an ideal history, the Gaussian shift bound and
\eqref{eq:mode-accuracy} bound the expected base TV error by
\begin{equation}\label{eq:descent-base-error}
C\delta2^{-J}\sqrt d+C(\delta2^{-J})^{m_0}\frac{2^{-J/2}\sqrt d}{1-\delta2^{-J}}.
\end{equation}
The second term uses Lemma~\ref{lem:descent-moments}. Choose fixed $C_1$
so $\delta2^{-J}\sqrt d\to0$ and fixed $C_2$ so the second term $\to0$;
replacing the base law under the ideal history changes TV by $o(1)$. By
\eqref{eq:conditional-descent} the fully ideal conditional algorithm
returns $\kappa_y$, and averaging over $y\sim\widetilde\pi$ gives $\pi_V$.
Combining the base comparison, \eqref{eq:descent-tv}, and the initial
replacement proves
\begin{equation}\label{eq:unnoised-output-error}
\norm{\mathcal L(X_{\rm out})-\pi_V}_{\TV}=o(1),
\end{equation}
uniformly over admissible $V$, the law being the original target without an
extra convolution.
\begin{theorem}[Subpolynomial upper bound]\label{thm:upper}
For every fixed $\varepsilon>0$ there is $C_\varepsilon<\infty$ such that,
for every integer $d\ge2$, a measurable algorithm using the exact
value-and-gradient oracle makes at most $C_\varepsilon d^\varepsilon$
queries on every execution and returns a sample within TV distance $1/10$
of $\pi_V$ for every $V$ satisfying \eqref{eq:original-class}; thus
$Q(d)\le C_\varepsilon d^\varepsilon$.
\end{theorem}
\begin{proof}
Given $\varepsilon>0$, choose $e\in(0,1)$ with $2e<\varepsilon$, then an
integer $K_0>2+1/e$, then $t>0$ as in \eqref{eq:t-choice}; fix the orders
in prescribed order and finally the chain, descent, and base constants.
These may depend on $\varepsilon$ but are fixed before any oracle call.
Every invocation of Theorem~\ref{thm:primitive} uses at most
$\polylog(d)\delta^{-t}$ gradient evaluations, each one original query by
\eqref{eq:proximal-primitive}, even at random descent centers. There are
$O(\delta^{-1}\log d)$ proximal calls, $O(\log d)$ descent calls, and
$O(\log d)$ base queries, giving the deterministic cap
\begin{equation}\label{eq:upper-query-count}
\polylog(d)\delta^{-(1+t)}=d^{e(1+t)+o(1)}.
\end{equation}
Query counts and recursion shapes are fixed by the scalars, independent of
oracle magnitudes or approximation events; all maps are measurable by the
primitive theorem and the explicit finite operations. Since $t<1$,
$e(1+t)<2e<\varepsilon$, and polylogarithmic factors fit
$C_\varepsilon d^\varepsilon$ for large $d$;
\eqref{eq:unnoised-output-error} gives TV error at most $1/10$ above a
uniform threshold.
For each of the finitely many smaller dimensions $d$, use independent
standard Gaussian proposals, accepting $X$ with probability
$\exp(-V(X)+|X|^2/2)$. Integrating the Hessian bounds from the known
minimizer gives $|x|^2/2\le V(x)\le|x|^2$, so the acceptance probability is
at most one with average at least $\E e^{-|X|^2/2}=2^{-d/2}$; the proposal
times acceptance is proportional to $e^{-V(x)}$, so an accepted proposal
has law $\pi_V$. Cap trials at $\lceil4\cdot2^{d/2}\rceil$; return the
first accepted or zero. Failure probability is at most $e^{-4}<1/10$, so
TV error at most $e^{-4}$, with at most the prescribed number of queries
(one per trial). Increasing $C_\varepsilon$ covers these finitely many
caps, proving the assertion for all $d\ge2$.
\end{proof}
This proves part (a) of Theorem~\ref{thm:main}. The orders depend on the
fixed $\varepsilon$, and the result bounds oracle calls, not computation.
The independent logarithmic lower bound is proved subsequently.
The lower bound already holds for centered Gaussian targets with known
minimizer, so the obstruction is present in the same oracle model. The
block-Krylov strategy of
Chewi, de Dios Pont, Li, Lu, and Narayanan (Section 5.3.3 and Theorem 45) is adapted, giving a complete
argument for the spectral interval $[1,2]$ and total-variation accuracy
$1/10$, covering arbitrary measurable adaptive algorithms and unrestricted
real computation.
\url{https://arxiv.org/pdf/2304.02599v2}
\begin{theorem}\label{thm:lower}
There are absolute constants $c>0$, $d_0$ with $Q(d)\ge c\log d$ for every
$d\ge d_0$, even restricting to $V(x)=x^TAx/2$, $I\preceq A\preceq2I$. The
query cap is on every execution, and the algorithm may use arbitrary
private randomness independent of its input.
\end{theorem}
For these potentials the gradient reply is $Ax$, and the value
$x^TAx/2$ is determined by the query and its reply, so an exact
value-and-gradient query is one query to $x\mapsto Ax$. Each potential is
in $\mathcal V_d$, has minimizer $0$, target $N(0,A^{-1})$.
Three steps: deferred revelation of a random rotation represents any
output in a Gaussian block-Krylov span (including outputs outside the span
of actual queries and replies); a cosine statistic is uniformly small on
this enlarged random span; the same statistic has a large bias under the
target Gaussian, separating the joint experiments in total variation. None
uses the upper bound.
An arbitrary algorithm's output is represented; the auxiliary Gaussians are
independent of one another but need not be independent of the input matrix.
\begin{lemma}[Deferred rotation]\label{lem:deferred-haar}
Let $\Lambda$ be fixed real diagonal, $O$ Haar on the orthogonal group,
$A=O\Lambda O^T$. For any measurable randomized algorithm, private
randomness independent of $O$, making at most $q$ adaptive queries to
$x\mapsto Ax$ with $2q+1<d$ and output $X\in\R^d$, there is an extension
carrying independent standard Gaussians $g_1,\ldots,g_{q+1}\in\R^d$ with
\begin{equation}\label{eq:krylov-span}
O^TX\in\operatorname{span}\{\Lambda^jg_i:1\le i\le q+1,\ 0\le j\le q\}
\quad\text{almost surely}.
\end{equation}
The joint law including the randomness and $O$ is unchanged; the $g_i$
need not be independent of $O$ or $A$.
\end{lemma}
\begin{proof}
Include all randomness, including at output, in a seed independent of $O$;
fix it and pad to exactly $q$ query slots.
Haar property: for a fixed restriction $T:E\to W$ of an orthogonal map, the
remaining restriction is uniform among isometries $E^\perp\to W^\perp$; a
specified unit vector's partner is uniform on the other sphere, and after
revealing the pair the map between smaller complements is again uniform.
This follows from invariance under orthogonal maps fixing $E,W$; iterating
also permits each direction to be a measurable function of the revealed
data. Work with regular conditional laws on finite-dimensional data after
fixing the seed, representing subspaces and partial isometries measurably by
Gram--Schmidt in order of appearance, discarding zero residuals, completing
bases by scanning standard coordinates (Borel on each rank stratum).
Initially $E=W=\{0\}$. After each reply reveal $T=O|_E$, $E$ the span of
eigenframe coordinates of all queries and replies so far, $W=OE$; the
invariant is that the unrevealed map is Haar on the complements given the
full latent history. Let $x$ be the next query, fixed given the history;
write $x=x_W+x_\perp$ on $W\oplus W^\perp$. The preimage of $x_W$ is
known. If $x_\perp\ne0$, put $w=x_\perp/|x_\perp|$, $v=O^Tw\in
E^\perp$, conditionally uniform on the complementary sphere. Let $h$ be a
standard Gaussian in $E$ and $a$ the norm of a standard Gaussian in
$\R^{d-\dim E}$, drawn independently of the unrevealed map given the
history, and define $g_i=h+av$; the Gaussian polar decomposition makes
$g_i$ a standard Gaussian in $\R^d$ given the history, and $a>0$ a.s. gives
$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 $g_i$ reveals
no further part of $O$ (its parallel component and complementary radius are
independent auxiliary variables, its normalized complementary component is
$v$), preserving the invariant. If $x_\perp=0$, draw a fresh independent
standard Gaussian $g_i$ and make no extension.
Either way the preimage $u=O^Tx$ is now known in $E+\operatorname{span}\{g_i\}$,
reply preimage $\Lambda u$. Project onto the enlarged eigenframe space; if
the residual is nonzero, reveal the image of its unit direction (other form
of the same rule) and extend once more. The reply is exactly $Ax$ and the
invariant persists.
Each $g_i$ has conditional law $N(0,I_d)$ given earlier vectors and history,
hence they are independent. For polynomial degrees, let $E_i$ be the
eigenframe space after reply $i$; induction gives
$E_i\subseteq\operatorname{span}\{\Lambda^jg_a:1\le a\le i,\ 0\le j\le i\}$,
since the new preimage is in $E_{i-1}+\operatorname{span}\{g_i\}$ and its
reply applies $\Lambda$ once, old directions staying in the larger span.
Treat the output as one more query with no reply: repeat the preimage
construction using $g_{q+1}$, giving
$O^TX\in E_q+\operatorname{span}\{g_{q+1}\}$, i.e.
\eqref{eq:krylov-span}. At most two dimensions exposed per query and one
for the output, so the dimension inequality keeps every complement
available. Started from the genuine Haar input and never altered it; the
added variables are measurable kernels from the original variables and
independent draws, so the construction is jointly measurable for fixed-seed
algorithms, and integrating the independent seed handles randomized ones.
\end{proof}
Henceforth set
\begin{equation}\label{eq:lower-spectrum}
\theta_k=\frac{(k-1/2)\pi}{d},\quad
\lambda_k=\frac{3+\cos\theta_k}{2},\quad
\Lambda=\diag(\lambda_1,\ldots,\lambda_d),\ 1\le k\le d,
\end{equation}
all eigenvalues strictly between $1$ and $2$. For $q\ge0$ put
$n=2q+1$, $p_k=\cos(n\theta_k)$, $D_p=\diag(p_1,\ldots,p_d)$; the frequency
$n$ exceeds twice the polynomial degree in \eqref{eq:krylov-span}, making
$u^TD_pu$ small for all $u$ in that random span.
\begin{lemma}\label{lem:lower-span-statistic}
Suppose $4q+1<2d$, let $g_1,\ldots,g_{q+1}$ be independent standard
Gaussians, $\mathcal K_q=\operatorname{span}\{\Lambda^jg_i:1\le i\le q+1,
0\le j\le q\}$, $\tau_d=40(q+1)^2/\sqrt d$. If $\tau_d\le1/2$, then with
probability at least $0.99$,
\begin{equation}\label{eq:lower-uniform-statistic}
|u^TD_pu|\le2\tau_d|u|^2\qquad\text{for every }u\in\mathcal K_q.
\end{equation}
\end{lemma}
\begin{proof}
With $m=(q+1)^2$, form the $d\times m$ matrix $B$ with columns indexed by
$(i,j)$, $1\le i\le q+1$, $0\le j\le q$:
$B_{k,ij}=d^{-1/2}(g_i)_kc_j\cos(j\theta_k)$, $c_0=1$, $c_j=\sqrt2$
($j\ge1$). The $\cos(j\theta)$ are degree-$j$ polynomials in $\cos\theta$
with nonzero leading coefficient, and $\cos\theta_k=2\lambda_k-3$, so the
range of $B$ is exactly $\mathcal K_q$.
Discrete cosine orthogonality gives
\begin{equation}\label{eq:lower-gram-means}
\E B^TB=I_m,\qquad \E B^TD_pB=0 .
\end{equation}
For the half-grid,
\begin{equation}\label{eq:lower-cosine-sum}
\sum_{k=1}^d\cos(\ell\theta_k)
=\frac{\sin(\ell\pi)}{2\sin(\ell\pi/(2d))}=0,\quad
1\le\ell<2d,\ \ell\in\mathbb Z,
\end{equation}
the denominator nonzero in range. Expand products of two cosines for the
first identity, and of three with frequencies $n,j,j'$ for the second (all
nonzero since $n>j+j'$, of absolute value at most $4q+1<2d$); different
Gaussian blocks have zero cross means.
Each Gram entry is a sum of independent row terms with coefficient
absolute value at most $2/d$ and product-variance at most $2$, so entry
variance at most $8/d$ and
\[
\E\bigl[\|B^TB-I_m\|_{\HS}^2+\|B^TD_pB\|_{\HS}^2\bigr]\le\frac{16m^2}{d}.
\]
Markov gives
$\P(\|B^TB-I_m\|_{\HS}^2+\|B^TD_pB\|_{\HS}^2>\tau_d^2)\le
\frac{16m^2/d}{1600m^2/d}=0.01$. Since $\|\cdot\|_{\op}\le\|\cdot\|_{\HS}$,
outside this event both operator norms are at most $40m/\sqrt d=\tau_d$.
Then for every $a\in\R^m$, $|Ba|^2\ge(1-\tau_d)|a|^2\ge\tfrac12|a|^2$ and
$|(Ba)^TD_p(Ba)|\le\tau_d|a|^2$, giving
\eqref{eq:lower-uniform-statistic} on the whole range of $B$.
\end{proof}
The uniformity matters: the output's coefficients in the block span may
depend arbitrarily on the whole transcript.
An ideal target has a detectable bias at the same frequency. Write
$\varrho=3-\sqrt8\in(0,1)$; the geometric-series formula gives the
absolutely convergent identity
\begin{equation}\label{eq:lower-fourier}
\frac{2}{3+\cos\theta}=\frac1{\sqrt2}\Bigl(1+2\sum_{j\ge1}(-\varrho)^j
\cos(j\theta)\Bigr),
\end{equation}
the parenthesis summing to $(1-\varrho^2)/(1+2\varrho\cos\theta+\varrho^2)$
with $1+\varrho^2=6\varrho$, $1-\varrho^2=4\sqrt2\,\varrho$.
\begin{lemma}\label{lem:lower-target-statistic}
Let $U\sim N(0,\Lambda^{-1})$ with $\Lambda$ from
\eqref{eq:lower-spectrum} and $n=2q+1<d$. Then
\begin{align}
\mu_{d,q}:=\E[U^TD_pU]&=\frac d{\sqrt2}(-\varrho)^n+R_{d,q},\quad
|R_{d,q}|\le\frac{\sqrt2\,d\,\varrho^{2d-n}}{1-\varrho},
\label{eq:lower-target-mean}\\
\Var(U^TD_pU)&\le2d,\quad \E|U|^2\le d,\quad \Var(|U|^2)\le2d.
\label{eq:lower-target-variance}
\end{align}
\end{lemma}
\begin{proof}
Independence gives $\mu_{d,q}=\sum_kp_k/\lambda_k$. Substitute
\eqref{eq:lower-fourier}; the constant term vanishes by
\eqref{eq:lower-cosine-sum}. For $1\le j<2d-n$,
$\cos(n\theta_k)\cos(j\theta_k)=\tfrac12\cos((n-j)\theta_k)
+\tfrac12\cos((n+j)\theta_k)$; if $j\ne n$ both frequencies are nonzero and
below $2d$ so their sums vanish, while if $j=n$ the constant term
contributes $d/2$ and frequency $2n<2d$ contributes zero. Thus only $j=n$
contributes there; each remaining cosine product has absolute sum at most
$d$, so the geometric tail gives \eqref{eq:lower-target-mean}. A centered
normal of variance $v$ has $\Var(U_k^2)=2v^2$; with $|p_k|\le1$,
$\lambda_k\ge1$, summing coordinate variances gives
\eqref{eq:lower-target-variance}.
\end{proof}
\begin{proof}[Proof of Theorem~\ref{thm:lower}]
Take $A=O\Lambda O^T$, $O$ Haar, $\Lambda$ as in
\eqref{eq:lower-spectrum}. The reduction identifies an exact
value-and-gradient algorithm on $V_A(x)=x^TAx/2$ with an algorithm querying
$x\mapsto Ax$, same cap and output law.
Fix any algorithm with at most $q$ queries on every execution,
$0\le q\le c_0\log d$, $c_0=1/(16\log(1/\varrho))$. All dimension
restrictions in Lemmas~\ref{lem:deferred-haar},
\ref{lem:lower-span-statistic}, and~\ref{lem:lower-target-statistic} hold
for 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, and define
\begin{equation}\label{eq:lower-separating-event}
\mathcal E_{d,q}=\Bigl\{(O,x): |O^Tx|^2\le2d,\ \
|(O^Tx)^TD_p(O^Tx)|>160(q+1)^2\sqrt d\Bigr\}.
\end{equation}
This event is measurable and determined by $(A,x)$: with $T_n$ the
polynomial $T_n(\cos\theta)=\cos(n\theta)$, $OD_pO^T=T_n(2A-3I)$.
Apply Lemma~\ref{lem:deferred-haar} to the output, then
Lemma~\ref{lem:lower-span-statistic} to its Gaussians. On the latter's
event the norm condition forces the quadratic form to be at most
$2\tau_d(2d)=160(q+1)^2\sqrt d$, so
\begin{equation}\label{eq:lower-algorithm-event}
\P_{\mathrm{alg}}(\mathcal E_{d,q})\le0.01 .
\end{equation}
For $\P_{\mathrm{target}}$, sample $O$ from the same prior then $X$ from
$N(0,O\Lambda^{-1}O^T)$; conditional on $O$, $U=O^TX\sim N(0,\Lambda^{-1})$.
The choice of $c_0$ gives $\varrho^{2q+1}\ge\varrho d^{-1/8}$ and
\[
\frac{d\varrho^{2q+1}}{(q+1)^2\sqrt d}
\ge\frac{\varrho d^{3/8}}{(1+c_0\log d)^2}\longrightarrow\infty .
\]
The remainder satisfies
$|R_{d,q}|/(d\varrho^n/\sqrt2)\le2\varrho^{2d-2n}/(1-\varrho)\to0$
uniformly since $n\le1+2c_0\log d$, hence is negligible; with
$T_{d,q}=160(q+1)^2\sqrt d$, $|\mu_{d,q}|>T_{d,q}$ holds for large $d$,
and by Chebyshev
\[
\P_{\mathrm{target}}(\mathcal E_{d,q}^c)
\le\frac2d+\frac{2d}{(|\mu_{d,q}|-T_{d,q})^2}=o(1)
\]
uniformly. Thus the two joint experiments' probabilities of
$\mathcal E_{d,q}$ differ by at least $0.99-o(1)$.
To compare with the guarantee on each fixed input, let $K_O$ be the
output law on $A=O\Lambda O^T$ and $G_O=N(0,A^{-1})$; for a joint event
$\mathcal E$ write $\mathcal E_O=\{x:(O,x)\in\mathcal E\}$. Both laws use
the same Haar prior, so
\begin{align*}
\bigl|\P_{\mathrm{alg}}(\mathcal E)-\P_{\mathrm{target}}(\mathcal E)\bigr|
&=\Bigl|\int[K_O(\mathcal E_O)-G_O(\mathcal E_O)]\,d\mathrm{Haar}(O)\Bigr|\\
&\le\int\|K_O-G_O\|_{\TV}\,d\mathrm{Haar}(O).
\end{align*}
The pointwise guarantee would bound this by $1/10$, and the event testing
the joint laws need not be computable by the algorithm. This is impossible
for \eqref{eq:lower-separating-event}. So no algorithm with
$q\le c_0\log d$ meets the accuracy for all inputs; decreasing the constant
proves the theorem.
\end{proof}
\end{document}