Proof. Proposition 6 states the sub-Gaussian bound for the forecast laws \(\tilde{\mu}_{XY}^{f,j}\), but its proof also controls the corresponding centered analysis variables \(\tilde{X}^{a,j}\) at each step of the same induction. Hence, for fixed observation paths \(y^{1:j}\), the analysis laws satisfy \(\left\|\tilde{\mu}^{a,j}_{X}-\mathbb{E}\tilde{\mu}^{a,j}_{X}\right\|_{2,\psi_2}\leq \tilde{C}_{\mathrm{subG}}\) for all \(1\le j\le J\), after increasing the constant if necessary, with a constant uniform over observation paths. Let us condition on \(Y^{1:J}\). Then \(\tilde{\mu}_{X,N}^{a,j}\) is an empirical measure of i.i.d.samples from the deterministic law \(\tilde{\mu}_X^{a,j}\), so Proposition 9 and a union bound over \(j=1,\ldots,J\) give \[\mathbb{P}\left(W_p(\tilde{\mu}_{X,N}^{a,j},\tilde{\mu}_X^{a,j})\le C\tilde{\gamma}_{k,N}^{p,n}\;\text{for all }1\le j\le J\,\middle|\,Y^{1:J}\right)\ge 1-\frac{1}{2k},\] after adjusting the constants. By the definition of \(\tilde{\gamma}_{k,N}^{p,n}\) and the elementary bound \(\max(a,b)\le a+b\), \[\tilde{\gamma}_{k,N}^{p,n}\le \gamma_{k,N}^{p,n}+\sqrt{\log k/N}.\] Since \(\sqrt{\log k/N}\lesssim_J (\log k)^{1+j/2}/\sqrt N\) for \(k\ge2\) and \(1\le j\le J\), the square-root term is absorbed by the interaction term. Integrating this conditional statement over \(Y^{1:J}\), intersecting it with the event from Theorem 1, and applying the triangle inequality yields the high-probability bound. The expectation bound follows similarly by conditioning on \(Y^{1:J}\), applying the \(k\)-free empirical-measure expectation rate \(\bar\gamma_N^{p,n}\) uniformly in the observation path, and combining with Corollary 1. The latter term is dominated by \(\bar\gamma_N^{p,n}\) up to constants for \(p\in[1,2]\). ◻

Proof. Pick \(K_5 =C_{\psi_2}\|X\|_{\psi_2}\) as in Proposition 4 such that Equation 27 holds. For any \(\lambda \in \mathbb{R}\), the moment generating function of \(\pi_{X|Y=y} \sim X^\prime\) is \(M_\lambda(y) \mathrel{\vcenter{:}}= \mathbb{E}\exp\left(\lambda X^\prime \right).\)

Bound mgf for a fixed value: Pick any \(|\lambda|\geq 1/K_5\). As a consequence of the disintegration theorem and sub-Gaussianity of \(X\), the following is true: \[\begin{align} \;\mathbb{E}M_\lambda(Y) &= \mathbb{E}\exp\left(\lambda X \right) \; \leq \;{\exp\left(K_5^2|\lambda|^2\right)}. \end{align}\] By Markov’s inequality, for all \(k\geq 1\), \(\mathbb{P}\big(M_\lambda(Y) \geq k\mathbb{E}(M_\lambda(Y)) \big) \leq k^{-1}.\) Since we have an upper bound on the expectation, we plug in \(k = \exp\left(cK_5^2|\lambda|^2\right)\) for some \(c \geq 0\) and conclude \[\mathbb{P}\left(M_\lambda(Y) \geq \exp\left(K_{5,c}^2|\lambda|^2\right) \right) \leq \exp\left( -cK_5^2|\lambda|^2\right)\] for any \(|\lambda| \geq 1/K_5\) and \(K_{5,c}\mathrel{\vcenter{:}}= \sqrt{1 + c}K_5\).

Bound on dyadic net: Define the set and the event \[\Lambda = \{\pm 2^{j}/K_{5} |j \in \mathbb{N}\}, \quad A_{c} = \left\{\exists \lambda \in \Lambda: M_\lambda(Y) \geq \exp\left(K_{5,c}^2|\lambda|^2\right) \right\}.\] Apply a union bound over the net, for any \(c>0\) \[\label{eq:prob95ac95bound} \mathbb{P}\left(A_{c} \right) \leq 2\sum\limits_{j = 0}^\infty \exp\left( - {c}2^{2j}\right) \;\leq\; 2\exp({-c})\left(1+\frac{1}{c\log 4}\right).\tag{28}\] The upper bound in the sum above is derived in Lemma 4 (Appendix 5.2).

Extending outside the dyadic net: Define the event \[B_{c} = \left\{ \pi_{X|Y=Y} \text{ is not sub-Gaussian with parameter } \tilde{K}_{5,c}\right\}\] for \(\tilde{K}_{5,c} = 2\sqrt{1+c}K_{5,c}\) and some choice of the regular conditional \(\pi_{X|Y}\). The parameter we refer to in the definition of \(B_c\) is the one as in Equation 27 of Proposition 4. The randomness in the event \(B_{c}\) comes from the randomness of the conditioning variable \(Y\). Therefore, equivalently, \[B_{c} = \left\{ \exists |\lambda| \geq 1/ \tilde{K}_{5, c} \text{ s.t. } M_{\lambda} (Y) > \exp\left({\tilde{K}_{5,c}^2|\lambda|^2}\right)\right\}.\] We will show that \(B_{c} \subseteq A_{c}\). Let \(|\lambda| \geq 1/ \tilde{K}_{5, c}\) be such that \[M_{\lambda} (Y) > \exp\left({\tilde{K}_{5,c}^2|\lambda|^2}\right).\] In what follows, we split into the cases (1) \(|\lambda| \geq 1/ K_{5}\), and (2) \(1/K_{5} \geq |\lambda| \geq 1/ \tilde{K}_{5, c}\).

• Case 1: Assume that \(|\lambda| \geq 1/ K_{5}\). By the mean value theorem, there are \(j \in \mathbb{N}, p \in [0, 1], s \in \{\pm 1\}\) such that we have the convex decomposition \[\lambda = p \cdot s2^{j}/K_5 + (1-p)\cdot s2^{j+1}/K_5.\] By log-convexity, this shows that \[\begin{align} K_{5, c}^2\left(\frac{2^{j+1}}{K_5}\right)^2& \leq \tilde{K}_{5,c}^2 |\lambda|^2 < \log M_{\lambda}(Y) \\ & \leq p \log M_{s2^{j}/K_5}(Y) + (1-p)\log M_{s2^{j+1}/K_5}(Y). \end{align}\] By the pigeonhole principle, this means that either \(\log M_{s2^{j}/K_5}(Y) > K_{5, c}^2\left(\frac{2^{j + 1}}{K_5}\right)^2\) or \(\log M_{s2^{j + 1}/K_5}(Y) > K_{5, c}^2\left(\frac{2^{j+1}}{K_5}\right)^2\), each of which implies that we are on \(A_{c}\).

• Case 2: Assume that \(1/K_{5} \geq |\lambda| \geq 1/ \tilde{K}_{5, c}\). Pick \(s= \text{sgn}(\lambda)\) and decompose convexly as \[\lambda = \left(1-|\lambda| K_5\right)\cdot 0 + \left(|\lambda| K_5\right)\cdot \frac{s}{K_5}.\] Log-convexity implies \[\begin{align} \tilde{K}_{5,c}^2 |\lambda|^2 < \log M_{\lambda}(Y) & \;\leq\; \left(1-|\lambda| K_5\right)\cdot\log M_0(Y) + \left(|\lambda| K_5\right) \cdot\log M_{\frac{s}{K_5}}(Y)\\ & \;=\; (|\lambda| K_5)\cdot\log M_{\frac{s}{K_5}}(Y). \end{align}\] This means, in particular, that \[\begin{align} 2(1 + c)= \frac{\tilde{K}_{5,c}}{K_5} \;\leq\; \frac{\tilde{K}_{5,c}^2 |\lambda|^2 }{|\lambda| K_5} \;<\; \log M_{\frac{s}{K_5}}(Y) \end{align}\] which is equivalent to \[M_{\frac{s}{K_5}}(Y)> \exp\left(2(1+c)\right) \geq \exp\left(K_{5,c}^2\left|\frac{s}{K_5}\right|^2\right)\] and so we are on the event \(A_{c}\).

This concludes the proof of \(B_{c} \subseteq A_{c}\). To conclude the entire argument, define any \(\Delta C>0\) sufficiently large, set \(C_{\mathrm{cond}}\mathrel{\vcenter{:}}= 2C_{\psi_2}^2(1+\Delta C)\), use equivalence of sub-Gaussian definitions (Proposition 4) and apply \(B_{c} \subseteq A_{c}\): \[\begin{align} \mathbb{P}\bigl(\|\pi_{X|Y=Y}\|_{\psi_2} > C_{\mathrm{cond}}(1 + c)\|X\|_{\psi_2} \bigr) & \leq \mathbb{P}\bigl(\|\pi_{X|Y=Y}\|_{\psi_2}\\ &\qquad > 2C_{\psi_2}^2(1 + c +\Delta C)\|X\|_{\psi_2} \bigr) \\ & \leq \mathbb{P}(B_{c + \Delta C}) \\ & \leq \mathbb{P}(A_{c + \Delta C}) \\ &\leq \exp(-c). \end{align}\] By inequality 28 , the last line holds for sufficiently large \(\Delta C>0\), independent of \(c\) and \(\|X\|_{\psi_2}\), completing the proof. ◻

Proof of Proposition 1. Let \(k\geq 2\). For \(j=1\), we interpret \(Y^{1:0}\) as the trivial \(\sigma\)-algebra, so that \(\mathbb{E}(Y^1\mid Y^{1:0})=\mathbb{E}Y^1\). Lemma 1 shows that the Euclidean norm of the centered process satisfies \[\left\|(X^j, Y^j) - \mathbb{E}(X^j, Y^j)\right\|_{2,\psi_2} \leq C_{\text{subG}}\] for \(1\leq j \leq J\).

Conditional expectation: In particular, \(\left\| Y^j - \mathbb{E}(Y^j)\right\|_2\) is sub-Gaussian with the same Orlicz-norm bound. Further, note that \[\|Y^j - \mathbb{E}(Y^j|Y^{1:j-1})\|_2 \leq \|Y^j - \mathbb{E}(Y^j)\|_2 + \|\mathbb{E}(Y^j|Y^{1:j-1})- \mathbb{E}(Y^j)\|_2.\] By Proposition 2, \(\|\mathbb{E}(Y^j|Y^{1:j-1})- \mathbb{E}(Y^j)\|_2\) is sub-Gaussian with the same Orlicz-norm bound. This shows \(\left\|Y^j - \mathbb{E}(Y^j|Y^{1:j-1})\right\|_{2,\psi_2} \leq 2 C_{\text{subG}}\) for all \(1 \leq j \leq J\). Appealing to the equivalent definitions of sub-Gaussianity in Proposition 4, we showed that \[\mathbb{P}\left( \|Y^j - \mathbb{E}(Y^j|Y^{1:j-1})\|_2 \geq c_1\sqrt{\log k}\right) \leq \frac{1}{3kJ}\] for a suitable choice of \(c_1 \mathrel{\vcenter{:}}= c_1(\sigma_X,\sigma_{\max}, \lambda_{\max}, L_\Psi, L_h, J)\) since \(C_{\text{subG}}\) is a constant depending only on \(\sigma_X,\sigma_{\max}, \lambda_{\max}, L_\Psi, L_h, J\) according to Lemma 1 .

Moments: Next, consider a version of the regular conditional \(\mu_{X Y}^{f, j}\) and put \[Z(y^{1:j-1}) \sim \mu_{X Y}^{f, j}(y^{1:j-1})\; \forall y^{1:j-1} \in (\mathbb{R}^m)^{j-1}\] where we make the dependence on the observations \(y^{1:j-1}\) explicit here. In particular, by disintegration \(Z(Y^{1:j-1})\) is the same in law as \((X^j,Y^j)\). By definition, for every \(q \geq 1\), \[\begin{align} \overline{M}_q\left( Z(y^{1:j-1})\right) & = \left(\mathbb{E}\left\|Z(y^{1:j-1})- \mathbb{E}\left(Z(y^{1:j-1}) \right) \right\|_2^q\right)^\frac{1}{q} \\ & \leq \underbrace{\left(\mathbb{E}\left\|Z(y^{1:j-1})- \mathbb{E}\left(X^j, Y^j \right) \right\|_2^q\right)^\frac{1}{q}}_{\text{I}}\\ &\quad +\underbrace{ \left\|\mathbb{E}(X^j, Y^j)- \mathbb{E}\left(Z(y^{1:j-1}) \right) \right\|_2}_{\text{II}}. \end{align}\] We proceed term-by-term.

Term I: By the disintegration theorem, we have \[\mathrm{Law}\left(\|Z(Y^{1:j-1})- \mathbb{E}(X^j, Y^j)\|_2\right) = \mathrm{Law}\left( \|(X^j, Y^j)- \mathbb{E}(X^j, Y^j)\|_2\right).\] Therefore, by Lemma 1, the scalar random variable \(\|Z(Y^{1:j-1})- \mathbb{E}(X^j, Y^j)\|_2\) is sub-Gaussian with Orlicz norm bounded by \(C_{\text{subG}}\). We now apply Proposition 5 to this scalar norm, with conditioning variable \(Y^{1:j-1}\). Thus, for all \(c \geq 0\) it holds that \[\left\|Z(y^{1:j-1})- \mathbb{E}(X^j, Y^j)\right\|_{2,\psi_2} \leq C_{\mathrm{cond}}(1+c) C_{\text{subG}}\] with probability \(1 - \exp\left( -c\right)\) in \(y^{1:j-1} = Y^{1:j-1}\). Using the moment definition of sub-Gaussianity in Proposition 4, we have that with the same probability and for all \(q\geq 1\), \[\left(\mathbb{E}\|Z(y^{1:j-1})- \mathbb{E}(X^j, Y^j)\|_2^q\right)^{\frac{1}{q}} \leq \sqrt{q} C_{\psi_2}C_{\mathrm{cond}}(1+c) C_{\text{subG}}.\] Reminding ourselves that \(k\geq 2\), put differently, there is a constant \[c_2 \mathrel{\vcenter{:}}= c_2(\sigma_X,\sigma_{\max}, \lambda_{\max}, L_\Psi, L_h, J)\] such that with probability \(1 - \frac{1}{3kJ}\) in \(y^{1:j-1}=Y^{1:j-1}\) \[\left(\mathbb{E}\|Z(y^{1:j-1})- \mathbb{E}(X^j, Y^j)\|_2^q\right)^{\frac{1}{q}} \leq c_2\sqrt{q} \log k.\] Term II: Further, by Proposition 2, \(\left\|\mathbb{E}(X^j,Y^j)- \mathbb{E}\left(Z(Y^{1:j-1})\right)\right\|_{2}\) is sub-Gaussian with norm \(C_{\text{subG}}\). Therefore, there is a constant \(c_3 = c_3(\sigma_X,\sigma_{\max}, \lambda_{\max}, L_\Psi, L_h, J)\) such that with probability \(1 - \frac{1}{3kJ}\) in \(y^{1:j-1} = Y^{1:j-1}\) it holds that \[\left\|\mathbb{E}(X^j, Y^j)- \mathbb{E}\left(Z(y^{1:j-1}) \right) \right\|_2 \leq c_3 \sqrt{\log k}.\] When \(j=1\), \(Y^{1:0}\) is trivial and \((\mathbb{R}^m)^0\) is a singleton; then \(Z(Y^{1:0})\sim\mathrm{Law}(X^1,Y^1)\), so term II is zero and the same estimates apply.

Union Bound: Taking a union bound, we showed that there is a constant \[C_M\mathrel{\vcenter{:}}= C_M(\sigma_X,\sigma_{\max}, \lambda_{\max}, L_\Psi, L_h, J, q)\] such that for all \(k \geq 2\) and \(1\leq q^\prime \leq q\), \[\begin{align} \mathbb{P}&\left( \overline{M}_{q^\prime}\left(\mu_{XY}^{f,j}\right) \leq C_M\log k \right.\\ \quad \text{and} \quad &\left.\|Y^j - \mathbb{E}(Y^j|Y^{1:j-1})\|_2^2 \leq C_M{\log k} \text{ for all }j \in \{1, \ldots, J\}\right) \geq 1 - \frac{1}{k}. \end{align}\] ◻

Proof. Inequality 30 does not depend on the particular coupling of \((X^j,Y^j)|Y^{1:j-1}\) and \((\tilde{X}_f^j,\tilde{Y}_f^j)|Y^{1:j-1}\), so we are free to choose whichever coupling we want, as long as the marginals are correct. The joint distribution of \(\{(X^j, Y^j)\}_{j\geq 1}\) is described through the dynamics in Equation 1 , namely \[\begin{align} X^{j}&=\Psi\left(X^{j-1}\right)+\xi^{j-1}, \qquad \xi^{j-1} \sim \nu_\xi^{j-1} \text{ independently} \\ Y^{j}& =h\left(X^{j}\right)+\eta^{j}, \qquad \eta^{j} \sim \nu_\eta^{j} \text{ independently} \end{align}\] for all \(j\geq 1\) and \(X^0 \sim \mu^0\). We couple \(\left\{(\tilde{X}^j_f, \tilde{Y}^j_f)\right\}_{j\geq 1}\) through \[\begin{align} \tilde{X}^j_f &= \Psi\left(\tilde{X}^{j-1}\right)+\tilde{\xi}^{j-1},\qquad \tilde{\xi}^{j-1} \sim \nu_\xi^{j-1} \text{ independently} \\ \tilde{Y}^{j}_f& =h\left(\tilde{X}^{j}_f\right)+\tilde{\eta}^{j}, \qquad \tilde{\eta}^{j} \sim \nu_\eta^{j} \text{ independently} \\ \tilde{X}^{j} &= T^{\tilde{\mu}^{f, j}_{XY}}_{Y^j} ( \tilde{X}^j_f, \tilde{Y}^j_f, \omega^j), \quad \omega^{j} \sim \kappa \text{ independently} \end{align}\] for all \(j\geq 1\) and \(\tilde{X}^0 \sim\mu^0\) independently.
Step I. We begin the proof by showing, inductively, that \[\Delta^j \mathrel{\vcenter{:}}= \|X^j - \tilde{X}_f^j\|_2\] is sub-Gaussian with deterministic parameter \(\bar t^j\). We will use this fact later to control \[\left\|\mathbb{E}\left(\tilde{Y}_f^j - Y^j \mid Y^{1:j-1}\right)\right\|_2 .\] For \(j = 1\), we have that \(\left\|\Delta^1\right\|_{\psi_2} \leq 2L_\Psi\sigma_X + 2 \sigma_{\max}=\bar t^1\). Let \(j\geq 2\) and assume \(\left\|\Delta^{j-1}\right\|_{\psi_2} \leq \bar t^{j-1}.\) Then, \[\begin{align} \Delta^j &= \left\|X^j -\tilde{X}^j_f\right\|_2 \;\leq\; L_\Psi\left\|\tilde{X}^{j-1}- X^{j-1}\right\|_2 +\left\|\xi^{j-1} \right\|_2 + \left\|\tilde{\xi}^{j-1}\right\|_2 . \end{align}\] Apply Equation 29 to bound \[\begin{align} \left\|\tilde{X}^{j-1}- X^{j-1}\right\|_2 &= \left\|T^{\tilde{\mu}^{f, j-1}_{XY}}_{Y^{j-1}} ( \tilde{X}^{j-1}_f, \tilde{Y}^{j-1}_f, \omega^{j-1}) - X^{j-1}\right\|_2 \\ &\leq \left\|T^{\tilde{\mu}^{f, j-1}_{XY}}_{Y^{j-1}} ( \tilde{X}^{j-1}_f, \tilde{Y}^{j-1}_f, \omega^{j-1}) - \tilde{X}^{j-1}_f\right\|_2 + \left\|X^{j-1} - \tilde{X}^{j-1}_f\right\|_2\\ & \leq L_y\left(1 + \text{Tr}\mathrm{Cov}\left( \tilde{\mu}^{f, j-1}_{XY}\right)^{e_y}\right) \left( \left\|\tilde{Y}^{j-1}_f - Y^{j-1} \right\|_2+ \left\|\omega^{j-1}\right\|_2\right)+ \Delta^{j-1}. \end{align}\] Step II. Further, we have \(\left\|\tilde{Y}^{j-1}_f - Y^{j-1} \right\|_2 \leq L_h\Delta^{j-1} + \left\|\tilde{\eta}^{j-1}\right\|_2 +\left\| \eta^{j-1} \right\|_2.\) Combining these estimates yields \[\begin{align} \left\|\Delta^j \right\|_{\psi_2} & \leq L_\Psi\Bigl( L_y\left(1 + \text{Tr}\mathrm{Cov}\left( \tilde{\mu}^{f, j-1}_{XY}\right)^{e_y}\right) \\ &\qquad\qquad\cdot \left( L_h\left\|\Delta^{j-1}\right\|_{\psi_2}+2\lambda_{\max}+\sigma_\kappa\right)\\ &\qquad\qquad + \left\|\Delta^{j-1}\right\|_{\psi_2}\Bigr)+2\sigma_{\max} \\ & \leq L_\Psi\Bigl( L_y\left(1 + \left(C_{\mathrm{tr}}\tilde{C}_{\mathrm{subG}}^2\right)^{e_y}\right) \\ &\qquad\qquad\cdot \left( L_h\bar t^{j-1}+2\lambda_{\max}+\sigma_\kappa\right) + \bar t^{j-1}\Bigr)+2\sigma_{\max} =\bar t^j\,, \end{align}\] and this completes the induction.

We have \[\begin{align} \left\|\mathbb{E}\left( \tilde{Y}^j_f - Y^j\left|Y^{1:j-1}\right.\right) \right\|_2& \leq \mathbb{E}\left( \left\|\tilde{Y}^j_f - Y^j\right\|_2\left|Y^{1:j-1}\right.\right) \\ &= \mathbb{E}\left( \left\|h\left(\tilde{X}^{j}_f\right)+\tilde{\eta}^{j} - h\left(X^{j}\right)-\eta^{j}\right\|_2\left|Y^{1:j-1}\right.\right) \\ &\leq \mathbb{E}\left( \left\|h\left(\tilde{X}^{j}_f\right) - h\left(X^{j}\right)\right\|_2\left|Y^{1:j-1}\right.\right) + 2\mathbb{E}\left( \left\|\eta^{j}\right\|_2\right) \\ & \leq L_h\mathbb{E}\left( \Delta^j\left|Y^{1:j-1}\right.\right)+ 2 \lambda_{\max}. \end{align}\] Here we used \(\mathbb{E}\|\eta^j\|_2\le \|\eta^j\|_{2,\psi_2}\le\lambda_{\max}\), which follows from Jensen’s inequality and the definition of \(\|\cdot\|_{2,\psi_2}\) in Equation 18 . Step I and Proposition 2 prove that \(\mathbb{E}\left( \Delta^j\left|Y^{1:j-1}\right.\right)\) is sub-Gaussian with parameter \(\bar t^j\). Further, Property 2 in Proposition 3 proves that \[\mathbb{P}\bigl(\mathbb{E}\left( \Delta^j\left|Y^{1:j-1}\right.\right)\ge \alpha\bigr)\leq 2\exp\!\bigl(-\alpha^2/K_2^2\bigr)\quad \text{ for all } \alpha\geq0\] for \(K_2 = C_{\psi_2}\bar t^j\). Therefore, for any \(k\geq 1\), with probability \(1 - \frac{1}{kJ}\), \[\mathbb{E}\left( \Delta^j\left|Y^{1:j-1}\right.\right)< C_{\psi_2}\bar t^j\sqrt{\log(2kJ)}\] and with probability \(1 - \frac{1}{k}\) \[\left\|\mathbb{E}\left( \tilde{Y}^j_f - Y^j\left|Y^{1:j-1}\right.\right) \right\|_2 \leq L_hC_{\psi_2}\bar t^j\sqrt{\log(2kJ)} +{2}\lambda_{\max} \text{ for all }j =1,\ldots, J.\] ◻

Proof. Let \(Y_i\mathrel{\vcenter{:}}= X_i-M\).

(i) By the monotonicity and centering properties of \(\psi_2\) norms (see [43]), for each \(v\in\mathbb{S}^{d-1}\) the marginal \(\langle v,Y_1\rangle\) is centered sub-Gaussian with parameter \(2K\). Therefore, \(\langle v,\overline{X}_N-M\rangle\) is sub-Gaussian with parameter \(2K/\sqrt{N}\). In particular, \(\|\overline{X}_N - M\|_{\psi_2} \leq \frac{2K}{\sqrt{N}}\). By Exercise 3.39 of [43], this implies \[\left\|\overline{X}_N-M\right\|_{2,\psi_2} \lesssim \sqrt{d/N}K.\] The claim then follows by Proposition 3 (Property 2), after choosing the universal constant \(c_1\) sufficiently large.

(ii) Because the \(Y_i\) are \(2K\)-sub-Gaussian and \(\widetilde{\Sigma}_N=\frac{1}{N}\sum_{i=1}^N (X_i-M)(X_i-M)^\top\), Theorem 6.5 of [42] gives, for all \(s\ge0\) and universal \(\tilde{c}_1, \tilde{c}_2, \tilde{c}_3\geq 0\), \[\mathbb{P}\!\left( \|\widetilde{\Sigma}_N-\Sigma\|_2 > \tilde{c}_14 K^2\left(\sqrt{\frac{d}{N}}+\frac{d}{N}+s\right) \right)\le \tilde{c}_2e^{-\tilde{c}_3N\min\{s,s^2\}}.\] Choosing \[s= \tilde{c}_4\left(\sqrt{u/N}+u/N\right)\] with a sufficiently large constant \(\tilde{c}_4\) gives \[\tilde{c}_2e^{-\tilde{c}_3N\min\{s,s^2\}} \leq e^{-u}.\] Therefore, with probability at least \(1 - e^{-u}\) \[\|\widetilde{\Sigma}_N-\Sigma\|_2 \le \tilde{c}_5K^2 \frac{d + u}{\sqrt{N}}\] for another universal constant \(\tilde{c}_5\). Finally, consider the event from (i). Since \(\Sigma_N=\widetilde{\Sigma}_N-(\overline{X}_N-M)(\overline{X}_N-M)^\top\), the decomposition \(\|\Sigma_N-\Sigma\|_2\le \|\widetilde{\Sigma}_N-\Sigma\|_2+\|\overline{X}_N-M\|_2^2\) concludes the proof.

(iii) Finally, \[\left|\text{Tr}(\widetilde{\Sigma}_N - \Sigma)\right| \leq d \left\|\widetilde{\Sigma}_N - \Sigma\right\|_2\] shows the first inequality. The other inequality follows from the same argument, utilizing the result from (ii). ◻

Proof. We will repeatedly use the following two standard bounds: if \(f\) is \(L\)-Lipschitz and \(Z\) has second moments \[\label{eq:two-sample} \overline{M}_2\big(\mathrm{Law}(f(Z))\big) \le L\overline{M}_2(\mathrm{Law}(Z))\tag{31}\] Further, for any joint random variable \((U,V) \in \mathbb{R}^n \times \mathbb{R}^m\), \[\label{eq:prod-moment-bound} \overline{M}_2(\mathrm{Law}(U,V))\le \overline{M}_2(\mathrm{Law}(U))+\overline{M}_2(\mathrm{Law}(V)) .\tag{32}\] Introduce the independent noise-pair and seed samples \[\varepsilon_\ell^{j-1}\mathrel{\vcenter{:}}= (\xi_\ell^{j-1},\eta_\ell^{j}) \sim \nu_\xi^{j-1} \otimes \nu_\eta^{j}, \qquad \omega_\ell^{j}\sim\kappa,\] and the empirical measures \[\begin{gather} \hat{\mu}^{a,0}_{X,N}=\tfrac1N\sum_{\ell=1}^N\delta_{x_{\ell}^{a,0}},\quad \hat{\mu}^{j-1}_{\varepsilon,N}=\tfrac1N\sum_{\ell=1}^N\delta_{\varepsilon_\ell^{j-1}},\\ \hat{\mu}^{f,j}_{\omega,N}=\tfrac1N\sum_{\ell=1}^N\delta_{\omega_\ell^{j}}, \quad \hat{\mu}^{j-1}_{X\varepsilon,N} = \tfrac1N\sum_{\ell=1}^N\delta_{(x_{\ell}^{a,j-1}, \varepsilon_\ell^{j-1})}. \end{gather}\] By Assumption 1, \(x_{\ell}^{a,0}\sim\mu^0\) are sub-Gaussian with parameter \(\sigma_X\), \(\xi_\ell^{j-1}\) and \(\eta_\ell^{j}\) are sub-Gaussian with parameters \(\sigma_{\max}\) and \(\lambda_{\max}\), and \(\omega_\ell^{j}\) has parameter \(\sigma_\kappa\). Further, note that for any two distributions \(\mu, \nu \in \mathcal{P}_2\left(\mathbb{R}^u\right)\), \(\overline{M}_2^2(\mu) = \text{Tr}\mathrm{Cov}(\mu) \leq |\text{Tr}\mathrm{Cov}(\mu) -\text{Tr}\mathrm{Cov}(\nu)| +\text{Tr}\mathrm{Cov}(\nu)\). Therefore, by (iii) in Lemma 2 and a union bound over \(j=1,\dots,J\), we have that the event \[\mathcal{E}_k\mathrel{\vcenter{:}}= \left\{ \max\left(\overline{M}_2(\hat{\mu}^{a,0}_{X,N}), \overline{M}_2(\hat{\mu}^{j-1}_{\varepsilon,N}), \overline{M}_2(\hat{\mu}^{f,j}_{\omega,N})\right)\le C_1 + C_2\sqrt{\frac{\log k}{\sqrt{N}}}\quad \forall j\le J \right\}\] has probability at least \(1-\frac{1}{k}\) for constants \(C_1, C_2\) depending only on \(\sigma_X\), \(\sigma_{\max}\), \(\lambda_{\max}\), \(\sigma_\kappa\), \(n\), \(m\), \(a\), \(J\). While the constant \(C_2\) bounds the concentration terms (e.g., \(|\text{Tr}\mathrm{Cov}(\hat{\mu}^{a,0}_{X,N}) - \text{Tr}\mathrm{Cov}(\mu^0)|\)), the constant \(C_1\) controls the constant terms (e.g., \(\text{Tr}\mathrm{Cov}( \mu^{0})\)). Next, we condition on \(\mathcal{E}_k\). Write \[d_{XY}^j\mathrel{\vcenter{:}}= \overline{M}_2(\hat{\mu}_{XY,N}^{f,j}), \qquad d_X^{j-1}\mathrel{\vcenter{:}}= \overline{M}_2(\hat{\mu}_{X,N}^{a,j-1}).\] Define the forecast map \(\Phi(x,\xi,\eta)\mathrel{\vcenter{:}}= \big(\Psi(x)+\xi, h(\Psi(x)+\xi)+\eta\big),\) and its Lipschitz constant inherited from Equation 1 : \(L_{PQ}\mathrel{\vcenter{:}}= (1+L_h)(1+L_\Psi).\) By the definitions in Equations 10 and 11 , \(\hat{\mu}_{XY,N}^{f,j}=\Phi_{\sharp}\hat{\mu}^{j-1}_{X\varepsilon,N},\) so by Equations 31 and 32 , \[\label{eq:forecast-moment} d_{XY}^j =\overline{M}_2\big(\Phi_{\sharp}\hat{\mu}^{j-1}_{X\varepsilon,N}\big) \;\le\; L_{PQ}\overline{M}_2(\hat{\mu}^{j-1}_{X\varepsilon,N}) \;\le\; L_{PQ}\left(d_X^{j-1}+C_1 +C_2\sqrt{\frac{\log k}{\sqrt{N}}}\right).\tag{33}\] Let \(\hat{\mu}^{f,j-1}_{XY\omega,N} =\tfrac1N\sum_{\ell=1}^N\delta_{(x_{\ell}^{f,j-1},y_{\ell}^{f,j-1},\omega_\ell^{j-1})}.\) Then \(\hat{\mu}_{X,N}^{a,j-1} =(T_{y^{j-1}}^{\hat{\mu}_{XY,N}^{f,j-1}})_{\sharp}\hat{\mu}^{f,j-1}_{XY\omega,N}.\) By Assumption 2 (1a), \[\mathrm{Lip}\Big(T_{y^{j-1}}^{\hat{\mu}_{XY,N}^{f,j-1}}\Big) \le C_{\mathrm{L}}\Big(1+\text{Tr}\mathrm{Cov}(\hat{\mu}_{XY,N}^{f,j-1})^{e_{\mathrm{L}}}\Big) = C_{\mathrm{L}}(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}).\] Another application of Equations 31 and 32 yields \[d_X^{j-1} \le C_{\mathrm{L}}\big(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}\big) \overline{M}_2(\hat{\mu}^{f,j-1}_{XY\omega,N}) \le C_{\mathrm{L}}\big(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}\big) \big(d_{XY}^{j-1}+\overline{M}_2(\hat{\mu}^{f,j-1}_{\omega,N})\big).\] We can simplify this further to \[\label{eq:analysis-moment} d_X^{j-1} \le C_{\mathrm{L}}\big(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}\big) \left(d_{XY}^{j-1}+C_1 +C_2\sqrt{\frac{\log k}{\sqrt{N}}}\right).\tag{34}\] Set \(A\mathrel{\vcenter{:}}= C_1 +C_2\sqrt{\frac{\log k}{\sqrt{N}}}\). Insert Equation 34 into Equation 33 to get, for all \(j\ge2\), \[\label{eq:main-recursion} d_{XY}^{j} \;\le\;L_{PQ}\Big( C_{\mathrm{L}}\big(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}\big)\big(d_{XY}^{j-1}+A\big) +A\Big).\tag{35}\] For \(j=1\), Equation 33 gives \[d_{XY}^{1}\leq 2 L_{PQ}A.\] Choosing sufficiently large \(\alpha\geq 1\) depending only on \(L_{PQ}\) and \(C_{\mathrm{L}}\), Equation 35 implies that, for all \(j\ge2\), \[\begin{align} 1 +d_{XY}^{j} &\le L_{PQ}\Big( C_{\mathrm{L}}\big(1+(d_{XY}^{j-1})^{2e_{\mathrm{L}}}\big)\big(d_{XY}^{j-1}+A\big) +A\Big) + 1 \\ & \leq \alpha\left(\big(1 +d_{XY}^{j-1}\big)^{2e_{\mathrm{L}} + 1} + \big(1 +d_{XY}^{j-1}\big)^{2e_{\mathrm{L}}} (1+A)\right) \end{align}\] and \(d_{XY}^{1}\le \alpha A\). Define \(\beta\mathrel{\vcenter{:}}= 2e_{\mathrm{L}} + 1\), \(r_1\mathrel{\vcenter{:}}=1\), \(r_j\mathrel{\vcenter{:}}=\beta r_{j-1}\) for \(j\ge 2\), and \(D_1\mathrel{\vcenter{:}}=\alpha\). Since \(\alpha\ge1\), we have \(1 +d_{XY}^{1} \le D_1 (1+A)^{r_1}\). Assume \(1 +d_{XY}^{j-1} \le D_{j-1} (1+A)^{r_{j-1}}\). Then \[\begin{align} 1 +d_{XY}^{j} &\le \alpha\left(\big(1 +d_{XY}^{j-1}\big)^{\beta} + \big(1 +d_{XY}^{j-1}\big)^{\beta-1} (1+A)\right)\\ &\le \alpha\left(D_{j-1}^\beta(1+A)^{\beta r_{j-1}} + D_{j-1}^{\beta-1} (1+A)^{(\beta-1)r_{j-1}+1}\right)\,. \end{align}\] Since \(r_{j-1}\ge 1\), we have \((\beta-1)r_{j-1}+1\le \beta r_{j-1}= r_j\). Hence, define \[D_j\mathrel{\vcenter{:}}= \alpha \left( D_{j-1}^\beta+D_{j-1}^{\beta-1}\right).\] Then \(1 +d_{XY}^{j} \le D_j (1+A)^{r_j}\), and so this estimate holds for any \(j\ge 1\) by induction. Define \(C_{\mathrm{ips}}\mathrel{\vcenter{:}}=\max_{1\le j \le J} D_j\). Then \(C_{\mathrm{ips}}\) depends only on \(L_\Psi,L_h,C_{\mathrm{L}},e_{\mathrm{L}},J\), and it holds \[1 +d^j_{XY} \leq C_{\mathrm{ips}} (1+ A)^{r_j}\] for \(r_j \mathrel{\vcenter{:}}= (1 +2 e_{\mathrm{L}})^{j-1}\) and all \(j\in\{1,\ldots J\}\). Since \(k\geq 2, N\geq 1\), we can redefine \(C_{\mathrm{ips}}\) such that \[d^j_{XY} \leq C_{\mathrm{ips}} \left(1+ \left(\frac{\log k}{\sqrt{N}}\right)^{r_j/2}\right)\] and introducing additional dependencies on \(\sigma_X,\sigma_{\max},\lambda_{\max},\sigma_\kappa,n,m,a\). ◻

Proof of Theorem 1. Whenever we say something depends on the allowed quantities, we mean it depends on \(\mathcal{Q}\). We write \(\lesssim\) if the constant in the bound depends only on these quantities. Define the empirical input measures (including the latent seed) as \[\hat{\mu}^{j}_{XY\omega,N}\mathrel{\vcenter{:}}= \frac{1}{N}\sum_{\ell=1}^N \delta_{(x_{\ell}^{f,j},y_{\ell}^{f,j},\omega_\ell^{j})}, \qquad \tilde{\mu}^{j}_{XY\omega,N}\mathrel{\vcenter{:}}= \frac{1}{N}\sum_{\ell=1}^N \delta_{(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})}.\] Let \((\tilde{X}^j_f, \tilde{Y}^j_f) \sim \tilde{\mu}_{XY}^{f,j}\) denote samples from the mean-field algorithm \(\{\tilde{\mu}_{XY}^{f, j}\}_{j=1}^J\) defined in Equation 20 , independent of all other random variables. We also abbreviate \[\delta_{k,N} = \frac{\log k}{\sqrt{N}},\] as this term will come up repeatedly throughout our argument. Since we use the synchronous coupling as defined in Equation 21 , note that \[\label{eq:forecast95coupling} \|(x_{\ell}^{f,j},y_{\ell}^{f,j})-(v_\ell^{f,j},o_\ell^{f,j})\|_2 \le L_{PQ}\|x_{\ell}^{a,j-1}-v_\ell^{a,j-1}\|_2,\tag{36}\] for \(L_{PQ}\mathrel{\vcenter{:}}= (1+L_h)(1+L_\Psi)\). Set \[\varepsilon^j\mathrel{\vcenter{:}}= \sqrt{\frac{1}{N}\sum\limits_{\ell= 1}^N\left\|x_{\ell}^{a,j} - v_\ell^{a,j}\right\|_2^2}.\] Then \(\varepsilon^0=0\) by initialization. Our argument will be inductive over \(j\). We start by fixing \(k\ge2\) and it will be helpful to condition on an event \(\mathcal{A}_k\) that we define as the intersection of:

• The dynamical-process moment bounds \[\overline{M}_{2}^2\left(\mu_{XY}^{f,j}\right) \lesssim (\log k)^2 \text{ and } \|Y^j - \mathbb{E}(Y^j|Y^{1:j-1})\|_2^2 \lesssim {\log k} \text{for all } j=1,\ldots,J,\] from Proposition 1.

• The mean-field stability bounds \[\left\|\mathbb{E}\left( \tilde{Y}^j_f - Y^j\left|Y^{1:j-1}\right.\right) \right\|_2\lesssim \sqrt{\log k}\text{ for all }j =1,\ldots, J,\] as in Proposition 7, using Proposition 6 to bound \(\text{Tr}\mathrm{Cov}\left(\tilde{\mu}_{X Y}^{f, j-1}\right)\).

• Defining \(\tilde{\mu}^{f,j}_{Y}\) as the \(Y\)-marginal of \(\tilde{\mu}^{f,j}_{XY}\), noting that \(o^{f,j}_\ell \sim \tilde{\mu}^{f,j}_{Y}\) i.i.d., and writing \(m_f^j\mathrel{\vcenter{:}}=\mathbb{E}(\tilde{Y}_f^j|Y^{1:j-1})\), \[\begin{align} \frac{1}{N}\sum\limits_{\ell = 1}^N\left\|o^{f,j}_\ell - m_f^j\right\|_2^2 &= \text{Tr}\left(\frac{1}{N}\sum\limits_{\ell = 1}^N \bigl(o^{f,j}_\ell - m_f^j\bigr) \bigl(o^{f,j}_\ell - m_f^j\bigr)^\top\right. \\ &\left.\qquad-\mathrm{Cov}(\tilde{Y}^{j}_f|Y^{1:j-1}) \right) +\text{Tr}\bigl(\mathrm{Cov}(\tilde{Y}^{j}_f|Y^{1:j-1})\bigr)\\ &\lesssim\delta_{k,N} +1 \qquad \text{for all } j=1,\ldots,J, \end{align}\] where Lemma 2 controls the first term and Proposition 6 the second.

• The moment bounds \(\overline{M}_2\big(\hat{\mu}_{XY,N}^{f,j}\big)\lesssim 1+ \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}/2}\) for all \(j=1,\ldots,J\), following from Proposition 8.

• Concentration of mean-field covariance and mean: for all \(j = 1, \ldots, J\) \[\begin{align} \|\mathbb{E}\tilde{\mu}^{f,j}_{XY,N}-\mathbb{E}\tilde{\mu}^{f,j}_{XY}\|_2 &\lesssim \sqrt{\frac{\log k}{N }} \,,\\ \|\mathrm{Cov}(\tilde{\mu}^{f,j}_{XY,N})-\mathrm{Cov}(\tilde{\mu}^{f,j}_{XY}) \|_2&\lesssim \delta_{k,N}\,, \end{align}\] using Proposition 6 and Lemma 2. We also include the corresponding direct concentration bound after applying \(g\) to the i.i.d.mean-field samples: \[\|\mathrm{Cov}(g_\sharp\tilde{\mu}^{f,j}_{XY,N})-\mathrm{Cov}(g_\sharp\tilde{\mu}^{f,j}_{XY}) \|_2\lesssim \delta_{k,N} \qquad \text{for all }j = 1,\ldots,J,\] which follows from the same concentration lemma applied to \(g_\sharp\tilde{\mu}^{f,j}_{XY}\), using the Lipschitz continuity of \(g\).

• Concentration of auxiliary noise: for all \(j = 1, ..., J\), \[\sum\limits_{\ell = 1}^N\frac{1}{N} \|\omega^j_\ell\|_2^2 \lesssim \delta_{k,N}+ 1,\] through a similar decomposition as for \(\frac{1}{N}\sum\limits_{\ell = 1}^N\left\|o^{f,j}_\ell - m_f^j\right\|_2^2\), relying on Lemma 2 and Proposition 4.

By a union bound and our standing assumptions, we can tune the confidence levels in the constituent events so that their intersection satisfies \(\mathbb{P}(\mathcal{A}_k)\;\ge\;1-\tfrac{1}{k}\) for all \(k\ge 2.\) In particular, although each marginal statement is stated with probability \(1-\tfrac{1}{k}\), we may assign them probabilities \(1-\tfrac{1}{C k}\) (for a fixed constant \(C\)) so that the total failure probability remains \(\le \tfrac{1}{k}\) and the resulting \(C\)-factors are absorbed into harmless constants. The lower bound \(2\) for \(k\) is chosen arbitrarily in all proofs simply to ensure \(\log k\) is lower bounded by some fixed positive constant. Now, let \(j\geq 1\) and use the triangle inequality to decompose \[\begin{align} \varepsilon^j &= \sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left\|T_{y^j}^{\hat{\mu}_{XY,N}^{f,j}}(x_{\ell}^{f,j},y_{\ell}^{f,j},\omega_\ell^{j}) - T_{y^j}^{\tilde{\mu}_{XY}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})\right\|_2^2} \\ &\leq \underbrace{\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left\|T_{y^j}^{\hat{\mu}_{XY,N}^{f,j}}(x_{\ell}^{f,j},y_{\ell}^{f,j},\omega_\ell^{j}) -T_{y^j}^{\hat{\mu}_{XY,N}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})\right\|_2^2}}_{\text{Term A}}\\ &+ \underbrace{\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left\|T_{y^j}^{\hat{\mu}_{XY,N}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j}) - T_{y^j}^{\tilde{\mu}_{XY}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})\right\|_2^2}}_{\text{Term B}} \end{align}\] Term A. Note that for \(j=1\), \((x_{\ell}^{f,1},y_{\ell}^{f,1}) =(v_\ell^{f,1},o_\ell^{f,1})\), and so Term A vanishes. For \(j\ge 2\), Assumption 2 (1a) and Equation 36 yield \[\begin{align} \text{Term A} & \le L_T^{\hat{\mu}_{XY,N}^{f,j}}{\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left\|(x_{\ell}^{f,j},y_{\ell}^{f,j},\omega_\ell^{j}) -(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})\right\|_2^2}}\\ & \leq L_T^{\hat{\mu}_{XY,N}^{f,j}}L_{PQ}\varepsilon^{j-1}. \end{align}\] On \(\mathcal{A}_k\), \[L_T^{\hat{\mu}_{XY,N}^{f,j}} =C_{\mathrm{L}}\bigl(1+\text{Tr}\mathrm{Cov}(\hat{\mu}_{XY,N}^{f,j})^{e_{\mathrm{L}}}\bigr) \lesssim 1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}e_\mathrm{L}}.\] Combining both inequalities, \[\label{eq:termA95final} \text{Term A} \lesssim \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}e_\mathrm{L}}\right)\varepsilon^{j-1}.\tag{37}\] Term B. Apply Assumption 2 (1b) with \(\nu=\hat{\mu}_{XY,N}^{f,j}\), \(\mu=\tilde{\mu}_{XY}^{f,j}\) \[\begin{align} \text{Term B} &= \sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left\|T_{y^j}^{\hat{\mu}_{XY,N}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j}) - T_{y^j}^{\tilde{\mu}_{XY}^{f,j}}(v_\ell^{f,j},o_\ell^{f,j},\omega_\ell^{j})\right\|_2^2} \\ & \lesssim \left(1 + \text{Tr}\mathrm{Cov}\left( \hat{\mu}_{XY,N}^{f,j}\right)^{e_{\mathrm{est},1}} + \text{Tr}\mathrm{Cov}(\tilde{\mu}_{XY}^{f,j})^{e_{\mathrm{est},2}} \right)\\ & \cdot\left\|\mathrm{Cov}(g_\sharp \hat{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp \tilde{\mu}_{XY}^{f,j})\right\|_2\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \left(\|\omega^j_\ell\|_2+ \|o^{f,j}_\ell - y^j\|_2\right)^2}\\ & \lesssim \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}e_{\mathrm{est},1}}\right)\left\|\mathrm{Cov}(g_\sharp \hat{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp \tilde{\mu}_{XY}^{f,j})\right\|_2\\ & \cdot\left(\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \|o^{f,j}_\ell - y^j\|_2^2}+ \sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \|\omega^j_\ell\|_2^2}\right) \end{align}\] We bound the covariance term and the innovation term separately. Note that \[\begin{align} \overline{M}_2^2(\tilde{\mu}_{XY,N}^{f,j}) &= \text{Tr}\left(\mathrm{Cov}(\tilde{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(\tilde{\mu}_{XY}^{f,j})\right) + \text{Tr}\mathrm{Cov}(\tilde{\mu}_{XY}^{f,j}) \\ & \lesssim 1 +\delta_{k,N} \end{align}\] and by inequality 36 , \(W_2\left(\hat{\mu}_{XY,N}^{f,j},\tilde{\mu}_{XY,N}^{f,j}\right) \lesssim \varepsilon^{j-1}\). Therefore, for the covariance term, we have by Lemma 5 (Appendix 5.4) \[\begin{align} \left\|\mathrm{Cov}(g_\sharp \hat{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp \tilde{\mu}_{XY}^{f,j})\right\|_2 &\leq \left\|\mathrm{Cov}(g_\sharp \hat{\mu}_{XY,N}^{f,j}) -\mathrm{Cov}(g_\sharp \tilde{\mu}_{XY,N}^{f,j}) \right\|_2\\ &\quad +\left\|\mathrm{Cov}(g_\sharp \tilde{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp \tilde{\mu}_{XY}^{f,j})\right\|_2 \\ &\lesssim \left(\overline{M}_2(g_\sharp \hat{\mu}_{XY,N}^{f,j}) + \overline{M}_2(g_\sharp \tilde{\mu}_{XY,N}^{f,j})\right)\\ &\quad\cdot W_2\left(g_\sharp \hat{\mu}_{XY,N}^{f,j}, g_\sharp \tilde{\mu}_{XY,N}^{f,j}\right)\\ &\quad+\left\|\mathrm{Cov}(g_\sharp \tilde{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp \tilde{\mu}_{XY}^{f,j})\right\|_2\\ &\lesssim \left(\overline{M}_2(\hat{\mu}_{XY,N}^{f,j}) +\overline{M}_2(\tilde{\mu}_{XY,N}^{f,j})\right)\\ &\quad\cdot W_2\left(\hat{\mu}_{XY,N}^{f,j},\tilde{\mu}_{XY,N}^{f,j}\right)\\ &\quad+\left\|\mathrm{Cov}(g_\sharp\tilde{\mu}_{XY,N}^{f,j}) - \mathrm{Cov}(g_\sharp\tilde{\mu}_{XY}^{f,j})\right\|_2\\ &\lesssim \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}/2}\right)\varepsilon^{j-1} + \delta_{k,N}\,, \end{align}\] where the last estimate holds on the event \(\mathcal{A}_k\). We can bound the innovation term as: \[\begin{align} \sqrt{\frac{1}{N}\sum\limits_{\ell = 1}^N\|o^{f,j}_\ell - y^j\|_2^2 } &\leq \sqrt{\frac{1}{N}\sum\limits_{\ell = 1}^N\|o^{f,j}_\ell - m_f^j\|_2^2} + \|m_f^j- \mathbb{E}(Y^{j}|Y^{1:j-1})\|_2 \\ &\quad +\|Y^j - \mathbb{E}(Y^{j}|Y^{1:j-1})\|_2 \\ &\lesssim \sqrt{\delta_{k,N}}+ 1+\sqrt{\log k} + \sqrt{\log k} \lesssim \sqrt{\log k}. \end{align}\] Finally, noting \(\sqrt{\sum\limits_{\ell = 1}^N\frac{1}{N} \|\omega^j_\ell\|_2^2}\lesssim \sqrt{\log k}\) on \(\mathcal{A}_k\), we can combine these estimates to \[\begin{align} \text{Term B} &\lesssim \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}e_{\mathrm{est},1}}\right)\left(\left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}/2}\right)\varepsilon^{j-1} + \delta_{k,N}\right) \sqrt{\log k} \\ & \lesssim \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}(e_{\mathrm{est},1} +\frac{1}{2})}\right)\sqrt{\log k}\varepsilon^{j-1} + \left(1 + \delta_{k,N}^{(1+2e_\mathrm{L})^{J-1}e_{\mathrm{est},1}}\right)\frac{(\log k)^{3/2}}{\sqrt{N}}. \end{align}\] Closing the recursion. Combining Term A and B yields \[\varepsilon^j \lesssim \left(1 + \delta_{k,N}^{e_{\mathrm{rec},1}}\right)\sqrt{\log k}\varepsilon^{j-1} + \left(1 + \delta_{k,N}^{e_{\mathrm{rec}, 2}}\right)\frac{(\log k)^{3/2}}{\sqrt{N}}\text{ for } j=1,\ldots,J,\qquad \varepsilon^0=0\] for \(e_{\mathrm{rec}, 1} \mathrel{\vcenter{:}}= (1+2e_\mathrm{L})^{J-1}\max\left(e_\mathrm{L}, e_{\mathrm{est},1} +\frac{1}{2}\right)\) and \(e_{\mathrm{rec}, 2} \mathrel{\vcenter{:}}= (1+2e_\mathrm{L})^{J-1}e_{\mathrm{est},1}\). Inducting this inequality over \(1\leq j\le J\) and defining \((q_1, p_1)= (3/2,e_{\mathrm{rec},2})\), \((q_{j+1}, p_{j+1}) = (q_j +\frac{1}{2}, p_j +e_{\mathrm{rec},1})\) shows we can bound this by \[\varepsilon^j\;\lesssim \left(1 + \delta_{k,N}^{p_j}\right)\frac{(\log k)^{q_j}}{\sqrt{N}}\,.\] For the induction step \(j\ge2\), if \(\varepsilon^{j-1}\lesssim \left(1 + \delta_{k,N}^{p_{j-1}}\right)\frac{(\log k)^{q_{j-1}}}{\sqrt{N}}\) holds, then we have \[\begin{align} \varepsilon^j & \lesssim \left(1 + \delta_{k,N}^{e_{\mathrm{rec},1}}\right)\sqrt{\log k}\varepsilon^{j-1} + \left(1 + \delta_{k,N}^{e_{\mathrm{rec}, 2}}\right)\frac{(\log k)^{3/2}}{\sqrt{N}}\\ & \lesssim \left(1 + \delta_{k,N}^{e_{\mathrm{rec},1}}\right)\sqrt{\log k}\left(1 + \delta_{k,N}^{p_{j-1}}\right)\frac{(\log k)^{q_{j-1}}}{\sqrt{N}} + \left(1 + \delta_{k,N}^{e_{\mathrm{rec}, 2}}\right)\frac{(\log k)^{3/2}}{\sqrt{N}}\\ & \lesssim\left(1 + \delta_{k,N}^{p_{j-1} + e_{\mathrm{rec},1}}\right)\frac{(\log k)^{q_{j-1} +\frac{1}{2}}}{\sqrt{N}} + \left(1 + \delta_{k,N}^{e_{\mathrm{rec}, 2}}\right)\frac{(\log k)^{3/2}}{\sqrt{N}}\\ & \lesssim\left(1 + \delta_{k,N}^{p_{j}}\right)\frac{(\log k)^{q_{j} }}{\sqrt{N}}\,. \end{align}\] Finally, the coupling \(\frac{1}{N}\sum_{\ell=1}^N\delta_{(x_{\ell}^{a,j},v_\ell^{a,j})}\) is an admissible coupling between \(\hat{\mu}_{X,N}^{a,j}\) and \(\tilde{\mu}_{X,N}^{a,j}\), and the parameters \(p_j, q_j\) can be explicitly expressed as \[\begin{align} p_j \mathrel{\vcenter{:}}=(1+2e_\mathrm{L})^{J-1}\left(e_{\mathrm{est},1} + (j-1)\max\left(e_\mathrm{L}, e_{\mathrm{est},1} +\frac{1}{2}\right)\right)\,,\qquad q_j=1+j/2\,. \end{align}\] Therefore, for all \(1\le j \le J\), \[W_2\bigl(\hat{\mu}_{X,N}^{a,j},\tilde{\mu}_{X,N}^{a,j}\bigr)\le \varepsilon^j \lesssim \left(1 + \left( \frac{\log k}{\sqrt{N}}\right)^{p_j}\right)\frac{(\log k)^{1 + j/2}}{\sqrt{N}},\] which gives the stated bound. ◻

Proof. We will show the three conditions required by Assumption 2.

(1a) By definition of \(T^{\text{K},\mu}_{y^\star}\), for all \(x,x^\prime \in \mathbb{R}^n\) and \(y,y^\prime \in \mathbb{R}^m\) it holds that \[\begin{align} &\left\| T^{\text{K},\mu}_{y^\star}(x,y) - T^{\text{K},\mu}_{y^\star}(x^\prime, y^\prime)\right\|_2\\ &\quad = \left\|(x - x^\prime) + \mathrm{Cov}\left(\mu\right)_{XY} \left(\mathrm{Cov}\left(h_\sharp\mu_X\right) + \Lambda\right)^{-1}(y^\prime - y)\right\|_2 \\ &\leq \left( 1 + \left\|\Lambda^{-1}\right\|_{2}\left\|\mathrm{Cov}\left(\mu\right)_{XY}\right\|_{2}\right)\left\|(x,y) - (x^\prime,y^\prime)\right\|_2 \\ &\leq \left( 1 + \left\|\Lambda^{-1}\right\|_{2}{\text{Tr}\mathrm{Cov}\left(\mu\right)}\right)\left\|(x,y) - (x^\prime,y^\prime)\right\|_2 \\ &\leq \left( 1 + \lambda_{\min}^{-1}\right)\left( 1 + {\text{Tr}\mathrm{Cov}\left(\mu\right)}\right)\left\|(x,y) - (x^\prime,y^\prime)\right\|_2, \end{align}\] where we use Equation 38 .

(1b) We abbreviate \(\Delta y = y^\star - y\) and set \[H_\rho\mathrel{\vcenter{:}}= \mathrm{Cov}\left(h_\sharp\rho_X\right), \qquad A_\rho\mathrel{\vcenter{:}}= H_\rho+\Lambda .\] Let \(x \in \mathbb{R}^n\), \(y, y^\star \in \mathbb{R}^m\), \(\nu, \mu \in \mathcal{P}_{2}(\mathbb{R}^n\times\mathbb{R}^m)\). Then, \[\begin{align} \left\| T^{\text{K},\nu}_{y^\star}(x,y) - T^{\text{K},\mu}_{y^\star}(x,y)\right\|_2 &= \left\|\mathrm{Cov}\left(\nu\right)_{XY}A_\nu^{-1}\Delta y - \mathrm{Cov}\left(\mu\right)_{XY}A_\mu^{-1}\Delta y \right\|_2 \\ &\leq \left\|\mathrm{Cov}\left(\nu\right)_{XY}A_\nu^{-1} - \mathrm{Cov}\left(\mu\right)_{XY}A_\mu^{-1}\right\| \left\|\Delta y \right\|_2 \\ &\leq \left\|\Lambda^{-1}\right\|_{2}\left\|\mathrm{Cov}\left(\nu\right)_{XY} - \mathrm{Cov}\left(\mu\right)_{XY}\right\|_{2}\left\|\Delta y \right\|_2 \\ & + \left\|\mathrm{Cov}(\nu)_{XY}\right\|_{2} \left\|A_\nu^{-1}-A_\mu^{-1}\right\|_{2} \left\|\Delta y \right\|_2 \\ &\leq \left\|\Lambda^{-1}\right\|_{2}\left\|\mathrm{Cov}\left(\nu\right) - \mathrm{Cov}\left(\mu\right)\right\|_{2}\left\|\Delta y \right\|_2 \\ & + \left\|\mathrm{Cov}(\nu)\right\|_{2} \left\|A_\nu^{-1}-A_\mu^{-1}\right\|_{2} \left\|\Delta y \right\|_2\\ &\leq \left\|\Lambda^{-1}\right\|_{2}\left\|\mathrm{Cov}\left(\nu\right) - \mathrm{Cov}\left(\mu\right)\right\|_{2} \left\|\Delta y \right\|_2 \\ & + \text{Tr}\mathrm{Cov}(\nu)\left\|\Lambda^{-1}\right\|_{2}^2 \left\| H_\nu - H_\mu\right\|_{2} \left\|\Delta y \right\|_2. \end{align}\] Stacking the covariance matrices together, i.e., inserting \(g(x,y) = ((x,y), h(x))\), verifies (1b).

(3) Finally, let \(\nu \in \mathcal{P}_2\left(\mathbb{R}^n\times\mathbb{R}^m\right), x\in\mathbb{R}^n,y\in\mathbb{R}^m, y^\star \in\mathbb{R}^m\). Then, \[\begin{align} \left\|T^{\text{K},\nu}_{y^\star}(x,y) - x \right\|_2 &= \left\| \mathrm{Cov}(\nu)_{XY}A_\nu^{-1}(y^\star - y)\right\|_2 \\ &\leq \lambda_{\min}^{-1}\left(1 + \text{Tr}\mathrm{Cov}(\nu)\right)\left\|y^\star - y\right\|_2. \end{align}\] ◻

Proof. By the sub-Gaussian assumption \[\mathbb{E}\exp\bigl(\| X-\mathbb{E}X\|_2^2/K^2\bigr) \leq 2< \infty.\] Since a translation does not change Wasserstein distances, we may assume \(\mathbb{E}X=0\) and thus \(\int e^{\|y\|_2^2/K^2}\mu(dy) < \infty.\) This is the exponential-moment assumption of Example 3.4 in [34], with \(\beta=2\ge p\) and constants depending only on \(d,p,K\). Let \(\mathcal{T}_p(\mu,\nu):=W_p(\mu,\nu)^p\). Example 3.4(a) with \(\beta=2\) gives, for all \(N\in\mathbb{N}\), \(x>0\), \[\mathbb{P}\bigl(\mathcal{T}_p(\mu,\mu_N)>x\bigr) \le C_1\exp\bigl(-c_1 N\phi(x)\bigr)\mathrm 1_{\{x\le1\}} + C_1\exp\bigl(-c_1 N x^{2/p}\bigr)\mathrm 1_{\{x\ge1\}},\] where the rate function \(\phi\) (depending on \(d,p\)) is \[\phi(x) = \begin{cases} x^2 & p>\frac{d}{2}\\ \bigl(x/\log(2+1/x)\bigr)^2 & p=\frac{d}{2}\\ x^{d/p} & 0<p<\frac{d}{2}. \end{cases}\] Plugging in \(x = t^p\) gives \[\mathbb{P}\bigl(W_p(\mu_N,\mu)>t\bigr) = \mathbb{P}\bigl(\mathcal{T}_p(\mu,\mu_N)>t^p\bigr) \le C_1\exp\bigl(-c_1 N\phi(t^p)\bigr) + C_1\exp\bigl(-c_1 N t^{2}\bigr).\] For \(t =A_1 \tilde{\gamma}_{k,N}^{p,d}\) the second term is upper bounded by \((2k)^{-1}\) for sufficiently large \(A_1 \geq 0\) since \(\tilde{\gamma}_{k,N}^{p,d} \geq \left(\log k/N \right)^{1/2}\) and \(\log(2k)\le 2\log k\) for \(k\ge2\). Similarly, plugging in \(t =A_1 \tilde{\gamma}_{k,N}^{p,d} \geq A_1\gamma_{k,N}^{p,d}\) for the first term and some elementary algebra in each of the three cases shows that it is smaller than \((2k)^{-1}\) for sufficiently large \(A_1 \geq 0\). The two bounds sum to the claimed failure probability \(k^{-1}\).

The expectation bound follows by the layer-cake representation. ◻

Proof of Proposition 10. Consider \(X^\prime\), an independent copy of \(X\). By repeatedly applying Jensen’s inequality and independence, we find \[\begin{align} \left\|f(X) - \mathbb{E}f(X)\right\|_{2,\psi_2}&= \inf \left\{t>0: \mathbb{E}\left(e^{ \left\|f(X) - \mathbb{E}f(X^\prime)\right\|_2^2 / t^{2}}\right) \leq 2\right\}\\ & \leq \inf \left\{t>0: \mathbb{E}\exp\left(\frac{\|f(X)-f(X^\prime)\|_2^2}{t^2}\right) \leq 2\right\}\\ & \leq \inf \left\{t>0: \mathbb{E}\left(e^{ \left\|X - X^\prime\right\|_2^2 / (t/L)^{2}}\right) \leq 2\right\}\\ & = \inf \left\{t>0: \mathbb{E}\left(e^{ \left\|(X- \mathbb{E}(X)) - (X^\prime - \mathbb{E}(X))\right\|_2^2 / (t/L)^{2}}\right) \leq 2\right\}\\ & \leq \inf \left\{t>0: \mathbb{E}\left(e^{ 2\left(\left\|X- \mathbb{E}(X)\right\|_2^2+ \left\|X^\prime - \mathbb{E}(X)\right\|_2^2\right) / (t/L)^{2}}\right) \leq 2\right\}\\ & = \inf \left\{t>0: \mathbb{E}\left(e^{ 2\left(\left\|(X- \mathbb{E}(X))\right\|_2^2\right) / (t/L)^{2}}\right)^2 \leq 2\right\}\\ & \leq \inf \left\{t>0: \mathbb{E}\left(e^{ 4\left\|X- \mathbb{E}(X)\right\|_2^2 / (t/L)^{2}}\right)\leq 2\right\}\\ &= \inf \left\{t>0: \mathbb{E}\left(e^{ \left\|X- \mathbb{E}(X)\right\|_2^2 / (t/2L)^{2}}\right)\leq 2\right\}\\ &= 2L \inf \left\{t>0: \mathbb{E}\left(e^{ \left\|X- \mathbb{E}(X)\right\|_2^2 / t^{2}}\right)\leq 2\right\}\\ &= 2L \left\|X - \mathbb{E}X\right\|_{2,\psi_2}. \end{align}\] ◻

Proof of Lemma 1. We prove by induction that \(\|X^j-\mathbb{E}(X^j)\|_2\) is sub-Gaussian.

For the base case, \(\left\|X^0\right\|_{2,\psi_2} \leq \sigma_X < \infty\) by Assumption 1. By \(\|X^0 - \mathbb{E}(X^0)\|_2 \leq \|X^0\|_2 + \mathbb{E}(\|X^0\|_2)\) we have \(\|X^0 - \mathbb{E}(X^0)\|_{2,\psi_2}\leq 2\sigma_X\) (see also Exercise 2.40e of [43] for tighter centering bounds). Let \(j\geq 0\) and \(\|X^{j} - \mathbb{E}(X^{j})\|_2\) be sub-Gaussian with norm bounded by a constant depending on the parameters in the statement. We have \(X^{j+1} = \Psi(X^j) + \xi^j\) and \(\|\Psi(X^j) - \mathbb{E}(\Psi(X^j))\|_2\) is sub-Gaussian by Proposition 10. Further, \[\|\Psi(X^j) + \xi^j - \mathbb{E}(\Psi(X^j) + \xi^j)\|_2 \leq \|\Psi(X^j) - \mathbb{E}(\Psi(X^j))\|_2 + \|\xi^j - \mathbb{E}(\xi^j)\|_2\] and by triangle inequality of the Orlicz norm the right-hand-side is sub-Gaussian. By sub-Gaussian summation properties, this also shows that the left-hand side \(\left\|X^{j+1} - \mathbb{E}(X^{j+1})\right\|_2\) is sub-Gaussian. This shows that all \(X^j, 1\leq j \leq J\) are sub-Gaussian with sub-Gaussian norm upper bounded by a parameter depending only on \(\sigma_X, \sigma_{\max}, L_\Psi, J\). \(\|Y^j - \mathbb{E}(Y^j)\|_2\) is sub-Gaussian by the same argument since \(Y^j = h(X^j) + \eta^j,\) \(h\) is Lipschitz, and \(\eta^j\) is sub-Gaussian. The final result follows with \[\left\|(X^j, Y^j) - \mathbb{E}(X^j, Y^j)\right\|_2 \leq \left\|X^j- \mathbb{E}(X^j)\right\|_2+ \left\|Y^j- \mathbb{E}(Y^j)\right\|_2.\] We emphasize that each of the cited results demonstrates the stability of sub-Gaussian behavior under operations such as addition, upper bounding, and Lipschitz transformations, with constants that depend solely on the parameters characterizing the sub-Gaussian and Lipschitz parameters. ◻

Proof of Proposition 2. By Jensen’s inequality, \[\begin{align} e^{\left\|\mathbb{E}(X|Y) / t\right\|_2^{2}} \le \mathbb{E}\left(e^{\|X / t\|_2^{2}}\,|\,Y\right) \end{align}\] conditional on \(Y\), and hence \[\begin{align} \left\|\mathbb{E}(X|Y)\right\|_{2,\psi_2} &= \inf \left\{t>0: \mathbb{E}\left(e^{\left\|\mathbb{E}(X|Y) / t\right\|_2^{2}}\right) \leq 2\right\} \\ &\leq \inf \left\{t>0: \mathbb{E}\left(e^{\|X / t\|_2^{2}}\right) \leq 2\right\} \\ &= \|X\|_{2,\psi_2} . \end{align}\] ◻

Proof. For \(c>0\), let \(f(x)=\exp({-c4^x})\). Since \(f\) is positive and decreasing on \([0,\infty)\), \[\sum_{j=0}^\infty f(j) \le f(0)+\int_0^\infty f(x)\,dx = \exp({-c})+\int_0^\infty \exp(-c4^x)\,dx.\] With \(t=4^x\) so \(dx=\frac{dt}{t\log 4}\), \[\int_0^\infty \exp(-c4^x)dx = \int_{1}^{\infty}\frac{\exp(-ct)}{t\log 4}dt = \frac{1}{\log 4}\int_{c}^{\infty}\frac{\exp(-u)}{u}du \le \frac{\exp(-c)}{c\log 4}.\] Hence \[\sum_{j=0}^\infty \exp(-c4^j) \le \exp({-c})\left(1+\frac{1}{c\log 4}\right).\] ◻

Proof of Proposition 4. Assume there exists a constant \(K_5\) such that for any \(|t|\geq 1/K_5\) it holds that \(\mathbb{E}(\exp(tX)) \leq \exp\left( t^2K_5^2\right).\) Fix \(x\geq 2K_5\). For any \(|\lambda| \geq1/K_5\), a Chernoff bound gives \[\mathbb{P}\bigl(X \geq x\bigr) \le \exp(-\lambda x) \mathbb{E}\left(\exp(\lambda X)\right) {\le} \exp\!\Bigl( -\lambda x + \lambda^{2}K_5^{2} \Bigr).\] Minimizing the exponent in \(\lambda\) gives \(\lambda^{\star}=x/(2K_5^{2})\), which is valid for the Chernoff bound since \(|\lambda^{\star}|=|x|/(2K_5^{2}) = x/(2K_5^{2}) \geq1/K_5\). Hence, \(\mathbb{P}\bigl(X \geq x\bigr) \le \exp\left(-\frac{x^{2}}{4K_5^{2}}\right).\) It follows by applying the same argument to \(-X\) and a union bound that for all \(x \geq 2K_5\), \(\mathbb{P}\bigl(|X | \geq x\bigr) \le 2\exp\!\left(-\frac{x^{2}}{4K_5^{2}}\right).\) For \(K_2 = \frac{2K_5}{\sqrt{\log 2}}\), we have that \[\mathbb{P}\bigl(|X | \geq x\bigr) \le 2\exp\!\left(-\frac{x^{2}}{K_2^2}\right)\] for all \(x\geq 0\): For \(x \geq 2K_5\) this follows from what we just showed, i.e., \(\mathbb{P}\bigl(|X | \geq x\bigr) \le 2\exp\left(-\frac{x^{2}}{4K_5^2}\right) \leq 2\exp\left(-\log 2\frac{x^{2}}{4K_5^2}\right) = 2\exp\left(-\frac{x^{2}}{K_2^2}\right)\), as \(\log 2 < 1\). For \(0 \leq x\leq 2K_5\), we have \(2\exp\left(-\frac{x^{2}}{K_2^2}\right) \geq 2\exp\left(-\log 2\frac{4K_5^2}{4K_5^2}\right) = 1\), so it is trivially true. This shows sub-Gaussianity by Proposition 3 with \(K_2 = \frac{2}{\sqrt{\log 2}}K_5\). Conversely, let \(X\) be sub-Gaussian and write \(K_1=\|X\|_{\psi_2}\). By Proposition 3, the tail bound in Property 2 holds for some \(K_2\le \tilde{C}_{\psi_2}K_1\). Then, the moment generating function exists for any \(t \in \mathbb{R}\). The case \(t=0\) is immediate, so let \(t > 0\) and apply the layer-cake representation as well as Property 2 in Proposition 3: \[\begin{align} \mathbb{E}\exp(tX) & = \int\limits_0^\infty \mathbb{P}\left(\exp(tX)\geq s\right)\mathrm{d}s \\ &\leq 1 + \int\limits_1^\infty \mathbb{P}\left(\exp(tX)\geq s\right)\mathrm{d}s\;=\; 1 + \int\limits_1^\infty \mathbb{P}\left(X\geq \frac{\log s }{t}\right)\mathrm{d}s \\ &= 1 + t\int\limits_0^\infty \exp(tu)\mathbb{P}\left(X\geq u\right)\mathrm{d}u \;\leq\; 1 + 2t\int\limits_{0}^\infty \exp(tu)\exp\left(-\frac{u^2}{K_2^2}\right)\mathrm{d}u\\ &\leq 1 + 2t\int\limits_{-\infty}^\infty \exp(tu)\exp\left(-\frac{u^2}{K_2^2}\right)\mathrm{d}u \;=\; 1 + 2\sqrt{\pi t^2K_{2}^2} \cdot \exp \left(\frac{t^{2} K_{2}^2}{4}\right) \\ \end{align}\] Now, note that the function \(\alpha x\exp(-\beta x^2)\) restricted to \(x\geq 0\) has its maximum value at \(\frac{\alpha}{\sqrt{2\beta}}\exp(-1/2)\). Choosing \(\alpha = 2 \sqrt{\pi}\) and \(\beta = 2\pi\) yields \(2 \sqrt{\pi} x \exp(-2\pi x^2) \leq 1\) for all \(x \geq 0\). Therefore, we can further bound, for \(x = K_2 t\), \[\begin{align} \mathbb{E}\exp(tX) &\leq 1 + 2\sqrt{\pi} x \exp\left(- 2\pi x^2\right) \cdot \exp \left( \frac{x^2}{4} + 2\pi x^2\right) \\ & \leq 1 + 1 \cdot \exp \left(\left(\frac{1}{4} + 2\pi\right){t^{2} K_{2}^2}\right) \\ &\leq \exp \left(\underbrace{\left(\frac{1}{4} + 2\pi\right)K_2^2}_{=:\tilde{K}_2^2}\cdot {t^{2}} + \log 2\right). \end{align}\] We can generalize this to all \(t\in \mathbb{R}\) by applying the same argument to \(-X\) which is sub-Gaussian with the same parameter to obtain \[\mathbb{E}\exp(tX)\leq \exp \left(\tilde{K}_2^2 {t^{2}} + \log 2\right) \forall t\in\mathbb{R}.\] Now, choose \(K_5 = 4 \tilde{K}_2\). Let \(|t|\geq \frac{1}{K_5}\). We have \[\begin{align} t^2 K_5^2 - \left({t^{2} \tilde{K}_{2}^2}+ \log2\right)& = t^2\tilde{K}_2^2\left( 16- 1\right)- \log2 \;\geq\; \frac{15}{16} - \log2 \;\geq\; 0. \end{align}\] In particular, \[\mathbb{E}\exp(tX) \leq \exp(K_5^2t^2), \quad \forall |t| \geq \frac{1}{K_5}\] and we are done with the reverse bound since \(K_5 = 4\sqrt{1/4+2\pi}K_2\). Thus, with \[C_{\psi_2}\mathrel{\vcenter{:}}= \max\left(1,\tilde{C}_{\psi_2},\frac{2\tilde{C}_{\psi_2}}{\sqrt{\log 2}},4\sqrt{1/4+2\pi}\tilde{C}_{\psi_2}\right)\,,\] the first direction gives \(\|X\|_{\psi_2}\le C_{\psi_2}K_5\) whenever 27 holds with parameter \(K_5\), while the reverse direction gives an admissible \(K_5\le C_{\psi_2}\|X\|_{\psi_2}\). We conclude by showing the monotonicity property. Assume Equation 27 holds and let \(\tilde{K}_5 > K_5\). For \(|t| \geq 1/K_5\) there is nothing to show, so let \(1/\tilde{K}_5 \leq |t| \leq 1/K_5\). Without loss of generality, \(t\geq 0\) since we can apply the same argument to \(-X\). Decompose convexly as \(t = (1-p) \cdot 0 + p\cdot 1/K_5\), so \(p=K_5t\). By log-convexity of the moment generating function, we have that \[\log \mathbb{E}(\exp(tX)) \leq p \log \mathbb{E}(\exp(X/K_5)) \leq p \log e = p.\] We are done since \(p=K_5t\le \tilde{K}_5t\le \tilde{K}_5^2t^2\), where the last inequality uses \(t\ge 1/\tilde{K}_5\). ◻

Proof of Proposition 6. Fix the observation path \(y^{1:J}\) throughout this proof. Construct the mean-field dynamics \(\tilde{Z}^j_f\mathrel{\vcenter{:}}= (\tilde{X}^j_f, \tilde{Y}^j_f) \sim \tilde{\mu}^{f,j}_{XY}\) at different times on a common probability space by drawing independently \(X^0\sim \mu^0\), \(\{\xi^{t}\}_{t \ge 0}\) with \(\xi^{t} \sim \nu_\xi^t\), \(\{\eta^{t}\}_{t \ge 1}\) with \(\eta^{t} \sim \nu_\eta^{t}\), \(\{\omega^{t}\}_{t \ge 1}\) with \(\omega^{t} \sim \kappa\), setting \(\tilde{X}^{a,0}\mathrel{\vcenter{:}}= X^0\), and defining recursively, for \(j \ge 1\), by \[\tilde{X}_f^{j} = \Psi(\tilde{X}^{a,j-1}) + \xi^{j-1}, \qquad \tilde{Y}^{f,j} = h(\tilde{X}_f^{j}) + \eta^{j}, \qquad \tilde{X}^{a,j} = T_{y^{j}}^{\tilde{\mu}_{XY}^{f,j}}\big(\tilde{Z}_f^{j}, \omega^{j}\big).\] Write \(s^j \mathrel{\vcenter{:}}= \big\|\tilde{Z}_f^{j} - \mathbb{E}\tilde{Z}_f^{j}\big\|_{2,\psi_2}\) and \(a^j \mathrel{\vcenter{:}}= \big\|\tilde{X}^{a,j} - \mathbb{E}\tilde{X}^{a,j}\big\|_{2,\psi_2}\). Let \[\Phi(x,\xi,\eta) \mathrel{\vcenter{:}}= \big(\Psi(x) + \xi, h(\Psi(x) + \xi) + \eta\big), \qquad L_{PQ} \mathrel{\vcenter{:}}= (1 + L_h)(1 + L_\Psi).\] We induct an upper bound \(d^j\) on \(s^j\) over \(j\geq 1\). By Assumption 1, \(\Phi\) is \(L_{PQ}\)-Lipschitz, and \[a^0 = \|\tilde{X}^{a,0}- \mathbb{E}\tilde{X}^{a,0}\|_{2,\psi_2} \le 2\sigma_X, \quad \|\xi^{j-1}\|_{2,\psi_2} \le \sigma_{\max}, \quad \|\eta^{j}\|_{2,\psi_2} \le \lambda_{\max}\] for all \(j \in \mathbb{N}_{\geq 1}\). Further, by Proposition 10 with \(f=\Phi\), \[\label{eq:s95and95a95equation} s^j \le 2 L_{PQ} \big(a^{j -1} +2 \sigma_{\max} + 2\lambda_{\max}\big).\tag{52}\] In particular, this implies that \[s^1 \le 4 L_{PQ} \big(\sigma_X + \sigma_{\max} + \lambda_{\max}\big) =: d^1.\] For the inductive step, we use that \(\mathrm{Tr}\mathrm{Cov}(\tilde{\mu}_{XY}^{f,j-1}) = \mathbb{E}\|\tilde{Z}_f^{j-1} - \mathbb{E}\tilde{Z}_f^{j-1}\|_2^2 \le 2( s^{j-1})^2\) so that Assumption 2 (1a) yields \[\mathrm{Lip}\big(T_{y^{j-1}}^{\tilde{\mu}_{XY}^{f,j-1}}\big) \le C_{\mathrm{L}} \big(1 + 2^{e_{\mathrm{L}}}( s^{j-1})^{2 e_{\mathrm{L}}}\big).\] Thus, by Proposition 10 \[\begin{align} a^{j -1} &= \left\|T_{y^{j-1}}^{\tilde{\mu}_{XY}^{f,j-1}}(\tilde{Z}_{j-1}^f, \omega^{j-1}) - \mathbb{E}\left(T_{y^{j-1}}^{\tilde{\mu}_{XY}^{f,j-1}}(\tilde{Z}_{j-1}^f, \omega^{j-1}) \right)\right\|_{2,\psi_2}\\ &\le 2 C_{\mathrm{L}} \big(1 + 2^{e_{\mathrm{L}}}( s^{j-1})^{2 e_{\mathrm{L}}}\big)\cdot \big(s^{j-1} + \sigma_\kappa\big). \end{align}\] Combining this with Equation 52 gives the recursion \[s^j \le A \big(1 + (s^{j -1})^{2 e_{\mathrm{L}}}\big)\big(s^{j -1} + \sigma_\kappa\big) + B\] for fixed constants \(A, B\geq 0\) depending only on \(L_{PQ}\), \(\sigma_{\max}, \lambda_{\max}, C_{\mathrm{L}},e_{\mathrm{L}},\) and \(\sigma_\kappa\) with initial bound \(s^1 \leq d^1\). Define for \(j \ge 2\), \[d^j \mathrel{\vcenter{:}}= A \big(1 + (d^{j-1})^{2 e_{\mathrm{L}}}\big)\big(d^{j-1} + \sigma_\kappa\big) + B.\] By monotonicity in \(d^{j-1}\), we have \(s^j \le d^j\) for all \(j \le J\). Since the recursion is at most polynomial, \(d^j < \infty\) for each finite \(J\). Setting \[\tilde{C}_{\mathrm{subG}} \mathrel{\vcenter{:}}= \max_{1 \le j \le J} d^j\] gives \(s^j \le \tilde{C}_{\mathrm{subG}}\) for all \(j \le J\), proving the claimed \(\psi_2\) control. ◻

Proof. Define \(\Delta\mathrel{\vcenter{:}}= \mathrm{Cov}(g_\sharp\nu)-\mathrm{Cov}(g_\sharp\mu)\). Write \(m_\nu\mathrel{\vcenter{:}}= \mathbb{E}g_\sharp \nu\) and \(m_\mu\mathrel{\vcenter{:}}= \mathbb{E}g_\sharp\mu\). For any coupling \(\pi\in\Pi(\nu,\mu)\), \[\Delta =\int\Big((g(z)-m_\nu)(g(z)-m_\nu)^\top-(g(z')-m_\mu)(g(z')-m_\mu)^\top\Big)d\pi(z,z').\] Using \(\|aa^\top-bb^\top\|_2\le(\|a\|_2+\|b\|_2)\|a-b\|_2\), we get \[\bigl\|\Delta\bigr\|_2 \le \int \big(\|g(z)-m_\nu\|_2+\|g(z')-m_\mu\|_2\big)\|g(z)-g(z')+m_\mu-m_\nu\|_2d\pi(z,z').\] Split the last norm and apply Cauchy–Schwarz: \[\begin{align} &\int \big(\|g(z)-m_\nu\|_2+\|g(z')-m_\mu\|_2\big)\|g(z)-g(z')\|_2d\pi(z,z') \\ \le& \Big(\mathbb{E}_\nu\|g-m_\nu\|_2^2+\mathbb{E}_\mu\|g-m_\mu\|_2^2\Big)^{\!1/2} \Big(\int\|g(z)-g(z')\|_2^2d\pi(z,z')\Big)^{1/2} \\ \le &\big(\overline{M}_2(g_\sharp\nu)+\overline{M}_2(g_\sharp\mu)\big)L_g\Big(\int\|z-z'\|_2^2d\pi(z,z')\Big)^{\!1/2}. \end{align}\] Noting that \[\|m_\mu-m_\nu\|_2\le L_g\Big(\int\|z-z'\|_2^2d\pi(z,z')\Big)^{1/2}.\] and combining the previous inequalities, \[\bigl\|\Delta\bigr\|_2 \le 2L_g\big(\overline{M}_2(g_\sharp\nu)+\overline{M}_2(g_\sharp\mu)\big) \Big(\int\|z-z'\|_2^2d\pi(z,z')\Big)^{1/2}.\] Optimizing over \(\pi\) yields the bound with \(W_2(\nu,\mu)\). ◻

Proof. For any \(z \in \mathbb{R}^n, y \in\mathbb{R}^m\), we will show that the equation \(S(x, y) = z\) has a unique solution. Since \(\alpha_k > 0\), we can rewrite component \(k\) of the inverse map as \[\begin{align} x_k &= \frac{1}{\alpha_k} \left(z_k - c_k - \sum_{i=1}^{\iota_k} \theta_k^i f_k^i(x_{1:k-1},y)\right). \end{align}\] The existence result for the inverse then follows by induction.

We now show the smoothness of the inverse map. For \(S = \mathcal{S}(\alpha, c, \theta)\) with \(\alpha \in \mathbb{R}_{\geq \alpha_{\min}}^n\), the Jacobian determinant is lower bounded by construction: \[|\det\nabla_x S(x,y)|\geq \alpha_{\min}^n.\] Then, the inverse function theorem guarantees that the inverse is continuously differentiable. ◻

Proof of Proposition [prop:existence95minimizer95sm95KL]. For any \(y\), the density of \(\left(S(\cdot, y)\right)^\sharp\gamma_{\text{G}}\) is given by \[\left(\left(S(\cdot, y)\right)^\sharp\gamma_{\text{G}}\right)(x) = {(2\pi)^{-n/2}} \exp\left(-\left\|S(x,y)\right\|_2^2/2\right) |\det\nabla_x S(x,y)|,\] where \(|\det\nabla_x S(x,y)| = |\prod_{k = 1}^n \alpha_k|\) is the Jacobian determinant of \(S(\cdot, y)\). Thus, the cross-entropy objective in Equation 46 is given by \[\begin{align} J(S) :=&-\int \mathrm{d}\mu(x,y) \log\left(\left(S(\cdot, y)\right)^\sharp\gamma_{\text{G}}\right)(x) \\ =& \int \mathrm{d}\mu(x,y) \left(\left\|S(x, y) \right\|_2^2/2 - \log|\det\nabla_x S(x,y)| \right) + \frac{n}{2}\log(2\pi) \\ =& \sum\limits_{k = 1}^n\int \mathrm{d}\mu(x,y) \left(|S_k(x, y)|^2/2 - \log|\partial_{x_k}S_k| \right) + \frac{n}{2}\log(2\pi). \end{align}\] For \(S = \mathcal{S}(\alpha, c, \theta)\), we can rewrite the \(k\)-th summand in terms of only the \(k\)-lower indexed parameters. Thus, the cross-entropy objective is separable into objectives \(J_k\) for each map component \(k=1,\dots,n\) given by: \[\begin{align} J_k\left(\alpha_k, \tilde{c}_k, \{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k}\right) &= \int \mathrm{d}\mu(x,y) \left(|S_k(x, y)|^2/2 - \log|\partial_{x_k}S_k| \right) \\ &= \frac{\alpha_k^2}{2}\int \mathrm{d}\mu(x,y)\left|x_k + \tilde{c}_k + \sum_{i=1}^{\iota_k} \tilde{\theta}_k^i f_k^i(x_{1:k-1},y)\right|^2 - \log\alpha_k, \end{align}\] where \(\tilde{\theta}^i_k = \frac{\theta^i_k}{\alpha_k}, \tilde{\theta}_k \mathrel{\vcenter{:}}= ( \tilde{\theta}_k^1, \ldots, \tilde{\theta}_k^{\iota_k})^\top\); see [49] and [51] for similar derivations of this fact.

It follows that for the map to have a unique minimizer, it is sufficient to show that for every \(k = 1,\ldots, n\), the function \(J_k\) constrained to \(\mathbb{R}_{\geq \alpha_{\min}} \times \mathbb{R}\times \mathbb{R}^{\iota_k}\) attains a minimum. \(J_k\) is of the form \[J_k\left(\alpha_k, \tilde{c}_k, \{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k}\right) = \frac{\alpha_k^2}{2}\tilde{J}_k(\tilde{c}_k, \{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k})- \log \alpha_k\] for \[\tilde{J}_k(\tilde{c}_k, \{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k}) =\int \mathrm{d}\mu(x,y) \left|x_k + \tilde{c}_k + \sum_{i=1}^{\iota_k} \tilde{\theta}_k^i f_k^i(x_{1:k-1},y)\right|^2.\] As \(J_k\) is strictly monotone in \(\tilde{J}_k\) and depends on \(\tilde{c}_k\) and \(\{\tilde{\theta}_k^i\}_{i=1,\ldots,\iota_k}\) only through \(\tilde{J}_k\), continuity allows us to take the minimum over \(\tilde{J}_k\), provided it exists, before taking the infimum over \(\alpha_k\). By the \(L_2\)-characterization of the expectation, for every fixed \(\{\tilde{\theta}_k^i\}_{i=1,\ldots,\iota_k}\), there is a unique minimizer in \(\tilde{c}_k\), given by \[\tilde{c}_k =- \int \mathrm{d}\mu(x,y) \left(x_k + \sum_{i=1}^{\iota_k} \tilde{\theta}_k^i f_k^i(x_{1:k-1},y)\right).\] Therefore, we must choose \(\{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k}\) to minimize \[\int \mathrm{d}\mu(x,y) \left| \overline{x}_k + \sum_{i=1}^{\iota_k} \tilde{\theta}_k^i \overline{f}_k^i(x_{1:k-1},y)\right|^2\] where \[\begin{align} \overline{x}_k &\mathrel{\vcenter{:}}= x_k - \int \mathrm{d}\mu(x^\prime,y^\prime)x^\prime_k \\ \overline{f}_k^i(x_{1:k-1},y) & \mathrel{\vcenter{:}}= f_k^i(x_{1:k-1},y) - \int \mathrm{d}\mu(x^\prime,y^\prime) f_k^i(x^\prime_{1:k-1},y^\prime). \end{align}\] We can rewrite this as the unconstrained quadratic program \[\inf_{\tilde{\theta}_k \in \mathbb{R}^{\iota_k}} \left(\mathrm{Cov}\left(X_k\right) + \tilde{\theta}_k^\top \mathrm{Cov}\left( f_k(X,Y) \right) \tilde{\theta}_k + 2\mathrm{Cov}\left( f_k(X,Y), X_k \right)^\top \tilde{\theta}_k\right).\] This quadratic program has a unique minimizer if and only if \(\mathrm{Cov}\left(f_k(X,Y) \right) \succ 0,\) which is given by \[\tilde{\theta}_k = - \mathrm{Cov}\left( f_k(X,Y) \right)^{-1} \mathrm{Cov}\left( f_k(X,Y), X_k \right).\] This shows that \(\tilde{J}_k\) has a unique minimizer that achieves \[\tilde{J}_k(\tilde{c}_k, \{\tilde{\theta}_k^i\}_{i = 1,\ldots, \iota_k}) = R_k(\mu) .\] Note that \(R_k(\mu) > 0\) by assumption. So far, we have shown that \[\inf\limits_{\alpha_k \in \mathbb{R}_{\geq \alpha_{\min}}, \tilde{c}_k\in \mathbb{R}, \tilde{\theta}_k \in \mathbb{R}^{\iota_k}}J_k\left(\alpha_k, \tilde{c}_k, \tilde{\theta}_k\right) = \inf\limits_{\alpha_k \in \mathbb{R}_{\geq \alpha_{\min}}}\frac{\alpha_k^2}{2} R_k(\mu)- \log \alpha_k\] and that this infimum is uniquely achieved in \(\tilde{c}_k, \tilde{\theta}_k\). The function \(Cx^2 - \log x\) is convex for every \(C > 0\) and therefore the infimum is also achieved in \(\alpha_k\), uniquely given by \[\label{eq:alpha95first95order} \alpha_k = \max\left(\frac{1}{\sqrt{R_k(\mu)}}, \alpha_{\min}\right).\tag{53}\]  ◻

Proof. Fix \(x, x^\prime\in\mathbb{R}^n, y, y', y^\star \in \mathbb{R}^m\) and let \(w = T^{\text{SM},\mu}_{y^\star}(x,y), w '= T^{\text{SM},\mu}_{y^\star}(x',y')\). By definition of Equation 48 , we have the recursion \[w_k = x_k - \theta_k^\top\left(f_k(x,y) -f_k(w,y^\star)\right)\] with \(\theta_k \mathrel{\vcenter{:}}= \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right)\) and \((X,Y) \sim \mu\). Therefore, \[\begin{align} |w_k - w_k'|^2 &= \left|x_k - x_k' - \theta_k^\top\left(f_k(x,y) - f_k(x',y') -(f_k(w,y^\star) - f_k(w',y^\star)) \right) \right|^2\\ &\leq 2 \left|x_k - x_k'\right|^2 +4 \|\theta_k\|_2^2\Bigl( \left\|f_k(x,y) - f_k(x',y')\right\|^2_2\\ &\qquad\qquad + \left\|f_k(w,y^\star) - f_k(w',y^\star) \right\|^2_2\Bigr)\\ &\leq 2 \left|x_k - x_k'\right|^2 +4 \|\theta_k\|_2^2L_f^2\Bigl( \left\|(x_{1:k-1},y) - (x_{1:k-1}',y')\right\|^2_2\\ &\qquad\qquad + \left\|w_{1:k-1} - w_{1:k-1}' \right\|^2_2\Bigr), \end{align}\] where the last inequality follows from \(f_k\) being \(L_f\)-Lipschitz. Let \(\Delta_{1:k}^2 := \|w_{1:k} - w_{1:k}'\|_2^2\) and \(\tilde{\Delta}^2 := \|(x,y) - (x',y')\|_2^2\). Then, another upper bound for the result above implies the recursion: \[\begin{align} \Delta_{1:k}^2 \leq \left(2 + 4 \|\theta_k\|_2^2L_f^2\right)\tilde{\Delta}^2 + \left(2 + 4 \|\theta_k\|_2^2L_f^2 \right)\Delta_{1:k-1}^2 \end{align}\] The solution to this recursion is \[\begin{align} \Delta_{1:k}^2 & \leq \tilde{\Delta}^2\sum\limits_{\ell = 1}^k \left(\prod\limits_{g =\ell}^{k}\left(2 + {4\left\| \theta_g\right\|_2^2 L_f^2} \right)\right) \\ & \leq \tilde{\Delta}^2 k \prod\limits_{g =1}^{k}\left(2 + {4\left\| \theta_g\right\|_2^2 L_f^2} \right). \end{align}\] Finally, we need to bound \(\|\theta_k\|_2^2\). By applying Cauchy-Schwarz and bounding the change in moments under the pushforward of Lipschitz transformations, we have \[\begin{align} \left\|\theta_k \right\|_2^2& = \left\|\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right)\right\|_2^2 \nonumber \\ & \leq \hat{\sigma}_f^{-2}\left\|\mathrm{Cov}\left( f_k( X, Y), X_k \right) \right\|_2^2 \nonumber \\ & \leq 4\hat{\sigma}_f^{-2}L_f^2 \text{Tr}\mathrm{Cov}(\mu)^2. \label{eq:theta95k95norm95bound} \end{align}\tag{54}\] Therefore, applying Jensen’s inequality repeatedly results in the following upper bound for the Lipschitz constant: \[\begin{align} \frac{\Delta_{1:n}}{\tilde{\Delta}} &\leq \sqrt{n} \sqrt{\prod\limits_{g = 1}^n\left(2 + 16\hat{\sigma}_f^{-2}L_f^4 \text{Tr}\mathrm{Cov}(\mu)^2\right)} \\ & \leq \sqrt{n} {\left(\sqrt{2} + 4\hat{\sigma}_f^{-1}L_f^{2} {\text{Tr}\mathrm{Cov}(\mu)}\right)^n} \\ & \leq \sqrt{n}2^{n-1} \left({2}^{n/2} + 4^n\hat{\sigma}_f^{-n}L_f^{2n} \text{Tr}\mathrm{Cov}(\mu)^{n}\right)\\ & \leq \sqrt{n}2^{n-1} \left({2}^{n/2} + 4^n\hat{\sigma}_f^{-n}L_f^{2n}\right)\left(1+ \text{Tr}\mathrm{Cov}(\mu)^{n}\right). \end{align}\] ◻

Proof. Let \[\theta_k \mathrel{\vcenter{:}}= \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right)\] for \((X, Y) \sim \mu\). Define \(w = T^{\text{SM},\mu}_{y^\star}(x,y)\). By definition from Equations 48 and 50 , for all \(k \geq 1\): \[\begin{align} \left|x_k - w_k\right|^2 & = \left| \theta_k^\top (f_k(x,y) - f_k(w, y^\star))\right|^2 \\ &\leq L_f^2\left\|\theta_k\right\|_2^2\left( \left\|x_{1:k-1} - w_{1:k-1}\right\|_2^2 +\|y - y^\star\|_2^2\right). \end{align}\] This shows inductively that \[\left\|x_{1:k} - w_{1:k}\right\|_2^2 \leq \left\|y - y^\star\right\|_2^2\sum\limits_{\ell = 1}^k \prod\limits_{g = \ell}^k (1 + L_f^2\left\|\theta_g\right\|_2^2).\] We can further bound this using elementary inequalities and Equation 54 as \[\begin{align} \left\|x - T^{\text{SM},\mu}_{y^\star}(x,y)\right\|_2^2 & \leq \left\|x_{1:n} - w_{1:n}\right\|_2^2 \\ &\leq n\left\|y - y^\star\right\|_2^2 \prod\limits_{g = 1}^n (1 + L_f^2\left\|\theta_g\right\|_2^2) \\ &\leq n\left\|y - y^\star\right\|_2^2 (1 + 4L_f^4\hat{\sigma}_f^{-2} \text{Tr}\mathrm{Cov}(\mu)^2)^n \\ &\leq n2^{n-1}\bigl(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n} \text{Tr}\mathrm{Cov}(\mu)^{2n}\bigr) \left\|y - y^\star\right\|_2^2 \end{align}\] ◻

Proof. By definition, \[A_\upsilon =\mathop{\mathrm{arg\,min}}_{\tilde{A} \succeq \upsilon I_{u}} f_{A}(\tilde{A}), \qquad f_A(\tilde{A}) \mathrel{\vcenter{:}}=\frac{1}{2}\left\|\tilde{A} - A\right\|_F^2.\] Let \(\langle\cdot,\cdot\rangle\) be the inner product inducing the Frobenius norm, i.e., \[\langle A, B\rangle \mathrel{\vcenter{:}}= \sum_{i,j=1}^u A_{ij}B_{ij}.\] This is a convex problem with \(f_A\) convex and differentiable. Moreover, the set of matrices \(\{\tilde{A} \in \mathbb{R}^{u \times u}: \tilde{A} \succeq \upsilon I_{u}\}\) where the optimization is performed is convex. Thus, by Equation 4.21 of [52], the KKT conditions give the following optimality criterion: \[\langle\nabla_{\tilde{A}}f_{A} (A_\upsilon), \tilde{A}-A_\upsilon\rangle = \langle A_\upsilon - A, \tilde{A}-A_\upsilon\rangle \geq 0 \text{ for all }\tilde{A}\succeq \upsilon I_u.\] Similarly, \[\langle B_\upsilon - B, \tilde{B}-B_\upsilon\rangle \geq 0 \text{ for all }\tilde{B}\succeq \upsilon I_u.\] In particular, \[\langle A_\upsilon - A, B_\upsilon-A_\upsilon\rangle \geq 0 ,\qquad \langle B_\upsilon - B, A_\upsilon-B_\upsilon\rangle \geq 0.\] Adding these together gives \[\langle A_\upsilon - B_\upsilon, B_\upsilon-A_\upsilon\rangle + \langle B - A, B_\upsilon-A_\upsilon\rangle \geq 0 \Leftrightarrow \langle B - A, B_\upsilon-A_\upsilon\rangle \geq \|A_\upsilon - B_\upsilon\|_F^2.\] The result then follows from Cauchy-Schwarz. ◻

Proof. Let \[\begin{align} \theta_k &\mathrel{\vcenter{:}}= \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right) \\ \theta_k' &\mathrel{\vcenter{:}}= \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right)^{-1} \mathrm{Cov}\left( f_k( X', Y'), X_k '\right) \end{align}\] for \((X, Y) \sim \mu, ( X', Y') \sim \nu\). Further, let \[w = T^{\text{SM},\mu}_{y^\star}(x,y), \qquad w' = T^{\text{SM},\nu}_{y^\star}(x,y).\] By definition from Equations@eq:eq:sm95approximate95conditioning95recursion and 50 , we have \[\begin{align} \left|w_k - w_k'\right|^2 & = \left|\theta_k^\top\left(f_k(x,y) -f_k(w,y^\star)\right) - \theta_k'^\top\left(f_k(x,y) -f_k(w',y^\star)\right) \right|^2\\ & = \left|\theta_k^\top\left(f_k(x,y) -f_k(x, y^\star) + f_k(x,y^\star)- f_k(w,y^\star)\right)\right. \\ &\left.- \theta_k'^\top\left(f_k(x,y) -f_k(x, y^\star) + f_k(x,y^\star) -f_k(w',y^\star)\right) \right|^2\\ & = \left|(\theta_k-\theta_k')^\top\left(f_k(x,y) -f_k(x, y^\star)\right)\right. \\ &\qquad\left. + (\theta_k -\theta_k')^\top(f_k(x,y^\star)- f_k(w,y^\star))\right. \\ & \left.- \theta_k'^\top\left(f_k(w,y^\star) -f_k(w',y^\star)\right) \right|^2\\ & \leq 3\Bigl(\|\theta_k -\theta_k'\|_2^2 \left(\|f_k(x,y) - f_k(x,y^\star)\|_2^2 + \|f_k(x,y^\star) - f_k(w,y^\star)\|_2^2 \right)\\ &\qquad + \|\theta_k'\|_2^2\|f_k(w,y^\star) - f_k(w',y^\star)\|_2^2\Bigr)\\ & \leq 3L_f^2\Bigl(\|\theta_k -\theta_k'\|_2^2 \left(\|y-y^\star\|_2^2 + \|x - w\|_2^2\right)\\ &\qquad + \|\theta_k'\|_2^2\|w_{1:k-1}-w'_{1:k-1}\|_2^2\Bigr). \end{align}\] Plugging in Lemma 8 and defining \(\Delta_{1:k}^2 = \|w_{1:k} - w_{1:k}'\|_2^2\), this yields \[\begin{align} \Delta_{1:k}^2 &\leq 3L_f^2\|\theta_k -\theta_k'\|_2^2 \left(1 + n2^{n-1}(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n} \text{Tr}\mathrm{Cov}(\mu)^{2n}) \right) \|y-y^\star\|_2^2 \\ &\quad+\left(1 + 3L_f^2\|\theta_k'\|_2^2\right)\Delta_{1:k-1}^2 \end{align}\] Further, applying Equation 54 to \(\theta_k'\) with \((X',Y')\sim\nu\) shows \[\begin{align} \Delta_{1:k}^2 &\leq 3L_f^2\|\theta_k -\theta_k'\|_2^2 \left(1 + n2^{n-1}(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n} \text{Tr}\mathrm{Cov}(\mu)^{2n}) \right) \|y-y^\star\|_2^2 \\ &\quad+\left(1 + 12L_f^4\hat{\sigma}_f^{-2} \text{Tr}\mathrm{Cov}(\nu)^2\right)\Delta_{1:k-1}^2. \end{align}\] Letting \[\begin{align} C_1& =3L_f^2\max_{1\leq k\leq n}\|\theta_k -\theta_k'\|_2^2 \left(1 + n2^{n-1}(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n} \text{Tr}\mathrm{Cov}(\mu)^{2n}) \right) \|y-y^\star\|_2^2 \\ C_2& = \left(1 + 12L_f^4\hat{\sigma}_f^{-2} \text{Tr}\mathrm{Cov}(\nu)^2\right), \end{align}\] this shows inductively that \[\Delta_{1:n}^2 \leq \sum\limits_{\ell = 1}^n C_2^{n-\ell} C_1 \leq n C_1 C_2^n.\] Now, we bound \[\begin{align} \|\theta_k -\theta_k'\|_2^2 &= \left\| \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right) \right.\\ &\left.- \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right)^{-1} \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2\\ &\leq 2\left\| \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X, Y), X_k \right) \right.\\ &-\left. \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2 \\ &+2\left\| \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1} \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right.\\ &-\left. \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right)^{-1} \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2 \\ &\leq 2\hat{\sigma}_f^{-2} \left\| \mathrm{Cov}\left( f_k( X, Y), X_k \right) - \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2 \\ &+2\left\|\mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2\left\|\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right)^{-1}\right.\\ &\left. \cdot \left(\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right) - \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right)\right) \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right)^{-1} \right\|_2^2 \\ &\leq 2\hat{\sigma}_f^{-2} \left\| \mathrm{Cov}\left( f_k( X, Y), X_k \right) - \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2 \\ &\quad+2\hat{\sigma}_f^{-4} \left\|\mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2\\ &\qquad\cdot \left\|\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right) - \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right) \right\|_2^2. \end{align}\] Set \[D_g\mathrel{\vcenter{:}}= \left\| \mathrm{Cov}\left( g( X, Y)\right) - \mathrm{Cov}\left( g( X', Y') \right)\right\|_2.\] Invoking Lemma 9, we find the following results: \[\begin{align} &\left\|\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right) - \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right) \right\|_2^2\\ &\qquad \leq \left\|\mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X, Y) \right) - \mathrm{Cov}_{\hat{\sigma}_f}\left( f_k( X', Y') \right) \right\|_F^2 \\ &\qquad \leq \left\|\mathrm{Cov}\left( f_k( X, Y) \right) - \mathrm{Cov}\left( f_k( X', Y') \right) \right\|_F^2 \\ &\qquad \leq bD_g^2.\\ &\left\| \mathrm{Cov}\left( f_k( X, Y), X_k \right) - \mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2\\ &\qquad \leq D_g^2. \end{align}\] Moreover, \[\left\|\mathrm{Cov}\left( f_k( X', Y'), X_k' \right)\right\|_2^2 \leq 2(L_f^2 + 1)^2\text{Tr}\mathrm{Cov}(\nu)^2.\] Thus, it follows that \[\begin{align} \|\theta_k -\theta_k'\|_2^2 &\leq C_3 \left(1 + \text{Tr}\mathrm{Cov}( \mu) + \text{Tr}\mathrm{Cov}(\nu)\right)^2 D_g^2 \end{align}\] for \(C_3 \mathrel{\vcenter{:}}= 2\hat{\sigma}_f^{-2}+4b(L_f^2 + 1)^2\hat{\sigma}_f^{-4}\). Plugging into \(\Delta_{1:n}^2\), we showed \[\begin{align} \Delta_{1:n}^2 &\leq n \left(1 + 12L_f^4\hat{\sigma}_f^{-2} \text{Tr}\mathrm{Cov}(\nu)^2\right)^n 3L_f^2\|\theta_k -\theta_k'\|_2^2\\ & \cdot \left(1 + n2^{n-1}(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n} \text{Tr}\mathrm{Cov}(\mu)^{2n}) \right)\|y-y^\star\|_2^2 \\ & \leq C_4\left(1 + \text{Tr}\mathrm{Cov}(\mu)^{2n}\right) \left(1 + \text{Tr}\mathrm{Cov}(\nu)^{2n}\right) \|y-y^\star\|_2^2\|\theta_k -\theta_k'\|_2^2\\ &\leq C_3C_4\left(1 + \text{Tr}\mathrm{Cov}(\mu)^{2n}\right) \left(1 + \text{Tr}\mathrm{Cov}(\nu)^{2n}\right)\\ &\quad\cdot \|y-y^\star\|_2^2 \left(1 + \text{Tr}\mathrm{Cov}( \mu) + \text{Tr}\mathrm{Cov}(\nu)\right)^2 \\ & \cdot D_g^2 \\ &\leq C_5 \|y-y^\star\|_2^2 \left(1 + \text{Tr}\mathrm{Cov}( \mu)^{2n+1} + \text{Tr}\mathrm{Cov}(\nu)^{2n+1}\right)^2\\ &\cdot D_g^2, \end{align}\] where we defined the intermediary constant \[\begin{align} C_4 &\mathrel{\vcenter{:}}= 3L_f^2n2^{n-1} \left(1 + 12^nL_f^{4n}\hat{\sigma}_f^{-2n} \right)\\ &\quad\cdot \left(1 + n2^{n-1}(1 + 4^nL_f^{4n}\hat{\sigma}_f^{-2n}) \right), \end{align}\] and \(C_5 = C_5(n,b,L_f,\hat{\sigma}_f)\) absorbs \(C_3C_4\) and the polynomial bound in \(\text{Tr}\mathrm{Cov}(\mu)\) and \(\text{Tr}\mathrm{Cov}(\nu)\). Taking the square root gives the final result. ◻