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

Abstract

Nonlinearity presents a significant challenge in developing quantum algorithms involving differential equations, prompting the exploration of various linearization techniques, including the well-known Carleman Linearization. Instead, this paper introduces the Koopman Spectral Linearization method tailored for nonlinear autonomous ordinary differential equations. This innovative linearization approach harnesses the interpolation methods and the Koopman Operator Theory to yield a lifted linear system. It promises to serve 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 ordinary differential equations.

1. Introduction

Differential equations find extensive applications in diverse fields such as fluid dynamics, biology, and finance. For instance, in fluid dynamics, the Navier–Stokes equations serve as fundamental tools for modeling fluid behavior [1]. In epidemiology, ordinary differential equations, such as those based on the Susceptible–Infectious–Recovered (SIR) model, are indispensable for analyzing the spread of infectious diseases [2]. Additionally, in finance, the Black–Scholes model relies on differential equations to facilitate option pricing [3]. However, computing solutions for high-dimensional differential equations presents a significant challenge. This challenge arises from the curse of dimensionality, wherein traditional numerical methods entail discretizing equations at grid points. As the dimension of equations increases, the number of grid points grows exponentially. Consequently, even basic linear algebra computations within these numerical algorithms become prohibitively complex.

A promising strategy for mitigating this complexity involves using quantum computing techniques since quantum computer leverages unique features such as quantum entanglement and superposition, enabling exponential speedups in specific computational tasks. As previous research has pointed out, replacing the linear systems solvers in classical numerical methods with their quantum counterparts, the Quantum Linear Systems Algorithm (QLSA), [4] can yield provably efficient quantum algorithms for solving linear ordinary differential equations (ODEs) [5,6,7,8] and linear Partial Differential Equations (PDEs) [7,9,10,11]. Additionally, quantum algorithms based on Hamiltonian simulation, rather than the QLSA, have also been theoretically shown to efficiently solve computational problems involving linear differential equations, whether for the standard Hamiltonian simulation or non-Hermitian Hamiltonian simulation [12].

However, despite a series of theoretical successes in developing quantum algorithms for linear differential equations, addressing nonlinear cases remains an open problem. Early attempts [13] to directly combine QLSA with nonlinear ODEs proved to have poor scaling with evolution time. Recently, Liu et al. [14] demonstrated significant improvements by introducing Carleman Linearization to certain nonlinear ODEs. Using QLSA on the linearized differential equations greatly enhanced the complexity bound with respect to the evolution time. This approach has motivated many researchers to consider first employing linearization techniques, such as Carleman Linearization [15,16,17], Koopman von Neumann linearization [18,19,20] and the level set method [21], to approximate or represent low-dimensional nonlinear differential equations with high-dimensional linear ones, and then using quantum algorithms designed for linear differential equations to solve these problems. Figure 1 illustrates the work flow of these various quantum algorithms for solving nonlinear ODEs.

As highlighted in reference [22], the Carleman matrix is, in fact, a transposed finite-section matrix approximation of the Koopman Operator’s generator on a polynomial basis. This specific connection and the weaker performance of Carleman Linearization in nonpolynomial cases motivate our research. 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 explicit equation to find an approximation of the principal Koopman modes and, correspondingly, Koopman eigenpairs. In particular, numerical schemes such as dynamic mode decomposition (DMD) [23,24,25] and its diverse variations, including Extended Dynamic Mode Decomposition (EDMD) [26], make use of the 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 its variants, methods for obtaining tractable representations have been further connected with deep neural networks [27,28]. 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 higher-dimensional 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 [29].

Diverging from the previously discussed data-driven methods, our approach constitutes a progressive evolution of the linearization technique. This distinction mainly arises from our method of constructing the infinitesimal generator for the approximated Koopman operator. Rather than relying on data, we utilize the values from the right-hand side function of the ODEs at various discrete spatial points, thereby necessitating an explicit equation. More precisely, our method integrates the differentiation matrix in the spectral method [30,31], yielding an explicit matrix approximation of the Koopman Operator’s generator. Subsequently, this matrix generates 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 spectral method for differential equations discussed in [32,33]. Therefore, like other spectral methods, our approach’s effectiveness and convergence rate are influenced by the choice of the basis function. Potential options for the elements of D include polynomials [34], Fourier modes [32], radial basis functions [35], and spectral elements [36]. The best choice of basis functions is likely dependent on the specific characteristics of the underlying dynamical system and the data sampling strategy. Specifically, we adapt the Chebyshev Differentiation Matrix for demonstration purposes. This study focuses on autonomous systems and potentially serves as an alternative approach applicable in situations where Carleman Linearization has been conventionally utilized.

In addition to applications in quantum computing, our methods also hold potential for applications in traditional computational contexts. The most conventional method of linearization involves using a first-order Taylor expansion of the system dynamics. This technique, while proven effective and widely utilized in areas such as Model Predictive Control (MPC) and related domains, is limited to applications within a very narrow vicinity around the expansion point. To address these limitations more comprehensively, a range of more sophisticated linearization techniques such as Carleman Linearization [37] and the Koopman Operator [38] have been utilized. These methods were originally applied in the fields of dynamical systems and control [22].

The paper is structured as follows: Section 2 presents the background knowledge and discusses 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.

2. Background

2.1. Carleman Linearization

Carleman Linearization provides a systematic methodology for deriving the corresponding lifted dynamics from an initial value problem defined as follows [37]:

$$\frac{d\mathit{x}}{dt}=B_1\mathit{x}+B_2{\mathit{x}}^{\otimes 2}+\cdots +B_k{\mathit{x}}^{\otimes k}\mathrm{with}{\mathit{x}}_0=\mathit{x}\left(t_0\right)$$

(1)

Here, the matrices $B_i\in {\mathcal{R}}^{d\times d^i}$ are time-independent, k is the degree of the polynomial ODE, and ${\mathit{x}}^{\otimes i}$ denotes the i-fold tensor product for each positive integer i. We define an infinite-dimensional vector $\mathit{y}={({\mathit{y}}_{\mathbf{1}},{\mathit{y}}_{\mathbf{2}},{\mathit{y}}_{\mathbf{3}},\dots .)}^T$ with ${\mathit{y}}_{\mathit{i}}={\mathit{x}}^{\otimes i}$, and formulate matrices $A_{i+j-1}^i$

$$\begin{array}{c}\hfill A_{i+j-1}^i=\sum _{v=1}^i{\otimes }_i{B}_j\end{array}$$

(2)

where ${\otimes }_{i}B=I\otimes I\otimes I\otimes \cdots \otimes B\otimes \cdots \otimes I$ with B at the $(d-i+1)$-th position. Further, Formula (2) leads to the following equation:

$$\begin{array}{c}\hfill \frac{d{\mathit{y}}_{\mathit{i}}}{dt}=\sum _{j=1}^k{A}_{i+j-1}^i{\mathit{y}}_{\mathit{i}+\mathit{j}-\mathbf{1}}\end{array}$$

(3)

Formally, this can be represented as an infinite-dimensional ordinary differential equation:

$$\begin{array}{cc}\hfill \frac{d\mathit{y}}{dt}& =\mathcal{A}y\mathrm{with}{\mathit{y}}_{\mathit{i}}\left(t_0\right)=\mathit{x}{\left(t_0\right)}^{\otimes i}\hfill \end{array}$$

(4)

where $\mathcal{A}$ is an infinite-dimensional matrix defined as:

$$\begin{array}{c}\hfill \mathcal{A}=\left[\begin{array}{ccccccccc} & A_1^1& A_2^1& A_3^1& \cdots & A_k^1& 0& 0& \cdots \\ & 0& A_2^2& A_3^2& \cdots & A_k^2& A_{k+1}^2& 0& \cdots \\ & 0& 0& A_3^3& \cdots & A_k^3& A_{k+1}^3& A_{k+2}^3& \cdots \\ & ⋮& ⋮& ⋮& & ⋮& ⋮& ⋮& \end{array}\right]\end{array}$$

(5)

In practice, the system is often truncated at a certain order N. Here is an illustrative example [37] of Carleman Linearization applied to the following system:

$$\begin{array}{cc}\hfill \frac{dx}{dt}& =x^2\hfill \end{array}$$

(6)

The Carleman Linearization procedure can derive a linear ODE-like

$$\frac{d}{dt}\left[\begin{array}{c}x\\ x^2\\ x^3\\ ⋮\end{array}\right]=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& \cdots \\ 0& 0& 2& 0& 0& \cdots \\ 0& 0& 0& 3& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮\end{array}\right]\left[\begin{array}{c}x\\ x^2\\ x^3\\ ⋮\end{array}\right]$$

(7)

2.2. Koopman Operator

Firstly, consider a d-dimension (i.e., $\mathit{x}\in X\subset {\mathbb{R}}^d$) dynamical system described by first-order autonomous order ordinary differential equations with $t\in [t_0,T]$:

$$\frac{d\mathit{x}\left(t\right)}{dt}=\mathit{f}\left(\mathit{x}\left(t\right)\right)$$

(8)

where $\mathit{x}={(x_1,x_2,\cdots ,x_d)}^T$ belongs to a space X, and this dynamics $\mathit{f}$ is also a d-dimensional vector valued function with $f={(f_1,f_2,\cdots ,f_d)}^T$. For our analysis, we introduce observables or measurement functions, represented by the function $g:X\to \mathbb{R}$; note that this is a function on $\mathcal{G}\left(X\right)$ (e.g., an ${\mathcal{L}}^2$ Hilbert space). And the flow map $F_t$ represents the evolution of the system dynamics as a mapping on X:

$$F_t:\mathit{x}\left(t_0\right)\to \mathit{x}(t+t_0)$$

(9)

The Koopman operator family, denoted as $\left\{{\mathcal{K}}_t\right\}$, is defined as:

$${\mathcal{K}}_t:g\left(\mathit{x}\left(t_0\right)\right)\to g\left(F_t\left(\mathit{x}\left(t_0\right)\right)\right)=g\left(\mathit{x}(t_0+t)\right)$$

(10)

Alternatively, we can express ${\mathcal{K}}_t$ as a composition of functions:

$${\mathcal{K}}_t:g\to g\circ F_t$$

(11)

It is crucial to note that the Koopman operator is linear, owing to the linearity of the function space $\mathcal{G}\left(X\right)$. Further, we can define the generator of the Koopman operator $\mathcal{K}$

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

(12)

The generator can be explicitly formulated as follows:

$$\mathcal{K}g=\frac{\partial g}{\partial t}=\nabla g\left(\mathit{x}\right)\frac{d\mathit{x}}{dt}=\nabla g\left(\mathit{x}\right)·f\left(\mathit{x}\right)$$

(13)

The transformation of representing it in the form of an operator can be further extended as follows:

$$\mathcal{K}=\nabla ·\mathit{f}\left(\mathit{x}\right)=\sum _{i=1}^d{f}_i\left(x\right)\frac{\partial }{\partial x_i}$$

(14)

As a trade-off, we convert a nonlinear mapping on the variable x into a linear operator on the variable g, which leads to an infinite-dimensional space.

In Koopman Spectral Theory, eigenvalues and eigenfunctions (or eigenpairs for brevity) are employed for practical numerical computations [24,39]. Suppose $(\lambda ,\varphi (\mathit{x}\left)\right)$ is an eigenpair for the generator $\mathcal{K}$. By definition,

$$\begin{array}{c}\hfill \mathcal{K}\varphi \left(\mathit{x}(t+t_0)\right)=\frac{d\varphi \left(\mathit{x}(t+t_0)\right)}{dt}=\lambda \varphi \left(\mathit{x}(t+t_0)\right)\end{array}$$

(15)

Based on the ODE in Equation (15), we can derive the following outcomes:

$$\varphi \left(\mathit{x}(t+t_0)\right)=\varphi \left(\mathit{x}\left(t_0\right)\right)e^{\lambda t}$$

(16)

Actually, $\varphi \left(\mathit{x}\right)$ is also an eigenfunction for ${\mathcal{K}}_t$; this can be observed from

$${\mathcal{K}}_t\varphi \left(\mathit{x}(t+t_0)\right)=\varphi \left(\mathit{x}(t_0+2t)\right)=e^{\lambda t}\varphi \left(\mathit{x}(t_0+t)\right)$$

(17)

Therefore, if we denote $e^{\lambda t}$ as $\mu$, then $(\mu ,\varphi )=(e^{\lambda t},\varphi )$ is an eigenpair for ${\mathcal{K}}_t$.

If $g\in \mathcal{G}\left(X\right)\subset Span\left\{{\varphi }_k\right\}$, then even the nonlinear observable can be represented as a linear combination of eigenfunctions. This can be formulated as $g={\sum }_j{c}_j{\varphi }_j$, where $c_j\in \mathbb{R}$ is the corresponding coefficient of g concerning ${\varphi }_j$. Further, as pointed out in [39], the evolution of the observable can also be represented as a linear combination of eigenpairs which is derived from:

$${\mathcal{K}}_{t}g\left(\mathit{x}(t+t_0)\right)={\mathcal{K}}_t\left(\sum _{j=1}^{\infty }{c}_j{\varphi }_j\left(\mathit{x}(t_0+t)\right)\right)=\sum _{j=1}^{\infty }{c}_j{e}^{{\lambda }_{j}t}{\varphi }_j\left(\mathit{x}(t_0+t)\right)=\sum _{j=1}^{\infty }{c}_j{\mu }_j{\varphi }_j\left(\mathit{x}(t_0+t)\right)$$

(18)

Hence, we can infer

$$g\left(\mathit{x}(t+t_0)\right)=\sum _{j=1}^{\infty }{c}_j{e}^{{\lambda }_{j}t}{\varphi }_j\left(\mathit{x}\left(t_0\right)\right)$$

(19)

Multi-dimensional observables are based in the same manner. Just consider vector-valued $c_j$ with the same dimension of the observable. $c_j$ is also known as the j-th Koopman mode [29].

3. Koopman Spectral Linearization Method

3.1. Construction of the Lifted Matrix

By adopting a construction of the matrix approximation from [30,31], the generator of the Koopman Operator can be approximated as follows.

Currently, when d = 1, for some eigenpairs $(\varphi ,\lambda )$ based on Equation (7),

$$\mathcal{K}\varphi =f·\frac{\partial }{\partial x}·\varphi$$

(20)

To derive a finite-dimensional approximation, we start from a polynomial interpolation of the eigenfunction $\varphi \left(x\right)$ on spatial discretized Gauss–Lobatto points ${\left\{{\xi }_i\right\}}_{i=1}^N$, which gives us

$$\varphi \left(x\right)\approx {\varphi }^N\left(x\right)=\sum _{i=1}^N{\varphi }^N\left({\xi }_i\right)L_i\left(x\right)$$

(21)

where $L_j$ are Lagrange polynomials such that $L_j\left({\xi }_i\right)={\delta }_{ij}$, where ${\delta }_{ij}$ is the delta function. Therefore, we obtain

$$K\left[\begin{array}{c}{\varphi }^N\left({\xi }_0\right)\\ \\ {\varphi }^N\left({\xi }_1\right)\\ \\ ⋮\\ \\ {\varphi }^N\left({\xi }_N\right)\end{array}\right]=\left[\begin{array}{cccc}f\left({\xi }_0\right)& & & \\ & f\left({\xi }_1\right)& & \\ & & \ddots & \\ & & & f\left({\xi }_N\right)\end{array}\right]D\left[\begin{array}{c}{\varphi }^N\left({\xi }_0\right)\\ \\ {\varphi }^N\left({\xi }_1\right)\\ \\ ⋮\\ \\ {\varphi }^N\left({\xi }_N\right)\end{array}\right]$$

(22)

Based on Equation (22), the finite representation of K is

$$K=diag\{f\left({\xi }_0\right),f\left({\xi }_1\right),\dots ,f\left({\xi }_N\right)\}D,$$

(23)

where $diag$ creates a diagonal matrix with the main diagonal elements. When $d=2$, let ${\left\{{\xi }_i\right\}}_{i=1}^N$ and ${\left\{{\eta }_j\right\}}_{j=1}^N$ be the Gauss–Lobatto points near $x_1$ and $x_2$. And denote $\Theta ={\left\{({\xi }_i,{\eta }_j)\right\}}_{i,j=1}^N$. Now, for every bivariate eigenfunction $\varphi (x_1,x_2)$, we have polynomial interpolation ${\varphi }^N$

$$\varphi (x_1,x_2)\approx {\varphi }^N(x_1,x_2)=\sum _{i=1}^N\sum _{j=1}^N{\varphi }^N({\xi }_i,{\eta }_j)L_i\left(x_1\right)L_j\left(x_2\right)$$

(24)

Hence, all function values on the collocation points can be formed as a matrix denoted as ${\varphi }^N(\mathbf{\Theta })$

$${\varphi }^N(\mathbf{\Theta })=\left[\begin{array}{cccc}{\varphi }^N({\xi }_1,{\eta }_1)& {\varphi }^N({\xi }_1,{\eta }_2)& \cdots & {\varphi }^N({\xi }_1,{\eta }_N)\\ {\varphi }^N({\xi }_2,{\eta }_1)& {\varphi }^N({\xi }_2,{\eta }_2)& \cdots & {\varphi }^N({\xi }_2,{\eta }_N)\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }^N({\xi }_N,{\eta }_1)& {\varphi }^N({\xi }_N,{\eta }_2)& \cdots & {\varphi }^N({\xi }_N,{\eta }_N)\end{array}\right]$$

(25)

Let $D_1,D_2$ be the differentiation matrices for $x_1$ and $x_2$, respectively, and $f_1(\mathbf{\Theta })$,$f_2(\mathbf{\Theta })$ be the matrices of $f_1$, $f_2$ evaluated at $({\xi }_i,{\eta }_j)$. And we denote $\mathcal{K}{\varphi }^N{(\mathbf{\Theta })}_{ij}=\mathcal{K}{\varphi }^N({\xi }_i,{\eta }_j)$. Then, the formula below can be derived based on (7):

$$\mathcal{K}{\varphi }^N(\mathbf{\Theta })=f_1(\mathbf{\Theta })\odot \left(D_1{\varphi }^N(\mathbf{\Theta })\right)+f_2(\mathbf{\Theta })\odot \left({\varphi }^N(\mathbf{\Theta })D_2^T\right)$$

(26)

In the equation above, ⊙ means the Hadamard product. We vectorize all the matrix operations as vector operations via $vec$, which flatten a matrix into a column vector; this process can be formulated as:

$$\begin{array}{cc}\hfill \mathcal{K}vec\left({\varphi }^N(\mathbf{\Theta })\right)& =vec\left(f_1(\mathbf{\Theta })\right)\odot \left((I\otimes D_1)vec\left({\varphi }^N(\mathbf{\Theta })\right)\right)+vec\left(f_2(\mathbf{\Theta })\right)\odot \left((D_2\otimes I)vec\left({\varphi }^N(\mathbf{\Theta })\right)\right)\hfill \\ \hfill & =[diag\left(vec\left(f_1(\mathbf{\Theta })\right)\right)(I\otimes D_1)+diag\left(vec\left(f_2(\mathbf{\Theta })\right)(D_2\otimes I)\right)]vec({\varphi }^N(\mathbf{\Theta })\hfill \\ \hfill & =Kvec({\varphi }^N(\mathbf{\Theta })\hfill \end{array}$$

(27)

Therefore, for $d=2$ cases, the construction of matrix K is given by

$$K=diag\left(vec\left(f_1(\mathbf{\Theta })\right)\right)(I\otimes D_1)+diag\left(vec\left(f_2(\mathbf{\Theta })\right)(D_2\otimes I)\right)$$

(28)

In general, for a d-dimensional case, consider Gauss–Lobatto points for each coordinate as ${\left\{{\xi }_{i_1}^1\right\}}_{i_1=1}^N$,${\left\{{\xi }_{i_2}^2\right\}}_{i_2=1}^N$,...,${\left\{{\xi }_{i_d}^d\right\}}_{i_d=1}^N$. And use $\Theta$ to denote the collection of all the collocation points, i.e., $\Theta =\{({\xi }_{i_1}^1,{\xi }_{i_2}^2,\dots ,{\xi }_{i_d}^d):i_1,i_2,\dots ,i_d,\mathrm{traverse}\mathrm{from}1\mathrm{to}\mathrm{N}\}$. And the eigenfunction $\varphi$ is approximated by

$$\varphi (x_1,\dots ,x_d)\approx {\varphi }^N(x_1,\dots ,x_d)=\sum _{i_1,\dots ,i_d=1}^N{\varphi }^N({\xi }_{i_1}^1,{\xi }_{i_2}^2,\dots ,{\xi }_{i_d}^d)L_{i_1}\left(x_1\right)\dots L_{i_d}\left(x_d\right)$$

(29)

Correspondingly, ${\varphi }^N(\mathbf{\Theta })$ is a tensor of eigenfunction values on collocation points, and $D_1,D_2,\dots D_d$ means differentiation matrices. With n-mode multiplication in tensor algebra, we can write

$$\mathcal{K}{\varphi }^N(\mathbf{\Theta })=\sum _{i=1}^d{f}_i(\mathbf{\Theta })\odot \left({\varphi }^N(\mathbf{\Theta }){\times }_i{D}_i\right)$$

(30)

Pursuing the analogous vectorization approach in Equation (27), we can reformulate Equation (30) as

$$\begin{array}{cc}\hfill \mathcal{K}vec\left({\varphi }^N(\mathbf{\Theta })\right)& =\{\sum _{i=1}^{d}vec\left(f_i(\mathbf{\Theta })\right)\odot \left({\otimes }_i{D}_i\right)\}vec\left({\varphi }^N(\mathbf{\Theta })\right)\hfill \\ \hfill & =\left\{\sum _{i=1}^{d}diag\left(vec\left(f_i(\mathbf{\Theta })\right)\right)\left({\otimes }_i{D}_i\right)\right\}vec\left({\varphi }^N(\mathbf{\Theta })\right)\hfill \end{array}$$

(31)

Consequently,

$$K=\sum _{i=1}^{d}diag\left(vec\left(f_i(\mathbf{\Theta })\right)\right)\left({\otimes }_i{D}_i\right)$$

(32)

where ${\otimes }_i{D}_i=I\otimes I\otimes I\otimes \cdots \otimes D_i\otimes \cdots \otimes I$ with $D_i$ at the $(d-i+1)$-position. To be more intuitive, here is an example when $d=3$:

$$\begin{array}{cc}\hfill K=& diag\left\{vec\left(F_1\right)\right\}I\otimes I\otimes D_1\hfill \\ \hfill +& diag\left\{vec\left(F_2\right)\right\}I\otimes D_2\otimes I\hfill \\ \hfill +& diag\left\{vec\left(F_3\right)\right\}{D}_3\otimes I\otimes I\hfill \end{array}$$

(33)

3.2. Solution of the Linear System

Now, let us 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

$$g\left(\mathit{x}\left(t\right)\right)\approx g^N\left(\mathit{x}\left(t\right)\right)=\sum _{j=1}^N{\widehat{c}}_j{\varphi }_j^N\left(\mathit{x}\left(t_0\right)\right)e^{{\widehat{\lambda }}_j}$$

(34)

When $d=1$, for a scalar observable g, consider Gauss–Lobatto points on a region with radius r around $x_0$, where $\Theta ={\left\{{\xi }_i\right\}}_{i=1}^N$ and ${\xi }_1<{\xi }_2<\cdots <{\xi }_N$. Here, we pick the odd number N such that ${\xi }_{(N+1)/2}=x_0=x\left(t_0\right)$. Thus, all Gauss–Lobatto points are on the interval $[x_0-r,x_0+r]$. Suppose $({\widehat{\lambda }}_j,v_j)$ is an eigenpair of matrix K. Vector $v_j$ are used to approximate eigenfunction ${\varphi }_j^N$ evaluated at the collocation points. We have ${\left(v_j\right)}_i={\varphi }_j^N\left({\xi }_i\right)\approx {\varphi }_j\left({\xi }_i\right)$. Hence,

$$g^N\left(x\left(t\right)\right)=\sum _{j=1}^N{\widehat{c}}_j{\varphi }_j^N\left(x\left(t_0\right)\right)e^{{\widehat{\lambda }}_{j}t}=\sum _{j=1}^N{\widehat{c}}_j{\varphi }_j^N\left({\xi }_{(N+1)/2}\right)e^{{\widehat{\lambda }}_{j}t}=\sum _{j=1}^N{\widehat{c}}_j{\left(v_j\right)}_{\frac{N+1}{2}}{e}^{{\widehat{\lambda }}_{j}t}$$

(35)

Now consider the N-dimensional linear system below

$$\frac{dy}{dt}=Ky{y}_0={(g(x_0-r),\dots ,g(x_0+r))}^T$$

(36)

We can write the analytical solution

$$y\left(t\right)=e^{Kt}{y}_0$$

(37)

Suppose that the eigendecomposition of $K=V\Lambda V^{-1}$, where each column of V is an eigenvector of K denoted as $v_j$, and $\Lambda$ contains all the eigenvalues. Further, by setting $t=0$ in Equation (35), we have

$$g_N\left(x_0\right)=\sum _{j=1}^N{\widehat{c}}_j{\varphi }_j^N\left(x_0\right)$$

(38)

This is also true for different initial values. Thus,

$$g_N\left({\xi }_i\right)=\sum _{j=1}^N{\widehat{c}}_j{\varphi }_j^N\left({\xi }_i\right)=\sum _{j=1}^N{\widehat{c}}_j{\left(v_j\right)}_i,i=1,2,\dots ,N$$

(39)

Thus, we have $V^{-1}{y}_0=c$ since $y\left(0\right)={(g^N\left({\xi }_1\right),g^N\left({\xi }_2\right),\dots ,g^N\left({\xi }_N\right))}^T$, where c is called the Koopman mode:

$$\begin{array}{cc}\hfill y\left(t\right)& =e^{Kt}{y}_0\hfill \\ \hfill & =e^{V\Lambda t{V}^{-1}}{y}_0\hfill \\ \hfill & =V{e}^{\Lambda t}{V}^{-1}{y}_0\hfill \\ \hfill & =\left[\begin{array}{c}{v}_1,v_2,\dots ,v_N\end{array}\right]\left[\begin{array}{cccc}{e}^{{\lambda }_{1}t}& & & \\ & e^{{\lambda }_{2}t}& & \\ & & \ddots & \\ & & & e^{{\lambda }_{N}t}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_N\end{array}\right]\hfill \\ \hfill & =\sum _{j=1}^N{c}_j{e}^{{\lambda }_{j}t}{v}_j\hfill \end{array}$$

(40)

And $g^N\left(x\left(t\right)\right)={\sum }_{j=1}^N{c}_j{e}^{{\lambda }_{j}t}{\left(v_j\right)}_{\frac{N}{2}}$.

Notice that $c_j$ are Koopman modes, and $e^{{\lambda }_{j}t}$ are eigenfunctions. In the general case, for the d-dimensional nonlinear system, suppose $\mathbf{\Theta }$ represents d-dimensional collocation points, and K is constructed as Equation (32). Now consider the linear system given below

$$\frac{d\mathit{y}\left(t\right)}{dt}=K\mathit{y}\mathit{y}\left(0\right)=vec\left(g(\mathbf{\Theta })\right).$$

(41)

Again, by vectorization, ${\varphi }_j\left({\mathit{x}}_{\mathbf{0}}\right)$ can be approximated by the “middle term” of tensor ${\varphi }^N(\mathbf{\Theta })$, which leads to

$$y\left(t\right)=\sum _{j=1}^{N^d}{c}_j{e}^{{\lambda }_{j}t}\left(v_j\right).$$

(42)

Only the middle term is needed for this long vector to compute the observable:

$$g_N\left(x\left(t\right)\right)=\sum _{j=1}^{N^d}{c}_j{e}^{{\lambda }_{j}t}{\left(v_j\right)}_{\frac{N^d+1}{2}},$$

(43)

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 we need to compute multi-dimensional 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\left(\mathit{x}\right)=\mathit{x}=(g_1\left(\mathit{x}\right),g_2\left(\mathit{x}\right),g_3\left(\mathit{x}\right))$, we would need to utilize three different initial conditions: $g_1(\mathbf{\Theta }),g_2(\mathbf{\Theta }),g_3(\mathbf{\Theta })$ 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 $\mathcal{A}$ shall be truncated at a certain level N, which yields the ${\mathcal{A}}_{\mathit{N}}$ matrix as follows:

$${\mathcal{A}}_{\mathit{N}}=\left[\begin{array}{ccccccccc} & A_1^1& A_2^1& A_3^1& \cdots & A_k^1& 0& 0& \cdots \\ & 0& A_2^2& A_3^2& \cdots & A_k^2& A_{k+1}^2& 0& \cdots \\ & 0& 0& A_3^3& \cdots & A_k^3& A_{k+1}^3& A_{k+2}^3& \cdots \\ & ⋮& ⋮& ⋮& & ⋮& ⋮& ⋮& \\ & 0& 0& 0& \dots .& 0& 0& 0& A_N^N\end{array}\right].$$

(44)

Then, the infinite linear system can be approximated by

$$\frac{dy}{dt}={\mathcal{A}}_{\mathit{N}}y$$

(45)

and the approximated solution of an original nonlinear system is exactly the first d elements of the solution for Equation (45).

4. Numerical Result

This section presents the performance comparison of Koopman Spectral Linearization against Carleman Linearization on five nonlinear dynamic models in Section 4.1. In each example, we visualized the curves of specific solution components and investigated the influence of truncation order N on accuracy. As N increases, both the precision of linearization and the computational effort increase correspondingly. In practical applications, the selection of parameter N should be based on the required accuracy and available computational resources. Since Carleman Linearization can only be applied to polynomial cases, all the following nonpolynomial examples are converted to polynomial form using higher-order Taylor expansions before computation. The reference solution is generated by The 4th order Runge–Kutta method (RK4) is used if no closed-form solution is available. We also aim to compare the loss of accuracy in the linearization step under the assumption of an ideal quantum computing environment, where a high-precision quantum linear solver or other quantum algorithms for linear ODEs are available. Hence, the solutions of the resulting linear ODEs are directly computed using the exp function in MATLAB 2022a. Next, in Section 4.2, we further compare the matrix size and computational cost (all the MATLAB codes can be downloaded at https://github.com/DongDongShiShi/Koopman-Comparison (accessed on 30 May 2024).

4.1. Linearize Dynamics with Koopman Spectral Linearization

4.1.1. Quadratic Model

We start with a simple nonlinear ODE, and the governing ODE is given by

$$\frac{dx}{dt}=x^2$$

We set $x\left(0\right)=0.08$ and T = 10 in this example. This quadratic model has a closed form solution $x\left(t\right)=\frac{1}{(1/x_0)-t}$. In the tests, we fix $r=0.03$. Figure 2 summarizes the results of the first comparison. Figure 2b shows the exponential convergence of our approach concerning the truncation order, which is similar to the conclusion in the conventional Carleman Linearization method. Further, our approach shows a faster convergence rate.

4.1.2. Cosine Square Model

This cosine square model is a synthetic model invented as a nonpolynomial case for our demonstrative purposes. The governing equation is given as

$$\frac{dx}{dt}={cos}^2\left(x\right).$$

And we set $x\left(0\right)=0.1$, $T=10$. Despite the nonlinear nature, we still have a closed-form solution $x\left(t\right)=arctan(-0.5t+tan(x_0\left)\right)$. We fix $r=0.3$ in the tests. It is important to note that as the system evolves over time, errors gradually accumulate. This is particularly evident when using lower-order polynomial interpolation (i.e., $N=3$). Specifically, in Figure 3a, beyond $t=8$, the curve begins to exhibit a more pronounced deviation from the standard solution. Figure 3b shows the exponential convergence of our approach concerning the truncation order, which is still similar to the conclusion in the previous example. Conversely, for this nonpolynomial model, the accuracy of Carleman Linearization is very low, and we did not observe any significant changes in accuracy across the different values of N.

4.1.3. Simple Pendulum Model

The simple pendulum is a well-studied model in science and engineering. The movement of the pendulum is described by the following second-order ordinary different equation:

$$\frac{d^2\theta }{dt}=-\frac{\mathrm{g}}{\mathrm{L}}sin\left(\theta \right),$$

where L is the length of the pendulum, $\theta$ 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=\theta$ and $x_2=\frac{d\theta }{dt}$. Also, we set $L=g=9.8$ for simplicity. Correspondingly,

$$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =x_2,\hfill \\ \hfill \frac{d{x}_2}{dt}& =-sin\left(x_1\right),\hfill \end{array}$$

and we set $x\left(0\right)={(0.1,0.1)}^T$, $T=5$. We fix $r=1$ in the tests.

Figure 4 presents the numerical result of a simple pendulum. Similar to the previous result, the exponential convergence of error concerning truncation order N is observed from Figure 4b. Although the accuracy of Carleman Linearization for this nonpolynomial model is significantly higher compared to its application on the cosine square model and even more accurate than Koopman Spectral Linearization at $N=3$ and $N=5$, the accuracy does not significantly change as N increases.

4.1.4. Lotka–Volterra Model

The Lotka–Volterra model, also known as the predator–prey model [40], 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: prey and predator populations. Here, we define as below:

$$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =1.1x_1-0.4x_1{x}_2,\hfill \\ \hfill \frac{d{x}_2}{dt}& =0.1x_1{x}_2-0.4x_2.\hfill \end{array}$$

We set $x\left(0\right)={(5,5)}^T$, $T=10$ and $r=3$.

Figure 5 presents the results of these tests for the Lotka–Volterra model. Again, errors are decreased exponentially with respect to the truncation order. For this polynomial model, the accuracy of Carleman Linearization improves with increasing N. However, regarding both numerical accuracy and convergence speed, Koopman Spectral Linearization shows a significant advantage.

4.1.5. Kraichnan–Orszag Model

The Kraichnan–Orszag model was presented in [41] for modeling fluid dynamics. This three-dimensional nonlinear model is given by

$$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =x_2{x}_3,\hfill \\ \hfill \frac{d{x}_2}{dt}& =x_1{x}_3,\hfill \\ \hfill \frac{d{x}_3}{dt}& =-2x_1{x}_2.\hfill \end{array}$$

We set $x\left(0\right)=(0.1,$−${0.2,0.3)}^T$, $T=5$ and $r=0.1$.

Figure 6 presents the results of the Kraichnan–Orszag model. Different from the previous example, a shorter time period is set due to the strong oscillations of the Kraichnan–Orszag model.

The truncation order testing figures illustrate the error as a function of the truncation order N and provide a comparison between Carleman Linearization and our method. The mistake of Carleman Linearization in polynomial cases (i.e., Section 4.1.1, Section 4.1.4 and Section 4.1.5) exhibits exponential decay, a phenomenon that has been substantiated by recent theoretical research [42,43]. Similarly, 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 Section 4.1.1, Section 4.1.4 and Section 4.1.5, 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 increasing. Moreover, when applied to nonpolynomial dynamical systems, our method showcases superior numerical accuracy and faster convergence rates. In Section 4.1.2 and Section 4.1.3, we present a comparative analysis of the error–truncation order relationship in nonpolynomial dynamical systems. Both nonpolynomial dynamics are subjected to 12th-order Taylor expansions to transform them into polynomial forms. The Carleman method exhibits significantly lower accuracy in this scenario, with notably slower convergence. This limitation arises because in nonpolynomial 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 nonpolynomial 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.

4.2. Matrix Size and Computational Cost Comparison

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\times N^d$. However, we can reuse the lifted matrix with different initial values for multi-dimensional dynamics. In the case of the Carleman Linearization procedure, the matrix size will be ${\sum }_{i=1}^N{d}^i\times {\sum }_{i=1}^N{d}^i$, and we can approximate the size of the Carleman lifted matrix as $O\left(d^N\right)$. 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 than the matrix in Carleman Linearization. To further illustrate this point, Figure 7 compares our approach with the size of the Carleman Linearization matrix.

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

Table 1. Error and time cost compared with Carleman Linearization ($N=9$).

	Error		Time (s)
	Carleman	Koopman	Carleman	Koopman
Quadratic	$4.19\times {10}^{-3}$	$2.56\times {10}^{-5}$	0.0147	0.0089
Lotka-Volterra	$1.49\times {10}^{-1}$	$6.39\times {10}^{-6}$	0.1525	0.0188
Kraichnan-Orszag	$8.80\times {10}^{-3}$	$1.99\times {10}^{-7}$	32.0546	0.1370
Cosine Square	1.26	$7.39\times {10}^{-4}$	0.0059	0.0049
Simple Pendulum	$1.20\times {10}^{-3}$	$7.02\times {10}^{-5}$	0.0368	0.0141

for a truncation order of N = 9. 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 s, being considerably slower than our approach.

On the other hand, 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. 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 and surpassing the precision of Carleman Linearization remains unattainable.

5. Discussion and Conclusions

The Koopman Spectral Linearization method explicitly uses the differentiation matrix to derive the lifted matrix of nonlinear autonomous dynamical systems. It provides a finite-dimensional representation for the generator of the Koopman Operator with a polynomial basis. Therefore, similar to Carleman Linearization, our approach offers a linearization method for nonlinear autonomous ordinary differential equations. 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 nonpolynomial. 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 has some limitations. It is essential 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. We employ interpolation within a local space of radius r around the initial point $x_0$ to approximate the eigenfunctions across in a local domain. However, since the spatial range of the system evolution during the process is not known in advance, it is very challenging to select the value of r without additional information; when theoretical guidance is unavailable, r is chosen empirically as a relatively small value (i.e., 0.1–0.5) to ensure a higher precision of approximation. If the spatial range of the system’s evolution is known beforehand, such as the case of fluid moving in a pipeline, r can then be selected based on spatial constraints. To further understand how this parameter r influences the accuracy, a comprehensive numerical analysis will be included in future work. In addition, the global interpolation method might overcome the limitation of the estimate r. Besides concerns regarding the parameter r, the degree of nonlinearity in physical systems can vary significantly, which might also influence the performance of our linearization method. Particularly for highly nonlinear systems such as turbulence, no linearization technique is capable of providing high-precision linear approximations over extended periods and large areas. In such instances, data-driven approaches may offer more satisfactory solutions, addressing the complexities inherent in these systems more effectively.

Regarding the matrix size as indicated in Section 4.2, the Koopman Spectral Method obtains a much smaller matrix with even higher accuracy compared with Carleman Linearization. However, the matrix size is exponential increased with respect to dimension d (i.e., $N^d$) since the tensor product rule is applied. Therefore, both Carleman Linearization (i.e., $O\left(d^N\right)$) and our approach might not be efficient for the relatively high-dimensional problems on classic computers. A possible way to overcome this limitation is to adapt sparse grid methods to construct collocation points, which has been found successful in [30]. As pointed out in the Koopman Operator Theory, Carleman Linearization and our approach rely on tensor product rule; therefore, they need further modification to deal with high-dimensional problems on classical computers. However, on quantum computers, tensor products may not be a challenge, which may lead to innovative designs of algorithms that are completely different from their counterpart on classical computers. This is the main reason we conduct these comparisons.

As previous sections noted, the underlying principle (Schrödinger equation) of quantum computers is linear. Consequently, there are fundamental difficulties in using quantum computing to solve nonlinear ODEs/PDEs. Currently, the linearization of differential equations appears to be an unavoidable step in addressing this challenge. In particular, Carleman Linearization as a solution strategy has garnered significant attention in recent research [17]. By reviewing the underlying theory of Carleman Linearization—Koopman Operator Theory—and integrating it with the concept of differentiation matrices in spectral methods, we propose the Koopman Spectral Linearization method. Our novel linearization approach might be useful when designing quantum algorithms for nonlinear ODEs/PDEs under certain situations.