Carleman linearization

Examples

The following examples illustrate the construction of Type A and Type B Carleman approximants in one and two dimensions. Only one equilibrium point is considered in each Type A construction.

One-dimensional example: Riccati equation

Consider the Riccati equation

${\displaystyle \dot{u}=-u+u^2,u(0)=c.}$

Its exact solution is

${\displaystyle u(t)=\frac{c}{c+(1-c)e^t}.}$

The equation has two equilibrium points, ${\displaystyle u=0}$ and ${\displaystyle u=1}$. The Type A approximant below is constructed only around ${\displaystyle u=0}$.

Type A approximant

Introduce the Carleman variables

${\displaystyle u_k(t)=u^k(t),k=1,2,\dots .}$

Differentiation gives the infinite linear sequence

${\displaystyle {\dot{u}}_k=-k{u}_k+k{u}_{k+1},k=1,2,\dots .}$

A second-order approximation is obtained by retaining ${\displaystyle u_1=u}$ and ${\displaystyle u_2=u^2}$ and discarding the term ${\displaystyle u_3}$. The resulting homogeneous system is

${\displaystyle \frac{d}{dt}\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{cc}-1& 1\\ 0& -2\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right),\left(\begin{array}{c}{u}_1(0)\\ u_2(0)\end{array}\right)=\left(\begin{array}{c}c\\ c^2\end{array}\right).}$

Solving the linear system gives the second-order Type A approximant

${\displaystyle {\tilde{u}}_{N=2}^A(t)=u_1(t)=c\left[1+c\left(1-e^{-t}\right)\right]e^{-t}.}$

Type B approximant

To construct the approximant around the initial condition, introduce

${\displaystyle u=w+c,w(0)=0.}$

The differential equation becomes

${\displaystyle \dot{w}=(2c-1)w+w^2+c(c-1).}$

Define

${\displaystyle \alpha =2c-1,\beta =c(c-1).}$

Retaining ${\displaystyle w}$ and ${\displaystyle w^2}$ gives the second-order non-homogeneous linear system

${\displaystyle \frac{d}{dt}\left(\begin{array}{c}w\\ w^2\end{array}\right)=\left(\begin{array}{cc}\alpha & 1\\ 2\beta & 2\alpha \end{array}\right)\left(\begin{array}{c}w\\ w^2\end{array}\right)+\left(\begin{array}{c}\beta \\ 0\end{array}\right),\left(\begin{array}{c}w(0)\\ w^2(0)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).}$

Writing

${\displaystyle M=\left(\begin{array}{cc}\alpha & 1\\ 2\beta & 2\alpha \end{array}\right),\mathbf{b}=\left(\begin{array}{c}\beta \\ 0\end{array}\right),}$

the solution is

${\displaystyle \left(\begin{array}{c}w(t)\\ w^2(t)\end{array}\right)=M^{-1}\left(e^{Mt}-I\right)\mathbf{b},}$

when ${\displaystyle M}$ is invertible. The second-order Type B approximant is therefore

${\displaystyle \begin{array}{l}{\tilde{u}}_{N=2}^B(t)=c+w(t)=\\ {\displaystyle c\{1+\frac{c-1}{3c(c-1)+1}[1-2c+\frac{1}{2\xi }(\left[(2c-1)+6c(c-1)+1\right]e^{{\lambda }_{-}t}}\\ +\left[(2c-1)-6c(c-1)-1\right]e^{{\lambda }_{+}t})]\}\end{array}}$

where ${\displaystyle \xi \equiv \sqrt{1+12c(c-1)}}$, and ${\displaystyle {\lambda }_{\pm }=[3(2c-1)\pm \xi ]/2}$ are the eigenvalues of ${\displaystyle M}$. Although this expression is algebraically more involved than the Type A approximant, centring the construction at ${\displaystyle u(0)=c}$ generally improves its local accuracy.

The plot compares the exact solutions with second- and third-order Type A approximants constructed around the two equilibrium points and with a second-order Type B approximant constructed around each initial condition.

Two-dimensional example: center

Consider the two-dimensional nonlinear system

${\displaystyle \dot{x}=-y-xy,\dot{y}=x+xy,}$

with initial condition

${\displaystyle x(0)=a,y(0)=b.}$

The system has an equilibrium at ${\displaystyle (0,0)}$, which is a center and is stable in the sense of Lyapunov, but not asymptotically stable. It also has a saddle point at ${\displaystyle (-1,-1)}$. For ${\displaystyle x,y>-1}$, its trajectories satisfy the first integral

${\displaystyle x-\ln (1+x)+y-\ln (1+y)=C,}$

where ${\displaystyle C}$ is constant along each trajectory.

Type A approximant

For the Type A construction, consider only the equilibrium point ${\displaystyle (0,0)}$. At second order, introduce the enlarged set of variables

${\displaystyle \mathbf{X}=\left(\begin{array}{c}y\\ y^2\\ x\\ xy\\ x^2\end{array}\right).}$

After discarding terms of degree greater than two, the Carleman system is

${\displaystyle \frac{d}{dt}\left(\begin{array}{c}y\\ y^2\\ x\\ xy\\ x^2\end{array}\right)=\left(\begin{array}{ccccc}0& 0& 1& 1& 0\\ 0& 0& 0& 2& 0\\ -1& 0& 0& -1& 0\\ 0& -1& 0& 0& 1\\ 0& 0& 0& -2& 0\end{array}\right)\left(\begin{array}{c}y\\ y^2\\ x\\ xy\\ x^2\end{array}\right),}$

with initial condition

${\displaystyle \mathbf{X}(0)=\left(\begin{array}{c}b\\ b^2\\ a\\ ab\\ a^2\end{array}\right).}$

Thus,

${\displaystyle \mathbf{X}(t)=e^{Mt}\mathbf{X}(0),}$

where ${\displaystyle M}$ is the displayed coefficient matrix. The second-order Type A approximants ${\displaystyle {\tilde{x}}_{N=2}^A(t)}$ and ${\displaystyle {\tilde{y}}_{N=2}^A(t)}$ are respectively the third and first components of ${\displaystyle \mathbf{X}(t)}$.

Type B approximant

To centre the approximation at the initial condition ${\displaystyle (a,b)}$, introduce

${\displaystyle w=x-a,z=y-b,w(0)=z(0)=0.}$

The transformed nonlinear system is

${\displaystyle \dot{w}=-bw-(1+a)z-wz-b(1+a),}$

${\displaystyle \dot{z}=(1+b)w+az+wz+a(1+b).}$

At first Carleman order, the nonlinear term ${\displaystyle wz}$ is discarded, yielding

${\displaystyle \frac{d}{dt}\left(\begin{array}{c}w\\ z\end{array}\right)=\left(\begin{array}{cc}-b& -(1+a)\\ 1+b& a\end{array}\right)\left(\begin{array}{c}w\\ z\end{array}\right)+\left(\begin{array}{c}-b(1+a)\\ a(1+b)\end{array}\right).}$

Let

${\displaystyle M=\left(\begin{array}{cc}-b& -(1+a)\\ 1+b& a\end{array}\right),\mathbf{h}=\left(\begin{array}{c}-b(1+a)\\ a(1+b)\end{array}\right).}$

The first-order Type B approximant is

${\displaystyle \left(\begin{array}{c}{\tilde{x}}_{N=1}^B(t)\\ {\tilde{y}}_{N=1}^B(t)\end{array}\right)=\left(\begin{array}{c}a\\ b\end{array}\right)+M^{-1}\left(e^{Mt}-I\right)\mathbf{h},}$

when ${\displaystyle M}$ is invertible. The eigenvalues of the coefficient matrix are

${\displaystyle {\lambda }_{\pm }=\frac{1}{2}\left[a-b\pm \sqrt{(a-b{)}^2-4(1+a+b)}\right].}$

Despite being only first order, the Type B approximation can have a numerical accuracy intermediate between the first- and second-order Type A approximants because it is centred directly on the trajectory's initial condition.

The phase portrait compares second-order Type A approximants constructed around the two critical points with first-order Type B approximants constructed around several initial conditions.

Application to epidemiology

The Carleman approximants for the SIR model yield an explicit algebraic estimator of the time-dependent effective reproduction number ${\displaystyle R_t}$ from prevalence data. This application illustrates how an analytic approximation of a nonlinear system can be used not only to approximate its trajectories, but also to obtain closed-form expressions for quantities of practical interest.

Consider the SIR model

${\displaystyle \dot{S}=-\beta SI,\dot{I}=\beta SI-\xi I,\dot{R}=\xi I,}$

where ${\displaystyle S}$, ${\displaystyle I}$, and ${\displaystyle R}$ are the susceptible, infectious, and removed population fractions, respectively. The transmission and removal rates are denoted by ${\displaystyle \beta }$ and ${\displaystyle \xi }$. The effective reproduction number is

${\displaystyle R_t=S(t)\frac{\beta }{\xi }.}$

To account explicitly for the approximately exponential behaviour of the infectious population, the change of variable

${\displaystyle I=e^Y}$

transforms the first two SIR equations into the system

${\displaystyle \dot{S}=-\beta S{e}^Y,\dot{Y}=\beta S-\xi .}$

This formulation is referred to as the SYR system. Its Carleman sequence can be constructed using nonlinear terms of the form ${\displaystyle S^i{e}^{jY}}$, with approximation order ${\displaystyle N=i+j}$.

Let

${\displaystyle S(0)=a,I(0)=b,q=\frac{b}{a},R_q=a\frac{\beta }{\xi },}$

where ${\displaystyle R_q}$ is the value of the effective reproduction number at the time chosen as the origin. At second order, the Carleman approximant for the logarithmic prevalence is

${\displaystyle Y_{N=2}(t)=\ln b+(R_q-1)\xi t+q{R}_q^2\left(1-\xi t-e^{-\xi t}\right).}$

Because this approximation is an explicit analytic function of the epidemiological parameters, it can be fitted directly to prevalence data. Assume that the removal rate ${\displaystyle \xi }$ is known and remains approximately constant within a short time interval, and that prevalence ${\displaystyle I}$ has been reconstructed from incidence data. For observations indexed by ${\displaystyle i=-T,\dots ,T}$, define

${\displaystyle {\tau }_i=\xi t_i,z_i=1-{\tau }_i-e^{-{\tau }_i},y_i=\ln I(t_i).}$

The second-order approximation is fitted over a time window containing ${\displaystyle 2T+1}$ observations by minimizing the residual sum of squares

${\displaystyle \mathrm{RSS}=\frac{1}{2}\sum _{i=-T}^T{\left[\ln b+(R_q-1){\tau }_i+q{R}_q^2{z}_i-y_i\right]}^2}$

with respect to the parameters ${\displaystyle \ln b}$, ${\displaystyle R_q}$, and ${\displaystyle q}$. The parameters are assumed to remain approximately constant within the selected time window. The resulting value of ${\displaystyle R_q}$ estimates the effective reproduction number ${\displaystyle R_t}$ at the central point of the interval.

Although the resulting variational equations are nonlinear, their structure allows them to be solved in closed form. Define the centred quantities

${\displaystyle {\tilde{\tau }}_i={\tau }_i-\overline{\tau },{\tilde{z}}_i=z_i-\overline{z},{\tilde{y}}_i=y_i-\overline{y},}$

where

${\displaystyle \overline{\tau }=\frac{1}{2T+1}\sum _{i=-T}^T{\tau }_i,\overline{z}=\frac{1}{2T+1}\sum _{i=-T}^T{z}_i,\overline{y}=\frac{1}{2T+1}\sum _{i=-T}^T{y}_i.}$

Minimization of the residual sum of squares then gives the algebraic estimator

${\displaystyle R_t=1-\frac{\left(\sum _{i=-T}^T{\tilde{y}}_i{z}_i\right)\left(\sum _{i=-T}^T{\tilde{z}}_i{\tau }_i\right)-\left(\sum _{i=-T}^T{\tilde{y}}_i{\tau }_i\right)\left(\sum _{i=-T}^T{\tilde{z}}_i{z}_i\right)}{\left(\sum _{i=-T}^T{\tilde{\tau }}_i{z}_i\right)\left(\sum _{i=-T}^T{\tilde{z}}_i{\tau }_i\right)-\left(\sum _{i=-T}^T{\tilde{\tau }}_i{\tau }_i\right)\left(\sum _{i=-T}^T{\tilde{z}}_i{z}_i\right)}.}$

Thus, once the prevalence and removal time scale have been specified, ${\displaystyle R_t}$ can be estimated from finite sums over the data. The expression has been tested using simulated epidemic waves with time-dependent transmission rates and SARS-CoV-2 incidence data from the Valencian Country, Spain, yielding estimates consistent with the corresponding reference values.

Control-theoretic formulation

In control theory and state estimation, Carleman linearization can be applied to nonlinear systems with external inputs or disturbances. One such formulation considers the system

${\displaystyle \dot{x}=f(x)+\sum _{j=1}^r{g}_j(x)d_j(t),}$

where ${\displaystyle x\in {\mathbb{R}}^n}$ is the state vector, ${\displaystyle f}$ and ${\displaystyle g_j}$ are analytic vector fields, and ${\displaystyle d_j(t)}$ denotes the ${\displaystyle j}$-th component of an external input or disturbance. Because the disturbances may depend explicitly on time, the system is not autonomous in general. When the functions ${\displaystyle d_j(t)}$ are treated as inputs, the lifted model obtained below is bilinear in the lifted state and the inputs.

Let ${\displaystyle x_0}$ be a nominal point. After introducing the shifted coordinate ${\displaystyle z=x-x_0}$, the nominal point is placed at the origin. For notational simplicity, the shifted variable is again denoted by ${\displaystyle x}$. The vector fields may then be approximated by finite Taylor expansions,

${\displaystyle f(x)\simeq \sum _{k=0}^{\eta }{A}_k{x}^{[k]},g_j(x)\simeq \sum _{k=0}^{\eta }{B}_{j,k}{x}^{[k]},}$

where ${\displaystyle \eta }$ is the Taylor truncation degree and ${\displaystyle A_k}$ and ${\displaystyle B_{j,k}}$ are the corresponding matricized Taylor coefficient tensors,

${\displaystyle A_k=\frac{1}{k!}{D}^{k}f(0),B_{j,k}=\frac{1}{k!}{D}^k{g}_j(0).}$

Here, ${\displaystyle x^{[k]}}$ denotes the ${\displaystyle k}$-th Kronecker power,

${\displaystyle x^{[k]}=\underset{k\text{ factors}}{\underbrace{x\otimes x\otimes \cdots \otimes x}},}$

with the convention ${\displaystyle x^{[0]}=1}$.

The time derivative of the ${\displaystyle i}$-th Kronecker power is

${\displaystyle \frac{d{x}^{[i]}}{dt}=\sum _{l=0}^{i-1}{x}^{[l]}\otimes \dot{x}\otimes x^{[i-1-l]}.}$

Substitution of the Taylor expansions generates Kronecker powers of progressively higher degree. To obtain a finite-dimensional approximation, the hierarchy is truncated at a maximum lifted degree ${\displaystyle N}$. Discarding terms whose total degree exceeds ${\displaystyle N}$ gives, for ${\displaystyle i=1,\dots ,N}$,

${\displaystyle \frac{d{x}^{[i]}}{dt}\simeq \sum _{k=0}^{min(\eta ,N-i+1)}{A}_{i,k}{x}^{[i+k-1]}+\sum _{j=1}^r\sum _{k=0}^{min(\eta ,N-i+1)}{B}_{j,i,k}{x}^{[i+k-1]}{d}_j(t),}$

where

${\displaystyle A_{i,k}=\sum _{l=0}^{i-1}{I}_n^{[l]}\otimes A_k\otimes I_n^{[i-1-l]},}$

and similarly,

${\displaystyle B_{j,i,k}=\sum _{l=0}^{i-1}{I}_n^{[l]}\otimes B_{j,k}\otimes I_n^{[i-1-l]}.}$

In these expressions, ${\displaystyle I_n}$ is the ${\displaystyle n}$-dimensional identity matrix, ${\displaystyle I_n^{[l]}}$ denotes its ${\displaystyle l}$-th Kronecker power, and ${\displaystyle I_n^{[0]}=1}$.

Introducing the truncated lifted state

${\displaystyle x_{\otimes }={\left[\begin{array}{cccc}{x}^T& {\left(x^{[2]}\right)}^T& \cdots & {\left(x^{[N]}\right)}^T\end{array}\right]}^T,}$

the system can be written in block-matrix form as

${\displaystyle {\dot{x}}_{\otimes }=\mathcal{A}{x}_{\otimes }+{\mathcal{a}}_0+\sum _{j=1}^r\left({\mathcal{B}}_j{x}_{\otimes }+{\mathcal{b}}_{j,0}\right)d_j(t)+r_N(x,d),}$

where ${\displaystyle \mathcal{A}}$ and ${\displaystyle {\mathcal{B}}_j}$ are assembled from the blocks ${\displaystyle A_{i,k}}$ and ${\displaystyle B_{j,i,k}}$, respectively. The vectors ${\displaystyle {\mathcal{a}}_0}$ and ${\displaystyle {\mathcal{b}}_{j,0}}$ contain the constant terms of the Taylor expansions, while ${\displaystyle r_N(x,d)}$ represents the residual produced by the Taylor and Carleman truncations.

If the disturbances ${\displaystyle d_j(t)}$ are prescribed functions of time and the residual is neglected, the truncated model is a finite-dimensional affine linear time-varying system. When the disturbances are regarded as independent inputs, it is a bilinear state–input system. This formulation has been used, for example, in moving-horizon estimation of nonlinear processes.

Convergence and limitations

Except in particular systems for which the Carleman sequence closes after finitely many steps, a finite-order Carleman system provides an approximation rather than an exact representation of the original nonlinear dynamics. The accuracy depends on the order of truncation, the point in phase space around which the approximant is constructed, the initial condition, and the time interval considered.

The dimension of the enlarged linear system grows rapidly with both the number of original variables and the approximation order. Consequently, the computation and analytic exponentiation of the coefficient matrix may become difficult at high orders. Type B approximants may offer better local accuracy, but require the solution of a non-homogeneous system and generally have more complicated algebraic expressions than Type A approximants.

Convergence and truncation-error estimates have been established for particular classes of nonlinear systems, but the applicable bounds and time intervals depend on assumptions about the vector field and the solution.

Because finite-order Carleman approximants are solutions of linear differential systems with constant coefficients, they consist of linear combinations of functions such as ${\displaystyle t^k{e}^{\lambda t}}$, ${\displaystyle t^k\sin (\omega t)}$, and ${\displaystyle t^k\cos (\omega t)}$. They therefore do not reproduce finite singularities of an exact nonlinear solution in the complex time plane. They may also fail to preserve global dynamical structures such as closed orbits; at higher orders, repeated eigenvalues can generate secular terms that prevent an approximate trajectory from closing.