Milan Korda
email: korda@laas.fr

Didier Henrion
email: henrion@laas.fr

Igor Mezić
email: mezic@engineering.ucsb.edu

Convex computation of extremal invariant measures of nonlinear dynamical systems and Markov processes

Keywords: Invariant measure, convex optimization, ergodic optimization, physical measure

published or not. The documents may come

Convex computation of extremal invariant measures of nonlinear dynamical systems and Markov processes

Introduction

We propose a convex-optimization-based method for approximation of invariant measures. The method is based on the observation that the set of all invariant measures associated to a deterministic nonlinear dynamical system or a stochastic Markov process is given by the set of solutions to a linear equation in the space of Borel measures. The problem of finding an invariant measure can therefore be formulated as the feasibility of an infinitedimensional linear programming problem (LP). Adding an objective functional to this LP allows one to target a particular invariant measure (or a class of invariant measures) such as the physical measure, ergodic measures, absolutely continuous measures, atomic measures etc. The formulation is flexible in the sense that whenever a variational characterization of a given class of measures is known, then it can be used within the proposed framework. The approach is functional analytic in nature, by and large devoid of geometric or topological considerations. The only underlying assumption is that the dynamics is polynomial. This assumption is made for computational convenience even though the approach is far more general, applicable to any algebra of functions closed under function composition (in discretetime) or differentiation (in continuous-time).

The infinite-dimensional LP in the space of Borel measures is subsequently approximated along the lines of the classical Lasserre hierarchy [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] using a sequence of finite-dimensional convex semidefinite programming problems (SDPs). The optimal values of the SDPs are proven to converge from below to the optimal value of the infinite-dimensional LP (for the analysis of the speed of convergence in a related setting of optimal control, see [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for optimal control problems[END_REF]). The outcome of the SDP is an approximate truncated moment sequence of the invariant measure targeted; this sequence is proven to converge weakly to the moment sequence of the target invariant measure, provided this measure is unique (otherwise every accumulation point of the sequence corresponds to an invariant measure).

As a secondary contribution we describe a numerical procedure to approximate the support and density of the invariant measure using the truncated moment sequence obtained from the SDP. For the former, we provide confidence intervals enclosing, in the limit, a prescribed portion of the support; this is achieved using the Christoffel polynomial, an interesting object constructed from the Christoffel-Darboux kernel, which has already been utilized for support approximation in machine learning applications (e.g., [START_REF] Lasserre | The empirical christoffel function in statistics and machine learning[END_REF][START_REF] Pauwels | Sorting out typicality with the inverse moment matrix sos polynomial[END_REF]). For the latter, we construct a sequence of absolutely continuous measures with polynomial densities converging weakly to the target measure.

An interesting by product of the approach presented is the possibility to certify the nonexistence of an invariant measure, which is particularly pertinent for Markov processes. Such certification boils down to proving the emptiness of a spectrahedron defining the feasible set of the SDPs solved.

Finally, we also describe a generalization of the proposed approach to compute eigenmeasures of the Perron-Frobenius operator corresponding to a given (possibly complex) eigenvalue, with the invariant measures being a special case corresponding to eigenvalue one. This work is a continuation of the movement to apply convex optimization-based techniques to nonconvex problems arising from dynamical systems theory and control. For example, the related problem of invariant set computation was addressed in [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] whereas [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF][START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF][START_REF] Korda | Controller design and value function approximation for nonlinear dynamical systems[END_REF] addressed optimal control; [START_REF] Ozay | Convex certificates for model (in) validation of switched affine systems with unknown switches[END_REF], [START_REF] Ozay | Set membership identification of switched linear systems with known number of subsystems[END_REF] adressed model validation and switching system identification, respectively. The problem addressed here, i.e., invariant measure computation, was also addressed by this approach in [START_REF] Henrion | Semidefinite characterisation of invariant measures for one-dimensional discrete dynamical systems[END_REF] in one spatial dimension; this work can therefore be seen as a generalization of [START_REF] Henrion | Semidefinite characterisation of invariant measures for one-dimensional discrete dynamical systems[END_REF] to multiple dimensions and with a far more detailed theoretical and computational analysis. In the concurrent work [START_REF] Magron | Semidefinite characterization of invariant measures for polynomial systems[END_REF] the authors are also applying the Lasserre hierarchy for approximately computing invariant measures for polynomial dynamical systems, but there is no convex functional to be minimized and the focus is on distinguishing measures with different regularity properties (singular vs absolutely continuous).

Let us also mention the optimization-based approaches to invariant measure computation [START_REF] Bollt | The path towards a longer life: On invariant sets and the escape time landscape[END_REF][START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. These approaches are based on non-convex optimization and therefore have to deal with its inherent difficulties such as the existence of suboptimal local minimizers, saddle points or degeneracy. Therefore, contrary to the proposed convex-optimization based approach, these works do not provide convergence guarantees, despite being built on interesting ideas and showing promising practical performance.

The presented framework, where a convex user-specified functional is minimized among all invariant measures, can be seen as a generalization of and a computational method to carry out the so called ergodic optimization [START_REF] Jenkinson | Ergodic optimization[END_REF][START_REF] Jenkinson | Ergodic optimization in dynamical systems[END_REF][START_REF] Bochi | Ergodic optimization of Birkhoff averages and Lyapunov exponents[END_REF], where linear functionals are optimized among invariant measures (therefore leading to ergodic measures as the optimizers since these are the extreme points of the set of all invariant measures, hence the name ergodic optimization).

The presented approach based on optimization over Borel measures has a convex dual as an optimization over continuous functions that can be approximated by polynomial sumof-squares. This line of research has been investigated independently for various problems from dynamical systems (e.g., [START_REF] Chernyshenko | Polynomial sum of squares in fluid dynamics: a review with a look ahead[END_REF][START_REF] Fantuzzi | Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization[END_REF][START_REF] Goluskin | Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system[END_REF]). Of particular relevance to this work is [START_REF] Tobasco | Optimal bounds and extremal trajectories for time averages in dynamical systems[END_REF] which is dual to our approach in the continuous-time setting and when the objective functional in our approach is restricted to be linear. A by-product of our work is therefore an asymptotic convergence guarantee for the bounds obtained by [START_REF] Tobasco | Optimal bounds and extremal trajectories for time averages in dynamical systems[END_REF], provided that strong duality holds.

The paper is organized as follows. Section 2 formally states the problem of invariant measure computation. Section 3 describes the moment hierarchies and applies them to the invariant measure computation problem. Section 4 describes the reconstruction of the invariant measure from its moments. Section 5 discusses several concrete invariant measures to be targeted via the choice of objective functional to be optimized. Section 6 describes an extension to continuous time systems and Section 7 to Markov processes (both in discrete and continuous time) as well as discusses how to certify non-existence of invariant measures. Section 8 extends the method to eigenmeasures of he Perron-Frobenius operator and Section 9 presents numerical examples.

Problem statement

For concreteness we present the approach for deterministic discrete-time dynamical systems. The case of stochastic Markov processes is treated in Section 7; the continuous time cases are treated in Section 6 and 7.1.

Consider therefore a deterministic discrete-time nonlinear dynamical system

x + = T (x), (1)
where x ∈ R n is the state, x + ∈ R n is the successor state and each of the n components of the mapping T : R n → R n is assumed to be a multivariate polynomial.

An invariant measure for the dynamical system (1) is any nonnegative Borel measure µ satisfying the relation

µ(T -1 (A)) = µ(A) (2)
for all Borel measurable A ⊂ R n . In this paper we restrict our attention to invariant measures with support included in some compact set X ⊂ R n . With this assumption, the relation (2) reduces to

X f • T dµ = X f dµ (3)
for all f ∈ C(X).

When a measure µ is supported on a compact set X, it follows from the Stone-Weierstrass Theorem that it is entirely characterized by its moment sequence y := (y α ) α∈N n ∈ R ∞ , where

y α = X x α dµ(x), (4)
with x α := x α 1 1 • . . . • x αn n and α ∈ N n running over all n-tuples of nonnegative integers and with R ∞ denoting the space of all real-valued sequences. In particular, for the choice

f (x) = x α , relation (3) becomes X T α (x) dµ(x) = X x α dµ(x) (5)
for all α ∈ N n , where T α (x) := T 1 (x) α 1 • . . . • T n (x) αn . Since T is a polynomial, ( 5) is a linear constraint on the moments that can be written as

A(y) = 0, (6)
where A : R ∞ → R ∞ is a linear operator.

We remark that (6) characterizes all invariant measures associated to [START_REF] Bochi | Ergodic optimization of Birkhoff averages and Lyapunov exponents[END_REF] with support in X. In order to single out one invariant measure of interest we propose to use optimization.

In particular, we propose to solve the infinite-dimensional convex optimization problem min y∈M (X)

F (y) s.t. A(y) = 0 y 0 = 1 (7)
where the minimization is w.r.t. a sequence y belonging to the convex cone

M (X) := {y ∈ R ∞ : y α = X x α dµ(x), α ∈ N n , µ ∈ M(X) + }
of moments of non-negative Borel measures on X, the objective functional F : R ∞ → R is convex, and the constraint y 0 = 1 is a normalization constraint enforcing that the measure is a probability measure.

We note that, since T is polynomial and hence continuous and X is compact, the Krylov-Bogolyubov Theorem ensures that there exists at least one invariant measure for (1) and hence the optimization problem ( 7) is always feasible.

Remark 1 (Role of the objective function) The role of F is to target or single out a specific invariant measure from the set of all invariant probability measures characterized by the constraints of [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. In principle, F can be any convex functional that facilities this. In particular, it can be extended-valued (i.e. equal to +∞), therefore encoding any constraints of interest such as µ being absolutely continuous or singular w.r.t. to a given measure. See Section 5 for concrete choices of F .

A typical example encountered in practice for the choice of F is

F (y) = |α|≤d (y α -z α ) 2 , (8)
where (z α ) |α|≤d is a given finite vector of moments of total degree no more than d. The moments z α can be estimates of the first few moments of the invariant measure that we wish to compute obtained, e.g., from observed data or by analytical reasoning (e.g., based on the symmetries of the problem). The optimization problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] then seeks among all invariant measures µ the one which minimizes the discrepancy between the first moments (y α ) |α|≤d of µ and the given moments (z α ) |α|≤d . Of course, the least-squares criterion in [START_REF] Goluskin | Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system[END_REF] can be replaced by other convex metrics measuring the discrepancy between two truncated moment sequences.

Moment hierarchy

In this section we describe a hierarchy of finite dimensional convex optimization problems approximating the infinite-dimensional problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] and prove that the solutions obtained from these approximations converge to a solution of [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. For this we assume that the objective function F (y) in ( 7) depends only on finitely many moments (y α ) |α|≤d .

The finite-dimensional approximations are derived from the so-called Lasserre hierarchy of approximations to the moment cone. In particular, we use semidefinite programming representable outer approximations to this cone and in addition we truncate the first equality constraint of (7) by imposing the linear constraint (5) only for f (x) = x α , |α| ≤ k, i.e., for all monomials of degree no more than k. By linearity of the constraint, this implies that the constraint is satisfied for all polynomials of degree no more than k. The degree k is called relaxation degree.

Before writing down the finite-dimensional approximation of ( 7), we first describe the construction of the finite-dimensional outer approximations to M (X).

Finite-dimensional approximations of the moment cone

Here we describe the semidefinite-programming representable outer approximation to M (X).

For this we assume that the compact set X is of the form1

X := {x ∈ R n | g i (x) ≥ 0, i = 1, . . . , n g } (9)
with g i being multivariate polynomials. Let us denote the unit polynomial by g 0 (x) := 1.

The outer approximation M d (X) of degree d, d even, is

M d (X) := {y ∈ R ∞ : M d (g i y) 0, i = 0, 1, . . . , n g }, (10)
where • 0 denotes positive semidefiniteness of a matrix and M d (g i y) are the so-called localizing moment matrices, to be defined below. The convex cone M d (X) is an outer approximation to M (X) in the sense that for any non-negative measure µ on X the moment vector y of µ belongs to M d (X).

The localizing moment matrices M d (g i y) are defined by

M d (g i y) = l d y (g i v d i v d i ), (11)
where

d i = (d -deg g i )/2 , v d (x) := (x α ) |α|≤d
and where the Riesz functional

d y : R[x] d → R is defined for any f = |α|≤d f α x α by l d y (f ) = |α|≤d f α y α . (12)
This functional mimicks integration with respect to a measure; in particular when y ∈ R ( n+d d ) is a truncated moment vector of a measure µ, then l d y (f ) = f dµ for any f ∈ R[x] d . In this case, the localizing matrices M d (g i y) are necessarily positive semidefinite, corresponding to the fact that g i p 2 dµ ≥ 0 for any polynomial p. Importantly, the following fundamental converse result states that if localizing moment matrices are positive semidefinite for all d ∈ N, then necessarily y is a moment vector of a nonnegative measure supported on X. In order for this to hold, the defining polynomials g i need to satisfy the so-called Archimedian condition:

Assumption 1 There exist polynomials σ i and a constant r ∈ R such that

r 2 -x x = ng i=0 σ i g i , (13)
where σ i = j p 2 i,j with p i,j polynomial, i.e., each σ i is a sum of squares of other polynomials.

Assumption 1 is an algebraic certificate of compactness of X because this assumption implies that X ⊂ {x | x x ≤ r 2 }. Since X is assumed compact, this is a non-restrictive condition because a redundant constraint of the form g i = r 2 -x x can always be added to the definition of X for a sufficiently large r in which case Assumption 1 is satisfied trivially.

Theorem 1 (Putinar [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]) Suppose that Assumption 1 holds and that y ∈ M d (X) for all d ≥ 0. Then there exists a unique non-negative measure µ on X such that (4) holds for all α ∈ N n .

Approximation of the infinite-dimensional convex problem

Now we are ready to write down the finite-dimensional approximation to [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. The first step in the approximation is to impose the equality constraint of ( 7) only for all monomials of total degree no more than k instead of for all continuous functions. That is, we impose,

X T α (x)dµ = X x α dµ = 0, |α| ≤ k, where T α (x) = T 1 (x) α 1 • . . . • T n (x)
αn . This is a set of k linear equations and since T is polynomial it can be re-written in terms of the truncated moment sequence of the measure µ of degree no more than

d k := k deg T = k max i=1,...,n deg T i .
In terms of the Riesz functional, this set of equalities becomes

l d k y (T α (x)) = l d k y (x α ), |α| ≤ k. ( 14
)
This set of equalities can be re-written in a matrix form as

A k y = 0, for some matrix A k ∈ R ( n+k k )×( n+d k d k )
where y is the moment sequence of µ.

The equality constraint X 1 dµ translates to y 0 = 1. The conic constraint y ∈ M (X) is replaced, according to the previous section, by the constraint y ∈ M d k (X). This leads to the following finite-dimensional relaxation of order k of the infinite-dimensional problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] min y∈R

( n+d k d k ) F (y) s.t. Ay = 0 y 0 = 1 M d k (g i y) 0, ∀ i = 0, 1, . . . , n g . (15)
In optimization problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF], a convex function is minimized over a convex semidefiniteprogramming representable set and hence ( 7) is a convex optimization problem. Provided that the objective functional F (y) is also semidefinite programming representable (e.g., it is of the form (8)), then the problem ( 7) is a semidefinite programming problem and hence can be readily solved by off-the-shelf software (e.g., MOSEK or SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a Matlab toolbox for optimization over symmetric cones[END_REF]). Importantly, the finite-dimensional relaxation (15) can be derived from the abstract form [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] and passed to a selected SDP solver automatically with the help of the modelling software Gloptipoly 3 [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] and Yalmip [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in Matlab[END_REF].

An immediate observation is that problem ( 15) is a relaxation of problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] in the sense that the moment sequence of any measure feasible in (7) truncated up to degree d k is feasible in [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. Therefore in particular for any k, the optimal value of ( 15) provides a lower bound on the optimal value of (7). In the following section we study the convergence of these lower bounds to the optimal value of ( 7) as as well as convergence of the minimizers of ( 15) to a minimizer of (7).

Convergence of approximations

In this section we prove convergence of the finite-dimensional approximations [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] to a solution to the infinite-dimensional optimization problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF].

Theorem 2 Suppose that Assumption 1 holds, that the function F is lower semi-continuous2

and let y k ∈ R ( n+d k d k ) denote an optimal solution to (15) and p the optimal value of (7). Then the following holds:

1. lim k→∞ F (y k ) = p 2.
There exists a subsequence (k i ) ∞ i=1 such that y k i converges pointwise to a moment sequence of an invariant measure µ attaining the minimum in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF].

3. In particular, if there is a unique invariant measure µ attaining the minimum in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF], then y k converges to the moment sequence of µ .

Proof:

The proof follows a standard argument (see, e.g., [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]). By Assumption 1, for every ᾱ ≥ 0 there exists a k 0 such that |y α | ≤ M |α| for any vector y satisfying the constraints of [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] for k ≥ k 0 and for any α ∈ N n satisfying |α| ≤ ᾱ, where M is the constant from ( 13). This statement implies that each component of y k is bounded for sufficiently large k. To see this, let f = σ 0 + i σ i g i be the right-hand-side polynomial from [START_REF] Jenkinson | Ergodic optimization[END_REF] and let k 0 be the smallest number

d ≥ 1 such that d deg T ≥ deg σ 0 /2 and d deg T ≥ (σ i -deg g i )/2 . Then necessarily d k y (σ 0 ) ≥ 0 and d k y (σ i g i ) ≥ 0 for any k ≥ k 0 .
Therefore by linearity we have

d k y (f ) ≥ 0. But since f = r 2 -x x by (13) we get d k y (r 2 ) ≥ d k y (x x) = n i=1 d k y (x i 2 ). Since d k y (r 2 ) = r 2 d k y (1
) = r 2 y 0 and y 0 = 1 by the second constraint of (15) and since

d k y (x i 2 ) ≥ 0, we conclude that d k y (x i 2 ) ∈ [0, r 2 ].
Proceeding recursively, applying the same reasoning to

x α f with α i even for all α ∈ N n satisfying |α| ≤ ᾱ := max{α ∈ N n | |α|/2 + (deg σ i )/2 ≤ (d k -deg g)/2
}, we conclude that all even moments y α , |α| ≤ ᾱ lie in [0, r 2|α| ]. Since even moments are on the diagonal of the matrix M d (y) 0 and since the off-diagonal elements of a positive semidefinite matrix are bounded in magnitude by the diagonal elements, the conclusion follows.

Having established that lim sup k→∞ |y k α | < ∞ for each α ∈ N n , it follows using a standard diagonal argument that we can extract a subsequence y k i satisfying, for each α ∈ N n , lim i→∞ y k i α = y α with y α ∈ R. To conclude the proof it remains to show that y is a moment sequence of a measure attaining the minimum in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. Using Theorem 1, it follows that y is a moment sequence of a non-negative measure µ on X since M d (g i y ) 0 by continuity of the mapping M → λ min (M ), where λ min (M ) denotes the minimum eigenvalue of a symmetric matrix (or equivalently by closedness of the cone of positive semidefinite matrices). In addition, y satisfies the equality constraints of ( 7) by continuity since each row of the matrix A has only finitely many non-zero elements. Therefore µ is an invariant measure. Finally, since ( 15) is a relaxation of [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] we have F (y k i ) ≤ p . By the lower semicontinuity of F we also have F (y ) ≤ lim i→∞ F (y k i ) ≤ p and hence necessarily F (y ) = p since µ is feasible in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] and therefore F (y ) ≥ p .

Reconstruction of measure from moments

In this section we show how the solutions to the finite dimensional relaxations [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] in the form of a truncated moment sequence can be used to approximately reconstruct the invariant measure. In particular we show how to approximate the support of the measure and how to construct a sequence of absolutely continuous measures converging weakly to the invariant measure. In this section, we assume that the optimal invariant measure is unique, which holds generically (in the Baire category sense) by [START_REF] Jenkinson | Ergodic optimization[END_REF]Theorem 3.2].

Assumption 2 (Unique invariant measure) Convex problem [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] has a unique solution denoted by µ .

Approximation of the support

In this section we show how the solutions to the finite dimensional relaxations (15) can be used to approximate the support of the invariant measure µ . The approximations constructed here aim at enclosing a certain prescribed portion of the support. Guaranteed outer approximates to the global attractor (on which certain invariant measures are supported) can be computed using the approach of [START_REF] Schlosser | Converging outer approximations to global attractors using semidefinite programming[END_REF].

In order to construct the approximations we utilize a certain polynomial constructed for a vector of moments of a given measure. Assume that we are given the sequence of moments y of a non-negative measure µ on X. Then from the truncated moments of degree up to d, d even, we define the Christoffel polynomial

q y d (x) = v d/2 (x) M d (y) -1 v d/2 (x), (16)
where v d/2 (•) is the basis vector of all monomials up to degree d/2 with the same ordering as the vector of moments y. The polynomial q y d is well defined as long the moment matrix M d (y) is invertible, which is satisfied if and only if µ is not supported on the zero level set of a polynomial of degree d/2 or less.

The sublevel sets of the polynomial q y d (x) have a remarkable property of approximating the shape of the support of the measure µ. In the real multivariate domain, this was observed recently for empirical measures (sums of Dirac masses) in [START_REF] Pauwels | Sorting out typicality with the inverse moment matrix sos polynomial[END_REF] and subsequently studied analytically for measures with certain regularity properties in [START_REF] Lasserre | The empirical christoffel function in statistics and machine learning[END_REF]. In the complex domain, the theory is far more developed; see, e.g., [START_REF] Gustafsson | Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction[END_REF] and references therein.

Here we use the following simple result which holds for arbitrary probability measures µ.

Lemma 1 Let µ be a probability measure on X with moment sequence y, let ∈ [0, 1) be given and let

γ = 1 1- n+d/2 n
and assume that M d (y) is invertible. Then

µ({x : q y d (x) ≤ γ }) ≥ . ( 17
)
Proof: We bound the complementary event:

µ({x : q y d (x) > γ }) = {x:q y d (x)>γ } 1 dµ = {x:γ -1 q y d (x)>1} 1 dµ ≤ γ -1 X q y d (x) dµ(x),
where we have used the fact that q y d is nonnegative. Using the definition of q y d we get

X q y d (x) dµ(x) = X v d/2 (x) M d (y) -1 v d/2 (x) dµ(x) = trace M d (y) -1 X v d/2 (x)v d/2 (x) dµ(x) = trace{M d (y) -1 M d (y)} = trace I ( n+d/2 n ) = n + d/2 n ,
where I k denotes the identity matrix of size k. Therefore

µ({x : q y d (x) > γ ≤ 1 γ n + d/2 n
and the result follows since µ({x : q y d (x) ≤ γ }) = 1 -µ({x : q y d (x) > γ }). Lemma 1 can be readily used to construct approximations to the support of µ in the form

x : q y k d (x) ≤ γ where y k ∈ R ( n+d k d k
) is a solution to the kth order relaxation [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. We remark that y k only approximates the moments of µ , with guaranteed convergence by Theorem 2, and hence the bound from Lemma 1 may be violated for finite k.

Weakly converging approximations

In this section we show how the to construct a sequence of absolutely continuous measures (w.r.t. the Lebesgue measure) converging weakly to the invariant measure. This is especially useful if in fact the invariant measure possesses a density with respect to the Lebesgue measure although the approach is general and always provides a sequence of signed measures with polynomial densities that converges weakly to the invariant measure. The idea is simple: given a vector y ∈ R ( n+d n ) (e.g., a truncated moment vector of a measure), we can always represent the Riesz functional d y

: R[x] d → R as d y (p) = X p(x)q(x)dx (18)
for some polynomial q ∈ R[x] d , provided that the set X has a nonempty interior. Indeed, by linearity it suffices to satisfy (18) for p(x) = x α , |α| ≤ d, which leads to a system of linear equations

M L d q = y, ( 19
)
where q is the coefficient vector of polynomial q in the monomial basis with the same ordering as the vector of moments y, and M L d is the moment matrix of the Lebesgue measure on X of degree 2d, i.e.,

M L d := X v d (x)v d (x) dx.
Provided that the interior of X is nonempty, matrix M L d is invertible and hence the linear system of equations ( 19) has a unique solution q = (M L d ) -1 y. This approach applied to the solutions of (15) leads to the following result: Theorem 3 Let Assumption 2 hold, suppose that the interior of X is nonempty, denote y k ∈ R ( n+d k n ) any solution to [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF], and let

q k := (M L d k ) -1 y k . ( 20
)
Then the signed measures with densities q k (x) := v k (x)q k ∈ R[x] k with respect to the Lebesgue measure, converge weakly star on X to the invariant measure µ .

Proof: Verifying weak star convergence means that lim k→∞ X f q k dx = X f dµ for all f ∈ C(X). Since X is compact, it is enough to verify this relationship for all f of the form f = x α , α ∈ N n , which forms a basis of the space of all polynomials, which is a dense subspace of C(X). By construction we have X x α q k = y k α for k ≥ |α| and hence lim k→∞ X

x α q k (x) dx = lim

k→∞ y k α = y α = X x α dµ (x)
by Theorem 2, part 3.

Theorem 3 says that the density approximations constructed from the solutions to (15) using [START_REF] Lasota | Chaos, fractals, and noise -Stochastic aspects of dynamics[END_REF] converge in the weak star topology to the invariant measure optimal in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. From a practical point of view, e.g., for the purpose of visualization, we recommend using polynomial densities with coefficients

(M L k ) -1 ȳk , (21)
where ȳk ∈ R ( n+k n ) is the vector of the first n+k n elements of y k , rather than the full vector

y k ∈ R ( n+d k n )
. This is because the invariance constraint ( 14) is imposed only for all monomials up to degree k and hence moments of degrees higher than k are less constrained in [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] and hence are likely to be less accurate approximations to the true moments.

Applications -choosing the objective function

In this section we list a several classes of measures that can be targeted through the choice of the objective function F of convex problem (7).

Physical measures

Here we describe how the proposed methodology can be used to compute the moments of physical measures. Let Assumption 2 hold so that there is a unique physical measure µ with support included in X. Therefore for Lebesgue almost every x ∈ X for which the trajectory of (1) originating from x stays in X we have for any f ∈ C(X)

X f dµ = lim N →∞ 1 N N i=1 f (T i (x)).
Selecting f (x) = x α , we can approximately compute the moments of µ as

y num α = 1 N N i=1 f (T i (x)) ≈ X x α dµ (22)
with some N 1 and α running over a selected subset of multiindices I. The idea is that the number of moments we (inaccurately) compute using ( 22) is very small and then we use the optimization problem [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] to compute a much larger number of moments of µ, just from the information contained in {y num α : α ∈ I}. This is achieved by setting the objective function in [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] to

F (y) = α∈I (y α -y num α ) 2 (23)
or any other metric measuring the discrepancy among moments.

Ergodic measures

In this section we describe how one can target ergodic measures through the choice of the objective function F . In particular, these can be use used to locate embedded unstable fixed points or periodic orbits. The starting point is the well-known fact that ergodic measures are precisely the extreme points of the set of all invariant probability measures (see, e.g., [START_REF] Phelps | Lectures on Choquet's theorem[END_REF]Proposition 12.4]; see also [START_REF] Cross | Approximating extreme points of infinite dimensional convex sets[END_REF] for general results on the behavior of extreme points under projection in infinite-dimensional vector spaces). Therefore, we can target ergodic measures by selecting a linear objective functional in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] since the minimum of a linear program is attained at an extreme point, provided the minimizer is unique. This leads to the following immediate result:

Theorem 4 Let F be a continuous linear functional such that Assumption 2 holds and let

y k ∈ R ( n+d k d k
) denote the optimal solution to [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. Then y k converges to a moment sequence of an ergodic measure µ, provided Assumption 1 holds.

Proof: The result follows from Theorem 2 and the fact that the minimizer is unique, therefore necessarily an extreme point of the feasible set of [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF].

A typical choice of the objective functional is F (y) = α c α y α with c having only finitely many non-zero elements (i.e., F is a linear combination of a finite number of moments).

Absolutely continuous measures

In this section we describe how to target measures µ absolutely continuous w.r.t. a given measure ν through the choice of the objective function F . In most practical applications the measure ν will be the Lebesgue measure. The absolutely continuity of µ w.r.t. ν is equivalent to the existence of a density ρ such that µ = ρdν. Assuming that ρ(x) ≤ γ < ∞ for ν-almost all x ∈ X, we can impose the absolute continuity constraint by choosing

F (µ) = I {µ | µ≤γν} ,
where I A is the extended-value characteristic of a set A (i.e., I A (x) = 0 if x ∈ A and I A = +∞ if x / ∈ A). Then F (y) is a convex function of y and its domain is the set of moment vectors y such that

M d (γz -y) 0 ∀d ∈ N, (24)
where z denotes the moment sequence of ν. This is added to the constraints of ( 7) and ( 15) (for a particular fixed d in the latter case).

We remark that, in general, there may be multiple absolutely continuous measures in which case the constraint (24) can be combined with an additional choice of the objective functional F in order to target a specific measure (e.g., the physical measure as in Section 5.1).

Singular measures

In this section we describe how to target singular measures through the choice of the objective function F . In this work, we focus on atomic measures only. Whether there exists and effective variational characterization of the elusive singular-continuous measures (e.g., the Cantor measure) remains an open problem, to the best of our knowledge.

Atomic measures are of interested because they may correspond to unstable periodic orbits or fixed points, which are difficult to obtain using simulation-based techniques. These structures are typically ergodic and therefore the methods of this section can be combined with those of Section 5.2.

We use the observation that if a measure µ consists of K atoms, then the associated moment matrix M d (y) is of rank at most K for all d and of rank exactly K for sufficiently large d. Therefore, in order to seek an invariant measure consisting of K atoms (e.g., a K-period orbit) one would in principle want to choose

F (y) = I {y | rank M d (y)=K} .
Unfortunately, such F is not convex and therefore we propose to use a convex relaxation

F (y) = I {y | trace M d (y)≤γ} ,
where γ ≥ 0 is a regularization parameter. This translates to the constraint trace M d (y) ≤ γ which is added to the constraints of ( 7) and [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. We note that since M d (y) is positive semidefinite, its trace coincides with its nuclear norm, which is a standard proxy for rank minimization [START_REF] Fazel | Matrix rank minimization with applications[END_REF] (since the convex hull of the rank is the trace on the unit ball of symmetric matrices).

Continuous time version

In this section we briefly outline how the presented approach extends to continuous time.

Assume therefore that we are dealing with the dynamical system of the form

ẋ = b(x), (25)
each component of the vector field b is assumed to be a multivariate polynomial. We are seeking a non-negative measure µ on X invariant under the flow of dynamical system [START_REF] Lasserre | The empirical christoffel function in statistics and machine learning[END_REF], which is equivalent to the condition

X gradf • b dµ = 0 (26)
for all f ∈ C 1 (X).

The infinite-dimensional convex optimization problem (7) then becomes

min µ∈M + (X) F (µ) s.t. X gradf • b dµ = 0 ∀f ∈ C(X) X dµ = 1. (27)
This optimization problem is then approximated by taking f = x α , α ∈ N n , and proceeding in exactly the same way as described in Section 3.2, leading to a finite-dimensional relaxation of the same form as [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF], with the same convergence results of Theorems 2 and 3.

Markov processes

In this section we describe a generalization to Markov processes evolving on the state-space X. We assume a Markov chain in the state-space form3

x k+1 = T (x k , w k ),
where (w k ) ∞ k=0 is a sequence of independent identically distributed random variables with values in a given set W and the mapping T is assumed to be a polynomial in (x, w). The distribution of the random variables w k is denoted by P w , i.e., for all Borel A ⊂ X, P w (A) is the probability that w k ∈ A.

The condition for a probability measure µ to be invariant then reads

X W f (T (x, w)) dP w (w) dµ(x) = X f (x) dµ(x) (28)
for all f ∈ C(X). This equation is linear in µ and can therefore be used in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] instead of the first equality constraint. The approach then proceeds along the steps of Section 3.2, i.e., we set f (x) = x α and enforce ( 28) for all such f with |α| ≤ k, leading to

X W T α (x, w) dP w (w) dµ(x) = X x α dµ, |α| ≤ k. ( 29
)
Since T is a polynomial in (x, w), T α can be written as

T α (x, w) = β,γ t α,β,γ x β w γ
for some coefficients t α,β,γ . Therefore,

W T α (x, w)dP w (w) = β,γ t α,β,γ x β W w γ dP w (w) = β,γ t α,β,γ x β m γ ,
where

m γ = W w γ dP w (w)
are the moments of w k , which are fixed numbers that can be either precomputed analytically or using sampling techniques. The equation ( 29) can therefore be re-written as

β,γ t α,β,γ m γ X x β dµ(x) = X x α dµ(x), dµ, |α| ≤ k,
or, in terms of the moments of µ (4),

β,γ t α,β,γ m γ y β = y α , |α| ≤ k.
This is a finite-dimensional system of linear equations of the form

Ay = 0.
Adding the normalization constraint y 0 = 1 and the positive-semidefiniteness constraints M d k (y) 0 and M d k (y, g i ) 0 leads to the optimization problem [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. The objective functional F (y) of ( 15) is again chosen in order to target a particular class of invariant measures. The same convergence guarantees of Theorems 2 and 3 hold.

Remark 2 (Uniqueness) It is interesting to note that, in the presence of randomness, it is much more common for a unique invariant measure to exist. For example, a sufficient condition for this is the recurrence of the Markov chain; see [START_REF] Meyn | Markov chains and stochastic stability[END_REF]Chapter 10] for more details.

for every Borel measurable set A ⊂ X. Here M c (X) stands for the vector space of all complex-valued measures on X. Given λ ∈ C, a complex-valued measure µ is an eigenmeasure of P if f • T dµ = λ f dµ [START_REF] Ozay | Convex certificates for model (in) validation of switched affine systems with unknown switches[END_REF] for all f ∈ C(X). Since the moment-based approach developed previously applies to nonnegative measures, we use the Jordan decomposition

µ = µ + R -µ - R + i(µ + I -µ - I ), where i = √ -1 is the imaginary unit and µ + R ∈ M (X), µ - R ∈ M (X), µ + I ∈ M (X), µ - I ∈ M (X). Similarly, we write λ = λ R + iλ I .
Then the condition ( 31) is equivalent to

(f • T -λ R f ) dµ + R -(f • T -λ R f ) dµ - R + λ I f dµ - I -λ I f dµ + I = 0, ( 32a
) (f • T -λ R f ) dµ + I -(f • T -λ R f ) dµ - I + λ I f dµ - R -λ I f dµ + R = 0 (32b)
for all f ∈ C(X). When expressed with f = x α , α ∈ N n , this is equivalent to

A(y + R , y - R , y + I , y - I ) =
y + R ,y - R ,y + I ,y - I F (y + R , y - R , y + I , y - I ) s.t. A(y + R , y - R , y + I , y - I ) = 0 (y + R ) 0 + (y - R ) 0 + (y + I ) 0 + (y - I ) 0 = 1 y + R ∈ M (X), y - R ∈ M (X), y + I ∈ M (X), y - I ∈ M (X), (33)
which is then approximated by a sequence of finite-dimensional SDPs in exactly the same fashion as described in Section 3.2, with the convergence results of Theorem 2 also holding in this setting.

9 Numerical examples

Logistic map

As our first example we consider the Logistic map:

x + = 2x 2 -1 on the set X = [-1, 1] = {x ∈ R : (x + 1)(1 -x) ≥ 0}.

Physical measure

First, we compute the moments of the unique physical measure µ on X. For this example, the density of the physical measure is given by [START_REF] Lasota | Chaos, fractals, and noise -Stochastic aspects of dynamics[END_REF] ρ

(x) = 1 π 1 √ 1 -x 2 .
We used the first moment

1 π 1 -1
x √ 1-x 2 dx = 0 as data input for the objective function of the form [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF], i.e., we set F (y) = (y 1 -0) 2 (note that instead of an exactly computed value we could have used an imprecise value of the moment from a simulation or determine this value based on symmetry without analytically integrating the density). Then we solve [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] with k ∈ {5, 10, 100} and compute a degree k polynomial approximation to the density using [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF]; as in [START_REF] Henrion | Semidefinite characterisation of invariant measures for one-dimensional discrete dynamical systems[END_REF], for numerical stability reasons, we work in the Chebyshev basis rather than the monomial basis (i.e., we express the constraint ( 14) using the Chebyshev basis polynomials instead of monomials x α and replace the monomial basis vectors v d i by the vectors of Chebyshev polynomials up to degree d i in ( 11)); see [START_REF] Henrion | Semidefinite characterisation of invariant measures for one-dimensional discrete dynamical systems[END_REF] for more details on the use of Chebyshev polynomials in this context. In Figure 1 we compare the true density and the polynomial approximations. We observe a very good fit even for low-degree approximations and oscillations of the sign of the approximation error, akin to classical results from approximation theory. The computed moments (transformed to the monomial basis) are compared in Table 1; we see a very good match.

Hénon map

Our second example is the Hénon map:

x

+ 1 = 1 -1.4x 2 1 + x 2 x + 2 = 0.3x 1 .
The set X is the box [-1.5, 1.5] × [-0.4, 0.4] expressed as

X = {x | (x 1 -1.5)(1.5 -x 1 ) ≥ 0, (x 2 -0.4)(0.4 -x 2 ) ≥ 0}.
The goal is to compute the moments of the physical measure µ on X. The available information is only the first moment y num (1,0) (i.e., the expectation of the first coordinate) of µ approximately computed using [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] with f (x) = x 1 . The objective function is F (y) = (y (1,0) -y num (1,0) ) 2 . Table 2 compares the moments returned by relaxation (15) of order k = 10 with moments computed numerically using [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]; we observe a good agreement. Figure 3 then compares the support approximations computed using using the Christoffel polynomial [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] with y = y k . Table 2: Hénon map: comparison of moments computed by the SDP relaxation (15) and using [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF].

Only the first moment corresponding to x 1 was given as data input to [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF].

Moments

x

1 x 2 x 2 1 x 1 x 2 x 2 2 x 3 1 x 2 1 x 2 x 1 x 2 2
x 3 2 SDP 0.2570 0.0771 0.5858 -0.0379 0.0527 0.2468 0.0131 -0.0140 0.0067 Numeric 0.2570 0.0771 0.5858 -0.0291 0.0527 0.2320 0.0510 -0.0174 0.0063

Lorenz system

Next, we consider the classical Lorenz system ẋ1 = 10(x

2 -x 1 ) ẋ2 = x 1 (28 -x 3 ) -x 2 ẋ3 = x 1 x 2 - 8 3 x 3
scaled by the linear coordinate transformation x = diag([1/25, 1/30, 1/50])x. The goal is to compute the moments of the physical measure µ which is supported on the Lorenz attractor.

Here we investigate the effect of the amount of information available for the computation in terms of the number of moments used in the objective function

F (y) = N i=1 (y i - y num i ) 2
, where y num i are approximations of the true moments computed using simulation. We investigate as well the effect of the accuracy of the approximate moments provided by varying the simulation length M . The accuracy is measured in terms of the percentage root mean square error on the sequence of the first 56 moments (i.e., moments up to degree 5) with respect the numerical approximations from a simulation of length 10 6 . Figure 2 shows the results for M = 10 2 and M = 10 6 with N ∈ {1, . . . , 10}. We observe a rather small impact of the accuracy of the moments provided (determined by M ); on the other hand, the number of moments provided N plays a significant role, resulting in an error of roughly 15% with one moment provided and error of roughly 0.05% with seven or more moments.

Stochastic processes

OrnsteinUhlenbeck To demonstrate the approach of Section 7.1, we first choose the classical OrnsteinUhlenbeck process given by dx(t) = -ax(t)dt + σdW (t), where W (t) is the Wiener process. For a > 0, this process has a unique invariant measure equal to N (0, σ 2 2θ ), i.e., the normal distribution with zero mean and standard deviation σ √ 2θ . In order to compute the moments, we solved (15) without the objective function (i.e., as a feasiblity problem). The results for a = 1 and σ = 1 and d = 10 are shown in Table 3; we observe a precise match. Nonlinear drift Next, we test the approach on a stochastic differential equation with a nonlinear drift dx(t) = -ax 3 (t)dt + σdW (t). For a > 0, this stochastic process is uniformly ergodic and hence admits a unique invariant measure [START_REF] Meyn | Markov chains and stochastic stability[END_REF]. Table 5 shows the results of solving (15) for d = 10 in comparison with Monte Carlo simulation using the Euler forward discretization and 10 5 sample trajectories. Detecting nonexistence of invariant measures In this section we demonstrate the ability of the approach to numerically prove the non-existence of an invariant measure as described in Section 7.1. For this, we consider the process dx(t) = -ax 3 (t)dt + (1 + x 2 (t))dW (t).

Clearly, the process does not admit an invariant measure for a < 0. The situation is more interesting for a > 0. We carried out computations for a ∈ {0, 1, . . . , 10}, solving the SDP relaxation [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]. The results indicated that no invariant measure exists for this range of values of a, although the degree to certify it (i.e., to render (15) infeasible) increases with a. This is intuitive since, vaguely speaking, the larger the value of a, the "more stable" the process is. Table 5 shows this dependence; it appears that the minimum degree obeys the linear relation d = 2a + 2, although we did not attempt to prove this analytically. It should be noted that the infeasibility of (15) was decided by the interior point solver SeDuMi; in particular we did not use rigorous certification tools such as [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF] and therefore these results should be understood as a strong numerical evidence of non-existence rather than a rigorous proof.

Conclusion

This work presented a convex-optimization-based method for computation of invariant measures for continuous and discrete time deterministic and stochastic systems. We described how to cast the problem of invariant set computation as an infinite-dimensional LP in the space of Borel measures and how to target particular invariant measures by the choice of an objective functional. We showed how this infinite-dimensional LP can be approximated by a sequence of finite-dimensional SDPs with a guaranteed asymptotic convergence and how the results of this SDP can be used for support approximation of the invariant measure using the Christoffel-Darboux kernel. Interesting by-products of the approach are a method to certify non-existence of invariant measures and a method to compute eigenmeasures of the Perron-Frobenius operator.

Future work should focus on improving the scalability of the approach by exploit sparsity or symmetries of the problem at hand or on developing a data-driven counterpart of the approach where the model is unknown and only finite collection of observations is available, in the spirit of [START_REF] Korda | Computing controlled invariant sets from data using convex optimization[END_REF].

  0, where y + R etc are the moment sequences of the respective measures and A is a linear operator. Coupled with the normalization constraint (y + R ) 0 + (y - R ) 0 + (y + I ) 0 + (y - I ) 0 = 1 and and objective functional F (y + R , y - R , y + I , y - I ), we arrive at an infinite dimensional linear programming problem min

Figure 1 :

 1 Figure 1: Logistic map: approximation of the density. Only the first moment was given as input to optimization problem (15).

Figure 2 :

 2 Figure 2: Approximation of the support of the invariant measure of the Hénon map, confidence regions of the Christoffel polynomial built from 20 moments : green 90 %, red 50 %, blue simulated trajectories.

x 1 x 2 Figure 3 :

 23 Figure 3: Approximation of the support of the invariant measure of the Hénon map, confidence regions of the Christoffel polynomial built from 20 moments : green 90 %, red 50 %, blue simulated trajectories.

Table 1 :

 1 Logistic map: comparison of moments computed by the SDP relaxation[START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] with the true moments computed analytically. Only the first moment was given as data input to[START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF].

	Moments	x	x 2	x 3	x 4	x 5	x 6	x 7	x 8	x 9	x 10
	SDP	0.0000 0.5000 0.0047 0.3750 0.0062 0.3131 0.0068 0.2746 0.0071 0.2477
	True	0.0000 0.5000 0.0000 0.3750 0.0000 0.3125 0.0000 0.2734 0.0000 0.2461

Table 3 :

 3 OrnsteinUhlenbeck process: comparison of moments computed by the SDP relaxation[START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF]

	with moments computed analytically.
	Moments x x 2 x 3	x 4	x 5	x 6	x 7	x 8	x 9	x 10

Table 4 :

 4 Nonlinear drift: comparison of moments computed by the SDP relaxation (15) with moments computed using a Monte Carlo simulation.

	Moments	x	x 2	x 3	x 4	x 5	x 6	x 7	x 8	x 9	x 10
	SDP	-0.0000 0.4758 0.0000 0.5000 -0.0000 0.7138 0.0000 1.2500 -0.0000 2.4982
	Mont Carlo 0.0001 0.4776 0.0011 0.4997 0.0039 0.7192 0.0130 1.2637 0.0413 2.5683

Table 5 :

 5 Detecting nonexistence of invariant measure: minimum degree required to render the SDP relaxation[START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF] infeasible versus the value of the drift coefficient a.

	drift coefficient a	0 1 2 3 4	5	6	7	8	9 10
	minimum degree for infeasiblity of (15) 2 4 6 8 10 12 14 16 18 20 22

A set of the form (9) is called basic semialgebraic; this class of sets is very rich, including balls, boxes, ellipsoids, discrete sets and various convex and non-convex shapes.

More precisely, F is assumed to be lower semi-continuous with respect to the product topology on the space of sequences R ∞ . This is in particular satisfied if F depends only on finitely many moments as, for example, in[START_REF] Goluskin | Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system[END_REF].

For a relation of this form of a Markov process to the one specified by the transition kernel, see, e.g.,[START_REF] Hernández-Lerma | Discrete-time Markov control processes: basic optimality criteria[END_REF].

Acknowledgment

This work benefited from discussions with Victor Magron. This research was supported in part by the ARO-MURI grant W911NF-17-1-0306. The research of M. Korda was also supported by the Swiss National Science Foundation under grant P2ELP2 165166.

Continuous-time Markov processes

The extension to continuous-time stochastic processes is straightforward. The invariance condition [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in Matlab[END_REF] is simply replaced by X Af dµ = 0, [START_REF] Ozay | Set membership identification of switched linear systems with known number of subsystems[END_REF] where A is the infinite-dimensional generator of the process. For concreteness, let us consider the stochastic differential equation

where b : R n → R n is the drift, σ : R n → R n×k the diffusion matrix and W t is a vector-valued Wiener process. Then we have

which is a polynomial provided that b(x) and σ(x) are polynomial in x and we set f (x) = x α . The invariance condition [START_REF] Ozay | Set membership identification of switched linear systems with known number of subsystems[END_REF] can therefore be expressed solely in terms of the moments of µ (4). The approach then proceeds in exactly the same fashion as described in Section 3.2, leading to a finite-dimensional relaxation of the same form as [START_REF] Junge | On the sighting of unicorns: A variational approach to computing invariant sets in dynamical systems[END_REF].

Certifying non-existence of an invariant measure Contrary to the deterministic case, stochastic processes driven by a Wiener process typically evolve on non-compact domains, thereby rendering the question of the existence of an invariant measure much more subtle. For example, the Wiener process itself does not admit an invariant measure but the Ornstein-Uhlenbeck process dx(t) = -ax(t)dt + σdW (t) does for any a > 0. Interestingly, the proposed approach provides a means to numerically certify that no invariant measure exists. Indeed, since the feasible set of the SDP relaxation (15) contains the truncated moment sequences of all invariant measures, proving the emptiness of this feasible set (which is a finite-dimensional spectrahedron) implies the non-existence of an invariant measure. Indication of such infeasibility is detected during the solution of SDP ( 15) by most existing solvers, although its rigorous certification is more involved [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF]. The proposed method thereby complements methods based on sum-of-squares programming that can be used to prove the existence of an invariant measure using Foster-Lyapunov conditions [START_REF] Meyn | Markov chains and stochastic stability[END_REF].

Eigenmeasures of Perron-Frobenius

In this section we briefly describe how the presented approach can be extended to computation of the eigenmeasures of the Perron-Frobenius operator [START_REF] Mezić | Comparison of systems with complex behavior[END_REF]. For concreteness we work with the discrete-time dynamics (1), although analogous results can be obtained for continuous time, following the developments of Section 6. The Perron-Frobenius operator P : M c (X) → M c (X) is defined by (Pµ)(A) = µ(T -1 (A))