# Data-driven Discovery of Invariant Measures

Jason J. Bramburger1, Giovanni Fantuzzi2

( 1Department of Mathematics and Statistics, Concordia University,
Montréal, QC, Canada
2Department of Mathematics, Friedrich–Alexander–Universität Erlangen–Nürnberg,
Germany )

###### Abstract

Invariant measures encode the long-time behaviour of a dynamical system. In
this work, we propose an optimization-based method to discover invariant
measures directly from data gathered from a system. Our method does not
require an explicit model for the dynamics and allows one to target specific
invariant measures, such as physical and ergodic measures. Moreover, it
applies to both deterministic and stochastic dynamics in either continuous or
discrete time. We provide convergence results and illustrate the performance
of our method on data from the logistic map and a stochastic double-well
system, for which invariant measures can be found by other means. We then use
our method to approximate the physical measure of the chaotic attractor of the
Rössler system, and we extract unstable periodic orbits embedded in this
attractor by identifying discrete-time periodic points of a suitably defined
Poincaré map. This final example is truly data-driven and shows that our
method can significantly outperform previous approaches based on model
identification.

##  1 Introduction

Dynamical system analysis traditionally takes one of two complementary
approaches. The first is to study individual system trajectories in order to
identify fundamental structures such as steady-states, periodic orbits,
heteroclinic connections, and chaotic attractors. The second approach is to
understand the statistical behaviour of the system by examining the evolution
of measures under the dynamics. Paramount in the study of measure-theoretic
properties of dynamical systems is the identification of _invariant_ measures,
which are preserved by the flow of the system. Invariant measures lie at the
centre of ergodic theory and provide information regarding the long-time
behaviour of a system. Indeed, ergodic theorems allow one to calculate
infinite-time averages as expectations with respect to invariant measures [7],
transferring one from the realm of dynamics to one of probability theory.

In this paper, we describe an optimization-based approach to extract invariant
measures from dynamic data. Such a problem goes back at least to the work of
Ulam [35], who proposed discretizing the phase space into a finite collection
of disjoint sets and using the data to approximate transition probabilities
between these sets [17, 21, 23, 11]. This method produces a square stochastic
matrix that approximates the Perron–Frobenius operator, that is, the linear
operator describing how probability measures evolve under the system dynamics
(see 2.5 below for a precise definition). By the Perron–Frobenius theorem,
this matrix has a positive eigenvector with eigenvalue 1 that approximates the
density of an invariant measure of the system. However, the computational
complexity of Ulam’s original method grows cubically with the dimension of the
phase space and, even though recent variations of the method exhibit quadratic
scaling [14], it quickly becomes prohibitive. Moreover, accurate
approximations of invariant measures require fine discretizations of the phase
space and a considerable amount of data.

An alternative approach to approximating invariant measures from data rests on
the fact that the Perron–Frobenius operator is the formal adjoint of the
Koopman operator. The Koopman operator is linear and acts on real-valued
functions, often called _observables_ , through composition with the flow of a
dynamical system. Its action can be approximated on finite collections of
observables, termed _dictionaries_ , using extended dynamic mode decomposition
(EDMD) [36]. This method produces a square matrix that approximates the
Koopman operator and any left eigenvector with eigenvalue close to 1
approximates the density of an invariant measure of the underlying dynamical
system [21, 23, 17]. Despite EDMD’s rigorous convergence guarantees [4, 25,
36], however, the existence of nonnegative eigenvectors with unit eigenvalues
cannot be guaranteed. Thus, contrary to Ulam’s method, one must either settle
for signed densities coming from eigenvalues near 1, or augment EDMD to
produce a column-wise stochastic approximate Koopman operator matrix [14].

In this work, we propose to approximate invariant measures by combining EDMD
with convex optimization. Precisely, we first employ a variation of the EDMD
method from [4] to produce a rectangular (rather than square) matrix
approximation for the action of the Koopman operator on a dictionary of
polynomial observables. Using this approximation, we then provide a data-
driven implementation of the convex optimization framework for studying
invariant measures developed in [24] for dynamical systems with known
polynomial governing equations (see also [2, 1] for a different model-based
implementation of the same framework, which allows for non-polynomial
governing equations at the expense of other assumptions on the system’s
invariant measures). This strategy has two major advantages compared both to
Ulam’s method and to alternative EDMD-based techniques. First, we can easily
target different invariant measures if more than one exists for the dynamical
systems generating the data. In particular, we can search for _ergodic_
invariant measures, which are extreme points of the set of all invariant
measures in the sense of convex analysis. Second, we are not limited to
discovering densities of absolutely continuous invariant measures. In
particular, we are able to discover singular measures supported on invariant
trajectories for the underlying dynamical system.

The ability of our approach to discover extreme singular invariant measures
from data opens up new possibilities for the analysis of dynamical systems.
For example, extreme singular invariant measures for chaotic discrete
dynamical systems are often supported on fixed points or on unstable periodic
orbits (UPOs) that densely fill the system’s attractor. Thus, our approach
provides a way of systematically extracting UPOs from chaotic maps. We
demonstrate this by extracting both short- and long-period UPOs for the
continuous-time chaotic Rössler system using only data from a one-dimensional
Poincaré section. Precisely, we first use our data-driven technique to extract
extreme atomic measures supported on UPOs of the Poincaré return map of period
up to 19. Then, we recover the corresponding continuous-time UPOs for the
Rössler system. This strategy outperforms data-driven methods that first
identify a model for the Poincaré map data [5, 6, 3], because small
inaccuracies in the identified model propagate to forward iterations and make
it difficult to discover long-period UPOs. There are also advantages compared
to related methods for extracting UPOs from approximate invariant measures in
continuous-time [29, 28], as these rely on the target UPO being well
approximated by the zero level set of a nonnegative polynomial. This is
usually not the case unless the polynomial degree is large, because
continuous-time UPOs are not algebraic curves in general. In contrast, the
corresponding UPO for the Poincaré map generates an atomic measure, whose
atoms can be represented exactly as the zeros of a nonnegative polynomial. The
Poincaré map formulation is therefore better suited to algebraic methods and
our data-driven approach allows one to implement it even when, as is typical,
the Poincaré map is not known.

To make this article self-contained, section 2 introduces invariant measures
and reviews in detail the optimization framework for their approximation
developed in [24]. Section 3 shows how to implement this framework in a data-
driven setting by using EDMD, while convergence results for dynamical systems
governed by unknown polynomial maps are proved in section 4. In these sections
we focus on discrete-time deterministic dynamics for simplicity, but we
comment on extensions to continuous-time and stochastic dynamics in section 5.
Section 6 demonstrates the power of our data-driven approach on a series of
examples. Final remarks and open questions are offered in section 7.

##  2 Invariant Measures

We begin by reviewing invariant measures for discrete-time dynamical systems,
as well as existing model-based optimization tools for their approximation. As
mentioned above, these tools are well known (see, e.g., [18, 9, 13, 27, 32,
24, 2, 1]). Here we closely follow [24] and refer readers to this work for
further details and convergence proofs.

###  2.1 Background and Definitions

Let $X\subset\mathbb{R}^{n}$ be a compact set and, given a continuous function
$f:X\to X$, consider the discrete-time dynamical system governed by

$x_{n+1}=f(x_{n}).$ (2.1)

A nonnegative Borel measure $\mu$ supported on $X$ is _invariant_ for the
mapping $f$ if

$\mu(f^{-1}(A))=\mu(A)$ (2.2)

for all Borel measurable sets $A\subseteq X$. Equivalently, $\mu$ is invariant
if

$\int_{X}[g\circ f-g]\mathrm{d}\mu=0$ (2.3)

for all functions $g\in C(X)$. We assume throughout that $\mu$ is in the space
$\Pr(X)$ of probability measures on $X$, meaning $\mu(X)=1$. The set of all
invariant probability measures for the mapping $f$ is convex and its extreme
points are the ergodic measures.

Condition 2.3 has a deep connection to the Koopman operator and its adjoint,
the Perron–Frobenius operator. The Koopman operator associated to $f$ is the
linear operator $\mathcal{K}:C(X)\to C(X)$ encoding the composition of
continuous functions with $f$, i.e.,

$\mathcal{K}g:=g\circ f.$ (2.4)

The Perron–Frobenius operator, instead, is a linear operator $\mathcal{P}$ on
the space of measures that associates a measure $\mu$ to its pushforward by
$f$, i.e.,

$(\mathcal{P}\mu)(A)=\mu(f^{-1}(A)).$ (2.5)

It should be clear that 2.3 can be written in terms of the Koopman operator as

$\int_{X}[\mathcal{K}g-g]\mathrm{d}\mu=0,$ (2.6)

while 2.2 can be expressed as $\mathcal{P}\mu=\mu$. In other words, invariant
measures are eigenfunctions of the Perron–Frobenius operator with eigenvalue
$1$.

###  2.2 Convex Computation of Invariant Measures

Invariant measures for a known polynomial map $f$ can be approximated using a
convex optimization framework developed in [24], which we now review in detail
as it forms the basis of our data-driven method.

Fix an integer $d\geq 1$. For this section, and this section only, we assume
the following.

######  Assumption 2.1.

The function $f:X\to X$ is a polynomial of total degree $\leq d$.

We also introduce the following assumption about the set $X$.

######  Assumption 2.2.

$X\subset\mathbb{R}^{n}$ is a compact set such that:

  1. 1.

$X=\\{x\in\mathbb{R}^{n}|\ \sigma_{1}(x)\geq 0,\ldots,\sigma_{J}(x)\geq 0\\}$
for finitely many polynomials $\sigma_{1},\dots,\sigma_{J}$ of degree $\leq
d$.

  2. 2.

There exists an integer $R$ and polynomials $s_{0},\ldots,s_{J}$ that are sums
of squares of other polynomials such that
$R^{2}-|x|^{2}=s_{0}(x)+s_{1}(x)\sigma_{1}(x)+\cdots+s_{J}(x)\sigma_{J}(x)$.

Sets satisfying the first condition in assumption 2.2 are called basic
semialgebraic. Examples include boxes and balls, but the definition is
sufficiently general to include more complicated shapes, both convex and non-
convex. The second condition in assumption 2.2, instead, is a technical
requirement for the convergence results stated in theorem 2.1 below. This
requirement is mild when $X$ is compact, because one can always add the
inequality $R^{2}-|x|^{2}\geq 0$ to the definition of $X$ with $R$
sufficiently large not to change the set.

####  2.2.1 Moment-based Characterization of Invariance

Since $X$ is compact, any measure $\mu$ on $X$ is uniquely determined by its
moments

$y_{\alpha}=\int_{X}x^{\alpha}\mathrm{d}\mu,$ (2.7)

where $x^{\alpha}=x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{n}}$ is the monomial
in $x$ with exponent
$\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n}$. We write
$|\alpha|=\alpha_{1}+\cdots+\alpha_{n}$ for the degree of $x^{\alpha}$. The
first step to study invariant measures via convex optimization is to encode
the invariance of a measure using linear constraints on its moments.

Since polynomials are dense in $C(X)$, we conclude from 2.3 that a measure
$\mu$ supported on $X$ is invariant if and only if

$\int_{X}\left[(f(x))^{\alpha}-x^{\alpha}\right]\mathrm{d}\mu=0$ (2.8)

for every $\alpha\in\mathbb{N}_{0}^{n}$. Assumption 2.1 implies that
$(f)^{\alpha}$ is a polynomial of degree $|\alpha|d$. Thus, for each exponent
$\alpha$, condition 2.8 provides a linear constraint involving the (finitely
many) moments of $\mu$ of degree up to $|\alpha|d$. Letting $\alpha$ vary over
$\mathbb{N}_{0}^{n}$, we conclude that $\mu$ is an invariant measure if and
only if the corresponding moment sequence
$y=\\{y_{\alpha}\\}_{\alpha\in\mathbb{N}_{0}^{n}}$ satisfies an infinite set
of linear equality constraints. We write such constraints as

$\mathcal{A}y=0,$ (2.9)

where $\mathcal{A}:\mathbb{R}^{\infty}\to\mathbb{R}^{\infty}$ is an infinite-
dimensional linear operator. We also require that

$y_{0}=1$ (2.10)

because we focus on probability measures.

####  2.2.2 Optimization over Invariant Measures

We have seen above that finding an invariant measure is tantamount to finding
a sequence $y\in\mathbb{R}^{\infty}$ whose elements are moments of a measure
and satisfy 2.9 and 2.10. These conditions have multiple solutions if the
dynamical system 2.1 admits multiple invariant measures. To try and single out
an invariant measure with specific properties, one can optimize a cost
function $F(y)$. Following [24], we let $F(y)$ be a convex and lower
semicontinuous function that depends only on the first $\binom{n+d}{d}$
entries of $y$, which correspond to the moments of degree up to $d$. This
leads to the convex optimization problem

$\min_{\begin{subarray}{c}y\in\mathbb{R}^{\infty}\\\
\mu\in\Pr(X)\end{subarray}}\ F(y)\quad\text{such
that}\quad\begin{cases}\mathcal{A}y=0,\\\ y_{0}=1,\\\
y_{\alpha}=\int_{X}x^{\alpha}\,\mathrm{d}\mu\quad\forall\alpha\in\mathbb{N}_{0}^{n}.\end{cases}$
(2.11)

Examples of useful cost functions are discussed in [24]. Here, we will
consider only cost functions that promote the identification of physical
measures and ergodic measures.

##### Physical Measures.

A measure $\mu$ is said to be _physical_ if

$\int_{X}g\,\mathrm{d}\mu=\lim_{N\to\infty}\frac{1}{N}\sum_{n=0}^{N-1}g(f^{n}(x))$
(2.12)

for Lebesgue-almost-every point $x\in X$ and every $g\in C(X)$. Taking
$g(x)=x^{\alpha}$ and fixed $N\gg 1$, the moments of a physical measure $\mu$
satisfy

$y_{\alpha}\approx\frac{1}{N}\sum_{n=0}^{N-1}g(f^{n}(x_{0}))=:\tilde{y}_{\alpha}$
(2.13)

for Lebesgue-almost-every point $x_{0}\in X$. Then, to identify $\mu$, one can
compute $\tilde{y}_{\alpha}$ for some set $\mathcal{I}$ of exponents $\alpha$
with degree $|\alpha|\leq d$ and solve 2.11 with cost function

$F(y)=\sum_{\alpha\in\mathcal{I}}(y_{\alpha}-\tilde{y}_{\alpha})^{2}.$ (2.14)

##### Ergodic Measures.

An invariant probability measure $\mu\in\Pr(X)$ for 2.1 is said to be
_ergodic_ if every invariant set has either full or zero measure, that is,
$f^{-1}(S)=S$ implies $\mu(S)\in\\{0,1\\}$. Ergodic measures are extreme
points of the set of all invariant probability measures of $f$. Since the
minima of linear programs are attained at extreme points, one can target an
ergodic measure by solving 2.11 with a linear cost function $F$.

####  2.2.3 Approximation via Semidefinite Programming

The minimization problem in 2.11 is convex and returns an optimal invariant
measure for a given cost function, but it is infinite-dimensional and thus
computationally intractable. We now explain how to approximate it with a
hierarchy of semidefinite programs (SDP), which are tractable optimization
problems with positive semidefiniteness constraints on matrices that depend
linearly on the optimization variables.

The first step is to replace the infinite sequence $y$ with a finite vector,
also denoted $y$ with a slight abuse of notation. Fix an integer $k\geq 1$.
Rather than the infinite set of equality constraints in 2.9, which are
obtained by imposing 2.8 for all exponents $\alpha\in\mathbb{N}_{0}^{n}$, we
enforce this condition only for the $\binom{n+k}{k}$ exponents with degree
$|\alpha|\leq k$. This gives $\binom{n+k}{k}$ constraints for the vector $y$
of $\binom{n+kd}{kd}$ moments of $\mu$ of degree up to $kd$, written compactly
as

$A_{k}y=0$ (2.15)

for an explicitly computable matrix
$A_{k}\in\mathbb{R}^{\binom{n+k}{k}\times\binom{n+dk}{dk}}$. Observe that
$A_{k}$ has a high-dimensional kernel when $k$ is large and $f$ is nonlinear
($d\geq 1$).

Next, we seek to identify elements in the kernel of $A_{k}$ that correspond to
truncated moment sequences of probability measures on $X$. Unfortunately, this
problem is intractable in general. Following a standard relaxation strategy
for moment problems (see [30, 10] and references therein), we therefore limit
ourselves to identifying kernel elements $y$ that make the _moment matrices_
$M(\sigma_{j}y)$, $j=0,1,\dots,J$, positive definite. Setting $\sigma_{0}=1$,
these matrices are defined by

$M(\sigma_{j}y)=\int_{X}\sigma_{j}(x)p_{\gamma_{j}}(x)p_{\gamma_{j}}(x)^{\intercal}\mathrm{d}\mu,\qquad
j=0,1,\dots,J,$ (2.16)

where $\gamma_{j}=\lfloor(dk-\deg\sigma_{j})/2\rfloor$ and $p_{\gamma}(x)$ is
the vector of monomials in $x$ up to degree $\gamma$. The square matrix
$M(\sigma_{j}y)$ is therefore of size
$\binom{n+\gamma_{j}}{\gamma_{j}}\times\binom{n+\gamma_{j}}{\gamma_{j}}$. We
also impose the normalization condition $y_{0}=1$ to arrive at the following
tractable _semidefinite program_ (SDP):

$\min_{y\in\mathbb{R}^{\binom{n+k}{k}}}\ F(y)\quad\text{such
that}\quad\begin{cases}A_{k}y=0\\\ y_{0}=1\\\ M(\sigma_{j}y)\succeq
0,&j=0,1,\dots,J.\end{cases}$ (2.17)

This SDP is a relaxation of 2.11 in the sense that its feasible set contains
all moment vectors of all invariant measures for the discrete-time system 2.1.
In general, however, it also contains elements that are not moment vectors, so
the optimal $y$ need not be a vector of moments of an invariant measure.
Nevertheless, one can prove that optimal solutions to 2.17 converge (up to a
subsequence) to moment sequences of optimal invariant measures for 2.11.

######  Theorem 2.1 ([24, Theorem 2]).

Suppose assumptions 2.1 and 2.2 hold. Let
$y^{k}\in\mathbb{R}^{\binom{n+kd}{kd}}$ be the optimal solution of 2.17 for
fixed $k\in\mathbb{N}$, and let $\mu$ be optimal for 2.11. There exists an
increasing sequence $\\{k_{j}\\}_{j\geq 1}\subset\mathbb{N}$ such that

$\lim_{j\to\infty}y^{k_{j}}_{\alpha}=\int_{X}x^{\alpha}\,\mathrm{d}\mu\qquad\forall\alpha\in\mathbb{N}_{0}^{n}.$

In particular, if 2.11 has a unique optimizer $\mu$, then
$y_{\alpha}^{k}\to\int x^{\alpha}\,\mathrm{d}\mu$ as $k\to\infty$.

Finally, we remark that the SDP 2.17 can be formulated and solved even when
the set $X$ does not satisfy the second condition in assumption 2.2 and, in
particular, when $X=\mathbb{R}^{n}$ is the full space. In such cases, however,
it is currently not know if the convergence results from theorem 2.1 always
apply. Fortunately, convergence if often observed in examples (see, e.g., [9,
13, 15, 16]).

####  2.2.4 Recovery of an Approximate Invariant Measure

Solutions of the SDP 2.17 for increasing $k$ recover, possibly after passing
to a subsequence, the moments of an optimal invariant measure $\mu$ for 2.11.
We now review how to approximate $\mu$ using a sequence of _signed_ measures
$\mu_{k}$ constructed from the optimal solutions $y^{k}$ of 2.17. For
notational ease, we assume we have already passed (without relabelling) to a
convergent subsequence of optimal solutions as in theorem 2.1.

First, fix a finite measure $\pi$ (not necessarily a probability measure)
whose support is contained in $X$ and such that the associated moment matrix

$M_{\pi}:=\int_{X}p_{k}(x)p_{k}(x)^{\intercal}\,\mathrm{d}\pi$ (2.18)

is computable for every integer $k$. We assume $\pi$ is supported on a full-
dimensional or non-algebraic subset of $X$, so as to ensure $M_{\pi}$ is
invertible. In the examples of section 6 we will take $\pi$ to be the Lebesgue
measure on $X$, but here we work with general $\pi$ to highlight that other
choices are possible.

Next, consider the degree-$k$ polynomial

$\rho_{k}(x)=p_{k}(x)^{\intercal}\,M_{\pi}^{-1}y^{k}$ (2.19)

and define a $\pi$-absolutely-continuous signed measure $\mu_{k}$ via

$\mu_{k}(A)=\int_{A}\rho_{k}(x)\,\mathrm{d}\pi\qquad\forall A\subset X.$
(2.20)

This measure satisfies $\mu_{k}(X)=1$ by construction. Note, however, that
$\mu_{k}$ is a _signed_ measure in general because $\rho_{k}$ need not be
nonnegative on the support of $\pi$. Nevertheless, a straightforward extension
of arguments in [24] shows that $\mu_{k}$ converges to an optimal invariant
measure $\mu$ for 2.11 in the weak-star topology of $\Pr(X)$, meaning that

$\lim_{k\to\infty}\int_{X}g(x)\rho_{k}(x)\,\mathrm{d}\pi=\int_{X}g(x)\,\mathrm{d}\mu$
(2.21)

for every continuous function $g:X\to\mathbb{R}$. In other words, as the
polynomial degree $k$ increases, the signed measures $\mu_{k}$ enable one to
approximate as accurately as desired the expectation with respect to the
invariant measure $\mu$ of any fixed function $g$. Moreover, 2.21 implies that
if $\mu$ is $\pi$-absolutely-continuous with density $\rho$, then the signed
densities $\rho_{k}$ converge to $\rho$ in the weak topology of $L^{1}(\pi)$.

##  3 Discovering Extremal Invariant Measures from Data

In the convex optimization approach to identifying invariant measures reviewed
in section 2, the key ingredient is the matrix $A_{k}$. This matrix is the
only component of the optimization problem 2.17 that encodes the system
dynamics and it is built by considering the action of the Koopman operator,
$\mathcal{K}g=g\circ f$, on monomials $g(x)=x^{\alpha}$ of degree
$|\alpha|\leq k$. Recent advances in data-driven analysis of dynamical systems
make it possible to approximate the action of the Koopman operator on finite-
dimensional spaces of functions, called _observables_ , using only data
gathered from the system [36, 23, 22]. Here, we show how to use these data-
driven methods together with the computational approach from section 2 to
discover invariant measures entirely from data. We henceforth assume we do not
know the map $f$ in 2.1, but we have a set of $m$ data snapshots
$\\{(x_{i},z_{i})\\}_{i=1}^{m}\subset X\times X$ where $z_{i}=f(x_{i})$. We do
not assume that $f$ is polynomial, but still require $X$ to be a basic
semialgebraic set according to the first condition in assumption 2.2.

Fix integers $k$ and $\ell$ with $1\leq k\leq\ell$, and let $p_{k}(x)$ and
$p_{\ell}(x)$ be the vectors of monomials in $x$ up to degrees $k$ and $\ell$,
respectively. Let $k_{x}=\binom{n+k}{k}$ and $\ell_{x}=\binom{n+\ell}{\ell}$
denote the size of these vectors, respectively, and consider the matrices

$\displaystyle P$ $\displaystyle=\begin{bmatrix}|&&|\\\
p_{\ell}(x_{1})&\cdots&p_{\ell}(x_{m})\\\
|&&|\end{bmatrix}\in\mathbb{R}^{\ell_{x}\times m}$ (3.1) and $\displaystyle Q$
$\displaystyle=\begin{bmatrix}|&&|\\\ p_{k}(z_{1})&\cdots&p_{k}(z_{m})\\\
|&&|\end{bmatrix}\in\mathbb{R}^{k_{x}\times m}.$ (3.2)

Extended dynamic mode decomposition (EDMD) seeks to approximate the action of
the Koopman operator on $\mathrm{span}\\{p_{k}(x)\\}$ using a linear operator
$\mathcal{K}_{k\ell}:\mathrm{span}\\{p_{k}(x)\\}\to\mathrm{span}\\{p_{\ell}(x)\\}$
called the approximate Koopman operator. The matrix $K$ representing this
operator is obtained by minimizing the Frobenius matrix norm $\|Q-KP\|_{F}$.
The optimal solution is

$K_{k\ell}:=QP^{\dagger}\in\mathbb{R}^{k_{x}\times\ell_{x}}$ (3.3)

where the superscript $\dagger$ denotes the Moore–Penrose pseudoinverse. The
approximate Koopman operator then acts on elements $g(x)=c\cdot p_{k}(x)$ in
$\mathrm{span}\\{p_{k}(x)\\}$ by

$\mathcal{K}_{k\ell}\,g(x):=c\cdot K_{k\ell}\,p_{\ell}(x).$ (3.4)

Next, let $\Theta_{k\ell}\in\mathbb{R}^{k_{x}\times\ell_{x}}$ be the matrix
that extracts the monomials in $p_{k}(x)$ from $p_{\ell}(x)$, meaning that
$p_{k}(x)=\Theta_{k\ell}p_{\ell}(x)$. Define the matrix

$L_{k\ell}:=K_{k\ell}-\Theta_{k\ell}\in\mathbb{R}^{k_{x}\times\ell_{x}}.$
(3.5)

The approximate Lie derivative $\mathcal{L}_{k\ell}$ is a linear operator from
$\mathrm{span}\\{p_{k}(x)\\}$ to $\mathrm{span}\\{p_{\ell}(x)\\}$ that acts on
any function $g(x)=c\cdot p_{k}(x)$ in $\mathrm{span}\\{p_{k}(x)\\}$ by

$\mathcal{L}_{k\ell}\,g(x)=c\cdot
L_{k\ell}\,p_{\ell}(x)=\mathcal{K}_{k\ell}\,g(x)-g(x).$ (3.6)

Using this definition, condition 2.6 for the invariance of a measure $\mu$ can
be approximated by

$\int_{X}\mathcal{L}_{k\ell}g\,\mathrm{d}\mu=0\qquad\forall
g\in\mathrm{span}\\{p_{k}(x)\\}.$ (3.7)

Note that this approximation is built purely from the data snapshots and
requires no knowledge of the map $f$ driving the dynamics.

Condition 3.7 is equivalent to requiring that the vector
$y\in\mathbb{R}^{\ell_{x}}$ of moments of $\mu$ with degree up to $k$
satisfies $L_{k\ell}\,y=0$. This finite set of constraints is a data-driven
approximation of the condition $A_{k}y=0$ in 2.17. This suggests that an
approximation to an invariant measure $\mu$ for 2.1 may be discovered from the
data snapshots by first solving the SDP

$\min_{y\in\mathbb{R}^{\ell_{x}}}\ F(y)\quad\text{such
that}\quad\begin{cases}L_{k\ell}\,y=0\\\ y_{0}=1\\\ M(\sigma_{j}y)\succeq
0,&j=0,1,\dots,J,\end{cases}$ (3.8)

and then applying the measure recovery strategy outlined in section 2.2.4.
Once again, observe that the SDP in 3.8 is constructed only using data
gathered from the system and without requiring any knowledge of the underlying
map $f$. In particular, one can construct it and solve it even when the
dynamics are not governed by polynomial $f$.

######  Remark 3.1.

Contrary to standard practice, we implement EDMD using two different
polynomial dictionaries $p_{k}$ and $p_{\ell}$. This is useful because the
space of degree-$k$ polynomials is not generally invariant under the action of
the Koopman operator, so $\mathcal{K}g$ for a polynomial $g$ of degree $k$ can
be better approximated using polynomials of degree $\ell\geq k$. If the data
snapshots are governed by a polynomial map of known degree $d\geq 1$, then it
would suffice to set $\ell=dk$ as in section 2. When the map $f$ is unknown,
however, one should implement the method with fixed $k$ and increasingly large
values of $\ell$ until good performance is observed. This can be expected
given theoretical analysis in [4], which under mild conditions guarantees the
convergence of the approximate Lie derivative $\mathcal{L}_{k\ell}$ to the
true Lie derivative $\mathcal{L}$ acting on $\mathrm{span}\\{p_{k}(x)\\}$ as
both the number $m$ of data snapshots and the degree $\ell$ tend to infinity.
More precise results for the case of unknown polynomial maps are presented in
section 4.

######  Remark 3.2.

Above, as well as in section 2, one can replace the monomial vectors
$p_{k}(x)$ and $p_{\ell}(x)$ with any bases for the spaces of polynomials in
$x$ of degree $k$ and $\ell$, respectively. Moments are then expressed with
respect to these bases. Choosing a good polynomial basis is key to ensure the
SDPs 2.17 and 3.8 are well conditioned. For dynamical systems in one
dimension, for instance, it is usually convenient to use the Chebyshev basis.
Unfortunately, we are not aware of a basis selection strategy that works well
in general.

##  4 Convergence for Polynomial Maps

We now investigate the convergence of the data-driven approach introduced in
the previous section to the model-based approach of [24] when the data points
$\\{(x_{i},z_{i})\\}_{i=1}^{m}$ are generated by a polynomial map $f$ of
degree $d$. We are interested in particular in clarifying how solutions of the
data-driven SDP 3.8 behave as the number $m$ of datapoints becomes infinite.

Loosely speaking, one expects that optimal solutions of our data-driven SDP
should converge to optimal solutions to the model-based SDP 2.17. However,
these two problems have a different number of variables unless the polynomial
degrees $k$ and $\ell$ satisfy $\ell=dk$, so they cannot be compared directly.
Instead, under suitable conditions, we will show that for every $\ell\geq dk$
the optimal solutions of the data-driven SDP 3.8 converge (up to subsequences)
to an optimal solution of the ”limit” problem

$\min_{y\in\mathbb{R}^{\ell_{x}}}\ F(y)\quad\text{such
that}\quad\begin{cases}A_{k\ell}\,y=0\\\ y_{0}=1\\\ M(\sigma_{j}y)\succeq
0,&j=0,1,\dots,J.\end{cases}$ (4.1)

Here $A_{k\ell}$ is the unique $k_{x}\times\ell_{x}$ matrix whose first
$\binom{n+dk}{dk}$ columns are the columns of the matrix $A_{k}$ in 2.17, and
whose remaining columns are zero. This new SDP has the same number of
variables as the data-driven SDP 3.8 and, for every $\ell\geq dk$, it enjoys
the same convergence properties as 2.17 as $k\to\infty$. Precisely, the
following extension of theorem 2.1 holds.

######  Theorem 4.1.

Suppose assumptions 2.1 and 2.2 hold and assume $\ell\geq dk$. Let
$y^{k}\in\mathbb{R}^{\ell_{x}}$ be the optimal solution of 4.1 for fixed
$k\in\mathbb{N}$, and let $\mu\in\Pr(X)$ be optimal for 2.11. There exists an
increasing sequence $\\{k_{j}\\}_{j\geq 1}\subset\mathbb{N}$ such that
$y^{k_{j}}_{\alpha}\to\int_{X}x^{\alpha}\,\mathrm{d}\mu$ as $j\to\infty$ for
all $\alpha\in\mathbb{N}_{0}^{n}$.

The proof of this statement is the same as that of [24, Theorem 2], so we omit
the details. Briefly, assumption 2.2 ensures that the sequence
$\\{y_{\alpha}^{k}\\}_{k\geq 1}$ is uniformly bounded in $k$ for every
$\alpha\in\mathbb{N}_{0}^{n}$, so by a diagonalization argument one can
extract a subsequence $\\{k_{j}\\}_{j\geq 1}$ and real numbers
$y_{\alpha}^{*}$ such that $y_{\alpha}^{k_{j}}\to y_{\alpha}^{*}$ for every
$\alpha\in\mathbb{N}_{0}^{n}$. Moreover, always by assumption 2.2, there
exists a probability measure $\mu\in\Pr(X)$ such that
$y_{\alpha}^{*}=\int_{X}x^{\alpha}\,\mathrm{d}\mu$ for every $\alpha$. The
equality constraints in 4.1 ensure that $\mu$ is invariant for the map $f$.

The rest of this section is dedicated to proving that the data-driven SDP 3.8
reduces almost surely to problem 4.1 as the number $m$ of data points
increases. (The meaning of the quantifier ”almost surely” is made precise in
the next subsection.) We start in section 4.1 by stating key assumptions and
background results that enable our convergence analysis. In section 4.2 we
establish that the data-driven SDP 3.8 is feasible almost surely if the number
$m$ of data points is large enough. We show that its optimal solutions
converge almost surely to those of 4.1 as $m\to\infty$ in section 4.3, and
conclude in section 4.4 with further comments on our analysis.

###  4.1 Preliminaries

In order to study how solutions of the data-driven SDP 3.8 behave with
increasing amount of data, we need to introduce some assumptions on the data
collection process. The first one imposes some technical restrictions on the
polynomial map $f:X\to X$ from which the data is collected.

######  Assumption 4.1.

The data points $\\{(x_{i},z_{i})\\}_{i=1}^{m}$ satisfy $z_{i}=f(x_{i})$ where
$f:X\to X$ is a polynomial of degree $d$ such that:

  1. 1.

For every integer $p$ and every choice of multi-indices
$\alpha_{1},\ldots,\alpha_{p}\in\mathbb{N}_{0}^{n}$, the polynomials
$[f(x)]^{\alpha_{1}}-x^{\alpha_{1}},\,\ldots,\,[f(x)]^{\alpha_{p}}-x^{\alpha_{p}}$
and $1$ are linearly independent.

  2. 2.

There exists an invariant probability measure $\mu\in\Pr(X)$ for $f$ whose
support is not an algebraic set, meaning that no nonzero polynomial vanishes
on it.

Next, we assume that the data points $x_{i}$ are generated by a random
sampling process. The sequence of datasets $\\{(x_{i},z_{i})\\}_{i=1}^{m}$ is
then also a random variable, and we say that a statement holds _almost surely_
as $m\to\infty$ if it holds for almost every sequence of increasingly large
datasets.

######  Assumption 4.2.

Each data point $x_{i}$ is an independent realization of a random variable
whose distribution is a probability measure $\nu\in\Pr(X)$.

The assumption of independence can in fact be relaxed: the points $x_{i}$ only
need to be sampled from $\nu\in\Pr(X)$ in such a way that, as $m\to\infty$,
the empirical averages of polynomial functions over the dataset converge
almost surely to the averages calculated with respect to the sampling measure
$\nu$ [4, Assumption 4.2]. This is true, for instance, when the datapoints
$x_{i}$ are collected from an ergodic trajectory, as we do for the
computational examples of Section 6.

Irrespective of how exactly the data is sampled, we assume that no nonzero
polynomial vanishes on the support of the sampling measure $\nu$.

######  Assumption 4.3.

If $g$ is a polynomial such that $\int_{X}g(x)\,\mathrm{d}\nu=0$, then
$g(x)\equiv 0$.

This assumption ensures that two polynomials are equal if they take the same
values on the support of $\nu$. In particular, if a sequence of polynomials
$\\{g_{j}\\}$ converges to a polynomial $g$ as $j\to\infty$ in the Lebesgue
space $L^{2}(\nu)$ of functions that are square integrable with respect to
$\nu$, then they converge pointwise on $\mathbb{R}^{n}$ and uniformly on every
compact set.

Finally, we recall the following version of a more general convergence result
for approximate Lie derivatives from [4], where $\mathcal{P}_{\ell}(u)$
denotes the $L^{2}(\nu)$-orthogonal projection of a function $u$ onto the
space of degree-$\ell$ polynomials. This is a key ingredient in our
convergence analysis.

######  Theorem 4.2 ([4, Theorem 4.1]).

For every degree-$k$ polynomial $g$ and every $\ell\geq k$, the approximate
Lie derivative $\mathcal{L}_{k\ell}g$ in 3.6 converges to
$\mathcal{P}_{\ell}(\mathcal{K}g-g)$ almost surely as $m\to\infty$.

We stress that this theorem only guarantees convergence as the number of
datapoints $m$ grows without bound, but provides no convergence rates.
Analysis [37] and computational evidence [4] reveal that if the data is
generated by random sampling, then approximation errors decay like
$1/\sqrt{m}$. However, this rate could be even slower if data is collected
from a single trajectory because time averages can be made to converge to
their ergodic value arbitrarily slowly [26].

###  4.2 Almost-sure Strict Feasibility

Having stated all of our assumptions, we now show that if the number $m$ of
data points and the polynomial degree $\ell$ are large enough, then the data-
driven SDP 3.8 is almost surely _strictly feasible_. This means that there
exists a vector $y$ satisfying the constraints of 3.8 for which the moment
matrices $M(\sigma_{j}y)$ are positive definite.

######  Proposition 4.1.

Suppose assumptions 4.1, LABEL:, 4.2, LABEL:, and 4.3 hold. For every $k$ and
every $\ell\geq dk$, there almost surely exists $m_{0}\in\mathbb{N}$ such that
the SDP 3.8 is strictly feasible for all $m\geq m_{0}$.

###### Proof.

Let $A_{k\ell}$ be the $k_{x}\times\ell_{x}$ matrix appearing in problem 4.1.
Consider the vector $\bar{y}\in\mathbb{R}^{\ell_{x}}$ of degree-$\ell$ moments
of the invariant measure $\mu$ in the second part of assumption 4.1. This
assumption ensures that $A_{k\ell}\bar{y}=0$, $\bar{y}_{0}=1$, and
$M(\sigma_{j}\bar{y})\succ 0$, so $\bar{y}$ is a strictly feasible point of
4.1.

Next, let $e\in\mathbb{R}^{\ell_{x}}$ be the vector such that $1=e\cdot
p_{\ell}(x)$ and set

$B=\begin{bmatrix}A_{k\ell}\\\
e^{\top}\end{bmatrix}\in\mathbb{R}^{(k_{x}+1)\times\ell_{x}}\qquad
B_{m}=\begin{bmatrix}L_{k\ell}\\\
e^{\top}\end{bmatrix}\in\mathbb{R}^{(k_{x}+1)\times\ell_{x}}\qquad
c=\begin{bmatrix}0\\\ 1\end{bmatrix}\in\mathbb{R}^{k_{x}+1},$ (4.2)

so the equality constraints of the SDPs 4.1 and 3.8 can be written as
$B_{m}y=c$ and $By=c$, respectively. Note that $B_{m}$ depends on the number
of data points because so does $L_{k\ell}$, while $B$ does not. For each $m$,
the vector

$y^{m}:=(I-B_{m}^{\dagger}B_{m})\bar{y}+B_{m}^{\dagger}c$ (4.3)

satisfies $B_{m}y^{m}=c$. We need to prove that, almost surely, there exists
an integer $m_{0}>0$ such that $M(\sigma_{j}y^{m})\succ 0$ when $m\geq m_{0}$,
meaning that $y^{m}$ is a strictly feasible point for 3.8.

Since the matrices $M(\sigma_{j}\bar{y})$ are positive definite and their
eigenvalues are continuous with respect to perturbations in $\bar{y}$, it is
enough to show that $y^{m}$ converges to $\bar{y}$ almost surely as
$m\to\infty$. To see this, we use the definition of $y^{m}$ and the identity
$c=B\bar{y}$ to estimate

$\displaystyle\|y^{m}-\bar{y}\|$
$\displaystyle=\|B_{m}^{\dagger}B_{m}\bar{y}-B_{m}^{\dagger}c\|$ (4.4)
$\displaystyle=\|B_{m}^{\dagger}\left(B_{m}\bar{y}-B\bar{y}\right)\|$
$\displaystyle\leq\|B_{m}^{\dagger}\|\|B_{m}-B\|\|\bar{y}\|.$

We now claim that $B_{m}$ converges almost surely to $B$ and that the norm
$\|B_{m}^{\dagger}\|$ is almost surely uniformly bounded in $m$, so
$\|y^{m}-\bar{y}\|$ tends to zero almost surely as desired.

For the first claim, fix an arbitrary degree-$k$ polynomial $g(x)=w\cdot
p_{k}(x)$, where $w\in\mathbb{R}^{k_{x}}$. By theorem 4.2, the approximate Lie
derivative $\mathcal{L}_{k\ell}g=w\cdot L_{k\ell}p_{\ell}(x)$ converges almost
surely to $\mathcal{P}_{\ell}(\mathcal{K}g-g)=\mathcal{P}_{\ell}(g\circ f-g)$
as $m$ is raised. Since $g\circ f-g=w\cdot A_{k\ell}p_{\ell}$ is a polynomial
of degree at most $dk\leq\ell$, we must have $\mathcal{P}_{\ell}(g\circ
f-g)=w\cdot A_{k\ell}p_{\ell}$ on the support of the data sampling measure
$\nu$. Then, assumption 4.3 guarantees that $w\cdot L_{k\ell}p_{\ell}(x)$
converges to $w\cdot A_{k\ell}p_{\ell}(x)$ pointwise in $x$. This is true for
all vectors $w\in\mathbb{R}^{k_{x}}$, so $L_{k\ell}$ must converge almost
surely to $A_{k\ell}$ as $m\to\infty$. This, in turn, proves our claim that
$B_{m}\to B$ almost surely as $m\to\infty$.

To see that $\|B_{m}^{\dagger}\|$ is almost surely uniformly bounded in $m$,
instead, note that the first condition in assumption 4.1 implies that the
matrix $B$ has full row rank. This, combined with the almost sure convergence
$B_{m}\to B$ and the sequential lower semicontinuity of the rank function,
ensures that $B_{m}$ has a full row rank almost surely for large enough $m$.
Since pseudoinversion is continuous on constant rank sequences [34], we obtain
that $B_{m}^{\dagger}\to B^{\dagger}$ almost surely as $m\to\infty$, so
$\|B_{m}^{\dagger}\|$ is almost surely uniformly bounded. ∎

###  4.3 Convergence with Infinite Data

We now establish our main convergence result: for every sufficiently large
$\ell$, the sequence of optimal solutions of the SDP 3.8 obtained for
increasing $m$ has almost surely a subsequence that converges to an optimal
solution of the model-based SDP 4.1. This result can be bootstrapped to the
convergence results from theorem 4.1 to conclude that if one takes $\ell\geq
dk$ and lets the parameters $m$ and $k$ tend to infinity, then our data-driven
approach computes moments of invariant measures for the polynomial map $f$
generating the data. Observe that the parameters $m$ and $k$ must be increased
in this precise order because one must obtain 4.1 as the infinite-data limit
of 3.8 before one can let $k\to\infty$ and apply theorem 4.1.

######  Theorem 4.3.

Suppose assumptions 2.2, LABEL:, 4.1, LABEL:, 4.2, LABEL:, and 4.3 hold. Let
the cost function $F$ in 3.8 and 4.1 be continuous. For every $k$ and every
$\ell\geq dk$, the sequence $\\{y^{m}\\}_{m\geq 1}$ of optimal solutions of
the SDP 3.8 almost surely has a subsequence that converges to an optimal
solution of the SDP 4.1 as $m\to\infty$.

###### Proof.

Let us first prove that SDP 4.1 has an optimal solution. The second condition
in assumption 4.1 guarantees the existence of a strictly feasible point $y$,
while assumption 2.2 ensures that the feasible set is compact. Every
minimizing sequence must then have a convergent subsequence, whose limit is an
optimizer by the continuity of the cost function $F$. The same arguments show
that the SDP 3.8 has an optimal solution $y^{m}$ when it is feasible, which is
true almost surely for large enough $m$ by proposition 4.1.

Consider now the sequence $\\{y^{m}\\}_{m\geq 1}\subset\mathbb{R}^{\ell_{x}}$
of optimal solutions of 3.8. Since each vector $y^{m}$ belongs to the set
$\\{y\in\mathbb{R}^{\ell_{x}}:\;M(\sigma_{j}y)\succeq 0\quad\forall
i=1,\ldots,J\\}$, which is compact by assumption 2.2, we can almost surely
extract a subsequence of $\\{y^{m}\\}_{m\geq 1}$ (not relabelled for
simplicity) that converges to a point $y^{*}$. We claim that $y^{*}$ is almost
surely an optimal solution of 4.1.

To prove this claim, we observe that $y^{*}$ satisfies the matrix inequalities
$M(\sigma_{j}y)\succeq 0$ for all $i=1,\ldots,J$ by construction. Moreover,
with $B$, $B_{m}$ and $c$ defined as in the proof of proposition 4.1, we have

$\displaystyle\|By^{*}-c\|$ $\displaystyle=\|By^{*}-B_{m}y^{m}\|$ (4.5)
$\displaystyle\leq\|By^{*}-By^{m}\|+\|(B-B_{m})y^{m}\|$
$\displaystyle\leq\|B\|\|y^{*}-y^{m}\|+\|B-B_{m}\|\|y^{m}\|$

because the points $y^{m}$ satisfy $B_{m}y^{m}=c$. The expression in the last
line tends to zero almost surely as $m\to\infty$ because $y^{m}\to y^{*}$ and
$B_{m}\to B$ almost surely. Thus, $By^{*}=c$ almost surely and we conclude
that $y^{*}$ is almost surely feasible for 4.1.

Next, let $y^{\text{opt}}$ and $y^{\text{f}}$ be an optimal and a strictly
feasible point for 4.1, respectively. Their existence is guaranteed, as
explained at the start of the proof. Given arbitrary $\epsilon>0$, let
$\delta\in(0,1)$ be small enough that the point

$y^{\epsilon}:=(1-\delta)y^{\text{opt}}+\delta y^{\text{f}}$ (4.6)

satisfies

$F(y^{\epsilon})\leq F(y^{\text{opt}})+\epsilon.$ (4.7)

It is not difficult to verify that $y^{\epsilon}$ is strictly feasible for
4.1. Moreover, arguing as in the proof of proposition 4.1 shows that

$y^{\epsilon,m}:=(I-B_{m}^{\dagger}B_{m})y^{\epsilon}+B_{m}^{\dagger}c$ (4.8)

converges almost surely to $y^{\epsilon}$ as $m\to\infty$ and is therefore
almost surely feasible for 3.8 for all large enough $m$. The continuity of
$F$, the optimality of $y^{\text{opt}}$ and $y^{m}$, and inequality 4.7 then
imply

$F(y^{\text{opt}})\leq
F(y^{*})=\lim_{m\to\infty}F(y^{m})\leq\lim_{m\to\infty}F(y^{\epsilon,m})=F(y^{\epsilon})\leq
F(y^{\text{opt}})+\epsilon.$ (4.9)

Since $\epsilon>0$ is arbitrary we conclude that $F(y^{\text{opt}})=F(y^{*})$,
so $y^{*}$ is optimal for 4.1. ∎

###  4.4 Further Discussion

We conclude this section by showing that, in general, the convergence of
approximate Lie derivative as guaranteed by theorem 4.2 is by itself not
sufficient to ensure that the data-driven SDP 3.8 converges to its exact
counterpart 4.1 as the number $m$ of data points grows.

Consider the map $f(x)=2x-x^{2}$ on the interval $X=[0,1]$, which we represent
using polynomial inequalities as $X=\\{x\in\mathbb{R}:\;x-x^{2}\geq 0\\}$. We
seek the invariant measure $\mu$ with the largest possible first moment,
$y_{1}=\int_{0}^{1}x\,\mathrm{d}\mu$. The solution is easily seen to be the
atomic measure supported at the fixed point $x=1$, for which $y_{1}=1$.

We now set up our SDP approach for polynomial degrees $k=1$ and $\ell=2$,
which satisfy $\ell\geq dk$ because the polynomial map $f$ we are considering
has degree $d=2$. Upon enforcing the constraint $y_{0}=1$ explicitly, for
these parameter choices the model-based SDP 4.1 reads

$\inf_{y\in\mathbb{R}^{2}}\ -y_{1}\quad\text{such
that}\quad\begin{cases}y_{1}-y_{2}=0,\\\\[4.30554pt] \begin{pmatrix}1&y_{1}\\\
y_{1}&y_{2}\end{pmatrix}\succeq 0,\\\\[4.30554pt] y_{1}-y_{2}\geq
0.\end{cases}$ (4.10)

The first constraint in this problem comes from enforcing the invariance
condition
$\int_{X}\mathcal{L}x\,\mathrm{d}\mu=\int_{X}(x-x^{2})\,\mathrm{d}\mu=0$,
which is the only one we must consider given our choice $k=1$. The other two
constraints are the moment matrices obtained for $\ell=2$ from the general
definition 2.16 with the polynomials $\sigma_{0}(x)=1$ and
$\sigma_{1}(x)=x-x^{2}$. The optimal solution of 4.10 is easily seen to be
$y=(1,1)$ (see the left panel of figure 1).

Figure 1: Feasible sets for the SDPs 4.10 (left) and 4.11 (right). In both
panels, the orange shaded region is the set of points satisfying the (matrix)
inequality constraints. Straight black lines indicate the linear spaces
$y_{1}-y_{2}=0$ (left panel) and $\alpha_{m}y_{1}-\beta_{m}y_{2}=0$ for
$\alpha_{m}/\beta_{m}>1$ (right panel). Optimal points in both panels are
marked by a blue dot.

The data-driven SDP 3.8, instead, reads

$\inf_{y\in\mathbb{R}^{2}}\ -y_{1}\quad\text{such
that}\quad\begin{cases}\alpha_{m}y_{1}-\beta_{m}y_{2}=0,\\\\[4.30554pt]
\begin{pmatrix}1&y_{1}\\\ y_{1}&y_{2}\end{pmatrix}\succeq 0,\\\\[4.30554pt]
y_{1}-y_{2}\geq 0,\end{cases}$ (4.11)

for some coefficients $\alpha_{m}$ and $\beta_{m}$ that depend on the size $m$
of the dataset used to approximate the Lie derivative $\mathcal{L}x=x-x^{2}$
with EDMD. While theorem 4.2 guarantees that $\alpha_{m},\,\beta_{m}\to 1$ as
$m\to\infty$, it gives no additional control on approximation errors, so _a
priori_ one cannot rule out the possibility that $\alpha_{m}/\beta_{m}>1$ for
all $m$. In this case, as illustrated in the right panel of figure 1, for
every $m$ the optimal solution of 4.11 would be $y^{m}=(0,0)$. Thus, $y^{m}$
would not converge to the optimal solution $y=(1,1)$ of 4.10 despite the
convergence of $\alpha_{m}$ and $\beta_{m}$.

This potential failure of convergence stems from the lack of strict
feasibility for the SDP 4.10. As shown in the left panel of figure 1, all
feasible points lie on the boundary of the set defined by the (matrix)
inequality constraints. Optimal solutions are therefore not robust to
arbitrary perturbations to the equality constraints. Such a scenario is not
possible in section 4.3, as the SDP 4.1 is strictly feasible by assumption
4.1(ii). This assumption fails for the present example because the only
invariant measures for the map $f(x)=2x-x^{2}$ on $X=[0,1]$ are the atomic
measures supported at the fixed points $x=0$ and $x=1$, which are an algebraic
set.

We must however point out that the possible failure of convergence in our
present example can be ruled out by an argument different from the proof in
section 4.3. Indeed, since $f$ is a quadratic polynomial map, it can be
recovered exactly by sampling it at $m\geq 3$ distinct points
$x_{1},\ldots,x_{m}\in[0,1]$. With such a dataset, EDMD recovers the _exact_
Lie derivative $\mathcal{L}x$, giving $\alpha_{m}=\beta_{m}=1$ in 4.11. Since
the datapoints $x_{1},\ldots,x_{m}$ are almost surely distinct when sampled
according to assumptions 4.2 and 4.3, the SDP 4.11 coincides with 4.10 almost
surely for every $m\geq 3$. The reasoning can in fact be extended to any
$n$-dimensional polynomial map of degree $d$ and any choice of the parameters
$k$ and $\ell\geq dk$: problem 3.8 coincides with 4.1 if (i) a subset of the
data points $x_{i}$ are _unisolvent_ for the space of degree-$\ell$
polynomials, and (ii) the map $f(x)$, or at least the Lie derivatives of
degree-$k$ polynomials required by the EDMD step of our approach, are sampled
exactly. This finer result, however, is specific to deterministic polynomial
maps. The proofs in sections 4.2 and 4.3, instead, extend to generalizations
of our SDP approach to continuous-time and/or stochastic dynamical systems,
discussed in the next section. The only requirements for this extension are
that the dynamics be governed by polynomial equations, and that approximate
Lie derivatives built with EDMD converge to their exact counterparts.

##  5 Extensions to Other Dynamical Systems

Our data-driven approach to discovering invariant measures can be applied to
dynamical systems beyond those governed by discrete maps. The key observation
is that the invariance condition in section 2 involves the _Lie derivative_

$\mathcal{L}g:=g\circ f-g$ (5.1)

of the test function $g$, which describes how $g$ changes along trajectories
of system 2.1. As such, Lie derivatives can be defined for a much broader
class of dynamical systems (see [4, Section 3]). Classical examples include:

  * •

Dynamical systems governed by well-posed autonomous ordinary differential
equations. Precisely, let $x(t)$ solve the differential equation
$\frac{d}{dt}x=f(x)$. Then, the Lie derivative of any continuously
differentiable function $g$ is given by

$\mathcal{L}g=f\cdot\nabla g.$ (5.2)

  * •

Discrete-time systems governed by stochastic maps. Precisely, let
$x_{n+1}=f(\omega_{n},x_{n})$ be a stochastic process where the map
$\omega\mapsto f(\omega,\cdot)$ is a random variable from a probability space
$(\Omega,\mathcal{F},\pi)$ into the space of maps from $X$ to itself. For such
a stochastic process, the Lie derivative acts on continuous functions $g$ by

$\mathcal{L}g(x)=\int_{\Omega}g(f(\omega,x))\mathrm{d}\pi(\omega)-g(x).$ (5.3)

  * •

Stochastic processes governed by stochastic differential equations. Precisely,
let $X_{t}$ be the solution of
$\mathrm{d}X_{t}=a(X_{t})\mathrm{d}t+b(X_{t})\mathrm{d}W_{t}$, where $W_{t}$
is a standard $k$-dimensional Wiener process and
$a:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and
$b:\mathbb{R}^{n}\to\mathbb{R}^{n\times k}$ are given Lipschitz functions.
Then, the Lie derivative of any twice-differentiable function $g$ is

$\mathcal{L}g=a\cdot\nabla g+\frac{1}{2}\langle
bb^{\intercal},\nabla^{2}g\rangle,$ (5.4)

where $\nabla^{2}g$ is the Hessian of $g$ and the angled brackets indicate the
Frobenius (element-wise) inner product.

In all of these cases, invariant measures are characterized by the condition
that

$\int_{X}\mathcal{L}g\,\mathrm{d}\mu=0$ (5.5)

for all functions $g$ for which the Lie derivative is defined. In addition,
EDMD with polynomial dictionaries can be used to approximate the Lie
derivative from data with guaranteed convergence as the amount of data and the
size of the polynomial dictionaries increase [4]. One can then construct and
solve the SDP 3.8, obtaining approximations of invariant measures even for
continuous-time or stochastic systems. The only difference compared to
discrete-time deterministic maps is that the data snapshots
$\\{(x_{i},z_{i})\\}_{i=1}^{m}$ must be such that, for a fixed timestep
$\tau>0$, $z_{i}$ is the state of the system at time $t+\tau$ when the state
at time $t$ is $x_{i}$. The continuous-time Lie derivative is then
approximated using the matrix

$L_{k\ell}=\frac{K_{k\ell}-\Theta_{k\ell}}{\tau},$ (5.6)

where $K=QP^{\dagger}$, the matrices $P$ and $Q$ are as in 3.1, and
$\Theta_{k\ell}$ is the rectangular matrix satisfying
$p_{k}(x)=\Theta_{k\ell}\,p_{\ell}(x)$. In the discrete-time stochastic
setting one simply fixes $\tau=1$.

##  6 Examples

We now illustrate the data-driven discovery of invariant measures on a series
of examples of increasing difficulty. The first two examples demonstrate our
approach on the logistic map (section 6.1) and on a stochastic double-well
system (section 6.2). Invariant measures for these examples can be easily
found by other methods and compared to our data-driven approximations. In
section 6.3, instead, we apply our method to approximate the physical measure
of the three-dimensional Rössler attractor and to extract unstable periodic
orbits embedded in it. We do this by identifying discrete-time periodic orbits
for a suitably defined Poincaré return map. This last example is purely data-
driven, in the sense that no analytical expression for the Poincaré return map
is available. For all examples, we implemented the SDP 3.8 using the MATLAB
toolbox YALMIP [31] and solved it using MOSEK [33]. Code to reproduce our
results is available at https://github.com/jbramburger/data-measures.

###  6.1 The Logistic Map

We begin by discovering invariant measures for the scalar logistic map,

$x_{n+1}=2x_{n}-1,$ (6.1)

on the interval $X=[-1,1]$. Our goal is to illustrate the performance of our
data-driven approach and highlight its ability to identify the system’s
physical measure as well as various ergodic measures that have an atomic
structure.

####  6.1.1 Physical Measures

The logistic map 6.1 is one of few nonlinear discrete dynamical systems for
which an exact invariant measure is known [19]. The density of this invariant
measure is

$\rho_{*}(x)=\frac{1}{\pi\sqrt{1-x^{2}}}.$ (6.2)

Here, we seek to ‘discover’ this density from data collected by simulating $m$
iterations of 6.1 starting from $x_{0}=0.25$. We averaged the ensuing
trajectory $\\{x_{n}\\}_{n=1}^{m}$ to obtain an approximate value
$\tilde{y}_{1}$ for the first moment $y_{1}$ of the physical invariant
measure, and we then implemented the SDP 3.8 with the objective function

$F(y)=(y_{1}-\tilde{y}_{1})^{2}.$ (6.3)

Observe that the value of $\tilde{y}_{1}$ changes as the amount of data is
increased and converges to the exact value $y_{1}=0$ as $m\to\infty$. For
numerical stability, we used ChebFun [12] to represent polynomials in the
Chebyshev basis $p_{k}(x)=\\{T_{0}(x),T_{1}(x),\dots,T_{k}(x)\\}$, rather than
the monomial basis. We also took $\ell=2k$, so
$p_{\ell}(x)=\\{T_{0}(x),T_{1}(x),\dots,T_{2k}(x)\\}$, but nearly identical
results are obtained with $\ell>2k$ because this brings no notable improvement
in the data-driven approximation of the Lie derivative of degree-$k$
polynomials [4, Section 5.3].

Optimal solutions to 3.8 were used to construct signed polynomial densities
$\rho_{k}$ for the underlying approximate invariant measure as explained in
section 2.2.4. We fixed the reference measure $\pi$ to be the Lebesgue measure
on $[-1,1]$ and the results for $k=5$, $10$ and $20$ are compared to the true
invariant density $\rho_{*}$ in figure 2. It is clear that the signed
densities $\rho_{k}$ approximate $\rho_{*}$ well, but they also exhibit a
highly oscillatory behaviour. Such oscillations are consistent with the fact
that, as both the number $m$ of data snapshots and the polynomial degree $k$
become infinite, $\rho_{k}$ is expected to converge to $\rho_{*}$ only weakly
in $L^{1}$ (cf. section 2.2.4), which does not imply pointwise convergence.
Nevertheless, the moments of the approximate densities converge to exact
values, so our approximations can be used to make meaningful statistical
predictions.

Figure 2: Discovered densities (blue, solid) for the logistic map 6.1 using
$m=10^{4}$ data points compared to the exact density 6.2 (red, dashed).

To corroborate this claim, we examined the $L^{1}$ norm difference between the
approximate and exact cumulative distribution functions,

$R_{k}(x)=\int_{-1}^{x}\rho_{k}(s)\mathrm{d}s\qquad\text{and}\qquad
R_{*}(x)=\int_{-1}^{x}\rho_{*}(s)\mathrm{d}{s}=\frac{1}{2}+\frac{1}{\pi}\arcsin(x).$
(6.4)

Indeed, the weak convergence of $\rho_{k}$ to $\rho_{*}$ is equivalent to the
pointwise convergence of $R_{k}$ to $R_{*}$ at every point of continuity in
$R_{*}$ [20, Theorem 13.23]. Since $R_{*}$ is continuous for all $x\in[-1,1]$
and $0\leq R_{k}(x)\leq 1$ for all $x$ and $k$, the dominated convergence
theorem guarantees that $R_{k}$ converges to $R_{*}$ in the $L^{1}$ norm. The
$L^{1}$ norm difference between $R_{k}$ and $R_{*}$ is reported in table 1 for
increasing polynomial degrees $k$ and increasing size $m$ of the training
data. The table also reports the $L^{1}$ norm differences obtained by solving
the SDP 2.17, which uses the exact Lie derivative, and the exact first moment
value $\tilde{y}_{1}=0$. For $m\geq 10^{4}$ the data-driven results closely
match those with exact knowledge of the system. The main limitation seems to
be approximating the exact invariant measure using semidefinite programming
with finite polynomial degree $k$. This is likely due to the singular
behaviour of the optimal density $\rho_{*}$ at the boundaries of the interval
$[-1,1]$, which cannot be approximated well using polynomials of low degree.

Table 1: $L^{1}$ norm difference between the cumulative distributions 6.4 of
the discovered and exact densities. Numerical procedures use $m$ consecutive
data points under the dynamics of the Logistic map seeded by the initial
condition $x_{0}=0.25$ with an objective to fit the approximate first moment
$\tilde{y}_{1}$ extracted from the data. The values $\tilde{y}_{1}$ are
reported in the final column.

Exact values are computed using the exact Lie derivative and first moment
$\tilde{y}_{1}=0$ in place of those that are estimated from data.  $k=5$
$k=10$ $k=15$ $k=20$ $k=25$ $\tilde{y}_{1}$ $m=10^{2}$ 0.04488 0.04401 0.04376
0.04336 0.04325 -0.04123 $m=10^{3}$ 0.02717 0.01602 0.01474 0.01440 0.01389
-0.01339 $m=10^{4}$ 0.02872 0.00955 0.00531 0.00372 0.00285 -0.00102
$m=10^{5}$ 0.02891 0.00959 0.00545 0.00371 0.00278 -0.00035 Exact 0.02902
0.00962 0.00556 0.00378 0.00280 0.0000

There are of course other ways to approximate the density of a physical
measures from data, such as building a histogram. For demonstration, we take
our density used in Table 1 with $m=10^{3}$ and $k=20$ and compare it to the
density approximated by the histogram constructed from the same data and
organized into 101 bins. The resulting histogram is plotted in Figure 3, where
it is compared to the exact density (6.2) and to the approximate density
obtained with our SDP method and the same data. To evaluate the accuracy of
the two approximate densities, we compare their moments against the exact
values. Relative errors of these predictions are presented in Table 2 for the
first five nonzero moments ($y_{2}$, $y_{4}$, $y_{6}$, $y_{8}$, and $y_{10}$),
and three higher-order moments ($y_{20}$, $y_{30}$, and $y_{40}$). Notice that
while the relative errors are comparable for the lower-order moments, the
difference in the methods is clearly observed in the approximation of the
higher-order moments. In particular, the histogram method provides a poor
approximation of the higher-order moments, nearly doubling the relative error
of the method of this manuscript. Of course, better accuracy can be achieved
by using more data, but our investigations reveal that our SDP-based approach
can achieve greater accuracy, particularly for the higher-order moments.

Figure 3: Left: Histogram approximation (blue) of the density of the physical measure (6.2) of the logistic map against the exact density (red, dashed). We use $m=10^{3}$ datapoints along a single trajectory starting with $x_{0}=0.25$ and organize data into 101 bins. Right: The approximated density using the SDP methods of this manuscript (blue) against the true density (red, dashed) using the same datapoints as the left panel and taking $k=20$.  Table 2: Relative errors in approximating the moments of (6.2) via the approximated densities from the semidefinite programming (SDP) method of Section 3 and from building a histogram (Histogram) from the data. Both methods use the same $m=10^{3}$ datapoints, while the histogram organizes the data into 101 bins covering the domain $[-1,1]$.  Moment | $y_{2}$ | $y_{4}$ | $y_{6}$ | $y_{8}$ | $y_{10}$ | $y_{20}$ | $y_{30}$ | $y_{40}$
---|---|---|---|---|---|---|---|---
SDP | 1.33% | 2.23% | 2.81% | 3.25% | 2.63% | 5.04% | 6.11% | 7.02%
Histogram | 1.51% | 2.62% | 3.49% | 4.26% | 4.99% | 8.35% | 11.78% | 15.46%

####  6.1.2 Ergodic Measures

The logistic map 6.1 has a chaotic attractor that is densely filled by an
infinite number of unstable periodic orbits, each of which is associated with
an ergodic measure with an atomic structure. As discussed in section 2.2.2,
moments of such ergodic measures can be approximated by solving the exact SDP
2.17 or its data-driven approximation 3.8 with a linear objective function.
Moreover, one can approximate these atomic measures using absolutely
continuous ones as explained in section 2.2.4. In contrast, a data-based
histogram cannot provide approximations of these atomic measures.

To illustrate this, we solved the SDP 3.8 for polynomial degrees $k=5$,$10$
and $20$ and the linear cost function $F(y)=y_{1}$. The data snapshots are as
for the discovery of physical measures above. Figure 4 shows the densities of
the Lebesgue-absolutely-continuous signed measures recovered from the SDP
solution. These densities become increasingly localized near $x=-0.5$, which
is a fixed point for the logistic map. This suggests that our approximations
are converging (in a weak-star sense) to the invariant atomic measure
supported on this fixed point. In fact, if the optimal solution $y$ of the SDP
3.8 is such that the moment matrix $M(y)$ satisfies a rank deficiency
condition [30, Theorem 5.20], then one can extract an atomic measure from $y$
[30, Section 6.7]. We refer interested readers to this work and the references
therein for the technical details. For our logistic map example, this
extraction procedure111Implemented in the function extractminimizers from the
modified version of YALMIP available at: https://github.com/aeroimperial-
optimization/aeroimperial-yalmip. returns at atomic measure supported at
$x\approx-0.5$ even for $k=5$, where the approximate density in figure 4 is
far from being localized.

Figure 4: Polynomial densities of Lebesgue-absolutely-continuous measures that
approximates the atomic invariant measure supported on the fixed point
$x=-0.5$ of the logistic map 6.1. Results are for $m=10^{4}$ data snapshots
and polynomial degrees $k=5$, $10$ and $20$.

In a similar way, by varying the linear optimization objective we are able to
discover periodic points of the logistic map. For example, fixing
$(m,k)=(10^{4},5)$ and taking the optimization objective to be $F(y)=y_{3}$ we
are able to extract from the optimal solution of the SDP 3.8 an atomic measure
supported on the points $x\approx-0.93969$, $0.17365$ and $0.76604$, which
approximate a period-3 orbit of the logistic map to four decimal places.

######  Remark 6.1.

While linear objective functions are always optimized by at least one ergodic
measure, our approach is not guaranteed to find it. Indeed, suppose two
ergodic measures optimize the same objective. Then, solutions of the SDPs 2.17
and 3.8 may converge to moments of any convex combination of these ergodic
measures, which remains optimal but is not ergodic. For instance, with
$F(y)=y_{5}$ in the above example we discover an atomic measure supported on
the fixed point $x=-0.5$ and on a period-4 orbit of the logistic map, which is
a convex combination of the ergodic measures corresponding to each of these
two invariant trajectories. However, in the absence of symmetry or other
specific structure in the data, we expect such examples not to be generic.

###  6.2 A Stochastic Double-Well System

For our next example we consider the stochastic differential equation

$\displaystyle\dot{x}_{1}$
$\displaystyle=-16(x_{1}+x_{2})^{3}+2x_{1}+6x_{2}+\sigma\xi_{1}(t),$ (6.5)
$\displaystyle\dot{x}_{2}$
$\displaystyle=-16(x_{1}+x_{2})^{3}+6x_{1}+2x_{2}+\sigma\xi_{2}(t),$

where $\xi_{1}(t)$ and $\xi_{2}(t)$ are independent Gaussian white noise
processes and $\sigma$ determines the noise intensity. The unique invariant
measure associated to this equation is absolutely continuous with respect to
the Lebesgue measure on $\mathbb{R}^{2}$ and has the density

$\rho(x_{1},x_{2})=C\exp\left\\{-\frac{[4(x_{1}+x_{2})^{2}-1]^{2}+4(x_{1}-x_{2})^{2}}{2\sigma^{2}}\right\\},$
(6.6)

where the constant $C$ is such that
$\int_{\mathbb{R}^{2}}\rho(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}=1$. This
density can be found by solving the stationary Fokker-Plank equation
associated to 6.5. As illustrated in the left-most panel of figure 5, $\rho$
exhibits two peaks at $(x_{1},x_{2})=(\frac{1}{4},\frac{1}{4})$ and
$(x_{1},x_{2})=(-\frac{1}{4},-\frac{1}{4})$, which are the stable fixed points
of the deterministic ($\sigma=0$) counterpart of 6.5. The peaks become
increasingly localized as $\sigma$ decreases. In this example, we use our
data-driven discovery approach to approximate this two-peak invariant density
for $\sigma=0.75$. This value gives two clear peaks, which are however
sufficiently ”spread out” for $\rho$ to be approximated well by polynomials of
relatively small degree. Smaller $\sigma$ could be considered at the cost of
working with larger polynomial degrees (hence, larger computational cost) than
those reported below.

Figure 5: Left-most panel: Invariant density $\rho$ for the stochastic system
in 6.5 for with $\sigma=0.75$. The formula for $\rho$ is given by 6.6. Other
panels: Approximations to $\rho$ obtained from the dataset $\mathcal{D}_{1}$
with (in order) with our SDP approach, with a data-based histogram, and with a
two-component Gaussian mixture fit to the data. In each panel, dotted white
lines mark the reflection symmetry axes ($x_{2}=x_{1}$ and $x_{2}=-x_{1}$) for
the exact invariant density.

We consider two sets of data snapshots $(x(t),x(t+\tau))$ in
$X=\mathbb{R}^{2}$, both collected by approximating solutions of 6.5 with the
Euler–Maruyama scheme and timestep $\tau=10^{-4}$. The first dataset, referred
to as $\mathcal{D}_{500}$ in what follows, was obtained from $500$
realizations of the stochastic process, integrated for $10^{4}$ timesteps
starting from $500$ initial conditions sampled from the uniform distribution
on the square $[-\frac{1}{2},\frac{1}{2}]$. The second dataset, denoted by
$\mathcal{D}_{1}$, was instead obtained from a single realization with
$5\times 10^{6}$ timesteps, again started from a random initial condition on
the square $[-\frac{1}{2},\frac{1}{2}]$. Datapoints $(x(t),x(t+\tau))$ not in
the unit box $[-1,1]^{2}$ were dropped from both datasets to ensure a
numerically stable evaluation of Chebyshev polynomials (see the next
paragraph). As a result, the dataset $\mathcal{D}_{500}$ contains
$m=4\,999\,759$ snapshots from $500$ process realizations, while the dataset
$\mathcal{D}_{1}$ has $m=4\,998\,925$ snapshots from a single, but longer,
process realization. The datasets are included in the code repository
accompanying this article.

For each dataset, we solved the SDP 3.8 for increasing polynomial degree $k$
and polynomial EDMD dictionaries of degrees $\ell=k$, $k+2$ and $k+4$. As
before, for numerical stability we represent polynomials using the Chebyshev
basis, rather than the monomial basis. We also found that the quality of
approximate Lie derivatives could be considerably improved if we computed the
matrix $L_{k\ell}$ with an iterative procedure similar to the SINDy method
[8]. Briefly, for each row of the matrix, we find the entry with largest
magnitude, set to zero entries with magnitude smaller than 0.1% of this value,
and re-optimize the others. The process can be repeated until the matrix
$L_{k\ell}$ stabilizes. Finally, the objective function for the SDP 3.8 was
set to $F(y)=-2-y_{20}-y_{02}$, which corresponds to looking for the invariant
measure that maximizes the stationary expectation of the observable function
$x_{1}^{2}+x_{2}^{2}$. Results for different objective functions (omitted for
brevity) are similar, but are not identical even though 6.5 admits a unique
invariant measure. This is expected because the constraints of 3.8 do not
characterize invariant measures exactly.

Figure 6: Data-driven approximations of the invariant density for the stochastic double-well system 6.5, recovered from the solution of the SDP 3.8. Different panels correspond to different choices of the polynomial degrees $k$ and $\ell$, which are indicated in each case. All plots are for the dataset $\mathcal{D}_{500}$.  Figure 7: Data-driven approximations of the invariant density for the stochastic double-well system 6.5, recovered from the solution of the SDP 3.8. Different panels correspond to different choices of the polynomial degrees $k$ and $\ell$, which are indicated in each case. All plots are for the dataset $\mathcal{D}_{1}$.  Table 3: Stationary expectations of bivariate Chebyshev polynomials $T_{\alpha}(x_{1})T_{\beta}(x_{2})$ for the stochastic system (6.5), predicted using (i) the approximate densities found with the SDP method of Section 3, (ii) a histogram, and (iii) a two-component Gaussian mixture model. Results are for the dataset $\mathcal{D}_{1}$ and all values $\alpha,\beta$ such that $\alpha+\beta\in\\{2,4,6\\}$. Errors are with respect to the exact expectation values $\int_{\mathbb{R}^{2}}T_{\alpha}(x_{1})T_{\beta}(x_{2})\rho(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}$, calculated using the exact invariant density $\rho$ from 6.6 and reported in the first column.  |  | Exact |  | SDP |  | Histogram |  | Gaussian Mixture
---|---|---|---|---|---|---|---|---
|  | Value |  | Value | Error |  | Value | Error |  | Value | Error
$T_{2}(x_{1})T_{0}(x_{2})$ |  | -0.825 |  | -0.816 | 1.1% |  | -0.817 | 1.0% |  | -0.796 | 3.6%
$T_{1}(x_{1})T_{1}(x_{2})$ |  | 0.017 |  | 0.019 | 9.5% |  | 0.019 | 10.4% |  | 0.015 | 14.7%
$T_{0}(x_{1})T_{2}(x_{2})$ |  | -0.825 |  | -0.840 | 1.8% |  | -0.841 | 1.9% |  | -0.842 | 2.0%
$T_{4}(x_{1})T_{0}(x_{2})$ |  | 0.453 |  | 0.431 | 4.9% |  | 0.432 | 4.6% |  | 0.390 | 13.8%
$T_{3}(x_{1})T_{1}(x_{2})$ |  | -0.049 |  | -0.054 | 10.0% |  | -0.054 | 10.3% |  | -0.043 | 11.6%
$T_{2}(x_{1})T_{2}(x_{2})$ |  | 0.668 |  | 0.673 | 0.8% |  | 0.674 | 1.0% |  | 0.658 | 1.4%
$T_{1}(x_{1})T_{3}(x_{2})$ |  | -0.049 |  | -0.051 | 5.3% |  | -0.052 | 5.7% |  | -0.041 | 16.8%
$T_{0}(x_{1})T_{4}(x_{2})$ |  | 0.453 |  | 0.488 | 7.7% |  | 0.490 | 8.1% |  | 0.500 | 10.4%
$T_{6}(x_{1})T_{0}(x_{2})$ |  | -0.141 |  | -0.118 | 16.3% |  | -0.118 | 16.1% |  | -0.087 | 38.1%
$T_{5}(x_{1})T_{1}(x_{2})$ |  | 0.066 |  | 0.073 | 9.8% |  | 0.072 | 9.2% |  | 0.058 | 12.3%
$T_{4}(x_{1})T_{2}(x_{2})$ |  | -0.337 |  | -0.326 | 3.0% |  | -0.328 | 2.5% |  | -0.296 | 12.0%
$T_{3}(x_{1})T_{3}(x_{2})$ |  | 0.137 |  | 0.146 | 6.4% |  | 0.146 | 6.2% |  | 0.118 | 13.8%
$T_{2}(x_{1})T_{4}(x_{2})$ |  | -0.337 |  | -0.362 | 7.5% |  | -0.364 | 8.2% |  | -0.364 | 8.3%
$T_{1}(x_{1})T_{5}(x_{2})$ |  | 0.066 |  | 0.068 | 2.2% |  | 0.068 | 2.5% |  | 0.054 | 18.5%
$T_{0}(x_{1})T_{6}(x_{2})$ |  | -0.141 |  | -0.172 | 22.0% |  | -0.174 | 23.4% |  | -0.200 | 42.0%

Figure 6 shows approximate signed invariant densities for 6.5 obtained with
the dataset $\mathcal{D}_{500}$. Results for the dataset $\mathcal{D}_{1}$ are
similar, so we report only a subset of them in figure 7. In both cases, the
approximate signed invariant densities were recovered from the optimal SDP
solutions as explained in section 2.2.4, taking the reference measure $\pi$ to
be the Lebesgue measure on $[-1,1]^{2}$. Just like the true invariant density
$\rho$ in 6.6, the signed densities localize near the points
$(\frac{1}{4},\frac{1}{4})$ and $(-\frac{1}{4},-\frac{1}{4})$ as the
polynomial degree $k$ is raised. For sufficiently large values of the
parameters $k$ and $\ell$, moreover, they approximately satisfy the reflection
symmetry of the exact invariant density $\rho$ with respect to the lines
$x_{2}=x_{1}$ and $x_{2}=-x_{1}$. The lack of exact symmetry in the figures
can be attributed to a lack of symmetry in our datasets because, as
illustrated by figure 5, the same issue arises also if we approximate $\rho$
by fitting a probability histogram or a two-component Gaussian mixture model
to the data. (This can be done with the MATLAB functions histogram2 and
fitgmdist, respectively, which we used with default parameters.) Incidentally,
as evidenced by table 3, the invariant densities produced by the SDP approach
and by fitting a histogram predict the stationary expectation of bivariate
Chebyshev polynomials $T_{\alpha}(x_{1})T_{\beta}(x_{2})$ with similar
accuracy. In contrast, the two-component Gaussian mixture model has much
larger prediction errors.

To conclude this example, we point out that the results for $\ell=k+2$ and
$\ell=k+4$ in figure 6 become almost indistinguishable as $k$ increases. This
is to be expected because the right-hand side of 6.5 is a cubic polynomial, so
5.4 and the convergence results in [4] guarantee that the choice $\ell=k+2$ is
large enough to recover exactly the Lie derivative of degree-$k$ polynomials
as the number of data snapshots tends to infinity.

###  6.3 The Rössler System

Having demonstrated our data-driven approach to discovering invariant measures
on two relatively simple examples, we now showcase its usefulness in providing
insight into the statistics of a chaotic dynamical system. Precisely, we
consider the chaotic Rössler system

$\begin{split}\dot{x}_{1}&=-x_{2}-x_{3},\\\ \dot{x}_{2}&=x_{1}+0.1x_{2},\\\
\dot{x}_{3}&=0.1+x_{3}(x_{1}-18),\end{split}$ (6.7)

and show that our data-driven approach can produce good approximations of the
statistics on the attractor, as well as enable the discovery of a large number
of unstable periodic orbits. In both cases, the dataset consists of a
numerically integrated trajectory of 6.7 starting from
$(x_{1},x_{2},x_{3})=(0,-20,0)$ at $t=0$, with data collected every
$\tau=0.005$ time units up to $t=1000$. We look for an invariant measure in
the set $X=[-30,30]\times[-30,30]\times[0,60]$, which contains the attractor.
For computations, variables are scaled so that $X$ is the unit cube
$[-1,1]^{3}$, but results are plotted in the original coordinates of 6.7.

####  6.3.1 Physical Measures and Attractor Statistics

First, we approximate the physical measure on the chaotic attractor by solving
3.8 with polynomial degrees $(k,\ell)=(14,15)$ and optimization objective

$F(y)=\frac{(y_{100}-\tilde{y}_{100})^{2}}{\tilde{y}_{100}^{2}}+\frac{(y_{010}-\tilde{y}_{010})^{2}}{\tilde{y}_{010}^{2}}+\frac{(y_{001}-\tilde{y}_{001})^{2}}{\tilde{y}_{001}^{2}}.$
(6.8)

The values $\tilde{y}_{100}$, $\tilde{y}_{010}$ and $\tilde{y}_{001}$ are the
approximate moments of $x_{1}$, $x_{2}$ and $x_{3}$, respectively, which we
obtained from the training data via 2.13. We used the method of section 2.2.4
with the reference measure $\pi$ fixed to the Lebesgue measure on $X$ in order
to recover an absolutely-continuous measure with signed polynomial density of
degree $k=14$ that approximates the attractor’s physical measure.

Figure 8: Simulated data of the Rössler system 6.7 coloured according to the
value of the discovered approximate density using $(k,\ell)=(14,15)$ and the
objective 6.8. Data on the chaotic attractor in $(x_{1},x_{2},x_{3})$-space is
on the left, while timeseries plots for each variable are on the right.

The left panel in figure 8 shows the training data points on the system’s
attractor, coloured according to the value of the discovered signed density.
The large upward swings in $x_{3}$ along the attractor are identified as rare
events, which is confirmed by the three timeseries plots on the right of
figure 8. Moreover, the highest density component is in the region $x_{1}\leq
0$, through which all trajectories seem to be funnelled as they circle around
the attractor. Finally, all moments of the discovered density up to degree 2
agree with time averages integrated up to $t=10^{6}$ to within a $13\%$ error,
showing that we can make accurate statistical predictions even for moments not
included in the objective function.

####  6.3.2 Discovery of Unstable Periodic Orbits

Next, we turn to the discovery of UPOs for the Rössler system through the
analysis of a Poincaré return map. This is a truly data-driven example because
the Poincaré return map for 6.7 does not have a closed-form analytical
expression. Yet, analysing it is convenient because its UPOs have an atomic
structure and, therefore, can more easily be detected within the polynomial-
based convex optimization framework of sections 2 and 3.

Numerical evidence suggests that the Rössler system has the advantageous
property that points on the attractor with $x_{1}=0$ also have $x_{3}=0$.
Therefore, to define the Poincaré return map we choose the Poincaré section to
be the plane $x_{1}=0$ with the crossing condition $\dot{x}_{1}>0$ and only
track the value of $x_{2}$ at each crossing. We then simulate 6.7 from time
$t=0$ to time $t=10^{6}$ starting from $(x_{1},x_{2},x_{3})=(0,-20,0)$ and
measure the values of $x_{2}$ every time the simulated trajectory crosses the
Poincaré section. This produces $811$ values $x_{2,i}$, $i=1,\ldots,811$,
which we organize in $m=810$ data snapshots $(x_{2,i},x_{2,i+1})$
corresponding to measurements of the Poincaré return map. These datapoints are
plotted in the top-left panel of figure 9, which reveals that the Poincaré map
is unimodal. However, the exact mapping is unknown and can only be estimated
from data.

Figure 9: Poincaré map data and UPOs for the Rössler system 6.7. Each panel
has two subfigures. _Left subfigure:_ Data points on the Poincaré section,
plotted as $(x_{2,i},x_{2,i+n})$ to represent measurements of the $n$th
iterate of the Poincaré map with (a) $n=1$, (b) $n=2$, (c) $n=3$, and (d)
$n=4$. Intersections with the diagonal $(x_{2},x_{2})$ correspond to
period-$n$ fixed points of the Poincaré map, which are coloured according to
their minimal period. _Right subfigure:_ UPOs embedded in the continuous-time
Rössler attractor corresponding to the period-$n$ UPOs of the Poincaré map.

The UPOs of the Rössler attractor are in one-to-one correspondence with the
UPOs of the Poincaré map. We thus attempt to extract the latter by discovering
invariant atomic measures from the data collected in the Poincaré section. For
accuracy, we solve the SDP 3.8 with polynomial degrees $(k,\ell)=(20,80)$,
meaning that we are approximating the Lie derivative of polynomials in $x_{2}$
up to degree $k=20$ using polynomials up to degree $\ell=80$. Loosely
speaking, this corresponds to assuming the Poincaré map is well approximated
by a quartic polynomial map. As before, to improve numerical stability we
represent polynomials using the Chebyshev basis, rather than the monomial
basis. Then, to extract the periodic points from the Poincaré map we run the
optimization procedure 3.8 with 1000 randomized linear objective functions. We
are able to extract periodic points of periods 1 through 19, with the low-
period points plotted in figure 9. To confirm the accuracy of our method, the
figure also shows these points plotted against the forward iterates of the
Poincaré map. Observe that, as expected for periodic points, they lie
precisely along the $(x_{2},x_{2})$-diagonal. Importantly, by counting the
intersections of the Poincaré map data with the diagonal, we conclude we have
extracted _all_ UPOs of the Poincaré map with period up to at least 4.

Figure 10: A sample of long period UPOs of the Róssler attractor. Plotted UPOs
intersect the Poincaré section (a) 8 times, (b) 13 times, and (c) 17 times.

The periodic points in the Poincaré section can be further used as initial
conditions for root-finding methods that compute the corresponding UPOs for
the continuous-time Rössler system (see [3, Section 2.3] for details). Figure
9 shows all UPOs that intersect the Poincaré section up to four times. Figure
10 provides a small sample of longer UPOs that we are able to extract via the
same process.

A valid alternative to our data-driven discovery of atomic invariant measures
for the Poincaré map would be to perform model identification on the Poincaré
section data [6, 5, 3] and then numerically extract the periodic points from
the discovered symbolic mapping. This can perform well for low-period UPOs,
but small errors in the model discovery process accumulate through iterations
of the discovered map, making it difficult (if not impossible) to identify
UPOs with long period [6]. The advantage of our data-driven approach is that
we do not rely on iterates of an approximate Poincaré map to identify its
periodic points, so model discovery errors do not accumulate in the same way.

##  7 Conclusions

In this article, we have described a convex optimization approach to
approximate invariant measures for dynamical systems entirely from data. More
precisely, we have shown that the model-based convex optimization framework
for the approximation of invariant measures from [24] can be implemented in an
entirely data-driven manner by replacing Lie derivatives with data-driven
approximations obtained with EDMD. Our approach can be applied to a broad
range of dynamical systems, including those governed by deterministic or
stochastic differential equations and maps. Moreover, and crucially from a
practical point of view, the implementation of our method is the same
irrespective of whether the underlying dynamics are deterministic or
stochastic.

Compared to alternative data-driven approaches, such as Ulam’s method or the
EDMD-based strategy in [21], our method has the advantage of enabling one to
target a particular invariant measure through the choice of a suitable
optimization objective. We have demonstrated through a series of examples that
one can successfully approximate the physical measure of a chaotic attractor,
as well as ergodic invariant measures associated to different periodic orbits
embedded within it. In particular, the success in discovering UPOs for the
Rössler system suggests that our computational approach can be extremely
useful whenever UPOs are of interest. One example of practical relevance is
satellite mission design, since orbits that shadow UPOs require minimal use of
propellant. Orbital dynamics models such as the reduced three-body system,
which has a four-dimensional state space but a Hamiltonian structure that
allows for a two-dimensional Poincaré section, are only marginally more
complicated than the Rössler system studied in section 6.3. We believe they
can be studied successfully with our method, although this remains to be
confirmed.

Finally, although we have demonstrated that our approach works well in
practice, questions about its theoretical convergence remain. Known
convergence results for EDMD approximations of Lie derivatives [4] provide
only a heuristic justification for the good performance of our approach.
Rigorous statements require proving that optimizers of the SDP 3.8 converge to
optimizers of 2.11 as the number of data snapshots $m$ and the polynomial
degrees $k$ and $\ell$ of the EDMD dictionaries increases. We have done this
in section 4 when the data is governed by polynomial maps satisfying suitable
technical conditions. For general nonlinear maps, it should still be
relatively straightforward to study the limit of infinite data ($m\to\infty$)
for fixed $k$ and $\ell$. It seems however much harder to take $\ell\to\infty$
for fixed $k$, as this requires studying the fine convergence properties of
EDMD-based polynomial approximations to non-polynomial Lie derivatives. The
main technical challenges involved in this analysis are absent from our study
of polynomial maps in section 4 because in this case the infinite-data limit
of the matrix $L_{k\ell}$ coincides with the matrix $A_{k\ell}$ in 4.1 for all
$\ell\geq dk$, so our data-driven approach recovers (a straightforward
generalization of) the model-based strategy from [24] for a finite, albeit
unknown, value $\ell$. For general dynamics governed by non-polynomial maps,
however, one is forced to take $\ell$ to infinity. Identifying which problem
3.8 tends to in this limit and studying how this limiting problem behaves for
increasing $k$ remain two important open questions.

## Acknowledgments

JB was supported by an NSERC Discovery Grant. GF is grateful to Mayur Lakshmi
for pointing out the relevance of the methods in this paper for orbital
mechanics.

## References

  * [1] C. J. Bose and R. Murray.  Minimum ‘energy’ approximations of invariant measures for nonsingular transformations.  Discrete Contin. Dyn. Syst., 14(3):597–615, 2006.
  * [2] C. J. Bose and R. Murray.  Duality and the computation of approximate invariant densities for nonsingular transformations.  SIAM J. Optim., 18(2):691–709, 2007.
  * [3] J. J. Bramburger, S. L. Brunton, and J. Nathan Kutz.  Deep learning of conjugate mappings.  Phys. D, 427:Paper No. 133008, 16, 2021.
  * [4] J. J. Bramburger and G. Fantuzzi.  Auxiliary functions as koopman observables: Data-driven analysis of dynamical systems via polynomial optimization.  Journal of Nonlinear Science, 34(1):8, 2024.
  * [5] J. J. Bramburger and J. N. Kutz.  Poincaré maps for multiscale physics discovery and nonlinear Floquet theory.  Phys. D, 408:132479, 12, 2020.
  * [6] J. J. Bramburger, J. N. Kutz, and S. L. Brunton.  Data-driven stabilization of periodic orbits.  IEEE Access, 9:43504–43521, 2021.
  * [7] M. Brin and G. Stuck.  Introduction to dynamical systems.  Cambridge University Press, Cambridge, 2002.
  * [8] S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [9] S. I. Chernyshenko, P. Goulart, D. Huang, and A. Papachristodoulou.  Polynomial sum of squares in fluid dynamics: a review with a look ahead.  Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 372(2020):20130350, 18, 2014.
  * [10] E. de Klerk and M. Laurent.  A survey of semidefinite programming approaches to the generalized problem of moments and their error analysis.  In World women in mathematics 2018, volume 20 of Assoc. Women Math. Ser., pages 17–56. Springer, Cham, 2019.
  * [11] M. Dellnitz and O. Junge.  On the approximation of complicated dynamical behavior.  SIAM J. Numer. Anal., 36(2):491–515, 1999.
  * [12] T. A. Driscoll, N. Hale, and L. N. Trefethen.  Chebfun guide, 2014.
  * [13] G. Fantuzzi, D. Goluskin, D. Huang, and S. I. Chernyshenko.  Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization.  SIAM J. Appl. Dyn. Syst., 15(4):1962–1988, 2016.
  * [14] S. Gerber, S. Olsson, F. Noé, and I. Horenko.  A scalable approach to the computation of invariant measures for high-dimensional markovian systems.  Scientific Reports, 8(1):1796, 2018.
  * [15] D. Goluskin.  Bounding averages rigorously using semidefinite programming: mean moments of the lorenz system.  Journal of Nonlinear Science, 28(2):621–651, 2018.
  * [16] D. Goluskin and G. Fantuzzi.  Bounds on mean energy in the Kuramoto-Sivashinsky equation computed using semidefinite programming.  Nonlinearity, 32(5):1705–1730, 2019.
  * [17] D. Goswami, E. Thackray, and D. A. Paley.  Constrained Ulam dynamic mode decomposition: approximation of the Perron-Frobenius operator for deterministic and stochastic systems.  IEEE Control Syst. Lett., 2(4):809–814, 2018.
  * [18] D. Henrion.  Semidefinite characterisation of invariant measures for one-dimensional discrete dynamical systems.  Kybernetika (Prague), 48(6):1089–1099, 2012.
  * [19] M. V. Jakobson.  Absolutely continuous invariant measures for one-parameter families of one-dimensional maps.  Comm. Math. Phys., 81(1):39–88, 1981.
  * [20] A. Klenke.  Probability theory: A comprehensive course.  Universitext. Springer, Cham, 2020.
  * [21] S. Klus, P. Koltai, and C. Schütte.  On the numerical approximation of the Perron-Frobenius and Koopman operator.  J. Comput. Dyn., 3(1):51–79, 2016.
  * [22] S. Klus, F. Nüske, S. Peitz, J.-H. Niemann, C. Clementi, and C. Schütte.  Data-driven approximation of the Koopman generator: Model reduction, system identification, and control.  Phys. D, 406:132416, 15, 2020.
  * [23] S. Klus and C. Schütte.  Towards tensor-based methods for the numerical approximation of the Perron-Frobenius and Koopman operator.  J. Comput. Dyn., 3(2):139–161, 2016.
  * [24] M. Korda, D. Henrion, and I. Mezić.  Convex computation of extremal invariant measures of nonlinear dynamical systems and Markov processes.  J. Nonlinear Sci., 31(1):Paper No. 14, 26, 2021.
  * [25] M. Korda and I. Mezić.  On convergence of extended dynamic mode decomposition to the Koopman operator.  J. Nonlinear Sci., 28(2):687–710, 2018.
  * [26] U. Krengel.  On the speed of convergence in the ergodic theorem.  Monatshefte für Mathematik, 86(1):3–6, 1978.
  * [27] J. Kuntz, M. Ottobre, G.-B. Stan, and M. Barahona.  Bounding stationary averages of polynomial diffusions via semidefinite programming.  SIAM J. Sci. Comput., 38(6):A3891–A3920, 2016.
  * [28] M. V. Lakshmi, G. Fantuzzi, S. I. Chernyshenko, and D. Lasagna.  Finding unstable periodic orbits: a hybrid approach with polynomial optimization.  Phys. D, 427:Paper No. 133009, 11, 2021.
  * [29] M. V. Lakshmi, G. Fantuzzi, J. D. Fernández-Caballero, Y. Hwang, and S. I. Chernyshenko.  Finding extremal periodic orbits with polynomial optimization, with application to a nine-mode model of shear flow.  SIAM J. Appl. Dyn. Syst., 19(2):763–787, 2020.
  * [30] M. Laurent.  Sums of squares, moment matrices and optimization over polynomials.  In Emerging applications of algebraic geometry, volume 149 of IMA Vol. Math. Appl., pages 157–270. Springer, New York, 2009.
  * [31] J. Löfberg.  YALMIP: A toolbox for modeling and optimization in MATLAB.  In 2004 IEEE Int. Conf. Robotics and Automation, pages 284–289. IEEE, 2004.
  * [32] V. Magron, M. Forets, and D. Henrion.  Semidefinite approximations of invariant measures for polynomial systems.  Discrete Contin. Dyn. Syst. Ser. B, 24(12):6745–6770, 2019.
  * [33] A. Mosek.  The MOSEK optimization toolbox for MATLAB manual, 2015.
  * [34] G. W. Stewart.  On the Continuity of the Generalized Inverse.  SIAM J. Appl. Math., 17(1):33–45, 1969.
  * [35] S. M. Ulam.  A collection of mathematical problems.  Interscience Tracts in Pure and Applied Mathematics, no. 8. Interscience Publishers, New York-London, 1960.
  * [36] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley.  A data-driven approximation of the Koopman operator: extending dynamic mode decomposition.  J. Nonlinear Sci., 25(6):1307–1346, 2015.
  * [37] C. Zhang and E. Zuazua.  A quantitative analysis of koopman operator methods for system identification and predictions.  Comptes Rendus. Mécanique, 351(S1):1–31, 2023.