{'fraction_non_alphanumeric': 0.04077174067243581, 'fraction_numerical': 0.011037026528585615, 'mean_word_length': 4.48149300155521, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Koopman operator of algorithms

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 made it possible to analyze 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 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.

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 [38] and the article by Chu [8] provide a broad overview of the relation between algorithms and dynamical systems. A lot of recent research focuses on the discrete steps in the solution to linear systems and the relation to continuous dynamics on manifolds of matrices [8,15,22,29]. In particular, algorithms for eigenvalue problems have been studied in this context, with the relation to integrable systems (Lax-pair formulation) [40]. Studying the connection between algorithms and dynamcal 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 [14]. 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 [35]. Another connection relates particle dynamics to algorithmic solutions [10]. 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 [32]. 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,31].
The work of Koopman and von Neumann on ergodic, mixing, and chaotic dynamical systems [19,41] revealed that that there is a canonical, linear operator associated to each system. The operator acts on complex-valued functions of the system state, with the result of evaluation 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 [25,7]. The Koopman operator and its adjoint, the Frobenius-Perron operator, are also known to provide optimal basis functions for uncertainty quantification [1,13].
In this paper, we discuss benefits of the operator viewpoint on algorithms, when considering them as dynamical systems. We demonstrate that a numerical approximation of the operator may be possible, even for systems with partially available or high-dimensional data. We then show how properties of the Koopman operator, such as ergodic quotients and spectral analysis, offer valuable insight into algorithmic behavior. Finally, even if the only available data is a set of randomly distributed, single-iteration pairs of states, the operator can be approximately constructed.
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. At the end of section 4, 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 algo-rithms' state spaces X ⊆ R d are smooth, k-dimensional Riemannian submanifolds embedded in Euclidean space of dimension d ∈ N, d ≥ k. The Riemannian metric on X is induced by the embedding. We define iterative algorithms as differentiable maps a : X → X, with iteration number n ∈ N. For a single iteration of the map a on a state x n , we write 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 → R k . Starting from a given initial state x 0 ∈ X, this vector field generates a curve c : R + → X in the state space through Equation (2.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 : R + ×X → X is called a (semi-)flow of the vector field v if, for all t, s ∈ R + and all x ∈ X, If we fix a time interval ∆t ∈ R + , the map S(∆t, ·) =: S ∆t (·) : X → 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.

Introduction to the Koopman operator.
Given a deterministic dynamical system with vector field v : X → R k , flow S t : X → X, initial condition x 0 ∈ X, and the family of Koopman operators K t indexed by t ∈ R is a family of linear operators acting on the function space F of observables g : X → C such that The choice of function space F crucially determines properties of 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 µ, denoted Starting from a point x0 ∈ R n , the function f : R n → R (or its absolute value, in case of root-finding algorithms) is to be minimized, where ∇f denotes its gradient, and Hxf the Hessian matrix at x. If f is vector-valued, Jxf denotes its Jacobian matrix at x. The symbol h ∈ R + denotes a small, positive constant.
Gradient Decenṫ Davidenko / Newton-Raphson method for optimization Davidenko / Newton-Raphson method for root finding Nesterov method The space L 2 (X, C, µ) is typically equipped with an inner product where the bar over g 2 (x) denotes complex conjugation. The set of all 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 (F, · F ) [12]. The map t → K t is continuous in the strong operator topology on a compact subset is continuous for every f ∈ F. Continuous (semi-)groups are then defined as follows. The family {K t } t∈R is a continuous (semi-)group of operators if, in the strong operator topology, Typically, a fixed t = t 0 is chosen, and the associated Koopman operator 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. For details about the operator in the stochastic setting, and the relation of the operator with its adjoint, the Perron-Frobenius operator, we refer to [37,18]. An important feature of the operator is the spectrum σ(K) ⊂ C, where we consider the atomic (or pure-point) part σ pp and the absolutely continuous part σ ac (see [7,27] for discussions of the singularly continuous part). We can use the spectral theory for linear operators [30,25] to write where P λ and E(λ) are projection operators mapping functions to their projection in the eigenspaces associated to λ. Equation (2.3) decomposes the operator K t at t = 1. Following functional calculus, the spectrum associated to other values of t can be obtained by exponentiation of λ, i.e. if λ ∈ σ(K) then λ t ∈ σ(K t ). For an eigenvalue λ ∈ σ pp (K), there are functions φ λ ∈ F such that These functions are called eigenfunctions associated to the eigenvalue λ. To predict the time evolution of a function g ∈ span{φ λ,k } ⊂ F such that g = k c k φ λ,k , c k ∈ C, we can thus write 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 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 K t is the infinitesimal generator A K , where The infinitesimal generator has the same eigenfunctions as any K t . For fixed t > 0 and eigenvalues λ t of K t , the eigenvalues ω of A K satisfy exp(tω) = λ t . The generator is important because it induces the following relation on the eigenfunctions φ λ : for all x ∈ R n where the eigenfunctions and the vector field are defined (see [3] for a discussion). Note that equations (2.5) and (2.6) reveal the relation between the generator A K and the vector field v. If we interpret v as a map from a point x ∈ X to the coefficients v i for basis vectors With this reformulation, the vector field is an object that acts on functions g : which is exactly what equation (2.5) describes. Also, formally, In this sense, the vector field v acts as the generator of the Koopman operator (semi)group.
The generator A 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 [23], 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., [34,33,16]) is becoming increasingly popular in the machine learning literature.

Illustrative example:
Koopman operator of the Euler method. We use a typical example in 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, and solution S t (x) = exp(−at)x. A trajectory S t (x 0 ) can be approximated with a numerical algorithm, given only the ODE (2.7) and an initial condition x 0 ∈ 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 ∆t > 0, the forward Euler method is given through the following iterative scheme, starting at x 0 for n = 0: For the ODE (2.7), we thus have the linear, discrete-time dynamical system We now choose the function space F = span{x k , x ∈ R, k ∈ N}, with the identity function g(x) = x as the generator of the basis. The Koopman operator K n : F → F applied to the generator is which shows that g is an eigenfunction of K n associated to the eigenvalue (1 − a∆t). On the space F, the spectrum of K consists of the eigenvalues ω k = (1−a∆t) k , k ∈ N, associated to the eigenfunctions φ k (x) = g(x) k = x k . It is easy to see that for a∆t ∈ [0, 2], K n is a contraction operator on F: for all k, |λ k | ≤ 1. This is a sufficient condition for the numerical stability of the Euler method in this region. Note that system (2.7) has eigenvalues λ k = exp(−ka) on F, and so the discrete system (2.9) defined by the Euler method is not conjugate to the continuous system for any ∆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 := 1/∆t (rounding to the next largest integer) and applying the formula lim n→∞ 1 − a n n = exp(−a).
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. In this section, 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.
We consider the state space of a given system-usually, a smooth manifold M embedded in Euclidean space R d -sampled with a finite number of points X = {x k ∈ M ⊂ R d |k = 1, . . . , N X }, N X ∈ N. The main idea of (Extended) Dynamic Mode Decomposition is to approximate the action of the operator on a function space F over M by choosing a finitedimensional subspace spanned by a finite set of real-valued functions of M (called a "dictionary", with "observables" as elements). For F, we typically use the space L 2 (M, C, µ), see equation (2.2). 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: 2. For a fixed t > 0, compute the action of the Koopman operator K t on the dictionary elements using the flow map S t : M → M of the given system: and construct the matrix 3. Approximate the operator K t through the matrix The choice of the dictionary D is crucial for a successfuly approximation of the operator, see [21,42,43] and references therein. In this paper, we employ thin-plate radial basis functions [9], where the point x is the center of the radial basis function and δ is a small, positive constant to extend the function to points y = x (here: δ = 10 −3 ). The centers x are uniformly distributed over the data set X, where the uniform distribution is constructed through the k-means algorithm, clustering the data in N D clusters and using the centers of the clusters as the centers for the basis functions.

Approximation of the infinitesimal generator.
After some of the eigenfunctionsφ and eigenvalues λ of K are available through their numerical approximations, we can represent the approximation of the infinitesimal generator A K for the group K t as a matrix L = V ln [Λ] ∆t V −1 (see [23]), where V contain the eigenvectors of the matrix representation of a particular K ∆t , obtained with EDMD from snapshot data (x, S ∆t (x)). For the coordinate functions x i ∈ F n and a finite step size ∆t, where Λ is the diagonal matrix of eigenvalues of 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 → R. The approximation of the vector field is thus Note that the approximation of the generator through the logarithm of the eigenvalues involves several challenges: numerical issues with eigenvalues of Λ 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 [23]. 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 Re(ln λ) < −2/∆t to zero.

Approximation of the continuous spectrum.
The EDMD algorithm is applicable to approximate the Koopman operator if its spectrum only consists of eigenvalues (pure point spectrum). Not many numerical approximations are available for operators with continuous spectrum. Here, we briefly describe the results from Korda, Putinar and Mezić [20] for unitary Koopman operators, where the spectrum is concentrated on the unit circle 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 T ⊂ C, For any observation function g in the domain F of K, the measure E defines a real-valued, positive measure µ g on T through for all Borel sets A ⊂ T. The moments m k , k ∈ Z of µ g are defined by where the last identity follows from the spectral theorem. The measure µ g depends on the choice of observable g, but under certain conditions on g, µ g fully determines the operator K. The conditions are satisfied if g is -cyclic, meaning that repeated applications of K to g span the function space F, see [20]. In this case, approximating the measure µ g not only reveals information about the specific observable g, but also about the underlying system (and its operator K). The density of µ 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 Define ψ N (z) = [1, z, z 2 , . . . , z N ] T , and write the Christoffel-Darboux kernel as with a semi-positive definite, Hermitian Toeplitz matrix M defined through Inverting the matrix M is non-trivial, see [20] for details. Instead of inverting M directly, the authors propose to invert a matrixM where the elements m k in M are replaced by elements Given the kernelK N (now usingM instead of M ), the continuous density ρ of µ g at a given angle θ ∈ T ∼ = [0, 2π) can be approximated through

Analysis of algorithms through their Koopman operator. This section contains illustrative examples of algorithm analysis through the Koopman operator.
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. This may be useful to accelerate 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. Different types of training (e.g. through mini-batches, dropout) will induce different biases, and the Koopman operator based surrogate may provide a path towards quantifying 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 [17]-a standard test function for optimization algorithms: f (x 1 , x 2 ) = (x 2 1 + x 2 − 11) 2 + (x 1 + x 2 2 − 7) 2 . 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 (3.1)), 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 (2.4)), using all 503 eigenvalues and eigenfunctions.
The eigenfunctions of K t at eigenvalue λ = 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 λ = 1. To approximate the ergodic decomposition, we consider the map

Continuous versions of discrete algorithms.
In recent years, a lot of research was focused on rederiving 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 continous system derived by Candès et al. [39] 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 whereẋ 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) = 1 2 x 2 , x ∈ 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 (3.1)) as well as the x, v, t coordinates and the constant function as a basis for the function space of observables on 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 (3.2) and (3.3), which should be similar to the vector field in equation (4.1). Figure 3 shows streamplots of the vector field (4.1) for three different values of t, with an example trajectory, compared to the vector field approximated by using the Koopman operator. The 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 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. 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, we employ the Müller-Brown potential and 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 R 100 parametrized by x, y), the potential is defined by The cost function C : R 100 → R on the embedding space is then defined through where we rotate the data by a random, unitary matrix U ∈ R 100×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 , . . . , x 100 and the constant function (i.e., the final dictionary has 726 elements). As in section 4.1, we can again accurately predict trajectories of gradient descent (see figure 4).

Koopman operators for the Newton-Raphson method.
In this section, we mostly use the conjugacy results from [28] 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  Figure 4. Top-Left: Data gathered by sampling the 100-dimensional state space with a Gaussian distribution, converging to the low-dimensional manifold in five iterations of gradient descent (end points marked in orange), and then recording one additional iteration (blue points). The coordinates x1 and x2 here are only the first two of all 100 present in the data set. Top-Right: Predicted trajectories of the algorithm (blue) that were computed after the construction of the operator, for ten initial conditions chosen at random from the initial data set (orange points). The bottom panels show predicted (blue, dots) and actual (orange, dashed) trajectories for the coordinates x1, x2, and x20, to demonstrate that all coordinate trajectories are reconstructed accurately.
for cubic polynomials. The Newton-Raphson method to find roots of a function f is defined as the iterative scheme where [Jf ] ij = ∂f i ∂x j is the Jacobian matrix of f . In the example we consider here, the function f : C → C is a complex polynomial of degree two, for which the Newton-Raphson method simplifies to Without loss of generality, we assume f has the form f (z) = az 2 + bz + d with a, b, d ∈ C, a = 0. Any polynomial f of degree two is conjugate to exactly one polynomial g c (z) = z 2 + c, c ∈ C (see [28]), where the conjugacy is given through a linear transformation h(z) = az + b 2 , such that The constant c is related to a, b, d through c = ad + b 2 − 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 In the Koopman operator picture, the map N c defines the dynamical system on discrete time, with an associated Koopman operator K Nc . An important result in [28] is that the map N c is conjugate to the polynomial g 0 (z) = z 2 through the map which is a Möbius transformation, and hence invertible on C ∪ {∞}. Thus, eigenfunctions of K Nc can be constructed through eigenfunctions of K g 0 , which are defined through Since g 0 (z) = z 2 , there are eigenfunctions φ k of the form φ k (z) = ln(|z|) k , associated to the eigenvalues λ k = 2 k . Through the given conjugacy h 0 , we have that are eigenfunctions for K Nc , associated to the eigenvalues λ k . Fig. 5 shows ψ 1 for c = 1, i.e.
for the polynomial f (z) = z 2 + 1. The roots at ±i are clearly visible as the points where the eigenfunction diverges.

An example with chaotic behavior.
Consider the function f (z) = z 2 + 1 from the previous section, with its two roots at ±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 so that z n+1 = N f (z n ) yields a new iterate, starting at z 0 ∈ C. For z 0 ∈ C/R, the limit points of the iterative scheme are the extremal points, where lim n→∞ z n = i Im(z 0 ) > 0, lim n→∞ z n = −i Im(z 0 ) < 0. For z 0 ∈ R, however, the scheme does not converge (see Fig. 6) and follows the chaotic recurrence (which can be solved explicitly in this case)  We employ the numerical method by Korda, Putinar and Mezić [20] (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 µ g for g(z) = 1 + exp(2πiz), with z ∈ R from a trajectory of the Newton-Raphson method starting at z 0 = 1/2. The point spectrum is σ pp = {1} with eigenvalue 1 associated to the constant function, as expected for a chaotic system. The absolutely continuous spectrum seems to be σ c = T/{1}, with no singularly continuous part. The trajectories were computed with the iteration of N f (equation (4.3)) with 16 digits of accuracy, and with the explicit formula (equation (4.4)) 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 ≈ 10 31,000 for the last point in the longest trajectory.

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 [28]. 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.

Conclusions.
We studied the Koopman operator of several dynamical systems that represent the action of algorithms. In a series of examples, we demonstrated that the Koopman operator provides 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 already 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.
flow map on the entire set X is not defined for n = 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 [11]. Eigenvalues with a real part larger than one indicate unstable behavior, indicating the existence of unstable nodes or saddles in the data set. . The spectral density on the unit circle T, associated to the observable g(z) = 1 + exp(2πiz), is approximated with single trajectories of Newton's method with length 10 3 , 10 4 , and 10 5 (rows) starting at z0 = 1/2 + 0i. The point spectrum with eigenvalue at 1 + 0i associated to the constant function and the continuous spectrum on T/{1} can be distinguished. The blue curves are obtained from iterating the map N f with 16 digits of accuracy, the orange curves by using the exact formula and 35,000 digits of accuracy.   Figure 11. At the beginning, the right steady state is not inside the data domain (orange rectangle, A). When both steady states are contained in the data, the maximum of the real part of all eigenvalues is close to one (B). As the available data set moves to the right, the left steady state is no longer available to the approximation (C). When the rectangle is moved even further (D), the saddle between the two steady states also disappears and only one, attracting steady state is contained in the data.