text:
\begin{document}

\maketitle

\begin{abstract}
Optimization is central to problems across machine learning, estimation, automatic control, and many other areas. Convex optimization problems are an important subclass, admitting global solutions and mature solution methods. Recent attention has focused on geodesic convexity, or g-convexity, which allows non-convex problems to be transformed into convex ones posed over Riemannian manifolds. The main contribution of this paper is to provide a bridge between g-convexity and contraction analysis of
nonlinear dynamical systems. Specifically, we show the equivalence between strong geodesic convexity and contraction of natural gradient flows. Given this result, existing tools for the analysis and synthesis of nonlinear contracting systems can be considered to both discover and efficiently optimize g-convex functions. In addition, the contraction perspective allows one to easily apply results to time-varying optimization settings, and to recursively build large optimization structures out of simpler elements.
\end{abstract}

\section{Introduction}

This paper considers geodesically-convex optimization \cite{RAPCSAK91} through the lens of nonlinear contraction analysis.
Many problems in learning, estimation, automatic control and other areas can be phrased in terms of optimization problems posed over Riemannian manifolds~\cite{Absil08,Sra_First_Order_Methods,Bonnabel12}. Transitioning these problems away from a conventional Euclidean formulation has multiple advantages. Nonlinear constraints that may define a manifold can be eliminated by considering unconstrained formulations on the manifold directly. Second, functions that are not convex in a Euclidean sense may be {\em geodesically} convex along shortest paths on the manifold~\cite{udriste}.

    This second advantage presents opportunity to reformulate non-convex optimization problems over $\mathbb{R}^n$ as convex ones over a Riemannian manifold. This reformulation is accomplished not through any explicit coordinate change, but simply by endowing $\Rn$ with a suitable Riemannian metric. Principled search over metrics can be accomplished computationally and represents additional generality beyond considering nonlinear changes of coordinates. In effect, endowing $\mathbb{R}^n$ with a Riemannian metric may be viewed as providing a differential change of coordinates~\cite{Slotine98}. While an explicit coordinate change can be very effective if it fundamentally derives from the structure or the physics of the system, a differential coordinate change is not necessarily integrable and thus is considerably more general.

    A similar approach is at the heart of nonlinear contraction analysis. Contraction theory~\cite{Slotine98} allows the stability of nonlinear non-autonomous systems to be characterized through linear time-varying dynamics describing the propagation of infinitesimally-small displacements along the systems' flow. The existence of a Riemannian metric that contracts these virtual displacements (i.e., elements in the tangent space) is necessary and sufficient for exponential convergence of any pair of trajectories. Contraction naturally yields methods for constructing stable systems of systems, including synchronization phenomena and consensus as well as other key building-blocks that allow the construction of large contracting systems out of simpler elements.

    The main contribution of this paper is to provide a precise link between the geodesic convexity of functions on $\Rn$ and the contraction of their gradient flows. Section~\ref{sec:Convex} introduces key definitions and our main theoretical result. Specifically, we analyze continuous time gradient descent on a manifold, which takes the form of the natural gradient~\cite{Amari98b} in coordinates. Analysis in continuous time optimization is limited, in part, by the fact that dynamics can be arbitrarily sped up to achieve any convergence rate without regards for how it may affect a discrete implementation. However, a continuous perspective has yielded insight on important phenomena such as in the analysis~\cite{su2014differential}, discrete implementations \cite{Zhang18}, and extensions \cite{wibisono2016variational} of Nesterov's accelerated gradient descent method \cite{nesterov1998introductory}. It has also enabled analysis of primal-dual algorithms~\cite{Kostya}, where an absolute time reference is obtained by introducing additional fast dynamics or delays using a singular perturbation framework.

Following the main result linking g-convexity and contraction in Section~\ref{sec:Convex}, Section \ref{GPD_opt} presents an extension for the analysis of primal-dual type dynamics that appear in mixed convex/concave saddle systems. Section \ref{sec:application} then discusses how these insights can be used to design larger optimization structures out of simpler ones by leveraging combination properties of contracting systems. Section \ref{sec:Conclusions} provides an outlook on potential future advances that may stem from these connections.

\section{Definitions and Main Result}
\label{sec:Convex}
We first recall basic definitions and facts on convex
and g-convex functions, and on contraction analysis
of nonlinear dynamical systems.  These definitions
provide the context for our main result, Theorem \ref{th_main}, that establishes equivalence between g-convexity and contraction of gradient flows.

\subsection{Convex Optimization}

\begin{definition}[Strong Convexity]
\label{def:strongConv}
A function $f \in \mathcal{C}^2( \mathbb{R}^n , \mathbb{R})$ is $\alpha$-strongly convex with $\alpha>0$ if its Hessian matrix $\pdd{f}{\x}{\x}$ satisfies the matrix inequality
\[
\pdd{f}{\x}{\x} \succeq \alpha\, \mathbf{I} \quad \quad \forall \x \in \mathbb{R}^n
\]

\end{definition}
As its name suggests, a function that is strongly convex is convex in the usual sense, while the converse is not always true. From a dynamic systems perspective, strong convexity provides exponential convergence of gradient flows:
\begin{proposition}
\label{prop:euc_grad_desc}
\ If a function $f \in \mathcal{C}^2( \mathbb{R}^n , \mathbb{R})$ is $\alpha$-strongly convex, then its gradient system
\begin{equation}
\dot{\x} = - \pd{f}{\x}
\label{eq:grad_descent}
\end{equation}
converges to the global minimum exponentially with rate $\alpha$.
\end{proposition}

\subsection{g-Convex Optimization}

Geodesic convexity \cite{RAPCSAK91} generalizes conventional notions of convexity to the case where the domain of a function is viewed as a Riemannian manifold. A special case occurs in geometric programming (GP)~\cite{Boyd_GeomProgram}. In GP, a non-convex problem over positive variables $\{x_i\}_{i=1}^N$ can be transformed into a convex problem by a change of variables $y_i = {\rm log}(x_i)$. Alternately GP can be formulated over the positive reals viewed as a Riemannian manifold by measuring differential length elements $ds$ in a relative sense
\begin{equation}
ds^2 = \sum_{i=1}^N \left( \frac{dx_i}{x_i} \right)^2 = \sum_{i=1}^N dy_i^2
\label{eq:GP_Gconvex}
\end{equation}
Geodesically-convex optimization generalizes this transformation strategy to a broader class of problems~\cite{Sra_First_Order_Methods}. However, beyond special cases (see, e.g., \cite{SRA_Conic_Optimization}), generative procedures remain lacking to formulate g-convex optimization problems or recognize g-convexity.

To introduce g-convexity more formally, consider a function $f:\mathbb{R}^n \rightarrow \mathbb{R}$ and a positive definite metric $\M(\x):\Rn \rightarrow \R^{n\times n}$. We note that geodesic convexity of $f$ is not an intrinsic property of the function itself, but rather is a property of $f$ defined on the Riemannian manifold $(\Rn, \M)$.
\begin{definition}[g-Strong Convexity \cite{udriste}]
A function $f \in \mathcal{C}^2( \mathbb{R}^n , \mathbb{R})$ is said to be geodesically $\alpha$-strongly convex (with $\alpha>0$) in a symmetric positive definite metric $\M(\x)$ if its Riemannian Hessian matrix $\H(\x)$ satisfies:
\begin{equation}
\H(\x) \succeq \alpha\, \M(\x)\quad \quad \forall \x \in \mathbb{R}^n
\label{eq:Hessian}
\end{equation}
The elements of the Riemannian Hessian are given as
\begin{equation}
H_{ij} = \pdd{f}{i}{j} - \Gamma_{ij}^k \, \pd{f}{k}
\label{eq:riemhessian}
\end{equation}
where $\pdd{f}{i}{j} = \frac{\partial^2 f}{\partial x_i \partial x_j}$ provide the elements of the conventional (Euclidean) Hessian and $\Gamma_{ij}^k$ denotes the Christoffel symbols of the second kind
\[
\Gamma_{ij}^m = \frac{1}{2} \sum_{k=1}^n \left[ M^{mk} \left( \pd{M_{ik}}{j} + \pd{M_{jk}}{i} - \pd{M_{ij}}{k} \right)  \right]
\]
with $M^{ij} = (\M^{-1})_{ij}$. The function $f$ is g-convex when \eqref{eq:Hessian} holds with $\alpha=0$.
\end{definition}

The Riemannian Hessian generalizes the notion of the Hessian from a Euclidean context and captures the curvature of $f$ along geodesics. Likewise, the natural gradient generalizes the notion of a Euclidean gradient to the Riemannian context in the following sense.

\begin{definition}[Natural Gradient \cite{Amari98b}]
Consider $\Rn$ equipped with a Riemannian metric $\M(\x)$. The {\em natural} gradient of the function $f \in \mathcal{C}^2(\Rn,\R)$ is the direction of steepest ascent on the manifold and is given by $ \ \M(\x)^{-1} \pd{f}{\x} \ $.

\end{definition}

\begin{remark} From a differential geometric viewpoint, the first covariant derivative of $f$ is a covector field  given in coordinates by $\pd{f}{\x}$, while
the natural gradient is a vector field given in coordinates by $\M(\x)^{-1} \pd{f}{\x}$ ~\cite{Amari98b}. In a Euclidean context, where $\M(\x)$ is identity, this distinction between covariant (covector) and contravariant (vector) representations of the gradient is immaterial.

Similarly, the Riemannian Hessian $\vH$ represents in coordinates the second covariant derivative of $f$.

\end{remark}

When $\M(\x)$ is the identity metric, geodesic $\alpha$-strong convexity naturally coincides with the definition of $\alpha$-strong convexity in Definition~\ref{def:strongConv}. Toward providing a parallel to Proposition \ref{prop:euc_grad_desc}, we consider stability analysis in the Riemannian context through the lens of contraction theory.

\begin{definition}[Contraction Metric~\cite{Slotine98}]
\label{def:contraction}
A system $\dot{\x} = \h(\x,t)$ is said to be contracting at rate $\alpha>0$ with respect to
a symmetric positive definite metric $\M(\x) : \Rn  \rightarrow \R^{n\times n}$,  if for all  $t \ge 0$ and all $\x \in \Rn $,
\begin{equation}
\dot{\M} + \A\T\, \M + \M\, \A \preceq -2 \alpha \M
\label{eq:contraction_condition}
\end{equation}
where $\A(\x,t) = \pd{\h}{\x}$ is the system Jacobian and $\dot{\M}= \sum_i \left( \pd{\M}{\xsub{i}} \right) h_i(\x)$.
The system is said to be semi-contracting with respect to $\M(\x)$ when \eqref{eq:contraction_condition} holds with $\alpha=0$.
\end{definition}

As a key property, given an $\alpha$-contracting system and an arbitrary pair of initial conditions $\x_1(0)$ and $\x_2(0)$, the solutions $\x_1(t)$ and $\x_2(t)$  converge to one another exponentially
\begin{equation}\label{d_M}
d_{\mathcal{M}}( \x_1(t), \x_2(t)) \ \le \ {\rm e}^{-\alpha t} \ d_{\mathcal{M}}(\x_1(0), \x_2(0)) \ \
\end{equation}
where $d_{\mathcal{M}}(\cdot,\cdot)$ denotes the geodesic distance on the Riemannian manifold $\mathcal{M} = (\Rn, \M(\x))$. This property can be shown by considering the evolution of differential displacements $\delta \x$, which describe the evolution of nearby trajectories and coincide with the notion of virtual displacements in Lagrangian mechanics. More precisely, letting $\x(t;\x_0, t_0)$ denote the solution of $\dot{\x} = \h(\x,t)$ from initial condition $\x(t_0) = \x_0$, differential displacements satisfy
\[
\delta \x(t) = \frac{ \partial \x(t; \x_0, t_0) }{\partial \x_0} \ \delta \x(t_0)
\]

Property~\eqref{d_M} follows from the evolution of the squared length of these differential displacements~\cite{Slotine98}, which verifies
\begin{equation}\label{ddt_z^2}
    \ddt (\delta \x^T \M \delta \x) \le - 2 \alpha (\delta \x^T \M \delta \x)
\end{equation}
Furthermore, if a system is $\alpha$-contracting in a metric $\M(\x)$ that satisfies $\M(\x)\succ \beta \mathbf{I}$  uniformly
for some constant $\beta>0$, then any two solutions verify
\[
\| \x_1(t) - \x_2(t) \| \ \le \ \frac{1}{\sqrt{\beta}} \  {\rm e}^{-\alpha t} d_{\mathcal{M}}(\x_1(0), \x_2(0))
\]

\begin{example}
Consider an $\alpha$-strongly convex function $f$ and its associated gradient descent system~\eqref{eq:grad_descent}. Since $f$ is strongly convex, it has a unique global minimum $\x^*$, which is a stationary point of \eqref{eq:grad_descent}. It can be verified that the gradient descent dynamics of $f$ are contracting in an identity metric $\M = {\bf I}$ with rate $\alpha$. Since geodesic distances are just Euclidean distances in this metric,~\eqref{d_M} immediately implies that
\[
\forall t \ge 0, \ \ \|\x(t) - \x^*\| \ \le \ {\rm e}^{-\alpha t} \ \| \x(0) - \x^* \|
\]
thus proving Proposition~\ref{prop:euc_grad_desc}.
\end{example}

The following theorem uses the natural gradient to provide a critical link between g-convexity and contraction analysis of nonlinear systems, mirroring Proposition~\ref{prop:euc_grad_desc} within the Riemannian context.

\begin{theorem}\label{th_main} Consider a function $f(\x,t) \in \mathcal{C}^2(\mathbb{R}^n\times \mathbb{R}, \mathbb{R})$, a symmetric positive definite metric $\M(\x): \mathbb{R}^n \rightarrow \mathbb{R}^{n \times n}$, and the natural gradient system \cite{Amari98b}

\begin{equation}
\dot{\x} = \h(\x,t) =  -\M(\x)^{-1}\, \pd{f(\x,t)}{\x} \label{eq:NatGrad}
\end{equation}
Then, $f$ is $\alpha$-strongly g-convex in the metric $\M$ for each $t$  if and only if \eqref{eq:NatGrad} is contracting with rate $\alpha$ in the metric $\M$. More specifically, the Riemannian Hessian verifies
\begin{equation}
\vH = -\frac{1}{2}\left( \M \left( \pd{\h}{\x}\right)  + \left( \pd{\h}{\x}\right)\T \M + \dot{\M} \right)
\label{eq:H}
\end{equation}
\end{theorem}

\begin{proof}
See Appendix 1 for a self-contained proof.
\end{proof}

\begin{remark}
The result can also be derived as a special case of contraction analysis for complex Hamilton-Jacobi dynamics~\cite{Loh09}. A reorganization of \eqref{eq:NatGrad} as
\[
\M(\x) \dot{\x} = - \pd{f(\x,t)}{\x}
\]
may be recognized as the generalized momentum being the negative covariant gradient within a Hamiltonian context.
\end{remark}

\begin{remark}
While the above theorem applies to $\alpha$-strong convexity, the link between the Riemannian Hessian and the contraction condition~\eqref{eq:contraction_condition} also provides immediate equivalence between g-convexity of a function and semi-contraction of its natural gradient dynamics.
\end{remark}

\begin{remark}
The above result provides an alternate
way to compute the geodesic Hessian $\H$, and, as expected, leaves it invariant when the metric $\M$ is scaled by a strictly positive constant. Because of
the structure of the natural gradient dynamics, scaling
$\M$ is akin to scaling time and implies inversely scaling the contraction rate $\alpha$, consistently with~\eqref{eq:Hessian}.\\
By contrast, note that given a \emph{fixed} dynamics $\h$, the contraction metric analyzing it can always
be arbitrarily scaled while leaving the contraction rate unchanged.
\end{remark}

\begin{remark}
\label{rem:robust}
\renewcommand{\bd}{{\mathbf d}}
Note that contraction also provides robustness. Consider perturbed dynamics $\dot{\x} = \h(\x,t) + \bd(t)$ with $\sqrt{ \bd(t)\T \M(\x) \bd(t) } < R$ uniformly. If the dynamics are contracting with rate $\lambda$, then all trajectories contract to a geodesic ball of radius $R/\lambda$ \cite{Slotine98}.
\end{remark}

Let us illustrate Theorem~\ref{th_main} using the classical nonconvex Rosenbrock function:
\begin{equation}
f(\x) = 100 (x_1^2- x_2)^2 + (x_1-1)^2
\label{eq:rosenbrock}
\end{equation}
This function has a unique global optimum at $\x^*=[1,1]\T$, which is located along a long, shallow, parabolic-shaped valley. As a result, it exhibits poor scaling and is frequently used as a test problem for optimization.

\begin{example}
Consider the Rosenbrock function \eqref{eq:rosenbrock} and the metric
\[
\M(\x) = \begin{bmatrix} 400 x_1^2 + 1& -200 x_1 \\
        -200 x_1 &      100
\end{bmatrix}
\]
The metric $\M(\x)$ satisfies ${\rm tr}(\M(\x))= 400 x_1^2 + 101>0$ and ${\rm det}(\M(\x)) = 100>0$, and thus $\M(\x)\succ0$. Note that $\M(\x)$ is not simply the Hessian of $f(\x)$. The natural gradient dynamics follows
\[
\dot{\x} = \h(\x)= - \M(\x)^{-1} (\pd{f}{\x}) = -2 \begin{bmatrix} x_1 -1 \\ x_1^2 - 2 x_1 +  x_2 \end{bmatrix}
\]
It can be verified algebraically that
\[
\M \left( \pd{\h}{\x} \right) + \left(\pd{\h}{\x}\right)\T \M + \dot{\M} = -4 \M
\]
which shows that natural gradient descent is contracting with rate $\alpha=2$. This implies that the natural gradient dynamics satisfy
\[
d_{\mathcal{M}}(\x(t), \x^*) \ \le \ {\rm e}^{-2t} \ d_{\mathcal{M}}(\x(0), \x^*)
\]
where $\x^* = [1,1]\T$. Equivalently, the Rosenbrock function is geodesically $\alpha$-strongly convex with $\alpha = 2$.
\end{example}

The Rosenbrock metric $\M(\x)$ can be viewed as following from a differential change of variables
\[
\delt{z} = \boldsymbol{\Theta}(\x) \delt{x} = \begin{bmatrix} 20 x_1 & -10 \\ 1 & 0 \end{bmatrix} \delt{x}
\]
where $\M=\boldsymbol{\Theta}\T\boldsymbol{\Theta}$ yields $\ \delt{x}\T\M\delt{x} = \delt{z}\T \delt{z}$. This differential change of variables is integrable, so that g-convexity of the Rosenbrock can be shown using the explicit nonlinear coordinate change $z_1 = 10 x^2_1 - 10 x_2$ and $z_2 = x_1-1$ that provides $f = z_1^2 + z_2^2$. Although in this case the coordinate transform is integrable,  the general freedom to consider differential changes of coordinates provides additional flexibility and generality to both contraction analysis and g-convexity.

\begin{example}
Geodesically-convex optimization can also be used to carry out manifold-constrained optimization in an unconstrained fashion via recasting problems over a Riemannian manifold directly \cite{Absil08,SRA_Conic_Optimization}.

Consider for instance optimization over the set $\mathbb{S}^n_+$ of $n\times n$ positive definite matrices, and specifically the problem of finding the Karcher mean of $m$ matrices $\A_i \in \mathbb{S}^n_+$ \cite{Sra_First_Order_Methods}, which minimizes the objective function
\[
f(\X) = \frac{1}{2}\sum_{i=1}^m\| {\rm log}(\A_i^{-1}\,\X) \|_F^2
\]
where $\|\A\|_F=\sqrt{{\rm tr}(\A\T \A)}$ denotes the Frobenius norm of a matrix $ \A$. The function $f(\X)$ is $m$-strongly convex~\cite{Sra_First_Order_Methods} on $\mathbb{S}^n_+ \ $in the metric which measures square differential displacements as
\begin{equation}
ds^2 = {\rm tr}\left( \left( \delta \X \ \X^{-1} \right)^2 \right)
\label{eq:PSD_Riem}
\end{equation}
This metric generalizes the GP case \eqref{eq:GP_Gconvex}, and coincides with the Hessian of the log barrier $-{\rm logdet}(\X)$ \cite{BoydVanberghe04}. The gradient of $f(\X)$ can be written
\[
\pd{f}{\X} = \sum_{i=1}^m {\rm log}(\A_i^{-1}\,\X )\,  \X^{-1}
\]
and accordingly the natural gradient can be shown to satisfy
\[
\sum_{i=1}^m \ \X \, {\rm log}(\A_i^{-1}\,\X) = \X \left( \pd{f}{\X} \right) \X
\]
From Theorem \ref{th_main},  any trajectory with arbitrary initial condition $\X(0) \in \mathbb{S}^n_+$ will remain within $\mathbb{S}^n_+$ under the natural gradient descent dynamics
\[
\dot{\X} = - \sum_{i=1}^m  \ \X \, {\rm log}(\A_i^{-1}\,\X)
\]
since (intuitively) the Riemannian metric \eqref{eq:PSD_Riem} makes positive semi-definite matrices an infinite distance away from any positive definite one, and contraction of the natural gradient dynamics ensures that geodesic distances decrease exponentially.
\end{example}

\begin{example}
An approximation to the Riemmanian distance of two positive definite (PD) matrices on the PD cone is given by the Bregman LogDet divergence
\[
d(\A||\X) = {\rm logdet}(\A^{-1} \X) + {\rm tr}(\X^{-1} \A)
\]
The metric is convex in its first argument, and can be shown to be geodesically convex in the second. We illustrate the connection with contraction to show this property. Note that
\[
\pd{d(\A||\X)}{\X} = \X^{-1} - \X^{-1} \A \X^{-1}
\]
so that the natural gradient descent dynamics are simply
\[
\dot{\X} = - \X \left[ \pd{d(\A||\X)}{\X} \right] \X = -\X + \A
\]
with differential dynamics
\[
\delta \dot{\X} = - \delta \X
\]
Considering the rate of change in length of these differential displacements
\[
\ddt {\rm tr}\left( (\X^{-1} \, \delta \X)^2 ) \right) = -2{\rm tr}(\X^{-1} \A (\X^{-1} \delta \X^{-1})^2 )
\]
and defining the differential change of variables $\delta \Z = \X^{-\frac{1}{2}}\, \delta \X \, \X^{-\frac{1}{2}}$, one has $\tr{ \delta \Z^2} = {\rm tr}\left( (\X^{-1} \, \delta \X)^2 ) \right) $ and
\[
\ddt \tr{ \delta \Z^2} = -2 \tr{ \delta \Z \, \X^{-\frac{1}{2}} \A \X^{-\frac{1}{2}} \,\delta \Z} <0
\]
for all $\delta \Z\ne \mathbf{0}$. Hence, considering only the second argument to LogDet divergence, its Riemannian Hessian is positive definite, thus proving g-convexity.
\end{example}

\newcommand{\p}{\mathbf{p}}
\newcommand{\q}{\mathbf{q}}
\begin{example}

More generally, given a convex set $\mathcal{X}$ and a convex function $f\in \mathcal{C}^3( \mathcal{X}, \mathcal{R})$ , the Bregman divergence of $f$ is given by
\[
d_f(\p||\q) = f(\p) - f(\q) - \pd{f}{\x}(\q)\T (\p-\q)
\]
Consider the metric induced by the Hessian of $f$, $M_{ij}(\x) = \pd{f}{ij}$. The Riemannian Hessian of the divergence $d_f$ in its second argument has elements
\[
H_{ij} = \pd{f}{ij}(\q) - \frac{1}{2} \sum_k \pd{f}{ijk}(\q) \cdot (p_k-q_k)
\]
or more compactly, using tensorial multiplication with the tensor of third derivatives,
\[
\H = \pd{f}{\x\x}(\q) - \frac{1}{2} \pd{f}{\x\x\x}(\q) \cdot(\p-\q)
\]
This general form can be applied to understand the previous result in a more general setting.
Let $\mathcal{X}$ be the interior of the probability simplex $\{ \x \in \Rn~|~\sum x_i = 1, x_i > 0\}$ , and consider the convex function $f(\x) = \sum_i x_i {\rm log}(x_i)$ . The associated Bregman divergence
\[
d_f(\p||\q) = \sum_i \left[{\log}\left( \frac{p_i}{q_i} \right) - p_i + q_i\right]
\]
is the Kullback-Leibler (KL) divergence between the distributions. It is g-convex over $\q$ in the diagonal metric $M_{ii}(\x) = 1/x_i$ , with diagonal Riemannian Hessian
\[
H_{ii} = \frac{p_i+q_i}{2 q_i^2}
\]
The previous g-convexity result, on the divergence between positive definite matrices, can be equivalently viewed as the KL divergence on the statistical manifold of zero-centered multivariate Gaussians~\cite{Calvo90}.
\end{example}

Note that in the autonomous case, equations
governing the differential displacement follow
\begin{equation}\label{delta_x}
\ddt \ \delta \x \ = \ \partial_{\bf x} \h(\x) \ \delta \x
\end{equation}
which has a similar structure to the time evolution of $\h(\x)$
\begin{equation}\label{ddt_h}
\ddt \ \h(\x) \ = \ \left( \partial_{\bf x} \h(\x)\right) \ \h(\x)
\end{equation}
Thus, for natural gradient descent of an $\alpha$-strong g-convex function $f(\x)$, the same algebra leading to~\eqref{ddt_z^2} also gives
\[
\ddt (\h^T \M \h) = -2 \h^T \H \h \le - 2 \alpha (\h^T \M \h)
\]
so that
\[
V(\x) = \h(\x)^T \M(\x) \h(\x) = (\partial_{\bf x} f)^T \M(\x)^{-1} \partial_{\bf x}f
\]
can be viewed as an exponentially converging
Lyapunov function, with global minimum $V=0$ at the unique
minimum of $f(\x)$. Of course, \eqref{delta_x} remains
valid for non-autonomous systems as well, while~\eqref{ddt_h} does not.

\subsection{Non-autonomous systems}

In our optimization context, the fact that contraction analysis is directly applicable to non-autonomous systems can be exploited in a variety of ways. As we shall detail later, a key aspect is that it allows feedback combinations or hierarchies of contracting modules
to be exploited to address more elaborate optimization problems or architectures. Also, it makes the construction of \emph{virtual} systems~\cite{wei} possible to potentially extend results beyond strict natural-gradient descent.

\begin{remark}  The natural gradient $\ \M^{-1}(\x)\ \pd{f}{\x}(\x,t)$ represents the direction of steepest ascent on the manifold at any given time. With this in mind, Remark~\ref{rem:robust} on robustness enables convergence analysis for natural gradient descent within time-varying optimization contexts~\cite{Kostya}. Let $\x^*(t)$ denote the optimum of a time-varying $\alpha$-strongly g-convex function. If $\sqrt{\dot{\x}^*(t) \M(\x) \dot{\x}^*(t)}<R$, then the natural gradient will track $\dot{\x}^*(t)$ with accuracy $R/\alpha$ after exponential transient.
\end{remark}

\subsection{Semi-Contracting Case}

Consider the case when the Riemannian Hessian is only positive semi-definite. The natural gradient dynamics are then semi-contracting, i.e., the geodesic distance between any two solutions is non-increasing. We still have, from~\eqref{eq:H},
\[
\ddt\left( \delta \x\T \M( \x(t) )\, \delta \x \right)  \ = \ - \ 2 \ \delta \x\T  \H(\x(t),t ) \, \delta \x \ \le \ 0
\]
In the case that $\M$ is uniformly positive definite and one trajectory of the system is known to be bounded, this property then ensures that all trajectories are bounded. Then, from Barbalat's lemma and the symmetry of $\H$,
\[\H(\x,t)\, \delta \x \ \rightarrow \ 0\ \ {\rm~as}\ \ t\rightarrow \infty
\]
Thus, the difference between the solutions of any two nearby trajectories must approach the nullspace of the Riemannian Hessian.

This property allows us to derive the following theorem from Rapcsak~\cite{RAPCSAK91} using purely dynamic systems tools.

\begin{theorem}
\label{thm:semi}
Consider a function $f(\x) \in \mathcal{C}^2(\mathbb{R}^n\times \R, \R)$ and a uniformly symmetric positive definite metric $\M(\x):\Rn \rightarrow \R^{n\times n}$. If the Riemannian Hessian of $f$ in the metric $\M$ is positive semi-definite, then any stationary point of $f$ is a global optimum.
\end{theorem}
\begin{proof}
Consider two stationary points $\x_1$ and $\x_2$ along with a smooth path $\bgamma_0(s) : [0,1]\rightarrow \Rn$ such that $\bgamma_0(0) = \x_1$ and $\bgamma_0(1) = \x_2$. Consider further a family of solutions to the natural gradient dynamics $\bgamma(s,t)$ such that $\bgamma(s,0) = \bgamma_0(s)$ and
\[
\frac{\rm d}{{\rm d} t} \bgamma(s,t) = \h( \bgamma(s,t) )
\]
Since the Riemannian Hessian of $f$ is positive semi-definite, it follows that these dynamics are semi-contracting and thus have bounded solutions.

It then follows that
\begin{equation}
\H( \bgamma(s,t) ) \, \pd{\bgamma}{s} \rightarrow 0
\label{eq:H_gamma_s_convergence}
\end{equation}
as $t \rightarrow \infty$. Denote the conventional gradient along the solution family as
$
\g(s,t) = \pd{f}{\x}( \bgamma(s,t) )\,.
$
For any fixed time, the gradient $\g(s,t)$ can be constructed through a path integral, along the solution family $\bgamma(s,t)$, according to
\[
\g(s,t) = \int_0^s \boldsymbol{\Phi}_{\A_t}(s,\sigma)\  \H(\bgamma(\sigma,t)) \ \pd{\bgamma}{s}(\sigma,t) \ {\rm d}\sigma
\]
where $\A_t(s)$ depends on the Christoffel symbols of the metric, and $\boldsymbol{\Phi}_{\A_t}(\cdot,\cdot)$ the state transition matrix associated with $\A_t(\cdot)$. Further detail is given in Appendix \ref{app:Semi} $-$ essentially, the state transition matrix accounts for parallel transport of changes to the gradient along the solution family.

Since $\bgamma(s,t)$ bounded, it follows that $\A_t(s)$ is bounded, and likewise so is $\boldsymbol{\Phi}_{\A_t}(\cdot, \cdot)$. Combining this result with \eqref{eq:H_gamma_s_convergence} it follows that:
\[
\g(s,t) \rightarrow \bzero {\rm~as~} t\rightarrow \infty,~ \forall s \in [0,1]
\]
Since the gradient converges to zero on the path $\bgamma$, and the path distance is bounded, it follows that $f(\x_1) = f(\x_2)$. Thus, any two local optima must necessarily be global optima.
\end{proof}

 In the context of machine learning, note that so-called residual networks \cite{He16} largely reduce combinatorial redundancies and thus the frequency of local minima.

\section{Geodesic Primal-Dual Optimization}\label{GPD_opt}

The Primal-Dual algorithm is widely used
in optimization to determine saddle points.
Saddle points generally occur when a system is simultaneously minimizing a function over one set of its
variables and maximizing it over another variable set. Saddle points appear in economics in the context of
market equilibrium~\cite{arrow1958studies}, and also appear naturally in constrained optimization~\cite{BoydVanberghe04}, where
Lagrange parameters play the role of dual variables.
When a function is strictly convex in a subset of its variables, and strictly concave in the remaining, gradient descent/ascent dynamics converge to a unique saddle
equilibrium~\cite{arrow1958studies,paganini}. Within the context of constrained optimization, these dynamics are known as the primal-dual dynamics.
Such dynamics play an important role e.g., in machine learning, for instance in the adversarial training of deep networks~\cite{madry2017towards}, in
reinforcement learning~\cite{Wang17} and actor-critic methods, and in support vector machine representations~\cite{kosaraju2018primal}. More generally, they are central to a large class of practical min-max problems, such as, e.g.,
nonlinear electrical networks modeled in terms of Brayton-Moser mixed potentials~\cite{Kostya,ortega2003power, cavanagh2018transient}.

Continuous-time convex primal-dual optimization is analyzed from a nonlinear contraction perspective in~\cite{Kostya}, building on a earlier result of~\cite{modular}. As we now show, Theorem~\ref{th_main} yields a natural extension
to geodesic primal-dual optimization, where convexity in terms of primal and dual variables
is replaced by g-convexity, thus broadening the above results to state-dependent metrics.

Consider a scalar cost function $\script{L}({\bf x},{\bf \llambda}, t)$,
possibly time-dependent, with $\script{L}$ g-strongly convex over $\x$ and g-strongly concave over $\llambda$ in metrics
$\M_\x(\x)$ and $\M_{\llambda}(\llambda)$ respectively. Consider the geodesic primal-dual dynamics, which
we define as
\begin{subequations} \label{eq:PD}
\begin{align} \label{eq:primal}
{\bf M_\x}({\bf x})\ \dot{\bf x} \ &= -\partial_{\bf x} \script{L} \\
{\bf M_\llambda}({\bf \llambda})\ \dot{\llambda} \ &= \ \partial_{\llambda}  \script{L}
\end{align}
\end{subequations}
Using the metrics ${\bf M_x}({\bf x})$ and ${\bf M_\llambda}({\llambda})$ extends the standard case \cite{Kostya}, where they would be replaced by constant, symmetric positive definite matrices.

\begin{theorem}\label{th_PD}
Consider a scalar cost function $\script{L}({\bf x},{\bf \llambda}, t)$,
with $\script{L}$ g-strongly convex over $\x$ and g-strongly concave over $\llambda$ in metrics
$\M_\x(\x)$ and $\M_{\llambda}(\llambda)$ respectively. Then, the geodesic primal-dual dynamics~\eqref{eq:PD} is globally contracting, in metric
\begin{equation}\label{GPD}
   {\bf M}({\bf x},{\bf \llambda})  = \begin{bmatrix}
     {\bf M_x}({\bf x}) & {\bf 0}\\
    {\bf 0} & {\bf M_\llambda}({\bf \llambda})
    \end{bmatrix}
\end{equation}
\end{theorem}

\begin{proof} \ Letting $\z = [{\x}\T, {\llambda}{}\T]\T$ and $\dot{\z}  = \f(\z,t)$ denote the overall system dynamics, the system Jacobian can be written
\begin{equation} \label{eq:PD_Jacobian}
{\bf A}(\x,\llambda, t) =  \pd{\f}{\z} =
\begin{bmatrix}
- \partial_{\bf x} ({\bf M_x}^{-1} \partial_{\bf x} \script{L}) & - \ {\bf M_\x}^{-1}\ \pdd{\script L}{\llambda}{\x}\\
   {\bf M_\llambda}^{-1} \ \pdd{\script L}{\x}{\llambda} & \partial_{\llambda} ({\bf M_\llambda}^{-1} \partial_{\llambda} \script{L})
    \end{bmatrix}
\end{equation}
\noindent so that, using Theorem~\ref{th_main},
\begin{equation}\nonumber
     {\bf M}\T {\bf A}  + \ {\bf A}\T {\bf M} \ + \ \dot{\bf M } \ = \ -2 \ \begin{bmatrix}
     {\bf H_x} & {\bf 0}\\
    {\bf 0} & - {\bf H_\llambda}
    \end{bmatrix} \ < \ {\bf 0}
    \end{equation} \end{proof}

\section{Applying Contraction Analysis Tools to Geodesic Optimization}
\label{sec:application}

Theorem \ref{th_main} immediately implies that existing contraction
tools can be directly applied in the context of geodesic
optimization. These include, in particular, combination
properties, virtual systems,

and stochastic versions.

\subsection{Combination Properties}

We now discuss how basic combination properties of
contracting systems~\cite{Slotine98,modular} can be exploited to build on the results above.
While these properties derive from simple matrix algebra and in principle could be proven directly from the definition of geodesic convexity, as we will see most
rely for their practical relevance on the flexibility

afforded by the contraction analysis point of view.

\subsubsection{Sum of g-convex}

If two functions $f_1(\x,t)$ and $f_2(\x,t)$ are g-convex in the same metric for each $t$, then their sum $f_1(\x,t) + f_2(\x,t)$ is g-convex in the same metric.

\begin{example}
\label{example:Parallel_Geo_Desc}
Consider a function $f_1(\x_1,\y_1,t)$ g-convex for each $t$ in a block diagonal metric ${\rm BlkDiag}(\M_{x_1}(\x_1),\M_y(\y_1))$ and a function $f_2(\x_2,\y_2,t)$ g-convex for each $t$ in a block diagonal metric ${\rm BlkDiag}(\M_{x_2}(\x_2), \M_y(\y_2))$.
Then, the function:
\[
f(\x_1, \x_2, \y,t) =  f_1(\x_1,\y,t)+f_2(\x_2,\y,t)
\]
is g-convex in metric ${\rm BlkDiag}(\M_{x_1},\M_{x_2},\M_y)$ for each $t$.

\end{example}

\subsubsection{Skew-Symmetric Feedback Coupling}

The system Jacobian~\eqref{eq:PD_Jacobian} for the primal-dual dynamics is a special case of a contraction property that holds more broadly.

Assume that a scalar function $f_1(\x_1,\x_2)$ is $\alpha_1$-strongly g-convex in $\x_1$ in a metric $\M_1(\x_1)$ for each fixed $\x_2$, and similarly that a scalar function $f_2(\x_1,\x_2)$ is $\alpha_2$-strongly g-convex in a metric $\M_2(\x_2)$ for each fixed $\x_1$. If $f_1$ and $f_2$ satisfy the scaled skew-symmetry property
\[
\pdd{f_1}{\x_1}{\x_2} \ = \ - \ k \ \left(\pdd{f_2}{\x_2}{\x_1} \right)\T
\]
where $k$ is some strictly positive constant, then the natural gradient dynamics
\begin{align}
\label{eq:Nash_dyn}
\dot{\x}_1 &= - \M_1(\x_1)^{-1} \ \pd{f_1(\x_1,\x_2)}{\x_1}  \\
\dot{\x}_2 &= - \M_2(\x_2)^{-1} \ \pd{f_2(\x_1,\x_2)}{\x_2} \nonumber
\end{align}
is contracting with rate $ \ \min(\alpha_1, \alpha_2) \ $ in metric $\M(\x_1,\x_2) = {\rm BlkDiag}(\M_1(\x_1), k \M_2(\x_2))$. Since the {\em overall} system is both contracting and autonomous, it
tends to a unique equilibrium~\cite{Slotine98} $(\x_1^*, \x_2^*)$ which satisfies the Nash-like conditions
\begin{align*}

\x_1^* &= {\rm argmin}_{\x_1} f_1(\x_1,\x_2^*) \\
\x_2^* &= {\rm argmin}_{\x_2} f_2(\x_1^*,\x_2)
\end{align*}

Using again the machine learning context as an example, such two-player game dynamics can occur in certain types of adversarial training.

The result extends to a game with an arbitrary number of players. Consider $n$ functions $\{f_i(\x_1,\ldots,\x_n)\}_{i=1}^n$ such that each $f_i$ is $\alpha_i$-strongly g-convex over $\x_i$ in a metric $\M_i(\x_i)$. If the functions satisfy the skew-symmetry conditions
\[
\pdd{f_i}{\x_i}{\x_j} = - k_j \left( \pdd{f_j}{\x_j}{\x_i} \right)\T
\]
for each $j>i$, then the suitable generalizations of \eqref{eq:Nash_dyn} result in a coupled system that is contracting with rate ${\rm min}(\alpha_1, \ldots, \alpha_n)$ in the metric $\M= {\rm BlkDiag}(\M_1, k_2 \M_2, \ldots, k_n \M_n)$. The overall system converges to a unique Nash-like equilibrium satisfying
\[
\x_1^* = {\rm argmin}_{\x_1} f_1(\x_1,\x_2^*,\ldots,\x_n^*) \\
\]
and a similar relation for each other player.

\subsubsection{Hierarchical Natural Gradient}
Consider a function $f_1(\x_1)$ $\alpha_1$-strongly g-convex in a metric $\M_1(\x_1)$, and a function $f_2(\x_1,\x_2)$ $\alpha_2$-strongly g-convex in a metric $\M_2(\x_2)$ for each given $\x_1$, where the $\x_i$ may have different dimensions. Then, the hierarchical natural gradient dynamics
\begin{align*}
\dot{\x}_1 &= - \M_1(\x_1)^{-1} \ \pd{f_1(\x_1)}{\x_1}  \\
\dot{\x}_2 &= - \M_2(\x_2)^{-1} \ \pd{f_2(\x_1,\x_2)}{\x_2}
\end{align*}
is contracting with rate $ \ \min(\alpha_1, \alpha_2) \ $ in metric $\M(\x_1,\x_2) = {\rm BlkDiag}(\M_1(\x_1), \M_2(\x_2))$, under
the mild assumption that the coupling Jacobian
is bounded~\cite{Slotine98}. Since the {\em overall} system is both contracting and autonomous, it
tends to a unique equilibrium~\cite{Slotine98} at rate $\min(\alpha_1, \alpha_2)$, and thus to the unique solution of
\begin{align*}
\pd{f_1(\x_1)}{\x_1} &= {\bf 0}\\
\pd{f_2(\x_1,\x_2)}{\x_2}  &= {\bf 0}
\end{align*}

By recursion, this structure can be chained an arbitrary number of times, or applied to any cascade or directed acyclic graph of natural gradient dynamics. Such hierarchical optimization may play a role, for instance, in backpropagation of natural gradients in machine learning, with all descents occurring concurrently rather than in sequence.

\begin{remark}

In large-scale optimization settings such as those appearing commonly in machine learning, natural gradient with a fully-dense metric can become intractable. In specific cases, such as natural gradient descent based on Fisher information~\cite{Amari98b},
computationally effective approximations have been derived~\cite{Martens15}.
In addition, the combination of simple (e.g., diagonal) metrics through hierarchical structures lends an opportunity to recover significant complexity at broad scale $-$ see, e.g., the hierarchical combination of scalar metrics to learn hierarchical representations of symbolic data in  \cite{Nickel17,Nickel18}. Such simpler metrics are also
well motivated in the context of positive or monotone
systems \cite{Rantzer15,Manchester17}.\end{remark}

\subsection{Virtual Systems}

The use of virtual contracting systems~\cite{wei,jouffroy,ManchesterCCM} allows guaranteed exponential convergence to a unique minimum to be extended to classes of dynamics which are not pure natural gradient.

For instance, it is common in optimization to adjust
the learning rate as the descent progresses.
Consider a natural gradient descent with the function $f(\x)$ $\alpha$-strongly g-convex in metric $\M(\x)$, and define the new system
\begin{equation}
\dot{\x} = - p(\x,t)\, \M(\x)^{-1}\, \pd{f}{\x}(\x)
\label{p(x,t)}
\end{equation}
where the smooth scalar function $p(\x,t)$ modulates the learning rate~\cite{Amari98b} and is uniformly positive definite,
\[ \exists \ p_{min} > 0 , \ \forall t \ge 0, \ \forall \x, \ \ \ \
p(\x,t) \ge p_{min}
\]
Let us show that this system tends exponentially to the minimum $\x^*$ of $f(\x)$.

Consider the auxiliary, \emph{virtual} system,
\begin{equation}
\dot{\y} = - \M(\y)^{-1}\, \pd{(\ p(\x,t)f(\y)\ )}{\y}
\label{virtual_p(x,t)}
\end{equation}
For this system, $p(\x(t) ,t)$ is an external, uniformly positive definite function of time, and thus
\[\pd{(\ p(\x,t)f(\y)\ )}{\y} \ = \ p(\x,t) \ \pd{f}{\y}(\y)
\]
so that the contraction of~\eqref{eq:NatGrad} with rate $\alpha$ implies the contraction of~\eqref{virtual_p(x,t)} with rate $ \ \alpha p_{min} \ $. Since both $\x(t)$ and $\x^*$ are particular solutions
of~\eqref{virtual_p(x,t)}, this implies in turn that
$\x(t)$ tends to $\x^*$ with rate $ \ \alpha p_{min} \ $.

Note that since we only assumed that $p(\x,t)$ is uniformly positive definite, in general the actual system~\eqref{p(x,t)} is not contracting with respect to the metric $\M(\x)$. The approach extends immediately to the primal-dual
context of Section~\ref{GPD_opt}.

\section{Conclusions}
\label{sec:Conclusions}

This paper has outlined a fundamental connection between g-strongly convex functions and contraction of their natural gradient dynamics. This observation sets the foundation for contraction results to be exploited for the construction and analysis of g-convex optimization for complex systems. A natural next step for the application of contraction in g-convex optimization is to design geodesic quorum sensing ~\cite{tabareau,russo} algorithms for synchronization \cite{bouvrie}, as well as other consensus mechanisms considering time-delays~\cite{Wei_delay,cloud}, which may serve as the basis for distributed and large-scale optimization techniques on Riemannian manifolds. Other future applications will consider stochastic gradient descent in the Riemannian setting~\cite{Bonnabel13} with quorum sensing extensions (as e.g. in~\cite{Lecun}) to address inherent noise~\cite{tabareau}. Such advances could have clear impact, e.g., in the context of machine learning, among other applications.

\appendices
\section{Proof of Theorem 1}

\begin{proof}
Recall that $\alpha$-strong geodesic convexity of $f(\x,t)$ in the metric $ \M(\x)$ (for each $t$) is equivalent to the Riemannian Hessian of $f$, denoted $\H(\x,t)$, satisfying:
\[
\H(\x,t) \succeq \alpha \M(\x)
\]
We now show that this property is exactly the same as contraction of the natural gradient dynamics \eqref{eq:NatGrad} in the metric $\M(\x)$. Specifically,
given $\h(\x,t)$ from~\eqref{eq:NatGrad}, and
defining
\[
\mathbf{Q} = \M \left(\pd{\h}{\x} \right) + \left( \pd{\h}{\x} \right)\T \M + \dot{\M}
\]
we show that $\mathbf{Q} = -2 \mathbf{H}$, thus proving the theorem.

In coordinates, entries of $\H$ are given by
\[
H_{ij} = \pdd{f}{\xsubi}{\xsubj} - \Gamma_{ij}^k \left( \pd{f}{\xsub{k}} \right)
\]
where $\Gamma_{ij}^k$ denotes the Christoffel symbol of the second kind
\[
\Gamma_{ij}^m = \frac{1}{2} M^{mk} \left( \pd{M_{ik}}{\xsubj} + \pd{M_{jk}}{\xsubi} - \pd{M_{ij}}{\xsubk} \right)
\]
and the usual Einstein summation convention is applied  (implying e.g. a sum over $k$ in the above).

Consider the partials of the natural gradient system,
\begin{align*}
\pd{h_k}{\xsub{j}} &= \pd{}{\xsub{j}}\left[ - M^{k\ell}\, (\pd{f}{\xsub{\ell}}) \right] \nonumber \\
			  &= -M^{k\ell} (\pdd{f}{\xsub{\ell}}{\xsub{j}})  + M^{kr} (\pd{M_{rs}}{\xsub{j} }) M^{s\ell} (\pd{f}{\xsub{\ell}})
\end{align*}
Using this result
\begin{align*}
Q_{ij} =& M_{ik} (\pd{h_k}{\xsub{j}}) + M_{jk} (\pd{h_k}{\xsub{i}}) - (\pd{M_{ij}}{\xsub{k}}) M^{k\ell} (\pd{f}{\xsub{\ell}}) \\
			   =& -M_{ik} M^{k\ell} (\pdd{f}{\xsub{j}}{\xsub{\ell}})  + M_{ik} M^{kr} (\pd{M_{rs}}{\xsub{j}}) M^{s\ell} (\pd{f}{\xsub{\ell}}) \\
			    &- M_{jk} M^{k\ell} (\pdd{f}{\xsub{i}}{\xsub{\ell}})  + M_{jk} M^{kr} (\pd{M_{rs}}{\xsub{i}}) M^{s\ell} (\pd{f}{\xsub{\ell}}) \\
			    &- (\pd{M_{ij}}{\xsub{k}}) M^{k\ell} (\pd{f}{\xsub{\ell}})
\end{align*}
Noting that $M_{ik} M^{kj} = \delta_{ij}$, with $\delta_{ij}$ the Kronecker delta,
\begin{align*}
Q_{ij} &= -2 \pdd{f}{\xsub{i}}{\xsub{j}} + M^{s\ell}\left(\pd{M_{is}}{\xsub{j}} + \pd{M_{js}}{\xsub{i}}  - \pd{M_{ij}}{\xsub{s}}  \right) \pd{f}{\xsub{\ell}} \nonumber \\
			   &= -2 \left(\pdd{f}{\xsub{i}}{\xsub{j}} - \Gamma^k_{ij} \, \pd{f}{\xsub{k}} \right) \\&= - 2 H_{ij} \nonumber
\end{align*}
\end{proof}

\section{Supplemental Derivation for Theorem~\ref{thm:semi}}
\label{app:Semi}

The appendix supplements the derivation of the gradient $\g(s,t)$ via a path integral along the solution family $\bgamma(\sigma,t)$ introduced in the proof of Theorem~\ref{thm:semi}.

Consider $\g(s,t) = \pd{f}{\x}(\bgamma(s,t))$ for a fixed time $t$. The gradient can be computed at any point through integration of the conventional Hessian along the path $\bgamma$ as
\begin{equation*}
g_i(s,t) =\int_0^s \pd{f}{ij}(\bgamma(\sigma,t)) \ \left[ \pd{\bgamma}{s}(\sigma,t)\right]_j {\rm d} \sigma
\end{equation*}
Applying \eqref{eq:riemhessian}, this integral can be written as
\[
g_i(s,t) =\int_0^s \left[ \Gamma_{ij}^k \, g_k \, \pd{\gamma_j}{s} + H_{ij}\, \pd{\gamma_j}{s} \right]_{\x = \bgamma(\sigma,t)} \,  {\rm d} \sigma
\]
Consider the matrix $\A_t(s)$ defined as:
\[
\left[ \A_t(s) \right]_{ik} = \sum_j \left[ \Gamma_{ij}^k \, \pd{\gamma_j}{s} \right]_{\x = \bgamma(s,t)}
\]
It then follows that:
\[
\frac{\partial}{ \partial s} \g(s,t) =  \A_t(s) \g(s,t) +  \H(\bgamma) \, \pd{\bgamma}{s}
\]

 Applying linear systems theory, for each fixed $t$, the solution to this differential equation is given by
\[
\g(s,t) = \int_0^s \boldsymbol{\Phi}_{\A_t}(s,\sigma)\  \H(\bgamma(\sigma,t)) \ \pd{\bgamma}{s}(\sigma,t) \ {\rm d}\sigma
\]
where $\boldsymbol{\Phi}_{\A_t}(\cdot, \cdot)$ denotes the state-transition matrix of $\A_t(\cdot)$ satisfying $\boldsymbol{\Phi}_{\A_t}(s_0,s_0) = \mathbf{I}$ and
\[
\frac{\rm d}{{\rm d} s} \boldsymbol{\Phi}_{\A_t}(s,s_0) = \A_t(s)\,\boldsymbol{\Phi}_{\A_t}(s,s_0)
\]

\bibliography{Geodesic.bbl}

\end{document}


TITLE: Semantic proofs to syntactic proofs
QUESTION [9 upvotes]: Given a first-order logic theory $T$ and and a formula $F$, suppose I have semantically proved that $T\vdash F$. That is, I have proved that any model $M$ of $T$ satisfies $F$ and I conclude by Gödel's completeness theorem.
Do I have a general algorithm to extract from the above a syntactic proof of $T\vdash F$, i.e. a finite sequence of formulas that respects inference rules, uses $T$ and finishes at $F$ ?
If no such algorithm exists, then did I really prove $T\vdash F$ ? The completeness theorem was just an example of how to indirectly prove that there exists a proof of $T\vdash F$, without explicitly giving this latter proof. What if my indirect proof uses an inaccessible cardinal, do I have to mention the awkward
$$ (\text{ZFC + Inaccessible cardinal})\; \vdash\; (T \vdash F) $$
And then this proof can also be indirect, so I might continue stacking theories to the left and it becomes a nightmare. Don't semantic or other indirect formal proofs somewhat defeat the purpose of formal logic, that we should be absolutely certain that the formal proofs exist and are correct?

REPLY [6 votes]: Keep in mind that proofs are "easily-recognizable finite strings of symbols" - precisely, we can effectively enumerate all proofs from a given theory. So we can always find a formal proof of $F$ from $T$ - if one exists - effectively by simply checking each $T$-proof in order until we find one which is a proof of $F$. This is unsatisfying, but is perfectly precise and effective.


TITLE: Gradient descent inside the expectation-maximization (EM) algorithm
QUESTION [2 upvotes]: I am feeling super uncertain about how much I can play around with the EM algorithm. Here is my question:
In the EM algorithm, during the M-step, one attempts to find a parameter value, $\theta$, that maximizes $Q$.
Sometimes the function of $Q$ with respect to that parameter value is continuous, differentiable and concave, but does not have a closed form solution for the stationary points and requires some form of numerical techniques, like gradient descent or Newton-Raphson method, to find the stationary points.
Am I correct to conclude that one does not break the convergences of the EM algorithm if you use these numerical methods to optimize $Q$ with respect to some parameter value $\theta$.
In fact, am I correct to conclude that one could use any optimization technique in the E step for all the parameters we wish to optimize for?

REPLY [3 votes]: Yes, you can use any optimization technique, including numerical, in the M (maximization) step.
In fact, you needn't even maximize; as long as the M step improves the objective function Q. The Generalized EM Algorithm, is described for example, in section 7 of "The EM Algorithm As a Lower Bound Optimization Technique", by Rave Harpaz and Robert Haralick. I have substituted your notation for the paper's notation in the below quote from that section.

In the M-step, the parameters $\theta$ were chosen as the value for which the
lower bound of the objective Q was maximized. While this ensures the
greatest increase in Q and subsequently the
log-likelihood, it is possible to relax the requirement of
maximization to one of simply increasing Q. This approach
to simply increase and not necessarily maximize is known as the
Generalized Expectation Maximization (GEM) algorithm.