Subpolynomial queries for log-concave sampling

Source and verification

This article explains family 139 of the October 6, 2026 OpenAI collection. Research conclusions are attributed to its manuscripts. The review covers statements, introductions, and selected proof mechanisms; it is not an independent verification of the complete proofs.

🏷️ Preliminaries: log-concavity and sampling accuracy

A probability density proportional to $e^{-V(x)}$ is log-concave when $V$ is convex. Here the potential is twice continuously differentiable on ${\mathbb{R}}^d$, and the stronger assumptions are

$$\begin{array}{c}V(0)=0,\nabla V(0)=0,\\ I_d⪯{\nabla }^{2}V(x)⪯2I_d\text{for all }x.\end{array}$$

The matrix inequalities bound every directional curvature between one and two. The minimizer is known to be the origin; it is not an unknown point that the algorithm must first locate. Taylor’s integral formula gives

$$\frac{1}{2}|x{|}^2\le V(x)\le |x{|}^2,$$

so the normalizing integral is finite. Write the target law as

$${\pi }_V(dx)=Z_V^{-1}{e}^{-V(x)}dx.$$

Total variation measures the largest discrepancy in probabilities of events:

$$\Vert \mu -\nu {\Vert }_{\mathrm{TV}}=\underset{B}{\sup }|\mu (B)-\nu (B)|.$$

The required accuracy is $1/10$. This is a distributional guarantee, stronger in a different direction than merely matching the mean or returning a point near the minimizer. The prerequisites are graduate probability, Gaussian conditioning, multivariable calculus, and elementary operator norms.

🏷️ Exact oracle complexity and exponent zero

A first-order oracle query at $x$ returns the exact pair $(V(x),\nabla V(x))$. Queries may be adaptive. The model allows random real variables and unrestricted computation between queries; all algorithms must be measurable, and their only information about $V$ comes through oracle replies.

Let $Q(d)$ be the smallest worst-case query cap sufficient to return a sample within total variation $1/10$ for every potential in this class. The cap applies to every execution, not just to the expected number of queries.

Manuscript result

For every $\epsilon >0$ there is a constant $C_{\epsilon }$ such that $Q(d)\le C_{\epsilon }{d}^{\epsilon }$ for all $d\ge 2$. Conversely, $Q(d)\ge c\log d$ for all sufficiently large $d$, even for arbitrary randomized adaptive algorithms in this exact oracle model. (OpenAI, 2026)

Thus the optimal dimension exponent is zero:

$$\underset{d\to \infty }{\mathrm{limsup}}\frac{\log Q(d)}{\log d}=0.$$

The lower bound rules out a dimension-independent query cap. The upper bound does not establish a polylogarithmic bound: a function such as $\exp (\sqrt{\log d})$ is also smaller than every fixed positive power of $d$, while exceeding every fixed power of $\log d$.

The constants and approximation orders can depend on $\epsilon$. No arithmetic-operation or finite-bit running-time bound follows. In particular, this result does not say that a practical implementation costs less than reading or writing a $d$-dimensional vector.

🏷️ A Gaussian example and the source of difficulty

If $V(x)=x^{T}Ax/2$, the target is $N(0,A^{-1})$, with $I⪯A⪯2I$. A gradient query returns $Ax$, and the value is already determined by $x$ and that reply. Thus the oracle becomes matrix-vector multiplication.

If $A$ were supplied explicitly, one could factor it and transform a standard Gaussian. In this oracle problem the matrix is hidden. An algorithm must acquire information through its chosen vectors. The lower bound already occurs in this Gaussian subclass, despite its known center and bounded condition number.

For the general upper bound, the challenge is to approximate nonlinear conditional laws without assuming high derivatives of $V$. Standard high-order expansions cannot simply assume regularity beyond the stated $C^2$ hypotheses.

🏷️ Proximal conditioning creates a small perturbation

Suppose $X\sim {\pi }_V$ and $Y=X+\sqrt{h}G$, where $G$ is an independent standard Gaussian. The conditional law of $X$ given $Y=y$ has density

$$\begin{array}{rl}{\kappa }_y(dx)& \propto e^{-U_y(x)}dx,\\ U_y(x)& =V(x)+\frac{|x-y{|}^2}{2h}.\end{array}$$

Writing $x=y+\sqrt{h}z$ puts it into the form

$$\begin{array}{rl}{\kappa }_y(dz)& \propto e^{-|z{|}^2/2-F_y(z)}dz,\\ F_y(z)& =V(y+\sqrt{h}z).\end{array}$$

Its perturbation has gradient Lipschitz constant at most $2h$. Indeed,

$$\nabla F_y(z)=\sqrt{h}\nabla V(y+\sqrt{h}z),$$

and the chain rule contributes another $\sqrt{h}$ to the Hessian. One query to the original gradient evaluates this transformed gradient.

Small $h$ therefore turns each conditional sampling task into a small perturbation of a standard Gaussian. The manuscript develops a reusable primitive for such laws, including rescaled and shifted variants. It returns a conditional sample with added Gaussian noise and an appropriately noised conditional gradient mean.

The noise is part of the primitive’s output law. It is needed for the recursive implementation and final smoothing estimates; it cannot simply be deleted while preserving the theorem.

🏷️ Gaussian smoothing and controlled high derivatives

For a near-Gaussian variable $Z$, interpolate with a standard Gaussian $N$ by

$$Y_{\rho }=\rho Z+\sqrt{1-{\rho }^2}N.$$

The velocity transporting its marginal laws can be written in terms of conditional gradient means. Those conditional means belong to the same auxiliary class, allowing recursive evaluation.

High-order approximation is useful only if each derivative does not cost another power of dimension. The manuscript bounds derivative tensors after every division of their coordinate slots into two nonempty groups, treating each division as a matrix. These operator-norm bounds control the contractions that occur when derivative expressions are composed.

Repeated Gaussian integration by parts then introduces polylogarithmic dimension losses, followed by one $\sqrt{d}$ factor for the final vector output. The higher derivatives come from Gaussian conditioning, not from an added assumption that the original potential has bounded derivatives of every order.

A second construction couples a smoothed position law to a Gaussian momentum while preserving their product law. A specially chosen skew-symmetric drift converts an integral of conditional gradient means into a Gaussian displacement centered at the desired mean. High-order quadrature approximates this transport and centering mechanism.

These tensor and transport estimates are the analytic input to the claimed complexity bound. The elementary proximal change of variables alone does not give a subpolynomial-query algorithm.

🏷️ Recursive implementation without an exact-mean oracle

Finite Picard iterations and quadrature express the required transports in terms of conditional means. Merely writing such an expression is insufficient: a conditional mean is not an available oracle query.

The implementation uses a pathwise translation identity. Shifting the center of a routine can be compensated by a specified shift of its Gaussian seeds, leaving every primitive query unchanged. A deterministic covariance adjustment then absorbs independent Gaussian noise in the center. This converts the noisy-center call into the intended known-center law.

The identity is finite-dimensional and exact. It avoids treating an approximately known center as if it were exact, or silently assuming that errors from nested calls are independent.

Two depths must be distinguished: the dependency depth of the finite expression in conditional means, and the height of the routines used to evaluate those means. The construction bounds both before selecting its quadrature mesh. Error targets and perturbation sizes are propagated through the recursion; choosing sufficiently high fixed approximation orders makes the resulting query exponent arbitrarily small.

🏷️ Contraction, total variation, and removal of the noise

Consider the ideal transition that samples $X^{\prime }\sim {\kappa }_y$ and returns $X^{\prime }+\sqrt{h}N$. It preserves the noised law

$$\tilde{\pi }={\pi }_V\ast N(0,h{I}_d).$$

This follows by viewing the transition as two successive conditional draws in the joint law of $(X,Y)$. The $y$-marginal is precisely $\tilde{\pi }$.

Let $W_2$ denote the minimum root mean square distance over couplings. The conditional potential has curvature at least $1+1/h$, and changing $y_1$ to $y_2$ changes its gradient by a constant vector of norm $|y_1-y_2|/h$. Synchronous comparison yields

$$W_2({\kappa }_{y_1},{\kappa }_{y_2})\le \frac{|y_1-y_2|}{1+h}.$$

Using the same subsequent Gaussian preserves this contraction. Approximately $h^{-1}\log d$ transitions suffice for the required contraction scale, provided local errors are small enough to survive accumulation.

The paper takes $h=d^{-e}/2$ for a small fixed $e>0$. With $\delta =d^{-e}$, sufficiently high order $K$ produces error factors of the form

$$(1+\log d{)}^C{\delta }^K\sqrt{d}.$$

Increasing $K$ offsets the final dimensional factor and the losses from summing local errors. The source also controls dependence on the current center and propagates second moments; a uniform estimate only at the origin would not suffice along the chain.

One extra Gaussian-smoothed step turns the accumulated $W_2$ estimate into total-variation control at the appropriate noise scale. But the result so far concerns $\tilde{\pi }$, not ${\pi }_V$.

To remove the added noise, the construction repeatedly conditions and halves the remaining variance. The terminal conditional law is sufficiently close to a Gaussian centered at its mode. A final Gaussian approximation completes the descent to the unnoised target. Total-variation data processing controls the initial error when this same conditional procedure is applied to an approximate noised sample.

🏷️ Why logarithmically many queries remain necessary

The lower bound takes a diagonal matrix $\Lambda$ with eigenvalues in $[1,2]$ and hides its eigenvectors with a random orthogonal rotation $O$, giving $A=O\Lambda O^T$.

Deferred revelation of this rotation shows that, after $q$ adaptive queries, even an output outside the observed query span admits the representation

$$\begin{array}{c}{O}^{T}X\in {\mathcal{K}}_q,\\ {\mathcal{K}}_q=\mathrm{span}\{{\Lambda }^j{g}_i{\}}_{i,j}.\end{array}$$

Here $1\le i\le q+1$ and $0\le j\le q$. The auxiliary $g_i$ are mutually independent standard Gaussian vectors in an extension of the experiment. They need not be independent of the hidden matrix. The extra vector accounts for a previously unobserved direction chosen at output, so the argument is not restricted to algorithms that return a linear combination of past queries.

A cosine statistic is small uniformly on this enlarged span when $q$ is too small. Its frequency is chosen beyond twice the polynomial degree accessible in the displayed representation. The same statistic has a detectable bias under $N(0,A^{-1})$. Comparing these two experiments yields the $c\log d$ lower bound for arbitrary measurable randomized adaptive algorithms.

The theorem therefore determines the optimal power of dimension, while leaving the finer growth rate between logarithmic and subpolynomial bounds unresolved.

Concentration of information for log-concave distributions develops complementary distributional background. The graph switch chain studies total-variation mixing in a finite combinatorial state space, where transition cost and mixing time have a different computational interpretation.

🏷️ Key steps and proof outlines

An exact identity replaces the unavailable mean oracle

The recursive sampler needs conditional means at centers known only with Gaussian noise. The implementation section supplies an exact finite-dimensional remedy. Let $\mathsf{Q}(x;z_1,\dots ,z_m)$ be a deterministic seeded routine at scale $r$, and suppose a fixed vector $v\in {\mathbb{R}}^m$ satisfies

$$\mathsf{Q}(x+\Delta ;z+v\otimes \Delta /r)=\mathsf{Q}(x;z)$$

for every center, displacement, and seed list, with exactly the same physical oracle queries. Here $(v\otimes a{)}_j=v_{j}a$. This is a pathwise identity, so it remains valid for correlated seeds.

Choose $D_s>0$ with $D_s\ge \Vert v{\Vert }_2$, put $n=r/(2D_s)$, and set

$$P={\left(I_m-\frac{v{v}^{\mathsf{T}}}{4D_s^2}\right)}^{1/2}.$$

Its eigenvalues lie between $\sqrt{3}/2$ and one. For independent standard Gaussian vectors $Z,G_1,\dots ,G_m$, the translation identity rewrites the noisy-center call as

$$\mathsf{Q}(x+nZ;PG)=\mathsf{Q}\left(x;PG-\frac{v\otimes Z}{2D_s}\right).$$

The seed array on the right is centered Gaussian with slot covariance $P{P}^{\mathsf{T}}+v{v}^{\mathsf{T}}/(4D_s^2)=I_m$. Its slots therefore have exactly the standalone independent standard Gaussian law. They may remain correlated with $Z$, but the rewritten center is deterministic, so that correlation is irrelevant. Without the covariance reduction $P$, translating the center back would add unwanted seed variance. (OpenAI, 2026)

Recursion must preserve identities before using probability

At each routine level, the compiler first proves the query-preserving shift property for arbitrary fixed seeds. It builds child shift vectors before their covariance corrections and the parent’s shift vector. A separate deterministic stability estimate controls changes of center even for transformed or correlated seed arrays. Only after Gaussian absorption restores a child’s full standalone seed law is its probabilistic accuracy theorem invoked.

This order separates three claims that otherwise look similar: equality of physical queries under a pathwise shift, stability under an imperfect center, and distributional accuracy under independent Gaussian seeds. The first two need no independence. The third does, and the covariance identity supplies it at the precise point of use. Bounded routine height and bounded shift-vector norms then make the number of seeds and oracle calls a deterministic cap rather than an expected cost.

Contraction and Gaussian smoothing solve different error problems

For the ideal proximal conditional kernel, strong convexity gives the contraction factor $(1+h{)}^{-1}$ in $W_2$. If each approximate step has $W_2$ error at most $\eta$, coupling and the triangle inequality give the schematic recurrence

$$e_{j+1}\le \frac{e_j}{1+h}+\eta ,e_j\le (1+h{)}^{-j}{e}_0+\frac{1+h}{h}\eta .$$

Thus local error must be small compared with $h$ times the desired global error. The actual source estimates also track the center’s moments, because its local bounds depend on where the chain is. High fixed approximation order makes powers of the perturbation scale overcome the final $\sqrt{d}$ factor and the accumulated $h^{-1}$ loss.

A $W_2$ bound by itself does not imply total variation control. Added Gaussian noise supplies that conversion. Conditional on a coupling $(X,Y)$, the relative entropy between $N(X,hI)$ and $N(Y,hI)$ is $\Vert X-Y{\Vert }^2/(2h)$. Pinsker’s inequality and averaging therefore bound the total variation of their mixtures by $\mathbb{E}\Vert X-Y\Vert /(2\sqrt{h})\le W_2/(2\sqrt{h})$. This explains the additional noise-scale accuracy reserve.

The final target requires a separate noise-removal argument

The transition preserves ${\pi }_V\ast N(0,hI)$, not ${\pi }_V$. The source’s successive conditional refinements reduce the added variance and finish with a Gaussian approximation at a sufficiently small scale. Applying the same refinement kernel to both the ideal and approximate noised laws cannot increase their total variation distance, by data processing. The refinement’s own error is then added to the initial discrepancy.

This separates three distinct obligations: contract toward the noised stationary law, convert the chain’s metric error to total variation, and remove the noise. The Gaussian subclass lower bound addresses yet another issue: how much information an adaptive exact oracle reveals about hidden eigendirections. Its deferred-rotation representation includes an extra output direction, so it is not restricted to algorithms returning a vector in the observed query span. Neither the upper nor lower argument supplies a finite-bit implementation bound in this unrestricted arithmetic query model.

🏷️ Manuscripts and verification scope

• Subpolynomial query complexity for well-conditioned log-concave sampling — pinned manuscript and companion source.

The family scope page describes selected formalization coverage. It selects the exponent-zero statement with its subpolynomial upper bounds and eventual logarithmic lower bound in the exact adaptive value-and-gradient oracle model. Arithmetic and bit costs are outside this scope. No Lean build, axiom audit, or manuscript-to-formal-statement equivalence check was performed for these notes.

OpenAI collection provides prerequisites and the reading sequence; coverage register records completed and pending families.

🐻  OpenAI 2026. Subpolynomial query complexity for well-conditioned log-concave sampling.