$\displaystyle\leq\frac{3}{2\alpha_{\mu}}\left(\frac{d_{\operatorname{\mathcal{X}}}}{2\lambda\eta}+\frac{M_{\mu}^{2}}{\lambda^{2}}+R_{\mu}^{2}\,\mathbb{E}_{X_{k-\ell+1}^{i}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}[\lVert
X_{k-\ell+1}^{i}\rVert^{2}]\right)$
$\displaystyle\leq\frac{3}{2\alpha_{\mu}}\left(\frac{d_{\operatorname{\mathcal{X}}}}{2\lambda\eta}+\frac{M_{\mu}^{2}}{\lambda^{2}}+R_{\mu}^{2}\left(\lVert
X_{k-\ell}^{i}\rVert^{2}\vee\mathfrak{s}^{\mu}\right)\right)$
$\displaystyle<\frac{3}{2\eta^{2}(L_{\mu}+\lambda
R_{\mu})^{2}}(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu}).$

Hence we arrive at the desired statement,

$\operatorname{\mathrm{KL}}(\tilde{\mu}_{k}^{i}\|\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\leq\left(1+\frac{3\exp(-(\ell-1)\alpha_{\mu}\lambda\eta)}{2\eta^{2}(L_{\mu}+\lambda
R_{\mu})^{2}}\right)(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu}).$

∎

The subsequent lemmas provide control over the Wasserstein distance between
pushforward distibutions. In particular, Lemma C.8 is the discrete analogue of
the $O(\beta_{t}/B_{t})$ time derivative bound obtained in the proof of
Proposition 3.3.

######  Lemma C.7.

For any two measures
$\mu^{(N)},\tilde{\mu}^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}}^{N})$
it holds that

$W_{2}(\Pi\mu^{(N)},\Pi\tilde{\mu}^{(N)})\leq\frac{1}{\sqrt{N}}W_{2}(\mu^{(N)},\tilde{\mu}^{(N)}).$

###### Proof.

Recall the dual formulation of $W_{2}$,

$W_{2}^{2}(\mu,\tilde{\mu})=\sup_{\phi,\psi}\bigg{\\{}\int\phi\mathop{}\\!\mathrm{d}\mu-\int\psi\mathop{}\\!\mathrm{d}\tilde{\mu}\;\bigg{|}\;\phi,\psi:\operatorname{\mathcal{X}}\to\operatorname{\mathbb{R}},\;\phi(x)-\psi(y)\leq\left\lVert
x-y\right\rVert^{2}\bigg{\\}}.$

Then for any pair of functions $\phi,\psi$ such that
$\phi(x)-\psi(y)\leq\left\lVert x-y\right\rVert^{2}$, the pullback functions
$\Pi^{*}\phi,\Pi^{*}\psi$ on $\operatorname{\mathcal{X}}^{N}$ satisfy

$\Pi^{*}\phi(\mathscr{X})-\Pi^{*}\psi(\mathscr{Y})=\frac{1}{N}\sum_{i=1}^{N}\phi(X^{i})-\psi(Y^{i})\leq\frac{1}{N}\sum_{i=1}^{N}\|X^{i}-Y^{i}\rVert^{2}=\frac{1}{N}\left\lVert\mathscr{X}-\mathscr{Y}\right\rVert_{L^{2}(\operatorname{\mathcal{X}}^{N})}^{2}.$

Therefore,

$\displaystyle\int_{\operatorname{\mathcal{X}}}\phi(x)\Pi\mu^{(N)}(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{X}}}\psi(x)\Pi\tilde{\mu}^{(N)}(\mathop{}\\!\mathrm{d}x)$
$\displaystyle=\int_{\operatorname{\mathcal{X}}^{N}}\Pi^{*}\phi(\mathscr{X})\mu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})-\int_{\operatorname{\mathcal{X}}^{N}}\Pi^{*}\psi(\mathscr{X})\tilde{\mu}^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})\leq\frac{1}{N}W_{2}^{2}(\mu^{(N)},\tilde{\mu}^{(N)}),$

which yields the assertion by taking the supremum over all permissible
$\phi,\psi$. ∎

######  Lemma C.8.

The projected 2-Wasserstein distance between
$\operatorname{\widehat{\mu}}_{k}^{(N)}$,
$\operatorname{\widehat{\mu}}_{k-1}^{(N)}$ is bounded as

$W_{2}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)},\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})\leq\frac{2M_{\mu}\beta_{k}}{\alpha_{\mu}\lambda
B_{k}}.$

###### Proof.

The proof is deferred to Section C.4. ∎

###  C.3 Proof of Proposition 3.6

We take $\ell=\ell^{\mu}=O(\eta^{-1}\log\eta^{-1})$ as defined in (9)
throughout the proof and only consider the case $k\geq 2\ell$ in Steps 1
through 4.

_Step 1._ We first look $\ell-1$ steps back to the past and control the
displacement of the proximal $\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}$
from the stationary state $\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)}$ of
the modified process via Lemma C.8, conditioning on the earlier history
$(\mathscr{X},\mathscr{Y})_{1:k-\ell}$.

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k-\ell+1:k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\qquad+M_{\mu}\sum_{j=1}^{\ell-1}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k-\ell+1:k-j}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[W_{1}(\Pi\operatorname{\widehat{\mu}}_{k-j-1}^{(N)},\Pi\operatorname{\widehat{\mu}}_{k-j}^{(N)})\right]$
$\displaystyle\leq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]+\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}\sum_{j=1}^{\ell-1}\frac{\beta_{k-j}}{B_{k-j}}.$

It is simple to further verify that

$\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}\sum_{j=1}^{\ell-1}\frac{\beta_{k-j}}{B_{k-j}}\leq\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}\frac{(r+1)(\ell-1)}{k-\ell+1}.$

_Step 2._ Next, we look back to the future and convert the expectation with
respect to $\mu_{\mathscr{X}_{k}}$ to the corresponding expectation for the
modified process. The incurred error can be bounded by utilizing Lemmas C.2,
C.3 and C.4 as

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\qquad-\mathbb{E}_{(\widetilde{\mathscr{X}},\widetilde{\mathscr{Y}})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\mu_{\widetilde{\mathscr{X}}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle=\mathbb{E}_{(\mathscr{X},\widetilde{\mathscr{X}},\mathscr{Y},\widetilde{\mathscr{Y}})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\bigg{[}\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\mu_{\widetilde{\mathscr{X}}_{k}})(\mathop{}\\!\mathrm{d}x)$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{X}}}\left(F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})-F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})\right)(\mu_{\widetilde{\mathscr{X}}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\bigg{]}$
$\displaystyle\leq\mathbb{E}_{(\mathscr{X},\widetilde{\mathscr{X}},\mathscr{Y},\widetilde{\mathscr{Y}})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\Bigg{[}M_{\mu}W_{1}(\mu_{\mathscr{X}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}})$
$\displaystyle\qquad+\frac{1}{N}\sum_{i=1}^{N}\left\lVert
F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})-F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})\right\rVert_{\operatorname{\mathrm{Lip}}}W_{1}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\Bigg{]}$
$\displaystyle\leq\frac{(r+1)M_{\mu}}{k-\ell+1}\mathfrak{w}_{\ell}^{\mu}$
$\displaystyle\qquad+\frac{(r+1)L_{\mu}}{k-\ell+1}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\mathbb{E}_{(\widetilde{X}_{k}^{i}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\Bigg{[}\bigg{(}\frac{2}{N}\sum_{i=1}^{N}\int_{\operatorname{\mathcal{X}}}\lVert\widetilde{X}_{k}^{i}-x\rVert^{2}\,\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)}(\mathop{}\\!\mathrm{d}x)\bigg{)}^{\frac{1}{2}}\Bigg{]}$
$\displaystyle\leq\frac{(r+1)M_{\mu}}{k-\ell+1}\mathfrak{w}_{\ell}^{\mu}+\frac{(r+1)L_{\mu}}{k-\ell+1}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\bigg{(}\frac{2}{N}\sum_{i=1}^{N}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{q}^{\mu}+2\mathfrak{s}^{\mu}\bigg{)}^{\frac{1}{2}}.$

_Step 3._ For the modified process, we apply a leave-one-out argument and
consider the expectation with respect to each particle $\widetilde{X}_{k}^{i}$
which is independent of $\widetilde{\mathscr{X}}_{k}^{-i}$,
$\widetilde{\mathscr{Y}}_{k}$ when conditioned on the stopped history
$(\mathscr{X},\mathscr{Y})_{1:k-\ell}$. That is,

$\displaystyle\mathbb{E}_{(\widetilde{\mathscr{X}},\widetilde{\mathscr{Y}})_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\mu_{\widetilde{\mathscr{X}}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}^{-i},\widetilde{\mathscr{Y}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\mathbb{E}_{\widetilde{X}_{k}^{i}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\delta_{\widetilde{X}_{k}^{i}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}^{-i},\widetilde{\mathscr{Y}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\mathbb{E}_{\widetilde{X}_{k}^{i}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}^{-i}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\delta_{\widetilde{X}_{k}^{i}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\qquad+\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\left\lVert
F(\mu_{\widetilde{\mathscr{X}}_{k}},\nu_{\widetilde{\mathscr{Y}}_{k}})-F(\mu_{\widetilde{\mathscr{X}}_{k}^{-i}},\nu_{\widetilde{\mathscr{Y}}_{k}})\right\rVert_{\operatorname{\mathrm{Lip}}}W_{1}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\right]$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}^{-i},\widetilde{\mathscr{Y}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\mathbb{E}_{\widetilde{X}_{k}^{i}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}^{-i}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\delta_{\widetilde{X}_{k}^{i}}-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\qquad+\frac{L_{\mu}}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[W_{1}(\mu_{\widetilde{\mathscr{X}}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}^{-i}})W_{1}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\right]$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}^{-i},\widetilde{\mathscr{Y}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\widetilde{\mathscr{X}}_{k}^{-i}},\nu_{\widetilde{\mathscr{Y}}_{k}})(\mu_{k}^{i}(\widetilde{X}_{k}^{i})-\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\qquad+\frac{L_{\mu}}{N}\sum_{i=1}^{N}\mathbb{E}_{\widetilde{\mathscr{X}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[W_{1}(\mu_{\widetilde{\mathscr{X}}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}^{-i}})W_{1}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\right]$
$\displaystyle\leq\frac{M_{\mu}}{N}\sum_{i=1}^{N}W_{1}(\mu_{k}^{i},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})$
$\displaystyle\qquad+\frac{L_{\mu}}{N}\sum_{i=1}^{N}\left(\mathbb{E}_{\widetilde{\mathscr{X}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[W_{2}^{2}(\mu_{\widetilde{\mathscr{X}}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}^{-i}})\right]\mathbb{E}_{\widetilde{\mathscr{X}}_{k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[W_{2}^{2}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})\right]\right)^{\frac{1}{2}}$
$\displaystyle\leq\frac{2M_{\mu}}{N}\sum_{i=1}^{N}\sqrt{\alpha_{\mu}^{-1}(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu})}$
$\displaystyle\qquad+\frac{2L_{\mu}}{N}\sum_{i=1}^{N}\Bigg{(}\frac{2\mathfrak{s}^{\mu}}{N}+\frac{1}{N(N-1)}\sum_{j\neq
i}\lVert X_{k-\ell}^{j}\rVert^{2}+\frac{1}{N}\lVert
X_{k-\ell}^{i}\rVert^{2}\Bigg{)}^{\frac{1}{2}}\\!\left(\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}$

by applying Lemma C.2, Lemma C.4 and Proposition C.6.

_Step 4._ Putting things together, we obtain the conditional bound

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{k-\ell+1:k}|(\mathscr{X},\mathscr{Y})_{1:k-\ell}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}\frac{(r+1)(\ell-1)}{k-\ell+1}$
$\displaystyle\qquad+\frac{(r+1)M_{\mu}}{k-\ell+1}\mathfrak{w}_{\ell}^{\mu}+\frac{(r+1)L_{\mu}}{k-\ell+1}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\bigg{(}\frac{2}{N}\sum_{i=1}^{N}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{q}^{\mu}+2\mathfrak{s}^{\mu}\bigg{)}^{\frac{1}{2}}$
$\displaystyle\qquad+\frac{2M_{\mu}}{N}\sum_{i=1}^{N}\sqrt{\alpha_{\mu}^{-1}(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu})}$
$\displaystyle\qquad+\frac{2L_{\mu}}{N}\sum_{i=1}^{N}\Bigg{(}\frac{2\mathfrak{s}^{\mu}}{N}+\frac{1}{N(N-1)}\sum_{j\neq
i}\lVert X_{k-\ell}^{j}\rVert^{2}+\frac{1}{N}\lVert
X_{k-\ell}^{i}\rVert^{2}\Bigg{)}^{\frac{1}{2}}\\!\left(\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}.$

Recalling $\mathbb{E}[\lVert X_{k-\ell}^{i}\rVert^{2}]\leq\mathbb{E}[\lVert
X_{1}^{i}\rVert^{2}]+\mathfrak{s}^{\mu}$ from Lemma C.2, taking the
expectation with respect to the history $(\mathscr{X},\mathscr{Y})_{1:k-\ell}$
finally gives

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\frac{r+1}{k-\ell+1}\left(\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}(\ell-1)+M_{\mu}\mathfrak{w}_{\ell}^{\mu}+L_{\mu}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\left(2\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+4\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}\right)$
$\displaystyle\qquad+2M_{\mu}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\frac{1}{\alpha_{\mu}N}\sum_{i=1}^{N}(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu})\right]^{\frac{1}{2}}$
$\displaystyle\qquad+\frac{L_{\mu}}{N^{\frac{3}{2}}}\sum_{i=1}^{N}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\frac{1}{N-1}\sum_{\smash[b]{j}=1}^{N}\lVert
X_{k-\ell}^{j}\rVert^{2}+\frac{2N-3}{N-1}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{q}^{\mu}+3\mathfrak{s}^{\mu}\right]$
$\displaystyle\leq\frac{r+1}{k-\ell+1}\left(\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}(\ell-1)+M_{\mu}\mathfrak{w}_{\ell}^{\mu}+L_{\mu}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\left(2\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+4\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}\right)$
$\displaystyle\qquad+2M_{\mu}\left(\frac{\mathfrak{K}^{\mu}\mathfrak{p}^{\mu}+\mathfrak{L}^{\mu}}{\alpha_{\mu}}\right)^{\frac{1}{2}}+\frac{2L_{\mu}}{\sqrt{N}}\left(3\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+6\mathfrak{s}^{\mu}\right)$
$\displaystyle\leq\frac{r+1}{k}C_{1}(\eta)+C_{2}\sqrt{\eta}+\frac{C_{3}}{\sqrt{N}},$

where the last bound holds if $k\geq 2\ell^{\mu}$. To be explicit,

$\displaystyle
C_{1}(\eta)=2\left(\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}(\ell-1)+M_{\mu}\mathfrak{w}_{\ell}^{\mu}+L_{\mu}\left(\mathfrak{w}_{\ell}^{\mu}+\mathfrak{w}_{\ell}^{\nu}\right)\left(2\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+4\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}\right),$
$\displaystyle C_{2}=2M_{\mu}\left(\frac{\bar{\eta}R_{\mu}^{2}(L_{\mu}+\lambda
R_{\mu})^{2}\mathfrak{p}^{\mu}}{\alpha_{\mu}^{2}}+\frac{(L_{\mu}+\lambda
R_{\mu})^{2}}{\alpha_{\mu}^{2}\lambda^{2}}\left(\bar{\eta}M_{\mu}^{2}+\lambda^{2}\bar{\eta}R_{\mu}^{2}\mathfrak{s}^{\mu}+\lambda
d_{\operatorname{\mathcal{X}}}\right)\right)^{\frac{1}{2}},$ $\displaystyle
C_{3}=2L_{\mu}\left(3\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+6\mathfrak{s}^{\mu}\right).$

The constants $C_{2},C_{3}$ can be taken to be polynomial and independent of
$\eta$ by substituting in the upper bound
$\bar{\eta}=\frac{r_{\mu}\lambda}{2(L_{\mu}+\lambda
R_{\mu})^{2}}\wedge\frac{r_{\mu}}{4\lambda R_{\mu}^{2}}$ in the expressions
for $\mathfrak{s}^{\mu},\mathfrak{K}^{\mu}/\eta,\mathfrak{L}^{\mu}/\eta$,
while $\ell^{\mu}=O(\eta^{-1}\log\eta^{-1})$. However, $C_{1}(\eta)$ contains
the dependency

$O(\mathfrak{w}_{\ell}^{\mu})=O\left(\eta^{-1}\exp(\ell
L_{\mu}\eta)\right)=O\Bigg{(}\frac{1}{\eta}\left(\frac{3}{2\eta^{2}(L_{\mu}^{2}+\lambda
R_{\mu}^{2})^{2}}\right)^{\frac{L_{\mu}}{\alpha_{\mu}\lambda}}\Bigg{)},$

which is a consequence of uniformly bounding the perturbation from the
gradient stopped process over a time period of $\ell$.

_Step 5._ For $k<2\ell$, proceeding similarly without converting to the
modified process gives

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\frac{M_{\mu}}{N}\sum_{i=1}^{N}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[W_{1}(\delta_{X_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})\right]$
$\displaystyle\qquad+\frac{L_{\mu}}{N}\sum_{i=1}^{N}\left(\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[W_{2}^{2}(\mu_{\widetilde{\mathscr{X}}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}^{-i}})\right]\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[W_{2}^{2}(\delta_{\widetilde{X}_{k}^{i}},\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})\right]\right)^{\frac{1}{2}}$
$\displaystyle\leq\frac{M_{\mu}}{N}\sum_{i=1}^{N}\left(\mathbb{E}[\lVert
X_{1}^{i}\rVert^{2}]+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}$
$\displaystyle\qquad+\frac{2L_{\mu}}{N}\sum_{i=1}^{N}\Bigg{(}\frac{2\mathfrak{s}^{\mu}}{N}+\frac{1}{N(N-1)}\sum_{k\neq
i}\mathbb{E}[\lVert X_{1}^{k}\rVert^{2}]+\frac{1}{N}\mathbb{E}[\lVert
X_{1}^{i}\rVert^{2}]\Bigg{)}^{\frac{1}{2}}\\!\left(\mathbb{E}[\lVert
X_{1}^{i}\rVert^{2}]+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}$
$\displaystyle\leq
M_{\mu}\sqrt{\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}}+\frac{2L_{\mu}(3\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+3\mathfrak{s}^{\mu})}{\sqrt{N}}$
$\displaystyle<\frac{C_{1}(\eta)}{2\ell}+\frac{C_{3}}{\sqrt{N}},$

where the final bound follows by noting
$\eta<\frac{r_{\mu}}{4L_{\mu}R_{\mu}}\leq\frac{1}{4L_{\mu}}$, hence by
expanding $(1+2\eta L_{\mu})^{\ell}$

$\displaystyle(\mathfrak{w}_{\ell}^{\mu})^{2}>\frac{1}{2L_{\mu}}\cdot\frac{M_{\mu}^{2}}{\eta^{2}L_{\mu}^{3}}\left(2\eta
L_{\mu}\ell+2\eta^{2}L_{\mu}^{2}\ell(\ell-1)\right)>\frac{M_{\mu}^{2}\ell^{2}}{L_{\mu}^{2}}$

and so

$C_{1}(\eta)>2L_{\mu}\mathfrak{w}_{\ell}^{\mu}\left(2\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+4\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}>2M_{\mu}\ell\left(\mathfrak{p}^{\mu}+\mathfrak{q}^{\mu}+\mathfrak{s}^{\mu}\right)^{\frac{1}{2}}.$

Thus the bound holds for all integers $k$. We conclude the proof by taking the
maximum with the corresponding quantities for $\nu$. ∎

###  C.4 Properties of Conjugate Functionals

We proceed to develop the $N$-particle lifted analogues
$J_{k}^{(N)},\widehat{J}_{k}^{(N)}$ of the conjugate functionals in the proof
of Theorem 3.4. In order to deal with time and particle discretization, we
will need a more precise characterization of their perturbative properties.
Many of the subsequent results do not follow from standard methods and
requires a careful synthesis of the discussion thus far.

######  Lemma C.9.

Given Lipschitz functions
$\zeta_{\mu}:\operatorname{\mathcal{X}}\to\operatorname{\mathbb{R}}$,
$\zeta_{\nu}:\operatorname{\mathcal{Y}}\to\operatorname{\mathbb{R}}$ and a
pair of $N$-particle probability measures
$\mu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}}^{N})$,
$\nu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{Y}}^{N})$
define the functional

$\displaystyle J_{k}^{(N)}(\mu^{(N)},\nu^{(N)}|\zeta^{\mu},\zeta^{\nu})$
$\displaystyle=-\int_{\operatorname{\mathcal{X}}^{N}}\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)\mu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})+\int_{\operatorname{\mathcal{Y}}^{N}}\int_{\operatorname{\mathcal{Y}}}\zeta^{\nu}(\nu_{\mathscr{Y}}-\rho^{\nu})(\mathop{}\\!\mathrm{d}y)\nu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{Y})$
$\displaystyle\qquad-\frac{\lambda
B_{k}}{N}\left(\operatorname{\mathrm{KL}}(\mu^{(N)}\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\nu^{(N)}\|\rho^{\nu\otimes N})\right).$

Then the maximum

$\widehat{J}_{k}^{(N)}(\zeta^{\mu},\zeta^{\nu})=\max_{\mu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}}^{N})}\max_{\nu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{Y}}^{N})}J_{k}^{(N)}(\mu^{(N)},\nu^{(N)}|\zeta^{\mu},\zeta^{\nu})$

exists for all $k\in\operatorname{\mathbb{N}}$ and is uniquely attained by the
pair of distributions

$\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})\propto\rho^{\mu\otimes
N}\exp\left(-\frac{N}{\lambda
B_{k}}\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}\mu_{\mathscr{X}}(\mathop{}\\!\mathrm{d}x)\right),\;\operatorname{\widehat{\nu}}_{k}^{(N)}(\zeta^{\nu})\propto\rho^{\nu\otimes
N}\exp\left(\frac{N}{\lambda
B_{k}}\int_{\operatorname{\mathcal{Y}}}\zeta^{\nu}\nu_{\mathscr{Y}}(\mathop{}\\!\mathrm{d}y)\right).$

###### Proof.

The proof is similar to Lemma B.3; we only check the first-order condition by
setting

$\displaystyle\frac{\delta
J_{k}^{(N)}}{\delta\mu^{(N)}}(\mu^{(N)})(\mathscr{X})=-\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)-\frac{\lambda
B_{k}}{N}\log\frac{\mu^{(N)}(\mathscr{X})}{\rho^{\mu\otimes
N}(\mathscr{X})}=\text{const.}$

∎

The $N$-particle proximal distributions
$\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\operatorname{\widehat{\nu}}_{k}^{(N)}(\zeta^{\nu})$,
despite being defined over the configuration spaces
$\operatorname{\mathcal{X}}^{N},\operatorname{\mathcal{Y}}^{N}$ also satisfy
the log-Sobolev inequality with the same constant as before due to the
tensorization property of entropy.

######  Lemma C.10 (product log-Sobolev inequality).

Suppose that $\zeta^{\mu}/B_{k},\zeta^{\nu}/B_{k}$ are
$M_{\mu},M_{\nu}$-Lipschitz, respectively. Then
$\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\operatorname{\widehat{\nu}}_{k}^{(N)}(\zeta^{\nu})$
satisfy the LSI on
$\operatorname{\mathcal{X}}^{N},\operatorname{\mathcal{Y}}^{N}$, with the same
constants $\alpha_{\mu},\alpha_{\nu}$ as in Proposition 3.2.

###### Proof.

We can write $\mu^{(N)}=\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})$
as the symmetric product distribution

$\mu^{(N)}(\mathscr{X})=\prod_{i=1}^{N}\mu^{i}(X^{i}),\quad\mu^{i}(X^{i})=\rho^{\mu}(X^{i})\exp\left(-\frac{\zeta^{\mu}(X^{i})}{\lambda
B_{k}}\right),\quad 1\leq i\leq N,$

where the marginals $\mu^{i}(X^{i})$ each satisfy the LSI with constant
$\alpha_{\mu}$ by Proposition 3.2. Also write $\mu^{-i}(X^{-i})=\prod_{j\neq
i}\mu^{i}(X^{i})$. For an appropriately integrable function $f$ on
$\operatorname{\mathcal{X}}^{N}$, denote by $f^{i}$ for the functions
$f^{i}(X^{i})=f(X^{1},\cdots,X^{i},\cdots,X^{N})$. Then by Proposition 2.2 of
Ledoux (1999),

$\operatorname{\mathrm{Ent}}_{\mu^{(N)}}(f^{2})\leq\sum_{i=1}^{N}\mathbb{E}_{\mu^{-i}}[\operatorname{\mathrm{Ent}}_{\mu^{i}}((f^{i})^{2})]\leq\sum_{i=1}^{N}\frac{2}{\alpha_{\mu}}\mathbb{E}_{\mu^{-i}}\mathbb{E}_{\mu^{i}}[\lVert\nabla
f^{i}\rVert^{2}]=\frac{2}{\alpha_{\mu}}\mathbb{E}_{\mu^{(N)}}[\left\lVert\nabla
f\right\rVert^{2}].$

∎

######  Lemma C.11.

The functional $\widehat{J}_{k}^{(N)}$ is convex in both arguments, and admits
functional derivatives at any $(\zeta^{\mu},\zeta^{\nu})$ which are given as

$\frac{\delta\widehat{J}_{k}^{(N)}}{\delta\zeta^{\mu}}(\zeta^{\mu},\zeta^{\nu})=-\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})+\rho^{\mu},\quad\frac{\delta\widehat{J}_{k}^{(N)}}{\delta\zeta^{\nu}}(\zeta^{\mu},\zeta^{\nu})=\Pi\operatorname{\widehat{\nu}}_{k}^{(N)}(\zeta^{\nu})-\rho^{\nu}.$

###### Proof.

Substituting
$\widehat{J}_{k}^{(N)}(\zeta^{\mu},\zeta^{\nu})=J_{k}^{(N)}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\nu})|\zeta^{\mu},\zeta^{\nu})$,

$\displaystyle\frac{\delta\widehat{J}_{k}^{(N)}}{\delta\zeta^{\mu}}(\zeta^{\mu},\zeta^{\nu})=-\frac{\delta}{\delta\zeta^{\mu}}\int_{\operatorname{\mathcal{X}}^{N}}\\!\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)\mu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})\bigg{|}_{\mu^{(N)}=\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})}$
$\displaystyle-\int_{\operatorname{\mathcal{X}}^{N}}\\!\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)\frac{\delta\operatorname{\widehat{\mu}}_{k}^{(N)}}{\delta\zeta^{\mu}}(\zeta^{\mu})(\mathop{}\\!\mathrm{d}\mathscr{X})-\frac{\lambda
B_{k}}{N}\int_{\operatorname{\mathcal{X}}^{N}}\left(\log\frac{\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})}{\rho^{\mu\otimes
N}}\right)\frac{\delta\operatorname{\widehat{\mu}}_{k}^{(N)}}{\delta\zeta^{\mu}}(\zeta^{\mu})(\mathop{}\\!\mathrm{d}\mathscr{X})$
$\displaystyle=\frac{\delta}{\delta\zeta^{\mu}}\left(-\int_{\operatorname{\mathcal{X}}^{N}}\frac{1}{N}\sum_{i=1}^{N}\zeta^{\mu}(X^{i})\mu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})+\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}\rho^{\mu}(\mathop{}\\!\mathrm{d}x)\right)\Bigg{|}_{\mu^{(N)}=\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})}$
$\displaystyle=-\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})+\rho^{\mu}.$

The integral over the configuration space measure
$\operatorname{\widehat{\mu}}_{k}^{(N)}$ therefore lifts the expectation with
respect to the discrete measure $\mu_{\mathscr{X}}$ to a differentiable
functional of $\zeta^{\mu}$, which in turn pushes forward
$\operatorname{\widehat{\mu}}_{k}^{(N)}$ onto the space
$\operatorname{\mathcal{X}}$. ∎

The following proposition is crucial to controlling the evolution of the
conjugate functional as well as the proximal distributions over time.

######  Proposition C.12.

Suppose $\zeta^{\mu}/B_{k},\tilde{\zeta}^{\mu}/B_{k}$ are $M_{\mu}$-Lipschitz
functions such that the difference $\zeta^{\mu}-\tilde{\zeta^{\mu}}$ is
$m_{\mu}$-Lipschitz for some $m_{\mu}>0$. Then the projected proximal
distributions satisfy

$W_{2}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))\leq\frac{m_{\mu}}{\alpha_{\mu}\lambda
B_{k}}.$

###### Proof.

Taking the first-order conditions

$\displaystyle-\int_{\operatorname{\mathcal{X}}}\zeta^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)-\frac{\lambda
B_{k}}{N}\log\frac{\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})}{\rho^{\mu\otimes
N}}=\text{const.},$
$\displaystyle-\int_{\operatorname{\mathcal{X}}}\tilde{\zeta}^{\mu}(\mu_{\mathscr{X}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)-\frac{\lambda
B_{k}}{N}\log\frac{\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})}{\rho^{\mu\otimes
N}}=\text{const.}$

Subtracting both sides and integrating over the difference
$\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})-\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})$,
we obtain

$\displaystyle-\int_{\operatorname{\mathcal{X}}^{N}}\int_{\operatorname{\mathcal{X}}}(\zeta^{\mu}-\tilde{\zeta}^{\mu})\mu_{\mathscr{X}}(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})-\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))(\mathop{}\\!\mathrm{d}\mathscr{X})$
(10) $\displaystyle=\frac{\lambda
B_{k}}{N}\int_{\operatorname{\mathcal{X}}^{N}}\log\frac{\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})}{\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})-\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))(\mathop{}\\!\mathrm{d}\mathscr{X}).$

Now the left-hand side of (10) can be bounded from above by

$\displaystyle-\int_{\operatorname{\mathcal{X}}^{N}}\int_{\operatorname{\mathcal{X}}}(\zeta^{\mu}-\tilde{\zeta}^{\mu})\mu_{\mathscr{X}}(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})-\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))(\mathop{}\\!\mathrm{d}\mathscr{X})$
$\displaystyle=-\int_{\operatorname{\mathcal{X}}}(\zeta^{\mu}-\tilde{\zeta}^{\mu})(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})-\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))(\mathop{}\\!\mathrm{d}x)$
$\displaystyle\leq
m_{\mu}W_{1}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))\leq
m_{\mu}W_{2}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})),$

while the right-hand side of (10) is bounded from below by

$\displaystyle\frac{\lambda
B_{k}}{N}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu})\|\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})\|\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}))\right)$
$\displaystyle\geq\frac{\alpha_{\mu}\lambda
B_{k}}{N}W_{2}^{2}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu}))$
$\displaystyle\geq\alpha_{\mu}\lambda
B_{k}W_{2}^{2}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\zeta^{\mu}),\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}(\tilde{\zeta}^{\mu})),$

where we have used Talagrand’s inequality from Lemma C.10 and the $W_{2}$
pushforward bound from Lemma C.7. Combining the two results yields the desired
statement. ∎

Denote the unnormalized aggregate derivatives as

$\delta_{k}^{\mu}=\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}}),\quad\delta_{k}^{\nu}=\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})$

so that
$\operatorname{\widehat{\mu}}_{k}^{(N)}=\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\mu})$,
$\operatorname{\widehat{\nu}}_{k}^{(N)}=\operatorname{\widehat{\nu}}_{k}^{(N)}(\delta_{k}^{\nu})$.
The functions $\delta_{k}^{\mu}/B_{k}$ and $\delta_{k}^{\nu}/B_{k}$ are
$M_{\mu}$\- and $M_{\nu}$-Lipschitz, respectively, due to Assumption 2. Lemma
C.11 and Proposition C.12 then allow us to quantify the change in
$\widehat{J}_{k}^{(N)}(\delta_{k}^{\mu},\delta_{k}^{\nu})$ as time progresses.

######  Lemma C.13.

We have the following one-step relation for $\widehat{J}_{k}^{(N)}$, $k\geq
2$:

$\displaystyle\widehat{J}_{k}^{(N)}(\delta_{k}^{\mu},\delta_{k}^{\nu})-\widehat{J}_{k-1}^{(N)}(\delta_{k-1}^{\mu},\delta_{k-1}^{\nu})$
$\displaystyle\leq\beta_{k}\int_{\operatorname{\mathcal{X}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(-\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}+\rho^{\mu})(\mathop{}\\!\mathrm{d}x)+\beta_{k}\int_{\operatorname{\mathcal{Y}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})(\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}-\rho^{\nu})(\mathop{}\\!\mathrm{d}y)$
$\displaystyle\qquad-\frac{\lambda\beta_{k}}{N}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{k}^{(N)}\\!\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{k}^{(N)}\\!\|\rho^{\nu\otimes
N})\right)+\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\frac{\beta_{k}^{2}}{2\lambda
B_{k-1}}.$

###### Proof.

By the maximality of $\widehat{J}_{k}^{(N)}$,

$\displaystyle\widehat{J}_{k}^{(N)}(\delta_{k}^{\mu},\delta_{k}^{\nu})=J_{k}^{(N)}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\mu}),\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\nu})|\delta_{k}^{\mu},\delta_{k}^{\nu})$
$\displaystyle=J_{k-1}^{(N)}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\mu}),\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\nu})|\delta_{k}^{\mu},\delta_{k}^{\nu})-\frac{\lambda\beta_{k}}{N}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{k}^{(N)}(\delta_{k}^{\mu})\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{k}^{(N)}(\delta_{k}^{\nu})\|\rho^{\nu\otimes
N})\right)$
$\displaystyle\leq\widehat{J}_{k-1}^{(N)}(\delta_{k}^{\mu},\delta_{k}^{\nu})-\frac{\lambda\beta_{k}}{N}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{k}^{(N)}\\!\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{k}^{(N)}\\!\|\rho^{\nu\otimes
N})\right).$

Further defining the interpolations

$\delta_{k}^{\mu}(s)=\delta_{k-1}^{\mu}+s(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})=\sum_{j=1}^{k-1}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})+s\beta_{k}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}}),\quad
0\leq s\leq 1$

and similarly for $\delta_{k}^{\nu}(s)$, we have

$\displaystyle\widehat{J}_{k-1}^{(N)}(\delta_{k}^{\mu},\delta_{k}^{\nu})-\widehat{J}_{k-1}^{(N)}(\delta_{k-1}^{\mu},\delta_{k-1}^{\nu})=\int_{0}^{1}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}s}\widehat{J}_{k-1}^{(N)}(\delta_{k}^{\mu}(s),\delta_{k}^{\nu}(s))\mathop{}\\!\mathrm{d}s$
$\displaystyle=\int_{0}^{1}\int_{\operatorname{\mathcal{X}}}(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})\frac{\delta\widehat{J}_{k-1}^{(N)}}{\delta\zeta^{\mu}}(\delta_{k}^{\mu}(s),\delta_{k}^{\nu}(s))(\mathop{}\\!\mathrm{d}x)+\int_{\operatorname{\mathcal{Y}}}(\delta_{k}^{\nu}-\delta_{k-1}^{\nu})\frac{\delta\widehat{J}_{k-1}^{(N)}}{\delta\zeta^{\nu}}(\delta_{k}^{\mu}(s),\delta_{k}^{\nu}(s))(\mathop{}\\!\mathrm{d}y)\mathop{}\\!\mathrm{d}s$
$\displaystyle=\int_{0}^{1}\int_{\operatorname{\mathcal{X}}}-(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\delta_{k}^{\mu}(s))(\mathop{}\\!\mathrm{d}x)+\int_{\operatorname{\mathcal{Y}}}(\delta_{k}^{\nu}-\delta_{k-1}^{\nu})\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}(\delta_{k}^{\nu}(s))(\mathop{}\\!\mathrm{d}y)\mathop{}\\!\mathrm{d}s$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{X}}}(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})\rho^{\mu}(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{Y}}}(\delta_{k}^{\nu}-\delta_{k-1}^{\nu})\rho^{\nu}(\mathop{}\\!\mathrm{d}y)$
$\displaystyle\leq\int_{0}^{1}\int_{\operatorname{\mathcal{X}}}-(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\delta_{k-1}^{\mu})(\mathop{}\\!\mathrm{d}x)+\int_{\operatorname{\mathcal{Y}}}(\delta_{k}^{\nu}-\delta_{k-1}^{\nu})\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}(\delta_{k-1}^{\nu})(\mathop{}\\!\mathrm{d}y)\mathop{}\\!\mathrm{d}s$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{X}}}(\delta_{k}^{\mu}-\delta_{k-1}^{\mu})\rho^{\mu}(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{Y}}}(\delta_{k}^{\nu}-\delta_{k-1}^{\nu})\rho^{\nu}(\mathop{}\\!\mathrm{d}y)$
$\displaystyle\qquad+\int_{0}^{1}M_{\mu}\beta_{k}W_{1}(\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\delta_{k}^{\mu}(s)),\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\delta_{k-1}^{\mu}))\mathop{}\\!\mathrm{d}s$
$\displaystyle\qquad+\int_{0}^{1}M_{\nu}\beta_{k}W_{1}(\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}(\delta_{k}^{\nu}(s)),\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}(\delta_{k-1}^{\nu}))\mathop{}\\!\mathrm{d}s$
$\displaystyle\leq\beta_{k}\int_{\operatorname{\mathcal{X}}}-\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\mathop{}\\!\mathrm{d}x)+\beta_{k}\int_{\operatorname{\mathcal{Y}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})\Pi\operatorname{\widehat{\nu}}_{k-1}^{(N)}(\mathop{}\\!\mathrm{d}y)$
$\displaystyle\qquad+\beta_{k}\int_{\operatorname{\mathcal{X}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})\rho^{\mu}(\mathop{}\\!\mathrm{d}x)-\beta_{k}\int_{\operatorname{\mathcal{Y}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})\rho^{\nu}(\mathop{}\\!\mathrm{d}y)+\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\frac{\beta_{k}^{2}}{2\lambda
B_{k-1}},$

where for the first inequality we used the fact that
$\delta_{k}^{\mu}-\delta_{k-1}^{\mu}$ is $M_{\mu}\beta_{k}$-Lipschitz, and for
the second we applied Proposition C.12 with $m_{\mu}=sM_{\mu}\beta_{k}$. ∎

We now give the promised proof of the pushforward evolution bound.

Proof of Lemma C.8. Note that
$\operatorname{\widehat{\mu}}_{k-1}^{(N)}=\operatorname{\widehat{\mu}}_{k-1}^{(N)}(\delta_{k-1}^{\mu})$
may also be written as

$\operatorname{\widehat{\mu}}_{k-1}^{(N)}=\operatorname{\widehat{\mu}}_{k}^{(N)}\Bigg{(}\frac{B_{k}}{B_{k-1}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})\Bigg{)}=\operatorname{\widehat{\mu}}_{k}^{(N)}\left(\frac{B_{k}}{B_{k-1}}\delta_{k-1}^{\mu}\right).$

Since $\delta_{k-1}^{\mu}/B_{k-1}$ is $M_{\mu}$-Lipschitz and

$\delta_{k}^{\mu}-\frac{B_{k}}{B_{k-1}}\delta_{k-1}^{\mu}=-\frac{\beta_{k}}{B_{k-1}}\sum_{j=1}^{k-1}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})+\beta_{k}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{k}},\nu_{\mathscr{Y}_{k}})$

is $2M_{\mu}\beta_{k}$-Lipschitz, by Proposition C.12 we obtain the bound

$W_{2}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)},\Pi\operatorname{\widehat{\mu}}_{k-1}^{(N)})\leq\frac{2M_{\mu}\beta_{k}}{\alpha_{\mu}\lambda
B_{k}}.$

∎

###  C.5 Proof of Theorem 3.7

_Step 1._ We first prove a convergent upper bound of the following surrogate
$\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})$
for the NI error of the average distributions. Note that the defining maximum
is lifted to the configuration space and the discrete empirical distributions
have been replaced with their proximal counterparts for measuring relative
entropy. While $\mathfrak{N}$ is not exactly the desired quantity, it arises
naturally from the discrete conjugate argument and helps to bound the expected
error.

$\displaystyle\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})$
$\displaystyle:=\\!\max_{\mu^{(N)},\nu^{(N)}}-\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\Pi\mu^{(N)},\nu_{\mathscr{Y}_{j}})-\frac{\lambda}{N}\operatorname{\mathrm{KL}}(\mu^{(N)}\|\rho^{\mu\otimes
N})+\frac{\lambda}{NB_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu\otimes
N})$
$\displaystyle\qquad+\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu_{\mathscr{X}_{j}},\Pi\nu^{(N)})-\frac{\lambda}{N}\operatorname{\mathrm{KL}}(\nu^{(N)}\|\rho^{\nu\otimes
N})+\frac{\lambda}{NB_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu\otimes
N})$
$\displaystyle\leq\max_{\mu^{(N)},\nu^{(N)}}-\int_{\operatorname{\mathcal{X}}^{N}}\int_{\operatorname{\mathcal{X}}}\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\mu_{\mathscr{X}}-\mu_{\mathscr{X}_{j}})(\mathop{}\\!\mathrm{d}x)\mu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{X})$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{Y}}^{N}}\int_{\operatorname{\mathcal{Y}}}\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\nu_{\mathscr{Y}}-\nu_{\mathscr{Y}_{j}})(\mathop{}\\!\mathrm{d}y)\nu^{(N)}(\mathop{}\\!\mathrm{d}\mathscr{Y})$
$\displaystyle\qquad-\frac{\lambda}{N}\left(\operatorname{\mathrm{KL}}(\mu^{(N)}\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\nu^{(N)}\|\rho^{\nu\otimes N})\right)$
$\displaystyle\qquad+\frac{\lambda}{NB_{k}}\sum_{j=1}^{k}\beta_{j}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu\otimes
N})\right)$
$\displaystyle=\frac{1}{B_{k}}\Bigg{[}\widehat{J}_{k}^{(N)}(\delta_{k}^{\nu},\delta_{k}^{\nu})+\frac{\lambda}{N}\sum_{j=1}^{k}\beta_{j}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu\otimes
N})\right)$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{X}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\mu_{\mathscr{X}_{j}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{Y}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\nu_{\mathscr{Y}_{j}}-\rho^{\nu})(\mathop{}\\!\mathrm{d}y)\Bigg{]},$

due to the convex-concavity of $\operatorname{\mathcal{L}}$. Recursively
applying Lemma C.13 then yields

$\displaystyle\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})$
$\displaystyle\leq\frac{1}{B_{k}}\Bigg{[}\sum_{j=2}^{k}\left(\widehat{J}_{j}^{(N)}(\delta_{j}^{\nu},\delta_{j}^{\nu})-\widehat{J}_{j-1}^{(N)}(\delta_{j-1}^{\nu},\delta_{j-1}^{\nu})\right)+\frac{1}{B_{k}}\widehat{J}_{1}^{(N)}(\delta_{1}^{\nu},\delta_{1}^{\nu})$
$\displaystyle\qquad+\int_{\operatorname{\mathcal{X}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\mu_{\mathscr{X}_{j}}-\rho^{\mu})(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{Y}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\nu_{\mathscr{Y}_{j}}-\rho^{\nu})(\mathop{}\\!\mathrm{d}y)$
$\displaystyle\qquad+\frac{\lambda}{N}\sum_{j=1}^{k}\beta_{j}\left(\operatorname{\mathrm{KL}}(\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu\otimes
N})+\operatorname{\mathrm{KL}}(\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu\otimes
N})\right)\Bigg{]}$
$\displaystyle\leq\frac{1}{B_{k}}\Bigg{[}\sum_{j=1}^{k}\beta_{j}\int_{\operatorname{\mathcal{X}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\mu_{\mathscr{X}_{j}}-\Pi\operatorname{\widehat{\mu}}_{j-1}^{(N)})(\mathop{}\\!\mathrm{d}x)$
$\displaystyle\qquad-\sum_{j=1}^{k}\beta_{j}\int_{\operatorname{\mathcal{Y}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\nu_{\mathscr{Y}_{j}}-\Pi\operatorname{\widehat{\nu}}_{j-1}^{(N)})(\mathop{}\\!\mathrm{d}y)+\frac{1}{2\lambda}\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\sum_{j=2}^{k}\frac{\beta_{j}^{2}}{B_{j-1}}\Bigg{]},$

where the initial term is substituted as
$\widehat{J}_{1}^{(N)}(\delta_{1}^{\nu},\delta_{1}^{\nu})=J_{1}^{(N)}(\operatorname{\widehat{\mu}}_{1}^{(N)},\operatorname{\widehat{\nu}}_{1}^{(N)}\\!|\delta_{1}^{\nu},\delta_{1}^{\nu})$
with the convention that
$\operatorname{\widehat{\mu}}_{0}^{(N)}=\operatorname{\widehat{\mu}}_{1}^{(N)},\operatorname{\widehat{\nu}}_{0}^{(N)}=\operatorname{\widehat{\nu}}_{1}^{(N)}$.
Now taking the expectation over the full history and applying Proposition 3.6,
we arrive at

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})\right]$
$\displaystyle\leq\frac{1}{B_{k}}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\sum_{j=1}^{k}\beta_{j}\int_{\operatorname{\mathcal{X}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\mu_{\mathscr{X}_{j}}-\Pi\operatorname{\widehat{\mu}}_{j-1}^{(N)})(\mathop{}\\!\mathrm{d}x)$
$\displaystyle\qquad-\sum_{j=1}^{k}\beta_{j}\int_{\operatorname{\mathcal{Y}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu_{\mathscr{X}_{j}},\nu_{\mathscr{Y}_{j}})(\nu_{\mathscr{Y}_{j}}-\Pi\operatorname{\widehat{\nu}}_{j-1}^{(N)})(\mathop{}\\!\mathrm{d}y)+\frac{1}{2\lambda}\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\sum_{j=2}^{k}\frac{\beta_{j}^{2}}{B_{j-1}}\Bigg{]}$
$\displaystyle\leq\frac{1}{B_{k}}\Bigg{[}2\sum_{j=1}^{k}\beta_{j}\left(\frac{r+1}{j}C_{1}(\eta)+C_{2}\sqrt{\eta}+\frac{C_{3}}{\sqrt{N}}\right)+\frac{1}{2\lambda}\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\sum_{j=2}^{k}\frac{\beta_{j}^{2}}{B_{j-1}}\Bigg{]}$
$\displaystyle\leq\left(\frac{(r+1)^{2}}{rk}+O(k^{-2})\right)\left(2C_{1}(\eta)+\frac{1}{2\lambda}\bigg{(}\frac{M_{\mu}^{2}}{\alpha_{\mu}}+\frac{M_{\nu}^{2}}{\alpha_{\nu}}\bigg{)}\right)+2C_{2}\sqrt{\eta}+\frac{2C_{3}}{\sqrt{N}}$
$\displaystyle\leq\left(\frac{(r+1)^{2}}{rk}+O(k^{-2})\right)\cdot\frac{9}{4}C_{1}(\eta)+2C_{2}\sqrt{\eta}+\frac{2C_{3}}{\sqrt{N}}$

by simply using $\ell>1$. For $r=0$, the last expression is replaced by the
exact bound

$\frac{1+\log
k}{k}\cdot\frac{9}{4}C_{1}(\eta)+2C_{2}\sqrt{\eta}+\frac{2C_{3}}{\sqrt{N}}.$

_Step 2._ We now control the NI error of the averaged pushforward proximal
distributions using $\mathfrak{N}$. In the defining maximum over
$\mu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}}^{N}),\nu^{(N)}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{Y}}^{N})$,
we may restrict to product distributions $\mu^{(N)}=\mu^{\otimes
N},\nu^{(N)}=\nu^{\otimes N}$ so that

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})\right]$
$\displaystyle\geq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\max_{\mu,\nu}-\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu,\nu_{\mathscr{Y}_{j}})-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu})$
$\displaystyle\qquad+\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu_{\mathscr{X}_{j}},\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu})\Bigg{]}$
$\displaystyle\geq\max_{\mu,\nu}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}-\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu,\nu_{\mathscr{Y}_{j}})-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu})$
$\displaystyle\qquad+\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu_{\mathscr{X}_{j}},\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu})\Bigg{]}$
$\displaystyle\geq\max_{\mu,\nu}-\operatorname{\mathcal{L}}(\mu,\mathbb{E}[\nu_{\overline{\mathscr{Y}}_{k}}])-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\lambda\operatorname{\mathrm{KL}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}]\|\rho^{\nu})$
$\displaystyle\qquad+\operatorname{\mathcal{L}}(\mathbb{E}[\mu_{\overline{\mathscr{X}}_{k}}],\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\lambda\operatorname{\mathrm{KL}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}]\|\rho^{\mu})$

by convex-concavity of $\operatorname{\mathcal{L}}$ as well as convexity of KL
divergence, where we have written

$\overline{\Pi}\operatorname{\widehat{\mu}}_{k}:=\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\Pi\operatorname{\widehat{\mu}}_{j}^{(N)},\quad\overline{\Pi}\operatorname{\widehat{\nu}}_{k}:=\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\Pi\operatorname{\widehat{\nu}}_{j}^{(N)}.$

Again by Proposition 3.6, this is further bounded as

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})\right]$
$\displaystyle\geq\max_{\mu,\nu}-\operatorname{\mathcal{L}}_{\lambda}(\mu,\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}])-\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\frac{1}{B_{k}}\int_{\operatorname{\mathcal{Y}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu,\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}])(\nu_{\mathscr{Y}_{j}}-\Pi\operatorname{\widehat{\nu}}_{j}^{(N)})(\mathop{}\\!\mathrm{d}y)\Bigg{]}$
$\displaystyle\qquad+\operatorname{\mathcal{L}}_{\lambda}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}],\nu)+\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\frac{1}{B_{k}}\int_{\operatorname{\mathcal{X}}}\sum_{j=1}^{k}\beta_{j}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}],\nu)(\mu_{\mathscr{X}_{j}}-\Pi\operatorname{\widehat{\mu}}_{j}^{(N)})(\mathop{}\\!\mathrm{d}x)\Bigg{]}$
$\displaystyle\geq\operatorname{\mathrm{NI}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}],\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}])-\left(\frac{(r+1)^{2}}{rk}+O(k^{-2})\right)\cdot
2C_{1}(\eta)-2C_{2}\sqrt{\eta}-\frac{2C_{3}}{\sqrt{N}},$

with the appropriate modification for $r=0$.

_Step 3._ Finally, we convert the above pushforward proximal bounds back to a
Wasserstein distance bound for the expected empirical measures. By Lemma 3.5
and Talagrand’s inequality for the MNE $(\mu^{*},\nu^{*})$,

$\displaystyle
W_{2}^{2}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}],\mu^{*})+W_{2}^{2}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}],\nu^{*})$
$\displaystyle\leq\frac{2}{\smash[b]{\alpha_{\mu}}}\vee\frac{2}{\alpha_{\nu}}\left(\operatorname{\mathrm{KL}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}]\|\mu^{*})+\operatorname{\mathrm{KL}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}]\|\nu^{*})\right)$
$\displaystyle\leq\frac{2}{\smash[b]{\alpha_{\mu}\lambda}}\vee\frac{2}{\alpha_{\nu}\lambda}\operatorname{\mathrm{NI}}(\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}],\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\nu}}_{k}])$
$\displaystyle\leq\frac{2}{\smash[b]{\alpha_{\mu}\lambda}}\vee\frac{2}{\alpha_{\nu}\lambda}\Bigg{[}\left(\frac{(r+1)^{2}}{rk}+O(k^{-2})\right)\cdot\frac{17}{4}C_{1}(\eta)+4C_{2}\sqrt{\eta}+\frac{4C_{3}}{\sqrt{N}}\Bigg{]}.$

Note also by Proposition 3.6 and Lemma C.8 that

$\displaystyle
M_{\mu}W_{1}(\mathbb{E}[\mu_{\overline{\mathscr{X}}_{k}}],\mathbb{E}[\overline{\Pi}\operatorname{\widehat{\mu}}_{k}])$
$\displaystyle=\sup_{\left\lVert
F\right\rVert_{\operatorname{\mathrm{Lip}}}\leq
M_{\mu}}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\int_{\operatorname{\mathcal{X}}}F(\mu_{\mathscr{X}_{j}}-\Pi\operatorname{\widehat{\mu}}_{j}^{(N)})(\mathop{}\\!\mathrm{d}x)\right]$
$\displaystyle\leq\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\left(\frac{r+1}{j}C_{1}(\eta)+C_{2}\sqrt{\eta}+\frac{C_{3}}{\sqrt{N}}\right)+\frac{1}{B_{k}}\sum_{j=2}^{k}\beta_{j}W_{1}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)},\Pi\operatorname{\widehat{\mu}}_{j-1}^{(N)})$
$\displaystyle\leq\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\left(\frac{r+1}{j}C_{1}(\eta)+C_{2}\sqrt{\eta}+\frac{C_{3}}{\sqrt{N}}\right)+\frac{2M_{\mu}}{\alpha_{\mu}\lambda
B_{k}}\sum_{j=2}^{k}\frac{\beta_{j}^{2}}{B_{j}}$
$\displaystyle\leq\left(\frac{(r+1)^{2}}{rk}+O(k^{-2})\right)\cdot\frac{3}{2}C_{1}(\eta)+C_{2}\sqrt{\eta}+\frac{C_{3}}{\sqrt{N}},$

so the square of this term can be ignored. Hence we can conclude that

$W_{1}^{2}(\mathbb{E}[\mu_{\overline{\mathscr{X}}_{k}}],\mu^{*})+W_{1}^{2}(\mathbb{E}[\nu_{\overline{\mathscr{Y}}_{k}}],\nu^{*})\leq\frac{(r+1)^{2}}{rk}\widetilde{C}_{1}(\eta)+\widetilde{C}_{2}\sqrt{\eta}+\frac{\widetilde{C}_{3}}{\sqrt{N}},$

again with the $1+\log k$ modification when $r=0$. ∎

###  C.6 Expected Wasserstein Distance

Theorem 3.7 gives error bounds for the expected distributions
$\mathbb{E}[\mu_{\overline{\mathscr{X}}_{k}}]$ and
$\mathbb{E}[\nu_{\overline{\mathscr{Y}}_{k}}]$. This quantifies a sort of bias
of the MFL-AG outputs, but does not tell us anything about the variance. Can
we similarly bound the expected distance
$\mathbb{E}[W_{1}(\mu_{\overline{\mathscr{X}}_{k}},\mu^{*})+W_{1}(\nu_{\overline{\mathscr{Y}}_{k}},\nu^{*})]$
of the empirical distributions to the MNE? The following fundamental fact
about Wasserstein distance tells us that this is impossible:

######  Theorem C.14 (Rate of convergence of the empirical measure, adapted
from Fournier & Guillin (2015), Theorem 1).

Let $X^{i}$ be independent samples drawn from
$\mu^{i}\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathbb{R}}^{d})$ for
each $i\in[N]$. If $d\geq 3$, the 1-Wasserstein distance between the empirical
measure $\mu_{\mathscr{X}}=\frac{1}{N}\sum_{i=1}^{N}\delta_{X^{i}}$ and the
underlying averaged measure $\mu=\frac{1}{N}\sum_{i=1}^{N}\mu^{i}$ is bounded
in expectation as

$\mathbb{E}[W_{1}(\mu_{\mathscr{X}},\mu)]\leq
C_{W}\sqrt{\mathfrak{m}_{2}(\mu)}\cdot N^{-1/d},$

where $\mathfrak{m}_{2}(\mu)$ is the raw second moment of $\mu$ and $C_{W}$ is
a universal constant. If $d=2$, the rate is $O(N^{-1/2}(\log N)^{2})$; if
$d=1$, the rate is $O(N^{-1/2})$. Furthermore, this rate is tight up to
constants.

###### Proof.

The original theorem only considers i.i.d. samples
$\mu^{1}=\cdots=\mu^{N}=\mu$ and omits the $W_{1}$ case for simplicity, so we
present the necessary modifications.

For a Borel subset $A\subset\operatorname{\mathbb{R}}^{d}$, the quantity
$N\mu_{\mathscr{X}}(A)$ is not distributed as $\text{Binomial}(N,\mu(A))$ but
as a sum of independent $\text{Bernoulli}(\mu^{i}(A))$ random variables.
Nonetheless, we obtain the same bound

$\displaystyle\mathbb{E}[|\mu_{\mathscr{X}}(A)-\mu(A)|]$
$\displaystyle\leq(\mathbb{E}[\mu_{\mathscr{X}}(A)]+\mu(A))\wedge\sqrt{\operatorname{Var}\mu_{\mathscr{X}}(A)}$
$\displaystyle\leq 2\mu(A)\wedge\sqrt{\mu(A)/N}.$

We now repeat the same arguments and substitute $p=1,q=2$ to arrive at the
following inequality,

$\mathbb{E}[W_{1}(\mu_{\mathscr{X}},\mu)]\leq
C\sqrt{\mathfrak{m}_{2}(\mu)}\cdot\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}2^{-m}(2^{-n}\wedge(2^{dm}/N)^{1/2})$

from which point we give explicit computations. Defining

$m_{N}=\left\lceil\frac{\log_{2}N}{d}\right\rceil,\quad
n_{m}=\left\lceil\frac{\log_{2}N-dm}{2}\right\rceil,$

we have for $d\geq 3$ that

$\displaystyle\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}2^{-m}(2^{-n}\wedge(2^{dm}/N)^{1/2})$
$\displaystyle=\sum_{m=0}^{m_{N}-1}2^{-m}n_{m}(2^{dm}/N)^{1/2}+\sum_{m=0}^{m_{N}-1}\sum_{n=n_{m}}^{\infty}2^{-m-n}+\sum_{m=m_{N}}^{\infty}\sum_{n=0}^{\infty}2^{-m-n}$
$\displaystyle\leq\frac{1}{2\sqrt{N}}\sum_{m=0}^{m_{N}-1}(dm_{N}-dm+2)2^{(d/2-1)m}+\sum_{m=0}^{m_{N}-1}2^{1-m-n_{m}}+2^{2-m_{N}}$
$\displaystyle\leq\frac{(2+d)2^{(d/2-1)(m_{N}+1)}}{(2^{d/2}-2)\sqrt{N}}+\frac{2^{2+(d/2-1)m_{N}}}{(2^{d/2}-2)\sqrt{N}}+2^{2-m_{N}}$
$\displaystyle=O(N^{-1/d}).$

When $d=2$, the rate is easily checked to be $N^{-1/2}(\log N)^{2}$. The tight
rate in one dimension is derived using different techniques in Bobkov & Ledoux
(2016), Section 3. ∎

That is, even in the ideal case where chaos does not propagate and the
particles are somehow i.i.d. sampled directly from the true distribution, the
expected Wasserstein distance will always be of order
$N^{-1/d_{\operatorname{\mathcal{X}}}\vee d_{\operatorname{\mathcal{Y}}}}$,
automatically incurring the curse of dimensionality. We emphasize that the
uniform law of large numbers and short-term perturbation methods developed
throughout Section C as well as the presentation of Theorem 3.7 have been
carefully designed to bypass this technicality.

Nevertheless, it is still possible to bound the expected Wasserstein distance
in a similar manner save for the unavoidable
$N^{-1/d_{\operatorname{\mathcal{X}}}\vee d_{\operatorname{\mathcal{Y}}}}$
dependency.222Of course, we may also simply run the algorithm multiple ($M$)
times and take the average of the outputs, which would also bypass the issue
and yield the standard $1/\sqrt{M}$ convergence. We first present a more
direct bound for the proximal gap.

######  Proposition C.15.

The following inequality holds for all $k$,

$\mathbb{E}\left[W_{1}(\mu_{\mathscr{X}_{k}},\Pi\operatorname{\widehat{\mu}}_{k}^{(N)})\right]\leq\frac{r+1}{k}C_{1}^{\prime}(\eta)+C_{2}^{\prime}\sqrt{\eta}+C_{3}^{\prime}N^{-1/d_{\operatorname{\mathcal{X}}}}.$

###### Proof.

The derivations are similar and more straightforward compared to the proof of
Proposition 3.6. We only look at $k\geq 2\ell$ and directly compare
$\mu_{\mathscr{X}_{k}}$ to $\mu_{\widetilde{\mathscr{X}}_{k}}$ using Lemma
C.3, $\mu_{\widetilde{\mathscr{X}}_{k}}$ to the expected modified distribution
using Theorem C.14 (recall that the modified particle trajectories
$\widetilde{X}_{k}^{i}$ are independent when conditioned on
$(\mathscr{X},\mathscr{Y})_{1:k-\ell}$), the expected modified distribution to
the stationary distribution $\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)}$
using Proposition C.6, and $\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)}$
back to $\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}$ using Lemma C.8.

$\displaystyle\mathbb{E}[W_{1}(\mu_{\mathscr{X}_{k}},\Pi\operatorname{\widehat{\mu}}_{k}^{(N)})]$
$\displaystyle\leq\mathbb{E}[W_{2}(\mu_{\mathscr{X}_{k}},\mu_{\widetilde{\mathscr{X}}_{k}})]+\mathbb{E}[W_{1}(\mu_{\widetilde{\mathscr{X}}_{k}},\mathbb{E}[\mu_{\widetilde{\mathscr{X}}_{k}}])]$
$\displaystyle\qquad+\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[W_{2}(\mu_{k}^{i},\Pi\operatorname{\widehat{\mu}}_{k-\ell}^{(N)})]+\sum_{j=0}^{\ell-1}\mathbb{E}[W_{2}(\Pi\operatorname{\widehat{\mu}}_{k-j-1}^{(N)},\Pi\operatorname{\widehat{\mu}}_{k-j}^{(N)})]$
$\displaystyle\leq\frac{r+1}{k-\ell+1}\mathfrak{w}_{\ell}^{\mu}+C_{W}\sqrt{\mathbb{E}[\mathfrak{m}_{2}(\mu_{\widetilde{\mathscr{X}}_{k-\ell}})]}\cdot
N^{-1/d_{\operatorname{\mathcal{X}}}}$
$\displaystyle\qquad+\frac{1}{N}\sum_{i=1}^{N}\sqrt{\frac{4}{\smash[b]{\alpha_{\mu}}}(\mathfrak{K}^{\mu}\lVert
X_{k-\ell}^{i}\rVert^{2}+\mathfrak{L}^{\mu})}+\frac{2M_{\mu}^{2}}{\alpha_{\mu}\lambda}\sum_{j=0}^{\ell-1}\frac{\beta_{k-j}}{B_{k-j}}$
$\displaystyle\leq\frac{r+1}{k-\ell+1}\mathfrak{w}_{\ell}^{\mu}+C_{W}\sqrt{\mathfrak{p}^{\mu}+\mathfrak{s}^{\mu}}\cdot
N^{-1/d_{\operatorname{\mathcal{X}}}}+\sqrt{\frac{4}{\alpha_{\mu}}(\mathfrak{K}^{\mu}(\mathfrak{p}^{\mu}+\mathfrak{s}^{\mu})+\mathfrak{L}^{\mu})}+\frac{2M_{\mu}}{\alpha_{\mu}\lambda}\frac{(r+1)\ell}{k-\ell+1}$
$\displaystyle=\frac{r+1}{k}C_{1}^{\prime}(\eta)+C_{2}^{\prime}\sqrt{\eta}+C_{3}^{\prime}N^{-1/d_{\operatorname{\mathcal{X}}}}.$

∎

We now give the desired bound for the expected Wasserstein distance to the
MNE. Note the effect of dimensionality compared to Theorem 3.7.

######  Theorem C.16 (Variance of discretized MFL-AG).

If $\eta\leq\bar{\eta}$ and $\beta_{k}=k^{r}$ with $r>0$, the MFL-AG discrete
update satisfies for all $K,N$,

$\mathbb{E}[W_{1}(\mu_{\overline{\mathscr{X}}_{K}},\mu^{*})]^{2}+\mathbb{E}[W_{1}(\nu_{\overline{\mathscr{Y}}_{K}},\nu^{*})]^{2}\leq\frac{(r+1)^{2}}{rK}\widetilde{C}_{1}(\eta)+\widetilde{C}_{2}\sqrt{\eta}+\widetilde{C}_{3}N^{-1/d_{\operatorname{\mathcal{X}}}\vee
d_{\operatorname{\mathcal{Y}}}}$

with similar constants as in Proposition 3.6. When $r=0$, the first term is
replaced by $O(\log K/K)$. If $d_{\operatorname{\mathcal{X}}}\vee
d_{\operatorname{\mathcal{Y}}}=2$, the third term is replaced by
$O(N^{-1/2}(\log N)^{2})$; if
$d_{\operatorname{\mathcal{X}}}=d_{\operatorname{\mathcal{Y}}}=1$, by
$O(N^{-1/2})$.

###### Proof.

Note that by convexity,

$\operatorname{\mathcal{L}}(\mu_{\mathscr{X}_{k}},\nu)-\operatorname{\mathcal{L}}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)},\nu)\geq\int_{\operatorname{\mathcal{X}}}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\Pi\operatorname{\widehat{\mu}}_{k}^{(N)},\nu)(\mu_{\mathscr{X}_{k}}-\Pi\operatorname{\widehat{\mu}}_{k}^{(N)})(\mathop{}\\!\mathrm{d}x)\geq-
M_{\mu}W_{1}(\mu_{\mathscr{X}_{k}},\Pi\operatorname{\widehat{\mu}}_{k}^{(N)}).$

We can modify Step 2 of Section C.5 using Proposition C.15 as follows.

$\displaystyle\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\mathfrak{N}(\mu_{\overline{\mathscr{X}}_{k}},\nu_{\overline{\mathscr{Y}}_{k}})\right]$
$\displaystyle\geq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\max_{\mu,\nu}-\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu,\nu_{\mathscr{Y}_{j}})-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu})$
$\displaystyle\qquad+\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu_{\mathscr{X}_{j}},\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu})\Bigg{]}$
$\displaystyle\geq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Bigg{[}\max_{\mu,\nu}-\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\mu,\Pi\operatorname{\widehat{\nu}}_{j}^{(N)})-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\nu}}_{j}^{(N)}\\!\|\rho^{\nu})$
$\displaystyle\qquad+\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\\!\operatorname{\mathcal{L}}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)},\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\frac{\lambda}{B_{k}}\sum_{j=1}^{k}\beta_{j}\operatorname{\mathrm{KL}}(\Pi\operatorname{\widehat{\mu}}_{j}^{(N)}\\!\|\rho^{\mu})$
$\displaystyle\qquad-\frac{M_{\mu}}{B_{k}}\sum_{j=1}^{k}\beta_{j}{W_{1}(\mu_{\mathscr{X}_{k}},\Pi\operatorname{\widehat{\mu}}_{k}^{(N)})}-\frac{M_{\nu}}{B_{k}}\sum_{j=1}^{k}\beta_{j}{W_{1}(\nu_{\mathscr{Y}_{k}},\Pi\operatorname{\widehat{\nu}}_{k}^{(N)})}\Bigg{]}$
$\displaystyle\geq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\Big{[}\max_{\mu,\nu}-\operatorname{\mathcal{L}}(\mu,\overline{\Pi}\operatorname{\widehat{\mu}}_{k})-\lambda\operatorname{\mathrm{KL}}(\mu\|\rho^{\mu})+\lambda\operatorname{\mathrm{KL}}(\overline{\Pi}\operatorname{\widehat{\nu}}_{k}\\!\|\rho^{\nu})$
$\displaystyle\qquad+\operatorname{\mathcal{L}}(\overline{\Pi}\operatorname{\widehat{\mu}}_{k},\nu)-\lambda\operatorname{\mathrm{KL}}(\nu\|\rho^{\nu})+\lambda\operatorname{\mathrm{KL}}(\overline{\Pi}\operatorname{\widehat{\mu}}_{k}\\!\|\rho^{\mu})\Big{]}$
$\displaystyle\qquad-\frac{M_{\mu}}{B_{k}}\sum_{j=1}^{k}\beta_{j}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}[{W_{1}(\mu_{\mathscr{X}_{k}},\Pi\operatorname{\widehat{\mu}}_{k}^{(N)})}]-\frac{M_{\nu}}{B_{k}}\sum_{j=1}^{k}\beta_{j}\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}[{W_{1}(\nu_{\mathscr{Y}_{k}},\Pi\operatorname{\widehat{\nu}}_{k}^{(N)})}]$
$\displaystyle\geq\mathbb{E}_{(\mathscr{X},\mathscr{Y})_{1:k}}\left[\operatorname{\mathrm{NI}}(\overline{\Pi}\operatorname{\widehat{\mu}}_{k},\overline{\Pi}\operatorname{\widehat{\nu}}_{k})\right]-\frac{M_{\mu}}{B_{k}}\sum_{j=1}^{k}\beta_{j}\left(\frac{r+1}{j}C_{1}^{\prime}(\eta)+C_{2}^{\prime}\sqrt{\eta}+C_{3}^{\prime}N^{-1/d_{\operatorname{\mathcal{X}}}\vee
d_{\operatorname{\mathcal{Y}}}}\right).$

Combining with Step 1 and Lemma 3.5 gives that

$\mathbb{E}\left[\operatorname{\mathrm{KL}}(\overline{\Pi}\operatorname{\widehat{\mu}}_{k}\\!\|\mu^{*})+\operatorname{\mathrm{KL}}(\overline{\Pi}\operatorname{\widehat{\nu}}_{k}\\!\|\nu^{*})\right]\leq\frac{(r+1)^{2}}{rk}C_{1}^{\prime\prime}(\eta)+C_{2}^{\prime\prime}\sqrt{\eta}+C_{3}^{\prime\prime}N^{-1/d_{\operatorname{\mathcal{X}}}\vee
d_{\operatorname{\mathcal{Y}}}}.$

Finally, we convert back to a Wasserstein distance bound by invoking
Talagrand’s inequality and Proposition C.15 again:

$\displaystyle\mathbb{E}[W_{1}(\mu_{\overline{\mathscr{X}}_{k}},\mu^{*})]^{2}\leq
2\bigg{(}\frac{1}{B_{k}}\sum_{j=1}^{k}\beta_{j}\mathbb{E}[W_{1}(\mu_{\mathscr{X}_{j}},\Pi\operatorname{\widehat{\mu}}_{j}^{(N)})]\bigg{)}^{2}+\frac{4}{\alpha_{\mu}}\mathbb{E}[\operatorname{\mathrm{KL}}(\overline{\Pi}\operatorname{\widehat{\mu}}_{k}\\!\|\mu^{*})].$

This concludes the proof. ∎

_Remark._ If we assume a higher degree of regularity so that all relevant
distributions have finite fourth moments, say, then Theorem C.14 actually
holds for the 2-Wasserstein metric. Theorem C.16 can then be stated in terms
of the 2-Wasserstein distance to the MNE, guaranteeing us slightly better
control over the error compared to Proposition 3.6 which only allows a $W_{1}$
formulation.

##  Appendix D Convergence Analysis of MFL-ABR

###  D.1 Inner Loop Convergence

The convergence of the decoupled inner loop is a simple consequence of the
convex analysis for single optimization (Nitanda et al., 2022a); we reproduce
the proof here for completeness.

######  Proposition D.1 (Convergence of MFL-ABR inner loop).

Under Assumptions 1 and 3,

$\operatorname{\mathrm{KL}}(\mu_{k,\tau}^{\dagger}\|\operatorname{\widehat{\mu}}_{k})\leq\frac{2C_{\mu}}{\lambda}\exp(-2\alpha\lambda\tau),\quad\operatorname{\mathrm{KL}}(\nu_{k,\tau}^{\dagger}\|\operatorname{\widehat{\nu}}_{k})\leq\frac{2C_{\nu}}{\lambda}\exp(-2\alpha\lambda\tau).$

###### Proof.

For any $0\leq t\leq\tau$, the KL gap converges as

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\operatorname{\mathrm{KL}}(\mu_{k,t}^{\dagger}\|\operatorname{\widehat{\mu}}_{k})$
$\displaystyle=\int_{\operatorname{\mathcal{X}}}\log\frac{\mu_{k,t}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\partial_{t}\mu_{k,t}^{\dagger}(\mathop{}\\!\mathrm{d}x)=\lambda\int_{\operatorname{\mathcal{X}}}\log\frac{\mu_{k,t}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\nabla_{x}\cdot\bigg{(}\mu_{k,t}^{\dagger}\nabla_{x}\log\frac{\mu_{k,t}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\bigg{)}(\mathop{}\\!\mathrm{d}x)$
$\displaystyle=-\lambda\int_{\operatorname{\mathcal{X}}}\bigg{\lVert}\nabla_{x}\log\frac{\mu_{k,t}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\bigg{\rVert}^{2}\mu_{k,t}^{\dagger}(\mathop{}\\!\mathrm{d}x)\leq-2\alpha\lambda\cdot\operatorname{\mathrm{KL}}(\mu_{k,t}^{\dagger}\|\operatorname{\widehat{\mu}}_{k})$

by substituting the Fokker-Planck equation for $\mu_{k,t}^{\dagger}$ and
applying the LSI for $\operatorname{\widehat{\mu}}_{k}$ via Theorem A.6.
Invoking Gronwall’s lemma and Lemma D.2 below for
$\operatorname{\widehat{\mu}}_{k}$ concludes the proof. ∎

The following result gives uniform bounds to control the magnitude of
perturbations.

######  Lemma D.2.

For any $w>0$, define the class

$\mathcal{F}_{w}^{\mu}:=\left\\{\mu\in\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}}):\left\lVert\log\frac{\mu}{\rho^{\mu}}\right\rVert_{\infty}\leq\frac{wC_{\mu}}{\lambda}\right\\}.$

Then under Assumption 3, the distribution
$\operatorname{\widehat{\mu}}_{k}\in\mathcal{F}_{2}^{\mu}$ and
$\mu_{k},\mu_{k,\tau}^{\dagger}\in\mathcal{F}_{4}^{\mu}$.

###### Proof.

For $\operatorname{\widehat{\mu}}_{k}$, the exponential term and the
normalizing integral

$\exp\left(-\frac{1}{\lambda}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{k},\nu_{k})\right),\quad
Z_{k}^{\mu}=\int_{\operatorname{\mathcal{X}}}\rho^{\mu}\exp\left(-\frac{1}{\lambda}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{k},\nu_{k})\right)\mathop{}\\!\mathrm{d}x$

are both bounded by $T^{\mu}/\lambda$, proving the assertion. For
$\mu_{k,\tau}^{\dagger}$, define the density ratio
$h_{t}=\mu_{k,t}^{\dagger}/\operatorname{\widehat{\mu}}_{k}$. The Fokker-
Planck equation for $\mu_{k,t}^{\dagger}$ reads

$\partial_{t}\mu_{k,t}^{\dagger}=\nabla_{x}\cdot\left(\mu_{k,t}^{\dagger}\nabla_{x}\left(\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{k},\nu_{k})+\lambda\nabla_{x}U^{\mu}\right)\right)+\lambda\Delta_{x}\mu_{k,t}^{\dagger}=\lambda\nabla_{x}\cdot\left(\mu_{k,t}^{\dagger}\nabla_{x}\log\frac{\mu_{k,t}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\right),$

so that the parabolic partial differential equation satisfied by $h_{t}$ is
derived as

$\displaystyle\partial_{t}h_{t}$
$\displaystyle=\operatorname{\widehat{\mu}}_{k}^{-1}\partial_{t}\mu_{k,t}^{\dagger}$
$\displaystyle=\lambda\operatorname{\widehat{\mu}}_{k}^{-1}\nabla_{x}\cdot\left(\operatorname{\widehat{\mu}}_{k}h_{t}\nabla_{x}\log
h_{t}\right)$
$\displaystyle=\lambda\nabla_{x}\log\operatorname{\widehat{\mu}}_{k}\cdot\nabla_{x}h_{t}+\lambda\Delta
h_{t}$
$\displaystyle=-\nabla_{x}\left(\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu_{k},\nu_{k})+\lambda\nabla_{x}U^{\mu}\right)\cdot\nabla_{x}h_{t}+\lambda\Delta
h_{t}$ $\displaystyle=\mathbf{L}^{\dagger}h_{t},$

where $\mathbf{L}^{\dagger}$ is the infinitesimal generator for the stochastic
process $X_{t}^{\dagger}$. Hence by the Feynman-Kac formula, we may write for
any $t\in[0,\tau]$

$h_{t}(x)=\mathbb{E}^{x}[h_{0}(X_{t}^{\dagger})]=\mathbb{E}^{x}\left[\frac{\rho^{\mu}}{\operatorname{\widehat{\mu}}_{k}}(X_{t})\right].$

Since
$\lVert\log(\operatorname{\widehat{\mu}}_{k}/\rho^{\mu})\rVert_{\infty}\leq
2C_{\mu}/\lambda$ as discussed above, we infer that $\left\lVert
h_{\tau}\right\rVert\leq 2C_{\mu}/\lambda$ and therefore

$\bigg{\lVert}\log\frac{\mu_{k,\tau}^{\dagger}}{\rho^{\mu}}\bigg{\rVert}_{\infty}\leq\bigg{\lVert}\log\frac{\mu_{k,\tau}^{\dagger}}{\operatorname{\widehat{\mu}}_{k}}\bigg{\rVert}_{\infty}+\bigg{\lVert}\log\frac{\operatorname{\widehat{\mu}}_{k}}{\rho^{\mu}}\bigg{\rVert}_{\infty}\leq\frac{4C_{\mu}}{\lambda},$

i.e. $\mu_{k,\tau}^{\dagger}\in\mathcal{F}_{\infty}^{\mu}$ for all $k$.
Finally, since $\mathcal{F}_{\infty}^{\mu}$ is closed under linear
combinations in $\operatorname{\mathcal{P}_{2}}(\operatorname{\mathcal{X}})$
we conclude that

$\mu_{k}=\beta\mu_{k,\tau}^{\dagger}+\beta(1-\beta)\mu_{k-1,\tau}^{\dagger}+\cdots\beta(1-\beta)^{k}\mu_{0,\tau}^{\dagger}\in\mathcal{F}_{\infty}^{\mu}.$

∎

###  D.2 Proof of Theorem 4.1

We perform one-step analysis of the outer loop by setting for $0\leq s\leq 1$

$\mu(s)=(1-\beta s)\mu_{k}+\beta s\mu_{k,\tau}^{\dagger},\quad\nu(s)=(1-\beta
s)\nu_{k}+\beta s\nu_{k,\tau}^{\dagger},$

so that $\mu(0)=\mu_{k}$, $\mu(1)=\mu_{k+1}$ and $\nu(0)=\nu_{k}$,
$\nu(1)=\nu_{k+1}$. We track the KL divergence to the interpolated proximal
distributions defined as

$\operatorname{\widehat{\mu}}(s)=\frac{\rho^{\mu}}{Z^{\mu}(s)}\exp\left(-\frac{1}{\lambda}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu(s),\nu(s))\right),\quad\operatorname{\widehat{\nu}}(s)=\frac{\rho^{\nu}}{Z^{\nu}(s)}\exp\left(\frac{1}{\lambda}\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\nu}(\mu(s),\nu(s))\right).$

Note that the second order bounds in Assumption 3 immediately imply the
following Lipschitz property in TV distance,

$\left\lVert\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu,\nu)-\frac{\delta\\!\operatorname{\mathcal{L}}}{\delta\mu}(\mu^{\prime},\nu^{\prime})\right\rVert_{\infty}\leq
2C_{\mu\mu}\operatorname{TV}(\mu,\mu^{\prime})+2C_{\mu\nu}\operatorname{TV}(\nu,\nu^{\prime}).$

Similarly to Lascu et al. (2023), Lemma A.2 we can then prove that

$\displaystyle\operatorname{TV}(\operatorname{\widehat{\mu}}_{k},\operatorname{\widehat{\mu}}(s))$
$\displaystyle\leq\frac{1}{2\lambda}\left(\exp\left(\frac{C_{\mu}}{\lambda}\right)+\exp\left(\frac{2C_{\mu}}{\lambda}\right)\right)\left(2C_{\mu\mu}\operatorname{TV}(\mu_{k},\mu(s))+2C_{\mu\nu}\operatorname{TV}(\nu_{k},\nu(s))\right)$
$\displaystyle\leq\beta s\mathfrak{t}^{\mu},$

where we have written

$\mathfrak{t}^{\mu}:=\frac{C_{\mu\mu}+C_{\mu\nu}}{\lambda}\left(\exp\left(\frac{C_{\mu}}{\lambda}\right)+\exp\left(\frac{2C_{\mu}}{\lambda}\right)\right).$

Also, $\operatorname{\widehat{\mu}}(s)\in\mathcal{F}_{2}^{\mu}$ and
$\mu(s)\in\mathcal{F}_{4}^{\mu}$ by Lemma D.2 which implies
$\lVert\log(\mu(s)/\operatorname{\widehat{\mu}}(s))\rVert_{\infty}\leq
6C_{\mu}/\lambda$.

Now the derivative of the KL gap of the max policy for any $0\leq s\leq 1$ is

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}s}\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))=\int_{\operatorname{\mathcal{X}}}\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}\partial_{s}\mu(s)(\mathop{}\\!\mathrm{d}x)-\int_{\operatorname{\mathcal{X}}}\partial_{s}\log\operatorname{\widehat{\mu}}(s)\mu(s)(\mathop{}\\!\mathrm{d}x).$

The first term can be decomposed as

$\displaystyle\int_{\operatorname{\mathcal{X}}}\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}\partial_{s}\mu(s)(\mathop{}\\!\mathrm{d}x)$
$\displaystyle=\beta\int_{\operatorname{\mathcal{X}}}\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}(\mu_{k,\tau}^{\dagger}-\mu_{k})(\mathop{}\\!\mathrm{d}x)$
$\displaystyle=\beta\int_{\operatorname{\mathcal{X}}}\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}(\operatorname{\widehat{\mu}}(s)-\mu(s)+\mu(s)-\mu_{k}+\mu_{k,\tau}^{\dagger}-\operatorname{\widehat{\mu}}_{k}+\operatorname{\widehat{\mu}}_{k}-\operatorname{\widehat{\mu}}(s))(\mathop{}\\!\mathrm{d}x)$
$\displaystyle\leq-\beta\left(\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))+\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))\right)+2\beta^{2}s\left\lVert\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}\right\rVert_{\infty}\operatorname{TV}(\mu_{k,\tau}^{\dagger},\mu_{k})$
$\displaystyle\qquad+\beta\left\lVert\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}\right\rVert_{\infty}\sqrt{2\operatorname{\mathrm{KL}}(\mu_{k,\tau}^{\dagger}\|\operatorname{\widehat{\mu}}_{k})}+2\beta\left\lVert\log\frac{\mu(s)}{\operatorname{\widehat{\mu}}(s)}\right\rVert_{\infty}\operatorname{TV}(\operatorname{\widehat{\mu}}_{k},\operatorname{\widehat{\mu}}(s))$
$\displaystyle\leq-\beta\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))+\frac{12\beta
C_{\mu}}{\lambda}\Bigg{(}2\beta
s(\mathfrak{t}^{\mu}+1)+\sqrt{\frac{C_{\mu}}{\lambda}}\exp(-\alpha\lambda\tau)\Bigg{)}.$

by Proposition D.1. For the second term, we may follow the derivations
presented in Section 3 of Lascu et al. (2023) with minimal modifications to
obtain

$\displaystyle-\int_{\operatorname{\mathcal{X}}}\partial_{s}\log\operatorname{\widehat{\mu}}(s)\mu(s)(\mathop{}\\!\mathrm{d}x)$
$\displaystyle=-\frac{\beta}{\lambda}\iint_{\operatorname{\mathcal{X}}\times\operatorname{\mathcal{X}}}\frac{\delta^{2}\\!\operatorname{\mathcal{L}}}{\delta\mu^{2}}(\mu(s),\nu(s),x,z)(\operatorname{\widehat{\mu}}(s)-\mu(s))(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\mu}}_{k}-\mu_{k})(\mathop{}\\!\mathrm{d}z)$
$\displaystyle\qquad+\frac{\beta}{\lambda}\iint_{\operatorname{\mathcal{X}}\times\operatorname{\mathcal{Y}}}\frac{\delta^{2}\\!\operatorname{\mathcal{L}}}{\delta\mu\delta\nu}(\mu(s),\nu(s),x,w)(\operatorname{\widehat{\mu}}(s)-\mu(s))(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\nu}}_{k}-\nu_{k})(\mathop{}\\!\mathrm{d}w).$

When $s=0$, the first integral is nonpositive due to convexity while the
second integral cancels out when adding with the corresponding term for the KL
gap of the max policy, which completes the argument in Lascu et al. (2023).
Hence the remaining error we must control is

$\displaystyle-\frac{\beta}{\lambda}\iint_{\operatorname{\mathcal{X}}\times\operatorname{\mathcal{X}}}\frac{\delta^{2}\\!\operatorname{\mathcal{L}}}{\delta\mu^{2}}(\mu(s),\nu(s),x,z)(\operatorname{\widehat{\mu}}(s)-\operatorname{\widehat{\mu}}_{k}+\mu_{k}-\mu(s))(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\mu}}_{k}-\mu_{k})(\mathop{}\\!\mathrm{d}z)$
$\displaystyle\qquad+\frac{\beta}{\lambda}\iint_{\operatorname{\mathcal{X}}\times\operatorname{\mathcal{Y}}}\frac{\delta^{2}\\!\operatorname{\mathcal{L}}}{\delta\mu\delta\nu}(\mu(s),\nu(s),x,w)(\operatorname{\widehat{\mu}}(s)-\operatorname{\widehat{\mu}}_{k}+\mu_{k}-\mu(s))(\mathop{}\\!\mathrm{d}x)(\operatorname{\widehat{\nu}}_{k}-\nu_{k})(\mathop{}\\!\mathrm{d}w)$
$\displaystyle\leq\frac{4\beta}{\lambda}(C_{\mu\mu}+C_{\mu\nu})(\operatorname{TV}(\operatorname{\widehat{\mu}}(s),\operatorname{\widehat{\mu}}_{k})+\operatorname{TV}(\mu_{k},\mu(s)))$
$\displaystyle\leq\frac{4\beta^{2}s}{\lambda}(C_{\mu\mu}+C_{\mu\nu})(\mathfrak{t}^{\mu}+1).$

Adding everything up, we obtain

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}s}\left(\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))+\operatorname{\mathrm{KL}}(\nu(s)\|\operatorname{\widehat{\nu}}(s))\right)$
$\displaystyle\leq-\beta\left(\operatorname{\mathrm{KL}}(\mu(s)\|\operatorname{\widehat{\mu}}(s))+\operatorname{\mathrm{KL}}(\nu(s)\|\operatorname{\widehat{\nu}}(s))\right)$
$\displaystyle\qquad+\frac{12\beta
C_{\mu}}{\lambda}\Bigg{(}2\beta(\mathfrak{t}^{\mu}+1)+\sqrt{\frac{C_{\mu}}{\lambda}}\exp(-\alpha\lambda\tau)\Bigg{)}+\frac{4\beta^{2}}{\lambda}(C_{\mu\mu}+C_{\mu\nu})(\mathfrak{t}^{\mu}+1)$
$\displaystyle\qquad+\frac{12\beta
C_{\nu}}{\lambda}\Bigg{(}2\beta(\mathfrak{t}^{\nu}+1)+\sqrt{\frac{C_{\nu}}{\lambda}}\exp(-\alpha\lambda\tau)\Bigg{)}+\frac{4\beta^{2}}{\lambda}(C_{\mu\nu}+C_{\nu\nu})(\mathfrak{t}^{\nu}+1).$

By applying Gronwall’s lemma over $s\in[0,1]$ and iterating over $k$, we
conclude that

$\displaystyle\operatorname{\mathrm{KL}}(\mu_{k}\|\operatorname{\widehat{\mu}}_{k})+\operatorname{\mathrm{KL}}(\nu_{k}\|\operatorname{\widehat{\nu}}_{k})$
$\displaystyle\leq\frac{2(C_{\mu}+C_{\nu})}{\lambda}\exp(-\beta
k)+\frac{12}{\lambda^{\frac{3}{2}}}\left(C_{\mu}^{\frac{3}{2}}+C_{\nu}^{\frac{3}{2}}\right)\exp(-\alpha\lambda\tau)$
$\displaystyle\qquad+\frac{4\beta}{\lambda}\left((6C_{\mu}+C_{\mu\mu}+C_{\mu\nu})(\mathfrak{t}^{\mu}+1)+(6C_{\nu}+C_{\mu\nu}+C_{\nu\nu})(\mathfrak{t}^{\nu}+1)\right).$

Finally, applying Lemma 3.4 of Lascu et al. (2023) yields the suboptimality
bound

$\displaystyle\operatorname{\mathrm{NI}}(\mu_{k},\nu_{k})\leq
2(C_{\mu}+C_{\nu})\exp(-\beta
k)+\frac{12}{\sqrt{\lambda}}\left(C_{\mu}^{\frac{3}{2}}+C_{\nu}^{\frac{3}{2}}\right)\exp(-\alpha\lambda\tau)+C\beta.$

Hence an $\epsilon$-MNE may be obtained in
$k=O(\frac{1}{\epsilon}\log\frac{1}{\epsilon})$ outer loop iterations by
taking $\beta=O(\epsilon)$ and $\tau=O(\log\frac{1}{\epsilon})$. ∎