# A Convex Optimization Approach to Learning Koopman Operators

M. Sznaier The author is with the Department of Electrical and Computer
Engineering, Northeastern University, Boston, MA 02115, email
<EMAIL_ADDRESS>This work was partially supported by NSF grants
CNS–1646121, CMMI–1638234, IIS–1814631 and ECCS 1808381 and AFOSR grant
FA9550-19-1-0005.

###### Abstract

Koopman operators provide tractable means of learning linear approximations of
non-linear dynamics. Many approaches have been proposed to find these
operators, typically based upon approximations using an a-priori fixed class
of models. However, choosing appropriate models and bounding the approximation
error is far from trivial. Motivated by these difficulties, in this paper we
propose an optimization based approach to learning Koopman operators from
data. Our results show that the Koopman operator, the associated Hilbert space
of observables and a suitable dictionary can be obtained by solving two rank-
constrained semi-definite programs (SDP). While in principle these problems
are NP-hard, the use of standard relaxations of rank leads to convex SDPs.

##  1 Introduction and motivation

Many scenarios involve predicting the output of an unknown non-linear system
based on past measurements and some a-priori information. Recently,
substantial interest has been devoted to the use of Koopman operator based
methods to solve this problem, as a tractable alternative to nonlinear
identification. An excellent introduction to the topic is given in [8], and
more recent references can be found in [7, 10]. Given a non-linear discrete
time system of the form:

$\begin{split}\boldsymbol{\xi}_{k+1}&=f(\boldsymbol{\xi}_{k})\;\text{where}\;\boldsymbol{\xi}_{k}=\begin{bmatrix}\mathbf{x}_{k-r+1}^{T}&\ldots&\mathbf{x}_{k}^{T}\end{bmatrix}^{T},\;\mathbf{x}_{j}\in\mathbb{R}^{n}\end{split}$
(1)

let $\mathbb{H}$ denote a Hilbert space of functions
$\boldsymbol{\psi}(\boldsymbol{\xi})\colon\mathbb{R}^{nr}\to\mathbb{R}^{mr}$(the
so called observables). The Koopman $\mathcal{K}$ operator acts on the
elements of $\mathbb{H}$, by propagating their values one step into the
future:

$(\mathcal{K}\circ\boldsymbol{\psi})(\boldsymbol{\xi}_{k})=(\boldsymbol{\psi}\circ
f)(\boldsymbol{\xi}_{k})=\boldsymbol{\psi}(\boldsymbol{\xi}_{k+1})$ (2)

$\mathcal{K}$ is a linear operator, albeit typically infinite dimensional.
When it has a countable set of eigenfunctions $\boldsymbol{\phi}_{i}(.)$ with
eigenvalues $\mu_{i}$, the observables $\boldsymbol{\psi}(.)$ can be
propagated as follows. Let
$\mathbf{a}=\begin{bmatrix}a_{1}\ldots\end{bmatrix}^{T}$ denote the
coordinates of $\boldsymbol{\psi}(.)$ in the basis spanned by
$\boldsymbol{\phi}(.)$, that is

$\boldsymbol{\psi}(.)=\sum
a_{i}\boldsymbol{\phi}_{i}(.)\doteq\boldsymbol{\Phi}(.)\mathbf{a},\;\text{where:
}\boldsymbol{\Phi}(.)=\begin{bmatrix}\boldsymbol{\phi}_{1}(.)\ldots\end{bmatrix}$

Then

$(\mathcal{K}\circ\boldsymbol{\psi})(.)=\sum
a_{i}\mu_{i}\boldsymbol{\phi}_{i}(.)=\boldsymbol{\Phi}(.)\mathbf{M}\mathbf{a},\text{where
$\mathbf{M}$ = diag($\mu_{i}$)}$

In particular, if the state
$\boldsymbol{\xi}\in\text{span}\\{\boldsymbol{\phi}_{i}\\}$, then
$\boldsymbol{\xi}_{k+1}=\boldsymbol{\Phi}(\boldsymbol{\xi}_{k})\mathbf{M}\mathbf{a}$.
While this approach leads for to linear representations of (1), identifying
the Koopman eigenfunctions from data is not trivial.

Extended Dynamical Mode Decomposition (EDMD) type approaches seek to identify
approximations to Koopman operators over a restricted subspace, defined by the
span of a given dictionary
$\mathcal{D}(.)\doteq\begin{bmatrix}\boldsymbol{\psi}_{1}(.)\ldots\boldsymbol{\psi}_{N}(.)\end{bmatrix}$.
In this subspace, the Koopman operator can then be approximated by a matrix
$\mathbf{K}\in\mathbb{R}^{N\times N}$ that propagates the coefficients of the
expansion, that is, for $\boldsymbol{\psi}(.)=\mathcal{D}(.)\mathbf{a}$, then
$(\mathcal{K}\circ\boldsymbol{\psi})(.)=\mathcal{D}(.)\mathbf{Ka}$. Typically,
given experimental data
$\mathbf{X}\doteq\begin{bmatrix}\boldsymbol{\xi}_{1}&\boldsymbol{\xi}_{2}&\ldots&\boldsymbol{\xi}_{T}\end{bmatrix}$,
$\mathbf{K}$ is found by minimizing the one-step prediction error over a set
of observables. Specifically, this approach considers $m$ observables
$\boldsymbol{\psi}^{(j)}(.)\doteq\mathcal{D}(.)\mathbf{a}_{j}$, each defined
by a coordinate vector $\mathbf{a}_{j}$, and solves:

$\mathbf{K}=\mathop{argmin}\limits_{K}\sum_{j=1}^{m}\sum_{k=1}^{T-1}\|[\mathbf{D}(\boldsymbol{\xi}_{k+1})-\mathbf{D}(\boldsymbol{\xi}_{k})\mathbf{K}]\mathbf{a}_{j}\|_{2}^{2}$
(3)

where $\mathbf{D}(\boldsymbol{\xi}_{k})$ is the matrix obtained by evaluating
the dictionary a the point $\boldsymbol{\xi}_{k}$. EDMD often works well, but
requires choosing a suitable dictionary, with the approximation error strongly
hinging on this choice. This approximation error can be reduced by considering
larger dictionaries, but this may lead to overfitting of the data and poor
generalization capabilities [10].

Deep learning motivated approaches use a neural network parameterized by a set
of weights $W$ as dictionary. The

Figure 1: Top: Finding Koopman operators via Semi-Definite Programs. The first
SDP (Section 3.3) finds the observables $\mathbf{y}_{k}$ corresponding to
given data $\mathbf{x}_{k}$, the Koopman operator $\mathbf{K}$, and the
Loewner matrices that encode the mapping $\mathbf{x}_{k}\to\mathbf{y}_{k}$.
The second SDP (Section 3.4) finds the inverse mapping
$\mathbf{y}_{k}\to\mathbf{x}_{k}$. Bottom: The pipeline to predict
$\mathbf{x}_{k+1}/\mathbf{x}_{k}\ldots,\mathbf{x}_{k-r}$ uses explicit
expressions for the predictions of the model $\mathbf{K}$. Thus, it only
requires $\mathcal{O}(r)$ operations.

Koopman operator $\mathbf{K}$ is found by alternatively minimizing the
prediction error over $W$ and $\mathbf{K}$. Alternating minimization methods
can get trapped in local minima. Further, the issue of which architectures are
best suited to represent dynamical systems is largely open. Recent work [7,
10] proposed encoder/decoder type architectures that map states
$\boldsymbol{\xi}$ to latent variables $\mathbf{y}$ and impose approximately
linear dynamics for the evolution of the latter. A salient feature of these
approaches is that the states $\boldsymbol{\xi}$ are no longer required to be
in the span of the Koopman eigenfuctions. As shown in [10], the use of a
nonlinear decoder to map $\mathbf{y}$ back to $\boldsymbol{\xi}$ (as opposed
to a linear one if $\boldsymbol{\xi}\in$ span { $\mathcal{D}$}) results in
substantially smaller dictionaries. Still, these methods require ad-hoc
parameter selection (dimension of the latent variables, order of the dynamics)
and, as before, can lead to local minima.

An alternative approach, HAVOK [2], rooted in Takens embedding theorem [12],
seeks to model the trajectories of (1) by considering a forced linear system,
whose dynamics are precisely the Koopman operator. The states and forcing term
are obtained from the singular value decomposition of a Hankel matrix
$\mathbf{H}_{\mathbf{x}}$, formed by delayed measurements of $\mathbf{x}_{k}$.
As shown in [2] this approach successfully recovers the trajectories of
nonlinear chaotic systems, as linear combinations of a given basis. However,
this linear reconstruction, combined with the difficulty of identifying the
linear dynamics from the svd of $\mathbf{H}_{\mathbf{x}}$ [2] can lead to high
order models (e.g. a 14th order model for the third order Lorentz system).

In this paper, motivated by [4, 13, 2, 7, 10], we propose an alternative,
convex optimization based, approach to the problem of data-driven
identification of Koopman operators. The philosophy, illustrated in Fig. 1,
uses delay coordinates, but, as in [7, 10] does not impose that the state of
the system belongs to span of the Koopman eigenfunctions. Rather, we identify
a manifold of latent variables where the dynamics are linear and map back to
state-space via a non-linear transformations. The problems of finding the
embedding manifold, the associated Koopman operators and the mapping back to
state-space are all recast as rank-constrained semi-definite programs (SDPs).
In turn, these can be relaxed to convex optimizations using the standard
weighted nuclear norm surrogate for rank. Advantages of the proposed approach
include:

  * •

A simple rank check allows for certifying that the solution to these convex
SDPs is indeed the Koopman operator underlying the given data.

  * •

Does not specify a priory the dimension of the embedding or the order of the
dynamics. Rather, both of these can be obtained from the solution to the SDPs.

  * •

Minimizing the order of the linear dynamics leads to simpler models than
competing methods.

  * •

In cases where the spectrum of the Koopman operator is not finite, it allows
for obtaining finite dimensional approximations with guaranteed approximation
error.

  * •

These SDPs have an underlying structure, chordal sparsity, that can be
exploited to substantially reduce computational complexity, leading to
algorithms that scale linearly with the number of data points.

The paper is organized as follows. In section 2 we formally state the problem
under consideration and summarize some needed results on rational
interpolation. Section 3 contains the main results of the paper. It shows that
a Hilbert space $\mathbb{H}$ of observables, its associated Koopman dictionary
and eigenfunctions, and the mapping back to state-space can be found by
solving rank-constrained SDPs. Section 4 illustrates the proposed approach
with some simple examples. Finally, Section 5 summarizes the paper and points
out to directions for extending its results.

##  2 Preliminaries

For ease of reference, next we summarize our notation and recall some results
on interpolation.

###  2.1 Notation

$|\mathcal{S}|$ |  cardinality of the set $\mathcal{S}$
---|---
$\mathbf{x},\mathbf{M}$ |  a vector in $\mathbb{R}^{n}$ (matrix in $\mathbb{R}^{n\times m}$)
$\otimes$ |  Matrix Kronecker product
$\mathbf{M}\succeq 0$ |  the matrix $\mathbf{M}$ is positive semidefinite.
$\|\mathbf{M}\|_{*}$ |  nuclear norm: $\|\mathbf{M}\|_{*}=\Sigma\text{ singular values of $\mathbf{M}$}$.
$\mathbf{H}^{m}_{y}$ |  Hankel matrix with $m$ columns associated with a vector sequence $\mathbf{y}(.)$, with block elements $(\mathbf{H}_{y}^{m})_{i,j}=\mathbf{y}_{i+j-1}$
svec$(\mathbf{M)}$ |  (column-wise) vectorization of the unique elements of a symmetric matrix $\mathbf{M}$.
smat$(\mathbf{v)}$ |  create a symmetric matrix $\mathbf{M}$ from the elements of $\mathbf{v}$ such that svec$(\mathbf{M})=\mathbf{v}$

###  2.2 Rational Interpolants and Loewner Matrices

Given $2n$ scalar pairs $(x_{i},y_{i})$, consider the problem of finding a
rational function
$g(x)\doteq\frac{\sum_{k=1}^{m}a_{k}x^{k}}{\sum_{k=1}^{m}b_{k}x^{k}}$ such
that $y_{i}=g(x_{i}),\;i=1,\ldots 2n$. Define the Loewner matrix

$\mathbf{L}\in\mathbb{R}^{n\times n}$ as :

$\mathbf{L}(x,y)=\begin{bmatrix}\frac{y_{1}-y_{n+1}}{x_{1}-x_{n+1}}&\frac{y_{1}-y_{n+2}}{x_{1}-x_{n+2}}&\ldots&\frac{y_{1}-y_{2n}}{x_{1}-x_{2n}}\\\
\frac{y_{2}-y_{n+1}}{x_{2}-x_{n+1}}&\frac{y_{2}-y_{n+2}}{x_{2}-x_{n+2}}&\ldots&\frac{y_{2}-y_{2n}}{x_{2}-x_{2n}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{y_{n}-y_{n+1}}{x_{n}-x_{n+1}}&\frac{y_{n}-y_{n+2}}{x_{n}-x_{n+2}}&\ldots&\frac{y_{n}-y_{2n}}{x_{n}-x_{2n}}\end{bmatrix}$
(4)

Then, there exists a rational function of order at most $m$ that interpolates
the given data points if and only if rank$(\mathbf{L})\leq m-1$ [1, 6].

###  2.3 Statement of the problem

Consider the nonlinear dynamical system:

$\begin{split}\mathbf{x}_{k+1}&=f(\mathbf{x}_{k},\ldots,\mathbf{x}_{k-r+1})\;\qquad\mathbf{x}_{j}\in\mathbb{R}^{n}\end{split}$
(5)

where both the dynamics $f(.)$ and its order $r$ are unknown. Our goal is to
identify its associated Koopman operator, over a suitable space of
observables, from experimental data $\mathbf{x}$. Specifically:

######  Problem 1.

Given a set of $N$ trajectories
$\\{\mathbf{x}_{k}^{(i)}\\}_{k=1}^{T_{i}},\;i=1,\ldots,N$,
$\mathbf{x}_{k}^{(i)}\in\mathbb{R}^{n}$, find a (functional) dictionary
$\mathcal{D}(.)$, a Hilbert space $\mathbb{H}$ of observables
$\boldsymbol{\psi}(.)$ of the form:

$\begin{split}&\boldsymbol{\psi}(\boldsymbol{\xi}_{k})\doteq\begin{bmatrix}\mathbf{y}_{k-r+1}^{T}\ldots\mathbf{y}_{k}^{T}\end{bmatrix}^{T}\in\text{span
$\\{\mathcal{D}(\boldsymbol{\xi}_{k})\\}$ with
$\mathbf{y}_{j}\in\mathbb{R}^{m}$}\\\ &\text{where
$\boldsymbol{\xi}_{k}\doteq\begin{bmatrix}\mathbf{x}_{k-r+1}\ldots\mathbf{x}_{k}\end{bmatrix}^{T}$}\end{split}$
(6)

and an operator $\mathcal{K}\colon\mathbb{H}\to\mathbb{H}$ such that
$(\mathcal{K}\circ\boldsymbol{\psi})(\boldsymbol{\xi}_{k})=\boldsymbol{\psi}(\boldsymbol{\xi}_{k+1})$.

Problem 1 is reminiscent of EDMD approaches. However, the main difference is
that here we seek to learn the dictionary $\mathcal{D}$ and the dimensions of
the space $\mathbb{H}$ directly from the data, rather than postulating a fixed
dictionary and dimension. Further, if Problem 1 has a solution, the resulting
operator $\mathcal{K}$ is indeed the exact Koopman operator in $\mathbb{H}$.

######  Remark 1.

As stated, Problem 1 is ill posed, since $\|\boldsymbol{\psi}(.)\|$ can be
arbitrarily small or large. To avoid this, and with an eye towards
reconstruction of $\boldsymbol{\xi}$ from $\boldsymbol{\psi}$, we will impose
the additional constraints:

$\begin{split}&\frac{1}{M_{\ell}(\boldsymbol{\xi}_{i})}\|\boldsymbol{\psi}(\boldsymbol{\xi}_{i})-\boldsymbol{\psi}(\boldsymbol{\xi}_{j})\|_{2}\leq\|\boldsymbol{\xi}_{i}-\boldsymbol{\xi}_{j}\|_{2}\leq
M_{u}(\boldsymbol{\xi}_{i})\|\boldsymbol{\psi}(\boldsymbol{\xi}_{i})-\boldsymbol{\psi}(\boldsymbol{\xi}_{j})\|_{2}\;\text{$\forall\boldsymbol{\xi}_{j}$
such that $\|\boldsymbol{\xi}_{i}-\boldsymbol{\xi}_{j}\|_{2}\leq\delta$}\\\
&\mathbf{y}(\boldsymbol{\xi}_{k})_{j}^{T}\mathbf{y}(\boldsymbol{\xi}_{k})_{j}=(\boldsymbol{\xi}_{k})_{j}^{T}(\boldsymbol{\xi}_{k})_{j}\;j=1,\ldots
r,\;\text{for all $k$ in a given set of ``anchor" points
$\mathbb{I}$}\end{split}$

where $\mathbf{y}_{j}$ denotes the $j^{\rm{th}}$ block component of
$\boldsymbol{\psi}(\boldsymbol{\xi}_{k})$ and the scalar $\delta$ and the set
of anchor points $\mathbb{I}$ are design hyperparameters. That is, we impose
that (a) the mapping $\phi\colon\boldsymbol{\xi}\to\boldsymbol{\psi}$ and its
inverse are locally Lipschitz continuous, with Lipschitz constants
$M_{\ell}(\boldsymbol{\xi}_{i})$ and $M_{u}(\boldsymbol{\xi}_{i})$; and (b)
the function $\boldsymbol{\psi}(.)$ is normalized to have components with
unity gain at some given “anchor” points.

##  3 Learning Koopman Operators via Semi Definite Optimization

In this section we present the main theoretical result of the paper: a
reformulation of Problem 1 as a rank minimization subject to a positive semi-
definite constraint. Since this problem is generically NP hard, we then
develop a tractable convex relaxation, along with optimality certificates.

###  3.1 Finding Koopman operators as a constrained rank minimization

Consider the following feasibility problem (in $\mathbf{y},r,m$):

######  Problem 2.

Given a set of $N$ trajectories
$\\{\mathbf{x}_{\ell}^{(i)}\\}_{\ell=1}^{T_{i}},\;i=1,\ldots N$,
$\mathbf{x}_{\ell}^{(i)}\in\mathbb{R}^{n}$, find scalars $r,m$ and $N$
trajectories $\mathbf{y}_{k}^{(i)}\in\mathbb{R}^{m},k=1,\ldots,T_{i}$, such
that the following holds:

$\displaystyle\text{rank $(\mathbf{H_{y}}^{(r+1)})\leq r$, where
$\mathbf{H_{y}}^{(r+1)}\doteq\begin{bmatrix}\mathbf{H}^{(r+1)}_{\mathbf{y}^{(1)}}\\\
\vdots\\\
\mathbf{H}^{(r+1)}_{\mathbf{y}^{(N)}}\end{bmatrix}$}\;\text{and}\;\mathbf{H}_{\mathbf{y}^{(i)}}^{(r+1)}\doteq\begin{bmatrix}\mathbf{y}_{1}^{(i)}&\mathbf{y}_{2}^{(i)}&\cdots&\mathbf{y}_{r+1}^{(i)}\\\
\mathbf{y}_{2}^{(i)}&\mathbf{y}_{3}^{(i)}&\cdots&\mathbf{y}_{r+2}^{(i)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{y}_{T_{i}-r}^{(i)}&\mathbf{y}_{T_{i}-r+1}^{(i)}&\cdots&\mathbf{y}_{T_{i}}^{(i)}\end{bmatrix}$
(7) $\displaystyle\begin{array}[]{l}\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}\leq
M_{u}(\mathbf{x}_{s})\|\mathbf{y}_{s}-\mathbf{y}_{t}\|_{2}\\\
\|\mathbf{y}_{s}-\mathbf{y}_{t}\|_{2}\leq
M_{\ell}(\mathbf{x}_{s})\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}\end{array}\Big{\\}}\text{$\forall(s,t)$
such that $\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}\leq\delta$}$ (10)
$\displaystyle\;\|\mathbf{y}_{s}\|_{2}=\|\mathbf{x}_{s}\|_{2}\;\text{for all
$s\in\mathbb{I}$}$ (11)

As shown next, the solution to Problem 1 (e.g the dictionary $\mathcal{D}$,
the embedding Hilbert space $\mathbb{H}$ and the associated Koopman operator)
can be constructed from any feasible solution to (7)-(11).

######  Theorem 1.

Let $(\mathbf{y}^{(i)}_{k},r,m)$ denote a feasible solution to (7)-(11) with
$\mathbf{y}_{k}^{(i)}\in\mathbb{R}^{m}$ and rank$(\mathbf{H_{y}})=r*\leq r$.
Let $\mathbf{H_{y}}^{(r^{*}+1)}$ and
$\mathcal{N}_{R}(\mathbf{H_{y}}^{(r^{*}+1)})$ denote the Hankel matrix
obtained by rearranging the elements of $\mathbf{H_{y}}^{(r+1)}$ into
$r^{*}+1$ columns, and its right null space, respectively. Note that by
construction rank$(\mathbf{H_{y}}^{(r^{*}+1)})=r^{*}$ and thus
$dim(\mathcal{N}_{R}(\mathbf{H_{y}}^{(r^{*}+1)}))\geq 1$. Consider a vector
$\mathbf{p}\in\mathcal{N}_{R}(\mathbf{H_{y}}^{(r^{*}+1}))$, of the form
$\mathbf{p}=\begin{bmatrix}a_{0}&\ldots&a_{(r^{*}-1)}&-1\end{bmatrix}^{T}$.
Let $\rho_{j},\;j=1,\ldots,r^{*}$ denote the roots of the polynomial
$\mathcal{P}(\rho)\doteq\rho^{r^{*}}-\sum_{i=0}^{r^{*}-1}a_{i}\rho^{i}$ and
define the $r^{*}$ vectors

$\mathbf{v}_{j}\doteq\begin{bmatrix}1&\rho_{j}&\rho_{j}^{2}&\ldots\rho_{j}^{r^{*}}\end{bmatrix}^{T}$

Finally, let $\mathbf{V}$ denote the Vandermonde matrix
$\mathbf{V}=\begin{bmatrix}\mathbf{v}_{1}&\mathbf{v}_{2}&\ldots&\mathbf{v}_{r^{*}}\end{bmatrix}$.
Then:

  1. 1.

The desired dictionary $\mathcal{D}(.)$ has the matrix representation
$\mathbf{D}=\mathbf{V}\otimes\mathbf{I}_{m}$.

  2. 2.

The Hilbert space $\mathbb{H}$ of observables is given by span($\mathcal{D}$),
with the usual inner product.

  3. 3.

The operator $\mathcal{K}\colon\mathbb{H}\to\mathbb{H}$ with the matrix
representation
$\boldsymbol{\Lambda}=\text{diag}(\rho_{i})\otimes\mathbf{I}_{m}$ in the basis
defined by the columns of $\mathbf{V}$ is the Koopman operator associated with
(1) in the space $\mathbb{H}$.

###### Proof.

Given in the Appendix. ∎

Theorem 1 provides the foundation for constructing the Koopman operator from
the solution of an optimization problem, but is of limited practical value,
due to several reasons: (i) It does not indicate how to find $m$, the
dimension of $\mathbf{y}_{k}$, or $r$, the “memory” of the system, and (ii) it
leads to a difficult, non-convex problem. Motivated by [4], next we show that
Problem 2 is equivalent to a SDP constrained rank-minimization. The starting
point is to consider the Kernel matrix with entries
$\mathbf{K}_{r,s}\doteq\mathbf{y}_{r}^{T}\mathbf{y}_{s}$, where
$\mathbf{y}_{r},\mathbf{y}_{s}$ denote the observables corresponding to points
$\mathbf{x}_{r},\mathbf{x}_{s}$ drawn from (not necessarily the same) training
trajectories. Let $\mathbf{y}_{s}^{(i)},\;s=1,\ldots T_{i}$ denote the
observables corresponding to the ith trajectory and define the
$(r+1)\times(r+1)$ Gram matrix
$\mathbf{G}^{(i)}\doteq(\mathbf{H}^{(r+1)}_{\mathbf{y}^{(i)}})^{T}\mathbf{H}^{(r+1)}_{\mathbf{y}^{(i)}}$.
The key observation is that both the entries of $\mathbf{G}^{(i)}$ and the
argument of the constraints (10)–(11) are _affine_ functions of entries of
$\mathbf{K}$, leading to the following result:

######  Theorem 2.

Define the family of Gram matrices:
$\mathbf{G}^{(i)}=(\mathbf{H}^{(r+1)}_{\mathbf{y}^{(i)}})^{T}\mathbf{H}^{(r+1)}_{\mathbf{y}^{(i)}}=\sum_{\ell=0}^{T_{i}-r}\mathbf{K}_{\ell,r}^{(i)}$
where

$\mathbf{K}_{\ell,r}^{(i)}=\begin{bmatrix}(\mathbf{y}_{\ell}^{(i)})^{T}\mathbf{y}_{\ell}^{(i)}&(\mathbf{y}_{\ell}^{(i)})^{T}\mathbf{y}_{\ell+1}^{(i)}&\cdots&(\mathbf{y}_{\ell}^{(i)})^{T}\mathbf{y}_{\ell+r}^{(i)}\\\
\vdots&\vdots&\ddots&\vdots\\\
(\mathbf{y}_{\ell+r}^{(i)})^{T}\mathbf{y}_{\ell}^{(i)}&(\mathbf{y}_{\ell+r}^{(i)})^{T}\mathbf{y}_{\ell+1}^{(i)}&\cdots&(\mathbf{y}_{\ell+r}^{(i)})^{T}\mathbf{y}_{\ell+r}^{(i)}\\\
\end{bmatrix}$

(note that $\mathbf{K}_{\ell,r}^{(i)}$ are submatrices of $\mathbf{K}$).
Consider the following rank minimization problem:

$\displaystyle r^{*}=$ $\displaystyle\min_{\mathbf{K}\succeq
0}\text{rank}(\mathbf{G}\doteq\begin{bmatrix}(\mathbf{G}^{(1)})^{T}&\ldots&(\mathbf{G}^{(N)})^{T}\end{bmatrix}^{T})\;\text{subject
to:}$ (12)
$\displaystyle\left.\begin{array}[]{l}\frac{1}{M^{2}_{u}(\mathbf{x}_{s})}\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}^{2}\leq
K_{s,s}-2K_{s,t}+K_{t,t}\\\ K_{s,s}-2K_{s,t}+K_{t,t}\leq
M^{2}_{\ell}(\mathbf{x}_{s})\|\mathbf{x}_{s}-\mathbf{x}_{t}\|^{2}_{2}\end{array}\right\\}\text{$\forall(s,t)$
such that $\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}\leq\delta$}$ (15)
$\displaystyle K_{s,s}=\|\mathbf{x}_{s}\|_{2}^{2}\;\text{for all
$s\in\mathbb{I}$}$ (16)

Denote by $\mathbf{K}^{(i)}$ the submatrix of $\mathbf{K}$ with entries
$(\mathbf{K}^{(i)})_{\ell,j}=(\mathbf{y}^{(i)}_{\ell})^{T}\mathbf{y}^{(i)}_{j}$,
and let $m=\max_{i}\left\\{\text{rank}(\mathbf{K}^{(i)})\right\\}$. Consider
the factorizations $(\mathbf{Y}^{(i)})^{T}\mathbf{Y}^{(i)}=\mathbf{K}^{(i)}$
with $\mathbf{Y}^{(i)}\in\mathbb{R}^{m\times T_{i}}$. Then, if $r*<r+1$, the
columns $\mathbf{y}_{k}^{(i)}$ of $\mathbf{Y}^{(i)}$ solve Problem 2.

###### Proof.

Given in the Appendix ∎

###  3.2 Adding a regularization

Theorems 1 indicates how to find the observables
$\boldsymbol{\psi}(.)\in\mathbb{H}$ by solving a constrained optimization
problem. Further, these constraints guarantee that the mapping
$\boldsymbol{\psi}(.)\colon\mathbb{R}^{rn}\to\mathbb{H}$ locally satisfies
some Lipschiz and gain constraints. However these constraints alone do not
guarantee that $\boldsymbol{\psi}(.)$ is not arbitrarily complex, or even has
the same functional form for all $\boldsymbol{\xi}$. These issues can
complicate the task of finding an explicit form for the mapping, if one is
needed. Next, we briefly indicate how to use additional degrees of freedom
available in the problem to guarantee that $\boldsymbol{\psi}(.)$ is the
simplest possible mapping, in a sense precisely defined below, and has the
same functional form for all $\boldsymbol{\xi}$.

Consider a point
$\boldsymbol{\xi}_{k}\doteq\begin{bmatrix}\mathbf{x}_{k-r+1}^{T}&\ldots&\mathbf{x}_{k}^{T}\end{bmatrix}^{T}$,
and for each (block) component $\mathbf{x}_{j}$, denote by
$\mathbb{N}_{\mathbf{x}_{j}}$ the indexes of its nearest neighbors. Let
$\mathbf{K}_{\mathbf{x}_{j}}$ be matrix with elements
$(\mathbf{K}_{\mathbf{x}_{j}})_{r,s}=\mathbf{x}^{T}_{r}\mathbf{x}_{s}$ for all
$r,s\in\left\\{j\cup\mathbb{N}_{\mathbf{x}_{j}}\right\\}$. Similarly, given
$\boldsymbol{\psi}(\boldsymbol{\xi}_{k})\doteq\begin{bmatrix}\mathbf{y}_{k-r+1}\ldots\mathbf{y}_{k}\end{bmatrix}^{T}$,
let $\mathbf{K}_{\mathbf{y}_{j}}$ be the submatrix of $\mathbf{K}$ with
elements
$(\mathbf{K}_{\mathbf{y}_{j}})_{r,s}=\mathbf{y}^{T}_{r}\mathbf{y}_{s}$ for all
$r,s\in\left\\{j\cup\mathbb{N}_{\mathbf{x}_{j}}\right\\}$. For ease of
notation, let
$\boldsymbol{\kappa}_{\mathbf{x}}^{(j)}\doteq\textbf{svec$(\mathbf{K}_{\mathbf{x}_{j}})$}\in\mathbb{R}^{q}$,
$\boldsymbol{\kappa}_{\mathbf{y}}^{(j)}\doteq\textbf{svec$(\mathbf{K}_{\mathbf{y}_{j}})$}\in\mathbb{R}^{q}$,
where
$q\doteq\frac{(|\mathbb{N}_{\mathbf{x}_{j}}|+1)(|\mathbb{N}_{\mathbf{x}_{j}}|+2)}{2}$.
Note that these vectors contain the unique elements of the matrices
$\mathbf{K}_{\mathbf{x}_{j}}$, $\mathbf{K}_{\mathbf{x}_{j}}$. Finally, let
$p=\lfloor\frac{q}{2}\rfloor$ and define the Loewner matrix

$\mathbf{L}_{\mathbf{x}_{j}}\doteq\begin{bmatrix}\frac{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{y}_{p+1}}}{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+1}}}&\frac{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+2}}}{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+2}}}&\ldots&\frac{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{y}_{q}}}{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{q}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{\boldsymbol{\kappa}_{\mathbf{y}_{p}}-\boldsymbol{\kappa}_{\mathbf{y}_{p+1}}}{\boldsymbol{\kappa}_{\mathbf{x}_{p}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+1}}}&&\ldots&\frac{\boldsymbol{\kappa}_{\mathbf{y}_{p}}-\boldsymbol{\kappa}_{\mathbf{y}_{q}}}{\boldsymbol{\kappa}_{\mathbf{x}_{p}}-\boldsymbol{\kappa}_{\mathbf{x}_{q}}}\\\
\end{bmatrix}$ (17)

where
$\boldsymbol{\kappa}_{\mathbf{x}_{i}},\boldsymbol{\kappa}_{\mathbf{y}_{i}}$
denote the $i^{\rm{th}}$ component of $\boldsymbol{\kappa}_{\mathbf{x}}^{(j)}$
and $\boldsymbol{\kappa}_{\mathbf{y}}^{(j)}$ respectively. From the results in
section 2.2, it follow that if rank$(\mathbf{L}_{\mathbf{x}_{j}})<p$, then
there exists a rational mapping of degree up to $p-1$ that maps the elements
of $\mathbf{K}_{\mathbf{x}_{j}}$ to those of $\mathbf{K}_{\mathbf{y}_{j}}$.
Further, the degree of this mapping can be minimized by minimizing the rank of
$\mathbf{L}_{\mathbf{x}_{j}}$ with respect to the variables
$\boldsymbol{\kappa}_{\mathbf{y}_{i}}$, leading (locally) to the lowest order
rational mapping $\boldsymbol{\psi}(\boldsymbol{\xi}_{k})$. If a global,
rather than local, rational mapping is desired, a similar idea can be using
involving all pairs $\mathbf{x},\mathbf{y}$, rather than just the nearest
neighbors of each point.

###  3.3 A Convex Relaxation

Theorem 2 allows for reducing Problem 1 to a constrained rank minimization
problem. However, this problem is still NP-hard. In order to obtain a
tractable relaxation, we will replace the objective (12) by
$\sum_{i=1}^{N}\text{rank}(\mathbf{G}^{(i)})$ and add a term of the form
$\lambda_{1}\sum_{j=1}^{T}\text{rank}(\mathbf{L}_{\mathbf{x}_{j}})$, where
$T=\sum T_{i}$ is the total number of points. Then, proceeding as in [9], we
will replace rank with a convex surrogate, a weighted nuclear norm, where the
weights are updates as each step of the algorithm. Finally, in order to handle
outliers, we will consider a “soft” version of (15)-(16), where these are
added to the objective as penalties. The complete algorithm is outlined in
Algorithm 1. It is worth noting that if the algorithm yields a solution
$\mathbf{G}$ with rank$(\mathbf{G})<r$, this certifies that $m$ is indeed the
Koopman operator. On the other hand, if the algorithm yields a solution
$\mathbf{G}$ with minimum singular value $\sigma_{\text{min}}$, then an
$r^{th}$ order approximate model $m_{r}$ can be obtained by performing PCA on
$\mathbf{G}$. In this case the approximation error is bounded (in the Hankel
norm sense) by $\sqrt{\sigma_{\text{min}}}$.

Algorithm 1 Reweighted $\|.\|_{*}$ based Koopman Identification

1:initialize:
$iter=0,\mathbf{W}_{0}=\mathbf{I};\mathbf{V}_{0}^{(j)}=\mathbf{I},j=1,\ldots,T$;
$\lambda_{1},\lambda_{2},\lambda_{3},\delta\leftarrow$ hyperparameters,
---
$\mathbb{NN}=\left\\{(s,t)\colon\|\mathbf{x}_{s}-\mathbf{x}_{t}\|_{2}\leq\delta\right\\}$,
$\sigma\leftarrow$small number, $r\leftarrow$ upper bound on system order.

2:Repeat: Solve  $\begin{array}[]{l}\min_{\mathbf{K}^{(i)}\succeq
0}\quad\|\mathbf{W}_{iter}\mathbf{G}\|_{*}+\lambda_{1}\sum_{j=1}^{T}\|\mathbf{V}^{(j)}_{iter}\mathbf{L}_{\mathbf{x}_{j}}\|_{*}+\lambda_{2}\sum_{s\in\mathbb{I}}(K_{s,s}-\|\mathbf{x}_{s}\|_{2}^{2})^{2}\\\
+\lambda_{3}\sum_{{r,t}\in\mathbb{NN}}\max\left\\{0,\frac{1}{M^{2}_{u}(\mathbf{x}_{s})}\|\mathbf{x}_{s}-\mathbf{x}_{t}\|^{2}_{2}-K_{s,s}+2K_{s,t}-K_{t,t}\right\\}\\\
+\lambda_{3}\sum_{{r,t}\in\mathbb{NN}}\max\left\\{0,K_{s,s}-2K_{s,t}+K_{t,t}-M^{2}_{\ell}(\mathbf{x}_{s})\|\mathbf{x}_{s}-\mathbf{x}_{t}\|^{2}_{2}\right\\}\end{array}$
Update
$\begin{array}[]{l}\mathbf{W}_{(iter+1)}=\left(\frac{\mathbf{G}+\sigma\mathbf{I}}{\|\mathbf{G}+\sigma\mathbf{I}\|}\right)^{-1},\;\mathbf{V}^{(j)}_{(iter+1)}=\left(\frac{\mathbf{L}_{\mathbf{x}_{j}}+\sigma\mathbf{I}}{\|\mathbf{L}_{\mathbf{x}_{j}}+\sigma\mathbf{I}\|}\right)^{-1},\;iter=iter+1\end{array}$

3:Until: rank$(\mathbf{G})<r$.

4:$[\mathbf{U}^{(i)},\mathbf{S}^{(i)},(\mathbf{U}^{(i)})^{T}]\leftarrow\text{svd}(\mathbf{K}^{(i)})$,
$\mathbf{S}^{(i)}\leftarrow\frac{\mathbf{S}^{(i)}}{\|\mathbf{S}^{(i)}\|}$,
$r_{k}^{(i)}\leftarrow\min r\colon\sum_{j=1}^{r}\mathbf{S}_{jj}^{(i)}\geq
0.99$

5:$\mathbf{Y}^{(i)}\leftarrow[\mathbf{U}^{(i)}(:,1:r_{k})]^{T}$

6:$[\mathbf{U}_{\mathbf{G}},\mathbf{R},\mathbf{V}_{\mathbf{G}}^{T}]\leftarrow\text{svd}(\mathbf{G})$,
$\mathbf{m}\leftarrow\mathbf{V}_{\mathbf{G}}(:,r+1)$

7:Output: embeddings $\mathbf{Y^{(i)}}$, model $\mathbf{m}$.

###  3.4 Mapping observables to states

The approach presented in Section 3 finds the observables
$\boldsymbol{\psi}(\boldsymbol{\xi}_{k})$ corresponding to a given trajectory
$\boldsymbol{\xi}_{k},\;k=1,\ldots T$. However, it does not explicitly provide
a method for mapping a given $\boldsymbol{\psi}(\boldsymbol{\xi})$, obtained
for instance by using the Koopman operator to propagate a trajectory in
observable space, back to the corresponding point $\boldsymbol{\xi}$ in state
space. Motivated by [11] we propose to find (pointwise) the mapping
$\boldsymbol{\psi}\to\boldsymbol{\xi}$ by locally approximating the mapping
between the embedded space and ambient space kernels,
$\mathbf{K}_{\mathbf{y}}$ and $\mathbf{K}_{\mathbf{x}}$, with a rational
function. Specifically, given a point $\mathbf{y}^{*}\in\mathbb{R}^{m}$, let
$\mathcal{N}_{\mathbf{y^{*}}}\doteq\left\\{\mathbf{y}_{k}\colon\|\mathbf{y}^{*}-\mathbf{y}_{k}\|_{2}^{2}\leq\delta\right\\}$
and denote by $\mathcal{X}$ its preimage. We propose to estimate
$\mathbf{x}^{*}$ by first finding ${K}_{\mathbf{x}_{i},\mathbf{x}^{*}}$, the
elements of $\mathbf{K}_{x}$ corresponding to
$\mathbf{x}_{i}^{T}\mathbf{x}^{*},\;\forall\mathbf{x}_{i}\in\mathcal{X}$ and
then finding $\mathbf{x}^{*}$ by factorizing $\mathbf{K}_{x}$. Note that, in
order to get a valid kernel compatible with the priors, the elements
${K}_{\mathbf{x}_{i},\mathbf{x}^{*}}$ should be such that the completed matrix
$\mathbf{K}_{\mathbf{x}}\succeq 0$, rank$(\mathbf{K})\leq n$, and the
constraints (10)-(11) are satisfied. As shown next, under the assumption that
the mapping $G\colon\mathbf{K}_{\mathbf{y}}\to\mathbf{K}_{\mathbf{x}}$ is
rational, then $\mathbf{x}^{*}$ can be found by solving a rank minimization
problem subject to semi-definite constraints.

Consider the Kernel matrices
$\mathbf{K}_{\mathbf{x}},\mathbf{K}_{\mathbf{y}}\in\mathbb{R}^{(|\mathcal{X}|+1)\times(|\mathcal{X}|+1)}$,
where the entries have been ordered so that the elements of the form
$\mathbf{y}_{i}^{T}\mathbf{y}^{*}$ and $\mathbf{x}_{i}^{T}\mathbf{x}^{*}$
appear in the first row and column. As before, for ease of notation, let
$\boldsymbol{\kappa}_{\mathbf{x}}=\textbf{svec$(\mathbf{K}_{\mathbf{x}})$}$,
$\boldsymbol{\kappa}_{\mathbf{y}}=\textbf{svec$(\mathbf{K}_{\mathbf{y}})$}$.
Note that
$\boldsymbol{\kappa}_{\mathbf{x}},\boldsymbol{\kappa}_{\mathbf{y}}\in\mathbb{R}^{q}$,
with $q\doteq\frac{(|\mathcal{X}|+1)(|\mathcal{X}|+2)}{2}$, and that all inner
products involving $\mathbf{y}^{*}$ and $\mathbf{x}^{*}$ appear in the first
$|\mathcal{X}|+1$ elements of $\boldsymbol{\kappa}_{\mathbf{y}}$ and
$\boldsymbol{\kappa}_{\mathbf{x}}$. Let $p=\lfloor\frac{q}{2\rfloor}$ and
consider the following rank minimization problem:

$\displaystyle\min_{\boldsymbol{\kappa}_{\mathbf{x}}}\text{rank
$(\mathbf{L})$}\text{subject to:}$ (18)
$\displaystyle\left.\begin{tabular}[]{l}$\kappa_{\mathbf{x}_{1}}-2\boldsymbol{\kappa}_{\mathbf{x}_{i}}+\boldsymbol{\kappa}_{\mathbf{x}_{j}}\leq\delta^{2}$\\\
$\boldsymbol{\kappa_{\mathbf{x}_{1}}}-2\boldsymbol{\kappa}_{\mathbf{x}_{i}}+\boldsymbol{\kappa}_{\mathbf{x}_{j}}\leq
M_{u}^{2}(\mathbf{x}_{i})(\kappa_{\mathbf{y}_{1}}-2\boldsymbol{\kappa}_{\mathbf{y}_{i}}+\boldsymbol{\kappa}_{\mathbf{y}_{j}})$\\\
$\boldsymbol{\kappa_{\mathbf{y}_{1}}}-2\boldsymbol{\kappa}_{\mathbf{y}_{i}}+\boldsymbol{\kappa}_{\mathbf{y}_{j}}\leq
M_{\ell}^{2}(\mathbf{x}_{i})(\kappa_{\mathbf{x}_{1}}-2\boldsymbol{\kappa}_{\mathbf{x}_{i}}+\boldsymbol{\kappa}_{\mathbf{x}_{j}})$\\\
\end{tabular}\right\\}\begin{tabular}[]{l}$i=2,\ldots,|\mathcal{X}|+1$\\\
$j=\frac{(2|\mathcal{X}|+4-i)(i-1)}{2}+1$\end{tabular}$ (24)
$\displaystyle\mathbf{K}_{\mathbf{x}}\doteq\textbf{smat}(\boldsymbol{\kappa}_{\mathbf{x}})\succeq
0,\;\text{rank}(\mathbf{K}_{\mathbf{x}})\leq
n,\boldsymbol{\kappa}_{\mathbf{x}_{1}}=\boldsymbol{\kappa}_{\mathbf{y}_{1}}$
(25)
$\displaystyle\mathbf{L}=\begin{bmatrix}\frac{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+1}}}{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{y}_{p+1}}}&\frac{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+2}}}{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{y}_{p+2}}}&\ldots&\frac{\boldsymbol{\kappa}_{\mathbf{x}_{1}}-\boldsymbol{\kappa}_{\mathbf{x}_{q}}}{\boldsymbol{\kappa}_{\mathbf{y}_{1}}-\boldsymbol{\kappa}_{\mathbf{y}_{q}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{\boldsymbol{\kappa}_{\mathbf{x}_{p}}-\boldsymbol{\kappa}_{\mathbf{x}_{p+1}}}{\boldsymbol{\kappa}_{\mathbf{y}_{p}}-\boldsymbol{\kappa}_{\mathbf{y}_{p+1}}}&&\ldots&\frac{\boldsymbol{\kappa}_{\mathbf{x}_{p}}-\boldsymbol{\kappa}_{\mathbf{x}_{q}}}{\boldsymbol{\kappa}_{\mathbf{y}_{p}}-\boldsymbol{\kappa}_{\mathbf{y}_{q}}}\\\
\end{bmatrix}$

######  Theorem 3.

Let $\boldsymbol{\kappa}^{*}_{\mathbf{x}},\mathbf{L}^{*}$ denote the solution
to (18)-(25). If rank$(\mathbf{L}^{*})<p$, then (i) there exist a rational
function $g(.)$ of degree at most $p$ such that
$g(\kappa_{\mathbf{y}_{i}})=\kappa_{\mathbf{x}_{i}}$; and (ii) the vector
$\mathbf{x}^{*}$ defined by the first row of $\mathbf{X}$, where
$\mathbf{X}^{T}\mathbf{X}=\mathbf{K}_{\mathbf{x}}$ satisfies constraints
(10)-(11) in Problem 2.

###### Proof.

Given in the Appendix ∎

Relaxing the rank in (18) and (25) to a weighed nuclear norm, leads to an
algorithm similar to Algorithm 1, based on solving a sequence of SDPs until
rank deficient matrices $\mathbf{L},\mathbf{K_{x}}$ are obtained.

##  4 Illustrative Examples

Example 1: Lorentz Attractor. In this example we consider the Lorentz chaotic
system:

$\displaystyle\dot{x}_{1}=\sigma(x_{2}-x_{1});\dot{x}_{2}=x_{1}(\rho-
x_{3})-x_{2};\dot{x}_{3}=x_{1}x_{2}-\beta x_{3};$ (26)

with parameters $\sigma=28,\rho=10,\beta=\frac{8}{3}$. We used 400 points of
the trajectory starting at $[-10.38\;-4.5366\;35.1640]^{T}$, uniformly sampled
every $0.0271$ seconds to find the embeddings, and matlab’s command ssest to
estimate an 7th order model. Fig 2(a) shows the training and one step ahead
reconstructed data, that is the results of applying the encoder/decoder
illustrated on the top of Fig. 1 to (i) train, (ii) project the training data,
(iii) perform a one step ahead prediction and (iv) lift back. Figure 2(b)
shows the predictions obtained using the pipeline at the bottom of Fig. 1, for
points not part of the training data. As shown there, the proposed pipeline is
indeed able to predict with reasonable accuracy the one step ahead value of
the trajectory, using a 7th order Koopman operator. For comparison, [2] uses a
14th order model.

Figure 2: Lorentz attractor: Left: one step ahead prediction of training data.
Right: one step ahead prediction of new data

Example 2: The Duffing Oscillator. Here we consider the system111The
conventional Duffing equation is a forced oscillator. Here we use the last two
equations to generate the forcing term $sin(t)$.:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2};\dot{x}_{2}$
$\displaystyle=-0.5x_{2}-x_{1}-x_{1}^{3}+0.42x_{3};\;\dot{x}_{3}$
$\displaystyle=x_{4};\;\dot{x}_{4}$ $\displaystyle=-x_{3}$ (27)

In this case, Algorithm 1 yielded an embedding $\mathbf{y}\in\mathbb{R}^{3}$.
We then used matlab’s command ssest to estimate a second order model for each
component of $\mathbf{y}$. Fig 3 (left) shows the one step ahead prediction of
the training data. The right panel in Fig. 3 shows the predictions obtained
using the pipeline at the bottom of Fig. 1, for points not part of the
training data. As before, the proposed pipeline successfully predicts the next
point in the trajectory.

Figure 3: Duffing oscillator one step ahead predictions of (left) training
data and (right) new data.

Example 3: predator-prey model. In this example we considered the predator-
prey model:

$\displaystyle\dot{x}_{1}$ $\displaystyle=-x_{1}+x_{1}x_{2}$ (28)
$\displaystyle\dot{x}_{2}$ $\displaystyle=x_{2}-x_{1}x_{2}$

We used 120 points from the trajectory in ambient space to find the
embeddings, and matlab’s command ssest to estimate an 8th order model. Fig
4(a) shows the training and reconstructed data, that is the results of
applying back to back the encoder/decoder illustrated on the top of Fig. 1.
Figure 4(b) shows the predictions obtained using the pipeline at the bottom of
Fig. 1, starting from an initial condition not part of the training data. As
shown there, the proposed pipeline is indeed able to predict with reasonable
accuracy the trajectory over an 80 steps horizon that encompasses all regions
visited by the trajectory.

|
---|---
Figure 4: Predator-prey example. Left: encoding (black) and decoding (red) the
training data. Right: one step ahead state prediction over an 80 step horizon,
blue: ground truth, red: predictions using the proposed pipeline, black:
initial conditions. Figure 5: Two different trajectories corresponding to
close initial conditions.

Example 4: Another chaotic system. Here we consider the system:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2}$ (29)
$\displaystyle\dot{x}_{2}$ $\displaystyle=-x_{1}^{5}-0.1x_{2}+x_{3}$
$\displaystyle\dot{x}_{3}$ $\displaystyle=x_{4}$ $\displaystyle\dot{x}_{4}$
$\displaystyle=-x_{3}$

It is well know that this system has extreme sensitivity to initial condition.
This effect is illustrated in Fig 5 showing two different trajectories
corresponding to the initial conditions $[2.0,3.0,0.0,6.0]$ and
$[2.1,3.0,0.0,6.0]$. We used 100 points from the trajectory in ambient space
to find the embeddings. In this case Algorithm 1 yielded an embedding
$\mathbf{y}\in\mathbb{R}^{3}$. We then used matlab’s command ssest to estimate
a second order model for each component of $\mathbf{y}$. Fig 6(left) shows the
training and reconstructed data. The right panel in Figure 6 shows the
predictions obtained using the pipeline at the bottom of Fig. 1, starting from
an initial condition not part of the training data. As before, the proposed
pipeline successfully predicts with reasonable accuracy the trajectory over a
100 steps horizon, in spite of the sensitivity of the system to initial
conditions noted above.

Figure 6: Example 4: Left: training data. Right: one step ahead predictions
over an 100 step horizon, red: ground truth, blue: predictions using the
proposed pipeline.

##  5 Conclusions

This paper proposes a convex optimization approach to learning Koopman
operators from data. The main idea is to use delay coordinates and nonlinear,
kernel based embeddings to recast the problem as a rank-constrained
optimization. In turn, this optimization can be relaxed to a tractable semi-
definite program. Salient features of this approach are its ability to certify
that the solution to this SDP indeed solves the original problem, and the fact
that neither the order of the embedding nor of the dynamics governing their
evolution need to be specified a-priori. Further, by seeking embeddings that
minimize the order of these dynamics, it leads to simpler models than those
obtain for instance by simply factoring the Hankel matrix of the observed
data. The effectiveness of the proposed technique was illustrated with two
examples that exhibit chaotic behavior. In principle the approach proposed
here requires solving a large SDP, and it is well known that SDPs have poor
scaling properties. However, as shown in the Appendix, the specific
optimization arising in this paper exhibits an underlying sparse structure
(chordal sparsity) than can be exploited to obtain algorithms whose complexity
scales linearly with the number of data points, when these SDPs are solved
using an ADMM based method such as the one proposed in [14]. This extension
along with an extension to piecewise linear dynamics on the manifold, is
currently being explored.

## References

  * [1] A. C. Antoulas and B. D. O. Anderson.  On the scalar rational interpolation problem.  IMA J. of Mathematical Control and Information, 3:61–88, 1986.
  * [2] Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, Eurika Kaiser, and J. Nathan Kutz.  Chaos as an intermittently forced linear system.  Nature Communications, 8(1):19, 2017.
  * [3] Jerome Dancis.  Positive semidefinite completions of partial hermitian matrices.  Linear Algebra and its Applications, 175:97 – 114, 1992.
  * [4] Fei Xiong, O. I. Camps, and M. Sznaier.  Low order dynamics embedding for high dimensional time series.  In 2011 International Conference on Computer Vision, pages 2368–2374, 2011.
  * [5] Robert Grone, Charles R Johnson, Eduardo M Sá, and Henry Wolkowicz.  Positive definite completions of partial hermitian matrices.  Linear Algebra Appl., 58:109–124, 1984.
  * [6] A. C. Ionita.  Lagrange Rational Interpolation and Its Applications to Approximation of Large-Scale Dynamical System.  PhD thesis, Rice University, 2013.
  * [7] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature Communications, 9(1), 2018.
  * [8] Igor Mezić.  Analysis of Fluid Flows via Spectral Properties of the Koopman Operator.  Annual Review of Fluid Mechanics, 45:357–378, January 2013.
  * [9] Karthik Mohan and Maryam Fazel.  Iterative reweighted algorithms for matrix rank minimization.  Journal of Machine Learning Research, 13(110):3441–3473, 2012.
  * [10] Samuel E. Otto and Clarence W. Rowley.  Linearly recurrent autoencoder networks for learning dynamics.  SIAM Journal on Applied Dynamical Systems, 18(1):558–593, 2019\.
  * [11] Sam T. Roweis and Lawrence K. Saul.  Nonlinear dimensionality reduction by locally linear embedding.  Science, 290(5500):2323–2326, 2000.
  * [12] Floris Takens.  Detecting strange attractors in turbulence.  In David Rand and Lai-Sang Young, editors, Dynamical Systems and Turbulence, Warwick 1980, pages 366–381, Berlin, Heidelberg, 1981. Springer Berlin Heidelberg.
  * [13] F. Xiong, Y. Cheng, O. Camps, M. Sznaier, and C. Lagoa.  Hankel based maximum margin classifiers: A connection between machine learning and wiener systems identification.  In 52nd IEEE Conference on Decision and Control, pages 6005–6010, 2013.
  * [14] Y. Zheng, G. Fantuzzi, A. Papachristodoulou, P. Goulart, and A. Wynn.  Chordal decomposition in operator-splitting methods for sparse semidefinite programs.  Mathematical Programming, 180(1):489–532, 2020.

##  Appendix A Technical Proofs

Proof of Theorem 1. For simplicity assume that the roots of
$\mathcal{P}(\rho)$ are simple. Begin by noting that the $r^{*}$ linearly
independent vectors $\mathbf{v}_{i}$ are in $\mathcal{N}_{R}(\mathbf{p}^{T})$,
the right null space of $\mathbf{p}^{T}$. Since
dim$(\mathcal{N}_{R}(\mathbf{p}^{T}))=r^{*}$, it follows that these vectors
form a basis of $\mathcal{N}_{R}(\mathbf{p}^{T})$, e.g.
$\mathcal{N}_{R}(\mathbf{p}^{T})=\text{span}(\mathbf{V})$. Next, let
${y}_{k,j}^{(i)}$ denote the $j^{th}$ component of $\mathbf{y}_{k}^{(i)}$. By
construction the vectors

$\boldsymbol{\upsilon}_{k,j}^{(i)}\doteq\begin{bmatrix}{y}_{k-r^{*},j}^{(i)}&{y}_{k-r^{*}+1,j}^{(i)}&\ldots&{y}_{k,j}^{(i)}\end{bmatrix}^{T}\in\mathcal{N}_{R}(\mathbf{p}^{T})$

Hence $\boldsymbol{\upsilon}_{k,j}^{(i)}\in\text{span}(\mathbf{V})$ and can be
written as
$\boldsymbol{\upsilon}_{k,j}^{(i)}=\mathbf{V}\mathbf{c}_{k,j}^{(i)}$.
Repeating this reasoning for each component of $\mathbf{y}_{k}^{(i)}$ leads
to:

$\boldsymbol{\psi^{(i)}_{k}}=\left(\mathbf{V}\otimes\mathbf{I}_{m}\right)\mathbf{c}_{k}^{(i)}\doteq\mathbf{D}\mathbf{c}_{k}^{(i)}$

where
$\mathbf{c}_{k}^{(i)}\doteq\sum_{j=1}^{m}\mathbf{c}_{k,j}^{(i)}\otimes\mathbf{e}_{j}$,
$\mathbf{e}_{j}\doteq\begin{bmatrix}0\ldots 1\ldots 0\end{bmatrix}^{T}$ and
where, for notational simplicity we use the shorthand
$\boldsymbol{\psi^{(i)}_{k}}\doteq\boldsymbol{\psi}(\boldsymbol{\xi}^{(i)}_{k})$.
Applying the same reasoning to each (block) component of
$\boldsymbol{\psi^{(i)}_{(k+1)}}$ yields:

$\displaystyle\mathbf{y}_{\ell+1}^{(i)}$
$\displaystyle=(\begin{bmatrix}\rho_{1}^{\ell+1}\ldots\rho_{r^{*}}^{\ell+1}\end{bmatrix}\otimes\mathbf{I}_{m})\mathbf{c}^{(i)}_{k},\;\ell=k-r^{*},\ldots,k$
$\displaystyle=\left((\begin{bmatrix}\rho_{1}^{\ell}\ldots\rho_{r^{*}}^{\ell}\end{bmatrix}\text{diag}(\rho_{i}))\otimes\mathbf{I}_{m}\right)\mathbf{c}^{(i)}_{k}$
$\displaystyle=(\begin{bmatrix}\rho_{1}^{\ell}\ldots\rho_{r^{*}}^{\ell}\end{bmatrix}\otimes\mathbf{I}_{m})\boldsymbol{\Lambda}\mathbf{c}^{(i)}_{k}$

where $\boldsymbol{\Lambda}\doteq\text{diag}(\rho_{i})\otimes\mathbf{I}_{m}$
and where we used the Kronecker’s product property
$(\mathbf{A}\mathbf{C})\otimes(\mathbf{B}\mathbf{D})=(\mathbf{A}\otimes\mathbf{B})(\mathbf{C}\otimes\mathbf{D})$.
Thus

$\displaystyle\boldsymbol{\psi^{(i)}_{k+1}}$
$\displaystyle=\mathbf{D}\boldsymbol{\Lambda}\mathbf{c}^{(i)}_{k}\doteq\mathbf{D}\mathbf{c}^{(i)}_{k+1}\Rightarrow\mathbf{c}^{(i)}_{k+1}=\boldsymbol{\Lambda}\mathbf{c}^{(i)}_{k}$
(30)

It follows that $\boldsymbol{\Lambda}$ propagates $\mathbf{c}^{(i)}_{k}$, the
coordinates of $\boldsymbol{\psi}^{(i)}(\boldsymbol{\xi}_{k})$. Hence the
eigenfunctions $\boldsymbol{\phi}_{j,\ell}(\boldsymbol{\xi})$ of the Koopman
operator have the form
$\boldsymbol{\phi}_{j,\ell}(\boldsymbol{\xi})\doteq\mathbf{v}_{j}\otimes\mathbf{e}_{\ell}$,
$j=1,\ldots,r^{*}$, $\ell=1\ldots m$. ∎

Proof of Theorem 2. By construction,
$(\mathbf{y}_{s}^{(i)})^{T}\mathbf{y}^{(i)}_{t}=\mathbf{K}_{s,t}^{(i)}$, hence
satisfaction of (15)-(16) implies satisfaction of (10)-(11). Consider now the
corresponding Hankel matrices $\mathbf{H}_{\mathbf{y}^{(i)}}$. From the
definitions of $\mathbf{G}^{(i)}$ and $\mathbf{K}_{\ell,r}^{(i)}$ it follows
that
$\mathbf{H}^{T}_{\mathbf{y}^{(i)}}\mathbf{H}_{\mathbf{y}^{(i)}}=\mathbf{G}^{(i)}$.
Since by construction rank($\mathbf{G}^{(i)})$ $\leq r^{*}$ then
rank$(\mathbf{H}_{\mathbf{y}^{(i)}})\leq r^{*}<r+1$ and hence (7) is also
satisfied. ∎

Proof of Theorem 3. Existence of $g(.)$ follows from the properties of the
Loewner matrix discussed in Section 2.2. Since by construction
$\mathbf{K}_{\mathbf{x}}\succeq 0$, it follows that it defines a valid Kernel
in the set $\mathcal{X}\cup\\{\mathbf{x}^{*}\\}$. Further, since
rank$(\mathbf{K_{x}})\leq n$ it can be factored as
$\mathbf{K}_{\mathbf{x}}=\mathbf{X}^{T}\mathbf{X}=\mathbf{K}_{\mathbf{x}}$,
where, again by construction,
$\mathbf{X}(i+1,:)=\mathbf{x}_{i},i=1:|\mathcal{X}|$. Let
$\mathbf{x}^{*}=\mathbf{X}(1,:)$. From the definitions of
$\boldsymbol{\kappa_{\mathbf{x}}},\boldsymbol{\kappa_{\mathbf{y}}}$ it follows
that

$\displaystyle(\mathbf{x}^{*})^{T}\mathbf{x}^{*}=\boldsymbol{\kappa_{\mathbf{x}_{1}}},\;(\mathbf{y}^{*})^{T}\mathbf{y}^{*}=\boldsymbol{\kappa_{\mathbf{y}_{1}}}$
$\displaystyle\mathbf{x}_{i}^{T}\mathbf{x}^{*}=\boldsymbol{\kappa_{\mathbf{x}_{i}}},\;\mathbf{y}_{i}^{T}\mathbf{y}^{*}=\boldsymbol{\kappa_{\mathbf{y}_{i}}}$
$\displaystyle\mathbf{x}_{i}^{T}\mathbf{x}_{i}=\boldsymbol{\kappa_{\mathbf{x}_{j}}},\;\mathbf{y}_{i}^{T}\mathbf{y}_{i}=\boldsymbol{\kappa_{\mathbf{y}_{j}}}$

where $j$ is defined in (24). Thus, the constraint (24) is simply a
restatement of (15) in terms of the elements of
$\boldsymbol{\kappa_{\mathbf{x}}},\boldsymbol{\kappa_{\mathbf{y}}}$. ∎

##  Appendix B Exploting Chordal Sparsity

###  B.1 Semi-Definite Programs and Rank Minimization Over Chordal Graphs

In this paper, we will reduce the problem of identifying Koopman operators to
a constrained rank minimization of the form

$\begin{split}&\min_{\mathbf{X}\succeq 0}\text{rank $(\mathbf{X})$ subject to
}\\\ &\text{Trace $(\mathbf{A}_{i}\mathbf{X})\leq b_{i},i=1,\ldots,n_{c}$,
$\mathbf{X},\mathbf{A}_{i}\in\mathbb{R}^{n\times n}$}\end{split}$ (31)

In the specific problems arising in this paper only a small number of entries
of $\mathbf{X}$ appear in the trace constraints, while the role of the other
entries is just to enforce that $\mathbf{X}\succeq 0$. Thus, as long as
existence of a _minimum rank_ PSD completion is guaranteed, these variables do
not have to be explicitly found, allowing for a substantial computational
complexity reduction. Specifically, to the optimization (31) one can associate
a graph $\mathcal{G}(\mathcal{V,E})$ with $n$ vertices in $\mathcal{V}$ and
edge set $\mathcal{E}$, where there is an edge between vertices $j$ and $\ell$
if the element $(j,\ell)$ of any of the matrices $\mathbf{A}_{i}$ is nonzero.
Given a graph $\mathcal{G}$, define the cone

$\mathbb{S}_{+}^{n}(\mathcal{E},?)\doteq\\{\mathbf{X}\in\mathbb{S}_{+}^{n}\colon\text{$\mathbf{X}_{i,j}$
given if $(i,j)\in\mathcal{E}$}\\}$

that is, the cone of matrices with entries fixed over the edges $\mathcal{E}$
than can be completed to be PSD. When the graph $\mathcal{G}$ is chordal, the
minimum rank over all possible matrix completions over this cone has an
explicit expression, given by Dancis’ Theorem:

######  Theorem 4 ([3]).

Let $\mathcal{G}(\mathcal{V},\mathcal{E})$ be a chordal graph with a set of
maximal cliques
$\\{\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{n_{c}}\\}$. Then, for
any $\mathbf{X}\in\mathbb{S}^{n}_{+}(\mathcal{E},?)$ there exist at least one
minimum rank PSD completion where

$\text{rank}(\mathbf{X})=\max_{1\leq k\leq
n_{c}}\text{rank}(\mathbf{E}_{\mathcal{C}_{k}}\mathbf{X}\mathbf{E}_{\mathcal{C}_{k}}^{T})$

where the $0/1$ matrix $\mathbf{E}_{\mathcal{C}_{k}}$ selects the variables of
$\mathbf{X}$ corresponding to edges in the clique $\mathcal{C}_{k}$.

In addition, the cone $\mathbb{S}_{+}^{n}(\mathcal{E},?)$ can be characterized
using the following result (Grone’s Theorem):

######  Theorem 5 ([5]).

Let $\mathcal{G}(\mathcal{V},\mathcal{E})$ be a chordal graph with a set of
maximal cliques
$\\{\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{n_{c}}\\}$. Then,
$\mathbf{X}\in\mathbb{S}^{n}_{+}(\mathcal{E},?)$ if and only if

$\mathbf{X}_{k}=\mathbf{E}_{\mathcal{C}_{k}}\mathbf{X}\mathbf{E}_{\mathcal{C}_{k}}^{T}\in\mathbb{S}_{+},\qquad
k=1,\,\ldots,\,n_{c}$

Combining the two theorems above leads to the following result.

######  Corollary 1.

The optimization (31) is equivalent to:

$\begin{split}&\min\sum_{k}\text{rank
($\mathbf{E}_{\mathcal{C}_{k}}\mathbf{X}\mathbf{E}_{\mathcal{C}_{k}}^{T}$)
subject to }\\\
&\mathbf{E}_{\mathcal{C}_{k}}\mathbf{X}\mathbf{E}_{\mathcal{C}_{k}}^{T}\succeq
0\\\ &\text{Trace $(\mathbf{A}_{i}\mathbf{X})\leq
b_{i},i=1,\ldots,n_{c}$}\end{split}$ (32)

Since rank minimization problems are generically NP-hard, a standard convex
relaxation is to replace rank by is convex envelope, trace [9]. This
substitution leads to a convex SDP that can be solved to $\epsilon$–optimality
in polynomial time using interior point (IP) methods. However, while
efficient, these methods have relatively poor scaling properties
($\mathcal{O}(n^{2}n_{c}^{2}+n^{3}n_{c})$). On the other hand, as we will show
in the sequel, the specific problem arising in this paper has chordal
sparsity. Hence, the use of Corollary 1 to decompose the objective into
$\sum_{k}\text{Trace}(\mathbf{E}_{\mathcal{C}_{k}}\mathbf{X}\mathbf{E}_{\mathcal{C}_{k}}^{T})$
leads to a SDP where each of the PSD constraints has (size of
$\mathcal{C}_{k}$) variables. When combined with an ADMM approach where the
cost of each iteration is $\mathcal{O}(\text{number of variables}^{3})$ [14],
using this decomposition leads to a reduction of
$\frac{1}{n_{c}}(\frac{n}{\text{max size of clique}})^{3}$ in computational
complexity.

###  B.2 Exploiting Chordal Sparsity

The approach outlined in Section 3.3 works well for small to medium sized
problems. However, its computational complexity grows roughly as
$\mathcal{O}(\text{(number of data points)}^{6})$222The number of free
variables in $\mathbf{K}$ is $(\text{number of data points)}^{2}$. Thus, if
the SDP is solved using an interior point method, computational complexity
scales roughly as $\mathcal{O}(\text{(number of variables)}^{6})$..
Fortunately, as we show next, the convex relaxation of (12)-(16) is endowed
with chordal sparsity. Hence the decomposition outlined in Section B.1,
combined with an ADMM based algorithm such as the one proposed in [14], can be
exploited to substantially reduce computational complexity. Let $r$ be an
upper bound of the optimal rank $r^{*}$. Note that the only elements of
$\mathbf{K}^{(i)}$ that appear explicitly in (12)-(16) are those of the form
$\mathbf{K}^{(i)}_{s,t}$ where either $|s-t|\leq r$ or
$\|\mathbf{x}_{s}-\mathbf{x}_{r}\|_{2}\leq\delta$. Let

$\begin{split}&\mathcal{T}_{\ell}\doteq\\{(s,t):\ell\leq s,t\leq r\\}\\\
&\mathcal{S}_{\ell}\doteq\\{(s,q):\text{$\ell\leq s\leq r$ and
$\|\mathbf{x}_{s}-\mathbf{x}_{q}\|_{2}\leq\delta$}\\}\\\
&\mathcal{E}_{\ell}\doteq\mathcal{T}_{\ell}\cup\mathcal{S}_{\ell}\end{split}$
(33)

To the optimization problem (12)-(16) we can associate a graph with cliques
$\mathcal{C}_{\ell}$ defined by the edge sets $\mathcal{E}_{\ell}$. From
Corollary 1 and the fact that only the variables in $\mathcal{C}_{\ell}$
appear in the objective (12), it follows that each constraint
$\mathbf{K}^{(i)}\succeq 0$ can be replaced by a collection of smaller
constraints of the form
$\mathbf{E}^{T}_{\mathcal{C}_{\ell}}\mathbf{K}^{(i)}\mathbf{E}_{\mathcal{C}_{\ell}}\succeq
0$, where the matrix $\mathbf{E}_{\mathcal{C}_{\ell}}$ selects the entries of
$\mathbf{K}^{(i)}$ corresponding to edges in $\mathcal{C}_{\ell}$. Assuming a
fixed number $n_{v}$ of spacial neighbours, the size of each clique
$\mathcal{C}_{\ell}$ is given by $|\mathcal{C}_{\ell}|=(r+1)(1+n_{v})$ and
each trajectory has $T_{i}$ cliques. It follows that the computational
complexity when using the clique-based decomposition roughly decreases by a
factor of $\frac{(1+r)^{3}(1+n_{v})^{3}}{T_{i}}$. It is worth noting that,
when using the clique decomposition, the overall computational complexity
increases as $T_{i}(1+r)^{3}(1+n_{v})^{3}$. This scaling is linear, rather
than polynomial, in the number of data points.