Mihaela Claudia Rosca

On discretisation drift and smoothness regularisation in neural network training
Department of Computer Science

The deep learning recipe of casting real-world problems as mathematical optimisation and tackling the optimisation by training deep neural networks using gradient-based optimisation has undoubtedly proven to be a fruitful one. The understanding behind why deep learning works, however, has lagged behind its practical significance. We aim to make steps towards an improved understanding of deep learning with a focus on optimisation and model regularisation.
We start by investigating gradient descent (GD), a discrete-time algorithm at the basis of most popular deep learning optimisation algorithms. Understanding the dynamics of GD has been hindered by the presence of discretisation drift, the numerical integration error between GD and its often studied continuous-time counterpart, the negative gradient flow (NGF).
To add to the toolkit available to study GD, we derive novel continuous-time flows that account for discretisation drift. Unlike the NGF, these new flows can be used to describe learning rate specific behaviours of GD, such as training instabilities observed in supervised learning and two-player games.
We then translate insights from continuous time into mitigation strategies for unstable GD dynamics, by constructing novel learning rate schedules and regularisers that do not require additional hyperparameters.
Like optimisation, smoothness regularisation is another pillar of deep learning's success with wide use in supervised learning and generative modelling.
Despite their individual significance, the interactions between smoothness regularisation and optimisation have yet to be explored.
We find that smoothness regularisation affects optimisation across multiple deep learning domains, and that incorporating smoothness regularisation in reinforcement learning leads to a performance boost that can be recovered using adaptions to  optimisation methods.
We end by showing
that optimisation can also affect smoothness, as discretisation drift can act as an implicit smoothness regulariser in neural network training.

Our work focuses on understanding and improving components of deep learning systems. As such, it provides two avenues for impact: first, by providing insights that can be useful to other researchers in the field and second, by helping downstream applications that use such systems.

Inside deep learning academic research, this work provides a set of novel tools and insights concerning two main building blocks of deep learning systems, namely optimisation and models. Our insights include novel continuous-time flows to analyse optimisation dynamics, as well as a rethinking of the effect of smoothness constraints imposed on deep learning models. We have accompanied theoretical results with experimental frameworks for analysis and validation, including approaches that can stabilise training and reduce hyperparameter sensitivity and thus computational costs.
Future theoretical research directions based on our line of work include specialising our theoretical results to specific classes of neural networks, further expanding the family of available continuous-time flows capturing optimisation dynamics, and finding new beneficial or detrimental implicit regularisation forces in deep learning optimisation.  Promising avenues for practical impact include new methods for training  deep learning models, and finding new approaches for model selection. Throughout this work, we provide insights across multiple deep learning domains, including supervised learning, two-player games, and reinforcement learning.  The work included here has been presented at multiple conferences and published in journals.

While we do not tackle specific applications directly, this work can have impact outside the research domain by improving the stability of deep learning models used in applications across multiple domains; examples include image classification and generation, as well as game play using reinforcement learning.

One of the biggest opportunities I have been provided with throughout my PhD studies was to have not one, but two incredible supervisors. My primary supervisor, Marc Deisenroth, has been an wonderful source of guidance and strength. Marc, you have raised my bar for rigour in thinking and writing, and have shown me that technical excellence does not come in spite of kindness, but due to kindness; for that, I will be eternally grateful. My second supervisor, Shakir Mohamed, has been my champion and supporter for many years. Shakir, working with you has been the privilege of a life-time; your breadth of technical knowledge is outstanding, and your careful and deliberate application of it through a meticulous choice of impactful research directions should serve as a guide to us all.

I have also been lucky to have an amazing set of collaborators throughout my PhD: Andriy Mnih, Arthur Gretton, Benoit Dherin, Chongli Qin, Claudia Clopath, David G.T. Barrett, Florin Gogianu, Lucian Busoniu, Michael Figurnov, Michael Munn, Razvan Pascanu, Theophane Weber, Tudor Berariu, and Yan Wu. Special thanks to Benoit Dherin, for many evening chats on everything to do with deep learning optimisation and regularisation; to Yan Wu, for a long term fruitful collaboration; to Theophane Weber for many chats on reinforcement learning; to Tudor Berariu and Florin Gogianu for being so much fun while doing research; and to Arthur Gretton for guidance early in my PhD.

This work was done as part of the DeepMind-UCL PhD program, and I am grateful to those organising it and championing this program. Thanks also Frederic Besse, and to Sarah Hodkinson and Claudia Pope for their support.

During my PhD I was lucky to write a book chapter on Implicit Generative models with  Balaji Lakshminarayanan and Shakir Mohamed, for  Kevin P. Murphy's second edition of `Machine Learning: a Probabilistic Perspective'. Beyond allowing me to share my enthusiasm for this area of machine learning, a great clarity of thinking followed the writing process, which helped fledge some of the new ideas presented in this work. I thus must thank Kevin, Balaji, and Shakir for this opportunity.

Thanks also go to my examination committee, Ferenc Huszár and Patrick Rebeschini, for their time and an incredible and thoughtful discussion on my work.

I also have many reasons to be grateful to those who have supported me personally during this time and beyond. I am here because my parents have prioritised my education and encouraged me throughout the way. I also had the most loving of grandmothers, and I am sure she would have loved to see me pursue this degree.

I have been fortunate to have a great set of friends, especially in the though times of coronavirus lockdowns. Niklas, you have been there for me since the times of the Imperial labs, and have always had faith in me; there are no words that can express how grateful I am. Fabio, you are a consistent, reliable, and wise friend; you are also raising my favourite humans, who are a great source of happiness for me. Many thanks also go to Ada, Dhruva, David,  George, James, Ira, Sascha, Siqi, Slemi, Victor, and many other friends, for the time spent together and the great memories.

CHAPTER: INTRODUCTION

Why has deep learning been so successful in recent decades at solving a myriad of problems, such as scientific discovery and applications [Jumper et al., 2021, Ravuri et al., 2021, Fawzi et al., 2022], language and image generation [Ramesh et al., 2022, Brown et al., 2020, Hoffmann et al., 2022, Sauer et al., 2023], and super level human game play in intricate games [Silver et al., 2017, Silver et al., 2017, Berner et al., 2019]?
A common answer to the question of why deep learning has been so successful is that it scales well with increases in amount of data and computational resources, both of which have drastically increased in recent years.
But this answer raises further questions, since many of the reasons why deep learning scales well have yet to be fully uncovered. Particularly, we still do not know why the recipe behind deep learning—using simple gradient-based algorithms to optimise deep neural networks—is so effective on large datasets.
One valid response to this observation is to state that it might not matter why deep learning works, but that it does and to continue exploring its benefits on a wide range of domains, as well as continue the trend of scaling up data, architectures, and compute.
This approach works, can lead to progress, and has been widely used in other sciences. This is not the approach we take here, however.
Throughout this thesis, we take the point of view that understanding why a system works is not only desirable, but necessary in order to ensure targeted research progress, by exploring the areas most likely to lead to an acceleration of new methods, increased resource efficiency, and safe deployment.
In taking this point of view, we follow the footsteps of a growing body of work aiming to uncover, examine, and understand deep learning systems [Arora et al., 2018, Barrett and Dherin, 2021, Smith et al., 2021, Arora et al., 2018, Mescheder et al., 2017, Nagarajan and Kolter, 2019, Nagarajan and Kolter, 2017, Keskar et al., 2017, Dinh et al., 2017, Jiang et al., 2019, Smith et al., 2018].

There are many avenues worthy of study inside the deep learning domain, ranging from investigating why certain model architectures work well [Li et al., 2018, Bjorck et al., 2018, Brock et al., 2021, Ding et al., 2022], to studying why certain types of generative models or learning principles outperform others [Sriperumbudur et al., 2009, Arjovsky et al., 2017, Fedus et al., 2018].
Here, we choose to focus on the study of optimisation and model regularisation in the form of smoothness regularisation. We are motivated by their generality and applicability across deep learning domains from supervised learning, probabilistic modelling, and reinforcement learning; their significant impact on training performance and test set and out of distribution generalisation [Hoffman et al., 2019, Miyato et al., 2018, Radford et al., 2015]; and the opportunity for analysis due to the vast gap between their aforementioned practical performance and the understanding of their mechanisms. Unanswered questions include: What causes instability in deep learning? How can instability be mitigated? Why is there not more instability? What are the implicit regularisation effects of popular optimisers? What are the effects of smoothness regularisation? What are the interactions between model regularisation and optimisation in deep learning? Which optimisation and regularisation effects are domain specific and which transfer across deep learning domains? These are some of the questions we will tackle in this thesis.

A motivating example. Generative Adversarial Networks (GANs) [Goodfellow et al., 2014] are a type of generative model that has led to many image generation breakthroughs [Brock et al., 2018, Sauer et al., 2023, Karras et al., 2019]. Many GAN variants have been proposed, some based on minimising different distributional divergences or distances [Nowozin et al., 2016, Arjovsky et al., 2017, Gulrajani et al., 2017, Bińkowski et al., 2018, Mao et al., 2017] while others use convergence analysis around a Nash equilibrium [Mescheder et al., 2017, Nagarajan and Kolter, 2017]; each of these variants are  theoretically appealing and have different properties worth studying. Yet, the biggest successes in GANs have consistently come from changes to optimisation and regularisation: the choice of optimiser [Radford et al., 2015], large batch sizes [Brock et al., 2018], and incorporating smoothness regularisation [Miyato et al., 2018] are the key ingredients of all successful GANs, while the loss functions used are often those closest to the original formulation [Brock et al., 2018, Miyato et al., 2018, Zhang et al., 2019, Karras et al., 2021]. Similar optimisation and regularisation techniques have also allowed the expansion of GANs from continuous input domains, such as images, to discrete input domains, such as text [de Masson d'Autume et al., 2019]. The GAN example showcases what opportunities are available by understanding and improving optimisation and smoothness regularisation, opportunities we aim to explore here.

§ OPTIMISATION IN DEEP LEARNING

Machine learning is the science of casting real world problems such as prediction, generation, and action into optimisation problems.
Thus, many machine learning problems can be formulated as

\begin{align}
\min_{\vtheta} E(\vtheta),
\end{align}

where $E$ is a loss function suitably chosen for the underlying task and $\vtheta \in \mathbb{R}^D$. For many machine learning tasks, $E$ is an expected value of an unknown data distribution $p^*(\vx)$; $E$ gets estimated  given a training dataset $\mathcal{D}$ of unbiased samples from $p^*(\vx)$: ${E(\vtheta) = \mathbb{E}_{p^*(\vx)} E(\vtheta; \vx) \approx \frac{1}{|\mathcal{D}|} \sum_{\vx_i \in \mathcal{D}} E(\vtheta; \vx_i)}$.
A distinction is made between the model and the objective function components of $E(\vtheta; \vx)$, with the model $f(\vtheta; \vx)$ depending on the parameters $\vtheta$, and the objective function  depending only on the model's output.
In deep learning, $f$ is a neural network, $\vtheta$ are its parameters, and $E$ is non-convex. Due to the high dimensionality of $\vtheta$ or large dataset size, many optimisation algorithms are too costly or inefficient to be used in deep learning, as they scale unfavourably with parameter or dataset size. Luckily, however, we can use the compositional structure of neural networks to compute the gradient $\nabla_{\vtheta} E \in \mathbb{R}^{D}$ using the backpropagation algorithm [Werbos, 1982, Rumelhart et al., 1986, LeCun et al., 1989]. From there, we can use first-order gradient-based iterative algorithms, such as gradient descent and its variants [45, Robbins and Monro, 1951, Kingma and Ba, 2015, Tieleman and Hinton, 2012].

Studying the behaviour of optimisation algorithms can be challenging, even for algorithms as seemingly simple as gradient descent. Typical questions of study include: Will the algorithm converge to a local minima? Does the algorithm get stuck in saddle points?  When does the algorithm fail to converge and when it does converge, how quickly does it converge? The complexity of these issues gets compounded when studying the optimisation of deep neural networks with millions to billions of parameters that do not satisfy commonly used assumptions such as convexity; this often creates a discrepancy between what theory analyses and what is used in practice.
Fundamentally, the deep learning optimisation community is interested in answering the questions: Why do first-order optimisers work so well in training neural networks
and How can we improve optimisation in deep learning?
Progress in this area has been made recently both through theoretical and empirical means of studying neural network optimisation dynamics [Du et al., 2019, Zou et al., 2020, Du et al., 2017, Cohen et al., 2021, Keskar et al., 2017, Lewkowycz et al., 2020].
Novel insights challenge commonly held assumptions, such as the belief that local minima and saddle points are a major challenge for neural network optimisation due to high parameter dimensionality [Du et al., 2019, Zou et al., 2020, Du et al., 2017, Elkabetz and Cohen, 2021].
There has also been a growing body of work showing the importance of the Hessian $\nabla_{\vtheta}^2 E$ in supervised learning optimisation, and specifically on the connection between the learning rate and the largest Hessian eigenvalue and observed instabilities in training [Cohen et al., 2021, Keskar et al., 2017, Lewkowycz et al., 2020]; we will examine the importance of the Hessian in a new light in Chapter <ref>.

Many analyses of discrete-time methods, such as gradient descent, take a continuous-time approach [Glendinning, 1994, Saxe et al., 2014, Nagarajan and Kolter, 2017, Lampinen and Ganguli, 2019, Arora et al., 2018, Advani et al., 2020, Elkabetz and Cohen, 2021, Vardi and Shamir, 2021, França et al., 2020, Balduzzi et al., 2018] as continuous-time proofs tend to be easier to construct [May, 1976, Elkabetz and Cohen, 2021].
The downside of taking a continuous-time and not discrete-time approach to optimisation analysis is that one does not directly analyse the update of interest; the discrepancy between the discrete-time trajectory and its continuous counterpart, which we will call discretisation drift, can lead to conclusions that do not transfer from continuous-time analysis to discrete-time algorithms [Yaida, 2018, Liu et al., 2021].
Discretisation drift is still not greatly understood and has only recently come into attention in deep learning [Kunin et al., 2021, Barrett and Dherin, 2021, Smith et al., 2021].
The incorporation of techniques from the numerical integration community that construct continuous-time flows accounting for discretisation drift
has shed light on the implicit regularisation effects of gradient descent and highlighted the importance of the learning rate in generalisation [Barrett and Dherin, 2021, Smith et al., 2021]. We will use the same numerical integration techniques in Chapter <ref>, where we find a new continuous-time flow that, unlike existing flows, can capture instabilities induced by gradient descent.
We then use insights from our continuous-time analysis to devise an automatic learning rate schedule that can trade-off training stability and generalisation performance.

[Single objective.]

[Zero-sum games.]

[Overview of optimisation in deep learning.]Overview of optimisation in deep learning. Many optimisation problems in deep learning can be seen as either single-objective problems, where all parameters are updated to minimise the same loss, or multi objective problems where where different sets of parameters are updated to minimise different objectives; when the number of objectives is two, these are called two-player games. A special case of two-player games are adversarial games such as zero-sum games, where the objective of one player is to minimise a function $E$, while the objective of the other player is to maximise the same function. This adversarial structure can lead to challenges in optimisation, with approaches that lead to converge to local minima in single objective optimisation fig:single_objective_intro resulting in cyclic trajectories in zero sum games fig:two_player_games_intro.

While many studies focus on optimisation dynamics in the single-objective case, such as supervised learning, strides have also been made towards understanding optimisation in two-player games [Mescheder et al., 2017, Balduzzi et al., 2018, Schäfer and Anandkumar, 2019, Nagarajan and Kolter, 2017, Qin et al., 2020]. Unlike in the single-objective setting, in two-player games not all parameters minimise the same objective, as each player has its own set of parameters, here denoted by $\vphi$ and $\vtheta$, and its own objective:

\begin{align}
&\min_{\vphi} E_{\vphi}(\vphi, \vtheta) \\
&\min_{\vtheta} E_{\vtheta}(\vphi, \vtheta).
\end{align}

Two-player games can be adversarial, where the two players have opposing objectives, and even zero sum, where one player's gain is another's loss: ${E_{\vphi}(\vtheta, \vphi) = - E_{\vtheta}(\vtheta, \vphi)}$.
We visualise a simple two-player zero-sum game in Figure <ref>, showing that adversarial dynamics can make it challenging to reach a local equilibrium compared to a single-objective counterpart, shown in Figure <ref>.

Interest in the analysis of two-player game optimisation has been fuelled by interest in GANs [Goodfellow et al., 2014], a generative model trained via an adversarial two player game.
While GANs have been a successful generative model, they have a reputation of being notoriously hard to train, and  this has served as a motivation for finding approaches to stabilise and improve the dynamics of adversarial two-player games.
Analysing the behaviour and convergence of two-player games has led to progress in understanding the sources of instability and divergence in games and has led to many a regulariser that can stabilise training [Qin et al., 2020, Mescheder et al., 2018, Nagarajan and Kolter, 2017, Balduzzi et al., 2018, Wang et al., 2019, Mazumdar et al., 2019]. We add to this body of work by further investigating sources of instability in two-player games, coming from discretisation drift rather than the original continuous-time system, in Chapter <ref>. By modelling discretisation drift, we devise explicit regularisers that stabilise training without requiring any additional hyperparameter sweep.

§ SMOOTHNESS W.R.T. INPUTS

Optimisation provides an approach to approximate an unknown decision surface required for prediction, generation, or action from an available dataset. That is, given an optimal predictor $f^*$, optimisation provides a recipe to find a parametrised function $f(\cdot;\vtheta)$ close to $f^*$ on the training dataset. But there are often many such optimisation solutions, each with different properties outside the training data.
How useful this approximation is in the underlying application depends on how well it performs on new data.
This is an essential desideratum of machine learning algorithms, and it is relevant for generalisation, continual learning [Parisi et al., 2019], and transfer learning [Zhuang et al., 2020].   Since we know from the no-free-lunch theorem [Shalev-Shwartz and Ben-David, 2014] that no machine learning algorithm will perform better than all other algorithms for any data distribution and prediction problem, to provide guarantees beyond the training data requires making assumptions about the data distribution, the structure of $f^*$, or both.
One approach to encode beliefs about $f^*$ into $f(\cdot;{\vtheta})$ is through what is often referred to in machine learning as an inductive bias [Mitchell, 1980][When we talk about smoothness with respect to inputs, we often write the predictor as a function of data, given parameters. When discussing optimisation and the focus is on parameters, we write the reverse for clarity.].

In accordance with Occam's razor [MacKay, 1991], we often want to encode a simplicity inductive bias: we are searching for the least complex functions that can explain the data. One approach to formalising simplicity is by measuring the effect changes in function input have on the function's output. If small changes in function input do not result in large changes in function output, we call that function smooth. A smoothness inductive bias encodes a preference for learning smooth functions: between two functions that fit the data equally well, the smoothest one should be preferred. We show an intuitive example highlighting the importance of smoothness in Figure <ref>. Figure <ref> exemplifies how too much smoothness can hurt training performance by reducing model capacity, while Figure <ref> shows how the lack of smoothness can lead to overfitting.

[Very smooth fit,
hurts training data fit.]

[Desired level of smoothness.]

[Too smooth fit,
hurts generalisation.]

[The importance of smoothness on the model fit.]The importance of smoothness with respect to inputs on the model fit: too much smoothness can hurt capacity fig:intro_too_smooth, while not enough smoothness can hurt generalisation fig:intro_too_unsmooth.

The formal definition of smoothness is often taken to be Lipschitz smoothness or related measures.
A function $f(\cdot; \vtheta): \mathcal{X} \rightarrow \mathcal{Y}$ is $K$-Lipschitz  if

\begin{align}
\norm{f(\vx_1; \vtheta) - f(\vx_2; \vtheta)}_{\mathcal{Y}} \le K \norm{\vx_1-\vx_2}_{\mathcal{X}}  \hspace{3em} \forall \vx_1, \vx_2 \in \mathcal{X},
\end{align}

where $K$ is known as the Lipschitz constant of function $f(\cdot; \vtheta)$. Throughout this manuscript, when we refer to the smoothness of a neural network, we refer to the smoothness w.r.t. data $\vx$, and not w.r.t. parameters $\vtheta$.

For many a machine learning task smoothness with respect to inputs constitutes a reasonable prior;  this is can be easily illustrated in the image domain, where the prediction or action that ought to be taken does not change when a few pixels in the inputs change: an image of a dog still encodes a dog if a few pixels are modified. Thus, it does not come as a surprise that despite having vastly different motivations and implementations,  many existing deep learning regularisation methods target a smoothness inductive bias, including $L_2$ regularisation, dropout [Srivastava et al., 2014], early stopping.
We will discuss in detail how these regularisation approaches are connected to smoothness in Chapter <ref>, formalise these connections in Chapter <ref>, as well as link smoothness to important phenomena in deep learning, such as double descent [Belkin, 2021, Belkin et al., 2019, Nakkiran et al., 2019].

As the importance of smoothness is beginning to be recognised by the deep learning community, methods that directly target smoothness regularisation have been developed [Miyato et al., 2018, Yoshida and Miyato, 2017, Hoffman et al., 2019]. These regularisation methods have been incorporated in supervised learning to help generalisation [Hoffman et al., 2019, Bartlett et al., 2017] and adversarial robustness [Sokolić et al., 2017, Novak et al., 2018, Cisse et al., 2017], as well as into generative models leading to improvements in performance and stability [Miyato et al., 2018, Vahdat and Kautz, 2020, Zhang et al., 2019]. We investigate the importance of smoothness regularisation and its effects in Chapter <ref> as well as incorporate it into a new problem domain, reinforcement learning, in Chapter <ref>.

§ INTERACTIONS BETWEEN OPTIMISATION AND SMOOTHNESS REGULARISATION

We have thus far discussed optimisation and smoothness regularisation as separate pillars of a deep learning system. Optimisation—concerned with changes in parameters—affects training performance, speed, and stability.
Smoothness—concerned with changes in inputs—can be a powerful inductive bias avoiding learning overly complex models.
What is often overlooked is the implicit bias optimisation has on the model class being learned.
While neural architectures and model regularisers define the class of functions we can represent with our models, optimisation methods determine which functions in the function class we can learn [Barrett and Dherin, 2021]. Likewise, different model architectures are easier to learn than others, with ResNets [He et al., 2016] providing a prime example of an architecture that aids optimisation [Li et al., 2018], followed by model normalisation techniques such as batch normalisation [Ioffe and Szegedy, 2015], which have been shown to primarily benefit optimisation [Bjorck et al., 2018].

While not previously studied, we show that smoothness regularisation and  optimisation heavily interact in deep learning.
We start by showing examples of this interaction in Chapter <ref>, where we observe that smoothness regularisation methods can interact with hyperparameter choices such as learning rates.
We then highlight how in some problem domains smoothness regularisation leads to increased performance by improving optimisation dynamics in Chapter <ref>. In Chapter <ref>, we show how the implicit regularisation induced by the discretisation drift of an optimisation method can lead to a smoothness inductive bias, and that the strengh of this inductive bias is dependent on optimsation hyperparameters such as learning rate and batch size.

Throughout this thesis, we will see examples of how optimisation can affect the learned model class and how the model smoothness can affect optimisation. While it is tempting to think of each component of a deep learning system individually, it can hinder progress and misguide our intuitions. We will use these observations to argue for a cohesive view between optimisation and regularisation in deep learning.

CHAPTER: OPTIMISATION IN DEEP LEARNING

In this chapter we describe the optimisation tools and methods we will require in the rest of this thesis.
When describing the commonly used algorithms in deep learning we often take a continuous-time approach.
Compared to a discrete-time approach, a continuous-time approach tends to be more amenable to proof construction [Elkabetz and Cohen, 2021], having many tools available, including stability analysis, and can be used to construct conserved quantities [Kunin et al., 2021]. Indeed, the ease of analysis of continuous-time approaches has lead to a general uptake into continuous-time methods in deep learning recently, with use in generative models, architectures, and optimisation [Chen et al., 2018, Chen et al., 2018, Kidger, 2022, Grathwohl et al., 2019, Song et al., 2021, Qin et al., 2020].

§ THE DIRECTION OF STEEPEST DESCENT

Once we have formulated our machine learning problem into an optimisation problem

\begin{align}
\min_{\vtheta \in \mathbb{R}^D} E(\vtheta),
\label{eq:opt_problem_2}
\end{align}

we have to choose an optimisation algorithm with which to solve the above problem.
Solving the problem in Eq (<ref>) in closed form is intractable for the problems we are interested in. Thus, the optimisation algorithms used are iterative, and aim to answer the question: given a current set of parameters $\vtheta$, what is a direction of descent of function $E$?
The simplest choice of the direction is given by the negative gradient $- \nabla_{\vtheta} E(\vtheta)$, known as the direction of steepest descent.
To see why, consider the continuous-time flow, often referred to as the gradient flow or negative gradient flow (NGF)

\begin{equation}
\dot{\vtheta} = - \nabla_{\vtheta} E.
\label{eq:ngf_first}
\end{equation}

If parameters $\vtheta$ follow the NGF in continuous-time, $E(\vtheta)$ will decrease until a stationary point $\nabla_{\vtheta} E = \mathbf{0}$ is reached, since

\begin{equation}
\frac{d E}{d t} = \frac{d \vtheta}{d t}^T \nabla_{\vtheta} E  =  - \left(\nabla_{\vtheta} E\right)^T  \nabla_{\vtheta} E = -  || \nabla_{\vtheta} E||^2 \le 0.
\label{eq:e_min_ngf}
\end{equation}

If $E$ is convex, following the NGF will find the global minimum.
If $E$ is not convex, following the NGF in continuous-time will reach a local minimum (a proof will be provided in Corollary <ref>). A point $\vtheta^* \in \mathbb{R}^D$ is a local minimum if there is a neighbourhood $\mathcal{V}$
around $\vtheta^*$ such that $\forall \vtheta' \in \mathcal{V} \hspace{1em} E(\vtheta') \le E(\vtheta^*)$; a more amenable characterisation is that local minima are stationary points $\nabla_{\vtheta} E(\vtheta^*) = \mathbf{0}$ with a positive semi-definite Hessian, i.e. $\nabla_{\vtheta}^2 E(\vtheta^*)$ has only non-negative eigenvalues. It is often convenient to also define strict local minima, i.e. local minima where eigenvalues are strictly positive.

We have already found two insights that would follow us throughout the rest of this thesis: the ease of analysis in continuous-time, and the importance of the negative gradient. Before exploring the optimisation algorithms using these insights, we briefly turn to the question of computing the gradient.

§.§ Computing the gradient

To compute the gradient $\nabla_{\vtheta} E$ we have to consider the specification of the loss function $E$.
In machine learning, $E$ often takes the form of an expectation or sum of expectations; if the expectation is under the data distribution $p^*(\vx)$ and we can write $E(\vtheta) = \mathbb{E}_{p^*(\vx)} l(\vtheta; \vx)$, we have

\begin{align}
\nabla_{\vtheta} E(\vtheta) = \nabla_{\vtheta} \mathbb{E}_{p^*(\vx)} l(\vtheta; \vx) =  \mathbb{E}_{p^*(\vx)} \nabla_{\vtheta} l(\vtheta; \vx).
\label{eq:swap}
\end{align}

We can now compute an unbiased estimate of the gradient via Monte Carlo estimation:

\begin{align}
\nabla_{\vtheta} E(\vtheta) = \mathbb{E}_{p^*(\vx)}\nabla_{\vtheta} l(\vtheta; \vx) \approx \frac{1}{N}\sum_{i=1}^{N}\nabla_{\vtheta} l( \vtheta; \vx_i), \hspace{1em} \vx_i \sim p^*(\vx).
\label{eq:mc}
\end{align}

In deep learning, the gradient $\nabla_{\vtheta} l( \vtheta; \vx_i)$ is often computed via backpropagation, i.e. using the compositional structure of neural networks to compute gradients in a recursive fashion via the chain rule.
If the entire dataset is used to approximate the gradient in Eq (<ref>), this corresponds to full-batch training; if an unbiased sample from the dataset is used, the training procedure is referred to as mini-batch training.

Additional challenges arise for objectives of the form $E(\vtheta) = \mathbb{E}_{p(\vx; \vtheta)} l(\vx)$, as is the case for some reinforcement learning and generative modelling objectives. In such situations, the change of expectation and gradient in Eq (<ref>) is no longer possible.
For certain family of distributions $p(\vx; \vtheta)$, however, gradient estimators are available.
In later chapters we will use the pathwise estimator, where the random variable $X$ with distribution $p(\vx; \vtheta)$ can be written as $X = g(Z; \vtheta)$, in which case we can use the change of variable formula for probability distributions:

\begin{align}
\mathbb{E}_{p(\vx; \vtheta)} l(\vx) = \mathbb{E}_{p(\vz)} l(g(Z; \vtheta)).
\end{align}

From here we can write

\begin{align}
\nabla_{\vtheta}\mathbb{E}_{p(\vx; \vtheta)} l(\vx) = \nabla_{\vtheta}\mathbb{E}_{p(\vz)} l(g(Z; \vtheta)) = \mathbb{E}_{p(\vz)} \nabla_{\vtheta} l(g(Z; \vtheta)),
\end{align}

which is again amenable to Monte Carlo estimation.
We refer to our general overview for more details on the pathwise estimator as well as other estimators [Mohamed et al., 2020].

§ ALGORITHMS: FROM GRADIENT DESCENT TO ADAM

We have seen how following the negative gradient in continuous-time leads to convergence to a local minimum.
Implementing infinitely small updates in continuous time on our discrete computers is not feasible, however, and thus we need a discrete-time algorithm for optimisation. Discretising continuous-time flows is a heavily studied topic with its own area of research in applied mathematics, numerical integration [Hairer et al., 2006].
Numerical integrators provide approaches to approximate the solution of the flow at time $h$ and initial conditions $\vtheta(0)$, which we denote as $\vtheta(h; \vtheta(0))$.

A simple discretisation method is Euler discretisation, which when applied to the NGF leads to

\begin{align}
\vtheta(h; \vtheta(0)) &= \vtheta(0) + \int_{0}^h \dot{\vtheta}(t) d t = \vtheta(0) - \int_{0}^h \nabla_{\vtheta} E(\vtheta(t)) d t \\
&\approx \vtheta(0) - \int_{0}^h \nabla_{\vtheta} E(\vtheta(0)) d t = \vtheta(0) - h \nabla_{\vtheta} E(\vtheta(0)). \label{eq:euler_int}
\end{align}

Setting $\vtheta(0) = \vtheta_{t-1}$ in the above equation leads to the familiar gradient descent update:

\begin{equation}
\vtheta_t = \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1}).
\end{equation}

By assuming that the vector field of the flow—the gradient $\nabla_{\vtheta} E$—does not change over the time interval $h$ we obtained a fast discrete-time algorithm.
However, due to the approximation in Eq (<ref>), following gradient descent is not the same as following the NGF: the numerical integration error leads to a difference between the gradient descent and NGF trajectories. Throughout this thesis we will refer to this difference in trajectories as discretisation drift (often also known as discretisation error).
Due to discretisation drift, there are no guarantees that gradient descent will reach a local minimum. Unlike the NGF,
gradient descent can diverge or converge to saddle points (stationary points that have at least one negative eigenvalue of the Hessian).

We visualise the effects of discretisation drift in Figure <ref>. Discretisation drift can speed up training (Figure <ref>), destabilise it (Figure <ref>) or lead to divergence  (Figure <ref>).
Importantly, when using gradient descent, the discrete counterpart of the NGF, for large $h$ the function $E$ can increase, leading to training instabilities.
Understanding when gradient descent leads to instabilities can be challenging, due to intricate dependencies between the learning rate and the shape of $E$; this gets compounded in the case of deep learning due to the very high dimensional nature of the parameter space. We will tackle this question later in this thesis.

Recent work has made great strides in understanding the behaviour of gradient descent in deep learning. Questions of study include whether gradient descent converges to local or global minima [Du et al., 2019, Jacot et al., 2018, Du et al., 2018, Bartlett et al., 2018], the prevalence and effect of saddle points [Dauphin et al., 2014, Du et al., 2017],
the interaction between neural architectures and optimisation [Li et al., 2018, Bjorck et al., 2018, Li et al., 2021], the connection between learning rates and Hessian eigenvalues and instabilities [Cohen et al., 2021],
the effect of the learning rate on generalisation [Barrett and Dherin, 2021, Li et al., 2019].

[The importance of the learning rate in gradient descent.]The importance of the learning rate in gradient descent: small learning rates lead to convergence to local minima but slow to converge fig:dd_intro_1, while larger learning rates can lead to instability fig:dd_intro_2 and divergence fig:dd_intro_3.

Rprop and RMSProp. We have derived gradient descent as a discretisation of the NGF; we used the NGF since it decreases the value of the objective function $E$. We note, however, that only the sign of the gradient for each parameter is relevant in order to achieve a descent direction locally. To see why, consider constants $c_i >0$ and the system:

\begin{equation}
\dot{\vtheta_i} = -  c_i \nabla_{\vtheta_i} E.
\end{equation}

This flow, like the NGF, also decreases $E$, since as $c_i >0$

\begin{equation}
\frac{d E}{d t} = \sum_i \frac{d E}{d\vtheta_i} \frac{d \vtheta_i}{d t} =  \sum_i  \frac{d E}{d\vtheta_i} (-  c_i \nabla_{\vtheta_i} E) = -  \sum_i  c_i (\nabla_{\vtheta_i} E)^2\le 0.
\label{eq:c_i_ode}
\end{equation}

This observation is the motivation for Rprop [Riedmiller and Braun, 1993], which instead of using  the gradient as the parameter update, uses only its sign: $sign(\nabla_{\vtheta_i} E) = \frac{\nabla_{\theta_i} E}{\sqrt{\left(\nabla_{\vtheta_i} E\right)^2}}$. We can think of Rprop as the discretisation of the following flow at each iteration:

\begin{equation}
\dot{\vtheta_i} = -  \frac{1}{|\nabla_{\vtheta_i} E(\vtheta_0)|} \nabla_{\vtheta_i} E .
\end{equation}

While this can solve issues with exploding gradients and an imbalance between the magnitude of the updates of different parameters, it can be crude: Rprop cannot distinguish between areas of space where gradient are positive, but very small or areas of space with very large gradients; the update in both of these cases would be the learning rate $h$. The solution proposed to this issue was RMSprop [Tieleman and Hinton, 2012], which instead of normalising the gradient of each parameter by $\left(\nabla_{\vtheta_i} E\right)^2$, it normalises the gradient by a moving average of $\left(\nabla_{\vtheta_i} E\right)^2$ obtained from previous iterations. By using element-wise operations, we can write the RMSprop update as

\begin{align}
\vv_t &= \beta \vv_{t-1} + (1-\beta) \nabla_{\vtheta} E(\vtheta_{t-1})^2; \hspace{5em}
\vtheta_t = \vtheta_{t-1} - h \frac{\nabla_{\vtheta} E(\vtheta_{t-1})}{\sqrt{\vv_t}},
\end{align}

where $\beta$ is a hyperparameter often set to values in $[0.9, 0.999]$.
Since the denominator is positive, the RMSprop update has the same sign as the gradient, but a different magnitude. The use of moving averages can dampen the effect of a large gradient, which would lead to instabilities in gradient descent, but still accounts for the the magnitude of the gradient, not only its sign, unlike Rprop.
RMSprop and its variants are often used in reinforcement learning [Mnih et al., 2015, Kearney et al., 2018, Hessel et al., 2018].

Momentum. Moving averages are an important staple of optimisation algorithms, beyond adjusting the magnitude of the local gradient, as we have seen with RMSprop. The idea behind momentum algorithms is to speed up or stabilise training by using previous iteration gradients.

To provide intuition regarding momentum, we note that along a trajectory towards a local minimum the sign of the gradient will be the same throughout that trajectory.
Let $\vtheta^*$  be a local minimum of $E(\vtheta)$ and $\vtheta_{t+1}$ and $\vtheta_t$ iterates along a trajectory towards $\vtheta^*$ and consider parameter index $i$ such that $sign(\vtheta_{t+1,i} - \vtheta^*_{i}) = sign(\vtheta_{t, i} - \vtheta^*_{i})$; this encodes the assumption that the trajectory goes towards a $\vtheta^*$ without escaping it in dimension $i$. But $sign(\vtheta_{t+1, i} - \vtheta^*_i) = sign(\nabla_{\vtheta} E (\vtheta_t)_i)$ as they are both the direction which minimises $E$. Thus, $sign(\nabla_{\vtheta} E (\vtheta_t)_i)=sign( \nabla_{\vtheta} E(\vtheta_{t+1})_i)$. This tells us that as long as a trajectory is not jumping over a local minimum in a direction, the sign of the gradient in that direction does not change.
Hence, one can speed up training by taking bigger steps in that direction by accounting for the previous updates; alternatively this is often framed as an approach to speed up training in areas of low curvature (where by definition the gradient does not change substantially) [Sutskever et al., 2013].
This intuition suggests the following algorithm [Polyak, 1964]:

\begin{align}
\vtheta_t = \vtheta_{t-1} - h \sum_{i=1}^t \beta^{i-1} \nabla_{\vtheta}  E(\vtheta_{t-i}),
\end{align}

where the update is given by a geometrical average of the gradients at previous iterations, and $\beta$ is the decay rate, a hyperparameter with $\beta < 1$.
We can rewrite the momentum algorithm in its commonly known form:

\begin{align}
\vv_t &= \beta \vv_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1}); \hspace{10em}
\vtheta_t = \vtheta_{t-1} + \vv_t .
\end{align}

Traditionally momentum has been seen as a discretisation of the second order flow

\begin{align}
{\ddot{\vtheta}} + c_v \dot{\vtheta} = - \nabla_{\vtheta} E (\vtheta),
\end{align}

from where the vector $\vv$ can be introduced:

\begin{align}
{\dot{\vv}} = - c_v \vv - \nabla_{\vtheta} E (\vtheta); \hspace{10em}
{\dot{\vtheta}} = \vv
\end{align}

which, when setting $c_v = \frac{1-\beta}{h}$ and Euler discretised with learning rates $h$ and $1$, leads to the discrete updates above. Momentum can also be seen as a discretisation of the first-order continuous-time flow, by combining implicit and explicit Euler discretisation [Kunin et al., 2021]:

\begin{equation}
\dot{\vtheta} = -\frac{1}{1 - \beta} \nabla_{\vtheta} E .
\label{eq:momentum_ode}
\end{equation}

Since $\beta < 1 \implies \frac{1}{1 - \beta} > 1$, this is in line with our previous intuition of speeding up training. As with gradient descent, however, this intuition does not account for the effect of large learning rates leading to discretisation errors; while following the flow in Eq (<ref>) always decreases $E$, following the trajectory of the momentum optimiser instead loses that guarantee.

[$\beta=0.9$ and small $h$.]

[$\beta=0.9$ and large $h$.]

[$\beta=0.5$ and large $h$.]

[The importance of the learning rate and decay rate in gradient descent with momentum.]The importance of the learning rate $h$ and decay rate $\beta$ in gradient descent with momentum. fig:mom1: momentum can speed up training but can also introduce instabilities when the gradient changes direction; the learning rate used is the same as in Figure <ref>, where gradient descent converged to the minimum without instabilities. fig:mom2 and fig:mom3: for learning rates where gradient descent diverges, as we observed in Figure <ref>, momentum can stabilise training by dampening the effect of local large gradient magnitudes and large learning rates; from this simple example one can already observe the intricacies of the interactions of decay rates and learning rates when momentum is used.

We show the effect of $\beta$ and $h$ on the optimisation trajectory in Figure <ref>.
While momentum can speed up training and reach a neighbourhood of a local minimum quicker than vanilla gradient descent, for large decay rates $\beta$ it can lead to large updates that move further away from the equilibrium. We show this effect
in Figure <ref>, where we observe that despite using a small learning rate, with a high decay rate $\beta$ momentum can oscillate around the minimum significantly before reaching it, while gradient descent with the same learning rate does not oscillate around the minimum (Figure <ref>).  This effect can be mitigated by using a smaller decay rate, as exemplified in Figure <ref>.
Dampening the strength of the current gradient update through momentum can also stabilise training, as we show in Figure <ref>. While for the same learning rate gradient descent diverges (Figure <ref>), momentum does not, since the large local gradients moving away from the equilibrium are dampened by the moving averages obtained from past gradients. For quadratic functions, this difference between gradient descent and momentum can be seen in the smallest learning rate which leads to divergence: while gradient descent diverges if $h > 2/\lambda_0$ where $\lambda_0$ is the largest eigenvalue of the Hessian $\nabla_{\vtheta}^2 E$, for momentum divergence occurs when $h > 2/\lambda_0 + 2\beta/\lambda_0\ge 2/\lambda_0$ since $\beta \ge 0$ [Cohen et al., 2021]. We will come back to the importance of the learning rate relative to the Hessian eigenvalues in later chapters.

Gradient descent with momentum is still very much used for deep learning optimisation and has been credited with increased generalisation performance compared to other optimisers [Gupta et al., 2021, Zhou et al., 2020], though recent work has shown that other well tuned methods can outperform momentum if well tuned [Kingma and Ba, 2015].

The Adam optimiser [Kingma and Ba, 2015], perhaps the most popular optimisation algorithm in deep learning, combines some of the previously discussed insights: it uses momentum and an RMSprop-like approach of dividing by a moving average of the square gradient (all operations below are element-wise):

\begin{align}
\vm_t &= \beta_1 \vm_{t-1} + (1 - \beta_1) \nabla_{\vtheta}E(\vtheta_{t-1})\\
\vv_t &= \beta_2 \vv_{t-1} + (1 - \beta_2) \nabla_{\vtheta}E(\vtheta_{t-1})^2\\
\vtheta_t &= \vtheta_{t-1} - h \frac{\frac{1}{1 - \beta_1^t} \vm_t}{\sqrt{\frac{1}{1 - \beta_2^t}\vv_t} + \epsilon},
\label{eq:last_adam}
\end{align}

where $\beta_1, \beta_2$ and $\epsilon$ are hyperparameters usually set to $\beta_2 = 0.999$ and $\beta_1 \in [0, 0.9]$ and $\epsilon \in [10^{-1}, 10^{-8}]$ with the default at $10^{-8}$.
The division in Eq (<ref>) is obtained using a bias correction argument for the moving averages. To see this, we write:

\begin{align}
\vm_t &= (1 - \beta_1) \sum_{i=1}^t \beta_1^{i-1} \nabla_{\vtheta}E(\vtheta_{t-i}).
\label{eq:adam_m_expanded}
\end{align}

The authors make the argument that if the gradients at iteration $i$ $\nabla_{\vtheta}E(\vtheta_{i})$ are drawn from the distribution $p_{g,i}$

\begin{align}
\mathbb{E}_{p_m} [\vm_t] &=   (1 - \beta_1) \sum_{i=1}^t \beta_1^{i-1} \mathbb{E}_{p_{g, i}} \left[\nabla_{\vtheta}E(\vtheta_{t-i})\right]\\
& = (1 - \beta_1^{t}) \mathbb{E}_{p_g} \left[\nabla_{\vtheta}E(\vtheta_{t-1})\right] + \vxi,
\end{align}

where $\vxi$ is a bias term accounting for the changes in the expected value of the gradients at different iterations; it will be $\mathbf{0}$ if the gradient mean does not change between iterations otherwise it is expected to be small due to the decaying contribution of older gradients.
Thus, to ensure the expected value of the moving average stays the same as the gradient and does not decrease steadily as training progresses, the division in Eq (<ref>) is used; a similar argument can be made for the second moment. This ensures that in expectation the update (if a small $\epsilon$ is used), is approximately of magnitude $h$ for each parameter. We note that while the assumption that the bias $\vxi$ is small does not hold for certain areas of training, such as when the gradient changes direction due to instability around an equilibrium, the above argument holds in areas of low curvature (i.e. areas of space where the gradient does not change substantially); this is something we will use later.

By investigating the commonly used Adam hyperparameters, $\beta_1 \in [0, 0.9]$ and $\beta_2 = 0.999$, we get intuition about its behaviour. Since $\beta_1$ is smaller than $\beta_2$, Adam is more sensitive to the sign of the gradient (controlled by $\beta_1$) than its magnitude (controlled by $\beta_2$). Figures <ref> and <ref> show that the exponential moving averages in the denominator of the Adam update can dampen large updates compared to momentum (the corresponding momentum figures are shown in Figure <ref> and <ref>), but oscillations around the minimum still occur.
When a low decay rate for the nominator is used, such as $\beta_1 = 0.5$, in conjunction with a large $h$, instabilities can occur if small gradients are followed by a large gradient, since $\vm_t$ will be large but $\vv_t$ will remain small due to $\beta_2 = 0.999$. We see this in Figure <ref>, which we can contrast with momentum for $\beta=\beta_1$ and the same learning rate $h$ in Figure <ref>: momentum exhibits stable training, and does not exit the area around the equilibrium. Relatedly, the use of moving averages in the denominator $\vv_t$ has been linked to the lack of convergence in Adam for simple online convex optimisation problems [Reddi et al., 2019].

[$\beta_1=0.9$, small $h$, $\epsilon=10^{-8}$.]

[$\beta_1=0.9$, large $h$, $\epsilon=10^{-8}$.]

[$\beta_1=0.5$, large $h$, $\epsilon=10^{-8}$.]

[$\beta_1=0.9$, small $h$, $\epsilon=10^{-2}$.]

[$\beta_1=0.9$, small $h$, $\epsilon=10^{-2}$.]

[$\beta_1=0.5$, large $h$, $\epsilon=10^{-2}$.]

[The effect of hyperparameters in Adam.]The importance of the Adam hyperparameters: learning rate $h$ and decay rate $\beta_1$ and $\epsilon$. We use the default $\beta_2 = 0.999$. Top row: the update normalisation in Adam can stabilise training but can also lead to instability compared to momentum; this row uses $\epsilon=10^{-8}$, the default value and the other hyperparameters are the same as in Figure <ref> which shows the corresponding momentum trajectories. Bottom row: while keeping other hyperparameters fixed, changing $\epsilon$ to $10^{-2}$ leads to vastly different behaviours and can stabilise training by avoiding dividing by small values.

The hyperparameter $\epsilon$ does not change the sign of the update, but it can change its magnitude. Specifically, a high $\epsilon$ decreases the magnitude of the parameter update. We show the effect of $\epsilon$ in Figure <ref>, where we show that using $\epsilon =10^{-2}$ (bottom row) can lead to more stable trajectories compared to $\epsilon =10^{-8}$ (top row).

Adam is a very popular the optimiser in deep learning, and has been attributed with progress in various fields including GANs [Radford et al., 2015], and training modern models such as transformers [Liu et al., 2020, Choi et al., 2019]. The $\epsilon$ hyperparameter in Adam has been noted to be important for certain machine learning applications [Choi et al., 2019], including reinforcement learning, where the $\epsilon$ used is large compared to supervised learning [Bellemare et al., 2017, Hessel et al., 2018].

§ CONTINUOUS-TIME METHODS

We now describe the continuous-time approaches we will use in the rest of this thesis.

§.§ A brief overview of stability analysis

Stability analysis is a tool for continuous-time systems that helps decide whether $\vtheta^* \in \mathbb{R}^D$ is a stable fixed point for a continuous-time flow

\begin{align}
\dot{\vtheta} = f(\vtheta),
\end{align}

where $f: \mathbb{R}^D \rightarrow \mathbb{R}^D$ is the vector field of the flow.
The most used definitions of stability are Lyapunov stability and asymptotic stability. Lyapunov stability requires

\begin{align}
\forall \epsilon > 0, \exists \delta > 0 \hspace{3em} \norm{\vtheta(0) - \vtheta^*} < \delta \implies \ \forall t \ge 0 \ \norm{\vtheta(t) - \vtheta^*} < \epsilon.
\end{align}

Asymptotic stability requires Lyapunov stability and additionally

\begin{align}
\exists \delta > 0 \hspace{3em} \norm{\vtheta(0) - \vtheta^*} < \delta \implies \lim_{t \rightarrow \infty} \norm{\vtheta(t) - \vtheta^*} = \mathbf{0}.
\end{align}

The rate of convergence of asymptotic stability is exponential if

\begin{align}
\exists \delta, C, \alpha > 0 \hspace{1em}   \norm{\vtheta(0) - \vtheta^*} < \delta \implies \forall t \ge 0 \ \norm{\vtheta(t) - \vtheta^*} \le C \norm{\vtheta(0) - \vtheta^*} e^{- \alpha t}.
\end{align}

Exponential asymptotic stability can be established using the following criterion:
(From [Hartman, 1960]) A fixed point $\vtheta^*$ with $ f(\vtheta^*) = \mathbf{0}$ where the  Jacobian of the vector field $\jacthetaf(\vtheta^*)$ has eigenvalues with strictly negative real part is a stable fixed point of the flow $\dot{\vtheta} = f(\vtheta)$, with asymptotic exponential convergence. In contrast, if a fixed point has a Jacobian with an eigenvalue with strictly positive real part, that is an unstable equilibrium and the flow will not converge to it.

Remark <ref> provides us with a useful test to establish convergence based on a fixed point's Jacobian eigenvalues.
The test is inconclusive for equilibria with non-strictly negative real part Jacobian eigenvalues, i.e. Jacobians where a subset of eigenvalues are have strictly negative real part while some have $0$ real part.
Throughout this thesis, when performing stability analysis of a continuous-time system we will be testing exponential asymptotic stability using the above criterion.
The NGF $ \dot{\vtheta} = - \nabla_{\vtheta} E$ is only locally attracted to local minima of $E$, and all strict local minima are exponentially attractive under the NGF.

This follows from the application of Remark <ref>. If a stationary point ${\nabla_{\vtheta} E(\vtheta^*) = \mathbf{0}}$ is attractive under the NGF, the eigenvalues of the Jacobian $\jacparam{\vtheta}{(-\nabla_{\vtheta}E)}\left(\vtheta^*\right)$ have strictly negative or 0 real part. Since $\jacparam{\vtheta}{(-\nabla_{\vtheta}E)}\left(\vtheta^*\right) = - \nabla_{\vtheta}^2 E$ is a symmetric matrix, this condition is equivalent to the eigenvalues of $\nabla_{\vtheta}^2 E$ being non-negative, the local minimum condition. Conversely, since strict local minima have strictly positive eigenvalues, the Jacobian has strictly negative eigenvalues at strict local minima, leading to exponential local stability under the NGF.

§.§ Backward error analysis

Backward error analysis (BEA) is a tool in numerical analysis developed to understand the discretisation error of numerical integrators. We now present an overview of how to use BEA in the context of the gradient descent update $\vtheta_t = \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1})$ with $\vtheta \in \mathbb{R}^D$; for a general overview see Hairer et al., 2006.
BEA provides a modified vector field

\begin{equation}
\tilde f_n(\vtheta) =  - \nabla_{\vtheta} E + h f_1(\vtheta) + \cdots + h^n f_n(\vtheta)
\end{equation}

by finding functions $f_1:\mathbb{R}^D \rightarrow \mathbb{R}^D$, ..., $f_n: \mathbb{R}^D \rightarrow \mathbb{R}^D$, such that the solution of the modified flow at order $n$, that is,

\begin{align}
\bm{{\dot{\vtheta}}} =  - \nabla_{\vtheta} E + h f_1(\vtheta) + \cdots + h^n f_n(\vtheta)
\label{eq:general_modified_vector_field}
\end{align}
follows the discrete dynamics of the gradient descent update with an error $\| \vtheta_t - \vtheta(h)\|$ of order $\mathcal O(h^{n+2})$, where $\vtheta(h)$ is the solution of the modified equation truncated at order $n$ at time $h$ and $\vtheta(0) = \vtheta_{t-1}$.
The full modified vector field with all orders ($n \rightarrow \infty$)

\begin{equation}
\tilde f(\vtheta) = - \nabla_{\vtheta} E + h f_1(\vtheta) + \cdots + h^n f_n(\vtheta) + \cdots
\label{eq:bea_series}
\end{equation}

is usually divergent and only forms an asymptotic expansion. What BEA provides is the Taylor expansion in $h$ of an unknown $h$-dependent vector field $f_h(\vtheta)$ developed at $h=0$:

\begin{equation}
\tilde f(\vtheta) = \textrm{Taylor}_{h=0} f_h(\vtheta).
\end{equation}

Thus, a strategy for finding $f_h$ is to find a series of the form in Eq (<ref>) via BEA and then find the function  $f_h$ such that its Taylor expansion in $h$ at 0 results in the found series.
Using this approach we can find the flow $ \dot{\vtheta} = f_h(\vtheta)$ which  describes the gradient descent step $\vtheta_t = \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1})$ exactly.

While flows obtained using BEA are constructed to approximate one gradient descent step, the same flows can be used over multiple gradient descent steps as shown in Section <ref> in the Appendix.

Barrett and Dherin, 2021 used this technique to find the $\mathcal{O}(h^2)$ correction term of gradient descent in the single-objective setting. They showed that for a model with parameters  $\vtheta$ and loss $E(\vtheta)$, optimised with gradient descent $\vtheta_t = \vtheta_{t-1} - h\nabla_{\vtheta} E(\vtheta)$, the first-order modified equation is

\begin{align}
\dot{\vtheta} = -\nabla_{\vtheta}E   -\frac{h}{2} \nabla_{\vtheta}^2 E \nabla_{\vtheta} E,
\label{eq:first_igr}
\end{align}

which can be written as

\begin{align}
\dot{\vtheta} =-\nabla_{\vtheta} \tilde E(\vtheta) = -\nabla_{\vtheta}\left(E(\vtheta) + \frac h4 \norm{\nabla_{\vtheta} E(\vtheta)}^2\right).
\end{align}

Thus, gradient descent can be seen as implicitly minimising the modified loss

\begin{align}
\tilde E(\vtheta) = E(\vtheta) + \frac h4 \norm{\nabla_{\vtheta} E(\vtheta)}^2.
\label{eq:sup_learning_igr}
\end{align}

This shows that when training models with gradient descent, there is an implicit regularisation effect, dependent on learning rate $h$, which biases learning towards paths with low gradient norms. The authors refer to this phenomenon as `implicit gradient regularisation'; we will thus refer to the flow in Eq (<ref>) as the IGR flow.

§.§ BEA proofs

The general structure of BEA proofs is as follows: start with a Taylor expansion in $h$ of the modified flow in Eq (<ref>); write each term in the Taylor expansion as a function of $\nabla_{\vtheta} E$ and the desired $f_i$ using the chain rule; group together terms of the same order in $h$ in the expansion; and identify $f_i$ such that all terms of $\mathcal{O}(h^p)$ are 0 for $p \ge 2$, as is the case in the gradient descent update. A formal overview of BEA proofs can be found in Section <ref> in the Appendix.

We now exemplify how to use BEA to find the IGR flow we showed in Eq (<ref>) [Barrett and Dherin, 2021].  Since we are only looking for the first correction term, we only need to find $f_1$.
We perform a Taylor expansion to find the value of $\vtheta(h)$ up to order $\mathcal O(h^{3})$ and then identify $f_1$ from that expression such that the error $\| \vtheta_t - \vtheta(h)\|$ is of order $\mathcal O(h^{3})$.
We have

\begin{align}
\vtheta(h) = \vtheta_{t-1} + h \vtheta^{(1)}(\vtheta_{t-1}) + \frac{h^2}{2} \vtheta^{(2)}(\vtheta_{t-1}) +  \mathcal{O}(h^3).
\end{align}

We know by the definition of the modified vector field in Eq (<ref>) that

\begin{align}
\vtheta^{(1)} = - \nabla_{\vtheta} E + h f_1({\vtheta}).
\end{align}

We can then use the chain rule to obtain

\begin{align}
\vtheta^{(2)} = \frac{d \left(- \nabla_{\vtheta} E + h f_1({\vtheta})\right)}{dt}  = - \frac{d \nabla_{\vtheta} E}{dt} + \mathcal{O}(h) = \nabla_{\vtheta}^2 E \nabla_{\vtheta} E + \mathcal{O}(h).
\end{align}

\begin{align}
\vtheta(h) = \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1}) + h^2 f_1(\vtheta_{t-1}) + \frac{h^2}{2}  \nabla_{\vtheta}^2 E (\vtheta_{t-1})\nabla_{\vtheta} E (\vtheta_{t-1})+  \mathcal{O}(h^3).
\end{align}

We can then write

\begin{align}
\vtheta_t - \vtheta(h) &= \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1}) - \vtheta(h) \\
&= h^2 f_1(\vtheta_{t-1}) + \frac{h^2}{2}  \nabla_{\vtheta}^2 E (\vtheta_{t-1})\nabla_{\vtheta} E (\vtheta_{t-1})+  \mathcal{O}(h^3).
\end{align}

For the error to be of order $\mathcal{O}(h^3)$ the terms of order $\mathcal{O}(h^2)$ have to be $\mathbf{0}$.
This entails

\begin{align}
f_1(\vtheta_{t-1}) =  -\frac{1}{2} \nabla_{\vtheta}^2 E(\vtheta_{t-1}) \nabla_{\vtheta} E(\vtheta_{t-1}),
\label{eq:igr_proof_value_iterate}
\end{align}

from which we conclude to $f_1 =  -\frac{1}{2} \nabla_{\vtheta}^2 E \nabla_{\vtheta} E $ leading to Eq (<ref>).

§ MULTIPLE OBJECTIVE OPTIMISATION

We have thus far discussed optimisers for single-objective problems, which encompasses many machine learning use cases, including supervised learning and certain formulations of reinforcement learning and generative modelling.
There are, however, use cases in machine learning beyond single-objective problems, including GANs and types of reinforcement learning, such as actor-critic frameworks.
These settings are often formulated via nested minimisation problems:

\begin{align}
\min_{\vtheta}   E_{\vtheta}(\arg \min_{\vphi} E_{\vphi}(\vphi, \vtheta), \vtheta).
\label{eq:nested_op}
\end{align}

Solving nested optimisation problems is computationally expensive, as it would require solving the inner optimisation problem for each parameter update of the outer optimisation problem. Thus, in deep learning practice nested optimisation problems are often translated into problems of the form

\begin{align}
&\min_{\vphi} E_{\vphi}(\vphi, \vtheta) \label{eq:prob_games_1} \\
&\min_{\vtheta} E_{\vtheta}(\vphi, \vtheta) \label{eq:prob_games_2}. \
\end{align}

Multi-objective optimisation problems of this kind are often cast as a game. Each player has a set of actions to take, in the form of an update for their parameters, $\vphi$ for the first player, and $\vtheta$ for the second. The loss functions are the negative of the player's respective payoff. As it is common in the literature, we will often adapt the game theory nomenclature in this setting;
this connection allows for a bridge between the machine learning and the game theory communities.

When translating problems of the form in Eq (<ref>) into optimisation procedures that preserve the nested structure of a game it is common to use one of the single-objective algorithms discussed Section <ref> to update the first player's parameters, followed by the use of the same algorithm to update the second player's parameters. This approach is often referred to as alternating updates; we summarise it in Algorithm <ref>. Optimisation approaches that further account for the nested game structure do exist [Metz et al., 2017], but they tend to be more computationally prohibitive and thus less used in practice.
In some settings, there is no nested structure or the nested structure is further discarded and both players are updated simultaneously, as we show in Algorithm <ref>.

Alternating updates

$k \ge 0$; the number of first player updates for each second player update
Choice of optimiser_routine($E$, current_value); e.g: $gd(E, \vtheta_t) = \vtheta_t - h \nabla_{\vtheta} E(\vtheta_t)$
Initialise $\vphi_0, \vtheta_0$
$t\gets 0$
$\vphi \gets \vphi_t$
$i$ in {1, ..., $k$}
$\vphi \gets$ optimiser_routine($E_{\vphi}(\cdot, \vtheta_t), \vphi)$
$\vphi_{t+1} \gets \vphi$
$\vtheta_{t +1} \gets$ optimiser_routine($E_{\vtheta}(\vphi_{t+1}, \cdot), \vtheta_t)$
$t\gets t+1$

Simultaneous updates

Initialise $\vphi_0, \vtheta_0$
Choice of optimiser_routine($E$, current_value); e.g: $gd(E, \vtheta_t) = \vtheta_t - h \nabla_{\vtheta} E(\vtheta_t)$
$\vphi_{t +1} \gets$ optimiser_routine($E_{\vphi}(\cdot, \vtheta_t), \vphi_t)$
$\vtheta_{t +1} \gets$ optimiser_routine($E_{\vtheta}(\vphi_t, \cdot), \vtheta_t)$
$t\gets t+1$

Simultaneous gradient descent updates derived for solving the problem in Eqs (<ref>) and (<ref>)—using gradient descent as the optimiser in Algorithm <ref>—can be seen as the Euler discretisation of the continuous-time flow

\begin{align}
&\dot{\vphi} = -\nabla_{\vphi} E_{\vphi} \label{eq:gen_games_ode1} \\
&\dot{\vtheta} = -\nabla_{\vtheta} E_{\vtheta} \label{eq:gen_games_ode2}.
\end{align}

Analysis of these game dynamics is more challenging than in the single objective case (shown in Eq (<ref>)), as we can no longer ascertain that following such flows minimises the corresponding loss functions since

\begin{align}
\frac{dE_{\vphi}}{d t} = \frac{d \vphi}{d t}^T \nabla_{\vphi} E_{\vphi} + \frac{d \vtheta}{d t}^T \nabla_{\vtheta} E_{\vphi}
= - \norm{\nabla_{\vphi} E_{\vphi}}^2 -  \nabla_{\vtheta} E_{\vtheta}^T \nabla_{\vtheta} E_{\vphi},
\end{align}

which need not be negative. Furthermore, it is unclear how to transfer the alternating structure that preserves player order from Algorithm <ref> into a continuous-time flow; we will address this in a later chapter. Nonetheless, the flows in Eqs (<ref>) and (<ref>) are often studied to understand the behaviour of gradient descent in nested optimisation and games [Nagarajan and Kolter, 2017, Balduzzi et al., 2018, Jin et al., 2020, Qin et al., 2020].

The challenges pertaining to analysing optimisation in games get further compounded when considering the notions of equilibria. The natural translation of local minimum from single-objective optimisation to games is a local Nash equilibrium. A local Nash equilibrium is a point in parameter space $(\vphi^*, \vtheta^*)$ such that no player has an incentive to deviate from $(\vphi^*, \vtheta^*)$ to another point in its local neighbourhood; that is $\vphi^*$ is a local minimum for $E_{\vphi}(\cdot, \vtheta^*)$ and $\vtheta^*$ is a local minimum for $E_{\vtheta}(\vphi^*, \cdot)$. Equivalently, $\nabla_{\vphi} E_{\vphi}(\vphi^*, \vtheta^*) = \mathbf{0}$ and $\nabla_{\vtheta} E_{\vtheta}(\vphi^*, \vtheta^*) = \mathbf{0}$ and $\nabla_{\vphi}^2
E_{\vphi}(\vphi^*, \vtheta^*)$ and $\nabla_{\vtheta}^2 E_{\vtheta}(\vphi^*, \vtheta^*)$ are positive semi-definite, i.e. do not have negative eigenvalues.

While in the single-objective case, local minima are the only exponentially attractive equilibria of the NGF (see Corollary <ref>),
this correspondence between desired equilibrium and attractive dynamics does not generally occur in games. By applying stability analysis (Remark <ref>), we know that the system in Eqs (<ref>)
and (<ref>) is attracted to equilibria for which the Jacobian

\begin{align}
\vJ =
\begin{bmatrix}
-\jactwoparam{\vphi}{\vphi}{E_{\vphi}} & - \jactwoparam{\vtheta}{\vphi}{E_{\vphi}} \\
- \jactwoparam{\vphi}{\vtheta}{E_{\vtheta}} & -\jactwoparam{\vtheta}{\vtheta}{E_{\vtheta}}
\end{bmatrix}
\end{align}

evaluated at $(\vphi^*, \vtheta^*)$ has only eigenvalues with strictly negative real part.
In contrast, a Nash equilibrium requires that the block diagonals of $\vJ$,  namely $-\nabla_{\vphi}^2 E_{\vphi}(\vphi^*, \vtheta^*)$ and $-\nabla_{\vtheta}^2 E_{\vtheta}(\vphi^*, \vtheta^*)$, have negative eigenvalues.
Since for an unconstrained matrix there is no relationship between its eigenvalues and the eigenvalues of its block diagonals, there is no inclusion relationship between Nash equilibria and locally attractive equilibria of the underlying continuous game dynamics that applies to all two-player games.

For zero-sum games, i.e. games where  $E_{\vphi} = - E_{\vtheta} = E$, a strict Nash equilibrium will be locally attractive when following the flow in Eqs (<ref>)
and (<ref>). To see why, we have to examine the eigenvalues of the flow's Jacobian

\begin{align}
\vJ =
\begin{bmatrix}
- \jactwoparam{\vphi}{\vphi}{E} & - \jactwoparam{\vtheta}{\vphi}{E} \\
\jactwoparam{\vphi}{\vtheta}{E} &  \jactwoparam{\vtheta}{\vtheta}{E}
\end{bmatrix}
\end{align}

evaluated at a strict Nash equilibrium of such game.
Due to the zero-sum formulation, the off-diagonal blocks of $\vJ$ satisfy $- \jactwoparam{\vtheta}{\vphi}{E} = - (\jactwoparam{\vphi}{\vtheta}{E})^T$, and we have

\begin{align}
\vJ + \vJ^T = \begin{bmatrix}
- \nabla_{\vphi}^2 E  & \mathbf{0} \\
\mathbf{0} &  \nabla_{\vtheta}^2 E
\end{bmatrix}.
\end{align}

Thus,  $\vJ + \vJ^T$ at the Nash equilibrium $(\vphi^*, \vtheta^*)$ is a symmetric block diagonal matrix with negative definite blocks, $-\nabla_{\vphi}^2 E(\vphi^*, \vtheta^*) $ and $\nabla_{\vtheta}^2 E(\vphi^*, \vtheta^*)$, and thus has negative eigenvalues. From here it follows that $\vJ$ has eigenvalues with strictly negative real part, and thus satisfies the condition of Remark <ref>.
We note that the reverse does not hold, as being locally attractive in this case is a weaker condition than that of a Nash equilibrium.
This observation together with the realisation that in deep learning certain games might not have Nash equilibria [Farnia and Ozdaglar, 2020], has lead the search for other equilibrium measures including local stability and Stackelberg equilibria [Fiez et al., 2020, Berard et al., 2020, Wang et al., 2019].

§ CONCLUSION

In this section we have introduced the main optimisation algorithms used in deep learning, together with their main challenges and the analytical tools we will employ  in the rest of this thesis.

CHAPTER: A NEW CONTINUOUS-TIME MODEL OF GRADIENT DESCENT

CHAPTER: CONTINUOUS TIME MODELS OF OPTIMISATION IN TWO-PLAYER GAMES

§ A COMMENT ON DIFFERENT LEARNING RATES

When considering different learning rates for the two players, we have thus far kept the consistency between physical time and learning rates. That is, we assumed the gradient descent updates are obtained by Euler discretisation with learning rates $\lrp h$ and $\lrt h$ of the flow

\begin{align}
\dot{\vphi} &=  f( \vphi, \vtheta)    \label{eq:old_two_lr_ode1} \\
\dot{\vtheta}  &= g( \vphi, \vtheta) \label{eq:old_two_lr_ode2}.
\end{align}

We then used BEA to fine modified flow such that after time $\lrp h$ and $\lrt h$ the error between the modified flow and the discrete updates is $\mathcal{O}(h^3)$.

Alternatively, we can consider the gradient descent updates as the result of Euler discretisation with learning rate $h$ of the flow

\begin{align}
\dot{\vphi} &=  \lrp f( \vphi, \vtheta)  \label{eq:new_two_lr_ode1}  \\
\dot{\vtheta}  &= \lrt g( \vphi, \vtheta) .
\label{eq:new_two_lr_ode2}
\end{align}

While in this approach the learning rate loses the connection with physical time for each player, it keeps the consistency of time between two players: after one iteration the same amount of time, $h$, has passed for both players. With this in mind, we will call this `the same-physical time between players' approach.

Considering gradient descent as a discretisation of Eqs (<ref>) and (<ref>) instead of Eqs (<ref>) and (<ref>) has implications for the analysis of the behaviour of gradient descent using continuous-time tools such as stability analysis, as the two systems can lead to different results based on the two learning rates. Using Eqs (<ref>) and (<ref>) can also lead to a different perspective from the loss minimisation view. For the game with $f = - \nabla_{\vphi} E$ and $g = \nabla_{\vtheta} E$, the second system leads to the interpretation that the first player is minimising function $\lrp E(\vphi, \vtheta)$ while the second player is maximising function $\lrt E(\vphi, \vtheta)$. This is different than the zero-sum formulation of the first player minimising $E(\vphi, \vtheta)$ while the second player is maximising it, resulting from our previous interpretation.
The same-physical time between players approach also presents the following challenge:
Given learning rates $\lrp h$ and $\lrt h$ for the two-players, $h$ cannot be identified exactly and has to be chosen by the user. The choice of $h$ leads to different flows describing the dynamics of the discrete updates. This is in contrast with the approach taken previously, which only depends on $\lrp$ and $\lrt$ through the learning rates $\lrp h$ and $\lrt h$.

Following the same-physical time between players approach also has implications for the use of BEA to find flows with error $\mathcal{O}(h^3)$ after one Euler update. Luckily, to find the modified flows using BEA under the same-physical time approach, we simply have to apply our results from Section <ref> for the system in Eqs (<ref>) and (<ref>) by using their corresponding vector fields and learning rate $h$. We do so, and obtain:

The discrete simultaneous Euler updates
in Eqs (<ref>) and (<ref>) follow the continuous system

\begin{align}
\dot{\vphi} &= \lrp f - \frac{h}{2} \left(\lrp^2 \jacphif f + \lrp \lrt \jacthetaf g \right) \\
\dot{\vtheta} &= \lrt g - \frac{h}{2} \left(\lrp \lrt \jacphig f + \lrt^2 \jacthetag g \right)
\end{align}

with an error of size $\mathcal O(h^3)$  after one update step.

The discrete alternating Euler updates in (<ref>) and (<ref>) follow the continuous system

\begin{align}
\dot{\vphi} &= \lrp f - \frac {h}{2} \left(\frac{\lrp^2}{m}\jacphif f + \lrp \lrt \jacthetaf g \right) \\
\dot{\vtheta} &= \lrt g - \frac{h}{2} \left(-  \lrt \lrp \jacphig f + \frac{\lrt^2}{k} \jacthetag g \right)
\end{align}

with an error of size $\mathcal O(h^3)$ after one update step.

We observe that both modified flows each player's vector field depends on the learning rate of the other player, unlike in our previous results where for simultaneous updates, the vector field for each player depended only on its learning rate.
When interpreting these results to obtain implicit regularisation effects as we have done previously in zero-sum and common-payoff games, the resulting coefficients of the implicit regularisers are different than those obtained in Section <ref>. For zero-sum games, we obtain the following corollaries:
In a zero-sum two-player differentiable game, simultaneous gradient descent updates—as described in Eqs (<ref>) and (<ref>)—follows a gradient flow given by the modified losses

\begin{align}
\tilde E_{\vphi}&= - \lrp E + \frac{\lrp^2 h}{4} \norm{\nabla_{\vphi} E}^2  -  \frac{\lrp \lrt h}{4} \norm{\nabla_{\vtheta} E}^2
\label{eq:min_min_simultaneous_zero_sum1-diff-lr} \\
\tilde E_{\vtheta}&= \lrt E -  \frac{\lrp \lrt h}{4} \norm{\nabla_{\vphi} E}^2 + \frac{\lrt^2 h}{4} \norm{\nabla_{\vtheta} E}^2
\label{eq:min_min_simultaneous_zero_sum2-diff-lr}
\end{align}

with an error of size $\mathcal O(h^3)$ after one update step.

In a zero-sum two-player differentiable game, alternating gradient descent—as described in Eqs (<ref>) and (<ref>)—follows a gradient flow given by the modified losses

\begin{align}
\tilde E_{\vphi}&= -\lrp E + \frac{\lrp^2 h}{4m} \norm{\nabla_{\vphi} E}^2  - \frac{\lrp \lrt h}{4} \norm{\nabla_{\theta} E}^2
\label{eq:min_min_alt_zero_sum1_diff_lr}
\\
\tilde E_{\vtheta}&= \lrt E + \frac{\lrp \lrt h}{4} \norm{\nabla_{\vphi} E}^2  + \frac{\lrt^2 h}{4k}\norm{\nabla_{\theta} E}^2
\label{eq:min_min_alt_zero_sum2_diff_lr}
\end{align}

with an error of size $\mathcal O(h^3)$ after one update step.

From Corollary <ref> one concludes for alternating updates that the second player's interaction always minimises the gradient norm of the first player, which again helps explain the empirically observed stability of alternating updates over simultaneous updates; this is in contrast with our previous interpretation which suggested that this occurs only for certain learning rate ratios (Corollary <ref>).

Simultaneous updates       First player ($\vphi$)     Second player ($\vtheta$)

Corollary <ref>      $\frac{\lrp^2 h^2}{4} \nabla_{\vphi} \norm{\nabla_{\vtheta} E}^2$      $\frac{\lrt^2 h^2}{4} \nabla_{\vtheta} \norm{\nabla_{\vphi} E}^2$

Corollary <ref>      $\frac{\lrp \lrt h^2}{4} \nabla_{\vphi} \norm{\nabla_{\vtheta} E}^2$      $\frac{\lrp \lrt h^2}{4} \nabla_{\vtheta} \norm{\nabla_{\vphi} E}^2$

Alternating updates       First player ($\vphi$)     Second player ($\vtheta$)

Corollary <ref>       $\frac{\lrp^2 h^2}{4} \nabla_{\vphi} \norm{\nabla_{\vtheta} E}^2$      $\left(\frac{\lrt^2 h}{4} - \frac{2 \lrp \lrt h^2}{4}\right) \nabla_{\vtheta} \norm{\nabla_{\vphi} E}^2$

Corollary <ref>       $\frac{\lrp \lrt h^2}{4} \nabla_{\vphi} \norm{\nabla_{\vtheta} E}^2$      $-\frac{\lrp \lrt h^2}{4} \nabla_{\vtheta} \norm{\nabla_{\vphi} E}^2$

[Contrasting the implicit regularisation effects of the continuous-time flows proposed in zero-sum games.]The strength of the interaction terms in DD under the different modified continuous-time flows we find for simultaneous and alternating gradient descent in zero-sum games. We obtained the different flows via different approaches of handling different learning rates for the two players, namely $\lrp h \ne \lrt h$. We do not write the self terms here as they are the same for both approaches. We write the implicit regularisation coefficients as observed in the continuous-time flow displacement for one gradient descent update, rather than the modified loss, in order to account for the different learning rates under the two interpretations (which are $\lrp h$ and $h$ for the first player, respectively). For simultaneous updates, the insight that the interaction terms for each player leads to a maximisation of the other player's gradient norm remains true regardless of interpretation, though the strength of the regularisation changes in the two interpretations. For alternating updates, the different interpretations suggest not only different regularisation strengths but also signs, with the result from Corollary <ref> showing that for alternating updates, the interaction term  of the second player always minimises the gradient norm of the first player, as opposed to the result in Corollary <ref>, for which this effect depends on the learning rate ratio $\lrp/\lrt$. Both results are consistent in explaining why alternating updates are more stable than simultaneous updates, since in both cases for alternating updates the second player's interaction term has a smaller coefficient when maximising the first player's gradient norm compared to that of simultaneous updates, or is minimising the first player's gradient norm.

If we would like to implement explicit regularisation methods to cancel the implicit regularisers obtained above using the same-physical time between players approach, we have to discretise the modified losses with learning rate $h$; this is opposed to using learning rates $\lrp h$ and $\lrt h$ to discretise the modified flows obtained in Section <ref>. For the self terms, no difference is obtained compared to the previous interpretation; under both approaches of tackling different learning rates, the strength of the self terms in the parameters update is proportional to $\lrp^2 h^2$ and $\lrt^2 h^2$, respectively. For interaction terms, however, the coefficients required in the parameter update are changed. We highlighted these changes in Table <ref>; these coefficients also show that while in order to find the modified flow in the same-physical time between player approach we require to choose $\lrp$, $\lrt$ and $h$, for explicit regularisation we only need the learning rates of the two players.

We thus can investigate the effect of explicit regularisation to cancel the interaction terms under the same-physical time approach in zero-sum games. As before, we do so by using the generator saturating loss in GAN training. We show results in Figures <ref> and <ref> and compare with our previous results from Section <ref>. While for the simultaneous updates using the same physical time approach leads to more hyperparameters with increased performance, no significant difference is obtained for alternating updates.

[An alternative interpretation to different learning rates leads to a new set of flows for simultaneous gradient descent updates.]Simultaneous updates in zero-sum games: cancelling the interaction terms under the different interpretations of tackling different learning rates (Corollaries <ref> and <ref>) does not result in a significant empirical difference. We denote the approach of Corollary <ref> as `same physical time', as per its motivation of using the same physical time between the two players.

[An alternative interpretation to different learning rates leads to a new set of flows for alternating gradient descent updates.]Alternating updates in zero-sum games: cancelling the interaction terms under the different interpretations of tackling different learning rates (Corollaries <ref> and <ref>) does not result in a significant empirical difference. We denote the approach of Corollary <ref> as `same physical time', as per its motivation of using the same physical time between the two players.

We conclude by noting that when studying two-player games trained with different learning rates, one can choose between two approaches of constructing corresponding continuous dynamics and modified losses. Both are valid, in that both lead to modified flows with the same error in learning rate, and can be used to understand implicit regularisation effects and the discrepancies between simultaneous and alternating Euler updates. Additional examination is needed to understand how these approaches explain continuous-time behaviour over a larger number of steps; we leave that as future work.

§ EXTENDING THE PF FOR TWO-PLAYER GAMES

In Chapter <ref>, we introduced the principal flow (PF) as a model of single-objective gradient descent. We now briefly show how the same approach can be generalised to simultaneous Euler updates in two-player games, which will allow us to capture discretisation drift effects beyond those of order $\mathcal{O}(h^3)$ we explored in the rest of this chapter. For simplicity, we use assume equal learning rates for the two-players, though one can use the `same-physical time approach' from the previous section to adapt the results to different learning rates.

The approach we took in the single-objective setting had two steps: first, we developed the principal series with BEA (Theorem <ref>) and second, we used the principal series given by BEA to find the PF based on the eigen-decomposition of the Hessian (Corollary <ref> and Theorem <ref>). The first step can be readily translated to games as it applies to any vector field, using  Theorem <ref> in the Appendix, which Theorem <ref> is a corollary of. For completeness, we reproduce Theorem <ref> here:

The modified flow given by BEA with an error of order $\mathcal{O}(h^{p+1})$ to the Euler update  ${\vpsi_t = \vpsi_{t-1} + h \gamesvf(\vpsi_{t-1})}$
has the form:

\begin{align}
\dot{\vpsi} = \sum_{n=0}^{p} \frac{(-1)^{n}}{n+1} {(\jacparam{\vpsi}{\gamesvf})}^n \gamesvf + \csecondorderfamv
\end{align}

where $\csecondorderfamv$ denotes a class of functions defined as a sum of individual terms each containing higher than first order derivatives applied to $\gamesvf$.

To use the above result in games, consider the two-player game with the original the dynamics in Eqs (<ref>)-(<ref>):

\begin{align}
\dot{\left[ \begin{array}{c}
\vphi \\
\vtheta
\end{array}\right]}
&=  \left[ \begin{array}{c}
f \\
\end{array}\right].
\end{align}

If we define

\begin{align}
\vpsi
=  \left[ \begin{array}{c}
\vphi \\
\vtheta
\end{array}\right]
\hspace{3em}
=  \left[ \begin{array}{c}
f(\vphi, \vtheta) \\
g(\vphi, \vtheta)
\end{array}\right]
\end{align}

and the negative of the Jacobian as

\begin{align}
\vH(\vphi, \vtheta) = -
\begin{bmatrix}
\jacparam{\vphi}{f} & \jacparam{\vtheta}{f} \\
\jacparam{\vphi}{g} & \jacparam{\vtheta}{g}
\end{bmatrix},
\end{align}

we can replace this choice vector field and $\vH$ in Theorem <ref>, and obtain:

The full series of BEA constructed from simultaneous Euler updates (Eqs <ref> and <ref>) with learning rate $h$ exactly is of the form

\begin{align}
\dot{\left[ \begin{array}{c}
\vphi \\
\vtheta
\end{array}\right]}
&=  \sum_{p=0}^{\infty} \frac{1}{p+1} h^p {\vH(\vphi, \vtheta)}^p \left[ \begin{array}{c}
f \\
\end{array}\right] + \csecondorderfam.
\end{align}

We can no longer use the eigen-decomposition of $\vH$, as unlike in the single-objective case ($f = - \nabla_{\vphi} E$ and $g = - \nabla_{\vtheta} E$)  we examined in Chapter <ref>, $\vH(\vphi, \vtheta)$ need not be symmetric and thus we can no longer conclude Corollary <ref>. We instead consider the Jordan normal form of $\vH$

\begin{align}
\vH(\vphi, \vtheta) = \vP^{-1} \vJ \vP,
\end{align}

where $\vJ$ is a block diagonal matrix. We then have

\begin{align}
{\vH(\vphi, \vtheta)}^p = \vP^{-1} \vJ^p \vP,
\end{align}

leading to

\begin{align}
\dot{\left[ \begin{array}{c}
\vphi \\
\vtheta
\end{array}\right]}
&=  \sum_{p=0}^{\infty} \frac{-1}{p+1} \vP^{-1} (h\vJ)^p \vP \left[ \begin{array}{c}
-f \\
\end{array}\right] + \csecondorderfam,
\end{align}

where we use the negative of the vector field in order to obtain consistency with the PF results in Chapter <ref>, and recover single-objective results if $f = - \nabla_{\vphi} E$ and $g = - \nabla_{\vtheta} E$.
We are interested in finding a form for the series

\begin{align}
\sum_{p=0}^{\infty} \frac{-1}{p+1} (h \vJ)^p ,
\end{align}

which if converges can be written as another block diagonal matrix,

\begin{align}
\frac{1}{h} \log(\vI - h \vJ) \vJ^{-1}.
\end{align}

This result leads to the game equivalent of the PF:
When expanding the PF to simultaneous Euler updates (Eqs <ref> and <ref>) with learning rate $h$ in two-player games, we obtain:

\begin{align}
\dot{\left[ \begin{array}{c}
\vphi \\
\vtheta
\end{array}\right]}
&= \frac{1}{h}  \vP^{-1}  \log(\vI - h\vJ) \vJ^{-1} \vP \left[ \begin{array}{c}
-f  \\
\end{array}\right] + \csecondorderfam,
\label{eq:gen_pf_games}
\end{align}

where $\vH(\vphi, \vtheta) = \vP^{-1} \vJ \vP$ is the Jordan normal form of $\vH$.

While the Jordan normal form of $\vH$ is not unique due to permutations of the order of the blocks in $\vJ$, $\vP^{-1}  \log(\vI - h\vJ) \vJ^{-1} \vP$ will be the same regardless of the block order.

§.§ Stability analysis

We now perform stability analysis on the modified flow obtained in Eq (<ref>), and hope it sheds light on the behaviour of Euler updates—and thus gradient descent—in games. We start by estimating the eigenvalues of the flow's Jacobian:

\begin{align}
\frac{1}{h} \vP^{-1}  \log(\vI - h \vJ) \vP.
\end{align}

As $\vP$ is invertible, the Jacobian's eigenvalues are the same as the eigenvalues of

\begin{align}
\frac{1}{h}\log(\vI - h \vJ).
\end{align}

Since $\log(\vI - h \vJ)$ is block diagonal, its eigenvalues are the eigenvalues of its blocks $\frac{1}{h}\log(\vI_{n_k} - h \vJ_k)$, where $ \vI_{n_k}$ is the identity matrix of dimension $n_k$  and  $\vJ_k \in \mathbb{C}^{n_k, n_k}$.
Each block $\vJ_k$ corresponds to a unique eigenvalue $\lambda_k$ of $\vH$ with duplicity $n_k$. Given

\begin{align}
\vZ_{n_k} = \begin{bmatrix}
0 & 1&  \hdots & 0 & 0 \\
\vdots & \vdots& \ddots & \ddots &  \vdots\\
0 & 0 & 0 & \hdots & 1 \\
0 & 0 & 0 & 0 & 0
\end{bmatrix},
\end{align}

a matrix with all $0$s apart from the one off diagonal, $\vJ_{k}$ can be written as

\begin{align}
\vJ_{k} = \lambda_k \vI_{n_k} + \vZ_{n_k}.
\end{align}

From here

\begin{align}
\log (\vI - h \vJ_{k}) = \log(\vI - h\lambda_k \vI_{n_k} - h \vZ_{n_k}),
\end{align}

and since $\vI - h \vJ_{k}$ is upper triangular, its eigenvalues are its diagonal entries, which are $1 - h \lambda_k$ and the eigenvalues of $\log (\vI - h \vJ_{k})$ are $\log(1 - h \lambda_k)$. Thus, to perform stability analysis we have to assess whether $\Re{[\frac{1}{h}\log(1 - h \lambda_k)]} < 0$ for all $\lambda_k$. Since $h>0$, this is equivalent to $\Re{[\log(1 - h \lambda_k)]} < 0$ for all $\lambda_k$. Unlike the single-objective case, $\lambda_k$ might be complex.
If we write $\lambda_k = x_k + i y_k$ we have:

\begin{align}
\Re[\log(1 - h (x_k + i y_k))] = \log \sqrt{(1 - h x_k)^2 + h^2 y_k^2},
\end{align}

which is negative if $(1 - h x_k)^2 + (h y_k)^2 < 1$.

We now contrast this condition to the exponential stability convergence condition of the original continuous-time game dynamics. Under the original continuous-time dynamics, for a fixed point to be attractive, the real part of the $\vH$'s eigenvalues would need to be strictly positive since $\vH$ is the negative of the flow's Jacobian, i.e. $x_k > 0$. When accounting for DD, however, the importance of the imaginary part $y_k$ becomes apparent.
Moreover, if we rewrite the condition  $(1 - h x_k)^2 + (h y_k)^2 < 1$ as a condition on the maximum learning rate to be convergent, we obtain

\begin{align}
(1 - h x_k)^2 + (h y_k)^2 < 1 \implies h < \frac{1}{x_k} \frac{2}{1 + \left(\frac{y_k}{x_k}\right)^2},
\end{align}

which is the condition Mescheder et al., 2017 (their Lemma 4) obtained from a discrete-time perspective. Thus in the game setting too, we recover existing fundamental results using a continuous-time perspective. We note that using this intuition Mescheder et al., 2017 derive the algorithm Consensus Optimisation (CO), which we empirically assessed in Section <ref>.

If $f = - \nabla_{\vphi} E$ and $g = - \nabla_{\vtheta} E$ we have a single objective, with $y_k = 0$ $\forall k$. In this case we recover the convergence conditions we are familiar with from the PF, which we derived in Section <ref>.

In the zero-sum case, where $f = - \nabla_{\vphi} E$ and $g = \nabla_{\vtheta} E$, we have seen in Section <ref> that strict local Nash equilibria are local stable fixed points under the original game dynamics. That is, at a strict local Nash equilibrium we have that $x_k > 0, \forall k$. Under the continuous-time dynamics we derived in this section, which account for DD, this is no longer true, since $x_k$ being positive does not lead to $(1 - h x_k)^2 + (h y_k)^2 < 1$. We thus obtain an analogous result in zero-sum games to that of the single-objective case in Chapter <ref>, where we saw that unlike the NGF, whether or not the PF is attracted to strict local minima depends on the learning rate $h$.

§.§ An empirical example and different learning rates

[Same learning rate.]

[Different learning rates.]

[A visualisation of a simple example of the extension of the PF to two-games.]A visualisation of a simple example of the extension of the PF to two-games. Here we observe that the PF better tracks the gradient descent update compared to the flow obtained from using the $\mathcal{O}(h^3)$ error, and the NGF.

The game PF describes the Euler discretisation of the system in Eq (<ref>) exactly (see Section <ref> for a full proof). We visualise this example in Figure <ref>, where we observe that the PF tracks the behaviour of the discrete updates better than the modified third-order correction previously derived.
We note that the Jordan normal form (here implemented using `sympy' [Meurer et al., 2017]) is not always numerically stable and the observed (small) difference between the discrete updates and the PF is due to numerical issues.

If different learning rates are used for the two-players, we can use the approach outlined in Section <ref> to construct vector fields that account for the each player's learning rate, and then apply the same analysis as above to the corresponding Jacobian to construct a flow and stability analysis conditions. We use this approach to visualise the same example from Eq (<ref>) in Figure <ref>, and observe that here too the PF tracks the Euler updates closely.

CHAPTER: FINDING NEW IMPLICIT REGULARISERS BY REVISITING BACKWARD ERROR ANALYSIS

In previous chapters, we used Backward Error Analysis (BEA) to construct flows that approximate a discrete optimiser update and shed light on optimisation dynamics, both for single-objectives and two-player games trained with gradient descent.
The current usage of BEA is not without limitations, however, since not all the vector fields of continuous-time flows obtained using BEA can be written as a gradient, hindering the construction of modified losses revealing implicit regularisers.  Implicit regularisers uncover quantities minimised by following discrete updates and thus can be used to improve performance and stability across problem domains, from supervised learning to two-player games such as Generative Adversarial Networks [Barrett and Dherin, 2021, Smith et al., 2021, Geiping et al., 2022, Schäfer et al., 2020]; we have seen an example in common-payoff and zero-sum games in the previous chapter. In this chapter, we provide a novel approach to use BEA, and show how our approach can be used to construct continuous-time flows with vector fields that can be written as gradients. We then use this to find previously unknown implicit regularisation effects, such as those induced by multiple stochastic gradient descent steps while accounting for the exact data batches used in the updates, and in generally differentiable two-player games.

§ REVISITING BEA

In Section <ref>, we described BEA and showed how Barrett and Dherin, 2021 use BEA to find the IGR flow

\begin{align}
\dot{\vtheta} = -\nabla_{\vtheta}E   -\frac{h}{2} \nabla_{\vtheta}^2 E \nabla_{\vtheta} E  = -\nabla_{\vtheta}\left(E(\vtheta) + \frac h4 \norm{\nabla_{\vtheta} E(\vtheta)}^2\right),
\end{align}

which has an error of $\mathcal{O}(h^3)$ after one gradient descent update ${\vtheta_t = \vtheta_{t-1} - h\nabla_{\vtheta} E(\vtheta)}$.
Gradient descent can thus be seen as implicitly minimising the modified loss $E(\vtheta) + \frac h4 \norm{\nabla_{\vtheta} E(\vtheta)}^2$. This showcases an implicit regularisation effect induced by the discretisation drift of gradient descent, dependent on learning rate $h$, which biases learning towards paths with low gradient norms.

We would like to go beyond full-batch gradient descent and model the implicit regularisation induced by the discretisation of one stochastic gradient descent (SGD). Given the SGD update $\vtheta_t = \vtheta_{t -1} - h \nabla_{\vtheta} E(\vtheta_{t -1}; \vX_{t})$—where we denoted $E(\vtheta_{t-1}; \vX) = \frac{1}{B}\sum_{i=1}^B E(\vtheta_{t-1}; \vx_i)$ as the average loss over batch $\vX$ and $B$ is the batch size—
we can use the IGR flow induced at this time step:

\begin{align}
\dot{\vtheta} = - \nabla_{\vtheta} E(\vtheta; \vX_{t}) - \frac h 4 \nabla_{\vtheta} \norm{\nabla_{\vtheta} E(\vtheta; \vX_{t})}^2.
\label{eq:igr_stochastic}
\end{align}

Since the vector field in Eq (<ref>) is the negative gradient of the modified loss ${E(\vtheta; \vX_{t}) +  \frac h 4 \norm{\nabla_{\vtheta} E(\vtheta; \vX_{t})}^2}$, this reveals a local implicit regularisation which minimises the gradient norm $ \norm{\nabla_{\vtheta} E(\vtheta; \vX_{t})}^2$. It is not immediately clear, however, how to combine the IGR flows obtained for each SGD update in order to model the combined effects of multiple SGD updates, each using a different batch. What are, if any, the implicit regularisation effects induced by two SGD steps

\begin{align}
&\vtheta_{t} = \vtheta_{t-1} - h \nabla_{\vtheta} E(\vtheta_{t-1}; \vX_{t}) \label{eq:sgd1} \\
&\vtheta_{t+1} = \vtheta_t - h \nabla_{\vtheta} E(\vtheta_t; \vX_{t+1})?
\label{eq:sgd2}
\end{align}

Smith et al., 2021 find a modified flow in expectation over the shuffling of batches in an epoch, and use it to find implicit regularisation effects specific to SGD and study the effect of batch sizes in SGD. Since their approach works in expectation over an epoch, however, it does not account for the implicit regularisation effects of a smaller number of SGD steps, or account for the exact data batches used in the updates; we return to their results in the next section.
We take a different approach, and introduce a novel way to find implicit regularisers in SGD by revisiting the BEA proof structure and the  assumptions made thus far when using BEA.
Our approach can be summarised as follows:

Given the discrete update $\vtheta_t = \vtheta_{t-1} - h \nabla_{\vtheta}E (\vtheta_{t-1})$,
BEA constructs $\dot{\vtheta} = \nabla_{\vtheta} E + h f_1(\vtheta) + \cdots + h^p f_p(\vtheta)$, such that $\norm{\vtheta(h; \vtheta_{t-1}) - \vtheta_t} \in \mathcal{O}(h^{p+2})$ for a choice of $p \in \mathbb{N}, p \ge 1$. This translates into a constraint on the value of $\vtheta(h;\vtheta_{t-1})$.
Thus, BEA asserts only what the value of the correction terms $f_i$ in the vector field of the modified flow is at $\vtheta_{t-1}$.
Given the constraints on $f_i(\vtheta_{t-1})$, we can choose  $f_i: \mathbb{R}^D \rightarrow \mathbb{R}^D$ in the vector field of the modified flow $\dot{\vtheta}$ to depend on the initial condition $\vtheta_{t-1}$.

For example, in the proof by construction of the IGR flow (Section <ref>), we obtained that if $\vtheta_t = \vtheta_{t-1} - h \nabla_{\vtheta}E (\vtheta_{t-1})$ and we want to find $\dot{\vtheta} = -\nabla_{\vtheta} E(\vtheta) + h f_1(\vtheta)$ such that $\norm{\vtheta(h; \vtheta_{t-1}) - \vtheta_t}$ is of order $\mathcal{O}(h^3)$, then $f_1({\color{red}{\vtheta_{t-1}}}) = - \frac{1}{2} \nabla_{\vtheta}^2 E({\color{red}{\vtheta_{t-1}}})\nabla_{\vtheta} E({\color{red}{\vtheta_{t-1}}})$ (Eq (<ref>)). From there, following Barrett and Dherin, 2021, we concluded that $f_1(\vtheta) = - \frac{1}{2} \nabla_{\vtheta}^2 E(\vtheta)\nabla_{\vtheta} E(\vtheta)$. But notice how BEA only sets a constraint on the value of the vector field at the initial point $\vtheta_{t-1}$. If we allow the modified vector field to depend on the initial condition, an equally valid choice for $f_1(\vtheta) = - \frac{1}{2} \nabla_{\vtheta}^2 E(\vtheta)\nabla_{\vtheta} E({\color{red}{\vtheta_{t-1}}})$ or $f_1(\vtheta) =  -\frac{1}{2} \nabla_{\vtheta}^2 E({\color{red}{\vtheta_{t-1}}})\nabla_{\vtheta} E(\vtheta)$. By construction, the above flows will also have an error of $\mathcal{O}(h^3)$ after one gradient descent step of learning rate $h$ with initial parameters $\vtheta_{t-1}$. The latter vector fields only describe the gradient descent update with initial parameters $\vtheta_{t-1}$ and thus they only apply to this specific gradient descent step, though as previously noted that is also the case with the IGR flow due to the dependence on the data batch—see Eq (<ref>). The advantage of constructing modified vector fields that depend on initial parameters lies in the ability to write modified losses when a modified vector field depending only on $\vtheta$ cannot be written as a gradient operator, as we shall see in the next sections.

This observation about BEA leads us to the following remarks:

There are multiple flows that lead to the same order in learning rate error after one discrete update. Many of these flows depend on the initial conditions of the system, i.e. the initial parameters of the discrete update. We visualise this approach in Figure <ref> and contrast it with the existing BEA interpretation.

[Standard interpretation.]
[every text node part/.style=align=center,inner sep=0,outer sep=0, scale=0.9][overlay]
(theta_t_minus_1) at (0,0);
(theta_t) at (2.2,-1);
(theta_t_plus_1) at (4.4, -0.5);

(draw) at ($(theta_t_minus_1) + (+0.3,-0.32)$) $\vtheta_{t-1}$;
(draw) at ($(theta_t) + (-0.2,-0.2)$) $\vtheta_{t}$;
(draw) at ($(theta_t_plus_1) + (-0.05,-0.35)$) $\vtheta_{t+1}$;

(first_time_transition) at ($(theta_t_minus_1) + (+0.55,1.6)$);
(second_time_transition) at ($(first_time_transition) + (3,0)$);

(draw) at (first_time_transition) $t -1 \rightarrow t $;
(draw) at (second_time_transition) $t \rightarrow t + 1$;

(cont_theta_t) at (2.2,0.8);
(cont_theta_t_plus_1) at (4.4, 1.05);

(mod_cont_theta_t) at (2.2, -0.6);
(mod_cont_theta_t_plus_1) at (4.4, -0.2);

[NavyBlue,thick,dashed] (theta_t_minus_1) – (theta_t);
[NavyBlue,thick,dashed] (theta_t) – (theta_t_plus_1);

[OliveGreen,thick]  (theta_t_minus_1) to[out=50,in=180] node[midway,above] $\dot{\vtheta}$ (cont_theta_t);
[OliveGreen,thick]  (theta_t) to[out=50,in=180] node[midway,above]$\dot{\vtheta}$ (cont_theta_t_plus_1);

[Plum,thick]  (theta_t_minus_1) to[out=50,in=180]  node[near end,above] $\dot{\tilde{\vtheta}}$ (mod_cont_theta_t);
[Plum,thick]  (theta_t) to[out=30,in=170] node[near end,above] $\dot{\tilde{\vtheta}}$ (mod_cont_theta_t_plus_1);

[black,thick,dashed] ($(theta_t) + (0,-0.5)$) – ($(cont_theta_t) + (0,0.7)$);

] (theta_t_plus_1) – (mod_cont_theta_t_plus_1)
node [pos=0.5,anchor=west,xshift=0.15cm] $\mathcal{O}(h^3)$;

] (theta_t_plus_1) – (cont_theta_t_plus_1)
node [pos=0.5,anchor=west,xshift=1.6cm] $\mathcal{O}(h^2)$;

[Novel interpretation.]
[every text node part/.style=align=center,inner sep=0,outer sep=0,scale=0.9][overlay]
(theta_t_minus_1) at (0,0);
(theta_t) at (2.2,-1);
(theta_t_plus_1) at (4.4, -0.5);

(draw) at ($(theta_t_minus_1) + (+0.3,-0.32)$) $\vtheta_{t-1}$;
(draw) at ($(theta_t) + (-0.2,-0.2)$) $\vtheta_{t}$;
(draw) at ($(theta_t_plus_1) + (-0.05,-0.35)$) $\vtheta_{t+1}$;

(first_time_transition) at ($(theta_t_minus_1) + (+0.55,1.6)$);
(second_time_transition) at ($(first_time_transition) + (3,0)$);

(draw) at (first_time_transition) $t -1 \rightarrow t $;
(draw) at (second_time_transition) $t \rightarrow t + 1$;

(cont_theta_t) at (2.2,0.8);
(cont_theta_t_plus_1) at (4.4, 1.05);

(mod_cont_theta_t) at (2.2, -0.6);
(mod_cont_theta_t_plus_1) at (4.4, -0.2);

[NavyBlue,thick,dashed] (theta_t_minus_1) – (theta_t);
[NavyBlue,thick,dashed] (theta_t) – (theta_t_plus_1);

[OliveGreen,thick]  (theta_t_minus_1) to[out=50,in=180] node[midway,above] $\dot{\vtheta}$ (cont_theta_t);
[OliveGreen,thick]  (theta_t) to[out=50,in=180] node[midway,above]$\dot{\vtheta}$ (cont_theta_t_plus_1);

[Plum,thick]  (theta_t_minus_1) to[out=50,in=180]  node[near end,above] $\dot{\tilde{\vtheta}}_{\color{red}{t}}$ (mod_cont_theta_t);
[Plum,thick]  (theta_t) to[out=30,in=170] node[near end,above] $\dot{\tilde{\vtheta}}_{{\color{red}{t + 1}}}$ (mod_cont_theta_t_plus_1);

[black,thick,dashed] ($(theta_t) + (0,-0.5)$) – ($(cont_theta_t) + (0,0.7)$);

] (theta_t_plus_1) – (mod_cont_theta_t_plus_1)
node [pos=0.5,anchor=west,xshift=0.15cm] $\mathcal{O}(h^3)$;

] (theta_t_plus_1) – (cont_theta_t_plus_1)
node [pos=0.5,anchor=west,xshift=1.6cm] $\mathcal{O}(h^2)$;

[Visualising the standard approach to backward error analysis alongside an approach which constructs a different flow per gradient descent iteration.]Visualising the standard approach to BEA alongside an approach which constructs a different flow per gradient descent iteration. In our previous use of BEA tkz:orig_bea, we constructed modified flows $\dot{\tilde{\vtheta}}$ to capture the discretisation error of gradient descent; these flows did not depend on the initial iteration parameters. In this chapter, we take the second approach tkz:novel_bea, allowing us to construct additional flows $\dot{\tilde{\vtheta}}_{\color{red}{t}}$, which depend on initial parameters, and showcase additional implicit regularisation effects.

Implicit in the choice of the IGR flow [Barrett and Dherin, 2021], as well as the flows we have introduced thus far in this work using BEA, namely the PF (Chapter <ref>) and the modified flows of two-player games (Chapter <ref>), make an implicit assumption: that we are looking for the modified flows that hold for every training iteration and do not depend on the initial parameters. This is however challenged already in the case of stochastic gradient descent, where each implementation of the flows requires dependence on data batches, as shown in Eq (<ref>).

§ IMPLICIT REGULARISATION IN MULTIPLE STOCHASTIC GRADIENT DESCENT STEPS

We now use the above observations to build modified flows that capture multiple gradient descent updates with an error of $\mathcal{O}(h^3)$ after  $n$ SGD steps. We are interested in modified flows
that can be used to construct modified losses by writing the vector field of the flow as the negative gradient of a function; this enables us to capture the implicit regularisation effects of taking multiple SGD steps.
We analyse $n$ SGD steps, starting at iteration $t$

\begin{align}
\vtheta_{t+\mu} = \vtheta_{t+\mu-1} - \nabla_{\vtheta} E(\vtheta_{t+\mu-1}; \vX^{t+\mu}), \hspace{3em} \mu \in \{0, \dots, n-1\}
\label{eq:sgd_multiple_updates_bea}
\end{align}

where for simplicity we denoted $E(\vtheta_{t-1}; \vX) = \frac{1}{B}\sum_{i=1}^B E(\vtheta_{t-1}; \vx_i)$, with elements $\vx_i$ forming the batch $\vX$.
We start with the following remark:

Denote $E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \}) = \frac{1}{n} \sum_{\mu=0}^{n-1} E(\vtheta; \vX^{t+\mu})$, i.e. the average loss obtained from the $n$ data batches. Then the trajectory obtained by taking $n$ steps of stochastic gradient descent follows the trajectory of minimising the loss in continuous-time with $\mathcal{O}(h^2)$ error:

\begin{align}
\tilde{E}(\vtheta) &= E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \}).
\end{align}

Perhaps surprisingly, the above shows that effects from mini-batches appear only at higher-order terms in learning rate $h$. Thus, to capture implicit regularisation effects that account for the mini-batches used in the SGD updates, we use BEA. We find a modified flow that describes the SGD update with an error of  $\mathcal{O}(h^3)$,  with a vector field that can be written as a negative gradient. We provide proofs and the flows that construct the regularisers in the Appendix <ref>.
We use the same steps used for our other BEA proofs: 1) we expand the $n$ discrete updates in Eq (<ref>) to write a relation between the parameters at time $t-1$ and $t+n-1$ up to second order in learning rate; 2) expand the changes in continuous-time up to $\mathcal{O}(h^2)$; and 3) match the $\mathcal{O}(h^2)$ terms to find the terms which quantify the drift.
The significant difference with our previous approaches is that we allow the vector field to depend on the initial parameters, and choose a vector field that can be written as a gradient in order to construct a modified loss.

Denote $E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \}) = \frac{1}{n} \sum_{\mu=0}^{n-1} E(\vtheta; \vX^{t+\mu})$, i.e. the average loss obtained from the $n$ data batches. Then the trajectory obtained by taking $n$ steps of stochastic gradient descent follows the trajectory of minimising the loss in continuous-time with $\mathcal{O}(h^3)$ error

\begin{align}
\tilde{E}(\vtheta) &= E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \})
+ \frac{n h}{4} \underbrace{\norm{ \nabla_{\vtheta}E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \})}^2}_{\text{full batch norm regularisation}} \label{eq:sgd_norm_min}\\
& \quad - \frac{h}{n}\sum_{\mu =1 }^{n-1} \underbrace{\left[ \nabla_{\vtheta} E(\vtheta; \vX^{t+\mu})^T \left(\sum_{\tau = 0}^{\mu-1}    \nabla_{\vtheta}E({\color{red}{\vtheta_{t  -1}}}; \vX^{t+\tau})\right)\right]}_{\text{mini-batch gradient alignment}}.
\label{eq:modified_sup_learning}
\end{align}

We thus find the implicit regularisation effects induced by $n$ steps of SGD, capturing the importance of exact batches used and their order in Eq (<ref>). We note that without making use the observations regarding BEA in the previous section, and thus without the parameters $\vtheta_{t -1}$, a modified loss could not have been constructed outside of the full-batch case where one can recover the IGR flow.
The following remark immediately follows by setting $n=2$ in Eq (<ref>):
When taking a second stochastic gradient descent step, there is an implicit regularisation term maximising the dot product between the gradient at the current step and the gradient at the previous step: $\nabla_{\vtheta} E^T \nabla_{\vtheta} E(\vtheta_{t-1})$. This can be achieved by aligning the direction of the gradients between the two iterations or increasing the gradient norm, but we note that increasing gradient norm is counter other implicit regularisation effect in Eq (<ref>).

We now compare this novel modified loss with the modified loss we obtained by ignoring stochasticity and assuming all $n$ updates have been done with a full-batch; this entails using the IGR loss (proof for multiple steps in Section <ref>)

\begin{align}
\tilde{E} = E(\vtheta;\{\vX^{t}, \dots, \vX^{t+n-1} \}) + \frac{h}{4} \|\nabla_{\vtheta} E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \}) \|^2 .
\label{eq:igr_flow_mul}
\end{align}

The above modified losses show that both full-batch gradient descent and SGD have a pressure to minimise the gradient norm  $\|\nabla_{\vtheta} E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \})\|$. SGD leads to an additional regularisation effect capturing the importance of the order in which mini-batches are presented in training: maximising the dot product between gradients computed at the current parameters given a batch and the gradients computed at the initial parameters for all batches presented before the given batch. While this can be achieved both by increasing the norm of the gradients or by aligning gradients with those at the initial iteration, we note that increasing the gradient norm is counter to the other regulariser induced by SGD, the gradient norm  minimisation effect shown in Eq (<ref>).

The important role of learning rates in BEA appears here too: for our approximations to hold $nh$ has to be sufficiently small. If we adjust the learning rate by the number of updates, i.e. if we set the learning rate for SGD equal to $h/n$, where $h$ is the learning rate used by full-batch gradient descent, we obtain the same implicit regularisation coefficient  for the gradient norm $\|\nabla_{\vtheta} E(\vtheta;\{\vX^{t}, \dots, \vX^{t+n-1} \}) \|$ minimisation as the IGR flow in Eq <ref>, and the main difference between the two modified losses is given by the mini-batch gradient alignment terms present in Eq (<ref>).

The number of updates, $n$, also plays an important role. While the mini-batch gradient regularisation term in Eq (<ref>) has a coefficient of $\frac{h}{n}$, there are $\frac{n(n-1)}{2}$ terms composing the term. Thus the magnitude of the mini-batch gradient regularisation term can grow with $n$, but its effects strongly depend on the distribution of the gradients computed at different batches. For example, if gradients $\nabla_{\vtheta} E(\vtheta_{t-1};\vX^{t +i})$ are normally distributed with the mean at the full-batch gradient, i.e. $\nabla_{\vtheta} E(\vtheta_{t-1};\vX^{t +i}) \sim \mathcal{N}(\nabla_{\vtheta} E(\vtheta_{t-1}), \sigma^2)$, as the number of updates grows the regularisation effect in Eq <ref> will result in a pressure to align mini-batch gradients with the full-batch gradient at the initial parameters, $\nabla_{\vtheta} E(\vtheta_{t-1})$. Since our results hold for multiple values of $n$, empirical assessments need to be made to understand the interplay between the number of updates and the strength of  mini-batch gradient regularisation on training.

The approach provided here is complementary to that of that of Smith et al., 2021, who obtain a relationship similar to Eq (<ref>) but in expectation—they describe expected value of the modified loss $\mathbb{E}_{\sigma} \left[E_{sgd}(\vtheta; \{\vX^{\sigma(t)}, \dots, \vX^{\sigma(t+n-1)} \}) \right]$ where the expectation is taken over all possible data batch shufflings $\sigma$, but not the elements in the batch.
As before we denote $E(\vtheta; \vX^{k})$ as the loss given by mini-batch $k$ (the equivalent to $\hat{C}_k$ in their notation, see their Eqs (1) and (2)), and write the modified loss they obtain in expectation as

\begin{align}
\mathbb{E_\sigma} \left[E_{sgd}(\vtheta) \right] &= E(\vtheta;\{\vX^{0}, \dots, \vX^{n-1} \}) + \frac{h}{4n} \sum_{k=0}^{n-1} \norm{\nabla_{\vtheta} E(\vtheta;\vX^k)}^2 \\&= E(\vtheta; \{\vX^{0}, \dots, \vX^{n-1} \}) + \frac{h}{4}\|\nabla_{\vtheta} E(\vtheta; \{\vX^{0}, \dots, \vX^{n-1} \}) \|^2 \\ &\hspace{1em} + \frac{h}{4n} \sum_{k=0}^{n-1} \norm{\nabla_{\vtheta} E(\vtheta;\vX^k) - \nabla_{\vtheta} E(\vtheta; \{\vX^{0}, \dots, \vX^{n-1} \}}^2.
\end{align}

Our proposed approach does not require working in expectations and accounts for the exact batches sampled from the dataset; thus it directly describes SGD as used in practice.
If we make the same assumptions as Smith et al., 2021 and take expectation over all possible batch shufflings in an epoch in Eq (<ref>) we obtain (proof in Section <ref>):

\begin{align}
\hspace{-0.7em}\mathbb{E}_{\sigma} \left[E_{sgd}(\vtheta) \right]
&=  E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \}) \\&\hspace{1em} + \frac{nh}{4}  \norm{\nabla_{\vtheta} E(\vtheta;\{\vX^{t}, \dots, \vX^{t+n-1} \})}^2
\\&\hspace{1em}- \frac{h}{2n} \nabla_{\vtheta} E(\vtheta; \{\vX^{t}, \dots, \vX^{t+n-1} \})^T \nabla_{\vtheta} E({\color{red}\vtheta_{t-1}}; \{\vX^{t}, \dots, \vX^{t+n-1} \})
\\&\hspace{1em}  -  \frac{h}{2n}\left[\sum_{k=0}^{n-1}  \nabla_{\vtheta} E(\vtheta;\vX^{t+k}) ^T \nabla_{\vtheta} E({\color{red}\vtheta_{t-1}}; \vX^{t+k})\right].
\end{align}

As expected, we obtain a different result than that of Smith et al., 2021: while the gradient norm minimisation is still present in our results, the per-batch regularisation in their formulation gets translated into a dot product term, where both the full-batch and per-batch gradients are regularised to be aligned with those at the beginning of the epoch.  Both approaches share the limitation that $nh$ needs to be suitably small for approximations to be relevant; we note that since we do not require $n$ to be the number of updates in an epoch—and we obtain interestingly regularisation effects for $n=2$, see Remark <ref>—this is less of an issue for our approach. Their approach has the advantage of finding an implicit regularisation effect that does not depend on the initial parameters. By depending on initial parameters, however, our approach does not require working in expectations and accounts for the exact batches used in the SGD updates. We hope this can be used to stabilise SGD over multiple steps, as we have done with a single GD step in Section <ref>, that it can be used in continual and transfer learning [Iman et al., 2022, Zhuang et al., 2020, Parisi et al., 2019], as well as understanding the effects of the order of examples on model optimisation in online learning.

§ IMPLICIT REGULARISATION IN GENERALLY  DIFFERENTIABLE TWO-PLAYER GAMES

In Chapter <ref>, we presented a framework for the quantification of discretisation drift in two-player games.  We found distinct modified flows that  describe simultaneous Euler updates and alternating Euler updates.
In both cases, we found corrections $f_1$ and $g_1$ to the original system

\begin{align}
\dot{\vphi} &=  f( \vphi, \vtheta) \\
\dot{\vtheta}  &= g( \vphi, \vtheta),
\end{align}

such that the modified continuous system

\begin{align}
\dot{\vphi} &=  f( \vphi, \vtheta)  + h f_1( \vphi, \vtheta) \label{eq:general_bae_formulation1}\\
\dot{\vtheta}  &= g( \vphi, \vtheta) + h g_1( \vphi, \vtheta), \label{eq:general_bae_formulation2}
\end{align}

follows the discrete steps of the method with a local error of order $\mathcal O(h^3)$.
More precisely, if $(\vphi_{t}, \vtheta_{t})$ denotes the discrete step of the method at time $t$ and $( \vphi(h),  \vtheta(h))$ corresponds to the continuous solution of the modified system above starting at $(\vphi_{t-1}, \vtheta_{t-1})$,  $\| \vphi_{t} -  \vphi( h) \| \textrm{ and }  \| \vtheta_{t} -  \vtheta( h) \|$
are of order $\mathcal O(h^3)$. In this section, we assume for simplicity that both players use the same learning rate $h$ and simultaneous updates, but the same arguments can be made when they use different learning rates or alternating updates.

Using this framework, we constructed modified loss functions in the case of zero-sum (Section <ref>) and common-payoff games (Section <ref>). However, using the aforementioned modified flows, we cannot always write the vector fields of the modified flows as a gradient for differentiable two-player games and thus we cannot construct modified losses. To see why, consider the case where we have two loss functions for the two players respectively $E_{\vphi}(\vphi, \vtheta): \mathbb{R}^m\times \mathbb{R}^n \rightarrow \mathbb{R}$ and $E_{\vtheta}(\vphi, \vtheta): \mathbb{R}^m\times \mathbb{R}^n \rightarrow \mathbb{R}$.
This leads to the update functions $f = -\nabla_{\vphi} E_{\vphi}$ and $g = -\nabla_{\vtheta} E_{\vtheta}$.
By using $f = -\nabla_{\vphi} E_{\vphi}$ and $g = -\nabla_{\vtheta} E_{\vtheta}$ in Theorem <ref>, we have that for simultaneous gradient descent

\begin{align}
f_1 &= - \frac{1}{2}  \jacphif f- \frac{1}{2} \jacthetaf g\\
&= \underbrace{- \frac{1}{4}  \nabla_{\vphi} \norm{ \nabla_{\vphi} E_{\vphi}}^2}_{\text{self term}} \underbrace{- \frac{1}{2} \jacparam{\vtheta}{\nabla_{\vphi} E_{\vphi}} \nabla_{\vtheta} E_{\vtheta}}_{\text{interaction term}}
\end{align}

and similarly

\begin{align}
g_1 &= \underbrace{- \frac{1}{4}  \nabla_{\vtheta} \norm{ \nabla_{\vtheta} E_{\vtheta}}^2}_{\text{self term}}  \underbrace{- \frac{1}{2}
\jacparam{\vphi}{\nabla_{\vtheta} E_{\vtheta}} \nabla_{\vphi} E_{\vphi}}_{\text{interaction term}}.
\end{align}

Thus, it is not always possible to write $f_1$ and $g_1$ as gradient functions, since the interaction terms—see Definition <ref> for the definition of self and interaction terms—cannot always be written as a gradient. Thus, the modified flows cannot be used to construct modified losses leading to implicit regularisers in generally differentiable two-player games, as we have done for zero-sum games in Chapter <ref>.

We now use the interpretation of BEA we provided in this chapter, to choose two other functions $f_1$ and $g_1$ and we will use them to construct another set of modified flows. These flows will still satisfy the value constraints required by the BEA proofs and by construction after one gradient descent update the error between the modified flows and the discrete updates is $\mathcal{O}(h^3)$. Indeed, what the BEA proofs provide for simultaneous gradient descent is (see Eq (<ref>) in the Appendix)

\begin{align}
f_1(\vphi_{t-1}, \vtheta_{t-1}) &= - \frac{1}{2}  \jacphif(\vphi_{t-1}, \vtheta_{t-1}) f(\vphi_{t-1}, \vtheta_{t-1}) - \frac{1}{2}  \jacthetaf(\vphi_{t-1}, \vtheta_{t-1}) g(\vphi_{t-1}, \vtheta_{t-1}) \\
g_1(\vphi_{t-1}, \vtheta_{t-1}) &= - \frac{1}{2} \jacphig(\vphi_{t-1}, \vtheta_{t-1}) f(\vphi_{t-1}, \vtheta_{t-1}) - \frac{1}{2}  \jacthetag(\vphi_{t-1}, \vtheta_{t-1}) g(\vphi_{t-1}, \vtheta_{t-1}).
\end{align}

From here, we can choose $f_1$ and $g_1$, now depending on the iteration number $t$:

\begin{align}
f_{1, t}(\vphi, \vtheta) &= - \frac{1}{2}  \jacphif(\vphi, \vtheta) f(\vphi, \vtheta) - \frac{1}{2}  \jacthetaf(\vphi, \vtheta) g({\color{red}\vphi_{t-1}, \vtheta_{t-1}}) \\
g_{1, t}(\vphi, \vtheta) &= - \frac{1}{2}  \jacphig(\vphi, \vtheta) f({\color{red}\vphi_{t-1}, \vtheta_{t-1}}) - \frac{1}{2}  \jacthetag(\vphi, \vtheta) g(\vphi, \vtheta).
\end{align}

Here, we treat $g(\vphi_{t-1}, \vtheta_{t-1})$ and $f(\vphi_{t-1}, \vtheta_{t-1})$ as constants. We replace $f = - \nabla_{\vphi} E_{\vphi}$ and $g = - \nabla_{\vtheta} E_{\vtheta}$  in $f_{1, t}$ and $g_{1, t}$ and write the drift terms as gradient functions:

\begin{align}
f_{1, t}(\vphi, \vtheta) &= - \frac{1}{2} \jacparam{\vphi}{\nabla_{\vphi} E_{\vphi}} \nabla_{\vphi} E_{\vphi}  - \frac{1}{2} \jacparam{\vtheta}{\nabla_{\vphi} E_{\vphi}} \nabla_{\vtheta} E_{\vtheta}(\vphi_{t-1}, \vtheta_{t-1}) \\
&= - \frac{1}{2} \nabla_{\vphi} \norm{\nabla_{\vphi}E_{\vphi}}^2 - \frac{1}{2} \nabla_{\vphi}\left(\nabla_{\vtheta}  E_{\vphi} ^T \nabla_{\vtheta} E_{\vtheta}(\vphi_{t-1}, \vtheta_{t-1}) \right) \\
&= - \nabla_{\vphi} \left(\frac{1}{2} \norm{\nabla_{\vphi}E_{\vphi}}^2 + \frac{1}{2} \nabla_{\vtheta}  E_{\vphi}^T  \nabla_{\vtheta} E_{\vtheta}(\vphi_{t-1}, \vtheta_{t-1}) \right) \\
g_{1, t}(\vphi, \vtheta) &= - \nabla_{\vtheta} \left(\frac{1}{2} \norm{\nabla_{\vtheta}E_{\vtheta}}^2 + \frac{1}{2} \nabla_{\vphi}  E_{\vtheta}^T \nabla_{\vphi} E_{\vphi}(\vphi_{t-1}, \vtheta_{t-1}) \right).
\end{align}

Replacing the above in the modified flows in Eqs (<ref>) and (<ref>)), we obtain

\begin{align}
\dot{\vphi} &=  f( \vphi, \vtheta)  + h f_1( \vphi, \vtheta) \\
&= - \nabla_{\vphi} \left(E_{\vphi} +  h \left(\frac{1}{2} \norm{\nabla_{\vphi}E_{\vphi}}^2 + \frac{1}{2} \nabla_{\vtheta}  E_{\vphi}^T \nabla_{\vtheta} E_{\vtheta}(\vphi_{t-1}, \vtheta_{t-1}) \right)   \right) \label{eq:disc_int_app} \\
\dot{\vtheta}  &= g( \vphi, \vtheta) + h g_1( \vphi, \vtheta) \\
&= - \nabla_{\vtheta} \left(E_{\vtheta} +  h \left(\frac{1}{2} \norm{\nabla_{\vtheta}E_{\vtheta}}^2 + \frac{1}{2} \nabla_{\vphi}  E_{\vtheta} ^T \nabla_{\vphi} E_{\vphi}(\vphi_{t-1}, \vtheta_{t-1}) \right)\right).
\end{align}

We can now write modified losses for each gradient descent iteration of a general two-player differentiable game (which depend on the iteration $t$), which describe the local trajectory of simultaneous gradient descent up to $\mathcal{O}(h^3)$:

\begin{align}
\tilde{E}_{\vphi, t} = E_{\vphi} +  h \left(\underbrace{\frac{1}{2} \norm{\nabla_{\vphi}E_{\vphi}}^2}_{\text{self term}} + \underbrace{\frac{1}{2} \nabla_{\vtheta}  E_{\vphi}^T \nabla_{\vtheta} E_{\vtheta}({\color{red}\vphi_{t-1}, \vtheta_{t-1}}) }_{\text{interaction term}} \right) \label{eq:int_term_d_first} \\
\tilde{E}_{\vtheta, t} = E_{\vtheta} +  h \left(\underbrace{\frac{1}{2} \norm{\nabla_{\vtheta}E_{\vtheta}}^2}_{\text{self term}} + \underbrace{\frac{1}{2}  \nabla_{\vphi}  E_{\vtheta}^T \nabla_{\vphi} E_{\vphi}({\color{red}\vphi_{t-1}, \vtheta_{t-1}}}_{\text{interaction term}})\right).
\end{align}

The self terms result in the implicit gradient regularisation found in supervised learning [Barrett and Dherin, 2021] and  zero-sum games (Chapter <ref>): each player has an incentive to minimise its own gradient norm. The interaction term for each player encourages minimising the dot product between the gradients of its loss with respect to the other player's parameters and the previous gradient update of the other player.
Consider the first player's interaction term, equal to $(-\nabla_{\vtheta} E_{\vphi}(\vphi_t, \cdot))^T (- \nabla_{\vtheta}E_{\vtheta}({\color{red}\vphi_{t-1}, \vtheta_{t-1}}))$, where $-E_{\vtheta}({\color{red}\vphi_{t-1}, \vtheta_{t-1}})$ is the previous update direction of $\vtheta$ aimed at minimising $E_{\vtheta}$.
Its implicit regularisation effect depends on the functional form of $E_{\vphi}$ and $E_{\vtheta}$: if $\nabla_{\vtheta} E_{\vtheta}$ and $\nabla_{\vtheta} E_{\vphi}$ have aligned directions by construction, then the implicit regularisation effect nudges the first player's update towards a point in space where the second player's update changes direction; the opposite is true if $\nabla_{\vtheta} E_{\vtheta}$ and $\nabla_{\vtheta} E_{\vphi}$ are misaligned by construction.

§.§ The effect of the interaction terms: a GAN example

We work through the regularisation effect of interaction terms for a pair of commonly used GAN losses that do not form a zero-sum game (using the generator non-saturating loss [Goodfellow et al., 2015]) and contrast it with the zero-sum case (the saturating loss); we have previously investigated these losses in Chapter <ref>. As before, we denote the first player, the discriminator, as $D$, parametetrised by $\vphi$, and the generator as $G$, parametrised by $\vtheta$. We denote the data distribution as $p^*(\vx)$ and the latent distribution $p(\vz)$.
Given the non-saturating GAN loss Goodfellow et al., 2014:

\begin{align}
E_{\vphi}(\vphi, \vtheta) &= \mathbb{E}_{p^*(\vx)} \log D(\vx; \vphi) + \mathbb{E}_{p(\vz)} \log (1 - D(G(\vz; \vtheta); \vphi)) \\
E_{\vtheta}(\vphi, \vtheta) &=  \mathbb{E}_{p(\vz)} - \log D(G(\vz; \vtheta); \vphi),
\end{align}

we can then write the interaction term $\nabla_{\vtheta}  E_{\vphi}^T \nabla_{\vtheta} \nabla_{\vtheta} E_{\vtheta}({\color{red}\vphi_{t-1}, \vtheta_{t-1}})$ in Eq (<ref>) at iteration $t$  as (derivation in Appendix <ref>)

\begin{align}
\frac{1}{B^2}\sum_{i,j=1}^B &c^{non-sat}_{i,j}  \nabla_{\vtheta}  D(G(\vz^i_t; \vtheta); \vphi)^T   \nabla_{\vtheta}  D(G(\vz^j_{t-1}; {\color{red}\vtheta_{t-1}}); {\color{red}\vphi_{t-1})}), \hspace{2em} \text{with} \\
&c^{non-sat}_{i,j} = \frac{1}{1 - D(G(\vz^i_t; \vtheta); \vphi)} \frac{1}{D(G(\vz^j_{t-1}; {\color{red}\vtheta_{t-1}}); {\color{red}\vphi_{t-1})}} \label{eq:c_i_j_non_sat},
\end{align}

where $\vz^i_t$ is the latent variable with index $i$ in the batch at time $t$ with batch size $B$.
Thus, the strength of the regularisation—$c^{non-sat}_{i,j}$—depends on how confident the discriminator is. In particular, this implicit regularisation encourages the discriminator update into a new set of parameters where the gradient $\nabla_{\vtheta}  D(G(\vz^i_t; \vtheta); \vphi)$ points away from the direction of  $\nabla_{\vtheta}  D(G(\vz^j_{t-1}; \vtheta_{t-1}); \vphi_{t-1})$ when $c_{i,j}$ is large. This occurs for $\vz^i_t$ where the discriminator is fooled by the generator—i.e. $1 - D(G(\vz_i; \vtheta); \vphi)$ is close to 0 and  $\frac{1}{1 - D(G(\vz^i_t; \vtheta); \vphi)}$ is large—and samples $\vz^j_{t-1}$ where the discriminator was correct at the previous iteration—$D(G(\vz^j_{t-1}; \vtheta_{t-1}); \vphi_{t-1})$ low and thus $\frac{1}{D(G(\vz^j_{t-1}; \vtheta_{t-1}); \vphi_{t-1})}$ is large. This can be seen as beneficial regularisation for the generator, by ensuring the update direction $- {\nabla_{\vtheta} E_{\vtheta} = \mathbb{E}_{p(\vz)} \frac{1}{D(G(\vz; \vtheta); \vphi)} \nabla_{\vtheta}  D(G(\vz; \vtheta); \vphi)}$ is adjusted accordingly to the discriminator's output. We note, however, that regularisation might not have a strong  effect, as gradients $\nabla_{\vtheta}  D(G(\vz; \vtheta); \vphi)$ that have a high weight in the generator's update are those where the discriminator is correct and classifies generated data as fake—i.e. $\frac{1}{D(G(\vz; \vtheta); \vphi)}$ is large—but there the factor $\frac{1}{1 - D(G(\vz; \vtheta); \vphi)}$ in the interaction term coefficient in Eq (<ref>) will be low. This is inline with our empirical results in Section <ref>, where we saw that for the non-saturating loss discretisation drift does not have a strong effect on performance.

We can contrast this with what regularisation we obtain when using the saturating loss [Goodfellow et al., 2014], where $E_{\vtheta}(\vphi, \vtheta) = - \mathbb{E}_{p(\vz)} \log (1 - D(G(\vz; \vtheta); \vphi))$. We used the saturating loss extensively in experiments in Section <ref>, and we have shown it performs poorly in comparison to the non-saturating loss, in-line with the literature [Goodfellow et al., 2014]. If we perform the same analysis as above for the saturating loss we obtain the following implicit regulariser

\begin{align}
\frac{1}{B^2} \sum_{i,j=1}^B &c^{sat}_{i,j} \nabla_{\vtheta}  D(G(\vz^i_t; \vtheta); \vphi)^T   \nabla_{\vtheta}  D(G(\vz^j_{t-1}; {\color{red}\vtheta_{t-1}}); {\color{red}\vphi_{t-1})}),  \hspace{3em} \text{with}\\
&c^{sat}_{i,j} = \frac{1}{1 - D(G(\vz^i_t; \vtheta); \vphi)} \frac{1}{1-D(G(\vz^j_{t-1}; {\color{red}\vtheta_{t-1}}); {\color{red}\vphi_{t-1})})} \label{eq:c_i_j_sat}.
\end{align}