Search

Transformers are extremely successful machine learning models whose mathematical properties remain poorly understood. Here, we rigorously characterize the behavior of transformers with hardmax self-attention and normalization sublayers as the number of layers tends to infinity. By viewing such transformers as discrete-time dynamical systems describing the evolution of points in a Euclidean space, and thanks to a geometric interpretation of the self-attention mechanism based on hyperplane separation, we show that the transformer inputs asymptotically converge to a clustered equilibrium determined by special points called leaders. We then leverage this theoretical understanding to solve sentiment analysis problems from language processing using a fully interpretable transformer model, which effectively captures `context' by clustering meaningless words around leader words carrying the most meaning. Finally, we outline remaining challenges to bridge the gap between the mathematical analysis of transformers and their real-life implementation., Comment: 23 pages, 10 figures, 1 table. Funded by the European Union (Horizon Europe MSCA project ModConFlex, grant number 101073558). Accompanying code available at: https://github.com/DCN-FAU-AvH/clustering-hardmax-transformers

We consider an internally heated fluid between parallel plates with fixed thermal fluxes. For a large class of heat sources that vary in the direction of gravity, we prove that $\langle\delta T \rangle_h \geq \sigma R^{-1/3} - \mu$, where $\langle\delta T \rangle_h$ is the average temperature difference between the bottom and top plates, $R$ is a `flux' Rayleigh number and the constants $\sigma,\mu >0$ depend on the geometric properties of the internal heating. This result implies that mean downward conduction (for which $\langle\delta T \rangle_h< 0$) is impossible for a range of Rayleigh numbers smaller than a critical value $R_0$. The bound demonstrates that $R_0$ depends on the heating distribution and can be made arbitrarily large by concentrating the heating near the bottom plate. However, for any given fixed heating profile of the class we consider, the corresponding value of $R_0$ is always finite. This points to a fundamental difference between internally heated convection and its limiting case of Rayleigh-B\'enard convection with fixed flux boundary conditions, for which $\langle\delta T \rangle_h$ is known to be positive for all $R$., Comment: 13 pages, 4 figures

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

We consider finite element approximations to the optimal constant for the Hardy inequality with exponent $p=2$ in bounded domains of dimension $n=1$ or $n \geq 3$. For finite element spaces of piecewise linear and continuous functions on a mesh of size $h$, we prove that the approximate Hardy constant converges to the optimal Hardy constant at a rate proportional to $1/| \log h |^2$. This result holds in dimension $n=1$, in any dimension $n \geq 3$ if the domain is the unit ball and the finite element discretization exploits the rotational symmetry of the problem, and in dimension $n=3$ for general finite element discretizations of the unit ball. In the first two cases, our estimates show excellent quantitative agreement with values of the discrete Hardy constant obtained computationally., Comment: Review: Significantly improved estimates compared to the original version (23 pages, 6 figures)

We describe a `discretize-then-relax' strategy to globally minimize integral functionals over functions $u$ in a Sobolev space subject to Dirichlet boundary conditions. The strategy applies whenever the integral functional depends polynomially on $u$ and its derivatives, even if it is nonconvex. The `discretize' step uses a bounded finite element scheme to approximate the integral minimization problem with a convergent hierarchy of polynomial optimization problems over a compact feasible set, indexed by the decreasing size $h$ of the finite element mesh. The `relax' step employs sparse moment-sum-of-squares relaxations to approximate each polynomial optimization problem with a hierarchy of convex semidefinite programs, indexed by an increasing relaxation order $\omega$. We prove that, as $\omega\to\infty$ and $h\to 0$, solutions of such semidefinite programs provide approximate minimizers that converge in a suitable sense (including in certain $L^p$ norms) to the global minimizer of the original integral functional if it is unique. We also report computational experiments showing that our numerical strategy works well even when technical conditions required by our theoretical analysis are not satisfied., Comment: Improved main results, added section on computational complexity. 27 pages, 9 figures

We present a flexible data-driven method for dynamical system analysis that does not require explicit model discovery. The method is rooted in well-established techniques for approximating the Koopman operator from data and is implemented as a semidefinite program that can be solved numerically. Furthermore, the method is agnostic of whether data is generated through a deterministic or stochastic process, so its implementation requires no prior adjustments by the user to accommodate these different scenarios. Rigorous convergence results justify the applicability of the method, while also extending and uniting similar results from across the literature. Examples on discovering Lyapunov functions, performing ergodic optimization, and bounding extrema over attractors for both deterministic and stochastic dynamics exemplify these convergence results and demonstrate the performance of the method., Comment: We have significantly expanded the presentation. This has improved the presentation and made the paper more readable. Comments welcome

We use Lyapunov-like functions and convex optimization to propagate uncertainty in the initial condition of nonlinear systems governed by ordinary differential equations. We consider the full nonlinear dynamics without approximation, producing rigorous bounds on the expected future value of a quantity of interest even when only limited statistics of the initial condition (e.g., mean and variance) are known. For dynamical systems evolving in compact sets, the best upper (lower) bound coincides with the largest (smallest) expectation among all initial state distributions consistent with the known statistics. For systems governed by polynomial equations and polynomial quantities of interest, one-sided estimates on the optimal bounds can be computed using tools from polynomial optimization and semidefinite programming. Moreover, these numerical bounds provably converge to the optimal ones in the compact case. We illustrate the approach on a van der Pol oscillator and on the Lorenz system in the chaotic regime., Comment: 6 pages, 2 figures, 2 tables. Accepted to American Control Conference 2023

Internally heated convection involves the transfer of heat by fluid motion between a distribution of sources and sinks. Focusing on the balanced case where the total heat added by the sources matches the heat taken away by the sinks, we obtain \emph{a priori} bounds on the minimum mean thermal dissipation $\langle |\nabla T|^2\rangle$ as a measure of the inefficiency of transport. In the advective limit, our bounds scale with the inverse mean kinetic energy of the flow. The constant in this scaling law depends on the source--sink distribution, as we explain both in a pair of examples involving oscillatory or concentrated heating and cooling, and via a general asymptotic variational principle for optimizing transport. Key to our analysis is the solution of a pure advection equation, which we do to find examples of extreme heat transfer by cellular and `pinching' flows. When the flow obeys a momentum equation, our bound is re-expressed in terms of a flux-based Rayleigh number $Ra$ yielding $\langle |\nabla T|^2\rangle\geq CRa^{-\alpha}$. The power $\alpha$ is $0, 2/3$ or $1$ depending on the arrangement of the sources and sinks relative to gravity., Comment: Minor revisions from review, figure added

We analyze two recently proposed methods to establish a priori lower bounds on the minimum of general integral variational problems. The methods, which involve either `occupation measures' or a `pointwise dual relaxation' procedure, are shown to produce the same lower bound under a coercivity hypothesis ensuring their strong duality. We then show by a minimax argument that the methods actually evaluate the minimum for classes of one-dimensional, scalar-valued, or convex multidimensional problems. For generic problems, however, these methods should fail to capture the minimum and produce non-sharp lower bounds. We demonstrate this using two examples, the first of which is one-dimensional and scalar-valued with a non-convex constraint, and the second of which is multidimensional and non-convex in a different way. The latter example emphasizes the existence in multiple dimensions of nonlinear constraints on gradient fields that are ignored by occupation measures, but are built into the finer theory of gradient Young measures., Comment: 24 pages, updated references, minor modifications to text

New bounds are proven on the mean vertical convective heat transport, $\overline{\langle wT \rangle}$, for uniform internally heated (IH) convection in the limit of infinite Prandtl number. For fluid in a horizontally-periodic layer between isothermal boundaries, we show that $\overline{\langle wT \rangle} \leq \frac12 - c R^{-2}$, where $R$ is a nondimensional `flux' Rayleigh number quantifying the strength of internal heating and $c = 216$. Then, $\overline{\langle wT \rangle} = 0$ corresponds to vertical heat transport by conduction alone, while $\overline{\langle wT \rangle} > 0$ represents the enhancement of vertical heat transport upwards due to convective motion. If, instead, the lower boundary is a thermal insulator, then we obtain $\overline{\langle wT \rangle} \leq \frac12 - c R^{-4}$, with $c\approx 0.0107$. This result implies that the Nusselt number $Nu$, defined as the ratio of the total-to-conductive heat transport, satisfies $Nu \lesssim R^{4}$. Both bounds are obtained by combining the background method with a minimum principle for the fluid's temperature and with Hardy--Rellich inequalities to exploit the link between the vertical velocity and temperature. In both cases, power-law dependence on $R$ improves the previously best-known bounds, which, although valid at both infinite and finite Prandtl numbers, approach the uniform bound exponentially with $R$., Comment: 28 pages, 5 figures

Motivated by the application of Lyapunov methods to partial differential equations (PDEs), we study functional inequalities of the form $f(I_1(u),\ldots,I_k(u))\geq 0$ where $f$ is a polynomial, $u$ is any function satisfying prescribed constraints, and $I_1(u),\ldots,I_k(u)$ are integral functionals whose integrands are polynomial in $u$, its derivatives, and the integration variable. We show that such functional inequalities can be strengthened into sufficient polynomial inequalities, which in principle can be checked via semidefinite programming using standard techniques for polynomial optimization. These sufficient conditions can be used also to optimize functionals with affine dependence on tunable parameters whilst ensuring their nonnegativity. Our approach relies on a measure-theoretic lifting of the original functional inequality, which extends both a recent moment relaxation strategy for PDE analysis and a dual approach to inequalities for integral functionals., Comment: 8 pages. Submitted to IFAC for possible publication

We obtain an analytical bound on the mean vertical convective heat flux $\langle w T \rangle$ between two parallel boundaries driven by uniform internal heating. We consider two configurations, one with both boundaries held at the same constant temperature, and the other one with a top boundary held at constant temperature and a perfectly insulating bottom boundary. For the first configuration, Arslan et al. (J. Fluid Mech. 919:A15, 2021) recently provided numerical evidence that Rayleigh-number-dependent corrections to the only known rigorous bound $\langle w T \rangle \leq 1/2$ may be provable if the classical background method is augmented with a minimum principle stating that the fluid's temperature is no smaller than that of the top boundary. Here, we confirm this fact rigorously for both configurations by proving bounds on $\langle wT \rangle$ that approach $1/2$ exponentially from below as the Rayleigh number is increased. The key to obtaining these bounds are inner boundary layers in the background fields with a particular inverse-power scaling, which can be controlled in the spectral constraint using Hardy and Rellich inequalities. These allow for qualitative improvements in the analysis not available to standard constructions., Comment: 20 pages, 3 figures

We present a method for finding lower bounds on the global infima of integral variational problems, wherein $\int_\Omega f(x,u(x),\nabla u(x)){\rm d}x$ is minimized over functions $u\colon\Omega\subset\mathbb{R}^n\to\mathbb{R}^m$ satisfying given equality or inequality constraints. Each constraint may be imposed over $\Omega$ or its boundary, either pointwise or in an integral sense. These global minimizations are generally non-convex and intractable. We formulate a particular convex maximization, here called the pointwise dual relaxation (PDR), whose supremum is a lower bound on the infimum of the original problem. The PDR can be derived by dualizing and relaxing the original problem; its constraints are pointwise equalities or inequalities over finite-dimensional sets, rather than over infinite-dimensional function spaces. When the original minimization can be specified by polynomial functions of $(x,u,\nabla u)$, the PDR can be further relaxed by replacing pointwise inequalities with polynomial sum-of-squares (SOS) conditions. The resulting SOS program is computationally tractable when the dimensions $m,n$ and number of constraints are not too large. The framework presented here generalizes an approach of Valmorbida, Ahmadi, and Papachristodoulou (IEEE Trans. Automat. Contr., 61:1649--1654, 2016). We prove that the optimal lower bound given by the PDR is sharp for several classes of problems, whose special cases include leading eigenvalues of Sturm-Liouville problems and optimal constants of Poincar\'e inequalities. For these same classes, we prove that SOS relaxations of the PDR converge to the sharp lower bound as polynomial degrees are increased. Convergence of SOS computations in practice is illustrated for several examples.

The background method is a widely used technique to bound mean properties of turbulent flows rigorously. This work reviews recent advances in the theoretical formulation and numerical implementation of the method. First, we describe how the background method can be formulated systematically within a broader "auxiliary function" framework for bounding mean quantities, and explain how symmetries of the flow and constraints such as maximum principles can be exploited. All ideas are presented in a general setting and are illustrated on Rayleigh-B\'enard convection between stress-free isothermal plates. Second, we review a semidefinite programming approach and a timestepping approach to optimizing bounds computationally, revealing that they are related to each other through convex duality and low-rank matrix factorization. Open questions and promising directions for further numerical analysis of the background method are also outlined., Comment: Revised version: 29 pages, 3 figures, 2 tables. Added simple example to Section 1 and complexity estimates to Section 4. Improved discussion in Section 4

Chordal and factor-width decomposition methods for semidefinite programming and polynomial optimization have recently enabled the analysis and control of large-scale linear systems and medium-scale nonlinear systems. Chordal decomposition exploits the sparsity of semidefinite matrices in a semidefinite program (SDP), in order to formulate an equivalent SDP with smaller semidefinite constraints that can be solved more efficiently. Factor-width decompositions, instead, relax or strengthen SDPs with dense semidefinite matrices into more tractable problems, trading feasibility or optimality for lower computational complexity. This article reviews recent advances in large-scale semidefinite and polynomial optimization enabled by these two types of decomposition, highlighting connections and differences between them. We also demonstrate that chordal and factor-width decompositions allow for significant computational savings on a range of classical problems from control theory, and on more recent problems from machine learning. Finally, we outline possible directions for future research that have the potential to facilitate the efficient optimization-based study of increasingly complex large-scale dynamical systems., Comment: 49 pages, 21 figures, 4 tables

We prove a new rigorous bound for the mean convective heat transport $\langle w T \rangle$, where $w$ and $T$ are the nondimensional vertical velocity and temperature, in internally heated convection between an insulating lower boundary and an upper boundary with a fixed heat flux. The quantity $\langle wT \rangle$ is equal to half the ratio of convective to conductive vertical heat transport, and also to $\frac12$ plus the mean temperature difference between the top and bottom boundaries. An analytical application of the background method based on the construction of a quadratic auxiliary function yields $\langle w T \rangle \leq \tfrac{1}{2}\big(\tfrac{1}{2}+ \tfrac{1}{\sqrt{3}} \big) - 1.6552\, R^{-\frac13}$ uniformly in the Prandtl number, where $R$ is the nondimensional control parameter measuring the strength of the internal heating. Numerical optimisation of the auxiliary function suggests that the asymptotic value of this bound and the $-1/3$ exponent are optimal within our bounding framework. This new result halves the best existing (uniform in $R$) bound (Goluskin 2016, Springer, Table 1.2) and its dependence on $R$ is consistent with previous conjectures and heuristic scaling arguments. Contrary to physical intuition, however, it does not rule out a mean heat transport larger than $\frac12$ at high $R$, which corresponds to the top boundary being hotter than the bottom one on average., Comment: 11 pages, 3 figures

The mean vertical heat transport $\langle wT \rangle$ in convection between isothermal plates driven by uniform internal heating is investigated by means of rigorous bounds. These are obtained as a function of the Rayleigh number $R$ by constructing feasible solutions to a convex variational problem, derived using a formulation of the classical background method in terms of quadratic auxiliary functions. When the fluid's temperature relative to the boundaries is allowed to be positive or negative, numerical solution of the variational problem shows that best previous bound $\langle wT \rangle \leq 1/2$ can only be improved up to finite $R$. Indeed, we demonstrate analytically that $ \langle wT \rangle \leq 2^{-21/5} R^{1/5}$ and therefore prove that $\langle wT\rangle< 1/2$ for $R < 65\,536$. However, if the minimum principle for temperature is invoked, which asserts that internal temperature is at least as large as the temperature of the isothermal boundaries, then numerically optimised bounds are strictly smaller than $1/2$ until at least $R =3.4\times 10^{5}$. While the computational results suggest that the best bound on $\langle wT\rangle$ approaches $1/2$ asymptotically from below as $R\rightarrow \infty$, we prove that typical analytical constructions cannot be used to prove this conjecture., Comment: 33 pages, 12 figures

We present a novel method to compute unstable periodic orbits (UPOs) that optimize the infinite-time average of a given quantity for polynomial ODE systems. The UPO search procedure relies on polynomial optimization to construct nonnegative polynomials whose sublevel sets approximately localize parts of the optimal UPO, and that can be used to implement a simple yet effective control strategy to reduce the UPO's instability. Precisely, we construct a family of controlled ODE systems parameterized by a scalar $k$ such that the original ODE system is recovered for $k = 0$, and such that the optimal orbit is less unstable, or even stabilized, for $k>0$. Periodic orbits for the controlled system can be more easily converged with traditional methods and numerical continuation in $k$ allows one to recover optimal UPOs for the original system. The effectiveness of this approach is illustrated on three low-dimensional ODE systems with chaotic dynamics., Comment: 17 pages, 6 figures

We prove decomposition theorems for sparse positive (semi)definite polynomial matrices that can be viewed as sparsity-exploiting versions of the Hilbert--Artin, Reznick, Putinar, and Putinar--Vasilescu Positivstellens\"atze. First, we establish that a polynomial matrix $P(x)$ with chordal sparsity is positive semidefinite for all $x\in \mathbb{R}^n$ if and only if there exists a sum-of-squares (SOS) polynomial $\sigma(x)$ such that $\sigma P$ is a sum of sparse SOS matrices. Second, we show that setting $\sigma(x)=(x_1^2 + \cdots + x_n^2)^\nu$ for some integer $\nu$ suffices if $P$ is homogeneous and positive definite globally. Third, we prove that if $P$ is positive definite on a compact semialgebraic set $\mathcal{K}=\{x:g_1(x)\geq 0,\ldots,g_m(x)\geq 0\}$ satisfying the Archimedean condition, then $P(x) = S_0(x) + g_1(x)S_1(x) + \cdots + g_m(x)S_m(x)$ for matrices $S_i(x)$ that are sums of sparse SOS matrices. Finally, if $\mathcal{K}$ is not compact or does not satisfy the Archimedean condition, we obtain a similar decomposition for $(x_1^2 + \ldots + x_n^2)^\nu P(x)$ with some integer $\nu\geq 0$ when $P$ and $g_1,\ldots,g_m$ are homogeneous of even degree. Using these results, we find sparse SOS representation theorems for polynomials that are quadratic and correlatively sparse in a subset of variables, and we construct new convergent hierarchies of sparsity-exploiting SOS reformulations for convex optimization problems with large and sparse polynomial matrix inequalities. Numerical examples demonstrate that these hierarchies can have a significantly lower computational complexity than traditional ones., Comment: 32 pages, 7 figures, 4 tables. Established sparsity-exploiting versions of Reznick and Putinar--Vasilescu Positivstellens\"atze, and updated the second numerical example. Code for our numerical experiments is available here: https://github.com/aeroimperial-optimization/sos-chordal-decomposition-pmi

We prove a new rigorous upper bound on the vertical heat transport for B\'enard-Marangoni convection of a two- or three-dimensional fluid layer with infinite Prandtl number. Precisely, for Marangoni number $Ma \gg 1$ the Nusselt number $Nu$ is bounded asymptotically by $Nu \lesssim Ma^{2/7}(\ln Ma)^{-1/7}$. Key to our proof are a background temperature field with a hyperbolic profile near the fluid's surface, and new estimates for the coupling between temperature and vertical velocity., Comment: 10 pages, 1 figure

We study a convex optimization framework for bounding extreme events in nonlinear dynamical systems governed by ordinary or partial differential equations (ODEs or PDEs). This framework bounds from above the largest value of an observable along trajectories that start from a chosen set and evolve over a finite or infinite time interval. The approach needs no explicit trajectories. Instead, it requires constructing suitably constrained auxiliary functions that depend on the state variables and possibly on time. Minimizing bounds over auxiliary functions is a convex problem dual to the non-convex maximization of the observable along trajectories. This duality is strong, meaning that auxiliary functions give arbitrarily sharp bounds, for sufficiently regular ODEs evolving over a finite time on a compact domain. When these conditions fail, strong duality may or may not hold; both situations are illustrated by examples. We also show that near-optimal auxiliary functions can be used to construct spacetime sets that localize trajectories leading to extreme events. Finally, in the case of polynomial ODEs and observables, we describe how polynomial auxiliary functions of fixed degree can be optimized numerically using polynomial optimization. The corresponding bounds become sharp as the polynomial degree is raised if strong duality and mild compactness assumptions hold. Analytical and computational ODE examples illustrate the construction of bounds and the identification of extreme trajectories, along with some limitations. As an analytical PDE example, we bound the maximum fractional enstrophy of solutions to the Burgers equation with fractional diffusion., Comment: Revised according to comments by reviewers. Added references and rearranged introduction, conclusions, and proofs. 38 pages, 7 figures, 4 tables, 4 appendices, 87 references

We prove weak duality between two recent convex relaxation methods for bounding the optimal value of a constrained variational problem in which the objective is an integral functional. The first approach, proposed by Valmorbida et al. (IEEE Trans. Automat. Control 61(6):1649--1654, 2016), replaces the variational problem with a convex program over sufficiently smooth functions, subject to pointwise non-negativity constraints. The second approach, discussed by Korda et al. (arXiv:1804.07565v1 [math.OC]), relaxes the variational problem into a convex program over scaled probability measures. We also prove that the duality between these infinite-dimensional convex programs is strong, meaning that their optimal values coincide, when the range and gradients of admissible functions in the variational problem are constrained to bounded sets. For variational problems with polynomial data, the optimal values of each convex relaxation can be approximated by solving weakly dual hierarchies of finite-dimensional semidefinite programs (SDPs). These are strongly dual under standard constraint qualification conditions irrespective of whether strong duality holds at the infinite-dimensional level. Thus, the two relaxation approaches are equivalent for the purposes of computations., Comment: 11 pages

Tobasco et al. [Physics Letters A, 382:382-386, 2018; see https://doi.org/10.1016/j.physleta.2017.12.023] recently suggested that trajectories of ODE systems that optimize the infinite-time average of a certain observable can be localized using sublevel sets of a function that arise when bounding such averages using so-called auxiliary functions. In this paper we demonstrate that this idea is viable and allows for the computation of extremal unstable periodic orbits (UPOs) for polynomial ODE systems. First, we prove that polynomial optimization is guaranteed to produce auxiliary functions that yield near-sharp bounds on time averages, which is required in order to localize the extremal orbit accurately. Second, we show that points inside the relevant sublevel sets can be computed efficiently through direct nonlinear optimization. Such points provide good initial conditions for UPO computations. As a proof of concept, we then combine these methods with a single-shooting Newton-Raphson algorithm to study extremal UPOs for a nine-dimensional model of sinusoidally forced shear flow. We discover three previously unknown families of UPOs, one of which simultaneously minimizes the mean energy dissipation rate and maximizes the mean perturbation energy relative to the laminar state for Reynolds numbers approximately between 81.24 and 125., Comment: 22 pages, 7 figures. v3: update discussion after Theorem 1, fixed typos

Optimization over non-negative polynomials is fundamental for nonlinear systems analysis and control. We investigate the relation between three tractable relaxations for optimizing over sparse non-negative polynomials: sparse sum-of-squares (SSOS) optimization, diagonally dominant sum-of-squares (DSOS) optimization, and scaled diagonally dominant sum-of-squares (SDSOS) optimization. We prove that the set of SSOS polynomials, an inner approximation of the cone of SOS polynomials, strictly contains the spaces of sparse DSOS/SDSOS polynomials. When applicable, therefore, SSOS optimization is less conservative than its DSOS/SDSOS counterparts. Numerical results for large-scale sparse polynomial optimization problems demonstrate this fact, and also that SSOS optimization can be faster than DSOS/SDSOS methods despite requiring the solution of semidefinite programs instead of less expensive linear/second-order cone programs., Comment: 9 pages, 3 figures

This paper introduces a notion of decomposition and completion of sum-of-squares (SOS) matrices. We show that a subset of sparse SOS matrices with chordal sparsity patterns can be equivalently decomposed into a sum of multiple SOS matrices that are nonzero only on a principal submatrix. Also, the completion of an SOS matrix is equivalent to a set of SOS conditions on its principal submatrices and a consistency condition on the Gram representation of the principal submatrices. These results are partial extensions of chordal decomposition and completion of scalar matrices to matrices with polynomial entries. We apply the SOS decomposition result to exploit sparsity in matrix-valued SOS programs. Numerical results demonstrate the high potential of this approach for solving large-scale sparse matrix-valued SOS programs., Comment: 6 pages

We present methods for bounding infinite-time averages in dynamical systems governed by nonlinear PDEs. The methods rely on auxiliary functionals, which are similar to Lyapunov functionals but satisfy different inequalities. The inequalities are enforced by requiring certain expressions to be sums of squares of polynomials, and the optimal choice of auxiliary functional is posed as a semidefinite program (SDP) that can be solved computationally. To formulate these SDPs we approximate the PDE by truncated systems of ODEs and proceed in one of two ways. The first approach is to compute bounds for the ODE systems, increasing the truncation order until bounds converge numerically. The second approach incorporates the ODE systems with analytical estimates on their deviation from the PDE, thereby using finite truncations to produce bounds for the full PDE. We apply both methods to the Kuramoto-Sivashinsky equation, where we compute upper bounds on the spatiotemporal average of energy by employing polynomial auxiliary functionals up to degree six. The first approach is used for most computations, but a subset of results are checked using the second approach, and the results agree to high precision. These bounds apply to all odd solutions of period $2\pi L$, where $L$ is varied. Sharp bounds are obtained for $L\le10$, and trends suggest that more expensive computations would yield sharp bounds at larger $L$ also. The bounds are known to be sharp (to within 0.1% numerical error) because they are saturated by the simplest nonzero steady states, which apparently have the largest mean energy among all odd solutions. Prior authors have conjectured that mean energy remains $O(1)$ for $L\gg1$ since no particular solutions with larger energy have been found. Our bounds constitute the first positive evidence for this conjecture, albeit up to finite $L$, and they offer some guidance for analytical proofs., Comment: 32 pages, 6 figures; v2: revisions after review, new proposition added

We prove the first rigorous bound on the heat transfer for three-dimensional Rayleigh-B\'enard convection of finite-Prandtl-number fluids between free-slip boundaries with an imposed heat flux. Using the auxiliary functional method with a quadratic functional, which is equivalent to the background method, we prove that the Nusselt number $Nu$ is bounded by $Nu \leq 0.5999 R^{1/3}$ uniformly in the Prandtl number, where $R$ is the Rayleigh number based on the imposed heat flux. In terms of the Rayleigh number based on the mean vertical temperature drop, $Ra$, we obtain $Nu \leq 0.4646 Ra^{1/2}$. The scaling with Rayleigh number is the same as that of bounds obtained with no-slip isothermal, free-slip isothermal, and no-slip fixed flux boundaries, and numerical optimization of the bound suggests that it cannot be improved within our bounding framework. Contrary to the two-dimensional case, therefore, the $Ra$-dependence of rigorous upper bounds on the heat transfer obtained with the background method for three-dimensional Rayleigh-B\'enard convection is insensitive to both the thermal and the velocity boundary conditions., Comment: 5 pages, 3 figures. Version 2: changed format to JFM, minor rewriting following comments by referees, added Appendix A

When sum-of-squares (SOS) programs are recast as semidefinite programs (SDPs) using the standard monomial basis, the constraint matrices in the SDP possess a structural property that we call \emph{partial orthogonality}. In this paper, we leverage partial orthogonality to develop a fast first-order method, based on the alternating direction method of multipliers (ADMM), for the solution of the homogeneous self-dual embedding of SDPs describing SOS programs. Precisely, we show how a "diagonal plus low rank" structure implied by partial orthogonality can be exploited to project efficiently the iterates of a recent ADMM algorithm for generic conic programs onto the set defined by the affine constraints of the SDP. The resulting algorithm, implemented as a new package in the solver CDCS, is tested on a range of large-scale SOS programs arising from constrained polynomial optimization problems and from Lyapunov stability analysis of polynomial dynamical systems. These numerical experiments demonstrate the effectiveness of our approach compared to common state-of-the-art solvers., Comment: Added the results for matrix-valued SOS programs and weighted SOS contraints. Codes are available through https://github.com/OxfordControl/CDCS

The vertical heat transfer in B\'enard-Marangoni convection of a fluid layer with infinite Prandtl number is studied by means of upper bounds on the Nusselt number $Nu$ as a function of the Marangoni number $Ma$. Using the background method for the temperature field, it has recently been proven by Hagstrom & Doering that $ Nu\leq 0.838\,Ma^{2/7}$. In this work we extend previous background method analysis to include balance parameters and derive a variational principle for the bound on $Nu$, expressed in terms of a scaled background field, that yields a better bound than Hagstrom & Doering's formulation at a given $Ma$. Using a piecewise-linear, monotonically decreasing profile we then show that $Nu \leq 0.803\,Ma^{2/7}$, lowering the previous prefactor by 4.2%. However, we also demonstrate that optimisation of the balance parameters does not affect the asymptotic scaling of the optimal bound achievable with Hagstrom & Doering's original formulation. We subsequently utilise convex optimisation to optimise the bound on $Nu$ over all admissible background fields, as well as over two smaller families of profiles constrained by monotonicity and convexity. The results show that $Nu \leq O(Ma^{2/7}(\ln Ma)^{-1/2})$ when the background field has a non-monotonic boundary layer near the surface, while a power-law bound with exponent 2/7 is optimal within the class of monotonic background fields. Further analysis of our upper-bounding principle reveals the role of non-monotonicity, and how it may be exploited in a rigorous mathematical argument., Comment: Revised version: 33 pages, 9 figures. Extended discussion in Sections 6 and 7. Fixed mistakes in bibliography. Fixed typos. JFM style with patch for author-year references with hyperref

We employ chordal decomposition to reformulate a large and sparse semidefinite program (SDP), either in primal or dual standard form, into an equivalent SDP with smaller positive semidefinite (PSD) constraints. In contrast to previous approaches, the decomposed SDP is suitable for the application of first-order operator-splitting methods, enabling the development of efficient and scalable algorithms. In particular, we apply the alternating direction method of multipliers (ADMM) to solve decomposed primal- and dual-standard-form SDPs. Each iteration of such ADMM algorithms requires a projection onto an affine subspace, and a set of projections onto small PSD cones that can be computed in parallel. We also formulate the homogeneous self-dual embedding (HSDE) of a primal-dual pair of decomposed SDPs, and extend a recent ADMM-based algorithm to exploit the structure of our HSDE. The resulting HSDE algorithm has the same leading-order computational cost as those for the primal or dual problems only, with the advantage of being able to identify infeasible problems and produce an infeasibility certificate. All algorithms are implemented in the open-source MATLAB solver CDCS. Numerical experiments on a range of large-scale SDPs demonstrate the computational advantages of the proposed methods compared to common state-of-the-art solvers., Comment: To appear at Math. Prog. 36 pages, 7 figures; Codes available from https://github.com/oxfordcontrol/CDCS (sparse conic solver: CDCS)

This paper introduces an efficient first-order method based on the alternating direction method of multipliers (ADMM) to solve semidefinite programs (SDPs) arising from sum-of-squares (SOS) programming. We exploit the sparsity of the \emph{coefficient matching conditions} when SOS programs are formulated in the usual monomial basis to reduce the computational cost of the ADMM algorithm. Each iteration of our algorithm requires one projection onto the positive semidefinite cone and the solution of multiple quadratic programs with closed-form solutions free of any matrix inversion. Our techniques are implemented in the open-source MATLAB solver SOSADMM. Numerical experiments on SOS problems arising from unconstrained polynomial minimization and from Lyapunov stability analysis for polynomial systems show speed-ups compared to the interior-point solver SeDuMi, and the first-order solver CDCS., Comment: 6 pages; Codes available from https://github.com/zhengy09/SOSADMM; accepted in the IEEE Control Systems Letters

Using the energy method we investigate the stability of pure conduction in Pearson's model for B\'enard-Marangoni convection in a layer of fluid at infinite Prandtl number. Upon extending the space of admissible perturbations to the conductive state, we find an exact solution to the energy stability variational problem for a range of thermal boundary conditions describing perfectly conducting, imperfectly conducting, and insulating boundaries. Our analysis extends and improves previous results, and shows that with the energy method global stability can be proven up to the linear instability threshold only when the top and bottom boundaries of the fluid layer are insulating. Contrary to the well-known Rayleigh-B\'enard convection setup, therefore, energy stability theory does not exclude the possibility of subcritical instabilities against finite-amplitude perturbations., Comment: 11 pages, 2 figures. Preprint submitted to the Journal of Fluid Mechanics. Version 2: minor text and notational changes, added a new appendix A, added detail to Section 3

We propose an efficient first-order method, based on the alternating direction method of multipliers (ADMM), to solve the homogeneous self-dual embedding problem for a primal-dual pair of semidefinite programs (SDPs) with chordal sparsity. Using a series of block eliminations, the per-iteration cost of our method is the same as applying a splitting method to the primal or dual alone. Moreover, our approach is more efficient than other first-order methods for generic sparse conic programs since we work with smaller semidefinite cones. In contrast to previous first-order methods that exploit chordal sparsity, our algorithm returns both primal and dual solutions when available, and a certificate of infeasibility otherwise. Our techniques are implemented in the open-source MATLAB solver CDCS. Numerical experiments on three sets of benchmark problems from the library SDPLIB show speed-ups compared to some common state-of-the-art software packages., Comment: 6 pages; Codes are available from https://github.com/giofantuzzi/CDCS/tree/developer (conic solver CDCS); accepted in the IFAC 2017

Many problems in control theory can be formulated as semidefinite programs (SDPs). For large-scale SDPs, it is important to exploit the inherent sparsity to improve the scalability. This paper develops efficient first-order methods to solve SDPs with chordal sparsity based on the alternating direction method of multipliers (ADMM). We show that chordal decomposition can be applied to either the primal or the dual standard form of a sparse SDP, resulting in scaled versions of ADMM algorithms with the same computational cost. Each iteration of our algorithms consists of a projection on the product of small positive semidefinite cones, followed by a projection on an affine set, both of which can be carried out efficiently. Our techniques are implemented in CDCS, an open source add-on to MATLAB. Numerical experiments on large-scale sparse problems in SDPLIB and random SDPs with block-arrow sparse patterns show speedups compared to some common state-of-the-art software packages., Comment: 6 pages. Codes available from http://github.com/giofantuzzi/CDCS . Accepted in the American Control Conference, 2017 http://ieeexplore.ieee.org/abstract/document/7963462/

We introduce a new technique to optimize a linear cost function subject to a one-dimensional affine homogeneous quadratic integral inequality, i.e., the requirement that a homogeneous quadratic integral functional, affine in the optimization variables, is non-negative over a space of functions defined by homogeneous boundary conditions. Such problems arise in stability analysis, input-to-state/output analysis, and control of many systems governed by partial differential equations (PDEs), in particular fluid dynamical systems. First, we derive outer approximations for the feasible set of a homogeneous quadratic integral inequality in terms of linear matrix inequalities (LMIs), and show that under mild assumptions a convergent, non-decreasing sequence of lower bounds for the optimal cost can be computed with a sequence of semidefinite programs (SDPs). Second, we obtain inner approximations in terms of LMIs and sum-of-squares constraints, so upper bounds for the optimal cost and strictly feasible points for the integral inequality can also be computed with SDPs. To aid the formulation and solution of our SDP relaxations, we implement our techniques in QUINOPT, an open-source add-on to YALMIP. We demonstrate our techniques by solving problems arising from the stability analysis of PDEs., Comment: 31 pages, 3 figures. Codes available from https://github.com/aeroimperial-optimization/QUINOPT/. Latest version: improved exposition, fixed proof of Theorem 1, added comparison to other method, added discussion on scalability, corrected typos

We describe methods for proving upper and lower bounds on infinite-time averages in deterministic dynamical systems and on stationary expectations in stochastic systems. The dynamics and the quantities to be bounded are assumed to be polynomial functions of the state variables. The methods are computer-assisted, using sum-of-squares polynomials to formulate sufficient conditions that can be checked by semidefinite programming. In the deterministic case, we seek tight bounds that apply to particular local attractors. An obstacle to proving such bounds is that they do not hold globally; they are generally violated by trajectories starting outside the local basin of attraction. We describe two closely related ways past this obstacle: one that requires knowing a subset of the basin of attraction, and another that considers the zero-noise limit of the corresponding stochastic system. The bounding methods are illustrated using the van der Pol oscillator. We bound deterministic averages on the attracting limit cycle above and below to within 1%, which requires a lower bound that does not hold for the unstable fixed point at the origin. We obtain similarly tight upper and lower bounds on stochastic expectations for a range of noise amplitudes. Limitations of our methods for certain types of deterministic systems are discussed, along with prospects for improvement., Comment: 25 pages; Added new Section 7.2; Added references; Corrected typos; Submitted to SIADS