On the Koopman operator of algorithms

TL;DR

This paper proposes a Koopman operator framework to analyze and accelerate numerical algorithms by viewing them as dynamical systems, enabling data-driven reduced representations and comparisons.

Contribution

It introduces a novel application of Koopman operator theory to numerical algorithms, facilitating their analysis, comparison, and potential acceleration through data-driven methods.

Findings

Koopman operator provides a unified framework for analyzing algorithms.

Basis functions adapted to data improve reduced representations.

Demonstrated on gradient descent, Nesterov, and Newton-Raphson algorithms.

Abstract

A systematic mathematical framework for the study of numerical algorithms would allow comparisons, facilitate conjugacy arguments, as well as enable the discovery of improved, accelerated, data-driven algorithms. Over the course of the last century, the Koopman operator has provided a mathematical framework for the study of dynamical systems, which facilitates conjugacy arguments and can provide efficient reduced descriptions. More recently, numerical approximations of the operator have enabled the analysis of a large number of deterministic and stochastic dynamical systems in a completely data-driven, essentially equation-free pipeline. Discrete or continuous time numerical algorithms (integrators, nonlinear equation solvers, optimization algorithms) are themselves dynamical systems. In this paper, we use this insight to leverage the Koopman operator framework in the data-driven study of such algorithms and discuss benefits for analysis and acceleration of numerical computation. For algorithms acting on high-dimensional spaces by quickly contracting them towards low-dimensional manifolds, we demonstrate how basis functions adapted to the data help to construct efficient reduced representations of the operator. Our illustrative examples include the gradient descent and Nesterov optimization algorithms, as well as the Newton-Raphson algorithm.

Tables

Table 1. Table 1: Algorithms formulated as dynamical systems, with their continuous (left) and discrete versions (right). Starting from a point x 0 ∈ ℝ n subscript 𝑥 0 superscript ℝ 𝑛 x_{0}\in\mathbb{R}^{n} , the function f : ℝ n → ℝ : 𝑓 → superscript ℝ 𝑛 ℝ f:\mathbb{R}^{n}\to\mathbb{R} (or its absolute value, in case of root-finding algorithms) is to be minimized, where ∇ f ∇ 𝑓 \nabla f denotes its gradient, and H x ​ f subscript 𝐻 𝑥 𝑓 H_{x}f the Hessian matrix at x 𝑥 x . If f 𝑓 f is vector-valued, J x ​ f subscript 𝐽 𝑥 𝑓 J_{x}f denotes its Jacobian matrix at x 𝑥 x . The symbol h ∈ ℝ + ℎ superscript ℝ h\in\mathbb{R}^{+} denotes a small, positive constant.

Gradient Decent
$\dot{x}=-\nabla f\left(x\right)$	$x_{n+1}=x_n-h\nabla f\left(x_n\right)$
Davidenko / Newton-Raphson method for optimization
$\left(H_{x}f\right)\dot{x}+\nabla f\left(x_0\right)=0$	$x_{n+1}=x_n-{\left(H_{x_n}f\right)}^{-1}\nabla f\left(x_n\right)$
Davidenko / Newton-Raphson method for root finding
$\left(J_{x}f\right)\dot{x}+f\left(x_0\right)=0$	$x_{n+1}=x_n-{\left(J_{x_n}f\right)}^{-1}f\left(x_n\right)$
Nesterov method
$\ddot{x}+\frac{r}{t}\dot{x}+\nabla f\left(x\right)=0,r\ge 3$	$x_{n+1}=y_n-h\nabla f\left(y_n\right)$
	$y_{n+1}=x_{n+1}+\frac{n}{n+3}\left(x_{n+1}-x_n\right)$

Equations

$$x_{n+1}=a(x_n),n\in \mathbb{N}.$$

$$\frac{d}{ds}c(s)=v(c(s)),c(0)=x_0.$$

$$S(t+s,x)=S(t,S(s,x)),{\frac{d}{dt}S(t,x)|}_{t=0}=v(x).$$

$${\frac{d}{dt}{S}_t(x)|}_{t=0}=v(x),S_0(x)=x_0,$$

$$[{\mathcal{K}}^{t}g](x):=(g\circ S_t)(x).$$

$$\mathcal{F}=L^2(X,\mathbb{C},\mu ):=\{g:X\to \mathbb{C},\text{s.t.}{\int }_X|g(x){|}^{2}d\mu (x)<\infty \}.$$

$$⟨g_1,g_2⟩:={\int }_X{g}_1(x)\overline{g_2(x)}d\mu (x),$$

$$t\mapsto {\mathcal{K}}^{t}f$$

$$\underset{t\to 0}{\lim }\Vert {\mathcal{K}}^t-I\Vert =0.$$

$$\mathcal{K}=\sum _k\underset{\in {\sigma }_{pp}}{\underbrace{{\lambda }_k}}{P}_{{\lambda }_k}+{\int }_{{\sigma }_{ac}}\lambda \text{d}E(\lambda ),$$

$$[{\mathcal{K}}^t{\varphi }_{\lambda }](x)=({\varphi }_{\lambda }\circ S_t)(x)={\lambda }^t{\varphi }_{\lambda }(x),t\in {\mathbb{R}}^{+}.$$

$${\mathcal{K}}^{t}g={\mathcal{K}}^t\sum _k{c}_k{\varphi }_{\lambda ,k}=\sum _k{c}_k{\mathcal{K}}^t{\varphi }_{\lambda ,k}=\sum _k{c}_k{\lambda }^t{\varphi }_{\lambda ,k}.$$

$$A_{\mathcal{K}}g:=\underset{t\to 0}{\lim }\frac{{\mathcal{K}}^{t}g-g}{t}$$

$${⟨\nabla {\varphi }_{\lambda }(x),v(x)⟩}_{{\mathbb{R}}^n}=\omega {\varphi }_{\lambda }(x)$$

$$v(x)=\sum _i^d{v}_i(x){\partial }_i.$$

$${\mathcal{K}}^t=\exp \left(t{A}_{\mathcal{K}}\right).$$

$$\mathbb{P}[x_{t+\tau }\in A|x_t=x]={\int }_{A}p(\tau ,x,y)dy$$

$$p(t,x,y)={\delta }_{\varphi (t,x)}(y).$$

$$K^{\Delta t}f(t,x)={\int }_{X}p(\Delta t,x,y)f(t,y)d\mu (y)=\mathbb{E}\left[f(t,x_{t+\Delta t})|x_t=x\right],$$

$$P^{\Delta t}q(t,x)={\int }_{X}p(\Delta t,y,x)q(t,y)d\mu (y).$$

$${⟨K^{\Delta t}f,g⟩}_{\mathcal{F}}={⟨f,P^{\Delta t}g⟩}_{\mathcal{F}}.$$

$$p(t,x,y)=p(t,y,x)$$

$$\text{d}{\text{x}}_t=-D\nabla V({\text{x}}_t)\text{d}t+\sqrt{2dD}\text{d}{\text{W}}_t,$$

$${\frac{d}{dt}{S}_t(x)|}_{t=0}=-ax,x\in \mathbb{R},t\ge 0,$$

$$x_{n+1}=x_n+\Delta t{\frac{d}{dt}{S}_t(x_n)|}_{t=0}.$$

$$x_{n+1}=x_n-a\Delta t{x}_n=(1-a\Delta t)x_n.$$

$$[{\mathcal{K}}^{n}g](x)=g((1-a\Delta t{)}^{n}x)=(1-a\Delta t{)}^{n}g(x),$$

$$\underset{n\to \infty }{\lim }{\left(1-\frac{a}{n}\right)}^n=\exp (-a).$$

$$G=D(X)=\left[d_1(X),d_2(X),\dots ,d_{N_D}(X)\right]\in {\mathbb{R}}^{N_X\times N_D}.$$

$$[{\mathcal{K}}^t{d}_k](x)=(d_k\circ S_t)(x),$$

Full text

{AMS}

47B33, 68W40, 37C10

1 Introduction

The relation between algorithms and dynamical systems has been studied for a long time. Due to the breadth of the topic, we cannot give an exhaustive review of the literature and will only briefly discuss some of the research areas important to this work. The book by Stuart and Humphries [42] and the article by Chu [8] provide a broad overview of the relation between algorithms and dynamical systems. Much recent research focuses on the discrete steps in the solution to linear systems and the relation to continuous dynamics on manifolds of matrices [8, 16, 25, 33]. In particular, algorithms for eigenvalue problems have been studied in this context, with the relation to integrable systems (Lax-pair formulation) [44]. Studying the connection between algorithms and dynamical systems led to many surprising results, such as the definition of systems that can “sort numbers” [4], or, more generally, solve optimization problems with constraints [15]. Chaotic dynamical systems that “paint”, i.e., have stationary distributions according to the distribution of paint in an image, have been studied in relation to Markov Chain Monte Carlo algorithms [39]. Another connection relates particle dynamics to algorithmic solutions [11]. In the neural network research community, the idea of a (deep) neural network representing an iterated map of a discrete dynamical system has gained attention, particularly with respect to the continuous generator of this discrete map [36]. In general, optimizing the same loss function through different algorithms such as gradient descent or “stochastic” gradient descent may result in different trajectories, and different local or even global minima may be found. Optimization algorithms modeling second-order dynamical systems (for example, Nesterov’s method, sometimes called the “heavy ball with friction” algorithm) are favoured over traditional gradient descent methods, because of their ability to overcome local minima [2, 35].

The work of Koopman and von Neumann on ergodic, mixing, and chaotic dynamical systems [20, 45] revealed that that there is a canonical, linear operator associated to each system. The operator acts on complex-valued functions of the system state, the result being evaluation of the function at a future time. Linearity of the operator makes it amenable to finite-dimensional matrix approximations on computers, which was heavily exploited in the last twenty years [29, 7]. The Koopman operator and its adjoint, the Frobenius-Perron operator, are also known to provide optimal basis functions for uncertainty quantification [1, 14].

In this paper, we discuss benefits of the operator viewpoint on algorithms, when considering them as dynamical systems. The main problem we are addressing with the new operator viewpoint is the data-driven analysis of complex algorithms acting on high-dimensional state spaces. The Koopman operator and related numerical approximation methods for it allow us to treat this challenge in a unified way, solving many pressing issues such as acceleration and formulation of data-driven surrogate models, discovery of (almost) invariant sets and regions of convergence, and high-dimensional state spaces through model order reduction. We demonstrate that a numerical approximation of the operator may be possible, even for systems with partially available or high-dimensional data. We also show how properties of the Koopman operator, such as ergodic quotients and spectral analysis, offer valuable insight into algorithmic behavior. Even if the only available data is a set of randomly distributed, single-iteration pairs of states, the operator (and thus a sense of the algorithm) can be approximately constructed. We demonstrate a particularly novel approach to algorithm analysis through the computation of the spectrum of the Koopman operator. Summarizing a complex, nonlinear algorithm in such a meaningful way–and in the language of linear algebra–may provide a missing link to more fundamental results. An example is the comparison of algorithms based on their Koopman spectrum, which may lead to a definition of metrics or distances between algorithms. Conjugacy arguments have also been used to analyze dynamical systems for a long time. The operator viewpoint on algorithms can bring this successful approach to algorithm analysis, comparison, and development, even in a data-driven setting [3].

The remainder of the paper is organized as follows. In section 2.1, we introduce the relation between algorithms and dynamical systems. In section 2.2, we describe the Koopman operator framework, followed by section 3 with the Extended Dynamic Mode Decomposition algorithm for numerical approximation of the operator. Acceleration and domain decomposition of algorithms is outlined in section 4.1. In section 4.2, we describe how to construct a continuous version of the iterative Nesterov algorithm in a data-driven way. Section 4.3 describes numerical approximations of the Koopman operator for algorithms that act on high-dimensional spaces, but move their state close to a low-dimensional manifold after a few iterations. We analyze the Newton-Raphson method for root finding of polynomials on the complex plane in section 4.4, with explicit construction of the spectrum and eigenfunctions, as well as a numerical analysis of chaotic behavior of the method on the real line. Appendix A contains some preliminary computational experiments on the impact (on the numerical construction of the Koopman operator) of only partial information about the state space (due to finite sample size). We conclude with a summary of the results and an outlook on data-driven algorithm analysis in section 5.

2 The Koopman operator of algorithms

The focus of this paper is on potential benefits of the data-driven study of algorithms as dynamical systems from the Koopman operator viewpoint. In this section, we describe the necessary mathematical framework for dynamical systems, state the algorithms that we will use as examples, and briefly introduce the Koopman operator.

2.1 Algorithms are dynamical systems

In this paper, we consider algorithms that act on their state either continuously or in discrete steps (iterations). We assume that the algorithms’ state spaces $X\subseteq\mathbb{R}^{d}$ are smooth, $k$-dimensional Riemannian submanifolds embedded in Euclidean space of dimension $d\in\mathbb{N}$, $d\geq k$. The Riemannian metric on $X$ is induced by the embedding. We define iterative algorithms as differentiable maps $a:X\to X$, with iteration number $n\in\mathbb{N}$. For a single iteration of the map $a$ on a state $x_{n}$, we write

[TABLE]

For some algorithms, it is useful to consider continuously evolving state variables. In this case, the algorithm is typically represented as a differentiable vector field $v:X\to\mathbb{R}^{k}$. Starting from a given initial state $x_{0}\in X$, this vector field generates a curve $c:\mathbb{R}^{+}\to X$ in the state space through

[TABLE]

Equation (1) describes the action of a continuous algorithm. It can be reformulated as a discrete map through the definition of a flow: The map ${S}:\mathbb{R}^{+}\times X\to X$ is called a (semi-)flow of the vector field $v$ if, for all $t,s\in\mathbb{R}^{+}$ and all $x\in X$,

[TABLE]

If we fix a time interval $\Delta t\in\mathbb{R}^{+}$, the map ${S}(\Delta t,\cdot)=:{S}_{\Delta t}(\cdot):X\to X$ is an iterative algorithm (a discrete-time dynamical system). This construction is always possible: however, it is important to note that in general it is not possible to find a continuous algorithm $v$ for a given iterative algorithm $a$. We will discuss the necessary properties of $a$ for a continuous version to exist, and show how it can be approximated in a data-driven procedure through the consideration of the Koopman operator. In this paper, we study the prototypical forms of algorithms listed in Table 1, formulated in continuous and discrete time.

2.2 Introduction to the Koopman operator

We start with a definition of the family of Koopman operators in the deterministic setting in section 2.2.1. Details about the operator in the stochastic setting, and the relation with its adjoint, the Frobenius–Perron operator, are discussed in section 2.2.2.

2.2.1 Deterministic setting

Given a deterministic dynamical system with vector field $v:X\to\mathbb{R}^{k}$, flow ${S}_{t}:X\to X$, initial condition $x_{0}\in X$, and

[TABLE]

the family of Koopman operators $\mathcal{K}^{t}$ indexed by $t\in\mathbb{R}$ is a family of linear operators acting on the function space $\mathcal{F}$ of observables $g:X\to\mathbb{C}$ such that

[TABLE]

The choice of function space $\mathcal{F}$ crucially determines the properties of $\mathcal{K}^{t}$. A typical choice is the space of complex-valued functions on the domain $X$ that are square-integrable with respect to a measure $\mu$, denoted

[TABLE]

The space $L^{2}(X,\mathbb{C},\mu)$ is typically equipped with an inner product

[TABLE]

where the bar over $g_{2}(x)$ denotes complex conjugation. The set of all $\mathcal{K}^{t}$ forms a continuous group w.r.t the strong operator topology. It is a semi-group if ${S}_{t}$ is not invertible for all $t$ (e.g., if the system approaches spatial infinity in finite time). The strong operator topology is defined through pointwise convergence on $(\mathcal{F},\|\cdot\|_{\mathcal{F}})$ [13]. The map $t\mapsto\mathcal{K}^{t}$ is continuous in the strong operator topology on a compact subset $B\subset\mathbb{R}$ if, for $t\in B$,

[TABLE]

is continuous for every $f\in\mathcal{F}$. Continuous (semi-)groups are then defined as follows. The family $\{\mathcal{K}^{t}\}_{t\in\mathbb{R}}$ is a continuous (semi-)group of operators if, in the strong operator topology,

[TABLE]

Typically, a fixed $t=t_{0}$ is chosen, and the associated Koopman operator $\mathcal{K}^{t_{0}}$ is studied. In the remaining text, if the specific choice of time is not important in the context, we will choose $t=1$ and drop the superscript $t$ to simplify notation. An important feature of the operator is the spectrum $\sigma(\mathcal{K})\subset\mathbb{C}$, where we consider the atomic (or pure-point) part $\sigma_{pp}$ and the absolutely continuous part $\sigma_{ac}$ (see [7, 31] for discussions of the singularly continuous part). We can use the spectral theory for linear operators [34, 29] to write

[TABLE]

where $P_{\lambda}$ and $E(\lambda)$ are projection operators mapping functions to their projection in the eigenspaces associated to $\lambda$. Equation (4) decomposes the operator $\mathcal{K}^{t}$ at $t=1$. Following functional calculus, the spectrum associated to other values of $t$ can be obtained by exponentiation of $\lambda$, i.e. if $\lambda\in\sigma(\mathcal{K})$ then $\lambda^{t}\in\sigma(\mathcal{K}^{t})$. For an eigenvalue $\lambda\in\sigma_{pp}(\mathcal{K})$, there are functions $\phi_{\lambda}\in\mathcal{F}$ such that

[TABLE]

These functions are called eigenfunctions associated to the eigenvalue $\lambda$. To predict the time evolution of a function $g\in\text{span}{\{\phi_{\lambda,k}\}}\subset\mathcal{F}$ such that $g=\sum_{k}c_{k}\phi_{\lambda,k},\ c_{k}\in\mathbb{C}$, we can thus write

[TABLE]

The coefficients $c_{k}$ are sometimes called “Koopman modes” for the observable $g$, since they are used to reconstruct it. For many dynamical systems, the span of the eigenfunctions is a rich subspace of $\mathcal{F}$, and many observables can be defined through their Koopman modes.

For continuous-time dynamical systems, another important object related to the Koopman operator group $\mathcal{K}^{t}$ is the infinitesimal generator $A_{\mathcal{K}}$, where

[TABLE]

for all $g\in D(A_{\mathcal{K}})\subset\mathcal{F}$. For continuous semi-groups, the domain $D(A_{\mathcal{K}})$ is dense in $\mathcal{F}$ [13]. The infinitesimal generator has the same eigenfunctions as any $\mathcal{K}^{t}$. For fixed $t>0$ and eigenvalues $\lambda^{t}$ of $\mathcal{K}^{t}$, the eigenvalues $\omega$ of $A_{\mathcal{K}}$ satisfy $\exp(t\omega)=\lambda^{t}$. The generator is important because it induces the following relation on the eigenfunctions $\phi_{\lambda}$:

[TABLE]

for all $x\in\mathbb{R}^{n}$ where the eigenfunctions and the vector field are defined (see [3] for a discussion). Note that equations (6) and (7) reveal the relation between the generator $A_{\mathcal{K}}$ and the vector field $v$. If we interpret $v$ as a map from a point $x\in X$ to the coefficients $v_{i}$ for basis vectors $\partial_{i}$ of the tangent space $T_{x}X$, then we can write

[TABLE]

With this reformulation, the vector field is an object that acts on functions $g:X\to\mathbb{R}$ from the left, i.e. $vg=\sum_{i}^{d}v_{i}\partial_{i}g$, which is exactly what equation (6) describes. Also, formally,

[TABLE]

In this sense, the vector field $v$ acts as the generator of the Koopman operator (semi)group. The generator $A_{\mathcal{K}}$ can be used to obtain an approximation of the vector field $v$ from only discrete-time samplings of the dynamical system. This was used to identify vector fields from data of dynamical systems [26], and we will use it to identify vector fields for discrete algorithms in section 4.2.

The identification of vector fields using the Koopman operator is a global approach, taking into account all the data to obtain the vector field at every point. A different, more localized idea of identifying continuous vector fields on discrete data (e.g., [38, 37, 17]) is becoming increasingly popular in the machine learning literature.

2.2.2 Stochastic setting

We provide a brief introduction to the operator definitions in the stochastic setting. The manuscript mostly focuses on deterministic algorithms, and thus we refer to [41, 9, 19] for more theoretical details and applications of the stochastic case. Note that some numerical algorithms used to approximate the Koopman operator in the deterministic setting can be analyzed using a stochastic approximation of the system that accounts for numerical errors [46, 9].

Different from the deterministic systems based on equation (2), in this section, we start with a time-homogeneous stochastic process $\{x_{t}\}_{t\in\mathbb{R}^{+}_{0}}$ defined on the space $X$ with a probability measure $\mathbb{P}$. In this setting, the formulation of the Koopman operator and its adjoint, the Frobenius-Perron operator, can be stated in terms of the transition density function [19].

Definition 2.1**.**

The transition density function $p:\mathbb{R}^{+}_{0}\times X\times X\to\mathbb{R}^{+}_{0}$ of a process $\{x_{t}\}_{t\in\mathbb{R}^{+}_{0}}$ is defined by

[TABLE]

*for all measurable sets $A\subseteq X$. *

Note that for a deterministic system with flow $\phi$, the transition density function collapses to a Dirac delta centered on the point $\phi(t,x)$:

[TABLE]

For $1\leq p\leq\infty$, the spaces $L^{p}(X,\mu)$ denote the $p$-integrable functions w.r.t. $\mu$ (see equation (3) for $p=2$). The Koopman operator and its adjoint can now be defined as follows.

Definition 2.2**.**

Let $q(t,\cdot)\in L^{1}(X,\mu)$ be a probability density and $f(t,\cdot)\in L^{\infty}(X,\mu)$ an observable of the system with flow $\phi$. For a given time step $\Delta t\in\mathbb{R}^{+}$, the Koopman operator $K^{\Delta t}:L^{\infty}(X,\mu)\to L^{\infty}(X,\mu)$ is defined by

[TABLE]

and the Frobenius–Perron operator $P^{\Delta t}:L^{1}(X,\mu)\to L^{1}(X,\mu)$ is defined by

[TABLE]

The Koopman and Frobenius–Perron operators are adjoint w.r.t the inner product on a function space $\mathcal{F}$ if for all $f,g\in\mathcal{F}$,

[TABLE]

If for a given stochastic dynamical system the transition function satisfies the detailed balance condition

[TABLE]

for all $t\in\mathbb{R}^{+}$ and all states $x,y\in X$, the Koopman operator is self-adjoint (and hence equal to the Frobenius–Perron operator) on $L^{1}(X)\cap L^{\infty}(X)$, which is a direct consequence of the definitions of the operators (equations (10–11)).

The Frobenius–Perron and Koopman operators are the solution operators of the forward (Fokker–Planck) and backward Kolmogorov equations, respectively [22, section 11]. That is why they are also often referred to as forward and backward (transfer) operators.

General Smoluchowski equations of a $d$-dimensional system are given by

[TABLE]

with potential $V:X\to\mathbb{R}$, diffusion coefficient $D$, and Wiener process $W_{t}$. They are time-reversible [28, 23], and so their transition density function satisfies the detailed balance condition. Therefore, the two transfer operators are identical (self-adjoint). The fact that this holds for equations of type (13) is extremely strong and useful, particularly for the analysis of molecular dynamics simulators. As many algorithms can be reformulated in this setting, too, these strong results can be carried over directly.

2.3 Illustrative example: Koopman operator of the Euler method

We use a typical example from the stability analysis of algorithms to familiarize the reader with the concept of Koopman operators of algorithms. Consider the following ordinary differential equation (ODE) with parameter $a>0$,

[TABLE]

and solution ${S}_{t}(x)=\exp(-at)x$. A trajectory ${S}_{t}(x_{0})$ can be approximated with a numerical algorithm, an initial value solver, given only the ODE (14) and an initial condition $x_{0}\in\mathbb{R}$. Here, we choose the forward Euler method, compute its Koopman operator, and show how stability analysis of the method is related to the spectrum of the operator. For a given step size $\Delta t>0$, the forward Euler method is given through the following iterative scheme, starting at $x_{0}$ for $n=0$:

[TABLE]

For the ODE (14), we thus have the linear, discrete-time dynamical system

[TABLE]

We now choose the function space $\mathcal{F}=\text{span}\{x^{k},x\in\mathbb{R},k\in\mathbb{N}\}$, with the identity function $g(x)=x$ as the generator of the basis. The Koopman operator $\mathcal{K}^{n}:\mathcal{F}\to\mathcal{F}$ applied to the generator is

[TABLE]

which shows that $g$ is an eigenfunction of $\mathcal{K}^{n}$ associated to the eigenvalue $(1-a\Delta t)$. On the space $\mathcal{F}$, the spectrum of $\mathcal{K}$ consists of the eigenvalues $\omega_{k}=(1-a\Delta t)^{k}$, $k\in\mathbb{N}$, associated to the eigenfunctions $\phi_{k}(x)=g(x)^{k}=x^{k}$. It is easy to see that for $a\Delta t\in[0,2]$, $\mathcal{K}^{n}$ is a contraction operator on $\mathcal{F}$: for all $k$, $|\lambda_{k}|\leq 1$. This is a sufficient condition for the numerical stability of the Euler method in this region. Note that system (14) has eigenvalues $\lambda_{k}=\exp(-ka)$ on $\mathcal{F}$, and so the discrete system (16) defined by the Euler method is not conjugate to the continuous system for any $\Delta t>0$ (for conjugacy, the eigenvalues would have to be the same). Consistency–and thus, convergence–of the algorithm can now also be interpreted in the spectral sense, by setting $n:=\lceil{1/\Delta t}\rceil$ (rounding to the next largest integer) and applying the formula

[TABLE]

After this pedagogical example, we show how Koopman operators can provide a useful framework to understand more general algorithms. In a more applied setting, “the algorithm” is usually a large, complex piece of software that is treated as a black box data generator (for example, an optimization procedure for a supply chain). Then, numerical approximations of the Koopman operator can be used. In the next section, we briefly describe Extended Dynamic Mode Decomposition as an example for such an approximation procedure from data.

3 Data-driven approximation of the Koopman operator

We briefly cover the approximation methods employed in the computational experiments of later sections. The methods include an approximation of the Koopman operator with pure-point spectrum (Extended Dynamic Mode Decomposition), extracting the infinitesimal generator of the family of Koopman operators, as well as an approximation algorithm in case the Koopman operator is unitary and has continuous spectrum.

3.1 Extended Dynamic Mode Decomposition

One extension of Dynamic Mode Decomposition (DMD) [40] to nonlinear systems with pure-point spectrum is Extended DMD [46, 48, 47] (EDMD). Many more numerical methods exist, such as Generalized Laplace Averaging [30, 27], EDMD with dictionary learning through neural networks [24], dictionary selection through $L^{1}$-optimization [5], etc.

We consider the state space of a given system–usually, a smooth manifold $M$ embedded in Euclidean space $\mathbb{R}^{d}$–sampled with a finite number of points $X=\{x_{k}\in M\subset\mathbb{R}^{d}|k=1,\dots,N_{X}\}$, $N_{X}\in\mathbb{N}$. The main idea of (Extended) Dynamic Mode Decomposition is to approximate the action of the operator on a function space $\mathcal{F}$ over $M$ by choosing a finite-dimensional subspace spanned by a finite set of real-valued functions of $M$ (called a “dictionary”, with “observables” as elements). For $\mathcal{F}$, we typically use the space $L^{2}(M,\mathbb{C},\mu)$, see equation (3). The EDMD algorithm computes the action of the Koopman operator on the dictionary at the points in $X$, and then approximates the Koopman operator by solving a least-squares problem:

Given a dictionary $D=\{d_{k}:M\to\mathbb{R}|k=1,\dots,N_{D}\}\subset\mathcal{F}$ with $N_{D}$ observables, construct

[TABLE] 2. 2.

For a fixed $t>0$, compute the action of the Koopman operator $\mathcal{K}^{t}$ on the dictionary elements using the flow map ${S}_{t}:M\to M$ of the given system:

[TABLE]

and construct the matrix

[TABLE] 3. 3.

Approximate the operator $\mathcal{K}^{t}$ through the matrix

[TABLE]

where $(G^{T}G)^{\dagger}$ denotes the pseudo-inverse of the matrix $(G^{T}G)$.

The choice of the dictionary $D$ is crucial for a successful approximation of the operator, see [24, 46, 47] and references therein. In this paper, we employ thin-plate radial basis functions [10],

[TABLE]

where the point $x$ is the center of the radial basis function and $\delta$ is a small, positive constant to extend the function to points $y=x$ (here: $\delta=10^{-3}$). The centers $x$ are uniformly distributed over the data set $X$, where the uniform distribution is constructed through an appropriate k-means algorithm, clustering the data in $N_{D}$ clusters and using the centers of the clusters as the centers for the basis functions.

3.2 Approximation of the infinitesimal generator

After some of the eigenfunctions $\hat{\phi}$ and eigenvalues $\lambda$ of $\mathcal{K}$ are available through their numerical approximations, we can represent the approximation of the infinitesimal generator $A_{\mathcal{K}}$ for the group $\mathcal{K}^{t}$ as a matrix $L=V\frac{\ln[\Lambda]}{\Delta t}V^{-1}$ (see [26]), where $V$ contain the eigenvectors of the matrix representation of a particular $\mathcal{K}^{\Delta t}$, obtained with EDMD from snapshot data $(x,{S}_{\Delta t}(x))$. For the coordinate functions $x_{i}\in\mathcal{F}$ and a finite step size $\Delta t$,

[TABLE]

where $\Lambda$ is the diagonal matrix of eigenvalues of $\mathcal{K}$, $\ln$ is the complex logarithm (where we pick the principal branch), and $C$ are the “Koopman modes” associated to the coordinate functions $x_{i}:X\to\mathbb{R}$. The approximation of the vector field $v$ generating the flow $S_{t}$ is thus

[TABLE]

Note that the approximation of the generator through the logarithm of the eigenvalues involves several challenges: numerical issues with eigenvalues of $\Lambda$ close to zero, and non-invertibility (or rather, non-uniqueness) of the complex logarithm. The issue with non-uniqueness is noted by Mauroy and Goncalves [26]. In the example where we construct $L$ (section 4.2), the spectrum of the Koopman operator lies in the unit disk, so the complex logarithm is unique when picking its principal branch. To mitigate numerical issues, we set all eigenvalues with $\text{Re}{(\ln\lambda)}<-2/\Delta t$ to zero.

3.3 Approximation of the continuous spectrum

The EDMD algorithm is applicable for approximating the Koopman operator if its spectrum only consists of eigenvalues (pure-point spectrum). Not many numerical approximations are available for operators with continuous spectra. Here, we briefly describe the results from Korda, Putinar and Mezić [21] for unitary Koopman operators, where the spectrum is concentrated on the unit circle $\mathbb{T}$ in the complex plane. We will use this approximation to analyze the Newton algorithm in section 4.4. If the Koopman operator of a system is unitary, we can write it as an integral over a projection operator-valued measure $E$ on the unit circle $\mathbb{T}\subset\mathbb{C}$,

[TABLE]

For any observation function ${g}$ in the domain $\mathcal{F}$ of $\mathcal{K}$, the measure $E$ defines a real-valued, positive measure $\mu_{g}$ on $\mathbb{T}$ through

[TABLE]

for all Borel sets $A\subset\mathbb{T}$. The moments $m_{k}$, $k\in\mathbb{Z}$ of $\mu_{g}$ are defined by

[TABLE]

where the last identity follows from the spectral theorem. The measure $\mu_{g}$ depends on the choice of observable $g$, but under certain conditions on $g$, $\mu_{g}$ fully determines the operator $\mathcal{K}$. The conditions are satisfied if $g$ is $\star$-cyclic, meaning that repeated applications of $\mathcal{K}$ to ${g}$ span the function space $\mathcal{F}$, see [21]. In this case, approximating the measure $\mu_{g}$ not only reveals information about the specific observable $g$, but also about the underlying system (and its operator $\mathcal{K}$). The density of $\mu_{g}$ on the unit circle can be visualized after a numerical approximation of the moments. If we have $N$ observations $y_{i}=g(x_{i})$ along a trajectory of an ergodic system, we can estimate the moments through

[TABLE]

Define $\psi_{N}(z)=[1,z,z^{2},\dots,z^{N}]^{T}$, and write the Christoffel–Darboux kernel as

[TABLE]

with a semi-positive definite, Hermitian Toeplitz matrix ${M}$ defined through

[TABLE]

Inverting the matrix $M$ is non-trivial, see [21] for details. Instead of inverting $M$ directly, the authors propose to invert a matrix $\tilde{M}$ where the elements $m_{k}$ in $M$ are replaced by elements

[TABLE]

Given the kernel $\tilde{K}_{N}$ (now using $\tilde{M}$ instead of $M$), the continuous density $\rho$ of $\mu_{g}$ at a given angle $\theta\in\mathbb{T}\cong[0,2\pi)$ can be approximated through

[TABLE]

4 Analysis of algorithms through their Koopman operator

This section contains illustrative examples of algorithm analysis through their Koopman operator. We demonstrate several aspects of our solution to the main challenge: how to analyze complex, nonlinear algorithms in a data-driven setting within a unified framework. In particular, we show

how to accelerate algorithms by constructing data-driven surrogate models, 2. 2.

approximations of infinitesimal generators for discrete-time algorithms, 3. 3.

approximations of the (asymptotic, long-time) operator in high-dimensional state spaces, 4. 4.

and a combination of an explicit operator analysis for the Newton-Raphson method with a data-driven approximation of its continuous spectrum.

4.1 Acceleration

We can employ the predictive capabilities of the Koopman operator framework to accelerate the application of algorithms. First, the Koopman operator eigenvalues, eigenfunctions, and modes are constructed in an offline phase, using EDMD (see section 3). Second, given a new initial condition close to the ones in the data set used for constructing the operator, we can use the approximation as a data-driven surrogate model for the original algorithm.

Such a surrogate model may be useful in accelerating the training of neural networks. There, the state space is often very high-dimensional, because each state consists of all weights and biases of the network. The typical optimization algorithm, stochastic gradient descent, employed to change the state in order to minimize a certain loss function, is based on the gradient descent algorithm we study in this section. The stochastic version can be analyzed similarly, using the setting described in section 2.2.2. Different types of training (e.g. through mini-batches, or dropout) will induce different biases, and the Koopman operator based surrogate may provide a path towards quantifying and comparing such inductive biases.

Figure 2 shows the setting used for the example in this section, where we study the Gradient Descent algorithm minimizing Himmelblau’s function [18]—a standard test function for optimization algorithms:

[TABLE]

This function has one local maximum and four local minima with three connecting saddle points (see figure 1, a). As a dictionary for EDMD, we employ 500 thin-plate radial basis functions (see equation (17)), chosen uniformly distributed on $M=[-4,4]^{2}$. In addition to the radial basis functions, we add the two coordinate functions $x_{1},x_{2}$ as well as the constant function with value $1$ to the dictionary.

Figure 2 shows four predicted trajectories using the approximated Koopman operator. To generate a trajectory, we follow the prediction steps outlined in section 2.2 (equation (5)), using all 503 eigenvalues and eigenfunctions.

The eigenfunctions of $\mathcal{K}^{t}$ at eigenvalue $\lambda=1$ are special, as they are preserved under the flow of the system—they do not decay, expand, or oscillate. This can be used to construct an ergodic decomposition [6, 7] of the state space, separating the basins of attraction. In the system discussed here, we expect to see four basins of attraction, corresponding to four attracting fixed points, and accordingly, at least four eigenfunctions associated to eigenvalue $\lambda=1$. To approximate the ergodic decomposition, we consider the map $M\ni p\mapsto[\phi_{1}(p),\dots,\phi_{4}(p)]^{T}\in\mathbb{C}^{4}$, i.e. from the states $M$ into the values of eigenfunctions associated to $\lambda=1$. In this space (shown projected to two dimensions in figure 2, b), there should be four clusters of values corresponding to the different values the eigenfunctions take on over the corresponding basins of attraction. Clustering the data using the k-means algorithm then allows us to distinguish the original points in $M$ (plot c in figure 2).

4.2 Continuous versions of discrete algorithms

In recent years, extensive research has focused on re-deriving established algorithms as discrete versions of continuous dynamical systems, in this way “explaining” the original, iterative algorithm. In the case of the Nesterov method, the appropriate continuous system derived by Candès et al. [43] is a time-dependent, second-order ordinary differential equation. This ODE can be (trivially) transformed to a first-order system of autonomous equations, given through

[TABLE]

where $\dot{x}$ now refers to a derivative with respect to a new time variable. In the example we discuss below, the function to minimize is $f(x)=\frac{1}{2}x^{2}$, $x\in\mathbb{R}$. As in the previous section, we use EDMD to approximate the eigensystem of the Koopman operator with 125 thin-plate radial basis functions (see equation (17)) as well as the $x,v,t$ coordinates and the constant function as a basis for the function space of observables on $\mathbb{R}^{3}$ ($(x,v,t)$-space). The operator approximation is performed with data $(z_{n},z_{n+1})$ from the discrete algorithm with step size $h=0.01$ (see table 1), where $z_{n}:=(x(n),(x(n)-x(n-1))/h,t(n))$. We then approximate the corresponding vector field on the three coordinates of $z$ through equations (18) and (19), which should be similar to the vector field in equation (20). Figure 3 shows streamlines of the vector field (20) for three different values of $t$, with a distinguished sample trajectory, compared to the vector field approximated by using the Koopman operator. This example illustrates that it is possible to extract a continuous version (an infinitesimal generator) of a discrete algorithm.

4.3 High-dimensional state spaces

Many algorithms are operating on points in high-dimensional state spaces that are parametrized through a large number of variables. In this case, sampling the full state space to obtain an accurate approximation of the Koopman operator for the algorithm is difficult, especially when the action of the algorithm on the state is nonlinear. However, if a few iterations of the algorithm quickly contract the whole state space towards a low-dimensional subset, then numerical approximations can be successful in reducing the overall computational effort. The EDMD algorithm does not need an explicit expression for this low-dimensional subset if, for example, radial basis functions are used for the dictionary.

As an example, we construct a high-dimensional optimization problem with a low-dimensional, attracting manifold, where gradient descent is used as an optimizer. On the manifold, starting from the Müller-Brown potential, we add quadratic terms to obtain asymptotically stable, attractive behavior in the directions transverse to the manifold. On the manifold (here, a two-dimensional subspace of $\mathbb{R}^{100}$ parametrized by $x,y$), the potential is defined by

[TABLE]

with

[TABLE]

The cost function $C:\mathbb{R}^{100}\to\mathbb{R}$ on the embedding space is then defined through

[TABLE]

where we rotate the data by a random, unitary matrix $U\in\mathbb{R}^{100\times 100}$ into the 100 dimensional space to demonstrate that the EDMD algorithm with radial basis functions does not depend on the chosen embedding coordinate system. Now, all coordinates have non-trivial dynamics, but the low-dimensional, attractive manifold is still present. To approximate the operator, we use 625 thin-plate radial basis functions with centers randomly distributed over the initial data set in the 100-dimensional space, as well as the 100 coordinates $x_{1},\dots,x_{100}$ and the constant function (i.e., the final dictionary has 726 elements). The initial data consists of 2500 data points, sampled in a 100-dimensional standard normal distribution around zero, and then evolved forward with five iterations of the algorithm to converge towards the low-dimensional structure. Note that this number of points is–by far–not enough to densely sample the high-dimensional space, but the attracting behavior of the algorithm towards the low-dimensional structure is nonetheless sufficient to numerically approximate the operator. We use radial basis functions for EDMD to circumvent the need for an explicit approximation of the low-dimensional space: the functions only depend on the distance between the points, not the ambient dimension. As in section 4.1, we can again accurately predict trajectories of gradient descent (see figure 4).

4.4 Koopman operators for the Newton-Raphson method

In this section, we mostly use the conjugacy results from [32] to construct eigenfunctions for the Koopman operator of the Newton-Raphson method for root finding. We only consider polynomial functions of degree two in this analysis, and at the end give a hint of more complicated dynamics present for cubic polynomials. The Newton-Raphson method to find roots of a function $f$ is defined as the iterative scheme

[TABLE]

where $[Jf]_{ij}=\frac{\partial f_{i}}{\partial x_{j}}$ is the Jacobian matrix of $f$. In the example we consider here, the function $f:\mathbb{C}\to\mathbb{C}$ is a complex polynomial of degree two, for which the Newton-Raphson method simplifies to

[TABLE]

Without loss of generality, we assume $f$ has the form $f(z)={a}z^{2}+bz+d$ with $a,b,d\in\mathbb{C}$, $a\neq 0$. Any polynomial $f$ of degree two is conjugate to exactly one polynomial $g_{c}(z)=z^{2}+c$, $c\in\mathbb{C}$ (see [32]), where the conjugacy is given through a linear transformation $h(z)=az+\frac{b}{2}$, such that

[TABLE]

The constant $c$ is related to $a,b,d$ through $c=ad+\frac{b}{2}-\frac{b^{2}}{4}$. The map $h$ scales and translates the argument so that $g_{c}$ has the desired form. The conjugacy of $f$ and $g_{c}$ implies that we only need to analyze the Newton-Raphson method applied to functions of the form $g_{c}$, with a given $c$. For this simpler form of $f$, the Newton-Raphson method becomes

[TABLE]

In the Koopman operator picture, the map $N_{c}$ defines the dynamical system in discrete time with an associated Koopman operator $\mathcal{K}_{N_{c}}$. An important result in [32] is that the map $N_{c}$ is conjugate to the polynomial $g_{0}(z)=z^{2}$ through the map

[TABLE]

which is a Möbius transformation, and hence invertible on ${\mathbb{C}\cup\{\infty\}}$. Thus, eigenfunctions of $\mathcal{K}_{N_{c}}$ can be constructed through eigenfunctions of $\mathcal{K}_{g_{0}}$, which are defined through

[TABLE]

Since $g_{0}(z)=z^{2}$, there are eigenfunctions $\phi_{k}$ of the form $\phi_{k}(z)=\ln(|z|)^{k}$, associated to the eigenvalues $\lambda_{k}=2^{k}$. Through the given conjugacy $h_{0}$, we have that

[TABLE]

are eigenfunctions for $\mathcal{K}_{N_{c}}$, associated to the eigenvalues $\lambda_{k}$. Fig. 5 shows $\psi_{1}$ for $c=1$, i.e. for the polynomial $f(z)=z^{2}+1$. The roots at $\pm i$ are clearly visible as the points where the eigenfunction diverges.

4.4.1 An example with chaotic behavior

Consider the function $f(z)=z^{2}+1$ from the previous section, with its two roots at $\pm i$. Interpreting the Newton-Raphson method as a discrete dynamical system, the real line is the basin boundary separating the two basins of attraction on the complex plane. When started exactly on the real line, the Newton-Raphson method exhibits chaotic behavior, which we analyze in this section. The Newton-Raphson method applied to $f$ is defined by the map

[TABLE]

so that $z_{n+1}=N_{f}(z_{n})$ yields a new iterate, starting at $z_{0}\in\mathbb{C}$. For $z_{0}\in\mathbb{C}/\mathbb{R}$, the limit points of the iterative scheme are the extremal points, where

[TABLE]

For $z_{0}\in\mathbb{R}$, however, the scheme does not converge (see Fig. 6) and gives rise to the chaotic recurrence (which can be solved explicitly in this case)

[TABLE]

A similar discrete system with explicit solution is given by the logistic map at parameter $r=4$, with recurrence relation $y_{n+1}=ry_{n}(1-y_{n})$ and solution $y_{n}=\sin(c_{0}\,2^{n})^{2}$, $c_{0}=\arcsin\left(\sqrt{y_{0}}\right)$. Accordingly, except for the constant function, the continuous eigenfunctions of the Koopman operator computed in the previous section are all zero (or infinity, for $k<0$) on the real line, see equation (21). For the map $N_{f}$, the real line is an invariant, repelling limit set of measure zero.

The results so far have been constructed explicitly. Addressing the main challenge of the manuscript, we now want to demonstrate that numerical algorithms can also be used to analyze the behavior of the Newton-Raphson algorithm. To this end, we employ the numerical method by Korda, Putinar and Mezić [21] (as described in section 3.3) to approximate the continuous spectrum of the Koopman operator associated to the map $N_{f}$. Figure 7 shows how different numbers of points (i.e. lengths of trajectories) approximate the measure $\mu_{g}$ for $g(z)=1+\exp(2\pi iz)$, with $z\in\mathbb{R}$ from a trajectory of the Newton-Raphson method starting at $z_{0}=1/2$. The point spectrum is $\sigma_{pp}=\{1\}$ with eigenvalue $1$ associated to the constant function, as expected for a chaotic system. The absolutely continuous spectrum appears to be $\sigma_{c}=\mathbb{T}/\{1\}$, with no singularly continuous part. The trajectories were computed with the iteration of $N_{f}$ (equation (22)) with 16 digits of accuracy, and with the explicit formula (equation (23)) with 35,000 digits of accuracy, using the Wolfram Mathematica 11 software. The high accuracy is necessary to exactly represent the number $2^{100,000}\approx 10^{31,000}$ for the last point in the longest trajectory.

4.4.2 Fractal eigenfunctions

For cubic polynomial functions, application of the Newton-Raphson method is known to exhibit an even richer structure compared to the quadratic polynomials described in the previous section. We only give numerical results, inspired by [32]. Figure 8 shows the number of iterations (in color) of the Newton-Raphson method until the derivative of the objective function is smaller than $0.01$. Figure 9 shows eight iterations of the Newton-Raphson method on a subset of points sampled on the complex plane, for a polynomial of degree two (left three columns), and degree three (right three columns). The convergence to the two (and three) roots is visible in the third and sixth column. The first/second and fourth/fifth columns show the real/imaginary parts of the final positions as functions on the initial positions, which makes the fractal structure for degree three (and higher) apparent.

5 Conclusions

Many modern algorithms (e.g. those associated with training neural networks) can difficult to analyze in a classical sense: they are complex, interconnected, can be deterministic or stochastic, act on high-dimensional states, and are often only accessible as input-output blax box systems. To address this challenge, we propose to employ the unifying view of the Koopman operator framework, with several numerical approximations readily available. In a series of examples, we demonstrated that the Koopman operator provides such a unifying view of many important issues in the analysis of algorithms: (spectral) convergence, algorithm acceleration, state space decomposition, high-dimensional state spaces, partial information, and generating processes for discrete algorithms. In the last section, we gave an outlook to chaotic behavior of algorithms and how the behavior of statistics of observables can be analyzed through the spectrum of the operator. We discussed the possibility of using this approach to accelerate algorithms in high-dimensional embedding spaces whose long-term dynamics lie on low-dimensional manifolds. Another possibility is to use this approach to create surrogates of complicated “black box” algorithms that are difficult to analyze mathematically; the data-driven surrogates obtained by sampling iterated algorithm states may provide useful insights in their nature and behavior. The Koopman operator has been shown to provide a convenient framework for constructing data-driven homeomorphisms between dynamical systems; it is an intriguing possibility that we can use this framework to realize homeomorphisms between different algorithms for solving the same problem.

Appendix A Computational experiments with partial information

In most applications, the data will only partially cover the state space of an algorithm. In this section, we explore the behavior of the approximated operator spectrum on such partial domains. The theoretical underpinning of this experiment is provided by the concept of “open eigenfunctions” [31], i.e. eigenfunctions defined only on parts of the state space. Not all the results from the computational experiments can be readily explained by the current theory, however, as we will detail below.

We study discrete gradient descent on the function

[TABLE]

Figure 10 illustrates the function values on $[-1,1]^{2}$ through contour lines (color), and the sample points used for the approximation of the operator (red area). When using gradient descent to minimize $f$, we obtain a (discrete) dynamical system with two attracting steady states, and one saddle point in between them (see figure 10). The spectrum obtained by EDMD for the given sample domain (red) is shown on the right: all eigenvalues are inside the unit disk, indicating that the system behavior is purely attracting. This result is reasonable, because in the limit of infinite applications of the Koopman operator, all observable functions $g\in\mathcal{F}$ of the data set can be expressed by two piecewise constant eigenfunctions $\phi_{1}^{(1)},\phi^{(2)}_{1}$ associated to the eigenvalue $\lambda=1$:

[TABLE]

where $x^{(s)}_{1}$ is the $x_{1}$ coordinate of the saddle point between the two steady states. The constants $a_{1},a_{2}$ are determined by the value of $g$ on the two attracting steady states and the four constants $c_{1}^{(1)},c_{1}^{(2)},c_{2}^{(1)},c_{2}^{(2)}\in\mathbb{C}$ associated to the two eigenfunctions. Equation (24) illustrates that, in the function space $\mathcal{F}$, there is a four-dimensional (real and imaginary parts of two complex numbers $a_{1},a_{2}$), attracting subspace that is the limit set of the system $g_{n+1}=\mathcal{K}g_{n}$.

If we shift the domain to the right (figure 11, C), the left attracting steady state leaves the sampling domain. This causes trajectories starting in the sampling domain to leave it after some time, which typically indicates a repelling set—explaining the numerically obtained eigenvalues outside the unit disk in parts A and C of figure 11. Note that since many trajectories leave the domain after a finite number $N$ of iterations (depending on the initial state $x_{n}$), the dynamical system is not defined for iteration numbers $n>N$. For points on the boundary of the domain, $N=0$, so the flow map at these points is only defined for iterations backwards in time. Since there is a saddle inside the domain, the system is also not defined for some states $x\in X$ on the boundary that would leave the domain backwards in time (iteration numbers $n<0$). This poses problems with the definition of the Koopman operator family $\{\mathcal{K}^{n}\}$, since it acts on functions defined on all of $X$—and for all $n\neq 0$, some elements leave $X$, so that the flow map on the entire set $X$ is not defined for $n\neq 0$. If all vectors on the boundary of the data domain pointed inward (outward), the flow would exist for all forward (backward) time, and the problem would be related to inflowing (overflowing) invariant manifolds [12]. Eigenvalues with a real part larger than one indicate unstable behavior, indicating the existence of unstable nodes or saddles in the data set.