Duality of convex relaxations for constrained variational problems

TL;DR

This paper establishes the weak duality and strong duality between two convex relaxation methods for constrained variational problems, showing their equivalence in computational applications, especially with polynomial data.

Contribution

It proves the duality relationship and equivalence of two recent convex relaxation approaches for variational problems, including their finite-dimensional SDP approximations.

Findings

Weak duality between the two convex relaxation methods.

Strong duality when admissible functions are bounded.

Equivalence of the relaxation approaches for polynomial data.

Abstract

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.

Equations

$$L[u]:={\int }_{X}l(x,u,Du)\mathrm{d}x$$

$$\displaystyle\mathcal{U}:=\big{\{}u\in W^{1,\infty}(X;\,\mathbb{R}^{m}):\;\textstyle\int_{X}f(x,u,Du)\,{\rm d}x$$

$$g(x,u,Du)$$

$$h_1(x,u)$$

$$⋮$$

$$h_s(x,u)$$

$$u\in Y$$

$$Du\in Z$$

$$L^{\ast }:=\underset{u\in \mathcal{U}}{\inf }L[u]$$

$$L^{\ast }=\sup \left\{b\in \mathbb{R}:{\int }_{X}l(x,u,Du)\mathrm{d}x\ge b\forall u\in \mathcal{U}\right\},$$

$$\begin{gathered} L^{\ast }=\underset{\begin{array}{c}u\in W^{1,\infty }(X;{\mathbb{R}}^m)\\ u(x)\in Y\text{ a.e.}\\ Du(x)\in Z\text{ a.e.}\end{array}}{\inf }\underset{\begin{array}{c}a\in \mathbb{R}\\ \tilde{\varphi }\in C(X)\\ {\tilde{\psi }}_i\in C(\partial X_i)\end{array}}{\sup }\{{\int }_X\left[l(x,u,Du)+af(x,u,Du)+\tilde{\varphi }(x)g(x,u,Du)\right]\mathrm{d}x \\ +\sum _{i=1}^s{\int }_{\partial X_i}{\tilde{\psi }}_i(x)h_i(x,u)\mathrm{d}\sigma \}, \end{gathered}$$

$$\tilde{\varphi }(x)=\varphi [x,u(x),Du(x)],$$

$${\tilde{\psi }}_i(x)={\psi }_i[x,u(x)],i=1,\dots ,s.$$

$${\int }_X\left\{{\nabla }_x\cdot \theta (x,u)+\mathrm{tr}[v{D}_u\theta (x,u)]\right\}\mathrm{d}x-{\int }_{\partial X}\theta (x,u)\cdot \hat{n}(x)\mathrm{d}\sigma =0\forall \theta \in C^1(\overline{X}\times \overline{Y};{\mathbb{R}}^n),$$

$$L^{\ast }=\underset{\begin{array}{c}u\in W^{1,\infty }(X;{\mathbb{R}}^m)\\ u(x)\in Y\text{ a.e.}\\ v(x)\in Z\text{ a.e.}\end{array}}{\inf }\underset{\begin{array}{c}\alpha \in \mathbb{R}\\ \varphi \in C(\Omega )\\ {\psi }_i\in C({\Gamma }_i)\\ \theta \in C^1(\overline{X}\times \overline{Y};{\mathbb{R}}^n)\end{array}}{\sup }\left\{{\int }_{X}F(x,u,v)\mathrm{d}x+\sum _{i=1}^s{\int }_{\partial X_i}{G}_i(x,u)\mathrm{d}\sigma \right\},$$

$$F(x,u,v):=$$

$$+\varphi (x,u,v)g(x,u,v)+{\nabla }_x\cdot \theta (x,u)+\mathrm{tr}\left[v{D}_u\theta (x,u)\right],$$

$$G_i(x,u):=$$

$$L^{\ast }\ge \underset{\begin{array}{c}a\in \mathbb{R}\\ \varphi \in C(\Omega )\\ {\psi }_i\in C({\Gamma }_i)\\ \theta \in C^1(\overline{X}\times \overline{Y};{\mathbb{R}}^n)\end{array}}{\sup }\left\{\left|X\right|\underset{(x,y,z)\in \Omega }{\inf }F(x,y,z)+\sum _{i=1}^s\left|\partial X_i\right|\underset{(x,y)\in {\Gamma }_i}{\inf }{G}_i(x,y)\right\}=:\mathbb{D}.$$

$$\mathbb{D}=\sup_{\begin{subarray}{c}a,b,c_{i}\in\mathbb{R}\\ \phi\in C(\Omega)\\ \psi_{i}\in C(\Gamma_{i})\\ \theta\in C^{1}(\overline{X}\times\overline{Y};\,\mathbb{R}^{n})\end{subarray}}\!\!\!\!\!\!\begin{aligned} \\ \bigg{\{}\left|X\right|b+\sum_{i=1}^{s}\left|\partial X_{i}\right|c_{i}:\quad F(x,y,z)-b&\geq 0\text{ on }\Omega,\\[-12.91663pt] G_{i}(x,y)-c_{i}&\geq 0\text{ on }\Gamma_{i},\quad i=1,\,\ldots,\,s\bigg{\}}.\end{aligned}$$

$$\Omega =\overline{X}\times \overline{Y}\times \overline{Z}=\{(x,y,z):p_1(x,y,z)\ge 0,\dots ,p_r(x,y,z)\ge 0\}.$$

$${\Gamma }_i:=\partial X_i\times \overline{Y}=\{(x,y):q_{i,1}(x,y)\ge 0,\dots ,q_{i,t_i}(x,y)\ge 0\}.$$

$$Q(\Omega ):=\{w\in \mathbb{R}[x,y,z]:w={\sigma }_0+{\sigma }_1{p}_1+\cdots +{\sigma }_r{p}_r,{\sigma }_0,\dots ,{\sigma }_r\in \Sigma [x,y,z]\},$$

$$Q({\Gamma }_i):=\{w\in \mathbb{R}[x,y]:w={\sigma }_0+{\sigma }_1{q}_{i,1}+\cdots +{\sigma }_s{q}_{i,t_i},{\sigma }_0,\dots ,{\sigma }_{t_i}\in \Sigma [x,y]\},$$

$$\mathbb{D}\geq\sup_{\begin{subarray}{c}a,b,c_{i}\in\mathbb{R}\\ \phi\in\mathbb{R}_{d}[x,y,z]\\ \psi_{i}\in\mathbb{R}_{d}[x,y]\\ \theta_{1},\ldots,\theta_{n}\in\mathbb{R}_{d}[x,y]\end{subarray}}\;\begin{aligned} \\ \big{\{}\left|X\right|b+\left|\partial X\right|c:\quad F(x,y,z)-b&\in Q(\Omega),\\[-1.07639pt] G_{i}(x,y)-c&\in Q(\Gamma_{i}),\quad i=1,\,\ldots,\,s\big{\}}.\end{aligned}$$

$$\begin{array}{rl}{\zeta }^u:X& \to \overline{X}\times \overline{Y}\times \overline{Z}\\ x& \mapsto (x,u(x),Du(x)),\end{array}\begin{array}{rl}{\chi }_i^u:\partial X_i& \to \partial X_i\times \overline{Y}\\ x& \mapsto (x,u(x)).\end{array}$$

$${\int }_X\eta (x,u,Du)\mathrm{d}x={\iiint }_{\Omega }\eta (x,y,z)\mathrm{d}{\mu }_u(x,y,z)=:⟨\eta ,{\mu }^u⟩,$$

$${\int }_{\partial X_i}\xi (x,u)\mathrm{d}\sigma ={\iint }_{{\Gamma }_i}\xi (x,y)\mathrm{d}{\nu }_i^u(x,y)=:⟨\xi ,{\nu }_i^u⟩,$$

$$L^{\ast }=\underset{{\mu }^u}{\inf }⟨l,{\mu }^u⟩.$$

$$\Vert {\mu }^u\Vert :=⟨1,{\mu }^u⟩=\left|X\right|,$$

$$\Vert {\nu }_i^u\Vert :=⟨1,{\nu }_i^u⟩=\left|\partial X_i\right|.$$

$$⟨f,{\mu }^u⟩=0.$$

$$⟨\varphi g,{\mu }^u⟩=0\forall \varphi \in C(\Omega ).$$

$$⟨{\psi }_i{h}_i,{\nu }_i^u⟩=0\forall {\psi }_i\in C({\Gamma }_i),i=1,\dots ,s.$$

$$⟨{\nabla }_x\cdot \theta +\mathrm{tr}(z{D}_y\theta ),{\mu }^u⟩-\sum _{i=1}^s⟨\theta \cdot \hat{n},{\nu }_i^u⟩=0\forall \theta \in C^1(\overline{X}\times \overline{Y};{\mathbb{R}}^n).$$

$$\mathcal{M}:=\{(\mu ,{\nu }_1,\dots ,{\nu }_s):\mu ,{\nu }_1,\dots ,{\nu }_s\ge 0\text{ and satisfy (refe:measure-normsa,b)–lx@crefcreftype refnume:measure-ibp}\}.$$

Full text

1. Introduction

Constrained variational problems arise in a large number of fields, including nonlinear elasticity, fluid mechanics, and control theory. In this paper, we consider a general class of variational problems relevant to all these application domains: minimize an integral functional,

[TABLE]

where $l\in C(\overline{X}\times\mathbb{R}^{n}\times\mathbb{R}^{m\times n})$ is given, over all functions $u$ that belong to the set

[TABLE]

In these expressions and throughout this work $X\subset\mathbb{R}^{n}$ is an open bounded domain, whose boundary $\partial X=\overline{X}\setminus X$ is Lipschitz and consists of $s$ smooth surfaces $\partial X_{1},\,\ldots,\,\partial X_{s}$ that are disjoint up to a set of zero surface measure. The Sobolev space $W^{1,\infty}(X;\,\mathbb{R}^{m})$ consists of all weakly differentiable and essentially bounded functions $u:X\to\mathbb{R}^{m}$ whose Jacobian matrix $Du=(\partial_{x_{j}}u_{i})_{i=1,\,\ldots,\,m}^{j=1,\,\ldots,\,n}$ is also essentially bounded. The functions $f,g\in C(\overline{X}\times\mathbb{R}^{n}\times\mathbb{R}^{m\times n})$ and $h_{1}\in C(\partial{X}_{1}\times\mathbb{R}^{n}),\,\ldots,\,h_{s}\in C(\partial{X}_{s}\times\mathbb{R}^{n})$ define integral, differential and boundary constraint on $u$. The sets $Y$ and $Z$, to which the range and gradients of admissible functions are restricted, may or may not be bounded and may coincide with the full spaces $\mathbb{R}^{m}$ and $\mathbb{R}^{m\times n}$. We assume that $\mathcal{U}$ is nonempty and $L[u]$ is bounded below on $\mathcal{U}$, but not that a minimizer exists.

Minimizing $L[u]$ over the set $\mathcal{U}$ analytically is often impossible. The minimum value

[TABLE]

is typically approximated numerically either by computing a minimizing sequence using direct optimization techniques, or by discretizing and solving the Euler–Lagrange partial differential equations (PDEs). Unless 1.3 is a convex minimization problem, meaning that $L[u]$ is a convex functional and $\mathcal{U}$ is a convex set, such methods return approximate local minimizers with no way of checking whether they are globally optimal. In general, therefore, one only obtains an upper bound on $L^{*}$. This works investigates complementary techniques that bound $L^{*}$ from below.

Two such approaches have been proposed recently, and are particularly interesting because the search for lower bounds on $L^{*}$ is posed as a convex optimization problem even when 1.3 is nonconvex. The first approach is to rewrite 1.3 as

[TABLE]

augment the integral inequality using Lagrange multipliers to enforce the constraints that define $\mathcal{U}$ and the differential relation between $u$ and $Du$, and replace the augmented inequality with stronger pointwise inequalities on $\overline{X}\times\overline{Y}\times\overline{Z}$ and $\partial{X}\times\overline{Y}$. What results is an infinite-dimensional convex maximization problem with inequality constraints for $b$ and the Lagrange multipliers. The constraints are polynomial inequalities on semialgebraic sets when the sets $X$, $Y$, $Z$ are semialgebraic, the functions $f$, $g$, $h_{1},\,\ldots,\,h_{s}$, and $l$ are polynomial, and the Lagrange multipliers are restricted to be polynomials of fixed degree. Replacing polynomial non-negativity with sum-of-squares (SOS) conditions enables one to maximize lower bounds on $L^{*}$ numerically by solving a hierarchy of semidefinite programs (SDPs), indexed by the degree of the Lagrange multipliers. This strategy was proposed by Valmorbida et al. [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c] and Ahmadi et al. [Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018] in the context of stability analysis, input-output analysis, and safety verification for dynamical systems governed by polynomial PDEs, but applies equally well to constrained variational problems.

The second strategy to bound $L^{*}$ from below using convex optimization is to reformulate 1.3 as a minimization problem over the set of occupation and boundary measures on $\mathcal{U}$. These are the images of the Lebesgue measure on $X$ and the surface measure on $\partial{X}$ under the map $x\mapsto(x,u(x),Du(x))$ as $u$ varies in $\mathcal{U}$. The set of such measures is nonconvex in general, but relaxing the minimization to a convex superset gives an infinite-dimensional convex problem whose optimal value is a lower bound on $L^{*}$. This, in turn, can be estimated from below using a standard hierarchy of moment-SDP relaxations when the sets $X$, $Y$, $Z$ are semialgebraic and $f,\,g,\,h_{1},\,\ldots,\,h_{s}$ and $l$ are polynomials. These ideas have been applied to analyze and control linear PDEs [Bertsimas2006, Magron2018], and have recently been extended to the nonlinear case by Korda et al. [Korda2018]. Similar methods have also been used to approximate solutions to hyperbolic PDEs [Marx2018].

This work demonstrates that the infinite-dimensional convex problems obtained with the two relaxation methods just described are weakly dual in the sense of convex duality. In addition, we prove that the duality is strong when $Y$ and $Z$ are bounded sets, so both methods yield the same lower bound on $L^{*}$ in this case. The two hierarchies of finite-dimensional SDPs obtained for variational problems with polynomial data are also weakly dual. Strong duality of the SDPs at each level of the hierarchy can be established using standard constraint qualification conditions, which often hold in practice and are independent of whether the infinite-dimensional formulations are strongly dual. Moreover, many popular algorithms for solving SDPs require strong duality to guarantee convergence and avoid poor numerical conditioning. Consequently, the relaxations of 1.3 proposed by [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018] and [Korda2018] are equivalent from the point of view of numerical computations.

We present these new results in section 4 after reviewing the methods of [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018] and [Korda2018] in sections 2 and 3, respectively. While most of the material contained there has appeared in the cited literature, our discussion slightly differs from previous works and, sometimes, extends them. In particular, section 2 gives a new description of the methods of [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018], which does not start with 1.4 and makes the hitherto unrecognized connection with the techniques of [Korda2018] evident. Further comments are offered in section 5.

2. Convex relaxation using Lagrange multipliers

We begin by deriving an infinite-dimensional convex program over Lagrange multipliers, largely based on ideas from [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018], whose feasible solutions prove lower bounds on $L^{*}$. For notational simplicity we will often write $\Omega:=\overline{X}\times\overline{Y}\times\overline{Z}$ and $\Gamma_{i}:=\partial X_{i}\times\overline{Y}$, where $\partial X_{i}$ is any of the $s$ smooth components of the boundary $\partial X$. We will also write $\left|X\right|=\int_{X}{\rm d}x$ for the volume of $X$ and $\left|\partial X_{i}\right|=\int_{\partial X_{i}}{\rm d}\sigma$ for the surface area of each smooth portion of its boundary, where ${\rm d}\sigma$ is the surface measure.

2.1. An infinite-dimensional convex program over Lagrange multipliers

As in the classical approach to solving 1.3 using the calculus of variations, observe that

[TABLE]

where $a$, $\tilde{\phi}$ and $\tilde{\psi}_{1},\,\ldots,\,\tilde{\psi}_{s}$ are Lagrange multipliers for the integral, differential and boundary constraints that define $\mathcal{U}$. The spatial structure of the optimal multipliers clearly depends on $u$ and its derivatives. To enforce this without having to prescribe $u$, we let the multipliers be explicit functions of both $u$ and $Du$. Precisely, without loss of generality we consider functions $\phi\in C(\Omega)$ and $\psi_{i}\in C(\Gamma_{i})$, $i=1,\,\ldots,\,s$, and let

[TABLE]

Another source of difficulty is the differential relation between $u$ and $Du$. To handle this, we introduce a slack function $v\in L^{\infty}(X;\,\mathbb{R}^{m\times n})$ and replace $Du$ by $v$, subject to the constraint $v=Du$. We can then impose the differential relation between $u$ and its derivatives using a special type of vector-valued Lagrange multiplier, without having to consider $Du$ explicitly. The next proposition, which is similar to Lemma 1.1 in [Korda2018], makes this precise.

Proposition 1**.**

Let $u\in W^{1,\infty}(X;\,\mathbb{R}^{m})$ and $v\in L^{\infty}(X;\,\mathbb{R}^{m\times n})$. Then, $v=Du$ a.e. on $X$ if and only if

[TABLE]

where $\nabla_{x}\cdot\theta(x,u)$ is the divergence of $\theta$ with respect to $x$, $D_{u}\theta=(\partial_{u_{j}}\theta_{i})_{i=1,\,\ldots,\,n}^{j=1,\,\ldots,\,m}$ is the Jacobian of $\theta$ with respect to $u$, $\operatorname{tr}(\cdot)$ is the trace of a square matrix, and $\hat{n}$ is the outward unit vector normal to the boundary.

Proof.

The “only if” part follows from the divergence theorem. To prove the “if” part, we proceed as in Lemma 1.1 of [Korda2018] and set $\theta=(u_{j}\tau(x)\delta_{ik})_{i=1,\,\ldots,\,n}$ in 2.3 for fixed $(j,k)\in\{1,\,\ldots,\,m\}\times\{1,\,\ldots,\,n\}$ and arbitrary $\tau\in C^{\infty}(\overline{X})$, where $\delta_{ik}$ is the usual Kronecker delta. Upon integrating by parts the term $\nabla_{x}\cdot\theta=u_{j}\partial_{x_{k}}\tau$ we obtain $\int_{X}(v_{jk}-\partial_{x_{k}}u_{j})\tau\,{\rm d}x=0$, which implies $v_{jk}=\partial_{x_{k}}u_{j}$ a.e. because $\tau$ is arbitrary. Repeating this argument for all pairs $(j,k)$ concludes the proof. ∎

Combining Proposition 1 with 2.1 and (2.2a,b) shows that

[TABLE]

where

[TABLE]

Solving (2.4) is clearly just as hard as solving the original variational problem (1.3). However, starting with (2.4) it is almost immediate to derive an infinite-dimensional convex program that proves a lower bound on $L^{*}$. First, we exchange the inf and sup at the expense of replacing equality with a lower bound. Second, we estimate the integrals of $F$ and $G_{i}$ in an elementary way using the constraints $u(x)\in Y$, $v(x)\in Z$ to obtain

[TABLE]

This is a convex program for $a,\,\phi,\,\psi_{1},\,\ldots,\,\psi_{s}$ and $\theta$ because $F$ and each $G_{i}$ depend affinely on them, so $(\alpha,\phi,\psi_{1},\,\ldots,\,\psi_{s},\theta)\mapsto\inf_{\Omega}F(x,y,z)$ and $(\alpha,\phi,\psi_{1},\,\ldots,\,\psi_{s},\theta)\mapsto\inf_{\Gamma_{i}}G_{i}(x,y)$ are convex functions. Convexity can be made more explicit by introducing slack variables $b$ and $c_{1},\,\ldots,\,c_{s}$ and rewriting

[TABLE]

The maximization problem on the righthand side is still hard to solve even with computer assistance. However, observe that any choice of $a,\,\phi,\,\psi_{1},\,\ldots,\,\psi_{s}$ and $\theta$ such that $F(x,y,z)$ and $G_{1}(x,y),\,\ldots,\,G_{s}(x,y)$ are bounded below on $\Omega$ and $\Gamma_{1},\,\ldots,\,\Gamma_{s}$ produces a lower bound on $L^{*}$. This makes it possible to prove suboptimal bounds analytically. Moreover, as we discuss next, in certain cases it is possible to optimize $a,\,\phi,\,\tilde{\psi}_{1},\,\ldots,\,\tilde{\psi}_{s}$ and $\theta$ numerically.

Remark 1*.*

To pass from 2.4 to 2.6 we have estimated $\int_{\partial X}G_{i}(x,u){\rm d}\sigma\geq\left|\partial X_{i}\right|\inf_{\Gamma_{i}}G_{i}(x,y)$ for each $i=1,\,\ldots,\,s$. As already observed in [Valmorbida2015b], we could improve these generic estimates by optimizing a lower bound on each $\int_{\partial X_{i}}G_{i}(x,u){\rm d}\sigma$ while optimizing $a$, $\phi$, $\psi_{i}$ and $\theta$. Indeed, minimizing $\int_{\partial X_{i}}G_{i}(x,u){\rm d}\sigma$ is a variational problem with affine dependence on $a$, $\phi$, $\psi_{i}$ and $\theta$ on an $(n-1)$-dimensional surface. Consequently, it can be relaxed into a convex program exactly as explained in this section if an explicit $(n-1)$-dimensional parametrization of $\partial X_{i}$ is available. This, in turn, requires estimates on $(n-2)$-dimensional integrals, and the procedure can be iterated until one is left with a one-dimensional problem. We do not pursue this approach here, but [Valmorbida2015b] gives a detailed discussion for square domains in $\mathbb{R}^{2}$.

2.2. Optimizing bounds by solving SDPs

Let us now restrict the attention to variational problems with polynomial data. Specifically, we assume that the functions $f$, $g$, $h_{i}$ and $l$ are polynomials and that the sets $X$, $Y$ and $Z$ are semialgebraic, i.e., they are defined by a finite number of polynomial equations and inequalities. This means that we can find polynomials $p_{1},\,\ldots,\,p_{r}$ such that

[TABLE]

Similarly, for each smooth portion $\partial X_{i}$ of the boundary there exist polynomials $q_{i,1},\,\ldots,\,q_{i,t_{i}}$ such that

[TABLE]

For simplicity, we assume that the outward unit vector $\hat{n}$ normal to each $\partial X_{i}$ is polynomial. This is true, for example, when $X$ is a polyhedral domain. Appendix A shows that some cases in which $\hat{n}$ is not polynomial can also be handled after a small modification of 2.7.

Given an integer $d$, let us restrict the optimization in 2.7 to degree-$d$ polynomials $\phi\in\mathbb{R}_{d}[x,y,z]$ and $\,\psi_{1},\,\ldots,\,\psi_{s},\,\theta_{1},\,\ldots,\,\theta_{n}\in\mathbb{R}_{d}[x,y]$. Then, the constraints are polynomial inequalities on semialgebraic sets and depend affinely on $a$, $b$, $c_{1},\,\ldots,\,c_{s}$ and the (finitely many) coefficients of $\phi$, $\psi_{1},\,\ldots,\,\psi_{s}$ and $\theta=(\theta_{1},\,\ldots,\,\theta_{n})$. These are NP-hard in general, but can be strengthened into tractable conditions by requiring that non-negative polynomials are representable as weighted sums of squares. More precisely, to the semialgebraic sets $\Omega$ and $\Gamma_{1},\ldots\Gamma_{s}$ we associate the sets of polynomials

[TABLE]

where $\mathbb{R}[x,y,z]$ (resp. $\mathbb{R}[x,y]$) is the space of polynomials in variables $x,y,z$ (resp. $x,y$) and $\Sigma[x,y,z]$ (resp. $\Sigma[x,y]$) is its subset of SOS polynomials. In other words, elements of $Q(\Omega)$ are weighted sums of $r+1$ SOS polynomials with weights $1,\,p_{1},\,\ldots,\,p_{r}$, and similarly for each $Q(\Gamma_{i})$. All polynomials in $Q(\Omega)$ and $Q(\Gamma_{i})$ are clearly non-negative on $\Omega$ and $\Gamma_{i}$, respectively, although the converse is not true in general. Then, we can replace the polynomial inequalities in 2.7 with weighted SOS constraints to obtain

[TABLE]

It is well known that optimization problems with weighted SOS constraints can be recast into SDPs (see, e.g., section 2.4.2 in [Lasserre2015]), and can therefore be solved using a variety of algorithms with polynomial-time complexity. In addition, while the bounds obtained with (2.11) for finite $d$ are typically strictly lower than $\mathbb{D}$, they converge to $\mathbb{D}$ as $d$ is raised provided that $\overline{X}$, $\overline{Y}$ and $\overline{Z}$ satisfy suitable compactness assumptions. The next proposition makes this statement precise and—as the results in section 4 imply—is the dual counterpart to Theorem 3 in [Korda2018] on the convergence of the moment-SDP relaxations described in the next section.

Proposition 2**.**

Suppose that $\overline{X}$, $\overline{Y}$ and $\overline{Z}$ are compact. Suppose also that there exist positive constants $C_{0},\,\ldots,\,C_{s}$ such that $C_{0}-\|x\|^{2}-\|y\|^{2}-\|z\|^{2}$ is in $Q(\Omega)$ and $C_{i}-\|x\|^{2}-\|y\|^{2}$ is in $Q(\Gamma_{i})$ for each $i=1,\,\ldots,\,s$. Then, the righthand side of 2.11 converges to $\mathbb{D}$ as $d\to\infty$.

Proof.

The proof follows a standard template in SOS optimization, which combines polynomial approximation of $C^{1}$ functions and their derivatives with Putinar’s Positivstellensatz [Putinar1993, Lemma 4.1] on the existence of weighted SOS representations for strictly positive polynomials on a class of compact semialgeraic sets that includes $\Omega$ and $\Gamma_{1},\,\ldots,\,\Gamma_{s}$ (see section 2.4.1 of [Lasserre2015] for more on this result). The details are left to the interested reader. ∎

3. Convex relaxation using measures

Korda et al. [Korda2018] proposed a different approach to bounding $L^{*}$ from below. The key idea is to relax a variational problem into an infinite-dimensional convex program over scaled probability measures. For variational problems with polynomial data, the optimization over measures can be replaced with SDPs that optimize finite sequences of their moments. Here we review this approach in the context of 1.3.

3.1. A convex program over non-negative measures

For each $u$ in the set $\mathcal{U}$ of admissible functions for 1.3, consider the maps

[TABLE]

Let $\mu^{u}:=\zeta^{u}\sharp{\rm d}x$ be the pushforward by $\zeta^{u}$ of the Lebesgue measure on $X$. Similarly, for each $i=1,\,\ldots,\,s$ let $\nu_{i}^{u}:=\chi_{i}^{u}\sharp{\rm d}\sigma$ be the pushforward by $\chi_{i}^{u}$ of the surface measure on $\partial X_{i}$. Following [Korda2018], we refer to $\mu^{u}$ and $\nu_{1}^{u},\,\ldots,\,\nu_{s}^{u}$ as the occupation and boundary measures of $u$. They are defined on $\Omega$ and $\Gamma_{1},\,\ldots,\,\Gamma_{s}$, respectively, and satisfy

[TABLE]

whenever $\eta$ and $\xi$ are such that the lefthand sides are well defined (see, e.g., [Bogachev2007, Theorem 3.6.1]). In particular, for each $u\in\mathcal{U}$ we have $L[u]=\langle l,\mu^{u}\rangle$ and we can rewrite 1.3 as a minimization problem over occupation measures:

[TABLE]

While the objective in this problem is linear in $\mu^{u}$, the set of occupation measures is generally not convex and 3.3 is no easier than 1.3. The strategy of Korda et al. [Korda2018] is to construct a convex set $\mathcal{M}$ of measures $(\mu,\nu_{1},\ldots,\nu_{s})$ that contains all occupation and boundary measures, and minimize $\langle l,\mu\rangle$ as a linear function over $\mathcal{M}$. Clearly, this is a convex program that yields a lower bound on $L^{*}$.

To construct $\mathcal{M}$, note that all occupation measures $\mu^{u}$ and boundary measures $\nu_{1}^{u},\,\ldots,\,\nu_{s}^{u}$ are non-negative, which we write as $\mu^{u},\nu_{1}^{u},\,\ldots,\,\nu_{s}^{u}\geq 0$. This follows immediately after fixing $\eta$ and $\xi$ in 3.2a and 3.2b to be the characteristic functions of any $E\subseteq X$ and $E\subseteq\partial X_{i}$, respectively. In particular, for $E=X$ and $E=\partial X_{i}$ we obtain

[TABLE]

This shows that occupation measures are scaled probability measures on $\Omega$ with mass $\left|X\right|$, while each boundary measure $\nu_{i}^{u}$ is a scaled probability measure on $\Gamma_{i}$ with mass $\left|\partial{X}_{i}\right|$.

Additional conditions on the occupation and boundary measures can be derived from the integral, differential and boundary constraints that define the set $\mathcal{U}$ of admissible functions for 1.3. Specifically, applying 3.2a to the integral constraint $\int_{X}f(x,u,Du)\,{\rm d}x=0$ gives

[TABLE]

In addition, multiplying the PDE constraint $g(x,u,Du)=0$ by $\phi(x,u,Du)$ for any $\phi\in C(\Omega)$, integrating the results over $X$, and using 3.2a we obtain

[TABLE]

Similarly, for each $i=1,\,\ldots,\,s$, multiplying the boundary constraint $h_{i}(x,u)=0$ by $\psi_{i}(x,u)$ for any $\psi_{i}\in C(\Gamma_{i})$, integrating over $\partial X_{i}$, and using 3.2b we conclude that

[TABLE]

The uncountably infinite sets of conditions 3.6 and 3.7 correspond to equations (14b,c) in [Korda2018]. In fact, 3.7 slightly generalizes (14c) in [Korda2018] because we do not sum over $i$, and we do not require that $\psi_{i}$ and $\psi_{j}$ match on $\partial X_{i}\cap\partial X_{j}$ whenever $\partial X_{i}$ is adjacent to $\partial X_{j}$.

The last set of constraints on occupation and boundary measures considered in [Korda2018] comes from the divergence theorem. Specifically, applying 3.2a and 3.2b to 2.3 with $v=Du$ after writing the boundary integral as the sum of integrals over $\partial X_{1},\,\ldots,\,\partial X_{s}$ gives, with notation analogous to Proposition 1,

[TABLE]

Combining all the above results we conclude that the occupation and boundary measures associated to feasible $u$ for 1.3 belong to the set

[TABLE]

It is easy to verify that $\mathcal{M}$ is convex, so in general not all its elements are tuples of occupation and boundary measures. Thus, minimizing $\langle l,\mu\rangle$ over $\mathcal{M}$ typically yields a lower bound on $L^{*}$:

[TABLE]

The minimization problem on the righthand side, which is a restatement of problem (16) from [Korda2018] in the present context, is an infinite-dimensional convex program. As for the lower bound in 2.7, solving this convex program is generally beyond the reach of analytical work. However, in stark constrast to 2.7, it does not suffice to find a feasible point for 3.10 because only the optimal value $\mathbb{P}$ is guaranteed to be a lower bound on $L^{*}$. Thus, one must either construct a minimizer (or minimizing sequence), or estimate $\mathbb{P}$ from below.

Remark 2*.*

As already discussed in [Korda2018], it is of great interest to determine conditions on the functions $f$, $g$, $h_{i}$ and the sets $X$, $Y$, $Z$ under which the set $\mathcal{M}$ is the closed convex hull of the set of occupation and boundary measures in the weak-$\ast$ topology. If this were the case, the inequality in 3.10 would in fact be an equality, and the convex relaxation described in this section would be tight. We will not consider this problem in this work.

3.2. Optimizing bounds using SDPs

As in section section 2.2, let us now restrict the attention to polynomial $f$, $g$, $h_{i}$ and $l$, and to semialgebraic sets $X$, $Y$ and $Z$. In this case, the measures $\mu,\nu_{1},\ldots,\nu_{s}$ are supported on semialgebraic sets and $\mathbb{P}$ can be estimated from below by solving finite-dimensional SDPs derived with the so-called moment-SDP relaxation procedure described in [Korda2018]. We only give a brief overview of this approach here, and refer the interested reader to section 4 in [Korda2018] for a full discussion.

Moment-SDP relaxations rely on two observations. The first is that one can relax the minimization in 3.10 by imposing 3.5, 3.6, 3.7 and 3.8 only over finitely many choices for $\phi$, $\psi_{i}$ and $\theta=(\theta_{1},\,\ldots,\,\theta_{n})$. The resulting problem is still infinite-dimensional, because the optimization variables are non-negative measures, but has a finite number of constraints. The second observation is that, when $f,\,g,\,h_{1},\,\ldots,\,h_{s}$ and $l$ are polynomials, the objective $\langle l,\mu\rangle$ in 3.10 is a finite linear combination of moments of $\mu$. Similarly, 3.5 and (3.4a,b) are linear equalities relating a finite number of moments of $\mu$ and $\nu_{1},\,\ldots,\,\nu_{s}$. More constraints on the moments can be obtained from 3.6, 3.7 and 3.8 by taking

[TABLE]

to be monomials of total degree no larger than some integer $d$ of choice. (In these expressions we have used standard multi-index notation, e.g., $x^{\alpha}=x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{x}}$ and $\left|\alpha\right|=\alpha_{1}+\cdots+\alpha_{n}$.) Thus, the minimization in 3.10 can be relaxed into a finite-dimensional minimization problem for finitely many moments of $\mu$ and $\nu_{1},\,\ldots,\,\nu_{s}$. This, in turn, can be relaxed into an SDP because sets of truncated sequences of moments of measures with semialgebraic supports admit semidefinite-representable outer approximations; see, for instance, [Korda2018, section 3.2] and the monographs [Laurent2009, Lasserre2015].

4. Duality

The functions $\phi$, $\psi_{1},\,\ldots,\,\psi_{s}$ and $\theta$ used in sections 2 and 3 are evidently very similar. This similarity is not an artefact of our notation, but stems from the fact that the convex program over Lagrange multipliers on the righthand side of 2.6 (equivalently, 2.7) is the Lagrangian dual of the convex program over measures on the righthand side of 3.10. This observation is made precise by the following result, which is our main contribution.

Theorem 1**.**

The convex programs on the righthand sides of 2.6 and 3.10 are weakly dual and $\mathbb{D}\leq\mathbb{P}$. The duality is strong, meaning that $\mathbb{P}=\mathbb{D}$, if the sets $Y$ and $Z$ are bounded.

Proof.

To prove weak duality, recall the definitions of $F$ and $G$ from (2.5a,b) and note that

[TABLE]

The last equality follows after observing that if $\{(x_{j},y_{j},z_{j})\}_{j\geq 1}\subset\Omega$ and $\{(x_{ij},y_{ij})\}_{j\geq 1}\subset\Gamma_{i}$ are minimizing sequences for $F$ and $G_{i}$, then the collections of scaled Dirac measures

[TABLE]

form a minimizing sequence for the inner infimum on the second line.

To prove strong duality when $Y$ and $Z$ are bounded observe that $\overline{Y}$ and $\overline{Z}$, hence $\Omega=\overline{X}\times\overline{Y}\times\overline{Z}$ and $\Gamma_{i}=\partial{X}_{i}\times\overline{Y}$, are compact. Then, an abstract minimax theorem due to Sion [Sion1958, Theorem 3.3] guarantees that equality is preserved when exchanging the inf and the sup in the second line of 4.1. To verify that the hypotheses of Sion’s theorem hold in our case, observe that the tuple $(a,\phi,\psi_{1},\ldots,\psi_{s},\theta)$ belongs to the product space

[TABLE]

We consider $M$ as a linear (hence, convex) space with the product topology generated by the usual norm topologies on $\mathbb{R}$, $C(\Omega)$, $C(\Gamma_{i})$ and $C^{1}(\overline{X}\times\overline{Y};\,\mathbb{R}^{n})$. Moreover, the space

[TABLE]

is the product of spaces of scaled probability measures on the compact sets $\Omega,\,\Gamma_{1},\,\ldots,\,\Gamma_{s}$ with mass $\left|X\right|,\,\left|\partial X_{1}\right|,\,\ldots,\,\left|\partial X_{s}\right|$, respectively. We consider $N$ in the product weak-$\ast$ topology, so it is a compact linear (hence, convex) space. Finally, the function

[TABLE]

is linear and continuous (hence, quasiconcave and upper semicontinuous) on $M$ for each $(\mu,\nu_{1},\ldots,\nu_{s})\in N$. Conversely, for each $(a,\phi,\psi_{1},\ldots,\psi_{s},\theta)\in M$ the function

[TABLE]

is linear and continuous (hence, quasiconvex and lower semicontinuous) on $N$. ∎

In light of the well known duality between the cones of weighted SOS polynomials and moment sequences [Laurent2009, Lasserre2015], Theorem 1 implies that the finite-dimensional SDP relaxations briefly described in sections 2.2 and 3.2 are also weakly dual when the same value of $d$ is taken in 2.11 and (3.11a,b,c). Strong duality at the level of SDPs holds under general constraint qualification conditions, such as Slater’s condition [Boyd2004, section 5.9.2], which are often satisfied in practice and can be verified independently of whether $\mathbb{P}=\mathbb{D}$. Moreover, strong duality is needed to guarantee convergence and good numerical performance of many primal-dual algorithms for semidefinite programming that solve the Karush–Kuhn–Tucker (KKT) optimality conditions. Since many commonly used SDP solvers implement such algorithms, from the point of view of numerical computations the two approaches to bounding $L^{*}$ described in sections 2 and 3 are equivalent.

5. Further comments

In this work we have demonstrated the duality between the two convex relaxation methods proposed by [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018] and [Korda2018] to bound from below the optimal value $L^{*}$ of a constrained variational problem. Precisely, the convex program over scaled probability measures formulated in [Korda2018] is weakly dual to the infinite-dimensional convex program over Lagrange multipliers that we have derived using ideas from [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018]. Furthermore, we have proven that the duality is strong when the range and gradients of admissible function in the original variational problem are constrained to bounded sets.

For problems with polynomial data, the hierarchies of SDPs obtained with the SOS and moment-SDP relaxations briefly described in sections 2.2 and 3.2 are also dual. Thus, they can be interpreted as extensions to polynomial variational problems of the well known dual hierarchies of SOS and moment-SDP relaxations for standard polynomial optimization problems. In contrast to standard polynomial optimization problems, however, there is currently no guarantee that the lower bounds on $L^{*}$ computed using these hierarchies can be made arbitrarily sharp. Numerical experiments from [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018] and [Korda2018] suggest that sharp bounds are sometimes possible, but the issue should be investigated in more depth both theoretically and computationally.

All our results can be extended to variational problems with more general choices for the set $\mathcal{U}$ in 1.2. For instance, the extension to problems with multiple PDE, boundary, and integral constraints is immediate. In addition, one can replace the Sobolev space $W^{1,\infty}(X;\,\mathbb{R}^{m})$ with $W^{1,p}(X;\,\mathbb{R}^{m})$, $1\leq p<\infty$, provided suitable conditions are imposed on the problem data $f,g,h_{i},l$ and the multipliers $\phi,\psi_{i},\theta$ to ensure that all integrals being considered are well defined. For instance, if $p=3$ and $g$ is quadratic in $u$ and $Du$, then $\phi$ should grow no faster than a linear function of $u$ and $Du$, so $\int_{X}\phi(x,u,Du)g(x,u,Du){\rm d}x$ is well defined. Finally, following [Korda2018] it is not difficult to adapt our discussion to variational problems with second-order semilinear PDE constraints and whose objective includes an integral over the boundary.

In the latter case, however, it does not seem possible to formulate finite-dimensional SDPs unless the unit vector normal to each smooth part of the boundary is either a polynomial, or a rational function with nonvanishing denominator. Fortunately, this is true for many domains encountered in applications, including polyhedral domains. To see the source of the difficulty, observe that when

[TABLE]

one only needs to redefine each $G_{i}$ in 2.4, 2.6 and 2.7 as

[TABLE]

If $l_{\rm b}$ is a given (nonzero) polynomial, then any nonpolynomial dependence of $\hat{n}$ cannot be absorbed by a judicious choice of $\psi_{i}$ as we have done in appendix A for the case $l_{\rm b}=0$. The only tractable situation is when $\hat{n}$ is a rational function with nonvanishing (hence, sign-definite) denominator. In this case, each inequality $G_{i}(x,u)-c_{i}\geq 0$ can be strengthened into a weighted SOS constraint after multiplying through by the denominator of $\hat{n}$. Similar considerations hold for the measure-theoretic approach of [Korda2018] and section 3.

Finally, we remark that although the convex relaxation methods of sections 2 and 3 are essentially equivalent for a large class of variational problems, each approach has unique advantages. We have already mentioned that if analytical bounds on $L^{*}$ are of interest, then it is more convenient to work with (2.6) because any suboptimal choice of Lagrange multipliers produces a valid bound, whereas the convex program over measures in 3.10 must be solved exactly. On the other hand, studying the relation between the feasible set $\mathcal{M}$ of 3.10 and the set of occupation and boundary measures may help to identify conditions under which arbitrarily sharp bounds on $L^{*}$ are possible. Moreover, for variational problems more general than 1.3 the formulation of a convex relaxation may be easier if Lagrange multipliers are used instead of measures, or viceversa. The former case includes constructing Lyapunov-like functionals for dynamical systems governed by PDEs, as originally done in [Valmorbida2014, Valmorbida2015a, Valmorbida2015b, Valmorbida2015c, Ahmadi2014, Ahmadi2015, Ahmadi2016, Ahmadi2017, Ahmadi2018], as well as parametric problems of the form

[TABLE]

where $\mathcal{P}$ is a convex set of parameters and the dependence of $l$ and $\mathcal{U}$ on $p$ is affine. Instead, it is simpler to use occupation and boundary measures when the objective $L[u]$ in 1.3 is not an integral functional, but can still be expressed as a function of finitely many moments that is either convex, or can be bounded from below using convex conditions. One such example is when $L[u]$ is a polynomial of integral functionals. All these extensions to the methods described in this work should be explored further in the future.

Appendix A SDPs with nonpolynomial unit normal vector

When the domain $X$ is a semialgebraic set, the smooth portions of its boundary $\partial X_{1},\,\ldots,\,\partial X_{s}$ are level sets of polynomials. Suppose that each $\partial X_{i}$ is defined by the equation $S_{i}(x)=0$, where $S_{i}$ is a polynomial such that $\|\nabla S_{i}\|\neq 0$ on $\partial X_{i}$. Then, the unit normal vector to $\partial X_{i}$ is given by

[TABLE]

When $\|\nabla S_{i}\|$ is not constant but polynomial, 2.7 can be relaxed into an SDP as described in section 2.2 after multiplying the inequality $G_{i}(x,y)-c_{i}\geq 0$ by $\|\nabla S_{i}\|$.

When $\|\nabla S_{i}\|$ is not polynomial, one can still formulate and SDP if 2.7 is modified as follows. Recall that the multiplier $\psi_{i}$, on which the function $G_{i}(x,u)$ in 2.4 depends, is arbitrary. Consequently, without loss of generality we can write

[TABLE]

for some continuous function $\hat{\psi}_{i}(x,u)$. Then, 2.4 can be rewritten as

[TABLE]

where $F$ is as in 2.5a and $\hat{G}_{i}(x,u):=\hat{\psi}_{i}(x,u)h_{i}(x,u)-\theta(x,u)\cdot\nabla S_{i}(x)$. Upon estimating the integrals we can replace 2.6 with

[TABLE]

and 2.7 with

[TABLE]

The constraints on the righthand side are polynomial inequalities if the problem data and $\phi$, $\hat{\psi}_{i}$, $\theta$ are polynomials, and can be strengthened into weighted SOS conditions as outlined in section 2.2. Then, lower bounds on $\hat{\mathbb{D}}$ (hence, $L^{*}$) can be computed by solving an SDP provided that all integrals $\int_{\partial X_{i}}\frac{{\rm d}\sigma}{\|\nabla S_{i}\|}$ can be computed analytically or approximated numerically.