$B=\\{z_{i}\in\zeta_{0}\\}=\\{\overline{z_{i}}\in\overline{\zeta_{0}}\\}$.

This leads to

$\displaystyle\quad\partial_{i_{1}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}-\partial_{i_{1}}F_{\overline{{\zeta_{0}}}}\cdots\partial_{i_{l}}F_{\overline{{\zeta_{0}}}}$
$\displaystyle=\partial_{i_{1}}F_{\zeta_{0}^{(i)}}(\partial_{i_{2}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}-\partial_{i_{2}}F_{\overline{\zeta_{0}}}\cdots\partial_{i_{l}}F_{\overline{\zeta_{0}}})$
$\displaystyle\quad+\frac{1}{n_{b}}\mathds{1}_{B}\left(\partial_{i_{1}}f(x,z_{i})\partial_{i_{2}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}-\partial_{i_{1}}f(x,\overline{z_{i}})\partial_{i_{2}}F_{\overline{\zeta_{0}}}\cdots\partial_{i_{l}}F_{\overline{\zeta_{0}}}\right),$

where we observe that the terms inside the second bracket are $n$-independent.
If we take the expectation $\mathbb{E}_{\mathbf{z}}$ on these terms, we will
be averaging over $\binom{n}{n_{b}}$ terms, and only
$\frac{n_{b}}{n}\binom{n}{n_{b}}=\binom{n-1}{n_{b}-1}$ of these terms will be
non-zero due to the indicator function $\mathds{1}_{B}$.

Hence if we let $\phi_{\zeta_{0}}$ be any $n$-independent function, we will
have

$\mathbb{E}_{\mathbf{z}}\frac{1}{n_{b}}\mathds{1}_{B}\phi_{\zeta_{0}}=\frac{1}{n_{b}}\frac{1}{n}\sum_{\zeta_{0}\subset\mathbf{z}}\mathds{1}_{B}\phi_{\zeta_{0}}=\frac{1}{n}\widetilde{\phi},$
(6.25)

where $\widetilde{\phi}$ is also $n$-independent.

We can continue taking the expansion to get

$\displaystyle\quad\partial_{i_{1}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}-\partial_{i_{1}}F_{\overline{{\zeta_{0}}}}\cdots\partial_{i_{l}}F_{\overline{{\zeta_{0}}}}$
$\displaystyle=\frac{1}{n_{b}}\mathds{1}_{B}\sum_{j=1}^{l}\partial_{i_{1}}F_{\zeta_{0}^{(i)}}\cdots\partial_{i_{j-1}}F_{\zeta_{0}^{(i)}}\partial_{i_{j}}f(x,z_{i})\partial_{i_{j+1}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}$
$\displaystyle\quad-\frac{1}{n_{b}}\mathds{1}_{B}\sum_{j=1}^{l}\partial_{i_{1}}F_{\zeta_{0}^{(i)}}\cdots\partial_{i_{j-1}}F_{\zeta_{0}^{(i)}}\partial_{i_{j}}f(x,\overline{z_{i}})\partial_{i_{j+1}}F_{\overline{\zeta_{0}}}\cdots\partial_{i_{l}}F_{\overline{\zeta_{0}}}$
$\displaystyle=:\frac{1}{n_{b}}\mathds{1}_{B}\phi^{\alpha}_{\zeta_{0},\overline{\zeta_{0}}}\,,$

where the sum follows from recursively applying the first step to the terms
$\partial_{i_{1}}F_{\zeta_{0}^{(i)}}(\partial_{i_{2}}F_{\zeta_{0}}\cdots\partial_{i_{l}}F_{\zeta_{0}}-\partial_{i_{2}}F_{\overline{\zeta_{0}}}\cdots\partial_{i_{l}}F_{\overline{\zeta_{0}}})$.

Computing the expectation using Equation (6.25), we get that

$\mathbb{E}_{\mathbf{z}}\frac{1}{n_{b}}\mathds{1}_{B}\phi^{\alpha}_{\zeta_{0},\overline{\zeta_{0}}}=\frac{1}{n}\phi^{\alpha}.$

This implies we can define the desired operator as

$\widehat{A}_{N}=\sum_{\alpha}c_{\alpha}\phi^{\alpha}\partial_{\alpha},$

hence completing the proof for the first step.

Step 2. Finding $\widehat{L}_{N}$

We will prove by induction over $N$ that there exists an operator
$\widehat{L}_{N}$ independent of $n$ such that
$\frac{1}{n}\widehat{L}_{N}=L_{N,\mathbf{z}}-L_{N,\overline{\mathbf{z}}}$.

For the $N=1$ case, it follows by computing

$L_{N,\mathbf{z}}-L_{N,\overline{\mathbf{z}}}=L_{\mathbf{z}}-L_{\overline{\mathbf{z}}}=\frac{1}{n}\sum_{j=1}^{d}-\partial_{j}(f(x,z_{i})-f(x,\overline{z_{i}}))\partial_{j}.$

Now assuming the case for $N-1$, we will prove the statement for case $N$. We
first recall the definition of $L_{N}$

$L_{N}=A_{N+1}+\sum_{l=1}^{N}\frac{B_{l}}{l!}\sum_{n_{1}+\cdots
n_{l+1}=n-l}L_{n_{1}}\cdots L_{n_{l}}A_{n_{l+1}+1}\,.$

Since
$A_{N,\mathbf{z}}-A_{N,\overline{\mathbf{z}}}=\frac{1}{n}\widehat{A}_{N}$, it
is sufficient to analyze

$\displaystyle\quad L_{n_{1},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}A_{n_{l+1}+1,\mathbf{z}}-L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}}A_{n_{l+1}+1,\overline{\mathbf{z}}}$
$\displaystyle=(L_{n_{1},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}-L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}})A_{n_{l+1}+1,\mathbf{z}}+L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}}\frac{1}{n}\widehat{A}_{n_{l+1}+1}$

Now observe that since $n_{1},\cdots,n_{l}<N$, we have
$\frac{1}{n}\widehat{L}_{n_{j}}=L_{n_{j},\mathbf{z}}-L_{n_{j},\overline{\mathbf{z}}}$
for all $j=1,\cdots,l$. This allows us to write

$\displaystyle\quad L_{n_{1},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}-L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}}$
$\displaystyle=(L_{n_{1},\mathbf{z}}-L_{n_{1},\overline{\mathbf{z}}})L_{n_{2},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}+L_{n_{1},\overline{\mathbf{z}}}(L_{n_{2},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}-L_{n_{2},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}})$
$\displaystyle=\frac{1}{n}\widehat{L}_{n_{1}}L_{n_{2},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}+L_{n_{1},\overline{\mathbf{z}}}(L_{n_{2},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}-L_{n_{2},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{l}L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{j-1},\overline{\mathbf{z}}}\widehat{L}_{n_{j}}L_{n_{j+1},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}$

Putting it together we have

$\displaystyle\quad\widehat{L}_{N}$
$\displaystyle=\widehat{A}_{N+1}+\sum_{l=1}^{N}\frac{B_{l}}{l!}\sum_{n_{1}+\cdots
n_{l+1}=n-l}\sum_{j=1}^{l}L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{j-1},\overline{\mathbf{z}}}\widehat{L}_{n_{j}}L_{n_{j+1},\mathbf{z}}\cdots
L_{n_{l},\mathbf{z}}A_{n_{l+1}+1,\mathbf{z}}$
$\displaystyle\quad+\sum_{l=1}^{N}\frac{B_{l}}{l!}\sum_{n_{1}+\cdots
n_{l+1}=n-l}L_{n_{1},\overline{\mathbf{z}}}\cdots
L_{n_{l},\overline{\mathbf{z}}}\widehat{A}_{n_{l+1}+1}\,.$

Since the adjoint operation is linear, we can compute the adjoint for each
operator separately. Hence we have the desired result

$\frac{1}{n}\widehat{L}^{*}_{N}=L_{N,\mathbf{z}}-L_{N,\overline{\mathbf{z}}}\,.$

Step 3. Providing the Bound

To complete the bound it is sufficient to observe

$\left\|\phi\partial_{\alpha}((L^{*}_{\ell,\mathbf{z}}-L^{*}_{\ell,\overline{\mathbf{z}}})\mu_{k-\ell,\overline{\mathbf{z}}})q_{\overline{z_{i}}}\right\|_{L^{2}\left(\rho_{\mathbf{z}^{(i)}}\right)}=\frac{1}{n}\left\|\phi\partial_{\alpha}(\widehat{L}_{N}^{*}\mu_{k-\ell,\overline{\mathbf{z}}})q_{\overline{z_{i}}}\right\|_{L^{2}\left(\rho_{\mathbf{z}^{(i)}}\right)}\leq\frac{C}{n},$

where the bound follows from Lemma 6.11.

∎

##  7 Runtime Complexity: Proof of Corollary 2.5

We will restate and prove the result here.

######  Corollary 7.1 (Runtime Complexity).

Suppose $\\{X_{k}\\}_{k\geq 0}$ is any discretization of Langevin diffusion
(1.3) admiting an approximate stationary distribution $\pi^{N}_{\mathbf{z}}$
of the type in Theorem 2.3. Then there exists a constant $C>0$ (depending on
$N$), such that for all

$\epsilon>0\,,\quad
0<\eta<\min\left\\{\frac{2m}{M^{2}}\,,C\epsilon^{1/N}\right\\}\,,\quad
n\geq\frac{C}{\epsilon(1-\eta)}\,,\quad
k\geq\frac{C}{\epsilon^{1/N}}\log\frac{1}{\epsilon}\,,$ (7.1)

we achieve the following expected generalization bound

$\left|\,\mathbb{E}\left[F(X_{k})-F_{\mathbf{z}}(X_{k})\right]\,\right|\leq\epsilon\,.$
(7.2)

###### Proof.

We will start by decomposing the generalization error via an approximation
step with $\pi^{N}_{\mathbf{z}}$

$\displaystyle\left|\,\mathbb{E}\left[F(X_{k})-F_{\mathbf{z}}(X_{k})\right]\,\right|$
(7.3) $\displaystyle\leq$
$\displaystyle\left|\,\mathbb{E}\left[F(X_{k})-\pi^{N}_{\mathbf{z}}(F)\right]\,\right|+\left|\,\mathbb{E}\left[\pi^{N}_{\mathbf{z}}(F)-\pi^{N}_{\mathbf{z}}(F_{\mathbf{z}})\right]\,\right|+\left|\,\mathbb{E}\left[\pi^{N}_{\mathbf{z}}(F_{\mathbf{z}})-F_{\mathbf{z}}(X_{k})\right]\,\right|$
$\displaystyle\leq$ $\displaystyle C(e^{-\lambda
k\eta/2}+\eta^{N})+\frac{C}{n(1-\eta)}\,,$

where we used the results of Theorems 2.3 and 2.4 and absorbed dependence on
$x,F,F_{\mathbf{z}}$ into the constant $C$.

It is then sufficient to confirm all three terms are of order $O(\epsilon)$.
Observe clearly choosing $\eta=O(\epsilon^{1/N})$ and
$n=\Omega(\epsilon^{-1})$ is sufficient for the latter two terms. We will then
observe that

$e^{-\lambda k\eta/2}=O(\epsilon)\iff
k=\Omega(\eta^{-1}\log(\epsilon^{-1}))\,,$ (7.4)

and substituting in the choice of $\eta^{-1}=\epsilon^{-1/N}$ gives us the
desired result.

∎

## References

  * Abdulle et al. (2012) Assyr Abdulle, David Cohen, Gilles Vilmart, and Konstantinos C Zygalakis.  High weak order methods for stochastic differential equations based on modified equations.  _SIAM Journal on Scientific Computing_ , 34(3):A1800–A1823, 2012.
  * Abdulle et al. (2014) Assyr Abdulle, Gilles Vilmart, and Konstantinos C Zygalakis.  High order numerical approximation of the invariant measure of ergodic sdes.  _SIAM Journal on Numerical Analysis_ , 52(4):1600–1622, 2014.
  * Anton (2017) Cristina Anton.  Error expansion for a symplectic scheme for stochastic hamiltonian systems.  In _International Conference on Applied Mathematics, Modeling and Computational Science_ , pages 567–577. Springer, 2017.
  * Anton (2019) Cristina Anton.  Weak backward error analysis for stochastic hamiltonian systems.  _BIT Numerical Mathematics_ , 59(3):613–646, 2019\.
  * Bakry et al. (2013) D. Bakry, I. Gentil, and M. LeDoux.  _Analysis and Geometry of Markov Diffusion Operators_.  Springer, 2013.  ISBN 9783319002286.
  * Bakry et al. (2008) Dominique Bakry, Franck Barthe, Patrick Cattiaux, and Arnaud Guillin.  A simple proof of the poincaré inequality for a large class of probability measures.  _Electronic Communications in Probability [electronic only]_ , 13:60–66, 2008.
  * Block et al. (2020) Adam Block, Youssef Mroueh, Alexander Rakhlin, and Jerret Ross.  Fast mixing of multi-scale langevin dynamics underthe manifold hypothesis.  _arXiv preprint arXiv:2006.11166_ , 2020.
  * Boumal et al. (2016) Nicolas Boumal, Vlad Voroninski, and Afonso Bandeira.  The non-convex burer-monteiro approach works on smooth semidefinite programs.  In _Advances in Neural Information Processing Systems_ , pages 2757–2765, 2016.
  * Bousquet and Elisseeff (2002) Olivier Bousquet and André Elisseeff.  Stability and generalization.  _J. Mach. Learn. Res._ , 2:499–526, March 2002.  ISSN 1532-4435.
  * Burer and Monteiro (2003) Samuel Burer and Renato DC Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Mathematical Programming_ , 95(2):329–357, 2003\.
  * Cao et al. (2020) Yu Cao, Jianfeng Lu, and Lihan Wang.  Complexity of randomized algorithms for underdamped langevin dynamics.  _arXiv preprint arXiv:2003.09906_ , 2020.
  * Cattiaux and Guillin (2021) Patrick Cattiaux and Arnaud Guillin.  Functional inequalities for perturbed measures with applications to log-concave measures and to some bayesian problems, 2021.
  * Chaudhari et al. (2019) Pratik Chaudhari, Anna Choromanska, Stefano Soatto, Yann LeCun, Carlo Baldassi, Christian Borgs, Jennifer Chayes, Levent Sagun, and Riccardo Zecchina.  Entropy-sgd: Biasing gradient descent into wide valleys.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2019(12):124018, 2019.
  * Cheng et al. (2018) Xiang Cheng, Niladri S Chatterji, Yasin Abbasi-Yadkori, Peter L Bartlett, and Michael I Jordan.  Sharp convergence rates for langevin dynamics in the nonconvex setting.  _arXiv preprint arXiv:1805.01648_ , 2018.
  * Dalalyan and Karagulyan (2019) Arnak S Dalalyan and Avetik Karagulyan.  User-friendly guarantees for the langevin monte carlo with inaccurate gradient.  _Stochastic Processes and their Applications_ , 129(12):5278–5311, 2019.
  * Debussche and Faou (2011) A. Debussche and E. Faou.  Weak backward error analysis for SDEs.  _ArXiv e-prints_ , May 2011.
  * Durmus and Moulines (2017) Alain Durmus and Eric Moulines.  Nonasymptotic convergence analysis for the unadjusted langevin algorithm.  _The Annals of Applied Probability_ , 27(3):1551–1587, 2017.
  * Erdogdu et al. (2018) M. A. Erdogdu, L. Mackey, and O. Shamir.  Global Non-convex Optimization with Discretized Diffusions.  _ArXiv e-prints_ , October 2018.
  * Erdogdu and Hosseinzadeh (2020) Murat A Erdogdu and Rasa Hosseinzadeh.  On the convergence of langevin monte carlo: The interplay between tail growth and smoothness.  _arXiv preprint arXiv:2005.13097_ , 2020.
  * Evans (2010) Lawrence C Evans.  Partial differential equations.  2010\.
  * Gelfand and Mitter (1991) Saul B Gelfand and Sanjoy K Mitter.  Recursive stochastic algorithms for global optimization in r^d.  _SIAM Journal on Control and Optimization_ , 29(5):999–1018, 1991.
  * Hairer et al. (2006) Ernst Hairer, Christian Lubich, and Gerhard Wanner.  _Geometric numerical integration: structure-preserving algorithms for ordinary differential equations_ , volume 31.  Springer Science & Business Media, 2006.
  * (23) Wolfram Research, Inc.  Mathematica, Version 12.0.  Champaign, IL, 2019.
  * Jordan et al. (1998) Richard Jordan, David Kinderlehrer, and Felix Otto.  The variational formulation of the fokker–planck equation.  _SIAM journal on mathematical analysis_ , 29(1):1–17, 1998.
  * Kopec (2013) M. Kopec.  Weak backward error analysis for overdamped Langevin processes.  _ArXiv e-prints_ , October 2013.
  * Kopec (2015) Marie Kopec.  Weak backward error analysis for langevin process.  _BIT Numerical Mathematics_ , 55(4):1057–1103, 2015.
  * Kuo (2006) H.H. Kuo.  _Introduction to Stochastic Integration_.  Universitext. Springer New York, 2006.  ISBN 9780387310572.
  * Laurent and Vilmart (2020) Adrien Laurent and Gilles Vilmart.  Exotic aromatic b-series for the study of long time integrators for a class of ergodic sdes.  _Mathematics of Computation_ , 89(321):169–202, 2020.
  * Li and Erdogdu (2020) Mufan Bill Li and Murat A Erdogdu.  Riemannian langevin algorithm for solving semidefinite programs.  _stat_ , 1050:21, 2020.
  * Li et al. (2019) Xuechen Li, Yi Wu, Lester Mackey, and Murat A Erdogdu.  Stochastic runge-kutta accelerates langevin monte carlo and beyond.  In _Advances in Neural Information Processing Systems_ , pages 7748–7760, 2019.
  * Mattingly et al. (2010) J.C. Mattingly, A.M. Stuart, and M. Tretyakov.  Convergence of numerical time-averaging and stationary measures via the poisson equation.  _SIAM Journal of Numerical Analysis_ , 48:552–577, 2010\.
  * Menz and Schlichting (2014) Georg Menz and André Schlichting.  Poincaré and logarithmic sobolev inequalities by decomposition of the energy landscape.  _The Annals of Probability_ , 42(5):1809–1884, 2014.
  * Pardoux and Veretennikov (2001) E. Pardoux and Yu. Veretennikov.  On the poisson equation and diffusion approximation. i.  _Ann. Probab._ , 29(3):1061–1085, 07 2001.  doi: 10.1214/aop/1015345596.
  * Pardoux et al. (2003) E Pardoux, A Yu Veretennikov, et al.  On poisson equation and diffusion approximation 2.  _The Annals of Probability_ , 31(3):1166–1192, 2003.
  * Pardoux and Răşcanu (2014) Etienne Pardoux and Aurel Răşcanu.  Stochastic differential equations.  In _Stochastic Differential Equations, Backward SDEs, Partial Differential Equations_ , pages 135–227. Springer, 2014.
  * Pavliotis (2014) Grigorios A Pavliotis.  _Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations_ , volume 60.  Springer, 2014.
  * Raginsky et al. (2017) Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky.  Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis.  In _Proceedings of the 30th Conference on Learning Theory, COLT 2017, Amsterdam, The Netherlands, 7-10 July 2017_ , pages 1674–1703, 2017\.
  * Shalev-Shwartz et al. (2010) Shai Shalev-Shwartz, Ohad Shamir, Nathan Srebro, and Karthik Sridharan.  Learnability, stability and uniform convergence.  _Journal of Machine Learning Research_ , 11:2635–2670, 2010\.
  * Shardlow (2006) Tony Shardlow.  Modified equations for stochastic differential equations.  _BIT Numerical Mathematics_ , 46(1):111–125, 2006\.
  * Talay and Tubaro (1990) Denis Talay and Luciano Tubaro.  Expansion of the global error for numerical schemes solving stochastic differential equations.  _Stochastic analysis and applications_ , 8(4):483–509, 1990.
  * Vempala and Wibisono (2019) Santosh S Vempala and Andre Wibisono.  Rapid convergence of the unadjusted langevin algorithm: Log-sobolev suffices.  _arXiv preprint arXiv:1903.08568_ , 2019.
  * Villani (2009) Cédric Villani.  Hypocoercivity.  _Memoirs of the American Mathematical Society_ , 202(950), 2009.
  * Wang et al. (2016) Lijin Wang, Jialin Hong, and Liying Sun.  Modified equations for weakly convergent stochastic symplectic schemes via their generating functions.  _BIT Numerical Mathematics_ , 56(3):1131–1162, 2016.
  * Wilkinson (1960) James H Wilkinson.  Error analysis of floating-point computation.  _Numerische Mathematik_ , 2(1):319–340, 1960\.
  * Xu et al. (2017) P. Xu, J. Chen, D. Zou, and Q. Gu.  Global Convergence of Langevin Dynamics Based Algorithms for Nonconvex Optimization.  _ArXiv e-prints_ , July 2017.

##  Appendix A Signed Measure Results

We start this section by mentioning that the approximate stationary measure
constructed using weak backward error analysis can be a signed measure.

######  Proposition A.1 (The Approximate Stationary Measure is Signed).

Consider the following simple Ornstein-Uhlenbeck process in $\mathbb{R}^{d}$

$dX(t)=-X(t)dt+\sqrt{\frac{2}{\beta}}dW(t),$

with $f(x,z)=\frac{1}{2}|x|^{2}$ for all $z\in\mathcal{Z}$. We have that
whenever $\frac{2}{d}<\eta$, the approximate density $\pi^{1}(x)$ is not
always positive, specifically

$\pi^{1}(0)=(1+\eta\mu_{1}(0))\rho(0)<0.$

###### Proof.

We first write down the generator of this process

$L\phi=\sum_{i=1}^{d}-x_{i}\partial_{i}\phi+\frac{1}{\beta}\partial_{ii}\phi,$

and therefore leading to the following operator $A_{2}$,

$\displaystyle A_{2}\phi$
$\displaystyle=\sum_{\mathbf{k}}\sum_{i=1}^{d}A^{\mathbf{k}}_{1}(x)\left[(-x_{i})\partial_{i}+\frac{1}{\beta}\partial_{ii}\right]\partial_{\mathbf{k}}\phi$
$\displaystyle=\sum_{\mathbf{k}}\sum_{i,j=1}^{d}x_{i}x_{j}\partial_{ij}\phi-\frac{2}{\beta}x_{i}\partial_{ijj}\phi+\frac{1}{\beta^{2}}\partial_{iijj}\phi.$

Following the definition of $L_{1}$ (5.3), we can compute

$L_{1}=A_{2}+B_{1}L_{0}A_{1}=A_{2}+B_{1}L^{2},$

where $B_{1}$ is the first Bernoulli number. Now observe that since $L$ is
self-adjoint with respect to $L^{2}(\rho)$, i.e. $L^{\star}=L$, we have

$L_{1}^{\star}1=A_{2}^{\star}+B_{1}(L^{\star})^{2}1=A_{2}^{\star}1.$

Here we will compute explicitly $A_{2}^{\star}1$, and we start by writing

$(A_{2}^{\star}1)\rho=\sum_{i,j=1}^{d}\partial_{ij}(x_{i}x_{j}\rho)+\frac{2}{\beta}\partial_{ijj}(x_{i}\rho)+\frac{1}{\beta^{2}}\partial_{iijj}\rho=:T_{1}+T_{2}+T_{3}.$

We first compute one derivative of $T_{3}$ to match $T_{2}$ to get

$T_{3}=\sum_{i,j=1}^{d}\partial_{ijj}\frac{1}{\beta}(-x_{i}\rho),$

this implies that we can combine $T_{2}$ and $T_{3}$

$T_{2}+T_{3}=\sum_{i,j=1}^{d}\partial_{ijj}\frac{1}{\beta}(x_{i}\rho)=\sum_{i,j=1}^{d}\partial_{ij}\left(\frac{\delta_{ij}}{\beta}\rho-
x_{i}x_{j}\rho\right),$

where $\delta_{ij}$ denotes the Kronecker delta, and we add the $T_{1}$ term
as well to get

$T_{1}+T_{2}+T_{3}=\sum_{i,j=1}^{d}\partial_{ij}\frac{\delta_{ij}}{\beta}\rho=\sum_{i=1}^{d}\partial_{ii}\frac{1}{\beta}\rho=\sum_{i=1}^{d}\partial_{i}\left(-x_{i}\rho\right)=\sum_{i=1}^{d}\beta
x_{i}^{2}\rho-d\rho.$

This implies we have the following PDE for $\mu_{1}$

$L\mu_{1}=-L_{1}^{\star}1\implies\sum_{i=1}^{d}-x_{i}\partial_{i}\mu_{1}+\frac{1}{\beta}\partial_{ii}\mu_{1}=d-\sum_{i=1}^{d}\beta
x_{i}^{2}.$

Since the equation has a unique solution that satisfies the integral
constraint $\int\mu_{1}d\rho=0$ (see Proposition 6.8), we can explicitly guess
the solution

$\mu_{1}=\frac{\beta}{2}|x|^{2}-c,$

where to satisfy the integral constraint, we must have

$c=\int\frac{\beta}{2}|x|^{2}\frac{1}{(2\pi/\beta)^{d/2}}\exp\left(\frac{-\beta|x|^{2}}{2}\right)dx=\frac{d}{2}.$

This implies

$\inf_{x}\mu_{1}=\mu_{1}(0)=-\frac{d}{2}.$

Finally this implies that whenever $\eta>\frac{2}{d}$, we have

$\pi^{1}(0)=\left(1-\eta\frac{d}{2}\right)\rho(0)<0.$

∎

Now that we know our approximate measure $\pi^{N}$ can be signed, we can no
longer define uniform stability with respect to an expectation over the random
algorithm. Instead we naturally extend the definition based on the integral
with respect to the signed measure instead.

######  Definition A.2 (Uniform Stability).

A collection of distributions $\\{\pi_{\mathbf{z}}\\}$ on $\mathbb{R}^{d}$
indexed by $\mathbf{z}\in\mathcal{Z}^{n}$ is said to be $\epsilon$-uniformly
stable if for all $\mathbf{z},\overline{\mathbf{z}}\in\mathcal{Z}^{n}$ with
only one differing coordinate

$\sup_{z\in\mathcal{Z}}\left|\pi_{\mathbf{z}}(f(\cdot,z))-\pi_{\overline{\mathbf{z}}}(f(\cdot,z))\right|\leq\epsilon.$
(A.1)

######  Proposition A.3 (Generalization).

Suppose the collection of distributions $\\{\pi_{\mathbf{z}}\\}$ is
$\epsilon$-uniformly stable, and that for all
$(\mathbf{z},z)\in\mathcal{Z}^{n+1}$, we also have $f(\cdot,z)\in
L^{1}(\pi_{\mathbf{z}})$. Then $\\{\pi_{\mathbf{z}}\\}$ has
$\epsilon$-expected generalization error, or more precisely

$\left|\,\mathbb{E}\left[\pi_{\mathbf{z}}(F_{\mathbf{z}})-\pi_{\mathbf{z}}(F)\right]\,\right|\leq\epsilon\,.$
(A.2)

###### Proof.

it is sufficient to realize that since $f(\cdot,z)$ is integrable with respect
to $\pi_{\mathbf{z}}$, all the usual manipulations are well defined. We will
include the proof for completeness as it is an unintuitive claim.

We will start by denoting $\mathbf{z}=(z_{1},\ldots,z_{n})$,
$\widehat{\mathbf{z}}=(\overline{z_{1}},\ldots,\overline{z_{n}})$, and also
the replaced one data set
$\widehat{\mathbf{z}}^{(i)}=(z_{1},\ldots,z_{i-1},\overline{z_{i}},z_{i+1},\ldots,z_{n})$.

With this we can write

$\displaystyle\mathbb{E}_{\mathbf{z}\sim\mathcal{D}^{n}}\int_{\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}f(x,z_{i})d\pi_{\mathbf{z}}(x)$
$\displaystyle=\mathbb{E}_{(\mathbf{z},\widehat{\mathbf{z}})\sim\mathcal{D}^{2n}}\int_{\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}f(x,\overline{z_{i}})d\pi_{\widehat{\mathbf{z}}^{(i)}}(x)$
$\displaystyle=\mathbb{E}_{(\mathbf{z},\widehat{\mathbf{z}})\sim\mathcal{D}^{2n}}\int_{\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}f(x,\overline{z_{i}})d\pi_{\mathbf{z}}(x)+\delta,$

where we define $\delta$ as

$\delta:=\mathbb{E}_{(\mathbf{z},\widehat{\mathbf{z}})\sim\mathcal{D}^{2n}}\int_{\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}f(x,\overline{z_{i}})d\pi_{\widehat{\mathbf{z}}^{(i)}}(x)-\int_{\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}f(x,\overline{z_{i}})d\pi_{\mathbf{z}}(x).$

Therefore, $\epsilon$-uniform stability implies the difference between the two
integrals is at most $\epsilon$, hence we have the generalization error is
$|\delta|\leq\epsilon$.

∎

We observe that uniform stability is a stronger notion than needed here. In
fact, we will only need a sense of stability in expectation for generalization
in expectation (Shalev-Shwartz et al., 2010). Regardless, we have shown that
if $\pi^{N}_{\mathbf{z}}$ is $\epsilon$-uniformly stable, then we have also
have that $|T_{2}|\leq\epsilon$ as desired.

##  Appendix B Additional Calculations

In this section, we will consider the Ornstein-Uhlenbeck process in
$\mathbb{R}$, or $d=1$

$dX(t)=-X(t)dt+\sqrt{2/\beta}\,dW(t),$

where $W(t)$ is the standard one-dimensional Brownian motion. This corresponds
to the loss function being $F(x)=\frac{1}{2}x^{2}$. This implies we have the
Gibbs distribution is a one-dimensional Gaussian

$\rho(x)=\frac{1}{\sqrt{2\pi/\beta}}\exp\left(-\frac{\beta x^{2}}{2}\right).$

###  B.1 Computing the Relevant Operators

The above definition leads us to the following simple generator

$L\phi=-x\partial\phi+\frac{1}{\beta}\partial^{2}\phi,$

and therefore leading to the operators

$A_{k}=\sum_{\ell=0}^{k}\binom{k}{\ell}(-x)^{\ell}\frac{1}{\beta^{k-\ell}}\partial^{2k-\ell}.$

We will also obtain the operators

$\displaystyle L_{0}$ $\displaystyle=L,$ $\displaystyle L_{1}$
$\displaystyle=A_{2}+B_{1}L_{0}A_{1}=A_{2}+B_{1}L^{2},$ $\displaystyle L_{2}$
$\displaystyle=A_{3}+B_{1}(L_{0}A_{2}+L_{1}A_{1})+B_{2}(L_{0}L_{0}A_{1}).$

Next we will compute several adjoint operators with respect to $L^{2}(\rho)$.

First we recall $L$ is self-adjoint, i.e. $L^{*}=L$.

Then we can compute for all $f,g\in
C^{\infty}_{\operatorname{pol}}(\mathbb{R})$

$\int\partial fg\rho dx=\int-f\partial g\rho-fg\partial\rho dx=\int
f(-\partial+\beta x)g\rho dx,$

which means $(\partial)^{*}=-\partial+\beta x$.

At this point, we will use Mathematica (Inc., ) to compute all adjoint
operators by composition of the first order adjoint operator. For example, we
can compute all of the following using a recursion scheme, where $a\in
C^{\infty}_{\operatorname{pol}}(\mathbb{R})$

$\displaystyle(a\partial)^{*}$ $\displaystyle=-a\partial-\partial a+\beta xa,$
$\displaystyle(\partial^{2})^{*}$ $\displaystyle=\partial^{2}-2\beta
x\partial+\beta^{2}x^{2}-\beta,$ $\displaystyle(a\partial^{2})^{*}$
$\displaystyle=a\partial^{2}+2(\partial a-\beta
xa)\partial+(\partial^{2}a-\beta a-2\beta x\partial a+\beta^{2}x^{2}a),$
$\displaystyle(\partial^{3})^{*}$ $\displaystyle=-\partial^{3}+3\beta
x\partial^{2}+(3\beta-3\beta^{2}x^{2})\partial+(\beta^{3}x^{3}-3\beta^{2}x).$

Next we will compute $A_{2}^{*},A_{3}^{*}$

$\displaystyle A_{2}^{*}$
$\displaystyle=\left(x^{2}\partial^{2}-\frac{2x}{\beta}\partial^{3}+\frac{1}{\beta^{2}}\partial^{4}\right)^{*}$
$\displaystyle=\frac{1}{\beta^{*}}\partial^{4}-\frac{2x}{\beta}\partial^{3}+x^{2}\partial^{2}-2x\partial+(\beta
x-1).$ $\displaystyle A_{3}^{*}$
$\displaystyle=\left(-x^{3}\partial^{3}+\frac{3x^{2}}{\beta}\partial^{4}+\frac{-3x}{\beta^{2}}\partial^{5}+\frac{1}{\beta^{3}}\partial^{6}\right)^{*}$
$\displaystyle=\frac{1}{\beta^{3}}\partial^{6}-\frac{3x}{\beta^{2}}\partial^{5}+\frac{3x^{2}}{\beta}\partial^{4}-\left(x^{3}+\frac{6x}{\beta}\right)\partial^{3}+\left(9x^{2}-\frac{9}{\beta}\right)\partial^{2}-(9x-3\beta
x^{3})\partial.$

Finally, we can conclude

$\displaystyle L_{1}^{*}$ $\displaystyle=A_{2}^{*}+B_{1}L^{2},$ $\displaystyle
L_{2}^{*}$ $\displaystyle=A_{3}^{*}+B_{1}(A_{2}^{*}L+LL_{1}^{*})+B_{2}L^{3}.$

###  B.2 Solving the PDEs

To compute the approximations for
$\pi^{2}=\rho(1+\eta\mu_{1}+\eta^{2}\mu_{2})$, we will need to solve the
following two PDEs

$L\mu_{1}=-L_{1}^{*}1,\quad L\mu_{2}=-L_{1}^{*}\mu_{1}-L_{2}^{*}1,$

with the natural constraint that $\int\mu_{1}d\rho=\int\mu_{2}d\rho=0$.

####  B.2.1 First PDE

From the previous section, we have that

$-L_{1}^{*}1=-(A_{2}^{*}+B_{1}L^{2})1=-A_{2}^{*}1=1-\beta x^{2}.$

This means we need to first solve

$-x\partial\mu_{1}+\frac{1}{\beta}\partial^{2}\mu_{1}=1-\beta x^{2}.$

Here we can guess the solution $\mu_{1}=\frac{\beta}{2}x^{2}+c$, which gives
us

$-x(\beta x)+\frac{1}{\beta}\beta=1-\beta x^{2},$

and the constant is naturally

$c=-\int\frac{\beta}{2}x^{2}d\rho=-\frac{\beta}{2}\frac{1}{\beta}=\frac{-1}{2},$

where we used the fact that the Gaussian variance is $1/\beta$.

####  B.2.2 Second PDE

We need to compute $-L_{1}^{*}\mu_{1}=-(A_{2}^{*}+B_{1}L^{2})\mu_{1}$,
starting with the first term $A_{2}^{*}\mu_{1}$

$\left(\frac{1}{\beta^{2}}\partial^{4}-\frac{2x}{\beta}\partial^{3}+x^{2}\partial^{2}-\frac{1}{\beta}\partial^{2}+\beta
x^{2}-1\right)\left(\frac{\beta}{2}x^{2}-\frac{1}{2}\right)=\frac{1}{2}(\beta
x^{2}-1).$

###  B.3 Toy Example

In this subsection, we consider a one-dimensional problem of optimizing a loss
function $F(x)=\frac{1}{2}x^{2}$ with deterministic gradient. The Langevin
update can be written as

$X_{k+1}=X_{k}-\eta X_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k},\quad X_{0}=x.$

where $\xi_{k}\sim N(0,1)$ are i.i.d. random variables.

In this case, the corresponding continuous time Langevin process is described
by the following SDE

$dX(t)=-X(t)dt+\sqrt{\frac{2}{\beta}}dW(t),\quad X(0)=x,$

where $W(t)$ is a standard Brownian motion.

Observe we can compute the time-marginal distributions of both processes.
Let’s start by observing that $X_{k}$ is a sum of i.i.d. normal random
variables, hence it must also be a normal random variable. We simply need to
estimate the mean and variance.

To compute the parameters, we start by rewriting the update rule as

$\displaystyle X_{k+1}$ $\displaystyle=X_{k}-\eta
X_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k}$
$\displaystyle=(1-\eta)X_{k}+\sqrt{2\eta}{\beta}\xi_{k}$
$\displaystyle=(1-\eta)\left((1-\eta)X_{k-1}+\sqrt{\frac{2\eta}{\beta}}\xi_{k-1}\right)+\sqrt{2\eta}{\beta}\xi_{k}$
$\displaystyle=(1-\eta)^{k+1}X_{0}+\sum_{\ell=0}^{k}\sqrt{2\eta}{\beta}(1-\eta)^{\ell}\xi_{\ell},$

this implies $\mathbb{E}X_{k}=(1-\eta)^{k+1}x$, and we also have the variance

$\mathbb{E}(X_{k}-\mathbb{E}X_{k})^{2}=\sum_{\ell=0}^{k}\frac{2\eta}{\beta}(1-\eta)^{2\ell}=\frac{2}{\beta(2-\eta)}\left[1-(1-\eta)^{2(k+1)}\right].$

It is well known that the continuous time SDE is the Ornstein-Uhlenbeck
process, which has the following solution (Kuo, 2006)

$X(t)=X(0)e^{-t}+\int_{0}^{t}e^{-(t-s)}dW(s),$

which is a Gaussian process. The mean and variance can also be computed as

$\displaystyle\mathbb{E}X(t)$ $\displaystyle=xe^{-t},$
$\displaystyle\mathbb{E}(X(t)-\mathbb{E}X(t))^{2}$
$\displaystyle=\int_{0}^{t}e^{-2(t-s)}ds=\frac{1}{\beta}(1-e^{-2t}).$

In the actual plot of Figure 2 (b), we used parameters $\beta=20,\eta=0.5$.
Additionally, the stationary distribution of the discrete Langevin algorithm
is found by a kernel density smoothing for $10,000$ steps of simulation of
$X_{k}$, and we smoothed using a normal density with parameter $\sigma=0.1$.

###  B.4 Plot for Figure 2 (a)

In the plot, we consider the forward Euler discretization of the ODE

$\frac{dy}{dt}=y^{2},\quad y(0)=1,$

which gives the following update

$y_{k+1}=y_{k}+\eta y_{k}^{2},\quad y_{0}=1.$

The true solution is given by

$y(t)=\frac{1}{1-t},$

which the modified equation is given by

$\frac{dv}{dt}=v^{2}-\eta v^{3}+\eta^{2}\frac{3}{2}v^{4}+\cdots,\quad v(0)=1.$

To generate the plot in Figure 2 (a), we used a step size $\eta=\frac{1}{6}$
for the forward Euler solver, and for the modified equation we approximated
the solution using forward Euler with step size $\eta^{\prime}=\frac{1}{600}$
so it is sufficiently close to the true solution.


# A study on the plasticity of neural networks

Tudor Berariu  Wojciech Czarnecki  Soham De  Jorg Bornschein  Samuel Smith
Razvan Pascanu  Claudia Clopath

###### Abstract

One aim shared by multiple settings, such as continual learning or transfer
learning, is to leverage previously acquired knowledge to converge faster on
the current task. Usually this is done through fine-tuning, where an implicit
assumption is that the network maintains its plasticity, meaning that the
performance it can reach on any given task is not affected negatively by
previously seen tasks. It has been observed recently that a pretrained model
on data from the same distribution as the one it is fine-tuned on might not
reach the same generalisation as a freshly initialised one. We build and
extend this observation, providing a hypothesis for the mechanics behind it.
We discuss the implication of losing plasticity for continual learning which
heavily relies on optimising pretrained models.

Machine Learning, ICML, Continual Learning, Generalisation

Machine Learning, ICML

##  1 Introduction

Continual learning is concerned with training on non-stationary data. In a
practical description, an agent learns a sequence of tasks, being restricted
to interact with only one at a time. There are several desiderata for a
successful continual learning algorithm. First, agents should not forget
previously acquired knowledge, unless capacity becomes an issue or
contradicting facts arrive. Second, such an algorithm should be able to
exploit structural similarity between tasks and exhibit accelerated learning.
Third, backward transfer should be possible whenever new knowledge helps
generalisation on previously learnt tasks. Fourth, successful continual
learning relies on an enduring capacity to acquire new knowledge, therefore
learning now should not impede performance on future tasks.

In this work we focus on plasticity, namely the ability of the model to keep
learning. There are different nuances of not being able to learn. A neural
network might lose the capacity to minimise the training loss for a new task.
For example, PackNet (Mallya & Lazebnik, 2017) eventually gets to a point
where all neurons are frozen and learning is not possible anymore. In the same
fashion, accumulating constraints in EWC (Kirkpatrick et al., 2017) might lead
to a strongly regularised objective that does not allow for the new task’s
loss to be minimised. Alternatively, learning might become less data
efficient, referred to as negative forward transfer, an effect often noticed
for regularisation based continual learning approaches. In such a situation
one might still be able to reduce training error to $0$ and obtain full
performance on the new task, is just learning is considerably slower. Lastly,
a third meaning and the one we are concerned with, is that while training
error can be reduced to zero, and irrespective to how fast the model learns,
the optimisation might lead to a poor minimum which achieves lower
generalisation performance.

We define as the generalisation gap the difference in performance a pretrained
model can obtain — e.g. one that had learnt a few tasks already — versus a
freshly initialised one, without constraining the number of updates. Note that
this is similar to the notion of intransigence proposed by (Chaudhry et al.,
2018), but instead making the comparison against a model trained only on the
new data, rather than a multi-task solution. Our focus is on understanding if
a generalisation gap exists, whether it is positive or negative. The transfer
learning dogma used to indicate a positive gap, arguing that pretraining on
large sets of data provides a good initialisation for a related target task.
Recently, (He et al., 2019) showed that this does not necessarily hold,
reporting state-of-the-art results with randomly initialised models, although
at worse sample complexity. (Ash & Adams, 2019) considered an extreme transfer
scenario, where an agent is pretrained on data from the same distribution as
the target task, and reported a negative generalisation gap. See Figure 1 for
our reproduction of this finding. We build on this result in this work, trying
to further expand the empirical evidence on which factor affect the
generalisation gap and take a first step towards understanding its root
causes.

Igl et al. extend the observation that data non-stationarity affects
asymptotic generalisation performance to reinforcement learning scenarios.
Although both Ash & Adams, and Igl et al. propose solutions to close the
generalisation gap, the reasons for its occurrence in the first place are
still unclear. We argue that it can have a considerable impact on how we
approach continual learning, and one should track to what extend it affects
the algorithms we have.

Figure 1: Our reproduction of the core experiment performed by (Ash & Adams,
2019). A ResNet-18 model is pretrained on half of the CIFAR 10 training data,
and then tuned on the full training set. It generalises worse than the model
trained from scratch.

##  2 Generalisation gap - Experiments

In this section we present our analysis of the generalisation gap, and detail
a series of experiments we argue are indicative of the aggravation of the
phenomenon in continual learning. We ask a couple of questions and provide
empirical evidence to answer them, such as: How much pretraining is too much?
Is this negative effect additive when pretraining consists of several stages?
Is there a way to leverage the pretrained parameters for faster tuning?

We start with the same setup as in (Ash & Adams, 2019) training deep residual
networks (He et al., 2016) to classify the CIFAR 10 data set. We use the
average test accuracy in the last 100 training epochs (see the green box in
Figure 1) to compare different setups. We mention below the relevant details
for each experiment, and we offer a full description of the empirical setup in
Appendix A.

#### Does the optimization algorithm affect the gap?

Different optimizers have particular advantages in escaping sub-optimal
regions. We reproduced the warm start experiment for a couple of different
optimisers using constant learning rates: Adam, RMSprop, SGD, and SGD with
momentum (Figures 2, 8, 9, and 10). Ash & Adams reported similar results. The
fact that the generalisation gap manifests in all cases supports the
observation that it is rather a problem with the quality of the local minima
than one of finding appropriate descending trajectories.

Figure 2: Average test accuracy in the last 100 epochs of tuning after
pretraining the model for different numbers of epochs using Adam with a
constant learning rate ($10^{-3}$) for both phases.

#### How many pretraining steps are needed to produce a generalisation gap?

We investigated how much does the gap depend on a large number of optimisation
steps in the pretraining stage. Would early stopping close the gap? In our
experiments, although a larger number of pretraining steps hurt generalisation
more in the tuning stage, just a few passes through the data are enough to
observe a gap (5-10 for Adam, which is even before reaching 100% accuracy). In
conclusion the gap is there before any reason to do early stopping. This is
consistent across all considered optimisers (Figures 2, 8, 9, and 10). A
similar experiment was reported in (Ash & Adams, 2019).

#### Is the gap still there when data distribution slides smoothly?

We tested whether a smooth transition from the pretraining subset to the full
training set would remove the generalisation gap. But, as Figure 3 shows, the
generalisation gap manifests even for a small number of epochs with biased
sampling. This might have profound implications in reinforcement learning
where the data distribution changes slowly during training, as the policy
collecting the data changes.

Figure 3: Models trained in a single stage where each example is individually
sampled with probability $p=1-\gamma^{50n/N}$ from the full training data, and
with probability $1-p$ from the pretrain set ($n$ is the current step, while
$N$ represents the total number of steps – the equivalent of 500 epochs). A
few more details in Section A.2.

Concluding from the last two sets of experiments, a transient bias in the data
distribution significantly impacts generalisation performance.

#### Do multiple pretrain stages and/or class ordering matter?

Continual learning is concerned with possibly unlimited changes in the data
distribution. It is natural to ask whether the loss in generalisation
performance observed as a consequence of a single pretraining stage is
aggravated when data is incrementally added in more steps. In order to answer
this question we divided the data set in multiple splits training the model in
stages. We show in Figure 4 that the final generalisation performance (the
test accuracy achieved in the last stage training on the full training set)
degrades with the number of splits.

To get even closer to the usual continual learning setup, we considered splits
of the training set having some level of class imbalance, therefore exhibiting
larger differences between the data distributions considered at consecutive
stages (Figure 4, right). We tested for splits ranging from class partitions
where each stage would bring data from one or more new classes to the uniform
subsampling from the training set considered so far (see Section A.3 for the
detailed methodology used to split the data set). We noticed that higher
discrepancies between training stages lead to worse generalisation.

Figure 4: The same model (ResNet18) was trained in multiple stages. All but
the last pretrain stages consisted of a number of steps proportional with the
number of examples and sufficient to reach 100% accuracy on train. Right: New
data for a particular stage has a ratio of examples drawn uniformly from the
training set, and the rest from classes designated for that stage (see Section
A.3 for details).

We argue that this observations point out a core difficulty of continual
learning. When saving data from the past is feasible, retraining models seems
a better strategy than using pretrained models.

Figure 5: Average performance on the test set for residual networks of various
depths and widths. See Section A.4 for details on the models’ architectures.

#### How do model width and/or the depth change the generalisation gap?

We investigated whether increasing the capacity of the model helps recovering
the generalisation performance of a randomly initialised model. We show in
Figure 5 that even for very deep and wide models there is a significant gap
between pretrained models and those trained from random initialisations.

Figure 6: Performance of models for which a subset of layers were reset after
pretraining. 1 represents the first convolution, 2-5 are the four modules, and
6 is the fully connected output layer.

#### Which pretrained parameters should be kept for tuning not to have a gap?

Knowing that tuning the whole model leads to poor generalisation performance
we ask what the best strategy is for taking advantage from the pretrained
model? We conducted a series of experiments in which we re-sample the
parameters of some layers from the same distribution used at initialisation.
In our tests with ResNet-18 models on CIFAR 10, resetting just a small subset
of the layers is not enough to fully recover the gap. Our experiments,
summarised in Figure 6, indicate that in order to close the gap the top part
of the model must be reinitialised. Therefore it might be advantageous to keep
the first $k$ layers (as it is usually done in transfer learning), but in our
experiments, $k$ is quite small. Moreover, it seems that there is no advantage
in terms of training speed to keep the pretrained layers (see Section A.5 for
details.)

##  3 A possible account for the generalisation gap: Two Phases of Learning

One plausible hypothesis for the occurrence of the generalisation gap stems
from the flat versus sharp minima view on generalisation (Hochreiter &
Schmidhuber, 1997). Precisely, local minima which exhibit low curvature, and
wide basins of attraction generalise better than sharp ones. This could be
motivated from an information theoretic perspective: flat minima require less
precision to be described (the minimum description length argument made by
(Hochreiter & Schmidhuber, 1997)); or by thinking about the stability around
that point: flat minima are affected less by perturbations in the inputs.
Although there is no formal definition for flatness, previous works proposed
quantities such as the largest eigenvalue of the Hessian (Keskar et al.,
2017), or the local entropy (Chaudhari et al., 2019) to gauge it.

If optimisation were to follow the gradient flow (the infinitesimally precise
path determined by the gradient), then the minimum it converges to would be
determined by the initial random initialisation. In practice optimisation
diverges from that path. This is due to the noise induced by the randomness in
the mini-batch approximation of the loss function, and by the amplitude of the
update step. As a consequence the training dynamics allegedly traverse two
phases. In the early exploration phase, the parameters “bounce” form the
vicinity of a critical point to another until they land in the basin of
attraction of a minimum that is wide enough to trap the optimisation. This
also implies that larger learning rates used in the initial stage of training
are important for generalisation, which is consistent with findings reported
in the literature (e.g. (Li et al., 2019)).

Once the parameters get stuck in the basin of attraction of some minumum,
training goes into a refinement phase, where parameters converge to the said
critical point. In this phase optimisation follows the gradient flow.

Several works bring supporting evidence for the _two phases of learning_
hypothesis. (Achille et al., 2019) also identifies two phases by analysing the
information stored in the weights. Their observations are consistent with the
sharp versus flat minima view as the the Fisher Information Matrix used to
measure connectivity is also indicative about curvature. (Golatkar et al.,
2019) reveal that regularisation has an impact only in the initial phase,
while (Gur-Ari et al., 2018) show that after an early regime gradients reside
in a small subspace that remains constant during training. More relevant works
are mention in Section B.

Building on this hypothesis, it is natural to ask whether the generalisation
gap can be explained by how pretraining affects the exploration phase of
learning.

Note that the amount of exploration that learning is able to do in this
initial stage is proportional to several factors, among which the more
important ones are the learning rate and the inherent noise in the updates.
The role of the learning rate is self-explanatory, it can be seen as scaling
the amount of noise. As for the noise itself, there are multiple sources: the
inherent noise in the data (labelling imprecisions, noise in the observations,
irrelevant features), noise induced by the optimiser of choice (e.g. SGD
introduces noise by relying on mini-batches), noise induced by the data
augmentation procedure, etc. The noise profile of gradients and parameters
updates is important, and known to be non-Gaussian (Simsekli et al., 2019),
hence it is harder to replicate in practice.

We make the conjecture that, _given that everything else stays the same, the
pretraining stage can considerably reduce the amount of noise in the gradients
during the tuning stage, which leads to weaker exploration and convergence to
narrower minima, inducing the generalization gap._

In particular, in the case of discriminative learning which has been the focus
of this work, estimators tend to quickly become robust to many directions of
variation in the data which are not relevant for the classification task. For
example, the model starts ignoring non-discriminative features such as
background patterns early in training. In fact this property is of great
importance to the recent success of neural networks. Many architectural
advances are more efficient in doing so leading to more robust models with
better generalisation properties (e.g. the translation-invariant convolutions
compared with fully-connected linear layers). Data augmentation plays a
similar role.

However, these irrelevant features do potentially play a role in the initial
stage of learning as a source of noise, forcing the optimisation process to
focus on wider minima. When the model is pretrained, it becomes insensitive to
some of these irrelevant features (or some easy to discern sources of noise).
While it is true that when moving to the tuning stage the problem changes —
and hence the loss surface is different and there is likely no relationship
between the two loss surfaces and their critical points (e.g. this is the case
in a continual learning setting) — the model will still be insensitive to some
direction of change (either in the input space or in the latent space). Even
if those direction of variations are not relevant for the new task, it does
mean that in the early stage of tuning there will be considerably less noise,
and hence potentially less exploration. This means that optimisation in the
tuning stage will converge to a narrower minimum, which will generalise less,
leading to the generalization gap.

To test this hypothesis we check whether increasing the learning rate – which
would magnify the remaining noise in the parameter updates– during the tuning
stage helps. Figure 7 shows that a 10x larger learning rate reduces the gap
substantially. The rest of the performance gap could be explained by the fact
that a high constant learning rate does not benefit from the _refinement_
phase. While further empirical evidence is required to validate this
hypothesis, this result is encouraging. Under the assumption that this is the
cause for the gap, one question to be asked is why does the pretraining phase
reduce gradient noise for the tuning phase. The answer might rest in two
observations: 1) the strong data overlap between the two stages, and 2) neural
networks tend to filter early on information in the input that is not
discriminative. If the non-discriminative dimensions are already filtered out
by pretraining, the tuning stage might not benefit from the noise they would
induce.

Figure 7: In this figure we show the performance reached by tuned models for
three specific constant learning rates used during pretraining (the $0.0018$
is the one that generalizes best). Each point is the average of 9 seeds. The
horizontal lines represent the performance achieved by randomly initialised
models. Therefore each distance between a circle and the horizontal lines
represents the gap for a particular pair of learning rates.

##  4 Conclusions

We build on (Ash & Adams, 2019) and study the generalisation gap induced by
pretraining the model on the same data distribution. We extend the original
results, by looking at robustness of this gap to smooth transition between
data distributions, multiple stages of pretraining, model size or resetting
parts of the pretrained model. We take a first step towards understanding this
phenomenon by asking whether it is related to the two phases of learning
hypothesis.

The existence of this generalisation gap suggest that continual learning might
be hurt by using compact models that get finetuned on multiple tasks. We argue
that tracking the generalisation gap represents a new facet of forward
transfer that has not generally been measured or tracked in the literature.

## References

  * Achille et al. (2019) Achille, A., Rovere, M., and Soatto, S.  Critical learning periods in deep networks.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=BkeStsCcKQ.
  * Ash & Adams (2019) Ash, J. T. and Adams, R. P.  On the difficulty of warm-starting neural network training.  _arXiv preprint arXiv:1910.08475_ , 2019.
  * Chaudhari et al. (2019) Chaudhari, P., Choromanska, A., Soatto, S., LeCun, Y., Baldassi, C., Borgs, C., Chayes, J., Sagun, L., and Zecchina, R.  Entropy-sgd: Biasing gradient descent into wide valleys.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2019(12):124018, 2019.
  * Chaudhry et al. (2018) Chaudhry, A., Dokania, P. K., Ajanthan, T., and Torr, P. H. S.  Riemannian walk for incremental learning: Understanding forgetting and intransigence.  In _ECCV_ , 2018.
  * Ghorbani et al. (2019) Ghorbani, B., Krishnan, S., and Xiao, Y.  An investigation into neural net optimization via hessian eigenvalue density.  _arXiv preprint arXiv:1901.10159_ , 2019.
  * Golatkar et al. (2019) Golatkar, A. S., Achille, A., and Soatto, S.  Time matters in regularizing deep networks: Weight decay and data augmentation affect early learning dynamics, matter little near convergence.  In _Advances in Neural Information Processing Systems_ , pp. 10677–10687, 2019.
  * Gur-Ari et al. (2018) Gur-Ari, G., Roberts, D. A., and Dyer, E.  Gradient descent happens in a tiny subspace.  _arXiv preprint arXiv:1812.04754_ , 2018.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2016.
  * He et al. (2019) He, K., Girshick, R., and Dollár, P.  Rethinking imagenet pre-training.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 4918–4927, 2019.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Flat minima.  _Neural Computation_ , 9(1):1–42, 1997.
  * Igl et al. (2021) Igl, M., Farquhar, G., Luketina, J., Böhmer, W., and Whiteson, S.  Transient non- stationarity and generalisation in deep reinforcement learning.  In _Proceedings of the International Conference on Learning Representations_. OpenReview, 2021.
  * Jastrzebski et al. (2020) Jastrzebski, S., Szymczak, M., Fort, S., Arpit, D., Tabor, J., Cho, K., and Geras, K.  The break-even point on optimization trajectories of deep neural networks.  _arXiv preprint arXiv:2002.09572_ , 2020.
  * Keskar et al. (2017) Keskar, N. S., Mudigere, D., Nocedal, J., Smelyanskiy, M., and Tang, P. T. P.  On large-batch training for deep learning: Generalization gap and sharp minima.  In _ICLR 2017 : International Conference on Learning Representations 2017_ , 2017.
  * Kirkpatrick et al. (2017) Kirkpatrick, J. N., Pascanu, R., Rabinowitz, N. C., Veness, J., and et. al.  Overcoming catastrophic forgetting in neural networks.  _Proceedings of the National Academy of Sciences of the United States of America_ , 114 13:3521–3526, 2017.
  * Li et al. (2019) Li, Y., Wei, C., and Ma, T.  Towards explaining the regularization effect of initial large learning rate in training neural networks.  In _Advances in Neural Information Processing Systems_ , pp. 11669–11680, 2019.
  * Mallya & Lazebnik (2017) Mallya, A. and Lazebnik, S.  Packnet: Adding multiple tasks to a single network by iterative pruning.  _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 7765–7773, 2017.
  * Simsekli et al. (2019) Simsekli, U., Sagun, L., and Gurbuzbalaban, M.  A tail-index analysis of stochastic gradient noise in deep neural networks.  In _International Conference on Machine Learning_ , pp. 5827–5837. PMLR, 2019.

##  Appendix A Experiments on CIFAR-10

We detail below the experimental details for the observations made in Section
2. If it’s not explicitly mention otherwise, a ResNet-18 model was trained on
CIFAR 10 with batches of 128 examples, using an Adam optimiser with a constant
learning rate of 0.001. The optimiser’s statistics were reset between the warm
up and the tuning phase. In the first stage (during warm up) the model was
trained for 350 epochs on half of data. In the second stage the model is
trained for 500 epochs on all training data. We report here the average test
performance over all seeds during the last 100 training epochs. In all plots
the error bars measure standard deviation.

Due to space constraints we don’t show learning curves, but if it’s not
otherwise specified it’s implied that accuracy is 100% on the data used in
training for both stages as in Figure 1.

###  A.1 Different optimisers

In addition to Figure 2 in Section 2 we show here how generalisation
performance varies with the number of pretraining epochs on half of data for
three additional optimisers: RMSprop in Figure 8, Stochastic Gradient Descent
(SGD) in Figure 9, and SGD with a constant momentum (0.9) in Figure 10.

Figure 8: Final performance after warming up the model for different numbers
of epochs using RMSProp with a constant learning rate for both phases. Figure
9: Final performance after warming up the model for different numbers of
epochs using SGD with a constant learning rate for both phases. Figure 10:
Final performance after warming up the model for different numbers of epochs
using SGD with a constant learning rate and momentum for both phases.

###  A.2 Smooth transition between distributions.

In the experiments presented in Figure 3 we trained models in a single stage
of 500 epochs. In this case we called an epoch a sequence of 390 update steps
(which is the equivalent of 1 pass through the training data with a batch size
of 128). Note that such an epoch is not a permutation of the data. Each
example from each batch is individually sampled with probability $p$ from the
full training set, and with probability $1-p$ from the pretrainig set.

###  A.3 Class imbalance in the multiple stage setup

Given (i) a training set $\mathcal{D}$ with examples from a set of classes
$\mathcal{C}$, (ii) a number of pretraining stages $n$ and (iii) a number
$0\leq r\leq 1$ (the ratio of data from all classes) we constructed $n$
subsets of $\mathcal{D}$: $\left\\{\mathcal{D}_{i}\right\\}_{1\leq i\leq n}$
to be used for optimisation during the $n$ pretraining stages. In doing this
we applied the following methodology:

  1. 1.

We created a partition of the all classes:
$\left\\{\mathcal{C}_{1},\ldots\mathcal{C}_{n+1}\right\\}$ such that
$\mathcal{C}_{i}\cap\mathcal{C}_{j}=\emptyset,\forall i,j$, and
$\mathcal{C}=\bigcup_{i=1}^{n+1}\mathcal{C}_{i}$.

  2. 2.

We randomly split the full training set $\mathcal{D}$ in two:
$\mathcal{D}_{c}$, $\mathcal{D}_{u}$ such that
$\frac{|\mathcal{D}_{u}|}{|\mathcal{D}|}=r$. (of course:
$\mathcal{D}_{u}\cap\mathcal{D}_{c}=\emptyset$, and
$\mathcal{D}_{u}\cup\mathcal{D}_{c}=\mathcal{D}$).

  3. 3.

We partition $\mathcal{D}_{u}$ into $n+1$ subsets:
$\left\\{D_{u,1},\ldots\mathcal{D}_{u,n+1}\right\\}$.

  4. 4.

We now define the data sets used to optimise the model in each stage (and
considering $\mathcal{D}_{0}=\emptyset$):

$\mathcal{D}_{i}=\mathcal{D}_{i-1}\cup\mathcal{D}_{u,i}\cup\left\\{\left(x,c\right)\in\mathcal{D}_{c}\mid
c\in\mathcal{C}_{i}\right\\}$

The $n+1$-th dataset $\mathcal{D}_{n+1}\equiv\mathcal{D}$ corresponds to the
final tuning phase on the full data set.

###  A.4 Residual networks of various depths and widths

In the experiment with residual neural networks of different widths and depths
we changed the architecture of ResNet-18 (He et al., 2016) as follows.

Apart from the first convolution and the fully connected layer at the output,
ResNet-18 consists of four modules, each made up of $d=2$ residual blocks with
the same number of output channels. Each module doubles the number of channels
and halves the height and the width of the feature maps. The first module
receives a volume with $w=64$ channels, the second operates on $2w=128$, and
so on.

In our experiments we uniformly changed the depth $d$ of the four modules,
and/or scaled the number of channels in all modules ($w,2w,4w,8w$).

Note that this is not the standard way in which people design deeper residual
architectures such as ResNet-32, or ResNet-55. Deeper ResNets increase the
depth of the modules differently, and use bottleneck blocks to avoid an
explosion in the number of parameters.

###  A.5 Resetting the layers of the model

In the experiments presented in Figure 6 from Section 2 we reset subsets of
the model’s parameters. We reset entire modules referring with 1 to the first
convolution, with numbers from 2 to 5 to the four modules (each consisting of
2 residual blocks), and naming 6 the last fully connected layer.

As Figure 6 shows, resetting the last 4, or 5 modules seems to recover the
original performance of a model trained from random parameters. Therefore we
asked whether keeping the pretrained parameters of the first 1 or 2 modules
comes with any advantage in terms of training speed. As Figure 11 shows, in
our setup there seems to be no benefit from preserving parameters from the
pretraining stage.

Figure 11: Here we show the learning curves for three models pretrained for
350 epochs on half of data. For two of them we keep the first 1 or 2 modules,
reinitialise the rest and tune for 500 epochs.

##  Appendix B Supporting evidence for the two phases of learning hyothesis

A couple of works identify critical differences between the early stage and
the late stage of training, offering empirical evidence for the two phases of
learning hypothesis.

(Achille et al., 2019) identifies an initial memorisation phase when data
information is absorbed into the network’s weights, followed by a
reorganisation stage where unimportant connections are pruned and information
decreases while being redistributed among layers for efficiency. Achille et
al. used the Fisher Information Matrix to approximate the amount of
information stored in the weights. The FIM is also a curvature matrix,
therefore the observed regimes support the view that learning changes basins
of attraction of different minima until it lands in one with low curvature,
corresponding to a flat minimum. Achille et al. also point out that if data
statistics change after the initial phase, the network would remain trapped in
the valley the memorisation phase guided it into.

(Golatkar et al., 2019) empirically shows that regularisation has an impact on
final generalisation performance only in the early stages of training.
Applying weight decay or data augmentation only after this initial phase, or
stopping regularisation after that point would not affect generalisation. The
experiments using data augmentation later in training offer additional
evidence for the generalisation gap – if one thinks about that stage as tuning
on more data from the same distribution.

(Gur-Ari et al., 2018) shows that after an early training stage the gradients
reside in a small subspace that remains constant for the rest of training.
This reiterates the importance of the data used in the first steps of
training.

(Li et al., 2019) shows that in overparametrized networks the volume of good
minima dominates the volume of poor minima and underlines the importance of a
high learning rate to land in the basin of attraction of a well generalising
minimum of the loss function. (Jastrzebski et al., 2020) extends the
observation with the importance of using batches to induce the noise needed to
escape poorly generalising minima. Precisely, Jastrzebski et al. point out
that the ratio between learning rate and batch size determines the flatness of
the minumum.

(Ghorbani et al., 2019) computes the full spectrum of the Hessian showing that
Jastrzebski et al.’s claims about smaller learning rates guiding the network
into sharper minima doesn’t hold empirically. A possible explanation is that
the network is already trapped around some minima, and the slow learning rate
just reaches an even flatter region closer to the critical point.