Lebesgue and gaussian measure of unions of basic semi-algebraic sets

TL;DR

This paper introduces a systematic numerical method using semidefinite programming to approximate the measure of unions of semi-algebraic sets with arbitrary precision, leveraging available moments of the measure.

Contribution

It develops a hierarchy of semidefinite programs that converges to the measure of unions of semi-algebraic sets, enabling precise approximation from both above and below.

Findings

Convergent hierarchy of semidefinite programs for measure approximation.

Approximation of moments of the measure restricted to semi-algebraic sets.

Method applicable to Lebesgue measure with compact sets.

Abstract

Given a finite Borel measure $\mu$ on R n and basic semi-algebraic sets $\Omega$\_i $\subset$ R n , i = 1,. .. , p, we provide a systematic numerical scheme to approximate as closely as desired $\mu$(\cup\_i $\Omega$\_i), when all moments of $\mu$ are available (and finite). More precisely , we provide a hierarchy of semidefinite programs whose associated sequence of optimal values is monotone and converges to the desired value from above. The same methodology applied to the complement R n \ (\cup\_i $\Omega$\_i) provides a monotone sequence that converges to the desired value from below. When $\mu$ is the Lebesgue measure we assume that $\Omega$ := \cup\_i $\Omega$\_i is compact and contained in a known box B and in this case the complement is taken to be B \ $\Omega$. In fact, not only $\mu$($\Omega$) but also every finite vector of moments of $\mu$\_$\Omega$ (the restriction of $\mu$ on $\Omega$) can be approximated as closely as desired, and so permits to approximate the integral on $\Omega$ of any given polynomial.

Tables

Table 1. Table 1 . Example 1 : Values of ρ ¯ 10 subscript ¯ 𝜌 10 \overline{\rho}_{10} , ρ ¯ 10 subscript ¯ 𝜌 10 \underline{\rho}_{10} , ρ ¯ 10 S ​ t ​ o ​ k ​ e ​ s superscript subscript ¯ 𝜌 10 𝑆 𝑡 𝑜 𝑘 𝑒 𝑠 \overline{\rho}_{10}^{Stokes} , ρ ¯ 10 S ​ t ​ o ​ k ​ e ​ s superscript subscript ¯ 𝜌 10 𝑆 𝑡 𝑜 𝑘 𝑒 𝑠 \underline{\rho}_{10}^{Stokes} , ϵ 10 subscript italic-ϵ 10 \epsilon_{10} and ϵ 10 S ​ t ​ o ​ k ​ e ​ s superscript subscript italic-ϵ 10 𝑆 𝑡 𝑜 𝑘 𝑒 𝑠 \epsilon_{10}^{Stokes}

	$\text{u}=(0,0)$	$\text{u}=(0.1,0.5)$	$\text{u}=(0.5,0.5)$
${\overline{\rho }}_{10}$	$1.9649$	$1.9554$	$1.9484$
${\underset{¯}{\rho }}_{10}$	$1.6129$	$1.5752$	$1.5369$
${ϵ}_{10}$	$18\%$	$19\%$	$21\%$
${\overline{\rho }}_{10}^{Stokes}$	$1.8571$	$1.8308$	$1.8156$
${\underset{¯}{\rho }}_{10}^{Stokes}$	$1.7948$	$1.7746$	$1.7618$
${ϵ}_{10}^{Stokes}$	$3\%$	$3\%$	$3\%$

Table 2. Table 2 . Example 2 : Bounds and relative gap for d = 9 𝑑 9 d=9

${\overline{\rho }}_9$	${\underset{¯}{\rho }}_9$	${ϵ}_9$	${\overline{\rho }}_9^{Stokes}$	${\underset{¯}{\rho }}_9^{Stokes}$	${ϵ}_9^{Stokes}$
$2.0038$	$1.8252$	$8.9\%$	$1.9347$	$1.9019$	$1.7\%$

Table 3. Table 3 . Example 3 : Bounds and relative gap for d = 9 𝑑 9 d=9

${\overline{\rho }}_9$	${\underset{¯}{\rho }}_9$	${ϵ}_9$	${\overline{\rho }}_9^{Stokes}$	${\underset{¯}{\rho }}_9^{Stokes}$	${ϵ}_9^{Stokes}$
$2.0046$	$1.8342$	$8\%$	$1.9542$	$1.9083$	$2\%$

Table 4. Table 4 . Example 4 : Bounds and relative gap for d = 6 𝑑 6 d=6

${\overline{\rho }}_6$	${\underset{¯}{\rho }}_6$	${ϵ}_6$	${\overline{\rho }}_6^{Stokes}$	${\underset{¯}{\rho }}_6^{Stokes}$	${ϵ}_6^{Stokes}$
$2.8222$	$2.3123$	$18\%$	$2.6856$	$2.5360$	$5.6\%$

Table 5. Table 5 . Example 5 : Bounds and relative gaps for d = 7 𝑑 7 d=7

${\overline{\rho }}_7$	${\underset{¯}{\rho }}_7$	${ϵ}_7$	${\overline{\rho }}_7^{Stokes}$	${\underset{¯}{\rho }}_7^{Stokes}$	${ϵ}_7^{Stokes}$
$2.8143$	$2.3494$	$17\%$	$2.6887$	$2.5338$	$6\%$

Equations

$$\mu (\bigcup _{i=1}^p{\mathbf{\Omega }}_i)=\sum _{j=1}^p(-1{)}^{j+1}\sum _{1\le i_1<...<i_j\le p}\mu ({\mathbf{\Omega }}_{i_1}\cap ...\cap {\mathbf{\Omega }}_{i_j}),$$

$$Q(g_1,\dots ,g_m):=\left\{\sum _{j=0}^m{\sigma }_j{g}_j:{\sigma }_j\in \Sigma [\mathbf{x}]\right\}.$$

$$f\left(=\sum _{\gamma }{f}_{\gamma }{\mathbf{x}}^{\gamma }\right)\mapsto L_{\mathbf{y}}(f):=\sum _{\gamma }{f}_{\gamma }{y}_{\gamma },$$

$${\mathbf{M}}_d(\mathbf{y})(\alpha ,\beta ):=L_{\mathbf{y}}({\mathbf{x}}^{\alpha +\beta })=y_{\alpha +\beta },\forall \alpha ,\beta \in {\mathbb{N}}_d^n.$$

$${\mathbf{M}}_d(g\mathbf{y})(\alpha ,\beta ):=L_{\mathbf{y}}(g(\mathbf{x}){\mathbf{x}}^{\alpha +\beta })=\sum _{\gamma }{g}_{\gamma }{y}_{\alpha +\beta +\gamma },\forall \alpha ,\beta \in {\mathbb{N}}_d^n.$$

$$\displaystyle\mathbf{P}:\quad f^{*}=\underset{\phi}{\mbox{sup }}\Big{\{}\int_{\mathbf{K}}f\,d\phi:\lambda\leq\mu;\mbox{ }\phi\in\mathcal{M}(\mathbf{K})_{+}\Big{\}}$$

$$\mathbf{K}=\{\mathbf{x}\in {\mathbb{R}}^n:g_j(\mathbf{x})\ge 0,j=1,\dots ,m\},$$

$$z_{\alpha }={\int }_{\mathbf{B}}{\mathbf{x}}^{\alpha }d\mu (\mathbf{x}),\alpha \in {\mathbb{N}}^n,$$

$$\mathbf{Q}_{d}:\begin{array}[]{rl}\rho_{d}=\displaystyle\sup_{\mathbf{y}}&\{\,L_{\mathbf{y}}(f):\\ \mbox{s.t.}&\mathbf{M}_{d}(\mathbf{y})\,\succeq 0;\>\mathbf{M}_{d}(\mathbf{z}-\mathbf{y})\succeq 0\\ &\mathbf{M}_{d-r_{j}}(g_{j}\,\mathbf{y})\succeq 0,j=1,\ldots,m\}.\end{array}$$

$$\mathbf{Q}^{*}_{d}:\quad\rho^{*}_{d}=\displaystyle\inf_{p\in\mathbb{R}[\mathbf{x}]_{2d}}\,\{\,\int_{\mathbf{B}}p\,d\mu:\>p-f\geq 0\mbox{ on $\mathbf{K}$};\quad p\in\Sigma[\mathbf{x}]_{d}\,\},$$

$$\mathbf{\Omega }:=\bigcup _{i=1}^p{\mathbf{\Omega }}_i\subset \mathbf{B}.$$

$$S_k:=\sum _{1\le i_1<...<i_k\le p}\mu ({\mathbf{\Omega }}_{i_1}\cap ...\cap {\mathbf{\Omega }}_{i_k}),k=1,\dots ,p.$$

$$\mu (\bigcup _{k=1}^p{\mathbf{\Omega }}_k)=\sum _{k=1}^p(-1{)}^{k+1}{S}_k,$$

$$\mu (\bigcup _{i=1}^p{\mathbf{\Omega }}_i)$$

$$\ge \sum _{j=1}^{2k}(-1{)}^{j+1}{S}_j$$

$$\underset{d\to \infty }{\lim }\left(\sum _{i=1}^p(-1{)}^{k+1}\sum _{1\le i_1<\dots <i_k\le p}{\rho }_d^{i_1,\dots ,i_k}\right)=\mu (\mathbf{\Omega }).$$

$$\underset{d\to \infty }{\lim }\left(\sum _{i=1}^p(-1{)}^{k+1}\sum _{1\le i_1<\dots <i_k\le p}{y}_{d,0}^{i_1,\dots ,i_k}\right)=\mu (\mathbf{\Omega }),$$

$${\underline{\omega }}_d\le \mu (\mathbf{\Omega })\le {\overline{\omega }}_d,d\in \mathbb{N};\mu (\mathbf{\Omega })=\underset{d\to \infty }{\lim }{\underline{\omega }}_d=\underset{d\to \infty }{\lim }{\overline{\omega }}_d.$$

$${\mu }_{\alpha }={\int }_{\mathbb{B}}{\text{x}}^{\alpha }d\mu (\text{x}),\alpha \in {\mathbb{N}}^n,$$

$$\mathbf{Q}:\quad f^{*}=\sup_{\phi_{1},\ldots,\phi_{p}}\Big{\{}\sum_{i=1}^{p}\int_{\mathbf{\Omega}_{i}}fd\phi_{i}:\sum_{i=1}^{p}\phi_{i}\leq\mu;\mbox{ }\phi_{i}\in\mathcal{M}(\mathbf{\Omega}_{i})_{+},\>i=1,\ldots,p\Big{\}}.$$

$$f^{\ast }=\underset{\ell \to \infty }{\lim }\sum _{i=1}^p\int fd{\varphi }_i^{k_{\ell }}=\sum _{i=1}^p\underset{\ell \to \infty }{\lim }\int fd{\varphi }_i^{k_{\ell }}=\sum _{i=1}^p\int fd{\varphi }_i^{\ast },$$

$$\mathbf{x}\mapsto {\theta }_i(\mathbf{x})=\frac{1}{|\{j\in \{1,\dots ,p\}:\mathbf{x}\in {\mathbf{\Omega }}_j\}|}{1}_{{\mathbf{\Omega }}_i}(\mathbf{x}),\mathbf{x}\in \mathbf{\Omega }.$$

$${\varphi }_i^{\ast }(C):={\int }_C{\theta }_i(\mathbf{x})d\mu (\mathbf{x}),\forall C\in \mathcal{B}({\mathbb{R}}^n).$$

$$z_{\alpha }:=\int {\mathbf{x}}^{\alpha }d\mu (\mathbf{x}),\alpha \in {\mathbb{N}}^n.$$

$$\mathbf{Q}_{d}:\quad\begin{array}[]{rl}\rho^{f}_{d}=\displaystyle\sup_{\mathbf{y}^{1},\ldots,\mathbf{y}^{d}}&\Big{\{}\,\displaystyle\sum_{i=1}^{p}L_{\mathbf{y}^{i}}(f)\\ \mbox{s.t.}&\mathbf{M}_{d}(\mathbf{z}-\sum_{i=1}^{p}\mathbf{y}^{i})\succeq 0;\>\mathbf{M}_{d}(\mathbf{y}^{i})\succeq 0,\quad i=1,\ldots,p\\ &\mathbf{M}_{d-r_{ij}}(g_{ij}\,\mathbf{y}^{i})\succeq 0,\quad j=1,\ldots,m_{i};\>i=1,\ldots,p\,\Big{\}}.\end{array}$$

$$\rho^{f}_{d}\,\downarrow\,f^{*}\,=\,\int_{\mathbf{\Omega}}f\,d\mu,\quad\mbox{as $d\rightarrow\infty$.}$$

$$\underset{d\to \infty }{\lim }\sum _{i=1}^p{y}_{\alpha }^{i,d}=z_{\alpha }^{\ast }={\int }_{\mathbf{\Omega }}{\mathbf{x}}^{\alpha }d\mu .$$

$${\mathbf{M}}_d(\mathbf{z}-{\mathbf{y}}^{i,d})⪰{\mathbf{M}}_d(\sum _{j\ne i}{\mathbf{y}}^j)⪰0,i=1,\dots ,n.$$

$$\underset{k\to \infty }{\lim }{y}_{\alpha }^{i,d_k}=y_{\alpha }^{i,\ast },\forall \alpha \in {\mathbb{N}}^n,\forall i\in \{1,..,p\}.$$

$$f^{\ast }\le {\rho }_{d_k}^f=\sum _{i=1}^p{L}_{{\mathbf{y}}^{i,d_k}}(f)\downarrow \sum _{i=1}^p{L}_{{\mathbf{y}}^{i,\ast }}(f)=\sum _{i=1}^p\int fd{\varphi }_i.$$

$$\underset{k\to \infty }{\lim }\sum _{i=1}^p{y}_{\alpha }^{i,d_k}=\sum _{i=1}^p{y}_{\alpha }^{i,\ast }=z_{\alpha }^{\ast }={\int }_{\mathbf{\Omega }}{\mathbf{x}}^{\alpha }d\mu (\mathbf{x}).$$

Full text

Lebesgue and Gaussian measure of unions of basic semi-algebraic sets

Jean B. Lasserre

LAAS-CNRS and Institute of Mathematics

University of Toulouse

LAAS, 7 avenue du Colonel Roche

31077 Toulouse Cédex 4, France

lasserre@laas.fr

 and 

Youssouf Emin

Ecole Polytechnique

91 128 Palaiseau Cedex, France

youssouf.emin@polytechnique.edu

Jean B Lasserre: 7 Avenue du Colonel Roche, BP 54200, 31031 Toulouse cedex 4, France.

Tel: +66561336415; Fax: +33561336936; Email: lasserre@laas.fr

Youssouf Emin: Ecole Polytechnique, 91 128 Palaiseau, Cedex, France

Email: youssouf.emin@polytechnique.edu

Abstract

Given a finite Borel measure $\mu$ on $\mathbb{R}^{n}$ and basic semi-algebraic sets $\mathbf{\Omega}_{i}\subset\mathbb{R}^{n}$, $i=1,\ldots,p$, we provide a systematic numerical scheme to approximate as closely as desired $\mu(\bigcup_{i}\mathbf{\Omega}_{i})$, when all moments of $\mu$ are available (and finite). More precisely, we provide a hierarchy of semidefinite programs whose associated sequence of optimal values is monotone and converges to the desired value from above. The same methodology applied to the complement $\mathbb{R}^{n}\setminus(\bigcup_{i}\mathbf{\Omega}_{i})$ provides a monotone sequence that converges to the desired value from below. When $\mu$ is the Lebesgue measure we assume that $\mathbf{\Omega}:=\bigcup_{i}\mathbf{\Omega}_{i}$ is compact and contained in a known box $\mathbf{B}:=[-a,a]^{n}$ and in this case the complement is taken to be $\mathbf{B}\setminus\mathbf{\Omega}$. In fact, not only $\mu(\mathbf{\Omega})$ but also every finite vector of moments of $\mu_{\mathbf{\Omega}}$ (the restriction of $\mu$ on $\mathbf{\Omega}$) can be approximated as closely as desired, and so permits to approximate the integral on $\mathbf{\Omega}$ of any given polynomial.

Keywords: Lebesgue and Gaussian measure; semi-algebraic sets; moment problem and sums of squares; semidefinite programming; convex optimization

MSC: 44A60 28A75 90C05 90C22

1. Introduction

Given a set $\mathbf{\Omega}\subset\mathbb{R}^{n}$ and a finite Borel measure $\mu$ on $\mathbb{R}^{n}$, computing $\mu(\mathbf{\Omega})$ is a very challenging problem. In fact even approximating the Lebesgue volume of a convex body $\mathbf{\Omega}\subset\mathbb{R}^{n}$ (e.g. a polytope) is difficult; see e.g. Bollobás [2] and Dyer and Frieze [8]. However, in the latter case some efficient (non deterministic) algorithms with probabilistic guarantees are available and for more details the interested reader is referred to e.g. Dyer et al. [9], Cousins and Vempala [3, 4] and the references therein.

In the non convex case no such algorithm is available and one is left with approximating $\mu(\mathbf{\Omega})$ with Monte Carlo (or Quasi-Monte-Carlo) type methods as described in e.g. Niederreiter [17]. That is, one first generates a sample of $N$ points in $\mathbf{B}$ following the distribution $\mu$ on $\mathbf{B}$ and then one counts the number $N_{\mathbf{K}}$ of points that fall into $\mathbf{\Omega}$. This realization of the random variable $N_{\mathbf{K}}/N$ provides an estimate of $\mu(\mathbf{\Omega})$ but by no means an upper bound or a lower bound on $\mu(\mathbf{\Omega})$. Of course this method is quite fast, especially is small dimension.

Yet, as $\mu(\mathbf{\Omega})$ is indeed very difficult to compute exactly, a less ambitious but still useful goal would be to provide upper and/or lower bounds on $\mu(\mathbf{\Omega})$. Even better, a converging sequence of upper (or lower) bounds would be highly desirable. This is the strategy proposed in Henrion et al. [10] when $\mathbf{\Omega}$ is a compact basic semi-algebraic set and $\mu$ is the Lebesgue measure. In [10] the authors have provided a (deterministic) numerical scheme which yields a monotone sequence of upper bounds converging to $\mu(\mathbf{\Omega})$. It consists of solving a hierarchy of semidefinite programs of increasing size. By repeating the procedure but now with the complement $\mathbf{B}\setminus\mathbf{\Omega}$, one also obtains a monotone sequence of lower bounds converging to $\mu(\mathbf{\Omega})$. However, even on typical $2$ or $3$-dimensional examples, the convergence was rather slow and the authors proposed a slight modification which turned out to be much more efficient; the convergence was much faster but unfortunately not monotone anymore.

Contribution

The purpose of this paper is to introduce a deterministic method to approximate (in principle as closely as desired) the measure $\mu(\mathbf{\Omega})$ of the union $\mathbf{\Omega}=\bigcup_{i}\mathbf{\Omega}_{i}$ of finitely many basic semi-algebraic set. The finite Borel measure $\mu$ is any measure whose all moments are finite, e.g., the Lebesgue measure when $\mathbf{\Omega}$ is compact, the Gausssian measure $d\mu=\exp(-\|\mathbf{x}\|^{2})d\mathbf{x}$ for non-compact set $\mathbf{\Omega}$.

The method is similar in spirit to the one in [10] for a compact basic semi-algebraic set and the one in [13] for computing Gaussian measures of basic closed semi-algebraic sets (not necessarily compact), but with two important novelties.

$\bullet$ In contrast to [10] and [13], we consider a finite union $\mathbf{\Omega}$ of (non disjoint) basic semi-algebraic sets, which complicates matters significantly.

$\bullet$ We include a technique to accelerate the convergence different from the one described in [10]. Indeed in contrast to [10], it has the highly desirable feature to maintain the monotone convergence to $\mu(\mathbf{\Omega})$ which is essential if one wishes to obtain upper and lower bounds. It consists of using moments constraints coming from a particular application of Stokes’ theorem.

In fact this numerical scheme allows to approximate not only $\mu(\mathbf{\Omega})$ but also any fixed finite sequence of moments of the measure $\mu_{\mathbf{\Omega}}$ (where $\mu_{\mathbf{\Omega}}$ is the restriction of $\mu$ to $\mathbf{\Omega}$).

Remark 1.1**.**

One might invoke the inclusion-exclusion principle which states that

[TABLE]

so that in principle it suffices to compute (or approximate) $\mu(\mathbf{\Omega}_{i_{1}}\cap...\cap\mathbf{\Omega}_{i_{j}})$ for all possible intersections of the $\mathbf{\Omega}_{j}$’s, e.g. by the approach of [10] or [13]. But this approach has two major drawbacks. First there are possibly $2^{p}$ such sets and secondly, to compute an upper bound one has to compute an upper bound for such intersections with an odd number of elementary sets $\mathbf{\Omega}_{i_{j}}$, and a lower bound for such intersections with an even number of elementary sets. The latter lower bound in turn is obtained by computing an upper bound for the complement. This makes the whole procedure tedious and complicated. Finally, Bonferroni’s inequalities also provide a (finite) sequence of upper and lower bounds on $\mu(\mathbf{\Omega})$ but computing those bounds involves sums similar to the right-hand-side of (1.1), hence with the same drawbacks just mentioned. Our proposed technique is direct with no partial computation on intersections of elementary sets $\mathbf{\Omega}_{i_{j}}$.

Of course, the technique described in this paper is computationally expensive. In particular, its applicability is limited by the performance of the state-of-the-art semidefinite solvers because the size of the semidefinite programs increases fast with the rank in the hierarchy. Therefore it makes its application limited to small dimensional problems ($n\leq 3,4$). For higher dimensions only a few steps in the hierarchy can be implemented and therefore only upper and lower bounds (possibly crude) are expected. But the reader should keep in mind that the problem is very difficult and to the best of our knowledge we are not aware of an algorithm (at least at this level of generality) which provides certified upper and lower bounds with such convergence properties (even for convex sets and in particular for non compact sets $\mathbf{\Omega}$). In our opinion this methodology should be viewed as complementary to (rather than competing with) probabilistic methods.

2. Notation, definitions and preliminary results

2.1. Notation and definitions

Let $\mathbb{R}[\mathbf{x}]$ be the ring of polynomials in the variables $\mathbf{x}=(x_{1},\ldots,x_{n})$. Denote by $\mathbb{R}[\mathbf{x}]_{d}\subset\mathbb{R}[\mathbf{x}]$ the vector space of polynomials of degree at most $d$, which has dimension $s(d):=\binom{n+d}{d}$, with e.g., the usual canonical basis $(\mathbf{x}^{\gamma})_{\gamma\in\mathbb{N}^{n}_{d}}$ of monomials, where $\mathbb{N}^{n}_{d}:=\{\gamma\in\mathbb{N}^{n}\,:\,|\gamma|\leq d\}$, $\mathbb{N}$ is the set of natural numbers including [math] and $|\gamma|:=\sum_{i=1}^{n}\gamma_{i}$. Also, denote by $\Sigma[\mathbf{x}]\subset\mathbb{R}[\mathbf{x}]$ (resp. $\Sigma[\mathbf{x}]_{d}\subset\mathbb{R}[\mathbf{x}]_{2d}$) the cone of sums of squares (s.o.s.) polynomials (resp. s.o.s. polynomials of degree at most $2d$). If $f\in\mathbb{R}[\mathbf{x}]_{d}$, we write $f(\mathbf{x})=\sum_{\gamma\in\mathbb{N}^{n}_{d}}f_{\gamma}\mathbf{x}^{\gamma}$ in the canonical basis and denote by $\boldsymbol{f}=(f_{\gamma})_{\gamma}\in\mathbb{R}^{s(d)}$ its vector of coefficients. Finally, let $S^{n}$ denote the space of $n\times n$ real symmetric matrices, with inner product $\langle\mathbf{A},\mathbf{B}\rangle={\rm trace}\,\mathbf{A}\mathbf{B}$. We use the notation $\mathbf{A}\succeq 0$ (resp. $\mathbf{A}\succ 0$) to denote that $\mathbf{A}$ is positive semidefinite (definite). With $g_{0}:=1$, the quadratic module $Q(g_{1},\ldots,g_{m})\subset\mathbb{R}[\mathbf{x}]$ generated by polynomials $g_{1},\ldots,g_{m}$, is defined by:

[TABLE]

Definition 2.1** (Archimedean assumption).**

The quadratic module $Q(g_{1},\ldots,g_{m})$ is Archimedean if there exists $M>0$ such that the quadratic polynomial $\mathbf{x}\mapsto g_{m+1}:=M-\|\mathbf{x}\|^{2}$ belongs to $Q(g_{1},\ldots,g_{m})$. Notice that $g_{m+1}\in Q(g_{1},\ldots,g_{m})$ is an algebraic certificate that the set $\mathbf{K}:=\{\mathbf{x}:g_{j}(\mathbf{x})\geq 0,\>j=1,\ldots,m\}$ is compact.

If the set $\mathbf{K}:\{\mathbf{x}:g_{j}(\mathbf{x})\geq 0,\>j=1,\ldots,m\}$ is compact then $\|\mathbf{x}\|^{2}\leq M$ for some $M>0$, and one may always include the redundant quadratic constraint $\theta(\mathbf{x}):=M-\|\mathbf{x}\|^{2}\geq 0$ in the definition of $\mathbf{K}$ without changing $\mathbf{K}$. Then the quadratic module $Q(g_{1},\ldots,g_{m+1})$ is Archimidean.

Moment and localizing matrix

With a real sequence $\mathbf{y}=(y_{\gamma})_{\gamma\in\mathbb{N}^{n}_{d}}$, one may associate the (Riesz) linear functional $L_{\mathbf{y}}:\mathbb{R}[\mathbf{x}]_{d}\to\mathbb{R}$ defined by

[TABLE]

Denote by $\mathbf{M}_{d}(\mathbf{y})$ the moment matrix associated with $\mathbf{y}$, the real symmetric matrix with rows and columns indexed in the basis of monomials $(\mathbf{x}^{\gamma})_{\gamma\in\mathbb{N}^{n}_{d}}$, and with entries:

[TABLE]

Next, given $g\in\mathbb{R}[\mathbf{x}]$, denote by $\mathbf{M}_{d}(g\,\mathbf{y})$ the localizing moment matrix associated with $\mathbf{y}$ and $g$, the real symmetric matrix with rows and columns indexed in the basis of monomials $(\mathbf{x}^{\gamma})_{\gamma\in\mathbb{N}^{n}_{d}}$, and with entries:

[TABLE]

If $\mathbf{y}=(y_{\gamma})_{\gamma\in\mathbb{N}^{n}}$ is the sequence of moments of some Borel measure $\mu$ on $\mathbb{R}^{n}$ then $\mathbf{M}_{d}(\mathbf{y})\succeq 0$ for all $d\in\mathbb{N}$. However the converse is not true in general and it is related to the well-known fact that there are positive polynomials that are not sums of squares. Similarly, if the support of $\mu$ is contained in $\{\mathbf{x}:g(\mathbf{x})\geq 0\}$ then $\mathbf{M}_{d}(g\,\mathbf{y})\succeq 0$ for all $d\in\mathbb{N}$. For more details the interested reader is referred to e.g. [14, Chapter 3].

Given a Borel set $\mathbf{\Omega}\subset\mathbb{R}^{n}$ let $\mathcal{M}(\mathbf{\Omega})$ be the space of finite signed Borel measures on $\mathbf{\Omega}$ and let $\mathcal{M}(\mathbf{\Omega})_{+}\subset\mathcal{M}(\mathbf{\Omega})$ be the convex cone of finite Borel measures on $\mathbf{\Omega}$.

2.2. The measure of a basic semi-algebraic set

Let $\mathbf{B},\mathbf{K}\subset\mathbb{R}^{n}$ with $\mathbf{B}\supset\mathbf{K}$ and let $\mu$ be a finite Borel measure whose support is $\mathbf{B}$. (Typically $\mu$ is the Lebesgue measure on a box $\mathbf{B}$ and one wishes to compute the Lebesgue volume ${\rm vol}(\mathbf{K})$; alternatively $\mathbf{B}=\mathbb{R}^{n}$, $\mu$ is the Gaussian measure $d\mu=\exp(-\|\mathbf{x}\|^{2})d\mathbf{x}$ and one wishes to compute $\mu(\mathbf{K})$.)

An infinite-dimensional linear program $\mathbf{P}$

Let $f\in\mathbb{R}[\mathbf{x}]$ be positive almost everywhere on $\mathbf{K}$ and consider the following infinite-dimensional LP problem :

[TABLE]

Theorem 2.2** ([10]).**

The measure $\phi^{*}=\mu_{\mathbf{K}}$ (the restriction of $\mu$ to $\mathbf{K}$) is the unique optimal solution of $\mathbf{P}$. In particular, if $f(\mathbf{x})=1$ for all $\mathbf{x}$, then $f^{*}=\mu(\mathbf{K})$.

Semidefinite relaxations

Of course problem $\mathbf{P}$ in (2.2) is infinite-dimensional and cannot be solved directly. However, when $\mathbf{K}$ is a basic semi-algebraic set then Theorem 2.2 can be further exploited. So given $(g_{j})_{j=1}^{m}\subset\mathbb{R}[\mathbf{x}]$, let $\mathbf{K}\subset\mathbb{R}^{n}$ be the basic semi-algebraic set

[TABLE]

assumed to nonempty and compact. Let $\mathbf{B}\supset\mathbf{K}$ and let $\mu$ be a finite Borel measure whose all moments $\mathbf{z}=(z_{\alpha})$ with

[TABLE]

are available in closed form or can be computed.

To approximate $f^{*}$ as closely as desired in [10] the authors propose to solve the following hierarchy $(\mathbf{Q}_{d})_{d\in\mathbb{N}}$ of semidefinite programs111A semidefinite program (SDP) is a conic convex optimization problem with a remarkable modeling power. It can be solved efficiently (in time polynomial in its input size) up to arbitrary precision fixed in advance; see e.g. Anjos and Lasserre [1] indexed by $d\in\mathbb{N}$:

[TABLE]

Observe that $\mathbf{Q}_{d}$ is a relaxation of $P$, and so $\rho_{d}\geq\mu(\mathbf{K})$ for all $d$. In addition, the sequence $(\rho_{d})_{d\in\mathbb{N}}$ is monotone non increasing. The dual of (2.3) is the semidefinite program:

[TABLE]

and by weak duality, $\rho_{d}\leq\rho^{*}_{d}\leq f^{*}$ for all $d$.

Theorem 2.3** ([10]).**

Assume that $Q(g_{1},\ldots,g_{m})$ is Archimedean. Then $\rho_{d}\to f^{*}$ as $d\to\infty$. If $\mathbf{K}$ has nonempty interior then $\rho^{*}_{d}=\rho_{d}$ and (2.4) has an optimal solution $p^{*}\in\mathbb{R}[\mathbf{x}]_{2d}$.

So when $f=1$, $(\rho_{d})_{d\in\mathbb{N}}$ provides us with a monotone sequence of upper bounds on $f^{*}=\mu(\mathbf{K})$. Unfortunately the convergence is rather slow as observed on several numerical examples. This is because in the dual (2.4) the optimal solution $p^{*}\in\mathbb{R}[\mathbf{x}]_{2d}$ tries to approximate from above (in $L_{1}(\mathbf{B},\mu)$) the discontinuous function $1_{\mathbf{K}}$, which implies an annoying Gibb’s phenomenon222The Gibbs’ phenomenon appears at a jump discontinuity when one approximates a piecewise $C^{1}$ function with a continuous function, e.g. by its Fourier series.. To remedy this problem the authors in [10] propose to use a polynomial $f$, nonnegative on $\mathbf{K}$ and which vanishes on $\partial\mathbf{K}$. In this case the convergence $\rho_{d}\to\int_{\mathbf{K}}f\,d\mu$ as $d\to\infty$ is still monotone and if $\mathbf{y}^{d}=(y^{d}_{\alpha})_{\alpha\in\mathbb{N}^{n}_{2d}}$ denotes an optimal solution of (2.3) then $y^{d}_{0}\to\mu(\mathbf{K})$ as $d\to\infty$. However, while faster than with $f=1$, the latter convergence of $y^{d}_{0}$ to $\mu(\mathbf{K})$ is not monotone anymore, a rather annoying feature which prevents from obtaining a non increasing sequence of upper bounds.

3. Main result

The context

Let $\mathbf{B}\subset\mathbb{R}^{n}$ be a box, and for every $i=1,...,p$, let $\mathbf{\Omega}_{i}:=\{\,\mathbf{x}\in\mathbb{R}^{n}:g_{ij}(x)\geq 0,j=1,\ldots,m_{i}\}$, for some polynomials $(g_{ij})\subset\mathbb{R}[\mathbf{x}]$. Assume that $\mathbf{B}$ has been chosen so as to satisfy:

[TABLE]

The goal is to provide a numerical scheme to approximate as closely as desired the Lebesgue volume $\mu(\mathbf{K})$. (We will see how to adapt the methodology to also approximate as closely as desired $\mu(\mathbf{K})$ when $\mathbf{K}$ is not necessarily compact and $\mu$ is a Gaussian measure.) One possible approach described below is to use the powerful inclusion-exclusion principle and/or the associated Bonferroni inequalities.

3.1. The inclusion-exclusion principle and Bonferroni Inequalities

Let :

[TABLE]

By the inclusion-exclusion principle,

[TABLE]

which allows us to work with intersections of the $\mathbf{\Omega}_{k}$’s only. In addition, the Bonferroni inequalities state that

[TABLE]

which provides sequences of (increasingly tighter) upper and lower bounds.

Therefore to compute $\mu(\mathbf{\Omega})$ we only have to compute the measure of the intersection $\Theta_{i_{1},\ldots,i_{k}}:=\displaystyle\bigcap_{j=1,\ldots,k}\mathbf{\Omega}_{i_{j}}$, for all $1\leq i_{1}<...<i_{k}\leq p$. Notice that there are $2^{p}$ such sets. As each $\Theta_{i_{1},\ldots,i_{k}}\subset\mathbf{B}$ is a compact basic semi-algebraic set, one may apply the methodology described in §2.2, to obtain a sequence $(\rho^{(i_{1},\ldots,i_{k})}_{d})_{d\in\mathbb{N}}$ which converges to $\mu(\Theta_{i_{1},\ldots,i_{k}})$ as $d\to\infty$, and therefore

[TABLE]

Notice that the convergence is not monotone non increasing even if one solves (2.3) with $f=1$ because we sum up negative and positive terms. To maintain the monotone convergence (when $f=1$) it suffices to compute a lower bound on the complement $\mathbf{B}\setminus\Theta_{i_{1},\ldots,i_{k}}$ when $k$ is even. However as already mentioned the convergence is expected to be rather slow.

To accelerate the convergence one may use $f=\prod_{j=1}^{k}\prod_{\ell=1}^{m_{i_{j}}}g_{i_{j}\ell}$ when one solves (2.3) with $\mathbf{\Omega}=\Theta_{i_{1},\ldots,i_{k}}$ as $f\geq 0$ on $\Theta_{i_{1},\ldots,i_{k}}$ and $f=0$ on $\partial\Theta_{i_{1},\ldots,i_{k}}$. But then the convergence

[TABLE]

(where $\mathbf{y}^{i_{1},\ldots,i_{k}}_{d}=y^{i_{1},\ldots,i_{k}}_{d,\alpha}$ is an optimal solution of (2.3) with $\mathbf{\Omega}=\Theta_{i_{1},\ldots,i_{k}}$) is not monotone anymore.

3.2. A direct approach

In this section we describe a direct approach with two distinguishing features:

• •

It does not use the inclusion-exclusion principle and the need to approximate $\mu(\bigcap_{j=1}^{k}\mathbf{\Omega}_{i_{j}})$ for all $2^{p}$ such sets.

• •

The convergence to $\mu(\mathbf{\Omega})$ (and also to $\mu(\mathbf{B}\setminus\mathbf{\Omega}))$ is monotone non increasing, that is, we can compute two sequences $(\overline{\omega}_{d})_{d\in\mathbb{N}}$ and $(\underline{\omega}_{d})_{d\in\mathbb{N}}$ such that:

[TABLE]

Recall that any finite number of moments

[TABLE]

are either available in closed-form or can be obtained numerically.

A multi infinite-dimensional linear program $\mathbf{Q}$

As in §2.2 we first introduce an infinite-dimensional LP problem $\mathbf{Q}$ whose unique optimal solution is the restriction of $\mu$ on $\mathbf{\Omega}$ (and whose dual has a clear interpretation).

Let $f\in\mathbb{R}[\mathbf{x}]$ be positive almost everywhere on $\mathbf{\Omega}$ and consider the following infinite-dimensional LP problem :

[TABLE]

Theorem 3.1**.**

Problem $\mathbf{Q}$ has an optimal solution $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ and every optimal solution satisfies $\sum_{i=1}^{p}\phi^{*}_{i}=\mu_{\mathbf{\Omega}}$, where $\mu_{\mathbf{\Omega}}$ is the restriction of $\mu$ to $\mathbf{\Omega}$.

Proof.

We first prove that $\mathbf{Q}$ has an optimal solution. The set $\Delta_{\mu}:=\{\phi\in\mathcal{M}(\mathbb{R}^{n})_{+}:\phi\leq\,\mu\}$ is weakly sequentially compact; see e.g. Dunford & Schwartz [7, Theorem 1, p. 305]. Therefore let $(\phi^{k}_{1},\ldots,\phi^{k}_{p})_{k\in\mathbb{N}}$ be a maximizing sequence of feasible solutions of $\mathbf{Q}$. There exists a subsequence $(k_{\ell})_{\ell\in\mathbb{N}}$ such that for every $i=1,\ldots,p$, $\phi^{k_{\ell}}_{i}\stackrel{{\scriptstyle w}}{{\to}}\phi^{*}_{i}$ for some $\phi^{*}_{i}\in\mathcal{M}(\mathbb{R}^{n})_{+}$. The above weak convergence and $\displaystyle\int_{\mathbf{\Omega}_{i}^{c}}\,d\phi_{i}^{k_{\ell}}=0$ implies $\displaystyle\int_{\mathbf{\Omega}_{i}^{c}}\,d\phi_{i}^{*}=0$, that is, $\phi^{*}_{i}\in\mathcal{M}(\mathbf{\Omega}_{i})_{+}$ for all $i=1,\ldots,p$. Weak convergence again implies $\sum_{i=1}\phi^{*}_{i}\leq\mu$ and so $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ is a feasible solution of $\mathbf{Q}$. Finally weak convergence also implies

[TABLE]

which proves that $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ is an optimal solution of $\mathbf{Q}$.

We next prove that $f^{*}=\int fd\mu_{\mathbf{\Omega}}$. Indeed, firstly observe that for every feasible solution $(\phi_{1},\ldots,\phi_{p})$ of $\mathbf{Q}$, $\sum_{i=1}^{p}\int fd\phi_{i}\leq\int_{\mathbf{\Omega}}fd\mu=\int fd\mu_{\mathbf{\Omega}}$. On the other hand, for every $i=1,\ldots,p$, denote by $\theta_{i}$ the measurable function defined on $\mathbf{\Omega}$ by :

[TABLE]

The (discontinuous) functions $(\theta_{i})_{i=1,\ldots,p}$ form a partition of unity subordinate to the open cover $\bigcup_{i}{\rm int}(\mathbf{\Omega}_{i})$. For every $i=1,\ldots,p$, let $\phi^{*}_{i}\in\mathcal{M}(\mathbf{\Omega}_{i})_{+}$ be the finite Borel measure defined by:

[TABLE]

Hence, $\sum_{i=1}^{p}\phi^{*}_{i}(C)=\int_{C}\sum_{i=1}^{p}\theta_{i}(\mathbf{x})d\mu(\mathbf{x})=\int_{C}1_{\mathbf{\Omega}}(\mathbf{x})(\mathbf{x})d\mu(\mathbf{x})=\mu_{\mathbf{\Omega}}(C)\leq\mu(C)$. Therefore $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ is a feasible solution of $\mathbf{Q}$ such that $\sum_{i=1}^{p}\phi^{*}_{i}=\mu_{\mathbf{\Omega}}$, and so $\sum_{i=1}^{p}\int fd\phi^{*}_{i}=f^{*}$, i.e., $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ is an optimal solution of $\mathbf{Q}$. In fact, every optimal solution $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ of $\mathbf{Q}$ satisfies $\sum_{i=1}^{p}\int fd\phi^{*}_{i}=f^{*}$, and therefore $\phi^{*}:=\sum_{i=1}^{p}\phi^{*}_{i}\in\mathcal{M}(\mathbf{\Omega})_{+}$ is an optimal solution of $\sup_{\phi}\{\int fd\phi:\phi\leq\mu;\>\phi\in\mathcal{M}(\mathbf{\Omega})_{+}\}$. By Theorem 2.2 this solution $\phi^{*}$ is unique, which yields the desired result. ∎

A hierarchy of semidefinite relaxations

Let $\mathbf{z}=(z_{\alpha})_{\alpha\in\mathbb{N}^{n}}$ be the sequence of all moments of $\mu$ , that is,

[TABLE]

Let $\mathbf{B}\subset\mathbb{R}^{n}$ be a box and $\mathbf{\Omega}\subset\mathbf{B}$ be a compact semi-algebraic as in (3.1). With no loss of generality and possibly after scaling, we may and will assume that $\mathbf{B}\subset[-1,1]^{n}$ and $\mu$ is a probability measure. Therefore $|z_{\alpha}|\leq 1$ for all $\alpha\in\mathbb{N}^{n}$.

Let $r_{ij}=\lceil{\rm deg}(g_{ij})/2\rceil$ and let $f\in\mathbb{R}[\mathbf{x}]$ be a given polynomial positive almost everywhere on $\mathbf{\Omega}$ (and define $r_{00}:=\lceil{\rm deg}(f)/2\rceil$). For $d\geq d_{0}:=\mbox{max}_{i,j}r_{ij}$, consider the following hierarchy of semidefinite programs $(\mathbf{Q}_{d})$ indexed by $d\in\mathbb{N}$ :

[TABLE]

Observe that $\rho^{f}_{d}\geq f^{*}$ for all $d\in\mathbb{N}$. Indeed, if $(\mathbf{z}^{1},\ldots,\mathbf{z}^{p})$ is the sequence of moments of an optimal solution $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ of $\mathbf{Q}$ in (3.2) then $(\mathbf{z}^{1},\ldots,\mathbf{z}^{p})$ is also a feasible solution of $\mathbf{Q}_{d}$.

Theorem 3.2**.**

Consider the semidefinite programs $(\mathbf{Q}_{d})$, $d\geq d_{0}$. Then :

(i) $\mathbf{Q}_{d}$ has an optimal solution and the associated sequence of optimal values $(\rho^{f}_{d})_{d\in\mathbb{N}}$ is monotone non increasing and converges to $f^{*}$, that is:

[TABLE]

(ii) Let $(\mathbf{y}^{1,d},\ldots,\mathbf{y}^{p,d})$ be an optimal solution of $\mathbf{Q}_{d}$. Then for each $\alpha\in\mathbb{N}^{n}$:

[TABLE]

and in particular $\displaystyle\lim_{d\to\infty}\sum_{i=1}^{p}y^{i,d}_{0}=\mu(\mathbf{\Omega})$.

Proof.

For a sequence $\mathbf{y}=(y_{\alpha})$, let $\tau_{d}(\mathbf{y})=\max_{i=1,\ldots,n}L_{\mathbf{y}}(x_{i}^{2d})$ and recall that if $\mathbf{M}_{d}(\mathbf{y})\succeq 0$ then $|y_{\alpha}|\leq\max[y_{0},\max_{i}L_{\mathbf{y}}(x_{i}^{2d})]$ for every $\alpha\in\mathbb{N}^{n}_{2d}$; see [14, Proposition 3.6]. Next, observe that from $\mathbf{M}_{d}(\mathbf{z}-\sum_{i=1}^{p}\mathbf{y}^{i,d})\succeq 0$ and $\mathbf{M}_{d}(\mathbf{y}^{i})\succeq 0$,

[TABLE]

Hence the diagonal elements $z_{2\alpha}-y^{i,d}_{2\alpha}$ are all nonnegative which in turn implies $\tau_{d}(\mathbf{y}^{i,d})\leq\tau_{d}(\mathbf{z})\leq 1$ for all $i=1,\ldots,n$. As $z_{0}=1$ then by [14, Proposition 3.6 ] $|y^{i,d}_{\alpha}|\leq 1$ for every $\alpha\in\mathbb{N}^{n}_{2d}$, and so the feasible set of semidefinite program $\mathbf{Q}_{d}$ is closed, bounded, hence compact, and therefore $\mathbf{Q}_{d}$ has an optimal solution.

Next, let $(\mathbf{y}^{1,d},\ldots,\mathbf{y}^{p,d})$ be an optimal solution of $\mathbf{Q}_{d}$ and by completing with zeros, make $(\mathbf{y}^{1,d},\ldots,\mathbf{y}^{p,d})$ an element of the unit ball of $(\ell_{\infty})^{p}$ (where $\ell_{\infty}$ is the Banach space of bounded sequences, equipped with the sup-norm). As $(\ell_{\infty})^{p}$ is the topological dual of $(\ell_{1})^{p}$, by the Banach-Alaoglu Theorem, there exists $(\mathbf{y}^{1,*},..,\mathbf{y}^{p,*})\in(\ell_{\infty})^{p}$ and a subsequence $\{d_{k}\}$ such that $(\mathbf{y}^{1,d_{k}},\ldots,\mathbf{y}^{p,d_{k}})\rightarrow(\mathbf{y}^{1,*},\ldots,\mathbf{y}^{p,*})$ as $k\rightarrow\infty$, for the weak $\star$ topology $\sigma((\ell_{\infty})^{p},(\ell_{1})^{p})$. In particular,

[TABLE]

Next let $d\in\mathbb{N}$ be fixed arbitrary. From the pointwise convergence (3.7) we also obtain $\mathbf{M}_{d}(\mathbf{y}^{i,*})\succeq 0$ and $\mathbf{M}_{d}(\mathbf{z}-\sum_{i=1}^{p}\mathbf{y}^{i,*})\succeq 0$ for every $i=1,\ldots,p$. Similary, $\mathbf{M}_{d-r_{ij}}(g_{ij}\mathbf{y}^{i,*})\succeq 0$ for every $i$ and $j$. As $d$ was arbitrary, by Putinar’s Positivistellensatz [18], $\mathbf{y}^{i,*}$ has a representing measure $\phi_{i}$ supported on $\mathbf{\Omega}_{i}$ for all $i=1,\ldots,p$, and $\sum_{i=1}^{p}\phi_{i}\leq\mu$. In particular from (3.7), as $k\to\infty$,

[TABLE]

Therefore $(\phi_{1},\ldots,\phi_{p})$ is admissible for problem $\mathbf{Q}$ with value $\sum_{i=1}^{p}\int fd\phi_{i}\geq f^{*}$, and so $(\phi_{1},\ldots,\phi_{p})$ is an optimal solution of $\mathbf{Q}$. Finally, by Theorem 3.1, $\sum_{i=1}^{p}\phi_{i}=\mu_{\mathbf{\Omega}}$. And so for each $\alpha\in\mathbb{N}^{n}$:

[TABLE]

As the converging subsequence $(d_{k})_{k\in\mathbb{N}}$ was arbitrary, it follows that in fact the whole sequence $(\sum_{i=1}^{p}y^{i,d}_{\alpha})_{d}$ converges to $z_{\alpha}$, for all $\alpha\in\mathbb{N}^{n}$, that is, (3.6) holds. ∎

The dual of $\mathbf{Q}_{d}$

Let $g_{i0}(\mathbf{x})=1$ for all $\mathbf{x}\in\mathbb{R}^{n}$, $i=1,\ldots,p$. The dual of the semidefinite program $\mathbf{Q}_{d}$ is the semidefinite program:

[TABLE]

Proposition 3.3**.**

Assume that for every $i=1,\ldots,p$, both $\mathbf{\Omega}_{i}$ and $\mathbf{B}\setminus\mathbf{\Omega}_{i}$ have nonempty interior. Then there is no duality gap between (3.5) and its dual (3.8), that is, $\rho^{f}_{d}=(\rho^{f}_{d})^{*}$ for all $d\geq d_{0}$. Moreover (3.8) has an optimal solution $(q^{*},(\sigma_{ij}^{*})$.

Proof.

Let $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ be the measures defined in (3.3) the proof of the Theorem 3.1 and let $\mathbf{y}^{i}_{d}$ be the sequence of moments up to degree $d$ of $\phi^{*}_{i}$, $i=1,\ldots,p$. As every $\mathbf{\Omega}_{i}$ has nonempty interior, then clearly $\mathbf{M}_{d}(\mathbf{y}^{i})\succ 0$ and $\mathbf{M}_{d-r_{ij}}(g_{ij}\mathbf{y}^{i})\succ 0$ for every $j=1,\ldots,m_{i}$ and $i=1,\ldots,p$. As $\mathbf{B}\setminus\mathbf{\Omega}$ also has nonempty interior then $\mathbf{M}_{d}(\mathbf{z}-\sum_{i=1}^{p}\mathbf{y}^{i})\succ 0$. Therefore Slater’s condition holds for $\mathbf{Q}_{d}$. In addition, the set of admissible solution of $\mathbf{Q}_{d}^{*}$ is nonempty (set $q=f$ and $\sigma_{ij}=0$ for all $i,j$), and therefore a standard result in conic convex optimization yields the desired result333In fact as the set of optimal solutions of (3.5) is compact, the absence of a duality gap between (3.5) and (3.8) also follows from [19] without the conditions ${\rm int}(\mathbf{\Omega}_{i})\neq\emptyset$ and ${\rm int}(\mathbf{B}\setminus\mathbf{\Omega}_{i})\neq\emptyset$..

∎

As in the case of a basic closed semi-algebraic set, when $f$ is the constant function $1$ the convergence $\rho^{f}_{d}\to f^{*}=\mu(\mathbf{\Omega})$ is monotone non increasing, a highly desirable feature. However in typical examples this convergence is rather slow. Again one may take for $f$ a function that is nonnegative on $\mathbf{\Omega}$ and which vanishes on $\partial\mathbf{\Omega}$. This accelerates the convergence both $\rho^{f}_{d}\to f^{*}$ and $\sum_{i}y^{id}_{0}\to\mu(\mathbf{\Omega})$ as $d\to\infty$, but if by construction the former is monotone non increasing, the latter is not monotone anymore, a rather annoying feature if the goal is to obtain a converging sequence of upper bounds. In the next section we describe a technique that allows to accelerate significantly the convergence $\sum_{i}y^{id}_{0}\to\mu(\mathbf{\Omega})$ as $d\to\infty$, while maintaining its monotone non increasing character.

3.3. Convergence improvement using Stokes’ formula

In this section we show how to improve significantly the monotone non increasing convergence of $\rho_{d}^{1}$ (i.e. $\rho^{f}_{d}$ with $f=1$) to $\mu(\mathbf{\Omega})$. To do this we will use Stokes’ theorem for integration and in the sequel, to avoid technicalities we assume that $\mathbf{\Omega}\subset\mathbb{R}^{n}$ is the closure of its interior, i.e., $\mathbf{\Omega}=\overline{{\rm int}(\mathbf{\Omega})}$. The basic idea is simple to express in informal terms.

Since we know in advance that $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ in (3.3) is an optimal solution of problem $\mathbf{Q}$, every additional information in terms of linear constraints on the moments of $\phi^{*}_{i}$ can be included in $\mathbf{Q}$ without changing its optimal value. BUT when included in the relaxation $\mathbf{Q}_{d}$ it will provide useful additional constraints that restrict the feasible set of $\mathbf{Q}_{d}$ and so make its optimal value necessarily smaller.

Suppose for the moment that $\mathbf{\Omega}$ is compact with smooth boundary, and assume that the measure $\mu$ has a density $h$ with respect to Lebesgue measure $d\mathbf{x}$ of the form $q(\mathbf{x})\exp(r(\mathbf{x}))1_{\mathbf{B}}(\mathbf{x})$ for some polynomial $r,q\in\mathbb{R}[\mathbf{x}]$. Let $X$ be a given vector field and $f\in\mathbb{R}[\mathbf{x}]$. Then Stokes’ theorem states:

[TABLE]

where $\vec{n}_{\mathbf{x}}$ is outward pointing normal at $\mathbf{x}\in\partial\mathbf{\Omega}$, and $\sigma$ is the $(n-1)$-dimensional Hausdorff measure on $\partial\mathbf{\Omega}$. In particular if $f$ vanishes on $\partial\mathbf{\Omega}$ and $\mathbf{X}=e_{k}\in\mathbb{R}^{n}$ (where $e_{k}(j)=\delta_{k=j}$) Stokes’ formula becomes

[TABLE]

To exploit (3.9) in our particular context where $\mathbf{\Omega}$ is defined in (2.2), let $g=\prod_{i=1}^{p}\prod_{j=1}^{m_{i}}g_{ij}$ and let $\mathbf{x}\mapsto f(\mathbf{x})=\mathbf{x}^{\alpha}\,g(\mathbf{x})q(\mathbf{x})$ with $\alpha\in\mathbb{N}^{n}$ arbitrary. Then $f$ vanishes on $\partial\mathbf{\Omega}$ and on $\partial\mathbf{\Omega}_{i_{1},\ldots,i_{s}}:=\mathbf{\Omega}_{i_{1}}\cap\cdots\cap\mathbf{\Omega}_{i_{s}}$ for all $1\leq i_{1}<\ldots<i_{s}\leq p$, $s=1,\ldots,p$. Hence by (3.9):

[TABLE]

for all $1\leq i_{1}<\ldots<i_{s}\leq p$, $s=1,\ldots,p$, where

[TABLE]

Recalling how $\phi^{*}_{i}$ is defined in (3.3), it can be written as

[TABLE]

where each $\phi^{*}_{i,i_{1},\ldots,i_{s}}$ is supported on $\mathbf{\Omega}_{i}\cap\mathbf{\Omega}_{i_{1},\ldots,i_{s}}$ and has a constant density w.r.t. $\mu$. Therefore, for every $i=1,\ldots,p$:

[TABLE]

Hence (3.12) provides additional useful information on the optimal solution $(\phi^{*}_{1},\ldots,\phi^{*}_{p})$ of $\mathbf{Q}$ defined in (3.3). Namely it translates into

[TABLE]

i.e., linear constraints on the moments of $\phi^{*}_{i}$, for every $i=1,\ldots,p$.

Plugging this additional linear constraints on the moments of $\phi^{*}_{i}$ into the relaxation $\mathbf{Q}_{d}$, yields the following new hierarchy of SDP-relaxation $(\mathbf{Q}_{d}^{{\rm stokes}})_{d\geq d_{0}}$:

[TABLE]

By construction $\rho^{1}_{d}\geq\rho^{{\rm Stokes}}_{d}\geq\mu(\mathbf{\Omega})$ holds for every $d\geq d_{0}$, and the analogue of Theorem 3.2 (with $f=1$) reads:

Theorem 3.4**.**

Consider the semidefinite programs $(\mathbf{Q}_{d}^{{\rm Stokes}})$, $d\geq d_{0}$, defined in (3.13). Then :

(i) $\mathbf{Q}_{d}^{{\rm Stokes}}$ has an optimal solution and the associated sequence of optimal values $(\rho^{{\rm Stokes}}_{d})_{d\in\mathbb{N}}$ is monotone non increasing and converges to $\mu(\mathbf{\Omega})$, that is:

[TABLE]

(ii) Let $(\mathbf{y}^{1,d},\ldots,\mathbf{y}^{p,d})$ be an optimal solution of $\mathbf{Q}^{{\rm Stokes}}_{d}$. Then for each $\alpha\in\mathbb{N}^{n}$:

[TABLE]

The proof being almost a verbatim copy of that of Theorem 3.2, is omitted.

The important feature of Theorem 3.4 is that we now have the monotone non increasing convergence $\rho^{{\rm Stokes}}_{d}\downarrow\mu(\mathbf{\Omega})$ (compare with (3.6) (with $\alpha=0$) in Theorem 3.2).

3.4. Gaussian measure of non compact sets $\mathbf{\Omega}$

So far Theorem 3.2 and Theorem 3.4 have been given for $\mu$ supported on a box $\mathbf{B}$, and so only for sets $\mathbf{\Omega}$ in (3.1) that are compact.

It turns out that for a Gaussian measure $\mu$ of (possibly non-compact) sets $\mathbf{\Omega}=\bigcup_{i}\mathbf{\Omega}_{i}$, Theorem 3.2 (resp. Theorem 3.4) is still valid with exactly the same statement and exactly the same semidefinite relaxations (3.5) (resp. (3.13)), except that now $\mathbf{z}=(z_{\alpha})$ is the vector of moments of the Gaussian measure $\mu$ (instead of the moments of the Lebesgue measure on $\mathbf{B}$ previously).

However in the gaussian case the proof of Theorem 3.2(i)-(ii) and Theorem 3.4(i)-(ii) uses quite different arguments (some already used in [13] for a basic semi-algebraic set). Indeed as $\mathbf{\Omega}$ is not necessarily compact :

• The uniform bound $\sup_{\alpha}|\mathbf{z}_{\alpha}|\leq 1$ is not valid any more for the relaxations $\mathbf{Q}_{d}$ and $\mathbf{Q}^{Stokes}_{d}$.

• One cannot invoke Putinar’s Positivstellensatz [18] any more.

• The standard version of Stokes’ theorem where $\mathbf{\Omega}$ is compact cannot be invoked anymore either.

The new arguments that we need are the following:

$\bullet$ A crucial fact is that $\mu$ satisfies Carleman’s condition

[TABLE]

Then a sequence $\mathbf{y}=(y_{\alpha})_{\alpha\in\mathbb{N}^{n}}$ such that $\mathbf{M}_{d}(\mathbf{y})\succeq 0$ for all $d\in\mathbb{N}$, and

[TABLE]

has a unique representing measure $\phi$ on $\mathbb{R}^{n}$ which is moment determinate; see for instance [14, Proposition 3.5, p. 60].

$\bullet$ If in addition $\mathbf{M}_{d}(h\,\mathbf{y})\succeq 0$ for all $d\in\mathbb{N}$ (where $h\in\mathbb{R}[\mathbf{x}]$), and as $\phi$ satisfies (3.15), then $h(\mathbf{x})\geq 0$ for all $\mathbf{x}$ in the support of $\phi$; see Lasserre [15]. This argument is used to show that $\phi$ is supported on $\mathbf{\Omega}$.

$\bullet$ To obtain a version of Stokes for non-compact set $\mathbf{\Omega}$ with boundary $\partial\mathbf{\Omega}$, we invoke a limiting argument that uses (the standard) Stokes’s theorem on the compact $\mathbf{\Omega}\cap\mathbf{B}(0,M)$ (where $\mathbf{B}(0,M)=\{\mathbf{x}:\|\mathbf{x}\|\leq M\}$). Letting $M\to\infty$ and using the Monotone and Bounded Convergence theorems yields the desired result. For more details the reader is referred to [13] where such arguments have been used in the case of a basic semi-algebraic set.

Finally it is worth emphasizing that this methodology also works for any measure $\mu$ that satisfies (3.15) (and whose moments are known or can be computed); an important spacial case is the exponential measure on the positive orthant $\mathbb{R}^{n}_{+}$.

Remark 3.5**.**

As mentioned above, in [13] the first author has already used Stokes’ formula to accelerate the convergence of a hierarchy of semidefinite relaxations to approximate the Gaussian measure $\mu(\mathbf{\Omega})$ of a basic semi-algebraic set $\mathbf{\Omega}$, not necessarily compact. The important and non trivial novelty here is that (i) $\mathbf{\Omega}=\bigcup_{i=1}^{p}\mathbf{\Omega}_{i}$ is now a union of basic semi-algebraic sets, and (ii) even if this complicates matters significantly, we are still able to work with measures $\phi_{i}$, each supported on $\mathbf{\Omega}_{i}$ (a basic semi-algebraic set). It turns out that $\mu(\mathbf{\Omega})=\sum_{i=1}^{p}\phi_{i}^{*}$ where is each $\phi^{*}_{i}$ has a piecewise constant density w.r.t. $\mu$ (constant on each of the possible intersections $\mathbf{\Omega}_{i}\cap\mathbf{\Omega}_{i_{1},\ldots,i_{p}}$). By using a family of polynomials that all vanish on the boundary of each $\mathbf{\Omega}_{i}\cap\mathbf{\Omega}_{i_{1},\ldots,i_{p}}$, we can exploit Stokes’ Theorem on each piece and sum up to obtain a family of linear constraints on the moments of $\phi^{*}_{i}$.

4. Numerical experiments and discussion

For illustration purposes we have applied the methodology on a few (simple) examples. We report some numerical experiments carried out in Matlab and GloptiPoly3 [11], a software package for manipulating and solving generalized problems of moments. The SDP problems were solved with SeDuMi 1.1R3.

4.1. Lebesgue volume of a union of two ellipsoids

We first consider a simple example of two ellipsoids in $\mathbb{R}^{2}$ where the exact value $\mu(\mathbf{\Omega})$ can be computed exactly so that we can compare with our upper bounds. So we want to compute the Lebesgue measure of $\mathbf{\Omega}=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ with $\mathbf{\Omega}_{1}=\{(x_{1},x_{2})\in\mathbb{R}^{2}:\frac{x_{1}^{2}}{4}+x_{2}^{2}\leq 1\}$ and $\mathbf{\Omega}_{2}=\{(x_{1},x_{2})\in\mathbb{R}^{2}:\frac{x_{2}^{2}}{4}+x_{1}^{2}\leq 1\}$. In this example we take $\mathbf{B}:=[-2,2]^{2}$.

The results are displayed in the Figure 1 with: in orange the approximation of the Lebesgue volume $\mu(\mathbf{\Omega})$ without using Stokes’ formulas, in red the approximation when using Stokes’ formulas and in blue the exact value of $\mu(\mathbf{\Omega})$.

We next consider a union of two ellipsoids in dimension $n=3$. Let $\mathbf{\Omega}_{1}=\{\mathbf{x}\in\mathbb{R}^{3}:x_{1}^{2}+4x_{2}^{2}+4x_{3}^{2}\leq 1\}$, $\mathbf{\Omega}_{2}=\{\mathbf{x}\in\mathbb{R}^{3}:x_{2}^{2}+4x_{1}^{2}+4x_{3}^{2}\leq 1\}$, $\mathbf{\Omega}=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ and $\mathbf{B}=[-1,1]^{3}$. Results are displayed in Figure 2. In both examples one can check that the convergence is much faster when using Stokes’ formula.

4.2. Lebesgue measure a union of three ellipsoids

We next consider a union of three ellipsoid in dimension $n=2$, with:

[TABLE]

[TABLE]

[TABLE]

and $\mathbf{B}=[-1,1]^{2}$. In Figure 3 we also compare our results with those obtained when using Bonferroni inequalities. In red the upper bounds obtained by solving $\textbf{Q}^{Stokes}$, in orange the lower bounds obtained by solving $\textbf{Q}^{Stokes}$ for the complement, and in blue the upper bounds obtained by using Bonferroni inequalities. (For a fair comparison, for each relaxation in Bonferroni case we also use appropriate Stokes’ constraints.)

4.3. Examples for the Gaussian measure

In this section we consider the Gaussian measure $d\mu=\exp(-\frac{\left\|\textbf{x}\right\|^{2}}{\sigma^{2}})d\textbf{x}$ with variance $\sigma^{2}=0.8$. For each example we have computed two upper-bounds and two lower-bounds for $\mu(\mathbf{\Omega})$. The first (resp. second) upper-bound $\overline{\rho}_{d}$ (resp. $\overline{\rho}^{Stokes}_{d}$) is obtained by solving the semidefinite relaxation $\textbf{Q}_{d}$ (resp. $\textbf{Q}^{Stokes}_{d}$). Similary, the lower-bounds $\underline{\rho}_{d}$ (resp. $\underline{\rho}^{Stokes}_{d}$) are obtained from upper bounds for the complement $\mathbb{R}^{n}\setminus\mathbf{\Omega}$. The respective relative error-gap are denoted by $\epsilon_{d}=\frac{\overline{\rho}_{d}-\underline{\rho}_{d}}{\overline{\rho}_{d}}$ and $\epsilon^{Stokes}_{d}=\frac{\overline{\rho}^{Stokes}_{d}-\underline{\rho}^{Stokes}_{d}}{\overline{\rho}^{Stokes}_{d}}$.

Example 1**.**

In this example $\mathbf{\Omega}$ is the union of two ellipsoids. Let $\mathbf{\Omega}:=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ with $\mathbf{\Omega}_{1}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{u})^{T}\textbf{A}_{1}(\textbf{x}-\textbf{u})\leq 1\}$ and $\mathbf{\Omega}_{2}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{v})^{T}\textbf{A}_{2}(\textbf{x}-\textbf{v})\leq 1\}$ for the values

[TABLE]

$\textbf{v}=(1,0)$,

[TABLE]

In this case $\rho:=\mu(\mathbf{\Omega}$ can be computed exactly and so we have displayed the values of the relative errors denoted by $\epsilon_{d}=\frac{\overline{\rho}_{d}-\underline{\rho}}{\overline{\rho}}$ and $\epsilon^{Stokes}_{d}=\frac{\overline{\rho}^{Stokes}_{d}-\underline{\rho}}{\overline{\rho}}$ respectively, depending on whether or not we have used Stokes’ formula. As one can see in Table 1 for a reasonable value $d=10$ the relative error (when using Stokes’ formula) is quite good. The respective behaviors are displayed in Figure 4.

Example 2**.**

Consider $\mathbf{\Omega}=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ with $\mathbf{\Omega}_{1}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{u})^{T}\textbf{A}_{1}(\textbf{x}-\textbf{u})\leq 1\}$ and $\mathbf{\Omega}_{2}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{v})^{T}\textbf{A}_{2}(\textbf{x}-\textbf{v})\leq 1\}$ for the values

$\textbf{u}=(0,0),$ $\textbf{v}=(-2,0)$,

[TABLE]

In this case $\mathbf{\Omega}$ is not a compact set as it is unbounded. The results for $d=9$ displayed in Table 2 show that a good value is already obtained when using Stokes’ formula. The respective behaviors of $\epsilon_{d}$ and $\epsilon^{Stokes}_{d}$ displayed in Figure 5 also show that using Stokes’ formula yields a significant improvement.

Example 3**.**

Consider $\mathbf{\Omega}=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ with $\mathbf{\Omega}_{1}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{u})^{T}\textbf{A}_{1}(\textbf{x}-\textbf{u})\leq 1\}$ and $\mathbf{\Omega}_{2}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{v})^{T}\textbf{A}_{2}(\textbf{x}-\textbf{v})\leq 1\}$ for the values

$\textbf{u}=(0,0),$ $\textbf{v}=(-2,0)$,

[TABLE]

Again $\mathbf{\Omega}$ is not compact. The results in Table 3 and the respective behaviors of $\epsilon_{d}$ and $\epsilon^{Stokes}_{d}$ displayed in Figure 6 confirm that using Stokes’ formula yields a significant improvement.

Example 4**.**

We next consider an example in dimension $n=3$. Let $\mathbf{\Omega}_{1}=\{\textbf{x}\in\mathbb{R}^{3}:(\textbf{x}-\textbf{u})^{T}\textbf{A}_{1}(\textbf{x}-\textbf{u})\leq 1\}$ and $\mathbf{\Omega}_{2}=\{\textbf{x}\in\mathbb{R}^{2}:(\textbf{x}-\textbf{v})^{T}\textbf{A}_{2}(\textbf{x}-\textbf{v})\leq 1\}$ for the values

$\textbf{u}=(0,0,0),$ $\textbf{v}=(-2,0,-1)$,

[TABLE]

Results for $d=6$ are displayed in Table 4 and the relative errors $\epsilon_{d}$ and $\epsilon^{Stokes}_{d}$ are displayed in Figure 7. The quality of results is comparable to that in Examples 2 and 3 for $d=6$.

Example 5**.**

Still in dimension $n=3$, let $\mathbf{\Omega}=\mathbf{\Omega}_{1}\cup\mathbf{\Omega}_{2}$ with $\mathbf{\Omega}_{1}=\{\textbf{x}\in\mathbb{R}^{3}:\textbf{x}^{T}e\leq 1\}$ and $\mathbf{\Omega}_{2}=\{\textbf{x}\in\mathbb{R}^{3}:\textbf{x}^{T}\textbf{A}\textbf{x}\ \leq 1\}$, where $e=(1,1,1)$ and

[TABLE]

The relative errors $\epsilon_{d}$ and $\epsilon^{Stokes}_{d}$ are displayed in Figure 8.

One can see that in all examples quite good approximations are obtained with relatively few moments (up to order $2d\leq 18$ for $n=2$ and $2d\leq 14$ for $n=3$) provided that we use the hierarchy (3.13) with the additional moments constraints induced by Stokes’ formula. The convergence of the hierarchy (3.5) (without those Stokes constraints) is indeed much slower.

For all the examples that we have treated, the (crucial) moment and localizing matrices involved in (3.5) and in (3.13) have been expressed in the canonical basis $(\mathbf{x}^{\alpha})_{\alpha\in\mathbb{N}^{n}}$ of monomials for simplicity and easyness of implementation of the SDP relaxations. But this choice is in fact the worst from a numerical point of view (numerical stability and robustness) which prevented us from solving (3.5) and (3.13) for $d\geq 7$ when $n=3$. It is very likely that the basis of orthonormal polynomials w.r.t. $\mu$ (Legendre for the Lebesgue measure $\mu$ on $[-1,1]$ and Hermite for the Gaussian measure $\mu$) is a much better (and recommended) choice. Such a more sophisticated implementation was beyond the scope of this paper.

Conclusion

In this paper we have provided a numerical scheme to approximate as closely as desired the measure $\mu(\mathbf{\Omega})$ of a finite union $\mathbf{\Omega}=\cup_{i=1}^{p}\mathbf{\Omega}_{i}$ of basic semi-algebraic sets (the case of a single basic semi-algebraic set was treated in [13])). Surprisingly, even though the case of a union of semi-algebraic sets complicates matters significantly we are still able to adapt the methodology developed in [13] and provide a monotone non-increasing (resp. non-decreasing) sequence of upper (resp. lower) bounds that converges to $\mu(\mathbf{\Omega})$ as the number of moments considered increases. In addition we are also able to use additional moment constraints induced by an appropriate application of Stokes’ Theorem which permits to improve significantly the convergence. In fact those additional moment constraints are crucial to obtain good bounds rapidly as they permit strongly attenuate a Gibbs’ phenomenon that otherwise appears.

Our current implementation could be significantly improved by using a basis for polynomials more appropriate than the usual canonical basis of monomials (the worst choice from a numerical stability point of view). For instance in doing so it should be possible to implement step $d=8,9$ of the hierarchy in dimension $n=3$, and step $d=7$ for $n=4$. As the convergence seems to be fast, each additional step of the hierarchy can yield a significant improvement.

The methodology was presented for the Lebesgue measure $\mu$ when $\mathbf{\Omega}$ is compact and the Gaussian measure for non-compact sets $\mathbf{\Omega}$, but in fact and remarkably, the same methodology works for any measure $\mu$ that satisfies Carleman’s condition and provided that all its moments are available (or can be computed easily).

Of course the methodology proposed in this paper is computationally expensive, especially when compared with Monte-Carlo type methods. But the latter provide only an estimate of $\mu(\mathbf{\Omega})$ and by no means an upper or lower bound on $\mu(\mathbf{\Omega})$ and therefore these two types of methods should be seen as complementary rather than competing. In its present form it is also limited to small dimension problems (typically $n\leq 3,4$) because since each upper (or lower) bound requires to solve a semidefinite program whose size increases fast in the hierarchy, one is limited by the current efficiency of state-of-the-art semidefinite solvers. However to the best of our knowledge this is the first method that provides a sequence of upper and lower bounds with strong asymptotic guarantees, at least at this level of generality.

Acknowledgement

Research funded by by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement ERC-ADG 666981 TAMING)”