# Koopman Operator Theory for Nonlinear Dynamic Modeling using Dynamic Mode
Decomposition

Gregory Snyder and Zhuoyuan Song *This work was supported by the U.S. National
Science Foundation under awards CISE/IIS-2024928 and OIA-2032522.G. Snyder and
Z. Song are with the Department of Mechanical Engineering, University of
Hawai‘i at Mānoa, Honolulu, HI, 96822 USA. E-mails: {snyderg,
<EMAIL_ADDRESS>
###### Abstract

The Koopman operator is a linear operator that describes the evolution of
scalar observables (i.e., measurement functions of the states) in an infinite-
dimensional Hilbert space. This operator theoretic point of view lifts the
dynamics of a finite-dimensional nonlinear system to an infinite-dimensional
function space where the evolution of the original system becomes linear. In
this paper, we provide a brief summary of the Koopman operator theorem for
nonlinear dynamics modeling and focus on analyzing several data-driven
implementations using dynamical mode decomposition (DMD) for autonomous and
controlled canonical problems. We apply the extended dynamic mode
decomposition (EDMD) to identify the leading Koopman eigenfunctions and
approximate a finite-dimensional representation of the discovered linear
dynamics. This allows us to apply linear control approaches towards nonlinear
systems without linearization approximations around fixed points. We can then
examine the fidelity of using a linear controller based on a Koopman operator
approximated system on under-actuated systems with basic maneuvers. We
demonstrate the effectiveness of this theory through numerical simulation on
two classic dynamical systems are used to show DMD methods of evaluating and
approximating the Koopman operator and its effectiveness at linearizing these
systems.

##  I Introduction

In this work we explore the application of the Koopman operator to find
finite-dimensional linear representation of nonlinear dynamical systems. The
Koopman operator is an infinite-dimensional operator that evolves the
functions that evolve the state of a dynamical system linearly. In this work,
the Koopman operator is computed using data-driven methods such as the dynamic
mode decomposition (DMD) algorithm. DMD uses linear measurements of the states
of a dynamical system to find the dominant modes of the underlying system from
data to reconstruct a linear representation of the evolution of the system.

The theory behind the Koopman operator was first postulated as an infinite-
dimensional linear operator for uncontrolled systems in the seminal work by
Koopman and Neumann[1], which recently has sparked an reemergence of interest
in this method of modeling systems in the past two decades. More recently,
methods of finding the Koopman operator have been looked into by using
discrete data-driven methods for fluid dynamics [2] through dimensional
reduction algorithms developed by Schmid et al. [3].

Brunton et al. explored the relationship between the Koopman operator and
explored multiple observable functions to form a Koopman subspace to develop a
Koopman operator [4]. In their work, they showed that the state matrix found
by approximating the Koopman operater for specific systems can be used to
generate control laws with linear quadratic regression (LQR). Further more,
Proctor et al. [5] expanded on their previous work on DMD to create a DMD
method with control to extract low-order control models from higher-
dimensional systems building the base DMD algorithm. This DMD variant with
control, called DMDc, is demonstrated to show positive results in the analysis
of infections disease data. Proctor et al. [6] introduced a method of finding
the Koopman operator of a system that takes into account system’s inputs and
control (KIC) based on prior work on DMDc. Work has also been done to improve
the accuracy of the Koopman operater by ’lifting’ the states of the system to
a set of observables such that a solution can be found where the data required
from DMD is limited, extending the DMD algorithm (EDMD) [7]. Korda et al. [8]
presented the use of lifting nonlinear dynamics on augmented states with EDMD
and explored finding control laws for systems using model predictive control
(MPC). There have also been work done to optimize the methods in which the
lifting of dynamics is done through the optimization of the dictionary of
functions used to lift the state variables through EDMD with dictionary
learning [9].

Applications of the Koopman operator have been explored in the field of
robotics, where the operator is used to develop closed-loop controllers for
pendulum systems [10]. Kaiser et al. [11] explores the application of the
Koopman operator theory in generating energy-based control using a DMD variant
with control, extended dynamic mode decomposition with control, and compared
it against another method of finding the Koopman operator called the Koopman
reduced order nonlinear identification and control (KRONIC). Conversly, the
Koopman operator is not a panacea to solving dynamical systems; a recent work
by Gonzalez et al. [12] aimed at demonstrating some of the weakness of the
Koopman operator and the DMD algorithms and provides alternate approaches to
reconstruction.

This paper will illustrate the mathematics behind the Koopman operator theory
and show how augmenting the states of a dynamical system can create accurate
linear models through this process. It then shows how the DMD algorithm can be
used to find a linear approximation of dynamical systems before exploring
variants of DMD to create more accurate approximations of the Koopman operator
when representing system dynamics and other methods of using DMD to generate
control of a system represented via the Koopman operator. Numerical
simulations of two classical dynamical systems, the inverted pendulum and
cart-pole systems, are then studied using the DMD algorithms.

##  II Preliminaries

This work pertains to the implementation of Koopman operator theory on
dynamical systems to propagate non-linear dynamical equations with a data-
driven approximation method. Conventionally, a non-linear dynamical system
consists of a set of states and a function or rule that governs how the states
propagate either forward in time or with respect to each other [13, 14]. This
can be described in a continuous and a discrete method.

For a generic continuous system

$\displaystyle\frac{d}{dt}\mathbf{x}(t)$
$\displaystyle=F(\mathbf{x}(t),t;\mu),$ (1)

$\mathbf{x}(t)\in\mathbb{R}^{n}$ is the vector holding the states of the
dynamical system at time $t$, $n$ is the number of states that define the
system, $\mu$ is a set of parameters for the system dynamics, and $F(\cdot)$
is the rule describing the evolution of the state in a continuous sense. These
continuous-time dynamics can also be modeled in a respective discrete-time
representation where the system can be evaluated at every finite time
interval, $\Delta t$, which can otherwise be viewed as
$\mathbf{x}_{k}=x(k\Delta t)$ with the subscript $k$. The evolution of a
dynamic system in a discrete-time flow map can be formally portrayed as

$\displaystyle\mathbf{x}_{k+1}$ $\displaystyle=f(\mathbf{x}_{k}),$ (2)

by collecting the states at time $t_{k}$, $k=1,2,...,m$, for $m$ time steps,
where $x_{k}$ is a $n$-dimensional column vector of system states and
$x_{k+1}$ is the states of the system in the time step following $x_{k}$ [6,
5].

These rules governing how the states of the system advance tend to be
nonlinear equations to best simulate most of the systems in practice. An
unfortunate quality of nonlinear systems is that their dynamics are difficult
to solve analytically. As a result, modern control practices tend to turn to
approximations in order to produce a high-fidelity controller for the system.
However, if a dynamical system can be expressed with by a linear rule, we can
then achieve more accurate predictions of how a system advances in time. We
intend to show this using the Koopman operator theory.

###  II-A Koopman Operator

The Koopman operator is a linear but infinite dimensional operator that can be
defined for an autonomous, discrete time, dynamical system. Unlike the
function $f$ shown earlier in (2), the Koopman operator advances the
measurements of the changes in dynamics over time[13, 6, 8]. Consider the
evolution of a dynamical system $\mathbf{x}_{k+1}=f(\mathbf{x}_{k})$ where the
rule $f$ maps the state space onto itself, i.e.,
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$; using Koopman operator theory we
define a different rule, $g:\mathbb{R}^{n}\rightarrow\mathbb{M}^{n_{y}}$,
where $n_{y}$ is the dimension of a nearly infinite column vector defining the
dimension of the observable of $\mathbf{x}$ at a given time step. This implies
that the Koopman operator is defined for all observables, meaning that the
Koopman operator $\mathcal{K}$ is also infinite-dimensional. We will discuss
how we can make this a more reasonable attribute later on in this discussion.
Here, $g$ is a real-valued, scalar, measurement function which is an element
of an infinite-dimensional Hilbert space called an observable. The Koopman
operator acts on this observable such that

$\displaystyle\mathcal{K}_{t}g=$ $\displaystyle g\circ F(\mathbf{x}(t)),$ (3)
$\displaystyle\mathcal{K}_{\Delta t}g(\mathbf{x}_{k})=$ $\displaystyle
g(f(\mathbf{x}_{k})).$ (4)

In (3), the Koopman operator evolves the observable with respect to time in
the continuous fashion while in (4) the Koopman operator evolves the discrete-
time dynamics with respect to $\Delta t$. More details on the connections
between these two representations can be found in [13, 6, 8, 11, 15].
Equations (3) and (4) allow us to define an analogue for continuous and
discrete-time dynamical systems respectively as

$\displaystyle\frac{d}{dt}g=$ $\displaystyle\mathcal{K}g,$ (5) $\displaystyle
g(\mathbf{x}_{k+1})=$ $\displaystyle\mathcal{K}_{\Delta t}g(\mathbf{x}_{k}).$
(6)

This work explores the application of the Koopman operator on discret-time
dynamical systems. Due to the linear nature of the Koopman operator, we can
perform an eigen decomposition of $\mathcal{K}$ such that

$\displaystyle\mathcal{K}\varphi_{j}(\mathbf{x}_{k})=\lambda_{j}\varphi_{j}(\mathbf{x}_{k}),$
(7)

where $\lambda$ and $\varphi$ are the Koopman eigenvalue and eigenvector
describing the evolution of the Koopman operator. Considering this, the
observable, $g(x)$, can be expanded as

$\displaystyle g(\mathbf{x})$
$\displaystyle=\begin{bmatrix}g_{1}(\mathbf{x})\\\ g_{2}(\mathbf{x})\\\
\vdots\\\ g_{i}(\mathbf{x})\end{bmatrix}$
$\displaystyle=\sum_{j=1}^{\infty}\varphi_{j}(\mathbf{x})\bf{v}_{j},$ (8)

where $\bf{v}$ is a coefficient called the Koopman mode associated with its
corresponding Koopman eigenvector. This allows us to consider Koopman modes as
a projection of the observable:

$\bf{v}_{j}=\begin{bmatrix}\langle\varphi_{j},g_{1}\rangle\\\
\langle\varphi_{j},g_{2}\rangle\\\ \vdots\\\
\langle\varphi_{j},g_{i}\rangle\end{bmatrix}.$ (9)

By combining (7) and (8), one can show the relationship between the
observable, $g$, and the Koopman operator, $\mathcal{K}$. However, for all
practical proposes, using an infinite-dimensional vector and operator is not
always feasible so an approximation of the Koopman operator, $\mathcal{K}$,
will be typically used and found from collected data on the system. This
approximation allows us to combine the terms in (7) and (8) to show the
relationship between the Koopman modes in $g$ and the Koopman eigenvalues in
$\mathcal{K}$

$\mathcal{K}g(\mathbf{x}_{k})\approx\mathcal{K}g(\mathbf{x}_{k})=g(f(\mathbf{x}_{k}))=\sum_{j=1}^{\infty}\lambda_{j}\varphi_{j}(\mathbf{x})\bf{v}_{j}.$
(10)

This shows that the Koopman operator is an iterative set of triples,
$\lambda_{j},\varphi_{j}$ and $\bf{v}_{j}$, all of which make up the Koopman
mode decomposition[13, 6, 8, 11, 16, 17].

###  II-B Dynamic Mode Decomposition (DMD)

DMD is a data-driven method of exploring the behavior of complex systems.
Using measurements from numerical simulations or in experiments, we can
attempt to find the most dominant dynamic characteristics of the system. DMD
does this by finding the dynamic modes or eigen modes and eigenvalues of the
proposed system.
DMD acts on the assumption that the state of a system is connected to the next
by

$\mathbf{x}_{k+1}=\mathbf{A}\mathbf{x}_{k},$ (11)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ and $A\in\mathbb{R}^{n\times n}$ and is
the matrix describing the evolution of the state in a continuous-time manner.
Simulated or experimental measurements for $\mathbf{x}_{k}$ are then collected
at regular time intervals of $\Delta t$ to become snapshots to be used in a
discrete time system. These snapshots are collected and stored in sequence
like the following:

$\displaystyle\bf{\mathbf{X}}=$
$\displaystyle\begin{bmatrix}\vline&&\vline&&&&\vline\\\
\mathbf{x}_{1}&&\mathbf{x}_{2}&&\ldots&&\mathbf{x}_{m-1}\\\
\vline&&\vline&&&&\vline\end{bmatrix},$ (12)
$\displaystyle\bf{\mathbf{X}^{\prime}}=$
$\displaystyle\begin{bmatrix}\vline&&\vline&&&&\vline\\\
\mathbf{x}_{2}&&\mathbf{x}_{3}&&\ldots&&\mathbf{x}_{m}\\\
\vline&&\vline&&&&\vline\end{bmatrix},$ (13)

where $\bf{\mathbf{X}^{\prime}}$ is the time-shifted snapshot of matrix
$\bf{\mathbf{X}}$ such that

$\displaystyle\bf{\mathbf{X}^{\prime}}\approx A\bf{\mathbf{X}}.$ (14)

There is no set number of snapshots required for DMD but a sufficiently large
number is needed for the application and is closely related to the
approximation of the Koopman operator. For DMD we are attempting to find the
$\mathbf{A}$ in (14) using the snapshots of the system’s states. $\mathbf{A}$
can be approximated by

$\displaystyle\mathbf{A}=$ $\displaystyle\mathbf{X}^{\prime}\mathbf{X}^{+},$
(15)

where + is the Moore-Penrose pseudoinverse. Considering that $A$ has a large
$n$ that normal calculation would be a large computational load we can perform
Singular Value Decomposition (SVD) on the snapshots to find the dominant
characteristics of the pseudoinverse of $\bf{\mathbf{X}}$

$\displaystyle\bf{\mathbf{X}}\approx U\Sigma V^{*},$ (16)

where $\mathbf{U}\in\mathbb{R}^{n\times r}$,
$\mathbf{\Sigma}\in\mathbb{R}^{r\times r}$ and
$\mathbf{V}\in\mathbb{R}^{n\times r}$ where $*$ denotes the conjugate
transpose. $r$ is the reduced rank of the SVD approximation of
$\bf{\mathbf{X}}$. From SVD we can rearrange the singular values and the
eigenvectors of $\bf{\mathbf{X}}$ with the forward snapshot to $A$ that
satisfies (14):

$\displaystyle\mathbf{A}\approx\bar{\mathbf{A}}=$
$\displaystyle\bf{\mathbf{X}}^{\prime}\tilde{V}\tilde{\Sigma}^{-1}\tilde{U}^{*}.$
(17)

However in practice since the dimensionality of $A$ is so large we have to
approximate it as

$\displaystyle\tilde{\mathbf{A}}=\mathbf{U}^{*}\mathbf{X}^{\prime}V\mathbf{\Sigma}^{-1},$
(18)

where $\tilde{A}$ is a $r$ rank linear model of the dynamical system such that

$\displaystyle\tilde{\mathbf{x}}_{k+1}=$
$\displaystyle\tilde{\mathbf{A}}\tilde{\mathbf{x}}_{k}.$ (19)

Here $\tilde{\mathbf{x}}:\in\mathbb{R}^{r\times 1}$ is also a reduced order of
the state $\mathbf{x}$ which can be reconstructed by

$\displaystyle\tilde{\mathbf{x}}:\mathbf{x}=\tilde{\mathbf{U}}\tilde{\mathbf{x}}.$
(20)

The next step is to find the spectral decomposition of $\tilde{\mathbf{A}}$:

$\displaystyle\tilde{\mathbf{A}}\mathbf{W}=\mathbf{W}\mathbf{\Lambda},$ (21)

where $\mathbf{\Lambda}$ and $\mathbf{W}$ are the corresponding eigenvalues
(DMD modes) and eigenvectors of the full rank system dynamics matrix allowing
us to recover the full state system dynamics in a computationally efficient
manner since $\tilde{\mathbf{A}}\in\mathbb{R}^{r\times r}$ where $r<<n$. We
can now reconstruct the eigendecomposition in (16) for $\mathbf{A}$ from
$\mathbf{W}$ and $\mathbf{\Lambda}$ with the corresponding eigen vectors given
by the columns of $\mathbf{\Phi}$ such that

$\displaystyle\mathbf{\Phi}=$
$\displaystyle\mathbf{X}^{\prime}\mathbf{V\Sigma}^{-1}\mathbf{W}.$ (22)

From $\mathbf{\Phi}$ we can no reconstruct our approximation of the time
dynamics of $\mathbf{x}(t)$ by projecting our approximations into a future
solution

$\displaystyle\mathbf{x}(t)\approx\sum_{k=1}^{r}\boldsymbol{\phi}_{k}\exp{(\omega_{k}t)}b_{k}=\mathbf{\Phi}\exp{(\mathbf{\Omega}t)}b,$
(23)

where $b_{k}$ is the initial amplitude of each DMD mode, $\boldsymbol{\phi}$
is the columns that make up $\mathbf{\Phi}$ and $\mathbf{\Omega}$ is the
diagonal matrix of the eigenvalues $\omega$ in
$\omega_{k}=\ln{(\lambda_{k})}/\Delta t$[14, 13, 15, 5, 18].

###  II-C Extended Dynamic Mode Decomposition

Extended DMD is almost the same algorithm as the standard DMD one, however the
method of which we deploy EDMD is by using the observables of the system to
create a dictionary to pass through the normal DMD algorithm.[7, 13, 14, 19]
By performing regression on this new augmented vector containing linear and
non-linear measurements we can make a different approximation of the original
system dynamics. This augmented vector is constructed in the following manner

$\mathbf{y}=\mathbf{\Theta}^{T}(\mathbf{x})=\begin{bmatrix}\boldsymbol{\theta}_{1}(\mathbf{x})\\\
\boldsymbol{\theta}_{2}(\mathbf{x})\\\ \vdots\\\
\boldsymbol{\theta}_{p}(\mathbf{x})\end{bmatrix},$ (24)

where $p$ is the rank of the augmented state such that $p>>n$. Here
$\mathbf{\Theta}$ is the collection of measurements of the system possibly
containing the original state of the system, $\mathbf{x}$, as well as
nonlinear measurements. Once $\mathbf{y}$ is found, two data matrices are
created in the same manner seen above in the DMD algorithm (12) and (13). From
here a best-fit linear operator $\mathbf{A}_{Y}$ is found that maps (12) onto
(13)

$\mathbf{A}_{Y}=\operatorname*{argmin}_{\mathbf{A}_{Y}}\left\|{\mathbf{Y}^{\prime}-\mathbf{A}_{Y}\mathbf{Y}}\right\|=\mathbb{Y^{\prime}Y^{+}}.$
(25)

This regression can than be written in terms of the original data matrices
$\mathbf{\Theta}^{T}(\mathbf{x})$ as

$\mathbf{A}_{Y}=\operatorname*{argmin}_{\mathbf{A}_{Y}}\left\|\mathbf{\Theta}^{T}(\mathbf{x}^{\prime})-\mathbf{A}_{Y}\mathbf{\Theta}^{T}(\mathbf{x})\right\|=\mathbf{\Theta}^{T}(\mathbf{x}^{\prime})(\mathbf{\Theta}^{T}(\mathbf{x}))^{+}.$
(26)

$\mathbf{A}_{Y}$ is then the basis upon which we can derive the Koopman
operator. However, $\mathbf{\Theta}$ may not necessarily span the same
subspace as the Koopman operator and may consist of different eigenvalues and
eigenvectors of the Koopman operator when is why verification and re-
validation techniques need to be used to show that the EDMD model is properly
fit to the actual system.

###  II-D Approximating Koopman Eigenfunctions from Data

Beginning with the observable matrices in (24), the Koopman eigenfunction can
be approximated as

$\varphi(\mathbf{x})\approx\sum_{k=1}^{p}\boldsymbol{\theta}_{k}(\mathbf{x})\boldsymbol{\xi}_{k}=\mathbf{\Theta}(\mathbf{x})\boldsymbol{\xi},$
(27)

where by using (24), we can find

$\begin{bmatrix}\lambda\varphi(\mathbf{x}_{1})\\\
\lambda\varphi(\mathbf{x}_{2})\\\ \vdots\\\ \lambda\varphi(\mathbf{x}_{m})\\\
\end{bmatrix}=\begin{bmatrix}\varphi(\mathbf{x}_{2})\\\
\varphi(\mathbf{x}_{3})\\\ \vdots\\\
\varphi(\mathbf{x}_{m+1})\par\end{bmatrix}.$ (28)

We start seeing a connection between the Koopman operator and the process we
followed to approximate a linear evolution seen in DMD. We can expand (27) as

$\displaystyle[\lambda\mathbf{\Theta}(\mathbf{X})-\mathbf{\Theta}(\mathbf{X}^{\prime})]\boldsymbol{\xi}=0.$
(29)

Afterwards, we reduce (29) using a best least-squares fit to get

$\displaystyle\lambda\boldsymbol{\xi}=\mathbf{\Theta}(\mathbf{X})^{+}\mathbf{\Theta}(\mathbf{X}^{\prime})\boldsymbol{\xi}.$
(30)

Now we can compare (26) and (30); (26) is the transpose of the latter so that
the left eigenvectors are now the right eigenvectors. Comparing to (15) we can
use the eigenvectors $\boldsymbol{\xi}$ of
$\mathbf{\Theta}^{+}\mathbf{\Theta}^{\prime}$ to find the coefficients of the
eigenfunction $\varphi(\mathbf{x})$ that represents the basis of
$\mathbf{\Theta}(\mathbf{x})$. We now can confirm that the predicted
eigenfunctions behave linearly by comparing them to the predicted dynamics in
(28). We do this to determine if the regression done in (30) yields proper
eigenvalues and eigenvectors that span the Koopman invariant subspace for the
system[7].

###  II-E Methodology

The dynamical system model and the methods we used to discover the Koopman
operator are presented in the section. As discussed in the previous section,
the method that will be used to determine a finite-dimensional approximation
of the Koopman operator will be through DMD and EDMD. This approach has been
used in [20, 6, 11]. As stated in the previous section, DMD computes the
transformation matrix $\mathbf{A}$ that shows the evolution of two measurement
snapshot pairs by finding the corresponding eigenvectors and eigenvalues that
satisfy

$\displaystyle\mathbf{A}\mathbf{v}_{j}=$
$\displaystyle\lambda_{j}\mathbf{v}_{j}.$ (31)

Now assuming that $\mathbf{A}$ contains a complete set of eigenvectors, each
measurement can be found by the eigenvectors of $\mathbf{A}$ like so:

$\displaystyle g(x_{k})=\sum_{j=1}^{n}c_{jk}\mathbf{v}_{j}.$ (32)

In the case that these snapshot or augmented pairs in the case of EDMD and
DMD, the DMD modes and eigenvalues of (32) correspond to the Koopman modes of
(10). Since the principle concept of the Koopman operator, in the realm of
data-driven identification, is to advance the states of the system by
advancing the dynamics linearly.

##  III Results and Analysis

In this work we will numerically evaluate controlled and uncontrolled systems
to collect sufficiently large sets of data, convert it into snapshot pairs to
run through the DMD and EDMD algorithms and use the decomposition to
reconstruct linear models of the data to compare against the simulated
systems.

###  III-A Sample Systems

The classical systems that were used in this work to approximate the Koopman
operator are the inverted pendulum and cart pole system [13]. These systems
are used due to their periodic and transient behaviors and dynamics are well
known.

####  III-A1 Inverted Pendulum

Figure 1:  Classical inverted pendulum system with mass $m$ = 1 kg and arm
length of $L$ = 2 m. In this model the mass is free to swing about the pivot
360 degrees and will be treated as an undamped system.

In the inverted pendulum system seen in figure 1 we will start with an initial
state of $[\pi/4,0]$ and in the controlled examples we attempted to right the
pendulum mass to a terminal state of $[0,0]$. In this simulation we only
controlled the $\dot{\theta}$ term with a performance index of $[0,10]$ and a
control cost of 1 for our LQR gain.

####  III-A2 Cart-Pole

Figure 2:  Classical cart pole system with pendulum mass $m$ = 1kg, cart mass
$M$ = 5 kg, and arm length $L$ = 2 m. Like the pendulum system in Fig. 1, the
end mass is free to swing about the pivot 360 degrees and both the cart and
arm pivot are friction-less and undamped.

In the cart-pole system seen in Fig. 2 we started with an initial state of
$[-1,0,\pi,0]$ and attempted to right the cart-pole system after a brief
translation to $[1,0,\pi,0]$. In this simulation we created a control law
using LQR with a performance index of $[5,10,0,0]$ and a control cost of 1.

###  III-B Implementation

Figure 3: DMD performed on a pendulum system with uncontrolled data. Figure 4:
DMD performed on a pendulum system with with LQR controlled data.

We can see in Figs. 3 and 4 that the reconstruction of the original state
space succeeded with minimal errors. The periodic nature of the pendulum
system can still be seen in the error analysis and it recedes in time. This is
most likely due to the immediate transient behavior of the reference data,
though the over all behavior appears to be accurately represented. We can see
consistent low error in the reconstructed controlled and uncontrolled pendulum
data.

Figure 5: EDMD performed on a pendulum system with uncontrolled data lifted by
a second-order polynomial basis. Figure 6: EDMD performed on a pendulum system
with with LQR controlled data lifted by a second-order polynomial basis.

In Figs. 5 and 6 there is the repeated behavior of the DMD reconstruction seen
in Figs. 3 and 4. In both sets of the EDMD reconstruction we can see similar
or reduced error in reconstruction.

Figure 7: DMD performed on a cart-pole system with uncontrolled data. Figure
8: DMD performed on a cart-pole system with with LQR controlled data.

In Figs. 7 and 8 we can see significant losses in reconstructing the behavior
of the uncontrolled cart-pole system. Figure 8 shows that some of the
characteristics of the reference data have been captured in the reconstruction
as the reference data in this process is smoother and has less periodic and
transient characteristics.

###  III-C Discussions

We can see in the reconstruction results that the inverted pendulum system is
a good candidate for the system reconstruction using the approximated Koopman
operator. We can see a marginal improvement in the reconstruction of the state
data of the controlled data when we improve the DMD algorithm to EDMD thought
the improvement yet no change in the uncontrolled data. This shows that the
periodic behavior is not improved by a polynomial lifting function.
Additionally when moving to a higher order polynomial (or Fourier; see Figs.
11 and 12) we can additional improvements in the reconstruction of the
original data when the full lifted states are used. Oddly enough, as seen in
Fig. 13, we see no further improvement in the overall fidelity of the
reconstruction of the state space using a second-order Fourier fitting.

Figure 9: EDMD performed on a pendulum system with uncontrolled data lifted by
a third-order polynomial basis with no truncation. Figure 10: EDMD performed
on a pendulum system with with LQR controlled data lifted by a third-order
polynomial basis with no truncation.

Moving on to the second system in this work, we can see significant errors in
the reconstruction process as we are adding more transient components into the
system dynamics with the extra degrees of freedom of the cart-pole system. As
shown in Fig. 7, DMD (and by extension, EDMD) did not capture the uncontrolled
behavior of the system in question however when the strong periodic behavior
is removed in the controlled data set we find an improved reconstruction of
the system however weaker that that shown in the inverted pendulum system.
Additionally we note that the initial response generated by LQR appears to be
missed by the reconstruction possibly showing that more gradual control may
yield better results in the reconstruction. Additionally we can see no further
improvement on the reconstruction of the system when we attempted to use the
EDMD algorithm in Figs. 15 and 16 using either a polynomial or Fourier basis
to lift the original state space.

Figure 11: EDMD performed on a pendulum system with uncontrolled data using a
first-order Fourier basis to lift the state data.  Figure 12: EDMD performed
on a pendulum system with LQR controlled data using a first-order Fourier
basis to lift the state data.  Figure 13: EDMD performed on a pendulum system
with uncontrolled data using a second-order Fourier basis to lift the state
data.  Figure 14: EDMD performed on a pendulum system with LQR controlled data
using a second-order Fourier basis to lift the state data.  Figure 15: EDMD
performed on a cart-pole system with LQR controlled data using a second-order
polynomial basis to lift the state data.  Figure 16: EDMD performed on a
cartpole system with LQR controlled data using a first-order Fourier basis to
lift the state data.

##  IV Conclusions

This work presented an introduction on the usage of the Koopman operator
theory approach to model simple, nonlinear dynamical systems with a linear
operator. In this process, we simulated the states of classical non-linear
systems and feed that data as a set of snapshots through the DMD algorithm,
and its variants, to reconstruct the state data using the Koopman operator.
Through the usage of these algorithms on an inverted pendulum and cart-pole
system we can approximate the Koopman operator to recreate high fidelity
linear approximations of these systems and provide increased accuracy by
augmenting the states of the original data. These methods showed weakness in
reconstructing systems that act with coupled states that behave in different
matters as seen in the cart-pole system. In future work, we intend to explore
alternate lifting functions to represent the state spaces of our systems and
explore alternative methods of finding the Koopman operator of systems with
multiple coupled states to improve performance to develop linear control laws.

## References

  * [1] B. O. Koopman and J. V. Neumann, “Dynamical systems of continuous spectra,” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 18, no. 3, pp. 255–263, 1932. [Online]. Available: http://www.jstor.org/stable/86259
  * [2] C. Rowley, I. Mezic, S. Bagheri, P. Schlatter, and D. Henningson, “Spectral analysis of nonlinear flows,” _Journal of Fluid Mechanics_ , vol. 641, pp. 115 – 127, 12 2009.
  * [3] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” _Journal of Fluid Mechanics_ , vol. 656, p. 5–28, 2010.
  * [4] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz, “Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control,” _PloS one_ , vol. 11, no. 2, p. e0150171, 2016.
  * [5] J. Proctor, S. Brunton, and J. Kutz, “Dynamic mode decomposition with control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 15, 09 2014\.
  * [6] ——, “Generalizing Koopman theory to allow for inputs and control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 17, 02 2016.
  * [7] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the Koopman operator: Extending dynamic mode decomposition,” _Journal of Nonlinear Science_ , vol. 25, no. 6, p. 1307–1346, Jun 2015. [Online]. Available: http://dx.doi.org/10.1007/s00332-015-9258-5
  * [8] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” _Automatica_ , vol. 93, p. 149–160, Jul 2018. [Online]. Available: http://dx.doi.org/10.1016/j.automatica.2018.03.046
  * [9] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis, “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 10, p. 103111, Oct 2017. [Online]. Available: http://dx.doi.org/10.1063/1.4993854
  * [10] I. Abraham, G. de la Torre, and T. Murphey, “Model-based control using Koopman operators,” _Robotics: Science and Systems XIII_ , Jul 2017. [Online]. Available: http://dx.doi.org/10.15607/RSS.2017.XIII.052
  * [11] E. Kaiser, J. Kutz, and S. Brunton, “Data-driven discovery of Koopman eigenfunctions for control,” _Machine Learning: Science and Technology_ , vol. 2, 07 2017.
  * [12] E. Gonzalez, M. Abudia, M. Jury, R. Kamalapurkar, and J. A. Rosenfeld, “Anti-Koopmanism,” 2021.
  * [13] S. Brunton and J. Kutz, _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_. Cambridge University Press, 2019. [Online]. Available: https://books.google.com/books?id=gNcRuQEACAAJ
  * [14] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_. Philadelphia, PA, USA: SIAM-Society for Industrial and Applied Mathematics, 2016.
  * [15] A. Mauroy, I. Mezic, and Y. Susuki, _The Koopman Operator in Systems and Control Concepts, Methodologies, and Applications: Concepts, Methodologies, and Applications_ , 01 2020.
  * [16] I. Mezic, “Spectral properties of dynamical systems, model reduction and decompositions,” _Nonlinear Dynamics_ , vol. 41, pp. 309–325, 08 2005.
  * [17] ——, “Analysis of fluid flows via spectral properties of the koopman operator,” _Annual Review of Fluid Mechanics_ , vol. 45, pp. 357–378, 01 2013.
  * [18] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. Nathan Kutz, “On dynamic mode decomposition: Theory and applications,” _Journal of Computational Dynamics_ , vol. 1, no. 2, p. 391–421, 2014. [Online]. Available: http://dx.doi.org/10.3934/jcd.2014.1.391
  * [19] M. O. Williams, C. W. Rowley, and I. G. Kevrekidis, “A kernel-based approach to data-driven Koopman spectral analysis,” 2015.
  * [20] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, p. 3932–3937, Mar 2016. [Online]. Available: http://dx.doi.org/10.1073/pnas.1517384113


# Quantization in Spiking Neural Networks

Bernhard A. Moser
Institute of Signal Processing
Johannes Kepler University of Linz
bernhard.moser@{scch.at,jku.at} Michael Lunglmayr
Institute of Signal Processing
Johannes Kepler University of Linz, Austria
<EMAIL_ADDRESS>double affiliation: Software Competence Center
Hagenberg (SCCH), 4232 Hagenberg, Austria

###### Abstract

In spiking neural networks (SNN), at each node, an incoming sequence of
weighted Dirac pulses is converted into an output sequence of weighted Dirac
pulses by a leaky-integrate-and-fire (LIF) neuron model based on spike
aggregation and thresholding. We show that this mapping can be understood as a
quantization operator and state a corresponding formula for the quantization
error by means of the Alexiewicz norm. This analysis has implications for
rethinking reinitialization in the LIF model, leading to the proposal of
reset-to-mod as a modulo-based reset variant.

_Keywords_ Leaky-Integrate-and-Fire (LIF) Neuron $\cdot$ Spiking Neural
Networks (SNN) $\cdot$ Re-Initialization $\cdot$ Quantization $\cdot$ Error
Propagation $\cdot$ Alexiewicz Norm

##  1 Introduction

Though its simplicity, the leaky integrate-and-fire (LIF) neuron model is
widely spread in neuromorphic computing and computational neuroscience
Gerstner et al. (2014); Nunes et al. (2022). In contrast to more biophysically
realistic models, such as the Hodgkin-Huxley model, LIF is a middle ground to
capture essential features of bio-inspired time-based information processing,
and simple enough to be applicable from the point of view of neuromorphic
engineering. In particular, in biology spikes show varying shapes and its
generation and re-initialization follows a complex dynamics as revealed by the
Hodgkin-Huxley differential equations. In contrast, LIF is based on the
following major idealizations. (1), it neglects shape information and
idealizes spikes as shapeless impulses, (2), it is based on the assumption
that the triggering process is realized by a comparison with a given
threshold-level, and, (3), the process of triggering and re-initialization
acts instantaneously.

In this paper we take these assumptions as axioms for a mathematical operator
that maps a sequence of weighted Dirac impulses to another one. We ask about
properties of this mapping and find that it can be understood as a
quantization operation in the space of spike trains. The standard model of
quantization of a single number is integer truncation. That is, a given real
number $x=n+r$, $n\in\mathbb{Z}$, $r\in(-1,1)$ is mapped to the integer $n$
and the quantization error is given by $r$. The resulting mapping
$q:\mathbb{R}\rightarrow\mathbb{Z}$ can also be considered from a geometric
point of view. Consider a grid of vertexes of a tessellation of polytopes
induced by the unit ball of the maximum norm, $\|.\|_{\infty}$. This way, a
given point $x$ is element of such a polytope $P$ and mapped to that vertex of
$P$ that is closest to the origin.

In an abstract sense, we will show that LIF can be understood in the same way.
Instead of the maximum norm, $\|.\|_{\infty}$, the space $\mathbb{S}$ of spike
trains is tessellated by means of the unit balls of the Alexiewicz norm
$\|.\|_{A,\alpha}$, which will be outlined next. Then we get the formula
$\|\mbox{LIF}_{\alpha,\vartheta}(\eta)-\eta\|_{A,\alpha}<\vartheta$, where
$\mbox{LIF}_{\alpha,\vartheta}$ denotes the LIF mapping based on the leaky
parameter $\alpha$, the threshold $\vartheta>0$ and the input spike train
$\eta\in\mathbb{S}$.

The paper is outlined as follows. In Section 2 we reformulate the LIF model,
taking different re-initialization variants into account. In Section 3 we
recall the Alexiewicz norm $\|.\|_{A}$, generalize it to leaky variants
$\|.\|_{A,\alpha}$ and state the main theorem. Section 4 presents evaluations
considering the discussed reinitialization variants, based on Python code
available at https://github.com/LunglmayrMoser/AlexSNN.

##  2 LIF Model and Re-Initialization Variants

The flow of information in a LIF model is outlined in Fig. 1.

Figure 1: Illustration of steps of information processing by a LIF model. The
weighted superposition of incoming spike trains can result in spike amplitudes
that can exceed the threshold many times.

Arriving from several input channels, spike events trigger the dynamics of the
membrane potential by a weighted linear superposition across the channels and
an integration with a leaky term in time. The LIF model idealizes the
superposition (and, also that of re-initialization) by an instantaneous
aggregation. Due to this building process, the incoming spike events can
equivalently be represented by a single channel of spikes with amplitudes
resulting from a weighted sum of simultaneous spikes across the input
channels. This way we have to take spike amplitudes of virtually arbitrary
magnitude into account. In the literature, this aspect is pointed out as a
limitation of the LIF model, since it deviates from biology in this respect
Gerstner et al. (2014). As an idealized mathematical model, however, it has
its justification, but then this effect must also be taken into account in the
re-initialization step. While in the context of LIF there are discussed
basically two modes of re-initialization (synonymously, reset), namely reset-
to-zero and reset-by-subtraction, in Eqn. (1) we will propose a third one,
reset-to-mod, to be mathematically consequent regarding the outlined issue of
instantaneous superposition. According to Eshraghian et al. (2021), reset-to-
zero means that the potential is reinitialized to zero after firing, while
reset-by-subtraction subtracts the $\vartheta$-potential $u_{\vartheta}$ from
the membrane’s potential that triggers the firing event. As a third variant we
introduce reset-to-mod, which can be understood as instantaneously cascaded
application of reset-by-subtraction according to the factor by which the
membrane’s potential exceeds the threshold which results in a modulo
computation.

For an input spike train $\eta_{in}(t)=\sum_{i}a_{i}\delta_{t_{i}}(t)$ the
mapping
$\sum_{i}b_{i}\,\delta_{s_{i}}=\mbox{LIF}_{\vartheta,\alpha}(\eta_{in})$ is
recursively given by $s_{i+1}=\inf\left\\{s\geq
s_{i}:\,\left|u_{\vartheta,\alpha}(s_{i},s)\right|\geq\vartheta\right\\},$
where

$u_{\vartheta,\alpha}(t_{\mbox{event}},t):=\int^{t}_{t_{\mbox{event}}}e^{-\alpha(\tau-
t_{\mbox{event}})}\left(\eta_{in}(\tau)-\mbox{discharge}(t_{\mbox{event}},\tau)\right)d\tau$

models the dynamic change of the neuron membrane’s potential after an input
spike event at time $t_{\mbox{event}}$. The process of triggering an output
spike is actually a charge-discharge event

$\mbox{discharge}(t_{i},\tau):=\left\\{\begin{array}[2]{lcl}a_{i}\delta_{t_{i}}(\tau)&\ldots&\mbox{for
{\it reset-to-zero}},\\\
\mbox{sgn}(a_{i})\,\vartheta\,\delta_{t_{i}}(\tau)&\ldots&\mbox{for {\it
reset-by-subtraction}},\\\
q(a_{i}/\vartheta)\,\vartheta\,\delta_{t_{i}}(\tau)&\ldots&\mbox{for {\it
reset-to-mod}}\end{array}\right.$ (1)

where $\mbox{sgn}(x)\in\\{-1,0,1\\}$ is the signum and
$q(x):=\mbox{sgn}(x)\max\\{k\in\mathbb{Z}:k\leq|x|\\}$ is the integer
truncation quantization function.

Figure 2: Example of input spike train with different LIF outputs depending on
the re-initialization variant.

The reset-by-subtraction mode can be understood as compensation event so that
the net voltage balance of the spiking event equals zero, i.e., in case of an
output spike with amplitude $\vartheta>0$ the membrane is actually discharged
by this amount. Accordingly, though not always made clear in the literature,
see for example Eshraghian et al. (2021), this assumption has the subtle
consequence that an increase of the membrane potential $u$ by multiples
$q(u/\vartheta)$ of the threshold level $\vartheta$ results in a discharge of
the membrane’s potential by the same amount, that is
$q(u/\vartheta)\vartheta$. See Fig. 2 for an illustration and example.

##  3 Alexiewicz Norm and Spike Train Quantization

For $\eta=\sum_{i}a_{i}\delta_{t_{i}}\in\mathbb{S}$ we introduce the measure
$\|\eta\|_{A,\alpha}:=\max_{n}\left|\sum_{j=1}^{n}a_{j}e^{-\alpha(t_{n}-t_{j})}\right|$,
which satisfies the axioms of a norm on the vector space $\mathbb{S}$. For
$\alpha=0$ we obtain the Alexiewicz norm Alexiewicz (1948), which is
topologically equivalent to the discrepancy norm Moser (2012). The Alexiewicz
norm reveals the quantization character of the LIF model.

######  Theorem 1 (reset-to-mod LIF as $\|.\|_{A,\alpha}$-Quantization)

Given a LIF neuron model with reset-to-mod, the threshold $\vartheta>0$, the
leaky parameter $\alpha\in[0,\infty]$ and the spike train $\eta\in\mathbb{S}$
with amplitudes $a_{i}\in\mathbb{R}$. Then,
$\mbox{LIF}_{\vartheta,\alpha}(\eta)$ is a $\vartheta$-quantization of $\eta$,
i.e., the resulting spike amplitudes are multiples of $\vartheta$, where the
quantization error is bounded by

$\|\mbox{LIF}_{\vartheta,\alpha}(\eta)-\eta\|_{A,\alpha}<\vartheta.$ (2)

### Proof.

Figure 3: Illustration of Eqn. (4). The red arrows indicate reset by reset-to-
mod.

First of all, we introduce the following operation $\oplus$, which is
associative and can be handled with like the usual addition if adjacent
elements $a_{i}$ from a spike train $\eta=\sum_{i}a_{i}\delta_{t_{i}}$ are
aggregated:

$a_{i}\oplus a_{i+1}:=e^{-\alpha(t_{i+1}-t_{i})}a_{i}+a_{i+1}.$ (3)

This way we get a simpler notation when aggregating convolutions, e.g.,

$a_{i}\oplus\ldots\oplus a_{j}=\sum_{k=i}^{j}e^{-\alpha(t_{j}-t_{k})}a_{k}.$

For the discrete version we re-define $a_{i_{k}}\oplus
a_{i_{k+1}}:=\beta^{(i_{k+1}-i_{k})}a_{i_{k}}+a_{i_{k+1}}$, if $i_{k}$ and
$i_{k+1}$ refer to adjacent spikes at time $i_{k}$, resp. $i_{k+1}$. Further,
we denote $q[z]:=\mbox{sgn}(z)\,\left\lfloor|z|\right\rfloor$ which is the
ordinary quantization due to integer truncation, e.g. $q[1.8]=1$,
$q[-1.8]=-1$, where
$\left\lfloor|z|\right\rfloor=\max\\{n\in\mathbb{N}_{0}:\,n\leq|z|\\}$.

After fixing notation let us consider a spike train
$\eta=\sum_{j}a_{i}\delta_{t_{i}}$. Without loss of generality we may assume
that $\vartheta=1$. We have to show that
$\|\mbox{LIF}_{1,\alpha}(\eta)-\eta\|_{A,\alpha}<1$, which is equivalent to
the discrete condition that $\forall
n:\max_{n}\left|\sum_{i=1}^{n}\hat{a}_{i}\right|<1$, where
$\eta-\mbox{LIF}_{1,\alpha}(\eta)=\sum_{i}\hat{a}_{i}\delta_{t_{i}}$. Set
$\hat{s}_{k}:=\hat{a}_{0}\oplus\cdots\oplus\hat{a}_{k}$. We have to show that
$\max_{k}|\hat{s}_{k}|<1$. The proof is based on induction and leads the
problem back to the standard quantization by truncation.

Suppose that at time $t_{i_{k-1}}$ after re-initialization by reset-to-mod we
get the residuum $\Delta_{i_{k-1}}$ as membrane potential that is the starting
point for the integration after $t_{i_{k-1}}$. Note that

$\mbox{LIF}_{1,\alpha}(\eta)|_{t=t_{k}}=q(\Delta_{i_{k-1}}\oplus
a_{i_{k-1}+1}\cdots\oplus a_{i_{k}})$

Then, as illustrated in Fig. 3 the residuum $\Delta_{i_{k}}$ at the next
triggering event $t_{i_{k}}$ is obtained by the equation

$\Delta_{i_{k}}=\Delta_{i_{k-1}}\oplus a_{i_{k-1}+1}\oplus\ldots\oplus
a_{i_{k}}-q[\Delta_{i_{k-1}}\oplus\ldots\oplus a_{i_{k}}].$ (4)

Note that due to the thresholding condition of LIF we have

$|\Delta_{i_{k}}\oplus a_{i_{k}+1}\oplus\ldots\oplus a_{j}|<1$ (5)

for $j\in\\{i_{k}+1,\ldots,i_{k+1}-1\\}$. For the $\oplus$-sums
$\hat{s}_{i_{k}}$ we have

$\hat{s}_{i_{k+1}}=\hat{s}_{i_{k}}\oplus a_{i_{k}+1}\cdots
a_{i_{k+1}-1}\oplus\left(a_{i_{k+1}}-q[\Delta_{i_{k}}\oplus
a_{i_{k}+1}\oplus\ldots\oplus a_{i_{k+1}}]\right).$ (6)

Note that $\hat{s}_{0}=\Delta_{i_{0}}=a_{0}-q[a_{0}]$, then for induction we
assume that up to index $k$ to have

$\hat{s}_{i_{k}}=\Delta_{i_{k}}.$ (7)

Now, using (7), Equation (6) gives

$\displaystyle\hat{s}_{i_{k+1}}$ $\displaystyle=$
$\displaystyle\Delta_{i_{k}}\oplus a_{i_{k}+1}\oplus\ldots\oplus
a_{i_{k+1}-1}\oplus\left(a_{i_{k+1}}-q[\Delta_{i_{k}}\oplus
a_{i_{k+1}}\oplus\ldots\oplus a_{i_{k+1}}]\right)$ (8) $\displaystyle=$
$\displaystyle\Delta_{i_{k}}\oplus a_{i_{k}+1}\oplus\ldots\oplus
a_{i_{k+1}-1}\oplus a_{i_{k+1}}-q[\Delta_{i_{k}}\oplus
a_{i_{k+1}}\oplus\ldots\oplus a_{i_{k+1}}],$ $\displaystyle=$
$\displaystyle\Delta_{i_{k+1}}$

proving (7), which together with (5) ends the proof showing that
$|\hat{s}_{k}|<1$ for all $k$. $\square$

##  4 Evaluations

We evaluate the distribution of the quantization error (2) by means of box
whisker diagrams for different re-initialization variants depending on the
number of spikes in the spike train and different distributions of amplitudes.
Fig. 4 shows the result for incoming spike amplitudes that are below
threshold. In this case, (2) holds for reset-to-mod and reset-by-subtraction.
As expected, both variants behave similarly. For reset-to-zero the bound (2)
only holds approximately for larger leaky parameter $\alpha$, see first row of
Fig. 4. For reset-to-mod and reset-by-subtraction, With increasing number of
spikes one can observe a concentration of measure effect, see, e.g., Vershynin
).

Figure 4:  Evaluation of (2) for reset-to-mod, reset-by-subtraction and reset-
to-zero (1.st/2nd/3rd column), based on spike trains with spike amplitudes in
$[-\vartheta,\vartheta]$ and $100$ runs. The 1st row refers to $\alpha=1$ and
the 2nd row to $\alpha=0.1$.

However, as shown in Fig. 5, if the incoming spike amplitudes are not below
threshold anymore Eqn. (2) only holds for reset-to-mod. In this case the
measure of concentration effect become more apparent for smaller leaky
parameter. For the other variants the quantization error increases in average
with the number of spikes, and the theoretical bound (2), proven for reset-to-
mod, does not hold anymore.

Figure 5: The same as in Fig. 4 but with spike amplitudes in
$[-3/2\vartheta,3/2\vartheta]$.

##  5 Conclusion

In this paper we provide a novel view on the leaky integrate-and-fire model as
quantization operator in the Alexiewicz norm. This analysis gives rise to
rethinking the re-initialization modes reset-to-zero and reset-by-subtraction
that are commonly used in the context of spiking neural networks. These re-
initialization modes only hold under restricted conditions while our proposed
variant reset-to-mod satisfies the derived quantization bound under general
conditions. This general quantization error formula leads to new error bounds
for LIF and SNNs, such as a quasi-isometry relation in analogy to threshold-
based sampling Moser (2017); Moser and Lunglmayr (2019). Examples can be found
in the github repository https://github.com/LunglmayrMoser/AlexSNN and Moser
and Lunglmayr (2023).

## Acknowledgements

This work was supported (1) by the ’University SAL Labs’ initiative of Silicon
Austria Labs (SAL) and its Austrian partner universities for applied
fundamental research for electronic based systems, (2) by Austrian ministries
BMK, BMDW, and the State of Upper-Austria in the frame of SCCH, part of the
COMET Programme managed by FFG, and (3) by the NeuroSoC project funded under
the Horizon Europe Grant Agreement number 101070634.

## References

  * Gerstner et al. [2014] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.  _Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition_.  Cambridge University Press, USA, 2014.  ISBN 1107635195.
  * Nunes et al. [2022] João D. Nunes, Marcelo Carvalho, Diogo Carneiro, and Jaime S. Cardoso.  Spiking neural networks: A survey.  _IEEE Access_ , 10:60738–60764, 2022.  doi:10.1109/ACCESS.2022.3179968.
  * Eshraghian et al. [2021] Jason K. Eshraghian, Max Ward, Emre Neftci, Xinxin Wang, Gregor Lenz, Girish Dwivedi, Mohammed Bennamoun, Doo Seok Jeong, and Wei D. Lu.  Training spiking neural networks using lessons from deep learning.  _arXiv_ , 2021.  doi:10.48550/ARXIV.2109.12894.  URL https://arxiv.org/abs/2109.12894.
  * Alexiewicz [1948] Andrzej Alexiewicz.  Linear functionals on Denjoy-integrable functions.  _Colloq. Math._ , 1:289–293, 1948.
  * Moser [2012] Bernhard A. Moser.  Geometric characterization of Weyl’s discrepancy norm in terms of its n-dimensional unit balls.  _Discret. Comput. Geom._ , 48(4):793–806, 2012\.  doi:10.1007/s00454-012-9454-0.  URL https://doi.org/10.1007/s00454-012-9454-0.
  * [6] Roman Vershynin.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Number 47 in Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.  ISBN 978-1-108-41519-4.
  * Moser [2017] Bernhard A. Moser.  Similarity recovery from threshold-based sampling under general conditions.  _IEEE Transactions on Signal Processing_ , 65(17):4645–4654, 2017.  doi:10.1109/TSP.2017.2712121.
  * Moser and Lunglmayr [2019] Bernhard A. Moser and Michael Lunglmayr.  On quasi-isometry of threshold-based sampling.  _IEEE Transactions on Signal Processing_ , 67(14):3832–3841, 2019.  doi:10.1109/TSP.2019.2919415.
  * Moser and Lunglmayr [2023] Bernhard A. Moser and Michael Lunglmayr.  Spiking neural networks in the alexiewicz topology: A new perspective on analysis and error bounds.  _arXiv preprint arXiv:2305.05772_ , 2023.  URL https://doi.org/10.48550/arXiv.2305.05772.