remarkRemark Approximating maps into manifolds with lower curvature boundsS.
Jacobsson et al.

# Approximating maps into manifolds with lower curvature bounds††thanks:
Funding: This project was funded by BOF project C16/21/002 by the Internal
Funds KU Leuven and FWO project G080822N. J. Van der Veken is additionally
supported by the Research Foundation–Flanders (FWO) and the Fonds de la
Recherche Scientifique (FNRS) under EOS Project G0I2222N. R. Vandebril is
additionally supported by the Fund for Scientific Research–Flanders (Belgium),
projects G0A9923N and G0B0123N.

Simon Jacobsson KU Leuven, Department of Computer Science, Celestijnenlaan
200A – box 2402, B-3000 Leuven, Belgium<EMAIL_ADDRESS>Raf
Vandebril KU Leuven, Department of Computer Science, Celestijnenlaan 200A –
box 2402, B-3000 Leuven, Belgium<EMAIL_ADDRESS>Joeri van der
Veken KU Leuven, Department of Mathematics, Celestijnenlaan 200B – box 2400,
B-3000, Leuven, Belgium<EMAIL_ADDRESS>Nick Vannieuwenhoven
KU Leuven, Department of Computer Science, Celestijnenlaan 200A – box 2402,
B-3000 Leuven, Belgium (); Leuven.AI – KU Leuven Institute for AI, B-3000
Leuven, Belgium<EMAIL_ADDRESS>
###### Abstract

Many interesting functions arising in applications map into Riemannian
manifolds. We present an algorithm, using the manifold exponential and
logarithm, for approximating such functions. Our approach extends
approximation techniques for functions into linear spaces so that we can upper
bound the forward error in terms of a lower bound on the manifold’s sectional
curvature. Furthermore, when the sectional curvature of a manifold is
nonnegative, such as for compact Lie groups, the error is guaranteed to not be
worse than in the linear case. We implement the algorithm in a Julia package
and apply it to two example problems from Krylov subspaces and dynamic low-
rank approximation, respectively. For these examples, the maps are confirmed
to be well approximated by our algorithm.

###### keywords:

Riemannian manifold, manifold exponential, manifold logarithm, function
approximation, constructive approximation.

15A69, 53Z99, 65D15

##  1 Introduction

Approximation theory is concerned with approximating functions with simpler
functions. Classic approximation theory studies schemes for approximating
real-valued functions. For example, univariate schemes like interpolation or
regression [23, 62], but also multivariate schemes like hyperbolic crosses,
sparse grids, and greedy approximation [61]. More recently, multivariate
function approximation methods based on tensor decompositions have been
extensively investigated [12, 13, 14, 39, 38, 32, 31, 30, 58, 57]. All of
these schemes are naturally extended to approximate vector-valued functions by
applying them component-wise.

Many interesting functions arising in applications have some further
structure, like matrix-valued functions having _low-rank_ or _orthogonal_
values. For instance, the present work was conceived to address difficulties
encountered in the context of reduced order models of a microstructure
materials engineering problem in [27]. However, naively applying
aforementioned classic approximation schemes fails to preserve such
structures. To remedy this, we introduce a new scheme to approximate maps
whose values are points on a _Riemannian manifold_. Our scheme works by
pulling back the approximation problem to the tangent space. In this way, the
problem is reduced again to the classic function approximation problem between
vector spaces. Consequently, our scheme can leverage all of the previously
mentioned classic techniques.

A pullback to a tangent space is a natural idea to solve problems on manifolds
in general. In the context of interpolation, some previous related work based
on this idea is for example Zimmermann [68] and Zimmermann and Bergmann [70],
who solve uni- and multivariate Hermite interpolation problems on manifolds.
Bergmann and Gousenbourger [11] and Gousenbourger, Massart, and Absil [33]
define manifold-valued generalizations of Bézier curves. While these schemes
could be used in some cases to approximate maps into manifolds, they were
designed as interpolation methods. We focus on _approximation_ rather than
interpolation.

###  1.1 Contribution

We propose to approximate a map $f\colon R\to M$, where $R$ may be an
arbitrary space and $M$ is a Riemannian manifold, by a natural, generally
applicable three-step template:

  1. 1.

Choose a point $p\in M$ from which to linearize the approximation problem.

  2. 2.

Pull $f$ back to the tangent space $T_{p}M$ using a _normal coordinate chart_
[44], and approximate the resulting map $g:={\log_{p}}\circ f\colon R\to
T_{p}M$ by $\widehat{g}\colon R\to T_{p}M$ using an arbitrary approximation
scheme for maps into vector spaces.

  3. 3.

Push $\widehat{g}$ forward with $\operatorname{exp}_{p}$ to yield
$\widehat{f}:={\operatorname{exp}_{p}}\circ\widehat{g}\colon R\to M$.

Alternatively, if $\operatorname{exp}_{p}$ and $\operatorname{log}_{p}$ are
not known explicitly, a _retraction_ [1, 15] may be used instead. Algorithm 1
is a detailed description of a concrete algorithm that results from applying
this template to maps with domain $R=[-1,1]^{m}$.

Figure 1 is a small example of Algorithm 1 approximating a map into the sphere
$S^{2}$, defined by stereographic projection of
$(x,y)\mapsto(x^{2}-y^{2},2xy)$.

(a) $9$ samples.

(b) $16$ samples.

(c) $25$ samples.

Figure 1: Using different number of samples, we approximate
$f\colon[-1,1]^{2}\to S^{2}$. Let $\Gamma$ be a $10\times 10$ grid on
$[-1,1]^{2}$. Solid lines show $f(\Gamma)$ and dashed lines show
$\widehat{f}(\Gamma)$.

While the proposed template is conceptually straightforward, the main
challenge is establishing a bound on the approximation error, as any
discrepancy between $g$ and $\widehat{g}$ will be propagated through the
exponential map. This requires understanding how the endpoint of a geodesic,
which is a solution of a particular initial value problem on a manifold,
changes as the initial conditions of this differential equation are varied.
Theorem 3.1 shows how _Toponogov’s theorem_ can be leveraged to derive upper
bounds on this perturbation.

###  1.2 Outline

Toponogov’s theorem is a well-known Riemannian comparison theorem. In Section
2, we recall some geometric preliminaries to be able to utilize it.

Section 3 describes the error analysis of our algorithm template, where our
main original contribution is Theorem 3.1. Our analysis assumes a lower bound
on the _sectional curvature_ of the manifold. For completeness and ease of
reference, we lists several concrete manifolds along with explicit bounds on
their curvature in Appendix A. We are especially interested in manifolds with
nonnegative curvature, as their error bound is guaranteed not to be worse than
in the familiar linear case.

Inspired by Dolgov, Kressner, and Strössner [58] and by Strössner, Sun, and
Kressner [57], we use _tensorized Chebyshev interpolation_ in step 2 of our
template to implement Algorithm 1. We discuss the rest of the implementation
details and recall error bounds for tensorized Chebyshev interpolation in
Section 4.

We present two example approximation problems inspired by linear algebra
applications in Section 5. In these examples, the approximation error is
verified to be bounded by Theorem 3.1. Our implementation is available as a
Julia package, ManiFactor.jl [40].

Finally, the paper is concluded with a summary of the main results and a
future outlook in Section 6.

## Acknowledgements

We thank Arne Bouillon for spotting an error in Eq. 10 in an earlier draft,
and Astrid Herremans for her expertise in approximation theory and her
suggestions about the figures.

## Notation

The following notation will be used throughout the paper:

  * •

$\mathrm{id}\colon V\to V$ is the _identity map_.

  * •

$\norm{\cdot}_{\infty}$ is the _infinity norm_. For arrays
$A\in\mathbb{R}^{N_{1}\times\dots\times N_{m}}$, it is defined as the element
with maximum absolute value. For functions $g\colon R\to\mathbb{R}$, it is
defined as $\sup_{x\in R}\absolutevalue{g(x)}$. If $O$ is an operator, we will
also use $\norm{O}_{\infty}$ to mean the operator norm with respect to the
infinity norm.

  * •

$\angle(u,v)$ is the _angle_ between two vectors $u$ and $v$ in an inner
product space $(V,\left<{\cdot},{\cdot}\right>)$ with induced norm
$\norm{\cdot}$. It is defined as
$\acos(\left<{u/\norm{u}},{v/\norm{v}}\right>)$.

  * •

For a Riemannian manifold $M$, $\left<{\cdot},{\cdot}\right>_{p}$ is the inner
product on the tangent space $T_{p}M$ at $p\in M$, $\norm{\cdot}_{p}$ is the
induced norm on $T_{p}M$, and $d_{M}(\cdot,\cdot)$ is the induced geodesic
distance on $M$. $\operatorname{exp}_{p}$ and $\operatorname{log}_{p}$ are the
exponential and logarithmic maps at $p$ respectively.

##  2 Geometric preliminaries

We now present a high-level overview of the key concepts from Riemannian
geometry that we will use. The technical details, which the present work does
not require, can be found in standard references, such as Lee [43, 44].

Manifolds are abundant in numerical analysis. Examples are sets of orthogonal
matrices, positive definite matrices, upper triangular matrices, tensors of
fixed rank, and linear subspaces of fixed dimension. Morally, a _manifold_ is
a space $M$ that locally can be identified with open subsets of
$\mathbb{R}^{n}$. Such an identification is called a _chart_. A detailed
introduction to manifolds can be found in Lee [43, chapters 1–3].

A _Riemannian manifold_ $M$ is a manifold with a notion of distance between
points. A _geodesic_ is a curve in $M$ that is locally distance-minimizing.
The _tangent space_ $T_{p}M$ to $M$ at a point $p\in M$ is the space of
velocity vectors of curves in $M$ through $p$. It is an $n$-dimensional vector
space. Each tangent vector is associated with a geodesic through $p$, and vice
versa. This gives a canonical chart around each point $p\in M$, called a
_normal coordinate chart_. It consists of the _manifold exponential_ and its
inverse, the _manifold logarithm_ , which identify a neighbourhood $S\subset
T_{p}M$ of $0$ with a neighbourhood $S^{\prime}\subset M$ of $p$. We denote
them by

(1) $\displaystyle\operatorname{exp}_{p}\colon S\to
S^{\prime},\qquad\text{and}\qquad\operatorname{log}_{p}\colon S^{\prime}\to
S.$

A detailed introduction to Riemannian geometry is given in Lee [44, chapters
1–8].

A Riemannian manifold $M$ is _geodesically convex_ if for any two points on
$M$, there is a unique distance-minimizing geodesic connecting them. Such a
geodesic is called _minimal_. A _geodesic triangle_ is a set of three minimal
geodesics connecting three points on $M$.

A _Lie group_ is a manifold equipped with a group operation. Specifically,
Riemannian Lie groups whose group operation is compatible with their metric
have many desirable properties. They thus form a main class of examples in
Appendix A. A concise introduction to Lie groups can be found in
Arvanitoyeorgos [6].

###  2.1 Sectional curvature

If $M$ is a two-dimensional Riemannian manifold, then we can define the
_Gaussian curvature_ $K$ at a point $p\in M$ in several ways. One way is to
define it by the angular defect of small triangles:

(2) $\displaystyle\alpha+\beta+\gamma-\pi=KA+\textrm{higher order terms},$

where $\alpha$, $\beta$, $\gamma$ are the angles of a geodesic triangle in $M$
at $p$ with area $A$.

If $M$ is not two-dimensional, consider a two-dimensional linear subspace
$\Pi$ of $T_{p}M$. Then $\Sigma=\operatorname{exp}_{p}(\Pi)$ is a two-
dimensional submanifold of $M$. The Gaussian curvature of $\Sigma$ at $p$ is
called the _sectional curvature_ of $M$ along $\Pi$.

A manifold with constant sectional curvature $H$ is called a _model manifold_.
It is uniquely determined by $H$ and its dimension.

###  2.2 Riemannian comparison theory

Riemannian comparison theory studies inequalities involving geometric
quantities. A good reference is Cheeger and Ebin [21]. Many comparison
theorems relate properties of geodesics to curvature bounds. _Toponogov’s
theorem_ is such a theorem, bounding side lengths of geodesic triangles
defined on manifolds with lower bounded sectional curvature. It is illustrated
in Fig. 2.

$\leadsto$$M$

$|$

$||$

$|||$$N$

$|$

$||$

$|||$ Figure 2: Proposition 2.1 gives conditions for which the length of the
opposite (dashed) side in $M$ is bounded by the length of the opposite
(dashed) side in $N$. Equal length geodesics and equal angles have the same
number of notches.

######  Proposition 2.1 (Toponogov’s theorem [21, Theorem 2.2]).

Let $M$ be a geodesically convex Riemannian manifold with sectional curvature
bounded from below by some constant $H$. Let $\gamma_{1}$,
$\gamma_{2}\colon[0,1]\to M$ be minimal geodesics such that
$\gamma_{1}(0)=\gamma_{2}(0)$. If $H>0$, also assume that
$d_{M}(\gamma_{1}(0),\gamma_{1}(1))$,
$d_{M}(\gamma_{2}(0),\gamma_{2}(1))\leq\frac{\pi}{\sqrt{H}}$. Consider a model
manifold $N$ of constant sectional curvature $H$ and let $\lambda_{1}$,
$\lambda_{2}\colon[0,1]\to N$ be geodesics such that
$\lambda_{1}(0)=\lambda_{2}(0)$ and the angle where they meet satisfies
$\angle(\lambda_{1},\lambda_{2})=\angle(\gamma_{1},\gamma_{2})$. Then

(3) $\displaystyle d_{M}(\gamma_{1}(1),\gamma_{2}(1))\leq
d_{N}(\lambda_{1}(1),\lambda_{2}(1)).$

The advantage of looking at geodesic triangles on the model manifold is that
there are spherical and hyperbolic versions of the trigonometric identities.
These identities are used in Section 3 to find explicit error bounds when
approximating maps into manifolds.

######  Lemma 2.2 (Reid [50, sections 3.2 and 3.10]).

Let $N$ be a manifold of constant sectional curvature $H$ and consider a
geodesic triangle on $N$ with side lengths $A$, $B$, $C$ and opposite angles
$a$, $b$, $c$ respectively. Then

(4) $\displaystyle\cos(C\sqrt{H})={}$
$\displaystyle\cos(A\sqrt{H})\cos(B\sqrt{H})+\sin(A\sqrt{H})\sin(B\sqrt{H})\cos
c$

if $H>0$, and

(5) $\displaystyle\cosh(C\sqrt{\absolutevalue{H}})={}$
$\displaystyle\cosh(A\sqrt{\absolutevalue{H}})\cosh(B\sqrt{\absolutevalue{H}})-\sinh(A\sqrt{\absolutevalue{H}})\sinh(B\sqrt{\absolutevalue{H}})\cos
c$

if $H<0$.

##  3 Error analysis

The distance between two points on a Riemannian manifold depends on the
geometry of that manifold. For example, the distance between two orthogonal
matrices $A$ and $B$ in the space of orthogonal matrices will always be
greater than or equal to the distance between $A$ and $B$ in the ambient space
of matrices. When approximating maps into Riemannian manifolds it is natural
to measure the approximation error _intrinsically_ on the manifold, rather
than in some _extrinsic_ ambient space in which the manifold could be
embedded. Moreover, for embedded manifolds with large codimension, such as
low-rank matrices or tensors, it can be inefficient or even infeasible to
measure the error in the ambient space. The next result thus bounds the
intrinsic error of the approximation template from Section 1.1.

######  Theorem 3.1.

Let $M$ be a geodesically convex Riemannian manifold with sectional curvature
bounded from below by some constant $H$ and let $f\colon R\to M$, where $R$ is
an arbitrary set such that its image fits in a single normal coordinate chart
$S\subset T_{p}M$ around $p\in M$. Assume that the elements in $S$ have norm
bounded by some constant $\sigma$. Let

(6) $\displaystyle g={\operatorname{log}_{p}}\circ f\colon R\to
S\quad\text{and}\quad\widehat{f}={\operatorname{exp}_{p}}\circ\widehat{g}\colon
R\to M,$

where $\widehat{g}\colon R\to S$ is an approximation of $g$ such that

(7) $\displaystyle\norm{g(x)-\widehat{g}(x)}_{p}\leq\epsilon$

for all $x\in R$. Then, the distance between $f(x)$ and $\widehat{f}(x)$ on
$M$ obeys

(8) $\displaystyle d_{M}(f(x),\widehat{f}(x))\leq{}$
$\displaystyle\epsilon\quad$ $\displaystyle\text{if }H\geq 0,$ (9)
$\displaystyle d_{M}(f(x),\widehat{f}(x))\leq{}$
$\displaystyle\epsilon+\frac{2}{\sqrt{\absolutevalue{H}}}\operatorname{arcsinh}\frac{\epsilon\sinh(\sigma\sqrt{\absolutevalue{H}})}{2\sigma}\quad$
$\displaystyle\text{if }H<0,$

for all $x\in R$.

######  Remark 3.2.

We can in turn upper bound Eq. 9 by

(10) $\displaystyle d_{M}(f(x),\widehat{f}(x))\leq{}$
$\displaystyle\epsilon+\frac{2}{\sqrt{\absolutevalue{H}}}\log\left(\frac{\epsilon\exp(\sigma\sqrt{\absolutevalue{H}})}{2\sigma}+1\right).$

The assumptions in Theorem 3.1 about geodesic convexity, bounded curvature,
and bounded normal coordinate chart may seem very specific, but we note that
they are always satisfied in some neighbourhood of $p$. To formulate explicit
bounds, we however need explicit expressions for $H$. Therefore, we list
several standard manifolds along with explicit lower bounds for their
sectional curvature in Appendix A.

######  Proof 3.3 (Proof of Theorem 3.1).

The proof can be summarized as applying Toponogov’s theorem to the geodesic
triangle $(f(x),\widehat{f}(x),p)$, as visualized in Fig. 3.

$\leadsto$$M$$p$$f(x)$$\widehat{f}(x)$

$|$

$||$

$|||$$N$$q$$h(x)$$\widehat{h}(x)$

$|$

$||$

$|||$ Figure 3: Geodesic triangles defined by $p$, $f$, and $\widehat{f}$. The
map from $M$ to $N$ is
$\operatorname{exp}_{N,q}\circ\operatorname{log}_{M,p}$.

Let $N$ be the manifold of constant curvature $H$ and let $q\in N$. Fix some
isometry $T_{p}M\sim T_{q}N$ and let $\operatorname{exp}_{N,q}\colon S\to N$
be the manifold exponential on $N$. Define $h={\operatorname{exp}_{N,q}}\circ
g$ and $\widehat{h}={\operatorname{exp}_{N,q}}\circ\widehat{g}$. By
Proposition 2.1,

(11) $\displaystyle d_{M}(f(x),\widehat{f}(x))\leq
d_{N}(h(x),\widehat{h}(x)).$

Now our task is to bound the right-hand side of this inequality.

For brevity, from now on we drop the argument $x$ of the functions $f$, $g$,
and $h$.

#### Case 1: $H\geq 0$

Without loss of generality, let $H=0$. Then $N=T_{p}M$ and so

(12) $\displaystyle d_{N}(h,\widehat{h})={}$
$\displaystyle\norm{g-\widehat{g}}_{p}\leq\epsilon.$

#### Case 2: $H<0$

Consider the geodesic triangle $(h(x),\widehat{h}(x),q)$. Define
$c(t)=\cosh(t\sqrt{\absolutevalue{H}})$ and
$s(t)=\sinh(t\sqrt{\absolutevalue{H}})$. Then Lemma 2.2 says that

(13) $\displaystyle c(d(h,\widehat{h}))={}$ $\displaystyle
c(\norm{g}_{p})c(\norm{\widehat{g}}_{p})-s(\norm{g}_{p})s(\norm{\widehat{g}}_{p})\cos{\left(\angle(g,\widehat{g})\right)}.$

Since ${\sinh}$ is positive on $\mathbb{R}^{+}$,

(14) $\displaystyle s(\norm{g}_{p})s(\norm{\widehat{g}}_{p})\geq 0.$

Hence using

(15)
$\displaystyle\cos{\left(\angle(g,\widehat{g})\right)}=\frac{\left<{g},{\widehat{g}}\right>_{p}}{\norm{g}_{p}\norm{\widehat{g}}_{p}}=\frac{\norm{g}_{p}^{2}+\norm{\widehat{g}}_{p}^{2}}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}-\frac{\norm{g-\widehat{g}}_{p}^{2}}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}\geq
1-\frac{\norm{g-\widehat{g}}_{p}^{2}}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}$

in Eq. 13 gives us

(16) $\displaystyle c(d(h,\widehat{h}))\leq{}$ $\displaystyle
c(\norm{g}_{p})c(\norm{\widehat{g}}_{p})-s(\norm{g}_{p})s(\norm{\widehat{g}}_{p})\left(1-\frac{\norm{g-\widehat{g}}_{p}^{2}}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}\right)$
(17) $\displaystyle={}$ $\displaystyle
c(\norm{g}_{p}-\norm{\widehat{g}}_{p})+\frac{\norm{g-\widehat{g}}_{p}^{2}s(\norm{g}_{p})s(\norm{\widehat{g}}_{p})}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}.$

Since ${\cosh}$ is even and increasing on $\mathbb{R}^{+}$, we may use the
reverse triangle inequality to see that

(18) $\displaystyle c(d(h,\widehat{h}))\leq{}$ $\displaystyle
c(\norm{g-\widehat{g}}_{p})+\frac{\norm{g-\widehat{g}}_{p}^{2}s(\norm{g}_{p})s(\norm{\widehat{g}}_{p})}{2\norm{g}_{p}\norm{\widehat{g}}_{p}}.$

Furthermore, noting that $\frac{\sinh{x}}{x}$ is increasing and positive on
$\mathbb{R}^{+}$, we have that

(19) $\displaystyle c(d(h,\widehat{h}))\leq{}$ $\displaystyle
c(\norm{g-\widehat{g}}_{p})+\frac{\norm{g-\widehat{g}}_{p}^{2}\,s(\sigma)^{2}}{2\sigma^{2}}.$

This is equivalent to

(20) $\displaystyle 1+2s\left(\frac{d(h,\widehat{h})}{2}\right)^{2}\leq{}$
$\displaystyle
1+2s\left(\frac{\norm{g-\widehat{g}}_{p}}{2}\right)^{2}+\frac{\norm{g-\widehat{g}}_{p}^{2}s(\sigma)^{2}}{2\sigma^{2}},$

so that

(21) $\displaystyle d(h,\widehat{h})\leq{}$ $\displaystyle
2{s}^{-1}\left(\sqrt{s\left(\frac{\norm{g-\widehat{g}}_{p}}{2}\right)^{2}+\frac{\absolutevalue{g-\widehat{g}}^{2}\,s(\sigma)^{2}}{4\sigma^{2}}}\right),$

and since $\operatorname{sinh}$, $\operatorname{arcsinh}$, $\cdot^{2}$, and
$\sqrt{\cdot}$ are all increasing on $\mathbb{R}^{+}$, we find

(22) $\displaystyle d(h,\widehat{h})\leq{}$ $\displaystyle
2{s}^{-1}\left(\sqrt{s\left(\frac{\epsilon}{2}\right)^{2}+\frac{\epsilon^{2}s(\sigma)^{2}}{4\sigma^{2}}}\right).$

Lastly we will use that concave functions are _subadditive_ ,
$\sqrt{a^{2}+b^{2}}\leq a+b$ and
$\operatorname{arcsinh}(\sinh{a}+\sinh{b})\leq a+b$ for nonnegative $a$ and
$b$. Thus

(23) $\displaystyle d(h,\widehat{h})\leq{}$
$\displaystyle\epsilon+2{s}^{-1}\left(\frac{\epsilon
s(\sigma)}{2\sigma}\right).$

This concludes the proof.

If the exponential or logarithmic map are not explicitly known for a manifold,
then we might still be able to proceed by using a _retraction_. This is an
injective map $r_{p}$ that approximates the exponential map at $p$. For a
precise definition, see Absil, Mahony, and Sepulchre [1, Definition 4.1.1].

Assume the same conditions as in Theorem 3.1. Moreover, assume the max error
between $r_{p}$ and $\operatorname{exp}_{p}$ is bounded by $\zeta$, assume the
max error between ${r}^{-1}_{p}$ and $\operatorname{log}_{p}$ is bounded by
$\eta$, and assume that the vector-valued approximation scheme
$\mathcal{I}\colon g\mapsto\widehat{g}$ is $\mathbb{R}$-linear and that its
operator infinity norm is $\Lambda$. Now, instead of using
$\widehat{f}=({\operatorname{exp}_{p}}\circ\mathcal{I}\circ{\operatorname{log}_{p}})(f)$,
we use $\widehat{f}=(r_{p}\circ\mathcal{I}\circ{r_{p}}^{-1})(f)$. It is then
straightforward to show the following bounds:

(24) $\displaystyle d_{M}(f(x),\widehat{f}(x))\leq{}$
$\displaystyle\epsilon+\Lambda\eta+\zeta$ if $H\geq 0$, (25) $\displaystyle
d_{M}(f(x),\widehat{f}(x))\leq{}$
$\displaystyle\epsilon+\Lambda\eta+\zeta+\frac{2}{\sqrt{\absolutevalue{H}}}\operatorname{arcsinh}\left(\frac{(\epsilon+\Lambda\eta)\sinh(\sigma\sqrt{H})}{2\sigma}\right)$
if $H<0$.

When $H\to 0^{-}$, Eq. 9 tends to

(26) $\displaystyle d_{M}(f,\widehat{f})\leq 2\epsilon,$

which shows that it is not a tight bound. To explain why, first note that Eq.
15 is an equality when $\norm{g(x)}_{p}=\norm{\widehat{g}(x)}_{p}$, but we
also used the reverse triangle inequality in Eq. 18 which is an equality when
$g$ and $\widehat{g}$ are collinear. Furthermore, in Eq. 23 we used
subadditivity for the square root, which is an equality only when one of the
terms is $0$. But for small curvatures and small $\epsilon$, the two terms in
Eq. 23 are approximately equal. One could say that we lost a factor $\sqrt{2}$
in each of these steps.

Even though the bound Eq. 9 is not tight, there is, in the following sense, no
better bound.

######  Proposition 3.4.

Recall the notation of Theorem 3.1. Assume that $M$ has constant negative
sectional curvature $H$ and let $\widehat{g}$ be an approximation to $g$ such
that the maximum error $\epsilon$ is attained at a point $x\in R$ with
$\norm{\widehat{g}(x)}_{p}=\norm{g(x)}_{p}=\sigma$. Then

(27) $\displaystyle
d_{M}(f(x),\widehat{f}(x))=\frac{2}{\sqrt{\absolutevalue{H}}}\operatorname{arcsinh}\left(\frac{\epsilon\sinh{(\sigma\sqrt{\absolutevalue{H}})}}{2\sigma}\right).$

######  Proof 3.5.

Consider again the proof of Theorem 3.1. If
$\norm{g}_{p}=\norm{\widehat{g}}_{p}$, then Eq. 15 is an equality, and we do
not need to use the reverse triangle inequality in Eq. 18, hence Eq. 20
reduces to

(28) $\displaystyle 1+2s\left(\frac{d_{M}(f,\widehat{f})}{2}\right)^{2}={}$
$\displaystyle 1+\frac{\epsilon^{2}s(\sigma)^{2}}{2\sigma^{2}}.$

This is equivalent to what we wanted to prove.

As an immediate consequence of Theorem 3.1 we can obtain bounds on the
propagation of a small perturbation of a tangent vector through the
exponential map $\mathrm{exp}_{p}$. The size of this output perturbation
relative to the input perturbation is quantified by the _condition number_.
Recall the definition of Rice’s condition number [51] of a map $\phi\colon
A\to B$ between metric spaces $(A,d_{A})$ and $(B,d_{B})$:

(29) $\displaystyle\kappa[\phi](x)=\limsup_{\begin{subarray}{c}x^{\prime}\in
A,\\\ d_{A}(x,x^{\prime})\to
0\end{subarray}}\frac{d_{B}(\phi(x),\phi(x^{\prime}))}{d_{A}(x,x^{\prime})}.$

Let $\phi=\mathrm{exp}_{p}\colon(R\to S)\to(R\to M)$, with the metrics

(30) $\displaystyle d_{R\to S}(g,\widehat{g})={}\sup_{x\in
R}\norm{g(x)-\widehat{g}(x)}_{p},\mbox{\quad and\quad}d_{R\to
M}(f,\widehat{f})={}\sup_{x\in R}d_{M}(f(x),\widehat{f}(x)).$

The condition number implies an asymptotically sharp error bound

(31) $\displaystyle d_{R\to M}\left(\mathrm{exp}_{p}\circ
g,\mathrm{exp}_{p}\circ\widehat{g}\right)\leq{}$
$\displaystyle\kappa[\mathrm{exp}_{p}](g)d_{R\to S}(g,\widehat{g})+o(d_{R\to
S}(g,\widehat{g}))$ (32) $\displaystyle={}$
$\displaystyle\kappa[\mathrm{exp}_{p}](g)\epsilon+o(\epsilon),$

We obtain the following corollary of Theorems 3.1 and 3.4.

######  Corollary 3.6.

Recall the notation of Theorem 3.1. The condition number $\kappa$ of the map
$g\mapsto f$ satisfies

(33) $\displaystyle 1\leq{}\kappa\leq{}$ $\displaystyle
1+\frac{\sinh(\sigma\sqrt{\absolutevalue{H}})}{\sigma\sqrt{\absolutevalue{H}}}\quad$
$\displaystyle\text{if }H<0,$ (34) $\displaystyle\kappa={}$ $\displaystyle
1\quad$ $\displaystyle\text{if }H\geq 0.$

Furthermore, if $M$ is a model manifold with constant negative sectional
curvature $H$, then it is also lower bounded by

(35)
$\displaystyle\frac{\sinh(\sigma\sqrt{\absolutevalue{H}})}{\sigma\sqrt{\absolutevalue{H}}}\leq\kappa.$

The condition number grows exponentially in both $\sigma$ and
$\sqrt{\absolutevalue{H}}$. For example, such innocent looking parameters as
$H=-10$ and $\sigma=10$ gives a condition number of almost $9\cdot 10^{11}$,
potentially losing about 12 digits of accuracy. In general, if we want to be
guaranteed a reasonable condition number, we should make sure that the inputs
have small norm: $\sigma\leq\frac{1}{\sqrt{\absolutevalue{H}}}$.

##  4 A concrete implementation

Henceforth, we want to approximate a function $f\colon[-1,1]^{m}\to M$, where
$M$ is a Riemannian manifold. In the notation of Section 3, $R=[-1,1]^{m}$. In
this section, we present a concrete realization of the approximation template
discussed in Section 1.1, using tensorized Chebyshev interpolation. It is
presented as Algorithm 1. We describe steps 1 and 2 of the template in detail.
The final step consists of applying the exponential map to the constructed
approximation, which requires no further discussion.

Algorithm 1 For approximating maps into manifolds.

0: $f\colon{[-1,1]}^{m}\to M$, number $N_{k}+1$ of Chebyshev nodes to use in
each direction $k=1$, …, $m$.

0: Approximant $\widehat{f}\colon[-1,1]^{m}\to M$ of $f$ satisfies the error
bound of Theorem 3.1.

1: Choose $p\in M$ as (an approximation of) the Karcher mean of a large number
of $f(x)$ for $x$ chosen from a uniform distribution on ${[-1,1]}^{m}$.

2: Let $g={\operatorname{log}_{p}}\circ f$.

3: for $k=1$, …, $m$ do

4: for $i_{k}=1$, …, $N_{k}+1$ do

5: Let $t_{i_{k}}=\cos\left\\{(2i_{k}-1)\pi/2(N_{k}+1)\right\\}$ be the
$i_{k}$th Chebyshev node.

6: for $j=1$, …, $n$ do

7: Let $G_{i_{1}\dots i_{m}j}=g_{j}(t_{i_{1}},\dots,t_{i_{m}})$.

8: end for

9: end for

10: end for

11: Let $G_{i_{1}\dots i_{m}j}=\sum_{a_{1},\dots,a_{m},b}C_{a_{1}\dots
a_{m}b}U^{(1)}_{i_{1}a_{1}}\dots U^{(m)}_{i_{m}a_{m}}V_{bj}$ be an ST-HOSVD of
$G$.

12: for $k=1$, …, $m$ do

13: Let $\widehat{h}^{(k)}_{a_{k}}$ be the degree-$N_{k}$ interpolating
polynomial satisfying
$\widehat{h}^{(k)}_{a_{k}}(t_{i_{k}})=U^{(k)}_{i_{k}a_{k}}$.

14: end for

15: for $j=1$, …, $n$ do

16: Let $\widehat{g}_{j}(x)=\sum_{a_{1},\dots,a_{m},b}C_{a_{1}\dots
a_{m}b}\widehat{h}^{(1)}_{a_{1}}(x_{1})\dots\widehat{h}^{(m)}_{a_{m}}(x_{m})V_{bj}$.

17: end for

18: Let $\widehat{f}={\operatorname{exp}_{p}}\circ\widehat{g}$.

19: return $\widehat{f}$.

###  4.1 Choosing from where to linearize

In practice, the choice of the point $p\in M$ in Theorem 3.1 affects the error
a lot. Because the radius $\sigma$ of the chart $S$ appears as an exponent in
the error bound Eq. 9, it is natural to try to minimize it. However, direct
minimization of $\sigma$ is a nontrivial Riemannian optimization problem, so
instead we consider an approximation. Consider a sequence
$(\xi_{i})_{i=1}^{N}$ of points in $R$, sampled independently from the same
distribution on $R$. Then, consider the total squared distance to $p$,

(36) $\displaystyle\sum_{i=1}^{N}d_{M}(f(\xi_{i}),p)^{2}.$

The idea is that Eq. 36 may be easier to minimize than $\sigma$ itself, while
they are correlated. A $p$ that minimizes Eq. 36 is called a _Riemannian
center of mass_ of $(f(\xi_{i}))_{i=1}^{N}$ [41]. An efficient way to
approximate a Riemannian center of mass is by successive geodesic
interpolation [20, 18, 19, 22], namely

(37) $\displaystyle m(p_{1})={}$ $\displaystyle p_{1},$ (38) $\displaystyle
m(p_{1},\dots,p_{N})={}$
$\displaystyle\operatorname{exp}_{p_{N}}\left(\frac{N-1}{N}\operatorname{log}_{p_{N}}m(p_{1},\dots,p_{N-1})\right).$

Note that, in the Euclidean case, $m(p_{1},\dots,p_{N})$ is the arithmetic
mean of $(p_{i})_{i=1}^{N}$, and is thus exactly the Riemannian center of
mass.

In our implementation, we choose $p$ as $m(f(\xi_{1}),\dots,f(\xi_{N}))$ with
$(\xi_{i})_{i=1}^{N}$ sampled from a uniform distribution on $R$.

###  4.2 Tensorized Chebyshev interpolation

For the second step in our implementation, we choose _tensorized Chebyshev
interpolation_.

Tensorizing a univariate approximation scheme is a standard approach for
multivariate approximation. Here, we first recall Schultz’s [53] approach,
specialized to the case of a linear approximation scheme. Let
$\mathcal{I}\colon([-1,1]\to\mathbb{R})\to([-1,1]\to\mathbb{R})$ be a _linear
evaluation-based scheme_ for approximating real univariate functions; that is,

  1. 1.

$\mathcal{I}$ is a linear map: $\mathcal{I}(\alpha f+\beta
g)=\alpha\mathcal{I}(f)+\beta\mathcal{I}(g)$ for all real univariate functions
$f$ and $g$ and for all real numbers $\alpha$ and $\beta$.

  2. 2.

$\mathcal{I}={\mathcal{F}}\circ{\mathcal{E}_{P}}$, where
$\mathcal{F}\colon\mathbb{R}^{N}\to([-1,1]\to\mathbb{R})$ is linear and
$\mathcal{E}_{P}\colon([-1,1]\to\mathbb{R})\to\mathbb{R}^{N}$ is the
_evaluation map_ that takes a function $f$ and maps it to its evaluations
$(f(x_{1}),\dots,f(x_{N}))$ in the sample points
$P=(x_{1},\dots,x_{N})\in[-1,1]^{N}$.

Then, the _tensor product_ $\mathcal{I}\otimes\dots\otimes\mathcal{I}$ is well
defined. Recall the definition of the tensor product of linear maps [36]. Let
$U_{1}$, …, $U_{m}$, $V_{1}$, …, $V_{m}$ be vector spaces. Given linear maps
$L_{1}\colon U_{1}\to V_{1}$, …, $L_{m}\colon U_{m}\to V_{m}$, their tensor
product is the unique linear map that satisfies

$\displaystyle L_{1}\otimes\cdots\otimes L_{m}\colon U_{1}\otimes\cdots\otimes
U_{m}$ $\displaystyle\to V_{1}\otimes\cdots\otimes V_{m},$ (39) $\displaystyle
u_{1}\otimes\cdots\otimes u_{m}$ $\displaystyle\mapsto
L_{1}(u_{1})\otimes\cdots\otimes L_{m}(u_{m}),$

for all pure tensors $u_{1}\otimes\cdots\otimes u_{m}$. The following result
bounds the error of a tensorized univariate approximation scheme.

######  Lemma 4.1 (Essentially, Schultz [53, Theorem 2.1]).

Let $g\colon[-1,1]^{m}\to\mathbb{R}$. For each $k=1$, …, $m$, let
$\mathcal{I}_{k}={\mathcal{F}_{k}}\circ{\mathcal{E}_{P_{k}}}$ be a linear
evaluation-based approximation scheme using the $N_{k}$ evaluation points
$P_{k}$ on $[-1,1]$. Consider $g$ as a map of the $k$th argument only, i.e.,
$h_{k}(t)=g(x_{1},\dots,x_{k-1},t,x_{k+1},\dots,x_{m})$ and assume that

(40)
$\displaystyle\absolutevalue{h_{k}(t)-\mathcal{I}_{k}(h_{k})(t)}\leq\epsilon_{k}$

for all $t\in[-1,1]$ and for all $x\in[-1,1]^{m}$. Let

(41) $\displaystyle
G=(\mathcal{E}_{P_{1}}\otimes\dots\otimes\mathcal{E}_{P_{m}})g\in\mathbb{R}^{N_{1}\times\dots\times
N_{m}}$

be the tensor obtained from evaluating $g$ on the product grid
$P_{1}\times\dots\times P_{m}$. Furthermore, let $\widehat{G}=G+E$ for some
$E$. If
$\widehat{g}=(\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m})\widehat{G}$,
then

(42) $\displaystyle\norm{g-\widehat{g}}_{\infty}\leq{}$
$\displaystyle\epsilon_{1}+\Lambda_{1}\epsilon_{2}+\dots+\Lambda_{1}\cdots\Lambda_{m-1}\epsilon_{m}+\Lambda_{1}\cdots\Lambda_{m}\norm{E}_{\infty},$

where $\Lambda_{k}=\norm{\mathcal{F}_{k}}_{\infty}$.

######  Proof 4.2.

We have

$\displaystyle\norm{\widehat{g}-g}_{\infty}={}$
$\displaystyle\norm{(\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m})(G+E)-g}_{\infty}$
$\displaystyle\leq{}$
$\displaystyle\norm{(\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m})G-g}_{\infty}+\norm{(\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m})E}_{\infty}$
(43) $\displaystyle\leq{}$
$\displaystyle\norm{(\mathcal{I}_{1}\otimes\dots\otimes\mathcal{I}_{m})g-g}_{\infty}+\Lambda_{1}\cdots\Lambda_{m}\norm{E}_{\infty}.$

The first term is bounded as in [53, Theorem 2.1], concluding the proof.

Recall that Chebyshev polynomial interpolation interpolates a function on
$[-1,1]$ in the _Chebyshev nodes_ $\cos((2k-1)\pi/2n)$, $k=1$, …, $n$.
Chebyshev interpolants have many nice properties. Specifically, their
approximation error is bounded, see, e.g., Trefethen [62, Theorems 7.2 and
8.2], and the operator norm of Chebyshev interpolation is also bounded, see,
e.g., Trefethen [62, Theorem 15.2]. Combining Lemma 4.1 with these bounds
gives a concrete error bound for tensorized Chebyshev interpolation.

######  Corollary 4.3.

Recall the notation of Lemma 4.1. If $g$ has absolutely continuous partial
derivatives along $x_{k}$ up to order $\nu_{k}-1$ and the total variation of
$g$’s partial derivative of order $\nu_{k}$ along $x_{k}$ is bounded by
$V_{k}$, then

$\displaystyle\norm{g-\widehat{g}}_{\infty}\leq{}$
$\displaystyle\frac{4V_{1}}{\pi\nu_{1}{(N_{1}-\nu_{1})}^{\nu_{1}}}+\frac{4\Lambda_{N_{1}}V_{2}}{\pi\nu_{2}{(N_{2}-\nu_{2})}^{\nu_{2}}}+\dots+\frac{4\Lambda_{N_{1}}\cdots\Lambda_{N_{m-1}}V_{m}}{\pi\nu_{m}{(N_{m}-\nu_{m})}^{\nu_{m}}}$
(44) $\displaystyle+\Lambda_{N_{1}}\cdots\Lambda_{N_{m}}\norm{E}_{\infty},$

where $\Lambda_{N_{k}}\leq\frac{2}{\pi}\log(N_{k}+1)+1$. Especially, if
$N_{1}=\dots=N_{m}=N$, $\nu_{1}=\dots=\nu_{m}=\nu$, and $V_{1}=\dots=V_{m}=V$,
then

(45) $\displaystyle\norm{g-\widehat{g}}_{\infty}\leq{}$
$\displaystyle\frac{4(\Lambda_{N}^{m}-1)V}{\pi\nu{(N-\nu)}^{\nu}(\Lambda_{N}-1)}+\Lambda_{N}^{m}\norm{E}_{\infty}.$

######  Corollary 4.4.

Recall the notation of Lemma 4.1. If $g$ is analytic in $x_{k}$ with an
analytic continuation inside the Bernstein ellipse
$(\rho_{k}\mathrm{e}^{i\theta}+{\rho_{k}}^{-1}\mathrm{e}^{-i\theta})/2$
bounded by $C_{k}$, then

$\displaystyle\norm{g-\widehat{g}}_{\infty}\leq{}$
$\displaystyle\frac{4C_{1}}{(\rho_{1}-1)\rho_{1}^{N_{1}}}+\frac{4\Lambda_{N_{1}}C_{2}}{(\rho_{2}-1)\rho_{2}^{N_{2}}}+\dots+\frac{4\Lambda_{N_{1}}\cdots\Lambda_{N_{m-1}}C_{m}}{(\rho_{m}-1)\rho_{m}^{N_{m}}}$
(46) $\displaystyle+\Lambda_{N_{1}}\cdots\Lambda_{N_{m}}\norm{E}_{\infty},$

where $\Lambda_{N_{k}}\leq\frac{2}{\pi}\log(N_{k}+1)+1$. Especially, if
$N_{1}=\dots=N_{m}=N$, $\rho_{1}=\dots=\rho=\rho$, and $C_{1}=\dots=C_{m}=C$,
then

(47) $\displaystyle\norm{g-\widehat{g}}_{\infty}\leq{}$
$\displaystyle\frac{4(\Lambda_{N}^{m}-1)C}{(\rho-1)\rho^{N}(\Lambda_{N}-1)}+\Lambda_{N}^{m}\norm{E}_{\infty}.$

Note that these bounds for tensorized Chebyshev approximation are slightly
better than those obtained by Mason [45, (21) and (22)]. In particular,
tensorized Chebyshev approximation yields quasi-optimal approximations of
continuous functions by bounded-degree polynomials [45].

The approximation $\widehat{g}$ in Lemma 4.1 is defined as applying the
multilinear map $\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m}$ to the
tensor of function evaluations $G$. Unfortunately, this suffers from the curse
of dimensionality because it requires access to all $N_{1}\cdots N_{m}$
entries of $G$. To mitigate it, a common assumption is that the function
evaluations $G$ can be well approximated by a tensor
$\widehat{G}\in\mathbb{R}^{N_{1}\times\dots\times
N_{m}}\sim\mathbb{R}^{N_{1}}\otimes\dots\otimes\mathbb{R}^{N_{m}}$ that admits
a data-sparse tensor decomposition.

Consider, for example, a _Tucker decomposition_ [63]
$\widehat{G}=(U^{(1)}\otimes\dots\otimes U^{(m)})C$, where
$C\in\mathbb{R}^{r_{1}}\otimes\dots\otimes\mathbb{R}^{r_{m}}$ and each
$U^{(k)}\colon\mathbb{R}^{r_{k}}\to\mathbb{R}^{N_{k}}$, $k=1$, …, $m$, is a
linear map, or, in coordinates, an $N_{k}\times r_{k}$ matrix. Discrete
tensors containing smooth function evaluations can be well approximated with a
Tucker decomposition with small $r_{k}$’s [54]. Given such a Tucker
decomposition of $G$, we can construct $\widehat{g}$ as

(48) $\displaystyle\widehat{g}(x)$
$\displaystyle=(\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m})((U^{(1)}\otimes\dots\otimes
U^{(m)})C)(x)$ (49)
$\displaystyle=\sum_{a_{1}=1}^{r_{1}}\cdots\sum_{a_{m}=1}^{r_{m}}C_{a_{1}\dots
a_{m}}\mathcal{F}_{1}(U^{(1)}_{:a_{1}})(x_{1})\cdots\mathcal{F}_{m}(U^{(m)}_{:a_{m}})(x_{m}),$

where $U^{(k)}_{:a}\in\mathbb{R}^{N_{k}}$ is the $a$th column of
$U^{(k)}\in\mathbb{R}^{N_{k}\times r_{k}}$. Each factor is approximated
independently using the chosen real univariate approximation method
$\mathcal{F}_{k}$.

Since $\mathcal{F}_{1}\otimes\dots\otimes\mathcal{F}_{m}$ is linear, we can
also efficiently construct $\widehat{g}$ whenever $\widehat{G}$ is a linear
combination of Tucker decompositions. This is called a block term
decomposition [25]. We can also impose additional structure on $C$, such as a
tree tensor network or _hierarchical Tucker decomposition_ [34]. Linear
combinations of such hierarchical Tucker decompositions are a special case of
structured block term decomposition [28]. Without going into further details,
we highlight that efficient construction of $\widehat{g}$ is possible for any
structured block term decomposition.

Evidently, the curse of dimensionality is only truly circumvented if a data-
sparse model can also be constructed from a _sparse_ sample of the tensor $G$,
i.e., from viewing only a few entries of $G$. In theory, this is possible for
any polynomial-based decomposition model with a number of samples that equals
the dimensionality of the model; see [16, 52] for precise statements. In
practice, compressed sensing approximation schemes have been proposed for
several tensor decomposition models; see, among others, [48, 8, 35, 3, 2, 56,
42, 55, 66].

The foregoing multivariate approximation scheme is naturally extended to
approximate vector-valued $g$ by applying it component-wise. We have chosen to
present Algorithm 1 with a basic orthogonal Tucker decomposition [63, 26] as
data-sparse model, computed with the sequentially truncated higher-order
singular value decomposition (ST-HOSVD) approximation algorithm [65] from the
full tensor.

##  5 Numerical experiments

To illustrate Algorithm 1 and to verify that the approximant it produces
satisfies the error bounds predicted by Theorem 3.1, we present some examples
inspired by applications in linear algebra. All experiments were run on an HP
EliteBook x360 830 G7 Notebook with 32 GiB memory and an Intel Core i7-10610U
CPU using Julia 1.9.3 and ManiFactor.jl 0.1.0 on Ubuntu 22.04.3 LTS.

###  5.1 Implementation details and Julia package

Algorithm 1 has been implemented in our Julia package, ManiFactor.jl, which is
available online along with the experiments from this section [40]. Our
package depends on the package Manifolds.jl by Axen, Baran, Bergmann, and
Rzecki [7]. Many of the exponential and logarithmic maps for the manifolds in
Appendix A are implemented there. Our package also depends on the package
ApproxFun.jl by Olver and Townsend [47], which is a convenient interface to
Chebyshev interpolation.

As Lemma 4.1 implies that any suitable univariate scheme $\mathcal{I}$ can be
tensorized, ManiFactor.jl can use any such $\mathcal{I}$. Interesting
alternatives to Chebyshev interpolation include Fourier-based schemes and
rational approximation schemes.

Similarly, ManiFactor.jl can use any tensor decomposition. Specifically, it
use incomplete decompositions such as _TT-cross_ , see Strössner, Sun, and
Kressner [58]. The default is an ST-HOSVD from TensorToolbox.jl by Lana Periša
[49].

###  5.2 An example from Krylov subspaces

In this example we approximate a map into the _Grassmannian manifold_
$\mathrm{Gr}(n,k)$ of $k$-dimensional linear subspaces of $\mathbb{R}^{n}$.
See Section A.2 for a discussion of its geometry. The example is inspired by
Zimmermann’s discussion [69, section 4.5] of Grassmannians in the context of
model order reduction, and Benner, Gugercin, and Willcox’s discussion [10] in
the context of Krylov subspace methods.

Consider the heat equation [5, section 9.8] on a finite interval with
endpoints held at a constant temperature $0$:

(50)
$\displaystyle{\partialderivative[2]{u}{x}}(x,t)={}a^{2}{\partialderivative{u}{t}}(x,t),\quad
u(0,t)={}0,\quad u(\pi,t)={}0,$

where $a$ is a diffusion constant. For $t>0$, $u$ is determined by the
_kernel_

(51) $\displaystyle
K(x,x^{\prime},t)=\frac{2}{\pi}\sum_{l=1}^{\infty}\sin(lx)\sin(lx^{\prime})\exp(-ta^{2}l^{2}).$

That is,

(52) $\displaystyle
u(x,t)=\int_{0}^{\pi}K(x,x^{\prime},t)u(x^{\prime},0)\differential{x^{\prime}}$

is an explicit expression for $u$ at time $t$ in terms of $u$ at time $0$. By
discretizing $K$ along $x$ and $x^{\prime}$, we may approximate Eq. 52 as a
matrix equation

(53) $\displaystyle A(t)y=b(t),$

where $A(t)$ is a $n\times n$ symmetric matrix with elements
$A_{ij}(t)=\frac{\pi}{n-1}K(i\frac{\pi}{n-1},j\frac{\pi}{n-1},t)$, $y$ is an
$n$-vector with elements $y_{i}=u(i\frac{\pi}{n-1},0)$, and $b(t)$ is an
$n$-vector with elements $b_{i}(t)=u(i\frac{\pi}{n-1},t)$. Given this setup,
consider the problem of finding a $y$ that solves Eq. 53 for a given $b$.

Krylov subspace methods solve high-dimensional linear equations $Ay=b$ by
looking for a solution in a _Krylov subspace_

(54)
$\displaystyle\mathcal{K}_{k}(A,v)=\operatorname{span}\set{v,Av,\dots,A^{k-1}v}.$

If we fix $k$ and $v\in\mathbb{R}^{n}$, then

(55) $\displaystyle f(t)=\mathcal{K}_{k}(A(t),v)$

is a map from $\mathbb{R}^{+}$ into $\mathrm{Gr}(n,k)$.

Rather than compute the Krylov subspace from its definition for each $t$, we
want to be able to quickly compute an approximate Krylov subspace. Because the
dimension, $k(n-k)$, of $\mathrm{Gr}(n,k)$, is much smaller for small $k$,
than the dimension, $n^{2}$, of the space that $A$ lives in, such an
approximation can be very profitable.

2610141822$undef$${10}^{-5}$${10}^{-10}$${10}^{-15}$$N$$(3+2\sqrt{2})^{-N}$Error
bound Eq. 8Measured error Figure 4: Approximation error for the map Eq. 55
into the Grassmannian $\mathrm{Gr}(n,k)$, compared against what is predicted
by Theorem 3.1, for different numbers of Chebyshev nodes $N$.

In this example, we approximate $f$ on $R=[1,2]$ with $a=1/2$, $n=100$, $k=3$,
and $v$ sampled from a uniform distribution on $[-1,1]^{n}$. These parameters
are chosen such that, for $t\in R$, and for a normalized $b(t)$ in the image
of $A(t)$, it is possible to find a $y$ with forward error on the order of
$10^{-3}$ using $\mathcal{K}_{k}(A(t),v)$. From there, it is straightforward
to extend the Krylov subspace further by Arnoldi iteration.

We now proceed as in Algorithm 1. Choose $p$ and let
$g={\operatorname{log}_{p}}\circ f$. Then consider the polynomial
$\widehat{g}$ that interpolates $g$ in $N+1$ Chebyshev nodes on $R$. Finally,
let $\widehat{f}={\operatorname{exp}_{p}}\circ{\widehat{g}}$.

Since both $f$ and the manifold logarithm are analytic, and because of
Corollary 4.4, we expect the error $\norm{g-\widehat{g}}_{p}$ to decay at
least exponentially in $N$. The actual expressions, however, seem too
complicated to analyze analytically, but since $f$ has a singularity at $t=0$
(the series Eq. 51 fails to converge) we guess that this is the first
singularity that you hit as you analytically extend $g$. In that case,
$\norm{g-\widehat{g}}_{p}$ decays as $(3+2\sqrt{2})^{-N}$, and since
$\mathrm{Gr}(n,k)$ has nonnegative curvature, Theorem 3.1 implies
$d(f,\widehat{f})\leq\norm{g-\widehat{g}}_{p}$.

Indeed, the numerical experiment presented in Fig. 4 supports this guess. In
this experiment, we restrict $f$ to a discrete subset $R^{\prime}\subset R$ of
$1000$ points, sampled independently from a uniform distribution on $R$. We
are thus able to use $\epsilon=\max_{x\in
R^{\prime}}\norm{g(x)-\widehat{g}(x)}_{p}$ in Eq. 8 to derive a bound for
$d(f(x),\widehat{f}(x))$. Figure 4 shows that this error bound is confirmed by
the experiment. After $N=20$, rounding causes the errors to plateau and the
error bound stops being meaningful. The full code to produce Fig. 4 is
available in the ManiFactor.jl package [40, Example3.jl]

###  5.3 An example from dynamic low-rank approximation

Expanding on a synthetic example by Ceruti and Christian [17], we consider

(56) $\displaystyle
A(x_{1},x_{2},x_{3})=\mathrm{e}^{x_{1}}\mathrm{e}^{x_{2}W_{1}}D{\left(\mathrm{e}^{x_{3}W_{2}}\right)}^{\\!\mathsf{T}},$

where $D$ is an $n\times n$ diagonal matrix with diagonal entries $1$,
$2^{-1}$, …, $2^{1-n}$, and $W_{1}$ and $W_{2}$ are skew-symmetric matrices,
with entries sampled independently from a uniform distribution on $[-1,1]$,
and then normalized. Recall that when $W$ is skew-symmetric, $\mathrm{e}^{W}$
is orthogonal. The best rank-$1$ approximation to $A$ is

(57) $\displaystyle
f(x)=\mathrm{e}^{x_{1}}\mathrm{e}^{x_{2}W_{1}}e_{1}{\left(\mathrm{e}^{x_{3}W_{2}}e_{1}\right)}^{\\!\mathsf{T}},$

where $e_{1}$ is the column vector $(1,0,\dots,0)$. Thus $f$ is a map from
$\mathbb{R}^{3}$ into the _Segre manifold_ of rank-$1$ matrices
$\operatorname{Seg}(\mathbb{R}^{n}\times\mathbb{R}^{n})$. See Section A.3 for
a discussion of its geometry. We approximate $f$ on $R=[-1,1]^{3}$ with
$n=100$.

Similarly to Section 5.2, we proceed as Algorithm 1, while restricting to a
discrete subset $R^{\prime}\subset R$ of $1000$ points sampled independently
from a uniform distribution on $R$. Again, because of Corollary 4.4, we expect
the error to decay exponentially in the number of Chebyshev nodes
$N=N_{1}=N_{2}=N_{3}$. By Eq. 67, we have a lower bound $H=-\mathrm{e}^{2}$
for the sectional curvature. Thus Theorem 3.1 implies that the approximation
error is bounded by Eq. 9.

The numerical experiment presented in Fig. 5 again confirms the error bound.
After $N=13$, rounding causes the errors to plateau and the error bound stops
being meaningful. The full code to product Fig. 5 is available in the
ManiFactor.jl package [40, Example4.jl].

2345678910111213141516$undef$${10}^{-5}$${10}^{-10}$${10}^{-15}$$N$Error bound
Eq. 9Measured error Figure 5: Approximation error for the map Eq. 57 into the
Segre manifold $\operatorname{Seg}(\mathbb{R}^{n}\times\mathbb{R}^{n})$,
compared against what is predicted by Theorem 3.1, for different numbers of
Chebyshev nodes $N$ in each variable.

###  5.4 Effects of truncating the tensor decomposition

In Section 4.2, we discussed how Algorithm 1 is implemented using a
_sequentially truncated_ HOSVD. We argued that working with sparse
representations is key to circumventing the curse of dimensionality. We thus
now illustrate how different levels of truncation affect the error. Consider
again the example from Section 5.3. Figure 6 shows the same numerical
experiment, but with the decomposition of $G$ truncated at different
multilinear ranks (mrank).

Without truncation, $G$ has $(2(n-1)+1)(N+1)^{3}$ entries since
$\operatorname{Seg}(\mathbb{R}^{n}\times\mathbb{R}^{n})$ has dimension
$2(n-1)+1$. Compare this to an ST-HOSVD $\widehat{G}$ of $G$ with multilinear
rank $(r_{0},r_{1},r_{2},r_{3})$, which only has
$2(n-1)r_{0}+(N+1)r_{1}+(N+1)r_{2}+(N+1)r_{3}+r_{0}r_{1}r_{2}r_{3}$ entries.
In Fig. 6, we see that multilinear rank $(15,2,8,8)$ is enough to reach
machine precision at $N=14$. In this case, $G$ has $668\,250$ entries, while
$\widehat{G}$ only has $5160$ entries. The multilinear ranks are chosen
automatically by TensorToolbox.jl to satisfy the error thresholds ${10}^{-6}$,
${10}^{-8}$, ${10}^{-10}$, and ${10}^{-13}$, respectively.

2345678910111213141516$undef$${10}^{-5}$${10}^{-10}$${10}^{-15}$$N$mrank
$(9,2,5,5)$mrank $(11,2,6,6)$mrank $(13,2,7,7)$mrank $(15,2,8,8)$ Figure 6:
Approximation error for a map into the Segre manifold
$\operatorname{Seg}(\mathbb{R}^{n}\times\mathbb{R}^{n})$, with the tensor of
samples is truncated at different multilinear ranks.

##  6 Conclusions

In conclusion, we find that the approximability of maps into a Riemannian
manifold is related to the curvature of the manifold. Moreover, for many
manifolds that are relevant for numerical analysis, especially matrix
manifolds, explicit curvature bounds allow us to derive explicit error bounds.

From here, several areas of further research are visible. First, one could
derive the geodesics and sectional curvature for more manifolds. For example
spaces of tensors with fixed rank or multilinear rank. Second, one could
develop more efficient and more optimal ways of choosing the point $p$ from
where the manifold is linearized. The current way, geodesic interpolation of
sample evaluations, is a heuristic to estimate the Riemannian center of mass
of those evaluations, which in turn is a heuristic to minimize $\sigma$ in
Theorem 3.1. Third, it could be possible to consider more than one tangent
space. Given several approximations of $f$, based around several points, what
is a good estimate of $f(x)$, and what are good ways of choosing such points?

##  Appendix A Manifolds with lower bounded sectional curvature

In this section, we list Riemannian manifolds with explicit lower bounds for
their sectional curvature $K$ and whose manifold exponential and logarithm are
known. For these manifolds, Theorem 3.1 gives an explicit bound for the error
on the manifold in terms of the error on the tangent space. While sectional
curvature is always bounded locally, for many manifolds that are relevant for
applications we can find explicit expressions for those local bounds. In many
cases, the sectional curvature is even bounded globally, so that we can choose
the same $H$ for every chart in Theorem 3.1.

###  A.1 Lie groups

Given an inner product on its Lie algebra, one can extend it to a unique
Riemannian metric on the Lie group by demanding invariance under left (or
right) translation. Because of this invariance, the set of sectional
curvatures at a point is independent of that point. Therefore, a lower bound
for the sectional curvature at one point will be a global lower bound. There
is an expression for the sectional curvature of a Lie group in terms of its
Lie bracket [21, Proposition 3.18 (3)]. We only note that, for matrix groups
whose Lie algebra has the Euclidean inner product, it follows from this
expression that111Using $\norm{[u,v]}\leq\sqrt{2}\norm{u}\norm{v}$, which also
implies $\norm{\operatorname{ad}_{u}}_{\mathrm{op}}\leq\sqrt{2}\norm{u}$.

(58) $\displaystyle-\frac{11}{2}\leq K.$

For example, $\mathrm{GL}(n)$ is the group of invertible $n\times n$ matrices.
Its metric is induced by the inner product,
$\left<{u},{v}\right>_{A}=\operatorname{tr}\left({\left({A}^{-1}u\right)}^{\\!\mathsf{T}}{A}^{-1}v\right)$,
where $u$, $v\in T_{A}\mathrm{GL}(n)=\mathbb{R}^{n\times n}$. Locally, all
(finite-dimensional) Lie groups are Lie subgroups of some $\mathrm{GL}(n)$.
Andruchow, Larotonda, Recht, and Varela [4] derive its geodesics.

The compact Lie groups and their products with vector spaces are precisely the
Lie groups that admit a _bi-invariant_ metric, i.e., a metric that is
invariant under both left and right translation. This is also equivalent to
the manifold exponential agreeing with the Lie group exponential. The
sectional curvature is positive for such manifolds [21, Corollary 3.19 and
Proposition 3.34].

The classification of compact Lie groups is a classic result. The real compact
Lie groups are [6, Theorem 2.17]

  * •

tori $S^{1}\times\dots\times S^{1}$,

  * •

orthogonal groups $\mathrm{O}(n)$ and special orthogonal groups
$\mathrm{SO}(n)$,

  * •

unitary groups $\mathrm{U}(n)$ and special unitary groups $\mathrm{SU}(n)$,

  * •

spin groups $\mathrm{Spin}(n)$,

  * •

symplectic groups $\mathrm{Sp}(n)$,

  * •

compact exceptional Lie groups $G_{2}$, $F_{4}$, $E_{6}$, and $E_{8}$,

  * •

products, disjoint unions, and finite covers of the above.

###  A.2 Homogeneous spaces

A _Riemannian homogeneous space_ is a Riemannian manifold equipped with some
transitive action that the metric is invariant under. Similarly to Lie groups,
the curvature is then the same at every point.

Any Riemannian homogeneous space $X$ is of the form $G/H$ where $G$ is the Lie
group that acts on $X$ and $H$ is a Lie subgroup. The _Stiefel manifold_

(59) $\displaystyle\mathrm{St}(n,k)=\mathrm{O}(n)\mathbin{/}\mathrm{O}(n-k)$

is the space of orthogonal $n\times k$ matrices. Zimmermann [67] considers
efficient geodesic computations on the Stiefel manifold. It has nonnegative
curvature.

The _Grassmannian manifold_

(60)
$\displaystyle\mathrm{Gr}(n,k)=\mathrm{O}(n)\mathbin{/}(\mathrm{O}(n-k)\times\mathrm{O}(k))$

is the space of $k$-dimensional linear subspaces of $\mathbb{R}^{n}$.
Bendokat, Zimmermann, and Absil [9] consider efficient geodesic computations
on the Grassmannian manifold. It also has nonnegative curvature.

For the Stiefel and Grassmannian manifolds, their metrics as Riemannian
homogeneous spaces are called _canonical metrics_ to distinguish them from
their metrics as embedded submanifolds of $\mathbb{R}^{n\times n}$.

$\sigma_{B}\colon A\mapsto BA{B}^{\\!\mathsf{T}}$ defines a transitive action
of $\mathrm{GL}(n)$ on the space $\mathcal{P}_{n}$ of _symmetric positive
definite $n\times n$ matrices_. The _affine invariant metric_ on
$\mathcal{P}_{n}$ is thus derived from the quotient

(61) $\displaystyle\mathcal{P}_{n}=\mathrm{GL}(n)\mathbin{/}\mathrm{O}(n).$

Dolcetti and Pertici [29, Proposition 3.4] compute its geodesics, while
Criscitiello and Boumal [24, Proposition I.1] show that its sectional
curvature satisfies

(62) $\displaystyle-\frac{1}{2}\leq K\leq 0.$

Similarly to $\mathcal{P}_{n}$, $\sigma_{B}$ defines a transitive action of
$\mathrm{GL}^{+}(n)$, the group of matrices with positive determinant, on the
space $\mathrm{S}_{+}(n,k)$ of _rank- $k$ positive semidefinite $n\times n$
matrices_. We have that

(63)
$\displaystyle\mathrm{S}_{+}(n,k)=\mathrm{GL}^{+}(n)\mathbin{/}\left(\left(\mathrm{SO}(k)\ltimes\mathbb{R}^{k\times(n-k)}\right)\times\mathrm{GL}^{+}(n-k)\right),$

where $\ltimes$ is the semidirect product. Vandereycken, Absil, and Vandewalle
[64] use a right-invariant metric on $\mathrm{GL}^{+}(n)$ to derive geodesics
on $\mathrm{S}_{+}(n,k)$.

We also mention another metric that can be put on the space of rank-$k$
positive semidefinite matrices. Let $\mathbb{R}^{n\times k}_{*}$ be the space
of full rank $n\times k$ matrices. Any rank-$k$ positive semidefinite matrix
can be written as $Y{Y}^{\\!\mathsf{T}}$, $Y\in\mathbb{R}^{n\times k}_{*}$.
Such a representation is unique up to right multiplication by an orthogonal
matrix. Hence, we define

(64) $\displaystyle\mathcal{S}_{+}(n,k)=\mathbb{R}^{n\times
k}_{*}\mathbin{/}\mathrm{O}(k),$

where $\mathbb{R}^{n\times k}_{*}$ is equipped with the Euclidean metric.
Although this manifold is not a Riemannian homogeneous space, a metric may
still be induced by demanding that the quotient map is a Riemannian
submersion. Massart and Absil [46] summarize its geometry. It has nonnegative
curvature.

###  A.3 Segre manifold

The _Segre manifold_ ,
$\mathrm{Seg}(\mathbb{R}^{N_{1}}\times\dots\times\mathbb{R}^{N_{m}})$, is the
space of rank-$1$ tensors in the product space
$\mathbb{R}^{N_{1}}\otimes\dots\otimes\mathbb{R}^{N_{m}}$ [37]. It is immersed
by

(65) $\displaystyle\operatorname{Seg}\colon\mathbb{R}^{+}\times
S^{N_{1}-1}\times\dots\times S^{N_{m}-1}$
$\displaystyle\to\mathbb{R}^{N_{1}}\otimes\dots\otimes\mathbb{R}^{N_{m}},$
(66) $\displaystyle(\lambda,x_{1},\dots,x_{m})$ $\displaystyle\mapsto\lambda
x_{1}\otimes\dots\otimes x_{m}.$

Swijsen, van der Veken, and Vannieuwenhoven [60] derives the exponential map
for this manifold in the metric induced by the immersion. Moreover, Swijsen
[59, Corollary 6.2.4] derive the manifold logarithm and shows that the
sectional curvature satisfies

(67) $\displaystyle-\frac{1}{\lambda^{2}}\leq K\leq 0.$

Note that in this case there is only a local lower bound, namely if the image
of $f$ fits in a chart with $\lambda$ lower bounded by $\lambda_{*}$, then we
may choose $H=-1/\lambda_{*}^{2}$.

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, vol. 78, Princeton University Press, 12 2008.
  * [2] E. Acar, D. M. Dunlavy, and T. G. Kolda, A scalable optimization approach for fitting canonical tensor decompositions, Journal of Chemometrics, 25 (2011), pp. 67–86.
  * [3] E. Acar, D. M. Dunlavy, T. G. Kolda, and M. Mørup, Scalable tensor factorizations for incomplete data, Chemometrics and Intelligent Laboratory Systems, 106 (2011), pp. 41–56.
  * [4] E. Andruchow, G. Larotonda, L. Recht, and A. Varela, The left invariant metric in the general linear group, Journal of Geometry and Physics, 86 (2014), pp. 241–257.
  * [5] G. B. Arfken and H.-J. Weber, Mathematical Methods for Physicists, Elsevier Academic, 6 ed., 2005.
  * [6] A. Arvanitogeorgos, An Introduction to Lie Groups and the Geometry of Homogeneous Spaces, American Mathematical Society, 2003.
  * [7] S. D. Axen, M. Baran, R. Bergmann, and K. Rzecki, Manifolds.jl: An extensible Julia framework for data analysis on manifolds, 2023.
  * [8] J. Ballani, L. Grasedyck, and M. Kluge, Black box approximation of tensors in hierarchical Tucker format, LAA, 438 (2013), pp. 639–657.
  * [9] T. Bendokat, R. Zimmermann, and P. A. Absil, A Grassmann manifold handbook: Basic geometry and computational aspects, 2023, 2011.13699.  arxiv.org.
  * [10] P. Benner, S. Gugercin, and K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review, 57 (2015), pp. 483–531, https://doi.org/10.1137/130932715.
  * [11] R. Bergmann and P.-Y. Gousenbourger, A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve, Frontiers in Applied Mathematics and Statistics, 4 (2018).
  * [12] G. Beylkin and M. J. Mohlenkamp, Numerical operator calculus in higher dimensions, Proceedings of the National Academy of Sciences of the United States of America, 99 (2002), pp. 10246–10251.
  * [13] , Algorithms for numerical analysis in high dimensions, SIAM Journal on Scientific Computing, 26 (2005), pp. 2133–2159.
  * [14] D. Bigoni, A. P. Engsig-Karup, and Y. M. Marzouk, Spectral tensor-train decomposition, SIAM Journal on Scientific Computing, 38 (2016), pp. A2405–A2439.
  * [15] N. Boumal, An Introduction to Optimization on Smooth Manifolds, Cambridge University Press, Mar. 2023.
  * [16] P. Breiding, F. Gesmundo, M. Michałek, and N. Vannieuwenhoven, Algebraic compressed sensing, Applied and Computational Harmonic Analysis, 65 (2023), pp. 374–406.
  * [17] G. Ceruti and L. Christian, An unconventional robust integrator for dynamical low-rank approximation, BIT Numerical Mathematics, 62 (2022).
  * [18] R. Chakraborty and B. C. Vemuri, Recursive Fréchet mean computation on the Grassmannian and its applications to computer vision, in 2015 IEEE International Conference on Computer Vision (ICCV), 2015, pp. 4229–4237.
  * [19] , Statistics on the Stiefel manifold: Theory and applications, The Annals of Statistics, 47 (2019), pp. 415 – 438.
  * [20] , Efficient recursive estimation of the Riemannian barycenter on the hypersphere and the special orthogonal group with applications, in Riemannian Geometric Statistics in Medical Image Analysis, X. Pennec, S. Sommer, and T. Fletcher, eds., Academic Press, 2020, pp. 273–297.
  * [21] J. Cheeger and D. G. Ebin, Comparison Theorems in Riemannian Geometry, AMS Chelsea Publishing, 2008.
  * [22] G. Cheng, J. Ho, H. Salehian, and B. C. Vemuri, Recursive Computation of the Fréchet Mean on Non-positively Curved Riemannian Manifolds with Applications, Springer International Publishing, Cham, 2016, pp. 21–43.
  * [23] O. Christensen and K. L. Christensen, Approximation Theory: from Taylor Polynomials to Wavelets, Applied and Numerical Harmonic Analysis, Birkhäuser, 2004.
  * [24] C. Criscitiello and N. Boumal, An accelerated first-order method for non-convex optimization on manifolds, 2021, 2008.02252.  arxiv.org.
  * [25] L. De Lathauwer, Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 1033–1066.
  * [26] L. De Lathauwer, B. De Moor, and J. Vandewalle, A multilinear singular value decomposition, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1253–1278.
  * [27] T. De Weer, N. Vannieuwenhoven, N. Lammens, and K. Meerbergen, The parametrized superelement approach for lattice joint modelling and simulation, Computational Mechanics, 70 (2022), pp. 451–475.
  * [28] N. Dewaele, P. Breiding, and N. Vannieuwenhoven, The condition number of many tensor decompositions is invariant under tucker compression, Numerical Algorithms, 94 (2023), pp. 1003–1029.
  * [29] A. Dolcetti and D. Pertici, Differential properties of spaces of symmetric real matrices, 2018, 1807.01113.  arxiv.org.
  * [30] S. Dolgov, D. Kressner, and C. Strössner, Functional tucker approximation using chebyshev interpolation, SIAM Journal on Scientific Computing, 43 (2021), pp. A2190–A2210, https://doi.org/10.1137/20M1356944.
  * [31] K. Glau, D. Kressner, and F. Statti, Low-rank tensor approximation for Chebyshev interpolation in parametric option pricing, SIAM Journal on Financial Mathematics, 11 (2020), pp. 897–927.
  * [32] A. Gorodetsky, S. Karaman, and Y. Marzouk, A continuous analogue of the tensor-train decomposition, Computer Methods in Applied Mechanics and Engineering, 347 (2019), pp. 59–84.
  * [33] P.-Y. Gousenbourger, E. Massart, and P.-A. Absil, Data fitting on manifolds with composite Bézier-like curves and blended cubic splines, Journal of Mathematical Imaging and Vision, 61 (2019), pp. 645–671.
  * [34] L. Grasedyck, Hierarchical singular value decomposition of tensors, SIAM Journal on Matrix Analysis and Applications, 31 (2010), pp. 2029–2054.
  * [35] L. Grasedyck, M. Kluge, and S. Krämer, Variants of alternating least squares tensor completion in the tensor train format, SIAM Journal on Scientific Computing, 37 (2015), pp. A2424–A2450.
  * [36] W. H. Greub, Multilinear algebra, Die Grundlehren der mathematischen Wissenschaften: 136, Springer, 2 ed., 1978.
  * [37] J. Harris, Algebraic Geometry, A First Course, vol. 133 of Graduate Text in Mathematics, Springer-Verlag, 1992.
  * [38] B. Hashemi and Y. Nakatsukasa, On the spectral problem for trivariate functions, BIT Numerical Mathematics, 58 (2018), pp. 981–1008.
  * [39] B. Hashemi and L. N. Trefethen, Chebfun in three dimensions, SIAM Journal on Scientific Computing, 39 (2017), pp. C341–C363.
  * [40] S. Jacobsson, ManiFactor.jl: Approximating maps into manifolds.  https://gitlab.kuleuven.be/numa/software/ManiFactor, 2024.
  * [41] H. Karcher, Riemannian center of mass and mollifier smoothing, Communications on Pure and Applied Mathematics, 30 (1977), pp. 509–541.
  * [42] D. Kressner, M. Steinlechner, and B. Vandereycken, Low-rank tensor completion by Riemannian optimization, BIT Numerical Mathematics, 54 (2014), pp. 447–468.
  * [43] J. M. Lee, Introduction to Smooth Manifolds, Graduate texts in mathematics: 218, Springer, 2013.
  * [44] , Introduction to Riemannian Manifolds., Graduate texts in mathematics: 176, Springer, 2018.
  * [45] J. Mason, Near-best multivariate approximation by Fourier series, Chebyshev series and Chebyshev interpolation, Journal of Approximation Theory, 28 (1980), pp. 349–358.
  * [46] E. Massart and P.-A. Absil, Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 171–198, https://doi.org/10.1137/18M1231389.
  * [47] S. Olver and A. Townsend, A practical framework for infinite-dimensional linear algebra, in Proceedings of the 1st Workshop for High Performance Technical Computing in Dynamic Languages – HPTCDL ‘14, IEEE, 2014.
  * [48] I. V. Oseledets and E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra and its Applications, 432 (2010), pp. 70–88.
  * [49] L. Periša, Tensortoolbox.jl: Julia package for tensors.  https://github.com/lanaperisa/TensorToolbox.jl, 2023.
  * [50] M. Reid and B. Szendroi, Geometry and Topology, Cambridge University Press, 2005.
  * [51] J. R. Rice, A theory of condition, SIAM Journal on Numerical Analysis, 3 (1966), pp. 287–310, https://doi.org/10.1137/0703023.
  * [52] Y. Rong, Y. Wang, and Z. Xu, Almost everywhere injectivity conditions for the matrix recovery problem, Applied and Computational Harmonic Analysis, 50 (2021), pp. 386–400.
  * [53] M. H. Schultz, $L^{\infty}$-multivariate approximation theory, SIAM Journal on Numerical Analysis, 6 (1969), pp. 161–183.
  * [54] T. Shi and A. Townsend, On the compressibility of tensors, SIAM Journal on Matrix Analysis and Applications, 42 (2021), pp. 275–298.
  * [55] M. Signoretto, Q. Tran Dinh, L. De Lathauwer, and J. A. K. Suykens, Learning with tensors: a framework based on convex optimization and spectral regularization, Machine Learning, 94 (2013), pp. 303–351.
  * [56] M. Steinlechner, Riemannian optimization for high-dimensional tensor completion, SIAM Journal on Scientific Computing, 38 (2016), pp. S461–S484.
  * [57] C. Strössner and D. Kressner, Low-rank tensor approximations for solving multimarginal optimal transport problems, SIAM Journal on Imaging Sciences, 16 (2023), pp. 169–191, https://doi.org/10.1137/22M1478355.
  * [58] C. Strössner, B. Sun, and D. Kressner, Approximation in the extended functional tensor train format, 2022, 2211.11338.  arxiv.org.
  * [59] L. Swijsen, Tensor decompositions and Riemannian optimization, PhD thesis, KU Leuven, 2022.
  * [60] L. Swijsen, J. van der Veken, and N. Vannieuwenhoven, Tensor completion using geodesics on Segre manifolds, Numerical Linear Algebra with Applications, 29 (2022).
  * [61] V. Temlyakov, Multivariate Approximation, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2018.
  * [62] L. N. Trefethen, Approximation Theory and Approximation Practice, Extended Edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2019, https://epubs.siam.org/doi/pdf/10.1137/1.9781611975949.
  * [63] L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31 (1966), pp. 279–311.
  * [64] B. Vandereycken, P.-A. Absil, and S. Vandewalle, A Riemannian geometry with complete geodesics for the set of positive semidefinite matrices of fixed rank, IMA Journal of Numerical Analysis, 33 (2012), pp. 481–514, https://academic.oup.com/imajna/article-pdf/33/2/481/1845277/drs006.pdf.
  * [65] N. Vannieuwenhoven, R. Vandebril, and K. Meerbergen, A new truncation strategy for the higher-order singular value decomposition, SIAM Journal on Scientific Computing, 34 (2012), pp. A1027–A1052, https://doi.org/10.1137/110836067.
  * [66] N. Vervliet, O. Debals, L. Sorber, and L. De Lathauwer, Breaking the curse of dimensionality using decompositions of incomplete tensors: Tensor-based scientific computing in big data analysis, IEEE Signal Processing Magazine, 31 (2014), pp. 71–79.
  * [67] R. Zimmermann, A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canonical metric, SIAM Journal on Matrix Analysis and Applications, 38 (2017), pp. 322–342, https://doi.org/10.1137/16M1074485.
  * [68] , Hermite interpolation and data processing errors on Riemannian matrix manifolds, SIAM Journal on Scientific Computing, 42 (2020), pp. A2593–A2619, https://doi.org/10.1137/19M1282878.
  * [69] , Manifold interpolation and model reduction, 2022, 1902.06502.
  * [70] R. Zimmermann and R. Bergmann, Multivariate Hermite interpolation of manifold-valued data, 2022, 2212.07281.  arxiv.org.


* (23) B. Carter  Proceedings of the Royal Society of London Series A 433, 45 (1991)

  * (24) B Carter and D. Langlois, Phys. Rev.D 51, 5855 (1995), arXiv: 9507058 [hep-th]

  * (25) B. Carter B and I.M. Khalatnikov  Phys. Rev. D 45, 4536 (1992)

  * (26) T.S. Olson and W. A. Hiscock W. A.  Phys. Rev. D 41 , 3687 (1990)

  * (27) L Gavassino  Stability and causality of Carter s multifluid theory, arXiv: 2202.06760 [gr-qc]

  * (28) U. Heinz and R. Snellings,  Collective Flow and Viscosity in Relativistic Heavy-Ion Collisions, Annu. Rev. Nucl. Part. Sci. 63, 123 (2013)

  * (29) L Gavassino and M. Antonelli  Unified Extended Irreversible Thermodynamics and the stability of relativistic theories for dissipation , arXiv: 2202.06760 [gr-qc]

  * (30) F. S. Bemfica, M. M. Disconzi, V, Hoang V, J. Noronha, and M Radosz M  Phys.Rev.Lett., 126, 222301, (2021)

  * (31) Gavassino L Antonelli and M, Haskell B Symmetric-hyperbolic quasi-hydrodynamics arXiv: 2207.14778 [gr-qc]

  * (32) Olson T. S. Annals of Physics 199, 18, (1990)

  * (33) I Müller and W. Weiss The Europ. Phys.Jour. H 37, 139, (2012)

  * (34) Eckart C 1940 Phys. Rev. 58, 919, (1940)

  * (35) L. Landau L and E.M. Lifshitz, 1958 Fluid Mechanics p. 500. Addison Wesley, Reading, Mass

  * (36) W.A. Hiscock and L. Lindblom Phys. Rev. D31, 725, (1985)

  * (37) W.A. Hiscock and L. Lindblom Phys. Rev. D35, 3723, (1985)

  * (38) I.S Liu , I Müller and T Ruggeri, Ann. Phys. 169 (1986)

  * (39) S Pennisi Proceedings of symposium of Kinetic Theory and Extended Thermodynamics, Eds. Müller I and Ruggeri T, Bologna (1987)

  * (40) R Geroch and L Lindblom, Phys.Rev. D41, 1855 (1990)

  * (41) Priou D, Phys. Rev. D 43, 1223 (1991)

  * (42) G.S. Denicol, H. Niemi, E.Molnar and D.H. Rischke  Phys. Rev. D 85, 114047 (2012)

  * (43) R. Baier, P. Romatschke, D. Thanh Son, A. O. Starinets, and M. A. Stephanov,  J. of High Energy Physics,100 (2008), arXiv:0712.2451 [hep-th].

  * (44) W. Florkowski and R. Ryblewski,  Phys. Rev. C 83, 034907 (2011)

  * (45) M. Alqahtani, M. Nopoush, and M. Strickland,  Prog. Part. Nucl. Phys. 101, 204 (2018)

  * (46) P. Kovtun,  J. Phys . A45, 473001, (2012), 1205.5040 [hep-th]

  * (47) P. Kovtun, Journal of High Energy Physics 34 (2019), arXiv:1907.08191 [hep-th]

  * (48) F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. D 98, 104064 (2018)

  * (49) F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. X 12, 021044, (2022)

  * (50) W Israel, in Relativistic Fluid Dynamics eds. A. Anile and Y. Choquet-Bruhat (Springer) (1989)

  * (51) R. Maartens, Lectures given at the Hanno Rund Workshop on Relativity (1996)

  * (52) J. F. Salazar and T. Zannias,International Journal of Modern Physics D 29, 2030010 (2020). For an extended version see: J. F. Salazar and T. Zannias (2021), arXiv:1904.04368v3

  * (53) J. F. Salazar, Progresos recientes en termodinamica irreversible de fluidos relativistas, Ph.D. Thesis (in Spanish), Instituto de Física y Matematicas, Univ. Michoacana Mx, (2021)

  * (54) W. Florkowski, M. P. Heller, and M. Spalins’ki, Reports on Progress in Physics  81, 046001 (2018), arXiv:1707.02282 [hep-ph].

  * (55) D. M. Aitken Relativistic Irreversible Thermodynamics, M. Sc. thesis, University of Alberta, Edmonton, Canada (1964)

  * (56) H. B. Callen Thermodynamics an Introduction to Thermostatistics, 2nd ed., Wiley, New York, (1985)

  * (57) W A Hiscock and L Lindblom Ann. Phys 151, 466 (1983)

  * (58) A L Garcia-Perciante, M E Rubio, O A Reula  Jour. Stat. Phys. 81, 246, (2020) [arXiv:1908.04445 [gr-qc]]. (2019)

  * (59) S. Diles,  Phys. Lett.B 825, 136850, (2022)

  * (60) M. P. Heller, A. Serantes, M. Spalinski, V. Svensson, and B. Withers, Phys. Rev. D 104, 066002, (2021)

  * (61) S Denicol, T Kodama, T Koide, and Ph Mota, Journal of Physics G: Nuclear and Particle Physics, 35, (11),115102, (2008)

  * (62) Shi Pu, Tomoi Koide, and Dirk H. Rischke,Phys. Rev. D, 81, 114039, (2010)

  * (63) S. Garcia-Colin A. Sandoval-Villalbazo, A.L. Garcia-Perciante Physica A 388, 3770, (2009).

  * (64) W A Hiscock W A and Lindblom L, Phys.Lett. A 131, 509 (1988)

  * (65) J Ehlers Continuum Mechanics and Thermodynamics in Relativity, Astrophysics and Cosmology, ed. W. Israel, Dordrecht, Reidel Pub. Com. (1973)

  * (66) H Andreasson, Living Rev. Relativ 5, 2005-2 ( 2005)

  * (67) O Sarbach and T Zannias Class. Quant. Grav. 31, 085013, (2014)

  * (68) R. O. Acuna-Cardenas, C. Gabarrete, O. Sarbach Gen. Rel. Grav., 54, 23, (2022) [arXiv:2106.09235]

  * (69) I.S Liu Arch. Rat.Mech.Anal, 46, (1972)

  * (70) J. F. Salazar and T. Zannias, Some Remarks on Relativistic Fluids of Divergence Type, Report IFM-UMSNH unpublished (2022).

  * (71) F. Becattini, L. Bucciantini, E. Grossi and L. Tinti Eur. Phys. J C 75, 191, (2015)

  * (72) G.S. Denicol, H. Niemi, E. Molnar, D.H. Rischke Phys. Rev D.85, 114047, (2012) [Erratum: Phys. Rev.D91,no.3, 039902 (2015)], 1202.4551.


11institutetext:  Department of Physics,
Indian Institute of Science Education and Research (IISER), Berhampur,
Postal Code-760010, Odisha, India.

# Energy Dependence of Strangeness Production in Heavy-ion Collision

Md Nasim 11<EMAIL_ADDRESS>
###### Abstract

An experimental overview of the energy dependence of strangeness production is
presented. The strange hadrons are considered a good probe to study the QCD
matter created in relativistic nucleus-nucleus collisions. The heavy-ion
experiments at SPS, RHIC, and LHC have recorded a wealth of data in proton-
proton, proton-nucleus and nucleus-nucleus collisions at several beam
energies. In this proceeding, I discuss the invariant yield and azimuthal
anisotropy measurement of strange hadrons in nucleus-nucleus collisions at
SPS, RHIC, and LHC.

##  1 Introduction

The relativistic heavy-ion collision provides a unique opportunity to study
the properties of QCD matter at various temperatures and densities. In order
to map out the phase diagram of the QCD matter, experimental programs started
in the early 1990s at the Brookhaven Alternating Gradient Synchrotron (AGS)
and the CERN Super Proton Synchrotron (SPS) followed by Relativistic Heavy Ion
Collider (RHIC) at Brookhaven and recently at Large Hadron Collider (LHC) at
CERN. These experiments allowed us to study the QCD matter over a large range
of baryonic chemical potential (or net-baryon density) by varying center-of-
mass energy( $\sqrt{s_{NN}}$) of the two colliding nucleus.
Strange hadron production in nucleus-nucleus collisions is believed to be an
important tool to study the dynamics of the produced QCD matter. Due to the
large production cross-section of strange quarks in the QGP, an enhancement in
the production of strange hadrons in QGP compared to non-QGP medium is
expected rafaleski . Due to low energy density, the matter produced in the
nucleon-nucleon (p+p) collisions is considered as a baseline for the non-QGP
medium. Production of strange hadrons is being measured in nucleus-nucleus
(A+A) collisions and compared with that in p+p collisions to look for any
deviation from the non-QGP baseline. Any significant deviation would hint at
the formation of deconfined matter of quarks and gluons in A+A collisions.
Among all strange hadrons, multi-strange hadrons (e.g. $\phi$, $\Omega$) are
expected to have relatively small hadronic interaction cross-sections and
therefore they can carry the information directly from the chemical freeze-out
stage with almost no distortion due to late-stage hadronic interaction
star_white ; shor . The measurement of production yield and azimuthal
anisotropy of multi-strange hadrons are considered as a useful probe to study
the hadronization and collective properties of the QGP medium.

##  2 Strange Hadron Production

###  2.1  Transverse momentum ($p_{T}$) integrated particle ratio

Test of thermal model:
Figure 1 (the left panel) shows $K/\pi$ ratios as a function of
$\sqrt{s_{NN}}$ in central A+A collision measured at AGS, SPS, and RHIC AGS_1
; AGS_2 ; AGS_3 ; SPS_1 ; SPS_2 ; RHIC_1 . The $K/\pi$ ratio is of interest,
as it reflects the strangeness content relative to entropy in heavy-ion
collisions. The $K^{+}/\pi^{+}$ and $K^{-}/\pi^{-}$ ratios are shown by red
and black symbols, respectively. The observed peak position in the energy
dependence of $K^{+}/\pi^{+}$ has been suggested to be a signature of a change
in degrees of freedom kpi_horn_1 ; kpi_horn_2 . The thermal model calculations
are shown as yellow band for $K^{+}/\pi^{+}$ and green band for
$K^{-}/\pi^{-}$ thermal_1 ; thermal_2 ; thermal_3 . It is shown that the
thermal model calculation explains the measured $K/\pi$ ratios fairly well.
Cannonical Ensemble at high-baryon density:
The $\phi/K^{-}$ ratios as a function of $\sqrt{s_{NN}}$ is presented in the
right panel of figure 1 HADES ; RHIC_3gev . The thermal model calculations,
based on the Grand Canonical Ensemble (GCE) and Canonical Ensemble (CE) for
strangeness with several different choices of strangeness correlation length
are shown by different line colors thermal_2 ; thermal_3 . We can see for
$\sqrt{s_{NN}}$ $>$ 5 GeV, both GCE and CE calculations explain the data,
however for $\sqrt{s_{NN}}$ $\leq$ 3.0 GeV GCE calculation under-predicts the
measured $\phi/K^{-}$ ratios. Calculation using Canonical Ensemble (CE) with
strangeness correlation length $\sim$ 3 fm can explain the data measured by
STAR and HADES experiment.
Evidence of hadronic rescattering in central A+A collision:
Figure 2 shows $K^{*0}/K$ and $\phi/K$ ratios as a function of multiplicity in
p+Pb and Pb+Pb collisions at 2.76 and 5.02 TeV measured by the ALICE
experiment LHC_5TeV . We see that $K^{*0}/K$ ratios decrease with increasing
multiplicity. Whereas the $\phi/K$ ratio is almost independent of
multiplicity. This observation can be understood considering the re-scattering
of daughter particles in the hadronic phase. The lifetime of $K^{*0}$
resonance is short ($\sim$ 4 fm/c) and therefore they decay inside the
fireball. The daughter particles of $K^{*0}$ undergo re-scattering in the
hadronic medium and their momentum may get changed. This would lead to a loss
in the reconstruction of the parent resonance. The lifetime of the $\phi$
meson is $\sim$ 42 fm/c. Because of the longer lifetime, the $\phi$ meson will
mostly decay outside the fireball and therefore its daughters will not have
much time to rescatter in the hadronic phase. We can also see from figure 2
that thermal model calculations, which do not include hadronic rescattering,
explain $\phi/K$ ratios and over-predict $K^{*0}/K$ at the central Pb+Pb
collisions. This is a clear piece of evidence of the hadronic rescattering
effect on short-lived resonance in central Pb+Pb collisions.

Figure 1: Left panel: The $K^{+}/\pi^{+}$ and $K^{-}/\pi^{-}$ ratios as a
function of $\sqrt{s_{NN}}$ in central A+A collision AGS_1 ; AGS_2 ; AGS_3 ;
SPS_1 ; SPS_2 ; RHIC_1 . Right panel: The $\phi/K^{-}$ ratios as a function of
$\sqrt{s_{NN}}$ HADES ; RHIC_3gev . Figure 2: The $K^{*0}/K$ and $\phi/K$
ratios as a function of multiplicity in p+Pb and Pb+Pb collisions at 2.76 and
5.02 TeV measured by the ALICE experiment LHC_5TeV .

###  2.2 Particle ratio as a function of $p_{T}$

Probing hadronization through baryon-to-meson ratio:
The baryon-to-meson ratios have been predicted to be sensitive to the parton
dynamics of the collision system. Figure 3 shows the $\bar{\Lambda}/K^{0}_{S}$
ratio as a function of $p_{T}$ for different centralities in Au+Au collisions
at $\sqrt{s_{NN}}$ = 7.7â54.4 GeV RHIC_BES_spectra . The enhancement of
$\bar{\Lambda}/K^{0}_{S}$ ratios at intermediate $p_{T}$ in central Au+Au
collisions compared to peripheral Au+Au at $\sqrt{s_{NN}}$ $\geq$19.6 GeV is
interpreted as a consequence of hadron formation through parton
recombination.We do not see any significant difference between central and
peripheral collisions for $\sqrt{s_{NN}}$ $\leq$ 11.5 GeV. This may indicates
formation of hadronic interaction dominated matter at $\sqrt{s_{NN}}$ $\leq$
11.5 GeV.

Figure 3: The $\bar{\Lambda}/K^{0}_{S}$ ratio as a function of $p_{T}$ for
different centralities in Au+Au collisions at $\sqrt{s_{NN}}$ = 7.7â54.4GeV
RHIC_BES_spectra .

###  2.3 Nuclear modification factor

Probing QGP through parton energy loss:
The nuclear modification factor (RCP) is defined as the ratio of the particles
yield in the central to peripheral collisions normalized by the number of
binary collisions (Nbin). The value of Nbin is calculated from the Monte Carlo
Glauber simulation. If RCP is equal to one, then the nucleus-nucleus collision
is simply the superposition of nucleon-nucleon collisions. Deviation of RCP
from the unity would imply contribution from the nuclear medium effects. The
measured RCP, which is less than one at high $p_{T}$, at top RHIC and LHC
energies indicates the energy loss of partons traversing through QGP medium.
Figure 4 show RCP of strange hadrons measured in Au+Au collisions at
$\sqrt{s_{NN}}$ = 7.7â54.4GeV RHIC_BES_spectra . From figure 4 one can see, at
the intermediate $p_{T}$ , RCP goes above unity with the decrease in beam
energy. This indicates that at lower beam energy the parton energy loss effect
could be less important.

Figure 4: The nuclear modification factor of strange hadrons measured in Au+Au
collisions at $\sqrt{s_{NN}}$ = 7.7â54.4GeV RHIC_BES_spectra .

##  3 Collective Flow of Strange Hadrons

Partonic collectivity at top RHIC and LHC energies:
Figure 5 shows the elliptic flow coefficient ($v_{2}$) as a function of
$p_{T}$ for $\pi$, $p$, $\phi$, and $\Omega$ for 0â80% centrality in Au+Au
collisions at $\sqrt{s_{NN}}$ = 200 GeV and for $\phi$, and $\Lambda$ in
30â40% central Pb+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV RHIC_200_v2 ;
LHC_5TeV_v2 . We see mass ordering in $v_{2}$ at low $p_{T}$ and particle type
dependence at the intermediate $p_{T}$. The observed mass ordering can be
explained by hydrodynamics models whereas particle type dependence at the
intermediate $p_{T}$ can be explained through the parton recombination model.
It is clear from figure 5 that the $v_{2}$ of hadrons consisting of only
strange quarks ($\phi$ and $\Omega$) is similar to that of non-strange hadrons
$\pi$ and $p$. As we know that the $\phi$ and $\Omega$ do not participate
strongly in the hadronic interactions unlike $\pi$ and $p$, therefore figure 5
suggests that the major part of collectivity is developed during the partonic
phase at top RHIC energy. This is also true for LHC energies as shown by the
ALICE experiment LHC_5TeV_v2 .
Turn-off of QGP and change of equation of states ?
The energy dependence of $\phi$ meson $v_{2}$ measured at intermediate $p_{T}$
is shown in figure 6 RHIC_BES_v2 . We find that the magnitude of $v_{2}$ of
$\phi$ is almost constant ($\sim$0.1) for $\sqrt{s_{NN}}$ $\geq$ 19.6 GeV.
However, for $\sqrt{s_{NN}}$ =11.5 and 7.7 GeV, magnitude of $\phi$ $v_{2}$ is
small or close zero within large uncertainties. Small $\phi$ $v_{2}$ is
considered as a signature of formation of hadronic interaction dominated
matter as $\phi$ mesons are expected to have small hadronic interaction cross
sections. The energy dependence of the slope of $\phi$ directed flow ($v_{1}$)
is shown in the right panel of figure 6 RHIC_BES_v1 . There are indications of
change in the sign of slope of $v_{1}$ of $\phi$ below $\sqrt{s_{NN}}$ = 10
GeV. This could be related to the change of equation of state at the low beam
energy. The uncertainties in the current measurement are large, upcoming BES-
II program at RHIC is expected to provide precision measurement of $\phi$
$v_{1}$ and $v_{2}$.

Figure 5: The elliptic flow as a function of $p_{T}$ for $\pi$, $p$, $\phi$,
and $\Omega$ for 0%â80% centrality in Au+Au collisions at $\sqrt{s_{NN}}$ =
200 GeV and for $\phi$, and $\Lambda$ in 30%â40% central Pb+Pb collisions
$\sqrt{s_{NN}}$ = 5.02 TeV RHIC_200_v2 ; LHC_5TeV_v2 .

Figure 6: The energy dependence of $\phi$ $v_{2}$ at intermediate $p_{T}$ and
slope of directed flow ($dv_{1}/dy$) measured by the STAR experiment
RHIC_BES_v2 ; RHIC_BES_v1 .

##  4 Summary

The energy dependence of strange hadron production in heavy-ion collisions is
presented. It is found that the thermal model explains the measured particle
ratios (e.g. $K/\pi$ and $\phi/K$) and its energy dependence. The ratio
$\phi/K$ is found to be a sensitive probe to study the strangeness canonical
suppression. We observed thermal model calculations over-predicts $K^{*0}/K$
at the central Pb+Pb collisions. Also, the $K^{*0}/K$ ratio decreases with
increasing event multiplicity. This clearly indicates the effect of hadronic
re-scattering on short-lived resonance in central A+A collisions. The measured
elliptic flow of strange hadrons found to be similar to that of non-strange
hadrons at top RHIC and LHC energies, indicating formation of partonic
collectivity.
The measure $\phi$ $v_{2}$, at intermediate $p_{T}$, is found to be close zero
within large uncertainties at $\sqrt{s_{NN}}$ =11.5 and 7.7 GeV. The $\phi$
$v_{1}$ slope changes sign below $\sqrt{s_{NN}}$ = 10 GeV. The nuclear
modification factor is greater that unity at the intermediate $p_{T}$ for
$\sqrt{s_{NN}}$ $\leq$ 19.6 GeV. The baryon-to-meson
($\bar{\Lambda}/K^{0}_{S}$) ratios as a function of $p_{T}$ is found to be
centrality independent for $\sqrt{s_{NN}}$ =11.5 and 7.7 GeV. These results
indicate the formation of hadronic interaction dominated matter below
$\sqrt{s_{NN}}$ = 20 GeV.

## References

  * (1) J. Rafelski and B. Muller, Phys. Rev. Lett. 48, 1066 (1982).
  * (2) J. Adam et al. (STAR), Nucl. Phys. A 757, 102 (2005).
  * (3) A. Shor, Phys. Rev. Lett. 54, 1122 (1985).
  * (4) L. Ahle et al. (E866, E917), Phys. Lett. B 476, 1 (2000).
  * (5) L. Ahle et al. (E-802, E-866), Phys. Rev. C 60, 044904 (1999).
  * (6) L. Ahle et al. (E866, E917), Phys. Lett. B 490, 53 (2000).
  * (7) S. Afanasiev et al. (NA49), Phys. Rev. C 66, 054902 (2002).
  * (8) C. Alt et al. (NA49), Phys. Rev. C 77, 024903 (2008).
  * (9) L. Adamczyk et al. (STAR), Phys. Rev. C 96, 044904 (2017) and sQM2021 proceedings.
  * (10) M. Gazdzicki and M. I. Gorenstein, Acta Phys. Polon. B 30, 2705 (1999).
  * (11) J. Cleymans, et al. Phys. Lett. B 615, 50 (2005).
  * (12) J. Randrup and J. Cleymans, Phys. Rev. C 74, 047901 (2006).
  * (13) A. Andronic, et al. Nucl. Phys. A 772, 167 (2006)
  * (14) J. Cleymans et al. Phys.Lett. B 603, 146 (2004).
  * (15) G. Agakishiev et al. (HADES) Phys. Rev. C 80, 025209 (2009)
  * (16) M. S. Abdallah et al. (STAR), arXiv:2108.00924 and sQM2021 proceedings.
  * (17) S. Acharya, et al. ( ALICE) PLB 802, 135225 (2020).
  * (18) J. Adam et al. (STAR), Phys. Rev. C 102, 34909 (2020) and sQM2021 proceedings.
  * (19) L. Adamczyk et al. (STAR), Phys. Rev. Lett. 116, 062301 (2016).
  * (20) S. Acharya, et al. ( ALICE) JHEP 09, 006 (2018).
  * (21) L. Adamczyk et al. (STAR), Phys. Rev. Lett. 110, 142301 (2013)
  * (22) L. Adamczyk et al. (STAR), Phys. Rev. Lett. 120, 62301 (2018) and sQM2021 proceedings.


| $\epsilon=0$ | $\epsilon=0.05$ | $\epsilon=0.1$
NPL-MMD | 0.0107 (0.00527) | 0.00889 (0.00699) | 0.0113 (0.00538)
WABC | 0.743 (0.0358) | 0.768 (0.0624) | 0.757 (0.0242)
NPL-WLL | 0.00510 (0.00293) | 0.945 (0.0572) | 3.57 (0.0967)
NPL-WAS | 0.00689 (0.00333) | 0.0189 (0.00627) | 0.0397 (0.0129)
MMD-ABC | 0.750 (0.105) | 0.755 (0.0451) | 0.760 (0.0411)
Table 3: Experiment results for the G-and-k distribution over 10 runs

Method | NMSE (std)
---|---
| $\epsilon=0$ | $\epsilon=0.05$ | $\epsilon=0.1$
NPL-MMD | 0.00791 (0.00524) | 0.0128 (0.0117) | 0.0593 (0.0255)
WABC | 0.00142 (0.000965) | 0.585 (0.0188) | 0.532 (0.0109)

Table 4: Experiment results for the Toggle Switch model over 5 runs

Method | NMSE (std)
---|---
NPL-MMD | 0.338 (0.277)
WABC | 13.99 (13.8)

##  Appendix D ADDITIONAL EXPERIMENTS

We provide some additional experiments for our method by considering a type of
misspecification other than a contamination model, comparison with MMD-Bayes
in [57] and exploring sensitivity to the DP prior (through $\alpha$ and
$\mathbb{F}$), Gaussian kernel length scale $l$ and truncation limit $T$.

###  D.1 Misspecified Gaussian Location Model

So far in our empirical experiments we have considered the contamination model
which is canonical for analysing misspecification in robust statistics, mostly
because of theoretical convenience. While the model is simple, methods that
are robust against contamination models usually fare well for more practically
relevant alternatives, too. Figure 9 illustrates this on a new numerical
example: We generate Cauchy-distributed data, but wrongly fit a Gaussian to
it. The plot shows the posterior marginals for the location parameter, and
vertically marks its true location. The MSEs over 20 runs was 0.0289 (NPL-
MMD), 48.2 (NPL-WLL) and 28.6 (Standard Bayes).

Figure 9: Marginal posterior distribution for mean of Normal location model
with Cauchy data.

###  D.2 Comparison to MMD-Bayes for the Gaussian Location Model

We further compare our method to MMD-Bayes (KernelScore) in [57] on an
$\varepsilon$-contaminated Gaussian location model with outliers at location
$z$, where the weight is chosen using the grid search as in Section 4.2 of
[57]. Figure 10 below, shows the marginal posterior distributions for Standard
Bayes, Kernel Score and the MMD Posterior bootstrap methods for an increasing
number of outliers and location parameter.

Figure 10: Contaminated Normal location model.

###  D.3 Sensitivity to Hyperparameters

In this section we empirically examine the sensitivity of the proposed method
to several hyperparameters for the G-and-k distribution model.

####  D.3.1 Sensitivity to the DP Prior

We first set $T=n$ in (2.3) and $\mathbb{F}$ to be the Normal distribution
with parameters equal to the mean and standard deviation of the observed data.
Here we examine the effect of the prior in two ways; first by altering the
hyperparameter $\alpha$ of the DP prior which characterizes how much certainty
we impose on $\mathbb{F}$ and second by choosing an increasingly ‘worse’ prior
$\mathbb{F}$, since a higher proportion of outliers leads to a worse empirical
estimates of the mean and standard deviation in the Normal prior.

We generate $B=2^{9}$ bootstrap samples $\theta_{1},\dots,\theta_{B}$ for
different values of $\alpha$ ranging in $[0.01,300]$ and take the mean
estimator $\hat{\theta}=\frac{1}{B}\sum_{b=1}^{B}\theta_{b}$ for each value of
$\alpha$. In Figure 11, we plot (a) the normalised mean squared error for each
estimator $\hat{\theta}$ for an increasing value of $\alpha$ and (b) samples
from the generator with parameter $\hat{\theta}$. We observe that in the well-
specified case, $\alpha$ has no significant effect in inference, however as we
would expect, a larger value of $\alpha$, in combination with a worse prior
centering measure increasingly affects the parameter inference.

(a) Normalised mean squared error of the obtained estimator for an increasing
value of $\alpha$, corresponding to a higher confidence in the prior centering
measure $\mathbb{F}$.

(b) Generator samples obtained from each estimator $\hat{\theta}$. Lighter
(resp. darker) curves correspond to a smaller (resp. larger) value of
$\alpha$. The dotted line curve corresponds to the observed dataset in the
well-specified case.

Figure 11: Sensitivity to DP prior.

####  D.3.2 Sensitivity to $T$

We further illustrate the effect of the truncation limit $T$ in the
approximation of the DP posterior measure in (2.3). We fix $\alpha=0.1$ and
set $\mathbb{F}$ as in section D.3.1. We generate $B=2^{9}$ bootstrap samples
$\theta_{1},\dots,\theta_{B}$ for different values of $T$ ranging in
$[10,5000]$ and take the mean estimator
$\hat{\theta}=\frac{1}{B}\sum_{b=1}^{B}\theta_{b}$ for each value of $T$.
Figure 12 shows (a) the normalised mean squared error for each estimator
$\hat{\theta}$ for an increasing value of $T$ and (b) samples from the
generator with parameter $\hat{\theta}$. We observe that in this example, the
method is not significantly sensitive to the choice of $T$ for all values of
$\epsilon$.

(a) Normalised mean squared error of the obtained estimator for an increasing
value of $T$.

(b) Generator samples obtained from each estimator $\hat{\theta}$. Lighter
(resp. darker) curves correspond to a smaller (resp. larger) value of $T$. The
dotted line curve corresponds to the observed dataset in the well-specified
case.

Figure 12: Sensitivity to parameter $T$, the truncation limit of the DP
posterior.

####  D.3.3 Sensitivity to the Hyperparameter of Gaussian Kernel

We lastly investigate the effect of the length scale $l$ of the Gaussian
kernel. We fix $\alpha=0.1$, $T=n$ and set $\mathbb{F}$ as in D.3.1. We
generate $B=2^{9}$ bootstrap samples $\theta_{1},\dots,\theta_{B}$ for
different values of the length scale $l$ of the Gaussian kernel ranging in
$[10^{-1},10^{2}]$ and take the mean estimator
$\hat{\theta}=\frac{1}{B}\sum_{b=1}^{B}\theta_{b}$ for each value of $l$.
Figure 13 shows (a) the normalised mean squared error for each estimator
$\hat{\theta}$ for an increasing value of $l$ in logarithmic scale and (b)
samples from the generator with parameter $\hat{\theta}$.

(a) Normalised mean squared error of the obtained estimator for an increasing
value of length scale (in logarithmic scale).

(b) Generator samples obtained from each estimator $\hat{\theta}$. Lighter
(resp. darker) curves correspond to a smaller (resp. larger) value of $l$. The
dotted line curve corresponds to the observed dataset in the well-specified
case.

Figure 13: Sensitivity to the length scale of the Gaussian kernel.

To get some more intuition we plot the MMD loss as a function of
$\theta_{3}=g$ around a neighborhood of the true parameter value
$\theta_{0}=1$. Since there is noise in our estimate of the MMD, it is
possible that we get a global minimum for some value of $\theta_{3}$ far away
from one. This is unlikely to happen for $l<1$ because there is a much bigger
dip near $\theta_{0}$ as can be seen in figure 14. Of course, this is just a
projection of the loss landscape during optimisation, however it gives some
intuition as to why a small choice of length scale in this model leads to
better results.

Figure 14: Estimate of the MMD loss as a function of $\theta_{3}$ in the
G-and-k distribution for different values of the length scale of the Gaussian
kernel.