# Curvature-Dependant Global Convergence Rates for Optimization on Manifolds
of Bounded Geometry

Mario Lezcano-Casado Mathematical Institute, University of Oxford, United
Kingdom.
Correspondence<EMAIL_ADDRESS>
###### Abstract

We give curvature-dependant convergence rates for the optimization of weakly
convex functions defined on a manifold of $1$-bounded geometry via Riemannian
gradient descent and via the _dynamic trivialization algorithm_. In order to
do this, we give a tighter bound on the norm of the Hessian of the Riemannian
exponential than the previously known. We compute these bounds explicitly for
some manifolds commonly used in the optimization literature such as the
special orthogonal group and the real Grassmannian. Along the way, we present
self-contained proofs of fully general bounds on the norm of the differential
of the exponential map and certain cosine inequalities on manifolds, which are
commonly used in optimization on manifolds.

##  1 Introduction

We are interested in approximating critical points of a—possibly non-
convex—function defined on a manifold

$\min_{x\in M}f(x).$

In particular, we are interested in giving global convergence bounds for such
problems. These problems form a particularly well-behaved subset of the field
of constrained optimization. Examples of the importance of this family of
problems are ubiquitous in engineering, statistics, and machine learning. For
example, they have found applications in image processing with the
Grassmannian for tracking subspaces; analysis of time series in deep learning
with the orthogonal group to avoid vanishing and exploding gradient problems;
Bayesian statistics and machine learning with the manifold of positive
definite (PD) matrices in kernel methods and metric learning; the space of
low-rank and fixed-rank matrices for low-rank models; and the hyperbolic space
for word embeddings, among many others.

The approach that we investigate in this paper is that of using the
exponential map to pullback the problem from the whole manifold to a fixed
tangent space, converting the problem into an unconstrained one over
$T_{p}M\cong\mathbb{R}^{n}$

$\min_{v\in T_{p}M}f(\exp_{p}(v)).$ (1)

For this to be a sensible method, the exponential map has to be surjective, so
we will assume $(M,\textsl{g})$ to be connected and complete throughout the
paper.

##### General pullbacks in constrained optimization.

The idea of pulling back a problem from a space with a complex geometry to a
simpler one via a particular family of maps lies at the core of the area of
constrained optimization. Mirror descent [47] and Lagrange multipliers are
examples of methods that allow simplifying a constrained optimization problem
into one in a simpler space. It is also a recurrent theme in the literature on
optimization on manifolds. We have the Burer-Monteiro method [12, 13] which
parametrizes the space of low-rank matrices via the multiplication of two tall
matrices. Other approaches via different factorizations have also been
explored, for example via an SVD factorization [58] or other factorizations
like the polar factorization [46]. Similar parametrizations have been explored
for the algebraic variety of matrices of rank at most $k$—low-rank
optimization—and manifolds such as positive semidefinite matrices of a fixed
rank [59, 45]. The idea behind these algorithms is that of parametrizing a
space in terms of simpler ones, effectively pulling back the problem to an
easier one and then mapping back the solution of the simpler problem to the
complex space via the parametrization map.

##### Pullbacks along the exponential in numerical analysis.

Instances of (1) can be found all throughout the area of numerical analysis.
For example, in the study of the center of mass in the manifold of PD matrices
[3, 4]; continuous dynamics on Lie groups, and their discretizations [44, 33,
32]; optimization on compact Lie groups [42]; and probabilistic methods on
manifolds for their use in machine learning [23]. More recently, these ideas
were unified into the framework of _static trivializations_ [41].

The main idea behind pulling back an optimization problem from a manifold to a
flat space is that of heavily simplifying the problem at hand, going from a
manifold which might have a complex topology and geometry, to one that is flat
and has trivial topology. Of course, that comes at a cost, as this change of
topology makes the new function $f\circ\exp_{p}$ have more critical points.
Luckily, for any differentiable manifold and any smooth metric, this set of
problematic points, called the _conjugate locus_ , has measure zero by a
theorem of Sard [53]. Since the conjugate locus is a subset of the _cut locus_
, it is far from the point $p$. This is because the cut locus can vaguely be
regarded as _the set that is opposite to the point $p$_. For example, on a
sphere, the cut locus of a point $p$ is its antipode and on a cylinder, the
cut locus of a point is the whole line opposite to that point. As such, the
cut locus is, in some sense, _as far from $p$ as possible_, and so is the
conjugate locus. At the expense of this technicality, one may heavily simplify
the problem of optimizing a function on a space with a non-trivial topology to
a Euclidean problem. We will expand on this technical point in Section 2.

##### Pullbacks along the exponential map: Trivializations.

If we pullback the problem to the tangent space at each iteration of the
algorithm along a retraction, we recover Riemannian Gradient Descent (rgd)
along that retraction [11]. It was then noted in [41] that one may change the
point $p$ in (1) to avoid converging to a point in the conjugate locus of $p$.
In fact, one may change the point according to an arbitrary stopping rule
giving the following meta-algorithm:

  1. 1.

Compute gradient steps of the function $f\circ\exp_{p_{i}}$ to obtain iterates
$v_{i,k}\in T_{p_{i}}M$ for $k=1,2,\dots$

  2. 2.

If we are close to a point in the conjugate locus of $p_{i}$, we set
$p_{i+1}=\exp_{p_{i}}(v_{i,k})$ and repeat

This translates to Algorithm 1 given a boolean stopping rule stop and it is
called the _dynamic trivializations framework_ [41].

Algorithm 1 Dynamic trivialization framework

1:A starting point $p_{0}\in M$, a boolean statement stop

2:for $i=0,\ldots$ do

3: $v_{i,0}=0$

4: for $k=1,\dots$ do

5: $v_{i,k}=v_{i,k-1}-\eta_{i,k-1}\nabla\lparen
f\circ\exp_{p_{i}}\rparen\lparen v_{i,k-1}\rparen$$\triangleright$ Optimize on
$T_{p_{i}}M$

6: if stop then$\triangleright$ Stopping condition

7: break

8: $p_{i+1}=\exp_{p_{i}}(v_{i,k})$ $\triangleright$ Update pullback point

Two possible stopping rules stand out over the others. If we let
$\texttt{stop}\equiv\texttt{False}$, then we recover the _static-
trivialization framework_ of pulling back the function to a fixed tangent
space as in (1). On the other hand, if $\texttt{stop}\equiv\texttt{True}$,
then the algorithm is exactly Riemannian gradient descent, as mentioned
before. The stopping rule of being too close to the conjugate locus could be
then written as

$\texttt{stop}\equiv\frac{\lVert\nabla\lparen f\circ\exp_{p_{i}}\rparen\lparen
v_{i,k}\rparen\rVert}{\lVert\nabla f\lparen
x_{i,k}\rparen\rVert}<\varepsilon.$

The ideas of pulling back the problem at the current point of the optimization
and then taking a few steps on the tangent space have been recently explored
in [16, 56] in the context of escaping saddle points for optimization on
manifolds.

##### Retractions with bounded second and third order derivatives.

To be be able to prove convergence of the dynamic trivialization framework, we
will give novel bounds on the norm of the Hessian of the exponential map in
terms of the curvature of the manifold. These bounds are one example of second
order bounds for a large family of retractions. These bounds lie at the heart
of the convergence results of Riemannian methods via retractions, where it is
often referred to as the _$L$ -smoothness of the retraction_. Examples of work
under this assumptions are adaptive methods on matrix manifolds [37], quasi-
Newton methods on manifolds [31], stochastic methods with variance reduction
[38], and general convergence of Newton methods in Riemannian manifolds [25],
among many others. The bounds developed in this paper give exact rates of
growth for the norm of the Hessian of the retraction in terms of the curvature
of the manifold when the retraction is the exponential map, which falls
exactly in the framework of these papers and many others.

Third order and higher order bounds are used when trying to converge to second
order optima and escaping saddle points [16] and when adapting Newton methods
to manifolds [1], and other higher order methods. The computations of higher
order bounds is analogous to that of second order bounds. In this paper, we
present all the necessary techniques to give $n$-th order bounds for a
manifold of $n$-bounded geometry (see Definition 4.3 for a definition),
although we do not perform these computations explicitly. Bounds of this
flavour were first developed in [19], where the rate of growth of the $n$-th
derivative of the Christoffel symbols in normal coordinates is estimated.

##### General assumptions.

In this paper, we look at the problem of proving convergence for the family of
dynamic trivialization methods all at once. To do so, it is enough to control
the iterates of the (Euclidean) gradient descent in the context of problem
(1). If the bounds do not have a dependency on $p$, then they will provide
guarantees for Algorithm 1. This is the road that we will take in this paper.

For some manifolds, such as the hyperbolic space or the space of symmetric
positive definite matrices (and more generally on Hadamard spaces), one may
find (geodesically) convex functions which can be optimized efficiently [5].
On the other hand, when the manifold is compact, as is the case of the sphere,
or when dealing with orthogonality constraints [18], there exists no convex
function besides the constant ones, so we inevitably have to deal with non-
convex problems.

To this end, we will consider a function $f$ which is $\alpha$-weakly convex
(i.e., with bounded Hessian in operator norm) and we ask ourselves whether
$f\circ\exp_{p}$ is $\alpha^{\prime}$-weakly convex for any $\alpha^{\prime}$.
The answer to this problem will be: Not in general. For example, for the
hyperbolic space, the exponential map grows as $\cosh(t)$, which has unbounded
first and second derivatives. For this reason, we will have to content
ourselves with proving tight bounds on the growth of the Hessian of
$f\circ\exp_{p}$ on a neighbourhood of a given radius.

When it comes to the assumptions on the manifold, we will look at the family
of manifolds of bounded first order geometry.

######  Definition 1.1 (Bounded first order geometry).

We say that a Riemannian manifold $(M,\textsl{g})$ has
$(\delta,\Delta,\Lambda)$-bounded geometry if it has uniformly bounded
injectivity radius $\operatorname{r_{inj}}>0$ and we have bounds on the
sectional curvature and the derivative of the curvature tensor $R$

$\delta\leq\sec\leq\Delta\qquad\lVert\nabla R\rVert<\Lambda.$

These manifolds have found many uses in the area of PDEs on manifolds and
geometric analysis. In particular, any compact manifold with any metric is of
bounded geometry. On the other hand, this family does not impose topological
restrictions on $M$ as any differentiable manifold admits a complete metric of
bounded geometry [27].

For some general introduction to this family of manifolds see the thesis [20].

###  1.1 Contributions

In the paper, we will implement the following strategy:

  1. 1.

Give first and second order bounds for $\exp_{p}$ that depend on the curvature
of the manifold

  2. 2.

Use these bounds to estimate the weak convexity of $f\circ\exp_{p}$ on a
bounded domain

  3. 3.

Use these weak convexity bounds to prove the convergence of the dynamic
trivialization framework

First order bounds for the Riemannian exponential are well known in geometry
under the name of _Rauch’s theorem_. We will present a self-contained proof of
them here in full generality in Section 3 as it is not easily found in the
literature. We also take the chance to prove a law of cosines in manifolds
that follows from this in Section 3.4. Particular cases of this law of cosines
have been used throughout the optimization literature, but it is not easy to
find a self-contained proof of the general case.

Second order bounds for the exponential are considerably more difficult to
prove. There have been some recent attempts at bounding these quantities in
the context of optimization on manifolds [56]. In particular, concurrent to
this work, some new bounds were given in [17]. The bounds that we present in
this work are tighter than those in these papers. Historically, Kaul was the
first to prove second order bounds in [39]. We borrow some ideas from this
paper, while greatly improving their bounds. In particular, these original
bounds were given in terms of a solution of a differential equation that
cannot be integrated. The bounds that we give in Section 4 are explicit in
terms of trigonometric functions.

######  Theorem 1.2 (Bounds on the Full Hessian).

Let $(M,\textsl{g})$ be a Riemannian manifold with
$(\delta,\Delta,\Lambda)$-bounded geometry. For a geodesic
$\gamma\colon[0,r]\to M$ with initial unit vector $v$,
$r<\pi_{\frac{\Delta+\delta}{2}}$, and any two vectors $w_{1},w_{2}\in
T_{p}M$, we have that

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w_{1},w_{2}\rparen\rVert\leq\mbox{{}\small$\displaystyle\frac{8}{3r^{2}}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}\lparen\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}\lparen
r\rparen\rparen\lVert w_{1}\rVert\lVert w_{2}\rVert.$

where

$\operatorname{sn}_{\delta}(t)\coloneqq\begin{cases}\frac{\sin(\sqrt{\delta}t)}{\sqrt{\delta}}\qquad&\text{if
}\delta>0\\\ t\qquad&\text{if }\delta=0\\\
\frac{\sinh(\sqrt{-\delta}t)}{\sqrt{-\delta}}\qquad&\text{if
}\delta<0\end{cases}$

and $\pi_{\kappa}\in(0,\infty]$ is the first positive zero of
$\operatorname{sn}_{\kappa}$. Furthermore, the radius
$\pi_{\frac{\Delta+\delta}{2}}$ is tight for $\operatorname{SO}\lparen
n\rparen$ with a bi-invariant metric.

We also give tighter bounds in Theorem 4.10, although we believe that this
more compact bound will be easier to handle in general applications.

We compute explicit bounds for some manifolds which find uses in optimization
in Section 4.6. For example, for the special orthogonal embedded in
$\mathbb{R}^{n\times n}$ group we get

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}(w_{1},w_{2})\rVert\leq\frac{r}{3}\lVert
w_{1}\rVert\lVert w_{2}\rVert\qquad{r<2\sqrt{2}\pi}$

and for the real Grassmannian with the canonical metric

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}(w_{1},w_{2})\rVert\leq\frac{8r}{3}\lVert
w_{1}\rVert\lVert w_{2}\rVert\qquad{r<\pi}.$

Once we have first and second order bounds, in Section 5, we bound the
differential and the Hessian of the exponential map using the chain rule and
Cauchy-Schwartz

$\lVert\operatorname{Hess}(f\circ\exp_{p})\rVert\leq\lVert\operatorname{Hess}f\rVert\lVert\mathrm{d}\exp_{p}\rVert^{2}+\lVert\mathrm{d}f\rVert\lVert\operatorname{Hess}\exp_{p}\rVert$

getting the following result.

######  Theorem 1.3 (Weak convexity of the pullback).

Let $(M,\textsl{g})$ be a connected and complete Riemannian manifold and let
$f$ be an $\alpha$-weakly convex function on it. Fix a point $p\in M$, and a
set $\overline{\mathcal{X}}\subseteq T_{p}M$ with
$\operatorname{diam}(\overline{\mathcal{X}})=r<\pi_{\frac{\Delta+\delta}{2}}$
that does not intersect the conjugate locus and has at least one critical
point of $f$ in it. Then the map $f\circ\exp_{p}$ is
$\widehat{\alpha}_{r}$-weakly convex with constant

$\widehat{\alpha}_{r}=\alpha(C_{1,r}+C_{2,r})$

for

$\displaystyle
C_{1,r}=\max\Big{\\{}1,\mbox{{}\small$\displaystyle\frac{\operatorname{sn}_{\delta}(r)^{2}}{r^{2}}$}\Big{\\}}$
$\displaystyle
C_{2,r}=\mbox{{}\small$\displaystyle\frac{8}{3}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}\lparen\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}\lparen
r\rparen\rparen.$

The quantity $\widehat{\alpha}_{r}$ depends on the size of $r$, so the
distortion produced by precomposing with the exponential map may not be
uniformly bounded if the manifold has negative curvature and is not compact—as
is the case with the hyperbolic space. On the other hand, we can give exact
bounds for this deformation so, in practice, we can assume an upper bound on
$r$ and then tune it as necessary. This approximation is common when proving
convergence in manifolds that depend on the curvature (see for example [9, 61,
54, 57, 2]).

Finally, we showcase how to use these bounds to prove convergence of static
and dynamic trivializations. This is a corollary from the previous results and
the convergence of gradient descent in $\mathbb{R}^{n}$.

##  2 Differential geometry: Conventions and notation

In this section, we stablish the notation used to refer to some recurrent
objects, such as the distance function, the segment domain, mixed derivatives
and pullbacks of connections. We also use it to recall some definitions from
differential geometry that are not that commonly seen in the area of
optimization, such as the definition of the Lie derivative and covariant
derivative of tensors.

##### The manifold

We will always work on a connected and complete Riemannian manifold
$(M,\textsl{g})$ of dimension $m\geq 2$ and at least $C^{4}$ regularity. We
will denote the sectional curvature on the plane defined by two vectors
$u,v\in T_{p}M$ as $\sec(u,v)$. When we write $\sec\leq\Delta$ for a constant
$\Delta\in\mathbb{R}$, we mean that for every $p\in M$, $u,v\in T_{p}M$,
$\sec(u,v)\leq\Delta$. By a _geodesic_ , we will always mean a unit speed
geodesic. We will implicitly identify $T_{\lparen p,v\rparen}(T_{p}M)\cong
T_{p}M$ for every $v\in T_{p}M$.

##### Distance function, segment domain, conjugate and cut locus

For a point $p\in M$, we define the _segment domain_ as

$\operatorname{seg}(p)\coloneqq\\{v\in T_{p}M\>|\>\exp_{p}(tv)\text{ is length
minimizing for }t\in[0,1]\\}$

and we denote its interior by $\overline{U}_{p}$. The set $\overline{U}_{p}$
is a star-shaped open neighbourhood of $0$ in $T_{p}M$. On this neighbourhood
the exponential map is a diffeomorphism. In particular, there exists an
inverse $\exp_{p}^{-1}$ and $\mathrm{d}\exp_{p}$ is full rank on
$\overline{U}_{p}$. We will also write
$U_{p}\coloneqq\exp_{p}(\overline{U}_{p})$. $U_{p}$ is sometimes referred as a
_normal neighborhood of $p$_, and its complement is called the _cut locus_
$\operatorname{cut}(p)\subseteq M$.

A conjugate point $\exp_{p}(v)=q$ of $p$ is one such at which
$\lparen\mathrm{d}\exp_{p}\rparen_{v}$ is not full rank. We denote the set of
all these points $\operatorname{conj}(p)$. We have that
$\operatorname{conj}(p)\subseteq\operatorname{cut}(p)$. By a theorem of Sard,
$\operatorname{conj}(p)$ has measure zero in $M$ [53]. Even more,
$\operatorname{cut}(p)$ has Hausdorff dimension at most $n-1$ [34]. As such,
$U_{p}$ is an open neighbourhood of $p$ that covers almost all the manifold.

We will write $r(x)=d(x,p)$ for the distance to a fixed point $p$. This
function is differentiable on $U_{p}\backslash\\{p\\}$, with gradient the unit
radial vector field emanating from $p$. In particular, we have that for a
geodesic $\gamma$ starting at $p$, $\nabla r|_{\gamma}=\dot{\gamma}$. The cut
locus of $p$ is exactly the set of points other than $p$ at which the distance
function $r$ is not differentiable. This will be import

##### Einstein convention

Whenever we refer to coordinates $\\{x^{i}\\}$ in $U_{p}$, we will always
assume that these are the normal coordinates given by $\exp_{p}^{-1}$ and
certain fixed frame on $T_{p}M$. We will denote
$\partial_{i}\coloneqq\frac{\partial}{\partial x^{i}}$ for short. We will use
Einstein’s summation convention that if an index appears as a super-index and
a sub-index in a formula, it means that we are summing over it. For example,
for a vector field $X$ in local coordinates we write

$X=\sum_{i=1}^{n}X^{i}\partial_{i}=X^{i}\partial_{i}.$

######  Remark.

Almost all the results in this paper can be developed working just on a
neighbourhood of a geodesic. As such, most of the times it will not be
necessary to work on a subset of $U_{p}$ but merely on a neighbourhood of a
geodesic at which on which the exponential is a diffeomorphism. We will make
sure of make this explicit in each result throughout the paper.

##### Connections and derivations

We will write $D_{X}f\coloneqq\mathrm{d}f(X)$ for the directional derivative
of a function along a vector field $X$ to disambiguate with the gradient of a
function, which we will denote by $\nabla f$ its gradient. We denote by
$\nabla$ the Levi-Civita connection on $U_{p}$ and by $\nabla^{\text{flat}}$
the pushforward of the flat connection on $T_{p}M$ along $\exp_{p}$. In the
same way that we do for the norms, we will abuse the notation and also denote
by $\nabla$ the associated connections defined by $\nabla$ in the associated
bundles. We recall that connections on tensor bundles are defined so that the
Leibnitz rule holds. For example, for the Hessian and three vector fields
$X,Y_{1},Y_{2}$,

$D_{X}\lparen\operatorname{Hess}f\lparen
Y_{1},Y_{2}\rparen\rparen=\lparen\nabla_{X}\operatorname{Hess}f\rparen\lparen
Y_{1},Y_{2}\rparen+\operatorname{Hess}f\lparen\nabla_{X}Y_{1},Y_{2}\rparen+\operatorname{Hess}f\lparen
Y_{1},\nabla_{X}Y_{2}\rparen$

so we define

$\lparen\nabla_{X}\operatorname{Hess}f\rparen\lparen
Y_{1},Y_{2}\rparen\coloneqq D_{X}\lparen\operatorname{Hess}f\lparen
Y_{1},Y_{2}\rparen\rparen-\operatorname{Hess}f\lparen\nabla_{X}Y_{1},Y_{2}\rparen-\operatorname{Hess}f\lparen
Y_{1},\nabla_{X}Y_{2}\rparen.$

We can do the same for a tensor $T$ of type $(0,s)$. We define its covariant
derivative $\nabla T$ as the tensor of type $(0,s+1)$ such that the Leibnitz
rule holds

$\lparen\nabla_{X}T\rparen(X_{1},\dots,X_{s})\coloneqq D_{X}\lparen T\lparen
X_{1},\dots X_{s}\rparen\rparen-\sum_{i=1}^{s}T\lparen
X_{1},\dots,\nabla_{X}X_{i},\dots,X_{s}\rparen.$

For example, if we have in local coordinates the differential of a function is
a $(0,1)$ tensor $\mathrm{d}f=f^{i}\mathrm{d}x_{i}$. We may compute its
Hessian in local coordinates as the $(0,2)$ tensor given by

$\nabla\mathrm{d}f=\nabla(f^{i}\mathrm{d}x_{i})=\mathrm{d}f^{i}\otimes\mathrm{d}x_{i}+f^{i}\nabla\mathrm{d}x_{i}$

where we have used that the connection on functions is just the differential,
by definition.

The Lie derivative of tensors is defined in the same way. For a $(0,s)$ tensor
$T$, the Lie derivative of $T$, $\mathcal{L}_{X}T$ in the direction of a
tensor $X$ is defined as the $(0,s)$ tensor

$\lparen\mathcal{L}_{X}T\rparen(X_{1},\dots,X_{s})\coloneqq D_{X}\lparen
T\lparen X_{1},\dots X_{s}\rparen\rparen-\sum_{i=1}^{s}T\lparen
X_{1},\dots,\mathcal{L}_{X}X_{i},\dots,X_{s}\rparen.$

##### Pullback connections

When dealing with a smooth curve $\gamma\colon[0,r]\to M$, we will sometimes
want to derive a vector field $X$ along it. We will write
$\nabla_{\partial_{t}}X$ for the _covariant derivative along $\gamma$_, that
is, the _pullback connection along $\gamma$_.

We will sometimes write the covariant derivative along $\gamma$ as
$\dot{X}\coloneqq\nabla_{\partial_{t}}X$. This comes from the fact that if we
choose a parallel frame $\\{e_{i}\\}$ along $\gamma$, that is
$\nabla_{\partial_{t}}e_{i}=0$, we have that by the Leibnitz rule

$\dot{X}=\nabla_{\partial_{t}}\lparen
X^{i}e_{i}\rparen=\nabla_{\partial_{t}}\lparen X^{i}\rparen
e_{i}+X^{i}\nabla_{\partial_{t}}\lparen e_{i}\rparen=\dot{X}^{i}e_{i}$

where $\dot{X}^{i}$ are just the regular derivatives of the coordinate
functions. With this notation, the usual equation for geodesics simply reads
$\ddot{\gamma}=0$.

If instead of a curve we have an embedded surface

$\\!\begin{aligned} c\colon[0,r]\times[-\varepsilon,\varepsilon]&\to M\\\
(t,s)&\mapsto c(t,s)\end{aligned}$

we will analogously write $\partial_{s}$ for the vector field in the direction
of the second component. Suppose that we have a vector field $J$ along
$\gamma(t)\coloneqq c(t,0)$ given by

$J(t)=\frac{\partial c}{\partial s}(t,0)=\mathrm{d}c(\partial_{s})|_{\gamma}$

where $\partial_{s}$ is the coordinate vector field on the second component of
$[0,r]\times[-\varepsilon,\varepsilon]$. We will abuse the notation and write
for a vector field $X$ along $\gamma$

$\nabla_{J}X\coloneqq\nabla_{\partial_{s}}X$

in the same way that we may write the equation for the geodesics as

$\nabla_{\dot{\gamma}}\dot{\gamma}\coloneqq\nabla_{\partial_{t}}\dot{\gamma}=0.$

Sometimes, we will also pullback a connection along the exponential map
$\exp_{p}$. For the distance function $r$ to $p$, we have the radial vector
field on $U_{p}\backslash\\{p\\}\subseteq M$ given by $\nabla r$. By Gauss’s
lemma, the pullback of this vector field to
$\overline{U}_{p}\backslash\\{0\\}\subseteq T_{p}M$ via the exponential map is
exactly the radial vector field on the tangent space $\partial_{r}$, that is

$\mathrm{d}\exp_{p}(\partial_{r})=\nabla r.$

With this notation, we can write the equation for the geodesics on the whole
$U_{p}\backslash\\{p\\}$ using the pullback connection along $\exp_{p}$ as

$\nabla_{\partial_{r}}\nabla r=0.$

##### Hessian and iterated Hessian

If we want to evaluate the $(0,2)$-Hessian twice on the same vector, we may do
so by simplifying this problem to just performing a Euclidean derivative. Let
$\gamma$ be the geodesic such that $\dot{\gamma}(0)=X$

$\lparen
f\circ\gamma\rparen^{\prime\prime}(0)=D_{\dot{\gamma}(0)}\langle\nabla
f,\dot{\gamma}\rangle=\langle\nabla_{\dot{\gamma}(0)}\nabla
f,\dot{\gamma}(0)\rangle=\nabla\mathrm{d}f(X,X)=\operatorname{Hess}f(X,X).$
(2)

We will use the notation $\nabla\mathrm{d}f$ to refer to the Hessian, given
that it generalizes to smooth maps between manifolds such as the exponential
map.

We also recall the definition of the iterated Hessian of a function as a
$(0,2)$ tensor. This is also best introduced in its $(1,1)$ form

$\operatorname{Hess}^{2}f(X)\coloneqq\lparen\operatorname{Hess}f\circ\operatorname{Hess}f\rparen(X)=\nabla_{\nabla_{X}\nabla
f}\nabla f.$

Its $(0,2)$ version is consequently defined as

$\operatorname{Hess}^{2}f(X,Y)\coloneqq\langle\operatorname{Hess}f(\operatorname{Hess}f(X)),Y\rangle=\operatorname{Hess}f(\operatorname{Hess}f(X),Y)=\langle\operatorname{Hess}f(X),\operatorname{Hess}f(Y)\rangle$
(3)

where we have used that the Hessian is symmetric. Through this formula it is
clear that $\operatorname{Hess}^{2}f$ is also symmetric.

##  3 First Order Bounds for the Exponential Map

In this section, we give bounds on the norm of the differential of the
exponential map of a manifold with bounded sectional curvature. In particular,
if the sectional curvature of any plane of $M$ is bounded above and below by
$\delta\leq\sec\leq\Delta$, for a ball $B_{p}(r)$, the bounds will be of the
form

$f_{\Delta}(r)\lVert
w\rVert\leq\lVert\lparen\mathrm{d}\exp_{p}\rparen_{v}(w)\rVert\leq
f_{\delta}(r)\lVert w\rVert\mathrlap{\qquad\forall w\in T_{p}M,\lVert
v\rVert\leq r}$

for $r\geq 0$ and suitable functions
$f_{\delta},f_{\Delta}\colon[0,r]\to\mathbb{R}^{+}$.

First order bounds for the exponential map have been known since the times of
Cartan [14], and are well known in areas such as comparison geometry or PDEs,
but they are not so known in the area of optimization on manifolds. The tight
bounds that we will present here are often referred to as _Rauch’s theorem_ as
he was the first to show these bounds in the positive curvature-case [51] and
later in the general case [52].

There are quite a few proof techniques of these theorems. There exist proofs
via the second variation formula of the energy [55], through bounds on the
norm via the Jacobi equation [35], by bounding Riccati equation on self-
adjoint operators and integrating the result [21], or by the study of the
distance function to a given point [28]. We choose to present here this fourth
approach.

The approach presented here was first introduced by Gromov in the context of
volume bounds, in what’s called now the Bishop–Gromov theorem. This approach
has the advantage of yielding upper and lower bounds for positive and negative
curvatures at the same time. Other approaches need different techniques for
the upper and lower bounds or they just give some weaker version of the
theorem with constraints on the sign of the curvature. This approach also has
a more geometric flavour as it accounts for bounding the principal curvatures
of geodesic balls on the manifold. The general strategy of this proof is to
decompose certain homogeneous second order differential equation into a
Riccati equation and a system of first order equations, bound the Riccati
equation, and then integrate the result to get the bounds on the differential
of the exponential. This approximation to Jacobi fields and parallel vector
fields is can be found in [15], and the approximation using curvature
equations is from [50]. Some of the proofs are either corrections or
simplifications of the original ones.

###  3.1 Jacobi fields

We begin by introducing what will be the main tool that we will just
throughout the paper: Jacobi fields. Jacobi fields are certain vector fields
along a given geodesic that describe the behavior of the differential of the
exponential map along this geodesic. We will see that they are the solution of
a certain differential equation, and we will use this equation to give the
bounds on the differential of the exponential.

Before defining what Jacobi fields are and deducing this differential equation
we will need a lemma. This lemma that can loosely be interpreted as _the Lie
derivative commutes with the differential of a smooth map_.

######  Lemma 3.1.

For a diffeomorphism $f$ and two vector fields $X,Y$ on $M$ we have that

$\mathrm{d}f([X,Y])=[\mathrm{d}f(X),\mathrm{d}f(Y)].$

###### Proof.

See, for example, [[]1\. Lemma 22]oneill1966fundamental. ∎

With this lemma in hand, we are ready to show that the differential of the
exponential satisfies certain differential equation.

######  Proposition 3.2 (First order Jacobi equation).

Let $(M,\textsl{g})$ be a complete Riemannian manifold. Let
$\gamma\colon[0,r]\to M$ be a geodesic in $U_{p}$ with initial unit vector
$v\coloneqq\dot{\gamma}(0)$ and a vector $w\in T_{p}M$ such that $w\perp v$.
Consider the following variation of $\gamma$ in the direction of $w$

$\\!\begin{aligned} c\colon[0,r]\times(-\varepsilon,\varepsilon)&\to M\\\
(t,s)&\mapsto\exp_{p}(t(v+sw))\end{aligned}$

where $\varepsilon$ is small enough so that $c$ defines an embedded surface.

Define the vector field $J$ along $\gamma$ as

$J(t)\coloneqq\mathrm{d}c(\partial_{s})|_{\gamma}=\frac{\partial c}{\partial
s}(t,0)=\lparen\mathrm{d}\exp_{p}\rparen_{tv}(tw).$

This vector field solves the following differential equation on $\gamma$

$\mathcal{L}_{\nabla r}J|_{\gamma}=0$

or equivalently, using the pullback connection

$\nabla_{\dot{\gamma}}J=\nabla_{J}\nabla r.$ (4)

###### Proof.

The fact that $\gamma$ is in $U_{p}$ is equivalent to saying that
$tv\in\overline{U}_{p}$ for $t\in[0,r]$. From this we get that the
$\varepsilon$ in the definition exists as $\overline{U}_{p}$ is open. Since
$c$ defines an embedded surface, $c$ is a diffeomorphism onto its image.

Since $\partial_{t},\partial_{s}$ are coordinate vector-fields of the surface
defined by $c$, we have that

$[\partial_{t},\partial_{s}]=0.$

Using that $c$ is a diffeomorphism, by Lemma 3.1 together with the fact that
the Lie derivative restricted to a submanifold is the Lie derivative of the
restrictions, we have that

$[\mathrm{d}c\lparen\partial_{t}\rparen,\mathrm{d}c\lparen\partial_{s}\rparen]=\mathrm{d}c\lparen[\partial_{t},\partial_{s}]\rparen=0.$

We get the Lie formulation of the Jacobi equation by noting that
$\mathrm{d}c(\partial_{t})|_{\gamma}=\nabla r|_{\gamma}$. The equation using
the pullback connection is just a reformulation of the Lie derivative using
that the Levi-Civita connection is torsion-free. ∎

######  Remark.

Note that, even though the differential equation involves a Lie bracket, we
just need to define $J$ along the flow of $\nabla r$, which is just a geodesic
$\gamma$. This equation can be extended into a PDE on the whole
$U_{p}\backslash\\{p\\}$, looking for vector fields on this domain such that

$\mathcal{L}_{\nabla r}W=0.$

One such a solution is clearly given by the radial vector field $W=\nabla r$.

To find other solutions, we imitate the construction that we did in
Proposition 3.2 and define $W$ on a sphere around $p$ and extend it radially
to the whole $U_{p}$ as

$W(x)=\mathrm{d}\exp_{p}\Big{\lparen}r(x)W\Big{\lparen}\mbox{{}\small$\displaystyle\frac{x}{r\lparen
x\rparen}$}\Big{\rparen}\Big{\rparen}.$

It is clear by Proposition 3.2 that this vector field solves the first order
Jacobi equation on the whole $U_{p}\backslash\\{p\\}$.

We will now use this equation to derive the linear version of the Jacobi
equation.

######  Proposition 3.3 (Jacobi equation).

Let $\gamma\colon[0,r]\to M$ be a geodesic with initial values $\gamma(0)=p$,
$\dot{\gamma}(0)=v$. Then $J(t)=\lparen\mathrm{d}\exp_{p}\rparen_{tv}(tw)$ is
a vector field along $\gamma$ that solves the following second order
homogeneous linear differential equation

$\ddot{J}+R(J,\dot{\gamma})\dot{\gamma}=0.$

We call this equation the _Jacobi equation_.

###### Proof.

If $w$ is parallel to $v$, $J$ takes the form $J(t)=t\lVert
w\rVert\dot{\gamma}(t)$ or $J(t)=-t\lVert\dot{\gamma}\rVert(t)$. These fields
solve the differential equation, since
$R(\dot{\gamma},\dot{\gamma})\dot{\gamma}=0$.

If $w$ is normal to $v$, we may differentiate (4) to get

$\nabla_{\partial_{t}}\nabla_{\partial_{t}}J=\nabla_{\partial_{t}}\nabla_{J}\nabla
r$

and using the definition of the curvature tensor and the fact that $[J,\nabla
r]=0$ we get

$\nabla_{\partial_{t}}\nabla_{\partial_{t}}J+R(J,\nabla r)\nabla r=0.$

We then get the Jacobi equation by restricting this equation to $\gamma$,
since $\nabla r\circ\gamma=\dot{\gamma}$.

Using an orthonormal parallel frame $\\{e_{i}\\}$ along $\gamma$ such that
$e_{1}=\dot{\gamma}$, and expressing a solution of this equation as
$J^{i}e_{i}$, where $J^{i}\colon[0,r]\to\mathbb{R}$, we see that, in fact,
this is a second order linear differential equation on $[0,r]$

$\ddot{J}^{i}+R_{j}^{i}J^{j}=0\mathrlap{\qquad i=1,\dots,n}$

with coefficients

$R_{j}^{i}=\begin{cases}0\qquad&\text{if }i=1\\\
R(e_{j},\dot{\gamma},\dot{\gamma},e_{i})\qquad&\text{if
}i=2,\dots,n.\end{cases}\qed$

######  Definition 3.4 (Jacobi field).

We say that a vector field $J$ along a geodesic $\gamma$ is a Jacobi field if
it satisfies the Jacobi equation.

######  Remark (Basic properties of Jacobi fields).

By elementary theory of differential equations, the Jacobi equation has $2n$
independent solutions, defined by the initial values $(J(0),\dot{J}(0))\in
T_{p}M\times T_{p}M$. By construction, we have found $n$ independent solutions
given by

$J(t)=\lparen\mathrm{d}\exp_{p}\rparen_{tv}(tw)\mathrlap{\qquad\forall w\in
T_{p}M.}$

These correspond to the initial values $J(0)=0$, $\dot{J}(0)=w$.
Geometrically, they correspond to a family of geodesics $c$ that fixes the
initial point, that is, $c(0,s)=p$.

It is also direct to show that $\dot{\gamma}$ is another solution to this
equation with initial values $J(0)=v$, $\dot{J}(0)=0$. The other $n-1$
independent solutions correspond to variations of the geodesic $\gamma$ that
do not fix the initial point $p$. These solutions will not be important for
our analysis.

It starts being clear now the importance of Jacobi fields. If we want to
control the norm of the differential of $\exp_{p}$ at a point
$rv\in\overline{U}_{p}$, for $r>0$, $\lVert v\rVert=1$, we may define the
Jacobi field along the geodesic $\gamma\colon[0,r]\to M$ with initial
conditions $(p,v)$ and we have that

$\lparen\mathrm{d}\exp_{p}\rparen_{rv}(w)=\frac{J(r)}{r}.$

In order to bound the norm of the differential of the exponential, we just
need to bound the norm of the solutions of the Jacobi equation.

If the vector $w$ is parallel to $v$, as we saw in the proof of Proposition
3.3, the Jacobi field takes the form $J(t)=t\lVert w\rVert\dot{\gamma}(t)$. As
such, in this direction we have an exact solution of the Jacobi equation and
we can compute its norm exactly as $\lVert J(t)\rVert=t\lVert w\rVert$.

If the vector $w$ is perpendicular to $\dot{\gamma}(0)$, then $J(t)$ is
perpendicular to $\dot{\gamma}(t)$ for every $t\in[0,r]$, as the equation for
$J^{1}$ would be given by

$\ddot{J}^{1}=0\mathrlap{\qquad(J^{1}(0),\dot{J}^{1}(0))=(0,0).}$

Geometrically, the previous proposition says that Jacobi fields just rotate
around the vector field defined by $\dot{\gamma}$. In symbols, this means that
if we split a Jacobi field along $\gamma$ into its radial and normal part as

$J^{{\dagger}}=\langle J,\dot{\gamma}\rangle\dot{\gamma}\qquad
J^{\perp}=J-J^{{\dagger}}$

if $\dot{J}^{{\dagger}}(0)=0$ then $J^{{\dagger}}(t)=0$ for every $t\in[0,r]$.

The results discussed in this remark are commonly known as the _Gauss’s
lemma_.

Jacobi fields are also very closely related to conjugate points and the
conjugate locus. A point $q=\exp_{p}(v)$ is conjugate to $p$ if the
exponential at $v$ is not full rank. Another way of looking at this definition
is by defining a point $q$ conjugate to $p$ if there exists a Jacobi field
connecting $p$ and $q$ such that it is zero at $p$ and $q$. We will see when
studying Rauch’s theorem that the estimates for Jacobi fields will work on
every point on the manifold but the conjugate locus, as these points will be
exactly the singularities of the maps that we work with. For a unit vector
$v$, we will denote by $r_{\operatorname{conj}}(v)\in(0,\infty]$ the smallest
number at which $q=\exp_{p}(r_{\operatorname{conj}}v)$ is conjugate to $p$.

###  3.2 Parallel vector fields

In the same way that Jacobi fields with initial condition $J(0)=0$ and
$\dot{J}(0)\perp\dot{\gamma}(0)$ are the vector fields such that
$\mathcal{L}_{\nabla r}J=0$, we have parallel vector fields.

######  Definition 3.5.

We say that $E$ is a parallel vector field on $U_{p}\backslash\\{p\\}$ if it
satisfies

$\nabla_{\nabla r}E=0.$

Parallel vector fields can be easily described along a geodesic as parallel
transporting a vector $w\in T_{\gamma(0)}M$ along it. Analogously, we can also
define a parallel vector field on the whole $U_{p}\backslash\\{p\\}$ by
specifying a vector field on a sphere $S_{\varepsilon}=\\{x\in
M\>|\>d(x,p)=\varepsilon\\}\subseteq U_{p}$ and extending it to the whole
$U_{p}$ parallel transporting this vector field along geodesics.

###  3.3 Rauch’s theorem

We now go back to the study of the norm of the differential of the
exponential. The strategy that we will follow will be that of bounding the
derivative of the norm of the exponential. For that end, let $J$ be a Jacobi
field along a geodesic $\gamma$ such that $J(0)=0$,
$\dot{J}(0)\perp\dot{\gamma}(0)$. As we already saw before, this vector field
takes the form

$J(t)=\lparen\mathrm{d}\exp_{p}\rparen_{t\dot{\gamma}(0)}\lparen
t\dot{J}(0)\rparen.$

The derivative of its norm is given by

$\frac{\mathrm{d}}{\mathrm{d}t}\lVert
J\rVert=\frac{\langle\dot{J},J\rangle}{\lVert J\rVert}.$

At first sight it looks like we have not achieved much, as we still have a
term involving the square of the norm of $J$ on the right hand side, but it
turns out that it can be rewritten in a convenient way using (4),

$\langle\dot{J},J\rangle=\langle\nabla_{\nabla
r}J,J\rangle=\langle\nabla_{J}\nabla r,J\rangle=\operatorname{Hess}r(J,J).$

Our plan will be to bound the Hessian of the distance to $p$ on vectors of
norm $1$ to get bounds on the log-derivative of the norm of the exponential
map and then integrate these to get bounds on $\lVert J\rVert$. To this end,
we start by computing formulas that relate Hessian of the distance function to
the curvature tensor.

A geometric interpretation of this approach comes after noting that the
Hessian of the distance function is exactly the second fundamental form (or
shape operator) of the distance function. As such, this quantity can be
interpreted geometrically as the variation of the curvature of geodesic
spheres $S_{t}(p)=\\{x\in M\>|\>d(x,p)=t\\}$ in the radial direction.

######  Proposition 3.6 (Radial Curvature Equations).

Let $(M,\textsl{g})$ be a Riemannian manifold. For a point $p\in M$, we have
on $U_{p}$

$\displaystyle\lparen\mathcal{L}_{\nabla r}\operatorname{Hess}r\rparen(X,Y)$
$\displaystyle-\operatorname{Hess}^{2}r(X,Y)=-R(X,\nabla r,\nabla r,Y)$ (5)
$\displaystyle\lparen\nabla_{\nabla r}\operatorname{Hess}r\rparen(X,Y)$
$\displaystyle+\operatorname{Hess}^{2}r(X,Y)=-R(X,\nabla r,\nabla r,Y).$ (6)

where $\operatorname{Hess}^{2}$ denotes the _iterated Hessian_ defined in (3).

###### Proof.

Since for a geodesic $\gamma$ we have that $\nabla r|_{\gamma}=\dot{\gamma}$,
the gradient of $r$ has $\nabla_{\nabla r}\nabla r=0$. From this we get

Equation (5) follows after expanding the first term as

$D_{\nabla
r}\lparen\operatorname{Hess}r(X,Y)\rparen=\lparen\mathcal{L}_{\nabla
r}\operatorname{Hess}r\rparen(X,Y)+\operatorname{Hess}r\lparen\mathcal{L}_{\nabla
r}X,Y\rparen+\operatorname{Hess}r\lparen X,\nabla_{\nabla r}Y-\nabla_{Y}\nabla
r\rparen$

and equation (6) follows after expanding it as

$D_{\nabla r}\lparen\operatorname{Hess}r(X,Y)\rparen=\lparen\nabla_{\nabla
r}\operatorname{Hess}r\rparen(X,Y)+\operatorname{Hess}r\lparen\nabla_{\nabla
r}X,Y\rparen+\operatorname{Hess}r\lparen X,\nabla_{\nabla r}Y\rparen.\qed$

######  Remark (A Riccati-type equation).

Note that the previous proposition holds for every vector field $Y$. As such,
we may write it in its $(1,1)$ form

$\lparen\nabla_{\nabla
r}\operatorname{Hess}r\rparen(X)+\operatorname{Hess}^{2}r(X)+R(X,\nabla
r)\nabla r=0.$

Denoting by $R_{2}(X)\coloneqq R(X,\nabla r)\nabla r$ the $(1,1)$ curvature
tensor in this equation, and the shape operator for the geodesic balls as
$S\coloneqq\operatorname{Hess}r$, this equation can be rewritten in the radial
direction as a Riccati equation on symmetric $(1,1)$ tensors

$S^{\prime}+S^{2}+R_{2}=0.$

When restricted to a geodesic, this can be seen as a differential equation on
matrices, where $R_{2}$ is self-adjoint with respect to g. As one does in one
dimension splitting a second order differential equation into a first order
equation and a Riccati equation, here we have split the Jacobi equation into
first order matrix equation—more formally, an equation on symmetric
endomorphisms along the tangent spaces of $\gamma$—and a first order equation
for a parallel vector field $E$ of the form

$\dot{E}=S(E).$

A careful analysis of this Riccati equation yields another particularly clean
proof of Rauch’s theorem, at the expense of the use of more abstract methods,
as presented in [22].

In contrast, we will use parallel and Jacobi vector fields to simplify these
matrix equations.

######  Proposition 3.7.

Let $\gamma$ be a geodesic, and let $J$ be a Jacobi field along it such that
$J(0)=0$, $\dot{J}\perp\dot{\gamma}(0)$. Let $E$ be a parallel vector field
along $\gamma$ that is normal to $\dot{\gamma}$. We have the following
differential equations along $\gamma$

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\lparen\operatorname{Hess}r(J,J)\rparen$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\lparen\operatorname{Hess}r(E,E)\rparen$
$\displaystyle+\operatorname{Hess}^{2}r(E,E)=-\sec(E,\dot{\gamma}).$ (7)

###### Proof.

We just evaluate (5) and (6) on $J,E$ and use that

$\mathcal{L}_{\nabla r}J=0\qquad\nabla_{\nabla r}E=0.$

To go from the Riemannian tensor to the sectional curvature we just recall
that if a Jacobi field or a parallel vector field are normal to $\dot{\gamma}$
at one point they are normal to $\dot{\gamma}$ at every point.

Along a geodesic we have that the radial vector field $\nabla r$ is just
$\dot{\gamma}$, and as such we have that $D_{\nabla
r}|_{\gamma}=D_{\dot{\gamma}}=\frac{\mathrm{d}}{\mathrm{d}t}$. ∎

It is now a bit clearer how Jacobi fields may be used to give lower bounds and
parallel vector fields to give upper bounds on the Hessian of the distance
function.

We will now define some generalized trigonometric functions, which will be
useful in the sequel.

######  Definition 3.8 (Generalized trigonometric functions).

For a constant $\kappa\in\mathbb{R}$, we define the generalized sine function
$\operatorname{sn}_{\kappa}$ as the solution to the differential equation

$\ddot{x}+\kappa x=0\qquad x(0)=0,\,\dot{x}(0)=1.$

In particular, we have

$\operatorname{sn}_{\kappa}(t)\coloneqq\begin{cases}\frac{\sin(\sqrt{\kappa}t)}{\sqrt{\kappa}}\qquad&\text{if
}\kappa>0\\\ t\qquad&\text{if }\kappa=0\\\
\frac{\sinh(\sqrt{-\kappa}t)}{\sqrt{-\kappa}}\qquad&\text{if
}\kappa<0\end{cases}$

with its first positive zero being at

$\pi_{\kappa}\coloneqq\begin{cases}\frac{\pi}{\sqrt{\kappa}}&\qquad\text{if
}\kappa>0\\\ \infty&\qquad\text{if }\kappa\leq 0\end{cases}$

We define the generalized cotangent $\operatorname{ct}_{\kappa}$ as

$\operatorname{ct}_{\kappa}(t)\coloneqq\frac{\operatorname{sn}^{\prime}_{\kappa}(t)}{\operatorname{sn}_{\kappa}(t)}=\begin{cases}\sqrt{\kappa}\cot(\sqrt{\kappa}t)&\mathrlap{\qquad\text{if
}\kappa>0}\\\ \frac{1}{t}&\mathrlap{\qquad\text{if }\kappa=0}\\\
\sqrt{-\kappa}\coth(\sqrt{-\kappa}t)&\mathrlap{\qquad\text{if
}\kappa<0}\end{cases}$

This function is smooth on $(0,\pi_{\kappa})$.

######  Remark.

Note that the equation $\ddot{x}+\kappa x=0$ is the associated second order
equation to the Riccati equation $\dot{x}+x^{2}+\kappa=0$. This is exactly the
Riccati equation that the shape operator for the distance function solves in
constant curvature, as pointed out before. Given that
$\operatorname{sn}_{\kappa}$ solves the second order equation,
$\operatorname{ct}_{\kappa}$ solves the Riccati equation.

From an analytical point of view, it is worth noting that the family of
functions $\operatorname{sn}_{\kappa}(t)$ is strictly increasing in $\kappa$
for $t\in(0,\pi_{\kappa})$. We also have that $\operatorname{sn}_{\kappa}$ has
a zero at $0$ for every $\kappa$ and another one at $\pi_{\kappa}$ for
$\kappa>0$, so

$\displaystyle\lim_{t\to 0}\operatorname{ct}_{\kappa}(t)$
$\displaystyle=\infty\phantom{-}\mathrlap{\qquad\kappa\in\mathbb{R}}$
$\displaystyle\lim_{t\to\pi_{\kappa}}\operatorname{ct}_{\kappa}(t)$
$\displaystyle=-\infty\mathrlap{\qquad\kappa>0.}$

All we need to do now is to relate the solutions of the Jacobi equation in
Riccati as given by the Hessian of the distance function in Proposition 3.7
with its solutions in constant curvature. To do so, we need the following
elementary comparison lemma for functions of real variable. It can be thought
of as a version of Sturm’s comparison theorem for Riccati equations.

######  Lemma 3.9 (Comparison Lemma for Riccati Equations).

Let $\rho_{1},\rho_{2}\colon(a,b)\to\mathbb{R}$ be two differentiable
functions such that

$\dot{\rho}_{1}+\rho_{1}^{2}\leq\dot{\rho}_{2}+\rho_{2}^{2}.$

If we have that $\rho_{1}(t_{0})\geq\rho_{2}(t_{0})$ for a $t_{0}\in(a,b)$,
then $\rho_{1}\geq\rho_{2}$ on $(a,t_{0}]$.

###### Proof.

Let $u=\rho_{1}-\rho_{2}$ and
$s=\dot{\rho}_{2}-\dot{\rho}_{1}+\rho_{2}^{2}-\rho_{1}^{2}$. We have by
hypothesis that $u(t_{0})\geq 0$ and $s\geq 0$ on $(a,b)$, and it is enough to
prove that $u\geq 0$ on $(a,t_{0}]$. We can write the differential inequality
as a differential equation in terms of $u$ and $s$ as

$\dot{u}=-(\rho_{1}+\rho_{2})u-s.$

and we can solve this differential equation using the method of variation of
parameters. We compute the general solution of the homogeneous system

$\dot{v}=-(\rho_{1}+\rho_{2})v$

as $v=Cv_{1}$ for $v_{1}\coloneqq e^{-\int\rho_{1}+\rho_{2}}>0$ and a constant
$C\in\mathbb{R}$.

To get a particular solution to the inhomogeneous equation, we let $u=\eta
v_{1}$ and, plugging it into the inhomogeneous equation, we see that the
function $\eta$ satisfies $\dot{\eta}=-sv_{1}^{-1}$. As $s$ and $v_{1}^{-1}$
are positive, we conclude that $\eta$ is decreasing.

Finally, we get the general solution by adding the particular solution and the
general solution to the homogeneous system $u=v_{1}(C+\eta)$. Since, by
hypothesis, $u(t_{0})\geq 0$, we get that $C+\eta(t)\geq 0$ for
$t\in(a,t_{0}]$ so that $u(t)\geq 0$ on $(a,t_{0}]$. ∎

######  Proposition 3.10 (Riccati Comparison Estimate).

Let $\rho\colon(0,b)\to\mathbb{R}$ be a differentiable function and fix a
$\kappa\in\mathbb{R}$.

  1. 1.

If $\dot{\rho}+\rho^{2}\leq-\kappa$ then
$\rho(t)\leq\operatorname{ct}_{\kappa}(t)$ on $(0,b)$. Furthermore,
$b\leq\pi_{\kappa}$

  2. 2.

If $\dot{\rho}+\rho^{2}\geq-\kappa$ and $\lim\limits_{t\to 0}\rho(t)=\infty$
then $\operatorname{ct}_{\kappa}(t)\leq\rho(t)$ on
$(0,\min\\{b,\pi_{\kappa}\\})$.

###### Proof.

For the first case, assume that $\rho\lparen
t_{0}\rparen>\operatorname{ct}_{\kappa}(t_{0})$ for some $t_{0}\in(0,b)$. By
continuity, we can choose an $\varepsilon>0$ such that

$\rho\lparen t_{0}\rparen\geq\operatorname{ct}_{\kappa}(t_{0}-\varepsilon).$

By Lemma 3.9, since $\operatorname{ct}_{\kappa}$ solves the differential
equation $\dot{x}+x^{2}+\kappa=0$, we have that this inequality holds for
$t\in(\varepsilon,t_{0})$, but the right hand side goes to infinity as
$t\to\varepsilon$. In this case we also have that $b\leq\pi_{\kappa}$, as the
comparison function goes to $-\infty$ as $t$ goes to $\pi_{\kappa}$.

Interchanging the roles of $\rho$ and $\operatorname{ct}_{\kappa}$ in the
previous proof we get the second inequality. In this case we need
$t\in(0,\min\lparen b,\pi_{\kappa}\rparen)$ for $\operatorname{ct}_{\kappa}$
to be well-defined as $\operatorname{ct}_{\kappa}\to-\infty$ when $t$ goes to
$\pi_{\kappa}$ for $\kappa>0$. ∎

We now have everything we need to prove Rauch’s theorem.

######  Theorem 3.11 (Rauch’s theorem).

Let $(M,\textsl{g})$ be a Riemannian manifold with bounded sectional curvature
$\delta\leq\sec\leq\Delta$. Fix a point $p\in M$ and define $r$ as the
distance to $p$. For any other point in a ball $x\in
B_{p}(\pi_{\Delta})\backslash\\{p\\}$, and any $w_{1},w_{2}\in T_{x}M$ we have

$\operatorname{ct}_{\Delta}(r(x))\langle
w^{\perp}_{1},w^{\perp}_{2}\rangle\leq\lparen\operatorname{Hess}r\rparen_{x}(w_{1},w_{2})\leq\operatorname{ct}_{\delta}(r(x))\langle
w^{\perp}_{1},w^{\perp}_{2}\rangle.$

The upper bound also holds for every $x\in
M\backslash\lparen\operatorname{cut}(p)\cup\\{p\\}\rparen$. These bounds are
tight on spaces of constant curvature.

###### Proof.

First, note that since $\operatorname{Hess}r|_{x}$ is a symmetric bilinear
form, it is diagonalizable and it will achieve its maximum and minimum at
certain eigenvectors. As such, we just have to prove the result for
$w_{1}=w_{2}$.

Consider a distance minimizing geodesic $\gamma\colon[0,r]\to M$ connecting
$p$ and $x$. For a radial vector $w=\lVert w\rVert\dot{\gamma}(r)\in T_{x}M$
we evaluate the Hessian on the vector field $\lVert w\rVert\dot{\gamma}$ using
that $\nabla r|_{\gamma}=\dot{\gamma}$

$\operatorname{Hess}r(\lVert w\rVert\dot{\gamma},\lVert
w\rVert\dot{\gamma})=\lVert
w\rVert^{2}\langle\nabla_{\dot{\gamma}}\dot{\gamma},\dot{\gamma}\rangle=0.$

So the Hessian of the distance function is zero in the radial direction.

By linearity, it is enough to prove the theorem for a unitary vector $w$
normal to $\dot{\gamma}(r)$.

For the rest of the bounds, we will take advantage of the simple form that the
curvature equations take when evaluated on Jacobi and parallel vector fields,
as computed in Proposition 3.7.

For the lower bound we consider a Jacobi field such that $J(r)=w$ (i.e.,
$J(0)=0$ and $\dot{J}(0)=\frac{w}{r}$) and define
$\rho(t)=\operatorname{Hess}r(\frac{J}{\lVert J\rVert},\frac{J}{\lVert
J\rVert})$. We can compute its derivative as

$\dot{\rho}=\frac{\mathrm{d}}{\mathrm{d}t}\frac{\operatorname{Hess}r(J,J)}{\langle
J,J\rangle}=\frac{\frac{\mathrm{d}}{\mathrm{d}t}\operatorname{Hess}r(J,J)\lVert
J\rVert^{2}-2\operatorname{Hess}r(J,J)^{2}}{\lVert J\rVert^{4}}.$

Plugging the formula for the derivative of the Hessian into
(LABEL:eq:curv1_jacobi) we get

$\dot{\rho}+2\rho^{2}-\operatorname{Hess}^{2}r\Big{\lparen}\frac{J}{\lVert
J\rVert},\frac{J}{\lVert J\rVert}\Big{\rparen}=-\sec(J,\dot{\gamma}).$

Using Cauchy-Schwarz, we can bound the iterated Hessian for any vector field
$X$ of norm $1$

$\operatorname{Hess}^{2}r(X,X)=\langle\nabla_{X}\nabla r,\nabla_{X}\nabla
r\rangle\geq\langle\nabla_{X}\nabla
r,X\rangle^{2}=\operatorname{Hess}r(X,X)^{2}$

and together with the upper bound on the sectional curvature we get the
expected differential inequality

$\dot{\rho}+\rho^{2}\geq-\Delta.$

Finally, since $\lVert J(0)\rVert=0$ and $\lVert\dot{J}(0)\rVert\neq 0$, we
have that $\lim\limits_{t\to 0}\rho(t)=\infty$, so we may use the second part
of Proposition 3.10 to finish the lower bound.

For the upper bound consider a parallel vector field $E$ such that $E(r)=w$
and define $\rho=\operatorname{Hess}r(E,E)$. Using (7) we get that

$\dot{\rho}+\rho^{2}=-\sec(E,\dot{\gamma})\leq-\delta$

and we finish by applying the first part of Proposition 3.10. ∎

######  Remark.

For the upper bound, we never used the fact that
$x\not\in\operatorname{cut}(p)$. We merely used that the geodesic $\gamma$
does not have any conjugate points on $[0,r]$, to be able to use Proposition
3.10, as if $\lVert J(t)\rVert=0$, then the Hessian of the distance function
at that point would be infinite.

As stated, the theorem is as general as it can be, as the distance function
$r$ is not differentiable on $\operatorname{cut}(p)$. On the other hand, if it
is stated in terms of Jacobi fields, one can further generalize it to Jacobi
fields along geodesics without conjugate points, since both the proof of the
theorem and that of Proposition 3.6 can be done in terms of a geodesic and
short variations of geodesics. We will use this more general version in the
following theorem.

######  Theorem 3.12 (First order bounds for the exponential map).

Let $(M,\textsl{g})$ be a Riemannian manifold with bounded sectional curvature
$\delta\leq\sec\leq\Delta$. Fix a point $p\in M$, for any unit vector $v\in
T_{p}M$ and an $r\in[0,\pi_{\Delta}]$,

$\min\Big{\\{}1,\frac{\operatorname{sn}_{\Delta}(r)}{r}\Big{\\}}\lVert
w\rVert\leq\lVert\lparen\mathrm{d}\exp_{p}\rparen_{rv}(w)\rVert\leq\max\Big{\\{}1,\frac{\operatorname{sn}_{\delta}(r)}{r}\Big{\\}}\lVert
w\rVert\mathrlap{\qquad\forall w\in T_{p}M.}$

The upper bound also holds for $r<r_{\operatorname{conj}}(v)$.

These bounds are tight on spaces of constant curvature.

###### Proof.

The bound by $1$ above and below comes from the Gauss’s lemma, as the
exponential is a radial isometry. From this, and using linearity, we just have
to prove the bound for a vector $w$ in the normal direction with $\lVert
w\rVert=1$.

All we have to do now is to transform the bounds on the Riccati equation on
bounds on the Jacobi equation. Define $\gamma$ as the geodesic starting at $p$
with $\dot{\gamma}(0)=\frac{v}{\lVert v\rVert}$, and choose a Jacobi field
along $\gamma$ such that $J(r)=w$, that is, choose $\dot{J}(0)=\frac{w}{r}$.
As we noted at the beginning of this section, we have that

$\operatorname{Hess}r\Big{\lparen}\frac{J}{\lVert J\rVert},\frac{J}{\lVert
J\rVert}\Big{\rparen}=\frac{\langle\dot{J},J\rangle}{\lVert
J\rVert^{2}}=\frac{\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}$}\lVert
J\rVert}{\lVert J\rVert}=\frac{\mathrm{d}}{\mathrm{d}t}\log\lparen\lVert
J\rVert\rparen$

so Rauch’s theorem may be rewritten on $(0,\pi_{\Delta})$ as

$\frac{\mathrm{d}}{\mathrm{d}t}\log\lparen\operatorname{sn}_{\Delta}(t)\rparen\leq\frac{\mathrm{d}}{\mathrm{d}t}\log\lparen\lVert
J\rVert\rparen\leq\frac{\mathrm{d}}{\mathrm{d}t}\log\lparen\operatorname{sn}_{\delta}(t)\rparen.$

Integrating, using that $\lVert
J(0)\rVert=\operatorname{sn}_{\Delta}(0)=\operatorname{sn}_{\delta}(0)=0$, and
computing the limit using l’Hopital, we get

$\operatorname{sn}_{\Delta}(t)\leq\lVert
J\rVert\leq\operatorname{sn}_{\delta}(t).$

and since $J(r)=\lparen\mathrm{d}\exp_{p}\rparen_{rv}(rw)$,

$\lVert\lparen\mathrm{d}\exp_{p}\rparen_{rv}(w)\rVert=\frac{\lVert
J(r)\rVert}{r}.$

As it was the case in the proof of Rauch’s theorem, we just used that $\gamma$
does not have conjugate points on $[0,r]$, rather than the stronger
$rv\in\overline{U}_{p}$. As such, just by noting that the proof does not rely
on the definition of $\operatorname{Hess}r$, and instead, can be carried out
in terms of Jacobi fields, we get the finer result for the upper bound. ∎

######  Example 3.13.

To see how the definition of the upper bounds are an improvement over just
considering $rv\in\overline{U}_{p}$, consider the flat torus. For this
manifold, $\overline{U}_{p}\subseteq T_{p}M$ is a square centered at $0$. On
the other hand, since for the flat torus we have that
$\operatorname{conj}(p)=\emptyset$, we get that the upper bound holds on all
$T_{p}M$ for every point $p$.

We now have as a corollary of the lower bounds the following strengthening of
the Cartan-Hadarmard theorem that was implicitly present in the statement of
the last theorem.

######  Theorem 3.14.

If $(M,\textsl{g})$ has $\sec\leq\Delta$, then the exponential map $\exp_{p}$
is not singular on $B_{p}\lparen\pi_{\Delta}\rparen$, so
$B_{p}\lparen\pi_{\Delta}\rparen\subseteq M\backslash\operatorname{conj}(p)$
and $\pi_{\Delta}\leq r_{\operatorname{conj}}(v)$ for every point $p\in M$ and
every unit vector $v\in T_{p}M$.

###  3.4 The law of cosines on manifolds

Rauch’s theorem has found countless applications in geometry, as those given
by Rauch himself, who proved with it a weak version of the sphere theorem, or
as used by Gromov or Karcher to give bounds the volume form of a Riemannian
manifold, or to control de size of the fundamental group of a given manifold.
Here we deviate slightly from the presentation to remark one that has been
found particularly useful in the context of optimization on manifolds.

######  Theorem 3.15 (Law of cosines).

Let $(M,\textsl{g})$ be a Riemannian manifold with bounded sectional curvature
$\delta\leq\sec\leq\Delta$. Let $p$ be a point on $M$ and let $x,y\in
B_{p}\lparen R\rparen$ for $R\leq\pi_{\Delta}$ such that the minimizing length
geodesic $\gamma$ that connects them also lies in this ball. Define the angle
$\alpha=\angle pxy$. We then have

$\displaystyle d(y,p)^{2}$
$\displaystyle\leq\zeta_{1,\delta}(R)d(x,y)^{2}+d(x,p)^{2}-2d(x,y)d(x,p)\cos(\alpha)$
$\displaystyle d(y,p)^{2}$
$\displaystyle\geq\zeta_{2,\Delta}(R)d(x,y)^{2}+d(x,p)^{2}-2d(x,y)d(x,p)\cos(\alpha)$

where

$\displaystyle\zeta_{1,\kappa}(r)$
$\displaystyle\coloneqq\max\\{1,r\operatorname{ct}_{\kappa}(r)\\}$
$\displaystyle=\begin{cases}1&\text{if }\kappa\geq 0\\\
\sqrt{-\kappa}r\coth(\sqrt{-\kappa}r)&\text{if }\kappa<0\end{cases}$
$\displaystyle\zeta_{2,\kappa}(r)$
$\displaystyle\coloneqq\min\\{1,r\operatorname{ct}_{\kappa}(r)\\}$
$\displaystyle=\begin{cases}\mathrlap{1}\hphantom{\sqrt{-\kappa}r\coth(\sqrt{-\kappa}r)}&\text{if
}\kappa\leq 0\\\ \sqrt{\kappa}r\cot(\sqrt{\kappa}r)&\text{if
}\kappa>0\end{cases}$

Furthermore, the first bound holds as long as there exists a distance
minimizing geodesic contained in $U_{p}$ that connects $x$ and $y$. These
bounds are tight on spaces of constant curvature.

###### Proof.

Define $f_{p}=\frac{1}{2}r^{2}$. We start by relating the Hessian of $f_{p}$
to the Hessian of $r$ using the chain rule for the differential and the
Leibnitz rule for the connection

$\operatorname{Hess}f_{p}=\nabla\lparen
r\mathrm{d}r\rparen=\mathrm{d}r\otimes\mathrm{d}r+r\operatorname{Hess}r.$

As such, for a radial vector $w$ of norm $1$, since
$\operatorname{Hess}r(w,w)=0$, this equation simplifies to

$\operatorname{Hess}f_{p}(w,w)=\mathrm{d}r(w)^{2}=\langle\nabla
r,w\rangle^{2}=1.$

For normal vectors we may use the bounds on the Hessian of $r$. It is now
evident that the functions $\zeta_{1,\delta}$ and $\zeta_{2,\Delta}$ are
defined to be upper and lower bound on the Hessian of $f_{p}$ acting on
vectors of norm $1$.

Let $\gamma\colon[0,t]\to M$ be a distance-minimizing geodesic such that
$\gamma(0)=x$, $\gamma(t)=y$ and $\gamma(s)\in U_{p}$ for every $s\in[0,t]$.
We will prove the first inequality, as the proof of the second one is
analogous. By the bounds above, using (2), we have that

$\operatorname{Hess}f_{p}(\dot{\gamma}(s),\dot{\gamma}(s))=(f_{p}\circ\gamma)^{\prime\prime}(s)\leq\zeta_{1,\delta}(R)$

for every $s\in[0,t]$, since $f_{p}$ is differentiable on $U_{p}$. Integrating
this inequality we get

$\int_{0}^{t}(f_{p}\circ\gamma)^{\prime\prime}(s)\,\mathrm{d}s=(f_{p}\circ\gamma)^{\prime}(t)-\langle\nabla
f_{p}(\gamma(0)),\dot{\gamma}(0)\rangle\leq\zeta_{1,\delta}(R)t.$

Integrating once again we have that

$(f_{p}\circ\gamma)(t)-(f_{p}\circ\gamma)(0)-\langle\nabla
f_{p}(\gamma(0)),\dot{\gamma}(0)\rangle
t\leq\zeta_{1,\delta}(R)\frac{t^{2}}{2}.$

which, after substituting the values of $f_{p}$ and $\gamma$, gives

$d(y,p)^{2}\leq\zeta_{1,\delta}(R)d(x,y)^{2}+d(x,p)^{2}+2d(x,y)\langle\nabla
f_{p}(x),\dot{\gamma}(0)\rangle.$

Finally, since $\nabla f_{p}(x)=r\nabla r=-\exp_{x}^{-1}(p)$ we get the
desired inequality. ∎

######  Remark (Obstructions to a global law of cosines).

One could ask whether it is possible to extend this local result on $U_{p}$ to
a global result for geodesic triangles whose sides are distance minimizing, as
one does in Toponogov’s theorem. This is, in general, not possible. Consider
the flat torus and a point $p$ in it and an equilateral triangle on a
generating circle. The sides of this triangle are distance minimizing of
length $x$—one third the circumference of the circle—and the angles of this
triangle will be $\alpha=\pi$. Since for the flat torus we have that
$\delta=0$, the first inequality would imply that $3x^{2}\leq 2x^{2}$.

The upper bound for the case $\delta<0$ has found multiples uses in the
optimization literature. Before starting, to be able to put what comes next in
context, we recall the definition of Hadamard manifolds.

######  Definition 3.16 (Hadamard manifold).

A Riemannian manifold is a Hadamard manifold if it is complete, simply
connected, and has non-positive sectional.

These manifolds are called Hadamard manifolds due to the Cartan-Hadamard
theorem.

######  Theorem 3.17 (Cartan-Hadamard).

A Hadamard manifold is diffeomorphic to $\mathbb{R}^{n}$.

###### Proof.

By Theorem 3.14, since $\sec\leq 0$, $\pi_{\Delta}=\infty$ and the
differential exponential map at a point $p\in M$ is always full rank. As such,
by the inverse function theorem, it is a covering map, and $\mathbb{R}^{n}$ is
the universal cover of $M$. Since $M$ is simply connected, $M$ is
diffeomorphic to its universal cover. ∎

It is now clear that Hadamard manifolds are a class of particularly well-
behaved manifolds, in that the $\operatorname{cut}(p)=\emptyset$ for every
$p\in M$, so $U_{p}=M$ and $\overline{U}_{p}=T_{p}M$, and all the theorems in
Section 3 hold for every $v\in T_{p}M$.

If one then considers a Hadamard manifold with sectional curvature bounded
below by $\delta<0$, then one may apply the upper bounds in Theorem 3.15 on
the whole manifold. This idea was heavily exploited in [10] to get convergence
rates for stochastic gradient descent, and in [61] to prove convergence both
stochastic and non-stochastic setting. In [24] they used again the lower
bounds whenever $\delta<0$ to prove convergence rates of certain step-size
schedulers. In [6], they use the law of cosines on manifolds of either non-
negative or non-positive curvature to prove convergence rates for subgradient
methods and proximal-point methods. In [1], the lower bounds are used to prove
convergence rates of certain second order method.

In the previous presentation, we have showed that the Hadamard restriction,
that is, being simply connected and with sectional curvature bounded above, is
not really necessary. In particular, if we have a positive bound on the
curvature, we can still apply these kind of results, at the expense of working
on certain neighbourhood of the initial point. The same happens if we remove
the condition of the manifold being simply connected.

##  4 Second Order Bounds for the Exponential Map

In this section, we give bounds on the Hessian of the differential of the
exponential map. These kind of bounds were first given in the paper [39]. The
proof can just be found in German, so we will give a self-contained
presentation. Our proof uses a comparison lemma developed by Kaul, and then
streamlines all the other technical tools, considerably simplifying the proof
and obtaining tighter explicit bounds.

###  4.1 The differential equation

We start by giving a technical remark on the nature of the Hessian of the
exponential map. The differential of the exponential for a point $(p,v)\in TM$
is a linear map of the form

$\lparen\mathrm{d}\exp_{p}\rparen_{v}\colon T_{p}M\to T_{\exp_{p}(v)}M.$

For this reason, it can be seen as a section of the bundle
$T^{\ast}_{p}M\otimes\exp_{p}^{\ast}(TM)$ over $T_{p}M$. We have connections
$\nabla^{\text{flat}}$ and $\nabla$ on $T^{\ast}_{p}M$ and
$\exp_{p}^{\ast}(TM)$ respectively. Any two connections induce a connection on
a tensor product bundle so that the Leibnitz rule holds. In the concrete
example of the exponential map, this Leibnitz rule for vector fields
$\overline{W}_{1},\overline{W}_{2}$ on $\overline{U}_{p}$, takes the form

$\nabla_{\mathrm{d}\exp_{p}(\overline{W}_{1})}\lparen\mathrm{d}\exp_{p}(\overline{W}_{2})\rparen=\lparen\nabla_{\overline{W}_{1}}\mathrm{d}\exp_{p}\rparen(\overline{W}_{2})+\mathrm{d}\exp_{p}\lparen\nabla^{\text{flat}}_{\overline{W}_{1}}\overline{W}_{2}\rparen.$

Since $\exp_{p}$ is a diffeomorphism between $\overline{U}_{p}$ and $U_{p}$,
writing $W_{i}=\mathrm{d}\exp_{p}(\overline{W}_{i})$ for the pushforward of a
vector from $\overline{U}_{p}$ to $U_{p}$, we see that the Hessian of the
exponential map can be put in terms of the Christoffel symbols in normal
coordinates as

$\lparen\nabla_{\overline{W}_{1}}\mathrm{d}\exp_{p}\rparen(\overline{W}_{2})=\nabla_{W_{1}}W_{2}-\nabla^{\text{flat}}_{\overline{W}_{1}}\overline{W}_{2}=\Gamma(W_{1},W_{2})=\Gamma_{i,j}^{k}W_{1}^{i}W_{2}^{j}\partial_{k}.$
(8)

In other words, the Hessian of the exponential map is exactly the Christoffel
symbols of the connection in normal coordinates at $p$ evaluated at the
pushforward of the vector along $\exp_{p}$. In particular, this shows that the
Hessian of the exponential map is a symmetric bilinear map of the form

$\lparen\nabla\mathrm{d}\exp_{p}\rparen_{v}\colon T_{p}M\times T_{p}M\to
T_{\exp_{p}(v)}M.$

As we did for the first order bounds, our strategy to bound this quantity will
be to deduce a differential equation for $\nabla\mathrm{d}\exp_{p}$ and then
bound the norm of the solutions of this equation. In this case, we will have
to differentiate a second time through a second variation of a geodesic. Given
that this second derivative will not be in the direction of the geodesic, it
will not be enough to have vector fields along the geodesic—we will have to
have them defined also in the direction in which we want to differentiate
them.

######  Proposition 4.1.

Let $(M,\textsl{g})$ be a Riemannian manifold and let $\gamma$ be a geodesic
with initial unit vector $v\in T_{p}M$. For two vectors $w_{1},w_{2}\in
T_{p}M$ perpendicular to $v$, the vector field along $\gamma$

$K(t)\coloneqq\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}(tw_{1},tw_{2})$

satisfies the following second order inhomogeneous linear differential
equation along $\gamma$

$\ddot{K}+R(K,\dot{\gamma})\dot{\gamma}+Y=0\mathrlap{\qquad
K(0)=0,\dot{K}(0)=0}$

where $Y$ is the vector field along $\gamma$ given by

$Y\coloneqq 2R\lparen J_{1},\dot{\gamma}\rparen\dot{J}_{2}+2R\lparen
J_{2},\dot{\gamma}\rparen\dot{J}_{1}+\lparen\nabla_{\dot{\gamma}}R\rparen\lparen
J_{2},\dot{\gamma}\rparen J_{1}+\lparen\nabla_{J_{2}}R\rparen\lparen
J_{1},\dot{\gamma}\rparen\dot{\gamma}$

and $J_{1},J_{2}$ are the Jacobi fields along $\gamma$ with initial conditions
$J_{i}(0)=0$, $\dot{J}_{i}(0)=w_{i}$.

###### Proof.

Define the geodesic variation

$\\!\begin{aligned}
c\colon[0,r]\times\lparen-\varepsilon,\varepsilon\rparen\times\lparen-\varepsilon,\varepsilon\rparen&\to
M\\\ (t,s_{1},s_{2})&\mapsto\exp_{p}\lparen t\lparen
v+s_{1}w_{1}+s_{2}w_{2}\rparen\rparen\end{aligned}$

Consider the first variation along $\gamma(t)=\exp_{p}(tv)$ in the direction
of $s_{1}$ given by the family of Jacobi fields

$\tilde{J}_{1}(t,s)=\left.\frac{\partial c}{\partial
s_{1}}\right|_{s_{1}=0}(t,s)=\lparen\mathrm{d}\exp_{p}\rparen_{t(v+sw_{2})}(tw_{1}).$

We denote by $J_{i}(t)$ the Jacobi field along $\gamma$ with initial
conditions $(0,w_{i})$. In particular, we have that
$J_{1}(t)=\tilde{J}_{1}(t,0)$, so $\tilde{J}_{1}$ is an extension of $J_{1}$
in the direction of $J_{2}$ so that $\nabla_{J_{2}}\tilde{J}_{1}$ is well
defined. Moreover, we have that $\nabla^{\text{flat}}_{J_{2}}\tilde{J}_{1}=0$
as $\tilde{J}_{1}$ is constant in the direction of $J_{2}$, that is, in the
direction of $stw_{2}$ for $s\in(-\varepsilon,\varepsilon)$ and a fixed $t$.
For this reason, the vector field $K$ is the second variation along $\gamma$
in the directions $w_{1},w_{2}$

$K(t)=\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}\lparen
tw_{1},tw_{2}\rparen=\nabla_{J_{2}}\tilde{J}_{1}=\left.\frac{\partial\tilde{J}_{1}}{\partial
s}\right|_{s=0}(t).$

Furthermore, $\tilde{J}_{1}$ satisfies the Jacobi equation in all its domain
with initial conditions $(0,\frac{1}{\lVert v+sw_{1}\rVert})$ so that

$\nabla_{\nabla r}\nabla_{\nabla r}\tilde{J}_{1}+R(\tilde{J}_{1},\nabla
r)\nabla r=0.$

We may then derive a differential equation for $K$ along $\gamma$ by
differentiating this equation in the direction of $J_{2}$

$\nabla_{J_{2}}\nabla_{\nabla r}\nabla_{\nabla
r}\tilde{J}_{1}+\nabla_{J_{2}}\lparen R(\tilde{J}_{1},\nabla r)\nabla
r\rparen=0.$ (9)

Note that, as we just are interested in deriving a differential equation for
$K$ along $\gamma$, we will use that $\nabla r|_{\gamma}=\dot{\gamma}$ and
$\tilde{J}_{1}=J_{1}$ whenever we have that a quantity just depends on the
values of the vector fields involved along $\gamma$, rather than in a
neighborhood around $\gamma$ in the direction of $J_{2}$.

We may put the first term of (9) in terms of $K$ by repeatedly using the
definition of the curvature tensor and the fact that since $[J_{2},\nabla
r]|_{\gamma}=0$, we have that $\nabla_{J_{2}}\nabla
r=\nabla_{\dot{\gamma}}J_{2}=\dot{J}_{2}$

$\displaystyle\nabla_{J_{2}}\nabla_{\nabla r}\nabla_{\nabla r}\tilde{J}_{1}$
$\displaystyle=\nabla_{\dot{\gamma}}\nabla_{J_{2}}\nabla_{\nabla
r}\tilde{J}_{1}+R(J_{2},\dot{\gamma})\nabla_{\dot{\gamma}}J_{1}$
$\displaystyle=\nabla_{\dot{\gamma}}\nabla_{\dot{\gamma}}K+\nabla_{\dot{\gamma}}\lparen
R\lparen J_{2},\dot{\gamma}\rparen
J_{1}\rparen+R(J_{2},\dot{\gamma})\dot{J}_{1}$
$\displaystyle=\nabla_{\dot{\gamma}}\nabla_{\dot{\gamma}}K+\lparen\nabla_{\dot{\gamma}}R\rparen\lparen
J_{2},\dot{\gamma}\rparen J_{1}+R\lparen\dot{J}_{2},\dot{\gamma}\rparen
J_{1}+2R\lparen J_{2},\dot{\gamma}\rparen\dot{J}_{1}.$

We can expand the second term as

$\nabla_{J_{2}}\lparen R(\tilde{J}_{1},\nabla r)\nabla
r\rparen=\lparen\nabla_{J_{2}}R\rparen\lparen
J_{1},\dot{\gamma}\rparen\dot{\gamma}+R\lparen
K,\dot{\gamma}\rparen\dot{\gamma}+R\lparen
J_{1},\dot{J}_{2}\rparen\dot{\gamma}+R\lparen
J_{1},\dot{\gamma}\rparen\dot{J}_{2}.$

Putting everything together and using the symmetries of the curvature tensor
and the first Bianchi identity we get the differential equation.

For the initial conditions, we have that

$K(0)=\lparen\nabla\mathrm{d}\exp_{p}\rparen_{0}(0,0)=0.$

Differentiating the Hessian in the direction of
$\dot{\gamma}=\frac{\mathrm{d}}{\mathrm{d}t}$ we have that

$\dot{K}(0)=\lparen\nabla\nabla\mathrm{d}\exp_{p}\rparen_{0}(v,0,0)+\lparen\nabla\mathrm{d}\exp_{p}\rparen_{0}(w_{1},0)+\lparen\nabla\mathrm{d}\exp_{p}\rparen_{0}(0,w_{2})=0.\qed$

###  4.2 Manifolds of bounded geometry

By Proposition 4.1, we have that the Hessian of the exponential map solves a
Jacobi-like differential equation with two main differences: It is
inhomogeneous and depends on the covariant derivative of the curvature tensor.

The inhomogeneity of the differential equation will stop us from simplifying
the differential equation into two first order differential equations, as we
did for the Jacobi equation. In contrast, this time we will have to deal
directly with the second order equation. Luckily, we have already developed
most of the tools to do so.

The fact that involves the covariant derivative of the curvature tensor is a
more intrinsic difference. Since the sectional curvature completely defines
the curvature tensor, bounds on the sectional curvature can be translated into
bounds on the norm of the curvature tensor. One question that naturally arises
is whether these $0$-th order bounds also give first order bounds. As one may
expect, this is not the case.

######  Example 4.2 (Manifold of bounded $0$-th order and unbounded $1$-st
order).

Consider a rotationally symmetric surface of the form
$\mathrm{d}r^{2}+\rho(r)^{2}\mathrm{d}\theta^{2}$ on the cylinder
$(0,1)\times\SS^{1}$ for a function $\rho>0$. These metrics have a
particularly simple formula for the sectional curvature

$\sec(\partial_{r},\partial_{\theta})=\frac{R(\partial_{r},\partial_{\theta},\partial_{\theta},\partial_{r})}{\rho^{2}}=-\frac{\ddot{\rho}}{\rho}.$

All we have to do now is to choose a positive function $\rho$ with bounded
second derivative and arbitrarily large third derivative, for example
$\rho(r)=2+r^{5}\sin(1/r)$. This metric has $\lvert\sec\rvert<12$, but

$\lparen\nabla_{\partial_{r}}R\rparen\lparen\partial_{r},\partial_{\theta},\partial_{\theta},\partial_{r}\rparen=D_{\partial_{r}}\lparen
R\lparen\partial_{r},\partial_{\theta},\partial_{\theta},\partial_{r}\rparen\rparen-2R\lparen\partial_{r},\nabla_{\partial_{r}}\partial_{\theta},\partial_{\theta},\partial_{r}\rparen=\dot{\rho}\ddot{\rho}-\rho{\mathop{\kern
0.0pt\rho}\limits^{\makebox[0.0pt][c]{\vbox
to-2.2pt{\kern-2.0pt\hbox{\scalebox{0.89}{.}\kern-0.5pt\scalebox{0.89}{.}\kern-0.5pt\scalebox{0.89}{.}}\vss}}}}$

which is unbounded as $r\to 0$.

This example motivates the following definition.

######  Definition 4.3 (Manifold of bounded geometry).

Let $(M,\textsl{g})$ be a Riemannian manifold. We say that a manifold has
$k$-bounded geometry if for every $i=0,\dots,k$ there exist constants
$C_{i}\geq 0$ such that

$\lVert\nabla^{i}R\rVert\leq C_{i}$

and the injectivity radius is uniformly bounded below by a positive constant

$\operatorname{r_{inj}}\coloneqq\inf_{p\in M}\operatorname{r_{inj}}(p)>0.$

######  Remark.

A manifold has $0$-bounded geometry if it has bounded sectional curvature
above and below and positive injectivity radius. As one may expect, given that
the $1$-bounded geometry is given by the Christoffel symbols, which are given
by directional derivatives of the metric g in normal coordinates, the
condition on the boundedness of the derivatives of the curvature tensor is
equivalent to asking for the entries of the metric tensor to be bounded in
normal coordinates. This was first proved in [19].

It is clear that any compact manifold with any given metric will be of
$k$-bounded geometry for every $k\in\mathbb{N}$. A natural question would be
whether this condition puts any constraint in the topology of the manifold.
This was answered in the negative in [27].

######  Theorem 4.4 (Greene).

Every differentiable manifold admits a complete metric of bounded geometry.

Given that we usually have access to upper and lower bounds on the sectional
curvature, and given the form that takes the covariant derivatives in
Proposition 4.1, we give the following definition, which is just a
generalisation of that given in the introduction.

######  Definition 4.5.

We say that a Riemannian manifold $(M,\textsl{g})$ has
$(\delta,\Delta,\Lambda)$-bounded geometry if $\operatorname{r_{inj}}>0$ and

$\displaystyle\delta\leq\sec\leq\Delta$
$\displaystyle\lVert\lparen\nabla_{u}R\rparen(w,u)w+\lparen\nabla_{w}R\rparen\lparen
w,u\rparen u\rVert\leq 2\Lambda\lVert u\rVert^{2}\lVert
w\rVert^{2}\mathrlap{\qquad\forall p\in M\quad u,w\in T_{p}M.}$

######  Remark.

If a manifold is of $(\delta,\Delta,\Lambda)$-bounded geometry according to
the definition in the introduction, it is of $(\delta,\Delta,\Lambda)$-bounded
geometry according to this definition, but not necessarily the other way
around. In the proofs, we will just use this definition, as it is all we need,
and it also heavily simplifies the computations when one wants to estimate the
actual constant $\Lambda$ for a given manifold.

######  Remark.

We will not use the fact that these manifolds have uniformly bounded
injectivity radius in this paper, but this fact turns out to be crucial when
proving other flavour of theorems in optimization of manifolds that involve
$\exp^{-1}_{p}$ or parallel transport along geodesics or a retraction.

###  4.3 Curvature bounds

We start by recalling the following lemma that gives a closed formula for the
curvature tensor of manifolds of constant sectional curvature. In its more
general form, it says that the Riemannian curvature tensor is completely
specified by the sectional curvature.

######  Lemma 4.6.

Let $(M,\textsl{g})$ be a Riemannian manifold of constant sectional curvature
$\kappa\in\mathbb{R}^{n}$. The curvature tensor takes the form

$R(x,y)z\coloneqq R_{\kappa}(x,y)z=\kappa\lparen\langle z,y\rangle x-\langle
z,x\rangle y\rparen\mathrlap{\qquad\forall x,y,z\in T_{p}M.}$

To be able to bound the terms in the differential equation, we will need some
estimates on the norm of the curvature tensor in terms of the sectional
curvature. These are a generalization of Berger’s lemma [7]. The first three
inequalities are announced—the last one with a worse constant—in [39] citing
for their proof a monograph by Karcher that was never published. We provide
original proofs for these inequalities. The third inequality appears in [36],
but with a different proof. The last inequality is entirely from [36].

######  Proposition 4.7.

Let $(M,\textsl{g})$ be a Riemannian manifold with bounded sectional curvature
$\delta\leq\sec\leq\Delta$ and define

$\varepsilon=\frac{\Delta-\delta}{2}\qquad\mu=\frac{\Delta+\delta}{2}\qquad
K=\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}.$

Given a point $p\in M$, we have the following inequalities

$\displaystyle\lvert R(x,y,y,w)-R_{\mu}(x,y,y,w)\rvert$
$\displaystyle\leq\varepsilon\lVert x\rVert\lVert y\rVert^{2}\lVert
w\rVert\mathrlap{\qquad\forall x,y,w\in T_{p}M}$ (10) $\displaystyle\lVert
R(x,y,z,w)-R_{\mu}(x,y,z,w)\rVert$
$\displaystyle\leq\frac{4}{3}\varepsilon\lVert x\rVert\lVert y\rVert\lVert
z\rVert\lVert w\rVert\mathrlap{\quad\forall x,y,z,w\in T_{p}M}$ (11)

The constants $1$ and $\frac{4}{3}$ are tight.

Furthermore, if $e\in T_{p}M$ is a unitary vector and we define
$R^{\perp}(x,y)z$ as the component of the curvature tensor perpendicular to
$e$, we have that for every $y,z$ perpendicular to $e$

$\displaystyle\lVert R^{\perp}(e,y)z\rVert$
$\displaystyle\leq\frac{4}{3}\varepsilon\lVert y\rVert\lVert z\rVert$ (12)
$\displaystyle\lVert R(e,y)z\rVert$ $\displaystyle\leq 2K.$ (13)

###### Proof.

We may assume that the vectors $x,y,z,w$ are of norm $1$ by linearity. For
(10), since the curvature tensor is skew-symmetric in the first two components
and the last two components, we may assume that $x,w$ are perpendicular to
$y$. Fix a vector $y$ and consider the bilinear form

$T_{y}(x,w)=R(x,y,y,w)-R_{\mu}(x,y,y,w).$

This application is symmetric, and as such, it attains its maximum at an
eigenvector $x_{1}$ of norm one orthogonal to $y$ so that for every $x$

$\lvert T_{y}(x,w)\rvert\leq\lvert
R(x_{1},y,y,x_{1})-R_{\mu}(x_{1},y,y,x_{1})\rvert=\lvert\sec(x_{1},y)-\mu\rvert\leq\varepsilon.$

For the second inequality, by linearity we can assume that all the vectors
have the same norm. By polarization we have that

$\displaystyle 6R(x,y,z,w)$ $\displaystyle=R(x,y+z,y+z,w)-R(x,y-z,y-z,w)$
$\displaystyle-R(y,x+z,x+z,w)+R(y,x-z,x-z,w).$

We then apply the triangle inequality to the expression for
$(R-R_{\mu})(x,y,z,w)$ together with (10) to get

$6\lvert(R-R_{\mu})(x,y,z,w)\rvert\leq\varepsilon\Big{\lparen}\lVert
x\rVert\lVert w\rVert\lparen\lVert y+z\rVert^{2}+\lVert
y-z\rVert^{2}\rparen+\lVert y\rVert\lVert w\rVert\lparen\lVert
x+z\rVert^{2}+\lVert x-z\rVert^{2}\rparen\Big{\rparen}$

and using the parallelogram law together with the fact that all the vectors
have the same norm, we get that

$\lVert x+y\rVert^{2}+\lVert x-y\rVert^{2}=2(\lVert x\rVert^{2}+\lVert
y\rVert^{2})=4\lVert x\rVert\lVert y\rVert$

and the second inequality follows.

These two inequalities are tight on $\mathbb{C}P^{n}$ seen as a real manifold
[36].

The bound on the normal part of the curvature tensor follows directly from
(11), since $R_{\mu}(e,y,z,w)=0$ whenever $y,z,w$ are perpendicular to $e$, so
choosing $w$ as the unitary vector in the direction of $R^{\perp}(e,y)z$ we
get that

$\lVert R^{\perp}(e,y)z\rVert=\lvert
R(e,y,z,w)\rvert\leq\frac{4}{3}\varepsilon\lVert y\rVert\lVert z\rVert.$

The last inequality is proved in [36]. ∎

###  4.4 A comparison lemma

We shall proceed in a very similar way to how we did in the case of the first
order equation. We will simplify the differential equation to one in one
dimension, and there, we will use a comparison theorem for functions of real
variable. The only difference is that, in this case, we will not be able to
simplify the computations to a Riccati equation. We start by proving the
second order version of the Riccati comparison lemma but for the Jacobi
equation.

######  Lemma 4.8 (Jacobi comparison lemma).

Fix a real number $\kappa\in\mathbb{R}$, and let
$f,g\colon[0,r]\to\mathbb{R}^{+}$ be two functions such that

$\ddot{f}+\kappa f\leq\ddot{g}+\kappa g\mathrlap{\qquad
f(0)=g(0),\dot{f}(0)\leq\dot{g}(0).}$

Then $f\leq g$ on $[0,\min\\{r,\pi_{\kappa}\\}]$.

###### Proof.

Let $h=g-f$ and $\zeta=\ddot{h}+\kappa h\geq 0$. Let us show that the solution
to the linear inhomogeneous initial value problem

$\ddot{h}+\kappa h=\zeta\mathrlap{\qquad h(0)=0,\ \dot{h}(0)\geq 0}$

is indeed positive on the given interval. Solving the equation, we find that
the solution is given by

$\displaystyle h(x)$
$\displaystyle=\operatorname{sn}_{\kappa}(x)\int_{0}^{x}\operatorname{sn}^{\prime}_{\kappa}(t)\zeta(t)\,\mathrm{d}t-\operatorname{sn}_{\kappa}^{\prime}(x)\int_{0}^{x}\operatorname{sn}_{\kappa}(t)\zeta(t)\,\mathrm{d}t+\dot{h}(0)\operatorname{sn}_{\kappa}(x)$
$\displaystyle=\int_{0}^{x}\operatorname{sn}_{\kappa}(x-t)\zeta(t)\,\mathrm{d}t+\dot{h}(0)\operatorname{sn}_{\kappa}(x).$

where we have used the trigonometric identity

$\operatorname{sn}_{\kappa}(x-t)=\operatorname{sn}_{\kappa}(x)\operatorname{sn}^{\prime}_{\kappa}(t)-\operatorname{sn}^{\prime}_{\kappa}(x)\operatorname{sn}_{\kappa}(t).$

From this we see that $h(x)\geq 0$ on $[0,\min\\{r,\pi_{\kappa}\\}]$ as
$\operatorname{sn}_{\kappa}$, $\zeta$ are positive on this interval. ∎

We can now show how to estimate the inhomogeneous differential equation for
the Hessian of the exponential, taking advantage of the fact that the initial
conditions are both zero.

######  Proposition 4.9 (Kaul, 1976).

Let $(M,\textsl{g})$ be a Riemannian manifold with bounded sectional curvature
$\delta\leq\sec\leq\Delta$. Let $\gamma\colon[0,r]\to M$ be a geodesic, and
let $X,Y$ be vector fields along $\gamma$ with $X,Y\perp\dot{\gamma}$ such
that

$\ddot{X}+R(X,\dot{\gamma})\dot{\gamma}=Y\mathrlap{\qquad
X(0)=0,\dot{X}(0)=0.}$

Assume that there exists a continuous $\eta$ function such that $\lVert
Y\rVert\leq\eta$ on $[0,r]$. Then, we have that $\lVert X\rVert\leq\rho$ on
$[0,\min\\{r,\pi_{\frac{\Delta+\delta}{2}}\\}]$, where $\rho$ is the solution
of

$\ddot{\rho}+\delta\rho=\eta\mathrlap{\qquad\rho(0)=0,\dot{\rho}(0)=0.}$

###### Proof.

We define again the quantities

$\varepsilon=\frac{\Delta-\delta}{2}\qquad\mu=\frac{\Delta+\delta}{2}.$

Fix a $t_{0}\in[0,r]$ and let $E$ be the parallel vector field along $\gamma$
such that $E(t_{0})=\frac{X(t_{0})}{\lVert X(t_{0})\rVert}$. The function
$f=\langle X,E\rangle$ satisfies the differential equation

where we have used Lemma 4.6 in the first equality and (10) for the bound.

Define $g$ as the solution to the differential equation

$\ddot{g}+\mu g=\lVert Y\rVert+\varepsilon\lVert X\rVert\mathrlap{\qquad
g(0)=\dot{g}(0)=0.}$

Using that $f(t_{0})=\lVert X(t_{0})\rVert$ and that $f(0)=\dot{f}(0)=0$, by
Lemma 4.8, we get that $\lVert X(t_{0})\rVert\leq g(t_{0})$, and since the
definition of $g$ does not depend on the chosen $t_{0}$, we may do this for
any $t_{0}\in[0,r]$ getting that $\lVert X\rVert\leq g$ on
$[0,\min\\{\pi_{\mu},r\\}]$. Now, Using that $\lVert X\rVert\leq g$, we get
that $g$ satisfies the inequality

$\ddot{g}+\delta g\leq\lVert Y\rVert\leq\eta.$

So, for the solution of

$\ddot{\rho}+\delta\rho=\eta\mathrlap{\qquad\rho(0)=0,\dot{\rho}(0)=0.}$

we have that $g\leq\eta$ on
$[0,\min\\{\pi_{\delta},r\\}]\supseteq[0,\min\\{\pi_{\mu},r\\}]$, getting the
result. ∎

###  4.5 A second order version of Rauch’s theorem

We are now ready to give bounds on the Hessian of the exponential map. We
first note that, since our goal is to bound the norm of the Hessian of
$f\circ\exp_{p}$, as this Hessian is symmetric, it attains its maximum at an
eigenvector. As such, we just need to bound the quantity

$\lparen\nabla\mathrm{d}\exp_{p}\rparen_{v}(w,w)\mathrlap{\qquad
v\in\overline{U}_{p},w\in T_{p}M.}$

For this reason, we will start by giving bounds on the diagonal of the Hessian
of the exponential map. We will then see that, since the Hessian is a
symmetric bilinear map, we can leverage these bounds to give bounds on the
full Hessian for any pair of vectors $w_{1},w_{2}\in T_{p}M$.

We give the second order bounds on the exponential map for a manifold of
$(\delta,\Delta,\Lambda)$-bounded geometry (cf. Definition 4.5).

######  Theorem 4.10 (Second order bounds for the exponential map).

Let $(M,\textsl{g})$ be a Riemannian manifold with
$(\delta,\Delta,\Lambda)$-bounded geometry. For a geodesic
$\gamma\colon[0,r]\to M$ with initial unit vector $v$, and a vector $w\in
T_{p}M$, we have that

  * •

If $w$ is radial to $\dot{\gamma}(0)$,

$\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}(w,w)=0.$

  * •

If $w$ is normal to $\dot{\gamma}(0)$, the radial part of the Hessian is
bounded as

$\Big{\lparen}\frac{1}{r}-\frac{\operatorname{sn}_{4\delta}(r)}{r^{2}}\Big{\rparen}\lVert
w\rVert^{2}\leq\langle\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen,\dot{\gamma}(r)\rangle\leq\Big{\lparen}\frac{1}{r}-\frac{\operatorname{sn}_{4\Delta}(r)}{r^{2}}\Big{\rparen}\lVert
w\rVert^{2}$

for $r<\pi_{\Delta}$ for the upper bound and $r<r_{\operatorname{conj}}(v)$
for the lower bound.

The normal part of the Hessian is bounded for
$r<\pi_{\frac{\Delta+\delta}{2}}$ as

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen^{\perp}\rVert\leq\rho(r)\lVert w\rVert^{2}$

where

$\rho(t)=\mbox{{}\small$\displaystyle\frac{8}{9r^{2}}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{t}{2}$}\big{\rparen}^{2}\big{\lparen}3\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{t}{2}$}\big{\rparen}^{2}+2\lparen\Delta-\delta\rparen\operatorname{sn}_{\delta}\lparen
t\rparen\big{\rparen}.$

Both bounds and their radii are tight in spaces of constant curvature.

Furthermore, the radius of convergence for the normal part is tight for
$\operatorname{SO}\lparen n\rparen$.

###### Proof.

As in Proposition 4.1, we write

$\displaystyle\tilde{J}(t,s)$
$\displaystyle=\lparen\mathrm{d}\exp_{p}\rparen_{t(v+sw)}\Big{\lparen}\frac{t}{r}w\Big{\rparen}$
$\displaystyle J(t)$ $\displaystyle=\tilde{J}(t,0)$ $\displaystyle K(t)$
$\displaystyle=\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}\Big{\lparen}\frac{t}{r}w,\frac{t}{r}w\Big{\rparen}.$

By linearity of the Hessian, it is enough to prove the result for $\lVert
w\rVert=1$.

If $w$ is radial, in Proposition 4.1 we have that $Y=0$, so $K$ is a solution
of the equation

$\ddot{K}+R(K,\dot{\gamma})\dot{\gamma}=0\mathrlap{\qquad
K(0)=0,\dot{K}(0)=0}$

and since it is a homogeneous second order linear equation with zero as the
initial condition, its solution is $K(t)=0$ for $t\in[0,r]$.

If $w$ is normal, we need to bound the quantity $\nabla_{J}\tilde{J}$. Note
that this is a vector field along $\gamma$, but to take the derivative in the
direction of $J$ we need to have $\tilde{J}$ defined in that direction as
well. We start by bounding its radial part

$\langle\nabla_{J}\tilde{J},\nabla r\rangle=D_{J}\langle\tilde{J},\nabla
r\rangle-\langle\tilde{J},\nabla_{J}\nabla
r\rangle=D_{J}\langle\tilde{J},\nabla r\rangle-\operatorname{Hess}r(J,J).$
(14)

We can compute the first term directly. By Gauss’s lemma we can simplify this
derivative to one in $T_{p}M$

$D_{J}\langle\tilde{J},\nabla
r\rangle(\gamma(t))=\frac{1}{r}\left.\frac{\mathrm{d}}{\mathrm{d}s}\right|_{s=0}\Big{\langle}\frac{t}{r}w,\frac{t\lparen
v+sw\rparen}{\lVert t\lparen
v+sw\rparen\rVert_{T_{p}M}}\Big{\rangle}_{T_{p}M}=\frac{t}{r^{2}}.$

Using the bounds on the Hessian of the distance function given in Theorem 3.11
we can bound the second term in (14). Evaluating these two quantities at $r$
and using the bounds on the size of the Jacobi fields together with the
trigonometric identity

$\operatorname{sn}_{\kappa}(t)\operatorname{sn}^{\prime}_{\kappa}(t)=\frac{\operatorname{sn}_{\kappa}(2t)}{2}=\operatorname{sn}_{4\kappa}(t)$
(15)

we get the bounds on the tangential part of the Hessian of the exponential.

Finally, for its normal part, consider the differential equation given by
Proposition 4.1. Since $\langle
K,\dot{\gamma}\rangle^{\prime}=\langle\dot{K},\dot{\gamma}\rangle$, the radial
(resp. normal) part of the derivative is the derivative of the radial (resp.
normal) part. For this reason, we have that

$\lparen
K^{\perp}\rparen^{\prime\prime}+R(K^{\perp},\dot{\gamma})\dot{\gamma}=-Y^{\perp}.$

Therefore, we just have to bound $\lVert Y^{\perp}\rVert$ to be able to use
Proposition 4.9 and finish.

We start by giving a bound on the norm of $\dot{J}$. Using that
$\dot{J}=\nabla_{J}\nabla r=\operatorname{Hess}r(J)$,

$\lVert\dot{J}\rVert\leq\lVert\operatorname{Hess}r\rVert\lVert J\rVert.$ (16)

We can then bound the norm of $Y^{\perp}$ as

$\displaystyle\lVert Y^{\perp}\rVert$
$\displaystyle\leq\lVert\lparen\nabla_{\dot{\gamma}}R\rparen\lparen
J,\dot{\gamma}\rparen J+\lparen\nabla_{J}R\rparen\lparen
J,\dot{\gamma}\rparen\dot{\gamma}\rVert+4\lVert R^{\perp}\lparen
J,\dot{\gamma}\rparen\dot{J}\rVert$ $\displaystyle\leq 2\Lambda\lVert
J\rVert^{2}+\frac{8(\Delta-\delta)}{3}\lVert J\rVert\lVert\dot{J}\rVert$
$\displaystyle\leq\Big{\lparen}2\Lambda+\frac{8(\Delta-\delta)}{3}\operatorname{ct}_{\delta}(t)\Big{\rparen}\lVert
J\rVert^{2}$
$\displaystyle\leq\Big{\lparen}2\Lambda+\frac{8(\Delta-\delta)}{3}\operatorname{ct}_{\delta}(t)\Big{\rparen}\frac{\operatorname{sn}_{\delta}(t)^{2}}{r^{2}}$
$\displaystyle=\frac{2}{r^{2}}\Big{\lparen}\Lambda\operatorname{sn}_{\delta}(t)^{2}+\frac{2(\Delta-\delta)}{3}\operatorname{sn}_{\delta}(2t)\Big{\rparen}.$

where we have used (12)—as since $J$ is perpendicular to $\dot{\gamma}$, so is
$\dot{J}$—to bound the normal part of the curvature tensor. We have also used
the first order bound on the differential of the exponential for perpendicular
initial conditions to bound the norm of the Jacobi fields and (16) to bound
their derivative. In the last equality we have used (15) again.

The result follows by noting that $\rho$ is the solution to the differential
equation

$\ddot{\rho}+\delta\rho=\eta\mathrlap{\qquad\rho(0)=\dot{\rho}(0)=0}$

where $\eta$ is the bound on $\lVert Y^{\perp}\rVert$ and applying Proposition
4.9.

We will see that the radius for the bound of the normal part is tight in the
case of $\operatorname{SO}\lparen n\rparen$ in Example 4.15. ∎

The first thing to note is that these bounds go to zero as $r$ tends to zero.
This is exactly what we expect, as the Christoffel symbols in normal
coordinates vanish at the origin.

Bounds on the whole Hessian can be easily deduced from those presented here
via polarization.

The bounds in this theorem provide a notable improvement compared to the best
bounds previously known (cf. [39]). For one, these bounds are tighter. We also
have that these bounds are explicit, compared to the previous bounds, which
were given in terms of a solution of a differential equation that did not have
an explicit integral. We also notably simplified the technical tools necessary
to get to these bounds.

We finish this section by giving a bound on the full Hessian of the
exponential.

######  Theorem 4.11 (Bounds on the Full Hessian).

Let $(M,\textsl{g})$ be a Riemannian manifold with
$(\delta,\Delta,\Lambda)$-bounded geometry. For a geodesic
$\gamma\colon[0,r]\to M$ with initial unit vector $v$,
$r<\pi_{\frac{\Delta+\delta}{2}}$, and any two vectors $w_{1},w_{2}\in
T_{p}M$, we have that

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w_{1},w_{2}\rparen\rVert\leq\mbox{{}\small$\displaystyle\frac{8}{3r^{2}}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}\lparen\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}\lparen
r\rparen\rparen\lVert w_{1}\rVert\lVert w_{2}\rVert.$

Furthermore, the radius $\pi_{\frac{\Delta+\delta}{2}}$ is tight for
$\operatorname{SO}\lparen n\rparen$.

###### Proof.

To bound the full Hessian we first see that it is just enough to bound its
diagonal part. For any symmetric bilinear form and any two vectors we have the
polarisation formula

$4\Phi(u,v)=\Phi(u+v,u+v)-\Phi(u-v,u-v).$

By linearity, we may assume that $\lVert u\rVert=\lVert v\rVert=1$. Taking
absolute values and applying the triangle inequality and Cauchy-Schwarz, we
get that

$4\lVert\Phi(u,v)\rVert\leq\lVert\Phi|_{\operatorname{diag}}\rVert\lparen\lVert
u+v\rVert^{2}+\lVert
u-v\rVert^{2}\rparen=4\lVert\Phi|_{\operatorname{diag}}\rVert$

where $\lVert\Phi|_{\operatorname{diag}}\rVert$ is the operator norm of the
application $u\mapsto\Phi(u,u)$. For this reason, it is enough to bound the
map $u\mapsto\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}\lparen u,u\rparen$.

Let $w$ be a vector normal to $v$. As we did in Theorem 4.10, we consider the
differential equation for

$K(t)=\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}\Big{\lparen}\frac{t}{r}w,\frac{t}{r}w\Big{\rparen}.$

This is exactly the same differential equation that we had for $K^{\perp}$,
only that rather than having to bound the normal part of the curvature tensor,
we bound the full curvature tensor. For that we use (13), and following with
the bounds as we did in Theorem 4.10 and solving the resulting equation we get
that for $w$ perpendicular to $v$, the norm of $K$ is bounded by the solution
of the differential equation

$\ddot{\rho}+\delta\rho=\frac{2}{r^{2}}\lparen\Lambda\operatorname{sn}_{\delta}(t)^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}(2r)\rparen\mathrlap{\qquad\rho(0)=\dot{\rho}(0)=0}$

which is solved by

$\rho(t)=\mbox{{}\small$\displaystyle\frac{8}{3r^{2}}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{t}{2}$}\big{\rparen}^{2}\lparen\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{t}{2}$}\big{\rparen}^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}\lparen
t\rparen\rparen.\qed$

######  Remark (Tighter bounds).

Another way to obtain bounds on the full Hessian of the exponential would be
to take the bounds from Theorem 4.10 for the normal and parallel part of the
Hessian and using the Cauchy-Schwarz inequality to get bounds of the form

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w_{1},w_{2}\rparen\rVert\leq\sqrt{\sigma^{2}+\rho^{2}}\lVert w_{1}\rVert\lVert
w_{2}\rVert.$

Where $\rho$ is the bound of the normal part defined in Theorem 4.10 and

$\sigma=\frac{1}{r^{2}}\max\\{\lvert
r-\operatorname{sn}_{4\delta}(r)\rvert,\lvert
r-\operatorname{sn}_{4\Delta}(r)\rvert\\}$

that is, $\sigma$ is the bound on the norm of the parallel part of the
Hessian.

This bound, although tighter in most specific examples, might be more
difficult to manipulate and lacks the simplicity of that presented in the
theorem above.

###  4.6 Concrete second order bounds

####  4.6.1 Constant curvature: Sphere, hyperbolic space, and Euclidean space

The simplest family to evaluate these bounds in is that of the spaces of
constant curvature such as the sphere, the hyperbolic plane, the Euclidean
space or the flat torus.

For a manifold of constant curvature $\kappa$, since the curvature is
constant, the derivative of the curvature tensor is zero. If we further assume
that their injectivity radius is positive111This is not really used to prove
the bounds, but it is part of the definition of bounded geometry., as is the
case of the hyperbolic space and the sphere, we have that they have
$(\kappa,\kappa,0)$-bounded geometry. Instantiating the bounds for these
spaces we see that

$\displaystyle\langle\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen,\dot{\gamma}(r)$
$\displaystyle=\Big{\lparen}\frac{1}{r}-\frac{\operatorname{sn}_{4\kappa}(r)}{r^{2}}\Big{\rparen}\lVert
w\rVert^{2}$
$\displaystyle\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen^{\perp}\rVert$ $\displaystyle=0$

for $r<\pi_{\kappa}$. As a first sanity check, we see that for the flat case
$\kappa=0$, the radial part is exactly equal to zero, as de whole Hessian of
the exponential map is everywhere zero in this case—the second derivative of
an affine function is zero.

To check that this solution is actually correct, we can solve the differential
equation in Proposition 4.1 for these spaces for

$K(t)\coloneqq\lparen\nabla\mathrm{d}\exp_{p}\rparen_{tv}(tw_{1},tw_{2}).$

Using Lemma 4.6, together with the fact that, since the sectional curvature is
constant, $\nabla R=0$, and the trigonometric identity

$\operatorname{sn}_{\kappa}(t)\operatorname{sn}^{\prime}_{\kappa}(t)=\frac{\operatorname{sn}_{\kappa}(2t)}{2}=\operatorname{sn}_{4\kappa}(t)$

we have that the differential equation for the constant curvature case is
given by

$\ddot{K}(t)+\kappa
K^{\perp}(t)=4\kappa\operatorname{sn}_{4\kappa}(t)\dot{\gamma}(t)\mathrlap{\qquad
K(0)=0,\ \dot{K}(0)=0}$

We can split this equation into its normal and radial part. The normal part is
clearly zero, since the right hand side is zero, and the remaining equation is
linear with zero as the initial condition. For the radial part, setting
$x=\langle K,\dot{\gamma}\rangle$ we get

$\ddot{x}(t)=4\kappa\operatorname{sn}_{4\kappa}(t)\mathrlap{\qquad x(0)=0,\
\dot{x}(0)=0.}$

Finally, $\frac{x(r)}{r^{2}}=\frac{r-\operatorname{sn}_{4\kappa}(r)}{r^{2}}$
is exactly the value announced before.

####  4.6.2 Locally symmetric spaces: Orthogonal group, Symmetric positive
definite matrices, and Grassmannian manifold

Locally symmetric spaces define a large family of particularly well-behaved
Riemannian manifolds. Examples of these spaces are the flat torus, the
orthogonal group (or any compact Lie group with a bi-invariant metric), the
space of symmetric positive definite matrices, the Grassmannian, the oriented
Grassmannian and the hyperbolic Grassmannian.222All these manifolds are to be
regarded as Riemannian manifolds with the metric inherited from their quotient
structure

We recall the algebraic definition of a locally symmetric space.

######  Definition 4.12.

A Riemannian manifold $(M,\textsl{g})$ is _locally symmetric_ if the curvature
tensor is covariantly constant, that is, $\nabla R=0$.

Locally symmetric spaces were introduced and studied by Cartan in 1926, who
also gave a complete classification of these in 1932.

The most notable examples of locally symmetric spaces are symmetric spaces
which are one of the most important families of real manifolds in Riemannian
geometry. These were intensively studied by Sigurður Helgason [29].

######  Definition 4.13.

A Riemannian manifold $(M,\textsl{g})$ is a _symmetric space_ if, for every
point $p\in M$ there exists an involutive isometry $\sigma_{p}$ that fixes
$p$, that is

$\sigma_{p}(p)=p\qquad\lparen\mathrm{d}\sigma_{p}\rparen_{0}=-\mathrm{Id}.$

As the name implies, one may prove that symmetric spaces are indeed locally
symmetric spaces.

For these manifolds we have the following result.

######  Proposition 4.14.

A symmetric space with bounded sectional curvature $\delta\leq\sec\leq\Delta$
is of $(\delta,\Delta,0)$-bounded geometry.

###### Proof.

Since it is a locally symmetric space, we have that $\nabla R=0$, so the first
order bound is clear. Since symmetric spaces are Riemannian homogeneous
spaces, for every two points there exists an isometry taking one to the other.
As such, the injectivity radius is constant throughout the manifold and thus,
positive. ∎

From this, we get that we just need to compute bounds on the sectional
curvature for these manifolds in order to give second order bounds for the
exponential map. Now, bounds on the sectional curvature of these manifolds are
well known. We give here some examples that are particularly useful in the
context of optimization.

######  Example 4.15 (The special orthogonal group).

The sectional curvature for a Lie group with a bi-invariant metric, after
identifying any pair of tangent vector with vectors in the Lie algebra, is
given for a pair of orthonormal vectors $X,Y\in\mathfrak{g}$

$\sec(X,Y)=\frac{1}{4}\lVert[X,Y]\rVert^{2}.$

In the special case of $\operatorname{SO}\lparen n\rparen$ for $n>2$, we have
that $\mathfrak{g}\cong\operatorname{Skew}\lparen n\rparen$. It is clear that
the sectional curvature is non-negative for any bi-invariant metric. In the
case of $\operatorname{SO}\lparen n\rparen$ this bound is tight.

Consider the bi-invariant metric given by the Frobenius norm—the metric
inherited from $\mathbb{R}^{n}$. For the upper bounds, if we work with a
matrix Lie group, it is enough to bound the norm of $XY-YX$ for matrices $X,Y$
of Frobenius norm $1$. In the general case, this inequality is called the
Böttcher-Wenzel inequality and it reads

$\lVert[X,Y]\rVert\leq 2\lVert X\rVert\lVert Y\rVert\mathrlap{\qquad\forall
X,Y\in\mathbb{R}^{n\times n}.}$

For a review of this inequality and a particularly clean proof see [43].

For the case of $\operatorname{SO}\lparen n\rparen$, that is, when $X,Y$ are
skew-symmetric, this inequality can be improved [8]

$\lVert[X,Y]\rVert\leq\lVert X\rVert\lVert Y\rVert\mathrlap{\qquad\forall
X,Y\in\operatorname{Skew}\lparen n\rparen.}$

For $n=3$ the constant can be further improved to $\frac{1}{2}$. These
constants are tight.

Wrapping all this together, we get that the bounds for
$\operatorname{SO}\lparen n\rparen$ with $n>2$ are given by

$\displaystyle 0\leq\langle\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen,\dot{\gamma}(r)$
$\displaystyle\leq\Big{\lparen}\frac{1}{r}-\frac{\sin\lparen
r\rparen}{r^{2}}\Big{\rparen}\lVert w\rVert^{2}$ $\displaystyle\qquad r<2\pi$
$\displaystyle\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen^{\perp}\rVert$ $\displaystyle\leq\frac{r}{9}\lVert w\rVert^{2}$
$\displaystyle\qquad r<2\sqrt{2}\pi.$

Note that the radius of definition of the normal part of the Hessian are
tight, as the conjugate radius of $\operatorname{SO}\lparen n\rparen$ is
exactly $2\sqrt{2}\pi$. This can be seeing by noting that the exponential of
matrices is not full rank at matrices with two eigenvalues that are $2\pi i$
apart.

In particular, we can give a bound on the full Hessian by bounding the
derivative of $\sqrt{\sigma^{2}+\rho^{2}}$ where $\sigma$ and $\rho$ are the
bounds on the tangential and normal part of the Hessian. This gives

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}(w_{1},w_{2})\rVert\leq\frac{2r}{9}\lVert
w_{1}\rVert\lVert w_{2}\rVert\qquad{r<2\pi}.$

We can see that Theorem 4.11 gives us a coarser bound, but on the other hand,
the bound is defined on a larger radius. In particular we get

$\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}(w_{1},w_{2})\rVert\leq\frac{r}{3}\lVert
w_{1}\rVert\lVert w_{2}\rVert\qquad{r<2\sqrt{2}\pi}.$

This radius is tight, as there are points that are $2\sqrt{2}\pi$ apart of any
given point which are conjugate to it. This can be seen by computing the
eigenvalues of the differential of the exponential map (see for example
[[]Theorem $D.2$]lezcano2019cheap).

Better bounds are possible by using a tighter version of the bounds on the
curvature tensor at the expense of having an uglier numeric constant. For
example, the $\frac{1}{3}$ constant can be improved this way to
$\frac{\sqrt{61}}{36}$.

These same ideas can be generalized to symmetric spaces. We will use this to
compute bounds for the Hessian of the exponential map in the Grassmannian.

######  Example 4.16 (Sectional Curvature of a Symmetric space).

Let $G/H$ be a symmetric space. Since a symmetric space is a Riemannian
homogeneous space, we just need to bound the sectional point at one point, as
the space has the same sectional curvature at every point. Denote by
$\mathfrak{m}=\mathfrak{h}^{\perp}\subseteq\mathfrak{g}$ the orthogonal
complement of the Lie algebra of $H$ with respect to the metric at the
identity in $G$. This complement to $\mathfrak{h}$ is not a Lie algebra itself
and in fact $[\mathfrak{m},\mathfrak{m}]\subseteq\mathfrak{h}$, as $G/H$ is a
symmetric space. This set $\mathfrak{m}$ may be identified isometrically with
the tangent space at $\pi(e)$, the projection of the identity element $e\in
G$, in other words, the map

$\lparen\mathrm{d}\pi\rparen_{e}|_{\mathfrak{m}}\colon\mathfrak{m}\to
T_{\pi(e)}G/H$

is a linear isometry. Via this identification we may treat vectors $X,Y\in
T_{\pi(e)}G/H$ as vectors
$\overline{X},\overline{Y}\in\mathfrak{m}\subseteq\mathfrak{g}$. Now, by
O’Neill’s formula [[]Chapter 7, Thm 47]oneill1966fundamental, the sectional
curvature for these manifolds has a particularly simple formula for
orthonormal vectors $X,Y\in T_{\pi(e)}G/H$

$\sec_{G/H}(X,Y)=\sec_{G}(\overline{X},\overline{Y})+\frac{3}{4}\lVert[\overline{X},\overline{Y}]\rVert^{2}$

where we have used that $[\mathfrak{m},\mathfrak{m}]\subseteq\mathfrak{h}$.

If $G$ is a Lie group with a bi-invariant metric we say that $G/H$ is a
_normal symmetric space_. For a normal metric, the sectional curvature
simplifies to

$\sec_{G/H}(X,Y)=\lVert[\overline{X},\overline{Y}]\rVert^{2}.$

Trough the study of symmetric spaces of non-compact type and their duality, it
is not difficult to prove that the sectional curvature on any symmetric space
is either the norm or minus the norm of the Lie bracket of a pair of vectors,
although we will not prove that as we will not need it.

######  Example 4.17 (The real Grassmannian).

The real Grassmannian as a symmetric space is given by the quotient
$\operatorname{Gr}\lparen n,k\rparen=\operatorname{SO}\lparen
n\rparen/\lparen\operatorname{SO}\lparen
k\rparen\times\operatorname{SO}\lparen n-k\rparen\rparen$, where the metric on
$\operatorname{SO}\lparen n\rparen$ is the bi-invariant metric generated by
the scalar product $\langle X,Y\rangle=\frac{1}{2}\operatorname{tr}\lparen
X^{\intercal}Y\rparen$ on the Lie algebra. For this symmetric space we have
that

$\displaystyle\mathfrak{h}=\operatorname{\mathfrak{so}}\lparen
k\rparen\oplus\operatorname{\mathfrak{so}}\lparen
n-k\rparen=\Big{\\{}\begin{pmatrix}B&0\\\
0&C\end{pmatrix}\>\Big{|}\>B\in\operatorname{Skew}\lparen
k\rparen,C\in\operatorname{Skew}\lparen n-k\rparen\Big{\\}}$
$\displaystyle\mathfrak{m}=\Big{\\{}\begin{pmatrix}0&A\\\
-A^{\intercal}&0\end{pmatrix}\>\Big{|}\>A\in\mathbb{R}^{n-k\times k}$

Note that the metric is chosen so that for $\overline{X}\in\mathfrak{m}$ we
have that $\lVert\overline{X}\rVert=\lVert X\rVert$, where
$X\in\mathbb{R}^{n-k\times k}$. This is so that this norm agrees with the
usual norm in the projective plane as $\operatorname{Gr}\lparen
n,1\rparen\cong\mathbb{R}P^{n}$.

Since the metric on $G=\operatorname{SO}\lparen n\rparen$ is bi-invariant, the
Grassmannian is a normal symmetric space, so the sectional curvature is given
by the norm of the commutator of elements in $\mathfrak{m}$ . Using the bound
on the norm of the Lie bracket for skew-symmetric matrices, we would get

$0\leq\sec(X,Y)\leq 4.$

As it can be seen by [[]Lemma 2.5]ge2014ddvv, these bounds are not tight. They
can be refined as announced in [[]Theorem 3a]wong1968sectional and proved in
[[]p.292]hildebrandt1980harmonic via a simple application of Cauchy-Schwarz as

$0\leq\sec(X,Y)\leq 2.$

Using these bounds, we get analogous bounds for the Hessian of the exponential
map for the Grassmannian,

$\displaystyle 0\leq\langle\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen,\dot{\gamma}(r)$
$\displaystyle\leq\Big{\lparen}\frac{1}{r}-\frac{\sin\lparen
2\sqrt{2}r\rparen}{2\sqrt{2}r^{2}}\Big{\rparen}\lVert w\rVert^{2}$
$\displaystyle\qquad r<\frac{\pi}{\sqrt{2}}$
$\displaystyle\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{rv}\lparen
w,w\rparen^{\perp}\rVert$ $\displaystyle\leq\frac{8r}{9}\lVert w\rVert^{2}$
$\displaystyle\qquad r<\pi.$

Note that for the real Grassmannian $\operatorname{r_{inj}}=\frac{\pi}{2}$
[40], so the radii of these equations should be enough for any practical
purposes. As in the case of $\operatorname{SO}\lparen n\rparen$, it is also
direct to get linear bounds on the full Hessian.

##  5 Convergence rates of dynamic trivializations

In this section, we use the bounds developed in the last two sections to prove
convergence of several descent algorithms on Riemannian manifolds.

###  5.1 Weakly Convex Optimization on Manifolds

We now have all the necessary tools to be able to complete the program stated
in the introduction. In particular, we can prove the convergence of the
dynamic trivialization framework for any stopping rule. We first recap the
results of the previous two sections in the following proposition.

######  Proposition 5.1.

Let $(M,\textsl{g})$ be a connected and complete Riemannian manifold of
$(\delta,\Delta,\Lambda)$-bounded geometry and let $f$ be an $\alpha$-weakly
convex function on it. Fix a point $p\in M$, and let $\mathcal{X}\subseteq
U_{p}$ be a subset of $M$ with
$\operatorname{diam}(\mathcal{X})=r\leq\pi_{\frac{\Delta+\delta}{2}}$, and at
least one critical point of $f$ in it. Denote the pullback of $\mathcal{X}$
under the exponential as
$\overline{\mathcal{X}}\coloneqq\exp^{-1}_{p}(\mathcal{X})\subseteq T_{p}M$.
Then, the map $f\circ\exp_{p}\colon\overline{\mathcal{X}}\to\mathbb{R}$ is
$\widehat{\alpha}_{r}$-weakly convex with constant

$\widehat{\alpha}_{r}=\alpha(C_{1,r}+C_{2,r})$

for

$\displaystyle
C_{1,r}=\max\Big{\\{}1,\mbox{{}\small$\displaystyle\frac{\operatorname{sn}_{\delta}(r)^{2}}{r^{2}}$}\Big{\\}}$
$\displaystyle
C_{2,r}=\mbox{{}\small$\displaystyle\frac{8}{3}$}\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}\lparen\Lambda\operatorname{sn}_{\delta}\big{\lparen}\mbox{{}\small$\displaystyle\frac{r}{2}$}\big{\rparen}^{2}+2\max\\{\lvert\Delta\rvert,\lvert\delta\rvert\\}\operatorname{sn}_{\delta}\lparen
r\rparen\rparen.$

###### Proof.

We will bound the Hessian of the map $f\circ\exp_{p}$ for an arbitrary $p\in
M$. By the Leibnitz rule, we have that

$\nabla\mathrm{d}\lparen
f\circ\exp_{p}\rparen=\nabla\lparen\mathrm{d}f\circ\mathrm{d}\exp_{p}\rparen=\nabla\mathrm{d}f\circ\mathrm{d}\exp_{p}+\mathrm{d}f\circ\nabla\mathrm{d}\exp_{p}.$

Taking norms, and since the left-hand side is a symmetric tensor, its maximum
value—i.e., its norm—it’s reached at a singular vector. As such, writing the
bound explicitly

$\displaystyle\lVert\nabla\mathrm{d}\lparen$ $\displaystyle
f\circ\exp_{p}\rparen\rVert_{\mathcal{X}}=$
$\displaystyle\max_{\begin{subarray}{c}v\in\exp^{-1}(\mathcal{X})\\\ w\in
T_{p}M,\ \lVert
w\rVert=1\end{subarray}}\Big{[}\lparen\nabla\mathrm{d}f\rparen_{\exp_{p}(v)}\lparen\lparen\mathrm{d}\exp_{p}\rparen_{v}(w),\lparen\mathrm{d}\exp_{p}\rparen_{v}(w)\rparen+\lparen\mathrm{d}f\rparen_{\exp_{p}(v)}\lparen\lparen\nabla\mathrm{d}\exp_{p}\rparen_{v}(w,w)\rparen\Big{]}.$

Since there exists a critical point of $f$ in $U_{p}$, and by the
$\alpha$-convexity, we have that $f$ is $\alpha r^{2}$-Lipschitz on
$\mathcal{X}$. Using this, the triangle inequality and Cauchy-Schwarz, we can
bound this quantity as

$\lVert\nabla\mathrm{d}\lparen
f\circ\exp_{p}\rparen\rVert_{\mathcal{X}}=\alpha\max_{\begin{subarray}{c}v\in\exp^{-1}(\mathcal{X})\\\
w\in T_{p}M,\ \lVert
w\rVert=1\end{subarray}}\lVert\lparen\mathrm{d}\exp_{p}\rparen_{v}(w)\rVert^{2}+\alpha
r^{2}\max_{\begin{subarray}{c}v\in\exp^{-1}(\mathcal{X})\\\ w\in T_{p}M,\
\lVert
w\rVert=1\end{subarray}}\lVert\lparen\nabla\mathrm{d}\exp_{p}\rparen_{v}(w,w)\rVert.$

From Theorems 3.12 and 4.11, we have that $C_{1,r},C_{2,r}$ are bounds for the
square of the norm of the differential and the norm Hessian of the Hessian of
the exponential respectively. ∎

After the heavy work of giving bounds on the Hessian of the pullback of a
function along the exponential map, we are in a position to prove convergence
rates for different instances of the dynamic trivialization framework in terms
of $\widehat{\alpha}_{r}$. We start with one of the simplest ones, namely
_static trivializations_. This is the algorithm that comes from choosing
$\texttt{stop}\equiv\texttt{False}$ in Algorithm 1. Equivalently, this is the
algorithm coming from solving

$\min_{v\in T_{p}M}f\lparen\exp_{p}(v)\rparen$

using gradient descent on $T_{p}M$.

The critical assumption in this result is that iterates remain bounded inside
a compact set $\mathcal{X}\subseteq T_{p}M$. This assumption is restrictive,
but it is standard in previous work [9, 61, 54, 57, 2].

######  Theorem 5.2 (Convergence of static trivializations).

Let $(M,\textsl{g})$ be a connected and complete Riemannian manifold of
$(\delta,\Delta,\Lambda)$-bounded geometry and let $f$ be an $\alpha$-weakly
convex function on it. Fix a point $p\in M$, and let $\mathcal{X}\subseteq
U_{p}$ be a subset of $M$ with $\operatorname{diam}(\mathcal{X})\leq r$, and
at least one critical point of $f$ in it. Consider Algorithm 1 with the
stopping rule $\texttt{stop}\equiv\texttt{False}$ and fixed step-size
$\eta_{i,k}=\mbox{{}\small$\displaystyle\frac{1}{\widehat{\alpha}_{r}}$}$
where $\widehat{\alpha}_{r}$ is as in Proposition 5.1. If all the iterates of
the method stay in $\mathcal{X}$, the method will find a point $v_{0,t}\in
T_{p}M$ such that $\lVert\nabla\lparen
f\circ\exp_{p}\rparen(v_{0,t})\rVert<\varepsilon$ in at most

$t\leq\Big{\lceil}\mbox{{}\small$\displaystyle\frac{2\widehat{\alpha}_{r}}{\varepsilon^{2}}$}\lparen
f(\exp_{p}(v_{0,0}))-f^{\ast}\rparen\Big{\rceil}$

steps, where $f^{\ast}$ is a lower-bound of $f$ on $\mathcal{X}$.

###### Proof.

We proved in Proposition 5.1 that the map $f\circ\exp_{p}$ is
$\widehat{\alpha}_{r}$-weakly convex on $\mathcal{X}$. This can be regarded as
a function on a Euclidean space, for which the convergence rate is well known
(see for example [48]). ∎

As we can see, the point $p\in M$ does not play any role in the proof of
Theorem 5.2, besides for the technical condition of the iterates being bounded
in $U_{p}$. Generalizing this assumption, we can prove at once the convergence
of the scheme of dynamic trivializations, for an arbitrary stopping rule.

######  Theorem 5.3 (Convergence of dynamic trivializations).

Let $(M,\textsl{g})$ be a connected and complete Riemannian manifold of
$(\delta,\Delta,\Lambda)$-bounded geometry and let $f$ be an $\alpha$-weakly
convex function on it. Assume that in the algorithm Algorithm 1 with an
arbitrary stopping rule stop, all the iterates $v_{i,k}$ are contained in sets
$\mathcal{X}_{i}\subseteq T_{p_{i}}M$ with
$\operatorname{diam}\lparen\mathcal{X}_{i}\rparen\leq r$, with at least one
critical point of $f$ in each of them. Then, for the choice
$\eta_{i,k}=\mbox{{}\small$\displaystyle\frac{1}{\widehat{\alpha}_{r}}$}$, the
algorithm will find a point $v_{i,k}\in T_{p}M$ such that $\lVert\nabla\lparen
f\circ\exp_{p_{i}}\rparen(v_{i,k})\rVert<\varepsilon$ in at most

$\Big{\lceil}\mbox{{}\small$\displaystyle\frac{2\widehat{\alpha}_{r}}{\varepsilon^{2}}$}\lparen
f(\exp_{p}(v_{0,0}))-f^{\ast}\rparen\Big{\rceil}$

steps, where $f^{\ast}$ is a lower-bound of $f$ on $\mathcal{X}$.

###### Proof.

Analogous to Theorem 5.2, as the proof does not depend on the pullback point
$p$. ∎

This theorem should can be considered as a building block to then be used in
specific examples to get actual convergence rates under weaker assumptions. In
plain words, this result asserts that if the algorithm is converging to a
critical point, it is fine to change the trivialization point, as the
exponential map will not distort the metric too much.

######  Remark (Practical considerations).

Of course, these bounds are worst-case bounds, but one can do better in
practice. As we have bounds for the distortion of the exponential map at every
point, we can choose a dynamic step-length that accounts for this. For
example, at a point with $\lVert v_{i,k}\rVert=s$, we could consider

$\eta_{i,k}=\frac{1}{\widehat{\alpha}_{s}}.$

In practice, this is not likely to give any improvement in a real-world
problem, as the weak-convexity constant of a function is itself an upper-bound
on the Hessian of the function itself, and the exponential map is not a well-
behaved map, as we have shown.

This is particularly true when using a dynamic stopping rule that controls how
much the norm of the gradient of the pullback deviates from the norm of the
gradient of the function. In this case, one could consider having a rule
similar to the one that we outlined in the introduction, but that does not
only account for the case in which the algorithm is converging to the cut
locus, which will mostly happen if the sectional curvature is positive, but
also accounts for the case when the gradient of the pullback is much larger
than that of the function, which might happen in cases of negative curvature

$\texttt{stop}\equiv\Big{\lparen}\frac{\lVert\nabla\lparen
f\circ\exp_{p_{i}}\rparen\lparen v_{i,k}\rparen\rVert}{\lVert\nabla f\lparen
x_{i,k}\rparen\rVert}<\varepsilon\Big{\rparen}\ \texttt{or}\
\Big{\lparen}\frac{\lVert\nabla\lparen f\circ\exp_{p_{i}}\rparen\lparen
v_{i,k}\rparen\rVert}{\lVert\nabla f\lparen
x_{i,k}\rparen\rVert}>\frac{1}{\varepsilon}\Big{\rparen}$

This rule can be read as _we change the trivialization point when we detect
that $\exp_{p}$ has deviated too much from being an isometry in the direction
of the gradient_. Of course, this rule can be adapted using the information
that one has a priori about the geometry of the manifold and, in particular,
its cut locus, conjugate locus or injectivity radius.

## Acknowledgements

We would like to thank Jaime Mendizabal and Prof. Vidit Nanda for checking the
proofs and giving very useful feedback on early versions of the work. We would
also like to thank Prof. Andrew Dancer, Prof. Raphael Hauser and Prof. Coralia
Cartis for useful feedback and pointers in the early stages of this project
and encouragement to follow the line of work that led to this final version.

This work was supported by an Oxford-James Martin Graduate Scholarship.

## References

  * [1] Naman Agarwal, Nicolas Boumal, Brian Bullins and Coralia Cartis “Adaptive regularization with cubics on manifolds” In _Mathematical Programming_ Springer, 2020
  * [2] Kwangjun Ahn and Suvrit Sra “From Nesterov’s Estimate Sequence to Riemannian Acceleration” In _arXiv preprint arXiv:2001.08876_ , 2020
  * [3] Vincent Arsigny, Olivier Commowick, Xavier Pennec and Nicholas Ayache “A log-Euclidean framework for statistics on diffeomorphisms” In _International Conference on Medical Image Computing and Computer-Assisted Intervention, MICCAI_ , 2006, pp. 924–931
  * [4] Vincent Arsigny, Pierre Fillard, Xavier Pennec and Nicholas Ayache “Geometric means in a novel vector space structure on symmetric positive-definite matrices” In _SIAM Journal on Matrix Analysis and Applications_ 29.1, 2006, pp. 328–347
  * [5] Miroslav Bačák “Convex analysis and optimization in Hadamard spaces” 22, De Gruyter Series in Nonlinear Analysis and Applications De Gruyter, Berlin, 2014 DOI: 10.1515/9783110361629
  * [6] Glaydston C. Bento, Orizon P. Ferreira and Jefferson G. Melo “Iteration-complexity of gradient, subgradient and proximal point methods on Riemannian manifolds” In _Journal of Optimization Theory and Applications_ 173.2, 2017, pp. 548–562 DOI: 10.1007/s10957-017-1093-4
  * [7] Marcel Berger “Sur quelques variétés riemanniennes suffisamment pincées” In _Bulletin de la Société Mathématique de France_ 88, 1960, pp. 57–71
  * [8] Anthony M. Bloch and Arieh Iserles “Commutators of skew-symmetric matrices” In _International Journal of Bifurcation and Chaos in Applied Sciences and Engineering_ 15.3, 2005, pp. 793–801 DOI: 10.1142/S0218127405012417
  * [9] Silvère Bonnabel “Stochastic gradient descent on Riemannian manifolds” In _IEEE Transactions on Automatic Control_ 58.9, 2013, pp. 2217–2229
  * [10] Silvère Bonnabel “Stochastic gradient descent on Riemannian manifolds” In _IEEE Transactions on Automatic Control_ 58.9, 2013, pp. 2217–2229
  * [11] Nicolas Boumal, Pierre-Antoine Absil and Coralia Cartis “Global rates of convergence for nonconvex optimization on manifolds” In _IMA Journal of Numerical Analysis_ 39.1, 2019, pp. 1–33
  * [12] Samuel Burer and Renato D.. Monteiro “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization” In _Mathematical Programming_ 95.2, Ser. B, 2003, pp. 329–357
  * [13] Samuel Burer and Renato D.. Monteiro “Local minima and convergence in low-rank semidefinite programming” In _Mathematical Programming_ 103.3, Ser. A, 2005, pp. 427–444
  * [14] Élie Cartan “Leçons sur la géométrie des espaces de Riemann” Gauthier-Villars, 1928
  * [15] Jeff Cheeger and David G. Ebin “Comparison theorems in Riemannian geometry” Revised reprint of the 1975 original AMS Chelsea Publishing, Providence, RI, 2008 DOI: 10.1090/chel/365
  * [16] Chris Criscitiello and Nicolas Boumal “Efficiently escaping saddle points on manifolds” In _Advances in Neural Information Processing Systems, NeurIPS_ , 2019, pp. 5985–5995
  * [17] Chris Criscitiello and Nicolas Boumal “An accelerated first-order method for non-convex optimization on manifolds” In _To Appear_ , 2020
  * [18] Alan Edelman, Tomás A. Arias and Steven Thomas Smith “The geometry of algorithms with orthogonality constraints” In _SIAM Journal on Matrix Analysis and Applications_ 20.2, 1998, pp. 303–353
  * [19] Jürgen Eichhorn “The boundedness of connection coefficients and their derivatives” In _Mathematische Nachrichten_ 152, 1991, pp. 145–158 DOI: 10.1002/mana.19911520113
  * [20] Jaap Eldering “Persistence of noncompact normally hyperbolic invariant manifolds in bounded geometry” In _Comptes Rendus Mathématique. Académie des Sciences. Paris_ 350.11-12, 2012, pp. 617–620 DOI: 10.1016/j.crma.2012.06.009
  * [21] Jost-Hinrich Eschenburg “Comparison Theorems in Riemannian Geometry”, Dipartimento di Matematica: Lecture notes series Università delli studi di Trento, Dipartimento di matematica, 1994
  * [22] Jost-Hinrich Eschenburg and Ernst Heintze “Comparison theory for Riccati equations” In _Manuscripta Mathematica_ 68.2, 1990, pp. 209–214
  * [23] Luca Falorsi, Pim Haan, Tim R. Davidson and Patrick Forré “Reparameterizing distributions on Lie groups” In _International Conference on Artificial Intelligence and Statistics, AISTATS_ , 2019, pp. 3244–3253
  * [24] O.. Ferreira, M.. Louzeiro and L.. Prudente “Iteration-complexity of the subgradient method on Riemannian manifolds with lower bounded curvature” In _Optimization_ 68.4, 2019, pp. 713–729
  * [25] O.. Ferreira and B.. Svaiter “Kantorovich’s theorem on Newton’s method in Riemannian manifolds” In _Journal of Complexity_ 18.1, 2002, pp. 304–329 DOI: 10.1006/jcom.2001.0582
  * [26] Jianquan Ge “DDVV-type inequality for skew-symmetric matrices and Simons-type inequality for Riemannian submersions” In _Advances in Mathematics_ 251, 2014, pp. 62–86 DOI: 10.1016/j.aim.2013.10.010
  * [27] R.. Greene “Complete metrics of bounded curvature on noncompact manifolds” In _Archiv der Mathematik_ 31.1, 1978, pp. 89–95 DOI: 10.1007/BF01226419
  * [28] Mikhael Gromov “Structures métriques pour les variétés riemanniennes”, Textes Mathématiques CEDIC, Paris, 1981
  * [29] Sigurður Helgason “Differential geometry, Lie groups, and symmetric spaces” Academic press, 1978
  * [30] S. Hildebrandt, J. Jost and K.-O. Widman “Harmonic mappings and minimal submanifolds” In _Inventiones Mathematicae_ 62.2, 1980, pp. 269–298 DOI: 10.1007/BF01389161
  * [31] Wen Huang, Kyle A. Gallivan and Pierre-Antoine Absil “A Broyden Class of Quasi-Newton Methods for Riemannian Optimization” In _SIAM Journal on Optimization_ 25.3, 2015, pp. 1660–1685 DOI: 10.1137/140955483
  * [32] Arieh Iserles, Hans Z. Munthe-Kaas, Syvert P. Nørsett and Antonella Zanna “Lie-group methods” In _Acta numerica_ 9, 2000, pp. 215–365
  * [33] Arieh Iserles and Syvert P. Nørsett “On the solution of linear differential equations in Lie groups” In _Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ 357.1754, 1999, pp. 983–1019
  * [34] Jin-ichi Itoh and Minoru Tanaka “The dimension of a cut locus on a smooth Riemannian manifold” In _Tohoku Mathematical Journal_ 50.4, 1998, pp. 571–575
  * [35] Jürgen Jost “Riemannian geometry and geometric analysis” Springer, 2017
  * [36] Hermann Karcher “A short proof of Berger’s curvature tensor estimates” In _Proceedings of the American Mathematical Society_ 26, 1970, pp. 642–644
  * [37] Hiroyuki Kasai, Pratik Jawanpuria and Bamdev Mishra “Riemannian adaptive stochastic gradient algorithms on matrix manifolds” In _International Conference on Machine Learning, ICML_ , 2019, pp. 3262–3271 URL: http://proceedings.mlr.press/v97/kasai19a.html
  * [38] Hiroyuki Kasai, Hiroyuki Sato and Bamdev Mishra “Riemannian stochastic recursive gradient algorithm with retraction and vector transport and its convergence analysis” In _International Conference on Machine Learning, ICML_ , 2018, pp. 2521–2529
  * [39] Helmut Kaul “Schranken für die Christoffelsymbole” In _Manuscripta Mathematica_ 19.3, 1976, pp. 261–273
  * [40] SE Kozlov “Geometry of real Grassmann manifolds. Part III” In _Journal of Mathematical Sciences_ 100.3 Springer, 2000, pp. 2254–2268
  * [41] Mario Lezcano-Casado “Trivializations for gradient-based optimization on manifolds” In _Advances in Neural Information Processing Systems, NeurIPS_ , 2019, pp. 9154–9164
  * [42] Mario Lezcano-Casado and David Martínez-Rubio “Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group” In _International Conference on Machine Learning, ICML_ , 2019, pp. 3794–3803