text:
\section{Introduction}
In recent years there has been an increasing number of fields where the underlying physical systems are represented by data. Increasing memory, processor speed, and the advent of distributed computing architectures have enabled us to handle and analyze data with more precision and to address learning problems for nonlinear complex systems at an unprecedented scale. Another advantage of data-driven analysis lies in the fact that for many naturally occurring complex systems and engineered systems with emergent phenomena, e.g. biological systems, inter-dependent critical infrastructure, social networks, financial systems, it may not always be possible to derive and analyze theoretical mathematical models of the underlying systems \cite{yeung2015global}. In such cases, one has to resort to data-driven techniques for understanding the behavior of such systems.

In the realm of dynamical systems, there is an increasing interest in transfer operator theoretic techniques for analysis and control \cite{yeung2018koopman,yeung2017learning, Dellnitz_Junge,Mezic2000,froyland_extracting,Junge_Osinga,Mezic_comparison,
Dellnitztransport,mezic2005spectral,Mehta_comparsion_cdc,Vaidya_TAC,
raghunathan2014optimal,susuki2011nonlinear,mezic_koopmanism,
mezic_koopman_stability,surana_observer}. The operator based techniques differ from the classical study of dynamical systems in the sense that instead of studying the dynamical system on the configuration manifold and the tangent and cotangent bundles associated with them, the evolution is studied on the space of functions \cite{johnson2018class} or on the space of measures. The big advantage of using the function space or the measure space instead of the configuration space is the fact that in those spaces, the evolution is linear, even if the actual system is nonlinear. Moreover, the operators, namely, Perron-Frobenius operator and the Koopman operator are positive Markov operators. These properties can be exploited to have probabilistic interpretations and can be used for various applications like the optimal placement of sensors and actuators \cite{optimal_placement_ECC, optimal_placement_JMAA} etc. However, the trade-off lies in the fact that the operators defined on the function space or the measure space are infinite dimensional.

Another advantage of the operator theoretic framework for analysis and control of dynamical systems is that the methods developed are data-centric and in many cases wholly data-driven. In particular, a finite-dimensional approximation of both Perron-Frobenius and Koopman operators can be constructed from time-series data obtained from experiments \cite{yeung2017learning}. Towards this goal various data-driven methods are proposed for the finite dimensional approximation of these operators \cite{dellnitz2002set, Mezic2000,DMD_schmitt,rowley2009spectral,EDMD_williams}, with Dynamic Mode Decomposition (DMD) and extended DMD being the ones which are used extensively. Recent works have also addressed the problem of computing these operators for systems with process and observation noise and for Random Dynamical Systems (RDS)  \cite{mezic_stochastic_koopman_spectrum,PhysRevE.96.033310, robust_DMD_ACC,robust_DMD_arxiv}. In \cite{mezic_stochastic_koopman_spectrum} the authors have provided a characterization of the spectrum and eigenfunctions of the Koopman operator for discrete and continuous time RDS, while in \cite{PhysRevE.96.033310}, the authors have provided an algorithm to compute the Koopman operator for systems with both process and observation noise. In \cite{robust_DMD_ACC,robust_DMD_arxiv} the authors used robust optimization-based techniques to compute the approximate Koopman operator for data sets of finite length and have shown that normal DMD or EDMD and subspace DMD \cite{PhysRevE.96.033310} lead to an unsatisfactory approximation of Koopman operator for data sets of finite length.

A different challenge that researchers often have to account for is the scenario when the data set is not only finite, but also small, that is sparse. Sparse data refers to
data sets with few data points and can have an immense effect on the ability to train the
Koopman operator into producing accurate predictions. In particular, in the case of sparse data the existing DMD and EDMD algorithms may lead to an ill-conditioned least-square problem.
In this paper, we address the problem of computation of Koopman operator when the data set is sparse. We append artificial data points to the sparse data set to enrich the data and use robust optimization-based techniques to obtain the approximate Koopman operator. The robust optimization problem is a min-max problem which can be approximated as a least squares problem with a regularization term. The regularization parameter imposes sparsity in the Koopman operator. Moreover, it prevents over-fitting of the data and hence can be used to design a data-driven predictor \cite{korda_mezic_predictor}.

The organization of the paper is as follows. In section \ref{section_transfer_operators} we provide the basics of transfer operators followed by a discussion of DMD and EDMD algorithms in section \ref{section_DMD}. In section \ref{section_sparse_Koopman} we present the main results of the paper and state the algorithm to construct the Koopman operator for sparse data. Design of the robust predictor is discussed in section \ref{section_predictor}, with simulation results in section \ref{section_simulation}. Finally we conclude the paper in \ref{section_conclusion}.

\section{Transfer Operators for Dynamical Systems}\label{section_transfer_operators}
Consider a discrete-time dynamical system
\begin{eqnarray}\label{system}
z_{t+1} = T(z_t)
\end{eqnarray}
where $T:Z\subset \mathbb{R}^N \to Z$ is assumed to be an invertible smooth diffeomorphism.  Associated with the dynamical system (\ref{system}) is the Borel-$\sigma$ algebra ${\cal B}(Z)$ on $Z$ and the vector space ${\cal M}(Z)$ of bounded complex valued measures on $X$. With this, two linear operators, namely, Perron-Frobenius (P-F) and Koopman operator, can be defined as follows \cite{Lasota} :
\begin{definition}[Perron-Frobenius Operator]
The Perrorn-Frobenius operator $\mathbb{P}:{\cal M}(Z)\to {\cal M}(Z)$ is given by
\[[\mathbb{P}\mu](A)=\int_{{\cal Z} }\delta_{T(z)}(A)d\mu(z)=\mu(T^{-1}(A))\]
$\delta_{T(z)}(A)$ is stochastic transition function which measure the probability that point $z$ will reach the set $A$ in one time step under the system mapping $T$.
\end{definition}

\begin{definition}[Invariant measures] Invariant measures are the fixed points of
the P-F operator $\mathbb{P}$ that are also probability measures. Let $\bar \mu$ be the invariant measure then, $\bar \mu$ satisfies
\[\mathbb{P}\bar \mu=\bar \mu.\]
\end{definition}

If the state space $Z$ is compact, it is known that the P-F operator admits at least one invariant measure.

\begin{definition} [Koopman Operator]
Given any $h\in\cal{F}$, $\mathbb{U}:{\cal F}\to {\cal F}$ is defined by
\[[\mathbb{U} h](z)=h(T(z))\]
where $\cal F$ is the space of function (observables) invariant under the action of the Koopman operator.
\end{definition}

Both the P-F operator and the Koopman operator are linear operators, even if the underlying system is non-linear. But while analysis is made tractable by linearity, the trade-off is that these operators are typically infinite dimensional. In particular, the P-F operator and Koopman operator often will lift a dynamical system from a finite-dimensional space to generate an infinite dimensional linear system in infinite dimensions.
\begin{properties}\label{property}
Following properties for the Koopman and Perron-Frobenius operators can be stated \cite{Lasota}.

\begin{enumerate}
\item [a).] For the Hilbert space ${\cal F}=L_2(Z,{\cal B}, \bar \mu)$
\begin{eqnarray*}
&&\parallel \mathbb{U}h\parallel^2=\int_Z |h(T(z))|^2d\bar \mu(z)
\nonumber\\&=&\int_Z | h(z)|^2 d\bar\mu(z)=\parallel h\parallel^2
\end{eqnarray*}
where $\bar \mu$ is an invariant measure. This implies that Koopman operator is unitary.

\item [b).] For any $h\geq 0$, $[\mathbb{U}h](z)\geq 0$ and hence Koopman is a positive operator.

\item [c).]For invertible system $T$, the P-F operator for the inverse system $T^{-1}:Z\to Z$ is given by $\mathbb{P}^*$ and $\mathbb{P}^*\mathbb{P}=\mathbb{P}\mathbb{P}^*=I$. Hence, the P-F operator is unitary.

\item [d).] If the P-F operator is defined to act on the space of densities i.e., $L_1(Z)$ and Koopman operator on space of $L_\infty(Z)$ functions, then it can be shown that the P-F and Koopman operators are dual to each other \footnote{with some abuse of notation we use the same notation for the P-F operator defined on the space of measure and densities.}
\begin{eqnarray*}
&&\left<\mathbb{U} f,g\right>=\int_Z [\mathbb{U} f](z)g(z)dx\nonumber\\&=&\int_Xf(y)g(T^{-1}(y))\left|\frac{dT^{-1}}{dy}\right|dy=\left<f,\mathbb{P} g\right>
\end{eqnarray*}
where $f\in L_{\infty}(Z)$ and $g\in L_1(Z)$ and the P-F operator on the space of densities $L_1(Z)$ is defined as follows
\[[\mathbb{P}g](z)=g(T^{-1}(z))|\frac{dT^{-1}(z)}{dz}|.\]

\item [e).] For $g(z)\geq 0$, $[\mathbb{P}g](z)\geq 0$.

\item [f).] Let $(Z,{\cal B},\mu)$ be the measure space where $\mu$ is a positive but not necessarily the invariant measure of $T:Z\to Z$, then the P-F operator $\mathbb{P}:L_1(Z,{\cal B},\mu)\to L_1(Z,{\cal B},\mu)$  satisfies  following property:

 \[\int_Z [\mathbb{P}g](z)d\mu(z)=\int_Z g(z)d\mu(x).\]\label{Markov_property}
\end{enumerate}
\end{properties}

\section{Dynamic Mode Decomposition (DMD) and Extended Dynamic Mode Decomposition (EDMD)}\label{section_DMD}
Dynamic Mode Decomposition was first proposed in \cite{DMD_schmitt} in the context of analysis of fluid flow analysis. DMD algorithm is a way of approximating the spectrum of the Koopman operator. DMD algorithm was improved upon (EDMD) in \cite{EDMD_williams}, where the authors used a choice of observables, called a \textit{dictionary}. In this improved algorithm, the Koopman operator is approximated as a linear map on the span of the finite set of dictionary functions. In this section, we briefly describe the EDMD algorithm for approximating the Koopman operator.

Consider snapshots of data set obtained from simulating a discrete time dynamical system $z\mapsto T(z)$ or from an experiment
\begin{eqnarray}
X_p = [x_1,x_2,\ldots,x_M],& X_f = [y_1,y_2,\ldots,y_M] \label{data}
\end{eqnarray}
where $x_i\in X$ and $y_i\in X$. The two pair of data sets are assumed to be two consecutive snapshots i.e., $y_i=T(x_i)$. Let $\mathcal{D}=
\{\psi_1,\psi_2,\ldots,\psi_K\}$ be the set of dictionary functions or observables, where $\psi : X \to \mathbb{C}$. Let ${\cal G}_{\cal D}$ denote the span of ${\cal D}$ such that ${\cal G}_{\cal D}\subset {\cal G}$, where ${\cal G} = L_2(X,{\cal B},\mu)$. The choice of dictionary functions are very crucial and it should be rich enough to approximate the leading eigenfunctions of Koopman operator. Define vector valued function $\mathbf{\Psi}:X\to \mathbb{C}^{K}$
\begin{equation}
\mathbf{\Psi}(\boldsymbol{x}):=\begin{bmatrix}\psi_1(x) & \psi_2(x) & \cdots & \psi_K(x)\end{bmatrix}
\end{equation}
In this application, $\mathbf{\Psi}$ is the mapping from physical space to feature space. Any function $\phi,\hat{\phi}\in \mathcal{G}_{\cal D}$ can be written as
\begin{eqnarray}
\phi = \sum_{k=1}^K a_k\psi_k=\boldsymbol{\Psi^T a},\quad \hat{\phi} = \sum_{k=1}^K \hat{a}_k\psi_k=\boldsymbol{\Psi^T \hat{a}}
\end{eqnarray}
for some set of coefficients $\boldsymbol{a},\boldsymbol{\hat{a}}\in \mathbb{C}^K$. Let \[ \hat{\phi}(x)=[\mathbb{U}\phi](x)+r,\]
where $r$ is a residual function that appears because $\mathcal{G}_{\cal D}$ is not necessarily invariant to the action of the Koopman operator. To find the optimal mapping which can minimize this residual, let $\bf K$ be the finite dimensional approximation of the Koopman operator. Then the  matrix $\bf K$ is obtained as a solution of least square problem as follows
\begin{equation}\label{edmd_op}
\min\limits_{\bf K}\parallel {\bf G}{\bf K}-{\bf A}\parallel_F
\end{equation}
\begin{eqnarray}\label{edmd1}
\begin{aligned}
& {\bf G}=\frac{1}{M}\sum_{m=1}^M \boldsymbol{\Psi}({x}_m)^\top \boldsymbol{\Psi}({x}_m)\\
& {\bf A}=\frac{1}{M}\sum_{m=1}^M \boldsymbol{\Psi}({x}_m)^\top \boldsymbol{\Psi}({y}_m),
\end{aligned}
\end{eqnarray}
with ${\bf K},{\bf G},{\bf A}\in\mathbb{C}^{K\times K}$. The optimization problem (\ref{edmd_op}) can be solved explicitly to obtain following solution for the matrix $\bf K$
\begin{eqnarray}
{\bf K}_{EDMD}={\bf G}^\dagger {\bf A}\label{EDMD_formula}
\end{eqnarray}
where ${\bf G}^{\dagger}$ is the psedoinverse of matrix $\bf G$.
DMD is a special case of EDMD algorithm with ${\bf \Psi}(x) = x$.

\section{Koopman Operator Construction for Sparse Data}\label{section_sparse_Koopman}
In many experiments, it often is the case that the data obtained is not rich enough to achieve a good enough approximation of Koopman operator. In this case, existing algorithms like DMD or EDMD fail to generate acceptable Koopman operators. In fact, in many instances, these algorithms lead to unstable eigenvalues, even though the underlying system is stable \cite{robust_DMD_ACC,robust_DMD_arxiv}. In this section, a robust optimization-based framework is presented for a better approximation of the Koopman eigenspectrum for sparse data.

\subsection{Enrichment of the existing dataset}
Let $\bar{X}_p = [{x}_1,{x}_2,\cdots , {x}_M]$ and $\bar{X}_f = [{y}_1,{y}_2,\cdots , {y}_M]$ be the snap-shots of data obtained from a discrete-time dynamical system $z\to T(z)$, where $z\in Z \subset \mathbb{R}^N$, $x_i\in Z$ and $y_i\in Z$. The two pair of data sets are assumed to be two consecutive snapshots i.e., $y_i=T(x_i)$. Consider ${x}_i+\delta x_i$, where $\parallel \delta x_i\parallel \leq \lambda_X$. Since, $T$ is a diffeoemorphism,
\begin{eqnarray}
T(x_i+\delta x_i)\approx T(x_i)+\frac{\partial T}{\partial x}\delta x_i = y_i + \delta y_i.
\end{eqnarray}
Since $\parallel \delta x_i\parallel \leq \lambda_X$, $x_i+\delta x_i\in B(x_i,\lambda_X)$, where
\begin{eqnarray*}
B(x_0,r) = \{x\in\mathbb{R}^N| \parallel x - x_0 \parallel \leq r\}.
\end{eqnarray*}
Again $T$ being a smooth diffeomorphism implies that $\frac{\partial T}{\partial x} \leq \lambda_T$ and hence $T(x_i+\delta x_i)\in B(y_i,\lambda_Y)$, where $\lambda_Y\leq \lambda_X\lambda_T$.

Hence, a \textit{sufficiently small perturbation} to any point $x_i$ will result in a small enough perturbation of $y_i$ and it is the boundedness of the small perturbations that is used to enrich the existing data set. In particular, to each observed data tuple $(x_i,y_i)$, we augment an extra data point $(x_i+\delta x_i, y_i+\delta y_i)$. Note that, for each data point more than one data point can be augmented, but for clarity, we will discuss the situation where only one extra artificial data point is augmented to each observed data point. Hence, an artificial data set
\begin{eqnarray}\label{enriched_dataset}
\begin{aligned}
{X}_p &= [{x}_1,\cdots , {x}_M,x_1+\delta x_1, \cdots , x_M + \delta x_M]\\
& = [{x}_1,\cdots ,{x}_{2M}]\\
{X}_f &= [{y}_1,\cdots , {y}_M,y_1+\delta y_1, \cdots , y_M + \delta y_M]\\
& = [{y}_1,\cdots ,{y}_{2M}]
\end{aligned}
\end{eqnarray}
is created with $2(M+1)$ data points. Here $x_{M+i} = x_i + \delta x_i$ and $y_{M+i} = y_i + \delta y_i$.

\subsection{Robust Optimization Formulation}

Let $X_p=[{x}_1,\cdots ,{x}_{2M}]$ and $X_f=[y_1,\cdots y_{2M}]$ be the enriched data set, as given in (\ref{enriched_dataset}). However, all the data points $(x_i,y_i)$ are not obtained from the system but has some error due to the artificial data points. The errors in the artificial data points can be thought of as measurement noise and the uncertainty acts as an adversary which tries to maximize the residual. Hence to obtain the Koopman operator $\bf K$ for uncertain data, a robust optimization problem can be formulated as the following $\min-\max$ optimization problem.

\begin{equation}\label{edmd_robust}
\min\limits_{\bf K}\max_{\delta\in \Delta}\parallel {\bf G}_\delta{\bf K}-{\bf A}_\delta\parallel_F=:\min\limits_{\bf K}\max_{\delta\in \Delta} {\cal J}({\bf K}, {\bf G}_\delta,{\bf A}_\delta)
\end{equation}
where
\begin{eqnarray}
&&{\bf G}_\delta=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({x}_i)\nonumber\\
&&{\bf A}_\delta=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({ y}_{i}),
\end{eqnarray}
with ${\bf K},{\bf G}_\delta,{\bf A}_\delta\in\mathbb{C}^{K\times K}$.

The robust optimization problem (\ref{edmd_robust}), is in general non-convex because the cost $\cal J$ may not be a convex function of $\delta$.

\begin{proposition}
The optimization problem (\ref{edmd_robust}) can be approximated as
\begin{equation}\label{edmd_robust_convex}
\min\limits_{\bf K}\max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{equation}
where $\cal U$ is a compact set in $\mathbb{R}^{K\times K}$.
\end{proposition}
\begin{proof}
From Taylor series expansion we have, ${\bf \Psi}(x_i+\delta x_i) = {\bf \Psi}(x_i) +  {\bf \Psi}'(x_i) \delta x_i+h.o.t.$, where ${\bf \Psi}'(x_i)$ is the first derivative of ${\bf \Psi}(x)$ at $x_i$. Hence,
\begin{eqnarray*}
{\bf G}_{\delta} &\approx& {\bf G} + \frac{1}{2M}\sum_{i=1}^{2M}{\bf \Psi}^\top (x_i)\delta x_i {\bf \Psi}'(x_i)\\
&=& {\bf G} +\delta {\bf G}
\end{eqnarray*}
where $\delta {\bf G} = \frac{1}{2M}\sum_{i=1}^{2M}{\bf \Psi}^\top (x_i)\delta x_i {\bf \Psi}'(x_i)$.

Moreover,
\begin{eqnarray*}
&&\parallel \delta {\bf G} \parallel_F = \parallel \frac{1}{2M}\sum_{i=1}^{2M}{\bf \Psi}^\top (x_i)\delta x_i {\bf \Psi}'(x_i) \parallel_F\\
&\leq& \frac{1}{2M}\sum_{i = 1}^{2M}\parallel {\bf \Psi}^\top (x_i)\delta x_i {\bf \Psi}'(x_i)\parallel_F \\
&\leq&\frac{1}{2M}\sum_{i=1}^{2M}\parallel {\bf \Psi}^\top (x_i)\parallel_F \cdot \parallel\delta x_i \parallel_F \cdot \parallel {\bf \Psi}'(x_i)\parallel_F
\end{eqnarray*}
Hence, $\delta {\bf G}$ belongs to a compact set ${\cal U}_1$. Similarly, one can show ${\bf A}_\delta \approx {\bf A} + \delta {\bf A}$ and $\delta {\bf A}$ belongs to a compact set ${\cal U}_2$. Letting ${\cal U}={\cal U}_1\cup{\cal U}_2$, proves the proposition.
\end{proof}
The above proposition allows us to compute the Koopman operator $\bf K$ as a solution of a robust optimization problem (\ref{edmd_robust_convex}). The optimization problem (\ref{edmd_robust_convex}) has interesting connections with optimization problem involving regularization. In particular, one has the following theorem.

\begin{theorem}
The optimization problem
\begin{equation}
\min\limits_{\bf K}\max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{equation}
is equivalent to the following optimization problem
\begin{eqnarray}\label{rob_eqv}
\min\limits_{\bf K}\parallel {\bf G}{\bf K}-{\bf A}\parallel_F+\lambda \parallel {\bf K}\parallel_F\label{regular}
\end{eqnarray}
\end{theorem}
\begin{proof}
The min-max optimization problem is
\begin{eqnarray}
{\cal J} = \min_{\bf K}\max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{eqnarray}
Fix ${\bf K}\in\mathbb{R}^{K\times K}$ and let
\begin{eqnarray*}
r = \max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{eqnarray*}
be the worst-case residual.
Let $\max\{\parallel u \parallel_F | u \in {\cal U}\}\leq \lambda$. Then,
\begin{eqnarray}\nonumber\label{opt_reg_leq}
&&r = \max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}} \parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F\\ \nonumber
&& \leq \max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}} \parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \parallel \delta {\bf G}{\bf K} - \delta {\bf A} \parallel_F\\
&& \leq \parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F + \lambda \parallel{\bf I}\parallel_F \\
&&\leq \parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F + \lambda K\end{eqnarray}
Again, let $\begin{pmatrix}
\delta {\bf G}^\star &, & \delta {\bf A}_i^\star
\end{pmatrix} = \frac{\lambda u}{\sqrt \parallel {\bf K}_i\parallel_2 + 1}\begin{pmatrix}
{\bf K}_i^\top &,& & -1
\end{pmatrix}$
where ${\bf A}_i$ and ${\bf K}_i$ are the $i^{th}$ columns of $\bf A$ and $\bf K$ respectively and
\begin{eqnarray*}
u = \begin{cases}
     \frac{{\bf G}{\bf K}_i-{\bf A}_i}{\parallel {\bf G}{\bf K}_i-{\bf A}_i \parallel_2},\textnormal{ if } {\bf G}{\bf K}_i\neq{\bf A}_i\\
     \textnormal{any unit norm vector otherwise.}
    \end{cases}
\end{eqnarray*}
Let
\begin{eqnarray*}
r_i = \max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}_i-({\bf A}_i+\delta {\bf A}_i)\parallel_F.
\end{eqnarray*}
Hence,
\begin{eqnarray}\nonumber
&&r_i\geq \parallel ({\bf G}{\bf K}_i-{\bf A}_i)+(\delta {\bf G}^\star {\bf K}_i -\delta {\bf A}_i^\star)\parallel_F\\ \nonumber
&& = \parallel ({\bf G}{\bf K}_i-{\bf A}_i) + \frac{\lambda({\bf G}{\bf K}_i-{\bf A}_i)}{\parallel {\bf G}{\bf K}_i-{\bf A}_i \parallel_2}({\bf K}_i^\top {\bf K}_i+1)\parallel_F\\
&&= \parallel {\bf G}{\bf K}_i-{\bf A}_i\parallel_F + \lambda \parallel {\bf K}_i\parallel_F + \lambda
\end{eqnarray}
Hence,
\begin{eqnarray}\label{opt_reg_geq}
r = \sum_{i=1}^K r_i \geq \parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F + \lambda K.
\end{eqnarray}

Finally, from (\ref{opt_reg_leq}) and (\ref{opt_reg_geq}) and minimizing over $\bf K$ we get
\begin{eqnarray*}
\min\limits_{\bf K}\max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{eqnarray*}
is equivalent to
\begin{eqnarray*}
\min\limits_{\bf K}\parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F + \lambda K.
\end{eqnarray*}
Moreover, since $\lambda$ and $K$ are positive constants,
\begin{eqnarray*}
\min\limits_{\bf K}\max_{\delta{\bf G},\delta{\bf A}\in  {\cal U}}\parallel ({\bf G}+\delta {\bf G}){\bf K}-({\bf A}+\delta {\bf A})\parallel_F
\end{eqnarray*}
is equivalent to
\begin{eqnarray*}
\min\limits_{\bf K}\parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F.
\end{eqnarray*}
\end{proof}

The above theorem allows computation of the approximate Koopman operator as a solution of an optimization problem with a regularization term. In particular, the approximate Koopman operator can be obtained as a solution of the following optimization problem
\begin{eqnarray}\label{robust_opt}
\parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F.
\end{eqnarray}
where
\begin{eqnarray}
\begin{aligned}
&{\bf G}=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({x}_i)\\
&{\bf A}=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({ y}_{i}).
\end{aligned}
\end{eqnarray}
\begin{algorithm}
\caption{Computation of Koopman Operator from Limited Data}
\begin{enumerate}
\item{To the existing data set $D = [x_1,\cdots, x_M]$ add new data points $\tilde{x}_i = x_i+\delta x_i$, where $\parallel \delta x_i\parallel \leq c$.}
\item{Form the enriched data set $\bar{D} = [x_1,\cdots, x_M, \tilde{x}_1,\cdots, \tilde{x}_M]$.}
\item{Form the sets $X_p = [x_1,\cdots, x_{M-1}, \tilde{x}_1,\cdots, \tilde{x}_{M-1}]$ and $X_f = [x_2,\cdots, x_M, \tilde{x}_2,\cdots, \tilde{x}_M]$.}
\item{Fix the dictionary functions ${\bf \Psi} = [\psi_1,\cdots , \psi_K]$.}
\item{Solve the optimization problem to obtain the approximate Koopman operator $\bf K$
\begin{eqnarray*}
\parallel {\bf G} {\bf K} - {\bf A} \parallel_F + \lambda \parallel {\bf K}\parallel_F.
\end{eqnarray*}
where
\begin{eqnarray*}
&&{\bf G}=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({x}_i)\nonumber\\
&&{\bf A}=\frac{1}{2M}\sum_{i=1}^{2M} \boldsymbol{\Psi}({ x}_i)^\top \boldsymbol{\Psi}({ y}_{i}).
\end{eqnarray*}
}
\end{enumerate}
\label{algo}
\end{algorithm}

\begin{remark}
Algorithm (\ref{algo}) describes the procedure of adding $M$ extra data points, but it should be noted that one can add more or less number of artificial data points.
\end{remark}

\section{Design of Robust Predictor}\label{section_predictor}
The Koopman operator generates a linear system in a higher dimensional space, even if the underlying system is linear. The linearity of the operator enables the design of linear predictors for nonlinear systems. The following is presented briefly for the self-containment of the paper and for details the readers are referred to \cite{korda_mezic_predictor}. Let $\{x_0,\ldots,x_M\}$ be the training data-set and $\bf K$  be the finite-dimensional approximation of the transfer Koopman operator obtained using algorithm ${\bf 1}$. Let $\bar x_0$ be the initial condition from which the future is to be predicted. The initial condition from state space is mapped to the feature space using the same choice of basis function used in the robust approximation of Koopman operator i.e., \[\bar x_0\implies {\bf \Psi}(\bar x_0)^\top=: {\bf z}\in \mathbb{R}^K.\] This initial condition is propagated using Koopman operator as \[{\bf z}_n={\bf K}^n{\bf z}.\]
The predicted trajectory in the state space is then obtained as
\[\bar x_n=C {\bf z}_n\]
where matrix $C$ is obtained as the solution of the following least squares problem
\begin{eqnarray}\label{C_pred}
\min_C\sum_{i = 1}^M \parallel x_i - C \boldsymbol \Psi (x_i)\parallel_2^2
\end{eqnarray}

\section{Simulations}\label{section_simulation}

In this section, we demonstrate the efficiency of the proposed algorithm on three different dynamical systems. In particular, we construct the Koopman operator for a linear system, a non-linear system and a system governed by a Partial Differential Equation (PDE).
\subsection{Network of Coupled Oscillators}
Consider a network of coupled linear oscillators given by
\begin{eqnarray}\label{coup_osc}
\ddot{\theta}_k &=& -\mathcal{L}_k\theta - d\dot{\theta}_k, \quad k = 1,\cdots , N
\end{eqnarray}
where $\theta_k$ is the angular position of the $k^{th}$ oscillator, $N$ is the number of oscillators, $\mathcal{L}_k$ is the $k^{th}$ row of the Laplacian $\mathcal{L}$ and $d$ is the damping coefficient. The Laplacian $\cal L$ is chosen such that the network is a ring network with 20 oscillators (Fig. \ref{fig_oscillator}).

\begin{figure}[htp!]
\centering
\includegraphics[scale=.6]{oscillator}
\caption{Ring network of 20 linear oscillators.}\label{fig_oscillator}
\end{figure}

In these sets of simulations, the damping coefficient $d$ has been assumed the same for all the oscillators and is set equal to $0.4$. Data for all the states were collected for 100-time steps, with sampling time $\delta t = 0.01$ seconds and since the system is linear, linear basis functions were used for computation of the Koopman operator.
\begin{figure}[htp!]
\centering
\subfigure[]{\includegraphics[scale=.21]{20_osc_15_training_points_eig}}
\subfigure[]{\includegraphics[scale=.21]{20_osc_15_training_points_eig_zoomed}}
\caption{(a) Eigenvalue obtained using normal DMD on original training data and Robust DMD on enriched data set. (b) Dominant eigenvalues.}\label{eig_osc}
\end{figure}
The first 15-time steps data was used for training the Koopman operators. Normal DMD on the 15 data points yields positive eigenvalues with a significant real part, as shown in Fig. \ref{eig_osc}. For the Robust identification of Koopman operator, the original data set was enriched by adding 30 artificial data points and Robust DMD formulation (algorithm \ref{algo}) yields a much better approximation of the eigenvalues for the original system. The eigenvalues obtained using normal DMD and Robust DMD are shown in Fig. \ref{eig_osc}, where in Fig. \ref{eig_osc}a the complete spectrum is plotted and in Fig. \ref{eig_osc}b the dominant eigenvalues are shown.

\begin{figure}[htp!]
\centering
\subfigure[]{\includegraphics[scale=.215]{err_in_x3}}
\subfigure[]{\includegraphics[scale=.215]{err_in_x4}}
\caption{(a) Errors in the prediction of position of oscillator 3. (b) Errors in the prediction of position of oscillator 4.}\label{err_osc}
\end{figure}
As mentioned earlier, data were obtained for 100 times steps and the first 15 time steps were used for training the Koopman operator. Koopman operators thus obtained was used to predict the next 45 time steps and was used to compare the error. The errors in the prediction of the positions of oscillators 3 and 4, using both normal DMD and Robust DMD, are shown in Fig. \ref{err_osc}a and Fig. \ref{err_osc}b respectively. It can be observed that Robust DMD formulation generates much smaller error compared to normal DMD. In fact, this was expected, since Robust DMD with enriched data-set approximates the eigenspectrum much better compared to normal DMD.

\subsection{Stuart-Landau Equation}
The nonlinear Stuart-Landau equation on a complex function $z(t) = r(t)\exp (i\theta(t))$ is given by
\begin{eqnarray}\label{stu_lan}
\dot{z} = (\mu + i\gamma)z - (1+ i\beta)|z|^2z ,
\end{eqnarray}
where $i$ is the imaginary unit. The solution of (\ref{stu_lan}) evolves on the limit cycle $|z| = \sqrt{\mu}$. Hence, the continuous time eigenvalues lie on the imaginary axis. The discretized version of (\ref{stu_lan}) is

{\small
\begin{eqnarray}\label{stu_lan_dis}
\begin{pmatrix}
r_{t+1}\\
\theta_{t+1}
\end{pmatrix} = \begin{pmatrix}
r_t + (\mu r_t -r_t^3)\delta t\\
\theta_t + (\gamma - \beta r_t^2)\delta t
\end{pmatrix}
\end{eqnarray}
}

The set of dictionary functions were chosen as
\begin{eqnarray}\label{stuart_dic}
{\bf \Psi}(\theta_t) = \begin{pmatrix}
e^{-10i\theta_t} & e^{-9i\theta_t} & \cdots & e^{9i\theta_t} & e^{10i\theta_t}
\end{pmatrix}
\end{eqnarray}
 and data was collected for 150 time steps, with $\delta t = 0.01$ and initial condition $(1,\pi)$.

\begin{figure}[htp!]
\centering
\subfigure[]{\includegraphics[scale=.2]{eig_30_training_step}}
\subfigure[]{\includegraphics[scale=.2]{eig_30_training_step_zoomed}}
\caption{(a) Eigenvalues obtained with EDMD on original data and Robust EDMD on enriched data set. (b) Dominant eigenvalues.}\label{eig_nonlinear}
\end{figure}

The first 30-time steps data were used as the training data for training the Koopman operator. An extra 30 artificial points were added to the obtained data set to form the enriched data set and this enriched data set was used to compute the eigenspectrum of the Koopman operator using Robust EDMD algorithm. The eigenvalues obtained using the dictionary functions given in (\ref{stuart_dic}), with normal EDMD and Robust EDMD with enriched data set is shown in Fig. \ref{eig_nonlinear}a. Fig. \ref{eig_nonlinear}b shows the dominant eigenvalues and it can be observed that Robust EDMD provides a better approximation of the original eigenspectrum. In particular, normal EDMD generates unstable eigenvalues.

\begin{figure}[htp!]
\centering
\includegraphics[scale=.35]{err_in_r_150_step_prediction}
\caption{Comparison of errors in prediction of $r$ using Robust EDMD and normal EDMD. The top figure shows the prediction error using Robust EDMD and the lower plot shows prediction error using normal EDMD.}\label{prediction_r}
\end{figure}

\begin{figure}[htp!]
\centering
\includegraphics[scale=.35]{err_in_theta_150_step_prediction}
\caption{Comparison of errors in prediction of $\theta$ using Robust EDMD and normal EDMD. The top figure shows the prediction error using Robust EDMD and the lower plot shows prediction error using normal EDMD.}\label{prediction_theta}
\end{figure}

Further, using the Koopman operators obtained using both normal EDMD and Robust EDMD, future values of both $r$ and $\theta$ was predicted for the next 70 time steps. The errors in the prediction of $r$ and $\theta$ are shown in Fig. \ref{prediction_r} and Fig. \ref{prediction_theta} respectively. In all the error plots, the errors are plotted against the actual values of $r$ and $\theta$ and it can be observed that the errors in prediction for both $r$ and $\theta$ with Robust EDMD are significantly smaller than the prediction errors using normal EDMD.

\subsection{Burger-Equation}
The third example considered in this paper is the Burger equation. Burger equation is a successful but simplified partial differential equation which describes the motion of viscous compressible fluids. The equation is of the form
\[\partial_t u(x,t)+u\partial_x u=k\partial_x^2u\]
where $u$ is the speed of the gas, $k$ is the kinematic viscosity, $x$ is the spatial coordinate and $t$ is time.

In the simulation, choosing $k=0.01$, we approximated the PDE solution using the Finite Difference method \cite{KUTLUAY1999251} with the initial condition $u(x,0)=sin(2\pi x)$ and Dirichet boundary condition $u(0,t)=u(1,t)=0$. Given the spatial and temporal ranges, $x\in[0,1],\;t\in[0,1]$, the discretizaion steps are chosen as $\Delta t=0.02$ and $\Delta x=1\times10^{-2}$. With the above set of conditions, the flow $u$ is shown in Fig. \ref{burger_flow}.
\begin{figure}[htp!]
\centering
\includegraphics[scale=.275]{PDE_data1}
\caption{Flow field of Burger equation.}\label{burger_flow}
\end{figure}

Since the space discretization was chosen as $\Delta x = 1\times 10^{-2}$, there are 100 state variables. For computing the Koopman operator, 8-time steps data were used. 40 extra data points were added to enrich the data set and the Robust Koopman operator was computed using the enriched data set. Koopman operator using normal DMD was also computed for comparing the errors in prediction. The errors in the prediction of 35 future time steps for $x_{40}$ and $x_{100}$ is shown in Fig. \ref{burger_error}(a) and Fig.. \ref{burger_error}(b) respectively. It can be seen that the error in prediction using Robust Koopman operator from the enriched data set is much smaller as compared to the normal DMD.

\begin{figure}[htp!]
\centering
\subfigure[]{\includegraphics[scale=.205]{err_x40_8_time_step_training}}
\subfigure[]{\includegraphics[scale=.2]{err_x100_8_time_step_training}}
\caption{(a) Errors in prediction of $x_{40}$. (b) Errors in prediction of $x_{100}$.}\label{burger_error}
\end{figure}

We further used different training size data for computing the Koopman operator and compared the mean square error in prediction of all the states. In particular, we used both Robust DMD approach and normal DMD to predict 35-time steps from $t=100$, with 7 different training size data, namely 5, 10, 15, 20, 25, 30 and 35-time steps. For each of the training size data, we appended the data set with artificial data points so that there are 40 data points in total. The mean square errors in the prediction of the states are shown in Fig. \ref{burger_mse}.

\begin{figure}[htp!]
\centering
\subfigure[]{\includegraphics[scale=.45]{MSE_robust}}
\subfigure[]{\includegraphics[scale=.45]{MSE_DMD}}
\caption{(a) Mean square error of prediction of all the states using Robust DMD aproach. (b) Mean square error of prediction of all the states using normal DMD.}\label{burger_mse}
\end{figure}

Fig. \ref{burger_mse}(a) shows the mean square error in prediction using the proposed approach and Fig. \ref{burger_mse}(b) shows the mean square error using normal DMD. It can be clearly seen that errors using the proposed method are much smaller (of the order of $10^2$). Another observation is that normal DMD is not much sensitive to small variations in training data size, whereas the proposed method is more sensitive to training data size.

\section{Conclusions}\label{section_conclusion}

In this paper, we addressed the problem of computation of Koopman operator from sparse time series data. In certain experimental applications, it may not be possible to obtain time series data which is rich enough to approximate the Koopman operator. We propose an algorithm to compute the Koopman operator for such sparse data. The intuition was based on exploiting the differentiability of the system mapping to append artificial data points to the sparse data set and using robust optimization-based techniques to approximate the Koopman eigenspectrum. The efficiency of the proposed method was also demonstrated on three different dynamical systems and the results obtained were compared to existing Dynamic Mode Decomposition and Extended Dynamic Mode Decomposition algorithms to establish the advantage of our proposed algorithm and in the future; we hope to investigate the performance of our approach on real experimental data sets.

\bibliographystyle{IEEEtran}


text:
\section{INTRODUCTION}
\label{sec:INTRODUCTION}

\begin{figure}[t!hbp]
\centering
\includegraphics[width =\linewidth]{figures/intro2.pdf}
\caption{
    We transform semantically annotated point sets provided by the SemanticKITTI \cite{Behley2019} data set into a multi-layer grid map representation.
    Using these grid maps we train a pixel-wise semantic labeling with five non-semantic input layers and the semantic grid map as training labels.
}
\label{fig:IntroGridMapLayer}
\end{figure}

A successful deployment of intelligent, mobile systems such as automated vehicles in complex urban scenarios demands an excellent scene understanding and a reliable interpretation of the surroundings.
To reach this goal the car has to be equipped with a range of complementary sensor systems and a sufficient amount of computational resources to extract the relevant information.
Semantic segmentation – assigning a semantic class to each single pixel of an image or point in a point-cloud – is a crucial task in the process of achieving a meaningful virtual representation of a vehicle's environment.
In the last years, deep convolutional neural networks have proven to be very successful in solving multiple computer vision tasks and significant progress has been made in terms of semantic segmentation.
Here, most research has been conducted in the vision domain since the majority of available data sets for semantic segmentation are image-based \cite{BDD100} \cite{Cordts2016}.
Recently, the first large-scale laser-based semantic data set, called SemanticKITTI, was published. This will allow researchers to focus on the semantic classification of LiDAR measurements.
The usage of laser sensors has several advantages: They provide a 360\si{\degree} surround-view of precise depth information, are robust against varying lighting conditions and are able to measure the reflectance of objects. %

\IEEEpubidadjcol
However, the processing and classification of high-resolution 3D point clouds requires high computational effort and most deep learning architectures are optimized to operate on well-structured input representaitions.
Considering this situation, we propose to encode the laser measurements into a multi-layer grid map representation.
In contrast to an unstructured point cloud, the grid map representation is well suited for sensor fusion applications \cite{Nuss2015} \cite{Richter2019}.
This enables us to leverage high-quality deep learning algorithms from the image domain due to the structural properties of the grid maps.
In contrast to the semantic segmentation of LiDAR-based range maps \cite{Milioto2019}, our approach predicts a semantic grid, i.e. a dense environment representation with an equidistant cell structure.

Our approach differs from previous publications by generating the semantic grids with labeled LiDAR scans.
The main contributions of this paper are the following: Firstly, we show that semantic segmentation of grid cells can be learned using sparse input data from LiDAR sensors.
We compare the usage of multi-layer input of different richness and evaluate the results in terms of accuracy and inference time.
Since a main advantage of our approach is the dense prediction of equidistant grid cells, we also propose a method which exploits the semantic information of up to 300 surrounding laser scans to form a dense semantic ground truth and enables up to evaluate the dense predictions of our models properly.

\section{Related Work} \label{sec:related_work}

Conceptually, the presented work is influenced by multiple research fields which we will introduce in the following subsections.

\subsection{Grid Mapping} \label{sec:related_work_grid_mapping}
2D occupancy grid maps were first introduced by Moravec \cite{Moravec1989} and Elfes \cite{Elfes1989}.
Due to the orthographic projection used during mapping, measurements are scale invariant which makes grid maps well-suited for sensor fusion tasks \cite{Nuss2015}.
When grid maps are modeled as regular grids, efficient operations such as unary operations and convolutions can be implemented which enables their use on GPUs \cite{Homm2010} and with machine learning applications (e.g. \cite{Wirges2018Aug}).
Homm et al. \cite{Homm2010} show that range sensor measurements can be mapped into polar grid coordinates and then efficiently transformed to a cartesian grid map.
In addition to occupancy, multi-layer grid maps also encode other information such as average reflected energy, minimum and maximum height of measurements and information on occluded areas obtained via ray-casting.
These different features enable the application of unsupervised methods (e.g. Connected Components Labeling) and supervised learning methods (e.g. convolutional object detectors \cite{Wirges2018Obj}) to grid maps.

\subsection{Benchmarks and Data Sets} \label{sec:related_work_data_sets}
The PASCAL VOC challenges \cite{Everingham2014} were the first benchmarks for the pixel-wise semantic classification of camera images.
Although its pixel-wise annotated images cover many scenarios, the amount of available on-road and traffic images is small.
Therefore, the KITTI Vision benchmark suite \cite{Geiger2013} was the first data set with benchmarks for automated driving applications.
It was gradually extended with benchmarks for other applications e.g. for visual odometry, scene flow estimation \cite{Menze2015} and object detection.
In 2018, a benchmark for semantic segmentation in camera images \cite{AbuAlhaija2018} was added which provides an alternative to the Cityscapes data set \cite{Cordts2016}.
Recently, Behley et al. \cite{Behley2019} presented a benchmark for the point-wise semantic classification of KITTI range sensor measurements.
The released part of the data set comprises a subset of the KITTI odometry benchmark sequences in which every point in more than 23000 lidar scans was annotated with one of 28 semantic classes. Some of these classes are further divided in \textit{moving} and \textit{non-moving} objects.

\subsection{Semantic Segmentation} \label{sec:related_work_semantic_segmentation}
Semantic segmentation in camera images refers to the task of assigning a semantic label to each image pixel.
Due to availability of large, annotated data sets (e.g. \cite{Everingham2014, AbuAlhaija2018}), Deep Learning approaches quickly became the state-of-the-art in semantic segmentation.
Significant improvements were made by Long et al. \cite{Long2015}, who proposed a fully convolutional network  for pixel-wise semantic segmentation which resembles an encoder / decoder structure with skip connections to combine strong semantic features with spatially dense features.
This idea of feature processing at different scales was later refined by Lin et al. \cite{Lin2017}.
Their Feature Pyramid Network (FPN) adds features via 1x1 convolutions from the top-down path to features from the bottom-up path and is a common structure for many semantic segmentation and object detection models.
Inspired by FPN, Chen et al. developed DeepLab (e.g. \cite{Chen2018}), a modular development suite.
It uses atrous convolutions, batch normalization and can use different feature extractors such as Xception \cite{Chollet2017} or MobileNet \cite{Howard2019}.
In addition to depth-wise separable convolutions, which are also used in Xception, recent MobileNet variants use blocks of narrow-wide-narrow convolutions - termed inverted residuals - with linear residuals and a general structure optimized for bounded memory and computation.

Before semantic range measurement labels were available, semantic segmentation in grid maps was performed by projecting camera image labels onto a ground plane, e.g. via inverse perspective mapping \cite{Erkent2018, Lu2019}.
These label transfer methods however do not perform well due to imperfect sensor calibration, parallax and a violated plane assumption for larger distances.
Due to the sparse mapping of range measurement end points, some grid map layers are sparsely occupied.
Jaritz et al. \cite{Jaritz2018} propose ways to resolve issues with sparse observations such as removing batch norm in the first layer and using a validity mask.
\section{DATA SET GENERATION}
\label{sec:DATASETGENERATION}

\begin{figure}[b]
    \centering
    \includegraphics[width =\linewidth]{figures/z_coord.pdf}
    \caption{
        Longitudinal ray geometry of a multi-layer LiDAR sensor.
        Ray casting from the sensor origin to point detections (red) determines the potentially occupied space above ground (green), which is registered as a single z-value in layer \textit{Min. observed height} (blue).
    }
    \label{fig:z_coord}
    \end{figure}

\begin{figure*}[thbp]
\centering
\includegraphics[width =0.95\textwidth]{figures/new_pipeline2.pdf}
\caption{
Each pose (blue) is separately processed by the framework resulting in five single-shot input layer and one sparse ground truth (sparse GT). Additionally, the framework determine the surrounding poses (green) within a certain radius ($r$). By assuming known poses, the framework fuses the static semantic information of all surrounding poses into one common grid map. This dense, static GT is finally enriched by the \textit{moving} labels of the sparse GT.
}
\label{fig:poles}
\end{figure*}

\subsection{Grid Map Framework}
Each semantically labeled point cloud is mapped into a multi-layer grid map containing the layers depicted in Fig.~\ref{fig:IntroGridMapLayer}.
The following sparse layers are used:
\begin{itemize}
    \item \textbf{Intensity}, containing the average overall reflection intensities in a given grid cell.
    \item \textbf{Min. detected height}, height of the lowest point detection.
    \item \textbf{Max. detected height}, height of the highest point detection.
\end{itemize}

In addition to the sparse layers, which can be directly deduced from the point detections, rays are cast between the sensor origin and the point detections to obtain the following dense layers:
\begin{itemize}
    \item \textbf{Observability}, number of transmissions in a cell.
    \item \textbf{Min. observed height}, height of the lowest laser beam traversing the cell, see Fig.~\ref{fig:z_coord}.
\end{itemize}

In order to facilitate parallel computation and account for geometric sensor characteristics, all layers are first computed in polar coordinates and subsequently remapped into a cartesian coordinate system utilizing the texture mapping described in \cite{Homm2010}.

\subsubsection{Grid Resolution and Perceiving Distance}
A 2:1 grid map aspect ratio is used, oriented in such a way that the ego vehicle appears to be driving to the right.
This aspect ratio and orientation corresponds more accurately with the spatial distribution of detections in urban scenarios than with a quadratic grid map.
Each grid cell represents a region of 10cm$\times10$cm and each image has a resolution of $1001\times501$ pixels.
This corresponds to a region of interest of 100m$\times50$m with the sensor located in its center.

\subsubsection{Split of Data Set}
As suggested in the SemanticKITTI paper, we select the sequences 0-7 and 9-10 with a total number of 19130 scans for the training set and use the sequence 8 for model evaluation.
In contrast to the authors in \cite{Behley2019} we evaluate on the full evaluation set of 4071 scans to get a meaningful evaluation.
We argue, that we require a larger set, as our method reduces the amount of labeled points.

\subsubsection{Label Set Definition}

In \cite{Behley2019}, it is suggested to use 19 of 28 classes are for training. In the course of this work, the label set comprising 19 SemanticKITTI classes is further reduced in order to counter false classification of classes with similar appearance in a grid map and under-representation of rare classes.
The classes \textit{car, truck} and \textit{other-vehicle} are labeled as \textit{vehicle}, the classes \textit{traffic-sign}, \textit{fence} and \textit{pole} as \textit{object}, the classes \textit{bicyclist} and \textit{motorcyclist} as \textit{rider}, the classes \textit{bicycle} and \textit{motorcycle} as \textit{two-wheel} and the classes \textit{other-ground} and \textit{parking} as \textit{other-ground}.
This leaves 12 relevant classes and an additional class called \textit{unlabeled}.
All pixels labeled with the latter are excluded from the training process by ignoring them in the loss function.

\subsubsection{Encoding of Ground Truth}
Next, we consider the label histogram for each grid cell individually and determine the cell class $k_\textit{i}$ by calculating a weighted arg max:
\begin{align}
k_{\mathrm{i}} = \underset{k \in K}{\mathrm{argmax}} \: w_\mathrm{k} \:  n_{\mathrm{i,k}},
\end{align}
where $n_{\mathrm{i,k}}$ is the total number of points of class $\mathrm{k}$  assigned to grid cell $\mathrm{i}$ and $w_\mathrm{k}$ is the class-specific weighting factor.

We choose the class-specific weighting factors as
\begin{align}
\label{eq:matching_costs}
w_\mathrm{k} =
\begin{cases}
0 & \text{for \textit{unlabeled} }\\
1 & \text{for \textit{others} }\\
5 & \text{for \textit{vehicle}, \textit{person}, \textit{rider}, \textit{two-wheel}}\\
\end{cases}
\text{.}
\end{align}
The weighting factor for traffic participants is higher than for other classes.
This is a compromise between favoring traffic participants over other classes, if there is evidence for their presence and preventing distortion of their appearance.
Since the weighting factor of \textit{unlabeled} is zero, the label is never chosen for a grid cell if any point is assigned to it.
If no point is assigned to a grid cell, it is labeled as \textit{unlabeled}.

The resulting label distribution is shown in Table \ref{tab::statistics}. Even after class-dependent weighting, the data set shows a significant under-representation of vulnerable road users in the data set.

\begin{table}[!htbp]
    \centering
        \caption{
    Class distribution of semantic grid maps. A comparison of the label distributions for the training and evaluation set for both the sparse and dense ground truth.}
    \label{tab::statistics}
    \begin{tabular}{lcccc}
        \hline
	{}   &  \multicolumn{4}{c}{label distribution}\\
	{}   &  \multicolumn{2}{c}{sparse}&  \multicolumn{2}{c}{dense}\\
     	\cmidrule(lr){2-3} \cmidrule(lr){4-5}
        label & {$\DS{train}$} & {$\DS{val}$} & {$\DS{train}$}& {$\DS{val}$}\\
        \hline
vehicle		&	0.275	\% & 	0.353	\% & 	0.844	\% & 	1.070	\%\\
person		&	0.004	\% & 	0.011	\% & 	0.019	\% & 	0.047	\%\\
two-wheel 	&	0.002	\% & 	0.006	\% & 	0.007	\% & 	0.011	\%\\
rider		&	0.001	\% & 	0.004	\% & 	0.001	\% & 	0.006	\%\\
road		&	2.347	\% & 	2.029	\% & 	14.842	\% & 	12.942	\%\\
sidewalk	&	1.549	\% & 	1.277	\% & 	9.883	\% & 	8.444	\%\\
other-ground&	0.222	\% & 	0.153	\% & 	1.421	\% & 	1.022	\%\\
building	&	0.472	\% & 	0.471	\% & 	1.739	\% & 	1.988	\%\\
object		&	0.395	\% & 	0.110	\% & 	1.101	\% & 	0.399	\%\\
vegetation	&	2.043	\% & 	2.016	\% & 	6.284	\% & 	6.205	\%\\
trunk		&	0.025	\% & 	0.046	\% & 	0.062	\% & 	0.122	\%\\
terrain		&	1.374	\% & 	2.312	\% & 	7.717	\% & 	16.264	\%\\
unlabeled	&	91.290	\% & 	91.211	\% & 	56.078	\% & 	51.480	\%\\
        \hline
    \end{tabular}
\end{table}

\subsection{Dense Ground Truth}
The highly accurate sensor poses from the SemanticKITTI data set are used to enhance the ground truth density for each time point.
First, we create a batch of point clouds for each scan of a sequence which contains all point clouds of this sequence that possibly contribute to the semantic ground truth corresponding to the given scan.
For this, all point clouds with a sensor distance $|\Delta p_x|$ smaller than twice the maximum scan range $r$ are considered.
This leads to 40 to 300 aggregated point clouds per scan depending on the current velocity of the vehicle.
Highly accurate SLAM poses are required, as the GNSS poses from KITTI lead to alignment errors in the aggregated points cloud which make the semantic ground truth inconsistent.

Since moving objects cannot be aggregated over time using this method, only information from the current scan is used for this.
A sparse label map, consisting of points from all moving objects in the current scan, is superimposed on the dense grid map, keeping the knowledge of moving objects from the current scan.

\begin{table*}[!htbp]
\centering
 \caption{Quantitative results on the full evaluation set.}\label{tab::exp_results}
\label{my-label}
\begin{tabular}{llr|*{11}{p{0.75cm}}cc}
\toprule
&Approach &\textbf{mIoU}&vehicle&person&two-wheel&rider&road&side-walk&other\par-gro.&building&object&vege-tation&trunk&terrain\\
\midrule
\multirow{9}{*}{\rotatebox{90}{Sparse Eval}}
&$\mathbf{m3L_{i}}$& \textbf{0.325} & 0.455 & 0.0 & 0.0 & 0.0 & 0.810 & 0.532 & 0.114 & 0.649 & 0.117 & 0.593 & 0.000 & 0.626 \\
&$\mathbf{x41_{i}}$& \textbf{0.328} & 0.444 & 0.0 & 0.0 & 0.0 & 0.812 & 0.549 & 0.123 & 0.655 & 0.134 & 0.605 & 0.000 & 0.610  \\
&$\mathbf{x65_{i}}$& \textbf{0.324} &0.431 & 0.0 & 0.0 & 0.0 & 0.818 & 0.556 & 0.163 & 0.657 & 0.108 & 0.593 & 0.000 & 0.560 \\
&$\mathbf{m3L_{id}}$&\textbf{ 0.364} & 0.634 & 0.0 & 0.0 & 0.0 & 0.835 & 0.556 & 0.141 & 0.662 & 0.138 & 0.670 & 0.043 & 0.685   \\
&$\mathbf{x41_{id}}$& \textbf{0.363} & 0.650 & 0.0 & 0.0 & 0.0 & 0.837 & 0.569 & 0.165 & 0.687 & 0.146 & 0.653 & 0.009 & 0.642 \\
&$\mathbf{x65_{id}}$& \textbf{0.373} & 0.675 & 0.0 & 0.0 & 0.0 & 0.846 & 0.584 & 0.181 & 0.703 & 0.156 & 0.651 & 0.041 & 0.637  \\
&$\mathbf{m3L_{ido}}$& \textbf{0.373} & 0.661 & 0.0 & 0.0 & 0.0 & 0.836 & 0.554 & 0.150 & 0.703 & 0.142 & 0.665 & 0.100 & 0.666  \\
&$\mathbf{x41_{ido}}$& \textbf{0.382} & 0.689 & 0.0 & 0.0 & 0.0 & 0.852 & 0.587 & 0.209 & 0.708 & 0.154 & 0.664 & 0.065 & 0.653   \\
&$\mathbf{x65_{ido}}$& \textbf{0.398} & 0.697 & 0.0 & 0.0 & 0.0 & 0.858 & 0.603 & 0.259 & 0.728 & 0.151 & 0.689 & 0.099 & 0.693  \\
\addlinespace
\multirow{9}{*}{\rotatebox{90}{Dense Eval}}
&$\mathbf{m3L_{i}}$& \textbf{0.243} & 0.213 & 0.0 & 0.0 & 0.0 & 0.777 & 0.424 & 0.079 & 0.398 & 0.066 & 0.390 & 0.000 & 0.569  \\
&$\mathbf{x41_{i}}$&\textbf{ 0.260 }& 0.253 & 0.0 & 0.0 & 0.0 & 0.794 & 0.456 & 0.095 & 0.456 & 0.084 & 0.395 & 0.000 & 0.591   \\
&$\mathbf{x65_{i}}$& \textbf{0.246 }& 0.238 & 0.0 & 0.0 & 0.0 & 0.802 & 0.454 & 0.091 & 0.424 & 0.058 & 0.360 & 0.000 & 0.525  \\
&$\mathbf{m3L_{id}}$& \textbf{0.291} &0.408 & 0.0 & 0.0 & 0.0 & 0.807 & 0.460 & 0.128 & 0.450 & 0.100 & 0.488 & 0.035 & 0.620   \\
&$\mathbf{x41_{id}}$& \textbf{0.298} & 0.408 & 0.0 & 0.0 & 0.0 & 0.818 & 0.481 & 0.133 & 0.515 & 0.106 & 0.472 & 0.005 & 0.634 \\
&$\mathbf{x65_{id}}$& \textbf{0.303} & 0.430 & 0.0 & 0.0 & 0.0 & 0.826 & 0.498 & 0.155 & 0.521 & 0.102 & 0.459 & 0.028 & 0.622  \\
&$\mathbf{m3L_{ido}}$& \textbf{0.303} & 0.430 & 0.0 & 0.0 & 0.0 & 0.819 & 0.467 & 0.124 & 0.511 & 0.101 & 0.477 & 0.081 & 0.630\\
&$\mathbf{x41_{ido}}$& \textbf{0.312} & 0.435 & 0.0 & 0.0 & 0.0 & 0.835 & 0.494 & 0.157 & 0.547 & 0.112 & 0.486 & 0.040 & 0.638 \\
&$\mathbf{x65_{ido}}$& \textbf{0.328} & 0.433 & 0.0 & 0.0 & 0.0 & 0.843 & 0.514 & 0.229 & 0.547 & 0.108 & 0.510 & 0.063 & 0.686 \\

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

\section{EXPERIMENTS}
\label{sec:EXPERIMENTS}

In this work, we present the results of nine different deep learning models, which are all based on the DeeplabV3 \cite{Chen2018} architecture. This section describes the general experimental setup and how the nine models differ.

\subsection{Experimental Setup}

\subsubsection{General}
In the training process, we purely consider sparse, single-shot laser scans - this holds for both input data and for the ground truth. For all experiments, we train on the full grid map resolution of $1001 \times 501$ pixels and use the full label map of 12 classes as shown in Table~\ref{tab::statistics}. Due to GPU limitations, we choose the batch size of the Xception models to two and trained the MobileNet model with a batch size of four. We trained all models for at least 200000 iterations. In the training process, all unlabeled ground truth cells are excluded from the optimization. The loss is scaled with respect to the share of unlabeled cells. %

\subsubsection{Data Augmentation}
In deep learning applications, augmentation methods avoid the network to overfit on the training data. Hence, they are considered to be crucial for the learning process as they increase the model's performance on the evaluation or test set \cite{Shorten2019}. In this work, we apply multiple augmentation methods to prevent the network from overfitting. The training data was randomly flipped  along the vertical axis and randomly scaled with a factor between $0.8$ and $1.2$. %
 \cite{Srivastava2014}.

\subsubsection{Initialization}
Generally, literature recommends to initialize a deep convolutional neural network instead of training it from scratch \cite{IanGoodfellowYoshuaBengio2015}. In this work, we use weights, pre-trained on ImageNet \cite{Krizhevsky2017}, which were released in the official DeepLab repository \cite{Chen2018}. As our input has up to five channels, we randomly initialize the weights of the very first layer. This is done for every number of input channels to ensure comparability between our experiments.

\subsection{Model Outline}
The models differ in the used backbone architecture and the amount of input channels which are also influencing the first network layers.

\subsubsection{Feature Extractor}
We consider three different feature extractors as backbone for the DeepLabV3 architecture. Two of them are based on the Xception model, a powerful, yet relative efficient network. It uses depthwise separable convolutions and is based on the empirically proven hypothesis of the inception modules, i.e. mapping the cross-channel correlation separately from the spatial correlations is effective and efficient. We refer to them as $\mathbf{x41}$ for the model with 41 layers and $\mathbf{x65}$ for the larger 65-layered network. Further, we chose the MobileNetV3-large \cite{Howard2019}, later called $\mathbf{m3L}$ because it is considered a state-of-the-art-network in terms of efficiency. %

\subsubsection{Definiton of Input Layer}
To fully evaluate the benefit of a multi-layer grid map representation, we choose three different input configurations, which are denoted by an index as follows:
\begin{itemize}
\item $\mathbf{i}$: Using \textit{intensity} as single layer input.
\item $\mathbf{id}$: Adding the height information, \textit{Min. detected height} and \textit{Max. detected height} to the  \textit{intensity} layer, resulting in three input channels.
\item $\mathbf{ido}$: Expanding the input size to five by adding the dense layers \textit{Observability} and \textit{Min. observed height} to the three layers above.
\end{itemize}

\section{EVALUATION}
\label{sec:EVALUATION}

\subsection{Metrics}
We choose the \textit{Intersection over Union} (IoU) \cite{Everingham2014}, also known as \textit{Jaccard index}, for the qualitative evaluation as it is a commonly used measure for semantic segmentation. The per-class IoU$_k$ is calculated by
\begin{equation}
\text{IoU}_k = \frac{T_{P_k} }{T_{P_k}  + F_{P_k}  +F_{N_k} },
\end{equation}
with k being one of 12 classes. Consequently, the  \textit{mean Intersection over Union}, mIoU, is calculated as follows:
\begin{equation}
\text{mIoU} = \frac{1}{|K|} \sum_{k \in K} \text{IoU}_k
\end{equation}
where $| K|$ is the the labelset's cardinality.

\subsection{Evaluation Modes}
Each experiment is evaluated using two different approaches. Note that the network always predicts a dense semantic grid map, independent of the evaluation method. For the first one, we take the sparse ground truth, which is generated by processing a single scan. This ground truth corresponds to the ground truths, which are used for the training process. We denote this mode as \textit{Sparse Eval} evaluation. The second method, referred to \textit{Dense Eval} evaluation, considers the dense ground truth. Instead of evaluating with the full dense grid map, we filter it first with the observability input layer, so that occluded areas are removed from the evaluation. A filtered prediction map can be seen in the middle of Fig.~\ref{fig:qualitative_results}.

\subsection{Experimental Results} \label{sec::EVAL::ExperimentalResults}
\begin{figure*}[!htbp]
    \renewcommand{\thesubfigure}{\arabic{subfigure}}
    \centering
    \begin{subfigure}{\textwidth}
    \parbox[c]{.03\linewidth}{\subcaption{}}
    \parbox[c]{\linewidth}{
        \includegraphics[width = 0.95\textwidth ]{figures/qualiResults1.pdf}
    }
    \end{subfigure}
    \\ \vspace*{1.0cm}
    \begin{subfigure}{\textwidth}
        \parbox[c]{.03\linewidth}{\subcaption{}}
        \parbox[c]{\linewidth}{
            \includegraphics[width = 0.95\textwidth ]{figures/qualiResults2.pdf}
        }
        \end{subfigure}
            \\
    \centering
    \caption{
        Qualitative results from the evaluation set, processed by $\mathbf{x65_{ido}}$\label{fig:augm_l1}.
        The network was trained on sparse, multi-layer input and sparse ground truth.
        We depict four out of five input layers in the left column - from top to down: \textit{intensity}, \textit{max. detected height}, \textit{observability} and \textit{min. observed height}. The middle column depicts the predictions as a sparse layer on top and a dense layer below, filtered by the observed area.
        The corresponding ground truth is depicted in the right column.
    }
    \label{fig:qualitative_results}
\end{figure*}
Table \ref{tab::exp_results} presents the mIoU scores obtained through the experiments. We observe high accuracy in the predictions of \textit{road}, \textit{sidewalk}, \textit{vegetation}, \textit{terrain} and \textit{building} up to an IoU of 85,8 \%. All of these classes are also well represented in the data set distribution (see Table \ref{tab::statistics}). In contrast, the three classes, which represent vulnerable road users cannot be learned at all. All of them are highly underrepresented, e.g. there are 2300 times more labels of the \textit{road} class than there are \textit{rider} labels. Another explanation might be that the grid map representation is not suitable for these classes. Furthermore, it can be observed that the Xception models outperform the MobileNet by a slight margin.

A closer look at the different input configurations reveals that a multi-layer input is highly beneficial for the task of semantic segmentation of grid maps. In the conducted experiments, both the additional height layers and the dense observation layers boost the IoU for each class. While some classes yield reasonable results with only the \textit{intensity} layer (e.g. \textit{road}, \textit{sidewalk}), the multi-layer approach seems to play a key role for others: For instance, the \textit{vehicle} class gains 20\% IoU in all experiments after adding the two height layers and \textit{trunk} seems to be unpredictable without the height information.

We now consider the results obtained in \textit{Dense Eval}. The IoUs do only slightly drop in comparison to \textit{Sparse Eval}, while at the same time a multitude of labels are available and considered in the evaluation. This means that many cells without corresponding detections in the input layers are still predicted correctly according to the dense ground truth. This evaluation shows that the grid map representation enables the dense prediction of semantic information, even if it was trained on sparse ground truth data.

Qualitatively, Fig.~\ref{fig:qualitative_results} shows the network's ability to accurately predict the general configuration of a road scene, the road course, drivable space and large, well-discernible objects such as vehicles or walls. Small or vertically aligned objects such as poles, trunks and especially pedestrians are predicted less accurately. The dense prediction from sparse input represents the scenery quite well with mostly smooth und continuous area shapes, which are similar to the dense ground truth.

\subsection{Inference Time}
We evaluated the processing times on an NVIDIA GeForce RTX 2080 Ti GPU with 11 GB graphics memory. The presented models are trained using the Tensorflow framework \cite{GoogleResearch2015} and executed in Python. Table \ref{tab::inference} summarizes the results. %

\begin{table}[h]
\centering
\caption {Comparison of model architectures with respect to inference time.}\label{tab::inference}
    \begin{tabular}{c@{\hspace{0.3cm}}ccc}
\toprule
     {}   &  \multicolumn{3}{c}{inference time in ms}\\
     \cmidrule(lr){2-4}
    {Input layers } & { $\mathbf{x41}$} & { $\mathbf{x65}$} & { $\mathbf{m3L}$}  \\ \midrule
    {$\mathbf{i}$} & $68,6$& $84,9$ &$32,7$\\
    {$\mathbf{id}$} & $69,2$& $85,7$ &$33.1$\\
    {$\mathbf{ido}$} & $70,0$& $87,3$ &$34,4$\\
    \bottomrule
\end{tabular}%
\end{table}

As shown in Section \ref{sec::EVAL::ExperimentalResults}, $\mathbf{m3L}$ yields classification results which are similar to the Xception networks. However, it is significantly faster than the $\mathbf{x41}$ and $\mathbf{x65}$ network, which makes it a promising network for mobile perception systems such as intelligent vehicles regarding our proposed classification task. We also see that the addition of more input layers has very little effect on the inference time, which can be explained by the fact that only the computational effort in the very first layer is influenced by the input's channel size.  %

\section{CONCLUSION}
\label{sec:CONCLUSION}
In this work, we presented an approach towards semantic grid map estimation based on sparse input data. Further we developed and presented a method to generate a dense semantic ground truth which can be used to refine and evaluate the network's output. By evaluating on sparse and dense ground truth, we demonstrated that our approach enables the learned prediction of dense labels given a sparse input. The experiments have shown that a multi-layer input boosts the performance regarding both sparse and dense prediction. While the algorithm achieves good classification results for most classes, it fails to predict vulnerable road users. This might be due to their under-representation in the created grid map data set.

As a next step, we will integrate the dense ground truth layer within the learning process to teach the model directly the mapping from sparse input to a dense representation. Furthermore, we aim to fuse sequential input data within the grid map framework to enrich the sparse input and enhance the semantic segmentation performance.

\section{Acknowledgment}
The  authors  thank  Daimler  AG  for  the  fruitful  collaboration and the support for this work. \printbibliography

\end{document}

redpajama_set_name: RedPajamaArXiv