Uniform-in-Time Particle Convergence
for Projected Gaussian SVGD

Abstract

Liu, Ghosal, Balasubramanian, and Pillai [12] conjectured uniform-in-time mean-field convergence for projected Gaussian Stein variational gradient descent (SVGD) beyond Gaussian targets. We give a quantitative affirmative answer for the continuous-time, exact-Gaussian-moment system with kernel $K_{1}(x,y)=x^{\top}y+1$ , assuming $\pi(dx)\propto e^{-V(x)}dx$ , $V\in C^{2}(\mathbb{R}^{d})$ , and $0<\alpha I\preceq\nabla^{2}V\preceq\beta I$ . For Gaussian i.i.d. initialization and every $N\geq d+1$ , the empirical particle law $\zeta_{N,t}$ and Gaussian mean-field law $\rho_{t}$ satisfy $\mathbb{E}[\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})]\leq Cr_{d}(N)$ . The asymptotic orders of $r_{d}(N)$ are $N^{-1}\log\log N$ , $N^{-1}(\log N)^{2}$ , and $N^{-2/d}$ in dimensions one, two, and at least three, respectively; each is sharp in particle number. Thus the time supremum can be placed inside the expectation without losing the Gaussian empirical-measure rate, even at the minimal full-rank sample size. We also obtain high-probability all-time bounds and a last-iterate guarantee to the reverse-KL optimal Gaussian at the same statistical scale. The key deterministic inequality separates initial moment mismatch from the Wasserstein discrepancy of the standardized cloud, using stability in square-root covariance coordinates and exact affine transport. At fixed $N$ , we characterize the discrete limiting cloud: its moments are optimal, but its preserved standardized shape leaves a strictly positive full-law error.

Note. This paper was generated entirely by AI, including MiMo, using an automated research pipeline developed by Chenghua Liu and Hanyu Li.

1 Introduction and prior conjecture

This section states the prior conjecture, explains the obstacle to uniform full-law control, and summarizes our resolution. Detailed definitions and results follow in Section 2.

Stein variational gradient descent (SVGD) transports particles by a kernelized deterministic velocity field [11, 10]. Finite-horizon propagation-of-chaos constants can grow with the horizon, and time-averaged bounds do not control a physical-time (fixed-time or last-iterate) law; uniform-in-time full-law control is therefore delicate. Our result concerns a continuous-time Gaussian-score projection with exact Gaussian expectations, not the ordinary unprojected SVGD particle system.

In the arXiv version of Liu et al. [12], Theorem 3.7 gives a Gaussian-target, pointwise-in-time expectation bound with a constant uniform in time under Gaussian i.i.d. initialization. For general targets, their Theorem 4.3 derives the exact-moment affine representation and is followed by a qualitative conjecture of uniform-in-time mean-field convergence. We formulate a stronger quantitative version, making the sampling model explicit.

Open problem 1.1 (Quantitative exact-moment conjecture).

For a possibly non-Gaussian target satisfying $V\in C^{2}(\mathbb{R}^{d})$ and $\alpha I\preceq\nabla^{2}V\preceq\beta I$ , initialize the projected $K_{1}(x,y)=x^{\top}y+1$ particle system with i.i.d. samples from a full-rank Gaussian law $\rho_{0}$ . At fixed dimension, can one place the time supremum inside the expectation, retain the sharp Gaussian empirical-measure rate, and obtain a physical-time last-iterate bound to the reverse-KL optimal Gaussian?

Exact moment closure alone does not answer this question. The empirical law retains a standardized shape that affine dynamics transports rather than Gaussianizes, while the moment vector field becomes sensitive near the positive-definite boundary. We resolve both issues under global strong convexity and a bounded continuous Hessian, without assuming a third derivative or a Hessian-Lipschitz modulus:

1.
a square-root-coordinate energy yields semiglobal all-time stability on every energy sublevel;
2.
a deterministic oracle separates moment mismatch from a time-independent standardized-shape discrepancy;
3.
Gaussian matching and log-Wishart estimates give the sharp expected $N$ -rate, a high-probability upper bound, and the endpoint $N=d+1$ ; and
4.
the exact finite-particle limit retains a positive shape floor.

A scalar counterexample delineates the scope of the first item: the semiglobal Euclidean Lipschitz estimate used here cannot be replaced by one initialization-uniform constant over the full positive-definite cone.

2 Model and main results

This section defines the projected moment and particle dynamics, then states the deterministic oracle, its probabilistic consequences, and the exact finite-particle limit. We close with the assumptions that delimit their scope.

Let $d\geq 1$ be fixed and let

\pi(dx)=Z_{V}^{-1}e^{-V(x)}\,dx,\qquad V\in C^{2}(\mathbb{R}^{d}),

where $Z_{V}<\infty$ by strong convexity and, for constants $0<\alpha\leq\beta<\infty$ ,

\alpha I\preceq\nabla^{2}V(x)\preceq\beta I,\qquad x\in\mathbb{R}^{d}.

(2.1)

We use $W_{2}$ for the quadratic Wasserstein distance on $\mathcal{P}_{2}(\mathbb{R}^{d})$ and $T_{\#}\nu$ for pushforward by $T$ . Vector norms are Euclidean; matrix norms carry subscripts $\mathrm{F}$ or $\mathrm{op}$ . For $q=(\mu,A)\in\mathbb{R}^{d}\times\Sym(d)$ , the product norm is $\left\lVert q\right\rVert^{2}=\left\lvert\mu\right\rvert^{2}+\left\lVert A\right\rVert_{\mathrm{F}}^{2}$ . For every $\Sigma\succ 0$ , $\Sigma^{1/2}$ denotes its principal symmetric square root, and

\mathcal{Q}=\mathbb{R}^{d}\times\mathbb{S}_{++}^{d},\qquad q=(\mu,A)=(\mu,\Sigma^{1/2}).

We also write $\gamma_{d}=\mathcal{N}(0,I_{d})$ , set $u_{+}=\max\{u,0\}$ and $\log_{+}x=(\log x)_{+}$ , and define $d_{\rm TV}(P,Q)=\sup_{B}\left\lvert P(B)-Q(B)\right\rvert$ . Unless indexed otherwise, $C,c>0$ may change from line to line and depend only on the fixed problem data in the relevant statement, never on $N, t$ , or a tail parameter.

For $(\mu,\Sigma)\in\mathbb{R}^{d}\times\mathbb{S}_{++}^{d}$ , define the exact Gaussian moments

m(\mu,\Sigma)=\mathbb{E}_{\mathcal{N}(\mu,\Sigma)}[\nabla V(X)],\qquad\Gamma(\mu,\Sigma)=\mathbb{E}_{\mathcal{N}(\mu,\Sigma)}[\nabla^{2}V(X)].

(2.2)

The Gaussian mean-field moment flow is


$\displaystyle\dot{\mu}$	$\displaystyle=(I-\Gamma\Sigma)\mu-(1+\left\lvert\mu\right\rvert^{2})m,$	(2.3a)
$\displaystyle\dot{\Sigma}$	$\displaystyle=2\Sigma-\Sigma(\Sigma\Gamma+\mu m^{\top})-(\Gamma\Sigma+m\mu^{\top})\Sigma,$	(2.3b)

where $m=m(\mu_{t},\Sigma_{t})$ and $\Gamma=\Gamma(\mu_{t},\Sigma_{t})$ . We write $\rho_{t}=\mathcal{N}(\mu_{t},\Sigma_{t})$ for the solution starting from fixed $\mu_{0}\in\mathbb{R}^{d}$ and $\Sigma_{0}\succ 0$ .

Given initial points $X_{1}(0),\dots,X_{N}(0)$ , set

\widehat{\mu}_{t}=\frac{1}{N}\sum_{i=1}^{N}X_{i}(t),\qquad\widehat{\Sigma}_{t}=\frac{1}{N}\sum_{i=1}^{N}(X_{i}(t)-\widehat{\mu}_{t})(X_{i}(t)-\widehat{\mu}_{t})^{\top}.

(2.4)

With

\widehat{m}_{t}=m(\widehat{\mu}_{t},\widehat{\Sigma}_{t}),\qquad\widehat{\Gamma}_{t}=\Gamma(\widehat{\mu}_{t},\widehat{\Sigma}_{t}),

the projected particle system is

\dot{X}_{i}(t)=C_{t}X_{i}(t)-\widehat{m}_{t},\qquad C_{t}=I-\widehat{\Gamma}_{t}\widehat{\Sigma}_{t}-\widehat{m}_{t}\widehat{\mu}_{t}^{\top}.

(2.5)

This is the Gaussian-score projection associated with $K_{1}$ ; it is not the ordinary, unprojected SVGD particle system. Define

\zeta_{N,t}=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{i}(t)}.

For Theorems 2.2 and 2.4, the initial points are drawn independently from $\mathcal{N}(\mu_{0},\Sigma_{0})$ . Under this Gaussian initialization, all particle paths and standardized quantities are assigned an arbitrary fixed value on the null event $\{\widehat{\Sigma}_{0}\not\succ 0\}$ . On its complement the paths and Gaussian moments are continuous in $t$ , hence so is $t\mapsto W_{2}^{2}(\zeta_{N,t},\rho_{t})$ . Consequently

\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})=\sup_{t\in\mathbb{Q}_{+}}W_{2}^{2}(\zeta_{N,t},\rho_{t}),

(2.6)

so every random supremum below is measurable.

For $Z\sim\gamma_{d}$ , use the constant-free energy representative

F(\mu,A)=\mathbb{E}V(\mu+AZ)-\log\det A.

(2.7)

Then $\KL(\mathcal{N}(\mu,A^{2})\|\pi)=F(\mu,A)+c_{\pi}$ , where $c_{\pi}=\log Z_{V}-\tfrac{d}{2}(1+\log(2\pi))$ . Let $q_{*}=(\mu_{*},A_{*})$ be the unique minimizer, put $\Sigma_{*}=A_{*}^{2}$ and $\rho_{*}=\mathcal{N}(\mu_{*},\Sigma_{*})$ , and define

\mathcal{K}_{R}=\{q:F(q)-F(q_{*})\leq R\}.

(2.8)

Section 4 proves that these sublevels are compactly contained in the positive covariance domain.

We first state the deterministic result from which the probabilistic all-time estimates are lifted. Given a deterministic cloud $x_{1},\ldots,x_{N}$ with empirical mean $\widehat{\mu}_{0}$ and covariance $\widehat{\Sigma}_{0}\succ 0$ , choose any whitening matrix $R_{0}$ satisfying $R_{0}\widehat{\Sigma}_{0}R_{0}^{\top}=I$ and set

z_{i}=R_{0}(x_{i}-\widehat{\mu}_{0}),\qquad\eta_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{z_{i}}.

(2.9)

The value of $W_{2}(\eta_{N},\gamma_{d})$ does not depend on the whitening: two whitening matrices differ by a left orthogonal factor, and $\gamma_{d}$ is rotation invariant.

Theorem 2.1 (Deterministic localized oracle).

Let $q_{0}=(\mu_{0},\Sigma_{0}^{1/2})$ , let $q_{t}=(\mu_{t},\Sigma_{t}^{1/2})$ be the solution of (2.3), set $\rho_{t}=\mathcal{N}(\mu_{t},\Sigma_{t})$ , and let $\zeta_{N,t}$ be the projected particle solution from an arbitrary deterministic cloud with $\widehat{\Sigma}_{0}\succ 0$ . Put $q_{N,0}=(\widehat{\mu}_{0},\widehat{\Sigma}_{0}^{1/2})$ . For every $0\leq R<\infty$ there are finite constants $S_{R}$ and $b_{R}$ , depending only on $R, d, V$ , such that, whenever $q_{0},q_{N,0}\in\mathcal{K}_{R}$ ,

\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\leq 2S_{R}^{2}\left\lVert q_{N,0}-q_{0}\right\rVert^{2}+2b_{R}^{2}W_{2}^{2}(\eta_{N},\gamma_{d}).

(2.10)

Both the mean-field and particle solutions are global under these conditions. No independence or distributional assumption on the initial cloud is required.

In particular, any sequence of full-rank deterministic or dependent clouds lying in a common sublevel and satisfying $\left\lVert q_{N,0}-q_{0}\right\rVert^{2}+W_{2}^{2}(\eta_{N},\gamma_{d})\to 0$ converges uniformly in physical time.

For an unambiguous statement at every admissible $N$ , let

r_{d}(N)=\begin{cases}N^{-1}\{1+\log\log(e^{e}+N)\},&d=1,\\ N^{-1}(1+\log N)^{2},&d=2,\\ N^{-2/d},&d\geq 3.\end{cases}

(2.11)

Thus, asymptotically, the first two rates are $N^{-1}\log\log N$ and $N^{-1}(\log N)^{2}$ , respectively. The regularization inside the logarithms only makes these rates finite for every admissible sample size.

We now specialize the oracle to Gaussian sampling, first in expectation.

Theorem 2.2 (Uniform-in-time particle convergence).

Assume (2.1), fix $d\geq 1$ , $\mu_{0}\in\mathbb{R}^{d}$ , and $\Sigma_{0}\in\mathbb{S}_{++}^{d}$ , let $N\geq d+1$ , and draw the initial particles i.i.d. from $\mathcal{N}(\mu_{0},\Sigma_{0})$ . Then $\widehat{\Sigma}_{0}\succ 0$ almost surely, the systems (2.3) and (2.5) have global solutions, and there is a finite constant $C$ , depending only on $d,V,\mu_{0},\Sigma_{0}$ , such that

\mathbb{E}\!\left[\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\right]\leq Cr_{d}(N).

(2.12)

In particular, the same bound holds with the supremum outside the expectation. The constant is independent of $N$ and $t$ .

The next proposition shows that the expected rate cannot be improved in its $N$ -dependence, even before the dynamics starts.

Proposition 2.3 (Sharpness in particle number).

Under the Gaussian i.i.d. initialization of Theorem 2.2, for every fixed $d,\mu_{0}$ , and $\Sigma_{0}\succ 0$ there is $c>0$ such that

\mathbb{E}\!\left[\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\right]\geq cr_{d}(N),\qquad N\geq d+1.

(2.13)

Thus the $N$ -dependence in Theorem 2.2 is optimal in every fixed dimension.

Proof.

At time zero, $\zeta_{N,0}$ is the empirical measure of $N$ i.i.d. samples from $\rho_{0}=\mathcal{N}(\mu_{0},\Sigma_{0})$ . Applying the invertible affine map $x\mapsto\mu_{0}+\Sigma_{0}^{1/2}x$ to the standard Gaussian empirical matching problem, with $\xi_{i}\overset{\mathrm{iid}}{\sim}\gamma_{d}$ and $\nu_{N}=N^{-1}\sum_{i}\delta_{\xi_{i}}$ , gives

\mathbb{E}W_{2}^{2}(\zeta_{N,0},\rho_{0})\geq\lambda_{\min}(\Sigma_{0})\mathbb{E}W_{2}^{2}(\nu_{N},\gamma_{d}).

The expectation on the right is bounded below by $c_{d}r_{d}(N)$ : in dimension one this is the two-sided Gaussian order-statistic estimate of Bobkov and Ledoux [5, Corollary 6.14] for $\mathbb{E}W_{2}^{2}$ ; in dimension two it follows from the nonstandard two-sample Gaussian matching lower bound of Talagrand [16, main result and display (1)] (if $\nu_{N}^{\prime}$ is an independent copy of $\nu_{N}$ , then $W_{2}^{2}(\nu_{N},\nu_{N}^{\prime})\leq 2W_{2}^{2}(\nu_{N},\gamma_{d})+2W_{2}^{2}(\nu_{N}^{\prime},\gamma_{d})$ ); and for $d\geq 3$ it is the Gaussian matching lower bound of Ledoux and Zhu [9, Theorem 1.1]. Since time zero is included in the supremum, (2.13) follows. Reducing the constant handles the finitely many nonasymptotic sample sizes. ∎

The same decomposition also yields exponential concentration around the dimension-dependent statistical scale.

Theorem 2.4 (High-probability all-time convergence).

Under the assumptions of Theorem 2.2, there are constants $C,c,c_{0}>0$ , depending only on the fixed problem data, such that for every $N\geq d+1$ and $1\leq x\leq c_{0}N$ , there is an event $H_{N}(x)$ with $\mathbb{P}(H_{N}(x)^{c})\leq Ce^{-cx}$ on which

\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\leq C\left(r_{d}(N)+\frac{x}{N}\right).

(2.14)

On $H_{N}(x)$ , simultaneously for every $t\geq 0$ ,

W_{2}^{2}(\zeta_{N,t},\rho_{*})\leq C\left\{r_{d}(N)+\frac{x}{N}+e^{-2ct}\right\}.

(2.15)

Equivalently, for $Ce^{-cN}\leq\delta\leq 1/2$ , the statistical term is $C\{r_{d}(N)+N^{-1}\log(C/\delta)\}$ with probability at least $1-\delta$ , after changing the fixed constants.

Combining all-time particle control with convergence of the mean-field flow gives a last-iterate guarantee rather than a time-averaged one.

Corollary 2.5 (Physical-time last iterate).

Under the assumptions of Theorem 2.2, there are constants $C,c>0$ , depending only on the fixed problem data, such that for all $t\geq 0$ ,

\mathbb{E}W_{2}^{2}(\zeta_{N,t},\rho_{*})\leq C\{r_{d}(N)+e^{-2ct}\}.

(2.16)

Consequently, once $t\geq(2c)^{-1}\log_{+}(1/r_{d}(N))$ , the error is at the statistical scale $O(r_{d}(N))$ .

The deterministic dynamics has an exact limit at every fixed particle number; the residual discrepancy is precisely an affine shape error.

Proposition 2.6 (Exact finite-particle limit and shape floor).

Let a deterministic full-rank cloud and mean-field initialization satisfy $q_{N,0},q_{0}\in\mathcal{K}_{R}$ for some $R>0$ . There exist $C_{R}<\infty$ , $\lambda_{R}>0$ , and an invertible matrix $B_{\infty}$ with $B_{\infty}B_{\infty}^{\top}=\Sigma_{*}$ such that, for

\zeta_{N,\infty}=(\mu_{*}+B_{\infty}\,\cdot)_{\#}\eta_{N},

(2.17)

one has

	$\displaystyle W_{2}^{2}(\zeta_{N,t},\zeta_{N,\infty})$	$\displaystyle\leq C_{R}e^{-2\lambda_{R}t},$		(2.18)
	$\displaystyle\lim_{t\to\infty}W_{2}(\zeta_{N,t},\rho_{t})$	$\displaystyle=W_{2}(\zeta_{N,\infty},\rho_{*}).$		(2.19)

Moreover,

\lambda_{\min}(\Sigma_{*})W_{2}^{2}(\eta_{N},\gamma_{d})\leq W_{2}^{2}(\zeta_{N,\infty},\rho_{*})\leq\lambda_{\max}(\Sigma_{*})W_{2}^{2}(\eta_{N},\gamma_{d}).

(2.20)

The constants $C_{R}$ and $\lambda_{R}$ depend only on $R, d, V$ ; the matrix $B_{\infty}$ may depend on the initial cloud and the whitening coordinates, while the law $\zeta_{N,\infty}$ does not. Thus the moments converge exactly to $(\mu_{*},\Sigma_{*})$ , while every full-rank finite cloud retains a strictly positive full-law error.

This shape floor persists for a fixed non-Gaussian sampling law, even as the number of particles grows.

Corollary 2.7 (Non-Gaussian initialization is not Gaussianized).

On one probability space, let $(X_{i})_{i\geq 1}$ be an infinite i.i.d. sequence from a non-Gaussian law $P_{0}$ having a density, positive-definite covariance $\Sigma_{0}$ , and finite second moment. For each $N\geq d+1$ , initialize the particle system with $X_{1},\ldots,X_{N}$ and define $\zeta_{N,\infty}$ by Proposition 2.6 on its almost-sure full-rank event. Write

\eta_{0}=\bigl(x\mapsto\Sigma_{0}^{-1/2}(x-\mathbb{E}_{P_{0}}X)\bigr)_{\#}P_{0}.

Put $D_{0}=W_{2}(\eta_{0},\gamma_{d})>0$ . Then, almost surely,

\lambda_{\min}(\Sigma_{*})D_{0}^{2}\leq\liminf_{N\to\infty}W_{2}^{2}(\zeta_{N,\infty},\rho_{*})\leq\limsup_{N\to\infty}W_{2}^{2}(\zeta_{N,\infty},\rho_{*})\leq\lambda_{\max}(\Sigma_{*})D_{0}^{2}.

(2.21)

In particular, the limiting error is bounded away from zero; if $\Sigma_{*}=\sigma_{*}^{2}I$ , it converges to $\sigma_{*}^{2}D_{0}^{2}$ . Thus Gaussianity is a structural consistency condition for a fixed i.i.d. initialization law in Theorems 2.2 and 2.4, not merely a convenient source of concentration. This does not preclude consistent deterministic or $N$ -dependent designs covered by the oracle.

Remark 2.8 (Scope).

The probabilistic theorems concern continuous time, exact evaluation of (2.2), fixed dimension, Gaussian i.i.d. initialization, and the projected system (2.5). The deterministic oracle does allow arbitrary full-rank initialization. These are not dimension-free results or theorems for ordinary SVGD, Euler discretization, or Monte Carlo estimators of $m$ and $\Gamma$ . Since every centered $N$ -point covariance has rank at most $N-1$ , the restriction $N\geq d+1$ is necessary for this unregularized positive-definite dynamics.

3 Proof overview

This section gives the mechanism behind the results and records where each technical step enters. Complete arguments, including the endpoint and rare event calculations, are developed in Sections 4–7.

In square-root coordinates $q=(\mu,A)$ , the energy in (2.7) has second variation

D^{2}F(\mu,A)[(u,H),(u,H)]=\mathbb{E}[(u+HZ)^{\top}\nabla^{2}V(\mu+AZ)(u+HZ)]+\tr(A^{-1}HA^{-1}H),

(3.1)

which is at least $\alpha(\left\lvert u\right\rvert^{2}+\left\lVert H\right\rVert_{\mathrm{F}}^{2})$ . Its sublevels are compactly contained in $\mathcal{Q}$ . The moment flow can be written as $\dot{q}=-P(q)\nabla F(q)$ with a positive-definite mobility $P(q)$ ; on each sublevel, $P(q)\succeq p_{R}I$ . This yields exponential convergence to $q_{*}$ . For two trajectories, a Gronwall bound controls the common compact sublevel until both enter a fixed neighborhood of $q_{*}$ . There the linearization is dissipative in the $D^{2}F(q_{*})$ -norm, and patching the two regimes proves the all-time stability estimate. Section 4 and Section 5 give the details.

The particle equation is affine. Its empirical moments solve the same closed ODE, and a fundamental matrix $L_{t}$ gives

X_{i}(t)=\widehat{\mu}_{t}+L_{t}\{X_{i}(0)-\widehat{\mu}_{0}\},\qquad\widehat{\Sigma}_{t}=L_{t}\widehat{\Sigma}_{0}L_{t}^{\top}.

(3.2)

After whitening the initial cloud, both the empirical law and its moment- matched Gaussian are pushforwards by the same affine map. Consequently,

\lambda_{\min}(\widehat{\Sigma}_{t})W_{2}^{2}(\eta_{N},\gamma_{d})\leq W_{2}^{2}\!\left(\zeta_{N,t},\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t})\right)\leq\lambda_{\max}(\widehat{\Sigma}_{t})W_{2}^{2}(\eta_{N},\gamma_{d}).

(3.3)

Combining this with semiglobal moment stability proves the deterministic oracle. Exponential moment convergence also makes the affine coefficient integrable, producing the limit and the two-sided shape floor in Proposition 2.6. These arguments are in Section 6.

For Gaussian initialization, concentration of the sample mean and covariance places the empirical initial state in a deterministic common sublevel, while Gaussian empirical matching controls $W_{2}(\eta_{N},\gamma_{d})$ . This proves the high-probability theorem directly from the oracle. For the expectation, the rare nearly singular event cannot be discarded; energy decrease reduces its contribution to a log-determinant moment. Bartlett’s decomposition gives a sum of $\log\chi_{k}^{2}$ variables whose second moments remain bounded even when $N=d+1$ . Cauchy–Schwarz then makes the rare-event contribution exponentially small. Section 7 contains the full localization and matching argument.

4 Reverse-KL geometry in square-root coordinates

This section identifies the gradient-flow structure behind the moment ODE. We first compute its Stein mobility and then use square-root coordinates to obtain strong convexity, compact sublevels, and quantitative convergence.

All matrix tangent and cotangent variables below lie in $\Sym(d)$ , equipped with the Frobenius inner product. The product inner product is

\left\langle(a,B),(u,H)\right\rangle=a^{\top}u+\tr(BH).

For $X\sim\mathcal{N}(\mu,\Sigma)$ , the exact reverse-KL energy is

\mathcal{E}(\mu,\Sigma)=\KL(\mathcal{N}(\mu,\Sigma)\|\pi)=\mathbb{E}V(X)-\frac{1}{2}\log\det\Sigma+c_{\pi}.

(4.1)

Lemma 4.1 (Energy gradient and Stein mobility).

For $m$ and $\Gamma$ from (2.2),

\nabla_{\mu}\mathcal{E}=m,\qquad\nabla_{\Sigma}\mathcal{E}=\frac{1}{2}(\Gamma-\Sigma^{-1}).

(4.2)

For $(a,B)\in\mathbb{R}^{d}\times\Sym(d)$ , define


$\displaystyle[\mathcal{M}_{\mu,\Sigma}(a,B)]_{\mu}$	$\displaystyle=2B\Sigma\mu+(1+\left\lvert\mu\right\rvert^{2})a,$	(4.3a)
$\displaystyle[\mathcal{M}_{\mu,\Sigma}(a,B)]_{\Sigma}$	$\displaystyle=\Sigma(2\Sigma B+\mu a^{\top})+(2B\Sigma+a\mu^{\top})\Sigma.$	(4.3b)

Then

\left\langle(a,B),\mathcal{M}_{\mu,\Sigma}(a,B)\right\rangle=\left\lVert 2B\Sigma+a\mu^{\top}\right\rVert_{\mathrm{F}}^{2}+\left\lvert a\right\rvert^{2}.

(4.4)

In particular, $\mathcal{M}_{\mu,\Sigma}$ is symmetric positive definite, and (2.3) is precisely

(\dot{\mu},\dot{\Sigma})=-\mathcal{M}_{\mu,\Sigma}\nabla\mathcal{E}(\mu,\Sigma).

(4.5)

Proof.

For $H\in\Sym(d)$ and all sufficiently small $\varepsilon$ , Gaussian differentiation (Price’s identity) gives

D_{\mu}\mathbb{E}V(X)[u]=u^{\top}m,\qquad\left.\frac{d}{d\varepsilon}\mathbb{E}_{\mathcal{N}(\mu,\Sigma+\varepsilon H)}V(X)\right|_{\varepsilon=0}=\frac{1}{2}\tr(\Gamma H).

The growth needed to differentiate under the integral follows from (2.1); differentiation of the log determinant then proves (4.2).

More generally, direct expansion gives the bilinear identity

\begin{split}\left\langle(a,B),\mathcal{M}_{\mu,\Sigma}(a^{\prime},B^{\prime})\right\rangle&=\left\langle 2B\Sigma+a\mu^{\top},2B^{\prime}\Sigma+a^{\prime}\mu^{\top}\right\rangle_{\mathrm{F}}+a^{\top}a^{\prime}.\end{split}

(4.6)

The right-hand side is unchanged when the primed and unprimed variables are interchanged, so $\mathcal{M}_{\mu,\Sigma}$ is self-adjoint. Taking the two arguments equal gives

\left\lvert a\right\rvert^{2}+\left\lvert\mu\right\rvert^{2}\left\lvert a\right\rvert^{2}+4a^{\top}B\Sigma\mu+4\tr(B\Sigma^{2}B),

which is (4.4); it vanishes only when $a=0$ and $B=0$ because $\Sigma\succ 0$ . Substituting $a=m$ and $B=(\Gamma-\Sigma^{-1})/2$ in $-\mathcal{M}(a,B)$ yields (2.3a) and (2.3b). ∎

Recall that $q=(\mu,A)$ with $A=\Sigma^{1/2}$ and that $F$ is defined by (2.7). Thus $\mathcal{E}(\mu,A^{2})=F(\mu,A)+c_{\pi}$ , so the two energies have identical derivatives in these coordinates.

The next lemma supplies the strong convexity and coercivity used throughout the stability analysis.

Lemma 4.2 (Regularity, strong convexity, and coercivity).

The function $F$ is $C^{2}$ on $\mathcal{Q}$ and, for $(u,H)\in\mathbb{R}^{d}\times\Sym(d)$ ,

$\displaystyle D^{2}F(\mu,A)[(u,H),(u,H)]$	$\displaystyle=\mathbb{E}\!\left[(u+HZ)^{\top}\nabla^{2}V(\mu+AZ)(u+HZ)\right]$
	$\displaystyle\quad+\tr(A^{-1}HA^{-1}H)$
	$\displaystyle\geq\alpha\bigl(\left\lvert u\right\rvert^{2}+\left\lVert H\right\rVert_{\mathrm{F}}^{2}\bigr).$	(4.7)

Moreover, $F$ has compact sublevel sets in $\mathcal{Q}$ and hence a unique minimizer $q_{*}=(\mu_{*},A_{*})$ .

Proof.

The Hessian bound implies $\left\lvert\nabla V(x)\right\rvert\leq C+\beta\left\lvert x\right\rvert$ and $\left\lvert V(x)\right\rvert\leq C(1+\left\lvert x\right\rvert^{2})$ . The first two derivatives of $V(\mu+AZ)$ are therefore dominated locally in $(\mu,A)$ by integrable polynomials in $Z$ . In particular,

DF(\mu,A)[u,H]=\mathbb{E}[\nabla V(\mu+AZ)^{\top}(u+HZ)]-\tr(A^{-1}H).

On a compact parameter neighborhood the integrands in the first and second variations are bounded respectively by $C(1+\left\lvert Z\right\rvert)^{2}$ and $C(1+\left\lvert Z\right\rvert)^{2}$ . Since $\nabla^{2}V$ is bounded and continuous, dominated convergence gives $F\in C^{2}$ and the displayed second-variation formula. The first term in (4.7) is at least

\alpha\mathbb{E}\left\lvert u+HZ\right\rvert^{2}=\alpha(\left\lvert u\right\rvert^{2}+\left\lVert H\right\rVert_{\mathrm{F}}^{2}),

and the log-determinant term is nonnegative because $\tr(A^{-1}HA^{-1}H)=\left\lVert A^{-1/2}HA^{-1/2}\right\rVert_{\mathrm{F}}^{2}$ .

Let $x_{V}$ be the unique minimizer of $V$ and set $c_{V}=V(x_{V})$ . Strong convexity gives, with $a_{1},\ldots,a_{d}$ the eigenvalues of $A$ ,

F(\mu,A)\geq c_{V}+\frac{\alpha}{2}\left\lvert\mu-x_{V}\right\rvert^{2}+\sum_{j=1}^{d}\left(\frac{\alpha}{2}a_{j}^{2}-\log a_{j}\right).

(4.8)

The right-hand side diverges if $\left\lvert\mu\right\rvert\to\infty$ , if an eigenvalue of $A$ tends to zero, or if $\left\lVert A\right\rVert_{\mathrm{op}}\to\infty$ . Thus every sublevel is compactly contained in $\mathcal{Q}$ . Existence follows by compactness and uniqueness by (4.7). Under the bijection $\Sigma=A^{2}$ , this minimizer corresponds to the unique reverse-KL optimal Gaussian $\rho_{*}=\mathcal{N}(\mu_{*},\Sigma_{*})$ , with $\Sigma_{*}=A_{*}^{2}$ . ∎

Let $\Phi(\mu,A)=(\mu,A^{2})$ . Its differential is

J_{q}(u,H)=(u,AH+HA).

(4.9)

The Sylvester map $H\mapsto AH+HA$ is invertible on $\Sym(d)$ when $A\succ 0$ . The chain rule and Lemma 4.1 therefore transform the moment flow into

\dot{q}=f(q)=-P(q)\nabla F(q),\qquad P(q)=J_{q}^{-1}\mathcal{M}_{\Phi(q)}J_{q}^{-\mathsf{T}}.

(4.10)

Here $P$ is smooth, symmetric, and positive definite. Lemma 4.2 gives $\nabla F\in C^{1}$ , and therefore $f=-P\nabla F\in C^{1}(\mathcal{Q})$ . This is the only local ODE regularity used below.

Proposition 4.3 (Sublevel ellipticity and convergence).

For $R\geq 0$ , recall

\mathcal{K}_{R}=\{q\in\mathcal{Q}:F(q)-F(q_{*})\leq R\}.

There are $M_{R},a_{R},b_{R}>0$ such that on $\mathcal{K}_{R}$ ,

\left\lvert\mu\right\rvert\leq M_{R},\qquad a_{R}I\preceq A\preceq b_{R}I.

(4.11)

For $R>0$ , one may take

P(q)\succeq p_{R}I,\qquad p_{R}=\frac{\min(1,4a_{R}^{4})}{2(1+M_{R}^{2})\max(1,4b_{R}^{2})}>0.

(4.12)

If $R>0$ , every solution of (4.10) starting in $\mathcal{K}_{R}$ is global, remains in $\mathcal{K}_{R}$ , and satisfies

	$\displaystyle F(q_{t})-F(q_{*})$	$\displaystyle\leq e^{-2\alpha p_{R}t}\{F(q_{0})-F(q_{*})\},$		(4.13)
	$\displaystyle\left\lVert q_{t}-q_{*}\right\rVert$	$\displaystyle\leq\sqrt{2R/\alpha}\,e^{-\alpha p_{R}t}.$		(4.14)

For $R=0$ , the only trajectory is $q_{t}=q_{*}$ , which is global and stationary.

Proof.

For $R=0$ , uniqueness gives $\mathcal{K}_{0}=\{q_{*}\}$ , so all claims follow with the stationary trajectory $q_{t}=q_{*}$ . Assume henceforth that $R>0$ . The bounds (4.11) follow from Lemma 4.2. To prove (4.12), put $s=\lambda_{\min}(\Sigma)$ and $U=2B\Sigma$ . Then

\left\lVert U\right\rVert_{\mathrm{F}}^{2}+\left\lvert a\right\rvert^{2}\leq 2(1+\left\lvert\mu\right\rvert^{2})\{\left\lVert U+a\mu^{\top}\right\rVert_{\mathrm{F}}^{2}+\left\lvert a\right\rvert^{2}\},

whereas

\left\lVert U\right\rVert_{\mathrm{F}}^{2}+\left\lvert a\right\rvert^{2}\geq\min(1,4s^{2})(\left\lVert B\right\rVert_{\mathrm{F}}^{2}+\left\lvert a\right\rvert^{2}).

Thus

\mathcal{M}_{\mu,\Sigma}\succeq\frac{\min(1,4s^{2})}{2(1+\left\lvert\mu\right\rvert^{2})}I.

Since $s\geq a_{R}^{2}$ and $\left\lVert J_{q}\right\rVert_{\mathrm{op}}^{2}\leq\max(1,4b_{R}^{2})$ on $\mathcal{K}_{R}$ , the congruence $P(q)=J_{q}^{-1}\mathcal{M}_{\Phi(q)}J_{q}^{-\mathsf{T}}$ gives (4.12): indeed, $\left\lVert J_{q}^{-\mathsf{T}}z\right\rVert\geq\left\lVert z\right\rVert/\left\lVert J_{q}\right\rVert_{\mathrm{op}}$ .

Strong convexity and the fact that the segment from $q$ to the interior minimizer $q_{*}$ remains in $\mathcal{Q}$ imply the Polyak–Lojasiewicz bound

\left\lVert\nabla F(q)\right\rVert^{2}\geq 2\alpha\{F(q)-F(q_{*})\}.

(4.15)

Along the flow,

\frac{d}{dt}\{F(q_{t})-F(q_{*})\}=-\nabla F(q_{t})^{\top}P(q_{t})\nabla F(q_{t})\leq-2\alpha p_{R}\{F(q_{t})-F(q_{*})\}.

This proves invariance and (4.13). Strong convexity also gives $F(q)-F(q_{*})\geq(\alpha/2)\left\lVert q-q_{*}\right\rVert^{2}$ , proving (4.14). Finally, a local solution trapped in the compact set $\mathcal{K}_{R}\Subset\mathcal{Q}$ cannot explode or reach the positive-definite boundary, so the standard continuation theorem gives global existence. ∎

5 Semiglobal all-time stability

The next lemma is the deterministic core of the argument. Its constant is uniform in time but is allowed to depend on the initial energy sublevel.

Lemma 5.1 (Semiglobal uniform stability).

For every $0\leq R<\infty$ there is $S_{R}<\infty$ such that any two solutions of (4.10) with $q_{0},\widetilde{q}_{0}\in\mathcal{K}_{R}$ satisfy

\sup_{t\geq 0}\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert\leq S_{R}\left\lVert q_{0}-\widetilde{q}_{0}\right\rVert.

(5.1)

If $R>0$ , there also exist $T_{R},D_{R},\kappa_{R}>0$ , depending only on $R, d, V$ , such that

\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert\leq D_{R}e^{-\kappa_{R}(t-T_{R})}\left\lVert q_{0}-\widetilde{q}_{0}\right\rVert,\qquad t\geq T_{R}.

Proof.

The argument has two regimes. Sublevel invariance first gives a common compact set and hence an early-time Gronwall bound; exponential energy decay then brings both trajectories into a neighborhood of $q_{*}$ where a fixed $H_{*}$ -norm is contractive. Patching the estimates yields the all-time constant and eventual exponential decay.

The case $R=0$ is immediate. Suppose $R>0$ and set

H_{*}=D^{2}F(q_{*}),\qquad P_{*}=P(q_{*}),\qquad\kappa_{*}=\lambda_{\min}(P_{*})\lambda_{\min}(H_{*})>0.

Since $\nabla F(q_{*})=0$ ,

Df(q_{*})=-P_{*}H_{*}.

In the fixed norm $\left\lVert v\right\rVert_{H_{*}}^{2}=v^{\top}H_{*}v$ , write $\cond(H_{*})=\lambda_{\max}(H_{*})/\lambda_{\min}(H_{*})$ . Then

v^{\top}H_{*}Df(q_{*})v=-(H_{*}v)^{\top}P_{*}(H_{*}v)\leq-\kappa_{*}\left\lVert v\right\rVert_{H_{*}}^{2}.

(5.2)

Continuity of $D f$ gives $\varrho>0$ such that the convex ellipsoid

\mathcal{U}=\{q:\left\lVert q-q_{*}\right\rVert_{H_{*}}<\varrho\}\Subset\mathcal{Q}

satisfies, for all $q\in\mathcal{U}$ and all $v$ ,

v^{\top}H_{*}Df(q)v\leq-\frac{\kappa_{*}}{2}\left\lVert v\right\rVert_{H_{*}}^{2}.

(5.3)

Indeed, continuity is used here to control the symmetric part of $H_{*}Df$ :

\sup_{\left\lVert v\right\rVert_{H_{*}}=1}\left\lvert v^{\top}H_{*}\{Df(q)-Df(q_{*})\}v\right\rvert\longrightarrow 0\qquad(q\to q_{*}).

For $e=q-q_{*}$ , the segment formula and $f(q_{*})=0$ give

f(q)-f(q_{*})=\int_{0}^{1}Df(q_{*}+\theta e)e\,d\theta.

Consequently, whenever the segment lies in $\mathcal{U}$ ,

\frac{d}{dt}\left\lVert q_{t}-q_{*}\right\rVert_{H_{*}}^{2}=2\int_{0}^{1}e_{t}^{\top}H_{*}Df(q_{*}+\theta e_{t})e_{t}\,d\theta\leq-\kappa_{*}\left\lVert e_{t}\right\rVert_{H_{*}}^{2}.

(5.4)

Since every concentric ellipsoid is convex, the usual first-exit argument applied to (5.4) shows that each closed concentric ellipsoid contained in $\mathcal{U}$ is forward invariant.

If $q_{t}$ and $\widetilde{q}_{t}$ lie in the same such ellipsoid, then the segment joining them also lies there and

f(q_{t})-f(\widetilde{q}_{t})=\int_{0}^{1}Df(\widetilde{q}_{t}+\theta(q_{t}-\widetilde{q}_{t}))(q_{t}-\widetilde{q}_{t})\,d\theta.

Thus

\frac{d}{dt}\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert_{H_{*}}^{2}\leq-\kappa_{*}\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert_{H_{*}}^{2}.

(5.5)

By (4.14), all trajectories starting in $\mathcal{K}_{R}$ enter the ellipsoid of radius $\varrho/2$ no later than

T_{R}=\frac{1}{2\alpha p_{R}}\log_{+}\!\left(\frac{8\lambda_{\max}(H_{*})R}{\alpha\varrho^{2}}\right).

(5.6)

Before $T_{R}$ , all trajectories remain in the convex compact set

\mathcal{C}_{R}=\{(\mu,A):\left\lvert\mu\right\rvert\leq M_{R},\,a_{R}I\preceq A\preceq b_{R}I\}\Subset\mathcal{Q}.

Let $\ell_{R}=\sup_{q\in\mathcal{C}_{R}}\left\lVert Df(q)\right\rVert_{\mathrm{op}}<\infty$ , where the operator norm is induced by the Euclidean–Frobenius product norm. The segment between the two trajectories remains in $\mathcal{C}_{R}$ , so Gronwall’s inequality gives

\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert\leq e^{\ell_{R}t}\left\lVert q_{0}-\widetilde{q}_{0}\right\rVert,\qquad 0\leq t\leq T_{R}.

(5.7)

After $T_{R}$ , both trajectories remain in the half-radius ellipsoid. Combining (5.5), norm equivalence, and (5.7) gives

\left\lVert q_{t}-\widetilde{q}_{t}\right\rVert\leq\sqrt{\cond(H_{*})}\,e^{-\kappa_{*}(t-T_{R})/2}e^{\ell_{R}T_{R}}\left\lVert q_{0}-\widetilde{q}_{0}\right\rVert,\qquad t\geq T_{R}.

Together with (5.7), this proves (5.1) and the claimed eventual exponential decay. In particular, the construction gives the explicit admissible choice

S_{R}=\sqrt{\cond(H_{*})}\,e^{\ell_{R}T_{R}}.

(5.8)

∎

Remark 5.2 (Dependence of the stability constant).

Finiteness of $S_{R}$ requires no global bound on $\nabla^{3}V$ : continuity of $D^{2}F$ on the relevant compact parameter sets suffices. A closed-form upper bound for $S_{R}$ , however, requires a quantitative modulus of continuity for $\nabla^{2}V$ on those sets. A global bound on $\nabla^{3}V$ is one sufficient, but not necessary, condition for such a modulus.

The synchronous coupling $\mu+AZ$ and $\widetilde{\mu}+\widetilde{A}Z$ gives the useful consequence

W_{2}^{2}(\mathcal{N}(\mu,A^{2}),\mathcal{N}(\widetilde{\mu},\widetilde{A}^{2}))\leq\left\lvert\mu-\widetilde{\mu}\right\rvert^{2}+\left\lVert A-\widetilde{A}\right\rVert_{\mathrm{F}}^{2}.

(5.9)

6 Exact affine decomposition and shape error

We next derive the moment closure and affine representation directly from the particle system. The resulting affine lift separates the evolving moment error from a time-independent discrepancy of the standardized shape.

Lemma 6.1 (Exact closure and affine transport).

If $\widehat{\Sigma}_{0}\succ 0$ , the particle system has a unique global solution. Its empirical moments solve (2.3) from $(\widehat{\mu}_{0},\widehat{\Sigma}_{0})$ and remain positive definite. If $L_{t}$ solves

\dot{L}_{t}=C_{t}L_{t},\qquad L_{0}=I,

(6.1)

then

	$\displaystyle X_{i}(t)$	$\displaystyle=\widehat{\mu}_{t}+L_{t}(X_{i}(0)-\widehat{\mu}_{0}),$		(6.2)
	$\displaystyle\widehat{\Sigma}_{t}$	$\displaystyle=L_{t}\widehat{\Sigma}_{0}L_{t}^{\top}.$		(6.3)

Proof.

We first construct the particle solution, thereby avoiding any continuation argument that presupposes positive definiteness of its empirical covariance. Let

q_{t}^{N}=(\mu_{t}^{N},A_{t}^{N}),\qquad\Sigma_{t}^{N}=(A_{t}^{N})^{2},

be the global solution of (4.10) from $(\widehat{\mu}_{0},\widehat{\Sigma}_{0}^{1/2})$ , whose existence follows from Proposition 4.3. Along this trajectory define

m_{t}^{N}=m(\mu_{t}^{N},\Sigma_{t}^{N}),\quad\Gamma_{t}^{N}=\Gamma(\mu_{t}^{N},\Sigma_{t}^{N}),\quad C_{t}^{N}=I-\Gamma_{t}^{N}\Sigma_{t}^{N}-m_{t}^{N}(\mu_{t}^{N})^{\top},

and let $L_{t}$ be the global solution of $\dot{L}_{t}=C_{t}^{N}L_{t}$ , $L_{0}=I$ . Now set

X_{i}(t)=\mu_{t}^{N}+L_{t}\{X_{i}(0)-\widehat{\mu}_{0}\}.

(6.4)

Its empirical mean is $\mu_{t}^{N}$ . Its empirical covariance

\widetilde{\Sigma}_{t}=L_{t}\widehat{\Sigma}_{0}L_{t}^{\top}

solves

\dot{\widetilde{\Sigma}}_{t}=C_{t}^{N}\widetilde{\Sigma}_{t}+\widetilde{\Sigma}_{t}(C_{t}^{N})^{\top},\qquad\widetilde{\Sigma}_{0}=\widehat{\Sigma}_{0}.

The covariance equation in (2.3) can equivalently be written

\dot{\Sigma}_{t}^{N}=C_{t}^{N}\Sigma_{t}^{N}+\Sigma_{t}^{N}(C_{t}^{N})^{\top}.

Uniqueness for this linear matrix equation therefore gives $\widetilde{\Sigma}_{t}=\Sigma_{t}^{N}$ . Likewise, the mean equation is $\dot{\mu}_{t}^{N}=C_{t}^{N}\mu_{t}^{N}-m_{t}^{N}$ . Differentiating (6.4) now verifies

\dot{X}_{i}(t)=C_{t}^{N}X_{i}(t)-m_{t}^{N},

so the constructed paths solve (2.5) with their own empirical moments.

For completeness, any local particle solution with positive-definite empirical covariance satisfies, by averaging and centering,

\dot{\widehat{\mu}}=C_{t}\widehat{\mu}-\widehat{m},\qquad\dot{\widehat{\Sigma}}=C_{t}\widehat{\Sigma}+\widehat{\Sigma}C_{t}^{\top}.

These are exactly (2.3); uniqueness of the $q$ -flow forces its moments to equal $(\mu_{t}^{N},\Sigma_{t}^{N})$ . Its centered particles then solve the same linear equation $\dot{Y}_{i}=C_{t}^{N}Y_{i}$ with the same initial data, proving uniqueness. Finally,

\det L_{t}=\exp\!\left(\int_{0}^{t}\tr C_{s}^{N}\,ds\right)\neq 0

for every finite $t$ , so $\Sigma_{t}^{N}=L_{t}\widehat{\Sigma}_{0}L_{t}^{\top}\succ 0$ . Since both the global $q$ -flow and its linear lift exist on every finite interval, this proves all claims, with $(\mu_{t}^{N},\Sigma_{t}^{N})=(\widehat{\mu}_{t},\widehat{\Sigma}_{t})$ . ∎

Henceforth write

q_{N,t}=(\widehat{\mu}_{t},\widehat{\Sigma}_{t}^{1/2})

(6.5)

for this empirical-moment trajectory.

Write the initial particles as

X_{i}(0)=\mu_{0}+\Sigma_{0}^{1/2}\xi_{i},\qquad\xi_{i}\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,I_{d}),

(6.6)

and define

\bar{\xi}=\frac{1}{N}\sum_{i=1}^{N}\xi_{i},\qquad S=\frac{1}{N}\sum_{i=1}^{N}(\xi_{i}-\bar{\xi})(\xi_{i}-\bar{\xi})^{\top}.

(6.7)

For $N\geq d+1$ , $S\succ 0$ almost surely. Let

z_{i}=S^{-1/2}(\xi_{i}-\bar{\xi}),\qquad\eta_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{z_{i}},\qquad\gamma_{d}=\mathcal{N}(0,I_{d}).

(6.8)

This is precisely (2.9) with the whitening $R_{0}=S^{-1/2}\Sigma_{0}^{-1/2}$ .

The following estimate relates the standardized shape to an ordinary Gaussian empirical measure without using inverse sample-covariance moments.

Lemma 6.2 (Standardized shape error).

Let $N\geq d+1$ , let $\xi_{1},\ldots,\xi_{N}$ be i.i.d. standard Gaussians, and define $\bar{\xi},S,z_{i},\eta_{N}$ , and $\gamma_{d}$ by (6.7)–(6.8). Then

\frac{1}{N}\sum_{i=1}^{N}\left\lvert z_{i}-\xi_{i}\right\rvert^{2}=\left\lVert I-S^{1/2}\right\rVert_{\mathrm{F}}^{2}+\left\lvert\bar{\xi}\right\rvert^{2}\leq\left\lVert I-S\right\rVert_{\mathrm{F}}^{2}+\left\lvert\bar{\xi}\right\rvert^{2}.

(6.9)

Moreover, there is $C_{d}<\infty$ , depending only on $d$ , such that

\mathbb{E}W_{2}^{2}(\eta_{N},\gamma_{d})\leq C_{d}r_{d}(N).

(6.10)

Proof.

Since $N^{-1}\sum_{i}(\xi_{i}-\bar{\xi})=0$ , the cross term vanishes in

z_{i}-\xi_{i}=(S^{-1/2}-I)(\xi_{i}-\bar{\xi})-\bar{\xi}.

Using the definition of $S$ gives

\frac{1}{N}\sum_{i}\left\lvert(S^{-1/2}-I)(\xi_{i}-\bar{\xi})\right\rvert^{2}=\tr\{(S^{-1/2}-I)S(S^{-1/2}-I)\}=\left\lVert I-S^{1/2}\right\rVert_{\mathrm{F}}^{2}.

The scalar inequality $\left\lvert 1-\sqrt{\lambda}\right\rvert\leq\left\lvert 1-\lambda\right\rvert$ for $\lambda\geq 0$ proves (6.9).

Let $\nu_{N}=N^{-1}\sum_{i}\delta_{\xi_{i}}$ . Coupling by index and using the squared triangle inequality yields

\mathbb{E}W_{2}^{2}(\eta_{N},\gamma_{d})\leq 2\mathbb{E}\{\left\lVert I-S\right\rVert_{\mathrm{F}}^{2}+\left\lvert\bar{\xi}\right\rvert^{2}\}+2\mathbb{E}W_{2}^{2}(\nu_{N},\gamma_{d}).

Because $N S$ is Wishart with $N-1$ degrees of freedom, the first expectation is explicit:

\mathbb{E}\left\lvert\bar{\xi}\right\rvert^{2}=\frac{d}{N},\qquad\mathbb{E}\left\lVert S-I\right\rVert_{\mathrm{F}}^{2}=\frac{(N-1)d(d+1)+d}{N^{2}}.

(6.11)

For the remaining empirical-matching term, the $d=1$ estimate is the $p=2$ Gaussian bound of Bobkov and Ledoux [5, Corollary 6.14]; the $d=2$ estimate is due to Ledoux [8, Theorem 1]; and for $d\geq 3$ the claimed rate follows from Ledoux and Zhu [9, Theorem 1.1, with = p 2 ]:

\mathbb{E}W_{2}^{2}(\nu_{N},\gamma_{d})\leq C_{d}r_{d}(N).

The finitely many sample sizes excluded by an asymptotic formulation of these estimates are absorbed by enlarging $C_{d}$ . Since $N^{-1}\leq C_{d}r_{d}(N)$ , this proves (6.10). ∎

To transfer the standardized discrepancy through the particle dynamics, we record a two-sided bound for every invertible affine map.

Lemma 6.3 (Two-sided affine lifting).

For any full-rank initial cloud, use the whitening and standardized law in (2.9). Then, for every $t\geq 0$ ,

\begin{split}\lambda_{\min}(\widehat{\Sigma}_{t})W_{2}^{2}(\eta_{N},\gamma_{d})&\leq W_{2}^{2}\!\left(\zeta_{N,t},\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t})\right)\\ &\leq\lambda_{\max}(\widehat{\Sigma}_{t})W_{2}^{2}(\eta_{N},\gamma_{d}).\end{split}

(6.12)

Proof.

Set $B_{t}=L_{t}R_{0}^{-1}$ . Then $B_{t}B_{t}^{\top}=L_{t}\widehat{\Sigma}_{0}L_{t}^{\top}=\widehat{\Sigma}_{t}$ , and (6.2) gives

\zeta_{N,t}=(x\mapsto\widehat{\mu}_{t}+B_{t}x)_{\#}\eta_{N}.

(6.13)

The same affine map pushes $\gamma_{d}$ to $\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t})$ . Pushing forward any coupling of $\eta_{N}$ and $\gamma_{d}$ multiplies its quadratic cost by at most $\left\lVert B_{t}\right\rVert_{\mathrm{op}}^{2}=\lambda_{\max}(\widehat{\Sigma}_{t})$ , proving the upper bound. Conversely, pull any coupling of the two image laws back through the invertible affine map and use $\left\lvert B_{t}u\right\rvert\geq\sigma_{\min}(B_{t})\left\lvert u\right\rvert$ . Taking the infimum over image couplings proves the lower bound because $\sigma_{\min}(B_{t})^{2}=\lambda_{\min}(\widehat{\Sigma}_{t})$ . ∎

Proof of Theorem 2.1.

Lemma 6.1 gives a unique global particle solution and identifies its moments with the $q$ -flow from $q_{N,0}$ . Since both initial moment states lie in $\mathcal{K}_{R}$ , Proposition 4.3 and Lemma 5.1 give

\sup_{t\geq 0}W_{2}^{2}\!\left(\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t}),\rho_{t}\right)\leq S_{R}^{2}\left\lVert q_{N,0}-q_{0}\right\rVert^{2},

(6.14)

where we used the synchronous Gaussian coupling (5.9). Invariance of $\mathcal{K}_{R}$ gives $\lambda_{\max}(\widehat{\Sigma}_{t})\leq b_{R}^{2}$ . The upper bound in Lemma 6.3, followed by the squared triangle inequality, now proves (2.10). ∎

Proof of Proposition 2.6.

The proof first uses exponential moment convergence to show that the affine coefficient is integrable. Its fundamental matrix then converges to an invertible limit, which identifies both the limiting empirical law and its two-sided shape floor.

Write

C(q)=I-\Gamma(\mu,A^{2})A^{2}-m(\mu,A^{2})\mu^{\top}.

The coordinate map $\Phi$ has invertible differential. Since $q_{*}$ is an interior minimizer, $\nabla\mathcal{E}(\mu_{*},\Sigma_{*})=0$ ; hence

m(\mu_{*},\Sigma_{*})=0,\qquad\Gamma(\mu_{*},\Sigma_{*})=\Sigma_{*}^{-1},\qquad C(q_{*})=0.

(6.15)

The map $C$ is Lipschitz on each $\mathcal{K}_{R}$ . For $m$ , this follows by synchronously coupling the Gaussians and using the $\beta$ -Lipschitzness of $\nabla V$ . For $\Gamma$ , boundedness of $\nabla^{2}V$ , Pinsker’s inequality, and the Gaussian KL formula give, for $q,\widetilde{q}\in\mathcal{K}_{R}$ ,

\left\lVert\Gamma(\mu,A^{2})-\Gamma(\widetilde{\mu},\widetilde{A}^{2})\right\rVert_{\mathrm{F}}\leq 2\sqrt{d}\,\beta\,d_{\rm TV}\!\left(\mathcal{N}(\mu,A^{2}),\mathcal{N}(\widetilde{\mu},\widetilde{A}^{2})\right)\leq C_{R}\left\lVert q-\widetilde{q}\right\rVert.

(6.16)

The last inequality holds because the Gaussian KL divergence is bounded by $C_{R}\left\lVert q-\widetilde{q}\right\rVert^{2}$ on the compact covariance box containing $\mathcal{K}_{R}$ . Set $\lambda_{R}=\alpha p_{R}$ . Proposition 4.3 and (6.15) therefore imply

\left\lVert C(q_{N,t})\right\rVert_{\mathrm{op}}\leq C_{R}e^{-\lambda_{R}t},\qquad\int_{0}^{\infty}\left\lVert C(q_{N,t})\right\rVert_{\mathrm{op}}\,dt<\infty.

(6.17)

The fundamental matrix from (6.1) consequently satisfies

\sup_{t\geq 0}\bigl(\left\lVert L_{t}\right\rVert_{\mathrm{op}}+\left\lVert L_{t}^{-1}\right\rVert_{\mathrm{op}}\bigr)<\infty.

Indeed, this follows by Gronwall for $L_{t}$ and for $\frac{d}{dt}L_{t}^{-1}=-L_{t}^{-1}C(q_{N,t})$ . Moreover,

L_{\infty}=I+\int_{0}^{\infty}C(q_{N,s})L_{s}\,ds

exists, is invertible, and $\left\lVert L_{t}-L_{\infty}\right\rVert_{\mathrm{F}}\leq C_{R}e^{-\lambda_{R}t}$ . Invertibility is explicit from

\sigma_{\min}(L_{\infty})=\lim_{t\to\infty}\sigma_{\min}(L_{t})\geq\left(\sup_{t\geq 0}\left\lVert L_{t}^{-1}\right\rVert_{\mathrm{op}}\right)^{-1}>0.

Set $B_{t}=L_{t}R_{0}^{-1}$ and $B_{\infty}=L_{\infty}R_{0}^{-1}$ . Because $q_{N,0}\in\mathcal{K}_{R}$ ,

\left\lVert R_{0}^{-1}\right\rVert_{\mathrm{op}}^{2}=\lambda_{\max}(\widehat{\Sigma}_{0})\leq b_{R}^{2},

so $\left\lVert B_{t}-B_{\infty}\right\rVert_{\mathrm{F}}\leq C_{R}e^{-\lambda_{R}t}$ . Passing to the limit in $B_{t}B_{t}^{\top}=\widehat{\Sigma}_{t}$ gives $B_{\infty}B_{\infty}^{\top}=\Sigma_{*}$ . Since the standardized cloud has zero mean and identity covariance, coupling particles by index yields

W_{2}^{2}(\zeta_{N,t},\zeta_{N,\infty})\leq\left\lvert\widehat{\mu}_{t}-\mu_{*}\right\rvert^{2}+\frac{1}{N}\sum_{i=1}^{N}\left\lvert(B_{t}-B_{\infty})z_{i}\right\rvert^{2}=\left\lvert\widehat{\mu}_{t}-\mu_{*}\right\rvert^{2}+\left\lVert B_{t}-B_{\infty}\right\rVert_{\mathrm{F}}^{2},

which proves (2.18). Metric continuity and $\rho_{t}\to\rho_{*}$ prove (2.19); applying the two-sided lifting argument to $B_{\infty}$ proves (2.20). Finally, $\eta_{N}$ is finitely atomic and $\gamma_{d}$ is non-atomic, so their Wasserstein distance is strictly positive. Equivalently, $O_{\infty}=\Sigma_{*}^{-1/2}B_{\infty}$ is orthogonal: the dynamics optimizes the Gaussian moments but can only rotate, not Gaussianize, the standardized initial shape. If $R_{0}^{\prime}=OR_{0}$ is another whitening, with $O$ orthogonal, then $\eta_{N}^{\prime}=O_{\#}\eta_{N}$ and $B_{\infty}^{\prime}=B_{\infty}O^{\top}$ . Hence the pushforward in (2.17), and therefore the limiting empirical law, is independent of the chosen whitening coordinates. ∎

Proof of Corollary 2.7.

Let $P_{N}=N^{-1}\sum_{i}\delta_{X_{i}(0)}$ and use the principal empirical whitening in (2.9). The empirical Wasserstein law of large numbers gives $W_{2}(P_{N},P_{0})\to 0$ almost surely; the sample mean and covariance converge almost surely as well. Hence $\widehat{\Sigma}_{0}^{-1/2}\to\Sigma_{0}^{-1/2}$ . Coupling corresponding sample points and then applying the fixed limiting affine map shows

W_{2}(\eta_{N},\eta_{0})\longrightarrow 0\qquad\text{almost surely}.

On the same probability-one event, every $N\geq d+1$ cloud is full rank and its moment state has finite energy, hence belongs to some sublevel $\mathcal{K}_{R}$ . The two-sided comparison below is independent of that cloud-dependent sublevel. Explicitly, the squared cost of replacing the empirical centering and whitening by their population values tends to zero by the sample second-moment law of large numbers, while the remaining term is at most $\left\lVert\Sigma_{0}^{-1/2}\right\rVert_{\mathrm{op}}W_{2}(P_{N},P_{0})$ . Taking the lower limit in the lower bound and the upper limit in the upper bound of (2.20) proves (2.21). Here $D_{0}>0$ because $\eta_{0}=\gamma_{d}$ would imply that $P_{0}$ itself is Gaussian. ∎

7 Localization and proof of the main theorem

This section upgrades the deterministic oracle to the expected and high-probability theorems. A fixed well-conditioned event supplies a common energy sublevel, while a log-Wishart estimate controls its rare complement.

Fix $\varepsilon_{0}=1/4$ and define the good event

G_{N}=\{\left\lvert\bar{\xi}\right\rvert\leq\varepsilon_{0},\ \left\lVert S-I\right\rVert_{\mathrm{op}}\leq\varepsilon_{0}\}.

(7.1)

Lemma 7.1 (Good-event control).

There are constants $C_{d},c_{d}>0$ such that, for all $N\geq d+1$ ,

\mathbb{P}(G_{N}^{c})\leq C_{d}e^{-c_{d}N}.

(7.2)

There is a deterministic $R<\infty$ , depending only on the fixed problem data, such that on $G_{N}$ both $q_{0}=(\mu_{0},\Sigma_{0}^{1/2})$ and $q_{N,0}=(\widehat{\mu}_{0},\widehat{\Sigma}_{0}^{1/2})$ belong to $\mathcal{K}_{R}$ . Furthermore, there is $C<\infty$ , depending only on $d,\mu_{0},\Sigma_{0}$ , such that

\mathbb{E}\!\left[\mathbf{1}_{G_{N}}\left\lVert q_{N,0}-q_{0}\right\rVert^{2}\right]\leq\frac{C}{N}.

(7.3)

Proof.

Since $\sqrt{N}\,\bar{\xi}\sim\mathcal{N}(0,I_{d})$ , the Gaussian mean has an exponential tail at the fixed threshold $\varepsilon_{0}$ . Furthermore,

S-I=\frac{1}{N}\sum_{i=1}^{N}(\xi_{i}\xi_{i}^{\top}-I)-\bar{\xi}\bar{\xi}^{\top}.

Thus

	$\displaystyle\mathbb{P}(\left\lVert S-I\right\rVert_{\mathrm{op}}>\varepsilon_{0})$	$\displaystyle\leq\mathbb{P}\!\left(\left\\|\frac{1}{N}\sum_{i=1}^{N}(\xi_{i}\xi_{i}^{\top}-I)\right\\|_{\mathrm{op}}>\frac{\varepsilon_{0}}{2}\right)$
		$\displaystyle\quad+\mathbb{P}(\left\lvert\bar{\xi}\right\rvert^{2}>\varepsilon_{0}/2).$

For each fixed pair of coordinates $(j,k)$ , the centered variable $\xi_{ij}\xi_{ik}-\delta_{jk}$ has a subexponential norm bounded by a constant independent of $i$ . Scalar Bernstein’s inequality, a union bound over the $d^{2}$ entries, and fixed-dimensional norm equivalence bound the first probability by $C_{d}e^{-c_{d}N}$ ; a Gaussian tail gives the same bound for the second. Together with the mean event in (7.1), and after enlarging $C_{d}$ for the finitely many small sample sizes, this proves (7.2).

On $G_{N}$ ,

\widehat{\mu}_{0}=\mu_{0}+\Sigma_{0}^{1/2}\bar{\xi},\qquad\widehat{\Sigma}_{0}=\Sigma_{0}^{1/2}S\Sigma_{0}^{1/2},

and

(1-\varepsilon_{0})\lambda_{\min}(\Sigma_{0})I\preceq\widehat{\Sigma}_{0}\preceq(1+\varepsilon_{0})\lambda_{\max}(\Sigma_{0})I.

Thus $q_{N,0}$ ranges over a fixed compact subset of $\mathcal{Q}$ ; continuity of $F$ places it, together with $q_{0}$ , in a deterministic sublevel $\mathcal{K}_{R}$ .

For positive definite $C, D$ , the identity

C-D=C^{1/2}X+XD^{1/2},\qquad X=C^{1/2}-D^{1/2},

and the Frobenius pairing give

\left\langle C-D,X\right\rangle_{\mathrm{F}}\geq\bigl(\sqrt{\lambda_{\min}(C)}+\sqrt{\lambda_{\min}(D)}\bigr)\left\lVert X\right\rVert_{\mathrm{F}}^{2}.

Cauchy–Schwarz therefore implies

\left\lVert C^{1/2}-D^{1/2}\right\rVert_{\mathrm{F}}\leq\frac{\left\lVert C-D\right\rVert_{\mathrm{F}}}{\sqrt{\lambda_{\min}(C)}+\sqrt{\lambda_{\min}(D)}}.

(7.4)

On $G_{N}$ the denominator has a deterministic positive lower bound. Since

	$\displaystyle\left\lvert\widehat{\mu}_{0}-\mu_{0}\right\rvert$	$\displaystyle\leq\left\lVert\Sigma_{0}^{1/2}\right\rVert_{\mathrm{op}}\left\lvert\bar{\xi}\right\rvert,$
	$\displaystyle\left\lVert\widehat{\Sigma}_{0}-\Sigma_{0}\right\rVert_{\mathrm{F}}$	$\displaystyle\leq\left\lVert\Sigma_{0}\right\rVert_{\mathrm{op}}\left\lVert S-I\right\rVert_{\mathrm{F}},$

the square-root bound (7.4) and the exact moments in (6.11) prove (7.3). ∎

Proof of Theorem 2.4.

We intersect concentration events for the empirical measure, mean, and covariance. These events place all empirical initial states in one deterministic sublevel, so the oracle applies with constants independent of $N$ and $x$ .

Let $\nu_{N}=N^{-1}\sum_{i}\delta_{\xi_{i}}$ . Matching corresponding atoms shows that the map

(\xi_{1},\ldots,\xi_{N})\longmapsto W_{2}(\nu_{N},\gamma_{d})

is $N^{-1/2}$ -Lipschitz on $\mathbb{R}^{dN}$ . Gaussian concentration and $(\mathbb{E}W_{2})^{2}\leq\mathbb{E}W_{2}^{2}\leq C_{d}r_{d}(N)$ therefore give

\mathbb{P}\!\left\{W_{2}^{2}(\nu_{N},\gamma_{d})>C\left(r_{d}(N)+\frac{x}{N}\right)\right\}\leq e^{-x},\qquad x\geq 1.

(7.5)

Indeed, outside an event of probability $e^{-x}$ , $W_{2}(\nu_{N},\gamma_{d})\leq\mathbb{E}W_{2}(\nu_{N},\gamma_{d})+\sqrt{2x/N}$ , and the squared triangle inequality gives (7.5).

Fixed-dimensional Gaussian concentration for $\bar{\xi}$ and scalar Bernstein concentration applied to every entry of $N^{-1}\sum_{i}(\xi_{i}\xi_{i}^{\top}-I)$ imply that, for $1\leq x\leq c_{0}N$ , outside an event of probability $Ce^{-cx}$ ,

\left\lvert\bar{\xi}\right\rvert^{2}+\left\lVert S-I\right\rVert_{\mathrm{F}}^{2}\leq C\frac{1+x}{N},\qquad\left\lVert S-I\right\rVert_{\mathrm{op}}\leq\frac{1}{2}.

(7.6)

Here $c_{0}>0$ is reduced, if needed, so that the second assertion follows from the usual Bernstein scale $C_{d}\{\sqrt{x/N}+x/N\}$ ; indeed, $x/N\leq c_{0}$ throughout the stated range. On the intersection of the good events underlying (7.5) and (7.6), the square-root estimate (7.4) gives

\left\lVert q_{N,0}-q_{0}\right\rVert^{2}\leq C\frac{1+x}{N},\qquad W_{2}^{2}(\eta_{N},\gamma_{d})\leq C\left(r_{d}(N)+\frac{x}{N}\right).

(7.7)

The second bound follows from the residual identity and coupling the atoms $z_{i}$ with $\xi_{i}$ . Moreover, all $q_{N,0}$ on these events, together with the fixed $q_{0}$ , lie in one deterministic energy sublevel $\mathcal{K}_{R_{\rm hp}}$ : the covariance eigenvalues stay in a fixed compact interval, and the sample means stay in a fixed ball because $x\leq c_{0}N$ . Here $R_{\rm hp}$ depends only on the fixed problem data and the chosen $c_{0}$ , not on $N$ or $x$ . Denote this intersection by $H_{N}(x)$ ; the preceding bounds give $\mathbb{P}(H_{N}(x)^{c})\leq Ce^{-cx}$ .

The deterministic oracle (2.10), the fact that $r_{d}(N)\geq N^{-1}$ , and a union bound prove (2.14). For the last iterate, Proposition 4.3 and the synchronous Gaussian coupling give the deterministic estimate

W_{2}^{2}(\rho_{t},\rho_{*})\leq\frac{2\{F(q_{0})-F(q_{*})\}}{\alpha}e^{-2\alpha p_{\overline{R}_{\rm mf}}t},\qquad\overline{R}_{\rm mf}=1+F(q_{0})-F(q_{*}).

(7.8)

Combining this with (2.14) by the squared triangle inequality proves (2.15). Finally take $x=c^{-1}\log(C/\delta)$ and rename the fixed constants to obtain the confidence-level formulation. ∎

On the good event, Lemmas 5.1 and 7.1, together with (5.9), imply

\mathbb{E}\!\left[\mathbf{1}_{G_{N}}\sup_{t\geq 0}W_{2}^{2}\!\left(\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t}),\rho_{t}\right)\right]\leq\frac{C}{N}.

(7.9)

Invariance of $\mathcal{K}_{R}$ also gives a deterministic upper bound for $\sup_{t}\lambda_{\max}(\widehat{\Sigma}_{t})$ on $G_{N}$ . Hence Lemmas 6.2 and 6.3 yield

\mathbb{E}\!\left[\mathbf{1}_{G_{N}}\sup_{t\geq 0}W_{2}^{2}\!\left(\zeta_{N,t},\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t})\right)\right]\leq Cr_{d}(N).

(7.10)

The complement of $G_{N}$ cannot simply be discarded: there the sample covariance can be arbitrarily close to singular. The following lemma is the required weighted tail estimate.

Lemma 7.2 (Weighted control outside the good event).

There are constants $C,c>0$ , depending only on $d,V,\mu_{0},\Sigma_{0}$ , such that

\mathbb{E}\!\left[\mathbf{1}_{G_{N}^{c}}\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\right]\leq Ce^{-cN},\qquad N\geq d+1.

(7.11)

Proof.

We first dominate the entire path by the initial energy, reduce that energy to polynomial moments and a log determinant, and then use a uniform log-Wishart second moment with Cauchy–Schwarz on $G_{N}^{c}$ .

For a probability law $\nu$ , write $M_{2}(\nu)=\int\left\lvert x\right\rvert^{2}\,d\nu(x)$ . The independent coupling through the origin gives

W_{2}^{2}(\nu,\widetilde{\nu})\leq 2M_{2}(\nu)+2M_{2}(\widetilde{\nu}).

(7.12)

For the particle cloud,

M_{2}(\zeta_{N,t})=\left\lvert\widehat{\mu}_{t}\right\rvert^{2}+\tr\widehat{\Sigma}_{t}=\left\lVert q_{N,t}\right\rVert^{2}.

Strong convexity of $F$ and energy decrease along the random moment flow give

\sup_{t\geq 0}M_{2}(\zeta_{N,t})\leq C\{1+F(q_{N,0})-F(q_{*})\}.

(7.13)

The upper Hessian bound on $V$ gives

F(q_{N,0})-F(q_{*})\leq C\{1+\left\lvert\widehat{\mu}_{0}\right\rvert^{2}+\tr\widehat{\Sigma}_{0}+(-\log\det\widehat{\Sigma}_{0})_{+}\}.

(7.14)

Indeed, Taylor’s theorem and the upper Hessian bound imply $V(x)\leq C(1+\left\lvert x\right\rvert^{2})$ , so

\mathbb{E}V(\widehat{\mu}_{0}+\widehat{\Sigma}_{0}^{1/2}Z)\leq C\{1+\left\lvert\widehat{\mu}_{0}\right\rvert^{2}+\tr\widehat{\Sigma}_{0}\};

the entropy term is $-\tfrac{1}{2}\log\det\widehat{\Sigma}_{0}$ , whose positive contribution is exactly one half of $(-\log\det\widehat{\Sigma}_{0})_{+}$ . Define

Z_{N}=1+\left\lvert\widehat{\mu}_{0}\right\rvert^{2}+\tr\widehat{\Sigma}_{0}+(-\log\det\widehat{\Sigma}_{0})_{+}.

(7.15)

All log-determinant quantities here and below are evaluated on the almost-sure full-rank event; on its null complement they inherit the fixed extension specified after (2.5).

We next prove

\sup_{N\geq d+1}\mathbb{E}Z_{N}^{2}<\infty.

(7.16)

The polynomial terms have uniform moments of every fixed order: the sample mean is Gaussian with uniformly bounded covariance, while $\tr\widehat{\Sigma}_{0}$ is bounded above by a fixed multiple of $N^{-1}\sum_{i}\left\lvert\xi_{i}\right\rvert^{2}$ . For the remaining term, Bartlett’s decomposition [1, Lemma 7.2.1] gives

NS\sim\operatorname{Wishart}_{d}(N-1,I),\qquad\det(NS)\overset{d}{=}\prod_{j=1}^{d}Y_{N-j},\quad Y_{k}\sim\chi_{k}^{2}\ \text{independently}.

(7.17)

For $s$ in a neighborhood of zero,

\mathbb{E}Y_{k}^{s}=2^{s}\frac{\Gamma(k/2+s)}{\Gamma(k/2)}.

Differentiating the cumulant generating function of $\log Y_{k}$ at zero gives the digamma and trigamma identities

\mathbb{E}\log Y_{k}=\psi(k/2)+\log 2,\qquad\operatorname{Var}(\log Y_{k})=\psi_{1}(k/2)

and hence the explicit formula

\mathbb{E}\left\lvert\log(Y_{k}/k)\right\rvert^{2}=\psi_{1}(k/2)+\{\psi(k/2)+\log 2-\log k\}^{2}.

Here $\psi_{1}(k/2)\leq\psi_{1}(1/2)$ , and the second term is uniformly bounded by the standard digamma asymptotic $\psi(x)=\log x+O(x^{-1})$ , together with the finitely many bounded small values of $k$ . Hence

\sup_{k\geq 1}\mathbb{E}\left\lvert\log(Y_{k}/k)\right\rvert^{2}<\infty.

(7.18)

Now

\log\det S=\sum_{j=1}^{d}\left\{\log\frac{Y_{N-j}}{N-j}+\log\frac{N-j}{N}\right\}.

For fixed $d$ and $N\geq d+1$ , the deterministic logarithms are uniformly bounded. Hence (7.18) and $(\sum_{j=1}^{d}a_{j})^{2}\leq d\sum_{j=1}^{d}a_{j}^{2}$ show that

\sup_{N\geq d+1}\mathbb{E}[(-\log\det S)_{+}^{2}]<\infty.

(7.19)

Since $\det\widehat{\Sigma}_{0}=\det(\Sigma_{0})\det S$ , the asserted uniform second-moment bound (7.16) follows.

Proposition 4.3 gives a deterministic uniform bound on $M_{2}(\rho_{t})$ . Combining (7.12)–(7.14) therefore gives the pathwise estimate

\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\leq CZ_{N}.

(7.20)

Cauchy–Schwarz, (7.16), and (7.2) now yield

\mathbb{E}\!\left[\mathbf{1}_{G_{N}^{c}}\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\right]\leq C\mathbb{P}(G_{N}^{c})^{1/2}(\mathbb{E}Z_{N}^{2})^{1/2}\leq Ce^{-cN/2}.

Renaming $c/2$ proves (7.11). ∎

The algebraic endpoint is independent of the dynamics and follows from an affine-span characterization.

Proposition 7.3 (Exact full-rank endpoint).

For any $N$ points in $\mathbb{R}^{d}$ , the centered empirical covariance has rank at most $N-1$ and is positive definite if and only if the points affinely span $\mathbb{R}^{d}$ . Consequently $N\geq d+1$ is necessary for the unregularized positive-definite dynamics, and it is sufficient almost surely for i.i.d. initialization from any distribution having a density.

Proof.

If $Y=[x_{1}-\bar{x},\ldots,x_{N}-\bar{x}]$ , then $Y\mathbf{1}_{N}=0$ and $\widehat{\Sigma}=N^{-1}YY^{\top}$ . Thus $\operatorname{rank}(\widehat{\Sigma})=\operatorname{rank}(Y)\leq N-1$ , and rank $d$ is equivalent to affine spanning. Under a density, affine dependence of any fixed $d+1$ samples is a Lebesgue-null event. ∎

Remark 7.4 (Logarithmic control at the endpoint).

At $N=d+1$ , the final Bartlett factor in (7.17) has one degree of freedom. The covariance is positive definite almost surely, but $\mathbb{E}\tr(S^{-1})=\infty$ , whereas $\mathbb{E}[(\log\chi_{1}^{2})^{2}]<\infty$ . This is why logarithmic-energy localization reaches the exact algebraic endpoint.

Proof of Theorem 2.2.

On $G_{N}$ , the squared triangle inequality with the intermediate Gaussian $\mathcal{N}(\widehat{\mu}_{t},\widehat{\Sigma}_{t})$ , followed by (7.9) and (7.10), gives

\mathbb{E}\!\left[\mathbf{1}_{G_{N}}\sup_{t\geq 0}W_{2}^{2}(\zeta_{N,t},\rho_{t})\right]\leq C\{N^{-1}+r_{d}(N)\}\leq Cr_{d}(N).

Lemma 7.2 controls $G_{N}^{c}$ . Since $e^{-cN}\leq Cr_{d}(N)$ for every fixed $d$ , the two event contributions prove (2.12). Almost-sure positive definiteness at time zero follows from Proposition 7.3; global existence follows from Lemma 6.1 and Proposition 4.3. ∎

Proof of Corollary 2.5.

Apply Proposition 4.3 to the deterministic initial condition. Since $q_{*}$ represents $\rho_{*}$ , (4.14) and (5.9) give

W_{2}^{2}(\rho_{t},\rho_{*})\leq Ce^{-2ct}

for a constant $c>0$ . The squared triangle inequality and Theorem 2.2 prove (2.16). The statistical-floor statement follows by balancing the two terms. ∎

8 Why a global Euclidean stability constant is impossible

The stability constant in Lemma 5.1 necessarily depends on an energy sublevel. The next proposition rules out a uniform replacement.

Proposition 8.1 (No global Euclidean Lipschitz constant for the scalar flow).

In dimension $d=1$ , there is a smooth, strongly convex potential with globally bounded Hessian for which no finite $S$ can satisfy

\sup_{t\geq 0}\left\lvert\phi_{t}(r_{0})-\phi_{t}(\widetilde{r}_{0})\right\rvert\leq S\left\lvert r_{0}-\widetilde{r}_{0}\right\rvert,\qquad r_{0},\widetilde{r}_{0}>0,

where $\phi_{t}(r_{0})$ is the standard-deviation component of the moment flow started from $(\mu_{0},\Sigma_{0})=(0,r_{0}^{2})$ .

Proof.

Fix $\kappa>a>0$ and $k>0$ , and take

V(x)=\frac{\kappa}{2}x^{2}+\frac{a}{k^{2}}\cos(kx).

(8.1)

Then $\kappa-a\leq V^{\prime\prime}\leq\kappa+a$ . The subspace $\mu=0$ is invariant. Writing $r=\sqrt{\Sigma}$ on this subspace gives

m=0,\qquad\Gamma(r)=\kappa-ae^{-k^{2}r^{2}/2},\qquad\dot{r}=g(r):=r\{1-r^{2}\Gamma(r)\}.

(8.2)

In fact,

\frac{d}{dr}\{r^{2}\Gamma(r)\}=2r\{\kappa-ae^{-k^{2}r^{2}/2}\}+ak^{2}r^{3}e^{-k^{2}r^{2}/2}>0,

because $\kappa>a$ . Thus $r^{2}\Gamma(r)$ is strictly increasing from zero to infinity and has a unique crossing $r_{*}>0$ of level one. In particular, $g>0$ on $(0,r_{*})$ . Fix $r_{c}\in(0,r_{*})$ . For $r_{0}\in(0,r_{c})$ , let $t_{\mathrm{hit}}(r_{0})$ be the time at which $\phi_{t}(r_{0})$ first reaches $r_{c}$ , and then freeze time at $t_{0}=t_{\mathrm{hit}}(r_{0})$ . This time is finite because $g$ is continuous and strictly positive on $[r_{0},r_{c}]$ . The variational identity for a scalar autonomous flow is

D_{r}\phi_{t}(r_{0})=\frac{g(\phi_{t}(r_{0}))}{g(r_{0})}.

Moreover, $g(r_{0})=r_{0}(1+o(1))$ as $r_{0}\downarrow 0$ . Therefore, with $t_{0}$ fixed while differentiating,

D_{r}\phi_{t_{0}}(r_{0})=\frac{g(r_{c})}{g(r_{0})}\sim\frac{g(r_{c})}{r_{0}}\longrightarrow\infty\qquad(r_{0}\downarrow 0).

(8.3)

Any common Lipschitz constant for all fixed-time maps $\phi_{t}$ would bound the left-hand side of (8.3), a contradiction. If the full moment flow had one initialization-uniform Euclidean Lipschitz constant, its restriction to the invariant subspace $\mu=0$ would have the same property, so the scalar contradiction also rules out that claim. Notice that we have not differentiated the composite $\phi_{t_{\mathrm{hit}}(r_{0})}(r_{0})=r_{c}$ , whose derivative is zero. ∎

This obstruction does not contradict Lemma 5.1: a fixed energy sublevel stays a positive distance away from $r=0$ . It explains why near-singular Wishart samples must be controlled probabilistically through energy rather than by a global derivative bound for the flow.

9 An explicit nonquadratic example

This section verifies the assumptions and identifies the variational optimum for a concrete nonquadratic target. It illustrates that the theorem is not confined to Gaussian potentials.

Example 9.1 (Cosine perturbation of a Gaussian target).

Let

V_{\varepsilon}(x)=\frac{1}{2}x^{2}+\varepsilon(1+\cos x),\qquad 0<\varepsilon\leq 0.1.

(9.1)

Then $1-\varepsilon\leq V_{\varepsilon}^{\prime\prime}(x)\leq 1+\varepsilon$ , so Theorem 2.2 applies. If $s>0$ denotes the variance of $\mathcal{N}(\mu,s)$ , then

m(\mu,s)=\mu-\varepsilon e^{-s/2}\sin\mu,\qquad\Gamma(\mu,s)=1-\varepsilon e^{-s/2}\cos\mu.

(9.2)

The moment flow is

\dot{\mu}=(1-\Gamma s)\mu-(1+\mu^{2})m,\qquad\dot{s}=2s(1-\Gamma s-\mu m).

(9.3)

Up to an additive constant, its reverse-KL energy is

E_{\varepsilon}(\mu,s)=\frac{1}{2}(\mu^{2}+s-\log s)+\varepsilon e^{-s/2}\cos\mu,

(9.4)

and

\begin{pmatrix}\dot{\mu}\\ \dot{s}\end{pmatrix}=-\begin{pmatrix}1+\mu^{2}&2s\mu\\ 2s\mu&4s^{2}\end{pmatrix}\nabla E_{\varepsilon}(\mu,s).

(9.5)

The mobility determinant is $4s^{2}>0$ . With $a_{0}=e^{-s/2}$ ,

\nabla^{2}E_{\varepsilon}=\begin{pmatrix}1-\varepsilon a_{0}\cos\mu&(\varepsilon a_{0}/2)\sin\mu\\ (\varepsilon a_{0}/2)\sin\mu&1/(2s^{2})+(\varepsilon a_{0}/4)\cos\mu\end{pmatrix}.

(9.6)

The upper-left entry is at least $1-\varepsilon$ , and its Schur complement is bounded below as follows:

		$\displaystyle\frac{1}{2s^{2}}+\frac{\varepsilon a_{0}}{4}\cos\mu-\frac{(\varepsilon a_{0}\sin\mu/2)^{2}}{1-\varepsilon a_{0}\cos\mu}$		(9.7)
		$\displaystyle\geq\frac{1}{2s^{2}}-\frac{\varepsilon e^{-s/2}}{4}-\frac{\varepsilon^{2}e^{-s}}{4(1-\varepsilon)}$
		$\displaystyle\geq\frac{1}{2s^{2}}\left(1-8e^{-2}\varepsilon-\frac{2e^{-2}\varepsilon^{2}}{1-\varepsilon}\right)>0.$

Indeed, $s^{2}e^{-s/2}\leq 16e^{-2}$ and $s^{2}e^{-s}\leq 4e^{-2}$ . Hence $E_{\varepsilon}$ is strictly convex and coercive, since the perturbation in (9.4) is bounded while $\mu^{2}+s-\log s$ diverges at every boundary of $\mathbb{R}\times(0,\infty)$ . Furthermore,

\partial_{\mu}E_{\varepsilon}(0,s)=0,\qquad\partial_{s}E_{\varepsilon}(0,s)=\frac{1}{2}\left(1-\frac{1}{s}-\varepsilon e^{-s/2}\right).

Let $h(s)=s(1-\varepsilon e^{-s/2})$ . Then

h^{\prime}(s)=1-\varepsilon e^{-s/2}+\frac{\varepsilon s}{2}e^{-s/2}>1-\varepsilon>0,\qquad h(0+)=0,\qquad h(\infty)=\infty.

Thus there is a unique $s_{*}>0$ with $h(s_{*})=1$ . The point $(0,s_{*})$ is stationary and, by strict convexity, is the unique minimizer; equivalently,

s_{*}\{1-\varepsilon e^{-s_{*}/2}\}=1.

(9.8)

This gives a fully worked, genuinely nonquadratic target with an explicit scalar characterization of the optimum.

10 Related work and positioning

This section distinguishes the projected exact-moment problem from ordinary SVGD and from time-averaged propagation of chaos. It also identifies the empirical-matching results used only as statistical inputs.

SVGD was introduced by Liu and Wang [11], and its gradient-flow and mean-field viewpoints were developed by Liu [10], Lu et al. [13], and Duncan et al. [6]. Nonasymptotic and finite-particle analyses include Korba et al. [7] and Shi and Mackey [14]. These works concern ordinary SVGD and do not themselves give uniform physical-time $W_{2}$ stability for the projected system (2.5).

The closest starting point is Liu et al. [12], whose exact-moment closure and affine representation lead to the prior conjecture discussed above. Banerjee et al. [3] prove finite-particle KSD and $W_{2}$ bounds for ordinary SVGD, including bilinear-plus-Matérn kernels, with long-time statements for time-averaged particle laws. The broad KSD and $W_{1}$ conclusions of Balasubramanian et al. [2], and their in-probability $W_{2}$ conclusion, concern time-averaged laws uniformly over deterministic averaging horizons, rather than physical-time last iterates; their finite-feature results have a different scope. Banerjee and Kim [4] obtain physical-time results for unprojected, Langevin-regularized SVGD with Brownian noise, whereas our system is deterministic and uses an exact Gaussian-score projection. Stein and Li [15] analyze generalized bilinear kernels for Gaussian targets. None of these results directly yields the fixed- $K_{1}$ , non-Gaussian, projected full-law theorem here.

The dimension-dependent statistical input comes from Gaussian empirical matching [5, 8, 16, 9]. Our contribution is not a new matching estimate: it shows that the nonlinear all-time dynamics adds no particle-number penalty to the sharp expected time-zero scale.

11 Implications and limitations

The deterministic oracle applies to each fixed full-rank cloud, with constants determined by a sublevel containing both initial moment states. A sequence of deterministic, dependent, or variance-reduced designs inherits a uniform guarantee when it lies in one common sublevel and its moment and standardized-shape errors vanish. Gaussian independence is used only to obtain the sharp expected $N$ -rate and the stated high-probability upper bound.

At fixed $N$ , the affine dynamics converges to a discrete law with moments $(\mu_{*},\Sigma_{*})$ but preserves standardized shape up to rotation. Thus moment convergence must not be confused with Gaussianization. The nonquadratic target in Section 9 shows that the theory is not limited to Gaussian potentials.

The restrictions are substantive. Constants depend on $d,V,\mu_{0},\Sigma_{0}$ and deteriorate near the positive-definite boundary. Section 8 proves, for a smooth scalar potential with bounded Hessian, that the Euclidean standard-deviation flow has no single initialization-uniform Lipschitz constant over $(0,\infty)$ . This does not exclude other metrics or weighted stability mechanisms. No dimension-free, Euler-discretization, or estimated-moment claim is made. Finally, $\rho_{*}$ is the reverse-KL Gaussian approximation, not the non-Gaussian target itself, so approximation bias is separate from particle error.

12 Conclusion

Square-root geometry and affine whitening turn the non-Gaussian projected Gaussian-SVGD problem into a stable moment flow plus a preserved shape term. This yields a deterministic all-time oracle, a sharp expected $N$ -rate, a high-probability upper bound at the corresponding scale, validity at the minimal count $N=d+1$ , and an exact description of the finite-particle limit. The result resolves the strengthened continuous exact-moment conjecture under global strong convexity and a bounded continuous Hessian, without a third-derivative or Hessian-Lipschitz assumption.

References

[1] T. W. Anderson. An Introduction to Multivariate Statistical Analysis. Third edition, Wiley-Interscience, Hoboken, NJ, 2003.
[2] K. Balasubramanian, S. Banerjee, and A. Korba. Uniform-in-time propagation-of-chaos for Stein variational gradient descent. arXiv:2607.00149, 2026.
[3] S. Banerjee, K. Balasubramanian, and P. Ghosal. Improved finite-particle convergence rates for Stein variational gradient descent. In International Conference on Learning Representations, 2025.
[4] S. Banerjee and D. Kim. Quantitative target convergence and uniform-in-time propagation of chaos for Langevin-regularized SVGD. arXiv:2608.28827, 2026.
[5] S. Bobkov and M. Ledoux. One-dimensional empirical measures, order statistics, and Kantorovich transport distances. Memoirs of the American Mathematical Society, 261(1259), 2019. doi:10.1090/memo/1259.
[6] A. Duncan, N. Nüsken, and L. Szpruch. On the geometry of Stein variational gradient descent. Journal of Machine Learning Research, 24(56):1–39, 2023.
[7] A. Korba, A. Salim, M. Arbel, G. Luise, and A. Gretton. A non-asymptotic analysis for Stein variational gradient descent. In Advances in Neural Information Processing Systems, volume 33, pages 4672–4682, 2020.
[8] M. Ledoux. On optimal matching of Gaussian samples. Journal of Mathematical Sciences, 238(4):495–522, 2019. doi:10.1007/s10958-019-04253-6.
[9] M. Ledoux and J.-X. Zhu. On optimal matching of Gaussian samples III. Probability and Mathematical Statistics, 41(2):237–265, 2021. doi:10.37190/0208-4147.41.2.3.
[10] Q. Liu. Stein variational gradient descent as gradient flow. In Advances in Neural Information Processing Systems, volume 30, pages 3115–3123, 2017.
[11] Q. Liu and D. Wang. Stein variational gradient descent: A general purpose Bayesian inference algorithm. In Advances in Neural Information Processing Systems, volume 29, pages 2378–2386, 2016.
[12] T. Liu, P. Ghosal, K. Balasubramanian, and N. S. Pillai. Towards understanding the dynamics of Gaussian–Stein variational gradient descent. In Advances in Neural Information Processing Systems, volume 36, pages 61234–61291, 2023. doi:10.52202/075280-2676; arXiv:2305.14076.
[13] J. Lu, Y. Lu, and J. Nolen. Scaling limit of the Stein variational gradient descent: The mean field regime. SIAM Journal on Mathematical Analysis, 51(2):648–671, 2019. doi:10.1137/18M1187611.
[14] J. Shi and L. W. Mackey. A finite-particle convergence rate for Stein variational gradient descent. In Advances in Neural Information Processing Systems, volume 36, pages 26831–26844, 2023. doi:10.52202/075280-1166.
[15] V. Stein and W. Li. Towards understanding accelerated Stein variational gradient flow: Analysis of generalized bilinear kernels for Gaussian target distributions. arXiv:2509.04008, version 2, 2025.
[16] M. Talagrand. Scaling and non-standard matching theorems. Comptes Rendus Mathématique, 356(6):692–695, 2018. doi:10.1016/j.crma.2018.04.018.

Uniform-in-Time Particle Convergence for Projected Gaussian SVGD

Abstract

1 Introduction and prior conjecture

Open problem 1.1 (Quantitative exact-moment conjecture).

2 Model and main results

Theorem 2.1 (Deterministic localized oracle).

Theorem 2.2 (Uniform-in-time particle convergence).

Proposition 2.3 (Sharpness in particle number).

Proof.

Theorem 2.4 (High-probability all-time convergence).

Corollary 2.5 (Physical-time last iterate).

Proposition 2.6 (Exact finite-particle limit and shape floor).

Corollary 2.7 (Non-Gaussian initialization is not Gaussianized).

Remark 2.8 (Scope).

3 Proof overview

4 Reverse-KL geometry in square-root coordinates

Lemma 4.1 (Energy gradient and Stein mobility).

Proof.

Lemma 4.2 (Regularity, strong convexity, and coercivity).

Proof.

Proposition 4.3 (Sublevel ellipticity and convergence).

Proof.

5 Semiglobal all-time stability

Lemma 5.1 (Semiglobal uniform stability).

Proof.

Remark 5.2 (Dependence of the stability constant).

6 Exact affine decomposition and shape error

Lemma 6.1 (Exact closure and affine transport).

Proof.

Lemma 6.2 (Standardized shape error).

Proof.

Lemma 6.3 (Two-sided affine lifting).

Proof.

Proof of Theorem 2.1.

Proof of Proposition 2.6.

Proof of Corollary 2.7.

7 Localization and proof of the main theorem

Lemma 7.1 (Good-event control).

Proof.

Proof of Theorem 2.4.

Lemma 7.2 (Weighted control outside the good event).

Proof.

Proposition 7.3 (Exact full-rank endpoint).

Proof.

Remark 7.4 (Logarithmic control at the endpoint).

Proof of Theorem 2.2.

Proof of Corollary 2.5.

8 Why a global Euclidean stability constant is impossible

Proposition 8.1 (No global Euclidean Lipschitz constant for the scalar flow).

Proof.

9 An explicit nonquadratic example

Example 9.1 (Cosine perturbation of a Gaussian target).

10 Related work and positioning

11 Implications and limitations

12 Conclusion

References

Uniform-in-Time Particle Convergence
for Projected Gaussian SVGD