\section{Introduction}
Black holes are peculiar, mysterious and probably the most exciting objects in our universe.
They arise mathematically from particular solutions of Einstein's equations of General Relativity (GR).
The first solutions describing black holes were found by Karl Schwarzschild in 1915
in his quest to describe stars by spherically symmetric and static solutions of the Einstein equations.
Shortly afterward, it was discovered that these solutions exhibit amazing physical properties,
the most interesting of which is that nothing, not even light, can escape from their interior.
In essence, this property defines a black hole.
Recently, astronomers have found strong evidence that black holes are ubiquitous in the
sense that they lie at the center of most galaxies.
These developments have brought black holes into the spotlight of astronomical research.
This has also inspired much interest in the study of theoretical properties of black holes.
In particular, the problem of stability of black holes is currently a question of high interest.
As Frolov and Novikov put it~\cite[page~143]{frolov+novikov},
this is ``one of the few truly outstanding problems that remain in the field of black
hole perturbations.''
We here report on recent work, partly carried out at the Center of Mathematical Sciences and
Applications at Harvard University, which shows that
rotating black holes are indeed stable under small perturbations (linear stability).

\section{The Kerr Black Hole}
In General Relativity, gravity is described geometrically in terms of
the curvature of space-time. Thus space-time is modelled mathematically
by a four-dimensional Lorentzian manifold~$(\mycal M, g)$ of signature~$(+ \ \!\! - \ \!\! - \ \! - )$
(for a more detailed introduction see~\cite{bull} or the textbooks on GR~\cite{adler, misner, wald, straumann2}). In GR, Newton's gravitational law is replaced by the Einstein equations
\[ R_{jk} - \frac{1}{2}\: R\, g_{jk} = 8 \pi \kappa\, T_{jk} \:, \]
where~$R_{jk}$ is the Ricci tensor, $R$ is scalar curvature, and~$\kappa$ denotes the gravitational constant.
Here~$T_{jk}$ is the energy-momentum tensor which describes the distribution of matter in space-time.

A rotating black hole is described by the celebrated Kerr geometry. It is
a solution of the vacuum Einstein equations discovered in 1963 by Roy Kerr.
In the so-called Boyer-Lindquist coordinates, the Kerr metric is given by
(see~\cite{chandra, oneill})
\[ ds^{2} = \frac{\Delta}{U} \,( dt-a\sin^{2}\vartheta
\,d\varphi)^{2}-U \left( \frac{dr^{2}}{\Delta}+d\vartheta^{2}
\right) -\frac{\sin^{2}\vartheta}{U} \Big( a\,dt-(r^{2}+a^{2})d\varphi
\Big)^{2}, \]
where
\[ U = r^{2}+a^{2}\cos^{2}\vartheta,\qquad \Delta = r^{2}-2Mr+a^{2} \:, \]
and the coordinates~$(t,r, \vartheta, \varphi)$ are in the range
\[ -\infty<t<\infty,\;\;\;
M+{\sqrt{M^{2}-a^{2}}}<r<\infty,\;\;\; 0<\vartheta<\pi,\;\;\;
0<\varphi<2\pi \:. \]
Here the parameters~$M$ and~$aM$ describe the mass and the angular
momentum of the black hole. It is easily verified that in the
case~$a=0$, one recovers the Schwarzschild metric
\[ ds^{2} = \left(1-\frac{2M}{r} \right) dt^{2}-
\left(1-\frac{2M}{r} \right)^{-1} dr^2 -r^{2}(d\theta^{2}+\sin^{2}\theta\,d\varphi^{2})\:. \]

For the Kerr metric to describe a physical black hole we need to assume
that~$M^{2}> a^{2}$, giving a bound for the angular momentum
relative to the mass. In this so-called {\em{non-extreme case}},
the hypersurface
\[ r = r_{1} := M+{\sqrt{M^{2}-a^{2}}} \]
defines the {\em{event horizon}} of the black hole.
For~$r<r_1$, the radial coordinate~$r$ becomes time, whereas~$t$ becomes a spatial coordinate.
Since time always propagates to the future, the event horizon can be regarded as the ``boundary of no escape.''
Here we shall only consider the region~$r > r_1$ outside the event horizon.
The coefficients of the Kerr metric are independent of~$t$ and~$\varphi$,
showing that the space-time geometry is {\em{stationary}} and {\em{axi-symmetric}}.
One of the key features of the Kerr geometry is the existence of an {\em
ergosphere}, a region which lies outside the event horizon and
in which the vector field $\frac{\partial}{\partial t}$
is {\em space-like}.

It is a major open problem whether the Kerr geometry is stable under small perturbations.
This question is of physical relevance because if black holes were unstable, they would
decay and should therefore most probably not be observable in our universe.
A possible scenario is that a rotating black hole loses more and more angular momentum
by emitting gravitational radiation. If this were the case, the Kerr black hole with~$a>0$ would be unstable,
because asymptotically for large times it would settle down to the
static and spherically symmetric Schwarzschild solution with~$a=0$.

When analyzing the stability of the Kerr solution, one must specify which type of
perturbations are considered. Perturbing in the class of vacuum solutions corresponds
to perturbations by {\em{gravitational radiation}}. But one can also consider solutions
with matter, for example by considering the Einstein-Maxwell equations or the Einstein
equations coupled to other matter fields.
For the stability problem, one considers the Cauchy problem for initial data
which is close to a Kerr solution (in a suitable weighted Sobolev norm).
Then the question is whether the resulting solution settles down to a Kerr solution
asymptotically for large time.
Due to the nonlinearity of the Einstein equations, this problem is very difficult.
As a first step, one can linearize the equations around a Kerr solution and
analyze the long-time behavior of the resulting linear wave equations. This
is referred to as {\em{linear stability problem}} for the Kerr black hole.

The problem of linear stability of the Schwarzschild black hole has a long history.
It goes back to the study by Regge-Wheeler~\cite{reggestab} who showed
that an integral norm of the perturbation of each angular mode is bounded uniformly in time.
Decay of these perturbations was first proved in~\cite{friedmanstab}.
For the Kerr black hole, linear stability has been an open problem for many years.
We now state our recent result on the linear stability of the Kerr black hole
and give an outline of the proof.

\section{The Teukolsky Equation and its Separation}
The problem of linear stability of black holes can be stated mathematically as the question of
whether solutions of massless linear wave equations in the Kerr geometry decay in time.
The different types of equations are characterized systematically in the Newman-Penrose
formalism by their spin~$s$, taking the possible values~$s=0, \frac{1}{2}, 1, \frac{3}{2}, 2, \ldots$.
From the physical point of view, the most interesting cases are $s=1$ (Maxwell field)
and~$s=2$ (gravitational waves).
A general framework for analyzing the equations of arbitrary spin
in the Kerr geometry is due to Teukolsky~\cite{teukolsky}, who
showed that the massless equations of any spin can be rewritten as a single wave equation for
a complex scalar field~$\phi$. The Teukolsky equation reads
\begin{align*}
&\bigg( \frac{\partial}{\partial r} \Delta \frac{\partial}{\partial r} - \frac{1}{\Delta}
\left\{ (r^2+a^2)\: \frac{\partial}{\partial t} + a\: \frac{\partial}{\partial \varphi} - (r-M) \,s \right\}^2
- 4 s \: (r+i a \cos \vartheta)\: \frac{\partial}{\partial t} \\
&\qquad + \frac{\partial}{\partial \cos \vartheta} \:\sin^2 \vartheta\: \frac{\partial}{\partial \cos \vartheta}
+ \frac{1}{\sin^2 \vartheta} \left\{ a \sin^2 \vartheta\: \frac{\partial}{\partial t}
+ \frac{\partial}{\partial \varphi} + i s \cos \vartheta \right\}^2 \bigg) \phi = 0 \:.
\end{align*}

We consider the Cauchy problem for the Teukolsky equation. Thus we seek a
solution~$\phi$ of the Teukolsky equation for given initial data
\[ \phi|_{t=0} = \phi_0 \qquad \text{and} \qquad \partial_t \phi|_{t=0} = \phi_1 \:. \]
Being a linear hyperbolic PDE, the Cauchy problem for the Teukolsky
equation has unique global solutions.
Also, taking smooth initial data, the solution is smooth for all time.
Our main interest is to show that solutions decay for large time.
In order to avoid specifying decay assumptions at the event horizon and at spatial infinity,
we restrict attention to compactly supported initial data outside the event horizon,
\begin{equation} \label{initsmooth}
\phi_0, \phi_1 \in C^\infty_0\big((r_1, \infty) \times S^2 \big) \:.
\end{equation}

Since the Kerr geometry is axisymmetric, the Teukolsky equation decouples into separate
equations for each azimuthal mode. Therefore, the solution of the Cauchy problem
is obtained by solving the Cauchy problem for each azimuthal mode and taking the sum of the resulting solutions.
With this in mind, we restrict attention to the Cauchy problem for a single azimuthal mode, i.e.
\begin{equation} \label{fouriermode}
\phi_0(r, \vartheta, \varphi) = e^{-i k \varphi}\: \phi_0^{(k)}(r, \vartheta)\:,\qquad
\phi_1(r, \vartheta, \varphi) = e^{-i k \varphi}\: \phi_1^{(k)}(r, \vartheta)
\end{equation}
for given~$k \in \mathbb{Z}/2$ (if~$s$ is half integer, then so is~$k$).
The main result of~\cite{stable} is stated as follows:
\begin{Thm} \label{thmdecay}
Consider a non-extreme Kerr black hole of mass~$M$ and angular momentum~$aM$
with~$M^2>a^2>0$. Then for any~$s \geq \frac{1}{2}$ and any~$k \in \mathbb{Z}/2$, the solution of the
Teukolsky equation with initial data of the form~\eqref{initsmooth} and~\eqref{fouriermode}
decays to zero in~$L^\infty_\text{\rm{loc}}((r_1, \infty) \times S^2)$.
\end{Thm} \noindent
This theorem establishes in the dynamical setting that the non-extreme Kerr black hole
is linearly stable.

Before outlining the proof of this theorem,
we remark that in the case~$a=0$
of the Schwarzschild geometry, the above result was already obtained in the paper~\cite{schdecay};
this was our starting point for attacking the problem with~$a>0$.

Restricting attention to a fixed azimuthal mode, the Teukolsky equation becomes
\begin{align*}
&\bigg( \frac{\partial}{\partial r} \Delta \frac{\partial}{\partial r} - \frac{1}{\Delta}
\left\{ (r^2+a^2)\: \frac{\partial}{\partial t} -i a k - (r-M) \,s \right\}^2
- 4 s \: (r+i a \cos \vartheta)\: \frac{\partial}{\partial t} \\
&\qquad + \frac{\partial}{\partial \cos \vartheta} \:\sin^2 \vartheta\: \frac{\partial}{\partial \cos \vartheta}
+ \frac{1}{\sin^2 \vartheta} \left\{ a \sin^2 \vartheta\: \frac{\partial}{\partial t}
-ik + i s \cos \vartheta \right\}^2 \bigg) \phi = 0 \:.
\end{align*}
This equation can be further separated. Namely, making the separation ansatz
\[ \phi(t,r, \vartheta) = \frac{1}{\sqrt{r^2+a^2}}\: e^{-i \omega t} \,X(r)\, Y(\vartheta) \:, \]
the Teukolsky equation gives rise to the coupled system of ODEs
\begin{equation} \label{coupled}
{\mathcal{R}}_\omega X(r) = -\lambda X(r) \:,\qquad
\mathcal{A}_\omega Y(\vartheta) = \lambda Y(\vartheta) \:,
\end{equation}
where~$\lambda$ is a separation constant, and the radial operator~${\mathcal{R}}_\omega$ as well as
the angular operator~$\mathcal{A}_\omega$ are given by
\begin{align}
{\mathcal{R}}_\omega &= -\frac{(r^2+a^2)^2}{\Delta} \left(\frac{\partial^2}{\partial u^2}
- \frac{ \partial_u^2 \sqrt{r^2+a^2}}{\sqrt{r^2+a^2}} \right) \notag \\
&\qquad + \frac{1}{\Delta} \Big( -i \omega \,(r^2+a^2) - iak - (r-M) \,s \Big)^2 - 4 i s r \omega + 4 k\, a \omega \label{Rop} \\
\mathcal{A}_\omega &= - \frac{\partial}{\partial \cos \vartheta} \:\sin^2 \vartheta\: \frac{\partial}{\partial \cos \vartheta}
+\frac{1}{\sin^2 \vartheta} \Big( -a \omega \sin^2 \vartheta +k - s \cos \vartheta \Big)^2 \:. \label{Aop}
\end{align}
Here the time separation reflects the fact that the Kerr metric is stationary.
However, as the Kerr metric is only axisymmetric (and not spherically symmetric),
the separation of the $\vartheta$-dependence does not correspond to
a symmetry of space-time. This separability goes back to the discovery by Carter
for the scalar wave equation~\cite{cartersep}.

We point out that, in contrast to the $\varphi$-dependence~\eqref{fouriermode},
the separation of the $\vartheta$-dependence cannot be performed for the initial data.
The reason is that the angular operator~$\mathcal{A}_\omega$ depends on~$\omega$.
Therefore, the decomposition into~$\vartheta$-modes makes it necessary to
restrict attention to fixed~$\omega$, meaning that the time dependence is that of a plane wave.
However, decomposing the wave into a superposition of plane waves~$e^{-i \omega t}$
would make it necessary to know the entire time evolution.

\section{Hamiltonian Formulation and Integral Representations}
In order to analyze the dynamics of the Teukolsky wave, it is useful to work with contour integrals
over the resolvent, as we now outline. In preparation, we write the Teukolsky equation in
Hamiltonian form. To this end, we introduce the two-component wave function
\[ \Psi = \sqrt{r^2+a^2} \begin{pmatrix} \phi \\ i \partial_t \phi \end{pmatrix} \]
and write the Teukolsky equation as
\begin{equation} \label{HPsi}
i \partial_t \Psi = H \Psi \:,
\end{equation}
where~$H$ is a second-order spatial differential operator.
We consider~$H$ as an operator on a Hilbert space~$\H$
with the domain
\[ \mathscr{D}(H) = C^\infty_0 \big( (r_1, \infty) \times S^2, \mathbb{C}^4 \big)\:. \]
We point out that the operator~$H$ is not symmetric on the Hilbert space~$\H$.
However, we choose the scalar product on~$\H$
as a suitable weighted Sobolev scalar product in such a way
that the operator~$H-H^*$ is bounded.

If~$H$ acted on a finite-dimensional vector space, the Cauchy problem for the
equation~\eqref{HPsi} with initial data~$\Psi_0$ could be solved with the Cauchy integral formula by
\begin{equation} \label{Cauchy}
\Psi(t) = -\frac{1}{2 \pi i} \ointctrclockwise_\Gamma e^{-i \omega t}\: \big(H-\omega)^{-1}\: \Psi_0\: d\omega\:,
\end{equation}
where~$\Gamma$ is a contour which encloses all eigenvalues of~$H$
(note that this formula holds for any matrix~$H$, even if it is not diagonalizable).
It turns out that in our infinite-dimensional setting, a similar formula holds.
Namely, using the fact that~$H-H^*$ is a bounded operator,
we prove that the resolvent~$R_\omega:= (H - \omega)^{-1}$
exists if~$\omega$ lies outside a strip enclosing the real axis (see~\cite[Lemma~4.1]{stable}):
\begin{Lemma} \label{lemmaresex}
For every~$\omega$ with
\[ |\im \omega| > c \:, \]
the resolvent~$R_\omega = (H-\omega)^{-1}$ exists and is bounded by
\[ \| R_\omega \| \leq \frac{1}{|\im \omega|-c}\:. \]
\end{Lemma} \noindent
When forming contour integrals, one must always make sure to stay outside
the strip~$|\im \omega| \leq c$, making it impossible to work with closed contours enclosing the spectrum.
However, we can work with unbounded contours in the following way (see~\cite[Corollary~5.3]{stable}):
\begin{Prp} For any integer~$p \geq 1$, the solution of the Cauchy problem for the Teukolsky equation
with initial data~$\Psi|_{t=0} = \Psi_0 \in \mathscr{D}(H)$ has the representation
\begin{equation} \label{propagatorp}
\Psi(t) = -\frac{1}{2 \pi i} \int_C e^{-i \omega t}\: \frac{1}{(\omega + 3 i c)^p} \;
\Big( R_\omega \,\big(H + 3 i c \big)^p \,\Psi_0  \Big)\: d\omega \:,
\end{equation}
where~$C$ is the contour
\begin{equation} \label{Cdef}
C = \big\{ \omega \:\big|\: \im \omega = 2c \big\} \cup
\big\{ \omega \:\big|\: \im \omega = - 2c \big\}
\end{equation}
with counter-clockwise orientation.
\end{Prp} \noindent
Here the factor~$(\omega + 3 i c)^{-p}$ gives suitable decay for large~$|\omega|$
and ensures that the integral converges in the Hilbert space~$\H$.

\section{A Spectral Decomposition of the Angular Teukolsky Operator}
The next task is to use the separation of variables in the integrand of
our contour integral representation~\eqref{propagatorp}.
Regarding the angular equation~\eqref{coupled}
as an eigenvalue equation, we are led to consider
the angular operator~$\mathcal{A}_\omega$ in~\eqref{Aop} as an operator on the
Hilbert space
\[ \H_k := L^2(S^2) \cap \{ e^{-i k \varphi}\: \Theta(\vartheta) \:|\: \Theta : (0, \pi) \rightarrow \mathbb{C} \} \]
with dense domain~$\mathscr{D}(\mathcal{A}_\omega) = C^\infty(S^2) \cap \H_k$.
Unfortunately, the parameter~$\omega$ is not real but lies on the contour~\eqref{Cdef}.
As a consequence, the operator~$\mathcal{A}_\omega$ in~\eqref{Aop} is not symmetric,
because its adjoint is given by
\[ \mathcal{A}_\omega^* = \mathcal{A}_{\overline{\omega}} \neq \mathcal{A}_\omega \:. \]
The operator~$\mathcal{A}_\omega$ is not even a normal operator,
making it impossible to apply the spectral theorem in Hilbert spaces.
Indeed, $\mathcal{A}_\omega$ does not need to be diagonalizable, because
there might be Jordan chains.
On the other hand, in order to make use of the separation of variables, we must decompose
the initial data into angular modes. This can be achieved by decomposing the
angular operator into invariant subspaces of bounded dimension, as is made
precise in the following theorem (see~\cite[Theorem~1.1]{tspectral}):
\begin{Thm} \label{thmmain}
Let~$U \subset \mathbb{C}$ be the strip
\[ |{\mbox{\rm{Im}}}\, \omega| < 3c \:. \]
Then there is a positive integer~$N$ and a family of bounded
linear operators~$Q^\omega_n$ on~$\H_k$
defined for all~$n \in \mathbb{N} \cup \{0\}$ and $\omega \in U$
with the following properties:
\begin{itemize}
\item[(i)] The image of the operator~$Q^\omega_0$ is an $N$-dimensional invariant
subspace of $\mathcal{A}_k$.
\item[(ii)] For every $n\geq1$, the image of the operator~$Q^\omega_n$ is an at most
two-dimensional invariant subspace of $\mathcal{A}_k$.
\item[(iii)] The $Q^\omega_n$ are uniformly bounded in~$\text{\rm{L}}(\H_k)$,
i.e. for all $n \in \mathbb{N} \cup \{0\}$ and~$\omega \in U$,
\[ \|Q^\omega_n\| \leq c_2 \]
for a suitable constant $c_2=c_2(s,k,c)$ (here~$\| \cdot \|$ denotes the $\sup$-norm on~$\H_k$).
\item[(iv)] The~$Q^\omega_n$ are idempotent and mutually orthogonal in the sense that
\[ Q^\omega_n\, Q^\omega_{n'} = \delta_{n,n'}\: Q^\omega_n \qquad \text{for all~$n,n' \in \mathbb{N} \cup \{0\}$}\:. \]
\item[(v)] The $Q^\omega_n$ are complete in the sense that for every~$\omega \in U$,
\begin{equation} \label{strongcomplete}
\sum_{n=0}^\infty Q^\omega_n = \mbox{\rm 1 \hspace{-1.05 em} 1}
\end{equation}
with strong convergence of the series.
\end{itemize}
\end{Thm} \noindent

\section{Invariant Disk Estimates for the Complex Riccati Equation}
In order to locate the spectrum of~$\mathcal{A}_\omega$, we use detailed ODE estimates.
The operators~$Q^\omega_n$ are then obtained similar to~\eqref{Cauchy}
as Cauchy integrals,
\[ Q^\omega_n := -\frac{1}{2 \pi i} \ointctrclockwise_{\Gamma_n} s_\lambda\: d\lambda
\:,\qquad n \in \mathbb{N}_0 \:, \]
where the contour~$\Gamma_n$ encloses the corresponding spectral points,
and~$s_\lambda=(\mathcal{A}_\omega - \lambda)^{-1}$ is the resolvent of the angular operator.
What makes the analysis achievable is the fact that~$\mathcal{A}_\omega$ is an ordinary
differential operator. Transforming the angular equation in~\eqref{coupled} into
Sturm-Liouville form
\begin{equation} \label{5ode}
\left( -\frac{d^2}{du^2} + V(u) \right) \phi = 0 \:,
\end{equation}
(where~$u=\vartheta$ and~$V \in C^\infty((0, \pi), \mathbb{C})$ is a complex potential),
the resolvent~$s_\lambda$ can be represented as an integral operator whose kernel
is given explicitly in terms of suitable fundamental solutions~$\phi^\text{\tiny{\rm{D}}}_L$ and~$\phi^\text{\tiny{\rm{D}}}_R$,
\begin{equation} \label{sldef}
s_\lambda(u,u') = \frac{1}{w(\phi^\text{\tiny{\rm{D}}}_L, \phi^\text{\tiny{\rm{D}}}_R)} \times \left\{
\begin{aligned} \phi^\text{\tiny{\rm{D}}}_L(u)\: \phi^\text{\tiny{\rm{D}}}_R(u') &\quad&& \text{if~$u \leq u'$} \\
\phi^\text{\tiny{\rm{D}}}_L(u')\: \phi^\text{\tiny{\rm{D}}}_R(u) &&& \text{if~$u' < u$}\:,
\end{aligned}  \right.
\end{equation}
where~$w(\phi^\text{\tiny{\rm{D}}}_L, \phi^\text{\tiny{\rm{D}}}_R)$ denotes the Wronskian.

The main task is to find good approximations for the solutions of the
Sturm-Liouville equation~\eqref{5ode} with rigorous error bounds
which must be uniform in the parameters~$\omega$ and~$\lambda$.
These approximations are obtained by ``glueing together'' suitable WKB, Airy and parabolic cylinder functions.
The needed properties of these special functions are derived in~\cite{special}.
In order to obtain error estimates, we combine several methods:
\begin{itemize}
\item[(a)] Osculating circle estimates (see~\cite[Section~6]{tspectral})
\item[(b)] The $T$-method (see~\cite[Section~3.2]{tinvariant})
\item[(c)] The $\kappa$-method (see~\cite[Section~3.3]{tinvariant})
\end{itemize}
The method~(a) is needed in order to separate the spectral points of~$\mathcal{A}_\omega$
(gap estimates).
The methods~(b) and~(c) are particular versions of {\em{invariant disk}} estimates
as derived for complex potentials in~\cite{invariant} (based on previous estimates for real potentials
in~\cite{angular} and~\cite{wdecay}).
These estimates are also needed for the analysis of the radial equation, see Section~\ref{secres} below.
We now explain the basic idea behind the invariant disk estimates.

Let~$\phi$ be a solution of the Sturm-Liouville equation~\eqref{5ode} with
a complex potential~$V$. Then the function~$y$ defined by
\[ y = \frac{\phi'}{\phi} \]
is a solution of the Riccati equation
\begin{equation} \label{riccati}
y'= V-y^2\:.
\end{equation}
Conversely, given a solution~$y$ of the Riccati equation,
a corresponding fundamental system for the Sturm-Liouville equation is obtained
by integration. With this in mind, it suffices to construct a particular approximate
solution~$\tilde{y}$ and to derive rigorous error estimates.
The invariant disk estimates are based on the observation that the Riccati
flow maps disks to disks (see~\cite[Sections~2 and~3]{invariant}).
In fact, denoting the center of the disk by~$m \in \mathbb{C}$ and its radius by~$R>0$, we get the
flow equations
\begin{align*}
R' &= -2 R \; {\mbox{\rm{Re}}}\, m \\
m' &= V - m^2 - R^2 \:.
\end{align*}
Clearly, this system of equations is as difficult to solve as the original Riccati equation~\eqref{riccati}.
But suppose that~$m$ is an approximate solution in the sense that
\begin{align*}
R' &= -2 R \; {\mbox{\rm{Re}}}\, m + \delta R \\
m' &= V - m^2 - R^2 \:+\: \delta m\:,
\end{align*}
with suitable error terms~$\delta m$ and~$\delta R$, then the Riccati flow
will remain inside the disk provided that its radius grows sufficiently fast, i.e.\
(see~\cite[Lemma~3.1]{invariant})
\[ \delta R \geq |\delta m|\:. \]
This is the starting point for the invariant disk method.
In order to reduce the number of free functions, it is useful to
solve the linear equations in the above system of ODEs by integration.
For more details we refer the reader to~\cite{invariant, tinvariant}.

\section{Separation of the Resolvent and Contour Deformations} \label{secres}
The next step is to use the spectral decomposition of the angular operator in Theorem~\ref{thmmain}
in the integral representation of the solution of the Cauchy problem.
More specifically, inserting~\eqref{strongcomplete} into~\eqref{propagatorp} gives
\begin{equation} \label{propagator2}
\Psi(t) = -\frac{1}{2 \pi i} \int_C \;\sum_{n=0}^\infty e^{-i \omega t}\: \frac{1}{(\omega + 3 i c)^p} \;
\Big( R_\omega \,Q^\omega_n \,\big(H + 3 i c \big)^p \,\Psi_0  \Big)\: d\omega \:.
\end{equation}
At this point, the operator product~$R_\omega Q^\omega_n$ can be expressed in
terms of solutions of the radial and angular ODEs~\eqref{coupled} which arise in
the separation of variables (see~\cite[Theorem~7.1]{stable}).
Namely, the operator~$Q^n_\omega$ maps onto an invariant subspace of~$\mathcal{A}_\omega$
of dimension at most~$N$, and it turns out that the operator product~$R_\omega \,Q^\omega_n$ leaves this
subspace invariant. Therefore, choosing a basis of this invariant subspace,
the PDE~$(H-\omega)R_\omega Q^n_\omega = Q^n_\omega$ can be rewritten as
a radial ODE involving matrices
of rank at most~$N$. The solution of this ODE can be expressed explicitly in terms
of the resolvent of the radial ODE. In order to compute this resolvent, it is useful to also transform the radial
ODE into Sturm-Liouville form~\eqref{5ode}.
To this end, we introduce the Regge-Wheeler coordinate~$u \in \mathbb{R}$ by
\[ \frac{du}{dr} = \frac{r^2+a^2}{\Delta} \:, \]
mapping the event horizon to $u=-\infty$. Then the radial ODE takes
again the form~\eqref{5ode}, but now with~$u$ defined on the whole real axis.
Thus the resolvent can be written as an integral operator with
kernel given in analogy to~\eqref{sldef} by
\[ s_\omega(u,v) = \frac{1}{w(\acute{\phi}, \grave{\phi})} \:\times\:
\left\{ \begin{array}{cl} \acute{\phi}(u)\, \grave{\phi}(v) & {\mbox{if~$v \geq u$}} \\[0.3em]
\grave{\phi}(u)\, \acute{\phi}(v) & {\mbox{if~$v < u$\:,}} \end{array} \right. \]
where~$\acute{\phi}$ and~$\grave{\phi}$ form a specific fundamental system
for the radial ODE.
The solutions~$\acute{\phi}$ and~$\grave{\phi}$ are constructed as
Jost solutions, using methods of one-dimensional scattering theory
(see~\cite{alfaro+regge} and~\cite[Section~6]{stable}, \cite[Section~3]{wdecay}).

The next step is to deform the contour in the integral representation~\eqref{propagator2}.
Standard arguments show that the integrand in~\eqref{propagator2} is holomorphic
on the resolvent set (i.e.\ for all~$\omega$ for which the resolvent~$R_\omega$ in~\eqref{propagatorp}
exists).
Thus the contour may be deformed as long as it does not cross singularities of the resolvent.
Therefore, it is crucial to show that the integrand in~\eqref{propagator2} is meromorphic
and to determine its pole structure.
Here we make essential use of Whiting's mode stability result~\cite{whiting}
which states, in our context, that every summand in~\eqref{propagator2} is holomorphic
off the real axis. In order to make use of this mode stability, we need to interchange the
integral in~\eqref{propagator2} with the infinite sum.
To this end, we derive estimates which show that the summands in~\eqref{propagator2}
decay for large in~$n$ uniformly in~$\omega$. Here we again use
ODE techniques, in the same spirit as described above for the angular equation
(see~\cite[Section~10]{stable}).
In this way, we can move the contour in the lower half plane arbitrarily close to the real axis. Moreover,
the contour in the upper half plane may be moved to infinity. We thus obtain
the integral representation (see~\cite[Corollary~10.4]{stable})
\[ \Psi(t) = -\frac{1}{2 \pi i} \sum_{n=0}^\infty\:\lim_{\varepsilon \searrow 0}
\int_{\mathbb{R} - i \varepsilon} \frac{e^{-i \omega t}}{(\omega + 3 i c)^p} \;
\Big( R_{\omega,n}\:Q_n^\omega \,\big(H + 3 i c \big)^p \,\Psi_0  \Big)\: d\omega \:. \]

The remaining issue is that the integrands in this representation might have
poles on the real axis.
These so-called {\em{radiant modes}} are ruled out by a causality argument
(see~\cite[Section~11]{stable}). We thus obtain the following result (see~\cite[Theorem~12.1]{stable}).
\begin{Thm} \label{thmrep}
For any~$k \in \mathbb{Z}/2$, there is a parameter~$p>0$ such that for any~$t<0$, the solution of the Cauchy problem
for the Teukolsky equation with initial data
\[ \Psi|_{t=0} = e^{-i k \varphi}\: \Psi_0^{(k)}(r, \vartheta) \qquad \text{with} \qquad
\Psi^{(k)}_0 \in C^\infty(\mathbb{R} \times S^2, \mathbb{C}^2) \]
has the integral representation
\begin{equation} \begin{split}
\Psi&(t,u,\vartheta, \varphi) \\
&= -\frac{1}{2 \pi i} \:e^{-i k \varphi}\: \sum_{n=0}^\infty \int_{-\infty}^\infty \frac{e^{-i \omega t}}{(\omega + 3 i c)^p} \;
\Big( R^-_{\omega,n} \:Q_n^\omega
\big(H + 3 i c \big)^p \,\Psi_0^{(k)}  \Big)(u, \vartheta)\: d\omega \:,
\end{split} \label{propfinal}
\end{equation}
where~$R^-_{\omega,n} \Psi := \lim_{\varepsilon \searrow 0} \big(R_{\omega-i \varepsilon,n} \Psi)$.
Moreover, the integrals in~\eqref{propfinal} all exist in the Lebesgue sense.
Furthermore, for every~$\varepsilon>0$ and~$u_\infty \in \mathbb{R}$, there is~$N$ such that
for all~$u<u_\infty$,
\begin{equation} \label{biges}
\sum_{n=N}^\infty \int_{-\infty}^\infty \bigg\| \frac{1}{(\omega + 3 i c)^p} \;
\Big( R^-_{\omega,n} \:Q_n^\omega
\,\big( H + 3 i c )^p \,\Psi_0^{(k)} \Big)(u) \bigg\|_{L^2(S^2)} \: d\omega < \varepsilon \:.
\end{equation}
\end{Thm}

\section{Proof of Decay} \label{secradial}
Theorem~\ref{thmdecay} is a direct consequence of the integral representation~\eqref{propfinal}
in Theorem~\ref{thmrep}. Namely, combining the estimate~\eqref{biges} with Sobolev methods,
one can make the contributions for large~$n$ pointwise arbitrarily small.
On the other hand, for each of the angular modes~$n=0,\ldots, N-1$,
the desired pointwise decay as~$t \rightarrow -\infty$ follows from the Riemann-Lebesgue lemma.
For details we refer to~\cite[Section~12]{stable}.

\section{Concluding Remarks} \label{secoutlook}
We first point out that the integral representation of Theorem~\ref{thmrep} is
a suitable starting point for a detailed analysis for the dynamics of the solutions of the Teukolsky equation.
In particular, one can study decay rates (similar as worked out for massive Dirac waves in~\cite{wdecay})
and derive uniform energy estimates outside the ergosphere (similar as for scalar waves in~\cite{sobolev}).
Moreover, using the methods in~\cite{penrose}, one could analyze superradiance phenomena
for wave packets in the time-dependent setting.

Clearly, the next challenge is to prove {\em{nonlinear stability}} of the Kerr geometry.
This will make it necessary to refine our results on the linear problem, for example
by deriving weighted Sobolev estimates and by analyzing the $k$-dependence of our estimates.
Moreover, it might be useful to combine our methods and results with microlocal techniques.

\vspace*{.5em} \noindent \thanks {{\em{Acknowledgments:}}
We would like to thank Shing-Tung Yau for his encouragement and support.
F.F.\ is grateful to the Center of Mathematical Sciences and Applications at
Harvard University for hospitality and support during his two visits from
September 15 until October 15, 2015, and again from March 15 until April 15, 2016.

\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
}
\providecommand{\href}[2]{#2}


\section{Introduction}

In a variety of neurological or musculoskeletal diseases, motor impairment leads to specific gait characteristics. In Parkinson's disease (PD), for example, the resulting symptoms of gait impairment are often characterised by shuffling steps, reduced stride length, or impaired gait initiation \cite{Klucken2013}. In advanced stages, PD patients can even experience spontaneous episodes of gait freezing. These symptoms can lead to a severe reduction in patient mobility and quality of life \cite{Ellis2011} and it is thus important to quantify and treat gait impairments as early as possible. However, the physician's rating of gait quality e.g. in terms of the Movement Disorder Society's revision of the Unified Parkinson Disease Rating Scale (MDS-UPDRS) \cite{Goetz2008} is often biased by subjectiveness.
Moreover, it requires long-time experience by the movement disorder specialist, is thus costly and restricted to the availability of the physician. With an ageing world population (\cite{UN2015,EU2005}) and a resulting increase in prevalence of chronic age-related disorders such as PD, this also requires a change in the healthcare system towards a new era of diagnostics and individualized treatment \cite{Pasluosta2015}.

In the past, this has led to the development of several electronic measurement systems with the aim to extract spatio-temporal gait parameters and aid the physician in an objective gait impairment scoring. These systems include for example camera based, stationary 3D optical tracking of markers affixed to the body \cite{Kressig2004} or electronic walkways with embedded pressure sensors \cite{Givon2009}. Optical tracking systems require time-consuming preparation of the subject, but enable three-dimensional, sub millimetre accurate tracking. Computerized walkways with embedded pressure sensors, however, require no preparation of the subject, but are limited in accuracy. Furthermore, they can only assess gait parameters related to ground contact, i.e. foot angles outside the ground plane can not be tracked. Both of these systems are expensive and need a laboratory environment for gait analysis which prohibits their use in everyday clinical practice or outpatient monitoring.

Unlike stationary systems, inertial sensors placed at the lower extremity enable mobile and unobtrusive gait analysis that is unconstrained by the laboratory setting \cite{Klucken2013}. This type of measurement system is inexpensive and requires little subject-preparation while at the same time enabling estimation of three-dimensional spatio-temporal gait parameters (\cite{Mariani2010,Rampp2015,Aminian2002,Salarian2004,Rebula2013,Trojaniello2014,Ferrari2015}). Although mobile systems are increasingly developed to support clinical diagnostics, monitor medical treatment or provide early indicators for gait impairment, they are still an open field of research \cite{Pasluosta2015}. One of the most crucial components of such systems is accurate estimation of stride-by-stride parameters. Therefore, the purpose of this work is accurate estimation of stride length based on stride-specific inertial sensor data captured at the subject's feet.

\section{Related Work}
There is a growing body of literature regarding stride length estimation from inertial sensors placed at the lower extremity of the human body. The methods covered generally divide into two classes: (Bio)mechanical model based approaches (\cite{Aminian2002,Salarian2004,Salarian2013}) and double-integration approaches (\cite{Veltink2003,Sabatini2005,Mariani2010,Rebula2013,Trojaniello2014,Rampp2015,Ferrari2015}).
\\
\paragraph{Model based approaches}
The model based approaches employ a double pendulum model for the swing phase and an inverse double pendulum for the stance phase. Three joints (hip and knees, ankle is considered fixed) and four segments (shank and thigh) are modelled. Integration of the respective equations of motion is driven by gyroscope data captured at the shank and thigh of the subjects (\cite{Aminian2002,Salarian2013,Salarian2004}). Salarian et al. \cite{Salarian2013} additionally try to reduce the number of sensing units either by duplicating data of the right thigh on the left thigh or by simplifying the model to a single-pendulum with two sensors at the shanks.
Apart from the amount of sensors needed on the patient, the models described in \cite{Aminian2002,Salarian2004} and \cite{Salarian2013} confine the movement-model to the sagittal plane and thus limit their applicability for analysis of impaired gait. Furthermore, specific subject characteristics (segment length etc.) are needed to scale the biomechanical models. One of the major benefits, however, is the intrinsic implementation of biomechanical constraints between the two legs.
\\
\paragraph{Double-integration approaches}
The vast majority of stride length estimation algorithms comes from the class of double-integration methods. This family of algorithms makes use of the stance phase as a zero-velocity update point (ZUPT) to re-initialize the integration process. Due to intrinsic measurement errors in currently available inertial sensors, this is necessary in order to minimize the integration drift \cite{Peruzzi2011}. For one gait cycle, the algorithm consists of sensor orientation estimation from gyroscope measurements in order to convert sensor readings from the local, constantly changing sensor frame to global coordinates. This is followed by gravity cancellation and double-integration of accelerometer readings. Due to aforementioned limitations in currently available mobile sensors, the integration result has to be corrected or de-drifted in a last processing step by enforcing several constrains (e.g. zero-velocity at the start and end of the stride, level floor, etc.). This class of methods is mostly evaluated per leg and biomechanical constrains are hard to implement. However, three-dimensional sensor trajectories can be reconstructed allowing a deeper analysis of gait characteristics or disease-specific alterations.

The contributions on the double-integration approach covered here (\cite{Veltink2003,Sabatini2005,Mariani2010,Rebula2013,Trojaniello2014,Rampp2015,Ferrari2015}) all employ accelerometer and gyroscope data captured at the subjects' feet, Trojaniello et al. \cite{Trojaniello2014} additionally use magnetometer information. Sampling rates are comparable and between 100 and 200 Hz for all studies. The contributions mainly differ in the choice of orientation estimation method. Rampp et al. \cite{Rampp2015} and Mariani et al. \cite{Mariani2010} apply a quaternion-based time integration of angular rates to estimate the sensor orientation, as proposed by Sabatini \cite{Sabatini2005}. Other contributions (\cite{Rebula2013,Trojaniello2014,Ferrari2015}) make use of a Kalman filter to obtain sensor orientations. Integration is mainly done directly, i.e. forward in time. Only Trojaniello et al. \cite{Trojaniello2014} use a direct and reverse integration scheme which spawns two integrations -- one from the start and one from the end of the current stride -- and computes a weighted mean between the two with weights depending on the distance to the next ZUPT. This is said to further resolve integration drifts \cite{Zok2004}.
De-drifting of integrated gravity-free acceleration signals is done in the spatial domain in \cite{Rebula2013,Rampp2015} and \cite{Mariani2010} with a linear de-drifting function.
Trojaniello et al. \cite{Trojaniello2014} use a frequency-based de-drifting approach to suppress the low-frequency drift component in the integrated signals. Another possibility is to combine all steps above in one Kalman filter  estimating orientation, velocity and position from angular rates and acceleration. This is implemented by Ferrari et al. \cite{Ferrari2015}.
\\
\paragraph{Clinical applicability}
Tab. \ref{tab:relatedwork} summarizes the dataset characteristics and results on all related work covered here.
Regarding the mean accuracies and precisions, the majority of methods are similar with accuracies within $1$ cm and precisions around $6 - 8$ cm. There are two exceptions: Rebula et al. \cite{Rebula2013} reach good precision but are outperformed by others in terms of accuracy. The method proposed by Trojaniello et al. \cite{Trojaniello2014} achieves the best accuracy and precision almost resolving the reference precision of $\pm 1.27\,\text{cm}$ \cite{Trojaniello2014}. They are the only ones employing a magnetometer which might greatly improve their orientation estimation. However, this is only true if indoor magnetic field distortions by e.g. ferroconcrete in the proximity of the gait acquisition track are small enough. In a general setting without prior knowledge about earth magnetic field disturbances due to construction materials, this out-rules the use of magnetometers for indoor gait data acquisitions with wearable sensors \cite{DeVries2009}.

This also applies to other studies listed in Tab. \ref{tab:relatedwork}.
Regarding the number and type of subjects, there are major differences in terms of variability within the datasets. Not listed are other quite important characteristics of the study populations, i.e. stride length distribution, diversity of gait alterations, precise age distribution, gender, study protocol or track length used during data acquisition. This is greatly threatening the comparability of methods and their clinical application due to non-representative study populations. It shows the urgent need of a unified evaluation of all stride length estimation methods on the same, clinically relevant benchmark dataset.

\newlength{\doubledaggerlength}
\settowidth{\doubledaggerlength}{\footnotesize$^{\dagger\dagger}$}
\def\hspace{\doubledaggerlength}\null{\hspace{\doubledaggerlength}\null}
\begin{table*}[!ht]
	\centering
	\caption{Overview of related work: Dataset characteristics vs. mean accuracies and precisions}
	\begin{tabular}{lrlrlrr}
		\toprule
		& Subjects & Diagnose &\# Strides & Reference &  Mean acc. $\pm$ prec. & rel. prec. \\\toprule
		Rebula et al. \cite{Rebula2013} 			& 9 \hspace{.2em}young 	& healthy		&5538	& MoCap	&$-1.2\pm3.7\text{ cm}^{\dagger\dagger}$ & $3.2~ \%^{~}$ \\
		Mariani et al. \cite{Mariani2010}			& 10 \hspace{.2em}young		& healthy		&482	& MoCap	&$2.4\pm7.5$ cm\hspace{\doubledaggerlength}\null& $6.8~ \%^{~}$\\\midrule
		Mariani et al. \cite{Mariani2010}	  		& 10 elderly 	& healthy		&492	& MoCap	&$0.7\pm6.1$ cm\hspace{\doubledaggerlength}\null& $6.1~ \%^{~}$\\
		Trojaniello et al. \cite{Trojaniello2014} 	& 10 elderly 	& healthy		&576	& GaitRite	&$0.0\pm 1.9$ cm\hspace{\doubledaggerlength}\null& $1.0~ \%^{~}$ \\
		Trojaniello et al. \cite{Trojaniello2014} 	& 10 elderly 	& Parkinson's	&532	& GaitRite	&$0.1\pm 1.9$ cm\hspace{\doubledaggerlength}\null& $2.0~ \%^{~}$ \\
		Rampp et al. \cite{Rampp2015} 				& 101 elderly 	& geriatric 	&1220	& GaitRite	&	 $-0.3\pm8.4$ cm\hspace{\doubledaggerlength}\null& $10.7~ \%^\dagger$\\
		Ferrari et al. \cite{Ferrari2015} 			& 14 elderly 	& Parkinson's 	&1314	& GaitRite	&	 $-0.2\pm7.0$ cm\hspace{\doubledaggerlength}\null& $5.0~ \%^{~}$\\
		Salarian et al. \cite{Salarian2013} 			& 15 elderly 	& healthy \& impaired 	&229	& MoCap	&	 $-0.8\pm6.6$ cm\hspace{\doubledaggerlength}\null& $6.0~ \%^{\dagger}$\\
		\bottomrule
		\multicolumn{7}{l}{
		\footnotesize $^\dagger$ ~Computed based on mean stride length reported.\hfill
		\footnotesize $^{\dagger\dagger}$ Computed after correspondence with the authors.}
	\end{tabular}
	\label{tab:relatedwork}
\end{table*}

Although there has been much progress on stride length estimation from inertial sensor data over the last years, a variety of clinically relevant phenomena still cannot be resolved due to the achieved measurement precision. These phenomena include for example the decline in stride length with age in healthy controls, the reduction in stride length with disease progression in e.g. PD patients or the effect of medication. Hollman et al. \cite{Hollman2011} report normative data on mean stride length in elderly, healthy controls. Between the age groups 70-75 years and 85+ years, stride length decreases by 20 cm in males and 14 cm in females in the mean values over 294 participants \cite{Hollman2011}. This translates to a reduction around 1-2 cm/year. Resolving this difference with mobile sensors requires the effect to be larger than typically twice the standard deviation of the error distribution (i.e. precision) obtained in the validation study of a stride length estimation method. In order to resolve the yearly decline with age, the precision of a mobile stride length estimation method would therefore need to be around 5 mm.
A second clinical application involves stride length estimation in PD patients. Hass et al. \cite{Hass2012} report normative values on mean stride length in PD patients over the course of the disease as measured by the Hoehn\&Yahr (H\&Y) scale \cite{Hoehn1967}. Between the groups H\&Y < 1.5 (mild) and H\&Y = 3-4 (severe), they observed a reduction in mean stride length of 11 cm for males and 12 cm for females based on data from 310 PD patients \cite{Hass2012}. Following the argument above, mobile stride length estimation methods would require precisions around 5 cm in order to resolve the effect of stride length reduction with disease progression even on the very coarse H\&Y scale ranging from 0-5. Lastly, Bryant et al. \cite{Bryant2011} investigated the change in stride length between the best possible (on) and worst (off) medication state in PD patients. Based on data from 21 PD patients, they found a stride length difference of 18 cm between the on- and off-state. These values help to get a feeling about the ranges of stride length effects between best possible to worst medication levels or mild to severe disease progression. Moreover, the clinician is probably interested in more minute changes in clinical scores or medication levels corresponding to smaller effects on the parameter stride length. To this date, mobile stride length estimation is not yet precise enough for this kind of tasks. From a clinical perspective, the problem thus stays challenging.
\\
\paragraph{Towards a new type of methodology}
Until now, the choice of method for stride length estimation has been based on biomechanical, physical or geometrical reasoning and these approaches have proven to work in technical validation settings. However, one could also aim to learn a regression function linking the raw sensor data of a stride directly with the corresponding stride length. Similar approaches were reported by Aminian et al. \cite{Aminian1995} in 1995 where a two-layer perceptron was used to estimate speed and incline of human walking. Due to the limits on computing power 20 years ago, the neural network structure presented by Aminian et al. \cite{Aminian1995} is very shallow with 10 input nodes, one hidden layer of five units and one output node. A 10 feature parametrization of the accelerometer signal based on physiological and statistical reasoning is used as input to the network.

With the recent advances in deep learning (\cite{LeCun2015,Krizhevsky2012,Hinton2006}) and user-friendly interfaces to such techniques like torch \cite{torch7} or google's TensorFlow library \cite{tensorflow2015-whitepaper}, the parametrization of the input data has become obsolete. By regressing against raw sensor data instead of parametrizations and by employing sufficiently deep architectures, the true potential of the neural network approach can be utilized. Based on a publicly available and clinically relevant benchmark dataset, we describe a method to learn a non-linear relationship between stride-specific sensor data and the corresponding stride length with a convolutional neural network. The network can be trained on big data sets \cite{Chen2014} while at the same time promising excellent individualisation capabilities to the patient at hand \cite{Stober2015}.
We evaluate the approach in a cross-validation and obtain mean accuracy and precision regarding stride length estimation. In order to assess the method's dependency on the data segment provided as input, this is done for three different stride definitions. These are either based on biomechanical events (heel-strike and mid-stance) or contain the relevant information (i.e. the swing phase where stride length is produced) but are not directly related to biomechanical events in the datastream.

The method we propose here adds a completely novel approach to the problem of stride length estimation from inertial sensor data. It is not subject to the zero-velocity assumption that the double-integration approaches need to re-initialize the integration process. Thereby, the current work enables spatial gait parameter estimation in situations where this assumption is violated, e.g. in patients with spasticity. Compared to the long processing chain employed by double-integration methods that results in amplification of errors and requires cumbersome, careful and error-prone tuning of counter measures such as de-drifting, the method proposed here is purely data-driven. There is no need to carefully keep track of and account for the accumulating error during processing.

Our main contributions are: 1) We show feasibility regarding estimation of spatial gait characteristics with deep convolutional neural networks. 2) With that change in methodology, we drop the zero-velocity phase assumption that constrains state-of-the-art double integration methods and limits their applicability.

\section{Methods}
\paragraph{Data Collection and Setup}
We use a benchmark dataset collected by Rampp et al. \cite{Rampp2015} that is publicly available at \url{https://www5.cs.fau.de/activitynet/benchmark-datasets/digital-biobank/} and briefly described here.
For data collection, the Shimmer2R sensor platform \cite{Burns2010} was used. It contains a 3d-accelerometer (range $\pm 6\,\text{g}$) and a 3d-gyroscope (range $\pm 500\, ^\circ$/s) and was attached laterally right below each ankle joint (see Fig. \ref{fig:shoes}). The same shoe model (adidas Duramo 3) was used by all subjects to avoid gait changes due to different shoe characteristics \cite{Menant2009}. Data was captured at $102.4$ Hz and a resolution of $12$ bit. Simultaneously, validation data was acquired with the well established, instrumented pressure mat GAITRite (\cite{McDonough2001,Webster2005}). The sensitive area of the pressure mat was $609.60~ \text{cm} \times 60.96~ \text{cm}$ with a spatial resolution of $\pm 1.27$ cm and a track width of $84$ cm.

\begin{figure}[!h]
		\centering
		\includegraphics[width=.6\columnwidth]{shoeSetup.pdf}
		\caption{Placement of the inertial sensor and axes definition.}
		\label{fig:shoes}
\end{figure}

In total, 116 geriatric inpatients were assessed with this setup at the Geriatrics Centre Erlangen (Waldkrankenhaus St. Marien, Erlangen, Germany). Written informed consent was obtained prior to the gait assessment in accordance with the ethical committee of the medical faculty at Friedrich-Alexander University Erlangen-N\"urnberg (Re.-No. 4208).

Patients performed a detailed geriatric assessment, described in detail by Rampp et al. \cite{Rampp2015}. For the scope of this manuscript, we focus on the free walking test over $10\,\text{m}$ at comfortable walking speed instrumented with the inertial sensors and the GAITRite system. Here, the number of strides acquired from each patient ranges from 7 -- 32 with a mean of 12 strides per patient. After excluding datasets from eight patients due to medical reasons, two due to inertial sensor malfunction and additional five due to measurement errors with the GAITRite system, a total of 101 datasets were left for training and evaluation of the stride length estimation method proposed here. In 54\% of the study population, gait disorders or fall proneness were diagnosed. The other top three diagnoses were heart rhythm disorder (70\%), arterial hypertension (69\%) and coronary artery disease (41\%) all of which are associated with gait and balance disorders \cite{Salzman2010}. Therefore, this dataset constitutes a clinically relevant study population both in terms of the number of subjects and the presence of unpredictable gait alterations.
\\
\paragraph{Preprocessing}
Before the inertial sensor data is fed to the convolutional neural network, we perform a series of preprocessing steps. These include extraction of annotated strides (a segmentation of the continuous signal into strides is already supplied by the dataset) and calibration from raw sensor readings to physical units \cite{Ferraris1995}. Due to different sensor mounting on the shoe, coordinate system transformations are needed to align sensor axes on left and right feet. Furthermore, the signals from accelerometer and gyroscope are normalised w.r.t. the respective sensor ranges and zero-padded to fixed length of 256 samples per stride to ensure equally scaled and fixed size input to the network.
\begin{figure}[h]
		\centering
		\includegraphics[width=.95\columnwidth]{inputSignal.pdf}
		\caption{Exemplary input signal for a stride defined from MS$\to$MS after preprocessing.}
		\label{fig:inputSignal}
\end{figure}

Since we will compare the algorithm on different stride definitions, we need to detect mid-stance (MS) and heel-strike (HS) events in the stride segmentation provided by the dataset and adjust the stride borders accordingly. Detection of these two events within the sensor data is done according to \cite{Rampp2015}. We call the stride definition provided by the dataset msDTW because of the multi-dimensional subsequence dynamic time warping they used for segmentation (for details see \cite{Barth2015,Rampp2015}).

An exemplary, pre-processed input signal for the network for one stride defined from MS$\to$MS is shown in Fig. \ref{fig:inputSignal}.\\

\paragraph{Network Architecture}

The neural network architecture we choose is a two layer convolutional network followed by one fully connected layer and a readout-layer (see Fig. \ref{fig:NetArchitecture}). The main theoretical motivation for this choice is the locality of features gained by this type of architecture \cite{LeCun2015}. The network is thus forced to take the true topology of the input data as multi-channelled, synchronized time series into account and maintain its temporal context. It only considers local connections and does not link	 information from temporally far away regions in the input signal.

\begin{figure}[!h]
\centering
\includegraphics[width=.9\columnwidth]{Architecture.pdf}
\caption{Architecture for the neural network regression approach presented in this paper: Two convolutional layers followed by max-pooling, one fully connected layer and a readout layer contracting the last hidden layer to the single target variable.}
\label{fig:NetArchitecture}
\end{figure}

The preprocessed inertial sensor data $x_{i,j}$ with $i=1\dots L_0$ samples and $j = 1\dots N_0$ channels serves as input to the network. In the first convolutional layer, feature maps are constructed from $x_{i,j}$ by convolution with $k=1\dots N_{1}$ kernels $\psi^{(1)}_{i,j,k}$ of size $i=1\dots L_1$ and $j\in[1,N_0]$. A bias term $b^{(1)}_k $ is added to each feature map before a ReLU activation function \cite{Glorot2011} produces the output activations of the first convolutional layer
\begin{equation}
	\vector c^{(1)}_{k} = \text{ReLU}\left(
						\sum_j\vector\psi^{(1)}_{j,k} * \vector x_{j} +b^{(1)}_{k}
						\right)
	\label{eq:outputNodesConv1}
\end{equation}
with $k\in[1,N_1]$ and $\text{ReLU}(z) = \max(0,z)$.

After each convolutional layer, we insert a max-pooling layer that downsamples the feature maps by taking the maximum in non-overlapping windows of length $r$. The downsampling factor thus is $1/r$. The output from the first max-pooling layer is called $o^{(1)}_{i,j}$ with $i\in [1,L_0/r]\text{ and } j\in[1,N_1]$.

In the second convolutional layer, the same approach is repeated with different kernels and biases. Given the input signal $o^{(1)}_{i,j}$, feature maps are constructed by convolution with $k=1\dots N_2$ kernels $\psi^{(2)}_{i,j,k}$ of length $i=1\dots L_2$ and $j\in[1,N_1]$. A bias term $b^{(2)}_k$ is added before a ReLU activation function produces the output nodes
\begin{equation}
	\vector c^{(2)}_{k} = {ReLU}\left(
						\sum_j\vector\psi^{(2)}_{j,k} * \vector o^{(1)}_j +b^{(2)}_{k}
						\right)
	\label{eq:outputNodesConv2}
\end{equation}
with $k\in[1,N_2]$. After max-pooling, the output $o^{(2)}_{i,j}$ is of size $i\in [1,L_0/r^2]\text{ and } j\in[1,N_2]$.
Similarly, up to $\log_{r} L_0$  convolutional layers followed by max-pooling could be stacked before the input signal is downsampled to size 1.

For the fully connected layer, the output from the last max-pooling layer $o^{(2)}_{i,j}$ is flattened to a one dimensional vector $\tilde o^{(2)}_i$ with $i\in[1, N_2\,L_0/r^2]$. This vector is multiplied with the  weights $W^\text{(fc)}_{i,j}$ with $i=1\dots(N_2\,L_0/r^2)$ and $j=1\dots N_\text{fc}$. A bias term $b^\text{(fc)}_j$ is added before a ReLU activation function produces the output of the fully connected layer
\begin{equation}
	o^\text{(fc)}_j = \text{ReLU}\left(\sum_i W^\text{(fc)}_{i,j} \; \tilde o^{(2)}_{i}  + b^\text{(fc)}_j\right)
	\label{eq:outputNodesHidden}
\end{equation}
with $j\in[1,N_\text{fc}]$.

Finally, a readout layer compresses the $N_\text{fc}$ output nodes of the hidden layer to the single target variable stride length. This is done by multiplication with a weight vector $W^\text{(ro)}_{i,1}$ with $i \in [1,N_\text{fc}]$. A final bias $b^\text{(ro)}$ is added to arrive at the estimate of the target variable
\begin{equation}
	y = \sum_i W^\text{(ro)}_{i,1} o^\text{(fc)}_i + b^\text{(ro)}
	\label{eq:target estimate}
\end{equation}
Similarly, a target vector $(y_1, \dots ,y_n)$ representing a set of stride features could be estimated by changing the dimensionality of weight and bias in the readout layer.

For our application of stride length estimation, the input data is of length $ L_0 = 256$ with $N_0 = 6$ channels. On the first layer of the network, we chose to learn $N_1=32$ filters of length $L_1=30$ samples, followed by $N_2=64$ and $L_2 = 15$ on the second.
Max-pooling is done in non-overlapping windows of length $r = 2$ and the fully connected layer has $N_\text{fc} =1024$ nodes. Given the downsampling factor and the sampling rate, the filter lengths both correspond to approximately $0.29$ s. The theoretical motivation for this choice is to keep the resulting receptive field size constant ($L_1 = r \, L_2 $) while increasing the amount of features maps ($N_2 > N_1$) with network depth.\\

\paragraph{Learning a Regression Function}

Given a training dataset, we now focus on determining the parameters involved in the previously described network architecture, i.e. kernels/weights and biases on each layer.

Training neural networks is commonly posed as an optimization problem regarding a scalar error function (implicitly) depending on the network parameters. This error implements a discrepancy measure between the predicted output and a ground truth reference on the training dataset or subsets thereof. Using back-propagation, weights and biases on all layers are changed with the aim to minimize the error. In practice, only random subsets of the training dataset, called mini-batches, are shown to the optimizer in one iteration of the training loop to speed up the learning phase (stochastic learning).

Given the relative error distribution
\begin{equation}
\epsilon_i = \frac{y_i-y_{i,\text{ref}}}{y_{i,\text{ref}}} \;\;\text{with }i\in[1,N_\text{batch}]
\label{eq:relativeError}
\end{equation}
on such a mini-batch of size $N_\text{batch}$, we define the discrepancy measure to minimize as
$E(\vector\epsilon) = \text{rmsq}(\vector\epsilon)$.
The term $\text{rmsq}(\vector\epsilon)$ is the root-mean-square-error on the mini-batch.

For optimization, we use Adam \cite{Kingma2015} which is a state-of-the-art optimization method for stochastic learning. It shows faster convergence than other stochastic optimization routines on benchmark datasets and we use default settings of $\alpha=1\text{e}^{-3},\beta_1=0.9,\beta_2 = 0.999 \text{ and } \epsilon=1\text{e}^{-8}$ (for details see \cite{Kingma2015}). All weights are initialized randomly by sampling a truncated normal distribution with standard deviation 0.1 and biases are initialized from 0.1. Training is done for a fixed number of 4000 iterations with a mini-batch size $N_\text{batch}=100$ strides.

As a measure against over-fitting we use dropout on the fully connected layer. Dropout effectively samples a large number of thinned architectures on the hidden layer by randomly dropping nodes and their connections during training. This technique has been shown to prevent over-fitting significantly in many use-cases and is superior to weight-regularisation methods \cite{Srivastava2014}. We use a fixed dropout probability of $p_\text{drop} = 0.5$, so every node on the fully connected layer has a 50/50 chance of being dropped during training. During testing, however, the full architecture is used and no connections are dropped.

The network is implemented and trained using google's TensorFlow library \cite{tensorflow2015-whitepaper}.
\\

\paragraph{Evaluation Scheme}
Evaluation of the proposed method is based on a 10-fold cross validation scheme. The strides from 101 patients on the dataset are sorted into training and test set depending on the patient identifier to ensure distinct splits of the dataset.
For each of the three stride definitions msDTW, HS$\to$HS and MS$\to$MS, the model is evaluated in such a cross-validation scheme and the stride length is estimated on the test set in each fold. The predictions from individual folds are then pooled to arrive at average statistics for the current stride definition.
Evaluation statistics include average accuracy $\pm$ precision which correspond to the mean and standard deviation of the (signed) error distribution. Additionally, we state precision relative to the mean stride length, average absolute accuracy $\pm$ precision as well as Spearman correlation coefficients between the reference system and our approach.

In order to assess the learning speed and performance of the convolutional neural network, we compute the training error and the corresponding precision over the training iterations for an exemplary 90/10\% (patient-wise) split of the dataset.

\section{Results}

\paragraph{Training}
On a mobile workstation equipped with an Intel Core i7 processor (8 cores), 16GB of RAM and no GPU acceleration, training takes around 20 min. per fold.

\begin{figure}[!h]
		\centering
		\includegraphics[width=.95\columnwidth]{trainingStats.pdf}
		\caption{Error $E(\vector\epsilon)$ and precision evaluated on the entire training set for an exemplary 90/10\% train/test split of the dataset during training.}
		\label{fig:training}
\end{figure}

For an exemplary fold, Fig. \ref{fig:training} shows the error as well as the precision on the entire training set over the iterations. The fixed number of 4000 iterations is sufficient to reach a stable regime of the error $E(\vector\epsilon)$ and precision evaluated on the entire training dataset. Moreover, we reach the reference precision after training which demonstrates the usability of the described approach for estimating stride length at high precision with a model tuned to the individual. Because we can not formally evaluate the performance of individualized models here due to the lack of data (6 strides on average per foot and patient), the results on the training set serve as a substitute for the individualized modelling case where prior information is available.\\

\FloatBarrier

\paragraph{Stride Length Estimation on Unseen Data}

Tab. \ref{tab:results} lists average statistics from the 10-fold cross validation for each of the three stride definitions. Bland-Altman plots in Fig. \ref{fig:bland-altman-plots} assess the agreement between the two measurements for all three stride definitions. Additionally, Fig. \ref{fig:bland-altman-plots} includes regression lines between the measurement error $y-y_\text{ref}$ and the measurement agreement $\frac{1}{2} (y+y_\text{ref}) $ in order to assess the dependency of the measurement error w.r.t. the stride length.

\begin{figure}[!h]
		\centering
		\includegraphics[width=.95\columnwidth]{CrossVal_MS-MS-strides_BlandAltman}\\\pagebreak
		\includegraphics[width=.95\columnwidth]{CrossVal_HS-HS-strides_BlandAltman}\\
		\includegraphics[width=.95\columnwidth]{CrossVal_msDTW-strides_BlandAltman}
		\caption{Bland-Altman plots from the evaluation of the proposed algorithm in a 10-fold cross validation for all three stride definitions. Mean accuracy and 95\% confidence intervals ($\pm 1.96~\sigma$) as well as a regression line between mean stride length and the measurement error are shown.}
		\label{fig:bland-altman-plots}
\end{figure}

\begin{table*}[!t]
	\centering
	\caption{Results of the deep convolutional neural network approach to stride length estimation for three stride definitions compared to the Gaitrite reference system. Listed are mean stride lengths for both systems, the error of measurement is reported as mean accuracy and precision, relative precision, mean absolute accuracy and precision and Spearman correlation (CC).}
	\begin{tabular}{lccrrcc}
		\toprule
		Stride definition & Our approach & GAITRite &Mean acc. $\pm$ prec. & Rel. prec.  & Mean abs. acc. $\pm$ prec. & CC \\\toprule
		\input{results_MS-MS.csv}\\
		\input{results_HS-HS.csv}\\
		\input{results_msDTW.csv}\\
		\bottomrule
	\end{tabular}

	\label{tab:results}
\end{table*}
The model yields best performance when strides are defined from MS$\to$MS with a mean accuracy and precision of $0.01\pm 5.37$ cm.
The range of precisions achieved for individual patients ranges between $\pm1.54$ cm in the best case (patient P114) and $\pm10.43$ cm in the worst case (patient P106).

Furthermore, the stride length estimation with the proposed approach is robust w.r.t. the stride definition. When trained and evaluated on the other two stride definitions (HS$\to$HS and msDTW), the mean precision differs only marginally from the one achieved on the MS$\to$MS strides. The mean accuracy, however, decreases to $-0.27$ cm and $-0.42$ cm for HS$\to$HS and msDTW strides respectively.

The regression lines that indicate the dependency of the measurement error on the stride length regime are comparable for all three stride definitions. The low stride length regime is on average overestimated while high stride lengths are on average underestimated. This reflects an imbalance found in the training data. Since a fixed track length was used for data acquisition, there are more examples from the low stride length regime on the dataset compared to high stride length strides.

\section{Discussion}
Overall, our results show that stride length estimation with deep convolutional neural networks based on stride-specific inertial sensor data is possible. This is the first major contributions of this work. We show feasibility and contribute a completely new type of algorithm to the field of stride length estimation. Although the achieved precision of 5.4 cm is not sufficient to track the subtle yearly decline in stride length in healthy individuals, it is sufficient to resolve stride length differences observed in PD patients \cite{Hass2012}.

The best result achieved with the method described here outperforms the approach by Rampp et al. \cite{Rampp2015} by $3.0$ cm (36\%) in precision and also gives a more balanced error distribution. We can directly compare these two approaches since they are evaluated on an identical dataset. This is not the case for the other methods found in literature (see Tab. \ref{tab:relatedwork}). Although our approach shows similar performance, the discrepancy in evaluation datasets and their characteristics limits comparability. This manifests the urgent need for a standardised and clinically relevant benchmark dataset regarding stride length estimation from inertial sensor data. Given the number of subjects and the presence of gait alterations in the study population used here, it would serve well as such a benchmark dataset. Furthermore, it is publicly available at \url{https://www5.cs.fau.de/activitynet/benchmark-datasets/digital-biobank/}.

The second main contribution of this work is the independence of stride definition. The mean accuracy and precision of the approach described is only marginally affected by the stride definition (Tab. \ref{tab:results}). This is particularly interesting since the different events in the datastream forming the stride definition can be determined with different degrees of exactness. The MS event, for example, is defined as the instant of lowest energy in all gyroscope axes during the stride (see \cite{Rampp2015} for details). In cases where the high frequency vibrations caused during heel-strike are not dampened enough to ensure a relatively stable sensor during the stance phase, this can lead to a temporal offset regarding individual MS events. The heel-strike detection, however, is defined as a minimum in the accelerometer x-axis (see \cite{Rampp2015} for details) and can thus be detected more accurately. Subsequently, the MS$\to$MS stride definition is not as exact as for example the HS$\to$HS stride definition. This results in more extreme outliers for the MS$\to$MS strides in the Bland-Altman plots (Fig. \ref{fig:bland-altman-plots}). The fact that this stride definition still yields best mean accuracy and precision is probably coincidence here.
\FloatBarrier

From a clinical perspective, the independence of stride definition opens up new possibilities regarding analysis of impaired gait. In clinical practice, the zero-velocity assumption is easily violated as in the case of spastic forms of gait alterations as e.g. experienced by post-stroke, multiple sclerosis or spastic paraplegia patients. Standard double-integration methods would fail in estimating spatial stride parameters in these clinically relevant situations. This is not the case with the approach presented here.

Lastly, the high precision on the training dataset shown in Fig. \ref{fig:training} promises excellent results in case of individualization. If we have a training/calibration dataset to tune the neural network to the patient at hand, i.e. by collecting one gait sequence instrumented with inertial sensors and the GAITRite, the proposed method could reach the precision of the reference. This would enable highly accurate stride-by-stride analysis of spatial gait parameters in daily clinical routine or outpatient monitoring at a level sufficient to address clinically relevant phenomena as for example the decline in stride length with disease progression in PD patients or medication effects. Certainly, this has to be formally evaluated in the future and is a possible extension of the presented method.

Furthermore, our methodology is not limited to stride length estimation. Targeting the underlying regression task at a different stride-specific gait characteristic e.g. stride width or timings of individual gait phases does not involve changing the underlying algorithms. The very same technology can be used to estimate different target parameters by simply changing or extending the reference annotation on the training dataset.

The main limitation of the proposed method is that the training of a regression function between stride-specific sensor data and the spatial parameter stride length implicitly depends on the training set. Therefore, the dataset used for training has to capture as much variability as possible for the problem at hand. Otherwise, application on a population that differs from the one used for training might result in lack of model-validity. This is the profound difference between the data-driven approach described here and a double-integration approach that is based on geometric and physical reasoning and does not encounter this kind of problem. However, given enough training data this problem can be solved. Additionally, preliminary experiments already show good generalisation to a PD population without any adaptation of the model trained on the geriatric population presented here. A second profound difference to double-integration approaches is the computational complexity and real-time capability. Although the trained model consists of convolutions and matrix-multiplications with a fixed set of kernels/weights, double-integration approaches will certainly be less complex. Depending on the (medical) application, however, real-time capability might not even be a necessity.

Future work includes a rigorous exploration of the parameter space (number of layers, kernels, weights, etc.) and introduction of a dynamic stop criterion for the optimization that depends on the current performance. Early stopping could also be investigated in this respect.  Additionally, the imbalance on the training data mentioned earlier needs to be addressed. Since limiting the number of strides per patient by random subsampling dramatically reduces the amount of training samples in the current dataset, a balanced training set could for example be acquired by means of simulation. Here, a biomechanical model could be used to sample the space of possible, three-dimensional foot trajectories and output the accelerometer and gyroscope readings for each stride as well as the corresponding target parameter of interest, e.g. stride length. This would also enable much richer datasets for training the network and thereby resolve the aforementioned limitation. Regarding the individualization of the proposed model to a specific patient, this work only shows promising substitute results. This has to be investigated in further detail. Finally, the application of deep convolutional neural networks to other spatial stride parameters, e.g. stride width, or the extraction of complete foot trajectories including foot orientation from stride-specific inertial sensor data will be investigated in future work.

In summary, we used convolutional neural networks to estimate stride length based on stride-specific inertial sensor data captured at the subject's feet. We provide technical validation of the proposed method on a publicly available and clinically relevant benchmark dataset. Furthermore, the proposed approach is robust w.r.t. the stride definition and thus not subject to the zero-velocity phase assumption. Consequently, the current work enables estimation of stride length in clinically relevant situations where the diversity of gait alterations easily violates the assumptions made by state-of-the-art methods.

\FloatBarrier
\section*{Acknowledgements}
This work was supported by the FAU Emerging Fields Initiative (EFIMoves). The authors would like to thank all participants of the study for their contributions.
\bibliographystyle{IEEEtran}