m#1 m⌊#1 ⌋ ✓ × mo #2  A: #1 #2  A: #1 m#1 mo #2  RA: #1 #2  RA: #1 mo #2  L:
#1 #2  L: #1 mo #2  GL: #1 #2  GL: #1 mo #2  JA: #1 #2  L: #1 m  #1  m#1

# Low-dimensional approximations of the conditional law of Volterra processes:
a non-positive curvature approach

Reza Arabpour111Vector Institute and Department of Mathematics, McMaster
University., Luca Galimberti222Department of Mathematics, King’s College
London., Anastasis Kratsios333Vector Institute and Department of Mathematics,
McMaster University., Giulia Livieri444Department of Statistics, London School
of Economics and Political Science., John Armstrong555Department of
Mathematics, King’s College London. Reza Arabpour11footnotemark: 1
<EMAIL_ADDRESS>
Vector Institute and Department of Mathematics, McMaster University John
Armstrong11footnotemark: 1<EMAIL_ADDRESS>
Department of Mathematics, King’s College London. Luca
Galimberti11footnotemark: 1 22footnotemark: 2<EMAIL_ADDRESS>
Department of Mathematics, King’s College London Anastasis
Kratsios11footnotemark: 1 22footnotemark: 2<EMAIL_ADDRESS>
Vector Institute and Department of Mathematics, McMaster University Giulia
Livieri11footnotemark: 1 22footnotemark: 2<EMAIL_ADDRESS>
Department of Statistics, London School of Economics and Political Science
Alphabetic ordering.Primary authors, equal contribution, alphabetical
ordering.Corresponding author.

###### Abstract

Predicting the conditional evolution of Volterra processes with stochastic
volatility is a crucial challenge in mathematical finance. While deep neural
network models offer promise in approximating the conditional law of such
processes, their effectiveness is hindered by the curse of dimensionality
caused by the infinite dimensionality and non-smooth nature of these problems.
To address this, we propose a two-step solution. Firstly, we develop a stable
dimension reduction technique, projecting the law of a reasonably broad class
of Volterra process onto a low-dimensional statistical manifold of non-
positive sectional curvature. Next, we introduce a sequentially deep learning
model tailored to the manifold’s geometry, which we show can approximate the
projected conditional law of the Volterra process. Our model leverages an
auxiliary hypernetwork to dynamically update its internal parameters, allowing
it to encode non-stationary dynamics of the Volterra process, and it can be
interpreted as a gating mechanism in a mixture of expert models where each
expert is specialized at a specific point in time. Our hypernetwork further
allows us to achieve approximation rates that would seemingly only be possible
with very large networks.

Keywords: Geometric Deep Learning, Measure-Valued Stochastic Processes, Non-
Positive Curvature, Barycenters, Universal Approximation, hypernetworks,
Mixture of Experts.

##  1 Introduction

The dynamic prediction of the conditional distribution
$\mathbb{P}[X_{t+1}|X_{[0:t]}=x_{[0:t]}]$ of a non-Markovian Volterra process
$X_{\cdot}$, conditioned on its realized path $x_{[0:t]}$ up to any given time
$t$, is a fundamental problem spanning the sciences, ranging from Bayesian
modelling (e.g., Bernardo and Smith (1994)) to mathematical finance (e.g.,
Back et al. (2004)). While there are several machine learning models for
learning a process conditional distribution on its historical paths; e.g.
Bishop (1994); Gauthier (2014); Mirza and Osindero (2014); Borovykh et al.
(2017); Chevyrev and Oberhauser (2022) and deep learning models which can
approximate signed measure-valued functions Chen and Chen (1995); Korolev
(2022); Benth et al. (2023); Cuchiero et al. (2023), the available
quantitative approximation bounds for measure-valued models; e.g. Kratsios
(2023); Acciaio et al. (2023) suggest that the measure-valued maps cannot be
efficiently approximated. This is due to two factors; firstly, they are
infinite-dimensional, meaning they suffer from extreme forms of the curse of
dimensionality; see Lanthaler and Stuart (2023) for a lower-bound in the
linear case. Secondly, most spaces of probability measures, e.g. Wasserstein
spaces, do not have any smooth or linear structure, which a deep learning
model can naturally leverage. To the best of the authors’ knowledge, there are
currently no available deep-learning models which can approximate the evolving
conditional distribution of most stochastic processes while also depending on
a computationally feasible number of parameters.

To address this issues, we present a two-step approach for dynamically
approximating the conditional evolution of a class of discrete-time stochastic
processes $X_{\cdot}$. By allowing for an irreducible dimension-reduction-type
error, our first step corrects the deficits of the geometries of most spaces
of probability measures by projecting
$\mathbb{P}[X_{t+1}|X_{[0:t]}=x_{[0:t]}]$ onto the (smooth) Riemannian
manifold $\mathcal{N}_{d}$ of non-singular $d$-dimensional Gaussian measures
equipped with a new Riemannian metric which has a non-positively curved
geometry when the $d\geq 2$. This can be contrasted with the standard
Fisher(-Rao) metric (Dowty, 2018), from information geometry, which is only
non-positively curved when $d=1$ (Skovgaard, 1984). Critically, we show that
(refer to Theorem 9), though our new projection shares an analogous
characterization as “distance” minimizer to information-theoretic projections,
e.g. the $I$-projection Csiszár and Matus (2003); Csiszár (1975), it preserves
the dependence of $\mathbb{P}[X_{t+1}|X_{[0:t]}=x_{[0:t]}]$ on the realized
path $x_{[0:t]}$ in a stable way as it is Lipschitz. Furthermore, the
dependence of the projected dynamical system (see below) on the historical
paths generated by $X_{\cdot}$ is found to be proportional to the persistence
of the memory of the stochastic process itself (Theorem 11).

The projection of the conditional distribution
$\mathbb{P}[X_{t+1}|X_{[0:t]}=x_{[0:t]}]$ onto $\mathcal{N}_{d}$, defines a
dynamical system on the Riemannian manifold $\mathcal{N}_{d}$. This long
memory dynamical system is then approximated by a dynamic version of the
geometric deep learning model of Kratsios (2023). We prove a quantitative
universal approximation theorem which shows that this system can be
approximated by our sequential deep learning model (Theorem 13). Our result
has the following properties: 1) it shows that any infinite-memory dynamical
systems between non-positively curved Riemannian manifolds, e.g. generalized
hyperbolic spaces, can be approximated provided that its memory fades. This is
similar to results in the reservoir computing literature Grigoryeva and Ortega
(2018, 2019). 2) It shows that the geometric deep learning model does not
suffer from the curse of dimensionality when the dynamical systems evolve
according to smooth dynamics; in particular, this happens for the projected
conditional law of $X_{\cdot}$ when its generator is “smooth” and has no
hidden random state. 3) In the static setting, i.e. when ignoring the flow of
time, our bounds are polynomially tighter than the non-Euclidean universal
approximation theorems of Kratsios (2023); Acciaio et al. (2023) and in the
dynamic setting, our bounds are exponentially tighter than the static setting.
4) Our constructive proof yields a federated algorithm where a sequence of
independent experts approximate the dynamical system at individual points in
time, after which an overarching hypernetwork Ha et al. (2017) is used to
synchronize them and create the recurrence without having to optimize
(backpropagate) in time. While our approximation theorems hold for dynamical
systems on any non-positively curved Riemannian manifold, we focus our
exposition on the statistical manifold $\mathcal{N}_{d}$. We empirically
verify the trainability of the proposed model and the role of each component
via an ablation study evaluating the dependence of the model on its
parameters.

All code can be found at https://github.com/arabporr/HyperNetwork .

##### Outline

The paper is organized as follows.

  1. 1.

Section 2 reviews notions from metric geometry required to formulate our main
results.

  2. 2.

Section 3 formalizes the setting and introduces our new perturbation of the
Fisher metric on $\mathcal{N}_{d}$, which has a well-behaved geometry.

  3. 3.

Section 4 contains our results on the Gaussian random projection of the
conditional distribution of Volterra processes.

  4. 4.

Section 5 contains universal approximation theorems for our model, first in
the static and then in the dynamic settings.

  5. 5.

Section 6 contains an ablation study of our model, which confirms our
theoretical results.

All proofs are deferred to the appendix, which also describes results from
differential geometry for which we have been unable to provide a proper
reference in the literature. The appendix also contains an example outlining
the problem with information projection for dimension reduction of the
conditional law of Volterra processes.

We assume the reader to be familiar with the terminology of stochastic
processes and optimal transport. Both notions are widely used in machine
learning, for example in denoising diffusion models Sohl-Dickstein et al.
(2015) and generative adversarial networks Arjovsky et al. (2017).

###  1.1 Related Work

Our results focus on the broad class of stochastic processes known as Volterra
processes, which represent a rich yet well-structured class of non-Markovian
stochastic differential equations (SDEs, henceforth) with latent stochastic
factors. Both the discrete and the continuous versions of stochastic Volterra
processes, and their generalizations Jaber et al. (2024), play a crucial role
in mathematical finance (e.g., (Jacquier et al., 2018; Abi Jaber et al., 2019;
Cuchiero and Teichmann, 2020; Bondi et al., 2024), reservoir computing (e.g.,
Grigoryeva and Ortega (2019); Gonon et al. (2022)), engineering (e.g., Shiki
et al. (2017)), and computational biology (e.g., Korenberg and Hunter (1996)).

Our projection of the conditional law of stochastic Volterra processes
exploits the existence and stability of barycenters of probability measures on
Riemannian manifolds of non-positive curvature, in the sense of Aleksandrov
(1951) and Rešetnyak (1960). Various notions of barycenters, e.g. Wasserstein
barycenters (Cuturi and Doucet, 2014; Srivastava et al., 2018; Puccetti et
al., 2020; Altschuler and Boix-Adserà, 2021; Kolesov et al., 2024a, b), have
been used in generative modelling as they provide a natural tool for averaging
distributions. The new element of our approach is to consider barycenters of
distributions, also called the Karcher/Fréchet mean, whose definition is
rooted in Riemannian (Grove and Karcher, 1973; Karcher, 1977b; Kim and Pass,
2020) and metric geometry (Sturm, 2002, 2003; Ohta, 2012; Mendel and Naor,
2013; Yokota, 2017). These provide a non-linear generalization of the
expectation of a random variable taking values in a non-linear space and are
commonly used in computational geometry (Bini and Iannazzo, 2013), shape
analysis (Le, 2001; Fletcher et al., 2004), geometric statistics (Lim and
Pálfia, 2020; Bhattacharya and Patrangenaru, 2003, 2005), and recently in
differential privacy (Utpala et al., 2023). We will use them to determine the
single most representative $d$-dimensional non-singular Gaussian measure,
amongst a family of random Gaussian measures in $\mathcal{N}_{d}$ in which we
encode the conditional law of the Volterra process.

This notion of barycenter, or intrinsic expectation, is only well-defined
globally for manifolds with non-positive curvature (Afsari, 2011). Examples
include the hyperbolic spaces used in hierarchical clustering (Chami et al.,
2020) or latent graph inference for GNNs (de Ocáriz Borde et al., 2022). Such
geometries appear for various statistical manifolds (Pinele et al., 2020; Le
Brigant et al., 2021) with the Fisher-Rao metric (Ay et al., 2017); notably
the one-dimensional non-degenerate Gaussian distributions (Atkinson and
Mitchell, 1981; Costa et al., 2015; Pinele et al., 2020). However, as shown in
(Skovgaard, 1984), the space $\mathcal{N}_{d}$ of multivariate non-degenerate
Gaussian distributions fails to have non-positive curvature in dimension $2$
or greater. Our perturbation of the standard information-theoretic geometry on
the statistical manifold $\mathcal{N}_{d}$ (see Proposition 7) identifies a
non-positively curved complete Riemannian geometry on $\mathcal{N}_{d}$ which
coincides with the standard information geometry on any submanifolds of
Gaussian distributions with fixed mean. This allows us to apply the intuition
from information geometry, but with the analytic benefit of a non-positive
metric geometry. Moreover our construction provides a closed-form expression
for the Riemannian distance: for typical statistical families, the Riemannian
distance in the Fisher-Rao metric can only be computed numerically.

Projecting the conditional distribution of the stochastic Volterra process
$X_{\cdot}$, conditioned on its realized path $x_{[0:t]}$ up to any time $t$,
down to $\mathcal{N}_{d}$ results in a (generalized) dynamical system between
finite-dimensional spaces. Since all the resulting spaces are finite-
dimensional and well-structured, one can reasonably hope that this system can
be approximated without the curse of dimensionality if the involved maps are
smooth enough, a feature not shared by infinite-dimensional approximation
problems (Lanthaler and Stuart, 2023; Acciaio et al., 2023). There is a well-
developed literature on the approximation of dynamical systems by recurrent
deep learning models such as reservoir computers (Grigoryeva and Ortega, 2018;
Gonon and Ortega, 2019; Grigoryeva and Ortega, 2019; Gonon and Ortega, 2021),
recurrent neural networks (Li et al., 2022a; Hutter et al., 2022), or
transformers Yun et al. (2020). However, the available universal approximation
results in the literature only apply to dynamical systems between linear input
spaces and systems whose dynamics do not change in time. We, therefore first
extend the static geometric deep learning of Kratsios (2023) to a
sequential/dynamic model capable of processing sequences of inputs and outputs
in any given appropriate pair of non-positively curved Riemannian manifolds,
e.g. on $\mathcal{N}_{d}$, and we then prove a universal approximation showing
that it can approximate most time-inhomogeneous dynamical systems between
these spaces, possibly having infinite but polynomially fading memory (Theorem
15).

Our estimates on the parameters and guarantees are novel even in the Euclidean
setting: many of the aforementioned approximation theorems for recurrent
neural networks (RNNs) are only qualitative and no rates are provided,
especially for low-regularity dynamical systems with slowly-fading infinite
memory. Even in the static setting, our results significantly improve the
rates for the available universal approximation theorems for deep learning
models between Riemannian manifolds Kratsios (2023); Acciaio et al. (2023).
Furthermore, in the static setting (Theorem 17), our approximate rates are
optimal since they match those of Yarotsky (2018); Yarotsky and Zhevnerchuk
(2020); Lu et al. (2021); Shen et al. (2022) in the special case where the
input space is Euclidean and the output space is the real line. In the static
setting, there are other geometric deep learning models defined between
Riemannian manifolds such as the geodesic convolutional model of (Masci et
al., 2015), but no approximation theory for these models has yet been
developed.

Mixture of experts (MoE) models such as DBRX MosaicAI (2024), Gemini (Google,
2024), Switch Transformers Fedus et al. (2022), Mixtral Jiang et al. (2024),
and many others (e.g. Saad et al. (2023); Chowdhury et al. (2023); Li et al.
(2024); Puigcerver et al. (2024); Saqur et al. (2024)) have taken a central
role in modern deep learning, due to the need to scale-up the model complexity
while maintaining a constant computational cost on the forward pass Kratsios
et al. (2024); Furuya et al. (2024). This is achieved via a gating mechanism
which routes any given input to one of a large number of “expert” neural
network models, which is then used to produce a prediction from that input.
Thus, only the gating network parameters and the selected “expert” neural
network are ever activated for that input. One can interpret our (hyper-
geometric network) as a mixture of infinitely many experts, and each of them
specializes in predicting at exactly one specific moment in time. The
hypernetwork in our HGN model acts as a gating mechanism which, given the
current point in time, routes the input to the corresponding expert at that
moment in time.

##  2 Geometric Background

We briefly overview some of the relevant notions from the geometry of non-
positively curved metric spaces, also called (global) NPC spaces or Hadamard
spaces. For further details we recommend the book by Bridson and Haefliger
(1999) as well as the results of Sturm (2003).

######  Definition 1 (NPC space)

A metric space $(N,d)$ is called global NPC space if it is complete and for
each $x_{0},x_{1}\in N$ there exists a point $y\in N$ with the property that
for all points $z\in N$

$d^{2}(z,y)\leq\frac{1}{2}d^{2}(z,x_{0})+\frac{1}{2}d^{2}(z,x_{1})-\frac{1}{4}d^{2}(x_{0},x_{1}).$

In particular, $(N,d)$ is contractible and hence simply connected.

A primary example of global NPC spaces is provided by Cartan-Hadamard
manifolds.

######  Proposition 2 (Manifolds (Proposition 3.1 in Sturm (2003))

Let $(N,g)$ be a Riemannian manifold and let $d$ be its induced Riemannian
distance. Then $(N,d)$ is a global NPC space if and only if it is complete,
simply connected and of non-positive (sectional) curvature.

We will describe some results on barycenters of probability measures on metric
spaces of NPC. First, we need to introduce the following notation; cfr. Sturm
(2003), Section 4. Let $(N,d)$ be a complete metric space. Denote by
$\mathcal{P}(N)$ the set of all probability measures $p$ on
$(N,\mathcal{B}(N))$ with separable support $\text{supp}(p)\subset N$. For
$1\leq\theta<\infty$, $\mathcal{P}^{\theta}(N)$ denotes the set of
$p\in\mathcal{P}(N)$ with $\int_{N}d^{\theta}(x,y)p(dy)<\infty$ for some
(hence all) $x\in N$, and $\mathcal{P}^{\infty}(N)$ denotes the set of all
$p\in\mathcal{P}(N)$ with bounded support. Denote by $\mathcal{W}_{1}$ the
$(L^{1}-)$ Wasserstein (or Kantorovich-Rubinstein) distance on
$\mathcal{P}^{1}(N)$. Let $(\Omega,\mathcal{F},\mathbb{P})$ be an arbitrary
probability space and $X:\Omega\rightarrow N$ a strongly666Given a measurable
space $(M,\mathcal{M})$ and $(N,d)$ a metric space, a map $f:M\rightarrow N$
is called strongly measurable iff it is $\mathcal{M}/\mathcal{B}(N)$
measurable and has separable range $f(M)$. measurable map. It defines a
probability measure $X_{*}\mathbb{P}\in\mathcal{P}(N)$, called the push
forward measure of $\mathbb{P}$ under $X$, by

$X_{*}\mathbb{P}[A]\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\mathbb{P}[X^{-1}(A)]=\mathbb{P}[\omega\in\Omega\,:\,X(\omega)\in
A],\quad\quad\forall A\in\mathcal{B}(N).$

In probabilistic language, a strongly measurable map $X:\Omega\rightarrow N$
is called $N$-valued random variable, the push forward measure
$X_{*}\mathbb{P}$ is called distribution of $X$ and denoted by
$\mathbb{P}_{X}$. In particular
$\mathbb{P}_{X}\in\mathcal{P}^{\theta}(N)\Leftrightarrow X\in
L^{\theta}(\Omega,N)$ where $L^{\theta}(\Omega,N)$ denotes the set of all
(strongly) measurable maps $f:\Omega\to N$ such that
$\int_{\Omega}d^{\theta}(x,f(\omega))\mathbb{P}(d\omega)<\infty$ for some (and
hence all) $x\in N$. We have the following

######  Proposition 3 (Existence of barycenters (Proposition 4.3 in Sturm
(2003)))

Let $(N,d)$ be a global NPC space and fix $y\in N$. For each
$q\in\mathcal{P}^{1}(N)$ there exists a unique point $z\in N$ which minimizes
the continuous function $z\mapsto\int_{N}[d^{2}(z,x)-d^{2}(y,x)]\,q(dx)$. This
point is independent of $y$: it is called barycenter (or, more precisely,
$d^{2}$-barycenter) of $q$ and denoted by

$\beta(q)=\underset{z\in
N}{\operatorname{argmin}}\int_{N}[d^{2}(z,x)-d^{2}(y,x)]\,q(dx).$

If $q\in\mathcal{P}^{2}(N)$, then $\beta(q)=\underset{z\in
N}{\operatorname{argmin}}\int_{N}d^{2}(z,x)\,q(dx)$.

As an example, if $q=\delta_{x_{0}}$, then trivially $\beta(q)=x_{0}$. For
$X\in L^{1}(\Omega,N)$ we define its expectation by

$\mathbb{E}[X]\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\underset{z\in
N}{\operatorname{argmin}}\,\mathbb{E}[d^{2}(z,X)-d^{2}(y,X)]=\beta(\mathbb{P}_{X}).$

That is, $\mathbb{E}[X]$ is the unique minimizer of the function
$z\mapsto\int_{N}[d^{2}(z,x)-d^{2}(y,x)]\mathbb{P}_{X}(\,dx)$ on $N$ (for each
fixed $y\in N$). We recall the following:

######  Theorem 4 (Fundamental contraction property; (Sturm, 2003, Theorem
6.3)))

For any given $p,q\,\in\mathcal{P}^{1}(N)$ it holds

$d(b(p),b(q))\leq\mathcal{W}_{1}(p,q).$ (2.1)

##  3 Problem Setting

In this paper, for a fixed $T\in\mathbb{N}_{+}$ we consider a discrete-time
stochastic process $X_{\cdot}=(X_{t})_{t=0}^{T}\subset\mathbb{R}^{d}$, which
evolves according to the following dynamics:

$X_{t+1}=X_{t}+\operatorname{Drift}(t,X_{[0:t]})+\operatorname{Diffusion}(t,X_{[0:t]},\mathbf{S}_{[0:t]})W_{t},\quad
t=0,\dots,T-1.$ (Volterra)

In the previous equation, $W_{\cdot}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(W_{t})_{t=0}^{T-1}$ is a Gaussian white-noise, i.e. an i.i.d.
collection of $\mathbb{R}^{d}$-valued standard normal random variables defined
on a probability space $(\Omega,\mathcal{F},\mathbb{P})$, and
$\mathbf{S}_{\cdot}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(S_{t})_{t=0}^{T-1}$ is a symmetric matrix-valued stochastic
process independent of $W_{\cdot}$, and $X_{0}=x_{0}\in\mathbb{R}^{d}$. In
this set-up, the direction and random fluctuations of the process $X_{\cdot}$,
depend on its realized path and the path of a latent (unobservable) stochastic
factor $\mathbf{S}_{\cdot}$

We focus on the case of discrete-time Volterra processes with non-singular
diffusion. Let $\operatorname{Sym}(d)$ (resp. $\operatorname{Sym}_{+}(d)$)
denote the set of $d\times d$ symmetric (resp. symmetric and positive
definite) matrices with real entries. We will call a function whose best
Lipschitz constant is at most $L$ an $L$-Lipschitz function. A Volterra
process with non-singular diffusion is a discrete-time stochastic process with
$\operatorname{Drift}:\mathbb{R}\times\mathbb{R}^{d(t+1)}\to\mathbb{R}^{d},$
and
$\operatorname{Diffusion}:\mathbb{R}\times\mathbb{R}^{d(t+1)}\times(\operatorname{Sym}(d))^{t+1}\to\operatorname{Sym}_{+}(d)$
satisfying:

$\displaystyle\operatorname{Drift}\big{(}t,x_{[0:t]}\big{)}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}$ $\displaystyle\sum_{r=0}^{t}\kappa(t,r)\,\mu(t,x_{r}),$ (3.1)
$\displaystyle\operatorname{Diffusion}\big{(}t,x_{[0:t]},s_{[0:t]}\big{)}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}$
$\displaystyle\exp\Biggl{(}\frac{1}{2}\,\sum_{r=0}^{t}\,\kappa(t,r)\,[\sigma(t,x_{r})+s_{r}]\Biggr{)},$

for all $t\in\mathbb{N}_{+}$, $x_{[0,t]}\in\mathbb{R}^{d(t+1)}$, and
$s_{[0,t]}\in(\operatorname{Sym}(d))^{t+1}$. Here $\exp$ denotes the matrix
exponential, $\mu:\mathbb{R}^{1+d}\rightarrow\mathbb{R}^{d}$,
$\sigma:\mathbb{R}^{1+d}\rightarrow\operatorname{Sym}(d)$ are
$L_{\mu}$-Lipschitz and $L_{\sigma}$-Lipschitz functions respectively and the
Volterra kernel $\kappa$ is a real-valued function satisyfing certain
properties which we will decribe below.

If the process $\mathbf{S}$ is $\mathbb{P}$-a.s. equal to $0$ then, we observe
that for each time $t\in\mathbb{N}_{+}$ and each path
$x_{[0:t]}\in\mathbb{R}^{(1+t)d}$ realized by the process $X_{\cdot}$ up to
time $t$, the distribution of $X_{t+1}$ conditioned on $X_{[0:t]}=x_{[0:t]}$
is Gaussian with mean $x_{t}+\operatorname{Drift}(t,x_{[0:t]})$ and covariance
equal to
$\operatorname{Diffusion}(t,x_{[0:t]},\boldsymbol{0})^{2}=\exp\big{(}\sum_{r=0}^{t}\,\kappa(t,r)\,\sigma(t,x_{r})\big{)}$.

The (discrete-time) Volterra kernel in Equation (3.1) determines the
persistence of $X_{\cdot}$’s memory on the distant past. A Volterra kernel is
defined to be a map $\kappa:\big{\\{}(t,r):\,t,r\in\mathbb{N}\mbox{ and }r\leq
t\big{\\}}\rightarrow[0,1]$ satisfying

$0<\sum_{r=0}^{t}\,\kappa(t,r)\leq 1\mbox{ and }\kappa(t,0)=0$

whenever $t>0$. We typically require that the Volterra kernel satisfies either
one of the following decay conditions

  1. i.

Exponential decay: For some $0<\alpha<1$ and $C>0$, $\kappa(T,r)\leq
C\,\alpha^{T-r}$ for all integers $0\leq r\leq T$;

  2. ii.

Polynomial decay: For some $\alpha<-1,C>0$, $\kappa(T,r)\leq C(T-r)^{\alpha}$
for all integers $0\leq r<T$.

We define the “tangential version” of the $\operatorname{Diffusion}$ map via
the matrix logarithm

$\operatorname{diffusion}(t,x_{[0:t]},s_{[0:t]})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\log\Big{(}\operatorname{Diffusion}(t,x_{[0:t]},s_{[0:t]})\Big{)}=\frac{1}{2}\,\sum_{r=0}^{t}\,\kappa(t,r)\,[\sigma(t,x_{r})+s_{r}],$

where the uniqueness of the choice of branch, in this case the logarithm, is
discussed in Appendix C.2. The intuition is that the inverse Riemannian
exponential map at $\mathcal{N}(0,I_{d})$ on the submanifold of
$\mathcal{N}_{d}$ of mean $0$ normal distributions coincides with the matrix
logarithm777See (3.10) in Proposition 7.. The stochastic diffusion factor
$s_{[0:t]}$ will be a realization of a $\operatorname{Sym}(d)$-valued
stochastic process $\mathbf{S}_{\cdot}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(\mathbf{S}_{t})_{t=0}^{T-1}$ on
$(\Omega,\mathcal{F},\mathbb{P})$. We require that the stochastic
$\operatorname{Sym}(d)$-valued process and the diffusive factor $\sigma$
satisfy the following conditions.

######  Assumption 5 (Uniformly bounded stochastic factor)

We assume there is a constant $R>0$ satisfying

$\sup_{t\leq T-1}\,\|\mathbf{S}_{t}\|_{F}\leq R,\quad\mathbb{P}-a.s.$ (3.2)

where $\|\cdot\|_{F}$ denotes the Frobenius norm of a matrix, and that there
is $M>0$ such that

$\sup_{t\leq T-1,\,x\in\mathbb{R}^{d}}\|\sigma(t,x)\|_{F}\leq M.$

###### Remark 6

Assumption 5 is required for the existence of a barycenter map. Though this
boundedness assumption is a bit strong, it can potentially be relaxed to an
exponentially decaying moment’s assumption in technical future work.

We now provide an example concerning discretization of stochastic delay
differential equations that fit into our framework.

######  Example 1 (Stochastic delay differential equations)

Considered a filtered probability space
$(\Omega,\mathcal{A},(\mathcal{A}_{u})_{u\in[0,\mathcal{T}]},\mathbb{P})$
satisfying the usual assumptions and carrying a standard one-dimensional
Wiener process $W(u)$ ($u$ denotes the time-variable here). Let
$f:\mathbb{R}\to\mathbb{R}$ and $g:\mathbb{R}\to\mathbb{R}_{>0}$ be bounded
and Lipschitz continuous functions. Fix a time delay $\tau>0$ and let
$\psi(u)$ be an $\mathcal{A}_{0}$-measurable $C([-\tau,0],\mathbb{R})$ valued
random variable satisfying
$\mathbb{E}[\left\lVert\psi\right\rVert_{\infty}^{2}]<\infty$. Consider the
following stochastic delay differential equation

$X(u)=X(0)+\int_{0}^{u}f(X(s-\tau))ds+\int_{0}^{u}g(X(s-\tau))dW(s)$

for $0\leq u\leq\mathcal{T}$ and with $X(u)=\psi(u)$ for $-\tau\leq u\leq 0$
(for further details see e.g. Buckwar (2000)). We consider a grid with a
uniform step $h$ on the interval $[0,\mathcal{T}]$ and
$h=\mathcal{T}/N,u_{t}=th,t=0,\dots,N$, and assume that there is an integer
$N_{\tau}$ such that $\tau=hN_{\tau}$. The Euler-Maruyama scheme is

$\tilde{X}_{t+1}=\tilde{X}_{t}+hf(\tilde{X}_{t-N_{\tau}})+g(\tilde{X}_{t-N_{\tau}})\Delta
W_{t+1},\quad 0\leq t\leq N-1$ (3.3)

where $\tilde{X}_{t-N_{\tau}}:=\psi(u_{t}-\tau)$ for $t-N_{\tau}\leq 0$, and
$\Delta W_{t+1}:=W((t+1)h)-W(nh)\sim N(0,h)$ are independent.

Consider an identically zero stochastic factor $\mathbf{S}$, setting
$\kappa(t,r)=0$ for $r\neq t-N_{\tau}\geq 0$ and $\kappa(t,r)=1$ for
$r=t-N_{\tau}\geq 0$ , $\mu(x):=hf(x)$ and $\sigma(x)=2h\log g(x)$, from which
we see that the discretization in (3.3) falls within our framework.
Nice applications involving stochastic differential are linked to, e.g.,
problems of optimal advertising under uncertainty for the introduction of a
new product to the market (Gozzi et al. (2005)), some infinite-dimensional
variants of the Black-Scholes equation for the fair price of an option
(Fuhrman et al. (2010) and Chang and Youree (2007)), stochastic control
problem with delay arising in a pension fund models (Federico (2011)).

###  3.1 A Globally NPC perturbation of the Fisher geometry for Gaussian
measures

Our target Riemannian (statistical) manifold of interest is the set
$\mathcal{N}_{d}$ of non-singular Gaussian measures of a given dimension $d$,
and formally defined below, endowed with a perturbation of the Fisher
Riemannian metric (Skovgaard (1984)), whose geometry we now review. Let
$\mathcal{N}_{d}$ be the collection of non-degenerate normal distributions of
fixed dimension $d$, i.e.

$\mathcal{N}_{d}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\\{\mathcal{N}_{d}(\mathfrak{m},\Sigma)\,:\,\mathfrak{m}\in\mathbb{R}^{d}\mbox{
and }\Sigma\in\operatorname{Sym}_{+}(d)\\},$ (3.4)

where, we recall, $\operatorname{Sym}_{+}(d)$ is the set of $d\times d$
positive definite symmetric matrices. Thus,
$\mathcal{N}_{d}(\mathfrak{m},\Sigma)$ has Lebesgue density given by

$p_{\mathfrak{m},\Sigma}(u)=\frac{e^{-\frac{1}{2}(u-{\mathfrak{m}})^{\top}\Sigma^{-1}(u-{\mathfrak{m}})}}{\sqrt{(2\pi)^{d}\det(\Sigma)}},\quad
u\in\mathbb{R}^{d}.$

Through the identification

$\varphi(\mathcal{N}_{d}(\mathfrak{m},\Sigma))=((\mathfrak{m}_{i})_{i=1}^{d},(\Sigma_{ij})_{i\leq
j}),$ (3.5)

where
$\mathfrak{m}=(\mathfrak{m}_{1},\dots,\mathfrak{m}_{d})^{\top},\Sigma=(\Sigma_{ij})_{i,j=1}^{d}$,
we see that $\mathcal{N}_{d}$ is isomorphic to an open subset
$\Theta\subset\mathbb{R}^{m}$, $m=(d^{2}+3d)/2$. Hence, $\mathcal{N}_{d}$ can
be considered as a differentiable manifold of dimension $m$, with
$(\mathcal{N}_{d},\varphi)$ as a global chart. We define a basis, in this
coordinate system, of the set of vector fields on $\mathcal{N}_{d}$

$\frac{\partial}{\partial\mathfrak{m}_{i}},\quad
i=1,\dots,d,\quad\frac{\partial}{\partial\Sigma_{ij}},\quad i,j=1,\dots
d,\,i\leq j.$

We identify these vector fields with the vectors and matrices

$\frac{\partial}{\partial\mathfrak{m}_{i}}\leftrightarrow
e_{i}\in\mathbb{R}^{d},\quad i=1,\dots,d,$
$\frac{\partial}{\partial\Sigma_{ij}}\leftrightarrow
E_{ij}\in\operatorname{Sym}_{+}(d),\quad i,j=1,\dots d,\,i\leq j,$

where $(e_{i})_{i}$ is the canonical basis of $\mathbb{R}^{d}$ and

$E_{ij}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\begin{cases}I_{(i,j)},\quad\quad\quad\quad i=j\\\
I_{(i,j)}+I_{(j,i)},\quad i\neq j,\end{cases}$

where $I_{(i,j)}$ denotes the $d\times d$ matrix whose $(i,j)$-$th$ entry is
equal to 1, while all the remaining entries are zero. Therefore, each smooth
vector field $A$ on $\mathcal{N}_{d}$ can be decomposed as

$A=\sum_{i=1}^{d}F_{i}e_{i}+\sum_{\overset{i,j=1}{i\leq j}}^{d}G_{ij}E_{ij}$

where $F_{i},G_{ij}:\mathcal{N}_{d}\to\mathbb{R}$ are smooth functions. We
shall call the two terms in this last equation the components of $A$ in the
$\mathfrak{m}$\- and $\Sigma$-direction, respectively. The space
$\mathcal{N}_{d}$ can be turned into a Riemannian manifold by using the Fisher
information as its Riemannian metric. More precisely, we define at an
arbitrary point $(\mathfrak{m},\Sigma)$ the following Riemannian metric
$\mathcal{I}^{F}$

$\begin{split}&\mathcal{I}^{F}\left(\frac{\partial}{\partial\mathfrak{m}_{i}},\frac{\partial}{\partial\mathfrak{m}_{j}}\right)\equiv\mathcal{I}^{F}(e_{i},e_{j})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\Sigma^{ij},\\\
&\mathcal{I}^{F}\left(\frac{\partial}{\partial\mathfrak{m}_{i}},\frac{\partial}{\partial\Sigma_{kl}}\right)\equiv\mathcal{I}^{F}(e_{i},E_{kl})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}0,\\\
&\mathcal{I}^{F}\left(\frac{\partial}{\partial\Sigma_{ij}},\frac{\partial}{\partial\Sigma_{kl}}\right)\equiv\mathcal{I}^{F}(E_{ij},E_{kl})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\frac{1}{2}\operatorname{tr}(\Sigma^{-1}E_{ij}\Sigma^{-1}E_{kl})=\Sigma^{il}\Sigma^{jk}+\Sigma^{ik}\Sigma^{jl},\end{split}$
(3.6)

for $i,j,k,l=1,\dots,d,\,i\leq j$ and $k\leq l$; where
$\Sigma^{-1}=(\Sigma^{ij})_{i,j=1}^{d}$ denotes the inverse of $\Sigma$, and
$\operatorname{tr}$ is the trace operator. From the last equation, by the
bilinearity of the metric $\mathcal{I}^{F}$, we see that

$\mathcal{I}^{F}(X,Y)=\frac{1}{2}\operatorname{tr}(\Sigma^{-1}X\Sigma^{-1}Y),\quad
X,Y\in\operatorname{Sym}(d).$ (3.7)

Though the Fisher metric is the canonical Riemannian metric on the smooth
manifold $\mathcal{N}_{d}$, in the sense of Chentsov’s Theorem (Csiszár, 1984;
Dowty, 2018) from information geometry, its geometry is not amenable to
intrinsically averaging in dimension above $1$. When $d>1$ the Fisher geometry
on $\mathcal{N}_{d}$ exhibits positive sectional curvature along tangent
planes in the directions spanned by
$\frac{\partial}{\partial\mathfrak{m}_{i}}$ and
$\frac{\partial}{\partial\mathfrak{m}_{j}}$ for distinct $i,j=1,\dots,d$ (see
(Skovgaard, 1984, Theorem 2.2)). These directions of positive sectional
curvature do not exist in univariate setting, since no such tangent plane
exist when $d=1$, and the only tangent plane at any point is generated by
$\frac{\partial}{\partial\mathfrak{m}_{1}}$ and
$\frac{\partial}{\partial\Sigma_{1,1}}$ whose sectional curvature equals
$-1/2$; this leads to the well-known identification of the Fisher geometry on
$\mathcal{N}_{1}$ with that of the hyperbolic plane.

As we have learnt above from the results of Sturm (2003), lack of globally NPC
on the Fisher geometry of $\mathcal{N}_{d}$ for $d>1$ makes it unsuitable to
intrinsic averaging. Furthermore, there are no known closed-form expressions
for the geodesic distance and for the Riemannian exponential map about any
point for Fisher geometry on $\mathcal{N}_{d}$, in dimension larger than $1$.
Both these limitations make the Fisher geometry on $\mathcal{N}_{d}$
incompatible with the available tools for probability theory on NPC spaces and
approximation theory between Riemannian manifold (Kratsios and Papon (2022)).
In light of this, we perturb the Fisher metric (3.6) as follows

$\begin{split}&\mathfrak{I}\left(\frac{\partial}{\partial\mathfrak{m}_{i}},\frac{\partial}{\partial\mathfrak{m}_{j}}\right)=\boldsymbol{\delta}_{ij},\\\
&\mathfrak{I}\left(\frac{\partial}{\partial\mathfrak{m}_{i}},\frac{\partial}{\partial\Sigma_{kl}}\right)=0,\\\
&\mathfrak{I}\left(\frac{\partial}{\partial\Sigma_{ij}},\frac{\partial}{\partial\Sigma_{kl}}\right)=\frac{1}{2}\operatorname{tr}(\Sigma^{-1}E_{ij}\Sigma^{-1}E_{kl})=\Sigma^{il}\Sigma^{jk}+\Sigma^{ik}\Sigma^{jl},\end{split}$
(3.8)

for $i,j,k,l=1,\dots,d,\,i\leq j$ and $k\leq l$, where $\boldsymbol{\delta}$
denotes the standard Riemannian metric of $\mathbb{R}^{d}$. Besides being a
non-positively curved perturbation of the geometry induced by Fisher metric,
which also has a closed-form distance function, the Riemannian metric
$\mathfrak{J}$ has a second simple geometric interpretation. Namely,
$(\mathcal{N}_{d},\mathfrak{J})$ is isometric as a Riemannian manifold to
$\mathbb{R}^{d}\times\operatorname{Sym}_{+}(d)$ equipped with the product
between the Euclidean metric $\boldsymbol{\delta}$ on $\mathbb{R}^{d}$ and the
well-studied affine-invariant metric on $\operatorname{Sym}_{+}(d)$, see Smith
(2005); Pennec et al. (2006); Schiratti et al. (2017); Brooks et al. (2019);
Criscitiello and Boumal (2022). Some of the key geometric properties of
$(\mathcal{N}_{d},\mathfrak{I})$ can now be recorded. Below,
$d_{\mathfrak{I}}$ will denote the distance function induced by
$\mathfrak{I}$.

######  Proposition 7 (The geometry of $(\mathcal{N}_{d},\mathfrak{I})$)

For any $d\geq 1$, $(\mathcal{N}_{d},d_{\mathfrak{I}})$ is a complete, path-
connected and simply connected, Riemannian manifold with sectional curvatures
in $[-1/2,0]$. Therefore, $(\mathcal{N}_{d},d_{\mathfrak{I}})$ is a global NPC
space, and for any
$\mathcal{N}_{d}(\mathfrak{m}_{0},\Sigma_{0}),\mathcal{N}_{d}(\mathfrak{m}_{1},\Sigma_{1})\in\mathcal{N}_{d}$,
$d_{\mathfrak{I}}$ admits the closed-form expression

$d_{\mathfrak{I}}\big{(}\mathcal{N}_{d}(\mathfrak{m}_{0},\Sigma_{0}),\mathcal{N}_{d}(\mathfrak{m}_{1},\Sigma_{1})\big{)}=\left[\left\lVert\mathfrak{m}_{0}-\mathfrak{m}_{1}\right\rVert^{2}_{2}+\frac{1}{2}\sum_{i=1}^{d}\,\ln(\lambda_{i})^{2}\right]^{1/2},$
(3.9)

where $0<\lambda_{1}\leq\dots\leq\lambda_{d}$ are the eigenvalues of
$\Sigma_{0}^{-1}\Sigma_{1}$.

For any $\mathfrak{m}\in\mathbb{R}^{d}$, $\Sigma\in\operatorname{Sym}_{+}(d)$,
and any $(\tilde{\mathfrak{m}},X)\in
T_{\mathcal{N}_{d}(\mathfrak{m},\Sigma)}\mathcal{N}_{d}\cong\mathbb{R}^{d}\times\operatorname{Sym}(d)$
the Riemannian exponential map is given by

$\operatorname{Exp}_{\mathcal{N}_{d}(\mathfrak{m},\Sigma)}(\tilde{\mathfrak{m}},X)=\mathcal{N}\left(\mathfrak{m}+\tilde{\mathfrak{m}},\Sigma^{1/2}\operatorname{exp}\big{(}\Sigma^{-1/2}X\Sigma^{-1/2}\big{)}\Sigma^{1/2}\right),$
(3.10)

where $T_{\mathcal{N}_{d}(\mathfrak{m},\Sigma)}\mathcal{N}_{d}$ denotes the
tangent space of $\mathcal{N}_{d}$ at the point
$\mathcal{N}_{d}(\mathfrak{m},\Sigma)$.

Proof  See Appendix D.1.1.

####  3.1.1 Comparison with other geometries on $\mathcal{N}_{d}$

We may further motivate the proposed Riemannian geometry on $\mathcal{N}_{d}$
by contrasting it with other standard geometries on the set $\mathcal{N}_{d}$,
occurring in optimal transport and information geometry. Our discussion is
summarized concisely in Table 1.

Table 1: Comparison of topological, metric, and computational properties
between different geometries of non-singular Gaussian measures from
information geometry and optimal transport.

Geometry |  Complete |  NPC (Barycenters) |  Closed-Form Dist. |  Isometric Copies $(\mathcal{N}_{d}^{\mathfrak{m}},\mathcal{I})$ | Reference
---|---|---|---|---|---
$(\mathcal{N}_{d},d_{\mathfrak{J}})$ |  $\good$ |  $\good$ |  $\good$ |  $\good$ | Equation (3.8)
$(\mathcal{N}_{d},d_{\mathcal{I}})$ |  $\good$ |  $\bad$ ($\good$ $d=1$) |  $\bad$ |  $\good$ | Skovgaard (1984)
$(\mathcal{N}_{d},\mathcal{W}_{2})$ |  $\bad$ |  ? |  $\good$ |  $\bad$ | Takatsu (2011)
$(N_{d},g_{SL})$ |  $\good$ |  $\good$ |  $\good$ |  $\bad$ | Lovrić et al. (2000)

In what follows, we linearly identify the spaces via the symmetrization map
$\operatorname{sym}$ defined by

$\operatorname{sym}:\mathbb{R}^{d(d+1)/2}\rightarrow\operatorname{Sym}(d),\quad
v\mapsto\begin{pmatrix}v_{1}&v_{2}&\dots&v_{d}\\\ v_{2}&v_{d+1}&&.\\\
\vdots&&\ddots&\vdots\\\ v_{d}&\dots&&v_{d(d+1)/2}\end{pmatrix}$ (3.11)

whose inverse is often denoted by $\operatorname{vec}$ (see (Petersen et al.,
2008, Section 10.2.2)). We observe that the linear identification is not an
isometry, but it is a bi-Lipschitz map, when $\operatorname{Sym}(d)$ is
equipped with Fröbenius inner-product
$\langle\Sigma_{1},\Sigma_{2}\rangle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{tr}(\Sigma_{1}\Sigma_{2})$ and
$\mathbb{R}^{d(d+1)/2}$.

Table 1 compares the Riemannian geometry $(\mathcal{N}_{d},\mathfrak{J})$ with
alternative geometries found in the literature. Specifically, it contrasts
$(\mathcal{N}_{d},\mathfrak{J})$ with the standard Fisher information geometry
as studied in Skovgaard (1984), the Riemannian geometry inducing the
$2$-Wasserstein metric ($\mathcal{W}_{2}$) on $\mathcal{N}_{d}$ as presented
in Takatsu (2011), and a non-positive curvature geometry on $\mathcal{N}_{d}$,
denoted $g_{SL}$, constructed in Lovrić et al. (2000) by identifying
$\mathcal{N}_{d}$ with a subset of a special linear group, modulo the action
of a special orthogonal group. The columns of Table 1 should be self-
explanatory apart from the penultimate column. Here we say that the geometry
admits isometric copies of $(\mathcal{N}_{d},\mathcal{I})$ if for each
$\mathfrak{m}\in\mathbb{R}^{d}$, the submanifold
$\mathcal{N}^{\mathfrak{m}}_{d}\subset\mathcal{N}_{d}$ whose points have mean
$\mathfrak{m}$ is isometric to $(\mathcal{N}_{d}^{\mathfrak{m}},\mathcal{I})$.

Finally, we remark that the the manifolds $(\mathcal{N}_{d},W_{2})$,
$(\mathcal{N}_{d},\mathfrak{I})$, and $(\mathcal{N}_{d},\mathcal{J})$ are
diffeomorphic; in particular, their topologies coincide. Nevertheless, we
emphasize that all three structures are different as Riemannian manifolds and
as metric spaces.

##  4 Gaussian random projections

The law of the process $X_{\cdot}$ in (Volterra), conditioned on any realized
path up to a given time $t\in\mathbb{N}$, is an arbitrary probability
distribution on $\mathbb{R}^{d}$ with a priori no favourable structure. For
instance, the space of probability distribution on $\mathbb{R}^{d}$, e.g. with
the Lévy-Prokhorov metric, see e.g. (Bogachev, 2007, Theorem 8.3.2), is
neither a finite-dimensional as a topological space nor is it a vector space
(even though it is convex). This lack of structure implies that its elements
cannot be approximated by simpler probability measures, such as finitely
supported measures, without facing the curse of dimensionality (see, e.g.
(Graf and Luschgy, 2000, Theorem 6.2)).

##### Intuitive idea

However, if we had access to the unobserved realized path of the “stochastic
factor process” $\mathbf{S}_{[0:t]}$, we could leverage the fact that
$X_{t+1}$ is Gaussian when conditioned on the realized path $X_{[0:t]}$ and
belongs to $\mathcal{N}_{d}$, the smooth manifold of non-singular
$d$-dimensional Gaussian measures. This motivates the family of probability
measures
$\\{\Lambda_{x_{[0:t]},s_{[0:t]}};\,x_{[0:t]}\in\mathbb{R}^{(1+t)d},\,s_{[0:t]}\in\operatorname{Sym}(d)^{1+t}\\}\subset\mathcal{N}_{d}$,
defined by

$\Lambda_{x_{[0:t]},s_{[0:t]}}[\cdot]\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\mathbb{P}[X_{t+1}\in\cdot|X_{[0:t]}=x_{[0:t]},\,\mathbf{S}_{[0:t]}=s_{[0:t]}\big{]},\quad
x_{[0:t]}\in\mathbb{R}^{(1+t)d},\,s_{[0:t]}\in\operatorname{Sym}(d)^{1+t}$
(4.1)

but since the process $\mathbf{S}_{\cdot}$ is latent, we need to reintroduce
its “random effects” into the family of probability measures
$\Lambda_{x_{[0:t]},s_{[0:t]}}$. Traditionally, this is achieved by
integrating against the law of $\mathbf{S}_{[0:t]}$. Informally, we instead
reinsert the observed path $\mathbf{S}_{[0:t]}$ as the conditional state
$s_{[0:t]}$, i.e.:

$\mathbb{Q}_{x_{[0:t]}}(\omega)\overset{\text{informal}}{=}\Lambda_{x_{[0:t]},\mathbf{S}_{[0:t]}(\omega)},\quad\omega\in\Omega.$
(4.2)

A rigorous construction of (4.2), formalized below in Definition 8, hinges on
our Riemannian geometry on $\mathcal{N}_{d}$ and it enjoys Lipschitz
stability, meaning that it is a Lipschitz function of path realized by the
process $X_{\cdot}$. In some cases, the dependence of $\mathbb{Q}_{x_{[0:t]}}$
on the realized path of $X_{\cdot}$ is even smooth, meaning that deep learning
model can approximate it with a feasible number of parameters. This is
generally not true with geometric deep learning models on non-smooth metric
spaces, see Acciaio et al. (2023). The results of Sturm (2003) guarantee that
there is a unique point $\Pi_{x_{[0:t]}}$ in $\mathcal{N}_{d}$ which best
represents the intrinsic average behaviour of $\mathbb{Q}_{x_{[0:t]}}$

$\Pi_{x_{[0:t]}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{argmin}_{u\in\mathcal{N}_{d}}\,\int_{v\in\mathcal{N}_{d}}d_{\mathfrak{J}}^{2}(u,v)\text{Law}(\mathbb{Q}_{x_{[0:t]}}\in
dv),$ (4.3)

where $d_{\mathfrak{J}}$ is the geodesic distance induced by our perturbation
of the Fisher Riemannian metric on $\mathcal{N}_{d}$. Unlike information
projections, the map defined by (4.3) enjoys the following stability
properties and can even be smooth in certain cases.

######  Informal Results 1 (Gaussian projections Theorem)

Under mild conditions, if the $\operatorname{Drift}$ and
$\operatorname{Diffusion}$ maps in (Volterra) are $L$-Lipschitz, then the map
$x_{[0:t]}\mapsto\Pi_{x_{[0:t]}}$ is well-defined, $\mathcal{O}(L)$-Lipschitz,
and it solves the following minimization problem

$\Pi_{x_{[0:t]}}=\underset{\mu\in\mathcal{N}_{d}}{\operatorname{argmin}}\,\mathbb{E}\big{[}d^{2}_{\mathfrak{J}}(\mu,\mathbb{Q}_{x_{[0:t]}})\big{]}$

Furthermore, if the latent factor process $\mathbf{S}_{\cdot}=0$ then

$\Pi_{x_{[0:t]}}=\mathbb{P}[X_{t+1}\in\cdot|X_{[0:t]}=x_{[0:t]}]$

for each time $t\in\mathbb{N}_{+}$ and every realized path
$x_{[0:t]}\in\mathbb{R}^{(1+t)d}$.

We now formalize our notion of projection of the conditional law of the
stochastic process (Volterra), called a Gaussian random projection. The
terminology should not be understood in the sense of the Johnson-Lindenstrauss
lemma, see (Matoušek, 2002, Theorem 12.2.1), but rather in the sense of Ohta
(2009). That is, as a Lipschitz measure-valued function with the property that
a certain class of inputs are fixed in certain sense (see Corollary 10).

###  4.1 Encoding the conditional distribution of a Volterra process in
$\mathcal{N}_{d}$

Our argument centers around the geometry of the regular conditional
distribution (RCD) function sending any $(x_{[0:t]},s_{[0:t]})$ to

$\mathbb{P}[X_{t+1}\in\cdot|X_{[0:t]}=x_{[0:t]},\mathbf{S}_{[0:t]}=s_{[0:t]}].$

Our first observation is to note that the RCD’s image lies in
$\mathcal{N}_{d}$.

The map
$\tilde{\varphi}:\mathcal{N}_{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d(d+1)/2}$,
defined by
$\mathcal{N}_{d}(\mathfrak{m},\Sigma)\stackrel{{\scriptstyle\tilde{\varphi}}}{{\mapsto}}(\mathfrak{m},\operatorname{vec}(\log(\Sigma)))$,
is a smooth chart. This is because it is the composition of the maps
$1_{\mathbb{R}^{d}}\times\operatorname{vec}$,
$\mathbb{R}^{d}\times\operatorname{Sym}(d)$ , $\log$ 888See (Hilgert and Neeb,
2012, Proposition 3.3.5). , $1_{\mathbb{R}^{d}}\times\log$ and $\varphi$,
defined in (3.5), each of which are diffeomorphisms. Thus,

$\displaystyle\psi:\mathbb{R}^{d(d+3)/2}$
$\displaystyle\rightarrow\mathcal{N}_{d}$ (4.4) $\displaystyle(u,v)$
$\displaystyle\mapsto\tilde{\varphi}^{-1}(u,v)=\mathcal{N}_{d}(u,\exp(\operatorname{sym}(v)))$

defines a smooth global chart on $\mathcal{N}_{d}$. For every
$t\in\mathbb{N}_{+}$, the maps

$\displaystyle\mathbb{R}^{(1+t)d}\ni x_{[0:t]}$
$\displaystyle\mapsto\operatorname{Drift}(t,x_{[0:t]})\in\mathbb{R}^{d}$
$\displaystyle\mathbb{R}^{(1+t)d}\times\operatorname{Sym}(d)^{(1+t)}\ni(x_{[0:t]},s_{[0:t]})$
$\displaystyle\mapsto\operatorname{Diffusion}(t,x_{[0:t]},s_{[0:t]})\in\operatorname{Sym}_{+}(d)$

are locally Lipschitz continuous. Thus, for every $t\in\mathbb{N}_{+}$, the
map

$\mathbb{R}^{(1+t)d}\times\mathbb{R}^{(1+t)d(d+1)/2}\ni(x_{[0:t]},y_{[0:t]})\mapsto\operatorname{Diffusion}(t,x_{[0:t]},\operatorname{sym}(y_{[0:t]}))\in\operatorname{Sym}_{+}(d),$

where $\operatorname{sym}$ is applied componentwise, is locally Lipschitz
continuous. Thus, for every $t\in\mathbb{N}_{+}$ the map
$\Psi_{t}:\mathbb{R}^{(1+t)d}\times\mathbb{R}^{(1+t)\,d(d+1)/2}\rightarrow\mathcal{N}_{d}$
defined for any $x_{[0:t]}\in\mathbb{R}^{(1+t)d}$ and
$y_{[0:t]}\in\mathbb{R}^{(1+t)d(d+1)/2}$ by

$\displaystyle\Psi_{t}(x_{[0:t]},y_{[0:t]})$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\psi\big{(}x_{t}+\operatorname{Drift}(t,x_{[0:t]}),\operatorname{vec}\circ\log\circ\operatorname{Diffusion}(t,x_{[0:t]},\operatorname{sym}(y_{[0:t]}))\big{)}$
(4.5)
$\displaystyle=\mathcal{N}_{d}\left(x_{t}+\operatorname{Drift}(t,x_{[0:t]}),\exp\left(\sum_{r=0}^{t}\kappa(t,r)\left[\,\sigma(t,x_{r})+\operatorname{sym}(y_{r})\right]\right)\right)$

must also be locally-Lipschitz when $\mathcal{N}_{d}$ is equipped with the
Riemannian distance function $d_{\mathfrak{J}}$. Note that
$\operatorname{Diffusion}(t,x_{[0:t]},\operatorname{sym}(s_{[0:t]}))$ is
positive definite. Therefore we can assume that its eigenvalues are such that
$0<\lambda_{1}\leq\dots\leq\lambda_{d}$. Furthermore, if $\mu$ and $\sigma$
are smooth, then each $\Psi_{t}$ is also smooth by composition. Consequently,
we may rigorously define the $\mathcal{N}_{d}$-valued stochastic process as
the composition of $x_{[0:t]}$ and $\mathbf{S}_{[0:t]}$ and the locally
Lipschitz continuous map $\Psi_{t}$, namely
$\Psi_{t}(x_{[0:t]},\operatorname{vec}\circ\mathbf{S}_{[0:t]})$.

We note that $\|e^{A}\|_{op}\leq e^{\|A\|_{op}}$ for any $d\times d$ matrix
$A$; hence, for some universal constants $c_{1},c_{2}>0$ depending only $d$,
$\|e^{A}\|_{F}\leq c_{1}e^{c_{2}\|A\|_{F}}$; therefore the boundedness of
$S_{\cdot}$ and the image of $\sigma$ in Assumption 5 imply that the law of
$\Psi_{t}(x_{[0:t]},\operatorname{vec}\circ\mathbf{S}_{[0:t]})$ is compactly
supported on $\mathcal{N}_{d}$. Therefore,
$\Psi_{t}(x_{[0:t]},\operatorname{vec}\circ\mathbf{S}_{[0:t]})$ is integrable
over the global NPC space $(\mathcal{N}_{d},d_{\mathfrak{J}})$, in the sense
that it belongs to $L^{1}(\Omega,\mathcal{N}_{d})$, and the barycenter of its
law is well-defined .

######  Definition 8 (Gaussian random projection of $X_{\cdot}$)

Given a Volterra process $X_{\cdot}$ in (3.1), we define the Gaussian Random
projection of $X_{\cdot}$ as the map, sending any
$x_{[0:t]}\in\mathbb{R}^{(1+t)d}$ to the measure
$\Pi_{x_{[0:t]}}\in\mathcal{N}_{d}$ defined by

$\displaystyle\Pi_{x_{[0:t]}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\beta\left(\operatorname{Law}\left(\Psi_{t}\left(x_{[0:t]},\operatorname{vec}(\mathbf{S}_{[0:t]})\right)\right)\right)\in\mathcal{N}_{d}$
(4.6)

where $\beta$ is given in Proposition 3. For the sake of brevity, we will
henceforth write $\Pi_{x_{[0:t]}}=\beta(\mathbb{Q}_{x_{[0:t]}})$, where
$\mathbb{Q}_{x_{[0:t]}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{Law}\left(\Psi_{t}(x_{[0:t]},\operatorname{vec}(\mathbf{S}_{[0:t]}))\right)$.

In particular, notice that if we have $\mathbf{S}_{[0:t]}=s_{[0:t]}$, then

$\Pi_{x_{[0:t]}}=\mathcal{N}_{d}(x_{t}+\operatorname{Drift}(t,x_{[0:t]}),\operatorname{Diffusion}(t,x_{[0:t]},s_{[0:t]})^{2})$

where
$\operatorname{Diffusion}(t,x_{[0:t]},s_{[0:t]})^{2}=\exp\left(\sum_{r=0}^{t}\kappa(t,r)\left[\,\sigma(t,x_{r})+s_{r}\right]\right)$.
We now have the following.

######  Theorem 9 (Optimality and Lipschitz Stability of the Gaussian Random
Projection)

Fix $t\in\mathbb{N}$ and a compact
$\mathcal{K}_{[0:t]}\subseteq\mathbb{R}^{(1+t)\,d}$. If (3.1) and Assumption 5
hold then, $\Pi:x_{[0:t]}\mapsto\beta(\mathbb{Q}_{x_{[0:t]}})$ is well-defined
and Lipschitz and

$\Pi_{x_{[0:t]}}=\underset{{u\in\mathcal{N}_{d}}}{\operatorname{argmin}}\,\int_{v\in\mathcal{N}_{d}}d_{\mathfrak{J}}^{2}(u,v)\mathbb{Q}_{x_{[0:t]}}(dv).$

Furthermore, the Lipschitz constant of $\Pi_{\cdot}$ depends only on $t$ and
$R$, as defined in Assumption 5.

Proof  As regards the Lipschitz part, we refer to Proposition E.5. The
remainder of the proof follows from the discussion near (4.6).

The smooth structure on $\mathcal{N}_{d}$ has the advantage of granting
meaning to the smoothness of the map $\Pi_{\cdot}$ when $\mu$ and $\sigma$ are
smooth and $\mbox{{S}}_{[0:t]}$ is almost surely constant. Smoothness is
beneficial from an approximation perspective since neural networks can
effectively approximate smooth Lipschitz functions between Riemannian
manifolds. This is demonstrated in Theorem 13 below. In contrast, even in the
Euclidean setting, non-smooth Lipschitz functions between Euclidean spaces on
simple, compact subsets cannot be accurately approximated by feasible neural
networks, as shown in (Yarotsky, 2017, Theorem 5).

The following Corollary provides a natural sufficient condition for the
smoothness of $\Pi_{\cdot}$. Namely, if the maps $\mu$ and $\sigma$, which
define the drift and diffusion maps of the Volterra process, are smooth, and
that $\mathbf{S}_{\cdot}$ is constantly $0$ then $\Pi_{\cdot}$ is also smooth.

######  Proposition 10 (Smoothness)

Suppose that $\mathbf{S}_{t}=0$ for every $t\in\mathbb{N}$, Assumption 5 holds
and that $\mu$ and $\sigma$ are of class $C^{\infty}$. Then for every
$t\in\mathbb{N}$, the following map is $C^{\infty}$

$\displaystyle\Pi_{\cdot}:\,(\mathbb{R}^{(1+t)d},\|\cdot\|)$
$\displaystyle\rightarrow(\mathcal{N}_{d},d_{\mathfrak{J}})$ $\displaystyle
x_{[0:t]}$ $\displaystyle\mapsto\Pi_{x_{[0:t]}}.$

Proof  See Appendix E.

Next, we study how the Volterra kernel dictates the rate at which the Gaussian
random projection forgets the distant history of any realized path of a
Volterra process.

###  4.2 Memory decay of $\Pi_{\cdot}$ given a decaying Volterra kernel

In order for the map $x_{[0:t]}\mapsto\Pi_{x_{[0:t]}}$ to be approximable by
any model which only processes the path realized by $X_{\cdot}$ up to the
current time, we need that the stochastic Volterra process does not exhibit a
persistent memory of its distant past. As one may expect, a stochastic
Volterra process exhibits no persistent long-term memory if its Volterra
kernel does not place large weights on long-past states realized by the
process $X_{\cdot}$. Similar phenomena have been observed in the RNN
literature when approximating linear state-space systems with exponentially
decaying impulse response in (Hutter et al., 2022) or in (Li et al., 2022b,
Theorem 11) when quantifying the approximability of linear functionals by
linear continuous-time RNNs. Our next result shows that if the Volterra
kernels decay at least polynomially, then one has a quantitative analogue of
the fading memory property in reservoir computing (Grigoryeva and Ortega,
2018, Definition 3), and of the approximable complexity property of (Acciaio
et al., 2023, Definition 4.9).

######  Theorem 11 (Memory Decay Rate of Gaussian Projection of Volterra
Process)

Suppose Assumption 5 holds and the Volterra-kernel $\kappa$ either satisfies
(i) or (ii). For every $T\in\mathbb{N}_{+}$, each compact
$\mathcal{K}_{[0:T]}\subseteq\mathbb{R}^{(1+T)d}$, there is a constant
$C_{M,R,\mu}>0$ such that for every pair of integers $\operatorname{{H}}$ and
$t$, with $0\leq\operatorname{{H}}<t\leq T$, there is a Lipschitz function
$\tilde{f}^{(t,\operatorname{{H}})}:(\mathbb{R}^{(1+\operatorname{{H}})d},\|\cdot\|)\rightarrow(\mathcal{N}_{d},d_{\mathfrak{J}})$
satisfying

$d_{\mathfrak{J}}\big{(}\Pi_{x_{[0:t]}},\tilde{f}^{(t,\operatorname{{H}})}(x_{[t-\operatorname{{H}}:t]})\big{)}\leq
C_{M,R,\mu}\,\operatorname{diam}(\mathcal{K}_{[0:T]})\,\begin{cases}\frac{C\alpha}{\alpha-1}\,(\alpha^{t}-\alpha^{\operatorname{{H}}})&\mbox{if
(i) holds}\\\
C\,(\operatorname{{H}}+1)^{\alpha}\,(t-\operatorname{{H}})&\mbox{if (ii)
holds}\end{cases}$

for every $x_{[0:T]}\in\mathcal{K}_{[0:T]}$.

Proof  See Appendix E

Our next objective is to approximate the map
$x_{[0:T]}\mapsto(\Pi_{x_{[0:t]}})_{t=0}^{T}$ while respecting the forward
flow of information in time. We will accomplish this by developing a more
general approximation theory for causal maps between Riemannian manifolds
which are global NPC spaces.

##  5 Universal approximation of manifold-valued processes

In this section, we introduce our sequential geometric deep learning model,
called the HGN and defined below in (Static). We then use it to prove a
universal approximation theorem guaranteeing that the map
$x_{[0:t]}\to\Pi_{x_{[0:t]}}$ can be causally (in time) approximated to
arbitrary precision if the memory of the Volterra process decays at-least
polynomially.

######  Informal Results 2 (HGNs can approximate projected RCDs)

For every time-horizon $T$, and every suitable process $X_{\cdot}$ defined by
(Volterra), for each $\varepsilon>0$ and every compact set of inputs
$\mathcal{K}_{[0:T]}\subseteq\mathbb{R}^{(1+T)\,d}$, there exists a HGN
$F:(\mathbb{R}^{d})^{T}\rightarrow\mathcal{M}^{T}$, as in Figure 2, satisfying

$\underset{x_{[0:T]\in\mathcal{K}_{[0:T]}}}{\max}\underset{t=T-\operatorname{{H}},\dots,T}{\max}\,\,d_{g}\big{(}\Pi_{x_{[0:t]}},F(x_{\cdot})_{t}\big{)}<\varepsilon.$

The number of parameters and structure of the HGN depends quantitatively on
$\varepsilon,\operatorname{{H}},\mathcal{K}_{[0:T]}$, and on the
$\operatorname{Drift}$ and $\operatorname{Diffusion}$, and on the Volterra
kernel $\kappa$, where $\operatorname{{H}}\in\mathbb{N}$ is the history of the
HGN.

The main result of this section (Theorem 17) is a significant generalization
of Informal Result 2, as it applies to infinite-memory dynamical systems
between any pair of well-behaved non-positively curved complete Riemannian
manifolds. We now present our main approximation theoretic results.

###  5.1 Standing Assumptions for Universal Approximation Results

In this section, we present our model and build on the theory of causal maps,
as formulated in Acciaio et al. (2023), between a source and target connected
Riemannian manifolds $(\mathcal{N},h)$ and $(\mathcal{M},g)$.

######  Assumption 12 (Regularity of input and output manifolds)

Consider two $C^{\infty}$ Riemannian manifolds $(\mathcal{N},h)$ and
$(\mathcal{M},g)$ of dimensions $D\in\mathbb{N}$ and $d\in\mathbb{N}$
respectively. We assume that:

  1. (i)

$(\mathcal{N},h)$ and $(\mathcal{M},g)$ are geodesically complete and simply
connected.

  2. (ii)

Their sectional curvatures are bounded in $[-\kappa,0]$ for some $\kappa>0$.

We denote the geodesic distance on these spaces respectively by $d_{h}$ and
$d_{g}$. Under Assumption 12, $(\mathcal{N},d_{h})$ and $(\mathcal{M},d_{g})$
are complete separable metrics spaces of NPC in the sense of Sturm (2003); see
Proposition 2.

###### Example 2

The Euclidean space $(\mathbb{R}^{d},\boldsymbol{\delta})$ satisfies
Assumption 12 with $\kappa=0$.

###### Example 3

By Proposition 7, $(\mathcal{N}_{d},h)$ satisfies Assumption 12 with
$\kappa=-1/2$.

###  5.2 Static Case

We begin by defining the static version of our main (dynamical) model, in
which time is not yet considered. This geometric deep-learning model extends
the geometric deep network of Kratsios and Papon (2022), thus we use the same
name. From the technical standpoint, our main approximation theorem for the
static case (Theorem 13) is a significant technical improvement of the main
result of Kratsios (2023) as it has optimal rates, matching the optimal bounds
of Shen et al. (2022) in the classical Euclidean setting when the target
function is Lipschitz and/or smooth; this is an exponential improvement of the
guarantees in Kratsios (2023) for the $\operatorname{ReLU}$ case. Moreover, we
obtain probabilistic approximation guarantees for Borel measurable and for
general functions being approximated on finite subsets of the input space. The
static case will serve as the first main technical tool for our main
approximation theorem for the dynamic case, (Theorem 17),

Figure 1: The GDN Model: GDN model process an input in
$x_{[-M:0]}\in\mathcal{N}^{1+\operatorname{{H}}}$, interpreted as sequential
points $x_{-M},\dots,x_{0}$ inputs in $\mathcal{M}$, in three steps: an
encoding, transformation, and decoding phase. First, it linearized (purple)
the inputs in $\mathcal{N}^{1+\operatorname{{H}}}$ along products of geodesics
emanating from a set of reference points $x_{0}^{\star},\dots,x^{\star}_{M}$
in $\mathcal{N}^{1+\operatorname{{H}}}$. It then transforms the linearized
features and maps them to a vector $v$ in the tangent space of $\mathcal{M}$
using a standard ReLU-MLP (yellow). In the decoding phase (green), the model
maps $v$ to a point $\hat{f}(x_{[-M:0]})$ on $\mathcal{M}$ by travelling
geodesics in $\mathcal{M}$ emanating from a reference point $y^{\star}$
therein with initial velocity $v$.

We now formalize the static version of our model, illustrated in Figure 1, and
extend the geometric deep network (GDN) of Kratsios and Papon (2022). Fix
$J,\operatorname{{H}}\in\mathbb{N}$, with $J>0$, and a multi-index
$[\mathbf{d}]\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(d_{0},\ldots,d_{J})$ with $d_{0}=(1+\operatorname{{H}})D$ and
$d_{J}=d$, and let $P([\mathbf{d}])\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(1+\operatorname{{H}})D+d+\sum_{j=0}^{J}\,d_{j+1}(1+d_{j})$.
Weights, biases, and base-point parameters are encoded as a unique parameter
$\theta\in\mathbb{R}^{P([\mathbf{d}])}$ with

$\theta\,\,\leftrightarrow\,\,\Big{(}(A^{(j)},b^{(j)})_{j=0}^{J},(\tilde{x}_{m})_{m=0}^{\operatorname{{H}}},\tilde{y}\Big{)}$
(5.1)

where for $j=0,\dots,J$, $A^{(j)}\in\mathbb{R}^{d_{j+1}\times d_{j}}$,
$b^{(j)}\in\mathbb{R}^{d_{j+1}}$, $\tilde{y}\in\mathbb{R}^{d}$, and
$\tilde{x}_{m}\in\mathbb{R}^{D}$ for $m=0,\dots,\operatorname{{H}}$. We
henceforth denote componentwise composition of the
$\operatorname{ReLU}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\max\\{0,\cdot\\}$ activation function with any vector
$u\in\mathbb{R}^{n}$ for any $n\in\mathbb{N}_{+}$, by
$\operatorname{ReLU}\bullet x\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\big{(}\operatorname{ReLU}(x_{i})\big{)}_{i=1}^{n}$. Given
$\operatorname{{H}}\in\mathbb{N}$, a GDN model is a (parameterized) map
$f_{\theta}:\mathbb{R}^{P([\mathbf{d}])}\times\mathcal{N}^{1+\operatorname{{H}}}\rightarrow\mathcal{M}$
sending any
$(\theta,x_{[0:\operatorname{{H}}]})\in\mathbb{R}^{P([\mathbf{d}])}\times\mathcal{N}^{1+\operatorname{{H}}}$
to

$\displaystyle
f(\theta,x_{[0:\operatorname{{H}}]})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}f_{\theta}(x_{[0:\operatorname{{H}}]})$
$\displaystyle=\operatorname{Exp}^{g}_{\bar{y}}\circ\iota^{-1}_{\mathcal{M},\bar{y}}(A^{(J)}x^{(J)}+b^{(J)})$
(Static) $\displaystyle x^{(j+1)}$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{ReLU}\bullet(A^{(j)}x^{(j)}+b^{(j)})\quad\text{for\,\,}j=0,\ldots,J-1,$
$\displaystyle x^{(0)}\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}$
$\displaystyle\big{(}\operatorname{Log}_{\bar{x}_{m}}^{h}\circ\iota_{\mathcal{N},\bar{x}_{m}}(x_{m})\big{)}_{m=0}^{\operatorname{{H}}},$

where the encoding/feature and decoding/readout maps are defined by the
product of Riemannian exponential map at pre-specified points and its inverse
up to identification with the relevant Euclidean spaces via choices of linear
isomorphisms
$\iota_{\mathcal{N},\bar{x_{m}}}:T_{\bar{x}_{m}}\mathcal{M}\to\mathbb{R}^{d}$
and $\iota_{\mathcal{M},y}:T_{y}\mathcal{M}\to\mathbb{R}^{d}$.

We encode the trainable parameters
$\bar{x}_{[0:\operatorname{{H}}]}=(\bar{x}_{m})_{m=0}^{\operatorname{{H}}}\in\mathcal{N}^{1+\operatorname{{H}}}$
and $\bar{y}\in\mathcal{M}$ vectorially via

$\displaystyle\bar{x}_{[0:\operatorname{{H}}]}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(\operatorname{Exp}_{x^{\star}}^{h}(\tilde{x}_{m}))_{m=0}^{\operatorname{{H}}}\,\,\,\mbox{
and }\,\,\,\bar{y}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{Exp}_{y^{\star}}^{g}(\tilde{y}),$ (5.2)

where
$\tilde{x}_{[0,\operatorname{{H}}]}\in\mathbb{R}^{(1+\operatorname{{H}})D}\cong(T_{x^{\star}}\mathcal{N})^{1+\operatorname{{H}}}$,
$\tilde{y}\in\mathbb{R}^{d}\cong T_{y^{\star}}\mathcal{M}$ for arbitrary, but
a-priori fixed, points $x^{\star}\in\mathcal{N}$ and
$y^{\star}\in\mathcal{M}$. The set of all GDNs with representation
(Static)-(5.2) is denoted by
$\mathcal{GDN}_{[\mathbf{d}],\operatorname{{H}}}$.

The following universal approximation theorem significantly improves on
(Kratsios and Papon, 2022, Theorem 9) in three ways. First, when the target
function is smooth, we obtain dimension-free approximation rates, unlike
(Kratsios and Papon, 2022, Theorem 9), which only gives cursed approximation
rates regardless of the regularity of the target function (on infinite sets).
Second, it exponentially improves on the loose rate given therein when
approximating $\alpha$-Hölder functions; indeed, when $(\mathcal{N},h)$ is a
Euclidean space and $(\mathcal{M},g)$ is the real line with Euclidean metric,
our result matches the optimal approximation rates of (Shen et al., 2022,
Theorem 2.4). Third, when the compact set on which the function is being
approximated is finite, we obtain approximation guarantees without assuming
any continuity of the target function using the fact that such functions have
smooth extensions.

######  Theorem 13 (Approximation capacity of GDNs for smooth and Hölder
concept classes)

Fix a depth parameter $J\in\mathbb{N}_{+}$. Suppose that $(\mathcal{N},h)$ and
$(\mathcal{M},g)$ satisfy Assumptions 12,
$\operatorname{{H}}\in\mathbb{N}_{+}$, and $0<\varepsilon\leq 1$. Fix
reference points $\bar{x}\in\mathcal{N}^{1+\operatorname{{H}}}$ and
$\bar{y}\in\mathcal{M}$ as in (Static). For every locally $\alpha$-Hölder map
(resp. $k$-times continuously differentiable, resp. arbitrary if $\mathcal{K}$
is finite) $f:\mathcal{N}^{\operatorname{{H}}+1}\rightarrow\mathcal{M}$, for
some $0<\alpha\leq 1$ (resp. $k\in\mathbb{N}_{+})$, and each non-empty compact
$\mathcal{K}_{[0:\operatorname{{H}}]}\subseteq\mathcal{N}^{\operatorname{{H}}+1}$
there is a GDN
$\hat{f}:\mathcal{N}^{\operatorname{{H}}+1}\rightarrow\mathcal{M}$ satisfying

$\sup_{x_{[0:\operatorname{{H}}]}\in\mathcal{K}_{[0:\operatorname{{H}}]}}\,d_{g}\big{(}f(x_{[0:\operatorname{{H}}]}),\hat{f}(x_{[0:\operatorname{{H}}]})\big{)}<\varepsilon.$

and depth and width as recorded in Table 10.

Proof  See Appendix F.

By randomizing, we obtain a global probabilistic universal approximation
theorem for general Borel functions. This result can be interpreted as a non-
Euclidean generalization of (Hornik et al., 1989, Theorem 2.3) where the
authors quantify the accuracy with which MLPs can approximate Borel functions
between Euclidean spaces, which are not necessarily in some $L^{p}$-space,
using a weaker notion than uniform convergence on compact sets.

Table 2: Complexity estimates for the GDN in Theorem 13 $C^{k}$-Times Cnt. Diff.  ($k\in\mathbb{N}_{+}$) |
---|---
Depth | $\tilde{\mathcal{O}}\big{(}D(J+k^{2}+\operatorname{{H}}d)\big{)}$
Width Parameter | $\tilde{\mathcal{O}}\big{(}J(\sqrt{D}/\varepsilon)^{(1+\operatorname{{H}})/(2k)}\big{)}$
Locally $\alpha$-Hölder ($0<\alpha\leq 1$) |
Depth | $\mathcal{O}(JD)$
Width Parameter | $\mathcal{O}\big{(}\operatorname{{H}}d\,V(J^{-2}\varepsilon)^{\frac{1}{\alpha(1+\operatorname{{H}})d}}\big{)}$
No. Reg. - Finite $N\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\mathcal{K}_{[0:\operatorname{{H}}]}$ |
Depth | $\mathcal{O}(N^{2}\,J)$
Width Parameter | $\mathcal{O}(\varepsilon^{-\frac{1}{2N\,J}}+N)$
Explicit expressions for all constants are given in Table 10.

######  Corollary 14 (GDN are universal approximators of Borel functions)

Let $(\mathcal{N},h)$ and $(\mathcal{M},g)$ be Riemannian manifolds satisfying
Assumptions 12. Let $\operatorname{{H}}\in\mathbb{N}_{+}$ and
$X_{0},\dots,X_{\operatorname{{H}}}$ are i.i.d. $\mathcal{N}$-valued random
variables defined on a common probability space
$(\mathcal{S},\mathcal{A},\mu)$. For any Borel function
$f:\mathcal{N}^{1+\operatorname{{H}}}\rightarrow\mathcal{M}$ and every
approximation error $0<\varepsilon\leq 1$, there is a GDN
$\hat{f}:\mathcal{N}^{1+\operatorname{{H}}}\rightarrow\mathcal{M}$ satisfying
the probabilistic bound

$\mu\big{(}d_{h}\big{(}f(X_{[0:\operatorname{{H}}]}),\hat{f}(X_{[0:\operatorname{{H}}]})\big{)}<\varepsilon\big{)}>1-\varepsilon.$

Proof  See Appendix F.

Next, we consider the dynamic version of our results, and our main
approximation theorem.

###  5.3 Dynamic/Sequential universality

Before presenting the dynamic universal approximation theorems for GDNs, we
formalize their architecture summarized graphically in Figure 2. Briefly, the
blue loop is an auxiliary feed-forward neural network defined on the GDN
model’s parameter space. This auxiliary feedforward neural network, called
hypernetwork (e.g., Von Oswald et al. (2019)), synchronizes the parameters of
several GDNs each of which independently/in parallel approximates the map
$x_{[0:T]}\mapsto(\Pi_{x_{[0:t]}})_{t=0}^{T}$ for a unique $t$. The
hypernetwork synchronizes these GDNs, memorizing the parameters of an optimal
GDN at time $t+1$ from those defining an optimal one at time $t$, for each
$t=0,\dots,T$. Since memorization requires exponentially fewer parameters than
approximation, then the HGN can approximate the causal map
$x_{[0:T]}\mapsto(\Pi_{x_{[0:t]}})_{t=0}^{T}$ using roughly the same number of
parameters as one GDN needed to approximate the map at any one step in time.
This allows us to obtain our main “dynamic” approximation theorem.

Figure 2: The HGN Model: The green layer encodes sequence segments in the
input manifold into distances relative to a reference/landmark point
$x^{\star}$ therein. These linearized features are then processed through a
ReLU MLP, illustrated by the yellow repeated applying fully-connected affine
(also called linear) layers interspersed with ReLU activation functions
orange. Finally, the purple decodes the vector $v$ generated by the
downsampled CNN into a manifold-valued prediction, by travelling along a
geodesic emanating from a reference/landmark point $y^{\star}$ therein with
initial velocity $v$.
The HGN Model: Applies the GDN model while iteratively updating its internal
parameters, at each time step, using an (blue) auxiliary ReLU network, called
a hypernetwork.

######  Definition 15 ((HGN) Hypergeometric Network)

A map $F:\mathcal{N}^{\mathbb{Z}}\rightarrow\mathcal{M}^{\mathbb{Z}}$ is
called a hyper-geometric network (HGN) if there are
$\operatorname{{H}},d^{\star},d,D\in\mathbb{N}_{+}$, a multi-index
$[\mathbf{d}]\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(d_{0}=(1+\operatorname{{H}})D,\dots,d_{J}=d)$, a hypernetwork
given by a ReLU MLP
$h:\mathbb{R}^{d^{\star}}\rightarrow\mathbb{R}^{d^{\star}}$, an initial
parameter $z\in\mathbb{R}^{d^{\star}}$, and a linear readout map
$\tilde{L}:\mathbb{R}^{d^{\star}}\rightarrow\mathbb{R}^{P([\mathbf{d}])}$
representing $F$ as

$\displaystyle F\big{(}x_{\cdot}\big{)}_{t}$
$\displaystyle=f_{\tilde{L}(z_{t})}\big{(}x_{[t-\operatorname{{H}}:t]}\big{)}$
$\displaystyle z_{t}$ $\displaystyle=\begin{cases}h(z_{t-1})&\mbox{ if }t>0\\\
z&\mbox{ if }t\leq 0\end{cases},$

for every $t\in\mathbb{Z}$ and each sequence
$x_{\cdot}\in\mathcal{N}^{\mathbb{Z}}$, where
$f_{\theta}\in\mathcal{GDN}_{[\mathbf{d}],\operatorname{{H}}}$. The quadruple
$(z,\tilde{L},h,[\mathbf{d}])$ is called a HGN with history $H$ and latent
dimension $d^{\star}$. The causal map $F$ is called the representation of the
HGN $(z,\tilde{L},h,[\mathbf{d}])$.

When clear from the context, we refer to the representation of a HGN simply as
a HGN.

We now formalize what it means for a map
$f:\mathcal{N}^{\mathbb{N}}\rightarrow\mathcal{M}^{\mathbb{N}}$ to flow
information forward in time while not overemphasizing the arbitrarily distant
past. When $\mathcal{N}$ and $\mathcal{M}$ are the Euclidean space, this
definition is a quantitative version of the fading memory property for
reservoir systems (e.g., Lukoševičius and Jaeger (2009); Gonon and Ortega
(2021)) and it is a variant of the approximable complexity property of Acciaio
et al. (2023) and the related notion introduced in (Galimberti et al., 2022,
Definition 9). We use
$C^{\lambda}_{\alpha}(\mathcal{N}^{\operatorname{{H}}},\mathcal{M})$ to denote
the space of $\alpha$-Hölder functions from $\mathcal{N}^{\operatorname{{H}}}$
to $\mathcal{M}$ whose $\alpha$-Hölder coefficients is at most
$\lambda$-Lipschitz. We use
$C^{k,\lambda}(\mathcal{N}^{\operatorname{{H}}},\mathcal{M})$ to denote the
space of all $C^{k}$-functions
$f:\mathcal{N}^{\operatorname{{H}}}\to\mathcal{M}$ for which there exists
$x_{1},\dots,x_{\operatorname{{H}}}\in\mathcal{N}^{\operatorname{{H}}}$ and
some $y\in\mathcal{M}$ such that

$\operatorname{Log}_{y}^{g}\circ
f\circ\Big{(}\prod_{m=1}^{\operatorname{{H}}}\operatorname{Exp}_{x_{m}}^{h}\Big{)}$

has at most $\lambda$-Lipschitz $k^{th}$ partial derivatives.

######  Definition 16 (Causal maps of finite virtual memory between Hadamard
manifolds)

A map $f\,:\,\mathcal{K}\rightarrow\mathcal{M}^{\mathbb{Z}}$ from a compact
subset $\mathcal{K}\subseteq\mathcal{N}^{\mathbb{Z}}$ is called a causal map
with virtual memory $r\geq 0$, if for every “compression level”
$\varepsilon>0$ and each “time-horizon” $T\in\mathbb{N}_{+}$ there is a
positive integer $\operatorname{{H}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{{H}}(\varepsilon,T,\mathcal{K})\in
O(\varepsilon^{-r})$ with $\operatorname{{H}}\leq T$ and functions
$f_{1},\dots,f_{T}\in C(\mathcal{N}^{\operatorname{{H}}},\mathcal{M})$
satisfying

$\max_{t\in[[T]]}\sup_{x\in\mathcal{K}}\,d_{g}\big{(}f(x)_{t},f_{t}(x_{(t-\operatorname{{H}}:t]})\big{)}<\varepsilon.$
(5.3)

We say that $f$ is $(r,k,\lambda)$-smooth if there are $k\in\mathbb{N}_{+}$
and $\lambda>0$ such that each $t\in\mathbb{N}_{+}$ the map $f_{t}$ belongs to
$C^{k,\lambda}(\mathcal{N}^{\operatorname{{H}}},\mathcal{M})$.
We say that $f$ is $(r,\alpha,\lambda)$-Hölder $k\in\mathbb{N}_{+}$ and
$0<\alpha\leq 1$ such that for each $t\in\mathbb{N}_{+}$ the map $f_{t}$
belongs to
$C_{\alpha}^{\lambda}(\mathcal{N}^{\operatorname{{H}}},\mathcal{M})$.

In order for causal maps to be uniformly approximable on compact sets, they
need to possess a basic level of continuity. Therefore, we focus on the
following two classes of “high regularity” and “low regularity” causal maps of
finite virtual memory as approximable target functions, respectively. These
definitions were adapted from the infinite-dimensional linear analogues in
Galimberti et al. (2022) to the finite-dimensional non-vectorial setting
considered herein.

We may now state our main universal approximation theorem, which shows that
HGN can approximate causal maps of finite virtual memory. Quantitative rates
are obtained when the target function is either smooth or Hölder in the above
sense. Our result is quantitative and depends only on the regularity and
memory of the target causal map being approximated, as well on the structure
of the compact path space over which the uniform spatio-temporal approximation
is guaranteed.

Table 3: Complexity estimates for the HGN in Theorem 17 Hyperparame. | Estimate
---|---
Depth | $\tilde{\mathcal{O}}(T^{3/2})$
Width | $\mathcal{O}\big{(}P([\bar{J}])+Q)\,T\big{)}$
N. Parameters | $\tilde{\mathcal{O}}\Big{(}T^{4}(P([{J}])+Q)^{3}\Big{)},$
Latent Dimension | $\Theta(P([\mathbf{d}]))$
The depth, width, and value of $P([\mathbf{d}])$ is as in Table 10 and as
determined by the regularity of $f$ (and the cardinality of
$\mathcal{K}_{[0:T]}$) in Theorem 13, with the value of $\varepsilon$ instead
set to be $\varepsilon/2$.

######  Theorem 17 (Universal approximation of causal maps between Hadamard
manifolds)

Under Assumption 12, let $\mathcal{K}$ be a non-empty compact subset of
$\mathcal{N}^{\mathbb{Z}}$, and let
$f:\mathcal{K}\rightarrow\mathcal{M}^{\mathbb{Z}}$ be a causal map of finite
virtual memory. Fix depth and memory hyperparameters $J,Q\in\mathbb{N}_{+}$
and $0<\delta<1$.
Suppose either that 1) $f$ is $(r,k,\lambda)$-smooth ($r,\lambda\geq 0$ and
$k\in\mathbb{N}_{+}$), or 2) $f$ is $(r,\alpha,\lambda)$-Hölder
($r,\lambda\geq 0$ and $\alpha\in(0,1]$) or 3) that each $\mathcal{K}_{t}$ is
finite for $t\in\mathbb{Z}$ and $f$ has a finite virtual memory $r\geq 0$.
For every “approximation error” $\varepsilon>0$ and each “time-horizon”
$T_{\delta,Q}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\lfloor\delta^{-Q}\rfloor$, there is an integer
$0\leq\operatorname{{H}}\lesssim\min\\{T_{\delta,Q},\varepsilon^{-r}\\}$, and
an HGN $F:\mathcal{N}^{\mathbb{Z}}\rightarrow\mathcal{M}^{\mathbb{Z}}$ with
history $\operatorname{{H}}$ and latent dimension $d^{\star}=Q+D$ satisfying
the uniform approximation guarantee

$\max_{t=0,\dots,T_{\delta,Q}}\,\sup_{x\in\mathcal{K}}\,\,d_{g}\big{(}f(x)_{t},F(x)_{t}\big{)}<\varepsilon.$

Moreover, the structure and number of parameters determining $\hat{h}$ are
recorded in Tables 2 and 3.

Proof  See Appendix G.

##  6 Ablation Study

We now complement our theoretic analysis of the Gaussian random projection and
the HGN with a numerical analysis of these new tools. The primary purpose of
this section is to show how such a pipeline can be implemented and to explain
the role of each component of our model and how each of these interacts with
the Volterra process $X_{\cdot}$. Since there are no available benchmarks for
approximating dynamical systems on Riemannian manifolds, which are guaranteed
to be universal, we instead perform an ablation study to better understand the
dependence of each of these various factors determining the process
$X_{\cdot}$. Additional details are provided in Appendix I.

##### Experiment Setup

Consider a family of i.i.d. random Bernoulli variables
$(B_{t})_{t=0}^{\infty}$ taking values in $\\{0,1\\}$ with equal probabilities
of each state. Fix a “randomness” parameter $\lambda\geq 0$ and define the
random matrices

$\mathbf{S}_{t}\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\lambda
B_{t}\cdot I_{d}$

where $t\in\mathbb{N}$. Fix a weight $w\in(0,1]$, Lipschitz functions
$\mu:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$,
$\varsigma:\mathbb{R}^{d}\rightarrow(0,2]$, and a $d\times d$ symmetric
positive-definite matrix $\sigma$. Consider the $d$-dimensional stochastic
process

$\displaystyle X_{t+1}$
$\displaystyle=X_{t}+\operatorname{Drift}(X_{[t-1:t]})+\operatorname{Diffusion}(X_{t},\mathbf{S}_{t})\,W_{t}$
(6.1) $\displaystyle\operatorname{Drift}(z,x)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}w\mu(x)+(1-w)\,\mu(z)$
$\displaystyle\operatorname{Diffusion}(x,s)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\varsigma(x)\cdot\sigma+s,$

for $t\in\mathbb{N}$, where $(W_{t})_{t=0}^{T}$ are i.i.d. $d$-dimensional
standard normal random variables independent of $(B_{t})_{t=0}^{T}$, and both
$X_{-1},X_{0}$ are $d$-dimensional standard normal random variables. The
diffusion component of the process $X_{t}$, conditionally on $X_{t}$, randomly
moves between $\varsigma(X_{t})\cdot\sigma$ and
$\varsigma(X_{t})\cdot\sigma+\lambda I_{d}$ with equal probabilities,
independently of the driving Brownian motion. For any $T\in\mathbb{N}_{+}$,
path $x_{[-1:T]}\in\mathbb{R}^{(2+T)d}$, and integer $0\leq t\leq T$, the
$\mathcal{N}_{d}$-valued random variable $\mathbb{Q}_{x_{[-1:t]}}$ is
distributed according to

$\displaystyle\mathbb{P}\big{[}\mathbb{Q}_{x_{[-1:t]}}=\mathcal{N}_{d}(x_{t}+\operatorname{Drift}(x_{[t-1:t]}),\varsigma(x_{t})^{2}\cdot\sigma^{2})\big{]}=$
$\displaystyle\frac{1}{2}$ (6.2)
$\displaystyle\mathbb{P}\big{[}\mathbb{Q}_{x_{[-1:t]}}=\mathcal{N}_{d}(x_{t}+\operatorname{Drift}(x_{[t-1:t]}),(\varsigma(x_{t})\cdot\sigma+\lambda
I_{d})^{2})\big{]}=$ $\displaystyle\frac{1}{2}.$

By (Sturm, 2003, Proposition 5.5) and the product Riemannian structure on
$(\mathcal{N}_{d},\mathfrak{J})$ we have that the barycenter of
$\operatorname{Law}(\mathbb{Q}_{x_{[0:t]}})$ is the Cartesian product of the
barycenters of its components, up to identification of
$(\mathcal{N}_{d},\mathfrak{J})$ with
$(\mathbb{R}^{d}\times\operatorname{Sym}_{+}(d),\boldsymbol{\delta}\oplus g)$.
Using, see (Bini and Iannazzo, 2013, page 1701), for the expression of the
barycenter between two-points in $(\operatorname{Sym}_{+}(d),g)$ we find that

$\beta(\mathbb{Q}_{x_{[-1:t]}})=\mathcal{N}_{d}\big{(}x_{t}+\operatorname{Drift}(x_{[t-1:t]}),\,\,\varsigma(x_{t})\cdot\sigma^{2}(\sigma^{-2}(\lambda
I_{d}+\varsigma(x_{t})\cdot\sigma)^{2})^{1/2}\big{)}.$ (6.3)

Next, we confirm that the HGN model can indeed approximate the map
$x_{[-1:t]}\to\beta(\mathbb{Q}_{x_{[-1:t]}})$ in practice. Furthermore, we
inspect the dependence of each of the components of our framework on the
parameters defining the Volterra process (6.1); namely, the drift $\mu$, the
diffusion $\sigma$, $\varsigma$, the “randomness of the stochastic factor
process” $\lambda>0$, and the effect of non-Markovianity $w$.

##### HGN Training Pipeline

The HGN model is trained as follows. First, we sample several
$N\in\mathbb{N}_{+}$ paths segments $\\{x^{(n)}_{[-1:T]}\\}_{n=1}^{N}$ up to
time $T\in\mathbb{N}_{+}$ and we train a GDN to predict
$y^{(n)}_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\beta(\mathbb{Q}_{x^{(n)}_{[-1:1]}})$ given each sampled path, by
minimizing the intrinsic mean squared error (IMSE)

$\ell_{1}(\theta)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{n=1}^{N}\,d_{g}(f_{\theta}(x_{[-1:1]}^{(n)}),y^{(n)}_{1})^{2}$

where $\theta$ parameterized a set of GDNs of pre-specified depth and width.
The IMSE is optimized using the native ADAM optimizer built into Pytorch until
a suitable GDN parameter $\theta_{1}$ is obtained.

Then, for every subsequent time $t$, each we train a GDN by rolling the
training window forward and minimizing the corresponding GDN

$\ell_{t}(\theta)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{n=1}^{N}\,d_{g}(f_{\theta}(x_{[-t-2:t]}^{(n)}),y^{(n)}_{t})^{2}$
(6.4)

where $y^{(n)}_{t}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\beta(\mathbb{Q}_{x^{(n)}_{[-t-2:t]}})$. To avoid instability due
to the several symmetries present in the parameter space of most MLPs, see
e.g. Ainsworth et al. (2023); Sharma et al. (2024), and thus of our GDNs, we
initialize the optimization of each GDN at time $t+1$ using the optimized
parameters $\theta_{t}$ obtained by minimizing $\ell_{t}$ at time $t$.
Additionally, this implicitly encodes a transfer-learning effect, whereby the
GDN responsible for predicting at time $t$ encodes the pre-trained structure
in previous times. We note that when training the first GDN, $20$ ADAM epochs
are used while subsequent GDNs are sequentially fine-tuned using $10$ ADAM
epochs.

Once we have trained each “expert” GDN $\\{f_{\theta_{t}}\\}_{t=0}^{T}$,
specialized only on approximating $\beta(\mathbb{Q}_{\cdot})$ at each time
$t$, the HGN can be trained by simply minimizing the “hyper-MSE”
$\ell_{\text{hyper}}$ in the common parameter space $\mathbb{R}^{p}$ of the
GDNs. Namely,

$\ell_{\text{hyper}}(\vartheta)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{t=0}^{T-1}\,\|h_{\vartheta}(\theta_{t})-\theta_{t+1}\|^{2},$

where $\vartheta$ encodes the parameters of a hypernetwork of fixed depth and
width. The HGN is then fully encoded into the pair $(\theta_{0},h)$.
Additional details, and pseudo-code is contained in Appendix I and our code
can be found at our Github page https://github.com/arabporr/HyperNetwork .

###  6.1 Ablation Results

We study the sensitivity of the HGN and GDN models to the principle
characteristics dictating the stochastic evolution of $X$. We subsequently
study the effect of encoding a large number GDN “experts” into a single HGN.

We fixed one base problem (a process with a specific set of parameters) and
$30$ additional variations used during our ablation study, each with a similar
set of parameters but with exactly one hyperparameter different (e.g. drift,
volatility, etc…) perturbed during each ablation. The training set consists of
the first $t=0,\dots,159$ time steps and the test set consists of the final
$160,\dots,200$ time steps of the process $X_{\cdot}$. In each result we
report $95\%$ empirical confidence intervals. All experiment details on the
computational resources used are in Appendix I.1.

##### Sensitivities to aspects of $X_{\cdot}$

We begin by ablating the sensitivity of the HGN and GDN models to: the
simplicity/complexity of the drift $(\mu)$, the level of randomness
$(\lambda)$ in the stochastic factor $\mathbf{S}_{\cdot}$, the effect of
large/small fluctuations $(\varsigma)$ in the diffusion, the dimension $(d)$
of the process $X_{\cdot}$, and the level of non-Markovianity/persistence of
memory $(w)$ of the process $X$. In each case, we report the intrinsic mean
squared error for the GDN and HGN models and the confidence intervals formed
from one standard deviation of the loss distribution about the (mean)
intrinsic mean squared error across all time steps in the test set.

Table 4: Drift Ablation: Sensitivity to the structure of the drift $(\mu)$ of
$X$.

$\mu$ | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
$\frac{1}{100}$ | $1.27\times{10^{-6}}$ | $[1.18,1.36]\times{10^{-6}}$ | $4.78\times{10^{-4}}$ | $[3.89,5.67]\times{10^{-4}}$
$\frac{1}{10}$ | $2.21\times{10^{-6}}$ | $[2.03,2.39]\times{10^{-6}}$ | $4.39\times{10^{-2}}$ | $[3.65,5.13]\times{10^{-2}}$
$\frac{1}{2}\,(\frac{1}{100}-x)$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$
$e^{-x}+\cos(\frac{x}{100})$ | $1.04\times{10^{-5}}$ | $[0.88,1.19]\times{10^{-5}}$ | $1.41\times{10^{+1}}$ | $[1.29,1.53]\times{10^{+1}}$

Table 4 shows that the HGN and GDN models can predict Volterra processes whose
drift is both simple, e.g. constant, or complicated, e.g. exhibiting
osculations $\cos(x/100)$ and decay such as $e^{-x}$. Nevertheless, as one
would expect, the more complicated drifts are more difficult to learn for both
models; as reflected by larger test set errors. Moreover, as the drift becomes
more complicated the gap between the test set performance of the HGN and GDN
models grows as the parameters of the GDN become increasingly difficult to
predict for the hypernetwork in the HGN model.

Table 5: Random Factor Ablation: Sensitivity to the randomness $(\lambda)$ in
the stochastic factor process $\mathbf{S}$.

$\lambda$ | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
0 | $3.47\times{10^{-7}}$ | $[3.36,3.58]\times{10^{-7}}$ | $3.49\times{10^{-7}}$ | $[3.41,3.58]\times{10^{-7}}$
0.1 | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$
0.25 | $2.32\times{10^{-5}}$ | $[2.19,2.46]\times{10^{-5}}$ | $7.67\times{10^{-4}}$ | $[6.87,8.46]\times{10^{-4}}$
0.5 | $8.94\times{10^{-5}}$ | $[8.25,9.63]\times{10^{-5}}$ | $4.22\times{10^{-3}}$ | $[3.75,4.69]\times{10^{-3}}$
0.75 | $2.25\times{10^{-4}}$ | $[2.07,2.44]\times{10^{-4}}$ | $1.34\times{10^{-2}}$ | $[1.21,1.47]\times{10^{-2}}$
0.9 | $3.30\times{10^{-4}}$ | $[3.00,3.60]\times{10^{-4}}$ | $1.69\times{10^{-2}}$ | $[1.55,1.83]\times{10^{-2}}$
1 | $4.30\times{10^{-4}}$ | $[3.91,4.69]\times{10^{-4}}$ | $2.12\times{10^{-2}}$ | $[1.94,2.30]\times{10^{-2}}$

Table 5 shows that all models have increasingly larger challenges when
predicting from processes with large levels of randomness $(\lambda)$ in the
stochastic factor $\mathbf{S}_{\cdot}$ influencing their diffusion component.
This is because the larger $\lambda$ is, the more spread out both states (6.2)
of the random variable $\mathbb{Q}_{x_{[-1:t]}}$ becomes and, consequentially,
the more information is lost when computing the intrinsic averaging using
$\beta$. As anticipated, highly random stochastic factors produce a larger gap
in the test set performance of the GDN and the HGN models. In contrast, less
randomness in the stochastic factor yields a smaller test set performance gap
between the GDN and HGN models.

Table 6: Dimension Ablation: Sensitivity to Dimension $(d)$ of the Volterra
Process $X$.

$d$ | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
2 | $3.84\times{10^{-6}}$ | $[3.32,4.36]\times{10^{-6}}$ | $4.27\times{10^{-5}}$ | $[3.83,4.71]\times{10^{-5}}$
5 | $4.44\times{10^{-6}}$ | $[4.14,4.75]\times{10^{-6}}$ | $9.56\times{10^{-5}}$ | $[0.88,1.04]\times{10^{-4}}$
10 | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$
20 | $1.77\times{10^{-3}}$ | $[1.73,1.81]\times{10^{-3}}$ | $1.76\times{10^{-3}}$ | $[1.72,1.80]\times{10^{-3}}$

Table 6 confirms the effect of dimensionality on the expressive power of the
HGN and GDN models, described in Theorems 13 and 17 and detailed in Table 2.
Importantly, the performance of the HGN consistently mirrors that of the GDN
model in dimensions $2$ to $20$. Since roughly the same number of parameters
are used in each case, then, naturally, the performance of both models is
higher in low dimensions than in higher dimensions.

Table 7: Non-Markovianity Ablation: Sensitivity to the persistence of memory
$(w)$ in $X$.

Memory | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
0 | $3.34\times{10^{-8}}$ | $[2.82,3.85]\times{10^{-8}}$ | $7.97\times{10^{-6}}$ | $[6.63,9.31]\times{10^{-6}}$
0.1 | $9.31\times{10^{-5}}$ | $[8.64,9.98]\times{10^{-5}}$ | $9.32\times{10^{-5}}$ | $[0.86,1.00]\times{10^{-4}}$
0.25 | $1.02\times{10^{-4}}$ | $[0.94,1.10]\times{10^{-4}}$ | $1.02\times{10^{-4}}$ | $[0.94,1.10]\times{10^{-4}}$
0.5 | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$

Both the HGN and GDN models perform nearly identically for all degrees of
memory persistence, from Markovianity to higher levels of non-Markovian
memory. This confirms that the hypernetwork can reliably predict GDN
parameters regardless of the degree of memory, as in Theorem 17.

Table 8: Diffusion Ablation: Sensitivity to the size of the fluctuations
$(\varsigma)$ in the diffusion of $X$.

$\varsigma$ | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
0.005 | $4.53\times{10^{-6}}$ | $[4.19,4.88]\times{10^{-6}}$ | $1.39\times{10^{-4}}$ | $[1.25,1.54]\times{10^{-4}}$
0.01 | $5.02\times{10^{-6}}$ | $[4.67,5.37]\times{10^{-6}}$ | $1.46\times{10^{-4}}$ | $[1.27,1.66]\times{10^{-4}}$
0.05 | $1.01\times{10^{-5}}$ | $[0.94,1.08]\times{10^{-5}}$ | $4.04\times{10^{-4}}$ | $[3.55,4.54]\times{10^{-4}}$
0.1 | $1.90\times{10^{-5}}$ | $[1.79,2.02]\times{10^{-5}}$ | $7.30\times{10^{-4}}$ | $[6.43,8.17]\times{10^{-4}}$
1 | $1.17\times{10^{-3}}$ | $[1.09,1.24]\times{10^{-3}}$ | $3.33\times{10^{-2}}$ | $[2.97,3.70]\times{10^{-2}}$
10 | $4.13\times{10^{-1}}$ | $[3.69,4.57]\times{10^{-1}}$ | $4.56\times{10^{+0}}$ | $[4.14,4.99]\times{10^{+0}}$
100 | $2.80\times{10^{+3}}$ | $[2.74,2.86]\times{10^{+3}}$ | $3.27\times{10^{+3}}$ | $[3.18,3.35]\times{10^{+3}}$
1000 | $3.10\times{10^{+5}}$ | $[3.01,3.18]\times{10^{+5}}$ | $3.15\times{10^{+5}}$ | $[3.07,3.23]\times{10^{+5}}$

Table 8 shows that both models can reliably predict regardless of the
diffusion component of the process $X_{\cdot}$ has large or small
fluctuations. As expected, the reliability of the HGN predictions deteriorates
when $\varsigma$ increases, as can be seen by an increase in the standard
deviation of the loss.

Table 9: Curvature Ablation: Sensitivity to the size of the fluctuations
$(\varsigma)$ in the diffusion of $X$.

$\varsigma$ | GDN Loss Mean | GDN Loss 95% C.I. | HGN Loss Mean | HGN Loss 95% C.I.
---|---|---|---|---
0.000001 | $2.34\times{10^{-3}}$ | $[2.30,2.38]\times{10^{-3}}$ | $2.32\times{10^{-3}}$ | $[2.28,2.36]\times{10^{-3}}$
0.0001 | $1.52\times{10^{-3}}$ | $[1.48,1.57]\times{10^{-3}}$ | $1.54\times{10^{-3}}$ | $[1.50,1.58]\times{10^{-3}}$
0.001 | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$ | $1.58\times{10^{-3}}$ | $[1.55,1.60]\times{10^{-3}}$

One may a priori expect that a very small variance hyperparameter $\varsigma$
(near $0$) would further improve the performance of the GDN and HGN models. As
illustrated by Table 9, this is not the case since Gaussian measures with
nearly singular covariance matrices (which happens when $\varsigma\approx 0$)
are near the (Gromov) boundary999This concept generalizes the circle/boundary
of the hyperbolic disc model used in the visualization in Figure 2.
Importantly, all points in this so-called Gromov boundary are infinitely far
from points genuinely inside the non-positively curved Riemannian manifold.
See (Bridson and Haefliger, 1999, Chapter II.8) for details. “at infinity” of
our non-positively curved Riemannian manifold
$(\mathcal{N}_{d},\mathfrak{J})$. The trouble here is that the Lipschitz
constant of the exponential map/layer based at a Gaussian measure far from the
boundary is very large (see Lemma C.1) which translates to constants in our
approximation guarantees in Theorem 13 (see Table 10). Consequentially,
significantly more parameters are required for the HGN to achieve a comparable
approximation accuracy when the target Gaussian measure is near the Gaussian
measure used as the base point of the exponential layer of the GDN (as in
Table 8). This shows that the constants in our main result concretely impact
practical implementations of the GDN and the HGN models.

The next set of ablation studies will continue to examine the efficacy of the
hypernetwork in encoding and predicting GDN parameters in the test set.

####  6.1.1 Ablation of the Hypernetwork Encoding

Theorem 17 guarantees that the hypernetwork can effectively encode a large
number of GDNs. In particular, doing so suggests that the HGN model can be
recursively rolled, allowing us to predict well into the future. Two questions
naturally arise: 1) In practice, is a hypernetwork encoding of a sequence of
GDN models legitimately trainable? 2) Does the hypernetwork continue to
generate well-performing GDN models out-of-sample in future times? This
section yields an affirmative yes to both of these questions; thus showing the
feasibility and reliability of the hypernetwork in the HGN model.

In this next experiment, we test this by comparing three different degrees of
hypernetwork encoding. These experiments are run with a subset of the same
configurations from the last section; which we annotate in each figure
caption. Each of the figures 3 and 4 plot the test set performance of each
model, which begins at time $t=160$ and ends at time $t=200$.

1) (GDN) no hypernetwork is used and only a different GDN “expert” is used to
generate predictions at any time, trained at time $t-1$. 2) (HGN $1$ step) at
every time $t$, the hypernetwork loads the predictions of the GDN at time
$t-1$ and uses them to predict a GDN, which is then used to predict at time
$t$. In this case, the hypernetwork component of the HGN is only ever used to
make one-step-ahead predictions. 3) (HGN) loads the GDN parameters at time
$160$ and then sequentially predicts the parameters at every subsequent time
$t$ using its predicted parameters at time $t-1$; up until the terminal time
$t=200$. All three models perform nearly identically, as illustrated by the
log-scale losses. This shows that the hypernetwork encoding of the GDN models
is as practically effective as in Theorem 17.

Figure 3: Situation I - Nearly Logarithmic Degradation of HGN Accuracy: The
HGN performance slowly (logarithmically) departs from that of the GDN as time
rolls forward. This is typically what is observed in most of our experiments.

Our experiments uncover two types of behaviours which the HGN can exhibit. In
the former case, seen in Figure 3, there is a small but growing gap between
the test set performance of the HGN and the GDN model which increases as time
flows forward. Furthermore, this gap is roughly the same for both the $1$-step
and recurrent HGN models.

Figure 4: Situation II - Nearly Perfect GDN Prediction by HGN: The HGN
continues to nearly perfectly predict the performance of the GDN as time rolls
forward. This occurs in a subset of experiments where the GDN parameters do
not change significantly between time steps.

Occasionally, the GDN training suffers from exploding gradients during
training and one can re-run the stochastic gradient descent algorithm when
this happens. We note that there is nothing apriori particular about instances
when this happens (see Appendix I.4).

##  7 Conclusion

We presented a framework for obtaining low-dimensional approximations of the
conditional distribution (non-Markovian) stochastic Volterra processes in
discrete time. First, we develop a tool, the Gaussian projection, for
projecting the condition of such processes down onto the smooth manifold
$\mathcal{N}_{d}$ of non-singular Gaussian measure with a perturbation
$\mathfrak{J}$ of its standard information geometry with favourable geometric
and computational properties (see Table 1). Like classical tools for dimension
reduction of probability measures (e.g. information projections) the Gaussian
projection is a projection-type optimization problem; however, unlike those
tools the Gaussian projection is a Lipschitz operation (Theorem 9 and
Proposition 10) and can even be a smooth map under additional conditions
(Proposition 10).

Using these insights, we then constructed a sequential geometric deep learning
model which is compatible with the non-positive curvature of
$(\mathcal{N}_{d},\mathfrak{J})$. We showed that this model is a universal
approximator of generalized dynamicals systems (Theorem 17), possibly with
long memory, and we obtained quantitative approximation rates therefor
(Theorem 13). In the static case our model reduces the GDN model of Kratsios
and Papon (2022); in this case, our main universal approximation theorem
yields improved rates (compared to (Kratsios and Papon, 2022, Theorem 9)) when
approximating both Hölder, as well as better rates when approximating
functions on finite set (compare Table 2 against (Kratsios and Papon, 2022,
Theorem 20)). We additionally obtain new probabilistic approximation
guarantees for any Borel measurable function in that setting (Corollary 14).

We numerically illustrated the HGN model showing its practical viability. We
conducted an ablation study, confirming our main theoretical results. Further,
our ablation study shows that certain geometric constants appearing in our
parametric complexity estimates (Table 2), depending on the sectional
curvature of the underlying space, are legitimately non-negligible in
practice.

### Future Research

In future work, we would like to extend our analysis to general stochastic
processes, going beyond the stochastic Volterra setting. We would also like to
explore the impact of projecting onto different information-like geometries
when approximation the conditional law of various processes and how to choose
such geometries if one has information on the structure of these processes.

##  8 Acknowledgment and Funding

AK acknowledges financial support from an NSERC Discovery Grant No.
RGPIN-2023-04482. AK also acknowledges that resources used in preparing this
research were provided, in part, by the Province of Ontario, the Government of
Canada through CIFAR, and companies sponsoring the Vector Institute
https://vectorinstitute.ai/partnerships/current-partners/.

The authors would like to extend a very special thanks to Deane Yang from
Courant Institute for his helpful references on Riemannian Centers of Mass.
They would also like to thank McKenzie Yuen-Kong Wang from McMaster University
for his references on the regularity of the exponential map. Both sets of
references significantly helped our derivations.

##  Appendix A Explicit Tables

This section records the versions of Tables 10, 11, and 12 with explicit
constants.

###  A.1 Static Case

Table 10: Detailed Complexity estimates for the GDN in Theorem 13.

$k$-Times Continuously Differentiable $f$ | $k\in\mathbb{N}_{+}$
---|---
Depth ($J)$ | $D\Big{(}1+18k^{2}(J+2)\log_{2}(4J)+2(1+\operatorname{{H}})d\Big{)}$
Width Parameter ($w$) | $C_{0}+C_{1}\Big{(}2+\Big{\lceil}J\Big{(}\varepsilon^{-1}C_{\bar{f}}D^{1/2}\Big{)}^{(1+\operatorname{{H}})/(2k)}\Big{\rceil}\Big{)}\log_{2}\Big{(}8\Big{\lceil}J\Big{(}\varepsilon^{-1}C_{\bar{f}}D^{1/2}\Big{)}^{(1+\operatorname{{H}})d/(2k)}\Big{\rceil}\Big{)}$
Locally $\alpha$-Hölder | $0<\alpha\leq 1$
Depth ($J)$ | $D\Big{(}11J+C_{2}\Big{)}$
Width Parameter ($w$) | $C_{0}+C_{1}\,\max\big{\\{}(1+\operatorname{{H}})d\big{\lfloor}\big{(}\big{\lceil}V\big{(}J^{-2}\,(C_{2}\varepsilon/C_{f})^{d/\alpha}\big{)}\big{\rceil}\big{)}^{1/((1+\operatorname{{H}})\,d)}\big{\rfloor},2+\big{\lceil}V\big{(}J^{-2}\,(C_{2}\varepsilon/C_{f})^{1/(\alpha(1+\operatorname{{H}})d)}\big{)}\big{\rceil}\big{\\}}$
No. Reg. - Finite $\mathcal{K}_{[0:\operatorname{{H}}]}$ | $N\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\\#\mathcal{K}_{[0:\operatorname{{H}}]}$
Depth ($J)$ | $D(28N^{2}J)$
Width Parameter ($w$) | $\lceil(\varepsilon^{-1}\,C)^{1/(2NJ)}\rceil-1\big{)}+2N-8$

##### Details about Constants in Table 10:

(1) When $f$ is of class $C^{k}$, the constants are
$C_{0}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}d(1+\operatorname{{H}})(D-1)$,
$C_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}17k^{d(1+\operatorname{{H}})+1}3^{d(1+\operatorname{{H}})}$ and

$C_{\bar{f}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\max_{i=1,\dots,D}\,\|(\rho^{-1}\circ f\circ\varphi^{-1}\circ
W)_{i}\|_{C^{k}([0,1]^{(1+\operatorname{{H}})\,d})}/\max\\{1,\operatorname{Lip}\big{(}\rho|f\circ\varphi(\overline{\operatorname{Ball}(\mathcal{K}_{[0:\operatorname{{H}}]},1})\big{)}\\}$

where $W$ is any bijective affine self-map on
$\mathbb{R}^{(1+\operatorname{{H}})\,d}$ mapping
$W(\overline{\operatorname{Ball}(\varphi(\mathcal{K}_{[0:\operatorname{{H}}]}),1)})$
into $[0,1]^{(1+\operatorname{{H}})\,d}$.
(2) When $f$ is $\alpha$-Hölder, the constants are
$C_{0}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(1+\operatorname{{H}})d(D-1)$,
$C_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}3^{(1+\operatorname{{H}})d+3}$,
$C_{2}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}19+2(1+\operatorname{{H}})d$,
$C_{3}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(131((1+\operatorname{{H}})dD)^{1/2})$,

$C_{f}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{Lip}_{\alpha}(\rho^{-1}\circ
f\circ\varphi^{-1}|\varphi(\overline{\operatorname{Ball}(\mathcal{K}_{[0:\operatorname{{H}}]},1)}),$

and $V:[0,\infty)\rightarrow[0,\infty)$ is the “special function” defined as
the inverse of $t\mapsto t^{2}\log_{3}(t+2)$.
(3) When $\mathcal{K}$ is finite,
$C=(1+\operatorname{{H}})d(D-1)+9\big{(}\operatorname{diam}(f(\mathcal{K}_{[0:\operatorname{{H}}]}))\,\Big{(}\frac{(1+\operatorname{{H}})d}{2(1+\operatorname{{H}})d+2}\Big{)}^{1/2}\,DN$.

###  A.2 Dynamic Case

Next, we provide detailed expressions for the constants in the dynamic case.

Table 11: Complexity estimates for the HGN in Theorem 17

$k$-Times Continuously Differentiable $f$ | $k\in\mathbb{N}_{+}$
---|---
Depth ($J)$ | $D\Big{(}1+18k^{2}(J+2)\log_{2}(4J)+2(1+\operatorname{{H}})d\Big{)}$
Width Parameter ($w$) | $C_{0}+C_{1}\Big{(}2+\Big{\lceil}J\Big{(}\varepsilon^{-1}C_{\bar{f}}D^{1/2}\Big{)}^{(1+\operatorname{{H}})/(2k)}\Big{\rceil}\Big{)}\log_{2}\Big{(}8\Big{\lceil}J\Big{(}\varepsilon^{-1}C_{\bar{f}}D^{1/2}\Big{)}^{(1+\operatorname{{H}})d/(2k)}\Big{\rceil}\Big{)}$
Locally $\alpha$-Hölder | $0<\alpha\leq 1$
Depth ($J)$ | $D\Big{(}11J+C_{2}\Big{)}$
Width Parameter ($w$) | $C_{0}+C_{1}\,\max\big{\\{}(1+\operatorname{{H}})d\big{\lfloor}\big{(}\big{\lceil}V\big{(}J^{-2}\,(C_{2}\varepsilon/C_{f})^{d/\alpha}\big{)}\big{\rceil}\big{)}^{1/((1+\operatorname{{H}})\,d)}\big{\rfloor},2+\big{\lceil}V\big{(}J^{-2}\,(C_{2}\varepsilon/C_{f})^{1/(\alpha(1+\operatorname{{H}})d)}\big{)}\big{\rceil}\big{\\}}$
No. Reg. - Finite $\mathcal{K}_{[0:\operatorname{{H}}]}$ | $N\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\\#\mathcal{K}_{[0:\operatorname{{H}}]}$
Depth ($J)$ | $D(28N^{2}J)$
Width Parameter ($w$) | $\lceil(\varepsilon^{-1}\,C)^{1/(2NJ)}\rceil-1\big{)}+2N-8$

##### Details about Constants in Table 17

(1) When $f$ is $(r,k,\lambda)$-smooth: the constants are
$C_{0}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}d(1+\operatorname{{H}})(D-1)$,
$C_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}17k^{d(1+\operatorname{{H}})+1}3^{d(1+\operatorname{{H}})}$ and
$C_{\bar{f}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\displaystyle\max_{\underset{i=1,\dots,D}{t=T-\operatorname{{H}},\dots,T}}\,\frac{\|(\rho^{-1}\circ
f_{t}\circ\varphi^{-1}\circ
W_{t})_{i}\|_{C^{k}([0,1]^{d})}}{\max\\{1,\operatorname{Lip}\big{(}\rho|f_{t}\circ\varphi(\mathcal{K}^{\star:T})\big{)}\\}}$;
where for each $t=\operatorname{{H}},\dots,T$,
$W_{t}:\mathbb{R}^{(1+\operatorname{{H}})d}\ni\frac{1}{C_{\bar{x}}\sqrt{M}\,2r}\,\sqrt{\frac{2d(1+\operatorname{{H}})+2}{d(1+\operatorname{{H}})}}\,(x-\tilde{x}_{t})\in\mathbb{R}^{(1+\operatorname{{H}})d}$
is affine.
(2) When $f$ is $(r,\alpha,\lambda)$-Hölder: the constants are
$C_{0}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(1+\operatorname{{H}})d(D-1)$,
$C_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}3^{(1+\operatorname{{H}})d+3}$,
$C_{2}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}19+2(1+\operatorname{{H}})d$,
$C_{3}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(131((1+\operatorname{{H}})dD)^{1/2})$,
$C_{f}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{Lip}_{\alpha}(\rho^{-1}\circ
f\circ\varphi^{-1}|\varphi(\overline{\operatorname{Ball}(\mathcal{K}_{[0:\operatorname{{H}}]},1)})$,
and $V:[0,\infty)\rightarrow[0,\infty)$ is the “special function” defined as
the inverse of $t\mapsto t^{2}\log_{3}(t+2)$.
(3) When $\mathcal{K}_{[0:T]}$ is finite: $C\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(1+\operatorname{{H}})d(D-1)+9\big{(}\operatorname{diam}(f(\mathcal{K}_{[0:\operatorname{{H}}]}))\,\Big{(}\frac{(1+\operatorname{{H}})d}{2(1+\operatorname{{H}})d+2}\Big{)}^{1/2}\,DN$
and $N\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\binom{\max_{t=0,\dots,T}\,\mathcal{K}_{t}+M-1}{M}$.

We complete this section by providing more resolution on the rates in Table 3,
with exact expressions wherever possible.

Table 12: Complexity estimates for the HGN in Theorem 17

Hyperparame. | Estimate
---|---
Depth | $\mathcal{O}\Biggr{(}T\biggl{(}1+\sqrt{T\log(T)}\,\Big{(}1+\frac{\log(2)}{\log(T)}\,\biggl{[}C+\frac{\Big{(}\log\big{(}T^{2}\,2^{1/2}\big{)}-\log(\delta)\Big{)}}{\log(2)}\biggr{]}_{+}\Big{)}\biggr{)}\Biggr{)}$
Width | $(P([\bar{J}])+Q)\,T+12$
N. Parameters | $\mathcal{O}\Bigl{(}T^{3}(P([\bar{J}])+Q)^{2}\,\left(1+(P([\bar{J}])+Q)\sqrt{T\log(T)}\,\left(1+\frac{\log(2)}{\log(T)}\,\left[C_{d}+\frac{\Big{(}\log\big{(}T^{2}\,2^{1/2}\big{)}-\log(\delta)\Big{)}}{\log(2)}\right]_{+}\right)\right)\Bigr{)},$
Latent Dimension |  $d^{\star}\geq\max\\{12,P([\mathbf{d}])+1\\}$ (free parameter)

The values of $(J)$ and $(w)$ are as in Table 10 and $P([\bar{J}])$ is
determined accordingly by Theorem 13, with the value of $\varepsilon$ set to
be $\varepsilon/2$. These depend only on the regularity of the map $f$ and the
cardinality of $\mathcal{K}_{[0:T]}$.

##  Appendix B A problem with information projections for regime-switching
processes

Consider the special case of (Volterra) where $X_{\cdot}$ evolves according to
the regime-switching type hidden Markov model, with dynamics given by

$X_{t}=\sum_{s=0}^{t}\,S_{s}\,W_{s},$

where $\\{\mathbf{S}_{s},W_{s}\\}_{s=0}^{T}$ are independent random variables
where the “regimes” are encoded by the identically distributed Bernoulli
$\\{0,1\\}$ random variables $S_{s}$ with equal probabilities and the
$(W_{s})_{s=0}^{T}$ are standard Gaussian random variables. Let
$x\in\mathbb{R}$ be a realized state of the process at time 0; the process
$X_{.}$ is defined on the measurable space
$(\mathbb{R}^{\mathbb{N}},\mathcal{B}(\mathbb{R}^{\mathbb{N}}))$ where
$\mathcal{B}(\mathbb{R}^{\mathbb{N}})$ is the Borel $\sigma$-algebra on
countable product $\mathbb{R}^{\mathbb{N}}$. These processes are such that the
corresponding RCD being projected does not admit a Lebesgue density. More
precisely, $\mathbb{P}[X_{1}\in\cdot|X_{0}=x]$ does not admit a well-defined
$I$-projection on the set of $d$-dimensional Gaussian probability measures
$\mathcal{N}_{d}$ with non-singular covariance, since for an $I$-projection to
exist and be unique, the measure $\mathbb{P}[X_{1}\in\cdot|X_{0}=x]$ must
admit a (Lebesgue) density.

##  Appendix C Auxiliary differential geometric results

This appendix contains the required differential geometric results on which
many of our proofs are built.

###  C.1 Results from Riemannian geometry

The following result is due to (Rauch, 1951), but we make use of the
formulation which is essentially found in (Karcher, 1977a, C.1 and the remark
following (A4.2)). We record the additional minor details.

In what follows, for any $d$-dimensional Riemannian manifold $(\mathcal{M},g)$
and any point $p\in\mathcal{M}$, we fix a linear isomorphism
$\iota_{p}:\mathbb{R}^{d}\rightarrow T_{p}\mathcal{M}$ which formalizes the
identification of the tangent space of $\mathcal{M}$ at $p$ with the
$d$-dimensional Euclidean space.

######  Lemma C.1 ((Rauch, 1951; Karcher, 1977a))

Let $(\mathcal{M},g)$ be a $d$-dimensional Cartan-Hadamard manifold with
sectional curvatures bounded between $[-\kappa,0]$ for some constant
$\kappa>0$. For every $\bar{r}>0$ and every $p\in\mathcal{M}$, we have that

$\operatorname{Lip}\left(\operatorname{Exp}_{p}\circ\iota_{p}|\overline{\operatorname{Ball}_{(\mathbb{R}^{d},\|\cdot\|)}(0,\bar{r}/c_{p})}\right)\leq
2c_{p}\,\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}},$

where
$\operatorname{Lip}\left(\operatorname{Exp}_{p}{\circ\iota_{p}}|\overline{\operatorname{Ball}_{(\mathbb{R}^{d},\|\cdot\|)}(0,\bar{r}/c_{p})}\right)$
denotes the best local Lipschitz constant of
$\operatorname{Exp}_{p}\circ\iota_{p}$ on the the closed ball about the origin
of radius $\bar{r}/c_{p}$ and $c_{p}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left|\iota_{p}\right|_{op}$ denotes the operator norm of
$\iota_{p}$.

Proof of Lemma C.1 First of all, let us recall the following general fact:
given $(\mathfrak{X},g)$ and $(\mathfrak{Y},h)$ Riemannian manifolds, with
$\mathfrak{X}$ connected, and $F:\mathfrak{X}\to\mathfrak{Y}$ a $C^{1}$ map,
we have that the following are equivalent:

  1. 1.

$d_{\mathfrak{Y}}(F(x_{1}),F(x_{2}))\leq Cd_{\mathfrak{X}}(x_{1},x_{2})$ with
$x_{1},x_{2}\in\mathfrak{X}$;

  2. 2.

$\left|dF_{x}\right|_{op}\leq C$ with $x\in\mathfrak{X}$;

where $\left|dF_{x}\right|_{op}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sup_{\left|V\right|_{x}=1}\left|dF_{x}(V)\right|_{F(x)}$ and
where $\left|\cdot\right|_{x}$ simply denotes $\sqrt{g(\cdot,\cdot)_{x}}$, and
similarly $\left|\cdot\right|_{F(x)}=\sqrt{h(\cdot,\cdot)_{F(x)}}$.

Fix $p\in\mathcal{M}$, and for any $w\in T_{p}\mathcal{M}$, identify
$T_{w}(T_{p}\mathcal{M})$ with $T_{p}\mathcal{M}$ itself.

Under our hypothesis, we can apply (Karcher, 1977a, C.1): therefore, for any
$A,v\in T_{p}\mathcal{M}$ with $\left|A\right|_{p}=\left|v\right|_{p}=1$ and
$A\bot v$ (i.e. orthgonal), and $r>0$ it holds

$\left|d(\operatorname{Exp}_{p})_{rv}(A)\right|_{\operatorname{Exp}_{p}(rv)}\leq\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r}$
(C.1)

where again $\left|\cdot\right|_{p}$ simply denotes
$\sqrt{g(\cdot,\cdot)_{p}}$, and similarly for
$\left|\cdot\right|_{\operatorname{Exp}_{p}(rv)}$.

Let us now show an analogous bound for the general case in which not
necessarily it holds $A\bot v$. In an obvious notation, write
$A=A^{\parallel}+A^{\bot}$, where $A^{\bot}$ is orthogonal to $v$ and
$A^{\parallel}$ is parallel to $v$. Then, by linearity of the differential, we
have, for arbitrary $\left|A\right|_{p}=\left|v\right|_{p}=1$:

$\begin{split}\left|d(\operatorname{Exp}_{p})_{rv}(A)\right|_{\operatorname{Exp}_{p}(rv)}&\leq\left|d(\operatorname{Exp}_{p})_{rv}(A^{\parallel})\right|_{\operatorname{Exp}_{p}(rv)}+\left|d(\operatorname{Exp}_{p})_{rv}(A^{\bot})\right|_{\operatorname{Exp}_{p}(rv)}\\\
&\leq\left|A^{\parallel}\right|_{p}\left|d(\operatorname{Exp}_{p})_{rv}(v)\right|_{\operatorname{Exp}_{p}(rv)}+\left|A^{\bot}\right|_{p}\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r},\end{split}$

where in the second line we have used (C.1). By (Karcher, 1977a, C.1) again,
we have that
$\left|d(\operatorname{Exp}_{p})_{rv}(v)\right|_{\operatorname{Exp}_{p}(rv)}=\left|v\right|_{p}=1$.
Besides, $1\leq\frac{\operatorname{sinh}(\xi)}{\xi}$ for any $\xi>0$. Hence,
we may write,

$\left|d(\operatorname{Exp}_{p})_{rv}(A)\right|_{\operatorname{Exp}_{p}(rv)}\leq\left[\left|A^{\parallel}\right|_{p}+\left|A^{\bot}\right|_{p}\right]\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r}.$

In conclusion, we have proved that for any $p\in\mathcal{M},r>0$ and
$\left|A\right|_{p}=\left|v\right|_{p}=1$, it holds

$\left|d(\operatorname{Exp}_{p})_{rv}(A)\right|_{\operatorname{Exp}_{p}(rv)}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r},$

i.e.

$\left|d(\operatorname{Exp}_{p})_{rv}\right|_{op}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}r)}{\sqrt{\kappa}r}.$

Since $(0,\infty)\ni\xi\mapsto\frac{\operatorname{sinh}(\xi)}{\xi}$ is
(strictly) increasing, given $\bar{r}>0$, we have
$\left|d(\operatorname{Exp}_{p})_{rv}\right|_{op}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}}$ for
$r\leq\bar{r}$, i.e.

$\left|d(\operatorname{Exp}_{p})_{w}\right|_{op}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}},\quad
p\in\mathcal{M},\,w\in T_{p}\mathcal{M},\,\left|w\right|_{p}\leq\bar{r}.$

Therefore, we conclude that

$d_{M}(\operatorname{Exp}_{p}(w_{1}),\operatorname{Exp}_{p}(w_{2}))\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}}\left|w_{1}-w_{2}\right|_{p},\quad
p\in\mathcal{M},\,w_{1},w_{2}\in
T_{p}\mathcal{M},\,\left|w_{1}\right|_{p},\left|w_{2}\right|_{p}\leq\bar{r}.$
(C.2)

Let $c_{p}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left|\iota_{p}\right|_{op}$ and suppose now that
$v_{1},v_{2}\in\overline{\operatorname{Ball}_{\mathbb{R}^{d},\|\cdot\|}(0,\bar{r}/c_{p})}$.
Then,
$\iota_{p}(v_{1}),\iota_{p}(v_{2})\in\overline{\operatorname{Ball}_{T_{p}\mathcal{M},g}(0,\bar{r})}$.
Whence, (C.2) implies that

$d_{M}\big{(}\operatorname{Exp}_{p}(\iota_{p}(v_{1})),\operatorname{Exp}_{p}(\iota_{p}(v_{2}))\big{)}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}}\left|\iota_{p}(v_{1})-\iota_{p}(v_{2})\right|_{p}\leq
2\frac{\operatorname{sinh}(\sqrt{\kappa}\bar{r})}{\sqrt{\kappa}\bar{r}}c_{p}\,\|v_{1}-v_{2}\|.$

This yields the conclusion.
Lemma C.1 and the Cartan-Hadamard Theorem ((Jost, 2017, Corollary 6.9.1)),
which states that the Riemannian exponential map is a global diffeomorphism
with $1$-Lipschitz inverse on a simply connected complete Riemannian manifold
of non-positive sectional curvature, imply local-regularity estimates for the
feature map $\varphi$, the readout map $\rho$, and their inverses; where
$\varphi$ and $\rho$ are defined as follows

$\displaystyle\rho(y)$ $\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{Exp}^{g}_{\bar{y}}\circ\iota^{-1}_{\mathcal{M},\bar{y}}(y)$
(C.3) $\displaystyle\varphi(x_{1},\dots,x_{m})$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\big{(}\operatorname{Log}_{\bar{x}_{m}}^{h}\circ\iota_{\mathcal{N},\bar{x}_{m}}(x_{m})\big{)}_{m=0}^{M},$

where $x_{1},\dots,x_{m}\in\mathbb{N}^{M}$ and $y\in\mathcal{M}$.

######  Lemma C.2 (Regularity of feature map ($\varphi$))

Suppose that $(\mathcal{N},h)$ satisfies Assumption 12, $M\in\mathbb{N}$,
$\bar{x}_{[0:\operatorname{{H}}]}\in\mathcal{N}^{\operatorname{{H}}}$, and
$\varphi$ is as in (C.3). Then
$\varphi:\mathcal{N}^{1+\operatorname{{H}}}\rightarrow\mathbb{R}^{(1+\operatorname{{H}})\,d}$
is

  1. (i)

smooth and $C_{\bar{x}}$-Lipschitz (globally) for some $C_{\bar{x}}>0$
depending only on $\bar{x}$ (see below),

  2. (ii)

$\varphi^{-1}$ exists, it is smooth, and for any
$\mathcal{K}_{[0:\operatorname{{H}}]}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\prod_{i=0}^{\operatorname{{H}}}\,\mathcal{K}_{i}\subseteq\mathbb{R}^{(1+\operatorname{{H}})\,d}$
where $\mathcal{K}_{0},\dots,\mathcal{K}_{\operatorname{{H}}}$ are non-empty
compact sets, we have

$\operatorname{Lip}(\varphi^{-1}|\mathcal{K}_{[0:\operatorname{{H}}]})\leq\max_{\tilde{m}=0,\dots,\operatorname{{H}}}c_{\bar{x}_{m}}\,\frac{\operatorname{sinh}(\kappa^{1/2}c_{\bar{x}_{m}}\,\|K_{m}\|)}{\kappa^{1/2}c_{\bar{x}_{m}}\,\|K_{m}\|},$

where $c_{\bar{x}_{m}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\|\iota_{\bar{x}_{m}}\|_{op}$ for $m=0,\dots,\operatorname{{H}}$.

Where for any non-empty subset $A\subseteq\mathbb{R}^{d}$ we define
$\|A\|\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\sup_{u\in A}\,\|u\|$
where $C_{\bar{x}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\max\limits_{m=0,\dots,\operatorname{{H}}}|\iota_{\bar{x}_{m}}^{-1}|_{op}$.

Proof of Lemma C.2 Since $\mathcal{N}$ is simply connected, complete, and has
non-positive sectional curvatures then the first statement in the Cartan-
Hadamard theorem, as formulated in (Jost, 2017, Corollary 6.9.1), guarantees
that each
$\operatorname{Log}^{h}_{\bar{x}_{m}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}(\operatorname{Exp}^{h}_{\bar{x}_{m}})^{-1}:\mathcal{N}\to
T_{\bar{x}_{m}}\mathcal{N}$, for $m=0,\dots,\operatorname{{H}}$ is well-
defined and smooth. Moreover, the second statement of (Jost, 2017, Corollary
6.9.1) implies that each $\operatorname{Log}^{h}_{\bar{x}_{m}}$ is a
$1$-Lipschitz map (globally on $(\mathcal{N},{d_{h}})$). Consequentially, for
every
$x_{[0:\operatorname{{H}}]},\tilde{x}_{[0:\operatorname{{H}}]}\in\mathcal{N}^{\operatorname{{H}}}$
it follows that

$\displaystyle\|\varphi(x_{[1:\operatorname{{H}}]})-\varphi(\tilde{x}_{[1:M]})\|=$
$\displaystyle\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,\big{\|}\operatorname{Log}_{\bar{x}_{m}}{\circ\iota^{-1}_{\bar{x}_{m}}}(x_{r})-\operatorname{Log}_{\bar{x}_{m}}{\circ\iota^{-1}_{\bar{x}_{m}}}(\tilde{x}_{r})\big{\|}^{2}_{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle{\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,[|\iota^{-1}_{\bar{x}_{m}}|\,d_{h}(x_{r},\tilde{x}_{r})]^{2}\biggr{)}^{1/2}}$
$\displaystyle\leq$
$\displaystyle{C_{\bar{x}}\,\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,d_{h}(x_{r},\tilde{x}_{r})^{2}\biggr{)}^{1/2}}$
$\displaystyle=$
$\displaystyle{C_{\bar{x}}}\,d_{h:2}(x_{[1:\operatorname{{H}}]},\tilde{x}_{[1:M]}).$

This yields claim (i), where $d_{h:2}$ is the $(1+\operatorname{{H}})$-fold
$2$-product metric on $\mathcal{N}^{1+\operatorname{{H}}}$ of the geodesic
distance $d_{h}$ on $(\mathcal{N},h)$.

Let
$u_{[0:\operatorname{{H}}]},\tilde{u}_{[0:\operatorname{{H}}]}\in(\mathbb{R}^{d})^{1+\operatorname{{H}}}\cong\mathbb{R}^{{(1+\operatorname{{H}})}\,d}$.
Under the sectional curvature lower-bound in Assumption 12 (iii), we may apply
Lemma C.1 to deduce that, for any $r_{{0}},\dots,r_{\operatorname{{H}}}>0$ and
each $m=0,\dots,\operatorname{{H}}$, the restriction of the map
$\operatorname{Exp}_{\bar{x}_{m}}$ to the closed Euclidean ball
$\overline{\operatorname{Ball}_{(\mathbb{R}^{d},\ell^{2}_{d})}(0,c_{\bar{x}_{m}}r_{m}/c_{\bar{x}_{m}})}$
is Lipschitz with Lipschitz constant at-most
$c_{\bar{x}_{m}}\,\frac{\operatorname{sinh}(\kappa^{1/2}c_{\bar{x}_{m}}r_{m})}{\kappa^{1/2}c_{\bar{x}_{m}}r_{m}}$.
Therefore, we have that

$\displaystyle
d_{h:2}(\varphi^{-1}(u_{[0:\operatorname{{H}}]}),\varphi^{-1}(\tilde{u}_{[0:\operatorname{{H}}]}))=$
$\displaystyle\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,d_{h}\big{(}\operatorname{Exp}_{\bar{x}_{m}}(u_{m}),\operatorname{Exp}_{\bar{x}_{m}}(\tilde{u}_{m})\big{)}^{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,\Big{(}c_{\bar{x}_{m}}\frac{\operatorname{sinh}(\kappa^{1/2}c_{\bar{x}_{m}}r_{m})}{\kappa^{1/2}c_{\bar{x}_{m}}r_{m}}\,\|u_{m}-\tilde{u}_{m}\|\Big{)}^{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,\biggl{[}\max_{\tilde{m}=0,\dots,\operatorname{{H}}}c_{\bar{x}_{m}}\,\Big{(}\frac{\operatorname{sinh}(\kappa^{1/2}r_{\tilde{m}})}{\kappa^{1/2}r_{\tilde{m}}}\Big{)}\biggr{]}^{2}\,\|u_{m}-\tilde{u}_{m}\|^{2}\biggr{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle\biggl{[}\max_{\tilde{m}=0,\dots,\operatorname{{H}}}c_{\bar{x}_{m}}\,\frac{\operatorname{sinh}(\kappa^{1/2}r_{\tilde{m}})}{\kappa^{1/2}r_{\tilde{m}}}\biggr{]}\,\biggl{(}\sum_{m=0}^{\operatorname{{H}}}\,\|u_{m}-\tilde{u}_{m}\|^{2}\biggr{)}^{1/2}$
$\displaystyle=$
$\displaystyle\biggl{[}\max_{\tilde{m}=0,\dots,\operatorname{{H}}}c_{\bar{x}_{m}}\,\frac{\operatorname{sinh}(\kappa^{1/2}r_{\tilde{m}})}{\kappa^{1/2}r_{\tilde{m}}}\biggr{]}\,\|u_{[0:\operatorname{{H}}]}-\tilde{u}_{[0:\operatorname{{H}}]}\|.$

Note that each
$K_{m}\subseteq\overline{\operatorname{Ball}_{(\mathbb{R}^{d},\ell^{2}_{d})}(0,\|K_{m}\|))}$,
for $m=0,\dots,\operatorname{{H}}$. Upon setting
$r_{m}\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}{\|K_{m}\|}$, for each
$m=0,\dots,\operatorname{{H}}$, we obtain the conclusion.

###### Remark C.3

Lemma C.2 shows that a poor choice of base-points $\bar{x}_{m}$ when defining
the feature map $\varphi$ will require a larger network to compensate for the
unnecessarily large Lipschitz constant of $\varphi^{-1}$.

Since $\rho$ is $\varphi^{-1}$ when interchanging the roles of $\mathcal{M}$
and $\mathcal{N}$ in Lemma (C.2), in the special case where $M=1$, then Lemma
C.2 directly implies the following.

######  Lemma C.4 (Regularity of readout map ($\rho$))

Suppose that $(\mathcal{M},g)$ satisfies Assumption 12 (i) and (iii). Fix a
$\bar{y}\in\mathcal{M}$ and define $\rho$ as in (C.3). Then,
$\rho:\mathbb{R}^{D}\rightarrow\mathcal{M}$ satisfies

  1. (i)

$\rho$ is smooth and on any non-empty compact $K\subseteq\mathbb{R}^{D}$ its
local-Lipschitz constant is

$\operatorname{Lip}(\rho|K)\leq
c_{\bar{y}}\frac{c_{\bar{y}}\operatorname{sinh}(\kappa^{1/2}\operatorname{diam}(K))}{c_{\bar{y}}\kappa^{1/2}\operatorname{diam}(K)},$

where $c_{\bar{y}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}2\,\|\iota_{\bar{y}}\|_{op}$.

  2. (ii)

$\rho^{-1}$ exists, it is smooth, and it is ${C_{p}}$-Lipschitz (globally)
where where $C_{\bar{p}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}|\iota_{p}^{-1}|_{op}$.

Proof  Point (i) follows from Lemma C.1, by setting
$\bar{r}=c_{\bar{y}}\,\operatorname{diam}(K)$ and noting that
$x\mapsto\operatorname{sinh}(\kappa^{1/2}x)/(\kappa^{1/2}x)$ is monotonically
increasing on $[0,\infty)$. Point (ii) follows from the second statement of
the Cartan-Hadamard Theorem (Jost, 2017, Corollary 6.9.1) and upon recalling
that we have linearly isomorphically identified the Hilbert spaces
$\iota_{p}:(T_{p}(\mathcal{M}),\|\cdot\|_{g_{p}})\to(\mathbb{R}^{D},\|\cdot\|)$
(which needn’t be an isometry); whence
$\operatorname{Lip}(\iota_{p}^{-1})=C_{\bar{p}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}|\iota_{p}^{-1}|_{op}>0$ (and need not equal to $1$).

###  C.2 Results from matrix theory

This appendix contains a short proof of some folklore results in matrix
theory.

Recall that, for a negative integer $t$ and an invertible matrix $\Sigma$, the
matrix $\Sigma^{t}$ is defined to be $(\Sigma^{-1})^{|t|}$. That is,
$\Sigma^{t}$ is the $t$-fold product of $\Sigma^{-1}$ with itself.

The next lemma justifies the definition of the real exponent of a symmetric
positive definite matrix.

######  Lemma C.5 (Powers of symmetric positive-definite matrices)

For each $\Sigma\in\operatorname{Sym}_{+}(d)$ and each $t\in\mathbb{R}$,
$\exp(t\log(\Sigma))$ is a well-defined element of
$\operatorname{Sym}_{+}(d)$. Furthermore, if $t\in\mathbb{Z}$ then

$\exp(t\log(\Sigma))=\Sigma^{t}.$

Proof of Lemma C.5 By (Hilgert and Neeb, 2012, Proposition 3.3.5) the
restriction of the matrix exponential map to $\operatorname{Sym}(d)$ is a
diffeomorphism onto the open subset $\operatorname{Sym}_{+}(d)$ of
$\operatorname{Sym}(d)$, and by definition, the matrix logarithm $\log$ is the
inverse of $\exp$ on $\operatorname{Sym}_{+}(d)$. Therefore,
$X\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\log(\Sigma)$ is a well-
defined element of $\operatorname{Sym}(d)$. Consequentially, $tX$ belongs to
$\operatorname{Sym}(d)$ for every $t\in\mathbb{R}$. Thus,
$\exp(tX)\in\operatorname{Sym}_{+}(d)$ nd the first claim follows.

Let now be $t\in\mathbb{Z}$ with $t\neq 0$ (otherwise the statement is
trivial). If $t>0$, then, since $X$ trivially commutes with itself, it easily
follows that

$\exp(tX)=\underbrace{\exp(X)\dots\exp(X)}_{\mbox{$t$-times}}=(\exp(X))^{t}.$

If on the other hand $t<0$, then, from we have just proved,
$\exp(tX)=\exp(-t(-X))=(\exp(-X))^{\left|t\right|}$. But since $\exp(X)$ is
invertible with $\exp(-X)=\exp(X)^{-1}$, from the observation at the beginning
of this section we infer $\exp(tX)=(\exp(X))^{t}$

So, we have that $\exp(tX)=\big{(}\exp(X)\big{)}^{t}$ for every
$t\in\mathbb{Z}$. Plugging in the definition of $X$ and using the fact that
$\log$ is the inverse of $\exp$ on $\operatorname{Sym}_{+}(d)$ yields the
conclusion, namely that

$\exp(t\log(\Sigma))=\big{(}\exp(\log(\Sigma))\big{)}^{t}=\Sigma^{t},$

concluding our proof.

##  Appendix D Proofs

This section contains the proofs of our main results.

###  D.1 Proof of Proposition 7

In order to prove Proposition 7, we need some preliminary results.

####  D.1.1 The global NPC geometry on the space of non-singular Gaussian
measures

The following result is well-known in the metric geometry literature; however,
parts of its proof were difficult to track down outside that literature, so we
record it here.

######  Lemma D.1 (The Geometry of $\operatorname{Sym}_{+}$)

For every $d\in\mathbb{N}$, $\operatorname{Sym}_{+}(d)$ is a
$d(d+1)/2$-dimensional Cartan-Hadamard manifold whose tangent spaces are
identifiable with $\operatorname{Sym}(d)$. The Riemannian metric is given by

$\hat{g}_{\Sigma}(X,Y)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{tr}(\Sigma^{-1}X\Sigma^{-1}Y),\quad\Sigma\in\operatorname{Sym}_{+}(d),\,X,Y\in\operatorname{Sym}(d).$

Moreover, the sectional curvatures are bounded in $[-1/2,0]$.

Proof of Lemma D.1 As discussed on (Bridson and Haefliger, 1999, page 314,
section 10.31), $\operatorname{Sym}_{+}(d)$ is a smooth manifold of dimension
$d(d+1)/2$ with tangent spaces identifiable with $\operatorname{Sym}(d)$ (see
(Petersen et al., 2008, Section 10.2.2)): a smooth global chart is provided by
$\operatorname{vec}:\operatorname{Sym}_{+}(d)\to\operatorname{vec}(\operatorname{Sym}_{+}(d))\subset\mathbb{R}^{d(d+1)/2}$
(recall that
$\operatorname{vec}:\mathbb{R}^{d(d+1)/2}\rightarrow\operatorname{Sym}(d)$ is
the inverse of the vector space isomorphism $\operatorname{sym}$ defined in
(3.11)); moreover $\operatorname{Sym}_{+}(d)$ can be endowed with the
Riemannian metric $\hat{g}_{\Sigma}(X,Y)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{tr}(\Sigma^{-1}X\Sigma^{-1}Y)$, where
$\Sigma\in\operatorname{Sym}_{+}(d)$ and $X,Y\in\operatorname{Sym}(d)$.

By Theorem (Bridson and Haefliger, 1999, 10.39),
$(\operatorname{Sym}_{+},d_{\hat{g}})$ is a $\operatorname{CAT}(0)$
space101010See (Bridson and Haefliger, 1999, Definition II.1.1). and since
$(\operatorname{Sym}_{+},\hat{g})$ is a Riemannian manifold then (Bridson and
Haefliger, 1999, Theorem II.1A.6)111111Due to Cartan (1988). implies that the
sectional curvatures of $(\operatorname{Sym}_{+},\hat{g})$ are all bounded-
above by $0$. On the other hand, (Criscitiello and Boumal, 2022, Proposition
I.1)121212In their extended appendix in the ArXiV version. we see that the
sectional curvatures in $(\operatorname{Sym}_{+},\hat{g})$ are bounded-below
by $-1/2$.

In what follows, for any $\mathfrak{m}\in\mathbb{R}^{d}$, we let
$\mathcal{N}_{d}^{\mathfrak{m}}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\\{\mathcal{N}_{d}(\mathfrak{m},\Sigma):\,\Sigma\in\operatorname{Sym}_{+}(d)\\}$.
Since $\mathcal{N}^{\mathfrak{m}}_{d}$ is a Riemannian submanifold of
dimension $d(d+1)/2$ of $(\mathcal{N}_{d},\mathcal{I}^{F})$ (with metric
pullbacked via the inclusion map $\iota$), by defining

$\Phi:\operatorname{Sym}_{+}(d)\to\mathcal{N}_{d}^{\mathfrak{m}},\quad\Sigma\mapsto\mathcal{N}(\mathfrak{m},\Sigma)$
(D.1)

it is easy to see that $\Phi$ is an isometry between
$(\operatorname{Sym}_{+}(d),\hat{g})$ and
$(\mathcal{N}^{\mathfrak{m}}_{d},\iota^{\ast}\mathcal{I}^{F})$: see also
(3.7). Lemma D.1 this implies the following.

######  Corollary D.2 (The geometry of $\mathcal{N}^{\mathfrak{m}}_{d}$)

For every $d\in\mathbb{N}$ and each $\mathfrak{m}\in\mathbb{R}^{d}$,
$\mathcal{N}^{\mathfrak{m}}_{d}$ is a $d(d+1)/2$-dimensional Cartan-Hadamard
manifold. Moreover, their sectional curvatures are bounded in $[-1/2,0]$. In