\section*{Introduction}
\subsection*{Lipschitz extension problem}
Let $X$, $Y$ be metric spaces.
Consider $X'\subset X$ and a Lipschitz map $f:X'\rightarrow Y$.
Can we extend $f$ to $F:X\rightarrow Y$ with the same constant $\mathrm{Lip}(F)=\mathrm{Lip}(f)$?
Failing that, can we bound the loss?
This potential ``loss'' can be encapsulated in a function $\mathcal{L}_{X,Y}$:
$$\begin{array}{rlcl} \mathcal{L}_{X,Y}: & \mathbb{R}^+ & \longrightarrow & \quad \mathbb{R}^+ \\ & C & \longmapsto &\displaystyle
\sup_{\substack{X'\subset X \\ f:X'\rightarrow Y \\ \mathrm{Lip}(f)\leq C}} \;
\inf_{\substack{F:X\rightarrow Y \\ F|_{X'}=f}} \;
\mathrm{Lip}(F).\end{array}$$
For example, maps to~$\mathbb{R}$, or more generally to a metric tree $T$, can always be extended without loss~\cite{m, las}: $\mathcal{L}_{X, \mathbb{R}}(C)=\mathcal{L}_{X, T}(C)=C$ for all $C\geq 0$.
Kirszbraun~\cite{k} proved that $\mathcal{L}_{X, Y}(C)=C$ when $X,Y$ are Euclidean spaces.

Recall that a Hadamard manifold is a complete, simply connected Riemannian manifold of nonpositive sectional curvature.
Lang and Schr\"oder~\cite{las}, extending work of Valentine~\cite{v} for the constant-curvature case, proved:
\begin{theoremA} \cite{las}
Let $\kappa_0, \kappa'_0<0$ be constants.
If $X,Y$ are Hadamard manifolds with $\kappa_X \geq \kappa_0$ and $\kappa_Y \leq \kappa'_0$, then $\mathcal{L}_{X,Y}(C)=C$ for all $C\geq \sqrt{\kappa_0/\kappa'_0}$.
\label{lasch}
\end{theoremA}
\subsubsection*{Main result}
Up to scaling, we may and always will assume $\kappa_0=\kappa'_0=-1$.
In that case, the above theorem also gives: $\mathcal{L}_{X,Y}(C)\leq 1$ when $C\leq 1$.
The goal of this note is to prove the following refinement:
\begin{theorem} \label{main}
For any $C<1$, $K\leq -1$ and $m \in \mathbb{N}$, there exists $C'<1$ such that
for any Hadamard manifolds $X,Y$ of dimension $\leq m$
satisfying $\kappa_X\geq -1 \geq \kappa_Y \geq K$,
one has $\mathcal{L}_{X,Y}(C)\leq C'$.
\end{theorem}
The $X=Y=\mathbb{H}^2$ case was conjectured in~\cite[App.\ C]{guk}, which put forward a strategy when $X'$ has bounded diameter.

\subsection*{About the method}
Our proof is based on the template of Lang and Schr\"oder's proof of Theorem~A, which we will recall in~\S\ref{prolasch} (slightly simplified, as~\cite{las} is set in the context of Alexandrov spaces).
The extra ingredients, which extend and uniformize arguments of~\cite{guk}, are based on the notion that under negative curvature, both in the small-scale limit (Euclidean geometry) and large-scale limit (real trees), loss-less extension is known to hold.
Thus, loss ($\mathcal{L}_{X,Y}(C)>C$) is in a sense a medium-range phenomenon, and can be controlled using a form of compactness and covering arguments.

When extending $f:X'\rightarrow Y$ to a single point $\xi \in X\smallsetminus X'$, we will see in~\S\ref{prolasch} that there is usually a natural ``optimal'' image $F(\xi)$, relative to the set $X'$ where the map is already defined.
Given a second point $\xi'$, we can then assign it an optimal image relative to $X'\cup \{\xi\}$, then pass to a third point $\xi''$ and so on, studying the loss incurred at each step.
One difficulty, which could cause the losses to pile up, is that the notion of ``optimal'', being relative to $X'\cup\{\xi, \xi', \dots\}$, changes as we go.

However, as pointed out in~\cite{las}, this difficulty \emph{disappears} when $Y$ is a metric tree: then, taking each $\xi\in X\smallsetminus X'$ to its optimal image (relative to $X'$ only) yields a globally Lipschitz map, with no loss.
This key feature, together with the fact that the curvature bounds force $Y$ coarsely to behave somewhat like a tree at large distances, is what allows us to prove Theorem~\ref{main}.
To patch together maps defined on different regions of $X$, we will use a standard interpolation procedure described in~\S\ref{interpol}.

\subsubsection*{Plan}
Section~\S\ref{prolasch} recalls the proof of Theorem~A; Section~\S\ref{secmain} proves Theorem~\ref{main}. Section~\S\ref{seconj} indulges in some speculation.

\subsection*{Notation}
Distances in metric spaces are all denoted $d$.

\noindent The open ball centered at $\xi$, of radius $r$, is written $\mathbb{B}_\xi(r)$. For a ball of unspecified center, we sometimes write $\mathbb{B}(r)$.

\noindent Given a point $\xi$ in a Hadamard manifold $X$, we write $\exp_\xi: \mathrm{T}_\xi(X) \rightarrow X$ the exponential map, and $\log_\xi$ its inverse.

\noindent Given $x,z\in X\smallsetminus \{\xi\}$, the notation $\widehat{x\xi z} \in [0,\pi]$ then refers to the angle between vectors $\log_\xi(x)$ and $\log_\xi(z)$, for the Euclidean metric on $\mathrm{T}_\xi(X)$.

\noindent The volume measure on $X$ is written $\mathrm{Vol}_X$.

\subsection*{Acknowledgements}
My interest for the Lipschitz extension problem in negative curvature started with Fanny Kassel and our work~\cite{guk}, which tackled an equivariant version. I am grateful to her and to Urs Lang for discussions related to this material.

\section{Proof of Theorem~A}\label{prolasch}
To build loss-less extensions, it is enough to do it one point $\xi\in X$ at a time: indeed, we can then repeat for a dense sequence $(\xi_n)_{n\in \mathbb{N}}$ of $X$, and pass to all of $X$ by continuity.

Let $X'\subset X$ and $f:X'\rightarrow Y$ be $C$-Lipschitz with
$$C\geq 1,$$
where $X,Y$ are Hadamard manifolds subject to curvature bounds $\kappa_X \geq -1 \geq \kappa_Y$.
We can restrict attention to $X'$ compact, nonempty. Consider $\xi \in X\smallsetminus X'$: the function defined by
$$\begin{array}{rcl} \varphi_\xi:  Y & \longrightarrow & \mathbb{R}^+  \\
    y  & \longmapsto & \displaystyle \max_{x\in X'} \, \frac{d(y,f(x))}{d(\xi, x)}
\end{array}$$
is proper and convex on $Y$, hence achieves a minimum
\begin{equation}\label{cxi} C_\xi:=\min \varphi_\xi = \varphi_\xi (\eta) \end{equation}
at some $\eta\in Y$ (in fact unique).
We can think of $\eta$ as an ``optimal candidate for $F(\xi)$'': Theorem~A will follow if we can prove
$$C_\xi \leq C.$$
If $C_\xi \leq 1$ we are done. If $C_\xi\geq 1$, define the compact set
\begin{equation}\label{xxi} X_\xi:=\Big \{x\in X', \frac{d(\eta, f(x))}{d(\xi, x)}=C_\xi \Big \}. \end{equation}

The exponential of any linear hyperplane $V \subset \mathrm{T}_\eta Y$ separates $Y$ into two half-spaces, each of which contains points of $f(X_\xi)$ in its closure: if not, we could push $\eta$ towards $f(X_\xi)$ (perpendicularly to $V$) to reduce $\varphi_\xi(\eta)$, contradicting minimality.
Hence, $\eta$ belongs to the convex hull of some points $y_i=f(x_i)$, $1\leq i \leq n$ where $x_i \in X_\xi$:
\begin{equation} \label{convhull} \sum_{i=0}^n \lambda_i \log_\eta (y_i) = 0_\eta \in \mathrm{T}_\eta Y \end{equation}
for some reals $\lambda_i>0$.
Since $x_i\in X_\xi$, the lengths $\ell_i:=d(\xi,x_i)= \Vert \log_\xi (x_i) \Vert$ satisfy $C_\xi \ell_i = d(\eta, y_i)= \Vert \log_\eta (y_i) \Vert$.
We can then write
\begin{align*} 0 &\leq \Big \Vert C_\xi \sum_{i=0}^n \lambda_i \log_\xi (x_i) \Big \Vert^2 - \Big \Vert \sum_{i=0}^n \lambda_i \log_\eta (y_i) \Big \Vert^2 \\
&= C_\xi^2 \sum_{i,j} \lambda_i \lambda_j \,  \ell_i \ell_j \, \big( \cos \widehat{x_i \xi x_j} - \cos \widehat{y_i \eta y_j} \big )
\end{align*}
hence at least one summand with $i\neq j$ is $\geq 0$, which happens if and only if $\widehat{x_i \xi x_j} \leq \widehat{y_i \eta y_j}$.
Hence, up to reindexing, we may assume
\begin{equation} \theta := \widehat{x_1 \xi x_2} \leq \widehat{y_1\eta y_2} =:\theta'. \label{cot} \end{equation}

Let $\mathcal{D}_\theta(\ell, \ell')$ denote the distance, in the hyperbolic plane $\mathbb{H}^2$, between the far ends of two segments of lengths $\ell, \ell'$ starting from a common vertex, an angle $\theta$ apart.
A~well-known trigonometric formula gives explicitly
\begin{equation} \label{trigo}
\mathcal{D}_\theta(\ell, \ell') = \mathrm{Arccosh} (\cosh \ell \, \cosh \ell' - \sinh \ell \, \sinh \ell' \, \cos \theta)
\end{equation}
but we will mostly use the following facts: the function $\mathcal{D}_\theta$ is convex in its two arguments, vanishes at $(0,0)$, and depends monotonically on $\theta$.
The Cartan-Alexandrov-Toponogov or CAT($-1$) comparison inequalities, whose interesting history is recounted in~\cite{cat}, say that
\begin{equation} \label{CAT}
\begin{array}{llll}
\mathrm{a.} & d(x_1, x_2) & \leq & \mathcal{D}_{\theta}(\ell_1, \ell_2) \\
\mathrm{b.} & d(y_1, y_2) & \geq & \mathcal{D}_{\theta'}(C_\xi \ell_1, C_\xi \ell_2)
\end{array} \end{equation}
due to the curvature bounds $\kappa_X \geq -1 \geq \kappa_Y$.
Therefore,

\begin{align}
 C\, d(x_1, x_2) & \geq  d(y_1, y_2) & \text{(Lipschitz bound)} \notag \\
& \geq \mathcal{D}_{\theta'}(C_\xi\ell_1, C_\xi\ell_2) & \text{by \eqref{CAT}}&\mathrm{.b} \notag \\
& \geq \mathcal{D}_{\theta}(C_\xi\ell_1, C_\xi\ell_2) & \text{by \eqref{cot}} \notag \\
& \geq  C_\xi \, \mathcal{D}_{\theta}(\ell_1, \ell_2) & \label{conv} \\
& \geq  C_\xi \, d(x_1, x_2). & \text{by \eqref{CAT}}&\mathrm{.a} \notag
\end{align}
where \eqref{conv} uses convexity of $\mathcal{D}_\theta$ and $C_\xi\geq 1$. Hence $C_\xi\leq C$ as desired,
proving Theorem~A. \qed

\section{Proof of Theorem~\ref{main}} \label{secmain}

\subsection{One-point extension}
We start by bounding the loss for extensions to a single point.

\begin{lemma} \label{1point}
For any $C<1$ there exists $C^*<1$ such that for any Hadamard manifolds $X,Y$ satisfying $\kappa_X \geq -1 \geq \kappa_Y$, any $X'\subset X$ and any $\xi\in X\smallsetminus X'$, every $C$-Lipschitz map $f:X'\rightarrow Y$ has a $C^*$-Lipschitz extension to $X'\sqcup \{\xi\}$.
\end{lemma}
\begin{proof}
Take $f, C, \xi$ as in the statement and define $C_\xi\geq 0$ (as well as $\eta \in Y$, $X_\xi \subset X'$, $y_i=f(x_i)\in f(X_\xi)$ and $\ell_i=d(\xi, x_i)$) as in the previous proof.
Theorem~A gives $C_\xi\leq 1$: let us bound $C_\xi$ away from $1$ in terms of $C$ alone.

Let $\Delta>0$ be such that
\begin{equation} \label{addloss} \mathcal{D}_{\pi/2}(\ell, \ell') \geq \ell + \ell' - \Delta ~\text{ for all }~ \ell, \ell'\geq 0 \end{equation}
(using~\eqref{trigo} one can show $\Delta=\log 2$ works).
Let $r>0$ be large enough that
\begin{equation} \label{C'} \widehat C:=C+{\Delta}/{r}<1. \end{equation}
We distinguish two cases.

$\bullet$ If $\ell_i \geq r$ for some index $i$, we use~\eqref{convhull} to find $j\neq i$ such that
\begin{equation} \label{obtus} \theta':= \widehat{y_i \eta y_j} \geq \pi/2 \end{equation}
and write:
\begin{align}
C\, d(x_i, x_j) & \geq  d(y_i, y_j) & \text{(Lipschitz bound)} \notag \\
& \geq \mathcal{D}_{\theta'}(C_\xi\ell_i, C_\xi\ell_j) & \text{by \eqref{CAT}}&\mathrm{.b} \notag \\
& \geq C_\xi (\ell_i + \ell_j) - \Delta & \text{by \eqref{addloss}--\eqref{obtus}} \notag \end{align}
hence $$C_\xi\leq
\widehat C$$
by \eqref{C'}, due to the triangle inequality $d(x_i, x_j) \leq \ell_i+\ell_j$ and $\ell_i \geq r$.

$\bullet$ If no such index $i$ exists, then we define $x_1, x_2\in X_\xi$ and $\theta\leq \theta' \in [0,\pi]$ as in the proof of Theorem~A and write, similar to~\eqref{conv}:
\begin{align}
C\, d(x_1, x_2) & \geq  d(y_1, y_2) & \text{(Lipschitz bound)} \notag \\
& \geq \mathcal{D}_{\theta'}(C_\xi\ell_1, C_\xi\ell_2) & \text{by \eqref{CAT}}&\mathrm{.b} \notag \\
& \geq \mathcal{D}_{\theta}(C_\xi\ell_1, C_\xi\ell_2) & \text{by \eqref{cot}} \notag \\
& \geq C'_\xi \, \mathcal{D}_{\theta}(\ell_1, \ell_2) & \text{(\emph{see \eqref{newc} below})}\label{pattes} \\
& \geq C'_\xi \, d(x_1, x_2) & \text{by \eqref{CAT}}&\mathrm{.a} \notag
\end{align}
where we use the new constant
\begin{equation} \label{newc} C'_\xi:= \frac{\sinh (C_\xi r)}{ \sinh (r)}. \end{equation}
Indeed, for a basepoint $o\in \mathbb{H}^2$, the differential of the exponential map $\mathrm{exp}_o:(\mathbb{R}^2,0) \rightarrow (\mathbb{H}^2,o)$ at
a point of the circle $\partial \mathbb{B}_0(\lambda)$
has principal values~$1$ radially and $\sinh (\lambda)$ along the circle --- this can be checked by differentiating~\eqref{trigo} near $(\ell, \ell', \theta)=(\lambda, \lambda, 0)$.
It follows that the radial map
$$ \begin{array}{rrcl} H: & (\mathbb{H}^2,o) & \longrightarrow & (\mathbb{H}^2,o) \\
& x & \longmapsto & \exp_{o} \big (C_\xi \log_o (x) \big ), \end{array} $$
defining a homothety of ratio $C_\xi$ on each line through $o$, satisfies
$$\mathrm{Lip}\left ( \left . H^{-1} \right |_{\mathbb{B}_o(C_\xi r)} \right ) = \frac{\sinh (r)}{\sinh(C_\xi r)}= \frac{1}{C'_\xi}$$
which means that step \eqref{pattes} holds (using $\ell_1, \ell_2\leq r$).
Therefore, $C'_\xi \leq C$.
Substituting in~\eqref{newc}, we find
\begin{equation} \label{arcsinh} C_\xi \leq \frac{\mathrm{Arcsinh}(C\sinh (r))}{r}<1. \end{equation}

In either case, we have bounded the Lipschitz constant $\mathrm{max}\{C,C_\xi\}$ (for the one-point extension $\xi \mapsto \eta$ of $f$) uniformly away from $1$.
\end{proof}

\subsection{Averaging maps} \label{interpol}
In curvature $\leq 0$, convex interpolation behaves well with respect to Lipschitz constants.
Namely, given $f_0, f_1 : X \rightarrow Y$, let $(f_t(x))_{t\in [0,1]}$ be the constant-speed parametrization of the geodesic segment $[f_0(x), f_1(x)]$, for all $x\in X$.
The ``barycenter'' maps $f_t:X\rightarrow Y$ thereby defined satisfy:
\emph{if $f_1$ agrees with $f_0$ on $X' \subset X$ then so does $f_t$.}

Moreover, for all $x,x'\in X$, if $(y_t)_{t\in [0,1]}$ denotes the constant-speed para\-metrization of the segment $[f_0 (x), f_1(x')]$, then
\begin{align*}
 d(f_t(x),f_t(x')) & \leq d(f_t(x),y_t) + d(y_t,f_t(x')) \\
 & \leq t\,  d(f_1(x),f_1(x')) + (1-t)\, d(f_0(x),f_0(x'))
\end{align*}
by CAT(0) comparison inequalities. It follows that
$$ \mathrm{Lip}(f_t) \leq t\, \mathrm{Lip}(f_1) + (1-t)\, \mathrm{Lip}(f_0).$$
We will simply use the \emph{notation} $$f_t =: t\, f_1 + (1-t)\, f_0.$$
We can also iterate the construction above, to define barycenters of $N$ maps: given maps $(f_i)_{i\geq 1}$, the maps $F_N=\sum_{i=1}^N \frac{1}{N} f_i$, defined inductively on $N$ by
$F_N : = \frac{1}{N} f_N + \frac{N-1}{N} \big ( \sum_{i=1}^{N-1} \frac{1}{N-1} f_i \big )$, inductively satisfy for any $Z\subset X$:
\begin{equation} \label{moyenne}
\mathrm{Lip}(F_N|_Z) \leq \sum_{i=1}^N \frac{1}{N} \mathrm{Lip}(f_i|_Z).
\end{equation}
When $N\geq 3$ this construction is not robust under permutation of the $f_i$; note however that symmetric constructions do exist~\cite{las}, which also satisfy a weakened form of associativity~\cite{guk}.

\subsection{Extensions to the whole space}
We now prove Theorem~\ref{main}.
Let $C<1$, $K\leq -1$, $m\in \mathbb{N}$ and Hadamard manifolds $X,Y$ be as in the theorem, and $C^*\in[C,1)$ be given by Lemma~\ref{1point}.

Let $f:X'\rightarrow Y$ be a $C$-Lipschitz map, where $X'\subset X$.
Again, we may assume $X'$ is compact.
By Lemma~\ref{1point}, we may consider a family of $C^*$-Lipschitz extensions $(f_\xi^*)_{\xi\in X}$ to $X'\cup\{\xi\}$, taking $\xi$ to its optimal candidate image.
We do allow $\xi \in X'$, in which case $f_\xi^*=f$.
Small balls in $X$ and $Y$ are uniformly $(1+o(1))$-bi-Lipschitz to Euclidean balls, by the curvature bounds $0\geq \kappa_X, \kappa_Y \geq K$
(in fact CAT-type inequalities~\eqref{CAT} show that this $o(1)$ tolerance is quadratic in the size of the balls).
By composition, loss-less extension in Euclidean geometry~\cite{k} implies that there exists $\varepsilon_0 \in (0,1)$ such that each $f_\xi^* \big |_{\mathbb{B}_\xi(\varepsilon_0) \cap (X'\cup \{\xi\})}$
has a $\sqrt{C^*}$-Lipschitz extension
\begin{equation} \label{patch} \widehat{f_\xi} : \mathbb{B}_\xi(\varepsilon_0) \longrightarrow Y. \end{equation}
Let $\varepsilon < \varepsilon_0$ be small enough, and $R>1$ large enough, that
\begin{equation} \label{buffer} \mathrm{(a)} \quad  \frac{C^* +  \varepsilon/\varepsilon_0}{1-  \varepsilon/\varepsilon_0}\leq 1
\quad\quad \text{and} \quad\quad \mathrm{(b)} \quad
\frac{(C^*+{\Delta}/{R})+2\varepsilon/R}{1-2\varepsilon/R} \leq 1\end{equation}
where $\Delta>0$ still satisfies~\eqref{addloss}.

\begin{lemma} \label{qtree}
Let $\xi, \xi' \in X$ be distance $\geq R$ apart.
Then $\mathrm{Lip}(G)\leq 1$ for $$G:= f \sqcup \widehat{f_\xi}\big |_{\mathbb{B}_\xi(\varepsilon)} \sqcup \widehat{f_{\xi'}} \big |_{\mathbb{B}_{\xi'}(\varepsilon)}.$$
\end{lemma}
\begin{proof}
Consider $x,x'\in X'\cup \mathbb{B}_\xi(\varepsilon) \cup \mathbb{B}_{\xi'}(\varepsilon)$. We distinguish several cases.

\noindent $\bullet$ (i) If $x,x'\in X'$ then $d(G(x), G(x'))=d(f(x), f(x'))\leq C d(x,x')$ because $f$ is $C$-Lipschitz.

\noindent $\bullet$ (ii) If $x,x'\in \mathbb{B}_\xi(\varepsilon)$ then by construction of $\widehat{f_\xi}$,
\begin{equation} \label{nest}
d(G(x), G(x'))=d(\widehat{f_\xi}(x), \widehat{f_\xi}(x'))\leq \sqrt{C^*} d(x,x').
\end{equation}

\noindent $\bullet$ (iii) If $x,x'\in \mathbb{B}_{\xi'}(\varepsilon)$, we do as in (ii), exchanging $\xi$ and $\xi'$.

\noindent $\bullet$ (iv) If $x\in \mathbb{B}_\xi(\varepsilon)$ and $x'\in X'$, we distinguish two cases:
if $x'\in X'\cap \mathbb{B}_\xi(\varepsilon_0)$, then~\eqref{nest} still applies. If not, then
we compute
$$ \frac{d(G(x'), G(x))}{d(x',x)}
\leq \frac{d(G(x'), G(\xi)) + d(G(\xi), G(x))}{d(x',\xi) - d(\xi, x)}
 \leq \frac{C^* d(x',\xi) + d(\xi, x)}{d(x',\xi) - d(\xi, x)} $$
which is $\leq 1$ by~\eqref{buffer}.a, since $d(\xi,x)\leq \varepsilon$ and $d(x',\xi)\geq \varepsilon_0$.

\noindent $\bullet$ (v) If $x\in \mathbb{B}_{\xi'}(\varepsilon)$ and $x'\in X'$, we do as in (iv), exchanging $\xi$ and $\xi'$.

\noindent $\bullet$ (vi) Up to exchanging $x$ and $x'$, the only remaining case is that $x\in \mathbb{B}_{\xi}(\varepsilon)$ and $x'\in \mathbb{B}_{\xi'}(\varepsilon)$.
It is only here that we will use the assumption $d(\xi, \xi')\geq R$.

We first treat the case $(x,x')=(\xi, \xi')$.
Recall from~\eqref{cxi} the optimal candidates $\eta=G(\xi)$ and $\eta'=G(\xi')$ and optimal constants
$C_\xi, C_{\xi'}<1$ used in the proofs of Theorem~A and Lemma~\ref{1point}.
By symmetry, we may assume
\begin{equation} \label{match3}C_{\xi'} \leq C_\xi  \end{equation}
and by definition of $C_{\xi'}$ we have
\begin{equation} \label{match} {d(\eta',f(z))} \leq C_{\xi'}\, {d(\xi', z)} ~\text{ for all }~z\in X'.\end{equation}
Recall also from~\eqref{xxi} the compact subset $X_\xi \subset X'$, satisfying
\begin{equation} \label{match2} {d(\eta,f(z))} = C_{\xi}\, {d(\xi, z)} ~\text{ for all }~z\in X_\xi. \end{equation}
By Lemma~\ref{1point} we know
\begin{equation} \label{match4}  C_\xi \leq C^* < 1. \end{equation}
Since $\eta$ lies by~\eqref{convhull} in the convex hull of $f(X_\xi)$, we can find $y_1=f(x_1) \in f(X_\xi)$ such that $\widehat{\eta' \eta y_1} \geq \frac{\pi}{2}$. Then,
\begin{align}
 d(\eta, \eta')  & \leq d(y_1,\eta')-d(y_1, \eta)+\Delta & \text{by \eqref{CAT}.b--\eqref{addloss}} \notag \\
 & \leq C_{\xi'} d(x_1,\xi')-C_\xi d(x_1, \xi) + \Delta & \text{by \eqref{match}--\eqref{match2}} \notag \\
 & \leq C_\xi \left ( d(x_1,\xi')-d(x_1, \xi) \right )+ \Delta & \text{by \eqref{match3}} \label{comput} \\
 & \leq C_\xi \, d(\xi', \xi) + \Delta & \text{(triangle inequality)} \notag \\
 & \leq C^* \, d(\xi', \xi) + \Delta & \text{by \eqref{match4}} & . \notag
\end{align}
 Since by assumption $d(\xi, \xi')\geq R$, it follows that
\begin{align}
{d(\eta, \eta')}/{d(\xi, \xi')} & \leq C^*+{\Delta}/{R}  \quad <1 & \text{by \eqref{buffer}.b.}  \label{buffer3} \end{align}
This
(\footnote{~For $Y$ a tree and $X$ a general metric space, a variant of the computation~\eqref{comput} holds with $\Delta=0$, and a variant of the argument in~\S\ref{prolasch} yields $C_\xi\leq C$.
Taking each $\xi, \xi' \in X \smallsetminus X'$ (independently) to its optimal image $\eta, \eta' \in Y$ therefore produces a global, loss-less extension of $f$: this was proved in~\cite[Th.~B]{las}, as alluded to in the Introduction.})
deals with the case $(x,x')=(\xi, \xi')$.

The general case of (vi) is now similar to (iv-v): we can compute
\begin{align*} \frac{d(G(x), G(x'))}{d(x,x')}
& \leq \frac{d(G(x), \eta) + d(\eta, \eta') + d(\eta', G(x'))}{-d(x,\xi) + d(\xi, \xi') - d(\xi', x')} \\
& \leq \frac{(C^* + \Delta/R) d(\xi,\xi') + 2\varepsilon }{d(\xi,\xi') - 2 \varepsilon} & \text{by~\eqref{buffer3}} \\
& = \frac{(C^* + \Delta/R) +{2\varepsilon}/{d(\xi,\xi')}}{1-{2\varepsilon}/{d(\xi,\xi')}} \quad \leq 1 & ~\text{by~\eqref{buffer}}&\mathrm{.b,} \end{align*}
using again $d(\xi, \xi')\geq R$.
Therefore, $\mathrm{Lip}(G)\leq 1$.
\end{proof}

To finish proving Theorem~\ref{main}, consider a maximal $\varepsilon$-sparse subset
$$\Xi=\big \{\xi_i \big \}_{i\in \mathbb{N}} \subset X.$$
This means that the closed balls $\overline{\mathbb{B}_{\xi_i}(\varepsilon)}$ cover $X$ but the $\mathbb{B}_{\xi_i}(\varepsilon/2)$ are pairwise disjoint (i.e.\ the $\xi_i \in \Xi$ are mutually $\geq \varepsilon$ apart).
For example, $\Xi$ can be constructed from a dense sequence $(x_\iota)_{\iota\in \mathbb{N}}$ of $X$ by setting $\xi_1:=x_1$ and letting inductively $\xi_i$ be the first $x_\iota$ lying outside $\mathbb{B}_{\xi_1} (\varepsilon) \cup \dots \cup \mathbb{B}_{\xi_{i-1}} (\varepsilon)$.

Since $0\geq \kappa_X \geq -1$, the volume of a ball in $X$ is bounded above (resp.\ below) by the volume of a ball of the same radius in hyperbolic space $\mathbb{H}=\mathbb{H}^{\dim(X)}$ (resp.\ in Euclidean space $\mathbb{E}=\mathbb{R}^{\dim(X)}$): indeed, CAT-type inequalities~\eqref{CAT} show that the \emph{Jacobians} of the exponential maps in $\mathbb{H}$, $X$, and $\mathbb{E}$ form, in that order, a weakly decreasing sequence.
Let $N \in \mathbb{N}$ satisfy
\begin{equation} \label{bins}
N \geq \frac{\mathrm{Vol}_{\mathbb{H}}(\mathbb{B}(R+\varepsilon/2))}{\mathrm{Vol}_{\mathbb{E}}(\mathbb{B}(\varepsilon/2))}.
\end{equation}
Each ball $\mathbb{B}_{\xi_i}(R)$ contains at most $N$ points of $\Xi$, because the $\varepsilon/2$-balls centered at those points are disjoint and contained in $\mathbb{B}_{\xi_i}(R+\varepsilon/2)$.
Therefore, we can find a partition of $\Xi$ into ``bins''
$$\Xi=\Xi_1 \sqcup \dots \sqcup \Xi_N$$
such that any distinct $\xi, \xi' \in \Xi_j$ satisfy $d(\xi, \xi')\geq R$:
for example, the $\Xi_j$ can be constructed inductively by putting $\xi_1$ in $\Xi_1$, and then dropping in turn each $\xi_i$ into any bin $\Xi_j$ disjoint from $\{\xi_1, \dots, \xi_{i-1} \} \cap \mathbb{B}_{\xi_i}(R)$.

Recall from~\eqref{patch} the $\sqrt{C^*}$-Lipschitz maps $\widehat{f_{\xi_i}}$ defined in $\varepsilon_0$-neighborhoods of the $\xi_i$.
For each $1\leq j \leq N$, define the map
$$F_j:= \Big ( \bigsqcup_{\xi\in \Xi_j} \widehat{f_\xi}|_{\mathbb{B}_\xi(\varepsilon)} \sqcup f \Big ) :\quad  \bigcup_{\xi \in \Xi_i} \mathbb{B}_\xi(\varepsilon) \cup X' \longrightarrow Y.$$
By Lemma~\ref{qtree}, since the Lipschitz property can be tested one pair of points at a time, we have in fact $\mathrm{Lip}(F_j) \leq 1$.
By Theorem~A, the $F_j$ admit $1$-Lipschitz extensions $\widehat{F}_j$ to $X$.
Finally we claim that
\begin{equation} \label{finalC'} F:=\sum_{j=1}^N \frac{1}{N} \widehat{F}_j \: : X \longrightarrow Y
~\text{ satisfies }~
 \mathrm{Lip}(F) \leq 1 - \frac{1-\sqrt{C^*}}{N} =:C' <1. \end{equation}
Indeed, this can be verified in restriction to each ball $\overline{\mathbb{B}_{\xi_i}(\varepsilon)}$ of the covering of~$X$: if $\xi_i$ falls in the bin $\Xi_j$, then on that ball $\widehat{F}_j$ is $\sqrt{C^*}$-Lipschitz by construction while all other $\widehat{F}_{j'}$ are $1$-Lipschitz; we conclude using~\eqref{moyenne}. \qed

\section{Conclusion} \label{seconj}
It seems natural to expect that the lower bound $K$ on curvature, and the upper bound $m$ on dimension, are not necessary in Theorem~\ref{main}.
\begin{conjecture}\label{conj}
For any $C<1$ there exists $C'\in (C,1)$ such that for any Hadamard manifolds $X,Y$ satisfying $\kappa_X\geq -1 \geq \kappa_Y$, every $C$-Lipschitz map from a subset of $X$ to $Y$ has a $C'$-Lipschitz extension to $X$:
$$\mathcal{L}_{X,Y}(C)\leq C' < 1.$$
\end{conjecture}
\noindent This statement should still hold if both the map and its extension are required to be equivariant under a given pair of actions on $X$ and $Y$: see~\cite{guk}

\smallskip

Loss does occur, i.e.\ $C'>C$ in general, as testified by many examples.
For instance, since $\ell \mapsto \mathcal{D}_{2\pi/3} (\ell, \ell)$ is strictly convex (see~\eqref{trigo}), a map $f$ that takes just the vertices of, say, a medium-sized equilateral triangle of $\mathbb{H}^2$ to the vertices of a smaller one, cannot be extended
without loss to the center of the triangle.
In such examples however, the ratio
${(1-C')}/{(1-C)}$
never seems to get very small.
Thus we propose the following strengthening:
\begin{conjecture} \label{conj2}
There exists a universal $\alpha \in (0,1)$ such that $\mathcal{L}_{X,Y}(C) \leq C^\alpha$
\end{conjecture}
Interestingly, this conjecture appears to be open even for $C$ close to $0$.
The article~\cite{lps} shows that $\mathcal{L}_{X,Y}(C)/C$ is bounded above (which for small $C$ is a stronger property), but only under some extra assumptions on the Hadamard manifold~$Y$, such as fixed dimension with pinched curvature.

As $C$ approaches $1$, bounds on the constant $C'$ extracted from our proof of Theorem~\ref{main} are not very stringent.
Fixing $K \leq -1$ and the dimension, we can estimate~\eqref{arcsinh} for $r=\frac{2\Delta}{1-C}$ to find that $1-C^*$ is on the order of $(1-C)^2$, yielding
$\varepsilon_0 \approx 1-C$, $\varepsilon \approx (1-C)^3$ and crucially $R\approx (1-C)^{-2}$ in~\eqref{buffer}.
In~\eqref{bins} this entails
$N\approx \mathrm{e}^{-(\Lambda+o(1))/(1-C)^2}$
for some $\Lambda>0$,
hence in~\eqref{finalC'} $$1-C'\approx \mathrm{e}^{-\frac{\Lambda+o(1)}{(1-C)^2}} \quad \text{as } C \rightarrow 1^-,$$
i.e.\ our upper bound $C'$ for $\mathcal{L}_{X,Y}(C)$ is a far cry from Conjecture~\ref{conj2}.


\section{ Introduction}

It is a matter of fact that  many interesting dynamic problems in science and
engineering are characterized by the presence of
 a variety of degrees of freedom with very different time scales.
As  important examples we can mention proteins~\cite{karplus76} and climate dynamics~\cite{donges09}:
we remind that  the time scale of vibration of covalent bonds is
$O(10^{-12} s) $, while the folding time for proteins may be of the order
of seconds;
in a similar way  the characteristic times of the  processes
involved in climate  vary from  seconds
 for the $3D$ turbulence,  to days for the atmosphere, to   $O(10^4 yr)$
 for the deep ocean currents and ice shield dynamics.
\\
Due to the multi-scale character of such kind of systems,
is not possible to perform  a direct  simulation
of  all the relevant  involved scales,
even with the support of modern supercomputers  and advanced
numerical algorithms.
These practical difficulties force us to reduce our ambitions;
a (non-trivial) possibility in this sense is to describe  the  "slow dynamics" in terms
of effective equations.
Using such an approach  one  has  both practical and conceptual advantages:
for instance, it is possible to decrease the computational effort, e.g. by reducing the number of equations and adopting a "large" $\Delta t$;
in addition, the effective equations are able to catch some general features and to reveal dominant ingredients which can remain hidden in the
detailed description~\cite{tnt_entropy}.
\\
Disappointingly enough, only in few cases it is possible to derive effective equations with a systematic approach: important examples are dilute gases \cite{dorfman72},
 harmonic chains~\cite{rubin60, zwanzig73} and the Markovian limit of Hamiltonian dynamics \cite{spohn80}.
\\
On the other hand, in the history of science there is a series
of clever practical approaches for the study of multi-scale problems
that do not rely on rigorous derivations, e.g. the averaging method in  celestial mechanics~\cite{mitropolsky67},
the Langevin equation for colloids~\cite{langevin},
the homogenization for partial differential equations~\cite{pavliotis08},
the Born-Oppenheimer "approximation"~\cite{hutter09} and
the Carr-Parrinello method~\cite{car85}.
\\
In order to give an idea of the general methodology
let us briefly remind a  well known example of effective model,
the advection-diffusion equation for a passive scalar $\rho({\bf x},t)$
(e.g. the concentration of a pollutant)
in an incompressible flow $(\nabla \cdot {\bf u}=0)$:
\begin{equation}
 \label{eq:adv-diff}
{\partial \rho \over \partial t}+({\bf u} \cdot \nabla) \rho=
D_0 \Delta \rho \,\, .
\end{equation}
Maxwell had the idea, now  supported
 by mathematics (under rather general conditions), to consider the solution of Eq.~\eqref{eq:adv-diff}
at large scale and asymptotically in time; in these limits one obtains the so called standard diffusion,
i.e. a  Fick's law holds of the form
\begin{equation}
 {\partial \Theta \over \partial t}=
\sum_{i,j} {\cal D}_{ij} {\partial^2 \Theta \over \partial x_i \partial x_j}
\end{equation}
where $\Theta$ is the spatial coarse graining of $\rho$, and ${\cal D}_{ij}$ is the
effective (eddy) diffusion tensor, depending (often in a non trivial way)
on $D_0$ and the field ${\bf u}$
(just for simplicity we considered the case $\langle{\bf u}\rangle=0$).
If some (rather general) conditions on  the field ${\bf u}$ are satisfied then
one can use a precise protocol to compute the tensor ${\cal D}_{ij}$~\cite{biferale95}.
\\
In this paper we review some important aspects of the problem of finding effective
equations for complex systems. In Section~\ref{sec:generalizing} we review the classical
Langevin equation and show how it can be extended to cases in which the system
obeys a Hamiltonian with a generalized form of the kinetic energy; Section~\ref{sec:testing}
is devoted to the discussion of a data-driven method that allows to test such generalization;
Section~\ref{sec:granular} shows how this method can be applied to experimental cases, and how it can be used to
infer coarse-grained models whose behavior nicely agree with that of the real system; then in Section~\ref{sec:bigdata}
we briefly comment on the lessons that we can learn from the problem of effective equations in physics,
and how such warnings can reveal useful also in the context of big data and machine learning. In Section~\ref{sec:conclusion}
we briefly sketch our conclusions.

\section{Generalizing the Langevin Equation}
\label{sec:generalizing}
In his celebrated paper about Brownian motion, Paul Langevin addressed the problem of properly describing the irregular behavior of pollen particles suspended in water~\cite{langevin}. Following Einstein, he assumed that
both the colloidal particle and the molecules of the fluid could be modeled as material points with masses $M$ and $m\ll M$ respectively. The motion of the heavy particle is due to the collisions with the molecules of the liquid, which are assumed to be uncorrelated. To account for the discontinuous action of the hitting molecules, Langevin relied upon the introduction of a stochastic term in the evolution equation of the colloid, namely a white Gaussian noise $\boldsymbol{\xi}(t)=(\xi_1(t),\xi_2(t),\xi_3(t))$ such that
\begin{equation}
\begin{aligned}
 \langle \xi_i (t) \rangle&=0\\
 \langle \xi_i (t) \xi_j (t') \rangle&=\delta(t-t')\delta_{ij}\,.
\end{aligned}
\end{equation}
He supposed that the impulsive force acting on the colloid was proportional to this discontinuous function $\vect{\xi}$.
On the other hand, he argued that the interaction with the fluid results into an average damping force acting on the colloidal particle, and proportional to its velocity (Stokes law). The combination of the above effects leads to the celebrated Langevin Equation~(LE)
\begin{equation}
\label{eq:LE}
 \dot{\vect{P}}=-\gamma \vect{P}+\sqrt{2 \gamma k_B T} \vect{\xi}
\end{equation}
which characterizes the evolution of the momentum $\vect{P}=M\dot{\vect{Q}}$ of the heavy colloid ($\vect{Q}$ being its position in the three-dimensional space). Here $\gamma$ is the friction term due to the interaction of the colloidal particle with the fluid, while $T$ is the temperature and $k_B$ is the Boltzmann constant.
The noise amplitude is determined by the Einstein relation, which relies on the fact that the interested particle is at thermal equilibrium with the fluid, so that equipartition theorem holds and
\begin{equation}
\label{eq:equipartition}
 \langle P^2/M \rangle =3 k_B T\,.
\end{equation}
\\
Eq.~\eqref{eq:LE} clearly shows that the Brownian motion is the result of the competing actions of a damping force and a thermal noise. It seems reasonable that such mechanism should hold, under appropriate modifications, also for Hamiltonian systems with non-standard, generalized forms of the kinetic energy $K(P)$ (e.g., non quadratic functions of the momenta). In these cases there would be no reason for the damping force to be proportional to the momentum, and the equipartition theorem could assume a formulation very different from Eq.~\eqref{eq:equipartition}.
\\
This problem has been addressed in~\cite{baldovin18} and, as it will be discussed in the following, it assumes particular conceptual relevance when Hamiltonian systems living in bounded phase-spaces are taken into account, so that the absolute temperature of the system can assume negative values.
\\
Let us consider the general case of a Hamiltonian system in the form
\begin{equation}
\label{eq:ham_gen}
  H(P,\{p_n\},Q,\{q_n\})=K(P)+\sum_n \tilde{K}(p_n)+ U(Q) + V(Q,\{q_n\})
 \end{equation}
 where $(P,Q)$ is a ``slow'' degree of freedom. For example, in a system with the usual quadratic kinetic energy it could represent a particle with a mass much higher than the others ($K(P)=P^2/2M$, $\tilde{K}(p_n)=p_n^2/2m$, $M\gg m$). $U(Q)$ is the external potential which the slow particle is subjected to, while $V(Q,\{q_n\})$ takes into account the interactions occurring among different degrees of freedom. For the sake of simplicity we consider here only Hamiltonian systems in one dimension, but all the results can be straightforwardly generalized to the multi-dimensional case.
 \\
 In what follows we will limit our analysis to Hamiltonians of the form~\eqref{eq:ham_gen}, in which the kinetic energy is the sum of single-particle contributions only depending on the momentum.
 \\
 The Hamilton equations describing the motion of the slow degree of freedom read:
 \begin{equation}
 \label{eq:motion}
  \begin{cases}
   \dot{Q}=\partial_PK(P)\\
   \dot{P}=-\partial_QU(Q)-\partial_QV(Q,\{q_n\})
  \end{cases}
 \end{equation}
 At this stage we introduce a first, strong hypothesis, in the spirit of the one done by Langevin: we suppose that the time-scale separation between the dynamics of the slow particle and the fast ones allows us to approximate the former through an effective stochastic equation. In other words, we assume that Eqs.~\eqref{eq:motion} can be rewritten as
 \begin{equation}
 \label{eq:kleinkramers}
  \begin{cases}
   \dot{Q}=\partial_PK(P)\\
   \dot{P}=-\partial_QU(Q)+\Gamma(P)+\sqrt{2D}\xi(t)\,.
  \end{cases}
 \end{equation}
  Here the term $\Gamma(P)$ can be seen as a generalization of the Stokes force, while $D$ is a constant which determines the amplitude of the noise. In this way we are ignoring the details of the interactions between the slow and the fast degrees of freedom. The possibility to perform such averaging procedure on rigorous mathematical grounds is a non-trivial, largely studied problem in the field of dynamical systems~\cite{givon04, kifer04}: the above approximation should be therefore viewed as an ansatz, whose validity needs to be checked \textit{a posteriori}.
 \\
 We aim at finding some kind of generalized Einstein relation to relate the constant $D$ to $\Gamma(P)$. Let us introduce the steady probability density $f(Q,P)$ of the considered degree of freedom (to be determined) and the corresponding steady probability currents:
\begin{equation}
\begin{cases}
 J_Q(Q,P)=f(Q,P)\partial_PK(P)\\
 J_P(Q,P)=-f(Q,P)\partial_QU(Q,P)+\Gamma(P)f(Q,P)-D\partial_P f(Q,P)\,.
\end{cases}
\end{equation}
In terms of the above quantities, the Fokker-Planck equation corresponding to Eq.~\eqref{eq:kleinkramers} reads:
\begin{equation}
\label{eq:fokkerplanck}
 \partial_QJ_Q(Q,P)+\partial_PJ_P(Q,P)=0\,.
\end{equation}
\\
We assume now that the system is in thermal equilibrium (which is the same hypothesis done by Einstein and Langevin when
exploiting the equipartition theorem~\eqref{eq:equipartition}). We require therefore that detailed balance is satisfied, i.e. that the irreversible part of $J_P$ vanishes, so that
\begin{equation}
\label{eq:gamma}
\Gamma(P)f(Q,P) -D\partial_Pf(Q,P)=0\,.
\end{equation}
Exploiting the factorization of the equilibrium distribution $f(Q,P) = f_Q(Q) f_P(P)$ and the fact that
\begin{equation}
 f_P(P)\propto e^{-\beta K(P)}
\end{equation}
where $\beta=(k_B T)^{-1}$, one easily finds from Eq.~\eqref{eq:gamma}:
\begin{equation}
\label{eq:einstein_gen}
   \Gamma(P)=-D \beta \partial_P K(P)\,.
\end{equation}
This last equation can be seen as a generalization of the Einstein relation to cases with non-quadratic kinetic energy. It tells that the Stokes law is always proportional to the velocity $\dot{Q}=\partial_P K(P)$, no matter what the form of the kinetic energy is, and that their ratio is fixed by $-D\beta$.
Let us stress that in the usual, Newtonian case $K(P)=P^2/2M$, the Einstein relation is exactly recovered, as it should.
\\
The above argument gives a relation between $D$ and $\Gamma$, but it is not sufficient to determine $\Gamma$ (or $D$) from the knowledge of the Hamiltonian. When the nature of the bath is specified one may try perturbative methods in the limit of large scale separation to derive all the parameters of the effective equation. An analytically tractable case has been discussed in Ref.~\cite{baldovin19}, where the thermal bath is constituted by a large number of Ising spins, which are kept at a fixed temperature by a Glauber dynamics, and the ``slow'' degree of freedom is an oscillator with generalized kinetic energy.
All the spins feel a magnetic field that depends on the position of the oscillator. In the limit in which the typical frequency of the oscillator $\omega_{osc}$ is much slower than the rate of the Glauber dynamics $r_{spin}$, a Chapman-Engsok expansion of the Fokker-Planck equation of the particle can be performed, for which the small parameter is given by
\begin{equation}
\epsilon=\frac{\omega_{osc}}{r_{spin}}\,.
\end{equation}
\\
The obtained  Langevin equation for the slow dynamics is of the form
\begin{equation}
\label{eq:prados_final}
\begin{aligned}
 \dot{Q}&= K'(P)\\
 \dot{P}&=-U'_R(Q) - \Gamma(P,Q) K'(P) +\sqrt{2 D(Q)} \eta
\end{aligned}
\end{equation}
where  $U_R(Q), \Gamma(P,Q)$ and  $D(Q)$ can be explicitly computed.
Remarkably, this result basically coincides with the one obtained with  the previous
phenomenological argument: Eq.~\eqref{eq:einstein_gen} still holds, with the only
difference that $D$ is now a function of $Q$.
\\
It can be verified that both equilibrium features (e.g. stationary probability density function), as well as non-equilibrium  ones
(e.g. correlations and relaxations) obtained with  the effective Langevin equation~\eqref{eq:prados_final}
are in perfect agreement  with the actual numerical  results from the complete system~\cite{baldovin19}.

\section{Testing the generalized LE}
\label{sec:testing}
In order to check the validity of the generalized form for the LE,
\begin{equation}
\label{eq:gl}
 \dot{P}= -\partial_QU(Q) -D \beta \partial_P K(P) + \sqrt{2D} \xi\,,
\end{equation}
we can perform computer simulations of a large, compound system in which both ``heavy'' and ``light'' particles are present, and compare the effective behavior of one slow particle with the stochastic description given by Eq.~\eqref{eq:gl}. Our strategy articulates into four steps:
\begin{enumerate}
 \item Design a suitable Hamiltonian system in which time-scale separation may be expected;
 \item Simulate a (deterministic) evolution of such system;
 \item Extrapolate \textit{a posteriori} the coefficients of the effective LE which approximates the dynamics to the best extent;
 \item Compare them to Eq.~\eqref{eq:gl}.
\end{enumerate}
First, we will carry out our program with a Hamiltonian system which is well-known to reproduce the Brownian motion in the thermodynamic limit, i.e. a harmonic chain with an heavy ``intruder'':
 \begin{equation}
 \label{eq:ham_osc}
 H=\frac{P^2}{2M} + \sum_{i=\pm 1, ..., \pm N} \frac{p_i^2}{2m} +  \frac{k}{2}\sum_{i=-N}^{N+1} (q_i-q_{i-1})^2, \quad \quad Q \equiv q_0\,.
  \end{equation}
  Here $(p_i,q_i),\,i=-N,...,-1,1,...,N$ are the canonical coordinates of the ``light'' particles, with equal masses $m$, while $(P,Q)$ are those of the heavy intruder of mass $M\gg m$; $k$ is the elastic constant. We consider fixed boundary conditions $q_{-N-1}=q_{N+1}=0$.
The above model and similar harmonic chains have been analytically studied since the 1960's and represent one of the few examples in which stochastic differential equations can be exactly derived starting from first principles~\cite{rubin60, turner60, ford-kac-mazur65, zwanzig73}.
\\
Hamiltonian~\eqref{eq:ham_osc} is integrable, so that the energy assigned to each normal mode at the beginning of the dynamical evolution is conserved; as a consequence, if the system is initialized in such a way that energy is shared among only few degrees of freedom, thermodynamic equilibrium will never be reached and the Langevin description~\eqref{eq:LE} will necessarily fail. If, conversely, the system starts at equilibrium, it can be rigorously shown that the dynamics of $P$ is approximated by a Markovian stochastic process, whose autocorrelation function reads
\begin{equation}
C(t) \simeq \exp\left(-\frac{2\sqrt{k m}}{M} t \right) + O(m/M)\,.
\end{equation}
\\
We numerically simulate Hamiltonian~\eqref{eq:ham_osc} with a standard velocity Verlet update, choosing the time-step in such a way that the relative fluctuations on the total energy are of order $O(10^{-5})$. We start from equilibrium initial conditions.
\\
Given a generic Langevin Equation
\begin{equation}
 \dot{P}=F(P)+\sqrt{2 D(P)}\xi
\end{equation}
the drift term $F(P)$ and the diffusivity $D(P)$ can be computed from the temporal evolution of $P$ using the definitions~\cite{gardiner85}
 \begin{subequations}
 \label{eq:d&d}
   \begin{align}
  F(P)&=\lim_{\Delta t \to 0}\, \frac{\langle \Delta P | P(t_0)=P \rangle}{\Delta t}\\
  D(P)&=\lim_{\Delta t \to 0}\, \frac{ \langle \Delta P^2 | P(t_0)=P \rangle}{2 \Delta t}\,.
  \end{align}
 \end{subequations}
 In other words, we can evaluate the Langevin coefficients for a given value $P$ of the variable by looking at the average behavior of the trajectory after it passes through $P$.
This approach has been used in several contexts ranging from physics to biology and finance~\cite{peinke11, peinke18}.
\\
Of course, the above limits of conditioned moments need to be evaluated with care. One has to be sure that the sampling rate is much higher than the typical frequencies of the dynamics, so that the quantities on the r.h.s. of Eq.~\eqref{eq:d&d} can be evaluated for time intervals $\Delta t$ smaller than any characteristic time of the evolution. In this case, in particular, one needs $\Delta t \ll M/\sqrt{k m}$.
\\
On the other hand, the evolution cannot be Markovian on all time-scales, because the original dynamics is deterministic, and it depends on the interactions with other degrees of freedom. This can be also understood by considering the velocity autocorrelation functions $C_D(t)$ and $C_L(t)$ in a deterministic and in a Langevin process respectively; for small times $t$, they can be expanded as:
\begin{subequations}
\begin{align}
 C_D(t)=1- {t^2 \over \tau_D^2}+O(t^3)\\
 C_L(t)=1- {t \over \tau_L}+O(t^2)\,.
\end{align}
\end{subequations}
Therefore it exists a minimal time-scale $\tau_M$ (sometimes called the Markov-Einstein time~\cite{peinke11}), such that the process can be considered Markovian only on time-scales much larger than $\tau_M$. Such threshold should be at least $O\Bigl({\tau_D^2 \over \tau_L}\Bigr)$, in order for the differences between $C_D$ and $C_L$ to be negligible.
\\
At a practical level, a good strategy consists in evaluating the quantities~\eqref{eq:d&d}a and ~\eqref{eq:d&d}b (for a fixed starting value $P$) as functions of the time interval $\Delta t$, then looking at their behavior for $\Delta t \le \tau_M$ (but still small with respect to the typical times of the evolution) and extrapolating the limit $\Delta t \to 0$ (Fig.~\ref{fig:d&d}).
\begin{figure}[!h]
  \includegraphics[keepaspectratio=true, width=0.48\linewidth]{dp_t-eps-converted-to.pdf}
  \includegraphics[keepaspectratio=true, width=0.48\linewidth]{dp2_t-eps-converted-to.pdf}
 \caption{Evaluation of the conditioned moments on the r.h.s. of Eq.~\eqref{eq:d&d}. We numerically compute such quantities as functions of $\Delta t$ (points), then we fit the curves with low-order polynomials (solid lines) and we consider the limits for $\Delta t \to 0$. Left: r.h.s. of Eq.~\eqref{eq:d&d}a, linear fit. Right: r.h.s. of Eq.~\eqref{eq:d&d}b, parabolic fit. Different colors and shapes of the points correspond to different values of the initial value of $P$. All fits are computed between $\Delta t=0.25$ and $\Delta t=1.5$. Parameters of the simulation  $M=200$, $m=1$, $k=2500$, $2N=2000$, $\beta\simeq1.0$. }\label{fig:d&d}
\end{figure}
In order to numerically evaluate the conditioned moments relative to an initial value, say, $P_0$, we have to wait until the trajectory passes through a (small) interval $(P_0-\Delta, P_0+\Delta)$, and then to look at its displacement $\Delta P$ after a time $\Delta t$. This is repeated many times, so that the averages in Eq.~\eqref{eq:d&d} are evaluated as temporal averages.
\\
The above procedure is performed for several values of the starting value of $P$, so that at the end we can appreciate the dependence of the drift and the diffusivity on $P$. The results are shown in Fig.~\ref{fig:fit_osc}.
\begin{figure}[!h]
\centering
   \includegraphics[keepaspectratio=true, width=0.7\linewidth]{fit_osc-eps-converted-to.pdf}
   \caption{Drift term $F(P)$ and diffusivity $D(P)$ of the process $P(t)$ as determined by the data-driven procedure discussed in the text. Red circles (black diamonds) represent the values obtained for the drift (diffusivity) from the limits~\eqref{eq:d&d}; solid lines are linear fits.}\label{fig:fit_osc}
\end{figure}
As expected, the drift linearly depends on the momentum and the diffusivity is constant. Relation~\eqref{eq:einstein_gen} is also verified.
In order to check that the reconstructed LE actually reproduces the behavior of the slow particle, we can do an additional check: we can compute the steady probability density and the autocorrelation function in this new coarse-grained dynamics and compare them to the original, deterministic evolution. In this simple case, since the dynamics is linear, we can determine such observables analytically once we know $F(P)$ and $D$; in more complicated cases one can rely on numerical simulations of the stochastic process. Fig.~\ref{fig:autoc_osc}(a) shows both quantities in the original and in the reconstructed dynamics: the agreement is quite good.
\begin{figure}[!h]
\centering
 \includegraphics[keepaspectratio=true, width=0.48\linewidth]{autoc_osc-eps-converted-to.pdf}
  \includegraphics[keepaspectratio=true, width=0.48\linewidth]{eta_autoc_log-eps-converted-to.pdf}
 \caption{Left: autocorrelation function of the velocity in the original dynamics (circles) and in the reconstructed stochastic process (solid line). In the inset, the steady probability density functions for the deterministic evolution (boxes) and for the stochastic one (solid red line) are compared. Right: autocorrelation function of the velocity (blue) compared to the one of the ``thermal noise'' defined by Eq.~\eqref{eq:zeta} (red).}\label{fig:autoc_osc}
\end{figure}
Finally, we have to check that {time-scale separation} hypothesis is valid, i.e. that the ``thermal noise''
\begin{equation}
\label{eq:zeta}
  \zeta(t)=\dot{P}(t)-F(P(t))
\end{equation}
 decorrelates much faster than  $P$. The autocorrelation functions of the two quantities are shown in Fig.~\ref{fig:autoc_osc}(b): the time-scale separation is evident.

\subsection{A case with also negative temperature}
According to Statistical Mechanics, the inverse temperature $\beta$
is defined  in terms of the microcanonical entropy $S(E)$
of the system described by $H(\mathbf{X})$:
\begin{equation}
\label{eq:beta}
 \beta=\frac{1}{k_B}\frac{\partial S}{\partial E},\quad \quad
 S(E)=k_B \ln\int d{\mathbf{X}}\, \delta(E-H(\mathbf{X})) \, .
\end{equation}
\\
From the above expression one realizes that
 $\beta$ becomes negative if $S(E)$ is a decreasing function. It is well known
 that this never happens for systems with the usual, quadratic kinetic energy~\cite{ramsey56},
 so that negative temperatures can be only observed for peculiar systems,
 typically living in bounded phase-spaces.
\\
Systems with negative temperature, however, are not mere curiosities:
among the many interesting physical cases we can mention
$2D$ inviscid hydrodynamics, including the point-vortices model discussed
by Onsager during the first StatPhys conference (Florence 1949)~\cite{onsager49}, systems of nuclear spins~\cite{ramsey56},
the discrete nonlinear Schr\"odinger equation~\cite{iubini13} and systems of cold atoms trapped in optical lattices~\cite{braun13}.
\\
Since in our derivations we never assumed that the temperature of the system is positive, it can be interesting to see if formula~\eqref{eq:einstein_gen} still holds for Hamiltonian models which can assume negative temperature.
We will consider the following form for the kinetic energy:
\begin{equation}
 K(p)=mc^2\left[1-\cos \left(\frac{p}{mc}\right)\right]\,.
\end{equation}
Here $c$ is a constant with the physical dimensions of a velocity, while $m$ can be seen as a generalized ``mass'': it is straightforward to verify that, at fixed velocity, both kinetic energy and momentum are proportional to $m$, as one would expect from additivity.
Momentum $p$ lives in the interval $[-mc \pi, mc \pi)$, so that it can be considered as an angular variable. If also the ``positions'' are chosen to live in a bounded space, we can expect to observe negative temperature at high values of the energy.
Apart from the constant, the above form of the kinetic energy resembles the one that has been observed in a famous experiment on cold atoms~\cite{braun13}, and it has been used in the definition of mechanical models for systems with negative temperature~\cite{hilbert14, cerino15, baldovin17, miceli19}. One of the reasons is that for small energies such term reduces to the usual, quadratic form $p^2/2m$, so that in the low (positive) temperature regime we recover the usual statistical properties.
\\
Also in this case, we want to study the effective motion of a ``heavy'' particle subjected to the action of a thermal bath of ``light'' molecules. We need that the intruder interacts simultaneously with many, uncorrelated, degrees of freedom, in such a way that it can be considered at thermal equilibrium with such reservoir. In the following, we will model the bath as a chain of ``oscillators'' with equal masses $m$:
 \begin{equation}
  H_{chain}=\sum_{i=1}^{N}mc^2\left[1-\cos \left(\frac{p_i}{mc}\right)\right]+\epsilon \sum_{i=1}^{N+1}[1-\cos(\theta_i-\theta_{i-1})]
 \end{equation}
where $(p_i,\theta_i)$ is the $i-$th pair of canonical coordinates.
We coupled then such chain to a ``heavy'' degree of freedom $(\Theta,P)$ of mass $M\gg m$ via a bounded potential, so that the total Hamiltonian reads
  \begin{equation}
  H=H_{chain}+Mc^2[1-\cos(P/cM)] + k \sum_{i=1}^{N/n} [1-\cos(\Theta-\theta_{i\cdot n})]\,.
 \end{equation}
Let us stress that the heavy particle is only coupled to those oscillators whose label is a multiple of the integer parameter $n\gg 1$.
\\
We can expect that the different inertia of the intruder makes its dynamics much slower than that of the other degrees of freedom, so that a time-scale separation should be observed. Moreover, our choice of the interaction potential should keep the particle at the same temperature of the thermal bath. Remarkably enough, in this case such temperature can be negative, due to the bound on the total phase space volume. If the generalized Einstein relation~\eqref{eq:einstein_gen} holds for $\beta<0$, it means that $F(P)$, in such regime, must be positive when $P>0$ and negative when $P<0$, contrary to what happens at positive temperature. In other words, the (generalized) Stokes force tends to give energy to the particle, instead of subtracting it. At first sight this behavior may appear very unphysical: it can be understood by remembering that a thermal bath at negative temperature increases its entropy by decreasing the internal energy, i.e. by releasing heat.
\\
We can apply to this new model the data-driven analysis discussed in the previous section.
\begin{figure}
 \centering
  \includegraphics[keepaspectratio=true, width=0.49\linewidth]{fit_pos-eps-converted-to.pdf}
  \includegraphics[keepaspectratio=true, width=0.49\linewidth]{fit_neg-eps-converted-to.pdf}
  \caption{Drift and diffusivity evaluated with the discussed data-driven approach (points) and corresponding fits (solid lines). Left: $\beta \simeq +0.11$. Right: $\beta \simeq -0.10$. Drift terms are fitted with sinusoidal functions, diffusivities with constant values. Parameters of the simulations: $M=8$, $m=1$, $k=0.5$, $N=600$, $n=15$. }\label{fig:posneg}
\end{figure}
 Fig.~\ref{fig:posneg} shows the drift and diffusivity terms as reconstructed by our approach. We considered two cases, one at positive and one at negative $\beta$.
 The drift term $F(P)$ is clearly proportional to $\partial_P K(P)$, as it should. The diffusivity is almost constant. Relation~\eqref{eq:einstein_gen} holds both at positive and at negative temperature. In particular, as it is clear from the figure, passing from positive to negative $\beta$ the proportionality factor between $F(P)$ and $K(P)$ changes its sign.
\begin{figure}
  \centering
  \includegraphics[keepaspectratio=true, width=0.48\linewidth]{autoc_pos-eps-converted-to.pdf}
   \includegraphics[keepaspectratio=true, width=0.48\linewidth]{autoc_neg-eps-converted-to.pdf}
   \caption{Autocorrelation functions of the velocity in the original dynamics (circles) and in the coarse-grained model (solid lines) at $\beta=2$ (left) and $\beta=-2$ (right). In the insets, probability density functions of $P$ in the two cases, both for the original dynamics (boxes) and for the reconstructed model (solid lines).} \label{fig:autoc}
\end{figure}
\\
Finally, we can look at the steady probability density functions and at the autocorrelations of $P$, in order to compare the ones computed with the coarse-grained model and those from the original deterministic system. As shown in Fig.~\ref{fig:autoc}, the agreement is quite good, both at positive and negative $\beta$.
\\
The fact that it is possible to write down a Langevin-like equation with $\beta<0$, which properly reproduces the behavior of a slow degree of freedom subjected to the action of a bath, sounds quite relevant in the context of the long-lasting debate about negative temperatures. Several Authors propose the adoption of a definition of entropy alternative to Eq.~\eqref{eq:beta}, the so-called ``volume entropy'' or ``Gibbs entropy'' $S_G$~\cite{dunkel14}. A statistical description based on $S_G$ is able to reproduce the classical results of thermodynamics \textit{exactly}, even if the number of degrees of freedom of the system is small~\cite{romero13, dunkel14, hilbert14, campisi15}: due to this remarkable property, $S_G$ may appear as the right mechanical analogue of the thermodynamic entropy. However, since $S_G$ is by definition a non-decreasing function of energy, $\beta$ would never be negative. We stress that with this alternative choice it would not be possible, in our opinion, to give a coherent description of the above results.

\subsection{ Some remarks on the method}
At a first glance it seems that everything is quite easy.
On the other hand one has to face several difficulties:
\begin{enumerate}
 \item Time-scales separation  does not imply Markovianity.
 \item It is true that, in some cases, adding variables can
lead to a Markovian process
(e.g. if colored noise is present); unfortunately
 a general method to find the ``right'' variables does not exist.
Such a difficulty has been expressed in a rather vivid way
in the caveat of Onsager and Machlup~\cite{onsager53}:
   {\it How do you know you have taken enough variables,
for it to be Markovian?}
\item The  procedure discussed in this Section cannot be applied ``blindly'',
and Markovianity should always be checked a posteriori.
\end{enumerate}
In order to give an example of the above troubles and
how mathematics can help to select/eliminate models
we discuss the case shown in Fig.~\ref{fig:contro}, which is obtained from a system with Hamiltonian
\begin{equation}
\label{eq:fakeH}
 H=\sum_{i=-N}^{N} [1-\cos(p_i/\sqrt{m_i})] + k\sum_{i=-N}^{N+1}[1-\cos(q_{i}-q_{i-1})]
\end{equation}
where $(p_0,q_0)\equiv(P,Q)$ and $m_i=m+\delta_{i,0}(M-m)$.
\iffalse
\begin{equation}
 m_i=
 \begin{cases}
  M \quad\mbox{ if } i=0\,,\\ m\quad\mbox{ otherwise. }
 \end{cases}
\end{equation}
\fi
\\
Also in this case one can expect that the presence of a large ``mass'' would lead to a scale-separation between the dynamics of the slow and those of the fast particles. However let us note that the additivity condition mentioned above is not satisfied in this case.
\begin{figure}
  \centering
  \includegraphics[keepaspectratio=true, width=0.5\linewidth]{autoc_contro-eps-converted-to.pdf}
   \caption{A case in which the described method fails: the autocorrelation function of the original system~\eqref{eq:fakeH} (circles) and that of the reconstructed model (solid red line) do not coincide. Similarly, in the inset, the original (boxes) and the reconstructed p.d.f. (solid red line) are different. Paramters: $M=200$, $m=1$, $k=2500$, $2N=2000$, $\beta=1.04$.} \label{fig:contro}
\end{figure}
In such a situation, in spite of
time-scales separation, it is impossible to
find a $1D$ Langevin equation for $P$.
The reason of such impossibility is  a property of the Fokker-Planck operator: in $1D$ equilibrium Langevin equation the autocorrelation function
cannot be negative~\cite{gardiner85}. Therefore if we want to build an effective Langevin equation,
it is necessary to introduce (at least) another variable.

\section{Data driven models: LE for a granular system}
\label{sec:granular}

Let us now discuss how to use the  procedure  for the building of effective
equations when just  experimental data are available,
in  a system  for which we have not  a well established
theoretical understanding.
\\
The analyzed time series is that of the angular velocity of a rotator suspended
 in a vibrofluidized granular medium, as discussed in Ref.~\cite{scalliet15}. In Fig. \ref{fig:scheme} we show a sketch of the experiment set-up. We recall that in a granular medium kinetic energy is not conserved and therefore Hamiltonian modelling discussed before is not applicable.
\begin{figure}
  \centering
  \includegraphics[keepaspectratio=true, width=0.48\linewidth]{E1-eps-converted-to.pdf}
      \caption{Sketch of the set-up for the experiment discussed in Ref.~\cite{scalliet15}. The system is composed of a large number of beads immersed in a cylindrical container, which is shaken from below with a fixed frequency. A blade is immersed in the granular gas and it is free to rotate along the vertical axis. Its angular position is measured by an encoder. The acceleration felt by the container is also measured.} \label{fig:scheme}
\end{figure}
\\
The granular medium is composed of $N$ spheres of diameter $d=4$ mm,
placed in a cylindrical container whose volume is $\sim$ 7300 times that of a single sphere.
The container is vertically shaken with a signal whose spectrum is approximately flat in a range
between  $f_{min}= 200$ Hz and $f_{max}= 400$ Hz.
A blade, a "massive tracer", with cross section $\sim 8 d \times  4 d$, is suspended into the granular
medium and rotates around a vertical axis.
Its angular velocity $\omega(t)$ and the traveled angle of rotation $\theta(t)=\int_0^t \omega(t')dt'$
are measured with a time-resolution of $2$ kHz.

\subsection{The Gas limit}

In the dilute limit, e.g. with $N=350$, and
 packing fraction  $\sim 5 \%$, the simplest possible scenario holds.
The blade performs  a standard rotational Brownian motion, and a $1D$
linear Langevin equation for $\omega$ is enough for a good description of the observed behavior.
We have
$$
{d \omega \over dt}= -{\omega \over \tau} + c \eta\,,
$$
where the parameters $\tau$ and $c$ can be easily obtained with the procedure
discussed in Sec.~\ref{sec:testing}, see Fig.~\ref{fig:gran_dilute}. This analysis has been done in Ref.~\cite{baldogran}.
\begin{figure}
  \centering
  \includegraphics[keepaspectratio=true, width=1.\linewidth]{Fig4-eps-converted-to.pdf}
     \caption{Velocity power density spectrum (left) and mean square displacement (right) for the angular position of the blade in the dilute-gas limit. Squares represent experimental data, solid lines are computed with the reconstructed model.} \label{fig:gran_dilute}
\end{figure}

\subsection{The Cold Liquid limit}
On the contrary, in the dense regime, e.g.
with $N=2600$ and  packing fraction  $ \sim 30 \%$, the
scenario of the usual standard rotational Brownian motion
fails, and a rich phenomenology appears (e.g., a superdiffusive behavior at long time-scales).
A single  equation for $\omega$ is not able to describe in a proper way the observed results.
It is necessary to introduce, at least,  a second variable.
These new variable, which  describes the slow behavior of the probe,
has been obtained by performing a running average with a Gaussian window function
$$
\theta_0(t)={1 \over \sqrt{2 \pi \sigma^2}} \int e^{-{(t-t')^2 \over 2 \sigma^2}} \theta(t') \, dt'
$$
the range of values in which $\sigma$ should be chosen is suggested by the shape of $S(f)$.
\\
The reconstructed Langevin equation for the variables
$$
 \theta \,, \,  \theta_0 \,, \,  \omega={d\theta \over dt} \,, \,  \omega_0={d\theta_0 \over dt}\,,
$$
discussed in Ref.~\cite{baldogran}, is able to reproduce the main statistical features, including the anomalous
diffusion, see Fig~\ref{fig:gran_dense}.
\\
\\
\begin{figure}
  \centering
  \includegraphics[keepaspectratio=true, width=1.\linewidth]{Fig6-eps-converted-to.pdf}
   \caption{Velocity power density spectrum (left) and mean square displacement (right) for the angular position of the blade in the cold-liquid limit. Squares represent experimental data, solid lines are computed with the reconstructed model. Different colors of the lines correspond to differnt values of the parameter $\sigma$ used for the analysis, as mentioned in the text.} \label{fig:gran_dense}
\end{figure}
\\
\\
\section{ About enthusiasm for Big Data,  Artificial Intelligence
and Machine Learning }
\label{sec:bigdata}
Recently the possibility of extracting knowledge by data mining (i.e.,
through the algorithmic analysis of large amounts of data) seems to
suggest the emergence of a fourth paradigm, a new scientific
methodology.
In particular a consistent fraction of the scientific community seems close
 to conclude that it is possible to
build models with  black box protocols.
\\
On the other hand, in our opinion, the past experience show that
for the problem of the prediction, there are rather severe
limitations for  such an approach.
For instance in Refs.~\cite{cecconi12,hosni18} the reader can find a detailed discussion
about the
good reasons (mainly due to  the Kac's lemma on the recurrence time) to be
skeptical about the data-centric enthusiasm supporting a general
philosophy starting from "raw data", without constructing modeling
hypotheses and, therefore, without theory.
\\
Let us discuss again how the ability of choosing the "right"
variables typically requires a conceptual abstraction which is key
to scientific discoveries.
\\
In the '80s, some researchers in the field of artificial intelligence
(AI) devised BACON, a computer program to automatize scientific
discoveries.  Apparently, BACON was able to "discover" the third
Kepler's law~\cite{bradshaw71,grabiner86}.
\\ Let us look at the details of the procedure used by
BACON, and by Kepler:
\begin{itemize}
 \item BACON used as input the numerical values
of distance from the Sun, $D$, and revolution period, $P$, of
planets. The program, then, discovered that $D^3$ is proportional to
$P^2$.
\item For Kepler the raw observables were not $D$ and $P$, but
a list of planetary positions seen from the Earth at different
times. In his discovery, Kepler chose the "right" variables $D$ and
$P$ as he was guided by strong beliefs in mathematical harmonies as
well as the (at that time) controversial theory of Copernicus.
\end{itemize}
Is it appropriate to claim that  AI methods can replace the traditional creative
approach to scientific discoveries~\cite{buchanan19}?

\section{Summary and conclusions}
\label{sec:conclusion}

Let us briefly summarize the main points discussed in this paper,
and try to fix some conclusions.
\begin{enumerate}
 \item LE can be generalized to cases with ``unusual'' kinetic energy, and
 in particular to systems which can show negative temperature.
\item In order to check its validity, one can simulate a  bath  and find the effective LE,
 this approach works in presence of:
 \begin{itemize}
  \item time-scale separation;
  \item Markovianity of the variable we are looking at.
 \end{itemize}
\item A procedure to find LE for the slow dynamics can be used in the treatment of data if
 one is able to identify the relevant variables.
\item The choice of the "good variables" does not follow from a mechanical protocol,
 it can be only suggested by intuition and/or a preliminary understanding of
 the main aspects of the phenomena under investigation.
\end{enumerate}
We conclude with a remark about the machine learning methods (MLM) and AI in the research.
It is matter of fact that in several cases MLM and AI had been able to succeed;
we believe that  community  have a great opportunity to give a contribute in the
understanding of  the range of applicability of MLM and AI.

\begin{acknowledgments}
 {We aknowldge fruitful collaboration and useful discussions with F. Borra, F. Cecconi, M. Cencini, M. Falcioni, A. Gnoli, A. Prados, A. Puglisi and A. Sarracino.}
\end{acknowledgments}

\section*{Bibliography}

\bibliographystyle{unsrt}


\section{Introduction} \label{sec:intro}

\input{1-KFJLT-guarantees-intro-arxiv-v2.tex}

\section{Other Related Work} \label{sec:related-work}

\input{2-KFJLT-guarantees-related-work-arxiv-v2.tex}

\section{Preliminaries} \label{sec:preliminaries}

\input{3-KFJLT-guarantees-prerequisites-arxiv-v2.tex}

\section{Main Results} \label{sec:main-results}

\input{4-KFJLT-guarantees-main-result-arxiv-v2.tex}

\section{Numerical Experiments} \label{eq:experiments}

\input{5-KFJLT-guarantees-experiments-arxiv-v2.tex}

\section{Conclusion} \label{sec:conclusion}

\input{6-KFJLT-guarantees-conclusions-arxiv-v2.tex}

\section*{Acknowledgments}

We would like to thank Tammy Kolda, Ruhui Jin and Rachel Ward for providing feedback on an early version of this work. We also thank the anonymous reviewers for their comments which helped improve the paper.

This material is based upon work supported by the National Science Foundation under Grant No.\ ECCS-1810314.

\bibliographystyle{plainnat}

\subsection{Experiment 1: Synthetic Data}

In this section, we present the results from an experiment which compares five different sketches when applied to random Kronecker vectors with three different random distributions. The five methods we compare are the following.
\begin{itemize}
	\item \textbf{Gaussian} sketch uses an unstructured $J \times \tilde{I}$ matrix with i.i.d.\ standard normal entries which are scaled by $1/\sqrt{J}$. This approach is not scalable, but interesting to use as a baseline in this experiment.
	\item \textbf{KFJLT} is the sketch discussed in this paper, and which is defined in Definition~\ref{def:KFJLT}, with the only difference that the uniform sampling of rows is done without replacement.
	\item \textbf{TRP} is the method proposed in \citep{sun2018}. As sub-matrices, we use matrices with i.i.d.\ standard normal entries of size $I_p \times J$ and rescale appropriately.
	\item \textbf{TensorSketch} is the method developed in \citep{pagh2013, pham2013, avron2014, diao2018}.
	\item \textbf{Sampling} is the method proposed by \citep{cheng2016}. It computes an estimate of the leverage scores for each row of the matrix to be sampled and uses these to compute a distribution which is used for sampling.
\end{itemize}
All of these methods, except the first, are specifically designed for sketching vectors with Kronecker structure. As input, we use random Kronecker vectors $\xbf = \kron_{p=1}^3 \xbf^{(p)}$, where each $\xbf^{(p)} \in \Rb^{16}$ has one of the following distributions.
\begin{itemize}
	\item Each $\xbf^{(p)}$ has i.i.d.\ standard normal entries.
	\item Each $\xbf^{(p)}$ is sparse with three nonzero elements, which are independent and normally distributed with mean zero and standard deviation 100. The positions of the three nonzero elements are drawn uniformly at random without replacement.
	\item Each $\xbf^{(p)}$ contains a single nonzero entry, which is chosen uniformly at random. This nonzero entry is equal to 100.
\end{itemize}
Sparse Kronecker vectors are interesting in many data science applications, and arise in decomposition of sparse tensors, for example.

In the experiment, we draw two random vectors $\xbf, \ybf \in \Rb^{4096}$, each of which is a Kronecker product of three smaller vectors of length $16$, drawn according to one of the three distributions above. For each of the five sketches $f$ and for some embedding dimension $J$, we then compute how well they preserve the distance between $\xbf$ and $\ybf$ by computing the quantity
\begin{equation} \label{eq:distortion-metric}
	\Big| \frac{\|f(\xbf)-f(\ybf)\|_2}{\|\xbf - \ybf\|_2} - 1 \Big|.
\end{equation}
For each of the three distributions and each embedding dimension
\begin{equation} \label{eq:J-range}
	J \in \{100, 200, \ldots, 1000\}
\end{equation}
we repeat this 1000 times and compute the mean, standard deviation and maximum of \eqref{eq:distortion-metric} over those 1000 trials.
For $J$ as in \eqref{eq:J-range}, applying one of the sketches to $\xbf$ and $\ybf$ reduces the number of entries in those vectors by between 76\% (for $J=1000$) and 98\% (for $J=100$).
Figures~\ref{fig:experiment1-normal}--\ref{fig:experiment1-large-single} present the results for each of the three distributions.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=1\textwidth]{experiment1-normal.pdf}
	\caption{Mean, standard deviation, and maximum of the quantity in \eqref{eq:distortion-metric} over 1000 trials when the test vectors are Kronecker products of vectors with i.i.d.\ standard normal entries.}
	\label{fig:experiment1-normal}
\end{figure}

\begin{figure}[ht!]
	\centering
	\includegraphics[width=1\textwidth]{experiment1-sparse.pdf}
	\caption{Mean, standard deviation, and maximum of the quantity in \eqref{eq:distortion-metric} over 1000 trials when the test vectors are Kronecker products of vectors with three nonzero elements which are independent and normally distributed with mean zero and standard deviation 100.}
	\label{fig:experiment1-sparse}
\end{figure}

\begin{figure}[ht!]
	\centering
	\includegraphics[width=1\textwidth]{experiment1-large-single.pdf}
	\caption{Mean, standard deviation, and maximum of the quantity in \eqref{eq:distortion-metric} over 1000 trials when the test vectors are Kronecker products of vectors containing a single nonzero entry equal to 100.}
	\label{fig:experiment1-large-single}
\end{figure}

No one method produces the best results for all three distributions. The leverage score sampling approach does very well on dense vectors, even outperforming the Gaussian sketch, but does less well on sparser inputs. Although TensorSketch has an impressive mean performance on the two sparser inputs, it sometimes produces high distortion rates on those inputs. On the sparser inputs, the KFJLT seems to strike the best balance between mean and worst case performance. TRP does poorly for all three distribution types.

\subsection{Experiment 2: MNIST Handwritten Digits}

In this experiment, we consider a subset of the MNIST Handwritten Digits dataset \citep{lecun1998}, which is a standard benchmark dataset in machine learning.\footnote{We downloaded the MNIST dataset using the scripts provided at \url{https://github.com/sunsided/mnist-matlab}.}
The dataset consists of images of handwritten digits between 0 and 9.
Each image is in gray scale and of size 28 by 28 pixels.
To make the image width and height powers of two, we pad the images with zeros so that their size is 32 by 32 pixels.
We arrange 100 images depicting fours into a tensor $\Xe \in \Rb^{32 \times 32 \times 100}$ and 100 images depicting nines into another tensor $\Ye \in \Rb^{32 \times 32 \times 100}$.
Handwritten fours and nines can look quite similar and can be difficult to distinguish, which is why we choose this particular pair of digits.
We then compute a rank-10 approximate CP decomposition of each tensor using \verb|cp_als| in the Tensor Toolbox for Matlab \citep{bader2006,bader2015}. These take the form
\begin{equation}
\begin{aligned}
	\hat{\Xe} &\defeq \sum_{r=1}^{10} \Abf^{(1)}_{:r} \circ \Abf^{(2)}_{:r} \circ \Abf^{(3)}_{:r} \approx \Xe, \\
	\hat{\Ye} &\defeq \sum_{r=1}^{10} \Bbf^{(1)}_{:r} \circ \Bbf^{(2)}_{:r} \circ \Bbf^{(3)}_{:r} \approx \Ye, \\
\end{aligned}
\end{equation}
where $\circ$ denotes outer product, and each $\Abf^{(1)}, \Bbf^{(1)}, \Abf^{(2)}, \Bbf^{(2)} \in \Rb^{32 \times 10}$ and $\Abf^{(3)}, \Bbf^{(3)} \in \Rb^{100 \times 10}$ are called factor matrices; see \citep{kolda2009} for further details on tensor decomposition.
Notice that the factor matrices require much less storage than the original tensors. Figure~\ref{fig:4-9} shows an example of a four and a nine in the MNIST dataset, and their corresponding approximations in the CP tensors.

\begin{figure}[ht!]
	\centering
	\includegraphics[width=1\textwidth]{4-9-plot.pdf}
	\caption{Example of a four and a nine with their corresponding approximations.}
	\label{fig:4-9}
\end{figure}

In this experiment, we apply the KFJLT, TRP, TensorSketch and sampling sketches to $\hat{\Xe}$ and $\hat{\Ye}$ to see how well they preserve the distance between these tensors.
The four sketches can be applied efficiently to $\hat{\Xe}$ and $\hat{\Ye}$ in their decomposed form.
To see this, define $\Xbf \defeq \Abf^{(1)} \odot \Abf^{(2)} \odot \Abf^{(3)} \in \Rb^{102400 \times 10}$ and $\Ybf \defeq \Bbf^{(1)} \odot \Bbf^{(2)} \odot \Bbf^{(3)} \in \Rb^{102400 \times 10}$, and let $\ubf \in \Rb^{20}$ denote a column vector with elements $\ubf_i = 1$ if $1 \leq i \leq 10$ and $\ubf_i = -1$ if $11 \leq i \leq 20$.
Then the following relation holds:
\begin{equation}
	\|\hat{\Xe} - \hat{\Ye}\|_\F = \| [ \Xbf, \, \Ybf ] \, \ubf \|_2,
\end{equation}
where $\|\cdot\|_\F$ denotes the tensor Frobenius norm.
Since the columns of $[ \Xbf, \, \Ybf ]$ are Kronecker products, each of the four sketches under consideration can be applied efficiently to this matrix.
The Gaussian sketch requires too much memory and is therefore not considered.

For any sketch $\Mbf$ with the property
\begin{equation} \label{eq:M-tensor-sketch}
	\| [ \Xbf, \, \Ybf ] \, \ubf \|_2 \approx \| \Mbf \, [ \Xbf, \, \Ybf ] \, \ubf \|_2,
\end{equation}
we may use $\| \Mbf \, [ \Xbf, \, \Ybf ] \, \ubf \|_2$ as an estimate for $\|\hat{\Xe} - \hat{\Ye}\|_\F$.
In the case of KFJLT, a guarantee of the form \eqref{eq:M-tensor-sketch} follows from Theorem~\ref{thm:subspace-embedding} when $J$ is large enough.
For each of the four sketches and some embedding dimension $J$, we compute the quantity
\begin{equation} \label{eq:distortion-metric-2}
\Big| \frac{\| \Mbf \, [ \Xbf, \, \Ybf ] \, \ubf \|_2}{\|\hat{\Xe} - \hat{\Ye}\|_\F} - 1 \Big|
\end{equation}
as a measure of performance.
For each embedding dimension
\begin{equation} \label{eq:J-range-2}
	J \in \{ 100, 200, \cdots, 5000 \}
\end{equation}
we repeat this 1000 times and compute the mean, standard deviation and maximum of \eqref{eq:distortion-metric-2} over those 1000 trials.
The pair $(\hat{\Xe}, \hat{\Ye})$ remains the same in all trials.
For $J$ as in \eqref{eq:J-range-2}, applying one of the sketches to $[ \Xbf, \, \Ybf ]$ reduces the number of rows by between 95\% (for $J = 5000$) and 99.9\% (for $J = 100$).
Figure~\ref{fig:experiment2} presents the results of the experiment.
\begin{figure}[ht!]
	\centering
	\includegraphics[width=1\textwidth]{experiment2.pdf}
	\caption{Mean, standard deviation, and maximum of the quantity in \eqref{eq:distortion-metric-2} over 1000 trials.}
	\label{fig:experiment2}
\end{figure}
All methods have similar performance.
The fact that the factor matrices are mostly dense (some rows are zero due to the padding) may explain why the occasional large errors for TensorSketch observed in Figure~\ref{fig:experiment1-large-single} are avoided.
Our results here indicate that the different sketches may have more similar performance on vectors like $[ \Xbf, \, \Ybf ] \, \ubf$ which have less structure than the Kronecker vectors in the synthetic experiment.
\citet{jin2019} made the related observation that KFJLT does a better job of embedding unstructured vectors than Kronecker structured ones; see Section~5.2 in their paper for further details.