Keble College Doctor of Philosophy Hilary 2021

# Geometric Optimisation on
Manifolds with Applications to
Deep Learning

Mario Lezcano Casado

###### Acknowledgements.

I would like to begin by thanking my supervisors, Prof. Raphael Hauser and
Prof. Vidit Nanda for their time, effort, and guidance throughout this thesis.
Thank you for giving me the flexibility to delve into whatever mad topic I
found interesting at the time. Discussions with Jaime Mendizabal and David
Martínez-Rubio have been truly excellent. Their suggestions, clarifications,
and probing questions pushed me to think in new ways and greatly improved the
thesis. My collaborators Adam Golinski and Tom Rainforth helped me to explore
the applications that this research has in other fields. Prof. Jose Manuel
Gamboa Mutuberría has been a constant source of support and encouragement all
these years. To all of them a warm thank you. A heartfelt thank you to my
parents and my brother for their unconditional support, in particular during
the final writing of this thesis. On the less academic side, I would like to
thank everybody in and around La Katacumbia. Borondiño, Juan, Adri, Álvaro,
Jorge, Recio, Miguel, David vinos, Sam, Felix, Steve, Adam, Sergio, Daniel,
María, Laura and everyone at $11$ Iffley road… You have made these years in
the small city of Oxford incredibly enjoyable. From the parties by the river,
to nights at Plush and Cellar, barbecues, swimming pool, squash, gigs at the
Tap Social and even a very special summer solstice in Stonehenge followed by a
weekend in Bristol that was not for the faint-hearted. In all fairness, a long
list of things to have done in the small town of Oxford. I must single out
Garima, Adri, and David, with whom I spent the best of lockdowns. Starting as
an odd night in London that looked straight from a B film, it was followed by
some of the best food I’ve had in my life—the arroz con leche, risotto, and
banana bread come to mind—Smash-up, Hollow Knight, Cuphead, the month that we
worked out… For these things and many more, thank you for making this surreal
year special. I would also like to thank Amigas y Conocidas, Roberto, Lorenzo,
Cámara, Charles, Visierdo, Inés and those floating around it such as my
beloved Juan. It is always fantastic to be back and be able to forget whatever
I’m working on and just have a great night out, whether in Burgos or in
Berlin. There are a few more that have made the last few years enjoyable—Irene
and Elyem from la Caixa, Mendi, Paco, David, Aguirre, Aitor and Adrià from
Comeis, Guillem and Sara, Debbie, Lina, Elena, Noémie, Manolo…—with whom I
have shared some truly amazing moments. To all of you: Thank you.

Contents

toc

##  Chapter 1 Introduction

###  1.1 Why this Thesis?

We start this thesis by justifying the theoretical and practical approach that
this thesis takes. We outline the general problems that we study on throughout
it. After this, we summarise the results in the thesis and its structure.

####  1.1.1 A practical reason

Optimisation has always been a field with vast applications in the real world,
bringing together the fields of analysis and constructive mathematics—and
subsequently computer science—in order to be able to approximate in a
principled way the minima of functions traditionally defined on some subset of
the Euclidean space.

A natural extension of this setting is that of optimisation on Riemannian
manifolds. Most of the interesting sets of matrices used in the context of
optimisation turn out to have a manifold structure. Optimisation on manifolds
is both theoretically and practically challenging due to the inherent
complexity of the objects involved. Even then, optimisation on matrix
manifolds has proven to be rather useful in many subfields of machine learning
and neural networks (nn). Examples of interesting matrix manifolds in the
context of gradient-based optimisation are the set of positive definite
matrices in Bayesian statistics as kernels [148], orthogonal matrices within
recurrent neural networks [8, 75, 113], nns with structured linear layers via
matrix factorisations such as the qr or the svd decomposition [183, 100], or
invertible matrices in normalising flows [19] and vaes [167].

In parallel, during the last $10$ years, neural networks have revolutionised
many fields such as speech recognition [82, 155], image recognition [106, 87],
and natural language processing [44, 143]. Initially, neural networks were a
concatenation of linear layers, but with time, these architectures have been
getting more and more complex with the introduction of architectures such as
attention mechanisms [15, 172], neural Turing machines [69], variational auto
encoders (vae) [101] and generative adversarial networks (gan) [68].

This trend of building larger models on more advanced hardware, such as gpus
and tpus has gone hand in hand with the development of large linear algebra
libraries such as PyTorch, Tensorflow, or jax. These libraries are similar to
Numpy or Scipy, but also come with full gpu and tpu support, an
autodifferentiation engine, support for automatic mixed-precision (amp), and
many other capabilities to help develop modern high-dimensional deep learning
models that scale to datasets with millions of datapoints.

At the time of writing, there is no library to put these two trends together
in a practical way. The current libraries that implement optimisation on
manifolds techniques are designed for researchers, not for practitioners. For
this reason, the thesis aims to fulfil the following two objectives:

###### GeoTorch: A library for optimisation on manifolds at scale

We design and implement a Python library to help the non-expert using all
these powerful tools in a way that is efficient, extensible, and simple to
incorporate into the workflow of the data scientist, practitioner, and applied
researcher. The algorithms implemented in this library have been designed with
usability and gpu efficiency in mind, and they can be added to any PyTorch
model with just one extra line of code.

###### Orthogonal rnns

We showcase the effectiveness of these tools on an application of optimisation
on manifolds in the setting of time series analysis. In this setting,
orthogonal and unitary optimisation is used to constraint and regularise
recurrent models and avoid _vanishing_ and _exploding gradient_ problems. The
algorithms designed for GeoTorch allow us to achieve state of the art results
in the standard tests for this family of models.

####  1.1.2 A theoretical reason

The current approach to optimisation on manifolds rarely makes use of the
abstract approach, favouring the embedded setting. This seems like a natural
thing to do, as the manifolds of interest in Riemannian optimisation are, by
definitions, manifolds embedded in $\mathbb{R}^{n\times k}$—i.e., matrix
manifolds. When one has a matrix manifold, it is natural to use the structure
of the manifold as a manifold embedded in a linear space and identify points
on the manifold, vectors, actions by a (matrix) Lie group, isometries and
similar objects as matrices or linear functions. This approach often puts
aside the entire abstract toolbox from modern differential geometry, like
fibred manifolds, group actions, principal connections, or Riemannian
homogeneous spaces, in favour of matrix manipulations.

Using matrices and algebraic manipulations rather than sections of the
different associated bundles and group actions often hides the actual
symmetries that make these manifolds tractable. It is these symmetries that
allow us to find closed-form solutions to some differential equations on the
manifold. Cartan made a similar point when talking about the abstract index
formalism:

> The great opportunities that the _absolute differential calculus_ of Ricci
> and Levi-Civita has given us must not prevent us from avoiding purely formal
> calculations, where the profusion111It is of historical interest to note
> that this expression was popularised in English by Spivak, who jokingly
> translated it as “debauch of indices” in the first volume of his treaty on
> differential geometry. The word “débauches” in French has the same meaning
> as “debauch” in English in almost any situation, being this one of the few
> exceptions in which it simply means “profusion” with no extra connotation.
> of indices masks an often very simple geometric reality. It is this reality
> that I set out to highlight here. (Original in French222Les services
> éminents qu’a rendus et que rendra encore le Calcul différentiel absolu de
> Ricci et Levi-Civita ne doivent pas nous empêcher d’éviter les calculs trop
> exclusivement formels, où les débauches d’indices masquent une réalité
> géométrique souvent très simple. C’est cette réalité que j’ai cherché à
> mettre partout en évidence.) [36]

Identifying geometric objects on the manifold with matrices simplifies the
background necessary to check whether a proof is true, but it often hides the
true reasons why the proof is true. Drawing similarities with type systems in
computer science, we can see that while it is true that dropping the types of
a correct program does not make it incorrect, it certainly makes it more
difficult to reason about. The same happens if we simply write points on the
manifold or sections on the tangent bundle of the manifold as matrices: It
hides the geometric meaning from these operations.

In this sense, this more abstract approach, whenever applicable, can be used
to explain the nature of some numerical computations. In particular, in the
field of constrained optimisation, when the set of constraints is a Riemannian
manifold with enough symmetries, we are in the best of scenarios as we can use
all the tool set from differential, Riemannian, and comparison geometry to
explain the behaviour of our algorithms. In this sense, this helps us in
getting closer to the goal already set by Hamming in his book on numerical
methods

> The purpose of computing is insight, not numbers. [73]

With these points in mind, in this thesis we will use the language of
differential geometry for studying the field of optimisation on manifolds. In
particular, we will explore the following two research directions:

###### Use abstract differential geometry to study optimisation problems over
matrices

We use the more abstract language of global differential geometry and group
actions to prove convergence bounds on certain optimisation problems on
manifolds, postponing the instantiation of the results on concrete manifolds
as much as possible.

In doing so, we follow the approach started by Smith in his PhD thesis [160],
which was written aimed for a numerical analysis audience in the influential
paper [51].

###### Use comparison geometry to study optimisation problems on Riemannian
manifolds

We use tools from comparison geometry to give bounds on quantities that are of
interest in optimisation problems. In particular, we build on the work of [98]
to give explicit bounds on the norm of the second derivative of the Riemannian
exponential.

###  1.2 Contributions and Outline

This thesis puts together three fields that are not usually found together.
These are global differential geometry, optimisation on manifolds, and deep
learning. Researchers in each of these fields approach problems lying at the
intersection with completely different tools, notation, and goals in mind. It
is for this reason that we include at the beginning of the thesis two
introductory chapters, one on classical optimisation, machine learning, and
their ties to deep learning (Chapter 2), and another one on differential
geometry from a global perspective (Chapter 3). All the results in these
chapters are well-known to researchers from the fields of optimisation,
machine learning and/or differential and Riemannian geometry, so it is only
the exposition and comments throughout it that are original.

In the last three chapters, we present the research carried by the author
during these last years. All the results in these chapters are original except
when specified otherwise. When a section is well-known to geometers, we
carefully announce so either at the beginning of each section or before each
result so that it is perfectly clear what is classic and what is new.

Chapter 4. Its main theme is optimisation on manifolds, and its secondary
theme is differential and Riemannian geometry. In this chapter, we study how
the theory of pullbacks along submersions is a fruitful one (Section 4.3).
This theory helps to simplify problems in optimisation from spaces with a
difficult geometry to simpler ones. We also show how to use the theory laid
out in Chapter 3 to prove second order bounds on the pullback of a function by
a Riemannian submersion (Theorems 4.3.1 and 4.3.2). We then go and look at the
special case when we pullback an optimisation problem to a Euclidean space via
a retraction. We call this the static trivialisation framework. We make
special emphasis on the case when we use the Riemannian exponential map to
pullback the problem to a tangent space (Section 4.4.2), as this will be the
main tool studied and used in Chapter 5 and implemented in Chapter 6. At the
end, we introduce the dynamic trivialisation framework (Algorithm 1), which
sees the Riemannian exponential as a map from $TM$ to $M$ and uses the vector
bundle structure of $TM$ to solve the shortcomings of static trivialisations.
These ideas were was presented at ICML 2019 and NeurIPS 2019 [113, 112].

Chapter 5. Its main theme is Riemannian and comparison geometry, and its
secondary theme is optimisation on manifolds. In this chapter, we give bounds
on $\lVert\operatorname{Hess}\lparen f\circ\exp_{p}\rparen\rVert$ in terms of
first and second order bounds on $f$ and the curvature of the Riemannian
manifold $(M,\textsl{g})$. We start the chapter by doing a literature review
of previous approaches to the subject, where we show how these bounds can be
readily used in a number of papers that use second order bounds on
retractions. We then offer a self-contained proof of Rauch’s theorem, as we
heavily rely on it throughout the second part of the chapter, and it is not
very well-known in the field of optimisation. We also look at its implications
in the area of convex optimisation on manifolds at the end of the section.
After this, we use these results and a few more to develop novel second order
bounds on the Riemannian exponential (Theorems 5.4.13 and 5.4.14). From the
first and second order bounds, we give conditional convergence rates for sgd
on the function $f\circ\exp_{p}$ (Theorem 5.6.2), and the equivalent result
for the dynamic trivialisation algorithm (Theorem 5.6.4). The second order
bounds and their application to get convergence rates for static and dynamic
trivialisations are the main theoretical results of the thesis. We finish this
chapter with a short characterisation result for retractions that are
Riemannian exponential maps for some connection in terms of sprays on the
second tangent bundle $TTM$ (Theorem 5.7.1).

Chapter 6. Its main themes are practical optimisation on manifolds and deep
learning. In this chapter, we explain how to implement the ideas of static and
dynamic trivialisations that were laid out in Chapter 4. The exposition in
this chapter is done by developing the example of optimisation with orthogonal
constraints. We then present a literature review of the current state of the
art of practical optimisation on manifolds. After this, we go on to describe
the implementation details of library GeoTorch that we have developed to
perform scalable optimisation on manifolds in large models (Section 6.2.2). We
finish this chapter and the thesis showing a number of experiments to showcase
the practical efficiency of dynamic trivialisations in the context of
orthogonal rnns(Sections 6.3 and 6.4).

GeoTorch is available under mit license at

https://github.com/Lezcano/geotorch

The ideas in this chapter have been submitted to core PyTorch and have been
accepted to be part of PyTorch 1.9.0. The work in Section 6.3.2 led to the
addition of the matrix exponential into core PyTorch by Nikita Vedeneev. This
implementation yielded a $\times 12$ speed-up over the previous fastest GPU
implementation.

##  Chapter 2 Machine Learning and Neural Networks

In this chapter, we aim to contextualise the field of optimisation on
manifolds and its applications to neural networks within the broader field of
machine learning. As we will see, the field of theoretical machine learning
draws from fields as disparate as probability in Banach spaces, optimisation
in Hilbert spaces, computer science, and engineering. This makes it difficult
to find consistent terminology for this field, where some concepts have been
rediscovered several times during the last $50$ years. As such, this chapter
has an introductory as well as a unifying character, giving a short
introduction to all these ideas, their development throughout history and the
kind of results that one encounters in each relevant field, while also
introducing the nomenclature used in the fields of machine learning and deep
learning, which will be used throughout the rest of the thesis. This is done
in the hope that readers coming from other fields are not left with a text
full of foreign terminology without a proper dictionary.

All the theoretical ideas in this chapter and many more can be found in the
books vapnik1995naturenesterov2004introductory.

###### Outline of the chapter.

We will start by giving a short motivation for the ideas that underlay the
current theoretical approach to offline machine learning, with the usual split
between generalisation and optimisation. Then, we will provide a historical
review of some classical results in the field of statistical learning theory
and optimisation in $\mathbb{R}^{n}$, to serve as an introduction and
contextualisation to the convergence results that we will give in the thesis.
After that, we will give specific examples of a popular family of models
called neural networks, and we will show how to implement them in the PyTorch
library. Throughout the course of this introduction to the practical side of
machine learning, we will make sure to always draw connections to the theory
outlined before. We hope this connection between theory and practice helps to
put in context the ideas that will be introduced in the thesis, where
practical remarks will often be interleaved with more theoretical ones to give
algorithms that are solid on both ends.

###  2.1 Offline Statistical Learning

Before introducing what neural networks are, we shall first present the
standard framework for offline statistical learning in its most basic setting.

####  2.1.1 Framework and terminology

###### Features and labels.

Let $\mathcal{X}$ and $\mathcal{Y}$ be two sets. We say that $\mathcal{X}$ is
the feature space and $\mathcal{Y}$ is the label space. For example, if we
wanted to model whether the income of the parents affects whether their
children end up studying at Oxford, we would have
$\mathcal{X}\cong\mathbb{R}_{\geq 0}$ as the set of all possible incomes (for
example in pounds), and $\mathcal{Y}\cong\\{0,1\\}$ which models whether
certain children of theirs studied at Oxford or not.111The use of $\cong$
instead of $=$ here has a bit of a pedantic nature. Note that the set of all
weights $\mathcal{X}$ _is not_ intrinsically equal to the abstract set
$\mathbb{R}_{\geq 0}$ but it is merely identified mathematically with it.
Hence, $\mathcal{X}$ and $\mathbb{R}_{\geq 0}$ are isomorphic, and not
equal.222The model and methodology presented here is an over-simplification
which would not be of much use in the real world. In a real model, one would
compare many other variables, together with a much more complex analysis than
the one that we will present here. In general, we could have $\mathcal{X}$ and
$\mathcal{Y}$ to be any pair of sets, discrete or continuous, without any
structure.

###### Distribution and training set.

We assume that there exists a random variable $(X,Y)$ on
$\mathcal{X}\times\mathcal{Y}$ that we want to study. We do not have direct
access to $(X,Y)$, but we can sample from it. One convenient way of thinking
of the random variable $(X,Y)$ is as thinking that there exists certain random
variable $X$ on $\mathcal{X}$ which induces a conditional one $Y\mid X$ on
$\mathcal{Y}$.

In the context of offline statistical learning, we assume that we have access
to $m$ i.i.d. samples
$(X_{1},Y_{1}),\dots,(X_{m},Y_{m})\in\mathcal{X}\times\mathcal{Y}$ from
$(X,Y)$. We call these samples the training set, and we assume $m$ to be fixed
a priori. This is where offline learning differs from online learning. In
online learning, we assume that we have access to a sampling oracle that gives
us i.i.d. samples of $(X,Y)$ on demand, so our algorithm may adapt as we get
more and more samples.333Often these $m$ samples would be split into training
set, test set and validation set. Although this practice is fundamental in
practice, we will not talk about it here to simplify this introduction.

In the example of the acceptance at Oxford, the distribution $(X,Y)$ would be
the actual unknown distribution, having certain distribution on the income
$\mathcal{X}$ and certain probability of being accepted at Oxford given
certain fixed income $Y\mid X$. The training set could be $m$ samples that we
collect via a survey.

###### Model and loss function.

In order to model $Y\mid X$, we consider a family of functions
$f_{\theta}\colon\mathcal{X}\to\mathcal{Y}$ indexed by some parameter
$\theta\in\Theta$. Usually $\Theta$ will be a Euclidean space $\mathbb{R}^{n}$
for some $n>0$, but in more general settings, like those that we will consider
in this thesis, $\Theta$ can be a manifold or a more general topological
space. The family of functions $\\{f_{\theta}\>|\>\theta\in\Theta\\}$ is
called the model. The last thing that we will need, is a loss function
$\ell\colon\mathcal{Y}\times\mathcal{Y}\to[0,\infty)$, which serves as a proxy
for a distance on $\mathcal{Y}$, although it need not be symmetric—and often
it is not.

Following the previous example of admissions to Oxford, we would have to model
the binary output $\mathcal{Y}=\\{0,1\\}$. To do so, one usually considers the
relaxation $\mathcal{Y}\cong\\{0,1\\}\subseteq[0,1]$. An example of a model in
this case could be what is called polynomial linear regression of degree $n$
for a fixed $n>0$ together with a mapping $\sigma\colon\mathbb{R}\to[0,1]$.
This model can be mathematically described setting
$\Theta\cong\mathbb{R}^{n+1}$ and

$\\!\begin{aligned}
f_{\theta}\colon\mathcal{X}\times\mathbb{R}^{n+1}&\to[0,1]\\\
(x,\theta)&\mapsto
f_{\theta}(x)=\sigma\Big{\lparen}\sum_{k=0}^{n}\theta_{k}x^{k}\Big{\rparen}\end{aligned}$

and we could choose $\operatorname{sigm}(x)\coloneqq\frac{1}{1+\exp(-x)}$ as
the mapping. This is called the sigmoid function.444The image of the sigmoid
function is $(0,1)$ rather than the full $[0,1]$, which seems in contradiction
with the fact that we are interested in modelling a random variable with
values on the ends of $[0,1]$. This is done due to convenience, as $(0,1)$ is
diffeomorphic to $\mathbb{R}$ (through the sigmoid function, for example),
while $[0,1]$ is not even homeomorphic to $\mathbb{R}$.

In classification models, it is customary to think of the map $f_{\theta}$
with image contained in $[0,1]$ as a random variable on $[0,1]$ modelling
$Y|X$. As such, the usual loss function is the cross-entropy between the
distribution induced by $f_{\theta}$ and the empirical distribution

$\ell(x,y)=\begin{cases*}\log(x)&\text{if $y=1$}\\\ \log(1-x)&\text{if
$y=0$}\end{cases*}$

Note that, in this case, $y$ is assumed to be in $\\{0,1\\}$, while $x$ is
assumed to be in $(0,1)$, which is the image of $f_{\theta}$.

######  Remark 2.1.1 (Modelling a problem).

The choice of the model and the loss function heavily depends on the problem
at hand, and currently, it is closer to an art than a science. For this
reason, we will not explore these further theoretically, as the scope of doing
so is limited. We will come back to them in the experiment sections, when we
implement certain models to benchmark the algorithms that we will propose
against those in the literature.

###### Expected and empirical risk.

With these objects in hand, we can define the expected risk or population risk
as

$r(\theta):=\operatorname{\mathbb{E}}[\ell\lparen f_{\theta}(X),Y\rparen].$

We would then like to approximate a solution to the expected risk problem

$\operatorname*{arg\,min}_{\theta\in\Theta}\operatorname{\mathbb{E}}[\ell\lparen
f_{\theta}(X),Y\rparen].$ (2.1)

Of course, we cannot directly compute $r(\theta)$ as we do not have direct
access to the distribution of $(X,Y)$, but what we can do is to estimate it
via the empirical risk

$R_{m}(\theta):=\frac{1}{m}\sum_{i=1}^{m}\ell\lparen f_{\theta}\lparen
X_{i}\rparen,Y_{i}\rparen.$

We can compute the empirical risk, as opposed to the population risk. It also
has the desirable property of being an unbiased consistent estimator of
$r(\theta)$. Using this estimator, we can define the empirical risk
minimisation objective (erm)

$\min_{\theta\in\Theta}\frac{1}{m}\sum_{i=1}^{m}\ell\lparen f_{\theta}\lparen
X_{i}\rparen,Y_{i}\rparen.$ (2.2)

Denote an expected optimal value for the expected risk problem and an
empirical optimal value for the empirical risk problem by

$\theta^{\ast}\in\operatorname*{arg\,min}_{\theta\in\Theta}r(\theta)\qquad\psi^{\ast}\in\operatorname*{arg\,min}_{\theta\in\Theta}R_{m}(\theta).$

In general, we will not be able to compute $\psi^{\ast}$ exactly. We will use
some optimisation method that gives a sequence of estimates $\psi_{t}$ given
by some optimisation algorithm. Putting all these quantities together, we have
the so-called statistics-optimisation trade-off.

######  Proposition 2.1.2 (Statistics-optimisation trade-off).

We have the following factorisation:

$r(\psi_{t})-r(\theta^{\ast})\leq\underbrace{\sup_{\theta\in\Theta}\lparen
r(\theta)-R_{m}(\theta)\rparen+\sup_{\theta\in\Theta}\lparen
R_{m}(\theta)-r(\theta)\rparen}_{\text{Statistical
term}}+\underbrace{R_{m}(\psi_{t})-R_{m}(\psi^{\ast})}_{\text{Optimisation
term}}.$ (2.3)

###### Proof.

Interpolating we have

$r(\psi_{t})-r(\theta^{\ast})=\lparen
r(\psi_{t})-R_{m}(\psi_{t})\rparen+\lparen
R_{m}(\psi_{t})-R_{m}(\psi^{\ast})\rparen+\lparen
R_{m}(\psi^{\ast})-R_{m}(\theta^{\ast})\rparen+\lparen
R_{m}(\theta^{\ast})-r(\theta^{\ast})\rparen.$

Noting that the term in the third parenthesis is non-positive and taking
suprema we get the result. ∎

The area of statistical learning theory studies how to give bounds for the
statistical term for different classes of models,555More generally,
statistical learning theory tries to bound $r(\psi)$ for an arbitrary point
$\psi\in\Theta$. The point does not need to come from erm, and the bounds may
be different to the ones presented here. That being said, this is by far the
most used approach. and tries to give bounds on the number of samples $m$
needed to make the statistical term smaller than certain given quantity
$\varepsilon$. Given that the samples $(X_{i},Y_{i})$ are stochastic, the
results in this field are given as results with high probability, that is, one
proves that the statistical term is less than $\varepsilon$ with probability
at least $1-\delta$, where $\varepsilon,\delta>0$ are some parameters given a
priori. The abstract idea behind these type of results is that of
concentration of measure, quantifying for example the rate of convergence of
certain central limit theorem. These usually come from fine high-dimensional
versions of Markov’s and Chebyshev’s inequalities and other tail inequalities.

On the other hand, the area of stochastic optimisation studies how to design
algorithms to approximate solutions to the problem in (2.2). Given an
algorithm, one proves the convergence of $\psi_{t}$ as a function of the
number of steps $t$ to a minimiser $\psi^{\ast}$. One often assumes some
notion of convexity or some degree of regularity of the derivatives of the
function being minimised—$R_{m}(\theta)$ in this case—to be able to prove this
kind of results.

####  2.1.2 A note on statistical learning theory

The statistical term is studied in the statistical learning theory, often also
referred to as generalisation theory, given that it studies how well an
empirical optimum—for example $R_{m}(\psi^{\ast})$ or any other
estimator—approximates $r(\theta^{\ast})$. In this approach, the quantity
$R_{m}(\theta)$ is treated as random variable that depends on
$(X_{1},Y_{1}),\dots,(X_{m},Y_{m})$.

###### Dimensionless bounds.

In practice, machine learning models are larger and larger. In particular, the
parameter space $\Theta$ can have a dimension in the order of the billions. As
such, the bounds that this theory hopes to get are dimensionless, that is,
independent of the dimension of $\Theta$.666Note that this need not hold on
the optimisation term. When computing bounds for the optimisation term later
on, we will often put them in terms of the Lipschitz constant. This Lipschitz
constant often depends on the dimension of the problem. For this reason, this
theory is based in the field of high dimensional probability. When one is able
to get dimensionless bounds, it tends to be the case that the proofs and
results can be adapted to work in Hilbert and Banach spaces, and as such, it
should be no surprise that this theory is itself based in the theory of
probability in Banach spaces. The reference book in this field is Ledoux and
Talagrand’s [110].

###### Modelling the generalisation term.

The approach to get bounds on the statistical term in (2.3) comes from
assuming a metric space $(\Theta,d)$ on the parameter space and, for a fixed
$m>0$, considering the stochastic process $X_{\theta}=r(\theta)-R_{m}(\theta)$
that depends on the random variables
$(X,Y),(X_{1},Y_{1}),\dots,(X_{m},Y_{m})$. The problem can be then translated
to bounding above and below the quantity

$\operatorname{\mathbb{E}}\sup_{\theta\in\Theta}X_{\theta}.$

Generalising the problem by forgetting the particular structure of
$X_{\theta}$, one looks to give upper and lower bounds on the expectation of a
general stochastic process $X_{\theta}$ on a metric space $(\Theta,d)$.

The first result in this direction was the chaining argument. It was first
developed by Kolmogorov in $1934$ but never published [40], and it was then
generalised by Dudley [50], who proved Dudley’s entropy integral under the
assumption that $X_{\theta}$ is subgaussian, giving the upper bound

$\operatorname{\mathbb{E}}\sup_{\theta\in\Theta}X_{\theta}\leq
12\int_{0}^{\infty}\sqrt{\log
N\lparen\Theta,d,\varepsilon\rparen}\,\mathrm{d}\varepsilon$

where $N(\Theta,d,\varepsilon)$ is the $\varepsilon$-covering number of the
metric space $(\Theta,d)$, that is, the minimum number of balls of radius
$\varepsilon$ needed to cover the whole $\Theta$.

Lower bounds on this quantity were given in the theory of Gaussian processes
by Slepian, Fernike, and Sudakov. These results assume $X_{\theta}$ to be a
Gaussian process and consider the natural metric induced by it in $\Theta$,
which is the $L^{2}$ metric

$d(\theta_{1},\theta_{2})=\lVert
X_{\theta_{1}}-X_{\theta_{2}}\rVert_{2}=\operatorname{\mathbb{E}}[\lvert
X_{\theta_{1}}-X_{\theta_{2}}\rvert^{2}]^{1/2}.$

Using this distance, a Gaussian process is, by definition, also subgaussian
under this metric. These natural metrics allow giving not only lower bounds,
but actually matching upper and lower bounds. This was done by Talagrand in
the celebrated majorising measure theorem
talagrand1987regularitytalagrand1992simpletalagrand1994constructions. This
result is quite technical, as it involves some fairly complex constructions on
metric spaces. Intuitively, these structures can be seen as a refinement of
the global quantity $N(\Theta,d,\varepsilon)$ to be rather a finer local one.

Going back to statistical learning theory, one of the most basic models of
learning to consider is the set

$\Theta=\\{(-\infty,x]\>|\>x\in\mathbb{R}\\}$

and $f_{\theta}$ to be the indicator function on the interval
$\theta\in\Theta$. A natural metric for these functions is the $L^{\infty}$
metric, as it is a metric that makes the empirical measure
$\mu_{m}=\sum_{i=1}^{m}\delta_{X_{i}}$ subgaussian. On the other hand, it is
direct to see that the covering numbers
$N(\Theta,\lVert-\rVert_{\infty},\varepsilon)=\infty$ for any $\varepsilon$.
As such, all the previous bounds become trivial for the case of infinite sets
of indicator functions. However, it is possible to give finite bounds for
these functions via a process called symmetrisation, which can be summarised
in the following chain of inequalities

$\operatorname{\mathbb{E}}\bigg{[}\sup_{\theta\in\Theta}\sum_{i=1}^{m}\lparen
f_{\theta}(X_{i})-\operatorname{\mathbb{E}}f_{\theta}\rparen\bigg{]}\leq
2\operatorname{\mathbb{E}}\bigg{[}\sup_{\theta\in\Theta}\sum_{i=1}^{m}\varepsilon_{k}f(X_{i})\bigg{]}\leq\sqrt{2\pi}\operatorname{\mathbb{E}}\bigg{[}\sup_{\theta\in\Theta}\sum_{i=1}^{m}g_{k}f(X_{i})\bigg{]}$

where $\varepsilon_{k}$ are i.i.d. symmetric Bernoulli and $g_{k}$ are i.i.d.
$\mathcal{N}(0,1)$. The quantities in the first and second bound are called
the Rademacher complexity of the set of functions $\Theta$ and the Gaussian
complexity of $\Theta$ respectively. These inequalities tell us that bounding
any of these two complexities is enough to give bounds on the empirical
process, and as such, to give a priori bounds on the generalisation of the
class of our class of models. Most of the theory in statistical machine
learning centres around giving bounds for these two complexity quantities
under different assumptions on the family $f_{\theta}$.

For a rather in-depth introduction to these and other topics from a
probability theory perspective, we strongly recommend the book by Ramon van
Handel [74]. For a book with more applications towards the machine learning
see [173]. There is also the classical book by Vapnik which lays the
foundations of the field of statistical learning [171].

To finish this note on the statistical term, we shall give an example of the
magnitude of this term for a simple family of classifiers.

######  Proposition 2.1.3.

Let $(X,Y)$ be a random variable on
$\mathcal{X}\times\mathcal{Y}\subseteq\mathbb{R}^{n}\times\mathbb{R}$. Let
$B_{n}(r)\subseteq\mathbb{R}^{n}$ be the ball centred at zero with respect to
the $\ell^{2}$ norm. Consider two bounded sets $\Theta\subseteq
B_{n}(r_{\Theta})$, $\mathcal{X}\subseteq B_{n}(r_{\mathcal{X}})$ for two
fixed $r_{\Theta},r_{\mathcal{X}}>0$. Consider the family of linear predictors
$f_{\theta}(x)=\langle\theta,x\rangle$ parametrised by $\Theta$. For a loss
function $L$-Lipschitz on the first variable, we can bound the statistical
term as

$\operatorname{\mathbb{E}}\sup_{\theta\in\Theta}\lparen
R_{m}(\theta)-r(\theta)\rparen=\operatorname{\mathbb{E}}\bigg{[}\sup_{\theta\in\Theta}\sum_{i=1}^{m}\ell\lparen\langle\theta,X_{i}\rangle,Y_{i}\rparen-\operatorname{\mathbb{E}}\ell\lparen\langle\theta,X\rangle,Y\rparen\bigg{]}\leq
2\frac{r_{\mathcal{X}}r_{\Theta}L}{\sqrt{m}}.$

###### Proof.

See for example [159, Lemma 26.9 and Lemma 26.10]. ∎

If one further assumes subgaussianity of the distribution of $(X,Y)$, it is
then possible to convert these bounds in expectation to bounds with high
probability, that is, one gets tail bounds on the distribution of
$R_{m}(\theta)$.

######  Remark 2.1.4 (Qualitative properties of the bound).

The first thing to note is that, as we mentioned before, this bound is
independent of the dimension $n$ of the parameter space $\Theta$ and the input
space $\mathcal{X}$.

Secondly, this bound tells us that, in the most basic case, the statistical
term in (2.3) scales as $\mathcal{O}\lparen\frac{1}{\sqrt{m}}\rparen$. It is
also possible to show that this bound is asymptotically tight for this family
of classifiers. This bound suggests that it is enough to optimise up to a
precision on the order of $\mathcal{O}\lparen\frac{1}{\sqrt{m}}\rparen$.

The asymptotics $\mathcal{O}\lparen\frac{1}{\sqrt{m}}\rparen$ do not appear
just in the case of linear classifiers, but also in many other architectures.
In essence, they stem from the convergence rates given by the central limit
theorem. These rates are asymptotically tight unless extra assumptions are
considered. In typical scenarios, under more restrictive hypotheses, the rate
of convergence can be improved to $\mathcal{O}\lparen\frac{1}{m}\rparen$,
which is usually referred to as fast rate.

The take home from these ideas is that it should be just necessary to bound
the optimisation term at a rate $\mathcal{O}\lparen\frac{1}{\sqrt{m}}\rparen$,
as this bound is already imposed by the statistical term in the statistics-
optimisation factorisation.

###  2.2 First Order Optimisation

####  2.2.1 The complexity model

In the rest of the thesis, we will just look at the optimisation term
introduced in the previous section. More generally, we will be interested in
studying algorithms to approximate a solution of problems of the form

$\min_{x\in\mathsf{X}}f(x)$

where $\mathsf{X}$ is a subset of $\mathbb{R}^{n}$.777The contributions in
this thesis will be in the more abstract setting where $\mathsf{X}$ is a
differentiable manifold not necessarily embedded in $\mathbb{R}^{n}$. We will
restrict the presentation in this introductory section to subsets of
$\mathbb{R}^{n}$ for the sake of simplicity. It will be in Chapter 4 where we
will give an introduction to the subject of optimisation on manifolds together
with a literature review of the field in Section 5.1.2, once we have given a
presentation of the geometric tools and notation in Chapter 3. Even though we
will not mention the manifold setting in this section, we will write most
definitions in a coordinate-independent fashion, so that they readily
generalise to the manifold setting without much effort. Note that the set
$\mathsf{X}$ is what we called $\Theta$ in the context of modelling and
machine learning.

As it is customary when developing a complexity theory for a class of
algorithms, to be able even talk about what an algorithm is without getting
lost in computational semantic considerations, it is necessary to introduce a
computational model. A computational model in the theory of optimisation
abstracts the idea of _what is a function_ , as a function in computer science
can be defined in many possible ways and a definition of a function may be as
concrete as to be language-dependant.

In optimisation, the computational model abstracts a function $f$ in terms of
its value and derivatives at its points.

######  Definition 2.2.1.

A $k$-th order oracle $\mathscr{O}$ for a class $C^{k}(\mathsf{X})$ function
$f\colon\mathsf{X}\to\mathbb{R}$ is a mapping that, given a point
$x\in\mathsf{X}$, returns the $k$-th order information of $f$:

$\mathscr{O}(x)=(f(x),\nabla f(x),\operatorname{Hess}f(x),\dots,f^{k)}(x)).$

An algorithm in this model may query the oracle for the value of the function
and its derivatives at one or more points. Then, the algorithm will process
this information together with the information from previous calls to output a
point $x_{t}\in\mathsf{X}$. This new point will be used to generate more
queries to the oracle and so on. We note that, in general, there is no reason
for which an algorithm should perform just one call to the oracle in every
step. In fact, some very important algorithms require of two calls or even a
logarithmic number of calls to the oracle in every step. On the other hand, in
this thesis we will mostly consider algorithms that perform one call to the
oracle in every step.

In this thesis, we will centre our attention on first order methods. These are
methods that just have access to the value of the function and its gradient.
This is the most common class of methods used in the context of neural
networks since if a model has $n$ parameters—$\dim\mathsf{X}=n$—the Hessian
would have $n^{2}$ entries, which is prohibitively large even for small
models. It is worth mentioning that, when implementing second order
optimisation, one almost never computes the full Hessian. Instead, algorithms
are specified in terms of _Hessian-vector products_ of the form $H^{-1}(v)$
where $v$ is a vector and $H$ is the Hessian of $f$—second order methods—or an
approximation of it from first order information—quasi-Newton methods. Even
then, computing Hessian-vector products tends to be quite costly when compared
with first order methods both in time and memory.

Due to memory limitations given the large models used in deep learning, we
will just look at first-order algorithms that do not depend on previous calls
to the oracle. An algorithm $\mathcal{A}$ of this form is just a function

$x_{t+1}=\mathcal{A}(x_{t},\nabla f(x_{t})).$

These algorithms are inherently local, as they just process the information at
the current point to make their next guess.

######  Remark 2.2.2 (Optimisation algorithms as discretisations of odes).

A more formal definition of an algorithm being local comes by interpreting
such an algorithm as a discretisation of a differential equation. For example,
the gradient step is the forward-Euler discretisation of the differential
equation

$\dot{x}=\nabla f(x).$

This is an autonomous differential equation which gives rise to a flow called
the gradient flow. Since this is an invariant equation, this view can be
generalised to Riemannian manifolds provided a suitable discretisation scheme.
This more general definition of local algorithms allows for a constant number
of calls to the oracle to be used in each step, rather than just one. The
perspective of looking at the algorithm $\mathcal{A}$ as an integrator has
yielded a unified way of looking at large families of optimisation algorithms
under a common lens
wibisono2016variationalwilson2018lyapunovfranca2020dissipative.

######  Remark 2.2.3 (A note on stochastic optimisation).

As we mentioned in the previous section, in machine learning, the objective
function often depends on an unknown random variable $\xi$ so that the
objective function takes the form888The random variable $\xi$ represents in
the offline setting the product random variable of the dataset
$(X_{1},Y_{1}),\dots,(X_{m},Y_{m})$.

$f(x)=\operatorname{\mathbb{E}}[f(x,\xi)].$

In the offline learning setting, the random variable $\xi$ has finite support
and represents the dataset so that $f(x)=R_{m}(x)$.

The field that studies functions with this structure is called stochastic
optimisation, and it was introduced in the seminal paper by Robbins and Monroe
[153]. Stochastic optimisation has become particularly popular in the last
decade due to its applications in machine learning.

The simplest stochastic algorithm approximates the gradient of
$\operatorname{\mathbb{E}}[f(x,\xi)]$ using $B$ i.i.d. samples
$\xi_{i}\sim\xi$ to form a stochastic approximation to the direction of
steepest descent of $f$

$x_{t+1}=x_{t}-\frac{\eta}{B}\sum_{i=1}^{B}\nabla f(x_{t},\xi_{i}).$ (2.4)

This is called stochastic gradient descent (sgd). The samples $\xi_{i}$ are
called the minibatch—or sometimes simply the batch—and the number $B$ is
called the minibatch size. Many other algorithms have been developed in recent
years for the stochastic optimisation in the online and offline setting,
yielding algorithms widely used in practice such as adagrad [49], and its
heuristic variations rmsprop [82] and adam [99]. In parallel, the field of
stochastic optimisation has grown in its own directions, for example by using
this framework to set online decision problems [29].

In the rest of the thesis, we will develop the theory in the non-stochastic
setting. On the other hand, the results that we will present are quite
general, so it should not be too difficult to generalise the ideas to the
stochastic setting. That being said, we will see in the experiments section
that the theory developed here works surprisingly well in practice in the
stochastic setting. We leave the theoretical development of these ideas open
for future research.

Having introduced the abstract computational model, we will spend the rest of
the section introducing and giving examples of the two main classes of
functions considered in the optimisation literature.

####  2.2.2 Convex optimisation

Convex optimisation studies the setting in which the function $f$ is convex
and the set of constraints $\mathsf{X}$ is a convex subset of
$\mathbb{R}^{n}$. This is the most classical scenario as there are many
interesting problems in practice that are either convex or can be relaxed into
a convex problem, while still reflecting most of the initial structure of the
problem.

The convexity hypothesis is a quite restrictive one, but at the same time, it
allows applying local algorithms to solve global problems over the whole space
$\mathsf{X}$. This comes from the fact that for a convex function on a convex
space, every local minimiser is a global minimiser. Furthermore, if the
function is strictly convex, there is a unique global minimiser. Convexity
(resp. strict convexity) is a local property when the function is of class
$C^{2}(\mathsf{X})$—it is equivalent to $\operatorname{Hess}f(x)\succeq 0$
(resp. $\operatorname{Hess}f(x)\succ 0$) at every point $x\in\mathsf{X}$—so it
makes tackling the global problem of optimisation by local methods feasible.

It would be rather naïve to try to give a comprehensive introduction to convex
optimisation in a couple of pages given the vast literature on the subject.
For example, when the function is linear or quadratic, one enters the fields
of linear programming, conic programming and semidefinite programming among
others, which study algorithms such as interior point methods and simplex
methods, each of which decomposes in many other subfamilies. To avoid falling
down that rabbit hole, we will just mention methods that can be applied in
practice to general non-linear functions, even if we are just able to prove
convergence for a subfamily of convex functions.

Furthermore, to simplify the exposition, we will assume that $f$ is at least
of class $C^{2}(\mathsf{X})$. Virtually all the theory in convex optimisation
can be reformulated without this hypothesis via the use of subgradients. On
the other hand, the $C^{2}(\mathsf{X})$ regularity assumption often clarifies
the exposition, and proofs in the $C^{2}(\mathsf{X})$ setting can often be
generalised to the convex non-differentiable setting via standard algebraic
arguments.

To exemplify the kind of results that one gets in this field, let us consider
the simplest of the algorithms: Gradient descent. Gradient descent
(gd)—sometimes called steepest descent—accounts for following the direction of
steepest descent for a time $\eta>0$. Starting at a point
$x_{0}\in\mathsf{X}$, the update rule is given by

$x_{t+1}=x_{t}-\eta\nabla f(x_{t}).$

The quantity $\eta$ is usually called the step-size, although in machine
learning is often called the learning-rate. Now, this begs the question: What
is a reasonable step-size?

In the same way that one does when giving bounds on the discretisation term of
an ode, to answer this question one has to put some assumptions on the
function. The simplest assumption is that of smooth convexity.

######  Definition 2.2.4.

A function $f\in C^{2}(\mathsf{X})$ is said to be of $\alpha$-bounded hessian
if

$-\alpha\mathrm{I}_{n}\preceq\operatorname{Hess}f(x)\preceq\alpha\mathrm{I}_{n}\mathrlap{\qquad\forall
x\in\mathsf{X}.}$

If the function $f$ is also convex (i.e., $0\preceq\operatorname{Hess}f(x)$),
$f$ is said to be $\alpha$-smooth.

In other words, a $C^{2}$ function is of $\alpha$-bounded Hessian if the
eigenvalues of its Hessian lie in the interval $[-\alpha,\alpha]$, and it is
$\alpha$-smooth if they lie in $[0,\alpha]$.

Expanding the function as a Taylor series along the straight line that
connects two points $x,y\in\mathsf{X}$ and bounding Lagrange’s remainder we
get the following bound around $x\in\mathsf{X}$ for a function of
$\alpha$-bounded Hessian

$f(y)\leq f(x)+\langle\nabla f(x),y-x\rangle+\frac{\alpha}{2}\lVert
x-y\rVert^{2}.$ (2.5)

Denote $\gamma(t)=x+t(y-x)$ for $t\in[0,1]$ the line connecting $x$ and $y$.
This bound can be interpreted geometrically as having an upper bound on
$f|_{\gamma}$ by a quadratic function tangent to $f$ at $x$ with leading term
$\frac{\alpha}{2}$. Note that this bound is everywhere sharp for the distance
function to a fixed point $f(x)=\frac{\alpha}{2}d\lparen x,x_{0}\rparen^{2}$
as this is a quadratic function with Hessian equal to
$\alpha\cdot\mathrm{Id}$. For this function, in the unconstrained case
$\mathsf{X}=\mathbb{R}^{n}$, it is possible to get to the minimum $x_{0}$ from
any point in one gradient descent step letting $\eta=\frac{1}{\alpha}$.
Algebraically, we can see all these ideas by setting
$y=x_{t+1}=x_{t}-\frac{1}{\alpha}\nabla f(x_{t})$ in (2.5) to get

$f(x_{t})-f(x_{t+1})\geq\frac{1}{2\alpha}\lVert\nabla f(x_{t})\rVert^{2}.$
(2.6)

Geometrically, this means that we improve the value of the function at least
by $\frac{1}{2\alpha}\lVert\nabla f(x_{t})\rVert^{2}$ after every step. On the
other hand, by convexity and Cauchy-Schwartz, one may bound below the norm of
the gradient by

$\frac{f(x_{t})-f(x^{\ast})}{\lVert x_{t}-x^{\ast}\rVert}\leq\lVert\nabla
f(x_{t})\rVert.$

Plugging this into (2.6), putting the resulting equation in terms of
$\delta_{t}=f(x_{t})-f(x^{\ast})$ and showing that $\lVert
x_{t}-x^{\ast}\rVert$ is decreasing in $t$, one gets the inequality

$\delta_{t+1}\leq\delta_{t}-\frac{1}{2\alpha\lVert
x_{0}-x^{\ast}\rVert^{2}}\delta^{2}_{t}.$

and bounding the general term of this recurrence, the following full
convergence result follows.

######  Theorem 2.2.5.

Assume that $\mathsf{X}=\mathbb{R}^{n}$ and $f$ is $\alpha$-smooth and has a
minimiser $x^{\ast}$. Then, gradient descent with step $\eta=\frac{1}{\alpha}$
starting at a point $x_{0}$ converges as

$f(x_{t})-f(x^{\ast})\leq\frac{2\alpha\lVert x_{0}-x^{\ast}\rVert^{2}}{t+4}.$

This result says that the optimisation term in the setting of convex
$\alpha$-smooth functions decays as $\mathcal{O}\lparen\frac{1}{t}\rparen$. A
much more impressive result is the one developed by Nesterov in $1983$.
Nesterov proposed an algorithm such that, under the same $\alpha$-smoothness
hypothesis, converged at a rate of $\mathcal{O}\lparen\frac{1}{t^{2}}\rparen$.
This was a major result, as a few years before it had been proven that this
rate was optimal for the set of convex $\alpha$-smooth functions. After that
result by Nesterov, many other algorithms that converge at this fast rate have
been developed. In general, an algorithm that converges faster than vanilla
gradient descent for a class of functions is called an accelerated algorithm.

This type of results gives us a worst-case analysis of the optimisation
oracle, as one may reword this as “gradient descent converges to a point with
a value at most $\varepsilon$ away of the minimum in at most
$\mathcal{O}(\frac{1}{\varepsilon})$ steps”.

Of course one can easily generalise this algorithm to one where
$\mathsf{X}\subseteq\mathbb{R}^{n}$ is a proper convex subset, provided that
we have an efficient way to project from $\mathbb{R}^{n}$ onto $\mathsf{X}$.
This is called projected gradient descent. In this case, a similar convergence
analysis gives the same rates as in the unconstrained case.

There are many other classes that are studied in the setting of convex
optimisation, such as $\alpha$-strongly convex functions
($\operatorname{Hess}f\succeq\alpha\mathrm{I}_{n}$ for $\alpha>0$), Lipschitz
functions and combinations of these. Another line of research is that of
studying convergence under weaker regularity assumptions than
differentiability. In this case, one of the most commonly known algorithm
given its versatility is Nemirovsky and Yudin’s mirror descent [137], which
may be applied in the more general setting of Banach spaces. In general, one
can find a vast zoo of algorithms with convergence guarantees in convex
optimisation. To name a few, by replacing the penalty function in mirror
descent by a non-linear function one gets proximal methods. If one can do
linear optimisation over the domain, one may use Frank-Wolfe. In coordinate
gradient descent and block-coordinate gradient descent one optimises one or
more of the coordinates of the function at a time, leaving the others fixed
during that iteration. In alternated direction method of multipliers (admm)
for constrained optimisation one tries to approximate a solution to equation
given by the Lagrange multipliers iteratively updating the multipliers… For a
theoretical introduction to all these methods and others see
nesterov2004introductoryboyd2004convex.

####  2.2.3 Non-convex optimisation

In the setting of non-convex optimisation, one drops completely the assumption
of $f$ being convex. On the one hand, this gives a much wider range of
applicability to the results developed by this theory. On the other, the
results that this theory is able to prove are much weaker. In particular, in
this theory, as the global properties given by convexity are dropped, one is
likely to encounter critical points and local minima of $f$, to the point
that, in general, not even finding but checking that a given point is a local
minimum for a non-convex function is NP-hard [133]. For this reason, in this
setting, one often has to content oneself with approximating critical points
of $f$, that is, points at which $\nabla f(x)=0$.

Of course, as it was the case with convex functions, if one wants to get
reasonable bounds, it is necessary to restrict the attention to a subset of
functions. Again, the class of functions of $\alpha$-bounded Hessian is
suitable for this task. For these functions, we get the following standard
result. We prove it here for completeness, as it will be the result that we
will generalise to manifold in Chapter 5 (cf., Theorem 5.6.4).

######  Theorem 2.2.6.

Let $f\colon\mathbb{R}^{n}\to\mathbb{R}$ be a function of $\alpha$-bounded
Hessian bounded below by $f^{\ast}\in\mathbb{R}$. Gradient descent on $f$,
with a step-size of $\frac{1}{\alpha}$ starting at a point $x_{0}$, converges
as

$\min_{0\leq k\leq t}\lVert\nabla f(x_{k})\rVert\leq\sqrt{\frac{2\alpha\lparen
f(x_{0})-f^{\ast}\rparen}{\lparen t+1\rparen}}.$

###### Proof.

By (2.6), we have that for a function of $\alpha$-bounded Hessian

$\frac{1}{2\alpha}\lVert\nabla f(x_{k})\rVert^{2}\leq f(x_{k})-f(x_{k+1}).$

Summing over all the indices and telescoping

$\frac{1}{2\alpha}\sum_{k=0}^{t}\lVert f(x_{k})\rVert^{2}\leq
f(x_{0})-f(x_{t+1})\leq f(x_{0})-f^{\ast}$

which yields the desired bound. ∎

In other words, if the function is of $\alpha$-bounded Hessian, gradient
descent with constant step-size finds a point with gradient
$\varepsilon$-close in norm to a critical point of $f$ in at most
$\mathcal{O}(\frac{1}{\varepsilon^{2}})$ steps, as opposed to the
$\mathcal{O}(\frac{1}{\varepsilon})$ needed to approximate the global
minimiser when the function is convex.

######  Remark 2.2.7 (Projected gradient descent on non-convex functions).

It is worth noting that the previous theorem assumes the problem to be
unconstrained, that is, $\mathsf{X}=\mathbb{R}^{n}$. One could hope to extend
this proof to work on constrained optimisation problems, as in the convex
case. A moment’s reflection shows that this is not easily achieved. In fact,
non-convex constrained optimisation is notoriously difficult, as the
projection does not respect the derivatives in any reasonable way, which
destroys the convergence of the algorithm. As such, these two elements do not
interact well together.

Even though the proof for gradient descent for non-convex objectives assures
that gradient descent will converge to a critical point, this critical point
might not be a local minimum, but perhaps a saddle point or a local maximum.
On the other hand, via an application of the stable manifold theorem from
dynamical systems, it is possible to prove that if the initial point $x_{0}$
is sampled from a smooth distribution on $\mathbb{R}^{n}$, gradient descent
will almost always converge in the limit to a local minimiser [111], although
it may take exponentially many steps in $\varepsilon$ [48]. These results
inspired many others that tried to estimate Hessian information to escape
saddle-points by perturbing the gradient [5].

In the last $5$ years, there has been renewed interest in the field of non-
convex optimisation due to the popularity of neural networks. Before that,
there had been very little progress since the years $1965$–$1975$ [138]. The
new lower-bounds developed for first order (and higher-order) methods on
several classes of non-convex functions—assuming certain first and second
order bounds—have been of special interest
cartis2010complexitycarmon2020lower. In parallel to this work, several
accelerated algorithms for non-convex optimisation [33] and algorithms in the
stochastic and non-stochastic setting for functions of the form
$f(x)=\frac{1}{n}\sum f_{i}(x)$, like the empirical risk. For a recent review
of these results see the papers allenzhu2018natasha2allenzhu2018neon2.

Other directions of research in this area are those coming from convex
relaxations, where one is able to find a convex problem that has the same
minimiser as the non-convex one, such as those coming from geometric
programming, many linear programming problems, and some eigenproblems such as
the svd decomposition. Many of these non-convex problems that happen to have a
fast way to solve them while being non-convex often have a hidden Riemannian
structure that makes them geodesically convex under some convenient metric. We
will explore in depth this kind of problems in Chapter 4. For now, we will
close this review section referencing what is, as the time of this writing,
the largest compilation to our knowledge of papers and reviews in the area of
provable non-convex optimisation

https://sunju.org/research/nonconvex/

The web-page is active as of $2020$, and contains lists of relevant papers
written in the last $5$–$10$ years in most of the sub-fields of non-convex
optimisation.

###  2.3 Neural Network Architectures

In this section, we introduce neural networks in the framework of statistical
models. We give some examples of different classes of neural networks, and we
show how to implement them in the Python library PyTorch. A good introductory
reference for all these ideas and many more is the book [67].

####  2.3.1 Neural networks

######  Definition 2.3.1.

Let $\Theta\subseteq\mathbb{R}^{n}$ be open and connected. A neural network
architecture is a mapping $f\colon\Theta\times\mathcal{X}\to\mathcal{Y}$ that
is differentiable in the first variable and that factorises into a composition
of simpler functions
$f^{(i)}\colon\Theta_{i}\times\mathbb{R}^{d_{i-1}}\to\mathbb{R}^{d_{i}}$ for
$i=1,\dots,p$ with $\mathcal{X}\cong\mathbb{R}^{d_{0}}$,
$\mathcal{Y}\subseteq\mathbb{R}^{d_{p}}$ and
$\Theta=\Theta_{1}\times\cdots\times\Theta_{p}$

$f_{\theta}=f_{\theta_{p}}^{(p)}\circ\dots\circ f_{\theta_{1}}^{(1)}.$

We will often omit the parameters of the functions and write
$f^{(i)}=f_{\theta_{i}}^{(i)}$. Each parametrised function $f^{(i)}$ is called
a layer. We say that $f^{(1)}$ is the input layer, $f^{(i)}$ for
$i=2,\dots,p-1$ are the hidden layers, and $f^{(p)}$ is the output
layer.999Sometimes the input and output layers are also called first and last
hidden layers. A network is said to be deep whenever it has three or more
layers. Deep learning is the area that studies deep neural networks.

Each of the layers in a neural network is designed to be a (parametrised)
transformation that can be computed very efficiently. In particular, most of
the layers used in practice are designed so that they can be heavily
parallelised, making their implementation very fast on a gpu. We will touch on
this in more detail later in this section.

The size of deep neural networks has seen an exponential explosion in size in
the last ten years, seeing models that have hundreds of layers and up to a
hundred billion parameters the current largest one [28].

######  Remark 2.3.2.

In most standard neural networks, $\Theta$ is just a Euclidean space. As we
did in the previous section, we will assume that $\Theta$ is isomorphic to a
Euclidean space in this section for simplicity, but the reader should bear in
mind that the direction of the rest of the thesis will be that of considering
$\Theta$ to be a connected manifold—often a matrix manifold.

####  2.3.2 Feedforward networks

The most classical neural network architecture is that of the feedforward
network. These are models where each layer is an affine transformation of its
input together with a coordinate-wise non-linear transformation.

######  Definition 2.3.3.

Let $\mathcal{X}\cong\mathbb{R}^{k}$, $\mathcal{Y}\cong\mathbb{R}^{d}$ and
$\Theta\cong\mathbb{R}^{d\times k}\times\mathbb{R}^{d}$. A linear layer of
width $d$ is the affine transformation defined by

$\\!\begin{aligned} f\colon\Theta\times\mathcal{X}&\to\mathcal{Y}\\\
((A,b),x)&\mapsto f_{A,b}(x)=Ax+b\end{aligned}$

We say that $A$ is the matrix of weights and $b$ is the bias.

Given a function $\sigma\colon\mathbb{R}\to\mathbb{R}$, we define a non-
linearity as the associated map $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$
that comes from applying the map $\sigma$ on vectors in $\mathbb{R}^{d}$
component-wise. We will denote both maps with the same symbol.

A feedforward network is a neural network where every layer is a linear layer
followed by a non-linearity, $f^{(i)}_{\theta_{i}}=\sigma\circ
f_{A_{i},b_{i}}$, $\theta_{i}=(A_{i},b_{i})$. Usually, the non-linearity in
the last layer is chosen to be the identity—i.e., there is no non-linearity in
the last layer.

###### Strengths, weaknesses, and limitations of feedforward networks

A feedforward network is one of the simplest neural network architectures used
in practice. It is sometimes also referred to as multi-layer perceptron (mlp).
The work of Cybenko [46] and Hornik [84] shows that feedforward networks are
universal approximators, in the sense that they are dense in the set of
continuous functions supported on a compact subset of $\mathbb{R}^{k}$. On the
other hand, an important remark is in order: While these functions may
approximate any continuous function arbitrarily well, the width of a shallow
network necessary to get an $\varepsilon$-approximator grows exponentially
fast in $\varepsilon$ [54]. It is this one of the reasons why feedforward
networks are often replaced in practice by more specialised layers that try to
encode domain-specific biases. Examples of these are cnns [62], attention
[15], neural machine translation [94] and transformer models [172] among many
others. More generally, the field of geometric deep learning studies how to
encode invariances present in the data so that the model structurally
preserves them, via the use of equivariant layers [43], graph neural networks
[158], and other architectures.

####  2.3.3 Recurrent neural networks (rnns)

The simplest non-trivial data topology is that of sequential data. This comes
up naturally when processing an input with variable length like audio samples,
time series, or sentences in a language. In these cases, every element of the
sequence is usually encoded into a vector of a fixed dimension, but the number
of these vectors is unknown a priori. Using computer science nomenclature,
these objects are encoded in the free monoid
$\lparen\mathbb{R}^{k}\rparen^{\ast}\coloneqq\oplus_{t=0}^{\infty}\lparen\mathbb{R}^{k}\rparen^{t}$.
Recurrent Neural Networks were designed to process this kind of sequential
data, constituting a map from $\lparen\mathbb{R}^{k}\rparen^{\ast}$ into a
finite-dimensional embedding space $\mathbb{R}^{d}$.

######  Definition 2.3.4.

Given a sequence of inputs $\lparen
x_{t}\rparen_{t=1}^{T}\in\lparen\mathbb{R}^{k}\rparen^{\ast}$, we define a
recurrent neural network (rnn) with hidden size $d>0$ as

$h_{t}=\sigma\lparen Bh_{t-1}+Cx_{t}\rparen\mathrlap{\qquad t=1,\dots,T}$

where $B\in\mathbb{R}^{d\times d}$ and $C\in\mathbb{R}^{d\times k}$ are
parameters, $\sigma$ is some fixed non-linearity, and $h_{0}\in\mathbb{R}^{d}$
is some fixed initial vector. $C$ is called the input kernel and $B$ is the
recurrent kernel. In most applications, the starting vector $h_{0}$ is set to
zero.

In classification tasks, where one needs an encoding of the whole sequence
into a vector, one often chooses either $h_{T}$, or $\frac{1}{T}\sum h_{t}$ as
the encoding of the whole sequence in $\mathbb{R}^{d}$. The problem with these
encodings is that they are biased towards the last elements of the sequence,
given that these were the elements that were encoded the last. For that
reason, one often uses a bi-directional rnn, which accounts for running two
separate rnns, one on $\lparen x_{t}\rparen_{t=1}^{T}$ and another one on the
reversed sequence $\lparen x_{T-t}\rparen_{t=1}^{T}$ and computes the
embedding as the sum of the hidden states of the two.

The intermediate vectors $h_{t}$ are usually regarded as the memory of the rnn
at time $t$. Each of the applications of the recurrent step is called a time-
step, and the dimension $t$ is called the time dimension.

rnns are a rather straightforward generalisation of feedforward neural
networks to sequences. On the other hand, they present some problems when the
matrix $B$ is not well-conditioned. Consider the degenerate case where
$\sigma=\mathrm{Id}$. If $B$ has eigenvalues of norm larger than $1$, the
gradient of $B$ will explode in some directions for long sequences, while if
$B$ has eigenvalues of norm less than $1$, the gradients of $B$ with respect
to the first time-steps will be close to zero, so the gradient will
incorporate almost no information of the initial elements of the sequence.
These problems are regarded as exploding gradient and vanishing gradient
problems [16]. We discuss them in more detail in Section 6.1.

The Long-Short Term Memory cell (lstm) [83] is a more complex recurrent
architecture that tries to heuristically alleviate the vanishing gradient and
exploding gradient problems via a system referred to as gating. The idea of
this heuristic is to _control the flow of the gradients along the lstm_, by
allowing to retain some information in the internal state, while giving
mechanisms to delete other information—often thought of as _forgetting_. The
gating mechanism is composed of $4$ vectors: input $i$, forget $f$, output,
and output modulation $g$. The lstm also separates the hidden state into an
internal state $c_{t}$—sometimes referred to as cell—and the hidden state
$h_{t}$. If we denote by $\operatorname{sigm}(t)=\frac{1}{1+e^{-t}}$ the
sigmoid function and by $\tanh$ the hyperbolic tangent—keeping in mind that
the range of these functions are $(0,1)$ and $(-1,1)$ respectively—the lstm is
defined by the equations

$\displaystyle\begin{pmatrix}i\\\ f\\\ o\\\
g\end{pmatrix}=\begin{pmatrix}\operatorname{sigm}\\\ \operatorname{sigm}\\\
\operatorname{sigm}\\\ \tanh\end{pmatrix}\Big{\lparen}W\begin{pmatrix}x_{t}\\\
h_{t-1}\end{pmatrix}\Big{\rparen}$ $\displaystyle
c_{t}=f\mathbin{\odot}c_{t-1}+i\mathbin{\odot}g\qquad
h_{t}=o\mathbin{\odot}\tanh(c_{t})$

where $\mathbin{\odot}$ denotes the component-wise multiplication and
$W\in\mathbb{R}^{4d\times(d+k)}$ is a matrix of parameters. There are other
popular variations of this recurrent layer, such as the Gated Recurrent Unit
(gru) [41], but we will spare the reader their description, as we will not use
them in the practical work in this thesis.

####  2.3.4 Why affine layers?

The layers that we have introduced until now are all variations of the same
idea: Applying an affine transformation to the inputs and then apply a
component-wise non-linearity. This is not a coincidence. This design comes
from engineering considerations arising from the implementation of neural
networks. Even though neural networks were first devised in the $80$’s [62],
they were partly forgotten until the first decade of the XXI century
ciresan2010deepkrizhevsky2012imagenet. Current interest has built upon two
large engineering efforts that made neural networks accessible to users
outside of academia. These two engineering advances are Graphics Process Units
(gpus) and auto-differentiation engines.

#####  2.3.4.1 Graphics Processing Units (gpus)

gpus have been around for many years now, mostly tied to the design and video
game industry. It has been just in the last $10$ years that we have seen the
field of gpu computing soar as a major tool in data analysis, allowing
processing large quantities of datapoints in parallel.

A coarse but fairly accurate description of gpus would come from thinking
about them as pieces of hardware that perform operations massively in parallel
much more quickly than a cpu. In general, they can apply the same operation in
parallel to many datapoints. While a cpu can execute many different operations
on individual pieces of data with a very low overhead between operations, a
gpu shines at executing the same operation in parallel on many datapoints. We
can think of a gpu as having a cpu with a very large number of cores, where
the cores are tied to execute the same operation.

If we look at a matrix-vector multiplication, for a matrix of size $d\times
k$, this can be decomposed into $d$ inner-products on $\mathbb{R}^{k}$, each
of which can be computed in parallel. Furthermore, each inner-product accounts
for adding $k$ different elements, which can be performed in $\log(d)$ time.

At the same time, when performing stochastic gradient descent as described in
(2.4), in every step of the algorithm we need to compute the value of the
function at $B$ datapoints, where $B$ is the minibatch size. This can also be
performed in parallel, taking again advantage of the fast parallel processing
given by gpus. This has allowed neural networks to be able to scale the
approximation of the minimisation problem to large datasets, having some of
them in the order of $10^{8}$ datapoints
deng2009imagenetkuznetsova2020openmahoney2006large.

###### A word on efficiency and gpus

gpus provide an efficient way to execute the same operation on many datapoints
in parallel. Each of these highly parallel operations is called a kernel.
Launching a kernel has a much larger overhead than executing an operation on a
cpu. For this reason, gpus shine when reducing the amount of kernels launched
during a program. We can think of executing a kernel as having a unitary fixed
cost—that of the overhead—when executed on very few points. For this reason,
computing a sequence of multiplications of sparse matrices can often be more
efficiently implemented on a cpu than on a gpu, due to the lower overhead per
operation of the cpu. gpus also suffer from a penalty coming from memory
latency, as each minibatch has to be moved from cpu to gpu, which can result
on a computational bottleneck in some applications. The take home from all
these ideas is that one has to be careful when computing the efficiency of an
algorithm implemented on a gpu, taking into account the size of the data one
is working with and how parallelisable are the operations that are being
executed. This latter point will be of particular importance when devising a
fast implementation of the methods proposed in this thesis in Section 6.3.2.

#####  2.3.4.2 Auto-differentiation: Forward and backward pass

In the offline learning setting, as we detailed in Section 2.1, we have access
to a dataset $\\{(X_{i},Y_{i})\\}$ of $m$ pairs of features and labels. As we
mentioned before, in deep learning, in many datasets we can have $m$ to be in
the order of the millions. Computing the gradient of the empirical risk

$R_{m}(\theta)=\frac{1}{m}\sum_{i=1}^{m}\ell\lparen f\lparen
X_{i},\theta\rparen,Y_{i}\rparen$

would be infeasible, as every step would require to evaluate the function $f$
on millions of points. As an alternative, stochastic gradient descent is used.
As we already mentioned, this algorithm chooses in every step a subset of $B$
elements $\\{(X_{a_{i}},Y_{a_{i}})\\}$ from the dataset to form a minibatch
and averages their losses to form

$\widehat{R}_{B}(\theta)\coloneqq\frac{1}{B}\sum_{i=1}^{B}\ell\lparen f\lparen
X_{a_{i}},\theta\rparen,Y_{a_{i}}\rparen$ (2.7)

performing the update

$\theta\leftarrow\theta-\eta\nabla\widehat{R}_{B}(\theta).$ (2.8)

In order to implement this algorithm in practice, it is necessary to be able
to compute the gradients of the function $\theta\mapsto\ell\lparen f\lparen
X_{i},\theta\rparen,Y_{i}\rparen$. It was in the seminal work [42] where it
was first documented that this naïve looking algorithm was able to perform
astonishingly well at optimising highly non-convex neural networks.
Furthermore, in this same work, they showed that this optimisation process
could be accelerated several orders of magnitude by implementing it on a gpu.
In this first work, the authors implement themselves the necessary gpu
kernels. This is a very difficult task, as programming on a gpu is notoriously
technical.

In parallel to this work, there had been a line of research that involved the
computation of the derivatives of functions automatically [176]. These are
called auto-differentiation engines or simply autodiff engines.

It was just then that Theano was born [152]. Theano was the first library
marrying auto-differentiation and gpus first released in 2007 as a library for
the programming language Lua. An autodiff engine helps the user define
differentiable functions: The user writes the function using methods provided
by the library and the library will automatically compute its derivative at a
given point. In some sense, it implements the chain rule—more formally, the
adjoint method—for a wide range of functions and constructions in a
programming language. In practice, this is a rather challenging problem.

####  2.3.5 An implementation of a feedforward network

Let us explain the rudiments of how all the ideas in this section come
together via an example written on the PyTorch library, as this will be the
library that we will use all throughout Chapter 6.101010A disclaimer to the
seasoned practitioner: This is intended to be a simple working example to
showcase the ideas in this chapter. It does not try to be a performant one by
any possible measure.

######  Example 2.3.5 (mnist classification with a feedforward network).

Consider a two layer feedforward neural network with
$\texttt{relu}(x):=\max(x,0)$ non-linearities designed to classify greyscale
images of handwritten digits. A dataset that provides this exact set-up is
that of mnist [109]. The mnist dataset consists of $60.000$ labelled images of
size $28\times 28$ represented by a floating-point number at every pixel
describing the intensity of the grey (see Figure 2.1). As such,
$\mathcal{X}=\mathbb{R}^{28\times 28}$, $\mathcal{Y}=\\{0,1,\dots,9\\}$.

Figure 2.1: Example of digits present in the mnist dataset.

The model will output a distribution on the set of $10$ elements by processing
the $10$-dimensional output of the last layer with a $\operatorname{softmax}$
non-linearity, which maps a $10$ dimensional vector into the $10$-dimensional
simplex by mapping its $k$-th coordinate to

$\operatorname{softmax}\lparen
x\rparen_{k}=\frac{e^{x_{k}}}{\sum_{i=0}^{9}e^{x_{i}}}.$

In practice we will output $\log\circ\operatorname{softmax}$ for numerical
stability reasons.

###### Loss

The loss function will be the negative log-likelihood at a given label

$\\!\begin{aligned}
\ell\colon\mathbb{R}^{10}\times\\{0,1,\dots,9\\}&\to\mathbb{R}\\\ \lparen
x,k\rparen&\mapsto-\log\lparen\operatorname{softmax}\lparen
x\rparen_{k}\rparen\end{aligned}$

###### Hyperparameters

We perform stochastic gradient descent with $\eta=10^{-3}$ and a minibatch
size of $128$. We will iterate over the whole dataset $80$ times. The hidden
size of the two layer network will be $256$.

###### Data processing

The dataset is normalised to have mean $0$ and standard deviation of $1$. This
is a common preprocessing step. The data will also be first flattened by the
model, that is, it will be converted from data in $\mathbb{R}^{28\times 28}$
into data in $\mathbb{R}^{784}$.111111By applying this transformation, we are
effectively _forgetting_ the $2$D structure of the data. This is far from
being an optimal way to process an image, but it is the simplest.

The procedure to fit this model to the dataset can be written in PyTorch 1.6
as shown in LABEL:lst:feedforward.

⬇

1import torch

2import torch.nn as nn

3import torch.nn.functional as F

4from torchvision import datasets, transforms

5

6# Hyperparameters

7epochs = 80

8batch_size = 128

9learning_rate = 0.001

10

11# Data processing

12# Centre the MNIST dataset, which has mean 0.1307 and standard deviation
0.3081

13transform = transforms.Compose([

14 transforms.ToTensor(),

15 transforms.Normalize((0.1307,), (0.3081,))

16 ])

17dataset = datasets.MNIST(’./data’, train=True, download=True,
transform=transform)

18loader = torch.utils.data.DataLoader(dataset, batch_size=batch_size,
shuffle=True)

19

20# Run the model on GPU if possible

21device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

22

23# Instantiate model and optimiser (Stochastic Gradient Descent)

24model = nn.Sequential(

25 nn.Flatten(),

26 nn.Linear(28*28, 256),

27 nn.ReLU(),

28 nn.Linear(256, 10),

29 nn.LogSoftmax(dim=1)

30 ).to(device)

31optim = torch.optim.SGD(model.parameters(), lr=learning_rate)

32

33# Training loop

34for epoch in range(epochs):

35 for idx, (batch_x, batch_y) in enumerate(loader):

36 batch_x, batch_y = batch_x.to(device), batch_y.to(device)

37

38 out = model(batch_x)

39 loss = F.nll_loss(out, batch_y)

40

41 optim.zero_grad()

42 loss.backward()

43 optim.step()

44 print("Train Epoch: {} [{}/{}]\tLoss: {:.6f}"

45 .format(epoch, idx*batch_size, len(dataset), loss))

Listing 1: PyTorch code implementing a two layer feedforward network to
classify mnist

Let us explain what is happening in this example.

###### Data processing

In lines 11 to 18, we download and normalise the dataset. In particular, we
download it automatically, and we use the PyTorch convenience functions to
apply a transform pipe to the points of this dataset. In this case, we are
transforming every image so that the whole dataset has zero mean and unitary
standard deviation. The values $0.1307$ and $0.3081$ are the empirical mean
and variance of the dataset and were computed a priori.

In this part we also wrap the dataset into a loader. The loader is an
abstraction that samples batches of size batch_size for us from the dataset.
We can iterate over it like a list, as we do in Line 35. We also set the
DataLoader to shuffle the data every time we start iterating over it.

In Line 21, we check whether the computer has a gpu available. If so, we will
be computing every operation in a gpu. If there is no gpu available, we fall
back to the cpu.

###### Model specification

In Line 24, we define a Sequential model. This model is one that splits
exactly as a sequential composition of functions, as defined in Definition
2.3.1. We start by flattening the input, that is, applying the linear
isomorphism $\mathbb{R}^{28\times 28}\cong\mathbb{R}^{784}$. After that, we
apply a linear transformation (cf., Definition 2.3.3) together with the
$\texttt{relu}(x)=\max\lparen x,0\rparen$ non-linearity, where the $\max$ is
taken coordinate-wise. Finally, in Line 30 we move the model to gpu if there
was one available. If not, this line does not do anything. Being able to
switch between cpu and gpu with a single line of code is one of the strengths
of these modern machine learning libraries.

In Line 31, we define the optimiser to be sgd with a learning rate of
$\eta=10^{-3}$.

###### Forward pass

The gist of the code starts in the training loop. The training loop iterates
over the whole dataset $80$ times—also called epochs—in batches of size $128$.
As we had mentioned, we can iterate over the loader, which returns in every
step a batch of size $128$ of images and labels. In particular, we have that
batch_x has dimensions $(128,28,28)$ while batch_y is a tensor of one
dimension of size $128$ of integers. As we did with the model, if there was an
available gpu, we move the data from cpu to gpu.

In Line 38, we compute the image under our model of every point of the
dataset. Using our previous notation, if we have a batch
$\\{(X_{i},Y_{i})\\}$, this is equivalent to computing
$\\{f_{\theta}(X_{i})\\}$. As such, the output out is of size $(128,10)$,
where every vector out[i] is the component-wise logarithm of a point on the
$10$-dimensional simplex.

In Line 39, we compute the average negative log-likelihood loss over the
minibatch using the logarithmic probabilities and the true labels. As such,
the variable loss is just a floating-point number encapsulated in a PyTorch
class, that is, a tensor of dimension $1$ holding the value of the empirical
risk of the minibatch $\widehat{R}_{B}(\theta)$ as defined in (2.7).

Lines 38 and 39 do much more heavy work than it looks at first sight. When
these lines are called, each operation in model, that is, every matrix
multiplication, every addition of a bias vector, every coordinate-wise non-
linearity, etc., is recorded into a graph called the computation graph as it
is being executed. This graph contains the dependencies between all the
variables in the program as a dependency tree, where the parameters of the
model are in the leaves and the inner nodes of this tree are operations
putting them together. The computation of $f_{\theta}(X_{i})$ and the creation
of the computation graph is what is often called the forward pass.

###### Backward pass

The computation of the gradients and the update of the parameters happens in
the last three lines of the loop. As we have seen, in Line 31 we have passed
the parameters of our model to the optimiser. Line 41 initialises the
optimiser, zeroing-out the gradients of all the parameters that it controls,
this is just a technical point. Line 42 computes the gradients traversing the
computation graph backwards using the chain rule—in this context often called
the adjoint method. Finally, Line 43 updates the variables with the update
step from the optimisation algorithm (cf., Equation 2.8).

##  Chapter 3 Differential and Riemannian Geometry

This section is intended as an introduction to some results and computations
in differential geometry and Riemannian geometry that will be useful later on.
We will skip the proofs of some constructions at the beginning for the sake of
the exposition, as they are quite technical, but we will give references to
what, in our opinion, are some particularly clean proofs of them. We will
always assume all the objects to be in the smooth category, unless stated
otherwise.

###### Outline of the chapter.

The problem motivating this whole chapter is that of computing geodesics on a
large class of manifolds—naturally reductive homogeneous spaces—as this will
be of interest in the rest of the thesis. The chapter is organised going from
more general objects to more particular ones, stopping at each layer looking
into what it adds to the global picture. We finish the chapter in Section 3.5
showing how to instantiate the constructions in the chapter in a number of
manifolds commonly used in optimisation.

Even though all the ideas in this chapter are well-known to geometers, we are
not aware of any source that contains all of them. It is for this reason that
we include at the beginning of each section a short comment on both historical
references and standard references. We also include the proofs necessary to
derive the equations for the geodesics in the studied spaces.

###  3.1 Fibred Manifolds

In this section we provide an introduction to the theory of submersions and
Riemannian submersions. This can be seen as a first approach to the theory of
Ehresmann connections on fibre bundles. The theory of connections on a fibred
manifold as distributions on the total space was first developed in [52],
while most of the formulae relating the geometry of the base space to the
geometry of the total space on a Riemannian submersion is due to O’Neill
[141]. A more detailed review of this subject can be found in [21, Chapter 9].

We start by recalling the definition of a submersion.

######  Definition 3.1.1.

Let $\pi\colon\overline{M}\to M$ be a map between manifolds. We say that $\pi$
is a submersion if $\lparen\mathrm{d}\pi\rparen_{p}\colon T_{p}\overline{M}\to
T_{\pi(p)}M$ is surjective for every $p\in\overline{M}$. In particular, if
$\pi$ is surjective, this is equivalent to saying that $\mathrm{d}\pi\colon
T\overline{M}\to TM$ is surjective. We say that $\overline{M}$ is the total
manifold, and $M$ is the base manifold. We will assume submersions to be
surjective unless stated otherwise.

The fibre over $x\in M$ is defined as $\overline{M}_{x}\coloneqq\pi^{-1}(x)$.
By the implicit function theorem, the fibres of a submersion are manifolds of
dimension $\dim\lparen\overline{M}\rparen-\dim\lparen M\rparen$. A pair of
manifolds together with a submersion between them is called a fibred manifold.

######  Notation.

Even though the notation may suggest otherwise, the manifold $\overline{M}$ is
completely unrelated to the manifold $M$. We choose this notation for fibred
manifold as it is possible to identify vector fields $X$ on $M$ with certain
class of vector fields $\overline{X}$ on $\overline{M}$. This notation makes
this identification easier to follow.

If all the fibres are diffeomorphic, then we get the concept of a space that
locally looks like a product of the base space and its fibre. This idea is
formalised through fibre bundles.

######  Definition 3.1.2.

Let $\pi\colon E\to M$ be a fibred manifold. Then $\pi\colon E\to M$ is a
fibre bundle over $M$ with fibre $F$ if for every $x\in M$ there is an open
neighbourhood $U$ and a map $\alpha\colon\pi^{-1}(U)\to F$ such that

$(\pi,\alpha)\colon\pi^{-1}(U)\to U\times F$

is a diffeomorphism. We say that $(\pi,\alpha)$ is a bundle chart. For every
$x\in M$ we have that $\alpha\mid_{x}\colon E_{x}\to F$ is a diffeomorphism.

All the fibred manifolds we will work with will be fibre bundles.

A key property of fibred manifolds is that they admit local smooth sections.
These are mainly used to show that objects defined on fibre bundles are
smooth.

######  Proposition 3.1.3 (Existence of local sections).

Let $\pi\colon\overline{M}\to M$ a fibred manifold. Then, for every
$p\in\overline{M}$, there exists a neighbourhood $U\subseteq\overline{M}$, and
a smooth function $\sigma\colon\pi(U)\to U$ such that
$\pi\circ\sigma=\mathrm{Id}_{\pi(U)}$. This is called a local section through
$p$.

###### Proof.

See for example [108, Ch. II Prop. 2.2] ∎

######  Remark 3.1.4 (Fibred spaces).

It is in fact possible to show via a simple application of the inverse
function theorem that submersions are characterised by the property of having
smooth sections. This allows us to generalise the concept of fibred manifolds
to the category of topological manifolds simply by swapping smooth sections
with continuous sections. This is of particular interest when studying
stratified spaces, such as low-rank matrices

$\mathbb{R}^{n\times k}_{\leq r}=\\{X\in\mathbb{R}^{n\times
k}\>|\>\operatorname{rank}(X)\leq r\\}.$

In this case, we have that this is not a manifold, but a collection of
manifolds—the collection of the manifolds of fixed-rank matrices (Section
3.5.4)—called strata, with some homogeneity properties. These spaces appear
naturally when considering decompositions such as the qr decomposition, which
is smooth on each of the strata of this topological manifold. Sadly, we will
not have time to explore this connection between topology, numerical analysis,
and optimisation on manifolds in this thesis.

######  Remark 3.1.5 (Global sections).

Global sections need not exist in general. In particular, we may not find in
general an immersed copy of the base space inside the total space. To see
this, consider the following circle bundle over the round sphere:

$U\SS^{2}\coloneqq\\{(x,v)\in T\SS^{2}\>|\>\lVert v\rVert=1\\}.$

It is not difficult to show that this is a fibre bundle over the sphere with
fibre $\SS^{1}$. On the other hand, by the hairy ball theorem, it has no
global sections, as a global section would give an everywhere non-zero vector
field over the $2$-sphere. This is a particular example of a general behaviour
of non-trivial principal bundles. Note that this is an example of a bundle
with no global section, but we can still find an immersed $\SS^{2}$ inside its
total space.

The study of the existence or lack of existence of global sections in fibred
spaces and fibre bundles is of paramount importance in the field of algebraic
topology. An introductory exposition to fibre bundles in this context can be
found in [162].

For a submersion, we define the vertical bundle as
$V\overline{M}\coloneqq\ker\mathrm{d}\pi\subseteq T\overline{M}$. This is a
vector subbundle of $T\overline{M}$.

If we also have a metric $\overline{\textsl{g}}$ on $\overline{M}$, we can
define the horizontal bundle $\mathcal{H}$ as the fibre-wise orthogonal
complement with respect to $\overline{\textsl{g}}$ of $V\overline{M}$, that
is,

$\mathcal{H}_{p}\coloneqq\lparen
V_{p}\overline{M}\rparen^{\bot}=\lparen\ker(\mathrm{d}\pi)_{p}\rparen^{\bot}=\\{X\in
T_{p}\overline{M}\>|\>\overline{\textsl{g}}(X,Y)=0,\,\forall Y\in
V_{p}\overline{M}\\}.$

These two vector bundles form a split
$T\overline{M}=V\overline{M}\oplus\mathcal{H}$. For an element $\zeta\in
T\overline{M}$ we denote its split into its horizontal and vertical components
by $\zeta=\zeta_{V}+\zeta_{\mathcal{H}}$.

Note that since $\lparen\mathrm{d}\pi\rparen_{p}$ is a surjective linear map,
$\lparen\mathrm{d}\pi\rparen_{p}\colon\mathcal{H}_{p}\to T_{\pi(p)}M$ is a
linear isomorphism of vector spaces.

In particular, $\mathcal{H}$ is a smooth choice of a complement of
$V\overline{M}$. Note that the only property that we have used of the metric
$\overline{\textsl{g}}$ is that it is non-degenerate, so that the complement
$\lparen V\overline{M}\rparen^{\perp}$ together with $V\overline{M}$ span all
of $T_{p}\overline{M}$. We can abstract these properties to define a general
connection on a fibred manifold.

######  Definition 3.1.6.

We say that a vector subbundle $\mathcal{H}$ of $T\overline{M}$ is an
Ehresmann connection on a fibred manifold $\pi\colon\overline{M}\to M$ if it
is a complement of the vertical bundle, that is,
$T\overline{M}=V\overline{M}\oplus\mathcal{H}$.

The split given by an Ehresmann connection makes
$\mathrm{d}\pi|_{\mathcal{H}}$ into a vector bundle morphism along $\pi$, in
the sense that the following diagram commutes

$\mathcal{H}$$TM$$\overline{M}$$M$$\mathrm{d}\pi|_{\mathcal{H}}$$\pi_{T\overline{M}}|_{\mathcal{H}}$$\pi_{TM}$$\pi$
(3.1)

and $\mathrm{d}\pi|_{\mathcal{H}}$ is a fibre-wise linear isomorphism.

######  Definition 3.1.7.

Let $\pi\colon\overline{M}\to M$ be a fibred manifold with an Ehresmann
connection $\mathcal{H}$. A vector field $X$ on $M$ can be uniquely identified
with a vector field $\overline{X}$ on $\overline{M}$ with values on
$\mathcal{H}$. We call this vector field the horizontal lift of $X$. In
symbols,

$\overline{X}_{p}\coloneqq\lparen\mathrm{d}\pi|_{\mathcal{H}_{p}}\rparen^{-1}\lparen
X_{\pi(p)}\rparen\mathrlap{\qquad\forall p\in\overline{M}.}$

This lift gives a projectable vector field, meaning that
$\mathrm{d}\pi(\overline{X})=X$ is well-defined. This construction can be
readily generalised to any tensor bundle of $TM$, being able to lift not only
vector fields, but also any kind of tensor. In technical words, a connection
induces connections on all the associated bundles.

######  Definition 3.1.8.

Let $(\overline{M},\overline{\textsl{g}}),(M,\textsl{g})$ be two Riemannian
manifolds and let $\pi\colon\overline{M}\to M$ be a submersion between them
with induced Ehresmann connection $\mathcal{H}$. We say that $\pi$ is a
Riemannian submersion if for every $p\in\overline{M}$

$\mathrm{d}\pi|_{\mathcal{H}_{p}}\colon\mathcal{H}_{p}\to T_{\pi(p)}M$

is a linear isometry. In other words,

$\overline{\textsl{g}}_{p}(\zeta_{1},\zeta_{2})=\textsl{g}_{\pi\lparen
p\rparen}\lparen\lparen\mathrm{d}\pi\rparen_{p}(\zeta_{1}),\lparen\mathrm{d}\pi\rparen_{p}(\zeta_{2})\rparen\mathrlap{\qquad\forall\zeta_{1},\zeta_{2}\in\mathcal{H}_{p}.}$

We may equivalently write this condition in terms of horizontal vectors as

$\overline{\textsl{g}}_{p}(\overline{u},\overline{v})=\textsl{g}_{\pi\lparen
p\rparen}\lparen u,v\rparen\mathrlap{\qquad\forall u,v\in T_{\pi(p)}M.}$

######  Remark 3.1.9 (A Riemannian submersion does not restrict to a local
isometry).

By the definition above, together with the existence of local sections of
$\pi$, one might be tempted to think that for a point $x\in M$ there exists a
section $s\colon\Sigma\to s(\Sigma)$ through $p\in\pi^{-1}(x)$ such that
$\pi|_{s(\Sigma)}$ is an isometry. This is not true as, by the local version
of Frobenius theorem [108, Ch. VI Thm. 1.1], $\mathcal{H}$ is only
integrable—there exists a submanifold $N\subseteq\overline{M}$ with
$T_{p}N=\mathcal{H}_{p}$ for $p\in N$—whenever it is involutive around $p$,
that is, it is closed under the Lie bracket. This means that we may only find
a section giving an immersed manifold $s(\Sigma)$ with tangent bundle
$\mathrm{d}s(\Sigma)\subseteq\mathcal{H}$ whenever $\mathcal{H}$ is
involutive. But $\mathcal{H}$ is often not involutive. In fact, the curvature
of $\mathcal{H}$ is defined as how much $\mathcal{H}$ deviates from being
involutive. Hence, we may only locally integrate $\mathcal{H}$ when the
connection is locally flat. In this case, $\overline{M}$ is locally a
Riemannian product of manifolds.

If we have a fibred manifold $\pi\colon\overline{M}\to M$ with an Ehresmann
connection $\mathcal{H}$ and a metric g on $M$, we may lift g to a metric
tensor on $\mathcal{H}$ by setting
$\overline{\textsl{g}}_{\mathcal{H}}\coloneqq\lparen\mathrm{d}\pi|_{\mathcal{H}}\rparen^{\ast}\lparen\textsl{g}\rparen$
where $\lparen-\rparen^{\ast}$ denotes the pullback

$\overline{\textsl{g}}_{\mathcal{H}_{p}}(u,v)\coloneqq\textsl{g}_{\pi(p)}\lparen\mathrm{d}\pi\lparen
u\rparen,\mathrm{d}\pi\lparen v\rparen\rparen\mathrlap{\qquad\forall
u,v\in\mathcal{H}_{p}.}$

Using the decomposition $T\overline{M}=V\overline{M}\oplus\mathcal{H}$, we can
choose a metric tensor $\overline{\textsl{g}}_{V}$ on $V\overline{M}$ and
define a metric on $\overline{M}$ declaring the horizontal and vertical
bundles to be orthogonal, that is,
$\overline{\textsl{g}}\coloneqq\overline{\textsl{g}}_{V}\oplus\overline{\textsl{g}}_{\mathcal{H}}$.
This metric makes $\pi\colon\overline{M}\to M$ into a Riemannian submersion.
We can summarise this discussion as follows:

######  Proposition 3.1.10 (Lifting a metric on a fibred manifold).

Let $\pi\colon\overline{M}\to M$ be a fibred manifold. Given a metric g on
$M$, an Ehresmann connection $\mathcal{H}$ and a metric tensor
$\overline{\textsl{g}}_{V}$ on $V\overline{M}$ we can form a metric that turns
$\pi$ into a Riemannian submersion by setting
$\overline{\textsl{g}}=\overline{\textsl{g}}_{V}\oplus\lparen\mathrm{d}\pi|_{\mathcal{H}}\rparen^{\ast}(\textsl{g})$.

Conversely, if we have a metric $\overline{\textsl{g}}$ on $\overline{M}$ that
is
projectable—$\overline{\textsl{g}}_{p_{1}}(\overline{u},\overline{v})=\overline{\textsl{g}}_{p_{2}}\lparen\overline{u},\overline{v}\rparen$
for every $p_{1},p_{2}\in\overline{M}_{x},\,u,v\in T_{x}M$—we may push it
forward to form a metric on $M$ such that $\pi$ is a Riemannian submersion.

We will show in the next sections how to construct metrics on $V\overline{M}$
in the particular case of principal bundles and homogeneous spaces, where the
fibre has the structure of a Lie group.

Riemannian submersions generalise the idea of an isometry between spaces of
the same dimension. In particular, they provide a strong link between the
geometry of $\overline{M}$ and $M$.

######  Proposition 3.1.11.

Let $\pi\colon\overline{M}\to M$ be a Riemannian submersion. We have that:

  1. 1.

A geodesic in $\overline{M}$ with horizontal initial conditions is horizontal
everywhere. $\pi$ takes horizontal geodesics on $\overline{M}$ to geodesics on
$M$. In particular, if $(\overline{M},\overline{\textsl{g}})$ is complete, so
is $(M,\textsl{g})$.

  2. 2.

Any geodesic on $M$ may be locally lifted to a geodesic on $\overline{M}$ with
horizontal initial conditions.

  3. 3.

$\pi$ is distance decreasing: $\overline{d}\lparen p,q\rparen\geq
d\lparen\pi\lparen p\rparen,\pi\lparen q\rparen\rparen$, $\forall
p,q\in\overline{M}$.

  4. 4.

[78] If $(\overline{M},\overline{\textsl{g}})$ is complete, any curve in $M$
may be lifted globally to a curve in $M$. In particular, any geodesic in
$(M,\textsl{g})$ may be lifted globally to a horizontal geodesic of
$(\overline{M},\overline{\textsl{g}})$.

  5. 5.

[141] For $X,Y$ vector fields on $M$ we have

$\displaystyle\nabla_{\overline{X}}\overline{Y}$
$\displaystyle=\overline{\nabla_{X}Y}+\mbox{{}\small$\displaystyle\frac{1}{2}$}[\overline{X},\overline{Y}]_{V}$
$\displaystyle[\overline{X},\overline{Y}]_{\mathcal{H}}$
$\displaystyle=\overline{[X,Y]}.$

###### Proof.

Item 1 can be found in [142, Ch. 7 Corol. 45]. Item 2 is a corollary of the
formula for the connection in Item 5. Item 3 is a direct corollary of Item 1.
Item 4 is proved in [78, Prop. $3.2$]. Item 5 is proved in [141, Lemma $3$ and
Lemma $1.2$] ∎

According to this proposition, if we have a Riemannian submersion, and we know
how to compute geodesics with horizontal initial conditions on the total
space, we can compute geodesic on the base space. Furthermore, since
$\mathrm{d}\pi$ is surjective, any geodesic on the base space is of this form.

Before finishing this section, we look at Riemannian submersions with totally
geodesic fibres, as these will be rather important in the sequel. We start by
defining a totally geodesic submanifold.

######  Definition 3.1.12.

An immersed submanifold $(F,\textsl{g}_{F})$ of a Riemannian manifold
$(M,\textsl{g}_{M})$ is totally geodesic if geodesics on $(F,\textsl{g}_{F})$
are also geodesics on $(M,\textsl{g}_{M})$.

######  Example 3.1.13.

On a round sphere maximum circles are geodesic submanifolds but any other
circle is not.

The deviation of an immersed submanifold from being a geodesic submanifold is
given by the shape tensor sometimes also called the second fundamental form.

######  Definition 3.1.14.

Given an immersed submanifold $(F,\textsl{g}_{F})$ of a Riemannian manifold
$(M,\textsl{g}_{M})$ with Levi-Civita connections $\nabla^{F},\nabla^{M}$
respectively, we define the shape tensor of $F$ for two vector fields $X,Y$ on
$F$ as the $(2,1)$ tensor

$II(X,Y)=\nabla^{F}_{X}Y-\nabla^{M}_{X}Y.$

A simple computation involving Gauss’ equation gives the following
characterisation of totally geodesic submanifolds.

######  Proposition 3.1.15.

A Riemannian submanifold $(F,\textsl{g}_{F})$ of $(M,\textsl{g})$ is totally
geodesic if and only if $II\equiv 0$.

###### Proof.

See [142, Ch. 4 Prop. 13]. ∎

Totally geodesic submanifolds will be important when we talk about Riemannian
submersions with totally geodesic fibres. Most of the Riemannian submersions
that we will encounter in the optimisation context will have this property.

In the next sections we will look at two problems: How to construct Riemannian
submersions and how to compute geodesics on some total spaces of these
submersions. Once we know how to do this, we may repeat this process forming a
chain of submersions

$M_{1}\xrightarrow{\pi_{1}}M_{2}\xrightarrow{\pi_{2}}M_{3}\xrightarrow{\pi_{3}}\dots,$

where, if we know how to compute geodesics on $M_{1}$, and every map is a
Riemannian submersion, we know how to compute geodesics on all the other
manifolds in the chain. This will be particularly useful for optimisation with
orthogonal constraints as shown in Section 3.5, where we will see that the
special orthogonal group, the Stiefel manifold, and the Grassmannian form such
a tower. We will also leverage this in Chapter 6 to implement optimisation on
different manifolds by implementing it on one manifold, and then pushing the
geodesics forward along Riemannian submersions.

###  3.2 Invariant Metrics and Principal Bundles

In this section we will study manifolds together with a group action. In
particular, we will be interested in studying conditions for the existence of
an invariant metric under this action. This will take us naturally to the
definition of principal bundles which, in turn, will provide us with a way to
construct Riemannian submersions. Most of the theory presented in this section
and the next one can be generalised to study connections on principal and
associated bundles. The presentation in this section differs from the
classical one in that we emphasise the Riemannian side of principal bundles as
spaces on which we can find a metric invariant under a Lie group action that
descends into the quotient of the manifold by the action.

Principal bundles were first proposed in full generality by Palais, who
introduced the study of proper actions on manifolds, as a generalisation of
actions by compact groups [145]. The results of this section can mostly be
found across the books kobayashi1963foundationskobayashi1969foundations,
written by two of the fathers of the subject.

We start by defining some notation and recalling some rudiments of Lie groups.

######  Notation.

Let $\mu\colon G\times M\to M$ be an action of a Lie group $G$ on a manifold
$M$. We will write $\mu_{g}(x)=\mu^{x}(g)=\mu(g,x)$ for the different partial
maps defined by it. In particular, $\mu_{g}$ is a diffeomorphism with inverse
$\mu_{g^{-1}}$.

We will denote left and right actions of a Lie group onto itself by $L_{g}$
and $R_{g}$, and the conjugacy action of a Lie group by $c_{g}\coloneqq
L_{g}\circ R_{g^{-1}}$.

For a Lie group $G,H$, we will denote their Lie algebras as
$\mathfrak{g},\mathfrak{h}$. Given a vector $X\in\mathfrak{g}$, we denote the
left-invariant vector field associated to $X$ by

${X}^{\texttt{L}}_{g}\coloneqq\lparen\mathrm{d}L_{g}\rparen_{e}(X).$

Left-invariant vector fields provide a linear isomorphism between $T_{g}G$ and
$\mathfrak{g}$ and they allow us to define left-invariant tensors on the whole
manifold by defining them just on the Lie algebra. For example, if we have an
inner product $\textsl{g}_{e}$ on the Lie algebra, we can define a left-
invariant metric as

${\textsl{g}}^{\texttt{L}}_{g}\lparen{X}^{\texttt{L}}_{g},{Y}^{\texttt{L}}_{g}\rparen\coloneqq\textsl{g}_{e}\lparen
X,Y\rparen\mathrlap{\qquad X,Y\in\mathfrak{g}.}$

We recall that the Lie exponential is defined as the integral of left-
invariant vector fields, that is, if $v\in\mathfrak{g}$,
$\operatorname{expm}(tX)=\gamma(t)$ we have that

$\dot{\gamma}(t)\coloneqq{X}^{\texttt{L}}_{\gamma(t)}.$

We will denote it by $\operatorname{expm}$, to differentiate it from the
Riemannian exponential map.

We denote the orbit of a point $x\in M$ as

$G(x)\coloneqq\\{\mu_{g}(x)\>|\>g\in G\\},$

and the stabiliser of $x$, also referred to as the isotropy group at $x$, as

$G_{x}\coloneqq\\{g\in G\>|\>\mu_{g}(x)=x\\}.$

Note that for continuous actions $G_{x}=\lparen\mu^{x}\rparen^{-1}(\\{x\\})$
is closed.

######  Remark 3.2.1 (Left and right actions).

In the sequel we will be using left and right actions. The use of one or the
other is merely due to convenience, since we may turn a left action $\mu_{g}$
into a right action by considering the map $g\mapsto\mu_{g^{-1}}$. We will
reserve right actions for actions we use to quotient by, and left actions to
talk about actions that act on this quotient.

We are interested in studying metrics that are invariant under the action of a
group, as these metrics will descend into the quotient of the manifold by the
group action.

######  Definition 3.2.2.

Given a manifold $M$ together with a right action $\mu$ by a Lie group $G$, we
say that a metric g on $M$ is $G$-invariant if

$\textsl{g}_{\mu_{g}\lparen
p\rparen}\lparen\mathrm{d}\mu_{g}(u),\mathrm{d}\mu_{g}(v)\rparen=\textsl{g}_{p}\lparen
u,v\rparen\mathrlap{\qquad\forall g\in G,p\in M,u,v\in T_{p}M,}$

in other words, the action $\mu_{g}$ is an isometry for every $g\in G$.
Furthermore, if the group of isometries $\mu_{G}=\\{\mu_{g}\>|\>g\in G\\}$ is
a closed subgroup of $\operatorname{Iso}(M,\textsl{g})$, we say that $G$ acts
by isometries.

######  Remark 3.2.3 (Effective actions).

Note that $\mu_{G}$ is isomorphic to $G$ if $\mu_{g}=\mathrm{Id}$ implies that
$g=e$. We call these effective actions. If an action is not effective, we can
define the ineffective kernel of the action $N=\bigcap_{x\in M}G_{x}$. This is
a normal subgroup of $G$, and in this setting $\mu_{G}$ is isomorphic to $G/N$
as Lie groups. In particular, the action of $G/N$ on the manifold $M$ is
equivalent to that of $G$, so in most situations we can safely assume that the
action is effective, although at times it is convenient to work with non-
effective actions, when doing explicit computations in concrete manifolds.

######  Remark 3.2.4 (Closed subgroups of isometries).

By a theorem of [135], the isometry group $\operatorname{Iso}(M,\textsl{g})$
is a Lie group. Hence, the requirement that $\mu_{G}$ is a closed subgroup of
$\operatorname{Iso}(M,\textsl{g})$ implies that $\mu_{G}$ is itself a Lie
group. We will see in shortly that having a closed subgroup of isometries will
allow us to _quotient_ by it—i.e., putting a unique smooth structure on the
space of orbits such that it makes the projection a submersion—which may not
be possible when the group is not closed.

The first kind of actions that we will consider are proper actions. An action
is said to be proper if the map from $G\times M$ to $M\times M$,
$(g,x)\mapsto(\mu_{g}(x),x)$ is proper, in the sense that the preimage of a
compact set is compact. Proper actions give us a way to construct
$G$-invariant metrics. Even though this result goes back to Palais, it was
recently shown that this metric can be chosen to be complete.

######  Theorem 3.2.5 ([95]).

Let $M$ be a manifold with a proper $G$-action $\mu$. Then, there exists a
complete metric g which makes $\mu_{G}$ into a closed subgroup of
$\operatorname{Iso}(M,\textsl{g})$.

The second kind of actions that we will look at are free actions. A free
action is one that has no non-trivial fixed points. In other words, if
$\mu_{g}(x)=x$ for an $x\in M$, then $g=e$. Another way of looking at this is
that every isotropy group is trivial. Free isometric actions are actions by
which we can quotient and obtain a smooth manifold.

######  Theorem 3.2.6 (Free slice theorem).

Let $G$ be a closed subgroup of $\operatorname{Iso}(M,\textsl{g})$ that acts
freely. Then, the set of orbits $M/G=\\{G(x)\>|\>x\in M\\}$ can be given a
unique manifold structure and a Riemannian metric such that $\pi\colon M\to
M/G$ is a Riemannian submersion.

###### Proof.

See [147, Theorem 5.6.21]. ∎

Combining these two theorems we have the following corollary.

######  Corollary 3.2.7.

Let $M$ be a manifold together with a free proper action. There exists a
$G$-invariant metric on $M$ that descends into a metric on $M/G$ which makes
$\pi\colon M\to M/G$ into a Riemannian submersion.

Proper actions are a reasonably general setting on which to study invariant
metrics.

######  Theorem 3.2.8 (Isometric actions are proper).

Let $(M,\textsl{g})$ be a Riemannian manifold and $G$ be a closed subgroup of
$\operatorname{Iso}(M,\textsl{g})$. Then, its action on $M$ is proper.

###### Proof.

See [4, Theorem 3.62]. ∎

######  Remark 3.2.9 (Isometric actions vs. effective and proper).

Theorem 3.2.8 gives a converse to Theorem 3.2.5. Between these two theorems,
we have an equivalence between effective and proper $G$-actions on a manifold
$M$ and closed subgroups $G$ of $\operatorname{Iso}\lparen
M,\textsl{g}\rparen$ where g is a complete metric.

######  Remark 3.2.10 (Non-free isometric actions).

In some situations, it is not strictly necessary for the action to be free but
it certainly makes things easier. If the action is not free, one still has an
open and dense subset of the manifold whose quotient by the action is a
manifold. This is part of the content of the principal orbit theorem, and the
start of the theory of Riemannian foliations [4, Thm. $3.82$].

Let us now define principal bundles, which are fibre bundles whose fibre is a
Lie group, together with a compatible fibre-preserving action.

######  Definition 3.2.11.

Let $\pi\colon P\to M$ be a fibre bundle with fibre a Lie group $G$. We say
that $P$ is a principal bundle if we have a free and transitive (on the
fibres) fibre-preserving right action $\mu\colon P\times G\to P$, and an atlas
such that the bundle charts $\alpha\colon\pi^{-1}(U)\to G$ are equivariant
with respect to the action,

$\alpha(\mu(p,g))=\alpha(p)g\mathrlap{\qquad\forall p\in\pi^{-1}(U),g\in G.}$

This makes the charts compatible with the action $\mu$ on $P$ and the action
of multiplication on the right on $G$.

We now show that principal bundles are the right abstraction to look at
metrics that are $G$-invariant by which we can quotient.

######  Theorem 3.2.12 (Construction of principal bundles).

Let $M$ be a manifold with a free and proper $G$-action. Then, $\pi\colon M\to
M/G$ has the structure of a principal bundle.

###### Proof.

See [4, Theorem 3.34]. A sketch of the proof goes as follows:

Since the action is free and proper, by Theorem 3.2.6 the map $\pi\colon M\to
M/G$ is a submersion, so every fibre $M_{x}\coloneqq\pi^{-1}(\\{x\\})$ is a
manifold. These fibres are, by definition, the orbits of the action and, since
the action is free, they are all isomorphic to $G$, so a principal bundle is a
fibre bundle with fibre $G$ and a fibre-preserving free and proper action.
Adapting the proof of Theorem 3.2.6 one proves that the bundle charts are
equivariant. ∎

And we have the following converse.

######  Proposition 3.2.13 (Characterisation of principal bundles).

The action $\mu$ of a principal bundle is free and proper.

###### Proof.

See [4, Proposition 3.33]. These properties follow from the fact that an
action of a Lie group on itself is free and proper, and using the equivariance
of the bundle charts. ∎

In particular, every principal bundle admits an invariant metric that descends
into the base space and makes $\pi$ into a Riemannian submersion.

Until now, we have argued that principal bundles are spaces that accept
$G$-invariant metrics, and that if we have a $G$-invariant metric on the total
space, it descends into a metric on the base space. We will now see some sort
of converse to this result. In particular, we will see how to lift a metric
from the base space into the total space such that the metric on $P$ is
$G$-invariant. In general, we described in Section 3.1 how, given a metric on
a base space and a connection on the total space, we may lift it to a metric
on the total space of a fibred manifold by declaring $\mathrm{d}\pi$ to be an
isometry and choosing a metric tensor on the vertical bundle. In the case of
principal bundles, the fact that the fibres look like Lie groups allows us to
simplify the process of choosing a metric tensor on $VP$ to choosing one on
the Lie algebra.

We start by defining a connection on a principal bundle.

######  Definition 3.2.14.

A connection on a principal $G$-bundle $\pi\colon P\to M$ (or simply a
principal connection) $\mathcal{H}\subseteq TP$ is a $G$-equivariant Ehresmann
connection, that is, a $G$-equivariant sub-bundle complementary to $VP$. In
symbols,

$\lparen\mathrm{d}\mu_{g}\rparen_{p}(\mathcal{H}_{p})=\mathcal{H}_{\mu(p,g)}\mathrlap{\qquad\forall
p\in P,g\in G.}$

It is direct to see that a $G$-invariant metric on $P$ induces a connection on
$P$ defining $\mathcal{H}_{p}\coloneqq\lparen V_{p}P\rparen^{\bot}$, and it is
exactly the $G$-invariance of the metric what gives the equivariance. In
particular, since principal bundles accept $G$-invariant metrics, every
principal bundle admits a principal connection. We may use these to lift a
metric from $M$ to $P$ making $\pi\colon P\to M$ into a Riemannian submersion.

######  Theorem 3.2.15 (Lifting a metric on a principal bundle [174]).

Let $\pi\colon P\to M$ be a principal $G$-bundle. Given a metric g on $M$, an
inner product $\textsl{g}_{e}$ on $\mathfrak{g}$ and a principal connection
$\mathcal{H}$, there exists a unique metric $\overline{\textsl{g}}$ on $P$
such that $\pi$ is a Riemannian submersion with totally geodesic fibres
isometric to $(G,{\textsl{g}_{e}}^{\texttt{L}})$—the left-invariant metric on
$G$ associated to $\textsl{g}_{e}$—and horizontal distribution $\mathcal{H}$.

###### Proof.

Given the principal connection $\mathcal{H}$, we can put a metric on the
horizontal space by pulling back the metric on $M$ as
$\textsl{g}_{\mathcal{H}}=\lparen\mathrm{d}\pi|_{\mathcal{H}}\rparen^{\ast}(\textsl{g})$.
We can then declare $\mathcal{H}$ and $VP$ to be orthogonal, and put a metric
on $VP$ by transferring the right-invariant metric on $G$ to $P$ via the
bundle charts. If $\alpha\colon\pi^{-1}(U)\to G$ is a bundle chart, the metric
is defined locally on $\pi^{-1}(U)$ as

$\overline{\textsl{g}}|_{\pi^{-1}(U)}=\alpha^{\ast}({\textsl{g}_{e}}^{\texttt{L}})\oplus\lparen\mathrm{d}\pi|_{\mathcal{H}\cap
T\pi^{-1}(U)}\rparen^{\ast}(\textsl{g}|_{U}).$

This is well-defined as, since the charts are right-equivariant, we have that
for any two bundle charts $\alpha_{U},\alpha_{V}$ with $U\cap V\neq\emptyset$,

$\alpha_{V}\circ\alpha_{U}^{-1}=L_{g}$

for some $g\in G$ (see e.g., [136, Lemma $4.2.1$]).

The fact that the fibres in this construction are totally geodesic
submanifolds of $P$ is proved in [174, Theorem $3.5$].

The uniqueness of this construction comes from the fact that the fibres are
totally geodesic, so they have to be orthogonal to the connection subbundle. ∎

######  Remark 3.2.16 (Associated bundles).

The construction above goes through in the more general setting of associated
bundles, as proved in the original paper. On the other hand, the simpler case
of principal bundles is enough for what we will need in the sequel.

Let us now define a family of vector fields which will be useful later.

######  Definition 3.2.17.

For a principal $G$-bundle $\pi\colon P\to M$ with action $\mu$ we define the
infinitesimal generators of the action as

$\\!\begin{aligned} \xi\colon P\times\mathfrak{g}&\to VP\\\
(p,v)&\mapsto\xi_{v}(p)\coloneqq\lparen\mathrm{d}\mu^{p}\rparen_{e}(v)\end{aligned}$

######  Example 3.2.18.

In the simplest case of the principal bundle $\pi\colon G\to\\{e\\}$,
$\xi_{v}(g)=\lparen\mathrm{d}L_{g}\rparen_{e}(v)$. In other words, $\xi_{v}$
are just left-invariant vector fields on $G$, that is,
$\xi_{v}={v}^{\texttt{L}}$.

The infinitesimal generators of the action take values on $VP$, since for a
$p\in P$,
$\lparen\mathrm{d}\pi\rparen_{p}\circ\lparen\mathrm{d}\mu^{p}\rparen_{e}=\mathrm{d}\lparen\pi\circ\mu^{p}\rparen_{e}=0$,
where we have used that the action $\mu^{p}$ is fibre-preserving, and thus
$\pi\circ\mu^{p}$ is constant. The infinitesimal generators provide with
bundle-isomorphism between $VP$ and $P\times\mathfrak{g}$, and hence, they
show that the vertical bundle is always trivial.

######  Proposition 3.2.19.

Let $\pi\colon P\to M$ be a principal $G$-bundle, then
$\lparen\mathrm{d}\mu^{p}\rparen_{e}\colon\mathfrak{g}\to V_{p}P$ is a linear
isomorphism, so that the vertical space factorises as $VP\cong
P\times\mathfrak{g}$. We also have that the infinitesimal generators induce a
Lie algebra homomorphism

$\xi_{[u,v]}=[\xi_{u},\xi_{v}]\mathrlap{\qquad\forall u,v\in\mathfrak{g}.}$

###### Proof.

For a $p\in P$, $V_{p}P$ can be identified with the tangent to the orbit of
$G$ at $p$. This orbit is diffeomorphic to $G$, as the action is free, so
$V_{p}P$ and $\mathfrak{g}$ have the same dimension. If we prove that
$\lparen\mathrm{d}\mu^{p}\rparen_{e}$ is injective, then it is a linear
isomorphism. Let $v\in\mathfrak{g}$ such that
$\lparen\mathrm{d}\mu^{p}\rparen_{e}(v)=0$ for $p\neq e$. This is equivalent
to saying that

$\left.\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}$}\right\rvert_{t=0}\mu(p,\operatorname{expm}(tv))=0.$

Thus, $p$ is locally a fixed point of the map $\mu^{\operatorname{expm}(tv)}$,
but since the action is free this only happens when $v=0$ so
$\lparen\mathrm{d}\mu^{p}\rparen_{e}$ is injective.

To prove that $\xi$ is an Lie algebra homomorphism, we have by equivariance
that

$\mu_{g}\circ\mu^{p}(\operatorname{expm}(tv))=\mu^{\mu_{g}(p)}(g^{-1}\operatorname{expm}(tv)g)\quad\forall
v\in\mathfrak{g},g\in G,p\in P$

and differentiating at $t=0$ we get the following infinitesimal equivariance
property:

$\lparen\mathrm{d}\mu_{g}\rparen_{p}(\xi_{v}(p))=\xi_{\operatorname{Ad}_{g^{-1}}(v)}(\mu_{g}(p))\mathrlap{\qquad\forall
g\in G,p\in P,v\in\mathfrak{g}.}$ (3.2)

Finally, letting $g=\operatorname{expm}(-tu)$, substituting $p$ with
$\mu\lparen\operatorname{expm}\lparen tu\rparen,p\rparen$ in (3.2), and
differentiating

$\left.\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}$}\right\rvert_{t=0}\lparen\mu_{\operatorname{expm}(-tu)}\rparen_{\ast}\lparen\xi_{v}\lparen\mu\lparen\operatorname{expm}(tu),p\rparen\rparen\rparen=\left.\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}$}\right\rvert_{t=0}\xi_{\operatorname{Ad}_{\operatorname{expm}\lparen
tu\rparen}(v)}(p)\mathrlap{\qquad\forall p\in P,u,v\in\mathfrak{g}.}$

The left-hand side is the definition of $[\xi_{u},\xi_{v}](p)$, while the
right-hand side is equal to $\xi_{[u,v]}(p)$, after using that the derivative
of a linear map is the map itself. ∎

######  Remark 3.2.20 (Lifting a metric on a principal $G$-bundle for $G$
compact).

In the case when $G$ is compact, the construction of the metric on $VP$ in
Theorem 3.2.15 may be greatly simplified by the use of infinitesimal
generators. As it was shown in Equation 3.2,
$\lparen\mathrm{d}\mu_{g}\rparen_{p}$ may be evaluated using infinitesimal
generators. In particular, it acts as the adjoint action on $\mathfrak{g}$. If
we have an $\operatorname{Ad}(G)$-invariant inner product on $\mathfrak{g}$,
we may then pull it back to a $G$-invariant metric tensor on $VP$. This metric
is much simpler to work with and compute, as it does not make use of the
bundle charts.

As we will see in Proposition 3.3.9, on compact Lie groups one may construct
$\operatorname{Ad}(G)$-invariant metrics, so we may carry this construction
whenever $G$ is compact. Luckily, we will just study the setting of
$G$-compact in the sequel, as all Riemannian homogeneous spaces are of this
form.

###  3.3 Homogeneous Spaces

We now come back to the computation of geodesics on different manifolds. The
discussion in the previous section tells us that, by Corollary 3.2.7, if one
has a $G$-invariant metric on a principal bundle, then this metric descends to
a metric on the quotient $P/G\cong M$, turning the projection into a
Riemannian submersion. Since a Riemannian submersion maps horizontal geodesics
into geodesics, if we know how to compute geodesics on $P$, then we can
compute geodesics on $M$ just by projecting the horizontal ones.

In this section, we will address the questions of how to implement these ideas
in concrete families of manifolds. In particular, we will look at families of
manifolds on which we can compute the geodesics easily, so that we can get the
geodesics on the quotient simply by projecting onto it.

The theory of (left-) invariant metrics on a space with a transitive action
was pioneered by [139], which was the first to realise that, if one has a
transitive action, it can be used to translate geometric problems on a
manifold to algebraic problems on a tangent space. An introduction to
reductive spaces in terms of representations can be found in [103, Ch. X Sec.
$1$ and $2$]. An early introduction to the subject of reductive homogeneous
spaces avoiding representations and delving deeper into the theory of the
canonical connection of a reductive space can be found in [104, Ch. $1$]. The
theory of reductive structures was first systematically studied by
Lichnerowicz in [114]. A concise and self-contained introduction to a number
of the ideas in this section and the next one can also be found in the
excellent text [39, Chapter $3$].

We start by looking at the simplest example. Let $H$ be a closed subgroup of a
Lie group $G$, consider the right action of $H$ on $G$ by right
multiplications. It is not difficult to see that this action is proper and
free, so we have our first specific family of principal bundles.

######  Proposition 3.3.1.

Let $H$ be a closed subgroup of a Lie group $G$ acting on $G$ via right
multiplications. Then, $\pi\colon G\to G/H$ is a principal $H$-bundle.

Manifolds that are a quotient of a Lie group by a closed subgroup are
particularly important and are called homogeneous spaces.

######  Definition 3.3.2.

A manifold $M$ is a $G$-homogeneous space if it is diffeomorphic to a manifold
of the form $G/H$ with $H$ a closed subgroup of a Lie group $G$.

The name of homogeneous spaces comes from the fact that these are exactly the
manifolds that admit transitive actions, as we will show in a second.

######  Definition 3.3.3.

An $G$-action is transitive on $M$ if for every $x,y\in M$, there exists a
$g\in G$ such that $\mu_{g}(x)=y$.

When $G$ is a Lie group acting on a smooth manifold, a transitive action
$\mu_{g}$ allows us to move any point to any other point. In particular, since
$\mu_{g}$ is a diffeomorphism (its inverse is $\mu_{g^{-1}}$), these manifolds
look the same from a smooth point of view at every point.

######  Proposition 3.3.4 (Construction and characterisation of homogeneous
spaces).

Let $H$ be a closed subgroup of $G$. The $G$-action by left multiplications
induced on the quotient $G/H$ is transitive, so $G/H$ is a $G$-homogeneous
space.

Conversely, let $M$ be a $G$-homogeneous space with action $\mu$, and fix a
point $x_{0}\in M$. Then, $M$ is diffeomorphic to the quotient of $G$ by the
stabiliser of $x_{0}$, $H\coloneqq G_{x_{0}}$ by right multiplications,
through the diffeomorphism

$\\!\begin{aligned} \widetilde{\mu}^{x_{0}}\colon G/H&\to M\\\
\pi(g)&\mapsto\mu^{x_{0}}(g)\end{aligned}.$

######  Remark 3.3.5 (Maps that descend to the quotient).

A map $\Phi\colon P\to P$ on a principal bundle map that commutes with the
action $\mu_{g}$ descends into a map on the quotient
$\widetilde{\Phi}\coloneqq\pi\circ\Phi$. We will use this in the setting of
the principal bundle of a homogeneous space $\pi\colon G\to G/H$, with $R_{h}$
for $h\in H$ as the bundle right-action and $\Phi=L_{g}$. The left-action then
descends into a left-action $\widetilde{L}_{g}$ on $G/H$ since left and right
multiplication on a Lie group always commute.

For a homogeneous space $M$, there might be different groups that act
transitively on it. Therefore, its representation as a quotient of Lie groups
may not be unique. To see this, we shall first introduce a large family of
examples for homogeneous spaces.

######  Definition 3.3.6.

A Riemannian manifold $(M,\textsl{g})$ is a Riemannian homogeneous space if
its isometry group $\operatorname{Iso}(M,\textsl{g})$ acts transitively on
$M$.

######  Example 3.3.7 (Non-unique quotient representation of a homogeneous
space).

Let $M=\mathbb{R}^{n}$ with the canonical metric, the full isometry group is
generated by rotations, reflections, and translations, but the subgroup of
translations acts transitively on $\mathbb{R}^{n}$ and is a closed subgroup of
$\operatorname{Iso}(\mathbb{R}^{n},\textsl{g}_{\mathbb{R}^{n}})$. In other
words, both the full isometry group and the subgroup of translations act
transitively on $\mathbb{R}^{n}$. Hence, if we denote the translations on
$\mathbb{R}^{n}$ by $\operatorname{T}\lparen n\rparen$, we have two
representations of $\mathbb{R}^{n}$ as
$\mathbb{R}^{n}\cong\operatorname{T}\lparen
n\rparen\cong\operatorname{Iso}\lparen\mathbb{R}^{n},\textsl{g}_{\mathbb{R}^{n}}\rparen/\operatorname{O}\lparen
n\rparen$, where we have used that the group of translations has a trivial
stabiliser.

Riemannian homogeneous spaces have a particularly simple geometry. In
particular, by definition, for any two points there exists an isometry that
transforms one into the other. Hence, quantities that are invariant under
isometries, such as the curvature tensor or the Levi-Civita connection, may be
computed at all points just by computing them at one point and translating
them using the isometries of the manifold. In other words, from the
perspective of Riemannian geometry, these manifolds look the same at every
point, so it is enough to study them at one of their points—or in a
neighbourhood of one point—to understand their geometry at any given point.

The stabilisers of Riemannian manifolds are always compact. This means that a
Riemannian homogeneous space can always be represented as $G/H$ with $H$
compact.

######  Proposition 3.3.8 (Necessary condition for Riemannian homogeneous
spaces).

Let $(M,\textsl{g})$ be a Riemannian manifold and let $G$ be a closed subgroup
of $\operatorname{Iso}(M,\textsl{g})$, then its stabiliser $H$ at a point
$x_{0}\in M$ is compact.

###### Proof.

By Theorem 3.2.8 an action by isometries is proper, and since
$H=\lparen\mu^{x_{0}}\rparen^{-1}(\\{x_{0}\\})$, $H$ is compact. ∎

The converse is also true. In order to prove it, we first need the following
important construction on Lie groups called the averaging trick.

######  Proposition 3.3.9 (Bi-invariant metric by a compact subgroup).

Let $G$ be a Lie group and $H$ a compact subgroup. Then $G$ admits a
left-$G$-invariant metric that is also right-$H$-invariant

###### Proof.

Let $\omega_{\mathfrak{h}}$ be a non-zero alternating multilinear map on
$\mathfrak{h}$ in $k=\dim\mathfrak{h}$ variables (i.e.,
$\omega_{\mathfrak{h}}\in\bigwedge^{k}\mathfrak{h}^{\ast}$). Since
$\dim\bigwedge^{k}\mathfrak{h}^{\ast}=1$, this form is unique up to a
multiplicative constant. We can then extend it to all of $H$ by pushing it
forward along $R_{h}$, defining $\lparen\omega_{H}\rparen_{h}\coloneqq\lparen
R_{h}\rparen_{\ast}\lparen\omega_{\mathfrak{h}}\rparen$. This makes
$\omega_{H}$ into a right-invariant top-dimensional form on $H$, and its
associated measure is called the right Haar measure on $H$.

Let $\langle\cdot,\cdot\rangle$ be any inner product on $\mathfrak{g}$. We can
turn this inner product into an $\operatorname{Ad}(H)$-invariant inner product
on $\mathfrak{g}$ by averaging it over the elements of the group using the
right Haar measure $\omega_{H}$

$\textsl{g}_{\mathfrak{g}}\lparen
u,v\rparen\coloneqq\int_{H}\langle\operatorname{Ad}_{h}(u),\operatorname{Ad}_{h}(v)\rangle\,\omega_{H}\mathrlap{\qquad
u,v\in\mathfrak{g}.}$

Note that this integral is well-defined since $H$ is compact. The
$\operatorname{Ad}(H)$-invariance follows from the right-invariance of
$\omega$

$\textsl{g}_{\mathfrak{g}}\lparen\operatorname{Ad}_{h^{\prime}}\lparen
u\rparen,\operatorname{Ad}_{h^{\prime}}\lparen
v\rparen\rparen=\int_{H}\langle\operatorname{Ad}_{hh^{\prime}}(u),\operatorname{Ad}_{hh^{\prime}}(v)\rangle\,\omega_{H}=\textsl{g}_{\mathfrak{g}}\lparen
u,v\rparen\mathrlap{\qquad\forall u,v\in\mathfrak{g},\forall h^{\prime}\in
H.}$

Finally, we can extend this inner product to a metric on all of $G$ by pushing
it forward along $L_{g}$, $\textsl{g}_{g}\coloneqq\lparen