text:
\begin{document}
\date{\bf Version of \today}
\maketitle
\begin{abstract} In this paper we make a clear relationship between the automorphic representations and the quantization through the Geometric Langlands Correspondence. We observe that the discrete series representation are realized in the sum of eigenspaces of Cartan generator, and then  present the automorphic representations in form of induced representations with inducing quantum bundle over a Riemann surface and then use the loop group representation construction to realize the automorphic representations. The Lanlands picture of automorphic representations is precised by using the Poisson summation formula.
\end{abstract}
\tableofcontents
\section{Introduction}
The representation theory of $\SL(2,\mathbb R)$ is well-known in the Bargman classification: every irreducible unitary/nonunitary representation is  unitarily/nonunitarily equivalent to one of representions in the five series:
\begin{enumerate}
\item the principal continuous series representations $(\pi_s,\mathcal P_s)$,
\item the discrete series representations $(\pi^\pm_k,\mathcal D_n), n\in \mathbb N, n \ne 0$,
\item the limits of discrete series representations $\mathcal (\pi^\pm_0,D_{\pm})$,
\item the complementary series representations $(\pi_s,\mathcal C_s), 0<s<1$,
\item the trivial one dimensional representation $1$,
\item the nonunitary finite dimensional representations $V_k$.
\end{enumerate}

Looking at a Fuchsian discrete subgroup $\Gamma$ of type I, i.e.
$$ \Gamma \subseteq \SL(2,\mathbb Z), \vol(\Gamma\backslash \SL(2,\mathbb R)) < +\infty$$
one decomposes
the cuspidal parabolic part $ L^2_{cusp}(\Gamma \backslash \SL(2,\mathbb R))$ of  $L^2(\Gamma \backslash \SL(2,\mathbb R))$, that is consisting of the so called {\it automorphic representations} in subspaces of automorphic forms, of which each irreducible component with {\it multiplicity} equal to the dimension of the space of modular forms on the upper Poincar\'e half-plane
$$\mathbb H = \{z\in \mathbb C | \Im(z) > 0\}.$$
There are a lot of studies concerning the automorphic forms and automorphic representations. Most of them realize the representations as some induced ones. Therefore some {\it clear geometric realization} of these representations should present some interests.

In this paper, we use the {\it ideas of geometric quantization} to realize the automorphic representations in form of some Fock representations of loop algebras, see Theorem \ref{mthm}, below.
In order to do this, we first present the action of the group $\SL(2,\mathbb R)$ in the induced representation of discrete series as the action of some loop algebra/group  in the heighest/lowest weight representations.

The representations are uniquely up to equivalence defined by the character, which are  defined as some distribution function.
Beside of the main goal we make also some clear presentation of automorphic representations in this context, see Theorems \ref{thm1}, \ref{thm2}.
In the theorem \ref{thm3} we expose the trace formula using both the spectral side and geometric side.

\section{Endoscopy groups for $\SL(2,\mathbb R)$}

We introduce in this section the basic notions and notations concerning $\SL(2,\mathbb R)$, many of which are folklore or well-known but we collect all together  in order to fix a consistent system of our notations.

The unimodular group $G=\SL(2,\mathbb R)$ is the matrix group
$$\SL(2,\mathbb R) = \left\{\left.g=\begin{pmatrix} a & b\\ c & d \end{pmatrix} \right| a,b,c,d \in \mathbb R, \det g = 1  \right\}.$$
The group has finite center $\bbZ/2\bbZ$.
It complexified group $G_\mathbb C = \SL(2,\mathbb C)$.
The unique maximal compact subgroup $K$ of $G$ is
$$K = \left\{\left. k(\theta) = \pm\begin{pmatrix} \cos \theta & \sin\theta\\ - \sin\theta & \cos\theta\end{pmatrix} \right| \theta\in [0,2\pi)  \right\}.$$
The group is simple with the only nontrivial parabolic subgroup, which is minimal and therefore is a unique split Borel subgroup
$$B = \left\{\left. b= \begin{pmatrix} a & b \\ 0 & d \end{pmatrix} \right| a,b,d\in \mathbb R, ad=1 \right\}.$$
The Borel group $B$ is decomposed into semidirect product of its unipotent radical $N$ and a maximal split torus $T\cong \mathbb R^*_+$ and a compact subgroup $M= \{\pm 1\}$
It is well-known the Cartan decomposition $G=B\rtimes K= BK$
$$\begin{pmatrix} a & b \\ c & d\end{pmatrix} = \begin{pmatrix} y^{1/2} & x\\ 0 & y^{-1/2} \end{pmatrix}\begin{pmatrix} \pm\cos\theta & \pm\sin\theta\\ \mp\sin\theta & \pm\cos\theta \end{pmatrix}$$
and it is easy to compute
$$y = \frac{1}{c^2 + d^2}$$
$$\cos\theta = \pm y^{-1/2}d,\mbox{ or } \theta = \arccos \frac{d}{\sqrt{c^2+d^2}}$$
$$\pm y^{1/2}\sin\theta \pm y^{-1/2}x \cos\theta = b \mbox{ or  } x = \pm \frac{(b-d)}{c},$$
and the Langlands decomposition of the Borel subgroup
$B=  M\ltimes(A \ltimes N) = MAN$.
The Lie algebra of $G=\SL(2,\mathbb R)$ is $\mathfrak g = \sl(2,\mathbb R) = \langle H,X,Y\rangle_\mathbb R$ where
$$H = \begin{pmatrix} 1 & 0\\ 0& -1\end{pmatrix},
X = \begin{pmatrix} 0 & 1\\ 0& 0\end{pmatrix},
Y = \begin{pmatrix} 0 & 0\\ 1& 0\end{pmatrix},$$
satisfying the Cartan commutation relations
$$[H,X]=2X, [H,Y]= -2Y, [X,Y] = H.$$
The Lie algebra of $A$ is $\mathfrak a = \langle H \rangle_\mathbb R$, the Lie algebra of $N$ is $\mathfrak n = \langle X \rangle_\mathbb R$.
The Lie algebra of $B$ is $$\mathfrak b = \mathfrak a \oplus \mathfrak n = \langle H,X\rangle_\mathbb R.$$
The complex Cartan subalgebra of $\mathfrak g$ is a complex subalgebra
$$\mathfrak h = \langle H \rangle_\mathbb C \subset \mathfrak g_\mathfrak C,$$ which is concided with it normalizer in $\mathfrak g_\mathbb C$. The corresponding subgroup of $G$ such that its Lie algebra is a Cartan subalgebra, is called a Cartan subgroup.
The root system of $(\mathfrak g,\mathfrak h)$ is of rank 1 and is
$$\Sigma = \{\pm\alpha\}, \alpha = (1,-1) \in \mathbb Z(1,-1)\subset \mathbb R(1,-1).$$ There is only one positive root $\alpha = (1,-1)$, which is simple. There is a compact Cartan group $T= K = \O(2)$ The  coroot vector is
$H_\alpha = (1,-1)$ and the split Cartan subgroup of $B$ is $H = \mathbb Z/2\mathbb Z \times \mathbb R^*_+ \cong \mathbb R^*$.

\begin{defn}
An endoscopy subgroup of $G=\SL(2,\mathbb R)$ is the connected component of identity in the centralizer of a regular semisimple element of $G$.
\end{defn}

\begin{prop} \label{endoscopy}
The only possible endoscopy subgroups of $G=\SL(2,\mathbb R)$ are itself or $\SO(2)$.
\end{prop}
\beginpf
The regular semisimple elements of $\SL(2,\mathbb R)$ are of the form $g = \diag(\lambda_1,\lambda_2)$, $\lambda_1\lambda_2 =1$. If $\lambda_1 \ne \lambda_2$, the centralizer of $g$ is the center $C(G) = \{\pm 1 \}$ of the group $\SL(2,\mathbb R)$. If $\lambda_1 = \lambda_2$ and they are real, the centralizer of $g$ is the group $\SL(2,\mathbb R)$ itself. If they are complex and their arguments are opposite, we have $g= \pm\diag(e^{i\theta},e^{-i\theta})$. In this case the connected component of identity of the centralizer is $\SO(2)$. The endoscopy groups  $\SL(2,\mathbb R)$ itself or the center $\{\pm 1\}$ are trivial and the only nontrivial endoscopy group is $\SO(2)$.
\endpf

\section{Automorphic Representations}
In this section we make clear the construction of automorphic representations.

The following lemma is well-known.
\begin{lem}
There is a one-to-one correspondence between any irreducible $2n-1$ dimensional representations of $\SO(3)$ and the $n$ dimensional representations of $\SO(2)$.
\end{lem}
\beginpf
There is a short exact sequence
\begin{equation} \CD
1 @>>> \{\pm I\} @>>> \SU(2) @>>> \SO(3) @>>> 1
\endCD\end{equation}
The characters of the $2n-1, n=1,2,\dots$ dimensional representation of $\SO(3)$ is
\begin{equation}
\chi_{2n-1}(k(\theta)) = \frac{\sin (2n-1)\theta}{\sin \theta}
\end{equation}
where $$k(\theta) = \begin{pmatrix}
\cos\theta & \sin\theta\\  -\sin\theta & \cos\theta \end{pmatrix}$$
\endpf

It is well-known the Iwasawa decomposition $ANK$: every element $g= \begin{pmatrix} a & b\\ c & d \end{pmatrix}$ has a unique decomposition of form
\begin{equation}\label{iwasawa}
\begin{pmatrix} a & b\\ c & d \end{pmatrix} =
\begin{pmatrix} y^{1/2} & 0\\  0 & y^{-1/2} \end{pmatrix}
\begin{pmatrix} 1 & y^{-1/2}x\\ 0 & 1 \end{pmatrix}
\begin{pmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{pmatrix}
\end{equation}
of $\SL_2(\mathbb R)$, where
\begin{equation}
N = \left\{ \begin{pmatrix} 1 & x\\ 0 & 1 \end{pmatrix} \right\} \end{equation}
 is the unipotent radical of the Borel subgroup,
\begin{equation}
A = \left\{ \begin{pmatrix} y^{1/2} & 0\\  0 & y^{-1/2} \end{pmatrix}\right\}
\end{equation}
is the split torus in the Cartan subgroup, and
\begin{equation}
K = \left\{\begin{pmatrix} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{pmatrix}\right\} \end{equation}
is the maximal compact subgroup.
\vskip .2cm

To each modular form $f\in\mathcal S_k(\Gamma)$ of weight $k$ on the Poincar\'e plane $\mathbb H = \SL(2,\mathbb R)/\SO(2)$, we associate the automorphic form $\varphi_f \in \mathcal A_k(\SL(2,\mathbb R))$
\begin{equation}
\varphi_f(\begin{pmatrix} a & b\\ c & d \end{pmatrix} ) = y^{k/2}e^{ik\theta} f(x+iy),
\end{equation}
where $x,y,\theta$ are from the Iwasawa decomposition (\ref{iwasawa}).

The discrete series representations are realized on the function on $L^2(\mathbb H,\mu_k)$, $\mu_k = y^k\frac{dxdy}{y^2}$ of weight $k$ modular form by the formula
\begin{equation}
\pi_k(\begin{pmatrix} a & b\\ c & d \end{pmatrix})f(z) = (cz+d)^{-k}f(z). \end{equation}
Denote by $\mathcal D_k$ the discrete series representation of weight $k$.
The cuspidal automorphic representations are realized in the space $L^2_{cusp}(\Gamma\backslash \SL(2,\mathbb R))$ of automorphic forms

\begin{thm}\label{thm1}
The set of interwining $\SL(2, \mathbb R)$ homomorphisms  from the set of discrete series representations to the set $L^2_{cusp}(\Gamma \backslash\SL(2, \mathbb R))$ of the  automorphic representations of
$\SL(2, \mathbb R)$ is equal to the set $S_k(\Gamma)$ of modular forms, i.e.
\begin{equation}
\Hom_{\SL(2, \mathbb R)}(\mathcal D_k, L^2_{cusp}(\Gamma \backslash \SL(2, \mathbb R)) = \mathcal S_k(\Gamma)
\end{equation}
\end{thm}
\beginpf
The theorem is well-known in literature, see for example A. Borel \cite{borel}. Starting from some intertwining operator $A\in \Hom_{\SL(2, \mathbb R)}(\mathcal D_k,L^2_{cusp}(\Gamma \backslash\SL(2, \mathbb R)) $ we may construct the $L$-series which is an element in $S_k(\Gamma)$; and conversely, starting from some a modular form $f\in S_k(\Gamma)$ we construct the corresponding $L$-series $L_f$. There exists a unique interwining operator $A$ such that the $L$-series of which is equal to $L_f$.
\endpf

\subsection{Geometric Langlands Correspondence}
The general Geometric Langlands Conjecture was stated by V. Drinfel'd and was then proven by E. Frenkel, D. Gaitsgory and Vilonen \cite{vilonenetc} and became the Geometric Langlands Correspondence (GLC).
We will start to specify the general GLC in our particular case of group $\SL_2(\mathbb R)$.

\begin{thm}[Geometric Langlands Correspondence]\label{thm2}
There is a bijection
\begin{equation}
[\pi_1(\Sigma), \SO(3)] \longleftrightarrow \mathcal A(\SL(2,\mathbb R))
\end{equation}
between the set  of equivalence classes of representation of the fundamental group $\pi_1(\Sigma)$ of the Riemannian surface $\Sigma = \Gamma\backslash\SL(2,\bbR)/\SO(2)$ in $\SO(3)$ and the set $\mathcal A(\SL(2, \mathbb R))$ of equivalence classes of automorphic representations of $\SL(2, \mathbb R)$.
\end{thm}
\beginpf The theorem was proven in the general context of a reductive group in the works of \cite{vilonenetc}. For the particular case of $\SL(2,\mathbb R)$ many things are simplified, what we want to point out here.

The main idea to prove the theorem is consisting of the following ingredients:
\begin{enumerate}
\item
It is well-known that every element of $\SO(3)$ is conjugate with some element of the maximal torus subgroup $\SO(2)$. Therefore the set of homomorphism from $\pi_1(\Sigma)$ into $\SO(3)$ is the same as the set of character of the Borel subgroup (minimal parabolic subgroup) from which the discrete series representations are induced.
\begin{lem}
There is one-to-one correspondencee between the conguacy classes of $\SO(3)$ in itself and the inducing character for the discrete series representations.
\end{lem}
\item
Moreover, the $\Gamma$ invariance condition is the same as the
 condition to be extended from the local character of some automorphic component. The following two lemmas are more or less known \cite{bailly}.
\begin{lem}\label{conjcl}
Every representation of $\pi_1(\Sigma)$ in $\SO(3)$ is defined by a system of conjugacy classes in $\SO(3)$, one per generator of $\pi_1(\Sigma)$
\end{lem}

\begin{lem}
Every system of conjugacy classes in the previous lemma \ref{conjcl} defines a unique modular form on $\mathbb H$ and hence a unique automorphic representation of $\SL_2(\mathbb R)$.
\end{lem}

\end{enumerate}
\endpf

\subsection{Geometric Quantization}
The idea of realizing the automorphic representations of reductive Lie groups was done in \cite{diep1}. In this section we show the concrete computation for the case of $\SL(2, \mathbb R)$.

\begin{thm}\label{mthm}
The automorphic representations are obtained from the quantization procedure of fields based on geometric Langlands correpondence.
\end{thm}
\beginpf
The discrete series representations can be realized through the geometric quantization as follows.
\begin{enumerate}
\item
The representation space of discrete series representation is consisting of square-intergrable holomorphic functions.
\begin{equation} \label{convseries}
f(z) = \sum_{n=0}^\infty c_nz^n; \quad \sum_{n=0}^\infty |c_n|^2 < \infty \end{equation}
\begin{lem}
The Hardy space of holomorphic quare-integrable functions can be realized as the expnential Fock space of the standard representation of $\SO(2)$ in $\mathbb C$.
\end{lem}
Indeed every module-square convergence series of type ref{convseries} can be express as some element
\begin{equation} \label{convseries}
f(z) = \sum_{n=0}^\infty n!c_n\frac{z^n}{n!}; \quad \sum_{n=0}^\infty |c_n|^2 < \infty \end{equation}
 in the exponential vector space
\begin{equation}\label{exp}
\EXP \mathbb C = \bigoplus_{n=0}^\infty \frac{C^{\otimes n}}{n!},\quad \mathbb C^{\otimes n} \cong \mathbb C
\end{equation}
\item
The Lie algebra $\sl_2(\mathbb R)$ with 3 generators $\sl_2(\mathbb R) = \langle H, X, Y \rangle_\mathbb R$ in the induced representations of discrete series act through the action of one-parameters subgroups
\begin{equation}
g_3(t) = \exp (tH) = \begin{pmatrix} e^t & 0\\ 0 & e^{-t} \end{pmatrix}
\end{equation}
$$g_1(t) = \exp(t(X-Y)) = \begin{pmatrix} \cos t & \sin t\\ -\sin t & \cos t \end{pmatrix}$$
$$g_2(t) = \exp(t(X+Y)) = \begin{pmatrix} \cosh t & \sinh t\\ \sinh t & \cosh t \end{pmatrix}$$
Under the representation $\pi^{\pm}_n$ we have
\begin{equation}
\pi^{\pm}_n(g_k(t)) = e^{it\hat U_k}
\end{equation}
and we have the relations
\begin{equation}
i\hat U_3 = \hat H = 2z\frac{\partial}{\partial z} +(n+1),
\end{equation}
$$i\hat U_1 = \hat X - \hat Y = -(1+z^2)\frac{\partial}{\partial z} - (n+1)z,$$
$$i\hat U_2 = \hat X + \hat Y = (1-z^2)\frac{\partial}{\partial z} - (n+1)z .$$

In this representation, the action of the element
\begin{equation}
\Delta = \frac{-1}{4}(\hat U_1^2 - \hat U_2^2 - \hat U_3^2 )=
\end{equation}

 $$=\frac{-1}{4}\left((\hat X - \hat Y)^2 - (\hat X + \hat Y)^2 - (\hat H)^2\right) = \frac{1}{4} (\hat H^2  + 4\hat X\hat Y)=$$
$$= (z -\bar z)^2\frac{\partial^2}{\partial z\partial \bar{z}}= -y^2\left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y} \right)$$

This action can be represented as the action of the loop algebra over the Riemann surface $\Sigma$ with values in $\SO(2)$, the elements of the loop algebra are presented in form of some formal/conformal Laurent series with values in the corresponding Lie algebra of form of connection appeared in the construction of induced representations, i.e.
\begin{equation} T(z) = \sum_{n=-\infty}^\infty c_nz^n, c_n\in \so(2), z\in\Sigma . \end{equation}
In our case the Lie algebra $\so(2)$ is one dimensional and we have all number coefficients $c_n, n\in \mathbb Z$.
\begin{lem}
The decompositon \ref{exp}  of $\EXP \mathbb C$ presents the weight decomposition of $\sl_2(\mathbb R)$ in which $H$ keeps each component, $X$ acting as creating operator and $Y$ is acting as some annihilating operator.
\end{lem}
\item The lowest weight representations of the loop algebras are realized through the lowest weight representations of the Virassoro algebra as follows. Let us consider the generators
\begin{equation} L_n = \int_\Sigma z^{n+1} T(z)dz \end{equation}
for any element $$T(.): \Sigma \to \SO(2)$$ from the loop algebra presentation are realized in the Fock space of the standard representation.
These generators satisfy the Virasoro algebra relations
\begin{equation}
[L_m,L_n] = (n-m)L_{m+n} + \delta_{n,-m}\frac{n(n^2-1)}{12}L_0,
\end{equation}
where in an irreducible representation, $Z=cL_0= cI$ is the central charge element.
\end{enumerate}
The proof of Theorem \ref{thm2} is therefore achieved.
\endpf

\section{Arthur-Selberg Trace Formula}
\subsection{Trace Formula}
Let us review the trace formula due to Selberg and J. Arthur.
Remind that by $\mathbb H = \SL(2,\mathbb R) \backslash \SO(2)$ denote the upper Poincar\'e half-plane, $$\mathbb H = \{ z= x+iy \in \mathbb C | \Im z = y > 0 \},$$
$\Gamma$ denote a Langlands type discrete subgroup of finite type with finite number of cusps $\kappa_1,\dots,\kappa_h$. Let $$\Gamma_0 =\left\{\left. \gamma = \begin{pmatrix}
1 & b\\ 0 & 1 \end{pmatrix} \right| b \in \mathbb Z  \right\} \subset \SL(2,\mathbb Z)$$ and $$\Gamma_i= \{\sigma \in \Gamma | \sigma \kappa_i = \kappa_i  \} \subset \SL(2,\mathbb Z)$$ and $\sigma_i\in \SL(2,\mathbb R), i=\overline{1,h}$ are such that  $\sigma_i\Gamma_i\sigma_i^{-1} = \Gamma_0$.
$\sigma_i\Gamma_i\sigma_i = \Gamma_0$,
$\mathfrak H = L^2(\Gamma\backslash G)$ denote the space of quare-integrable functions on $G$ on which there is a natural rugular representation $R$ of $G$ by formula
$$[R(\begin{pmatrix} a & b \\ c & d \end{pmatrix})f](z) = f(\frac{az+b}{cz+d}), z \in \mathbb H.$$ In particular, the unipotent radical $N= \left\{\left. \begin{pmatrix} 1 & x\\ 0 & 1 \end{pmatrix} \right| x \in \mathbb R  \right\}$ of $B$ is acting by the translation on varable $z$ to $z +x$. For any funciton $\psi: N \backslash \mathbb H \to \mathbb C$ of variable $y$ decreasing fast enough when $y$ approaches to 0 or $\infty$,  one defines the incomplete $\theta$-series
$$\theta_{t,\psi}(z) = \sum_{\sigma\in \Gamma_i\backslash \Gamma}  \psi(\sigma_i^{-1}\sigma z)$$
which is certainly of class $L^2(N\cap \Gamma\backslash \mathbb H)$. Denote by $$\Theta = \langle \theta_{t,\psi} | \forall \psi, t \rangle\subset L^2(\Gamma \backslash H)$$ the space of incomplete $\theta$-series

It is well-known that the orthogonal complement $\Theta^\perp$ of $\Theta$ in $L^2(\Gamma \backslash \mathbb H)$ is isomorphic to the space $\mathfrak H_0$ of cuspidal parabolic forms or in other words, of automorphic forms with zero Fourier constant  terms of  automorphic representations, $\mathfrak H = \Theta \oplus \mathfrak H_0$.

Consider the Hecke algebra $\mathcal H(\SL(2,\mathbb R))$ of all convolution Hecke operators of form as follows. Let $F : \mathbb H \to \mathbb C$ be a function $K$-invariant with respect to transformations of form $z \mapsto \gamma z$, for all $\gamma \in K$. Such a function is uniquely defined by a function on $K \backslash G \slash K$, or a so called spherical function $F$ on $G$ which is left and right invariant by $K$. Under convolution these functions provide the Hecke algebra of Hecke operators by convolution with functions in representations. All the Hecke operators have kernel as follows.
For any automorphic function $f$. $f(\gamma z) = f(z), \forall \gamma \in \Gamma$,
$$(F*f)(z) = \int_\mathbb H F({g'}^{-1}g)f(g')dg' = \int_\mathbb H k(z,z')f(z')dz'= \int_{\Gamma\backslash \mathbb H} \sum_{\gamma\in \Gamma} k(z,\gamma z')f(z')dz',$$
$$= \int_{\Gamma\backslash \mathbb H} K(z,z')f(z')dz',$$
 where $z = gi, z' g'i, i = \sqrt{-1}$, $$K(z,z')= \sum_{\gamma\in \Gamma}k(z,\gamma z').$$ Denote the sum of kernel over cups  by
$$H(z,z') = \sum_{i=1}^h H_i(z,z') = \sum_{i=1}^h\sum_{\gamma_i\in \Gamma_i\backslash \Gamma} \int_{-\infty}^{+\infty} K(z,\gamma_i n(x)\gamma_i^{-1}\gamma z')dx$$
The Hecke operator with kernel $K(z,z')$ has the same spectrum as the operator with kernel $K^*(z,z') = K(z,z') - H(z,z')$. The kernel $K^*(z,z')$ are bounded and the fundamental domain $\Gamma\backslash \mathbb H$ has finite volume, therefore the kernels $K^*(z,z')$ are of class $L^2$ on $\mathcal D \times \mathcal D, \mathcal D = \Gamma\backslash \mathbb H$, and all the Hecke operators are compact operators. The Hecke operators keeps each irreducible components of $\Theta^\perp$ invariant and therefore are scalar on each automorphic representation.
On each irreducible component, the Laplace operators has also a fixed eigenvalue
$$\Delta f = \lambda f, \lambda = \frac{s(s-1)}{4}, \Delta = -y^2\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}  \right)$$

One deduces therefore the theorem of spectral decomposition for the discrete part of the regular representation.

\begin{thm}[Spectral decomposition] In the induced representation space of $\ind_B^G\chi_{\lambda,\varepsilon}$, choose $$H_n = \left\{f \in H | \pi(\begin{pmatrix} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{pmatrix})f = e^{in\theta}f \right\}.$$. They are all of dimension 1 and $$H = \bigoplus_{n\in \mathbb Z} H_n.$$
The discrete part $R|_{L^2_{cusp}(\Gamma\backslash \SL(2,\mathbb R)}$ of the regular representation can be decomposed as the sum of the discrete series representations $\pi_n^\pm$ in spaces $$D_{s+1}^+=\bigoplus_{n\equiv\varepsilon \mod 2, n \geq m} H_n$$ or $$D_{s+1}^-=\bigoplus_{n\equiv\varepsilon \mod 2, n \leq -m}H_n,$$ $s\in \mathbb Z, s>0$ and $s+1 \equiv \varepsilon \mod 2$ and there exists $m\in \mathbb Z, m=s+1, m>0$, induced from $\chi_{i\lambda,\varepsilon} = |a|^{i\lambda}(\sgn a)^\varepsilon$ and limits of principal series representations $\pi_0^\pm$ in $\mathcal D_1^+$ or $\mathcal D_1^-$ as two component of the representation $\pi^{0,1} = \ind_B^G\chi_{0,1}$,   induced from the character $\chi_{0,1}$.
\end{thm}
\begin{rem}
In the spaces $\bigoplus_{-m<n<m} H_n$ of dimension $2m-1$ the finite dimensional representations $V_m$ are realized.
\end{rem}

\begin{cor}
For any function $\varphi$ of class $C^\infty_c(G)$, the operator $\pi^\pm_n(\varphi)$ is of trace class and is a distribution denoted by $\Theta^\pm_n$ (following Harish-Chandra) which are uniquely defined by their restriction to the maximal compact subgroup $K = \SO(2)$ and
$$\tr R(\varphi) = \sum_{n \in \mathbb Z, n\geq 0, \pm} m(\pi_n^\pm)\Theta^\pm_n(\varphi) ,$$ with multiplicities $m(\pi^\pm_n)$.
\end{cor}

\subsection{Stable Trace Formula}

The Galois group $\Gal(\mathbb C/\mathbb R) = \mathbb Z_2$ of the complex field $\mathbb C$ is acting on the discrete series representation by character $\kappa(\sigma) = \pm 1$. Therefore the sum of characters can be rewrite as some sum over stable classes of characters.
$$\tr R(f) = \sum_{n=1}^\infty \sum_{\varepsilon = \pm 1} (\Theta_n^+(f) - \Theta_n^-(f)).$$
\begin{rem}The stable trace formula is uniquely defined by its restriction to the maximal compact subgroup $K=\SO(2)$ and is
$$S\Theta_n = \frac{\Theta_{n}^+ - \Theta_n^-}{e^{i\theta} - e^{-i\theta}} = -2i\sin\theta(e^{in\theta}-e^{-in\theta}).$$
\end{rem}

\section{Endoscopy}
\subsection{Stable orbital integral}
Let us remind that the \textit{orbital integral} is defined as
$$\mathcal{O}_\gamma(f) = \int_{G_\gamma\backslash G} f(x^{-1}\gamma x)d\dot x$$ The main idea of computation of endoscopy transfer was explained in \cite{labesse}, we make it here in more detail to clarify some points.

The complex Weyl group is isomorphic to $\mathfrak S_2$ while the real Weyl group
is isomorphic to $\mathfrak S_1$ . The set of conjugacy classes inside a strongly regular
stable elliptic conjugacy class is in bijection with the pointed set
$\mathfrak S_2 /\mathfrak S_1= \mathfrak S_2$
 that can be viewed as a sub-pointed-set of the group
$\mathfrak E(\mathbb R, T, G) = (Z_2 )$
We shall denote by $\mathfrak K(\mathbb R, T, G)\cong \bbZ_2$ its Pontryagin dual.

Consider $\kappa \ne 1$ in $\mathfrak K(\mathbb R, T, G)$ such that $\kappa(H) = -1$.
 Such a $\kappa$ is unique.
The endoscopic group $H$ one associates to $\kappa$ is isomorphic to
$\SO(2)$
the positive root of $\mathfrak h$ in $H$ (for a compatible order) being $\alpha = \rho$

Let $f_\mu$ be a pseudo-coefficient for the discrete series representation $\pi_\mu$ then the \textit{$\kappa$-orbital integral} of a  regular element $\gamma$
in $T (R)$ is given by
$$\mathcal O^\kappa_\gamma(f_\mu) =\int_{G_\gamma\backslash G} \kappa(x) f_\mu(x^{-1}\gamma x) d\dot x = \sum_{\sign(w) =1} \kappa(w) \Theta^G_\mu(\gamma^{-1}_w)  = \sum_{\sign(w) =1} \kappa(w)\Theta_{w\mu}(\gamma^{-1}),$$ because there is a natural bijection between the left coset classes and the right coset classes.

\subsection{Endoscopic transfer}
We make details for the guides of J.-P. Labesse \cite{labesse}.

\textsl{The simplest case} is the case when $\gamma = \diag(a, a^{-1})$. In this case, because of Iwasawa decomposition $x=auk$, and the $K$-bivariance,  the orbital integral is
$$\mathcal O_\gamma(f) = \int_{G_\gamma\backslash G} f(x^{-1}\gamma x)dx = \int_U f(u^{-1}\gamma u)du$$ $$=\int_\mathbb R f(\begin{pmatrix} 1 & x\\ 0 & 1\end{pmatrix}^{-1}\begin{pmatrix} a & 0 \\ 0 & a^{-1} \end{pmatrix} \begin{pmatrix} 1 & x\\ 0 & 1\end{pmatrix})dx$$
$$=\int_\mathbb R f(\begin{pmatrix} 1 & -x\\ 0 & 1\end{pmatrix}\begin{pmatrix} a & 0 \\ 0 & a^{-1} \end{pmatrix} \begin{pmatrix} 1 & x\\ 0 & 1\end{pmatrix})dx$$ $$= \int_\bbR f(\begin{pmatrix} a & (a-a^{-1})x\\ 0 & a^{-1}\end{pmatrix})dx = |a-a^{-1}|^{-1} \mathcal O_\gamma(f).$$
The integral is abosolutely and uniformly convergent and therefore is smooth function of $a\in \bbR^*_+$. Therefore the function $$f^H(\gamma) = \Delta(\gamma)\mathcal O_\gamma(f), \quad \Delta(\gamma) = |a- a^{-1}|$$ is a smooth function on the endoscopic group $H= \bbR^*$.

\textsl{The second case} is the case where $\gamma = k_\theta = \begin{pmatrix} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{pmatrix}$. We have again, $x= auk$ and
$$\mathcal O_{k(\theta)}(f) = \int_{G_{k(\theta)}\backslash G} f(k^{-1}u^{-1}a^{-1}k(\theta)auk)dx\frac{dy}{|y|}d\theta$$
$$ = \int_{G_{k(\theta)}\backslash G} f(u^{-1}a^{-1}k(\theta)au)dx\frac{dy}{|y|}d\theta$$
$$ = \int_{G_{k(\theta)}\backslash G} f(\begin{pmatrix} 1 & -x \\ 0 & 1 \end{pmatrix}\begin{pmatrix} a^{-1} & 0\\ 0 & a \end{pmatrix} \begin{pmatrix} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{pmatrix}\begin{pmatrix} a^{-1} & 0\\ 0 & a \end{pmatrix}\begin{pmatrix} 1 & x \\ 0 & 1 \end{pmatrix})dx\frac{dy}{|y|}d\theta$$
$$= \int_1^\infty f(\begin{pmatrix}\cos\theta & t\sin\theta\\ -t^{-1}\sin\theta & \cos\theta \end{pmatrix})|t-t^{-1}|\frac{dt}{t}$$ $$= \int_0^{+\infty}\sgn(t-1) f(\begin{pmatrix}\cos\theta & t\sin\theta\\ -t^{-1}\sin\theta & \cos\theta \end{pmatrix})dt.$$ When $f$ is an element of the Hecke algebra, i.e. $f$ is of class $C^\infty_0(G)$ and is $K$-bi-invariant, the integral is converging absolutely and uniformly. Therefore the result is a function $F(\sin\theta)$.
The function $f$ has compact support, then the integral is well convergent at $+\infty$. At the another point $0$, we develope the function $F$ into the Tayor-Lagrange of the first order with respect to $\lambda = \sin\theta \to 0$
$$F(\lambda) = A(\lambda) + \lambda B(\lambda),$$ where $A(\lambda) = F(0)$ and $B(\lambda)$ is the error-correction term $F'(\tau)$ at some intermediate value $\tau, 0 \leq \tau \leq t$.
Remark that
$$ \begin{pmatrix} \sqrt{1-\lambda^2} & t\lambda\\ -t^{-1}\lambda & \sqrt{1-\lambda^2} \end{pmatrix}=\begin{pmatrix} t^{1/2}& 0\\ 0 &t^{-1/2} \end{pmatrix} \begin{pmatrix} \sqrt{1-\lambda^2} & \lambda\\ \lambda & \sqrt{1-\lambda^2} \end{pmatrix}\begin{pmatrix} t^{-1/2} & 0\\ 0 &t^{1/2} \end{pmatrix} $$ we have
$$B= \frac{dF(\tau)}{d\lambda} =\frac{d}{d\lambda} \left. \int_0^{+\infty} \sgn(t-1)f(\begin{pmatrix} \sqrt{1-\lambda^2} & t\lambda\\ -t^{-1}\lambda & \sqrt{1-\lambda^2} \end{pmatrix})dt \right|_{t=\tau}$$
$$=  \int_0^{+\infty} \sgn(t-1)g(\begin{pmatrix} \sqrt{1-\lambda^2} & t\lambda\\ -t^{-1}\lambda & \sqrt{1-\lambda^2} \end{pmatrix})\frac{dt}{t}, $$
where $g\in C^\infty_c(N)$
and $g(\lambda) \cong O(-t^{-1}\lambda)^{-1}.$ $B$ is of logarithmic growth and $$B(\lambda) \cong \ln(|\lambda|^{-1})g(1)$$ up to constant term,  and therefore is contimuous.
$$A = F(0) = |\lambda|^{-1} \int_0^\infty f(\begin{pmatrix}1 & \sgn(\lambda)u\\ 0 & 1 \end{pmatrix}) du - 2f(I_2) + o(\lambda)$$
Hence the functions
$$G(\lambda) = |\lambda|(F(\lambda) + F(\lambda)),$$
$$H(\lambda) = \lambda(F(\lambda) - F(-\lambda))$$ have the Fourier decomposition
$$G(\lambda) = \sum_{n=0}^N (a_n|\lambda|^{-1} + b_n)\lambda^{2n} + o(\lambda^{2N})$$
$$H(\lambda) = \sum_{n=0}^N h_n\lambda^{2n} + o(\lambda^{2N})$$
Summarizing the discussion, we have that in the case of $\gamma = k(\theta)$, there exists also a continuous function $f^H$ such that $$f^H(\gamma) = \Delta(\gamma) (\mathcal O_\gamma(f) - \mathcal O_{w\gamma}(f))=\Delta(k(\theta))\mathcal {SO}_\gamma(f), $$ where $\Delta(k(\theta)) = -2i\sin\theta$.

\begin{thm}
There is a natural function $\varepsilon : \Pi \to \pm 1$ such that in the Grothendieck group of discrete series representation ring, $$\sigma_G = \sum_{\pi\in \Pi} \varepsilon(\pi)\pi,$$ the map $\sigma \mapsto \sigma_G$ is dual to the map of geometric transfer, that for any $f$  on $G$, there is a unique $f^H$ on $H$
$$\tr \sigma_G(f) = \tr\sigma(f^H).$$
\end{thm}
\textsc{Proof.}
There is a natural bijection $\Pi_\mu \cong \mathfrak D(\mathbb R,H,G)$, we get a pairing
$$\langle .,.\rangle : \Pi_\mu \times \mathfrak k(\mathbb R,H,G)\to \mathbb C.$$ Therefore we have
$$\tr \Sigma_\nu(f^H) = \sum_{\pi\in \Pi_\Sigma} \langle s,\pi\rangle \tr \pi(f).$$
\hfill$\Box$

Suppose given a complete set of endoscopic groups $H = \mathbb S^1 \times \mathbb S^1 \times \{\pm 1\}$ or $\SL(2,\mathbb R) \times \{\pm 1\}$.  For each group, there is a natural inclusion
$$\eta: {}^LH \hookrightarrow {}^LG$$

Let $\varphi: DW_\bbR  \to {}^LG$ be the Langlands parameter, i.e. a homomorphism from the Weil-Deligne group $DW_\bbR = W_\bbR\ltimes \bbR^*_+$ the Langlands dual group,
$\bbS_\varphi$ be  the set of conjugacy classes of Langlands parameters modulo the connected component of identity map. For any $s\in
\bbS_\varphi$,
$\check{H}_s = \Cent(s,\check G)^\circ$ the connected component of the centralizer of $s\in\bbS_\varphi$ we have $\check{H}_s$ is conjugate with $H$.
Following the D. Shelstad pairing $$\langle s, \pi\rangle : \bbS_\varphi \times \Pi(\varphi) \to \bbC$$
$$\varepsilon(\pi) = c(s)\langle s,\pi\rangle.$$
Therefore, the relation
$$\sum_{\sigma\in \Sigma_s} \tr \sigma(f^H) = \sum_{\pi\in\Pi} \varepsilon(\pi) \tr \pi(f)$$
can be rewritten as
$$\widetilde{\Sigma}_s(f^H) = \sum_{s\in\Pi} \langle s,\pi\rangle \tr\pi(f)$$
and
$$\widetilde{\Sigma}_s(f^H) = c(s)^{-1}\sum_{\sigma\in\widetilde{\Sigma}_s} \tr \sigma(f^H) .$$
We arrive, finally to the result
\begin{thm}
$$\tr \pi(f) = \frac{1}{\#\bbS_\varphi}\sum_{s\in\bbS_\varphi} \langle s,\pi\rangle\widetilde{\Sigma}_s(\check{f}^H).$$
\end{thm}

\section{Poisson Summation Formula}
In the Langlands picture of the trace formula, the trace of the restriction of the regular representation on the cuspidal parabolic part is the coincidence of of the spectral side and the geometric side.
\begin{equation}
\sum_{\pi} m(\pi) \hat{f}(\pi) = \sum_{\gamma\in \Gamma\cap H} a^G_\gamma \hat{f}(\gamma) \end{equation}
Let us do this in  more details.

\subsection{Geometric side of the trace formula}
\begin{thm} \label{thm61}
The trace formula for the regular representation of $\SL_2(\mathbb R)$ in the space of cusp forms is deccomposed into the sum of traces of automorphic representations with finite multiplicities is transfered into the modified Poisson summation formula
\begin{equation}
\sum_{\gamma\in \Gamma\cap H} \varepsilon(\gamma)\mathcal O_\gamma(f) = \sum_{\gamma\in\Gamma\cap H} \varepsilon(\gamma)\vol(\Gamma\cap H) \int_{H\backslash G} f(x^{-1}\gamma x)dx
\end{equation}

\end{thm}
\beginpf
It is easy to see that the restriction of the Laplace operator $\Delta$ on the Cartan subgroup $H$ is elliptic and therefore the Cauchy problem for the other variables has a unique solution. The solution is the trace formula for the cuspidal parabolic part of the regular representation.
\begin{equation} \tr R(f)|_{ L^2_{cusp}(\Gamma\backslash G)} = \sum_{\gamma\in \Gamma\cap H} \varepsilon(\gamma)\mathcal O_\gamma(f)  \end{equation}
From another side  we have
\begin{equation}\tr R(f)|_{ L^2_{cusp}(\Gamma\backslash G)}=
\sum_{\gamma\in\Gamma\cap H}  \varepsilon(\gamma)\vol(\Gamma\cap H) \int_{H\backslash G} f(x^{-1}\gamma x)dx$$ $$=\sum_{\gamma\in \Gamma\cap H} \varepsilon(\gamma)\vol(\Gamma\cap H)\mathcal {SO}_\gamma(f^H).
\end{equation}
\endpf

\subsection{Spectral side of the trace formula}
The following result is well-known, see e.g. \cite{ggps}
\begin{thm}[Gelfand - Graev - Piateski-Shapiro] \label{GGPS}
For any compactly supported function $f\in C^\infty_0(\SL_2(\mathbb R))$ the operator $R(f)|_{L^2_{cusp}(\Gamma\backslash \SL_2(\mathbb R))}$ is of trace class and each irreducible component is of finite multiplicity
\begin{equation}
R(f)|_{L^2_{cusp}(\Gamma\backslash \SL_2(\mathbb R))} = \sum_{\pi\in \mathcal A(\SL_2(\mathbb R))} m(\pi) \pi(f),
\end{equation}
where $m(\pi) = \dim_{\mathbb C}\Hom_{ \SL_2(\mathbb R)}(\mathcal D_k, L^2_{cusp}(\Gamma\backslash \SL_2(\mathbb R)))$.
\end{thm}

\subsection{Poisson Summation formula}
Therefore, following the Poisson summation formula we have the equlity of the both sides.  Let us denote by $\chi_k$ the character of $\SO(2)$ that induce  the discrete series representation $\mathcal D_k$ of  $\SL_2(\mathbb R)$. In our case the Cartan subgroup is $\SO(2)$ and we have the ordinary
\begin{lem} \label{PS}
De note the universal covering of $H= \SO(2)$ by $\tilde H = \Spin(2)$
\begin{equation}
\sum_{\gamma\in \Gamma\cap \tilde H}\delta(x+\gamma) = \sum_{(\varepsilon,k)\in\mathbb Z_2 \times \mathbb Z=\hat H}  (\sign x)^\varepsilon\chi_k(x)
\end{equation}
\end{lem}
\beginpf It is the same as ordinary Poisson summation, lifted to the universal covering $\tilde H$ of $\SO(2)$.
\endpf

\begin{thm} \label{thm3}
The trace $\tr R(f)$ of the restriction of the regular representation on the cuspidal parabolic part $ L^2_{cusp}(G)$ is computed by the formulas
\begin{equation} \sum_{\pi\in \mathcal{A}(G)} m(\pi) \pi(f) = \sum_{\gamma\in \Gamma\cap \tilde H\atop \varepsilon(\gamma) = \pm 1}\varepsilon(\gamma) \vol(\Gamma\cap \tilde H)\mathcal O_\gamma(f)=\sum_{\gamma\in \Gamma\cap H} \vol(\Gamma\cap H)\mathcal {SO}_\gamma(f^H)\end{equation}
\end{thm}
\beginpf
The proof a combination of the previous theorems \ref{thm61}, \ref{GGPS} and \ref{PS} therefore is complete.
\endpf


TITLE: determine two-term outer, inner and uniform expansions
QUESTION [0 upvotes]: Consider the equation;
determine the two-term outer, inner and uniform expansions assuming that $0<\epsilon<<1$
$$\epsilon \frac{d^2y}{dx^2}+\frac{dy}{dx}+y=0, \hspace{5mm} y(0)=0,\hspace{2mm}y(1)=1$$
I have worked out the exact solution of why but don't know if that applies here?
I have also found by setting $\epsilon=0$ that the first order outer layer term is $\exp(1-x)$. I know how to work out the first term for the inner layer aswell, but could someone explain how to work out the second terms.

REPLY [2 votes]: To get the first term, you took the smallest powers of $\epsilon$. To get the second terms, take the second smallest terms. I'll go through the outer solution.
Let $y=y_0+\epsilon y_1+O(\epsilon^2)$, so
$$ \epsilon y_0''+\epsilon^2 y_1'+y_0'+\epsilon y_1'+y_0+\epsilon y_1=0. $$
The smallest power of $\epsilon$ is $0$, so the leading order equation is
$$y_0'+y_0=0,$$
the next order equation is $O(\epsilon)$,
$$y_0''+y_1'+y_1=0,$$
and you can continue this as far as you want.
Once you know $y_0$, you can use it to find $y_1$, then use that to find $y_2$ and so on.


text:
\begin{document}

\maketitle

\begin{center}
$*$: authors share first authorship, order selected at random. \\
$**$: authors share senior authorship, order selected at random.
\end{center}

\begin{abstract}
Recurrent neural networks (RNNs) have been successfully applied to a variety of problems involving sequential data, but their optimization is sensitive to parameter initialization, architecture, and optimizer hyperparameters.

Considering RNNs as dynamical systems, a natural way to capture stability, i.e., the growth and decay over long iterates, are the Lyapunov Exponents (LEs), which form the Lyapunov spectrum.

The LEs have a bearing on stability of RNN training dynamics because forward propagation of information is related to the backward propagation of error gradients.

LEs measure the asymptotic rates of expansion and contraction of nonlinear system trajectories, and generalize stability analysis to the time-varying attractors structuring the non-autonomous dynamics of data-driven RNNs.

As a tool to understand and exploit stability of training dynamics, the Lyapunov spectrum fills an existing gap between prescriptive mathematical approaches of limited scope and computationally-expensive empirical approaches.

To leverage this tool, we implement an efficient way to compute LEs for RNNs during training, discuss the aspects specific to standard RNN architectures driven by typical sequential datasets, and show that the Lyapunov spectrum can serve as a robust readout of training stability across hyperparameters.

With this exposition-oriented contribution, we hope to draw attention to this under-studied, but theoretically grounded tool for understanding training stability in RNNs.
\end{abstract}

\vspace*{-1mm}
\section{Introduction}
\vspace*{-1mm}
The propagation of error gradients in deep learning leads to the study of recursive compositions and their stability~\cite{Bengio1994}.

Vanishing and exploding gradients arise from long products of Jacobians of the hidden state dynamics whose norm exponentially grows or decays, hindering training~\cite{Pascanu}.

To mitigate this sensitivity, much effort has been made to mathematically understand the link between model parameters and the the eigen- and singular-value spectra of these long products~\cite{Chen2018,Pennington2017,Poole2016}.

For architectures used in practice, this approach appears to have limited a scope~\cite{Yang2017,zheng2020r} likely due to spectra having non-trivial properties reflecting complicated long-time dependencies within the trajectory.

Fortunately, the theory of dynamical systems has formulated this problem for general dynamical systems, i.e., for arbitrary architectures.

The Lyapunov spectrum (LS) formed by the Lyapunov exponents (LEs) is precisely the measurement that characterizes rates of expansion and contraction of nonlinear system dynamics.

While there are results using LEs in machine learning~\cite{Legenstein:2007th,Pennington:2018tg,Laurent:2016vs}, most only look at the largest exponent, the sign of which indicates the presence or absence of chaos. Additionally, to our knowledge, current machine learning uses of LEs are for the autonomous case,  i.e. , analyzing RNN dynamics in the absence of inputs, a regime at odds with the most common use of RNNs as input-processing systems.

Thus, there is likely much to be gained from both the study of other features of the LS and from including input-driven regimes into LEs analysis. Indeed, a concurrent work coming from theoretical neuroscience takes up this approach \cite{Engelken} and outlines the ongoing work using LEs to study models of neural networks in the brain. However, the machinery to compute the LS and the theory surrounding its interpretation is still relatively unknown in the machine learning community.

In this exposition-style paper, we aim to fill this gap by presenting an overview of LEs in the context of RNNs, and discuss their usefulness for studying training dynamics. We present encouraging results which suggest that LEs could serve in efficient hyperparameter tuning by allowing early detection of performance.

To compute LEs, we present a novel algorithm for various RNN architectures, taking advantage of modern deep learning environments.

Finally, we highlight key directions of research that have the potential to leverage LEs for improvement of RNNs robustness and performance.

\vspace*{-1mm}
\section{Motivation and Definitions}
\vspace*{-1mm}
Here we develop the problem of spectral constraints for robust gradient propagation. For transparent exposition, in this section we will consider the "vanilla" RNN defined as
\begin{align}
    \mathbf{o}_t=\mathbf{W}\mathbf{h}_t,\;\mathbf{h}_{t}&=\phi(\mathbf{a}_t),\; \mathbf{a}_t=\mathbf{V}\mathbf{h}_{t-1}+\mathbf{U}\mathbf{x}_{t}+\mathbf{b}\;,\label{eq:vanilla}
\end{align}
where $\mathbf{V}$ is the recurrent weight matrix that couples elements of the hidden state vector $\mathbf{h}_t\in\mathbb{R}^{N}$, $\mathbf{U}$ projects the input, $\mathbf{x}_t$, into the network, $\mathbf{b}$ is a constant bias vector, and $\phi$ applies a nonlinear scalar transformation element-wise. The readout weights, $\mathbf{W}$, output the activity into the output variable, $\mathbf{o}_t$. The loss over $T$ iterates is $L=\sum_{t=1}^{T}L_t$, with $L_t=f(\mathbf{y}_t,\hat{\mathbf{y}}_t)$, with $f$ some scalar loss function (e.g., cross-entropy loss,  $f(\mathbf{y},\hat{\mathbf{y}})={\mathbf{y}}\cdot\log\hat{\mathbf{y}}$), $\hat{\mathbf{y}}_t$ the prediction (e.g., $\hat{\mathbf{y}_t}=\mathrm{softmax}(\mathbf{o}_t)$), and $\mathbf{y}_t$ is a one-hot binary target vector. The parameters, $\Theta=(\mathbf{V},\mathbf{U},\mathbf{b},\mathbf{W})$, are learned by following the gradient of the loss, {\it e.g.} in the space of recurrent weights $\mathbf{V}$,
\begin{align}
    \nabla_{\mathbf{V}}L=\sum_{t=1}^T\sum_{i=1}^N \frac{\partial L}{\partial h_{t,i}}\nabla_\mathbf{V}h_{t,i}=\sum_{t=1}^T \textrm{diag}(\phi'(\mathbf{a}_t))\nabla_{\mathbf{h}_t}L \;\mathbf{h}_{t-1}^\top\;,
\end{align}
where $\textrm{diag}(\mathbf{x})$ is the diagonal matrix formed by the vector $\mathbf{x}$ and $\phi'$ is the derivative of $\phi$ and $^\top$ denotes transpose. Here,
\begin{align}
    \nabla_{\mathbf{h}_t}L=\sum_{s=t}^T\left(\prod_{r=t+1}^s \mathbf{J}^\top_r\right)\mathbf{W}^\top\nabla_{\mathbf{o}_s}L\;,\label{eq:gradloss}
\end{align}
where $\nabla_{\mathbf{o}_s}L$ is some simple expression ({\it e.g.} $\hat{\mathbf{y}}-\mathbf{y}_t$ for cross-entropy loss) and  $\mathbf{J}_t=\frac{\partial \mathbf{h}_{t}}{\partial \mathbf{h}_{t-1}}$ is the Jacobian of the hidden state dynamics,
\begin{align}
    \mathbf{J}_{t}&=\textrm{diag}\left(\phi'(\mathbf{a}_t)\right)\mathbf{V}\;.
\end{align}
$\mathbf{J}_t$ varies in time with $\mathbf{x}_t$ and $\mathbf{h}_{t-1}$ via $\mathbf{a}_t$ and so is treated as a random matrix with ensemble properties arising from the specified input statistics and the emergent hidden state statistics.

Gradients can vanish or explode with $T$ according to the singular value spectra of the products of Jacobians appearing in Eq.~\eqref{eq:gradloss}~\cite{Bengio1994}. This fact motivated the study of hyperparameter constraints that control the average of the latter over the input distribution. Specifically, the scale parameter of Gaussian or orthogonal initializations of $\mathbf{V}$ is chosen such that the condition number (ratio of largest to smallest singular value) of $\mathbf{J}_t$ is bounded around 1 for all $t$ \cite{Pascanu}, or such that the first and second moments of the distribution of squared singular values (averaged over inputs) of the Jacobian products are chosen near 1 and 0, respectively according to dynamical isometry~\cite{Pennington2017}. Parameter conditions for the latter have been derived for i.d.d. Gaussian input. The approach has been extended to RNNs under the assumption of untied weights with good empirical correspondence in vanilla and minimalRNN~\cite{Chen2018}, but larger discrepancy in LSTMs~\cite{Gilboa}, suggesting there are other contributions to stability in more complex, state-of-the-art architectures~\cite{Yang2017,zheng2020r}.

Better understanding how robust gradient propagation emerges in these less idealized settings that include high-dimensional and temporally correlated input sequences demands a more general theory. These Jacobian products appearing in Eq.~\eqref{eq:gradloss} can each be expressed as the transpose of a forward sequence, $\prod_{r=t+1}^s \mathbf{J}^\top_r=\left( \mathbf{J}_{s}\cdots\mathbf{J}_{t+1}\right)^\top$. Note that the backward time gradient dynamics shares properties with the forward time hidden state dynamics. This is an important insight that suggests the use of advanced tools from dynamical systems theory may help better understand gradient propagation.
\begin{figure}[!t]
\centering{}
\vspace*{-1mm}
\includegraphics[width= \linewidth]{fig1_spectra.png}

\caption{{\it Distinguishing stability regimes of forward dynamics using Lyapunov spectra}. (a-c) Example unit activations from simulations computed from 10 initial conditions and 2 input sequences ($\mathbf{x}_t$ (light) and $\mathbf{x}_t'$ (dark)) for variance of weight matrix $\mathbf{V}$ set to gives 3 examples spanning the 3 qualitatively distinct stability regimes based on the sign of $\lambda_{max}$: stability, $\lambda_{max}<0$ (a; blue,  $\sigma^2_\mathbf{V}=1/500$); marginal stability, $\lambda_{max}\approx0$ (b; gray, $\sigma^2_\mathbf{V}=100$); and chaos $\lambda_{max}>0$ (c; green,  $\sigma^2_\mathbf{V}=500$).
(d) Lyapunov spectra computed from (a-c) (same colors). (e) First 20 exponents from (d). (vanilla RNN with $\phi=\tanh$; $\mathbf{V}$ orthogonal;  $\mathbf{U}=\mathbb{I}$; 2 realizations (light/dark) of $\mathbf{x}$ Gaussian, $\mathbf{\Lambda}_\mathbf{x}=\sigma^2_\mathbf{x}\mathbb{I}$ with $\sigma^2_\mathbf{x}=0.6$;  10 Gaussian initial conditions, $\mathbf{\Lambda}_{\mathbf{h}_0}=\sigma^2_{\mathbf{h}_0}\mathbb{I}$ with $\sigma^2_{\mathbf{h}_0}=1$).}
\label{fig:untrained}
\vspace*{-2mm}
\end{figure}
\section{Background}
Here we describe the stochastic Lyapunov exponents from the ergodic theory of non-autonomous dynamical systems and outline their connection to the conditions that support gradient-based learning in RNNs. Jacobians are linear maps that evolve state perturbations forward along a trajectory, $\mathbf{u}_{t}=\mathbf{J}_{t}\mathbf{u}_{t-1}$, {\it e.g.} for some initial perturbation $\mathbf{h}_0(\epsilon)=\mathbf{h}_0+\epsilon \mathbf{u}_0$ with $\epsilon\ll 1$ and a given vector $\mathbf{u}_0$. Thus, perturbations can vanish or explode under the hidden state dynamics according to the singular value spectrum of $\mathbf{T}_t=\mathbf{J}_t\cdots \mathbf{J}_1$. The linear stability of the dynamics is obtained taking $\epsilon\to0$ and then $t\to\infty$. In particular, the Lyapunov exponents, $\{\lambda_i\}_{i=1}^N$, are the exponential growth rates associated with the singular values of $\mathbf{T}_t^{1/t}$ for $t\to\infty$, {\it i.e.} the logarithms of the eigenvalues of $\lim_{t\to\infty}\left(\mathbf{T}_t^\top\mathbf{T}_t\right)^\frac{1}{2t}$. As a property of the stationary dynamics, the Lyapunov exponents are independent of initial condition for ergodic systems, i.e. those with only one or the same type of attractor. If the maximum Lyapunov exponent $\lambda_{\textrm{max}}$ is positive, the stationary dynamics is chaotic and small perturbations explode, otherwise it is stable and small perturbations vanish. The theory was initially developed for autonomous dynamical systems, where stable dynamics implies limit-cycle or fixed-point attractors. How the shape of the Lyapunov spectrum varies with network model parameters has provided unprecedented insight into autonomous ( evolving in the absence of inputs) neural network models in theoretical neuroscience \cite{Monteforte2010,Engelken,PuelmaTouzel2015a}.

RNN dynamics are non-autonomous because the hidden state dynamics $\mathbf{h}_t$ are driven by inputs $\mathbf{x}_t$. The theory of {\it random dynamical systems} generalizes stability analyses to non-autonomous dynamics driven by an input sequence sampled from a stationary distribution \cite{arnold2013random}. Analytical results typically employ uncorrelated Gaussian inputs, but the framework is expected to apply to a wider range of well-behaved input statistics. This includes those with finite, low-order moments and finite correlation times like character streams from written language and sensor data from motion capture systems. The time course of the driven dynamics depends on the specific input realization but, critically, the theory guarantees that the stationary dynamics for all input realizations share the same stability properties, which will in general depend on the input distribution (e.g., its variance). Stability here is quantified using the same definition of Lyapunov exponents in the autonomous case (now called {\it stochastic Lyapunov exponents}; we hereon will omit stochastic for brevity). In Fig.~\ref{fig:untrained}, we show vanilla RNN dynamics in stable, marginally stable, and chaotic regimes, as measured by the sign of the largest LE. For stable driven dynamics, the stationary activity on the time-dependent attractor (called a random sink) is independent of initial condition. This holds true also in the marginally stable case, $\lambda_{max}\approx0$, where in addition there are directions along which the sizes of perturbations (error gradients) are maintained over many iterates. For chaotic driven dynamics, the activity variance over initial conditions fluctuates in time in a complex way. Understanding stability properties in this non-autonomous setting is at the frontier of analysis of neural network dynamics in both theoretical neuroscience~\cite{Engelken,Lajoie2013,Hennequin2012} and machine learning~\cite{NIPS2019_9513,Liu}.

The practical calculation of the Lyapunov spectra (see Sec.~\ref{sec:praccalc}) is fast and offers more robust numerical behaviour and generality over singular value decomposition. Together, this suggests the Lyapunov spectra as a useful diagnostic for RNN sequence learning.

There are a handful of features of the Lyapunov spectrum that dictate gradient propagation and thus, influence training speed and performance. The maximum LE determines the linear stability. The mean exponent determines the rate of contraction of full volume elements. The LE variance measures heterogeneity in stability across different directions and can reflect the conditioning of the product of many Jacobians. We use the first two in the analysis of our experiments later in the paper.

We remind readers that much of the theory regarding high-dimensional dynamics is derived in the so-called thermodynamic limit $N\to\infty$. Here, the number of exponents diverges and the spectrum over $i/N$ becomes stationary (so called extensive dynamics), i.e., it is insensitive to $N$ so long as $N$ is large enough (e.g. \cite{Monteforte2010,Engelken}). Here, self-averaging effects can take hold and enable accurate analytical results based on Gaussian assumptions justified by central limit theorem arguments. `Large enough' $N$ is typically in the hundreds (though this should be checked for any given application) and so these results can be useful for studying RNNs.

Finally, we note that obtaining the Lyapunov exponents from their definition relates to the singular values of $\mathbf{T}_t^{1/t}$ for large $t$, and thus is numerically impractical. The standard approach \cite{Benettin1980} is to exploit the fact that $m$-dimensional volume elements grow at a rate $\lambda^{(m)}=\sum_{i=1}^m\lambda_i$ and so the desired rates,  $\lambda_1=\lambda^{(1)}$, $\lambda_2=\lambda^{(1)}-\lambda_1$, ... arise from volumes obtained by projecting orthogonal to subspaces of increasing dimension. These rates are a direct output of computationally efficient orthonormalization procedures such as QR-decomposition. Next, we discuss the practical aspects of computing the Lyapunov exponents in data-driven RNNs.
\vspace*{-1mm}
\section{Lyapunov Spectrum Estimation for RNNs}\label{sec:praccalc}
\vspace*{-1mm}
We extend this framework of Lyapunov exponents to recurrent neural networks by calculating the asymptotic trajectories of the hidden states of the networks when driven by the same signal.

\subsection{Algorithm description}
We adopt the well-established standard algorithm for computing the Lyapunov spectrum~\cite{Benettin1980,Dieci:1995gp}.
We choose the driving signal to be sampled from fixed-length sequences of the test set.
For each input sequence in a batch, a matrix \textbf{Q} is initialized as the identity, and the hidden states, \textbf{h}, are initialized as zero. To track the expansion and contraction of \textbf{Q} at each step, the Jacobian of the hidden states is calculated and then applied to a set of vectors of \textbf{Q}. The expansion of each vector is calculated and the matrix \textbf{Q} is updated using the QR decomposition at each step. If $r_i^t$ is the expansion of the $i^{th}$ vector at time step $t$ -- corresponding to the $i^{th}$ diagonal element of \textbf{R} in the QR decomposition-- then the $i^{th}$ Lyapunov exponent $\lambda_i$ resulting from an input signal of length $T$ is given by:
\begin{equation}\label{eq:Lyap_sum}
    \lambda_i = \frac{1}{T}\sum_{t=1}^T {\log}(r_i^t)
\end{equation}

The Lyapunov exponents resulting from each input $x^j$ in the batch of input sequences are calculated in parallel and then averaged.
For our experiments, the Lyapunov exponents were calculated over 100 time steps with 10 different input sequences. The mean of the 10 resulting Lyapunov spectra is reported as the spectrum. An example calculation of the Lyapunov spectrum is shown in Fig.~\ref{fig:Convergence}.

The algorithm for computing the Lyapunov exponents is described in Algorithm \ref{algo: lyap}.

\begin{algorithm}[H]
    \SetAlgoLined
    \For{$\mathbf{x}^j$ in \text{Batch}}{
    initialize \textbf{h}, \textbf{Q}\\
    \For{t = 1 $\rightarrow$ T}{
    \textbf{h} $\leftarrow$ f($\textbf{h}$, $\mathbf{x}_t^j$)\\
    J $\leftarrow \frac{d\textbf{f}}{d\textbf{h}}$\\
    Q $\leftarrow$ J$\cdot$Q\\
    Q, R $\leftarrow$ $qr$(Q)\\
    $\gamma_i$ += log(R$_{ii}$)
    }
    $\lambda^j_i = \gamma^j_i/T$
    }
    $\lambda_i$ = mean$_j(\lambda_i^j)$
    \caption{Lyapunov Exponents Calculation}
    \label{algo: lyap}
\end{algorithm}

\begin{figure}[t!]
	\centering
\vspace*{-1mm}
	\begin{subfigure}{.4\linewidth}
	   \includegraphics[width= .8\linewidth]{LE_Convergence.png}
        \caption{}
    \end{subfigure}
    \begin{subfigure}{.4\linewidth}
	   \includegraphics[width= .8\linewidth]{LE_Errors.png}
	   	\caption{}
	\end{subfigure}
\caption{{\it Convergence of Lyapunov spectrum estimator}. (a) Lyapunov exponent estimation as a function of time, as shown in Equation \eqref{eq:Lyap_sum}, where $T$ is the number of iterations. Each line represents a different exponent in the spectrum. By 100 iterations, it is clear that each exponent is changing very little as a function of iteration number. (b) The mean Lyapunov spectrum over the set of input sequences. Standard deviation is shown by the red bars.}
\label{fig:Convergence}
\vspace*{-2mm}
\end{figure}
\vspace*{-1mm}
\section{Experiments}
\vspace*{-1mm}
To illustrate the computation and application of Lyapunov exponents, we applied the approach to two tasks with distinct task constraints: character prediction from sentences, where performance depends on capturing low-dimensional and long temporal correlations, and signal prediction from motion capture data, where performance relies also on signal correlations across the many dimensions of the input.

\subsection{Task details} \label{sec: tasks}
For character prediction, we use Leo Tolstoy's \textit{War and Peace} as the data and follow the character-level language modelling task outlined by Karpathy et al. in \cite{karpathy2015visualizing}. For signal prediction, we use the CMU motion-capture dataset and pre-process it using the procedure outlined in \cite{Li_2018_CVPR}. Both tasks require predicting the next step in a sequence given the preceding input sequence.
For the character prediction task (CharRNN), we use input sequences of 100 characters. For the motion capture task (CMU Mocap), we use input sequences of 25 time steps.

For each task, we use a single-layer LSTM architecture with weight parameters initialized uniformly on $[-p, p]$, where $p$ is referred to as the initialization parameter. For further details on the implementation of both tasks, see Supplementary Materials and in addition the Gihub repository\footnote{https://github.com/shlizee/lyapunov-hyperopt} for available code and ongoing work.

\subsection{Algorithm convergence properties}
We find that the general shape of the spectrum is reached early in training. In Fig. \ref{fig:Inits}, we use the CharRNN task as an example to show that the spectrum rapidly changes in the first few epochs of training, quickly converging after a small number of training epochs to a spectrum near that of the trained network. The reduced mean-squared difference shows a $\mathcal{O}(t^{-1})$ convergence with learning epoch $t$, while the mean difference shows this convergence is from below. This was true across the range of initialization parameters tested.
\begin{figure}[!t]
\centering
\vspace*{-1mm}
	\begin{subfigure}{.19\linewidth}
	   \centering
	   \caption{}
	   \includegraphics[width= \linewidth]{Spectra_epoch0.png}
    \end{subfigure}
    \begin{subfigure}{.19\linewidth}
	   \centering
	   \caption{}
	   \includegraphics[width= \linewidth]{Spectra_epoch4.png}

	\end{subfigure}
	\begin{subfigure}{.19\linewidth}
	   \centering
	   \caption{}
	   \includegraphics[width= \linewidth]{Spectra_epoch15.png}
	\end{subfigure}
	\begin{subfigure}{.19\linewidth}
	   \caption{}
	   \centering
	   \includegraphics[width= \linewidth]{Spectra_MSELoss.png}
	\end{subfigure}
	\begin{subfigure}{.19\linewidth}
	   \caption{}
	   \centering
	   \includegraphics[width= \linewidth]{Spectra_AverageDiff.png}
	\end{subfigure}
\caption{{\it Evolution of the Lyapunov spectra over training}. For values of initialization parameter given in legend in (a) ($p = 0.04$ (Blue), $p = 0.08$ (Orange), $p=0.12$ (Green), $p=0.16$ (Red), $p=0.2$ (Purple), $p=0.5$ (Brown)). The spectra before training (a) and after four (b) and fifteen (c) epochs. The reduced mean-squared error between successive epochs (d) shows that the each spectrum converges quickly during training, with $t^{-1}$ plotted as a dashed black line for reference. Similarly, the difference between the mean of the Lyapunov exponents at the fifteenth epoch and the mean at previous epochs (e) rapidly approaches zero throughout training.}
\label{fig:Inits}
\vspace*{-1mm}
\end{figure}
\subsection{Performance efficiency relative to training}
The calculation of the Lyapunov spectrum is very efficient relative to the computation required for training, as long as the user has an analytical expression for the Jacobian of the recurrent layer (as we describe the expressions for GRU and LSTM in supplementary materials). The average time required for calculating the spectrum for a network by evolving 10 input sequences for 100 time steps relative to the average time required to train those networks on the same device is shown in Table~\ref{table: speed} for the two tasks we consider (c.f. section~\ref{sec: tasks}). The mean number of training units (epochs or iterations) needed to reach the minimum validation loss is calculated for each task, and the total training time is compared to the Lyapunov spectra calculation time. For both tasks, the time required to calculate the spectrum was equivalent to about $4\%$ of the total training time.

\begin{table}[h!]
  \caption{{\it LE calculation rate relative to training}}
  \label{table: speed}
  \centering
  \begin{tabular}{lrrr}
    \toprule
    \multicolumn{1}{c}{} &Training Units& LE Calculation & Fraction of Total \\
    Example  & Needed & Time (training units) & \multicolumn{1}{c}{Training Time}\\
    \midrule
    CharRNN    & $26 \pm 10$ epochs & $1.04\pm .04$ epochs & \multicolumn{1}{c}{$0.04\pm 0.02$}  \\
    CMU Mocap     & $2.4 \pm .5 \times 10^3$ iters & $104\pm 3$ iters &  \multicolumn{1}{c}{$0.04\pm 0.01$} \\

    \bottomrule
\end{tabular}
\end{table}

Further improvements to the computation time of the exponents can be made by making slight modifications to Alg.~\ref{algo: lyap}. First, we can "warm-up" the hidden states \textbf{h} and \textbf{Q} onto the attractor of the dynamical system without the computational cost of the QR decomposition, allowing a faster convergence of the spectrum. A further improvement can be made by not taking the QR decomposition every time step during the calculation of the exponents, but instead every $T_{ON}$ steps. Since the QR step is the most expensive computation in the algorithm, the increase in speed over orthonormalizing every step is approximately $T_{ON}$. However, since increasing this interval leads to greater expansion and contraction of the vectors of \textbf{Q} before orthonormalization steps, this can lead to a spurious plateau at higher indices due to the accumulation of rounding errors. However, if one cares only about the first few exponents, this effect is negligible for reasonable selections of $T_{ON}$ (see supplementary materials for details about warmup and effect of increasing $T_{ON}$).

\vspace*{-3mm}
\subsection{Lyapunov spectrum as a robust readout of training stability}
In general, the dependence of the Lyapunov spectrum on hyperparameters is tangled. To direct our exploration of this dependence and how it relates to performance, we used the task constraints to guide our interpretation of spectra behavior of trained networks from randomly sampled hyperparameters.

For CharRNN, we hypothesized from existing work that to satisfy the constraint of propagating a scalar signal over many iterates, it would be sufficient to have a single exponent approaching 0 with the other exponents more negative~\cite{Henaff2015}. Focusing here on hyperparameters, we uniformly sampled a fixed range of log-dropout and learning rate while keeping the initialization parameter fixed. In Fig.~\ref{fig:performance}(a,b), we observe the spectra and indeed find a correlation between maximum Lyapunov exponent and validation loss. Networks with a maximum LE closer to zero indeed performed better. Of these 60 spectra, we selected a smaller subset of 6 that spanned the range of the gross shapes. Interestingly, they also spanned the range of losses, suggesting there is a consistent signal about the loss encoded in the complex, yet systematic variation of the maximum Lyapunov exponent with hyperparameters. The gross shape of these spectra correlated less with the loss, supporting our hypothesis. A straightforward interpretation of the signal we have observed is that the slower the activity in the mode associated with this maximum LE, the longer gradient information can propagate backward in time. This would give the dynamics more predictive power and thus better performance, but it remains to be verified.
\begin{figure}[t!]
\vspace*{-1mm}
	\begin{subfigure}{.49\linewidth}
	   \centering
	   \includegraphics[width= \linewidth]{CharLE_Summary.png}
    \end{subfigure}
    \begin{subfigure}{.49\linewidth}
	   \centering
	   \includegraphics[width= \linewidth]{SeqLE_Summary.png}
	\end{subfigure}
\caption{{\it Correlation between Lyapunov spectra and performance}. Shown are spectra (a,c) and loss dependence on spectra (b,d) for networks trained on CharRNN (a,b) and CMU MoCap (c,d). Those in (a,b) are for 60 random combinations of learning rate and log-dropout rate, while those in (c,d) are for 30 values of the initialization parameter, $p$. The colors in (a,b) and (c,d) are 6 samples selected by eye to roughly span the range of gross spectra shapes seen in (a) and (c), respectively. All other spectra are shown in gray.
}
\label{fig:performance}
\vspace*{-3mm}
\end{figure}
Next, we focused on the motion capture prediction task. Here, the higher dimensional data means that the changes that training induces in the shape of the Lyapunov spectra will likely be more complicated. We hypothesized that nevertheless there would be more information in the gross shape of the spectra of the dynamics of trained networks and less in the maximum LE. Here, we vary the initialization parameter from 0.05 to 0.5 while keeping the dropout and learning rates fixed. Indeed, we find (Fig.~\ref{fig:performance}c,d) consistent variation of the gross shape with loss. We have quantified the variation of loss with gross shape using the mean Lyapunov exponent, though the spectra vary with higher order features as well.

The plots of Fig.~\ref{fig:performance} illustrate that distinct properties of the LE spectrum correlate with validation performance, depending upon task structure.

For the CharRNN task, the maximum exponent being larger (closer to 0), corresponds to better performance. For the Mocap learning task, the larger mean Lyapunov exponent correspond to better performance.
Please see  Supplementary Materials for other spectrum statistics of all networks tested.

Given our observation that the spectrum changes most rapidly early in training, we believe the most strongly-correlated metrics for a task could allow an alternative way to asses performance of a network early in training, thus serving as a hyper-parameter tuning tool.

\vspace*{-1mm}
\section{Conclusion \& Discussion}
\vspace*{-1mm}
In this paper, we presented an exposition and example application of Lyapunov exponents for understanding training stability in RNNs. We motivated them as a natural quantity related to stability of dynamics and useful in a natural complementary approach to existing mathematical approaches for understanding training stability focused on the singular value spectrum. We adapted the standard algorithm for their computation to RNNs, and showed how it could be made more efficient in standard machine learning development environments. We demonstrated that even the basic implementation runs almost 2 orders of magnitude faster than typical training time (see~\cite{Engelken} for further precision and efficiency considerations) and that it converges relatively quickly with learning time. The latter implies it could be useful as a readout of performance early in training. To test this application, we studied it in two tasks with distinct constraints. In both we find interpretable variation with performance reflecting these distinct constraints.

There are a few points to be made about LEs in RNNs. First, one should consider carefully the degrees of freedom in the architecture under consideration. We have considered the hidden state dynamics of an LSTM such that there are as many Lyapunov exponents as there are units. LSTMs have gating variables that can also carry their own dynamics and could be considered as degrees of freedom, adding another $N$ exponents to the spectrum. While we were able to obtain evidence from the spectra of hidden states alone, recent work \cite{Can2020,Engelken} suggests that studying the stability in the subspace of gating variables can also be informative. A second point is that one should also consider the role of the input weights $\mathbf{U}$ (c.f. Eq. \ref{eq:vanilla}). The strength of inputs (e.g. input signal variance) that drive high-dimensional dynamics has long been known to have a stabilizing effect~\cite{Molgedey1992,Lajoie2013,Schuecker2018}. Thus, the input statistics of the task's data can play a role that is modulated by $\mathbf{U}$, making the latter a target for gradient-based algorithms aiming to decrease loss. Due to limited space, we have not analyzed these weights here, but believe it is important to do so in applications.
Finally, we raise perhaps the most glaring adaptation: that of approximating this asymptotic quantity to settings of finite-length sequences. We proposed to use averaging over inputs, in addition to averaging over initial conditions, as a way to take advantage of batch tensor mechanics to achieve faster LE estimates. This is in contrast to long-time averages used in the LE definition, although our method is justified under assumptions of ergodicity and short transient dynamics.
Analyzing the influence of input batch size on this precision is a topic of future research.

The task analysis we have performed suggests a use case for the Lyapunov spectrum as having features that serve as an early readout of performance useful to speed up hyperparameter search. We looked the maximum and mean Lyapunov exponent as two features highlighted by hypotheses we made based on our knowledge of the task constraints. Presumably, there are others relevant here and in other tasks. The search for generic features that do not require knowledge of task constraints (e.g. dynamical isometry) is important, as is demonstrating how useful these features are in practise with complex sequence data. One limitation is how quickly they converge with learning. For example, as a composite quantity, the mean LE converges rather quickly, while the maximum LE is significantly slower to converge (see supplemental material). In \cite{Engelken}, it is shown that loss function/learning rule combinations can also sculpt the spectra and alter its statistics, demonstrating that these are more sources of variation to understand. Interestingly, Full Force~\cite{Depasquale2018}, an alternative to Back-propagation-through-time, has a qualitatively similar spectra (same max, min, and mean LE) but with much lower variance, presumably a desirable feature for generalization. This method also demands more information about performance than just a gradient, suggesting the interesting hypothesis that the precision with which the features of the Lyapunov spectra can be sculpted can scale with the amount of information provided in training. Finally, recent work has given theoretical grounding to why LSTMs avoid vanishing gradients~\cite{Can2020}. Of course this was the reason for the design of LSTM so it only recapitulates the original design intuition. Looking forward however as more complex architectures are developed with less interpretable design, this approach based on LEs can still provide novel insight into why they work.

We close with a short discussion of open problems as this is a prominent area in understanding network dynamics at the intersection of machine learning and computational neuroscience. Having supporting analytical results is essential for robust control of complex problems. To this end, understanding how and where the assumption of untied weights breaks down and how forward propagation differs from backward propagation will be important in extending analytical results on spectral constraints. The Lyapunov spectrum can guide this work. Last, we believe insightful parallels in theoretical neuroscience work. For example, Lyapunov spectrum have been derived for additional degrees of freedom in the unit dynamics, different connectivity ensembles, and more. Also, a discrepancy between the loss of linear stability and the onset of chaotic dynamics in driven systems must be understood~\cite{Schuecker2018}. Making these connections explicit will serve both fields.

\medskip

config: arxiv


TITLE: Is there a geometric projection for every complex function?
QUESTION [2 upvotes]: I was wondering about the best way to visualize complex functions. As they're $$ {\mathbb R}^2 \rightarrow {\mathbb R}^2\ ,$$ I think best way are complex plane image/grid transforms like they used in the Dimensions movie (part 6) or here. Now, is there is a geometric surface (which you could plot with 3dplot) for every complex function which when projected renders the grid transform? (For $$z \rightarrow\ z + k/z $$ this would probably have to cut or overlap itself.) Also what is the further relationship between geometry and complex numbers?
And maybe someone knows what software they used for the animations. Is there any sw along these lines able to visualize your own algebras (vs. functions), i.e. not just $(ac - bd), (ad + bc)$ for multiplication.

REPLY [1 votes]: It seems you have software available that produces 3D-plots from parametric represenations of surfaces, as
$$(\phi,\theta)\mapsto\bigl(\cos \phi\cos \theta,\ \sin \phi\cos \theta,\ \sin \theta\bigr)\qquad\Bigl(0\leq \phi\leq 2\pi, \ -{\pi\over2}\leq \theta\leq {\pi\over2}\Bigr)$$
(the unit sphere). You can use this software to get "grid plots" of complex analytic functions $f$ as well. When such an $f$ is given as an expression $f(z)$ write it out separating real and imaginary parts:
$$f(x+iy)=u(x,y)+i v(x,y)\ ,$$
where now $u$ and $v$ are real functions of the real variables $x$ and $y$. Then let your software produce the 3D-plot corresponding to the parametric representation
$$(x,y)\mapsto\bigl(u(x,y),\ v(x,y),\ 0\bigr)\ .$$
In the drawing  $\bigl($choose the viewpoint at $(0,0,\infty)\bigr)$ you will see the images of the gridlines $x={\rm const.}$ and $y={\rm const.}$ Similarly you could arrange things such that you see the images of circles $|z|={\rm const.}$ and halflines ${\rm arg}(z)={\rm const.}$


TITLE: Which functions have a compositional inverse which is the antiderivative?
QUESTION [2 upvotes]: Which functions $f$ of one variable have a compositional inverse which is the antiderivative of $f$?

REPLY [0 votes]: Following up Arthur's suggestion and assuming a form $a x^{b}$ for $f$:
Let $g(x) = f^{-1}(x)$ and let $g'(x) = f(x)$
$$
f = a x^b \Rightarrow g = \frac{a}{b+1} x^{b+1}
$$
I accept I've ignored a possible integration constant for $g$ above.
$$
f (g (x)) = x
\Rightarrow
f (\frac{a}{b+1} x^{b+1}) = x
$$
$$
a \left(
\frac{a}{b+1} x^{b+1}
\right)^{b}
=
x
$$
$$
\frac{a^{b+1}}{(b+1)^{b}}
x^{b(b+1)}
=
x
$$
Which can be solved with two assumptions:
$$
b(b+1) = 1
\Rightarrow
b^{2} + b - 1 = 0
$$
and
$$
\frac{a^{b+1}}{(b+1)^{b}}
=1
\Rightarrow
a = (b+1)^{\frac{b}{b+1}}
$$
The first gives
$$
b =
\frac{-1 \pm \sqrt{5}}{2}
$$
which can be used to obtain $a$.
So this form for $f$ seems solvable.
I wonder what other function forms are possible.


text:
\section{System Model}
\label{sec:sysmodel}
\subsection{System Dynamics}

We consider the linear dynamical system,
\begin{align}
{\widetilde{X}}_{t+1} &= A {\widetilde{X}}_{t} + B U_t + w_t,
& Y_{t} = C {\widetilde{X}}_{t} + v_t, \label{eq:sys_model}
\end{align}
{where ${\widetilde{X}}_t \in \mathbb{R}^n$ is the state of the system at time {$t \in \mathbb{N}$}, $U_t \in \mathbb{R}^m$ is the input to the system at time $t$, $w_t \in \mathbb{R}^n$ is the process noise, $Y_t$ {are the} system observations, and $v_t \in \mathbb{R}^p$ is the observation noise.}
We denote $\widetilde{X}_1^T$ by ${\widetilde{X}}$, ${U_{1}^{T-1}}$ by $U$ and $w_1^{T-1}$ by $w$.
Based on the initial state ${\widetilde{X}}_1$ and target state ${\widetilde{X}}_T$, the controller computes a sequence of inputs that moves the state from initial state ${\widetilde{X}}_1$ to the target state ${\widetilde{X}}_T$ in $T$ time instances.
{We assume that the system uses the obsevations $Y_{1}^{T}$ to optimally estimate the states ${\widetilde{X}}$. The optimal estimates of $\widetilde{X}$ made by the system are denoted by $X$ -- in the case of {\it perfect observation}, {\it i.e.}, noiseless and observable systems, then $X = \widetilde{X}$.}

 \subsection{Communication and Attacker/Defender Models}
At each time instance the system (Alice) transmits information about its state {estimate} to a legitimate receiver, which is referred to as Bob, via a noiseless link.
This situation occurs for example when Bob is remotely monitoring the execution of the system as in {Supervisory Control And Data Acquisition (SCADA)} systems or in the remote operation of drones.

{\noindent \textbf{Attacker Model.}}
A malicious receiver, referred to as Eve, is assumed to eavesdrop on the communication between the system and Bob and is able to receive all transmitted signals. The goal of Eve is to make an estimate that is as close to $X$ as possible: since Bob receives $X$ and makes control decisions with this information, Eve is interested in $X$. {We assume that Eve knows: 1) the encoding/decoding functions used by Alice and Bob, 2) the dynamical system and 3) the controller design. This information automatically implies knowing the prior probability distributions on the input and state vectors. With this information set, we assume that Eve uses the most adversarial eavesdropping strategy: one which minimizes our performance metric (see Section~\ref{sec::bob_eve}).} Eve is assumed to be passive: she does not actively communicate but is interested in learning the system's states from $t=1$ to $T$.\\

{\noindent \textbf{Defender Model.} We assume that Alice and Bob have a shared {$k$-bit} key $K$ which they use to encode/decode the transmitted messages. For a given encoding/decoding function, the assumed Eve adopts the most adversarial eavesdropping strategy (from the perspective of our chosen performance metric). Therefore, we assume that Alice/Bob attempt to design their encoding/decoding functions which optimize this worst-case performance. We elaborate more on that in Section~\ref{sec::bob_eve}.}

 \subsection{Inputs and States Random Process Model}
We assume that both receivers are only aware of the system model, the matrices $A, B, C$ and the statistics of noises. Therefore, from the perspective of the receivers, the input and output sequences have random distributions which depend on $A, B, C$ and the statistics of the noise.
In addition to the process noise $w$, the joint distribution $f(X,U,w)$ depends on {\it i)} the initial and target states, {\it ii)} the control law of the system and {{\it iii)} the state estimation process}. So, even in noiseless systems, $X$ and $U$ possess inherent randomness from a receiver's perspective due to its lack of knowledge about the initial and target states.

\subsection{Encoding Model}
{The system encodes and transmits packets $Z_1^T$ to ensure that Bob is able to accurately receive $X_1^T$, the optimal estimates of the system. To do so,}
the system transmits a packet $Z_t$ at each time step $t$.

{In this work, we use} light-weight memoryless encryption schemes. The $t$-th transmitted packet {is} a function of only the current state estimate and the shared keys, thus,
$Z_t \coloneqq \mathcal{E}_t (X_t, K)$, where  $\mathcal{E}_t$ is the encoding function used at time $t$. We will denote $Z_1^T$ by $Z$.\\
\subsection{Bob/Eve Models of Decoding}
\label{sec::bob_eve}
Bob {noiselessly receives the} transmitted packets from the {system, and decodes them using the shared key. Then, using the decoded information, it generates an estimate of the state of the system at times $t \in [T]$.}

{We require that Bob's estimate is as accurate as Alice's. If we assume that, at time $t \in [T]$, Bob's decoding function is $\Gamma_t\left( Z_1^t,K \right)$, then the previous condition is satisfied by ensuring that $\Gamma_t\left( Z_1^t,K \right) = X_t$ for all $t \in [T]$.}

Similarly, Eve also {receives} all transmissions from the system. {However,} unlike Bob, she does not have the key $K$.

Therefore, Eve's estimate of $X_t$ is $\hat{X}_{t} \coloneqq \Set{\phi}_t \left(Z_1^T\right),  t \in [T]$,
where $\phi_t$ is the decoding function used by Eve at time $t$.

\subsection{Distortion Metrics} We consider a distortion-based security metric which captures how far an estimate is from the actual value. In particular, {our analysis is based on the} Euclidean distance as our distance metric. {However, our analysis can be extended to any $p$-norm,} since other norms are just a constant factor away, \textit{i.e.}, $\|X\|_p \leq n^{\frac{1}{p} - \frac{1}{q}} \|X\|_q$. We assess the performance of Eve as how far its estimate $\hat{X}$, is from Alice's estimate $X$. Formally, for a given time instance $t$ and a transmitted codeword $Z_1^T$, we define the following quantity,
\begin{align}
D(t,Z_1^T) & \coloneqq \mathbb{E}_{X_t|Z_1^T} \left\Vert X_t - \hat{X}_t \right\Vert^2 \stackrel{(a)}{=} \text{tr}\left( R_{X_t|Z_1^T} \right), \label{DE_single}
\end{align}

where \eqref{DE_single} captures the distortion incurred by Eve while estimating $X_t$ for transmitted symbols $Z_1^T$.

Equality in (a) follows because the best (minimizing) estimates of Eve at time $t$ are, $
 \hat{X}_t  = \Set{\phi}_t \left(Z_1^T\right) = \mathbb{E} \left[X_t|Z_1^T \right].$

Note that Bob is required to successfully estimate $X_t$ knowing ${Z_1^t}$ and the key. Therefore, for a given realization of the key, the encoding function can only map one $X_t$ and that key realization to each value of $Z_1^T$. Therefore Eve realizes {that only trajectories from a {particular subset} can be the true trajectory} for a given $Z_1^T$: those are the ones which correspond to each key realization. Therefore, the expectation in \eqref{DE_single} is in fact taken over the randomness in the key {taking into account} {posterior} probabilities {given} $Z_1^T$.
{If Eve does not have observations, the expectation is taken} over $X_t$ with prior distribution and we get $D(t,Z_1^T) = \text{tr}(R_{X_t})$.

As $D(t,Z_1^T)$ is a function of time $t$ and the transmitted sequence $Z_1^T$, we consider two overall distortion metrics: the ``average case" distortion (denoted by $D_E$) where we take expectation over all possible $Z_1^T$ and average out over time; and the ``worst case" distortion (denoted by $D_W$) where we take minimum over all possible $Z_1^T$ and time instances.

\begin{align}
&\begin{array}{ll}
\text{Average} \\ \text{Distortion}
\end{array} - \:\:
D_E \coloneqq \mathbb{E}_{Z_1^T} \left[ \frac{1}{T} \sum\limits_{t=1}^T D(t, Z_1^T) \right]  \label{eq:average_case}\\

&\begin{array}{ll}
\text{Worst Case} \\ \text{Distortion}
\end{array} - \:\:
D_W \coloneqq \min\limits_{Z_1^T} \left[ \min\limits_{t\in [T]} D(t, Z_1^T) \right] \label{eq:worst_case} .
\end{align}

{It is worth to note that the definitions of $D_E$ and $D_W$ in~\eqref{eq:average_case} and~\eqref{eq:worst_case} imply that Eve's state estimation must be associated to a time instance. In other words,
making a random/constant estimate of the state hoping that it matches the actual state at some time will lead to high distortion values.}

Further, $D_W$ can be defined even when there is no prior distribution on  $X_1^T$. However, {to provide} a baseline comparison {with} the case when the adversary has no observations, we assume that $X_1^T$ always have a known prior distribution.

\subsection{Design Goals} Our goal is to choose the encoding and decoding functions, $\mathcal{E}_t$ and $\gamma_t$, so that Bob can decode loselessly while the distortion is maximized for Eve's estimate. In addition, we seek to achieve this with the minimum amount of shared keys $K$. In absence of any observation by Eve, these distortions will be,
\begin{align*}
D_E^{\text{max}} &=  \frac{1}{T} \sum\limits_{t=1}^T \text{tr}(R_{X_t}), & D_W^{\text{max}} & = \min\limits_{t \in [T]}\text{tr}(R_{X_t}).
\end{align*}
These will serve as upper bounds as,
\begin{align}
D_E &\!=\!\frac{1}{T} \mathbb{E}_{Z_1^T}\! \sum\limits_{t=1}^{T}\!\text{tr}(R_{X_t|Z_1^T})   \stackrel{(a)}{\leq} \frac{1}{T} \sum\limits_{t=1}^{T} \text{tr}(R_{X_t}) = D_E^{\text{max}},\label{eq::UpperBound} \\
D_W &\!= \min\limits_{Z_1^T}  \min\limits_{t\in [T]} \text{tr}(R_{X_t|Z_1^T}) \leq   \min\limits_{t\in [T]}\mathbb{E}_{Z_1^T} \left[ \text{tr}(R_{X_t|Z_1^T})\right] \nonumber \\
& \stackrel{(b)}{\leq} \min\limits_{t\in [T]} \text{tr}(R_{X_t}) = D_W^{\text{max}} \label{eq::UpperBoundDW} ,
 \end{align}
{where (a) and (b) follows from $\mathbb{E}_{Z_1^T} \left[  \text{tr}\left( R_{X_t} | Z_1^T  \right) \right]    \leq \text{tr}\left( R_{X_t} \right)$ which follows from the law of total variance.}


TITLE: Limit of sequence $\lim_{n\to\infty}\frac{1+(\sqrt{n}+1)^{3}+2\sqrt{n}}{n+\sin(n)}$
QUESTION [4 upvotes]: This is no homework. It's another task of a sample exam and I'd like to know how to solve it.

Find the limit of $$\lim_{n\to
 \infty}\frac{1+(\sqrt{n}+1)^{3}+2\sqrt{n}}{n+\sin(n)}$$

Both numerator and denominator go towards $\infty$ for $n \rightarrow \infty$. So I have tried using L'Hôpitals rule because we got $\frac{\infty}{\infty}$.
Here I had the first question in mind, shall I simplify everything before I derivate? (In most cases I think that would be very useful...?)
I have tried both ways (simplify before derivate, derivate and then simplify) and both ended up in $\frac{\infty}{\infty}$ even after several derivations.
Here is how I simplified it:
$$\lim_{n\rightarrow \infty}\frac{n\sqrt{n}+3n+5\sqrt{n}+2}{n+\sin(n)}$$
Then I differentiated the numerator and denominator:
$f'(x) = \frac{3}{2}\sqrt{n}+3+\frac{5}{2\sqrt{n}}$, $g'(x) = 1+\cos(n)$
giving
$$\lim_{n\rightarrow \infty}\frac{\frac{3}{2}\sqrt{n}+3+\frac{5}{2\sqrt{n}}}
{1+\cos(n)}$$
Doesn't really help me...

REPLY [9 votes]: Expand the cube and divide top and bottom by $n$ to get $$\begin{align*}\lim_{n \to \infty} \frac{n\sqrt{n} +3n + 5\sqrt{n} + 2}{n + \sin n} &= \lim_{n \to \infty} \frac{\sqrt{n} + 3 + \frac{5}{\sqrt{n}} + \frac{2}{n}}{1 + \frac{\sin n}{n}} \end{align*}$$ from which you can see the numerator diverges off to $\infty$ whilst the denominator tends to $1$. So the limit diverges off to infinity as whole.
To see that the denominator tends to $1$, it suffices to note that $\frac{\sin n}{n} \to 0$ this is demonstrated via $-1 \leq \sin n \leq 1 \Rightarrow -\frac{1}{n} \leq \frac{\sin n}{n} \leq \frac{1}{n}$ and then squeezing.

Intuitively, this is to be expected since the highest power of $n$ in the numerator $n^{3/2}$ grows faster than the highest power of $n$ in the denominator $n^1$ whilst $|\sin n|$ is bounded by $1$. So the numerator grows too wildly for the limit to be finite.


TITLE: A question about proof of Ternary Goldbach Conjecture.
QUESTION [1 upvotes]: Let's recall:
Von Mangoldt function $\Lambda$ is a following function:
$$\Lambda: \mathbb{N}\rightarrow \mathbb{R}$$
$$
\Lambda(n) = \left\{ \begin{array}{ll}
\log p & \textrm{if $n = p^k$ for some prime $p$}\\
0 & \textrm{otherwise}\\
\end{array} \right.
$$
Ternary Goldbach Conjecture states that:
$$W(N):=\sum _{k_1+k_2+k_3 = N}\chi(k_1)\chi(k_2)\chi(k_3)>0$$
For odd numbers greater than $5$. Here $\chi$ is characteristic function for prime numbers.
But here is the problem. In proof of this concjecture we see a quite different expression.
$$G(N) = \sum_{k_1+k_2+k_3 = N}\Lambda(k_1)\Lambda(k_2)\Lambda(k_3)$$
This sum is proved to be sufficiently big.
My question is:

How are $W(N)$ and $G(N)$ connected? How to estimate from below the first expression ($W(N)$) using just $G(N)$?

Thank you in advance.

REPLY [4 votes]: In analytic number theory, when one studies primes, it is customary to study sums involving $\Lambda$ instead of $\chi$. Note that $\Lambda$ is essentially supported on primes. For any function $f$, we have
$$
\sum_{n\leq x} f(n)\Lambda(n)  = \sum_{p\leq x} f(p)\log p + \sum_{\substack{p^k\leq x\\ k\geq 2}} f\left(p^k\right)\log p.
$$
If $f$ is small, say bounded by 1, the second sum is $O(\sqrt{x}\log x)$, which is quite small if the main term is expected to be of size $x$. This is a general theme that occurs when studying sums over primes, namely that the sum over prime powers (with the power at least 2) doesn't give a significant contribution to the estimate. Thus summing $f$ against $\Lambda$ is essentially summing $f$ over primes, with the additional weight $\log p$ (indeed, sums such as the one above are often regarded in the literature as simply "sums over primes"). To address your question specifically, we have
$$
G(N) = \sum_{p_1+p_2+p_3=N} (\log p_1)(\log p_2)(\log p_3) + O\Big(\sum_{\substack{p_1^k+p_2^k+p_3^k=N\\ k\geq 2}} (\log p_1)(\log p_2)(\log p_3)\Big)
$$
By a trivial estimation, the error term is $O\left(N^{3/2}(\log N)^3\right)$. Since the main term is known to be of size $N^2$, it follows that for $N$ sufficiently large, the sign of $G(N)$ is not affected by the error term, and so $G(N) > 0$ for $N$ sufficiently large. Note that if $G(N) > 0$, we certainly must have $W(N) > 0$ as well. That is, one does not bound $W(N)$ using $G(N)$ directly, but rather infers that $W(N) > 0$ from the fact that $G(N) > 0$, as both are non-negative sums over primes.