Convex computation of maximal Lyapunov exponents

The upper bound 6 holds for each bounded $x(t)$ and differentiable $V$, so one can maximize the left-hand side over initial conditions $x_{0}$ and minimize the right-hand side over $V$. Let $\mathcal{B}_{0}$ be a set of initial conditions whose trajectories remain bounded, uniformly in $t\geq 0$ but not necessarily uniformly in $x_{0}$. Suppose all of these trajectories eventually remain in a set $\mathcal{B}$, which may be unbounded. Then 6 can be minimized over $x_{0}\in\mathcal{B}_{0}$ and maximized over $V\in\mathcal{C}^{1}(\mathcal{B})$ to obtain

Here $\mathcal{C}^{1}(\mathcal{B})$ is defined as the restriction of $\mathcal{C}^{1}(\mathbb{R}^{n})$ functions to a domain $\mathcal{B}$ (possibly of lower dimension), which is equivalent to the set of functions mapping $\mathcal{B}$ to $\mathbb{R}$ that admit a differentiable extension to an open neighborhood of $\mathcal{B}$. In 7 one can always choose $\mathcal{B}=\mathbb{R}^{n}$ but often it is possible to choose smaller $\mathcal{B}$, and this may decrease the upper bound. If $\mathcal{B}_{0}$ lies in the basin of a local attractor, for instance, one can choose any $\mathcal{B}$ that contains the attractor. If $\mathcal{B}_{0}$ is forward invariant, meaning no trajectories leave the set forward in time, then one can choose $\mathcal{B}=\mathcal{B}_{0}$.

When $\mathcal{B}$ is forward invariant and also compact, and one chooses $\mathcal{B}_{0}=\mathcal{B}$, then 7 is an equality [62],

The equality 8 can be seen as a statement of strong convex duality: the right-hand side can be restated with the inner maximization over invariant measures, and the left-hand side can be stated in the same way with the order of minimization and maximization reversed. For a sketch of the proof of 8, see expression (17) in [62]. The duality result 8 has analogous precedents in the context of discrete-time dynamics, for which we refer to the discussion and references in [6, 35]. In [62] it is assumed that $\mathcal{B}$ is a full-dimensional region, but their proof works more generally for any set $\mathcal{B}$ that is compact and forward invariant [61]. This generality is needed in what follows, where we apply 8 to sets defined by Cartesian products involving the sphere $\mathbb{S}^{n-1}$.

Although the infimum on the right-hand sides of 7 and 8 may not be attained, there exist $V$ giving arbitrarily sharp upper bounds on the left-hand supremum. However, it can be difficult to evaluate the inner supremum on the right-hand side for a given $V$, and optimizing over $V$ is still harder. An important exception, in which computationally tractable relaxations of the right-hand minimax problem are possible, is the case where the ODE right-hand side $f(x)$ and the quantity $\Phi(x)$ are polynomials.

The minimax problem on the right-hand sides of 7 and 8 can be rewritten as a minimization subject to a pointwise inequality constraint,

where

In the case where the given $f$ and $\Phi$ are polynomial in the components of $x$, if $V$ is chosen to be polynomial also, then so is $S$. Even for polynomial $S$ the right-hand minimization in 9 is often intractable; deciding nonnegativity of a multivariate polynomial on a given set has NP-hard computational complexity in general [47]. However, there is a standard way to enforce the nonnegativity constraint on $S$ using SOS conditions that are more tractable.

The simplest way to enforce nonnegativity of $S$ using SOS constraints is to require that $S$ admits an SOS representation, meaning that there exist a finite number of polynomials $q_{k}(x)$ such that $S=\sum_{k=1}^{K}q_{k}^{2}$. In other words, $S$ belongs to the set $\Sigma_{n}$ of $n$-variate SOS polynomials, which is a strict subset (when $n\geq 2$) of the set of globally nonnegative polynomials. In the corresponding SOS relaxation of the right-hand side of 9, the set $\mathcal{C}^{1}(\mathcal{B})$ is replaced by its subset $\mathbb{R}[x]$ of real polynomials in $x$, and the constraint is strengthened to $S\in\Sigma_{n}$. In some applications [24, 26, 38, 50] this simple SOS relaxation has given an upper bound not only equal to the left-hand side of 9 but also to the maximal time average $\sup_{x_{0}\in\mathcal{B}_{0}}\overline{\Phi}$ that one seeks to bound. In general, however, enforcing $S\in\Sigma_{n}$ can strictly increase the right-hand infimum in 9. This would be the case in our application to LEs below, where nonnegativity of a polynomial must be enforced only when one of its vector arguments has unit length. Thus a constraint weaker than $S\in\Sigma_{n}$ is needed.

When $\mathcal{B}\subsetneq\mathbb{R}^{n}$ there is a standard way to formulate SOS constraints that imply nonnegativity on $\mathcal{B}$ but not on all of $\mathbb{R}^{n}$, assuming that $\mathcal{B}$ is a semialgebraic set. That is, we assume $\mathcal{B}$ can be specified by a finite number of polynomial inequalities and/or equalities:

A sufficient condition for $S$ to be nonnegative on $\mathcal{B}$ is that there exist polynomials $\rho_{j}\in\mathbb{R}[x]$ and SOS polynomials $\sigma_{i}\in\Sigma_{n}$ such that

When $x\in\mathcal{B}$, the first sum in 12 is nonnegative and the second vanishes, so $S(x)\geq 0$. When $x\notin\mathcal{B}$, either sign of $S$ is possible. Relaxing the right-hand minimization in 9 by restricting to $V\in\mathbb{R}[x]$ and strengthening its constraint to 12 gives a possibly larger infimum,

It is proved in [38] that 13 is an equality if the semialgebraic specification 11 of $\mathcal{B}$ satisfies a technical condition called the Archimedean property [40, definition 3.4]. Every set with an Archimedean specification is compact, and every compact semialgebraic set’s specification can be made Archimedean, if it isn’t already, by adding a constraint $g_{i}(x)=R^{2}-|x|^{2}\geq 0$ with $R$ large enough so that the constraint is redundant.

If a compact set $\mathcal{B}$ is specified in an Archimedean way so that 13 is an equality, and if $\mathcal{B}$ is forward invariant so that the equality 8 holds, then the maximal time average we seek to bound is equal to the infimum on the right-hand side of 13 [38, theorem 1]:

Theorem 4 in [36] implies an equality similar to 14 but where the infimum subject to SOS constraints is replaced by its convex dual – a maximization subject to constraints on moments of invariant measures.

One more relaxation of the right-hand minimization in 14 is needed to obtain a computationally tractable problem: one chooses finite polynomial spaces from which to seek $V$ and each $\sigma_{i}$ and $\rho_{j}$. This yields an SOS program—a convex optimization problem that has SOS constraints and whose tunable parameters appear only affinely in these constraints and in the optimization objective [55]. (Note that the tunable parameters in 14 are $B$ and the coefficients of $V$, which appear affinely in $S$, along with the coefficients of each $\sigma_{i}$ and $\rho_{j}$.) The resulting SOS program will be computationally tractable when the dimensions of the ODE state space and of the vector spaces of tunable polynomials are not too large. A sequence of SOS programs with increasing computational cost can be defined by choosing successively higher-dimensional spaces over which to optimize the tunable polynomials. The infima of these SOS programs converge from above to the right-hand side of 14. In cases where the equality 14 holds, this means that SOS programs can give arbitrarily sharp upper bounds on the left-hand maximal time average.

To bound the time average on the right-hand side of 16 using the framework of section 2.1, we need evolution equations for the $(x,z)$ state space. The ODE for $z$ is found by taking the time derivative of the definition $z(t)=y(t)/|y(t)|$ and applying the $y$ ODE 2. Coupling this with the $x$ ODE 1 gives

where

Recall that the framework of section 2.1 was described for an ODE with state space $\mathbb{R}^{n}$, right-hand side $f(x)$, and initial set $\mathcal{B}_{0}$ whose trajectories eventually remain in $\mathcal{B}$. To bound the maximal LE of the $x$ ODE, we apply the same framework to an ODE with state space $\mathbb{R}^{n}\times\mathbb{S}^{n-1}$, right-hand side as in 17, and initial set $\mathcal{B}_{0}\times\mathbb{S}^{n-1}$ whose trajectories eventually remain in $\mathcal{B}\times\mathbb{S}^{n-1}$. We choose $\Phi(x,z)=z^{\sf T}Df(x)z$, so that $\mu=\overline{\Phi}$. Then the left-hand side of 7 becomes the maximal LE $\mu^{*}_{\mathcal{B}_{0}}$ defined in 4, and the right-hand side of 7 becomes the upper bound that we seek. This upper bound on $\mu^{*}_{\mathcal{B}_{0}}$ is stated in the first part of proposition 1 below, which requires no further proof. When $\mathcal{B}$ is compact and forward invariant for the $x$ ODE, so that $\mathcal{B}\times\mathbb{S}^{n-1}$ is compact and forward invariant for 17, the equality 8 proved in [62] can be applied in the present context to give the second part of the proposition. The assumption $f\in\mathcal{C}^{2}$ may be stronger than needed but suffices to ensure $(f,\ell)\in\mathcal{C}^{1}$, in which case the result of [62] is directly applicable to the ODE with right-hand side $(f,\ell)$.

In the proposition, $\nabla_{x}V$ and $\nabla_{z}V$ denote gradients of the auxiliary function $V(x,z)$ with respect to $x$ and $z$, respectively. In the first part of the proposition, neither $\mathcal{B}_{0}$ nor $\mathcal{B}$ is necessarily bounded—for instance, 19 holds with $\mathcal{B}_{0}=\mathcal{B}=\mathbb{R}^{n}$ if each trajectory $x(t)$ is bounded forward in time. Note that the right-hand sides of 19 and 20 are identical, but 20 requires stronger assumptions on $\mathcal{B}$ along with choosing $\mathcal{B}_{0}=\mathcal{B}$ on the left-hand side. A result analogous to 20 in the context of linear systems with stochastic forcing is proved in [37].

If the $x$ ODE 1 is polynomial, then so are the $(x,z)$ ODE 17 and the quantity $\Phi=z^{\sf T}Df(x)z$, and we let $V(x,z)$ be polynomial also. The inner supremum on the right-hand sides of 19 and 20 is bounded above by a constant $B$ if and only if the polynomial

is nonnegative for all $(x,z)\in\mathcal{B}\times\mathbb{S}^{n-1}$. Note that the $P(x,z)\geq 0$ constraint in the context of the $(x,z)$ ODE 17 is the analogue of the $S(x)\geq 0$ constraint appearing in 9 in the context of the $x$ ODE alone. Just as $S(x)\geq 0$ on $\mathcal{B}$ if the SOS conditions in 12 hold, $P\geq 0$ on $\mathcal{B}\times\mathbb{S}^{n-1}$ if there exist polynomials $\rho_{j}(x,z)$ and SOS polynomials $\sigma_{i}(x,z)$ such that

where $h_{0}(z)=1-|z|^{2}$ so that it vanishes on $\mathbb{S}^{n-1}$, and the conditions $g_{i}(x)\geq 0$ and $h_{j}(x)=0$ for $i,j\geq 1$ specify $\mathcal{B}$ as in 11. Note that we allow non-compact $\mathcal{B}$, including the case $\mathcal{B}=\mathbb{R}^{n}$ wherein 22 should be interpreted as the single SOS constraint $P-\rho_{0}h_{0}\in\Sigma_{2n}$.

The SOS relaxation of the right-hand minimization in 19 and 20 is the minimization of the constant $B$ that appears in $P(x,z)$, where the polynomials $V,\sigma_{i},\rho_{j}$ are optimized subject to 22. The infimum of the relaxed problem is an upper bound on the right-hand infimum in 19, which gives the first part of proposition 2 below with no further proof needed. If the specification of $\mathcal{B}$ is Archimedean, then adding the constraint $h_{0}(z)=0$ gives an Archimedean specification of $\mathcal{B}\times\mathbb{S}^{n-1}$. In this case the SOS relaxation does not change the value of the right-hand infimum in 19 and 20, as explained after 13. The second part of proposition 2 thus follows from the second part of proposition 1.

The right-hand sides of 23 and 24 are identical, but 24 requires stronger assumptions on $\mathcal{B}$ along with choosing $\mathcal{B}_{0}=\mathcal{B}$ on the left-hand side. The minimization on the right-hand sides of 23 and 24 is over $V(x,z)$, $\rho_{j}(x,z)$, and $\sigma_{i}(x,z)$ in the infinite-dimensional space $\mathbb{R}[x,z]$ of $2n$-variate polynomials. Restricting each of these tunable polynomials to a finite-dimensional subspace—i.e., choosing a finite basis with tunable coefficients—gives an SOS program that can be solved computationally when $n$, $I$, $J$, and the polynomial subspace dimensions are not too large. As required for an SOS program, all tunable coefficients appear affinely in the SOS expressions since $B$ and $V$ appear affinely in the definition 21 of $P$. Although restricting the tunable polynomials $V,\sigma_{i},\rho_{j}$ to finite subspaces may increase the infimum on the right-hand sides of 23 and 24, one can always enlarge these subspaces and solve another SOS program with greater computational cost. In this way, a sequence of increasingly expensive SOS programs will give a nonincreasing sequence of infima that converge to the right-hand infimum in 23 and 24. Thus, when the equality 24 holds, SOS computations can give arbitrarily sharp upper bounds on the maximal LE $\mu^{*}_{\mathcal{B}}$, at least in theory. The computational examples presented below in section 3 show that convergence of upper bounds to the maximal LE indeed can be achieved in practice.

An ODE is said to be symmetric under a transformation if solutions are mapped to solutions. Let $\Lambda\in GL(n)$, where $GL(n)$ denotes the group of invertible linear transformations on $\mathbb{R}^{n}$. An ODE $\frac{\mathrm{d}}{\mathrm{d}t}x=f(x)$ on $\mathbb{R}^{n}$ is symmetric under $\Lambda$ if and only if $f$ is equivariant under $\Lambda$, meaning $f(\Lambda x)=\Lambda f(x)$ for all $x\in\mathbb{R}^{n}$. A function $\Phi$ is invariant under $\Lambda$ if $\Phi(\Lambda x)=\Phi(x)$ for all $x\in\mathbb{R}^{n}$. The invariant transformations of any function form a group, as do its equivariant transformations. A function is said to be $\mathcal{G}$-invariant (resp. $\mathcal{G}$-equivariant) for a group $\mathcal{G}$ if it is invariant (resp. equivariant) under all $\Lambda\in\mathcal{G}$. Turning particularly to the $(x,z)$ dynamics 17, we seek a symmetry group for which the right-hand side of the ODE is equivariant, and the function $\Phi=z^{\sf T}Df(x)z$ is invariant, under all transformations in the group.

The transformation $(x,z)\mapsto(x,-z)$ is a symmetry both of the $(x,z)$ ODE 17 and of the function $\Phi(x,z)=z^{\sf T}Df(x)z$, reflecting the fact that the growth or decay rate of a vector in the tangent space is unaffected by reversing its direction. The equivariance of the ODE right-hand side under this transformation amounts to the oddness $\ell(x,-z)=-\ell(x,z)$, where $\ell(x,z)$ is defined as in 18, and the invariance $\Phi(x,-z)=\Phi(x,z)$ is clear.

Further symmetries of the $(x,z)$ dynamics are guaranteed only if the $x$ ODE alone is symmetric under orthogonal transformations, common examples of which include sign changes or rotations. In particular, if the $x$ dynamics are symmetric under $\Lambda\in O(n)$, where $O(n)$ denotes the group of orthogonal linear transformations on $\mathbb{R}^{n}$, then the $(x,z)$ dynamics are symmetric under $(x,z)\mapsto(\Lambda x,\Lambda z)$, and so is the function $\Phi(x,z)$. To state this precisely in the following proposition, for each group $\mathcal{G}\subset O(n)$ we define the corresponding group $\mathcal{G}^{\prime}\in O(2n)$ as

If an ODE, its domain, and a time-averaged quantity share a compact linear group of symmetries, then the upper bound on this time average given by the right-hand infimum in 7 and 8 is unchanged if the auxiliary function $V$ is further constrained to be invariant under the same symmetries. In particular, any $V$ that yields an upper bound on a time average can be used to construct a symmetrized $V$ yielding the same bound. This is proved in [26, Proposition A.1] in the context of a finite cyclic symmetry group. Proposition 6 in the appendix generalizes the result to all compact groups of linear transformations—i.e., all groups of orthogonal transformations, along with similarity transformations 28 of such groups.

We apply these symmetry results in the context of bounding LEs by the infimum on the right-hand sides of 19 and 20, in which case the pertinent symmetry group $\mathcal{G}^{\prime}$ of the $(x,z)$ dynamics is given by proposition 3 above. Provided that the domain $\mathcal{B}$ is $\mathcal{G}$-invariant, the following proposition guarantees that corresponding symmetries can be imposed on $V(x,z)$ in the right-hand minimization of 19 and 20 without changing the resulting upper bound.

Proposition 4, which allows invariance to be imposed on $V$ in the minimization that appears on the right-hand sides of 19 and 20, has an analogue for the minimization’s SOS relaxation that appears on the right-hand sides of 23 and 24. This is stated by proposition 5 below, whose assumptions require that the $\mathcal{G}$-invariant set $\mathcal{B}$ is specified in terms of $\mathcal{G}$-invariant polynomials, and whose conclusions allow invariance to be imposed not only on $V$ but also on any other tunable polynomials $\sigma_{i}$ and $\rho_{j}$. We do not include a proof because the result is a direct application of proposition 7 in the same way that proposition 4 is an application of proposition 6.