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

First suppose \(H\) is smooth and has bounded derivatives of orders at least two. For a smooth compactly supported \(h\), let \(P_th(x)=\mathbb E[h(X_t)\mid X_0=x]\). Differentiating the synchronous flow in the initial state gives bounded first and second derivatives on bounded time intervals. Itô’s formula and the semigroup identity yield \[\partial_tP_th=P_t\mathcal Ah=\mathcal AP_th, \qquad \mathcal A=\Delta-\nabla H\cdot\nabla.\] The equality with \(\mathcal AP_th\) also follows by differentiating \(P_{t+s}h=P_s(P_th)\) at \(s=0\). Its at-most-linear growth is integrable against \(\mu\). Integration by parts gives \(\int\mathcal AP_th\,d\mu=0\), hence invariance. Coupling to an invariant starting point and using contraction shows \(P_th\to\mu h\) in \(L^2(\mu)\). Therefore \[\operatorname{Var}_\mu h=2\int_0^\infty \int|\nabla P_th|^2\,d\mu\,dt.\] The flow derivative has operator norm at most \(e^{-mt}\). Jensen’s inequality and invariance thus bound the inner integral by \(e^{-2mt}\int|\nabla h|^2d\mu\), proving the Poincaré inequality.

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

For (6), start the two diffusions with their respective invariant marginals, with any initial coupling, and use common Brownian motion. Writing \(a(t)=(\mathbb E|X_t-Y_t|^2)^{1/2}\) and \(\beta=(\mathbb E_\eta|\nabla H-\nabla\bar H|^2)^{1/2}\), strong monotonicity and Cauchy–Schwarz give, in the usual upper-derivative sense, \[a'(t)\le-ma(t)+\beta.\] The marginals remain fixed, so \(W_2(\mu,\eta)\le a(t)\). Letting \(t\) tend to infinity proves the claim. ◻

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

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

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

Define the symmetric Appell tensor polynomials by \[P_j(s)= \left.\nabla_\theta^j \frac{\exp(\theta\cdot s)}{\mathbb E\exp(\theta\cdot S)} \right|_{\theta=0}.\] For a coefficient tensor \(T\) of order \(j\), let \(Q_T(s)=T:P_j(s)\), where the colon denotes contraction on all slots. Differentiating the generating function in \(s\) multiplies it by copies of \(\theta\). Consequently \[\mathbb E\nabla^h Q_T(S)=0\quad(0\le h<j),\qquad \nabla^j Q_T=j!\operatorname{Sym}(T).\] Here \(\operatorname{Sym}(T)\) is the average over all permutations of the slots; hence \(|\operatorname{Sym}(T)|_F\le|T|_F\). Apply Poincaré to each component of \(\nabla^h Q_T\) and sum. Iterating from \(h=0\) to \(h=j-1\) yields \[ \mathbb E|T:P_j(S)|^2\le \beta^j(j!)^2|T|_F^2. \tag{10}\] Only first derivatives in the displayed Poincaré inequality are used at each application.

In particular, if \(T\) and \(V\) have orders \(j,k\ge1\), then \[ \left|\mathbb E\bigl[(T:P_j(S))(V:P_k(S))\bigr]\right| \le\beta^{(j+k)/2}j!k!\,|T|_F|V|_F. \tag{11}\] This is the direct covariance estimate. It concerns arbitrary coefficient tensors, rather than products of one-dimensional test directions, which is why it controls a matrix split.

To pass from covariance to cumulants while preserving a specified split, consider the two-variable generating function \[\mathcal G(\theta,\omega)= \frac{\mathbb E\exp((\theta+\omega)\cdot S)} {\mathbb E\exp(\theta\cdot S)\,\mathbb E\exp(\omega\cdot S)}.\] Its derivative of orders \(j\ge1\) in \(\theta\) and \(k\ge1\) in \(\omega\) at \((0,0)\) is the cross-moment tensor \(\mathbb 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 (11).

Both \(\mathcal G(\theta,0)\) and \(\mathcal G(0,\omega)\) are identically one. Thus every nonzero term of \(\mathcal G-1\) contains variables from both groups. In the formal identity \[\log\mathcal G =\sum_{m\ge1}\frac{(-1)^{m+1}}m(\mathcal G-1)^m,\] only \(m\le n\) can contribute at total order \(n\). Every differentiated product is a tensor product of the preceding cross-moment operators, up to permutations within the two sides. Its norm is therefore at most the product of their norms. To bound the total coefficient, distribute the \(n\) labeled derivative slots among at most \(n\) ordered nonempty factors. There are at most \(n^{n+1}\) such assignments; the product of the factorials within any assignment is at most \(n!\). Enlarging an absolute constant bounds their total by \(C^n n^{Cn}\).

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

Proof. The conditional potential in (9) has Hessian between \((a^{-1}-l)I\) and \((a^{-1}+l)I\). By Lemma 3, its Poincaré constant is at most \(Ca\). Under this conditional law the map \[z\longmapsto S(z)=(z,f(z)/L)\] has derivative norm at most \(\sqrt2\), since \(f\) is \(L\)-Lipschitz. Pulling back an ambient test function shows that the possibly singular law of \(S\) satisfies Poincaré with constant \(Ca\). It has exponential moments of every linear functional: \(f\) has at most linear growth and the conditional law has Gaussian tails. Lemma 5 therefore applies. This use of repeated Poincaré does not differentiate \(f\) repeatedly. Each application pulls back a new ambient polynomial and uses only \(Df\).

In the conditional density, the linear tilt parameter for \(z\) is \(\rho y/a\). Repeated differentiation of an expectation with respect to a linear tilt gives joint cumulants. In coordinates this says \[ (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). \tag{15}\] It can also be checked directly by differentiating the logarithm of the conditional moment-generating function with one auxiliary parameter coupled to \(f_i\). Restrict the cumulant of \(S\) to one \(f/L\) coordinate and \(q\) ordinary coordinates. Lemma 5 gives a factor \(L a^{(q+1)/2}\), and the prefactor \(a^{-q}\) gives (12).

All the conditional integrals used here are smooth in \((\rho,y)\) for \(\rho<1\). This follows by differentiating their densities: on compact parameter sets, every derivative inserts a polynomial with an integrable Gaussian majorant. Formula (15) gives the smooth extension at zero without division by \(\rho\). For later integrations by parts we record a slightly stronger regularity fact. On each compact interval in \(\rho<1\), conditional moments of any fixed order grow at most polynomially in \(y\). Indeed, the conditional mode satisfies \[z_0=\rho y-ar g(x+rz_0),\] so \(|z_0|\le C(|y|+r|g(x)|)\). The upper and lower conditional Hessian bounds compare the density and its normalizing integral around this mode to Gaussian integrals. This proves the claimed moment growth. The same argument applies to time derivatives of the conditional integrals. Dimension-dependent constants are harmless for this regularity justification; the quantitative estimates above remain dimension-uniform.

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

Proof. For the random interpolation used to define \(\xi_\rho\), \[\frac{d}{d\rho}(\rho Z+\sqrt a\,N) =\frac{Z-\rho Y}{a}.\] The conditional mean of this derivative is \(-rM_\rho(Y)\) by the conditional integration-by-parts identity above. Thus, for compactly supported smooth \(\phi\), \[\frac{d}{d\rho}\int\phi\,d\xi_\rho =-r\int M_\rho\cdot\nabla\phi\,d\xi_\rho.\] On the other hand, (12) with \(q=1\) gives \(r\|DM_\rho\|_{\mathrm{op}}\le Cl\). The velocity has at most linear growth, so its smooth ordinary differential equation has a global-in-space flow and inverse on every compact interval in \(\rho<1\). To identify its pushforward, fix a smooth compactly supported terminal test and transport it backward along this flow. The transported test has zero material derivative. Its supports remain in a bounded set on the interval, by the linear-growth bound. The last display then shows its expectation under \(\xi_\rho\) is constant, which proves the pushforward assertion. ◻

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

For the higher identity, coordinate differentiation gives \[\mathcal D\partial_I M_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 (17) and multiply by \(\rho^{-q}\). When all derivatives pass inside the adjoint divergence, the result is \(\nabla^*U_{q+1}\). In the other terms they hit the score \(y+r\rho M_\rho\). A derivative hitting its linear part contributes \(qU_q/\rho\), canceled by differentiating the prefactor \(\rho^{-q}\). The nonlinear score terms and the material/spatial commutator terms are identical after scaling: each is \[r\rho\sum_{\varnothing\ne S\subseteq I} (U_{|S|})_{k,S}(U_{q-|S|+1})_{i,k,I\setminus S}.\] This proves (18). Each contraction retains at least one free slot in each tensor.

Finally differentiate the adjoint identity against a transported scalar test, componentwise in the other slots. The terms involving \(\mathcal D\phi\) cancel. The gradient commutator is \[\mathcal D\partial_k\phi-\partial_k\mathcal D\phi =r(\partial_kM_j)\partial_j\phi.\] Moving it to the other side by the adjoint relation gives (19). All coefficients in the final formulas are nonsingular at zero, so the extensions follow from Proposition 6. ◻

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

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

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

It remains to consider \(E=\nabla^*V\), where \(V\) has weight \(k_V\) and \(E\) has at least two slots. Commute the \(j\) derivatives past the divergence. Apart from \[T=\nabla^*(\nabla^jV),\] the terms are single contractions of \(\nabla^{j-h}V\) with an order-\(h\) derivative of the score, \(1\le h\le j\). Both factors retain slots. The score derivative is bounded pointwise by a constant of the form \(C^h h^{Ch}(\sqrt a)^{-(h-1)}\), by (13)–(14). Lemma 9, the inductive hypothesis, and Leibniz multiplicities therefore give the required bound for all these terms, even with exponent \(k_V+j-1\) in place of the larger allowed exponent \(k_V+j+1\). It remains to bound \(T\).

Proof. First take a compactly supported smooth vector field, so that the output is scalar. Integrating by parts twice gives \[\int(\nabla^*V)^2\,d\mu =\int\sum_{i,j} \left[(\partial_iV_j)(\partial_jV_i) +V_i(\partial_{ij}H)V_j\right]d\mu.\] The first term is bounded in absolute value by \(|\nabla V|_F^2\), and the second by \(\|\nabla^2H\|_{\mathrm{op}}|V|^2\). For tensor output, apply the scalar identity to each fixed collection of remaining indices and sum. Taking a square root proves (23) for compactly supported fields. For general \(V\), multiply by smooth cutoffs \(\chi_R\) equal to one on the ball of radius \(R\), with \(|\nabla\chi_R|\le C/R\). Then \(\chi_RV\to V\) and \(\nabla(\chi_RV)\to\nabla V\) in \(L^2(\mu)\), while their adjoint divergences converge pointwise. Fatou’s lemma proves the stated estimate. All fields to which we apply it have polynomially growing derivatives and hence satisfy the required integrability. ◻

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

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

Proof. Differentiate the transport equation with respect to its starting position. Gronwall’s inequality gives a bounded first derivative and an \(O(l)\) bound for its difference from \(I_d\). For order \(k>1\), the chain rule has terms consisting of an order-\(j\) derivative of the velocity and derivatives of the state of orders \(i_1,\ldots,i_j\ge1\), where \(i_1+\cdots+i_j=k\). By (12), the outer tensor has every split norm at most \(C_j lR^{-(j-1)}\). Attach the inner tensors one at a time. Each attachment is a single contraction retaining the outer output and the inner derivative slots, so Lemma 9 applies without a dimension factor. Bound a first state derivative by a constant and use the induction hypothesis for all higher state derivatives. Since \(l\le1\), the total power of \(R^{-1}\) is at most \[(j-1)+\sum_{\nu=1}^j(i_\nu-1)=k-1,\] and each forcing term retains the outer factor \(l\). The term involving the unknown order-\(k\) state derivative is linear, with coefficient bounded by \(Cl\), and this derivative starts at zero. Gronwall proves the assertion. The same argument with absolute interval lengths applies to transport in the reverse direction.

For the second assertion, write \(C_\rho=rDM_\rho(Y_\rho)=r\rho U_1(\rho,Y_\rho)\), so that \(B_\rho'=B_\rho C_\rho\). The product rule expresses \(B_\rho^{(m)}\) as a finite sum of ordered products of \(B_\rho\) and \(k\ge1\) factors \(C_\rho^{(j_1)},\ldots,C_\rho^{(j_k)}\), with \(j_1+\cdots+j_k=m-k\). The zero-order jet bound and (25) imply, for every fixed moment order, \[\big\|\|C_\rho^{(j)}\|_{\mathrm{op}}\big\|_{L^p} \le \operatorname{polylog}(d)\,l(1-\rho^2)^{-j}.\] Here the path estimates apply because inverse transport sends \(y\sim\xi_*\) to \(Y_\rho(y)\sim\xi_\rho\). Apply Hölder with sufficiently high fixed moments to each product. Using \(l^k\le l\) and \((1-\rho^2)^{-(m-k)}\le(1-\rho^2)^{-m}\) gives the claimed conservative exponent. Moments below one, if desired, follow from the \(L^1\) bound. ◻

Proof. Fix an output coordinate \(i\) and a smooth compactly supported scalar function \(\phi\). The transported test \(\phi\circ T_\rho\) has zero material derivative: it is constant along each probability-transport path. Differentiate its product with \(M_{\rho,i}\) under the transported laws. Lemma 8 gives \(\mathcal DM_\rho=\nabla_{\xi_\rho}^*U_1\), so integration by parts and the chain rule yield \[\frac{d}{d\rho}\mathbb E_{\xi_\rho}[M_{\rho,i}\,\phi\circ T_\rho] =\mathbb E_{\xi_\rho}\sum_{j,k}(U_1)_{ik}(DT_\rho)_{jk} (\partial_j\phi)\circ T_\rho.\] Smooth transport on the compact time interval and Gaussian tails justify differentiation; alternatively the test can first be cut off and the resulting identities passed to the limit. Push the expectation forward to \(y\sim\xi_*\). Since \(U_1\) is symmetric, its coefficient of \(\partial_j\phi(y)\) is \[\sum_k (U_1)_{ik}(\rho,Y_\rho(y))(B_\rho(y))_{jk} =(B_\rho(y)U_1(\rho,Y_\rho(y)))_{ji}.\] The terminal expectation is \(\mathbb E_{\xi_*}[M_{*,i}\phi]\). At zero, \(M_0=u\) and \((T_0)_\#\xi_0=\xi_*\), so the initial expectation is \(u_i\mathbb E_{\xi_*}\phi\). Integrating the displayed derivative in \(\rho\) proves \[\mathbb E_{\xi_*}[(M_{*,i}-u_i)\phi] =\mathbb E_{\xi_*}\sum_jK_{ij}\partial_j\phi.\] This proves the adjoint-divergence identity weakly, and smoothness makes it a pointwise identity.

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

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

Write the product density as \[p(y,g)=(2\pi)^{-d/2}\exp(-H_*(y)-|g|^2/2).\] Differentiating in the position and momentum variables separately, with the row convention of the kernel lemma, gives \[\operatorname{div}_y(pv_Y) =\frac{p}{s}\sum_i g_i(\nabla_{\xi_*}^*K)_i =\frac{p}{s}\,g\cdot(M_*-u),\qquad \operatorname{div}_g(pv_G) =-\frac{p}{s}\,g\cdot(M_*-u).\] Thus \(\operatorname{div}(pv)=0\). The Jacobian formula for the complete flow now yields, pointwise, \[\frac{d}{dt}\{p(\Phi_t(z))\det D\Phi_t(z)\} =\operatorname{div}(pv)(\Phi_t(z))\det D\Phi_t(z)=0.\] The determinant starts at one and stays positive. Change of variables therefore proves that \(\Phi_t\) preserves the probability density \(p\).

For later differentiation it is useful to record the same cancellation as an operator identity. Relative to the product law, the Lie derivative is \[ \mathcal L=\nabla^*A\nabla,\qquad A=\begin{pmatrix}0&K^T/s\\-K/s&0\end{pmatrix}. \tag{30}\] Indeed the second-order part vanishes by skew symmetry. The first-order coefficients are \(-K^Tg/s\) in position and \(\nabla_{\xi_*}^*K/s=(M_*-u)/s\) in momentum, exactly as in (28).

Finally, integrating the momentum equation over \([0,1]\) gives the pathwise identity \[sG_0+\int_0^1M_*(Y_t)\,dt=u+sG_1.\] Stationarity makes \(G_1\) a standard Gaussian, proving (29). ◻

Proof. Put \(\phi_k=\mathcal L^k\phi_0\). We prove, for \(j\ge0\) and \(k\ge1\), the stronger estimate \[ \|\nabla^j\phi_k\|_{F,2} \le C_{j,k}a_0\sqrt d\,\alpha^kR^{-(j+2k)}. \tag{32}\] For \(k=0\) the same bound is assumed only for \(j\ge1\), which suffices: the iteration differentiates \(\phi_{k-1}\) before multiplying by \(A\). All \(L^2\) norms below are for the invariant density.

Set \(V=A\nabla\phi_{k-1}\). Leibniz’s rule and the induction hypothesis give, for each \(h\ge0\), \[\|\nabla^h V\|_{F,2} \le C_{h,k}a_0\sqrt d\,\alpha^k R^{-(h+1+2(k-1))}.\] To see why no dimension factor appears, in each coefficient tensor \(\nabla^vA\) isolate the slot contracted against the gradient of \(\phi_{k-1}\). Its one-slot operator norm multiplies the Frobenius norm of the other factor. Summing the finitely many Leibniz terms preserves the displayed estimate.

Now commute \(j\) derivatives past the adjoint divergence in \(\phi_k=\nabla^*V\). The main term is \(\nabla^*(\nabla^jV)\). The Frobenius estimate (23) bounds its \(L^2\) norm by a constant times \(\|\nabla^jV\|_{F,2}+\|\nabla^{j+1}V\|_{F,2}\). These two terms have powers \(R^{-(j+2k-1)}\) and \(R^{-(j+2k)}\), respectively. Each remaining commutator contracts an order-\(v\) derivative of the score, \(1\le v\le j\), with \(\nabla^{j-v}V\). The coefficient’s one-slot norm again suffices; the total exponent of \(R^{-1}\) is \[(v-1)+(j-v+1+2(k-1))=j+2k-2.\] All these terms are bounded by (32) because \(0<R\le1\). This closes the induction, including \(j=0\) for \(k\ge1\). Gaussian tails and polynomial growth justify the integrations by parts through cutoffs.

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

Proof. Write \(Y_\rho\) for the exact path from the given random input and \(Q_\rho=M_\rho(Y_\rho)\). The marginal identity in Proposition 7 and (24) 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 \(Q\) of degree \(m-1\) at \(a_c\). Minkowski’s inequality applied to the integral remainder gives, at each \(v\in I\), \[ \|Q_v-P(v)\|_2 \le \frac1{(m-1)!}\int_{a_c}^v (v-z)^{m-1}\|Q_z^{(m)}\|_2\,dz \le C_mB_dL\sqrt d\,h_0^m. \tag{38}\] The last step uses \(\ell_c\le Ch_0(1-z)\) on that cell and \(1-z^2\ge1-z\). The same estimate holds at the interpolation nodes. Because interpolation reproduces \(P\) and its basis has bounded absolute sum, \(\|Q_v-\mathcal I_cQ(v)\|_2\) obeys the same bound. Integrate this pointwise estimate over the cells, whose total length is at most one. This proves the required quadrature error, for partial integrals and for either orientation. No supremum of the random path has been taken.

Let \(e_k=\max_i\|y_i^{[k]}-Y_{\rho_i}\|_2\). The global bound \(\operatorname{Lip}(M_\rho)\le CL\), together with the absolute row sums, gives \[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\). Integrating the velocity shows \(e_0\le CrD\). Since \(C_ml<1/2\) for large \(d\), summing the geometric recurrence proves (37). The norm here is the maximum of the nodewise RMS errors; it is not the RMS of their maximum.

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

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

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

For the first cell, subtract the Taylor polynomial at zero of degree \(m-1\). Its derivative at zero agrees with the exact zero derivative. The second derivative of the remainder, and that of the Hermite interpolant of the remainder, are bounded in RMS by \(C_m\operatorname{polylog}(d)l\sqrt d\,\ell_c^{m-2}\). Both first derivatives vanish at zero. If \(E\) is the difference between the exact function and its Hermite interpolant on this cell, then \[\|E'(\rho)\|_2 \le\int_0^\rho\|E''(v)\|_2\,dv \le C_mB_dl\sqrt d\,\rho\ell_c^{m-2}.\] Consequently the integral of \(E'(\rho)/\rho\) over the first cell has RMS norm at most \(C_mB_dl\sqrt d\,\ell_c^{m-1}\). This is within the stated, more conservative \(h_0^{m-2}\) bound.

On an ordinary cell the absolute coefficient sum of the interpolant derivative is at most \(C_m/\ell_c\). Integrating and using \(\ell_c\ge c'h_0(1-\rho)\) gives a total bounded by \[C_mh_0^{-1}\int_{\rho_{\mathrm{first}}}^{b_*} \frac{d\rho}{\rho(1-\rho)}.\] This gives the logarithms in (44). On the first cell, scale to the unit interval. The polynomial derivative vanishes at zero and is bounded by distance to zero times the second derivative, so the integral of its absolute value divided by that distance is bounded. Scaling back gives a weight sum \(C_m/\ell_c\le C_mh_0^{-1}\). This proves the claim. ◻

Proof. Combine Lemmas 19 and 20, and divide by \(rs\). The value errors are amplified only by \(C_m\psi^{-1}\) and the weight sum (44). The exact stencil truncation has factor \(\psi^{m-1}R^{-2m}=\delta^{2tm-4t}\) and the same probability-quadrature weight sum. Crucially, \[\frac{l}{r\sigma}=\frac L\sigma,\qquad \frac{l^{N_p}rD}{r\sigma} =l^{N_p-1}r\frac L\sigma D.\] Thus neither division introduces a negative power of the potentially small scale \(r\). The resulting estimate is (47).

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

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

Proof. First compare the exact flow to polynomial integration of its exact velocity. Lemma 17, applied to the state with one additional derivative, bounds this quadrature error by \[C_m\sqrt d\,\bigl((L/s)R^{-2}\bigr)^{m+1}.\] The approximation error of \(\widehat v\) at every exact trajectory node is bounded by Proposition 21, by stationarity. Previous state errors propagate with its global Lipschitz constant and the absolute row sums. The true centering velocity need not be globally Lipschitz in the joint variables; it is the approximate velocity’s bound that is used in this recurrence.

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

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

Proof. For large \(d\), probability and harmonic Picard contractions \(C_ml\le\operatorname{polylog}(d)\delta^b\) can be bounded by \(\delta^{1/3}\). The only additional inverse powers in kernel evaluation are \(\delta^{-4t-w}\), besides the divisions already canceled in Proposition 21. The centering iteration also contracts by \(\delta^{1/3}\). The probability iteration remainder is at most \(C\delta^{N/3}D\); after the kernel weights the iteration remainder is at most \(C\delta^{N/3-4t-w}D\). Both exponents exceed \(J\) with a fixed positive margin. For the truncation terms it suffices to check the exponents \[wm,\quad w(m-2),\quad w(m-1)-4t,\quad 2tm-4t-w,\quad (m+1)/3.\] The first four come respectively from probability quadrature, the zero-endpoint rule, its amplified probability-value error, and the harmonic and angle errors. For the last one use \((L/s)R^{-2}\le\delta^{1/3}\) for large \(d\) in both outer quadratures, dividing the mean-output error by \(\sigma\) and using \(L/\sigma\le\delta^{1/3}\). Every displayed exponent exceeds \(J\) by (50). The outer Picard remainder contracts by the same \(\delta^{1/3}\) factor. Thus fixed constants and polylogarithms, including those depending on \(m\), fit within the strict exponent margins as \(d\to\infty\).

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

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

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

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

Centers and their coefficients. Consider a mean evaluation inside any of these formulas. Its physical center is \(x+r\rho v\) when its normalized position is \(v\), by (9). Substitute its direct starting states: a backward probability state starts at the current centering position; a harmonic state starts at that backward position and the current momentum; and a forward probability state starts at the harmonic position. Stop each substitution at the mean terms just described. Only the free terms pass through this chain. Thus their anchor coefficients are exactly \((\rho,0)\) or \((\rho\cos z,\rho\sin z)\). For sampling there is only the first probability computation, so the coefficient is \((\rho,0)\). Multiplication by \(r\rho\), with \(\rho\le1\), bounds the center coefficient sum in the mean case by \[B_d(r/\sigma)\delta^{-4t-w}+C_mr^2.\] The term \(C_mr^2\) accounts for the backward, harmonic, and forward integral corrections. In the sampling case the bound is simply \(C_mr^2\). Enlarging these positive upper bounds to (54) is convenient for the implementation. The mean output is the single quadrature in (51), so its coefficient sum is at most \(C_m\) and its means all have probability time \(b_*\). For sampling, dividing a probability correction by \(\rho_e\ge2/\sqrt5\) gives output coefficient sum at most \(C_mr\). These observations also give the stated final direct terms.

The scale at probability time \(\rho\) is \(r\sqrt{1-\rho^2}\). In the mean computation \(\rho\le b_*\), giving \(\tau\ge R\), including equality for the output means. In the sampling computation, \[\sqrt{1-\rho_e^2}=\frac\eta{\sqrt{4+\eta^2}} \ge\frac R{\sqrt5}.\] This proves every assertion about the anchors and scales.

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

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

Ideal-law errors. Couple the numerical mean formula and the exact centering path integral by the same ideal anchors. Corollary 23 bounds their RMS difference by \(B_d\sigma\delta^JD(x)\), while Proposition 16 identifies the exact output as \(N(u,\sigma^2I_d/4)\). This coupling proves (55). For sampling, the exact probability transport from a standard Gaussian has law \(\xi_{\rho_e}\). After division by \(\rho_e\) this is the law of \[Z+\frac{\sqrt{1-\rho_e^2}}{\rho_e}G =Z+\frac\eta2G, \qquad Z\sim\nu_{x,r},\quad G\sim N(0,I_d)\text{ independent}.\] The same-anchor RMS estimate and \(\rho_e\ge2/\sqrt5\) prove (56). ◻

Proof. The eigenvalues of \(P\) lie in \([\sqrt3/2,1]\). The translation identity rewrites the first output as \[\mathsf Q\left(x;PG-\frac{v\otimes Z}{2D_s}\right).\] We suppress the tensor product with \(I_d\) on slot-space matrices. The seed array in this expression is centered Gaussian and has slot covariance \[PP^T+\frac{vv^T}{4D_s^2}=I_m.\] Its slots are therefore independent standard Gaussian vectors. They may remain correlated with \(Z\); that is harmless because the rewritten call is at the deterministic center \(x\). ◻

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

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

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

Choose \(w>0\) with \[C_3w<t/2,\qquad w<t,\qquad c_1w<t,\] and only then choose the fixed interpolation orders needed in Theorem 24. For sufficiently large \(d\), the last shift bound is smaller than \(D_s=\delta^{-t}\). This is a valid bottom-up induction: it proves each child’s bound before using its inverse transformation in the parent. The count bound is also at most \(\operatorname{polylog}(d)\delta^{-t}\).

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

Proof. For the standard Gaussian seed \(Z\), the primitive moment bounds give \[\|g(x+rZ)-u\|_{L^2} \le CL\bigl(\sqrt d+r|g(x)|\bigr)\le CLD(x).\] Since \(L/\sigma\le\operatorname{polylog}(d)\delta^b\) and \(p\le b\), this is within the mean error budget; couple the added Gaussian identically.

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

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

For a nonfinal term, let \(r\tau_{\rm prev}\) be the scale of the node whose center is currently being supplied. Its expression is evaluated with budget \(r\tau_{\rm prev}/(2D_s)\), irrespective of the tier at which it occurs. Consequently \[\sigma_i=\frac{r\tau_{\rm prev}}{4D_sA}, \qquad \frac{L_i}{\sigma_i} \le \frac{CD_s\lambda A}{R}.\] The coefficient estimates give \(\lambda A\le\operatorname{polylog}(d)\delta^{b-4t-c_1w}\) for an \(\mathsf M\) parent. For an \(\mathsf S\) parent, \(\lambda A=C_1l\le\operatorname{polylog}(d)\delta^b\), since \(\eta\le1\). Using \(D_s/R=\delta^{-2t}\) and \(c_1w<t\) costs less than the allowed \(10t\) decrease in \(b\). Also \(l_i\le l\). For the position-anchor \(\mathsf S\) call of \(\mathsf M\), \(l/\eta\le l/R\), which has the same allowance. This proves the size conditions by the bounded routine height. All accumulated fixed and polylogarithmic factors are still one polylogarithmic factor.

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

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

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

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

For a nonterminal \(\mathsf M\), first couple its approximate \(\mathsf S\) position anchor to an ideal anchor, independently of all remaining seed blocks. The child labels are \((b-10t,p-1/4)\). By induction and Lemma 30, replacing this anchor costs at most \[\operatorname{polylog}(d)L\delta^{p-1/4}D(x).\] After division by the mean error unit \(\sigma\delta^pD(x)\), the cost is at most \(\operatorname{polylog}(d)\delta^{b-1/4}\), hence is bounded for large \(d\). We may now assume that the position anchor has law \(b_*Z+RN\), where \(Z\sim\nu_{x,r}\) and \(N\) is standard Gaussian, and that the momentum anchor is an independent standard Gaussian. These are exactly the input laws in Theorem 24.

Let \(x_i\) be the centers of the deterministic exact-mean circuit under these anchors, and put \(D_i=\sqrt d+|g(x_i)|\). Their RMS values over anchors obey \[ \max_i\|D_i\|_{L^2}\le\operatorname{polylog}(d)D(x). \tag{66}\] Indeed, both position and momentum anchors have RMS at most \(CD(x)\). Use \(|g(x_i)|\le|g(x)|+\lambda|x_i-x|\) and the circuit’s direct terms. An exact primitive mean at a preceding center has magnitude at most \(CD_j\). Their contribution to this maximum-of-RMS recurrence is at most \(C\lambda A\) times its previous maximum. This is small by (65); the direct terms cost at most \(CD(x)\). Induction through the node ordering proves (66) without a maximum inside the expectation.

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

Different term occurrences have disjoint blocks, so their outputs are independent conditional on the anchors. Product couplings compare them to independent noised exact means. Together with the independent residual in (63), the ideal sum has exactly covariance \(n^2I_d\). We conclude that \[ \epsilon_E\le\sum_i|a_i| \left[\operatorname{polylog}(d)\lambda\epsilon_{E_i} +\operatorname{polylog}(d)\sigma_i\delta^{p_i}D_i\right]. \tag{67}\] Notice the order of this argument: inaccurate supplied centers use only pathwise stability, and standalone accuracy is invoked only after Gaussian absorption. No assertion is made that the actual algorithm’s intermediate states have the ideal input laws.

Take RMS in the anchors and expand (67) through the finite expression tree. At its root, direct-child labels have \(p_i=p\) and \[\sum_i|a_i|\sigma_i \le A_f\frac{n}{2A_f}=n/2.\] Thus their intrinsic errors fit the requested budget. At tier \(j\ge2\), absolute coefficient sums and Lipschitz factors contribute at most \[A_f\bigl(\operatorname{polylog}(d)\lambda A\bigr)^{j-1}.\] All such children have the single accuracy-label loss \(p_i=p-1/4\), not one loss per expression tier. Their noise units are at most \(C\sigma\) in an \(\mathsf M\) parent and \(C/r\) in an \(\mathsf S\) parent by Lemma 29; these units combine with \(A_f=C_1\) or \(C_1r\), respectively. Already the first factor \(\operatorname{polylog}(d)\lambda A\le\operatorname{polylog}(d)\delta^{1/3}\) pays for the loss \(\delta^{-1/4}\), leaving a positive power \(\delta^{1/12}\). The bounded expression depth allows all further tiers and polylogarithmic factors to be summed. Leaf affine expressions have zero error. This proves the required final evaluation error.

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

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

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

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

Proof. Lemma 3, applied at the mode zero, gives \(\mathbb E_{\pi_V}|X|^2\le d\). Hence the root mean square under \(\widetilde\pi\) is at most \(C\sqrt d\). Let \(\mu_k\) be the actual law after \(k\) steps, and set \(a_k=W_2(\mu_k,\widetilde\pi)\). Then \[\left(\int |y|^2\,\mu_k(dy)\right)^{1/2} \le a_k+C\sqrt d.\] All these moments are finite. Indeed they are finite at \(k=0\), and the local \(W_2\) bound in (70), together with the conditional moment bounds for \(P_h\), propagates finite second moments to the next step.

Take the root mean square of the first local error in (70). Conditional coupling and Lemma 32 give \[a_{k+1}\le \left(\frac1{1+h}+2hB_d\delta^{K_0}\right)a_k +C\sqrt h\,B_d\delta^{K_0}\sqrt d.\] For all large \(d\), the coefficient of \(a_k\) is at most \(1-ch\) for an absolute \(c>0\). Since \(a_0\le C\sqrt d\), summing this geometric recurrence yields \[ a_k\le C\sqrt d\,e^{-chk} +Ch^{-1/2}B_d\delta^{K_0}\sqrt d. \tag{71}\] The second term is \(o(\sqrt h)\): after division by \(\sqrt h\), it is bounded by a constant times \(B_d\delta^{K_0-1}\sqrt d\), which tends to zero because \(e(K_0-1)>1/2\). At \(k=\lceil Ch^{-1}\log d\rceil\), a sufficiently large fixed \(C\) makes the first term \(o(\sqrt h)\) as well.

Use one more transition to obtain TV control. The mixture bound in Lemma 4 gives \[\lVert\mu_k\widehat P_h-\mu_kP_h\rVert_{\mathrm{TV}} \le B_d\delta^{K_0} \bigl[\sqrt d+2\sqrt h\,(a_k+C\sqrt d)\bigr]=o(1).\] For the remaining comparison, the conditional \(x'\)-laws from \(\mu_k\) and \(\widetilde\pi\) can be coupled with root mean square at most \(a_k/(1+h)\). Their common independent \(\sqrt h\)-Gaussian then gives, by (8), \[\lVert\mu_kP_h-\widetilde\pi P_h\rVert_{\mathrm{TV}} \le\frac{a_k}{(1+h)\sqrt h}=o(1).\] Invariance of \(\widetilde\pi\) completes the proof. ◻

Proof. The initial estimate follows from \(|g_y(0)|\le2\sqrt h\,|y|\) and the moment bound for \(\widetilde\pi\). Conditional on an ideal history, the next observation has the law of \(Z+N\) at center \(x_j\) and scale \(r_j=2^{-j/2}\). Equation (7) therefore gives \[\left(\mathbb E[|Y_j|^2\mid x_j,y]\right)^{1/2} \le C(\sqrt d+r_j|g_y(x_j)|).\] The center update and \(\operatorname{Lip}(g_y)\le\delta\) imply, by the conditional and then unconditional \(L^2\) triangle inequalities, \[b_{j+1}\le(1+C\delta r_j^2)b_j+C\delta r_j\sqrt d.\] Both \(\sum_{j\ge0}r_j\) and \(\sum_{j\ge0}r_j^2\) are bounded by absolute constants. In particular the products of \(1+C\delta r_j^2\) are bounded by \(\exp(C\delta\sum_jr_j^2)\). Expanding the recurrence proves the stated uniform bound. ◻

Proof. Translate by \(z_0\). Up to an additive constant, the potential of \(w=z-z_0\) is \(|w|^2/2+J(w)\), where \[J(w)=z_0\cdot w+F(x+r(z_0+w))-F(x+rz_0).\] The mode equation gives \(\nabla J(0)=0\), and \(\operatorname{Lip}(\nabla J)\le l\). In particular \(J(0)=0\) and \(|J(w)|\le l|w|^2/2\). For \(0\le s\le1\), the potential \(|w|^2/2+sJ(w)\) has mode zero and Hessian between \((1-l)I_d\) and \((1+l)I_d\). Its law \(\mu_s\) satisfies, by Lemma 3, \[\operatorname{Var}_{\mu_s}J \le\frac1{1-l}\mathbb E_{\mu_s}|\nabla J|^2 \le\frac{l^2}{1-l}\mathbb E_{\mu_s}|w|^2 \le\frac{l^2d}{(1-l)^2}\le4l^2d.\] The same quadratic bounds justify differentiating \[\Phi(s)=\log\int_{\mathbb R^d}e^{-|w|^2/2-sJ(w)}\,dw\] twice under the integral, and give \(\Phi''(s)=\operatorname{Var}_{\mu_s}J\). Since \(\mu_0=N(0,I_d)\), \[\operatorname{KL}(\mu_0\Vert\mu_1) =\Phi(1)-\Phi(0)-\Phi'(0) =\int_0^1(1-s)\Phi''(s)\,ds\le2l^2d.\] Pinsker’s inequality, followed by translation back by \(z_0\), proves the assertion. ◻

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

Every invocation of Theorem 2 uses at most \(\operatorname{polylog}(d)\delta^{-t}\) evaluations of the gradient oracle \(g_y\). By (69), each evaluation is one original query, including calls at random descent centers. The cap is uniform in those centers. There are \(O(\delta^{-1}\log d)\) proximal calls, \(O(\log d)\) descent calls, and \(O(\log d)\) additional gradient queries in the base rule. Summing gives the deterministic cap \[ \operatorname{polylog}(d)\delta^{-(1+t)}=d^{e(1+t)+o(1)}. \tag{78}\] The query counts and recursion shapes of this construction are fixed by the scalar parameters; they do not depend on oracle magnitudes or approximation events. All maps are measurable by the primitive theorem and the explicit finite operations above. Since \(t<1\), we have \(e(1+t)<2e<\varepsilon\). The fixed polylogarithmic factors fit \(C_\varepsilon d^\varepsilon\) for all sufficiently large \(d\). Equation (77) gives TV error at most \(1/10\) above a uniform dimension threshold.

For each of the finitely many smaller dimensions, use independent standard Gaussian proposals. Accept a proposal \(X\) with probability \[\exp\bigl(-V(X)+|X|^2/2\bigr).\] Integrating the Hessian bounds from the known minimizer gives \(|x|^2/2\le V(x)\le|x|^2\). The acceptance probability is therefore at most one, and its average is at least \(\mathbb Ee^{-|X|^2/2}=2^{-d/2}\). The proposal density times its acceptance probability is proportional to \(e^{-V(x)}\), so an accepted proposal has law \(\pi_V\).

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

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

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

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

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

Let \(x\) be the next query, which is fixed given that history. Write \(x=x_W+x_\perp\) according to \(W\oplus W^\perp\). The preimage of \(x_W\) is known. If \(x_\perp\ne0\), put \[w=\frac{x_\perp}{|x_\perp|},\qquad v=O^Tw\in E^\perp.\] Conditionally on the history, \(v\) is uniform on the complementary unit sphere. Let \(h\) be a standard Gaussian vector in \(E\), and let \(a\) have the distribution of the Euclidean norm of a standard Gaussian in \(\mathbb R^{d-\dim E}\). Draw \(h\) and \(a\) independently of one another and of the unrevealed map, conditional on the history, and define \[g_i=h+av.\] The Gaussian polar decomposition shows that \(g_i\) is a standard Gaussian in \(\mathbb R^d\) conditional on the entire previous history. Moreover, \(a>0\) almost surely, so \(v=(g_i-h)/a\in E+\operatorname{span}\{g_i\}\). Reveal \(Ov=w\) and extend \(T\) to \(E+\operatorname{span}\{v\}\). Conditioning on the full vector \(g_i\) does not reveal any further part of \(O\): its parallel component \(h\) and complementary radius \(a\) are independent auxiliary variables, and its normalized complementary component is precisely \(v\). Thus the conditional Haar invariant persists. If \(x_\perp=0\), instead draw a fresh independent standard Gaussian \(g_i\) and make no extension at this step.

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

Every \(g_i\) has conditional law \(N(0,I_d)\) given all earlier Gaussian vectors and the earlier history. Consequently these vectors are independent. To check the polynomial degrees, let \(E_i\) be the revealed eigenframe space after reply \(i\). Induction gives \[E_i\subseteq\operatorname{span} \{\Lambda^j g_a:1\le a\le i,\ 0\le j\le i\}.\] Indeed, the new query preimage is in \(E_{i-1}+\operatorname{span}\{g_i\}\), and adjoining its reply applies \(\Lambda\) once. All old directions remain in the indicated larger span.

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

Proof. Let \(m=(q+1)^2\) and form a \(d\times m\) matrix \(B\), whose columns are indexed by \((i,j)\) with \(1\le i\le q+1\) and \(0\le j\le q\), by \[B_{k,ij}=d^{-1/2}(g_i)_k c_j\cos(j\theta_k), \qquad c_0=1,\quad c_j=\sqrt2\quad(j\ge1).\] The functions \(\cos(j\theta)\) are degree-\(j\) polynomials in \(\cos\theta\) with nonzero leading coefficients. Since \(\cos\theta_k=2\lambda_k-3\), the range of \(B\) is precisely \(\mathcal K_q\).

Discrete cosine orthogonality gives \[ \mathbb EB^TB=I_m,\qquad \mathbb EB^TD_pB=0. \tag{82}\] For completeness, the equally spaced half-grid satisfies \[ \sum_{k=1}^d\cos(\ell\theta_k) =\frac{\sin(\ell\pi)}{2\sin(\ell\pi/(2d))}=0, \qquad 1\le\ell<2d,\quad \ell\in\mathbb Z. \tag{83}\] Here the finite cosine sum follows by summing a geometric progression, and its denominator is nonzero in the indicated range. Expand products of two cosines for the first identity in (82). For the second, expand the product of three cosines with frequencies \(n,j,j'\). Its frequencies are nonzero because \(n>j+j'\), and have absolute value at most \(4q+1<2d\). Different Gaussian blocks have zero cross means in both cases.

An entry of either Gram matrix is a sum of independent row terms. The deterministic coefficient in each term has absolute value at most \(2/d\), and the variance of a product of two standard Gaussian entries is at most \(2\). Thus every entry has variance at most \(8/d\), and \[\mathbb E\bigl[\|B^TB-I_m\|_{\mathrm F}^2+ \|B^TD_pB\|_{\mathrm F}^2\bigr] \le\frac{16m^2}{d}.\] By Markov’s inequality, \[\mathbb P\!\left(\|B^TB-I_m\|_{\mathrm F}^2+ \|B^TD_pB\|_{\mathrm F}^2>\tau_d^2\right) \le\frac{16m^2/d}{1600m^2/d}=0.01.\] Since \(\|C\|_{\mathrm{op}}\le\|C\|_{\mathrm F}\), outside this event both operator norms are at most \(40m/\sqrt d=\tau_d\). Whenever these bounds hold, for every \(a\in\mathbb R^m\), \[|Ba|^2\ge(1-\tau_d)|a|^2\ge\tfrac12|a|^2, \qquad |(Ba)^TD_p(Ba)|\le\tau_d|a|^2.\] These inequalities imply (81) on the entire range of \(B\). ◻

Proof. Independence of the coordinates gives \(\mu_{d,q}=\sum_k p_k/\lambda_k\). Substitute (84). The constant term vanishes by (83). For \(1\le j<2d-n\), use \[\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 have nonzero absolute value below \(2d\), so their sums vanish. If \(j=n\), the constant term contributes \(d/2\) and the frequency \(2n<2d\) contributes zero. Thus only \(j=n\) contributes in this range. Each remaining cosine product has absolute sum at most \(d\), so summing the geometric tail proves (85). Finally, a centered normal variable with variance \(v\) satisfies \(\operatorname{Var}(U_k^2)=2v^2\). Since \(|p_k|\le1\) and \(\lambda_k\ge1\), summing these coordinate variances proves (86). ◻

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

Fix any algorithm using at most \(q\) queries on every execution, where \[0\le q\le c_0\log d,\qquad c_0=\frac1{16\log(1/\varrho)}.\] All dimension restrictions in Lemmas 38, 39, and 40 hold for sufficiently large \(d\), uniformly over these \(q\). Let \(\mathbb P_{\mathrm{alg}}\) be the joint law of the Haar rotation and the algorithm’s output on this input. Define the event on pairs \((O,x)\) \[ \mathcal E_{d,q}= \left\{(O,x): |O^Tx|^2\le2d,\quad \left|(O^Tx)^TD_p(O^Tx)\right|> 160(q+1)^2\sqrt d\right\}. \tag{87}\] This event is measurable. It is also determined by \((A,x)\): if \(T_n\) is the polynomial defined by \(T_n(\cos\theta)=\cos(n\theta)\), then \(OD_pO^T=T_n(2A-3I)\). We retain \(O\) in the notation to use eigenframe coordinates throughout the comparison.

Apply Lemma 38 to the algorithm’s output and then Lemma 39 to its auxiliary Gaussian vectors. On the latter lemma’s event, the norm condition in (87) forces the quadratic form to be at most \(2\tau_d(2d)=160(q+1)^2\sqrt d\). Hence \[ \mathbb P_{\mathrm{alg}}(\mathcal E_{d,q})\le0.01. \tag{88}\]

For the target joint law \(\mathbb P_{\mathrm{target}}\), sample \(O\) from the same prior and then \(X\) from \(N(0,O\Lambda^{-1}O^T)\). Conditional on \(O\), the vector \(U=O^TX\) has law \(N(0,\Lambda^{-1})\). Our choice of \(c_0\) gives \[\varrho^{2q+1}\ge\varrho d^{-1/8},\qquad \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 in (85) satisfies \[\frac{|R_{d,q}|}{d\varrho^n/\sqrt2} \le\frac{2\varrho^{2d-2n}}{1-\varrho}\longrightarrow0\] uniformly, since \(n\le1+2c_0\log d\). Thus it is negligible relative to the leading term. In particular, with \(T_{d,q}=160(q+1)^2\sqrt d\), we have \(|\mu_{d,q}|>T_{d,q}\) for large \(d\), and \[\mathbb P_{\mathrm{target}}(\mathcal E_{d,q}^c) \le\frac2d+ \frac{2d}{(|\mu_{d,q}|-T_{d,q})^2}=o(1)\] uniformly, by Lemma 40 and Chebyshev’s inequality. Thus the probabilities of \(\mathcal E_{d,q}\) under the two joint experiments differ by at least \(0.99-o(1)\).

To compare these joint laws with the required guarantee on each fixed input, let \(K_O\) denote the algorithm’s output law on \(A=O\Lambda O^T\) and let \(G_O=N(0,A^{-1})\). For a measurable joint event \(\mathcal E\), write \(\mathcal E_O=\{x:(O,x)\in\mathcal E\}\). Both laws use the same Haar prior, so \[\begin{align*} \left|\mathbb P_{\mathrm{alg}}(\mathcal E) -\mathbb P_{\mathrm{target}}(\mathcal E)\right| &=\left|\int\bigl[K_O(\mathcal E_O)-G_O(\mathcal E_O)\bigr] \,d\operatorname{Haar}(O)\right|\\ &\le\int\|K_O-G_O\|_{\mathrm{TV}}\,d\operatorname{Haar}(O). \end{align*}\] The pointwise sampling guarantee would bound the last integral by \(1/10\). The event used to disprove that guarantee is a test of the joint laws; the algorithm need not know or compute it. This is impossible for (87). Therefore no such algorithm with \(q\le c_0\log d\) can meet the required accuracy for all inputs. Decreasing the constant if necessary proves the theorem. ◻