\section{Introduction}
\label{sec:introduction}
This chapter is a revision of the lecture and the related hand-out
which I presented to the École de Physique des Houches on July 26-31
2020. Due to its strictly tutorial character, I shall not give
detailed primary references. I shall mention some references that
expand and support specific points, and I shall add some in a final
section to point to further developments.

I aim to review the basics of a peculiar setting for Information
Geometry (IG) in the sense of \cite{amari|nagaoka:2000}. This setting
has the following peculiarities.

\begin{itemize} \item It is non-parametric and infinite-dimensional.
\item It provides an affine manifold modeled on a Banach space in the sense of \cite{lang:1995}, the Banach space being an Orlicz space as defined \cite[Ch. 8]{adams|fournier:2003}.
\item  It focuses on a particular expression of the tangent bundle, called Statistical Bundle (ST).
\item It allows for the use of (weakly) differentiable densities when the reference measure is Gaussian as in \cite[Ch. V]{malliavin:1995}.
\end{itemize}

A previous tutorial paper \cite{pistone:2020NPCS} presents this
non-parametric construction in the case of finite state space. A
general presentation of the finite case is in
\cite[Ch. 2]{Ay|Jost|Le|Schwachhofer:2017IGbook}. Here, I will focus
on the preliminaries of the infinite state space case. A comprehensive
modern introduction to the whole topic is found in
\cite{nielsen|2020-entropy-1100}.

There are many other successful presentations of IG which are indeed
non-parametric. I will give a few relevant references in the concluding
section. Any useful presentation should explain and
include its historical development elements
below $\bm 1$ to $\bm 5$.

\noindent $\bm 1.$ The starting point to consider is the now classical
work of R. Fisher, see, for example,
\cite[Ch. 4]{efron|hastie:2016}. A \emph{regular statistical model} is
a mapping from a set of parameters $\Theta$ which is an open domain of
$\mathbb{R}^d$ to probability densities on a given measured sample space
$\mathcal P(X,\mathcal X,\mu)$, $\theta \mapsto p(\theta)$, such that
the following computation is feasible. If $f$ is a given random
variable, one wants to compute the variation of the expectation
(assuming its existence), as
\begin{multline*}
  \pderivby {\theta_i} \expectat {p(\theta)} f = \pderivby {\theta_i}
  \int f\ p(\theta) \ d\mu = \pderivby {\theta_i} \scalarat \mu f
  {p(\theta)} = \\
  \scalarat \mu f {\pderivby {\theta_i} p(\theta)} = \scalarat \mu f
  {\frac{\pderivby {\theta_i} p(\theta)}{p(\theta)} p(\theta)} = \scalarat {p(\theta)} f
  {\pderivby {\theta_i} \log p(\theta)} = \\
  \scalarat {p(\theta)} {f - \expectat {p(\theta)} f }
  {\pderivby {\theta_i} \log p(\theta)} \ .
\end{multline*}

The random vector with components $\pderivby {\theta_i} \log
p(\theta)$ is the \emph{Fisher score} of the model al $\theta$. Its
expected value with respect to $p(\theta)$ is 0 and its
variance matrix with respect to $p(\theta)$ is the \emph{Fisher
information matrix},
\begin{equation}\label{eq:fisher}
  I(\theta) = \left[ \expectat {p(\theta)} {\pderivby {\theta_i} \log
      p(\theta) \pderivby {\theta_j} \log p(\theta)} \right]_{ij} =
  \left[ \int \frac{\pderivby {\theta_i}  p(\theta) \pderivby
      {\theta_j}  p(\theta)}{p(\theta)} \ d\mu \right]_{ij} \ .
\end{equation}

The computations above shows many peculiar  features:
\begin{itemize} \item  One computes the variation of the expectation
  as a function of the statistical model. Moreover, a moving inner
  product $\theta \mapsto \scalarat {p(\theta)} \cdot \cdot$ appears naturally.
\item The score $\partial_i \log p(\theta)$ represents the velocity of
  variation of the statistical model, while $f - \expectat {p(\theta)}
  f$ represents the gradient of the expectation function.
\item  The velocity at $\theta$ lives in the space of random variables
  centred at $p(\theta)$.
\item The information matrix provides squared norms and scalar
  products of the velocities in the moving inner product.
\end{itemize}

\noindent$\bm 2.$ One explanation of the Fisher computations results
from the assumption of an exponential model,
\begin{equation*}
  p(\theta) = \mathrm{e}^{\sum_i \theta_i u_i - \kappa(\theta)} \ ,
\end{equation*}
where $u_i$ are the sufficient statistics and $\kappa(\theta)$ is the
cumulant of $\sum_i \theta_i u_i$, see, for example,
\cite[Ch. 5]{efron|hastie:2016}. In such a case, the score is the
centered sufficient statistics,
\begin{equation*}
  \partial _i \log p(\theta) = u_i - \partial_i
  \kappa(\theta) = u_i - \expectat {p(\theta)} {u_i} \ ,
\end{equation*}
and the information matrix equals the Hessian of the convex function
$\theta \mapsto \kappa(\theta)$,
\begin{equation*}
  I(\theta) = \left[\partial_i \partial_j \kappa(\theta)\right]_{ij} \ .
\end{equation*}
See \cite[Ch. 2]{brown:86} for a full account of analytic properties
of exponential families.

\noindent$\bm 3.$ C. R. Rao has been the first statistician to remark
that the Fisher information matrix is positive definite and smooth in
an adequately defined regular model. Hence it defines a Riemannian
metric on the space of parameters. Moreover, he provided an embedding
argument for the resulting manifold. What possibly inspired him in his
construction was Differential Geometry use in Physics. The original
presentation of IG in \cite{amari|nagaoka:2000} can be considered the
full unfolding of this approach.

The mapping
\begin{equation*}
  \Theta \ni \theta \mapsto 2 \sqrt{p(\theta)} = P(\theta)
\end{equation*}
maps the parameters' space $\Theta$ into the $L^2(\mu)$-sphere of
radius $2$. The vectors
\begin{equation*}
  \partial_i P(\theta) = \partial_i 2 \sqrt{p(\theta)} =
  \frac{\partial_i p(\theta)}{\sqrt{p(\theta)}}
\end{equation*}
are in the tangent space at $P(\theta)$ of the sphere, and the inner
product between tangent vector is
\begin{equation*}
  \int \partial_i P(\theta) \partial_j P(\theta) \ d\mu = \int
  \frac{\partial_i p(\theta) \partial_j(x,\theta)}{p(x;\theta)} \ d\mu \ ,
\end{equation*}
that is, the $(i,j)$ element of the Fisher information matrix.

The Rao's computations above reproduce all the metric structure of the
Fisher computations but in one point. That is, now the velocity is not
expressed by the logarithmic derivative $\partial_i \log p(\theta) =
\partial_i p(\theta) / p(\theta)$, but it is expressed by $\partial_i
2 \sqrt{p(\theta)} = \partial_i p(\theta) / \sqrt{p(\theta)}$.

It is possible to make sense, at least formally, of the apparent
contradiction by considering that there are here three different
expressions of the same object.
\begin{itemize}
\item $T\mathcal P_>$ is the tangent bundle of the set of positive
  densities $\mathcal P_>$; the tangent vectors $\dot p$ satisfy
  $\int \dot p \ d\mu = 0$.
\item $TS_2$ is the tangent bundle of the $L^2(\mu)$ sphere $S_2$.
\item $S \mathcal P_>$ is the Fisher's statistical bundle consisting
  of all couples $(p.u)$ such that $p \in \mathcal P_>$ and
  $\expectat p u = 0$.
\end{itemize}

The Rao's embedding $p \mapsto 2 \sqrt p$ provides the identification
of $T\mathcal P$ with $TS_2$. In fact, the computation of the tangent
mapping of $s \colon P \mapsto P^2/4 = p$ gives $ds(P)[\dot P] = \sqrt
p \dot P$.

The identification of $TS_2$ with
$S\mathcal P$ is provided by
\begin{equation*}
  TS_2 \ni (P,\dot P) \mapsto \left(\frac14 P^2, 2 \frac {\dot P}
    P\right) = (p,u) \in S\mathcal P \ .
\end{equation*}
In fact,
\begin{gather*}
  \int p \ d\mu = \frac14 \int P^2 \ d\mu = 1 \ ,\\
  \int u \ p \ d\mu = \frac 12 \int \frac {\dot P}{P} P^2 \ d\mu = \frac
  12 \int P \dot P \ d\mu = 0 \ ,\\
  \int u_1 u_2 \ p \ d\mu = \int \frac{\dot P_1}{P} \frac{\dot P_2}{P}
  P^2 \ d\mu = \int \dot P_1 \dot P_2 \ d\mu \ .
\end{gather*}

A large part of the literature in IG uses the expressions
$T\mathcal P$ and $TS_2$. Still, my own choice is to use
$S\mathcal P$.  It fits well with the statistical picture and the
exponential representation of strictly positive densities.

The choice of the exponential expression and Fisher's score might seem
arbitrary, but the following argument shows is not, see
\cite[Ch. 3]{Ay|Jost|Le|Schwachhofer:2017IGbook}. Assume
$t \mapsto \mu(t)$ is a one-dimensional smooth model of probability
measures and assume the mapping $t \mapsto \mu(A;t)$ is smooth for all
measurable $A$. Fix a value $\bar t$ of the parameter $t$. If a
measurable set $A$ is a zero set for $\mu(\bar t)$, then
$t \mapsto \mu(A,t)$ is minimum at $t = \bar t$ then the derivative is
zero at $\bar t$, namely $\dot \mu(\bar t)=0$. Il follows that the
measure $\dot \mu(\bar t)$ is absolutely continuous with respect to
$\mu(\bar t)$. The resulting logarithmic derivative
$d \dot \mu(\bar t) / d \mu(\bar t)$ is clearly a generalization of
the Fisher's score.

\noindent$\bm 4.$ The Riemannian approach by C. R. Rao can lead to a
more in-depth study of the second-order properties of the manifold,
namely the Levi-Civita connection and the curvature. Several authors,
notably S-I. Amari, B. Efron, Ph. Dawid, S. Lauritzen have promoted a more general perspective in studying the geometry of statistical
models. In modern terminology, a
statistical manifold consists of a metric and a couple of flat
connections,  in duality for the given inner product. This
set-up nicely solves the divide between Fisher's and Rao's
approach by producing a unified theory. Moreover, the specific type of
affine manifold relevant for IG is a Hessian manifold. That
is, all its structure depends on a master convex functional.

\noindent$\bm 5.$ Both theoretical and applied research have recently
shown interest in a particular type of non-parametric statistical
models. Namely, models where the real space $\mathbb{R}^n$ is a model for
the sample space, the reference measure is either the Lebesgue measure
or the Gaussian measure and the densities are required to have some
level of smoothness.

For example, there is the statistical estimation method based on
Hyv\"arinen's divergence,
\begin{equation} \label{eq:hyvarinen}
  \KH p q = \frac12 \int \avalof{\nabla \log p(x) - \nabla \log q(x)}^2 \ p(x) \ dx \ ,
\end{equation}
see \cite[\S 13.6.2]{amari:2016},
where $p,q$ are positive probability densities of the $n$-dimensional
Lebesgue space. It is assumed that the $\log$-densities are smooth and
the integral exists. If we define the Otto's inner product by
\begin{equation} \label{eq:otto}
  \llangle f,g \rrangle _ p = \int \nabla f(x) \cdot \nabla g(x) \ p(x) \ dx \ ,
\end{equation}
where $p$ is a probability density and $f,g$ are smooth random
variables such that $\expectat p f = \expectat p g = 0$, then the
development of the square in the Hyv\"arinen divergence produces the
term
\begin{equation*}
  \llangle \log p - \expectat p {\log p} , \log q - \expectat p {\log q}
  \rrangle _ p \ .
\end{equation*}
Notice that the equation above with $p=q$ reminds of the Fisher
information of \cref{eq:fisher}, with parameter derivative replaced by
spatial derivatives. In particular, parameter and spatial derivatives
coincide in the case of translation models, that is, when
$p(x,\theta) = p(x+\theta)$.

In the Gaussian case, the celebrated log-Sobolev inequality, see
\cite[Ch. 5]{ledoux:2001-concentration}, can be written as
\begin{multline*}
  \operatorname{Ent}(p) = \int p(x) \log p(x) \ \gamma(x) dx \leq \\ 2
  \int \avalof {\nabla \sqrt {p(x)}} ^2 \ \gamma(x)dx = \frac 12
\llangle \log p - \operatorname{Ent}(p) , \log p - \operatorname{Ent}(p)
  \rrangle _ p \ , \end{multline*}
where $\gamma$ is the standard Gaussian density, $p$ is a density with
respect to $\gamma$ and $\avalof \cdot$ is the Euclidean norm.

In the following, I shall focus on the exponential representation of
positive densities $p = \mathrm{e}^{u - K(u)}$. The reference measure is
specialised to be $\mu(dx) = \gamma(x) dx$ where
$\gamma$ is the standard Gaussian density of $\mathbb{R}^n$. The
sufficient statistics $u$ is assumed to belong to an exponential
Orlicz space, to be defined next, and such that
$\gaussint {u(x)} x = 0$. The normalizing constant $u \mapsto K(u)$ is
a convex function defined on the exponential Orlicz space. The interior of
the proper domain of $K$ is the parameter set of a maximal
exponential Gaussian space model. Here, maximal means that the
model contains all possible finite dimensional exponential families. I
will aim to sketch a theory in which all concerns of items 1 to 5 meet
a solution of a sort.

\section{Orlicz spaces}
\label{sec:stat-grad-model}
First, let us review a few elements of the theory of Orlicz spaces and
fix a convenient notation. I do not not aim to full generality,
cf. \cite[Ch. 8]{adams|fournier:2003}. If $\phi \in C[0,+\infty[$ satisfies:
\begin{enumerate}
\item  $\phi(0)=0$,
\item $\phi$ is strictly increasing, and
\item $\lim_{u \to +\infty} \phi(u) = +\infty$,
\end{enumerate}
then its primitive function
\begin{equation*}
  \Phi(x) = \int_0^x \phi(u) \ du \ , \quad x \geq 0 \ ,
\end{equation*}
is strictly convex. The function $\Phi$ is extended to $\mathbb{R}$ by
symmetry, $\Phi(x) = \Phi(\avalof x)$, and is called \emph{Young
  function}.

The inverse function $\psi = \phi^{-1}$ has the same properties 1) to 3) as $\phi$, so that its primitive
\begin{equation*}
  \Psi(y) = \int_0^y \psi(v) \ dv \ , \quad y \geq 0 \ ,
\end{equation*}
is again a Young function. The couple $(\Phi, \Psi)$, is a couple of \emph{conjugate} Young functions. The relation is symmetric and we write both $\Psi=\Phi_*$ and $\Phi = \Psi_*$. The Young inequality holds true,
\begin{equation*}
  \Phi(x) + \Psi(y) \geq xy \ , \quad x,y \geq 0 \ ,
\end{equation*}
and the Legendre equality holds true ,
\begin{equation*}
  \Phi(x) + \Psi(\phi(x)) = x \phi(x) \ , \quad x \geq 0 \ .
\end{equation*}

Here are the specific cases I am going to use. The sub-2 index
denotes the 2nd Taylor remainder.
\begin{align}
&\Phi_\alpha(x) = \frac {x^\alpha} \alpha \ , \quad \Psi_\beta(y) =
                \frac {y^\beta}\beta \ , \quad \alpha, \beta > 1 \  , \quad \frac1\alpha+\frac1\beta = 1 \ ;   \label{eq:cases-1} \\
  &\exp_2(x) = \mathrm{e}^x - 1 - x \ , \quad (\exp_2)_*(y) = (1+y)\log(1+y) - y \ ; \label{eq:cases-2} \\
  &\cosh_2(x) = \cosh x - 1 \ , \quad (\cosh_2)_*(y) = \int_0^y \sinh^{-1}(v) \ dv \ ; \label{eq:cases-3} \\
&\operatorname{gauss}_2(x) = \expof{\frac12x^2}-1 \ . \label{eq:cases-4}
\end{align}

Given a Young function $\Phi$ and a probability measure $\mu$, the
\emph{Orlicz space} $\orliczof \Phi \mu$ is the Banach space whose closed
unit ball is $\setof{f \in L^0(\mu)}{\int \Phi(\avalof f) \ d\mu \leq
  1}$. This defines the \emph{Luxemburg norm}, characterized by
\begin{equation*}
  \normat {\orliczof \Phi \mu} f \leq \rho \quad \text{if, and only if,} \quad \int \Phi(\rho^{-1} \avalof f) \ d\mu   \leq 1 \ .
\end{equation*}

Because of the Young inequality, it holds
\begin{equation*}
  \int \avalof{uv} \ d\mu \leq \int \Phi(\avalof u) \ d\mu + \int \Phi_*(\avalof v) \ d\mu \ .
\end{equation*}
This provides a separating duality $\scalarat \mu u v = \int uv \ d\mu$ of $\orliczof {\Phi} \mu$ and $\orliczof {\Phi_*} \mu$ such that
\begin{equation*}
  \scalarat \mu u v \leq 2 \normat {\orliczof{\Phi} \mu} u  \normat {\orliczof {\Phi_*} \mu} v \ .
\end{equation*}
From the conjugation between $\Phi$ and $\Psi$, an equivalent norm can
be defined, namely, the \emph{Orlicz norm}
\begin{equation*}
  \normat {{\orliczof {\Phi} \mu }^*} f = \sup \setof{\scalarat \mu f
    g}{\normat {\orliczof {\Psi} \mu} f \leq 1} \ .
\end{equation*}

The domination relation between Young functions imply continuous
injection properties for the corresponding Orlicz spaces. We wll say
that $\Phi_2$ \emph{eventually dominates} $\Phi_1$, written
$\Phi_1 \prec \Phi_2$, if there is a constant $k$ such that
$\Phi_1(x) \leq \Phi_2(kx)$ for all $x$ larger than some $\bar x$. As,
in our case, $\mu$ is a probability measure, the continuous embedding
$\orliczof {\Phi_2} \mu \to \orliczof {\Phi_1} \mu$ holds if, and only
if, $\Phi_1 \prec \Phi_2$. If $\Phi_1 \prec \Phi_2$, then
$(\Phi_2)_* \prec (\Phi_1)_*$.

A special case occurs when there exists a function $C$ such that
$\Phi(ax) \leq C(a) \Phi(x)$ for all $a \geq 0$, which is the case, for
example, for a power function and in the case of the functions
$(\exp_2)_*$ and $(\cosh-1)_*$. In such a case, the dual couple is a
couple of reflexive Banach spaces, and bounded functions are a dense
set. I will return to this important topic below.

The spaces corresponding to power case
\eqref{eq:cases-1} coincides with the  ordinary
Lebesgue spaces. The norm are related by
\begin{equation*}
  \normat {\orliczof {\Phi_\alpha} \mu } f = \alpha^{1/\alpha} \normat
  {L^\alpha(\mu)} f \ .
\end{equation*}
With reference to our examples \eqref{eq:cases-2} and
\eqref{eq:cases-3}, we see that $\exp_2$ and $(\cosh-1)$ are
equivalent. They both are eventually dominated by
$\operatorname{gauss}_2$ \eqref{eq:cases-4} and eventually dominate
all powers \eqref{eq:cases-1}.  The cases \eqref{eq:cases-2} and
\eqref{eq:cases-3} provide isomorphic B-spaces
$\orliczof {(\cosh-1)} \mu \leftrightarrow \orliczof {\exp_2} \mu$
which are of special interest for us as they provide the model spaces
for our non-parametric version of IG, see
\cref{sec:exponential-bundle} below.

Clearly, a function belongs to the space $\orliczof {\cosh_2} \mu$ if,
and only if, its \emph{moment generating function
$\lambda \mapsto \int \mathrm{e}^{\lambda f}$ is finite in a neighborhood
of 0}. In turn, this implies that the moment generating function is
analytic at 0, see \cite[Ch. 2]{brown:86}.

The same property is equivalent to a large deviation inequality, see
\cite[Ch. 2]{wainwright:2019-HDS}. \emph{A function $f$ belongs to
  $\orliczof {\cosh-1} \mu$ if, and only if, it is
  \emph{sub-exponential}, that is, there exist constants $C_1,C_2 > 0$
  such that}
\begin{equation*}
  \mu (\avalof{f} \geq t) \leq C_1 \expof{-C_2 t} \ , \quad t \geq 0 \ .
\end{equation*}

Let us check the equivalence above. If $\normat {\orliczof {\cosh_2}
  \mu} f = \rho$, then $\int \mathrm{e}^{\rho^{-1} \avalof f} \ d\mu \leq 4$. It
  follows that
  \begin{equation*}
    \mu(\avalof f > t) =
    \mu\left(\mathrm{e}^{\rho^{-1}\avalof f} >
      \mathrm{e}^{\rho^{-1}t}\right) \leq
    \left(\int \mathrm{e}^{\rho^{-1}\avalof f} \ d\mu\right)
    \mathrm{e}^{- \rho^{-1}t} \leq
    4 \mathrm{e}^{- \rho^{-1}t} \ .
  \end{equation*}
The sub-exponential inequality holds with $C_1=4$ and $C_2 = \normat
{\orliczof {\cosh_2} \mu} f ^{-1}$. Conversely, for all $\lambda > 0$,
\begin{equation*}
  \int \mathrm{e}^{\lambda f} \ d\mu \leq \int_1^\infty \mu\left(\mathrm{e}^{\lambda
      f^+}> t\right) \ dt \leq C_1 \int_0^\infty \mathrm{e}^{-(C_2
    \lambda^{-1} - 1)s} \ ds \ .
\end{equation*}
The right-hand side is finite if $\lambda < C_2$ and the same
bound holds for $-f$.

A sub-exponential random variable is of particular interest in
applications because they admit an explicit exponential bound in the
Law of Large Numbers. Another class of interest consists of
the \emph{sub-Gaussian} random variables, that is, those
random variables whose square is sub-exponential.

The theory of sub-exponential random variables provides an
\emph{equivalent norm for the space} $\orliczof {\cosh_2} \mu$, namely
the norm
\begin{equation*}
  f \mapsto \sup_k \left((2k)!^{-1} \int f^{2k} \ d\mu
  \right)^{1/2k} = \left\bracevert f \right\bracevert _ {\cosh_2} \ .
\end{equation*}
See \cite{buldygin|kozachenko:2000} or
\cite{siri|trivellato:2020-SPL}. Let us prove the equivalence. If
$\normat {\orliczof {\cosh_2} \mu} f \leq 1$, then
\begin{equation*}
  1 \geq \int \cosh_2 f \ d\mu \geq \frac1{(2k)!} \int f^{2k} \ d\mu
    \quad \text{for all $k = 1,2,\dots$,}
\end{equation*}
so that $1 \geq  \left\bracevert f \right\bracevert _
{\cosh_2}$. Conversely, if the latter inequality holds, then
\begin{equation*}
  \int \cosh_2 (f/\sqrt 2) \ d\mu = \sum_{k=1}^\infty \frac1{(2k)!}
  \int f^{2k} \ d\mu \left(\frac12\right)^k \leq 1 \ ,
\end{equation*}
so that $\normat {\orliczof {\cosh_2} \mu} f \leq \sqrt 2$.

It is covenient to introduce a further notation. For each Young
function $\Phi$, the function $\overline \Phi(x) = \Phi(x^2)$ is again
a Young function such that
$\normat {\orliczof {\overline \Phi} \mu} f \leq \lambda$ if, and only
if, $\normat {\orliczof \Phi \mu} {\avalof f ^2} \leq
\lambda^2$. \emph{We will denote the resulting space by}
$\orliczpof \Phi 2 \mu$. For example, $\operatorname{gauss}_2$ and
$\overline{\cosh-1}$ are $\prec$-equivalent , hence the isomorphisn
$\orliczof {\operatorname{gauss}_2} \mu \leftrightarrow \orliczpof
{(\cosh-1)} 2 \mu$. As an application of this notation, consider that
for each increasing convex $\Phi$ it holds
$\Phi(fg) \leq \Phi((f^2+g^2)/2) \leq (\Phi(f^2) + \Phi(g^2))/2$. It
follows that when the $\orliczpof \Phi 2 \mu$-norm of $f$ and of $g$
is bounded by one, the $\orliczof \Phi \mu$-norm of $f$, $g$, and
$fg$, are all bounded by one. The space $\orliczof {\cosh_2} \mu$ has
a continuous injection in the Fr\'echet space
$L^{\infty-0}(\mu) = \cap_{\alpha>1} L^\alpha(\mu)$, which is an
algebra. When we need the product, we can either assume the factor are
both sub-Gaussian, or, move up the functional frawork to the
intersection of the Lebesgue spaces.

Let us now discuss specific issues of the Gaussian exponential Orlicz
spaces $\orliczof {\cosh_2} \gamma$, $\gamma$ the
standard $n$-variate Gaussian density. \emph{Dominated convergence
  does not hold in this space}. The squared-norm function
$f(x) = \avalof x ^ 2$ belongs to the Gaussian exponential Orlicz
space $\orliczof {\cosh_2} \gamma$ because
\begin{equation*}
  \int \cosh_2(\lambda f(x)) \ \gamma(x)dx < \infty \quad
  \text{for all} \quad \lambda < 1/2 \ .
\end{equation*}
The sequence $f_N(x) = f(x)(\avalof x \leq N)$ converges to $f$
point-wise and in all $L^\alpha(\gamma)$, $1 \leq \alpha < \infty$.
However, the convergence does not hold in the Gaussian exponential
Orlicz space. In fact, for all $\lambda \geq 1/2$,
\begin{multline*}
 \int  \cosh_2(\lambda (f(x)-f_N(x))) \ \gamma(x)dx =
 \int_{\avalof x > N} \cosh_2 (\lambda f(x)) \ \gamma(x)dx =
 \infty \ ,
\end{multline*}
but convergence would imply
\begin{equation*}
  \limsup_{N \to \infty} \int  \cosh_2(\lambda (f(x)-f_N(x))) \ \gamma(x)dx \leq 1 \quad \text{for all $\lambda >0$}\ .
\end{equation*}

The closure in $\orliczof {\cosh_2} \gamma$ of the vector space
of bounded functions is called \emph{Orlicz class} and it is denoted
by $M_{\cosh_2}(\gamma)$. One can prove that
$f \in M_{\cosh_2}(\gamma)$ if, and only if. the moment
generating function
$\lambda \mapsto \int \mathrm{e}^{\lambda f(x)} \ \gamma(x)dx$ is
finite for all $\lambda$, see \cite{pistone:2018}. An example is
$f(x) = x$. Bounded confergence holds in the Orlicz class. Assume
$f \in M_{\cosh_2}(\gamma)$ and consider the sequence
$f_N(x) = (\avalof x \leq N) f(x)$. Now,
\begin{equation*}
  \int \cosh_2(\lambda(f(x) - f_N(x)) \ \gamma(x) dx =
  \int_{\avalof x \geq N} \cosh_2(\lambda f(x)) \ \gamma(x) dz
\end{equation*}
converges to 0 as $N \to \infty$.

\section{Calculus of the Gaussian space}
\label{sec:calc-gauss-space}
I will review here a few simple facts about the analysis of the
Gaussian space, the so-called Malliavin's calculus, see
\cite[Ch.~V]{malliavin:1995}.

Let us denote by $C^k_\text{poly}(\mathbb{R}^n)$, $k = 0,1,\dots$, the vector space of functions which are differentiable up to order $k$ and which are bounded, together with all derivatives, by a polynomial. This class of functions is dense in $L^2(\gamma)$. For each couple $f,g \in \cpoly 1 n$, we have
\begin{equation*}
  \gaussint {f(x) \ \partial_i g(x)} x =   \gaussint {\delta_i f(x) \ g(x)} x \ ,
\end{equation*}
where the divergence operator $\delta_i$ is defined by $\delta_i f(x) = x_i f(x) - \partial_i f(x)$. Multidimensional notations will be used, for example,
\begin{equation*}
  \gaussint {\nabla f(x) \cdot \nabla g(x)} x = \gaussint {f(x) \ \delta \cdot \nabla g(x)} x \ , \quad f,g \in \cpoly 2 n \ ,
\end{equation*}
with $\delta \cdot \nabla g(x) = x \cdot \nabla g(x) - \Delta g(x)$.

For example, in this notation, the Hyv\"arinen divergence of \cref{eq:hyvarinen} with $P = p\cdot \gamma$, $Q = q \cdot \gamma$, and $p,q \in \cpoly 2 n$, becomes
\begin{equation*}
  \frac12 \int \avalof {\nabla \log p(x) - \nabla \log q(x)} ^ 2 p(x)
  \gamma(x)dx \ ,
\end{equation*}
while the Otto's inner product of \cref{eq:otto} becomes, with $P = p
\cdot \gamma$ and $f,g,p  \in \cpoly 2 n$, gives
\begin{equation*}
  \gaussint {\nabla f(x) \cdot \nabla g(x) \ p(x)} x  = \gaussint { f(x) \delta \cdot \nabla (g(x) p(x))} x \ .
\end{equation*}

Hermite polynomials $H_\alpha = \delta^\alpha 1$ provide an orthogonal
basis for $L^2(\gamma)$ such that
$\partial_i H_\alpha = \alpha_i H_{\alpha - e_i}$. Fourier expansion
in Hermite polynomials provides proof of the closure of both operator
$\partial_i$ and $\delta_i$ on a domain which is an Hilbert subspace
of $L^2(\gamma)$. Moreover, the closure of $\partial_i$ is the
translation operator's infinitesimal generator.

\section{Exponential statistical bundle}
\label{sec:exponential-bundle}
In this section, I will very briefly review and slightly generalize my
construction of the statistical manifold as a Banach manifold modeled
on the exponential Orlicz space $\orliczof {\cosh_2}
\gamma$. For a detailed presentation, see
\cite{pistone:2013GSI} and \cite{pistone:2018}. The general set-up is
specialized to the Gaussian space.

The support of the manifold is the \emph{maximal exponential model}
$\maxexpat \gamma$ consisting of probability densities on
$(\mathbb{R}^n,\gamma)$ of the form
\begin{equation*}
  q = \expof{u - K_1(u)} \ , \quad u \in B_1 = \orliczof {\cosh_2}
  \gamma \cap \setof{u}{\int u(x) \ \gamma(x)dx = 0} \ .
\end{equation*}
The quantity $K_1(u) = \log \int \mathrm{e}^{u(x)} \ \gamma(x) dx$
is the unique normalising constant ($\log$- partition function) of
$\mathrm{e}^u$ and is assumed to be finite. A further restriction is
needed to avoid the border of the set $\set{K_1(u)<\infty}$. We can
easily prove that the
mapping $K_1 \colon \orliczof {\cosh_2} \gamma \to \mathbb{R}$ is
convex and that the topological interior of its proper domain in
non-empty because it contains the
open unit ball of $B_1$. We restrict the model to all $u$ in such
domain. The mapping
\begin{equation*}
  s_1 \colon \maxexpat \gamma \ni q \mapsto u = \log q - \expectat \gamma {\log q}
  \in B_1
\end{equation*}
provides a global chart to the manifold. If we chose to express the
velocity of a curve $t \mapsto \maxexpat \gamma$ by the
Fishers score $\derivby t \log q(t)$, then the tangent bundle of the
manifold is expressed by the \emph{statistical bundle}
$S \maxexpat \gamma$ consisting of all the couples $(q,v)$ such
that $q$ is a density of the maximal exponential model and $v$ is a
$q$-centered random variable in the exponential Orlicz space.

The following statement is crucial to prove consistency in infinite
dimensions. It shows that the statistical bundle fibers are isomorphic
as Banach spaces.

\emph{For all $p, q \in \maxexpat \gamma$ it holds $q = \mathrm{e}^{u - K_p(u)} \cdot p$, where $u \in \orliczof {(\cosh-1)} \gamma$, $\expectat p u = 0$, and $u$ belongs to the interior of the proper domain of the convex function $K_p$. This property is equivalent to any of the following:
\begin{enumerate}\label{portmanteaux}
\item\label{portmanteaux-1} $p$ and $q$ are connected by an open exponential arc;
\item\label{portmanteaux-2} $\orliczof {\cosh_2}
  p = \orliczof {\cosh_2}
  q$ and the norms are equivalent;
\item\label{portmanteaux-3}
  $p/q \in \cup_{a > 1} L^a(q)$ and $q/p \in \cup_{a>1} L^a(p)$.
\end{enumerate}}
See
\cite{santacroce|siri|trivellato:2016,santacroce|siri|trivellato:2018}
for a detailed proof. The following argument is a generalization of the the proof of
\cref{portmanteaux-1}$\Rightarrow$\cref{portmanteaux-2}. Let $F$ be
logarithmically convex on $\mathbb{R}$ and such that $\Phi = F - 1$ is a
Young function. For example, the assumption holds for both
$F(x)=\cosh x$ and $F(x) = \mathrm{e}^{x^2/2}$. For all real $A$ and $B$,
the function
\begin{equation*}
  \mathbb{R}^2 \ni  (\lambda,t) \mapsto F(\lambda A) \mathrm{e}^{tB} = \expof{\log F(\lambda
    A) + tB}
\end{equation*}
is convex and so is
\begin{equation*}
  C(\lambda,t) = \int F(\lambda f(x)) \mathrm{e}^{tu(x)} \
  \gamma(x)dx \ ,
\end{equation*}
where $f \in \orliczof \Phi \gamma$ with and
$u \in \orliczof {\cosh_2} \gamma$ with
$\int u(x) \ \gamma(x)dx = 0$. Without restriction of
generality, assume $\normat {\orliczof \Phi \gamma} f
=1$. Let us derive two marginal inequalities. First, for $t=0$, the
definition of Luxemburg norm gives
\begin{equation*}
  C(\lambda,0) = \int F(\lambda f) \ \gamma(x) dx
  \leq 2 \ , \quad -1 \leq \lambda \leq 1 \ .
\end{equation*}
Second, for $\lambda = 0$, consider
$K_1(tu) = \log \int \mathrm{e}^{tu} \ \gamma(x)dx$,
where $t$ belongs to an an open interval $I$ containing $[0,1]$ and
such that $K_1(tu) < + \infty$. It follows that
\begin{equation*}
  C(0,t) = \int \mathrm{e}^{tu} \ \gamma(x)dx = \mathrm{e}^{K_1(tu)} <
  + \infty \ .
\end{equation*}

Choose a $t>1$ in $I$ and consider the convex combination
\begin{equation*}
  \left(\frac{t-1}t,1\right) = \frac{t-1}t (1,0)+\frac1t (0,t)
\end{equation*}
and the inequality
\begin{equation*}
  C\left(\frac{t-1}t,1\right) \leq \frac{t-1}t C(1,0)+\frac1t C(0,t)
  \leq 2 \frac{t-1}t + \frac1t \mathrm{e}^{K_1(tu)} \ .
\end{equation*}
Now,
\begin{multline*}
  \int \Phi\left(\frac{t-1}t f(x)\right) \mathrm{e}^{u(x)-K_1(u)}
  \gamma(x)dx = \\ \int F\left(\frac{t-1}t f(x)\right)
  \mathrm{e}^{u(x)-K_1(u)} \gamma(x)dx - 1 = \\
  \mathrm{e}^{-K_1(u)} C\left(\frac{t-1}t,1\right) - 1 \leq
   \mathrm{e}^{-K_1(u)} \left(2\frac{t-1}t + \frac1t
     \mathrm{e}^{K_1(tu)}\right) - 1 \end{multline*}

 As the right-hand-side is finite, we have proved that
 $f \in \orliczof \Phi p$ for $p = \mathrm{e}^{u - K_1(u)}$. Conversely, a
 similar argument shows the other implication. We have proved that all
 Orlicz spaces $\orliczof \Phi p$, $p \in \maxexpat \gamma$ are
 equal. In turn, equality of spaces implies the equivalence norms. It
 is possible to derive explicit bounds by choosing a $t$ such that the
 right-hand-side is smaller or equal to 1, see
 \cite{siri|trivellato:2020-SPL}.

\section{Gaussian Orlicz-Sobolev spaces}
In this final section, I plan to extend the exponential statistical
bundle's construction to allow for differentiable densities in the
Gaussian space. Namely, I suggest a model with sub-exponential random
variables of the Gaussian space whose weak derivatives are
sub-Gaussian random variables. The presentation focusses on functional
analytic properties, see \cite{adams|fournier:2003} and
\cite{brezis:2011FASSPDE} as a general reference.

\emph{As a model of the statistical manifold, let us define the
  function space
\begin{equation}
  \label{eq:sobolev+}
  W^1 \orliczpof {\cosh_2} {1,2} \gamma =  \setof{f\in \orliczpof
    {\cosh_2} 1 \gamma}{\partial_i f
    \in \orliczpof {\cosh_2} 2 \gamma, i=1,\dots,n} \ ,
\end{equation}
where the weak partial derivative $\partial_j f$ exists if
\begin{equation*}
  \scalarat \gamma {\partial_j f}{\phi} = \scalarat \gamma
  {f}{\delta_j \phi} \quad \text{for all} \quad \phi \in C^1_0(\mathbb{R}^n) \ .
\end{equation*}}

The above definition of weak derivative coincides this the usual
definition of derivative in the sense of Schwartz distributions
because $\phi \leftrightarrow \phi \cdot \gamma$ is a bijection of
$C_0^\infty(\mathbb{R}^n)$ and
\begin{equation*}
  \scalarat \gamma f {\delta_i \phi} = - \int f(x) \pderivby
  {x_i}
  \left(\phi(x)\gamma(x)\right) \ dx \ .
\end{equation*}

Weak derivatives do not provide tools for non-linear
computations. Hence we want to recall the relation between translation
and weak derivative, the weak version of Calculus's fundamental
theorem. Let be given a locally integrable real mapping
$G \in L^1_\text{loc}(\mathbb{R})$, and assume there exists a locally
integrable function $G'$ which is the weak derivative of $G$, that is,
\begin{equation}\label{eq:dG-weak}
  \int G(x) \phi'(x) \ dx = - \int G'(x) \phi(x) \ dx \ , \quad \phi
  \in C_0^\infty(\mathbb{R}) \ .
\end{equation}
Define the translation $\tau_h G$, $h\in\mathbb{R}$, by $t_h G(x) =
G(x-h)$, $h\in\mathbb{R}$. It follows immediatly $\tau_h G \in
L^1_\text{loc}(\mathbb{R})$ and
\begin{align*}
  \int (\tau_{-h} G(x) - G(x)) \phi(x) \ dx
  &= \int G(x+h) \phi(x) \
    dx - \int G(x) \phi(x) \
    dx \\
  &= \int G(x) ( \phi(x-h) - \phi(x)) \ dx \\
  &= \int G(x) \left. \phi(x-sh) \right|_{s=0}^{s=1} \ dx \\
  &= - h \int G(x) \int_0^1 \phi'(x-sh)\ ds \ dx \\
  &= - h \int_0^1 \int G(x) \phi'(x-sh) \ dx \ ds\\
  &= h \int_0^1  \int G'(x) \phi'(x-sh) \ dx \ ds = \\
  &= h \int_0^1  \int G'(x+sh) \phi'(x) \ dx \ ds = \\
  &= \int \left(h \int_0^1 G'(x+sh) \ ds\right) \phi(x) \ dx
\end{align*}
As $\phi$ is any function in $C^\infty_0(\mathbb{R})$, we have proved that
\begin{equation}\label{eq:delta-G}
  \tau_{-h} G - G = h \int_0^1 \tau_{-sh} G' \ ds = h G' + h \int_0^1 (\tau_{-sh} G' - G') \ ds
\end{equation}
in $L^1_\text{loc}(\mathbb{R})$. In particular, if $G'$ is bounded by a
constant $K$, then $G$ is almost surely $K$-Lipschitz, $\avalof{G(x-h)
  - G(x)} \leq K \avalof h$.

Conversely, if \cref{eq:delta-G} holds, then \cref{eq:dG-weak} holds,
see, for example, \cite[Lemma~8.1-2]{brezis:2011FASSPDE}.

The argument above extends to $n$-variate functions $f \in
W^{1,1}_\text{loc}(\mathbb{R}^n)$, that is, $f, \partial_i f
\in L^1_\text{loc}(\mathbb{R}^n)$, $i = 1,\dots,n$, by considering, for
each $h \in \mathbb{R}^n$, the univariate function $t \mapsto \tau_{th} f$
defined by $\tau_{th} f(x) = f(x-th)$. We obtain
\begin{equation}\label{eq:delta-f}
  \tau_{-th} f - f = t \nabla f \cdot h + t \int_0^1 (\tau_{-sth}
  \nabla f - \nabla f) \cdot h \ ds \ ,
\end{equation}
where the equality holds in $L^1_\text{loc}(\mathbb{R}^n)$. The same
equality holds in all function space whose elements are locally
integrable.

The result about the functional differentiability of translations is
the following. \emph{For all
  $f \in W^1 \orliczpof {\cosh_2} {1,2} \gamma$ the following
  first increment equation holds,
\begin{equation*
  (\tau_{-th} f - f) - t \nabla f \cdot h = t \int_0^1 (\tau_{-sth}
  \nabla f - \nabla f) \cdot h \ ds \ ,
\end{equation*}
  and differentiability holds in $L^{\infty-0}(\gamma) = \cap_{\alpha>1} L^\alpha(\gamma)$.
} In fact, the translations are continous in all $L^\alpha(\gamma)$.

It is possible to exted the previous result to a property of the
derivative of the composite function $G\circ f$.  The increment of the
composition expands as
\begin{multline}\label{eq:delta-Gf}
  G(f(x + th)) - G(f(x)) = G(f(x) + (f(x+th) - f(x)) = \\
  (f(x+th) - f(x)) G'(f(x)) + \\ (f(x+th) - f(x)) \int_0^1 (G'(f(x) + s
  (f(x+th)-f(x))) - G'(f(x))) \ ds \ ,
\end{multline}
and the weak derivative of the composite function exists if $G'$ is
bounded and $f \in W^1 \orliczpof {\cosh_2} {1,2}
\gamma$. However, differentiability holds in
$L^{\infty-0}(\gamma)$.

An interesting example of application is the neuron of a neural
network \cite[Ch. 18]{efron|hastie:2016}. If
$f_1, \dots, f_k \in W^1 \orliczpof {\cosh_2} {1,2} \gamma$,
$G$ is an activation function with linear growth, for example, the
linear rectifier
$G(x) = x^+$, and $a_i, w_{ij}, b_i$ are given constants, then
\begin{equation*}
  \sum_{i=1}^h a_i G\left(\sum_{j=1}^k w_{ij} f_j - b_i\right) \in W^1 \orliczpof {\cosh_2} {1,2}
\gamma \ .
\end{equation*}

A crucial feature of the space defined in \cref{eq:sobolev+} is
the fact that \emph{each element of $W^1 \orliczpof {\cosh_2} {1,2}
  \gamma$ has a continuous version}. The embedding
\begin{equation*}
  \orliczpof {\cosh_2} 1 \gamma, \orliczpof {\cosh_2} 2
\gamma \subset \cap_{\alpha \geq 1} L^\alpha(\gamma)
\end{equation*}
allows to use the standard Sobolev inequalities to our
case. $W_{\text{loc}}^{1,\alpha}(\mathbb{R}^n)$ denotes the space of
funcions whose restriction to each open ball
$B_\rho = \setof{x \in \mathbb{R}^n}{\avalof x < \rho}$ is
$\alpha$-integrable, together with all weak partial
derivatives. $C^\lambda(\overline B_\rho)$ denotes $\lambda$-H\" older
functions on the closed ball.

\emph{\begin{enumerate}
    \item
  The following restriction and embedding hold true and are continuous,
  \begin{equation*}
    W^1 \orliczpof {\cosh_2} {1,2} \gamma \to
    W^{1,\alpha}(\overline B_\rho) \subset C^\lambda(\overline
    B_\rho)\ , \quad \rho > 0 \ , \quad 0 < \lambda < 1 \ .
  \end{equation*}
\item
  The following inclusions hold true and are continuous:
  \begin{equation*}W^1 \orliczpof {\cosh_2} {1,2} \gamma \subset \cap_{\alpha
    \geq 1} W^{1,\alpha}_{\textnormal{loc}} \cap \orliczpof
  {\cosh_2} 1 \gamma \subset
  C(\mathbb{R}^n) \cap \orliczpof {\cosh_2} 1 \gamma \ ,
\end{equation*}
where the space of continuous functions $C(\mathbb{R})$ is endowed with
the uniform convergence on compact sets.
\end{enumerate}
}

The embedding are easily verified.  If
$f \in \orliczpof {\cosh-1} 1 \gamma$, then, for all
$k \in \mathbb{N}$, the inequalities $x^{2k}/(2k)! \leq \cosh_2(x)$ and
$(2\pi)^{-n/2} \mathrm{e}^{-\rho^2/2} \leq \gamma(x)$ for
$x \in B_\rho$ imply the inequality
 \begin{multline*}
   \frac{(2\pi)^{-n/2}\mathrm{e}^{-\rho^2/2}}{(2k)!} \int_{B_\rho}
   \left(\frac{f(x)}{\normat {\orliczpof {\cosh_2} 1 \gamma} f}\right)^{2k} \
   dx \leq \\ \int \cosh_2\left(\frac{f(x)}{\normat
       {\orliczpof {\cosh_2} 1 \gamma} f}\right) \gamma(x) \ dx \leq 1 \ ,
 \end{multline*}
so that
\begin{equation*}
  \normat {L^{2k}(B_\rho)} f ^{2k} \leq (2\pi)^{n/2} (2k)! \mathrm{e}^{\rho^2/2}
  \normat {\orliczpof {\cosh_2} 1 \gamma} f \ .
\end{equation*}
A similar argument applies to the weak partial derivatives. Now we can
use the Sobolev embedding theorem, see
\cite[Th. 4.12]{adams|fournier:2003}.

Let us conclude by explicitly reviewing the main properties of our
space.

\emph{\begin{enumerate}
\item The space $W^1 \orliczpof {\cosh_2} {1,2}
\gamma$ contains the constants and all polinomial up to order
2.
\item Each element has a continuous version.
\item If $G : \mathbb{R} \to \mathbb{R}$ is the primitive of a bounded
  function, then $G \circ f \in W^1 \orliczpof {\cosh_2} {1,2} \gamma$.
\item $\min(f,g), \max(f,g) \in W^1 \orliczpof {\cosh_2} {1,2} \gamma$.
\end{enumerate}
}

Moreover, our space is a Banach space: \emph{The mapping
  \begin{equation*}
  W^1 \orliczpof \Phi {1,2} \gamma \ni f \mapsto \normat
  {\orliczpof \Phi 1 \gamma} f + \sum_{i=1}^n \normat {\orliczpof \Phi 2
    \gamma} {\partial_i f}
  \end{equation*}
is a complete norm and thus defines a Banach space.}
The argument is a standard one in Functional Analysis. The weak gradient $\nabla$ is a closed operator from $\orliczpof \Phi
1 \gamma \to \left(\orliczpof \Phi 2 \gamma\right)^n$,
that is, the graph of $\nabla$ is closed in $\orliczpof \Phi
1 \gamma \times \left(\orliczpof \Phi 2
  \gamma\right)^n$. In fact, given a converging sequence in the
graph, say $f_n \to f$ and $\partial_i f
\to f^i$, it holds
\begin{equation*}
  \scalarat \gamma {\partial_i f} {\phi} = \scalarat \gamma {f}
  {\delta_i \phi} = \lim_n \scalarat \gamma {f_n} {\delta_i \phi}
  = \lim_n \scalarat \gamma
  {\partial_i f} {\phi} = \scalarat \gamma {f^i} {\phi} \ .
\end{equation*}
The identification of the space $W^1 \orliczpof \Phi {1,2}
\gamma$ with the graph of $\nabla$ provides a complete
norm.

Finally, all exponential Orlicz spaces are isomorphic in a
maximal exponential model. This, in turn implies the isomorphism
\begin{equation*}
  W^1 \orliczpof {\cosh_2} {1,2}
\gamma \leftrightarrow   W^1 \orliczpof {\cosh_2} {1,2}
{p \cdot \gamma} \quad p \in \maxexpat \gamma \ .
\end{equation*}
It follows that this space is suitable as a model of the fibers of a
statistical bundle.

\section{Selected bibliography}
\label{sec:select-bibl}
General monograph on IG are \cite{kass|vos:1997},
\cite{amari|nagaoka:2000}, and
\cite{Ay|Jost|Le|Schwachhofer:2017IGbook}.  The type of analytic
framework I suggest to use in IG is that of Banach manifolds as in
\cite{lang:1995}. A first version of this project has been developed
in
\cite{pistone|sempi:95,gibilisco|pistone:98,pistone|rogantin:99,cena:2002}. In
a series of papers \cite{pistone:2013Entropy},
\cite{lods|pistone:2015}, \cite{pistone:2017-GSI2017}, and
\cite{pistone:2018-IGAIA-IV}, we have explored a version of the
non-parametric Information Geometry (IG) for smooth densities on
$\mathbb{R}^n$. Especially, we have considered the IG associated to Orlicz
spaces on the Gaussian space. The analysis of the Gaussian space is
discussed, for example, in \cite{malliavin:1997}, and
\cite{nourdin|peccati:2012}.

This set-up provides a way to construct a statistical manifold
modelled on functional spaces of smooth densities, but other modelling
options are in fact available, for example the global analysis methods
of \cite{kriegl|michor:1997} and the approach based on deformed
exponentials \cite{newton:2019LNCS}. For the two examples involving
derivatives, see \cite{hyvarinen:2005}, \cite{lods|pistone:2015},
\cite{otto:2001}, and \cite{lott:2008calculations}.

 We have worked with a restricted type of Young functions. See the
 more general cases in \cite[Ch.~II]{musielak:1983} and
 \cite[Ch.~VII]{adams|fournier:2003}. There is a large literature
 about sub-exponential random variables, for example,
 \cite{buldygin|kozachenko:2000}, \cite{vershynin:2018-HDP},
 \cite{wainwright:2019-HDS}. Application to IG is discussed in
 \cite{siri|trivellato:2020-SPL}. The need to control the product of
 two random variables in $\orliczof {(\cosh-1)} \mu$ appears, notably,
 in the study of the covariant derivatives of the statistical bundle,
 see \cite{gibilisco|pistone:98}, \cite{lott:2008calculations},
 \cite{pistone:2018Lagrange}, and
 \cite{chirco|malago|pistone:2020-2009.09431}.

See \cite{pistone:2013GSI}, and \cite{pistone:2018-IGAIA-IV} for all
details about the contruction of the statistical bundle that are
missing here. In particular the isomorphism of fibers is discussed in
detail in \cite{santacroce|siri|trivellato:2016} and
\cite{santacroce|siri|trivellato:2018}. For Malliavin's calculus, see
\cite{malliavin:1997} and \cite{malliavin:1995}. General references
about the functional background used in the construction of the
Gaussian Orlicz-Sobolev space are \cite{adams|fournier:2003}, and
\cite{brezis:2011FASSPDE}.

\subsection*{acknowledgements}The author acknowledges the support by de Castro
  Statistics, Collegio Carlo Alberto, Turin, Italy. He is a member of
  GNAMPA-INDAM. The author thanks an anonymous reviewer and L. Malag\`o for their helpful comments.

\bibliographystyle{splncs04}


\section{Finite-State Machine}

A Finite-State Machine (FSM) represents a computation as a finite number of states
labeled $[0, N)$ and transitions between those states.  A FSM starts at an
initial state---represented as {\sethlcolor{green}\hl{green}} in
this paper---and performs a series of transitions following the directed edges
until it reaches a final state---represented as {\sethlcolor{red}\hl{red}} in
this paper.  Transitions in a FSM can be labeled with symbols (e.g. a character
in a string).  In general, a Finite-State Transducer (FST) contains two symbols
on each edge.  One symbol represents what is \emph{read} from the FST's input,
and the other represents what is \emph{written} to the output
(\cref{fig:fst-hello-world}).  In the special case that both symbols are the
same along all edges, a FST can also be called a Finite-State Acceptor (FSA)
(\cref{fig:fsa-hello}).  Furthermore, edges in a FST can also be weighted.
This allows for a FST to score different paths that accept a given string.
If an input can be mapped to multiple outputs, then each output will have a
different weight assigned (\cref{fig:wfst-hello-world-troll}).

In the remainder of this paper, we will assume the reader is familiar with the theory and
algorithms behind FSTs.  Instead, we will focus on {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace and how it makes working with
FSTs from Python easy.  The examples in this paper are intended to demonstrate
how easy it is to get started constructing simple models rather than
demonstrating ``real'' models that work with real scenarios.

\begin{figure}[h]
\centering
\includegraphics[width=.8\textwidth]{notebook-demo/hello.png}
\caption{Finite-State Acceptor (FSA) for the string \emph{``hello.''}  State \#0, is
  colored green as it is the initial state of the FSA.  State \#5 is
  marked red as it is a terminal state.  Strings that do not exactly match
  \emph{``hello''} will not reach the final state (red), thus are not accepted.}
\label{fig:fsa-hello}
\end{figure}

\begin{figure}[h]
  \centering
  \hspace{-10mm}
  \includegraphics[width=.8\textwidth,trim=120 0 100 450,clip]{notebook-demo/hello-world-transducer.png}
  \caption{Finite-State Transducer (FST) that transforms the string \emph{``hello''}
    $\to$ \emph{``world.''}  During each transition, a character of the string
    \emph{``hello''} is read from the input string---represented on the left side of the
    colon.  On the right hand side of the colon, the output ``world'' is generated one character at a time as each arc is transversed.  When the
    characters are the same on both the input and output side (e.g. between
    states $3 \to 4$), the colon is omitted for readability.}
  \label{fig:fst-hello-world}
\end{figure}

\begin{figure}[h]
  \centering
  \includegraphics[width=.8\textwidth,trim=120 0 80 555,clip]{notebook-demo/weighted-hello-world-two-paths.png}

  \caption{Weighted Finite-State Transducer where \emph{``hello''} is transduced
    non-deterministically into both the words \emph{``world''} and \emph{``troll.''}  Each
    output word is weighted with the product of the weights along a given path.
    Here, \emph{``hello''} is given the weight $1*1*1*2*1*3 = 6$, and \emph{``troll''} is
    given $1*1*1*2*2*3 = 12$.  The weight $1$ is omitted from the figure for
    readability.  If there exist multiple paths that generate the same output,
    the weight would be summed across the different paths.}
  \label{fig:wfst-hello-world-troll}
\end{figure}

\section{{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace: Getting Started}
{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace was initially developed for teaching advanced FST techniques in a classroom
setting~\cite{seq2class_assignment2}.  This has motivated {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace's simple interface,
which makes it quick to get started.  {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace includes sensible defaults all
accessible from Python, while not eliminating power user features of OpenFST.

\subsection{One Command Install}
{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace is installable with a single command.  {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace includes and automatically
compiles OpenFST, requiring no additional steps from a user\footnote{Note: The
  install command takes around 10 minutes to run without printing any indication
  it is running to the terminal.}.
\begin{lstlisting}
pip install 'git+https://github.com/matthewfl/openfst-wrapper.git'
\end{lstlisting}

\subsection{First Simple FSAs}
Once {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace is installed, it makes the first steps of interacting with a FSTs from
Python easy.  Here, a FSA is constructed for the string \emph{``hello.''}  The FSA is
constructed using {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace's default semiring---more on that later in \cref{sec:weighted-fsts}.

\begin{figure}[h]
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/hello2.png}
  \caption{The simplest way to construct a FSA for the string \emph{``hello.''}  When
    a FST is returned from a Jupyter notebook cell, it is drawn into the
    notebook.  All figures in this paper are from {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace's Jupyter drawings.
    Here we are exemplifying how any Python iterable (such as a string) is
    converted into the labels on a FSA.}
\label{fig:fsa-hello-code}
\end{figure}

{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace exposes all of OpenFST's
operations\footnote{OpenFST included operations:
  \url{http://www.openfst.org/twiki/bin/view/FST/FstQuickTour\#Available\%20FST\%20Operations}.\label{footnote:openfst_docs}}
for interacting and manipulating a FST.  Whenever an instance of the
\texttt{FST} class is returned from a cell in a Jyputer notebook, the FST is
automatically drawn into the notebook.  This makes {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace ideal for learning how
FSTs work, visualizing how the different operations transform a FST without
requiring any additional steps to draw the FST (shown in
\cref{fig:union-determinize}).  {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace, being a Python package, includes
documentation on all of its methods (adapted from OpenFST's documentation).
This documentation can be easily accessed in an interactive Python environment
using Python's \texttt{help()} method (\cref{fig:docs-included}).

\begin{figure}[h]
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/union-determinize.png}
  \caption{Simple FSA constructed by taking the union of the FSA for the string
    \emph{``hello''} and \emph{``help.''}  In cell \#3, the FSA is referenced using Pythons
    underscore variable (\texttt{\_}) to reference the output of the previous
    computation.  The operations \texttt{remove\_epsilon} and
    \texttt{determinize} are chained together here to create a FSA where the
    prefixes \emph{``hel''} is merged.}
  \label{fig:union-determinize}
\end{figure}

\begin{figure}[h]
  \centering
  \begin{subfigure}[b]{.62\textwidth}
    \includegraphics[width=\textwidth,trim=0 200 400 0,clip]{notebook-demo/help-docs-included.png}
    \caption{Using Python's \texttt{help()} function to access documentation for a
      method.  The documentation is largely adapted from OpenFST online
      documentation$^{\text{\ref{footnote:openfst_docs}}}$.}
  \end{subfigure}%
  \hfill
  \begin{subfigure}[b]{.35\textwidth}
    \includegraphics[width=\textwidth,trim=0 200 150 0,clip]{notebook-demo/tab-complete.png}
    \caption{Tab competition in a Jupyter notebook to list the available methods
      on the FST class.}
  \end{subfigure}
  \caption{The documentation on {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace methods is included as Python docstrings
    and can also be easily auto-completed using any
    Python IDE.}
  \label{fig:docs-included}
\end{figure}

\section{Building Transducers}
FST can also be built using the \texttt{add\_state} and \texttt{add\_arc}
methods in addition to the simplified constructor we have demonstated so far.  These methods
makes it possible to build more complicated FSTs which may have loops or may define
input and output labels and weights on all of the edges.  The main methods for
constructing a FST are as follows, and are shown in \cref{fig:build-fst}.

\begin{itemize}
\item \texttt{add\_state()$\to$state\_id}: Modifies the FST in place to add a new state.

\item \texttt{add\_arc(from\_state, to\_state, weight, input\_label, output\_label)$\to$unit}:\linebreak Modifies the FST in place to add a new arc between two states.  If \texttt{from\_state} and \texttt{to\_state} are the same, then this will create a self-loop.  The \texttt{weight} is cast to the FST's semiring automatically.  By default \texttt{weight} is set to the semiring's identity element of $\bar{1}$.  \texttt{input\_label} and \texttt{output\_label} are 64 bit integers.  The value $0$ is used to represent the special $\epsilon$ character in the FST, which indicates that no symbol is read/written when traversing an arc.  If a single character of a string is passed, then it will be automatically converted into an integer using it character code (the method \texttt{chr} in Python).

\item \texttt{set\_initial\_state(state\_id)$\to$unit}: Modifies the FST to set its initial state.  A FST can only have a single initial state, and a FST without an initial state will not be usable.

\item \texttt{set\_final\_weight(state\_id, weight)$\to$unit}: Modifies the FST to set the final state weight.  There can be multiple final states, and a well defined FST should have at least one final state.
\end{itemize}

\begin{figure}[h]
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/hello-world-transducer.png}

  \caption{In Jupyter cell \#4, an empty FST is first constructed using
    \texttt{FST()} with the default semiring.  The edges are added one at a time
    using \texttt{add\_arc}, where most edges have different labels on the input
    and output side.  }
  \label{fig:build-fst}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=.8\textwidth]{notebook-demo/hello-world-compose.png}
\caption{In cell \#5, the FST is composed with the FSA from
    \cref{fig:union-determinize} that was the union of the word \emph{``hello''} and
    \emph{``help.''}  Only the word \emph{``hello''} is transduced as there is no path that
    accepts the word \emph{``help.''}  The \texttt{project('output')} method turns the
    FST into a FSA using the labels on the output side.}
\end{figure}

\section{Weighted FSTs} \label{sec:weighted-fsts}

Weighted FSTs allow for scoring a given sequence rather than just accepting or
transducing an input.  This is done by assigning a weight to every edge in the
FST.  The weights on a FST must all be instances of the same semiring, meaning
that they share a multiplication and addition operator as well as a multiplicative and
additive identity $\langle +, *, \bar{0}, \bar{1} \rangle$.
To determine the weight for a given input/output sequence, the weights along a
given path are multiplied together.  When there are multiple paths that accept
an input or transduce a given input/output pair, the weight will be summed
across all paths.  Changing the semiring and defining different actions for
multiplication and addition can allow the same algorithm operating on the FST to
compute different results---such as in \cref{fig:cast-fst}.

{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace includes the standard semirings.  This includes:
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L15}{Python values},
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L85}{real numbers},
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L112}{min},
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L140}{max},
 \href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L168}{tropical}, and
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/semirings.py#L178}{boolean}.
Each of these semirings is implemented in approximately 20 lines of Python each.
Custom semirings can be defined by extending the
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/__init__.py#L16}{\texttt{AbstractSemiringWeight}}
class.

{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace's default semiring---that we have been using so far---is the boolean
semiring.  In the boolean semiring, all weights along the edges are $\bar{1}$
(as the $\bar{0}$ weights are omitted).  The boolean semiring is given
special precedence in {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace as the boolean semiring will be automatically cast
to another semiring as necessary.  All operations in {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace which involve two or
more FSTs---such as compose or union---require that all of its arguments
implement the same semiring.  The boolean semiring's automatic casting is useful
as it allows us to construct FSTs and compose them with weighted FSTs without
too much concern for which semiring is in use, such as in
\cref{fig:compose-aaa}.  {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace also provides the method
\texttt{lift(semiring\_class, cast\_method)$\to$FST} to convert a FST from one
semiring to another, such as in \cref{fig:cast-fst}.

\begin{figure}[H]
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/score-double-a.png}
  \caption{Here we define a weighted FST using the \texttt{RealSemiringWeight}
    $\langle +, * \rangle$.  This FST transduces between the alphabet
    $\{ \text{a}, \text{b} \}$ and itself.  Both symbols can be read unchanged
    (as they transverse the edge from state $0 \to 0$, which has weight 1).  If
    the FST encounters two \emph{``a''}s in a row, then it is transduced
    non-deterministically into the letter \emph{``b''} with a weight of $0.5$.  The
    weight $0.5$ comes from the fact that it multiplies the weights along the
    edges from state $0\to1$ and $1\to0$ when reading the sequence \emph{``aa.''}}
  \label{fig:double-a}
\end{figure}

\begin{figure}[H]
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/compose-double-a.png}
  \caption{Here the FST defined in \cref{fig:double-a} is composed on its input
    side with the FSA which represents \emph{``aaa''} (constructed in the same way as
    \cref{fig:fsa-hello-code}).  The \texttt{FST('aaa')} is using the default
    semiring, so it is automatically cast to the RealSemiringWeight which was
    used to define \texttt{fst}.  This FST has three different paths which can
    generate outputs \emph{``aaa,''} \emph{``ba,''} or \emph{``ab.''}  The weight associated with
    each path is the result of multiplying the weights from both FSTs which are
    composed together.  The \texttt{FST('aaa')} only has a single path.
    The FST from \cref{fig:double-a} scores the \emph{``aa''}$\to$\emph{``b''}
    transduction with weight $0.5$.}
  \label{fig:compose-aaa}
\end{figure}

\begin{figure}[h]
  \centering
  \begin{subfigure}[t]{.48\textwidth}
    \includegraphics[width=\textwidth]{notebook-demo/minplus-lift.png}
  \end{subfigure}
  \begin{subfigure}[t]{.48\textwidth}
    \includegraphics[width=\textwidth]{notebook-demo/maxplus-lift.png}
  \end{subfigure}
  \caption{The semiring is cast to the min-plus and max-plus semirings and then
    we use the shortest path algorithm to identify the shortest path.  In these
    semirings, the multiplication ($\odot$) operator is summing the length of the path
    and the addition ($\oplus$) operator is defined as min/max.  This allows us to find the shortest
    path using the ``$\le$'' operator where it is defined as
    $a \le b \iff a \oplus b = a$ for idempotent semirings.}
  \label{fig:cast-fst}
\end{figure}

\subsection{Custom Weighted Semirings}

All semirings in {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace are defined in pure Python by extending the
\href{https://github.com/matthewfl/openfst-wrapper/blob/dea1c93741b0c3be3dcedaa5d6f9158ce4126453/mfst/__init__.py#L16}{\texttt{AbstractSemiringWeight}}
class.  This means that we can easily define a custom semiring that makes use
of other libraries in Python.  For example, we may be interested in learning the
weights on a FST.  One way in which this can be done is to create a featurized
semiring, where we are tracking the weights associated with every feature
rather than the weight of an edge itself (\cref{sec:basic-featurized}).

Another approach that we can use for learning the weights on a FST is to
leverage a dynamic graph neural network frameworks, such as PyTorch.
To construct a graph using standard neural network forward propagation, we simply use PyTorch's operations whenever we perform a semiring operation.
In \cref{fig:pytorch-semiring} we construct a custom
semiring that wraps a PyTorch tensors as a weight.  We can implement the
\texttt{\_\_add\_\_} and \texttt{\_\_mul\_\_} methods with anything which
results in a semiring.  For example, if we were to implement the min
semiring, then we would use \texttt{torch.min} for the \texttt{\_\_add\_\_}
method and use \texttt{torch.add} for the \texttt{\_\_mul\_\_} method.  For this example, however, we have instead chosen to exemplify how we could
perform any neural operation.  Here we are using a forward pass through
a linear layer followed by a sigmoid operation.  Note, the
``semiring'' which we define in \cref{fig:pytorch-semiring} is not a \emph{real}
semiring as it does not follow the required distributivity and identity elements
of a semiring\footnote{\url{https://en.wikipedia.org/wiki/Semiring\#Definition}}.
While {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace will still run, the outputs of an invalid ``semiring'' is unlikely
to be stable or work in practice.

\begin{figure}[h]
  %
  \centering
  \includegraphics[width=\textwidth]{notebook-demo/pytorch-semiring-define.png}
  \caption{A custom semiring where weights are 10-element PyTorch tensors.  The
    multiplication and addition operations of the semiring are calling PyTorch Neural
    Network operators.  Note: This example is intended only to illustate
    that PyTorch tensors can be passed through OpenFST as weights.  The
    ``semiring'' defined here is not, in fact, a ``mathematically correct''
    semiring, and so actually using this as written here is ill-advised.}
  \label{fig:pytorch-semiring}
\end{figure}

\begin{figure}
  \centering
  \begin{subfigure}[t]{.45\textwidth}
    \includegraphics[width=\textwidth,clip,trim=0 0 2 0]{notebook-demo/pytorch-compose-e.png}
    \caption{Composition with the default BooleanSemiring automatically casts to
      the identity elements of the TorchSemiring and then composes the two FSTs.}
  \end{subfigure}%
  \hfill
  \begin{subfigure}[t]{.45\textwidth}
    \vspace{-22mm}
    \includegraphics[width=\textwidth]{notebook-demo/pytorch-compose-sum_paths.png}
    \vspace{7mm}
    \caption{\texttt{sum\_paths} sums all of the paths in the FST and returns
      the resulting weight element in the semiring.}
  \end{subfigure}
  \caption{\texttt{torch\_fst} manipulated through {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace's built-in operations.}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=.8\textwidth]{notebook-demo/pytorch-prob.png}

  \caption{Here we exemplify how the \texttt{torch\_fst} could be used with
    PyTorch's backpropagation.  First, the numerator for the path that we are
    interested in is selected using \texttt{compose} and then
    \texttt{sum\_paths} to get the value for those specific paths.  We then get
    the normalizing constant (denominator) for all of the paths in the FST.  The
    variable \texttt{numerator} and \texttt{denominator} are both PyTorch
    tensors.  As a result, we can use PyTorch methods such as \texttt{torch.sum} and
    \texttt{torch.exp} on these variables.  As with a standard PyTorch tensor,
    we can use \texttt{tensor.backwards()} to compute the gradient.  In this
    example, we have included an optimizer, though that would be included in
    general.  Again, as mentioned above, this example is intended to illustrate how
    PyTorch tensors can be passed around through a FST as a semiring weight,
    not to illustrate something that is expected to work well.}
\end{figure}

\section{Implementation}

{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace is a thin wrapper around OpenFST's C++ interface.  OpenFST~\cite{openfst}
has the ability to define custom semiring by defining custom C++
semiring
class\footnote{\url{http://www.openfst.org/twiki/bin/view/FST/FstAdvancedUsage\#Weights}}.
{\usefont{U}{Sanremo}{xl}{n} M}FST\xspace defines a custom semiring which wraps an opaque Python object.  The Python
C++ bridge makes use of PyBind11~\cite{pybind11}, which makes passing around
Python objects in C++ and exposing C++ methods to Python largely transparent.
Anytime an operation is performed by OpenFST on the semiring (such as addition
or multiplication), the call is redirected to the relevant method defined in Python.
This also means that we can easily define our own semiring from Python by
implementing the relevant methods.  This means that all FST operations
are handled by OpenFST, hence matured and well-tuned over the years.

The images that are drawn by {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace make use of Jyputer's
\texttt{\_repr\_html\_}
method\footnote{\url{https://iPython.readthedocs.io/en/stable/config/integrating.html\#rich-display}}.
By defining a\linebreak
\texttt{\_repr\_html\_} method, this allows a class to define how it is rendered in a Jupyter
notebook by returning a string of HTML code that is used instead of the standard
textual output.  {\usefont{U}{Sanremo}{xl}{n} M}FST\xspace makes use of this method to generate HTML code that
utilizes D3.js~\cite{d3js} and Dagre-d3~\cite{dagre} to render the diagrams in
the browser.