Time average bounds via Sum of Squares

9 February 2021 | R. Rosa

Motivation

In many real-world problems, one is interested in estimating certain key quantities related to the problem. For instance, in fluid flows, quantities of interest involve kinetic energy, enstrophy, drag coefficient, energy dissipation rate, and so on. In other applications, one might be interested in mechanical stress, chemical concentration, infected population, pharmaceutical dosage, etc.

Many such problems can be resonably modeled by differential equations, which may, however, exibit complicate, perhaps chaotic dynamics. In those situations, the instantaneous value of certain quantities vary unpredictably in time, but very often their mean value is reasonably steady.

This mean value can be considered in different ways, e.g. as time average, as ensemble average, or as spatial average, and are thus more ameanable to analysis. This article considers ways to estimate time and ensemble averages of certain quantities.

Problem setting

If a model for the problem is available in the form of an ordinary differential equation

dudt=F(u), \frac{\textrm{d} u}{\textrm{d}t} = F(u),

where F:XXF:X\rightarrow X is some locally Lipschitz function acting in some finite-dimensional space X=RnX=\mathbb{R}^n, nNn\in\mathbb{N}, then, for each u0Xu_0\in X, there exists a unique solution u=u(t)u=u(t) with u(0)=u0u(0)=u_0. If we assume all solutions are defined globally in the futures, we obtain a continuous semigroup {S(t)}t0\{S(t)\}_{t\geq 0} acting on XX, with given by S(t)u0=u(t)S(t)u_0=u(t).

Given a function ϕ:XR\phi:X\rightarrow \mathbb{R} representing some "real" quantity we want to measure, the asymptotic superior limit of the time-average of ϕ(u(t))\phi(u(t)) is given (and denoted) by

ϕˉ(u0)=lim supT1T0Tϕ(u(t)) dt. \bar\phi(u_0) = \limsup_{T\rightarrow} \frac{1}{T} \int_0^T \phi(u(t))\;\textrm{d}t.

The idea here is that we would like to find an upper bound for ϕˉ(u0)\bar\phi(u_0), for all possible initial conditions u0Bu_0\in B, in some subset of interest BXB\subset X.

Assuming that ϕ\phi is continuous, that BB is positively invariant (i.e. u(0)Bu(0)\in B implies u(t)Bu(t)\in B, t0\forall t\geq 0.), and BB is compact, then tϕ(u(t))t\mapsto \phi(u(t)) is bounded on t0t\geq 0, and the superior limit above is uniformily bounded in u0Bu_0\in B.

The problem, then, is to find the best possible bound C for supu0Bϕˉ(u0)\sup_{u_0\in B}\bar\phi(u_0).

Overview of some analytic and numerical methods

One way to bound ϕ(u0)\phi(u_0) is to do what is called energy-type estimates, which amounts to multiplying the equation by appropriate terms aiming to obtain inequalities that eventually lead to an estimate for ϕ(u0)\phi(u_0). Or by variational estimates, introducing an auxiliary function, within a special class of functions, and performing some minimization with respect to the auxiliary function.

Numerically, we can use a Monte-Carlo method and simulate the evolution of the equation in a computer, with randomly chosen initial conditions, and look for the time average over sufficiently long time intervals (sufficiently in the sense, for instance, that the finite-time average does not change much – according to a given error – by increasing the averaging time). Or one can vary the auxiliary function in the special class of functions and look for the best estimate.

More recently, however, a novel method is being used, which is a sort of variational technique, but in a different perspective and leading to an efficient numerical approach. There are some aspects I would like to discuss concerning this method:

  1. It can be seen as a convex minimization problem;

  2. When FF and ϕ\phi are polynomials, the minimization problem can be relaxed to a P-complete problem by looking for a Sum of Squares (SoS) representation of an appropriate term, at the cost of obtaining a larger bound, but which is often near the optimal value;

  3. The original convex minimization problem can be recast into a minimax problem and be showed to indeed yield the optimal estimate;

  4. The optimality result above has been proved first in the finite-dimensional case and, more recently, myself and a co-author extended it to dissipative evolutionary partial differential equations.

The Sum of Squares (SoS) problem

Looking for an expression for a nonnegative multivariate real polynomial p=p(x)p=p(x) as a Sum of Squares (SoS) of other polynomials, i.e. p(x)=j=1kpj(x)2p(x) = \sum_{j=1}^k p_j(x)^2 for other polynomials pj=pj(x)p_j=p_j(x), is not a new problem. In 1885, the 21-year old Minkowski made his inaugural dissertation on quadratic polynomials conjecturing that there must exist homogeneous, nonnegative real polynomials of degree mm in nn variables, for arbitrary m,n>2m, n > 2, which are not sums of squares of other homogeneous real polynomials. Hilbert initially attacked Minkowski's claim, but by the end of the presentation Hilbert was convinced that this might indeed be true at least starting with n=3n=3. Three years later, at the age of 26, Hilbert himself proved that the claim is not true for n=3,m=4n=3, m=4, but that for n3,m6n\geq 3, m\geq 6 or for n4,m4n\geq 4, m\geq 4, the set of nonnegative polynomials of degree mm in nn variables is indeed strictly larger than the set of sum of squares of polynomials.

Further work on the subject led him to formulate the 17th problem in his list of 23 problems presented in 1900: must every nonnegative homogenous polynomial be expressed as a sum of squares of rational functions?

Hilbert used tools from classical algebraic geometry at that time, without given explicit examples for the problem addressing Minkowski's dissertation. Explicit examples of homogenous polynomials which are not sum of squares of other polynomials were only given in the second half of the 20th century. One famous example is that of Motzkin:

f(x,y,z)=z6+x4y2+x2y43x2y2z2.{\displaystyle f(x,y,z)=z^{6}+x^{4}y^{2}+x^{2}y^{4}-3x^{2}y^{2}z^{2}.}

Hilbert's 17th problem was solved in the affimartive by Artin, in 1926. For further historical accounts related to Hilbert's 17th problem, see e.g. Reznick (2000).

More recently, a number of methods to actually test and find whether a given multivariate nonnegative polynomial is a sum of squares of polynomials have been devised (e.g. Shor 1980s, 1990s, Choi, Lam, Reznik 1990s). Then, Parrilo (2000) presented in his PhD thesis, and in subsequent articles (e.g. Parrilo (2003)), several applications to differential equations, including the search for Lyapunov functions and control strategies. By the early 2000s, a number of MATLAB toolbox solvers were already available.

Applications to local stability of PDEs and, in particular to 2D fluid flows were given, respectively by Papachristodoulou and Peet (2006) and Yu, Kashima, Imura (2008).

Finally we get to the articles related to the main motivation for this pots, which is that of globally estimating quantities related to the problem at hand and the global stability of the model.

The first article which seems to exploit the technique of Sum of Squares to the global analysis of PDEs seem to be that of Goulart and Chernyshenko (2012), which considered, in particular, the global stability of fluid flows. This was soon followed by a number of works by various authors: Fantuzzi, Goluskin, Doering, Goulart, Chernyshenko, Huang, Papachristodoulou (2010s), among others (see e.g. Chernyshenko, Goulart, Huang, and Papachristodoulou). These culminated with the work of Tobasco, Goluskin, and Doering (2018) showing that the convex optimization problem can be written as a minimax problem, which can then be proved to yield an optimal result for the estimate of the asymptotic time averages mentioned earlier. In turn, this gives the expectation that relaxing the problem to use the sum of squares approach yields sharp bounds, close to the optimal one.

The convex minimization problem

Now we begin to directly address the points raised above. Let us go back to the setting described earlier and see how the convex optimization problem appears.

One key aspect is to realize that, given any continuously differentiable function V:XRV:X\rightarrow \mathbb{R}, it follows, by the chain rule and integration, that the asymptotic time average of F(u)V(u)F(u)\cdot \nabla V(u) satisfies

0TF(u)V(u) dt=0Tu˙V(u) dt=0TddtV(F(u)) dt=V(u(T))V(u(0)). \begin{aligned} \int_0^T F(u)\cdot \nabla V(u) \;\textrm{d}t & = \int_0^T \dot u\cdot \nabla V(u) \;\textrm{d}t \\ & = \int_0^T \frac{\textrm{d}}{\textrm{d}t} V(F(u)) \;\textrm{d}t \\ & = V(u(T)) - V(u(0)). \end{aligned}

Since BB is assumed to be compact and positively invariant, then V(u(t))V(u(t)) is uniformly bounded in t0t\geq 0, so that

1T0TF(u)V(u)=V(u(T))V(u(0))T0, \frac{1}{T} \int_0^T F(u)\cdot \nabla V(u) = \frac{V(u(T)) - V(u(0))}{T} \rightarrow 0,

as TT\rightarrow \infty. Hence, using that a bar denotes limit superior of the time-averages, we write

FV=0. \overline{F\cdot\nabla V} = 0.

Since (6) actually holds for the limit itself, not only the superior limit, then we may add it to ϕˉ\bar\phi and have that

ϕˉCϕ+FVC, \bar\phi \leq C \Leftrightarrow \overline{\phi + F\cdot\nabla V} \leq C,

on BB, for arbitrary continuously differentiable function VV.

Since the above holds for arbitrary such VV and since the aim is to obtain the best possible CC, in the sense of being the smallest possible one, we can write the problem of finding such bound CC for ϕˉ\bar\phi as the minimization problem

supϕˉinf(C,V)R×C1, ϕ+(FV)CC, \sup \bar\phi \leq \inf_{(C,V)\in \mathbb{R}\times\mathcal{C}^1, \;\overline{\phi + (F\cdot \nabla V)} \leq C} C,

However, if we had to check whether the time averages ϕ+(FV)\overline{\phi + (F\cdot \nabla V)} are smaller than or equal to CC, for every initial condition u0u_0 in BB, we would actually have much more work than simply checking whether ϕˉC\bar\phi \leq C. So, the idea is to require the much stronger condition that ϕ(u)+F(u)V(u)C\phi(u) + F(u)\cdot \nabla V(u) \leq C, for every point uu in BB. Notice this is not a dynamic condition. We are not solving any differential equation. The point uu is an arbitrary point in BB. It is a pointwise bound, that certainly implies the time-average bound for any solution tu(t)t\mapsto u(t) starting at the positively invariant set BB.

It may seem, at first, that, by requiring this stronger condition, we end up with a much worse bound. However, it turns out that the minimization process somehow compensates for that and end up yielding an optimal bound just like we would obtain by requiring only that the time-average be smaller than or equal to CC. This magic is taken care of by the inclusion of the auxiliary function VV, which is sometimes called the reservoir function. Notice that the time-average FV\overline{F\cdot \nabla V} vanishes, but when considering FVF\cdot \nabla V for arbitrary points in BB, this term, for suitable VV, can be negative to compensate when ϕ\phi is large, and it is allowed to be positive, when ϕ\phi is small, such that at the end we find a relatively small bound CC.

Notice we don't expect FVF\cdot \nabla V to be negative all the time, otherwise VV would be like a Lyapunov function, or a La Salle-type function, and the solutions would converge to the invariant set included in the set {V=0}\{V=0\}. Some systems do have such a function, but this is not expected to exist in more complicate problems.

Now, by requiring that ϕ+FVC\phi + F\cdot \nabla V\leq C holds pointwise in all BB, instead of only along the time average of the trajectories u(t)u(t), we arrive at the following minimization problem:

supϕˉinf(C,V)R×C1, ϕ+(FV)CC. \sup \bar\phi \leq \inf_{(C,V)\in \mathbb{R}\times\mathcal{C}^1, \;\phi + (F\cdot \nabla V) \leq C} C.

We may rewrite this as

supϕˉinf(C,V)R×C1, SC,V(u)0,uBC, \sup \bar\phi \leq \inf_{(C,V)\in \mathbb{R}\times\mathcal{C}^1, \;S_{C,V}(u)\geq 0, \forall u\in B} C,

where SC,V(u)=Cϕ(FV)S_{C,V}(u) = C - \phi - (F\cdot \nabla V). This is a convex minimization problem, since the objetive function (C,V)C(C,V)\mapsto C is linear, and the minimization is sought after within the set SC,V0S_{C,V}\geq 0, which is convex since (C,V)SC,V(u)(C,V)\mapsto S_{C,V}(u) is linear and the half plane is convex.

Relaxing the minimization problem to a Sum of Squares minimization

The minimization problem (10) can be NP-hard to compute. However, when the differential equation term FF and the quantity of interest ϕ\phi are polynomials, the minimization problem can be relaxed to a P-complete convex minimization problem by restricting VV to be a polynomial of a given order, or some special type of polynomial, and requiring that the polynomial SC,VS_{C,V} be SoS, which certainly implies the condition that it be nonnegative. That might not yield an optimal bound, but it's been show to yield pretty sharp estimates for a number of equations.

This formulation takes the precise form

supϕˉinf(C,V)R×Pm(X), SC,V(u)=SoSC, \sup \bar\phi \leq \inf_{(C,V)\in \mathbb{R}\times\mathrm{P}_m(X), \;S_{C,V}(u) = \texttt{SoS}} C,

where above we denote by Pm(X)\mathrm{P}_m(X) the set of real polynomials on XX.

This problem can be regarded as a semidefinite programming. It is similar to linear programming, but in which the first orthant xi0x_i\geq 0 is replace by the cone of positive semidefinite matrices S0S\succeq 0. More precisely, we may start with the primal problem:

Minimize LS subject to AiS=bi and S0, \text{Minimize } L\cdot S \text{ subject to } A_i\cdot S = b_i \text{ and } S \succeq 0,

where bRMb\in \mathbb{R}^M is a given vector, LL and AiA_i are given symmetric N×NN\times N matrices, and SS is the decision variable, also assume to be symmetric. The dot product for matrices is element-wise, i.e. AiA=(Ai)jkAjk=0A_i\cdot A = (A_i)_{jk} A_{jk} = 0, and S0S\succeq 0 means that SS is positive semidefinite, i.e. ξSξ0\xi\cdot S \xi \geq 0, for every ξRN\xi\in \mathbb{R}^N.

The minization problem above has the dual formulation

Maximize bη subject to i=1MηiAiL, \text{Maximize } b\cdot \eta \text{ subject to } \sum_{i=1}^M \eta_iA_i\preceq L,

where ηRM\eta\in\mathbb{R}^M. Any solution of the dual problem is a lower bound for the primal problem, and, conversely, any solution of the primal problem yields an upper bound for the dual problem. In fact, this follows from

LSbη=LSi=1MηiAiS=(Li=1MηiAi)S0. L\cdot S - b\cdot \eta = L\cdot S - \sum_{i=1}^M \eta_iA_i\cdot S = (L - \sum_{i=1}^M \eta_iA_i)\cdot S \geq 0.

Thus,

bηLS b\cdot \eta \leq L\cdot S

for any feasible η\eta and SS in each problem.

The question now is how to frame the Sum of Squares problem into a semidefinite programming one. As described in Parrilo (2003), it is possible to write the sum of squares problem in either the primal form or the dual form. In theory, they are mathematically equivalent, but one formulation may be numerically more efficient than the other, depending on the dimension of the problem. For the sake of illustration, we describe below how to arrive at the primal problem.

So, suppose a multivariate real polynomial p(x)p(x), xRnx\in \mathbb{R}^n, of degree mm is given. It is easy to argue that, for p(x)p(x) to have any chance of being a sum of squares, or just nonnegative, the degree mm of pp has to be even, say m=2dm=2d. It is also not difficult to argue that it can be written in the form

p(x)=m(x)Sm(x), p(x) = m(x)\cdot Sm(x),

for a symmetric matrix S=(Sjk)S=(S_{jk}), where

ξ=ξ(x)=(1,x1,,xn,x1x1,,x1d,x1d1x2,,x1xnd1,xnd) \xi = \xi(x)=(1, x_1, \ldots, x_n, x_1x_1, \ldots, x_1^d, x_1^{d-1}x_2, \ldots, x_1x_n^{d-1}, x_n^d)

is the vector of all monomials in xx of degree up to d=m/2d=m/2. The dimension of the space for ξ(x)\xi(x) is N=(n+dd)N=\left(\begin{matrix} n + d \\ d \end{matrix}\right).

For example, consider the polynomial

p(x,y)=x44x3y6x2y24xy3+y4, p(x,y)= x^4 - 4x^3y - 6x^2y^2 -4xy^3 + y^4,

in (x,y)R2(x,y)\in \mathbb{R}^2, which we know is SoS since it is precisely (xy)4(x-y)^4. Then, with ξ(x,y)=(1,x,y,x2,xy,y2)\xi(x,y)=(1, x, y, x^2, xy, y^2), we can take

S=[000000000000000000000120000262000021]. S = \left[ \begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -2 & 0 \\ 0 & 0 & 0 & -2 & -6 & -2 \\ 0 & 0 & 0 & 0 & -2 & 1 \end{matrix} \right].

Since the elements of mm are not algebraically independent (e.g. x2y2=(x2)(y2)=(xy)(xy))x^2y^2 = (x^2)(y^2) = (xy)(xy)), such SS is usually not unique. For instance, we can also take

S=[000000000000000000000123000202000321], S = \left[ \begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -2 & -3 \\ 0 & 0 & 0 & -2 & 0 & -2 \\ 0 & 0 & 0 & -3 & -2 & 1 \end{matrix} \right],

or any convex combinations of the two.

Back to the general case (16), if there exists a symmetric matrix SS which is positive semidefinite, then it can be diagonalizable with the elements diid_{ii} in the diagonal being all non-negative, i.e.

ξSξ=ξ(Q1SQ)ξ=(Qξ)D(Qξ)=ζDζ=diiζi2=i(diiζi)2, \xi\cdot S\xi = \xi\cdot (Q^{-1}S Q)\xi = (Q\xi)\cdot D(Q\xi) = \zeta \cdot D \zeta = \sum d_{ii}\zeta_i^2 = \sum_i (\sqrt{d_{ii}}\zeta_i)^2,

where ζ(x)=Qξ(x)\zeta(x) = Q\xi(x) is a vector of polynomials in xx. This yields that pp is a SoS.

Hence, the problem becomes to find a symmetric positive semidefinite matrix SS satisfying (16). The polynomials p(x)p(x) and ξ(x)Sξ(x)\xi(x)\cdot S\xi(x) are equal if, and only if, their coefficients are equal, which is a linear problem for SS, with dimension M=N2=(n+dd)×(n+dd)M=N^2 = \left(\begin{matrix} n + d \\ d \end{matrix}\right) \times \left(\begin{matrix} n + d \\ d \end{matrix}\right). If we define the coeficients of p=p(x)p=p(x) by bib_i and those of ξSξ\xi\cdot S\xi by AiSA_i\cdot S, i=1,,Mi=1,\ldots, M, then the problem becomes to find a symmetrix matrix SS such that

S0,AiS=bi. S \succeq 0, \qquad A_i S = b_i.

If we further ask SS to minimize the quantity LSL\cdot S for some desirable symmetric matrix LL, then we end up with the primal semidefinite programming problem for SS.

The minimax formulation

The convex minimization problem (10) can easily be rewritten in the minimax form

supu0Bϕˉ(u0)minVC1(B)maxuB{ϕ(u)+F(u)V(u)}. \sup_{u_0\in B} \bar\phi(u_0) \leq \min_{V\in \mathcal{C}^1(B)} \max_{u\in B} \left\{\phi(u) + F(u) \cdot \nabla V(u)\right\}.

With this formulation in mind, Tobasco, Goluskin and Doering (2018) gave a beautiful proof that the bound is actually optimal, and that the supremum at the left hand side above is achieved!:

maxu0Bϕˉ(u0)=minVC1(B)maxuB{ϕ(u)+F(u)V(u)}. \max_{u_0\in B} \bar\phi(u_0) = \min_{V\in \mathcal{C}^1(B)} \max_{u\in B} \left\{\phi(u) + F(u) \cdot \nabla V(u)\right\}.

The proof uses Ergodic Theory and a minimax principle. In a future opportunity we will go through its proof, as well as to detail the extension done to the infinite-dimensional setting, which is briefly discussed next.

Extension to the infinite-dimensional setting

The proof in the finite dimensional case uses a few conditions that are delicate to extend to the infinite dimensional case:

  1. The positively invariant set BB has to be compact;

  2. The quantity of interest ϕ\phi has to be a continuous function on the phase space XX;

  3. Borel probability mesure are Lagrangian invariant if and only if they are Eulerian invariant.

By Lagrangian invariant we mean the classical invariant condition μ(S(t)1E)=μ(S(t)E)\mu(S(t)^{-1}E) = \mu(S(t)E), for any Borel set EXE\subset X, where μ\mu is the Borel probability measure in question and {S(t)}t0\{S(t)\}_{t\geq 0} is the semigroup generated by the equation. By Eulerian invariant we mean that μ\mu has to satisfy XF(u)V(u) dμ(u)=0\int_X F(u)\cdot \nabla V(u) \;\textrm{d}\mu(u) = 0, for all VC1(X)V\in \mathcal{C}^1(X).

The assumption that XX be finite-dimensional is not a requirement per se, but it makes the above conditions hold in more generality. For instance, it suffices to have BB closed and bounded to have it compact. And this compactness is needed both for the passage from time average to ensemble average (i.e. average with respect to the invariant measure) and for the minimax principle.

Concerning the assumption of continuity of ϕ\phi, it is not a big deal in finite dimensions, but it is quite restrictive for partial differential equations. For instance, if the phase space is L2(Ω)L^2(\Omega), one cannot consider ϕ(u)\phi(u) involving derivatives of uu. Even if we attempt to use extensions of the minimax principle, they require upper-semicontinuity of ϕ\phi, so even quantities like ϕ(u)=jxju2\phi(u) = \sum_j |\partial_{x_j}u|^2 would not work as is.

But at least for the case of a continuous quantity ϕ\phi in the infinite-dimensional (e.g. kinetic energy on L2L^2), one can go around the requirement of BB being compact by considering dissipative systems which possess a compact attracting set.

The remaining delicate condition is the equivalence between Lagrangian and Eulerian invariance, which is by no means trivial in the infinite-dimensional case. In fact, I know of only two equations for which this has been proved: the two-dimensional Navier-Stokes equations and a globally modified Navier-Stokes equations obtained by truncating the nonlinear term. However, it is my belief that the key tool is simply that it be possible to approximate the system (any solution) by a right-invertible semigroup (e.g. Galerkin approximation or a hyperbolic/wave-type approximation) and exploit the usual a~priori estimates. It is an open field to prove this for other systems or to come up with an easily-applicable general statement.

It should be said that even the notion of Eulerian invariance needs to be relaxed to working for special types of functions VV, which we call cylindrical test functionals. They are at the core of the notion of statistical solution.

As in the finite-dimensional case, we leave further details about the result in infinite dimensions to a future post. Meanwhile, the details can be found in Rosa and Temam (arxiv 2020)

Selected References: