Koopman Spectral Linearization vs. Carleman Linearization: A Computational Comparison Study

Nonlinearity presents a significant challenge in problems involving dynamical systems, prompting the exploration of various linearization techniques, including the well-known Carleman Linearization. In this paper, we introduce the Koopman Spectral Linearization method tailored for nonlinear autonomous dynamical systems. This innovative linearization approach harnesses the Chebyshev differentiation matrix and the Koopman Operator to yield a lifted linear system. It holds the promise of serving as an alternative approach that can be employed in scenarios where Carleman linearization is traditionally applied. Numerical experiments demonstrate the effectiveness of this linearization approach for several commonly used nonlinear dynamical systems.


Introduction
Nonlinear dynamical systems find extensive applications in both science and engineering.However, analyzing and computing solutions for nonlinear systems can pose a formidable challenge.A natural strategy for mitigating this complexity involves approximating a nonlinear system through a technique known as linearization.The most conventional method of linearization entails employing a first-order Taylor expansion of the system dynamics.While this technique has proven effective and is widely used in Model Predictive Control (MPC) and similar domains, it remains confined to a very narrow vicinity around the expansion point.To comprehensively address these limitations, a series of more sophisticated techniques and theories have emerged.Approaches such as Carleman linearization and the Koopman operator are rooted in the concept of lifting nonlinear systems to infinite-dimensional linear systems.These global linearization techniques and theories were developed by scholars like Torsten Carleman [8], Bernard Koopman [12], and John Von Neumann [13] in the 1930s.In practice, by employing finite section approximation, we can attain a finite yet higher-dimensional linear system.This empowers us to capture the characteristics of the original nonlinear system over a longer time period compared with the first-order Taylor expansion approach.Within the control system community and quantum computing community, the application of Carleman linearization concepts has yielded a range of notable accomplishments across various studies.For instance, in [24], the author emphasizes the applicability of Carleman Linearization in designing nonlinear controllers and state estimators.The work in [18] harnessed Carleman approximation to establish a correlation between the lifted system and the domain of attraction of the original nonlinear system.Noteworthy recent endeavors, including [10], skillfully harnessed Carleman linearization to implement efficient model predictive control mechanisms for nonlinear systems.The application of Carleman linearization concepts was extended to [1], where they were adeptly utilized for tasks such as state estimation and the formulation of feedback control laws.Moreover, [2,4] exploited the inherent structure of the elevated system to introduce a practical methodology for quantitatively addressing and solving the Hamilton-Jacobi-Bellman equation through a meticulously crafted iterative approach.Furthermore, in the realm of quantum computing, Carleman linearization has also found prominence.Leveraging its ability to transform nonlinear dynamical systems into linear forms, researchers have employed Carleman linearization to design efficient quantum algorithms, with provable expoential speed up with respect to dimension compared with best possible classic algorithms [17,11,14].
As highlighted in reference [21], the Carleman matrix is, in fact, a transposed finite section matrix approximation of the Koopman Operator Generator in the polynomial basis.Despite this connection, modern research on the Koopman Operator and Carleman Linearization follows distinct paths.Investigations into the Koopman Operator predominantly adopt a data-driven approach, as they use spatiotemporal data without requiring the explicitly equation to find an approximation of the principal Koopman modes and correspondingly Koopman eigenpairs.Notably numerical schemes like Dynamic Mode Decomposition (DMD) [26,25,5] and its diverse variations, including Extended Dynamic Mode Decomposition (EDMD) [29], make use of data collected from time snapshots of a specified dynamical system.These approaches employ Singular Value Decomposition (SVD) and other matrix decomposition techniques to reveal the temporal and spatial spectral attributes of the system.Expanding beyond DMD, EDMD and their variants, methods for obtaining tractable representations have been further connected with deep neural networks.[20,19] Two specific neural network architectures have been explored in this context.One is based on the autoencoder principle, employing a low-dimensional latent space to enhance the interpretability of outcomes.The other architecture involves a higherdimensional latent space, which often yields improved results when dealing with dynamical systems featuring continuous spectra.A comprehensive review of the theoretical framework and algorithms related to Data-Driven Koopman Operators is provided by [7].
Diverging from the previously discussed data-driven methods, our approach constitutes a progressive evolution of the Linearization Technique.This distinction arises from our specific need for an analytical equation.More precisely, our method integrates the Chebyshev Differentiation Matrix in the spectral method [28,27], yielding a clear matrix approximation of the Koopman Operator's generator.Subsequently, this matrix is employed to generate an explicit higher-dimensional linear equation that closely approximates the original nonlinear system.In this context, the identical approximation matrix was employed across various scalar observable configurations, with only the initial value being altered.This innovative approach was motivated by the numerical ordinary differential equation (ODE) algorithm discussed in [16,15].This study centers around autonomous systems and potentially serves as a viable alternative approach applicable in situations where Carleman linearization has conventionally been utilized.
The paper is structured as follows: section 2 presents the background knowledge, followed by the discussion of the Koopman Spectral Linearization Method in section 3.In section 4, we provide the numerical results, and the subsequent sections delve into the discussion and conclusions in section 5.

Carleman Linearization
Carleman Linearization provides a systematic methodology for deriving the corresponding lifted dynamics from an initial value problem defined as follows [9]: Here the matrices B i ∈ R d×d i are time-independent, k is the degree of the polynomial ODE and x ⊗i denotes the j-fold tensor product for each positive integer j.We define an infinite-dimensional vector y = (y 1 , y 2 , y 3 , ....) T with y i = x ⊗i , and formulate matrices A i i+j−1 Where ⊗ i B = I ⊗ I ⊗ I ⊗ ... ⊗ B ⊗ ... ⊗ I with B at (d − i + 1)-th position.Further, formula (2) lead to the following equation Formally, this can be represented as an infinite-dimensional ordinary differential equation: where A is an infinite-dimensional matrix defined as: In practice, the system is often truncated at a certain order N.Here is an illustrative example [7] of Carleman Linearization applied to the following system: Carleman Linearization procedure can derive a linear ode like

Koopman Operator
Firstly consider d-dimension (i.e.x ∈ X ⊂ R d ) dynamical system described by a first order autonomous order ordinary differential equations with t ∈ [t 0 , T ].
where x = (x 1 , x 2 , ..., x d ) T belongs to a space X and this dynamics f is also an d dimensional vector valued function with (f 1 , f 2 , ..., f d ) T serve as each coordinate.For our analysis, we introduce observables or measurement functions, represented by the function g : X → R, note this is a function on G(X) (i.e.L 2 Hilbert space).And the flow map F t represents the evolution of the system dynamics as a mapping on X.
The Koopman operator family, denoted as {K t }, is defined as: Alternatively, we can express K t as a composition of functions: It's crucial to note that the Koopman operator is linear, owing to the linearity of the function space G(X).Further we can define the generator of the Koopman operator K The generator can be explicitly formulated as follows: The transformation of representing it in the form of an operator can be further extended as follows: As a trade-off, we convert a nonlinear mapping on the variable x into a linear operator on the variable g, which, in turn, leads to an infinite-dimensional space.
In Koopman Spectral Theory, eigenvalues and eigenfunctions (or eigenpairs for brevity) are employed for practical numerical computations [22,25].Suppose (λ, ϕ(x)) is an eigenpair for the generator K.By definition Based on ode in equation ( 15), we can derive the following outcomes: Actually, ϕ(x) is also an eigen-function for K t ,this can be observed from Therefore if we denote e λt as µ than (µ, ϕ) = (e λt , ϕ) is an eigenpair for K t .
If g ∈ G(X) ⊂ Span{ϕ k }, then even nonlinear observable can be represented as a linear combination of eigenfunctions.This can be formulated as g = j c j ϕ j , where c j ∈ R is the corresponding coefficient of g with respect to ϕ j .Further as pointed in [22], the evolution of the observable can also be represented as a linear combination of eigenpairs which is derived from: Hence, we can infer Multi-dimensional observables are based on same manner.Just consider vector-valued c j with same dimension of the observable.c j is also known as j-th Koopman mode [7].

Construction of the lifted matrix
By adopting construction of the matrix approximation from [16] and [15], the generator of Koopman Operator can be approximated as following: Currently, when d=1, for some eigen-pairs (ϕ, λ) based on equation ( 7) To derive a finite dimensional approximation, we starting from polynomial interpolation of the eigenfunction ϕ(x) on spatial discretized Guass-Labatto points {ξ i } N i=1 .Which give us with L j are Lagrange Polynomials such that L j (ξ i ) = δ ij where δ ij is the delta function.Therefore we obtained . . .
Based on equation (22), the finite representation of K is When d = 2, Let {ξ i } N i=1 ,and{η j } N j=1 be the Guass-Lobatto Points of x 1 and x 2 .And denote Θ = {(ξ i , η j )} N i,j=1 Now for every bivariate eigenfunction ϕ(x 1 , x 2 ) we have polynomial interpolation Hence all function values on the collocation points can be formed as matrix denoted as ϕ N (Θ) Let D 1 , D 2 be the differentiation matrices for x 1 and x 2 respectively, and f 1 (Θ),f 2 (Θ) be the matrices of f 1 f 2 evaluated at (ξ i , η j ).And we denote Kϕ N (Θ) ij = Kϕ N (ξ i , η j ) Then the formula below can be derived based on (7) We vectorize all the matrix operations as vector operations, this process can be formulated as: Therefore, for d = 2 cases, the construction of matrix K is given by In general for a d-dimensional case, consider Guass-Labotta points for each coordinate as And use Θ denotes the collection of all the collocation points i.e.Θ = {(ξ 1 i1 , ξ 2 i2 , ..., ξ d i d ) : i 1 , i 2 , ..., i d , traverse from 1 to N}.And the eigenfunction ϕ is approximated by Corresponingly, ϕ N (Θ) is a tensor of eigenfunction values on collocation points.and D 1 , D 2 , ...D d means differentiation matrices.With n-mode multiplication in tensor algebra, we can write Pursuing the analogous vectorization approach in equation ( 27), we can reformulate equation (30) as Consequently, Where -position, to be more intuitive, here is an example when d = 3.

Solution of the linear system
Now, let's examine the solution of the linear system to provide further insight into why it serves as an approximation of the original nonlinear one.Referring back to equation ( 19), we have.
When d = 1, for a scalar observable g.Consider Guass-Lobatto points on a region with radius r around x 0 , where Θ = {ξ i } N i=1 and ξ 1 < ξ 2 < ... < ξ N .Here we pick odd number N such that ξ (N +1)/2 = x 0 = x(t 0 ).Thus all Guass-Lobatto Points are on the interval [x 0 − r, x 0 + r].Suppose ( λj , v j ) is an eigen-pair of matrix K. Vector v j are used to approximate eigenfunction ϕ N j evaluated at the collocation points.we have (v j ) i = ϕ N j (ξ i ) ≈ ϕ j (ξ i ).Hence Now consider the N-dimensional linear system below We can write the analytical solution y(t) = e Kt y 0 (37) Suppose the eigen-decompostion of K = V ΛV −1 , with each columns of V are eigenvectors of K, denoted as v j , and Λ contain all the eigenvalues.Further by setting t = 0 in equation ( 35), we have This is also true for different initial value thus Thus we have V −1 y 0 = c since y(0) = (g N (ξ 1 ), g N (ξ 2 ), ..., g N (ξ N ) T , where c is called the Koopman mode.
e λ1t e λ2t . . . And . Notice c j are Koopman modes, e λj t are eigenfunctions.In general case, for d-dimensional nonlinear system, suppose Θ represents d-dimensional collocation points, and K is constructed as equation (32).Now consider the linear system given below Again by vectorization, ϕ j (x 0 ) can be approximated by the "middle term" of tensor ϕ N (Θ) which leads to And to compute the observable, only the middle term is needed for this long vector.
which is exactly the solution of (41).Unlike Carleman Linearization, which captures the information of all variables within a single linear ordinary differential equation, our linearization approach is primarily tailored for scalar observables.When the need arises to compute multi-dimensional observables or distinct scalar observables, we can still employ the previously constructed lifted matrix K.However, this requires transforming the initial conditions accordingly.For instance, when computing an observable g(x) = x = (g 1 (x), g 2 (x), g 3 (x)),we would need to utilize three different initial conditions: g 1 (Θ), g 2 (Θ), g 3 (Θ) This adaptation allows us to extend the applicability of our approach to scenarios involving diverse observable functions since K can be reused.
As for the solution of Carleman Linearization, the matrix A shall be truncated at certain level N which yields A N matrix as following: Then the infinite linear system can be approximated by and the approximated solution of original nonlinear system is exactly the first d elements of the solution for equation (45).

Numerical Results
In this section, we present the performance of Koopman Spectral Linearization on five nonlinear dynamic models in subsection 4.1.In each example, we investigate the influence of truncation order N and radius r on accuracy.The reference solution is generated by RK4 if there is no closed-form solution available.To be more precise, in two or three-dimensional cases, we test the influence of the radius in each direction separately.Next, in subsection 4.2, we compare the accuracy and computational efficiency involved with Carleman Linearization.These comparisons include errors against truncation order, matrix size, and computational time.

Quadratic Model
The quadratic model is possible simplest nonlinear ode just designed for demonstration purpose.The governing ODE is given by To be more clarified, we set x(0) = 0.08 and T=10 in this example.This quadratic model has a closed form solution x(t) = .Test 1(a) shows that the exponential convergence of our approach with respect to truncation order, which is similar to the conclusion in conventional Carleman Linearization Method.Test 1(b) illustrates that the the error shows a "bowl shape" with respect to the radius r i.e. r cannot be too large or too small.shows that the exponential convergence of our approach with respect to truncation order, which is still similar to the conclusion in conventional Carleman Linearization Method.Further, test 2(b) illustrates that the the error shows a "V shape" with respect to the radius r (i.e.r cannot be too large or too small as in the first example).

Simple Pendulum Model
The simple pendulum is well studied model in science and engineering.The movement of the pendulum is described by a second order ordinary different equation, Where L is the length of the pendulum, θ is the displacement angle, and the parameter g is the gravity acceleration.This second order equation can be converted to a two dimensional first order ode system.As a notation we define x 1 = θ and x 2 = dθ dt .Also, we set L = g = 9.8 for simplicity.Correspondingly, And we set x(0) = (0.1, 0.1) T , T = 10.The three experiments are under the following parameters settings: • Test of solution: r = (1, 1) T ; • Test of Truncation Order N : r = (1, 1) T ; • Test of radius r: N=9

Lotka-Volterra Model
The Lotka-Volterra model, also known as the predator-prey model [6], is a mathematical model used to describe the interactions between predators and prey in an ecosystem.The model describes the interaction between two fundamental populations: the prey population and the predator population.Here we defined as below We set x(0) = (5, 5) T and T = 10.The parameters used in three different tests are as follows: • Test of solution: r = (3, 3); • Test of Truncation Order N : r = (3, 3) T ; • Test of radius r: N = 13  Again error decrease exponentially with respect to truncation order and "V shape" radius test results was observed.In particular, different from previous example, the results obtained between two different dimensions exhibit a consistent pattern, but the precision of the numerical outcomes for the second dimension is consistently higher in any given test.

Compared with Carleman Linearization
Following the construction procedure introduced in section 3, we apply the tensor product rule to construct collocation points.Consequently, the approximation matrix of the Koopman generator has a size of N d × N d .However, for multi-dimensional dynamics, we can reuse the lifted matrix with different initial values.In the case of the Carleman linearization procedure, the matrix size will be N i=1 d i × N i=1 d i , and we can approximate the size of the Carleman lifted matrix as O(d N ).Therefore, when considering the matrix size of the derived linear system with a relatively small truncation order N , the matrix in our approach can be significantly smaller compared to the matrix in Carleman linearization.To further illustrate this point, 1 presents a matrix size comparison between our approach and Carleman linearization.While our approach generates d linear systems with the same matrix but different initial values, which may not be as convenient as Carleman linearization, where only one linear system is produced, the significant difference in matrix sizes results in a notable increase in computational cost.As an illustrative example, we present the running time and error in 2 for a truncation order of N =9.It is evident that the computational time required for Carleman Linearization is slower than our approach, especially when d = 2 or 3.In particular, for the Kraichnan-Orszag Model (d=3), the Carleman linearization procedure took 32.0546 seconds, which is considerably slower than our approach.
Figure 6 illustrates the error as a function of the truncation order N and provides a comparison between Carleman linearization and our method.The error of Carleman linearization exhibits exponential decay, a phenomenon that has been substantiated by recent theoretical research [9,3].Our method leverages the Chebyshev node interpolation approach, which also enjoys exponential convergence guarantees and demonstrates similar exponential convergence in numerical experiments.For all the polynomial examples in 6(a), 6(b), and 6(c), it is evident that our approach exhibits higher accuracy at the same truncation order, and the error decreases more rapidly as the truncation order is reduced.
Moreover, when applied to non-polynomial dynamical systems, our method showcases superior numerical accuracy and faster convergence rates.In Figure 7, we present a comparative analysis of the error-truncation order relationship in the context of non-polynomial dynamical systems.Both nonpolynomial dynamics are subjected to 12th-order Taylor expansions to transform them into polynomial forms.It is evident that the Carleman method exhibits significantly lower accuracy in this scenario, with notably slower convergence.This limitation arises from the fact that, in non-polynomial cases, the relevant part of the Taylor expansion is confined to orders lower than the truncation order, resulting in substantial errors when employing relatively small expansion orders, as demonstrated in our study.Consequently, our method demonstrates heightened precision in non-polynomial cases.Even in the case of the 'simple pendulum' where N =3, our method initially exhibits lower accuracy.However, as N increases to 9, the errors obtained by our method rapidly become dramatically smaller than those produced by the Carleman method.With regard to limitations, our method's capacity to yield numerical results with smaller matrices, reduced computational time, and increased accuracy, all under the same truncation order, hinges upon a critical prerequisite-namely, the selection of an appropriate parameter r, a requirement not present in the Carleman Linearization method.As demonstrated in the results of the parameter r testing discussed in subsection 4.1, deviating r from the optimal setting point or range, either by making it excessively large or excessively small, results in gradually increasing errors.Furthermore, when only the time-domain evolution differential equations of the dynamical system are available in advance, without access to spatial evolution information of the state, our method cannot capture the state's spatial evolution.Therefore, unless supplementary information is accessible regarding the spatial evolution of the dynamical system, allowing for a reasonable estimation of r, surpassing the precision of Carleman Linearization remains unattainable.

Conclusion and Discussion
The Koopman Spectral Linearization method uses the Chebyshev Differentiation Matrix to explicitly derive the lifted matrix of nonlinear autonomous dynamical systems.It provides a finitedimensional representation for the generator of the Koopman Operator with a polynomial basis.Therefore, similar to Carleman Linearization, our approach provides a linearization method for nonlinear dynamical systems.In each numerical experiment presented, our approach exhibits exponential convergence with respect to truncation order N , just like Carleman Linearization.Under the same truncation order, Koopman Spectral Linearization tends to exhibit significantly higher accuracy associated with lower time cost compared to Carleman Linearization, especially when the dynamics are not polynomial.Therefore, it is more suitable as an alternative linearization method where high-accuracy approximations are needed and f is non-polynomial.Different from Carleman Linearization, which describes the evolution of the state itself, our approach is used to describe a scalar smooth observable, which is more flexible.
However, despite its advantages, the Koopman Spectral Linearization method is not without its limitations.It's important to consider certain drawbacks when evaluating its applicability since a reasonable estimation of radius r is necessary to obtain an accurate approximation of the eigenfunctions.Therefore if no additional information like range of states are given, its difficult to provide an accurate approximation.To further understand how this parameter r influence the accuracy, a comprehensive numerical analysis will be included in the future work.Besides, the global interpolation method might overcome the limitation of estimate r.
Regarding the matrix size as indicated in subsection 4.2, Koopman Spectral Method obtained a much smaller matrix with a even higher accuracy compared with Carleman Linearization.However the matrix size is exponential increase with respect to dimension d (i.e.N d ) since tensor product rule is applied.Therefore both Carleman Linearization (i.e.O(d N )) and our approach might be not efficiency for relatively high-dimensional problem.A possible way to overcome this limitation is to adapt sparse grid methods to construct collocation points which has been founded success in [15].As pointed in Koopman Operator theory, either Carleman Linearization and our approach relies on tensor product rule, therefore cannot deal with high dimension problems.By combining with sparse grid technique, our approach is possible to address much higher-dimensional problems [15].
Furthermore, an interesting potential application about our approach is to address efficient quantum algorithm for nonlinear ODE by combination of quantum linear system algorithm to solve the linearized system [17].Especially our approach only require the middle element of the output vector which might significantly improve the process of extract quantum information.

1 ( 1 / 11 Figure 1
Figure1summarizes these results in plots 1(a) and 1(b).Test 1(a) shows that the exponential convergence of our approach with respect to truncation order, which is similar to the conclusion in conventional Carleman Linearization Method.Test 1(b) illustrates that the the error shows a "bowl shape" with respect to the radius r i.e. r cannot be too large or too small.
(a) Test of Truncation Order (b) Test of Radius

Figure 1 : 3 • 9 Figure 2
Figure 1: Quadratic Model: (a) Testing the influence of truncation order N (b) Testing the influence of radius r (a) Test of Truncation Order (b) Test of Radius

Figure 2 :
Figure 2: Cosine Square Model: (a) Testing Truncation Order N ; (b) Testing the influence of radius r

Figure 3 :Figure 3
Figure 3: Simple Pendulum Model: (a) Testing Solution; (b) Testing Truncation Order N ; (c) Testing the influence of radius r

Figure 4 :Figure 4
Figure 4: Lotka-Volterra Model: (a) Testing Solution; (b) Testing Truncation Order N ; (c) Testing the influence of radius r

( a )
Test of Truncation Order (b) Test of Radius r 1 (c) Test of Radius r 2 (d) Test of Radius r 3

Figure 7 :
Figure 7: Error-Truncation Order compared with Carleman Linearization non-polynomial case (a) Cosine Square Model; (b) Simple Pendulum Model