Towards tensor-based methods for the numerical approximation of the Perron-Frobenius and Koopman operator

The global behavior of dynamical systems can be studied by analyzing the eigenvalues and corresponding eigenfunctions of linear operators associated with the system. Two important operators which are frequently used to gain insight into the system's behavior are the Perron-Frobenius operator and the Koopman operator. Due to the curse of dimensionality, computing the eigenfunctions of high-dimensional systems is in general infeasible. We will propose a tensor-based reformulation of two numerical methods for computing finite-dimensional approximations of the aforementioned infinite-dimensional operators, namely Ulam's method and Extended Dynamic Mode Decomposition (EDMD). The aim of the tensor formulation is to approximate the eigenfunctions by low-rank tensors, potentially resulting in a significant reduction of the time and memory required to solve the resulting eigenvalue problems, provided that such a low-rank tensor decomposition exists. Typically, not all variables of a high-dimensional dynamical system contribute equally to the system's behavior, often the dynamics can be decomposed into slow and fast processes, which is also reflected in the eigenfunctions. Thus, the weak coupling between different variables might be approximated by low-rank tensor cores. We will illustrate the efficiency of the tensor-based formulation of Ulam's method and EDMD using simple stochastic differential equations.


Introduction
The Perron-Frobenius operator and the Koopman operator enable the analysis of the global behavior of dynamical systems. Eigenfunctions of these operators can be used to extract the dominant dynamics, to detect almost invariant sets, or to decompose the system into fast and slow processes [1,2,3,4,5,6,7]. Assume the state space of your system is R d and you want to discretize each direction using n grid points (or boxes if Ulam's method is used), then overall n d values need to be stored. Even for d = 10 and n = 10, more than 70 gigabyte of storage space would be required, whereas typical systems might have hundreds or thousands of dimensions and naturally also require a more fine-grained discretization. This so-called curse of dimensionality can be overcome by using tensor formats which compress the data and store only the information that is relevant for the reconstruction. In general, only an approximation of the original data can be retrieved. Approximating the objects under consideration by sums of low-rank tensor products has become a powerful approach for tackling high-dimensional problems [8] and in many physically significant problems near-linear complexity can be achieved since the separation rank depends only weakly on the dimension [9]. For high-dimensional systems exhibiting multiscale behavior, it might be possible to represent the weak coupling between different variables by low-rank tensor cores. The leading eigenfunctions of the Koopman operator, for instance, are typically almost constant for the fast variables of the system and depend mainly on the slowly changing variables (see [6]). Thus, using tensor-based algorithms could reduce the amount of time and memory required to compute and store eigenfunctions significantly. In this way, analyzing high-dimensional systems that could not be tackled using standard methods might become feasible.
Tensors, in our sense, are just multidimensional arrays as shown in Figure 1.
Here and in what follows, standard vectors will be denoted by lower-case letters, e.g. v, matrices by upper-case letters, e.g. A, and tensors by the corresponding bold symbols, e.g. x. It is important to note that tensors are typically not explicitly given -for example by observed data -, but only implicitly as solutions of systems of linear or nonlinear equations or eigenvalue problems [10]. Thus, numerical methods that operate directly on tensor approximations need to be developed since the full tensors cannot be stored or handled anymore in practice.
Over the last years, low-rank tensor approximation approaches have become increasingly popular in the scientific computing community and are now becoming a standard tool to cope with large-scale problems that could not be handled before by standard numerical methods. An overview of different low-rank tensor approximation approaches can be found in [10].
In this paper, we will show that the use of low-rank tensor approximation schemes potentially enables the computation of eigenfunctions of high-dimensional systems. The aim of this paper, however, is not to show that our approach is more efficient -the tensor algorithms are mainly implemented in Matlab, a comparison with highly optimized numerical libraries implemented in C or C++ would not lead to meaningful results, developing high performance libraries for large-scale tensor problems is a separate task -, but to derive a tensor-based reformulation of existing methods and to show equivalency so that the theory available for the conventional matrix-vector based formulation can be carried over to multi-dimensional arrays. Furthermore, tensors could enable low-rank approximations of the Perron-Frobenius and Koopman operator as well as their eigenfunctions. One of the main future goals is to combine low-rank tensor decomposition techniques and the splitting of the dynamics into fast and slow processes. In [6], it has been shown that such a splitting of a multi-scale system exists and can be exploited to extract the slow dynamics. Another open problem is the generation of the low-rank approximations of the operators. Currently, the canonical tensor format representations are converted to the tensor-train format, which is time-consuming and in general leads to high ranks. Ideally, low-rank TT approximations should be directly generated from the given data.
We will start by introducing standard methods such as Ulam's method and the recently developed Extended Dynamic Mode Decomposition (EDMD) to approximate the eigenfunctions of the Perron-Frobenius operator and the Koopman operator in Section 2. Then, in Section 3, different tensor formats will be described. In Section 4, we will reformulate Ulam's method and EDMD as tensor-based methods. Section 5 contains a brief summary of simple power iteration schemes for the resulting tensor-based (generalized) eigenvalue problems. Simple examples which illustrate the proposed approaches are shown in Section 6. In Section 7, we will conclude with a short summary and possible future work.

Perron-Frobenius and Koopman operator approximation
In this section, we will briefly introduce the Perron-Frobenius operator P and the Koopman operator K as well as numerical methods to compute finite-dimensional approximations, namely Ulam's method and EDMD. The main difference between Ulam's method and EDMD is that the former uses indicator functions 1 for a given box discretization of the domain while the latter allows arbitrary ansatz functions such as monomials, Hermite polynomials, trigonometric functions, or radial basis functions. Although EDMD was primarily developed for the approximation of the Koopman operator, it can be used to compute eigenfunctions of the Perron-Frobenius operator as well [12]. Analogously, Ulam's method can also be used to compute eigenfunctions of the Koopman operator.

Perron-Frobenius and Koopman operator
Let S : X → X be a dynamical system defined on a domain X , for example X ⊆ R d . Then the Perron-Frobenius operator or transfer operator P is defined by for all f, g ∈ F, where F is an appropriately defined function space and • denotes function composition. We assume in what follows that F = L 2 (X ). The aim is to compute eigenfunctions of the Perron-Frobenius operator, given by The eigenfunction ϕ 1 corresponding to λ 1 = 1 is the invariant density of the system, i.e. Pϕ 1 = ϕ 1 . The magnitude of the second largest eigenvalue λ 2 can be interpreted as the rate at which initial densities converge to the invariant density (for more details and assumptions about the dynamical system, see e.g. [6] and references therein). More generally, the leading eigenvalues of the Perron-Frobenius operator close to one correspond to the slowly converging transients of the system. The Koopman operator K, on the other hand, is defined by and acts on functions f : X → C, f ∈ F. Correspondingly, the stochastic Koopman operator is defined by denotes the expected value with respect to the probability measure underlying S(x). We will only introduce the required notation and focus mainly on discrete-time dynamical systems, for more details on the Koopman operator and its properties, we refer to [4,5,13]. While the Perron-Frobenius operator describes the evolution of densities, the Koopman operator describes the evolution of observables, which could be measurements or sensor probes [4]. Instead of analyzing orbits {x, S(x), S 2 (x), . . . } of the dynamical system, we now analyze the measurements {f (x), f (S(x)), f (S 2 (x)), . . . } at these points. The Koopman operator K is the adjoint of the Perron-Frobenius operator P and thus an infinite-dimensional but linear operator. A finite-dimensional approximation (computed using generalized Galerkin methods) of this operator captures the dynamics of a nonlinear dynamical system without necessitating a linearization around a fixed point [4,5]. We are again interested in eigenfunctions of the operator, given by Let f : X → R be an observable of the system that can be written as a linear combination of the linearly independent eigenfunctions ϕ i , i.e.
Analogously, for vector-valued functions F = [f 1 , . . . , f n ] T , we obtain where v i = [c i,1 , . . . , c i,n ] T . These vectors v i corresponding to the eigenfunctions ϕ i are called Koopman modes.

Ulam's method
A frequently used method to compute an approximation of the Perron-Frobenius operator is Ulam's method, see e.g. [2,14,15,6]. First, the state space X is covered by a finite number of disjoint boxes {B 1 , . . . , B k }. Let 1 Bi be the indicator function for box B i , i.e.
Then a finite-dimensional approximation of the operator can be obtained as follows: Using definition (1) leads to This relationship can be represented by a matrixP = (p ij ) ∈ R k×k witĥ Here, in order to avoid confusion with the approximation P on a tensor space introduced below, we denote the matrix representationP instead of P . The denominator m(B i ) normalizes the entriesp ij so thatP becomes a rowstochastic matrix and defines a finite Markov chain. The left eigenvector corresponding to the eigenvalue λ 1 = 1 approximates the invariant measure of the Perron-Frobenius operator P.
The entriesp ij of the matrixP represent the probabilities of points being mapped from box B i to box B j by the dynamical system S. These entries can be estimated by randomly choosing a large number of test points x ( The eigenfunctions of the Perron-Frobenius operator are then approximated by the left eigenvectors of the matrix P , the eigenfunctions of the Koopman operator by the right eigenvectors.

Extended dynamic mode decomposition
An approximation of the Koopman operator, the Koopman eigenvalues, eigenfunctions, and eigenmodes can be computed using EDMD. The method requires data, i.e. a set of values x i and the corresponding y i = S(x i ) values, i = 1, . . . , m, written in matrix form as X = x 1 x 2 · · · x m and Y = y 1 y 2 · · · y m , and additionally a set of ansatz functions or observables The vectors x i are used as collocation points to approximate the integrals required for the approximation of the Koopman operator. Let Ψ : X → R k , be the vector of all ansatz functions, then K can be approximated by a matrixK ∈ R k×k , witĥ where + denotes the pseudoinverse. The matricesÂ,Ĝ ∈ R k×k are defined aŝ As before, we use theˆsymbol to distinguish the matrices from the tensor approximations that will be introduced in Section 4. An approximation of the eigenfunction ϕ i of the Koopman operator K is then given by where ξ i is the i-th left eigenvector of the matrixK T . Alternatively, the generalized eigenvalue problem can be solved, provided thatĜ is regular. To compute eigenfunctions of the Perron-Frobenius operator, the corresponding eigenvalue problem ξ iÂ T = λ i ξ iĜ needs to be solved. Note that this formulation is similar to the variational approach to compute eigenfunctions of transfer operators of reversible processes presented in [16,17]. For more details, we refer the reader to [12].

Tensor formats
Several different tensor formats have been developed in the past, e.g. the canonical format, the Tucker format, and the tensor-train format. In this section, we will briefly introduce tensors and the required notation. The overall goal is to rewrite the methods presented in the previous section as tensor-based methods and to take advantage of lowrank tensor approximations and the fact that the dynamics of high-dimensional systems can often be decomposed.

Full format
A tensor in full format is simply a multidimensional array v ∈ R k1×···×k d . (A variation of this format is the sparse format which stores only the nonzero entries and is used, for example, in the sparse grid approach [18].) The entries of a tensor v are indexed by Addition and subtraction are trivially defined element-wise. Multiplication of a tensor v by a scalar c ∈ R is naturally generalized Matrix-vector multiplication is defined as follows: Given a linear operator A defined on a tensor space R k1×···×k d , with or in shorthand notation, using multi-indices, and the outer product v ⊗ w ∈ R k1×···×k d ×k1×···×k d as a tensor with entries The outer product v ⊗ w can be regarded as a linear map that acts on tensors R k1×···×k d . Often it is required or convenient to rewrite a tensor as a vector. The vectorization of a tensor, denoted vec(v), where vec : R k1×···×k d → R k1···k d , reorders the entries of v into one column vector. For v ∈ R 2×3×2 , for example, For our purposes, we will mainly be interested in eigenvalue problems of the form Av = λv or Av = λBv.
The treatment of tensors in the full format often leads to storage problems. Thus, different formats have been developed to overcome this problem. Instead of working with the full format, we will use compressed formats such as the r-term or TT format for numerical computations in order to minimize computational costs as well as storage requirements. Except for very particular examples, it is impossible to compress the data without any compression error [18]. Typically, the tensor representation is just an approximation of the original data. Below, we will describe different compressed tensor formats, the introduction is based on [18].

Canonical format
A tensor space is given by . . , V d are vector spaces defined over the same field K, typically R or C. An elementary tensor is defined to be a product of the form This format is also called r-term format or CP format. That is, instead of trying to store the tensor in the dense format, one only considers tensors that can be written as products of the form If the best approximation of this form is not good enough, then the natural extension is to consider [19,9]. Defining Example 3.1. Let us consider a system of the where V is the energy landscape associated with the system. The invariant density of this process is given by µ(x) = 1 Z e −βV (x) , where β and Z are constants [20]. Assume that the state space is three-dimensional and that the potential function can be written as Thus, storing the invariant density for a grid with k grid points in each direction in the full format would require an array of size k 3 while storing it in the canonical tensor format would require only a tensor of rank 1 and thus an array of size 3k.
Given tensors in the r-term format, basic operations are defined as follows [18]: • Matrix-vector multiplication: Since these operations increase the rank, truncation is typically required to approximate the resulting tensor v with rank r v by a tensorṽ with a lower rank rṽ by either fixing the rank rṽ or by fixing ε such that v −ṽ < ε.

TT format
Another frequently used tensor format is the TT format -TT now stands for tensor train instead of the former tree tensor -, which can be obtained by successive singular value decompositions. This is a special case of a more general hierarchical format and has been introduced in quantum physics under the name Matrix Product States (MPS), see [18] for details. A tensor v ∈ R k1×···×k d is decomposed into d component tensors v i of at most order three (the first and last are of order two and are often, for the sake of simplicity, considered as tensors of order three with "boundary condition" ρ 0 = ρ d = 1). That is, the entries of v are given by For fixed indices, the component tensors of rank three can be regarded as matrices, which leads to a more compact and justifies the original name MPS. Here, the numbers ρ i are the ranks of the TT tensor, resulting in a rank vector ρ = [ρ 0 , ρ 1 , . . . , ρ d ] which determines the complexity of the representation. The main advantages of the TT format are its stability from an algorithmic point of view and reasonable computational costs, provided that the ranks of the tensors are small [8]. The basic operations such as addition and matrix-vector multiplication are more complex than in the canonical format and can be found, for example, in [21]. Converting a tensor from the canonical format to the TT format is trivial, but the TT representation requires more memory. Numerical toolboxes for the TT decomposition of tensors and several algorithms for solving linear systems of equations are available online, see e.g. [22].

Comparison
Complexity-wise, the canonical format would be the ideal candidate for representing tensors since the number of required parameters depends only linearly on the dimension d, the rank r, and the sizes of the individual vector spaces. It turned out, however, that solving even simple problems using the canonical format is hard in practice due to redundancies and instabilities which can lead to numerical problems [8]. The main advantage of the TT format is its structural simplicity, higher-order tensors are reduced to d products of tensors of at most order three. Similar approaches have been known in quantum physics for a long time, the rigorous mathematical analysis, however, is still work in progress (see [8] and references therein). For our purposes, we will rely on the TT format and the TT toolbox developed by Oseledets et al. [22] and implement simple power iteration schemes to solve the resulting eigenvalue problems as we will show in Section 5. Other tensor formats, however, might be advantageous as well for the analysis of the Perron-Frobenius and Koopman operator. This should be investigated further in the future. One drawback of the TT format is that the decomposition depends on the ordering of the dimensions and thus results in different tensor ranks for different orderings.

Tensor-based approximation
In this section, we will present a tensor-based reformulation of Ulam's method and EDMD and show that these methods are equivalent to the corresponding vector-based counterparts. The new formulation enables the use of the low-rank tensor approximation approaches described in the previous section. That is, variables can be approximated with different degrees of accuracy.

Reformulation of Ulam's method
Given a dynamical system S : X → X , X ⊂ R d , define a box B that contains X , i.e.
Using again standard multi-index notation, we will denote i = (i 1 , . . . , i d ). Equipped with the mapping each multi-index i corresponds to a numberî ∈ {1, . . . ,k}. This induces a canonical numbering of the boxes Bî and the entries of tensors x ∈ R k1×···×k d such that x i = xî, where x = vec(x). Thus, with the aid of Ulam's method we could now generate the finite-dimensional representation of the Perron-Frobenius operatorP ∈ Rk ×k as described in Section 2. Our goal, however, is to approximate the operator by a tensor P ∈ R k1×···×k d ×k1×···×k d . Note that the indicator function for the box B i can be written as That is, each d-dimensional indicator function 1 B i (x) is now written as a product of d one-dimensional indicator  Figure 2(a). Thus, using Ulam's method, we would obtain 9 indicator functions {1 B1 , . . . , 1 B9 } and the matrixP that approximates the Perron-Frobenius operator P would be a row-stochastic (9×9)-matrix. An example of such a matrix is shown in Figure 2(b), the underlying dynamical system is not relevant here. The goal now is to rewrite this matrix using tensors (cf. Example 4.3). (a) for µ = 1, 2 and i µ = 1, 2, 3, the ansatz functions for the box discretization can be written as The product formulation of the indicator functions naturally leads to a tensor approximation P of the Perron-Frobenius operator P. Let Q µ : R d → R be the projection onto the µ-th component of a vector, i.e. Q µ (x) = x µ . Then we define  (2) and (7) For d = 1, the multi-index i is mapped toî = i 1 and j toĵ = j 1 by (6). Furthermore,k = k 1 and Then, we obtain by induction Here, the matrices and vectors with the superscripts (i d+1 , j d+1 ) and (j d+1 ), respectively, are obtained by fixing the corresponding indices, that is, these matrices and vectors are lower-dimensional slices of the corresponding higherdimensional objects. Since the entries of v (j d+1 ) are indexed by (j 1 , . . . , j d ) →ĵ = 1 + d µ=1 µ−1 ν=1 k ν (j µ − 1), the entries of the larger vector v -obtained by stacking the vectors v (j d+1 ) ∈ R k1···k d -are indexed by Analogously, left eigenvectors of P correspond to left eigenvectors ofP . The implementation of the tensorbased formulation is straightforward since only index computations for intervals are required, a numbering of the d-dimensional boxes is not needed anymore. Let T be the set of all test points and ind : R d → N d the functions that returns the corresponding multi-index i for a point x ∈ R d so that x µ ∈ I iµ µ , µ = 1, . . . , d. Then, Ulam's method can simply be expressed as: n end for In the standard formulation, the rows of the matrixP sum up to one. Correspondingly, the sum of all entries of each subtensor of P with fixed multi-index i is one, i.e.  Instead of working with the full format, the matrix P can also be expressed directly using the canonical tensor format. Assume that a test point x is mapped from box B i to box B j , where i = (i 1 , . . . , i d ) and j = (j 1 , . . . , j d ) are again multi-indices. Now let e iµ µ ∈ R kµ be the i µ -th unit vector of size k µ and let Then the elementary tensor e i ⊗ e j ∈ R k1×···×k d ×k1×···×k d describes the mapping of this point. Furthermore, let ind be again the function that returns the multi-index of the box that contains the point x and T the set of all test points, then the matrix P can be represented as a sum of elementary tensors of this form, i.e. That is, the number of elementary tensors iskn, i.e. the number of boxes multiplied by the number of test points per box, and thus potentially too large to store (unless sparse tensor formats are used). However, we will store only a low-rank approximation to reduce the required storage space. Note that this elementary tensor representation can also be easily converted into the TT format for numerical computations using the TT toolbox.
The question now is whether the tensor representation offers advantages over the standard formulation of Ulam's method. The goal is to approximate the eigenfunctions of the Perron-Frobenius operator or Koopman operator using low-rank tensors, reducing the computational cost as well as the memory consumption. Before we present numerical results, let us also rewrite EDMD in tensor form.

Reformulation of EDMD
Instead of writing Ψ as a vector of functions Ψ = [ψ 1 , ψ 2 , . . . , ψk] T , we now write Ψ as a tensor of functions. We start by selecting basis functions for each dimension separately. Let where i = (i 1 , . . . , i d ) is a multi-index. Thus, That is, we have againk = That is, Such a tensor basis is often used for high-dimensional problems, see also [23]. Typical basis functions are monomials, Hermite polynomials, or trigonometric functions. In the standard formulation, all basis functions are enumerated and rewritten in vector form (3). The difference here is that the tensor form will be preserved. We will use again (6) as a mapping from multi-index to single index when required. Now A, G ∈ R k1×···×k d ×k1×···×k d can be constructed as follows: The entries are again -as in the standard EDMD formulation -approximated using a collocation approach. EDMD computes the entries as shown in (4) or in short form, using the outer product, which in turn results in a generalized eigenvalue problem of the form ξA = λξG.
Note that the eigenvalue problem is the same as (5) in the standard case. For the sake of simplicity, we are omitting the index i here.
Proposition 4.5. Provided that the basis functions can be written in tensor product form (9), Proof. We just show that the entries of A andÂ as well as the entries of G andĜ are identical, the rest -the equivalency of the matrix-vector and tensor products -follows from Proposition 4.2. Assuming that the basis can be written in the product form, we obtain from (4) Instead of storing the dense matrices A and G, we can again directly represent these matrices using the canonical tensor format. The basis was chosen in such a way that Ψ(x) can be written as (10) it follows that A and G can be written as sums of m elementary tensors. As before, we are not storing the full-rank tensor, but only low-rank approximations.
The eigentensors ξ ∈ R k1×···×k d of the generalized eigenvalue problem can then be used to approximate the eigenfunctions of the Perron-Frobenius operator or Koopman operator: Let ξ be a left eigentensor, then approximates an eigenfunction of the Koopman operator. Analogously, if ξ is a right eigentensor of the generalized eigenvalue problem -observe that G[i 1 , . . . , i d , j 1 , . . . , j d ] = G[j 1 , . . . , j d , i 1 , . . . , i d ] -, then ϕ(x) is an approximation of the corresponding eigenfunction of the Perron-Frobenius operator, see also [12].
To compute the dominant eigenfunctions of the Koopman operator or Perron-Frobenius operator, we will use simple power iteration schemes outlined in the next section. General purpose eigenvalue solvers for nonsymmetric generalized eigenvalue problems are, to our knowledge, not part of the tensor libraries yet. Solvers for symmetric (non-generalized) eigenvalue problems already exist and are part of the TT toolbox [22].

Eigenvalue problems
In Section 2 and Section 4, we have shown that in order to compute the eigenfunctions of the Perron-Frobenius operator and Koopman operator, respectively, using either Ulam's method or EDMD, we need to solve standard eigenvalue problems or generalized eigenvalue problems. For the reformulated version of these methods, we have to develop the required numerical algorithms to solve the resulting tensor-based eigenvalue problems. At the time of writing, we are not aware of any tensor toolbox containing numerical methods for nonsymmetric generalized eigenvalue problems. Methods for the computation of eigenvectors of symmetric positive definite matrices in the TT format have been proposed in [8], where the eigenvalue problem is rewritten as a (Rayleigh quotient based) minimization problem which is then solved using the Alternating Linear Scheme (ALS). In practice, these methods have recently also been successfully used for nonsymmetric problems, although convergence has not been shown yet [24].
Suitable methods for eigenvalue problems can be subdivided into two main categories as explained in [10] (see also references therein, e.g. [9]): The first category of methods is based on combining classical iterative algorithms with low-rank truncation after each step, the second is based on a reformulation as an optimization problem, where admissible solutions are constrained to the set of low-rank tensors. In this section, we will describe a generalization of simple power iteration methods -belonging to the first category -to tensor-based eigenvalue problems. For a detailed description of general power iteration methods, we refer to [25]. Power iteration and inverse power iteration for tensors have also been proposed in [9]. The main difference between the standard algorithms and the tensorbased counterpart is that for the latter truncation is used to keep the ranks of the tensors low. It is important that the iteration moves from the initial state to the final state without creating intermediate solutions with an excessive rank [9].

Power iteration methods for standard eigenvalue problems
In what follows, let T denote the truncation of a tensor. Then instead of a classical iteration scheme of the form x k+1 = F (x k ), we simply obtain x k+1 = T (F (x k )). For an eigenvalue problem of the form Av = λv, given an initial guess v 0 for the dominant eigenvector, the power iteration algorithm computes: The iteration converges to an eigenvector associated with the largest eigenvalue λ 1 of the truncated operator A if the eigenvalue is simple and the initial guess v (0) has a component in the direction of the corresponding dominant eigenvector v 1 [25]. Even if the initial guess does not have a component in the direction of v 1 , rounding errors typically ensure that this direction will be picked up during the iteration. The rate of convergence depends on the ratio between the second-largest and largest eigenvalue λ 2 /λ 1 . The main advantage of this method is that it requires only matrix-vector multiplications and can thus easily be used for tensor eigenvalue problems.
A modification of this algorithm to compute eigenvectors corresponding to any eigenvalue is the inverse power iteration with shift, which -assuming that A − θI is nonsingular -can be written in the form: end for The parameter θ is called shift and the iteration converges to the eigenvalue closest to θ. This method is just the standard power iteration applied to the matrix (A − θI) −1 . Here, the linear solver computes a low-rank approximation of the solution so that truncation is not required.

Power iteration methods for generalized eigenvalue problems
Given matrices A and B, the generalized eigenvalue problem is given by Av = λBv. In this case, the power iteration method also requires the solution of a linear system of equations. Thus, we can also directly apply the inverse power iteration, where the resulting systems of linear equations are again solved with ALS: In order to keep the ranks of the intermediate solutions low, we also approximate the matrices P (Ulam's method) or A and G (EDMD) by low-rank tensors. That is, the initial matrices are converted toP,Ã, andG with a lower rank since the rank of these matrices can initially be very high. We are typically only interested in the general behavior of the eigenfunctions, a highly accurate representation of the eigenfunctions is often not needed.

Examples
All examples presented within this section have been implemented in Matlab using -for the sake of efficiencymex-functions to integrate the SDEs. All tensor computations were carried out with the TT toolbox [22]. For the eigenvector computations, we used our implementation of the simple power iteration methods described in Section 5.

2-dimensional double well problem
Let us start with a simple 2-dimensional example, a stochastic differential equation of the form where W 1 and W 2 are two independent standard Wiener processes. Here, the potential is given by see Figure 3 (cf. [12]). Furthermore, we set σ = 0.7. Note that in this case the potential can be written as . In order to analyze the tensor-based methods, we rotate the potential by an angle α and obtain The two independent Wiener processes W 1 and W 2 are rotated accordingly. We would expect that the eigenfunctions of systems with small α can be accurately approximated by low-rank tensors, whereas systems with a larger value of α require higher ranks since the dynamics are not aligned with the axes anymore. The second eigenfunctions of the Perron-Frobenius operator for the systems with potential V and different values of α computed using the tensor-based version of Ulam's method are shown in Figure 4. We chose α = 0, α = π/12, α = π/6, and α = π/4. The domain X = [−2, 2] 2 was subdivided into 50 × 50 equally sized boxes. That is, P ∈ R 50×50×50×50 . For each box, 100 randomly chosen test points were generated. The Euler-Maruyama method with a step size h = 10 −3 was used for the numerical integration, where one evaluation of S corresponds to 10, 000 integration steps, that is, the integration interval is [0, 10]. The shift parameter θ of the power iteration method was set to a value slightly smaller than 1. Figure 5 illustrates how the tensor approximation, depending on the rank, successively picks up the information about the shape of the eigenfunction and generates more and more accurate representations. For α = π/4, a tensor of rank 1 cannot represent the minimum and maximum  simultaneously since this would lead to two additional peaks in the lower left and upper right corner. Here, the first pair of singular vectors represents the maximum, the second pair of singular vectors the minimum. Additionally, we computed the eigenfunctions with the standard version of Ulam's method to evaluate the accuracy of the approximation and compared it with the results obtained by using the new tensor-based formulation. Figure 6 shows the influence of the truncation of the operator as well as the influence of the truncation of the resulting eigenfunctions. Here, in order to analyze the accuracy, we also compare the first eigenfunction with the analytically computed invariant density. Since we are computing eigenfunctions of the Perron-Frobenius operator associated with a stochastic differential equation, the results depend strongly on the number of test points chosen for each box. The higher the number of test points per box, the smoother the eigenfunction approximation. Thus, in this case, the smoother low-rank solutions can counterintuitively lead to better approximations of the true eigenfunctions. The high ranks are mainly required to resolve the numerical noise introduced by the coarse approximation of the operator. This can be seen, for example, in Figure 6b. Decreasing the rank initially reduces the error -the truncation of the operator results in smoother eigenfunctions -until the shape of the eigenfunction cannot be described by a low-rank approximation anymore and the error increases. For α = 0, the x 1 and x 2 dynamics are independent and a low-rank approximation is sufficient. Furthermore, the results illustrate that for a fixed-rank approximation, the error is smaller when the system's dynamics are aligned with the axes.
This example shows that in order to be able to approximate eigenfunctions by low-rank tensors, the dynamics of the system should be aligned with the axes chosen, although even if the dynamics are not aligned, the tensor format might be advantageous. In general, the dynamics are unknown a priori and not necessarily aligned with the axes, but for higher-dimensional systems it is often possible to decompose a system into slow and fast subsystems. Not all variables of a system might be equally important to describe the system's behavior. The intuition would be that certain subsystems require less information and that the tensor approximation automatically captures the relevant dynamics, using high ranks only when necessary.

3-dimensional triple well problem
Let us consider a more complex 3-dimensional example with the potential function This is a potential taken from [20], augmented by a third dimension. We subdivided the domain X = [−2, 2] × [−1, 2] × [−2, 2] into 20 × 20 × 20 boxes of the same size. For each box, we randomly generated 1000 test points. The resulting finite-dimensional approximation is then a tensor P ∈ R 20×20×20×20×20×20 . Figure 7 shows a scatter plot of the first three eigenfunctions, which were computed using the inverse power iteration described in Section 5. The first eigenfunction clearly shows the three expected regions with high probabilities corresponding to the minima of the potential V . The second eigenfunction separates the two deeper wells at (−1, 0, 0) and (1, 0, 0), the third eigenfunction separates these wells from the third shallower well at (0, 1.5, 0). The differences between the eigenfunctions computed using the conventional and the tensor-based version of Ulam's method are negligible, the average difference between the first eigenvector v 1 and the tensor v 1 is of the order of 10 −6 , which is mainly due to the less accurate power iteration method applied to the tensor eigenvalue problem. For the sake of comparison, we computed the eigenfunctions of the Koopman operator using EDMD. We chose basis functions D = {x i1 1 x i2 2 x i3 3 , i 1 , i 2 , i 3 = 0, . . . , 5}. Hence, A, G ∈ R 6×6×6×6×6×6 . For higher-order monomials, the resulting matrices are ill-conditioned and the eigenfunctions cannot be computed accurately anymore. The results are shown in Figure 8. As before, the second eigenfunction separates the two deeper wells. The third eigenfunction separates these wells from the third shallower well and is close to zero for the regions around the deep wells. The trivial eigenfunction of the Koopman operator corresponding to λ 1 = 1 is not plotted here since it is constant and does not contain relevant information about the system. Note that the eigenvalues are slightly different due to the different set of basis functions used for EDMD.

Conclusion
We have reformulated the problems of computing finite-dimensional approximations of the Perron-Frobenius and Koopman operator in a different format using tensors instead of vectors. The matrices P (if Ulam's method is used) or A and G (if EDMD is used) can now either be assembled in the dense tensor format or directly in the canonical tensor format -which can then easily be converted into the TT format -, enabling low-rank approximations of the aforementioned operators. The next step is to systematically develop the numerical methods required to efficiently solve the resulting nonsymmetric generalized tensor eigenvalue problems and also to store and handle these tensors minimizing memory requirements so that even high-dimensional problems can be solved. First results obtained by applying simple algorithms such as power iteration methods are promising and show that the approaches presented within this paper might be able to tackle high-dimensional problems. Currently, several toolboxes for tensor-based problems are under development. Once these toolboxes contain methods for solving nonsymmetric generalized eigenvalue problems, the proposed approaches can be implemented easily, potentially facilitating the computation of meta-stable sets or almost invariant sets of dynamical systems that could not be handled before due to the a) b) c) d) Figure 6: (Top) Letv 1 denote the first eigenvector ofP , v 1 the (vectorized) eigentensor of the truncated tensor representation P, and µ inv the analytically computed invariant density. The error here is defined by e = 1 k v 1 −v 1 2 and e = 1 k v 1 − µ inv 2 , respectively, where k is the number of boxes. a) Difference between v 1 andv 1 depending on the rank of P. b) Difference between v 1 and µ inv . (Bottom) Influence of the truncation of the first eigentensor v 1 of the full tensor representation P on the accuracy. c) Difference between the truncated eigentensor v 1 and v 1 . d) Difference between the truncated eigentensor v 1 and µ inv . The dashed lines show the error for the full-rank approximation which is almost identical for the different values of α.
curse of dimensionality. We demonstrated the tensor-based version of Ulam's method and EDMD using two-and three-dimensional problems mainly because the results can be easily validated and visualized.
Future work also includes determining which tensor format is suited best for our purposes. Currently, one of the main bottlenecks is the simulation required to obtain the data. Even simple molecular dynamics simulations, for instance, might easily take several days, but in order to capture the behavior of the system, long trajectories or a large number of short simulations with different initial conditions are required. This huge amount of data must then be processed. Thus, also the construction of the matrices P or A and G is time-consuming. The number of boxes or basis functions required to represent each variable of the system accurately is in general unknown a priori. If we are, for instance, only interested in the leading eigenfunctions of the Koopman operator, the fast variables of the system are typically almost constant and require less information to be captured. Starting with a set of only a few basis functions for each unknown, which is then, if needed, augmented adaptively based on the system's behavior would greatly improve the efficiency. Adaptive methods combined with (sparse) tensor approaches might be able to tackle high-dimensional systems and diminish the curse of dimensionality. Furthermore, a detailed numerical analysis of the efficiency and accuracy of the proposed algorithms would help understand the limitations and find Figure 7: Scatter plot of the first three eigenfunctions of the Perron-Frobenius operator for the triple well system. Only entries whose absolute value is larger than a given threshold are plotted, entries close to zero are omitted. Another open problem is the -depending on the number of dimensions d -typically extremely large condition number of the matrices A or G if EDMD with, for instance, monomials are used to compute the eigenfunctions of the Perron-Frobenius or Koopman operator. Hence, the resulting eigenvalue problems cannot be solved accurately anymore for high-dimensional systems. A detailed understanding and numerical analysis of different basis functions might help mitigate this problem. Radial-basis functions or other more locally defined functions could lead to better results.