Exact Optimization via Sums of Nonnegative Circuits and Sums of AM/GM Exponentials

TL;DR

This paper introduces hybrid numeric-symbolic algorithms for exact polynomial optimization using sums of nonnegative circuits (SONC) and sums of AM/GM exponentials (SAGE), with proven termination guarantees and experimental validation.

Contribution

It presents novel hybrid algorithms for exact polynomial optimization based on SONC and SAGE decompositions, including a decision algorithm with polynomial complexity.

Findings

Algorithms successfully compute exact decompositions.

The decision algorithm terminates within polynomial time.

Experimental results compare the effectiveness of the two methods.

Abstract

We provide two hybrid numeric-symbolic optimization algorithms, computing exact sums of nonnegative circuits (SONC) and sums of arithmetic-geometric-exponentials (SAGE) decompositions. Moreover, we provide a hybrid numeric-symbolic decision algorithm for polynomials lying in the interior of the SAGE cone. Each framework, inspired by previous contributions of Parrilo and Peyrl, is a rounding-projection procedure. For a polynomial lying in the interior of the SAGE cone, we prove that the decision algorithm terminates within a number of arithmetic operations, which is polynomial in the number of terms of the input, and linear in the distance to the boundary of the cone. We also provide experimental comparisons regarding the implementation of the two optimization algorithms.

Tables

Table 1. Table 1. Dependency of the average bit size and the average running time of optsage , with the number of variables, for fixed values of degree d 𝑑 d and number of terms t 𝑡 t ; For n = 8 𝑛 8 n=8 we observe a drastic drop both in running time and bit size.

$d=8$, $t=20$
$n$	bit size	time
2	27723	5.82
3	23572	4.99
4	22965	4.83
8	5678	1.12
10	1749	0.35

Table 2. Table 2. Bit size (upper part) and running time (lower part) of optsage in dependency of the degree d 𝑑 d and the number of terms t 𝑡 t for up to 4 variables; A “ × \times ” indicates, that we do not have instances with these parameters in our data set.

$t\setminus d$	6	8	10	18	20	26	28
6	912 0.24	1000 0.26	1002 0.26	1170 0.28	1014 0.26	955 0.28	900 0.26
9	2731 0.66	2673 0.65	2808 0.70	2890 0.68	2621 0.61	3166 0.82	2471 0.62
12	5599 1.30	6054 1.40	5449 1.27	5747 1.27	5478 1.21	6007 1.53	5027 1.18
20	9990 2.26	24078 5.08	20985 4.40	21364 4.43	19324 3.96	24096 5.23	17210 3.59
24	$\times$	36301 7.62	33414 6.99	37080 7.49	29266 5.87	37618 7.87	28090 5.43
30	$\times$	57744 11.90	56354 11.44	61564 12.57	48622 9.32	59975 12.76	55000 10.80
50	$\times$	$\times$	180971 36.11	174464 34.64	146218 27.80	196511 38.19	183598 38.36

Table 3. Table 3. Comparison of running time and bit size of the certificates between optsonc and optsage ; optsonc runs faster and has significantly smaller certificates than optsage .

$t$	bit size SONC	bit size SAGE	time SONC	time SAGE
6	432	1005	0.06	0.26
9	806	2696	0.19	0.66
12	1261	5568	0.37	1.29
20	2592	19203	0.64	4.00
24	3826	32543	0.97	6.66
30	5029	53160	1.34	10.58
50	10622	167971	3.95	32.78

Equations

$$\mathbf{\beta }=\sum _{j=1}^r{\lambda }_j\mathbf{\alpha }(j)\text{ with }{\lambda }_j>0\text{ and }\sum _{j=1}^r{\lambda }_j=1.$$

$$\displaystyle{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\Theta_{p}}\ =\ \prod_{j=1}^{r}\left(\frac{b_{\bm{\alpha}(j)}}{\lambda_{j}}\right)^{\lambda_{j}}.$$

$$\begin{array}{rl}\sum _{\mathbf{\alpha }\in \mathrm{MoSq}\left(p\right)}{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}\cdot \mathbf{\alpha }& =\mathbf{\beta }\\ \sum _{\mathbf{\alpha }\in \mathrm{MoSq}\left(p\right)}{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}& =1\\ {\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}& \ge 0\text{for all }\mathbf{\alpha }\in \mathrm{MoSq}\left(p\right)\end{array}$$

$$\begin{array}{rl}{p}_{\mathrm{SONC}}=\underset{\mathbf{X}}{\min }& \sum _{\mathbf{\beta }\in \mathrm{NoSq}\left(p\right)}{X}_{\mathbf{\beta },0}\\ \text{s.t.}& \sum _{\mathbf{\beta }\in \mathrm{NoSq}\left(p\right)}{X}_{\mathbf{\beta },\mathbf{\alpha }}\le b_{\mathbf{\alpha }},\mathbf{\alpha }\in \mathrm{MoSq}\left(p\right),\mathbf{\alpha }\ne 0,\\ & \prod _{\mathbf{\alpha }\in {\mathrm{Cov}}^{\mathbf{\beta }}}{\left(\frac{X_{\mathbf{\beta },\mathbf{\alpha }}}{{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}}\right)}^{{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}}=-b_{\mathbf{\beta }},\mathbf{\beta }\in \mathrm{NoSq}\left(p\right),\\ & X_{\mathbf{\beta },\mathbf{\alpha }}\ge 0,\mathbf{\alpha }\in \mathrm{MoSq}\left(p\right),\mathbf{\beta }\in \mathrm{NoSq}\left(p\right).\end{array}$$

$$\displaystyle{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathcal{C}_{\text{AGE}}}:=\left\{(\bm{c},\beta)\in\mathbb{R}_{+}^{t}\times\mathbb{R}:\text{There exists }\bm{\nu}\in\mathbb{R}_{+}^{t}\text{ with }D(\bm{\nu},e\bm{c})\leq\beta\,,\sum_{j=1}^{t}\bm{\alpha}(j)\nu_{j}=(\bm{1}\cdot\bm{\nu})\,\bm{\alpha}(0)\right\}\,.$$

$$\begin{array}{r}\sum _{j=1}^t{\mathbf{c}}^{(j)}=\mathbf{b},\sum _{i=1}^t\mathbf{\alpha }(i){\mathbf{\nu }}_i^{(j)}=0,-1\cdot {\mathbf{\nu }}_{\setminus j}^{(j)}={\nu }_j^{(j)},\\ {\mathbf{c}}_{\setminus j}^{(j)},{\mathbf{\nu }}_{\setminus j}^{(j)}\ge 0,D\left({\mathbf{\nu }}_{\setminus j}^{(j)},e{\mathbf{c}}_{\setminus j}^{(j)}\right)\le c_j^{(j)},j=1,\dots ,t.\end{array}$$

$$f_{\text{SAGE}}=\sup \{C\in \mathbb{R}:f-C\in {\mathcal{C}}_{\text{SAGE}}\}.$$

$$X_{\mathbf{\beta },0}\ge {\lambda }_0^{\mathbf{\beta }}\cdot {\left(-b_{\mathbf{\beta }}\cdot \prod _{\mathbf{\alpha }\in {\mathrm{Cov}}^{\mathbf{\beta }}}{\left(\frac{{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}}{X_{\mathbf{\beta },\mathbf{\alpha }}}\right)}^{{\lambda }_{\mathbf{\alpha }}^{\mathbf{\beta }}}\right)}^{\frac{1}{{\lambda }_0^{\mathbf{\beta }}}}.$$

$$\begin{array}{r}\sum _{j=1}^t{\mathbf{c}}^{(j)}=\mathbf{b},\sum _{i=1}^t\mathbf{\alpha }(i){\mathbf{\nu }}_i^{(j)}=0,-1\cdot {\mathbf{\nu }}_{\setminus j}^{(j)}={\nu }_j^{(j)},\\ {\mathbf{c}}_{\setminus j}^{(j)},{\mathbf{\nu }}_{\setminus j}^{(j)}>0,D\left({\mathbf{\nu }}_{\setminus j}^{(j)},e{\mathbf{c}}_{\setminus j}^{(j)}\right)<c_j^{(j)},j=1,\dots ,t.\end{array}$$

$$O(({\tau }_{\text{SAGE}}(f)+\tau +t\log t)\cdot t^7\log t)$$

$$D\left({\mathbf{\nu }}_{\setminus j}^{(j)},e\left({\mathbf{c}}_{\setminus j}^{(j)}+\delta 1\right)\right)+\frac{\epsilon }{2}\le D\left({\mathbf{\nu }}_{\setminus j}^{(j)},e{\mathbf{c}}_{\setminus j}^{(j)}\right)\le {\tilde{c}}_j^{(j)}-\frac{\epsilon }{2}=c_j^{(j)}.$$

$$t$$

$$\le d^{n-1}(1+d)\le 2d^n,$$

$$(\nu +\delta )\log (\nu +\delta )\le 0,$$

$$(\nu +\delta )\log (1+\delta /\nu )\le 2\delta ,$$

$$\delta \log \left(\frac{\nu }{ec}\right)\le \delta \max \left\{0,\nu \log \left(\frac{\nu }{ec}\right)\right\}.$$

$$D\left({\tilde{\mathbf{\nu }}}_{\setminus j}^{(j)}+\delta 1,e{\tilde{\mathbf{c}}}_{\setminus j}^{(j)}\right)=\sum _{i\ne j}({\tilde{\nu }}_i^{(j)}+\delta )\log \frac{{\tilde{\nu }}_i^{(j)}+\delta }{e{\tilde{c}}_i^{(j)}}.$$

$$({\tilde{\nu }}_i^{(j)}+\delta )\log \frac{{\tilde{\nu }}_i^{(j)}+\delta }{e{\tilde{c}}_i^{(j)}}$$

$$\le 2\delta \log (e{\epsilon }_{\mathbf{c}}).$$

$$({\tilde{\nu }}_i^{(j)}+\delta )$$

$$\le {\tilde{\nu }}_i^{(j)}\log \frac{{\tilde{\nu }}_i^{(j)}}{e{\tilde{c}}_i^{(j)}}+\delta \log \frac{{\tilde{\nu }}_i^{(j)}}{e{\tilde{c}}_i^{(j)}}+({\tilde{\nu }}_i^{(j)}+\delta )\log \left(1+\frac{\delta }{{\tilde{\nu }}_i^{(j)}}\right)$$

$$\le {\tilde{\nu }}_i^{(j)}\log \frac{{\tilde{\nu }}_i^{(j)}}{e{\tilde{c}}_i^{(j)}}+2\delta \log (e{\epsilon }_{\mathbf{c}})+2\delta ,$$

$$({\tilde{\nu }}_i^{(j)}+\delta )\log \frac{{\tilde{\nu }}_i^{(j)}+\delta }{e{\tilde{c}}_i^{(j)}}\le {\tilde{\nu }}_i^{(j)}\log \frac{{\tilde{\nu }}_i^{(j)}}{e{\tilde{c}}_i^{(j)}}+\delta \max \left\{0,{\tilde{\nu }}_i^{(j)}\log \frac{{\tilde{\nu }}_i^{(j)}}{e{\tilde{c}}_i^{(j)}}\right\}+2\delta ,$$

$$D\left({\tilde{\mathbf{\nu }}}_{\setminus j}^{(j)}+\delta 1,e{\tilde{\mathbf{c}}}_{\setminus j}^{(j)}\right)$$

$$\le {\tilde{c}}_j^{(j)}+2\delta {\tilde{c}}_j^{(j)}-(1+2\delta )\epsilon +2\delta t|\log (e{\epsilon }_{\mathbf{c}})|$$

$$\le {\tilde{c}}_j^{(j)}-(1+2\delta )\epsilon +2\delta (|b_j|+t|\log (e{\epsilon }_{\mathbf{c}})|),$$

$$f(\mathbf{x})=$$

$$-89x_1^2{x}_2^1{x}_3^2-19x_1^2{x}_2^2{x}_3^1+74x_1^2{x}_2^2{x}_3^2+268x_1^6{x}_3^2.$$

$$\displaystyle\begin{aligned} \bm{\nu}^{(2)}=&\ \left(\frac{1494563}{131072},-\frac{1494563}{65536},0,\frac{1494563}{131072},0,0,0,0,0,0,0,0\right)\\ f_{2}=&\ \frac{7050}{5161}-\frac{47048170074075847768}{2063044386600414657}x_{2}^{2}+\frac{51364821929347737990}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{1284035788022670567885}{13490375489389305583}x_{2}^{4}+\\ &\frac{83062810507300624}{22200197261503152705103}x_{1}^{1}x_{2}^{1}x_{3}^{2}+\frac{16345503627618}{597428155608650885}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{151893109261090080}{12048313292258664398599}x_{1}^{2}x_{3}^{4}+\\ &\frac{6887262866990276056}{3068799030848037634611567}x_{1}^{2}x_{2}^{1}x_{3}^{1}+\frac{584984347065145672}{102911615960676281785795}x_{1}^{2}x_{2}^{1}x_{3}^{2}+\frac{35180881401784}{5417152094682478713}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{910066980468240}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{7052856072901897195}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \bm{\nu}^{(5)}=&\ \biggl{(}\frac{2763713}{131072},\frac{1404885}{65536},\frac{40783}{131072},\frac{2618913}{131072},-\frac{1033291}{8192},\frac{46873}{65536},\frac{8086635}{131072},\frac{2825}{16384},\frac{16731}{65536},\frac{5689}{32768},\frac{4177}{16384},\frac{3431}{65536}\biggr{)}\\ f_{5}=&\ \frac{4895}{6532}+\frac{397583449413593513}{62516496563648929}x_{2}^{2}+\frac{4282065321967986445380570}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{5340068537046034265265}{107923003915114444664}x_{2}^{4}-\\ &\frac{797551428730610634080}{6323041088437240873}x_{1}^{1}x_{2}^{1}x_{3}^{2}+\frac{624998409184656940}{119485631121730177}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{2505598243831220279021280}{12048313292258664398599}x_{1}^{2}x_{3}^{4}+\\ &\frac{317731662232293486548800}{3068799030848037634611567}x_{1}^{2}x_{2}^{1}x_{3}^{1}+\frac{7986264919496237716952}{20582323192135256357159}x_{1}^{2}x_{2}^{1}x_{3}^{2}+\frac{4894573585100619175}{16251456284047436139}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{62000736112963044488}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{98753640749442855613080}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \bm{\nu}^{(6)}=&\ \biggl{(}\frac{24235}{65536},\frac{11469}{16384},\frac{731111}{65536},\frac{413881}{65536},\frac{75499}{65536},-\frac{2262417}{65536},\frac{224607}{65536},\frac{48937}{65536},\frac{41999}{32768},\frac{5895}{4096},\frac{459131}{65536},\frac{30411}{32768}\biggr{)}\\ f_{6}=&\ \frac{11}{6124}+\frac{175751812330814062}{6189133159801243971}x_{2}^{2}+\frac{20993411202657969962132460}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{230789345356723546855}{107923003915114444664}x_{2}^{4}+\\ &\frac{3497123471444427592320}{22200197261503152705103}x_{1}^{1}x_{2}^{1}x_{3}^{2}-\frac{20624391427747375598}{597428155608650885}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{19031996317336046221320}{12048313292258664398599}x_{1}^{2}x_{3}^{4}+\\ &\frac{188140322704886960215808}{3068799030848037634611567}x_{1}^{2}x_{2}^{1}x_{3}^{1}+\frac{27412383525960721370432}{102911615960676281785795}x_{1}^{2}x_{2}^{1}x_{3}^{2}+\frac{616440189731061280}{1805717364894159571}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{232955118417446202960}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{239384185705575737209071}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \bm{\nu}^{(8)}=&\ \biggl{(}\frac{5736955}{131072},\frac{1718597}{65536},\frac{2677}{65536},\frac{283517}{16384},\frac{298687}{65536},\frac{5897}{32768},\frac{1205677}{131072},-\frac{1398845}{8192},\frac{338517}{65536},\frac{488327}{32768},\frac{37365}{8192},\frac{1469765}{32768}\biggr{)}\\ f_{8}=&\ \frac{16553}{6389}+\frac{80194631513313539404}{6189133159801243971}x_{2}^{2}+\frac{936097366115717531112852}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{962838001240718755020}{13490375489389305583}x_{2}^{4}+\\ &\frac{168516279078159280107840}{22200197261503152705103}x_{1}^{1}x_{2}^{1}x_{3}^{2}+\frac{261920871973041930}{119485631121730177}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{622190516331822168107520}{12048313292258664398599}x_{1}^{2}x_{3}^{4}-\\ &\frac{174673447091258679729208576}{1022933010282679211537189}x_{1}^{2}x_{2}^{1}x_{3}^{1}+\frac{1345600783114821896667952}{102911615960676281785795}x_{1}^{2}x_{2}^{1}x_{3}^{2}+\frac{233241950370990311440}{5417152094682478713}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{1847380121666164588230}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{140920372909036228167589650}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \end{aligned}$$

$$\displaystyle\begin{aligned} \bm{\nu}^{(9)}=&\ \biggl{(}\frac{942879}{131072},\frac{66451}{8192},\frac{13687}{131072},\frac{291521}{32768},\frac{487245}{65536},\frac{29821}{65536},\frac{2150473}{65536},\frac{252191}{32768},-\frac{3460417}{32768},\frac{170395}{16384},\frac{548667}{65536},\frac{925733}{65536}\biggr{)}\\ f_{9}=&\ \frac{181}{1076}+\frac{3266564780864805056}{2063044386600414657}x_{2}^{2}+\frac{945568751423587451936415}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{195552867021926712375}{13490375489389305583}x_{2}^{4}+\\ &\frac{108598341793757685149280}{22200197261503152705103}x_{1}^{1}x_{2}^{1}x_{3}^{2}+\frac{261624143096996295}{119485631121730177}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{876811143650234639492646}{12048313292258664398599}x_{1}^{2}x_{3}^{4}+\\ &\frac{9330362356036127696945600}{3068799030848037634611567}x_{1}^{2}x_{2}^{1}x_{3}^{1}-\frac{10867831413671203786608456}{102911615960676281785795}x_{1}^{2}x_{2}^{1}x_{3}^{2}+\frac{192909634965734635520}{16251456284047436139}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{1339551331036614124520}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{17532007073614874704137330}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \bm{\nu}^{(10)}=&\ \biggl{(}\frac{732091}{262144},\frac{167933}{32768},\frac{14541}{131072},\frac{7773365}{262144},\frac{1327}{512},\frac{29683}{65536},\frac{272641}{131072},\frac{293905}{65536},\frac{214965}{65536},-\frac{4886959}{65536},\frac{144047}{32768},\frac{1284635}{65536}\biggr{)}\\ f_{10}=&\ \frac{157}{2744}+\frac{5424537902241772912}{6189133159801243971}x_{2}^{2}+\frac{880139855950527034906380}{176335433100237218699413}x_{2}^{2}x_{3}^{6}+\frac{571069367996398273530}{13490375489389305583}x_{2}^{4}+\\ &\frac{33169145578648933133280}{22200197261503152705103}x_{1}^{1}x_{2}^{1}x_{3}^{2}+\frac{228162076394260360}{119485631121730177}x_{1}^{1}x_{2}^{2}x_{3}^{3}+\frac{48697840759706172793616}{12048313292258664398599}x_{1}^{2}x_{3}^{4}+\\ &\frac{1587820308255205569373114}{1022933010282679211537189}x_{1}^{2}x_{2}^{1}x_{3}^{1}+\frac{59150503389680767180777}{20582323192135256357159}x_{1}^{2}x_{2}^{1}x_{3}^{2}-\frac{403951932103643744176}{5417152094682478713}x_{1}^{2}x_{2}^{2}x_{3}^{1}+\\ &\frac{616254425060782282800}{55380305977850685287}x_{1}^{2}x_{2}^{2}x_{3}^{2}+\frac{21315777567587124730738350}{672038441901303429467107}x_{1}^{6}x_{3}^{2}\\ \end{aligned}$$

Full text

For a polynomial lying in the interior of the SAGE cone, we prove that the decision algorithm terminates within a number of arithmetic operations, which is polynomial in the number of terms of the input, and linear in the distance to the boundary of the cone. We also provide experimental comparisons regarding the implementation of the two optimization algorithms.

1. Introduction

In this paper, we focus on certifying the output of polynomial optimization problems in a rigorous way. Finding the minimal value of a given polynomial in $n$ variables under polynomial constraints is known to be NP-hard in general [Lau09]. The related problem of deciding nonnegativity of a polynomial under polynomial constraints is co-NP hard; see e.g., [BCSS12].

This decision problem can be solved with the Cylindrical Algebraic Decomposition algorithm [Col75], which runs in time doubly exponential in $n$ and polynomial in the maximal total degree $d$ of the input functions. Further improved algorithms [GV88, BPR98, BGHP05], relying on critical point methods, allow to decide nonnegativity in singly exponential time in $n$. More generally, the complement of the problem lies in the existential theory of the reals, which can be solved in polynomial space and single exponential time [Ren88]. Safe validation of optimization problem results is mandatory for guaranted evaluation of mathematical functions [CHJL11], certified roundoff error bounds [MCD17] or computer assisted proofs [MAGW15, HAB*+*17].

In contrast to the algorithms mentioned above, several numerical frameworks have been developed in the last two decades. In the unconstrained case, one way to ensure nonnegativity of a given polynomial $f$ is to decompose $f$ as a * sum of squares (SOS)* of polynomials [Par00, Las01], which provides a certificate that $f$ is nonnegative over the reals. An SOS decomposition can be computed by solving a semidefinite program (SDP) of size $\binom{n+d}{n}$. In the constrained case, certificates can be provided by the prominent moment-SOS hierarchy, also called Lasserre’s hierarchy [Las01, Las10]. Each relaxation is solved with a semidefinite programming solver, implemented in finite-precision arithmetic, whose output is an approximate certificate. A drawback of these methods is that the size $\binom{n+d}{n}$ of the SDP matrices blows up when the degree $d$ and number of variables $n$ increases.

For larger values of $n$, a remedy consists of exploiting a potential sparsity/symmetry pattern arising in the input polynomials. A sparse version of Lasserre’s hierarchy has been developed in [WKKM06, Las06] when the objective function can be written as a sum of polynomials, each of them involving a small number of variables. The framework from [RTAL13] allows to take into account the symmetries of the polynomial optimization problem. One can also rely on the so-called bounded degree SOS hierarchy (BSOS) [LTY17]. In this hierarchy, positive polynomials are represented as a sum of two terms. The first term is an SOS polynomial of degree fixed in advance, while the second one belongs to the set of Krivine-Stengle representations, that is, is a finite combination of positive scalar weights and cross-products of the polynomials defining the set of constraints. This allows to handle larger instances than with the standard SOS hierarchy. The sparse variant of the BSOS hierarchy [WLT18] can handle even bigger problems, under the same sparsity pattern assumptions than the ones used for the sparse SOS hierarchy. If the support, i.e., number of monomial terms, of the polynomials is small in comparison to $\binom{n+d}{n}$, alternative relaxations based on geometric programming (GP) [DPZ67], and, more generally, relative entropy programming (REP), potentially allow to obtain lower bounds in a more efficient way than SDP relaxations. Both GP and REP are (equivalent to) convex optimization problems over the exponential cone. These alternative relaxations also provide the possibility to obtain answers when the SDP relaxations cannot be implemented because their size is too large for state-of-the art SDP solvers.

A first class of alternative certificates is given by sums of nonnegative circuit (SONC) polynomials. A circuit polynomial is a polynomial with support containing only monomial squares, at the exception of at most one term, whose exponent is a strict convex combination of the other exponents. In [IdW16], the authors derive a necessary and sufficient condition to prove that a given circuit polynomial is nonnegative. When the input polynomial has a more general support, a first attempt is given in [GM12, GM13] to compute lower bounds while relying on GP. This approach is generalized in [DIdW16] to compute SONC certificates when the set of constraints is defined as a finite conjunction of polynomial inequalities. In [DIdW17] the authors provide a bounded degree hierarchy, which can be computed via relative entropy programming. In the recent contribution [SdW18], the second and the third author develop an algorithm computing SONC certificates for sparse unconstrained polynomials with arbitrary support, together with a software library [SdW19], called POEM (Effective Methods in Polynomial Optimization). Although this framework yields a very efficient way to obtain a lower bound for a given polynomial, a drawback is that it currently remains unclear whether the number of circuits involved in a SONC relaxation is exponential in the number of terms of this polynomial. A second class of alternative certificates is given by sums of arithmetic-geometric-mean-exponentials (SAGE) polynomials. An AGE polynomial refers to a signomial, i.e., a weighted sum of exponentials composed with linear functionals of the variables, which is globally nonnegative with at most one negative coefficient. The framework from [CS16] derives a hierarchy of convex relaxations providing a sequence of increasing lower bounds for the optimal value of signomial programs (SP). For an input polynomial belonging to the SAGE cone, one can compute a SAGE decomposition by solving an REP, involving linear and relative entropy functions. Furthermore, it is shown in [MCW18, Theorem 20] that the cones of SAGE and SONC polynomials are related through their equivalence in terms of extreme rays. Namely, the extreme rays of the SAGE cone are supported on either a single coordinate or a set of coordinates inducing a simplicial circuit (a circuit with $\ell$ elements containing $\ell-1$ extreme points). Hence, both cones contain the same polynomials.

However, these alternative schemes share the same certification issues than the ones based on SDP relaxations. GP/REP/SP solvers rely on interior-point algorithms, implemented in finite-precision. Thus, they output only approximate certificates. In the unconstrained case, the rounding-projection procedure by Peryl and Parrilo [PP08] allows to compute a weighted rational SOS decompositon of a polynomial $f$ of degree $d=2k$, belonging to the interior of the SOS cone. In the “rounding” step, the algorithm computes an approximate Gram matrix of $f$, i.e., a matrix $\widetilde{\bm{G}}$ such that $f\simeq\bm{v}_{k}^{T}\widetilde{\bm{G}}\bm{v}_{k}$, where $\bm{v}_{k}$ is the vector of all monomials of degree at most $k$, then rounds $\widetilde{\bm{G}}$ in the space of rational matrices. In the “projection” step, the algorithm performs an orthogonal projection to obtain a matrix $\bm{G}$, such that $f=\bm{v}_{k}^{T}\bm{G}\bm{v}_{k}$. With sufficient precision digits, there always exists an SDP matrix fulfilling the above equality, yielding an (exact) Gram matrix associated to $f$. The last step to retrive an exact weighted SOS decomposition for $f$ consists of performing an exact $\bm{L}\bm{D}\bm{L}^{T}$ procedure [GL96, § 4.1]. Another framework [MD18a, MED21], provides a hybrid numeric-symbolic framework, which computes exact SOS decompositions under the same assumptions. This is based on a “perturbation-compensation” algorithm. In the “perturbation” step, one considers an arbitrary small perturbation of the input polynomial, and computes an approximate SOS decomposition with an SDP solver. In the “compensation” step, one relies on the perturbation terms to recover an exact SOS decomposition. By comparison with the rounding-projection procedure, this algorithm perturbates the input and provides an approximate $\bm{L}\bm{D}\bm{L}^{T}$ decomposition of the approximate Gram matrix, instead of a projection. It is shown in [MED21] that both procedures have a boolean running time, which is singly exponential in $n$ and polynomial in $d$. Practical experiments emphasize that the bit size of the SOS outputs obtained with the rounding-projection algorithm is often larger than the one obtained with the perturbation-compensation algorithm. The perturbation-compensation algorithm is inspired from prior work [MDS18], focusing on weighted SOS decompositions for nonnegative univariate polynomials, where the algorithm from [CHJL11] is formalized and analyzed. In the unconstrained case, the framework from [MED21] also extends the perturbation-compensation and rounding-projection algorithms to compute exact Polya and Hilbert-Artin’s representations, respectively for positive definite forms and nonnegative polynomials, yielding decompositions into SOS of rational functions, under the assumption that the numerator belongs to the interior of the SOS cone. In the constrained case, further algorithms are proposed to derive exact rational Putinar’s representations for positive polynomials over basic compact semialgebraic sets. All corresponding algorithms are integrated in the RealCertify [MD18b] Maple library.

The motivation of the present work is to improve the scalability of these existing certification frameworks, especially for large-size problems, which are currently out of reach when relying on SOS-based methods.

Contributions

In this paper, we provide a hybrid numeric-symbolic framework, in a similar spirit as [PP08], to certify exactly lower bounds obtained after computing SONC/SAGE decompositions with GP/REP relaxations. The resulting rounding-projection algorithms allow to handle unconstrained polynomial problems with such exact rational decompositions. Our first contributions, given in section 3, are two procedures, called optsonc and optsage, providing exact rational SONC and SAGE decompositions, respectively. These two algorithms allow to certify exactly lower bounds of unconstrained polynomials. Our framework is inspired from [PP08], first by rounding the output of a given GP/REP relaxation into rational numbers, next by performing an appropriate scaling of these numbers to obtain a solution satisfying exactly the (in)-equality constraints of the relaxation. We present another rounding-projection procedure called intsage in section 4 to handle the case of polynomials belonging to the interior of the SAGE cone. When the input is an $n$-variate polynomial of degree $d$ with $t$ monomial terms, integer coefficients of maximum bit size $\tau$, and with a distance to the SAGE cone of bit size upper bounded by $\tau_{\text{SAGE}}$, we prove that Algorithm intsage outputs SAGE decompositions within $\operatorname{\mathcal{O}}{((\tau_{\text{SAGE}}+\tau+t\log t)\cdot t^{7}\log t)}$ arithmetic operations. This is in contrast with the decision algorithm intsos from [MED21], which certifies nonnegativity of polynomials lying in the interior of the SOS cone in boolean time $\tau^{2}\cdot d^{d^{\operatorname{\mathcal{O}}{(n)}}}$. The two optimization algorithms optsonc and optsage are available within the POEM software library. In section 5, we provide experimental comparisons of these two algorithms.

Acknowledgements

Victor Magron benefited from the support of the FMJH Program PGMO (EPICS project) and EDF, Thales, Orange et Criteo. Timo de Wolff and Henning Seidler are supported by the DFG grant WO 2206/1-1. The authors would like to specially acknowledge the help of Riley Murray for providing insights about the barrier complexity of relative entropy programming.

2. Preliminaries

Let {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{Z}} be the set of integers and let {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{R}}, {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{R}_{\geq 0}} and {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{R}_{>0}} be the set of real, nonnegative real and positive real numbers, respectively. With {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{Q}} being the set of rational numbers, one defines similarly {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{Q}_{>0}}, {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{Q}_{\geq 0}}. The bit size of $i\in\mathbb{Z}$ is denoted by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\tau(i)}:=\lfloor\log_{2}(|i|)\rfloor+1 with {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\tau(0)}:=1. Given $i\in\mathbb{Z}$ and $j\in\mathbb{Z}\backslash\{0\}$ with gcd$(i,j)=1$, we define {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\tau(i/j)}:=\max\{\tau(i),\tau(j)\}. For two mappings $g,h:\mathbb{N}^{l}\to\mathbb{R}$, we use the notation $g(\bm{v})\in\operatorname{\mathcal{O}}{(h(\bm{v}))}$ to state the existence of $i\in\mathbb{N}$ such that $g(\bm{v})\leq ih(\bm{v})$, for all $\bm{v}\in\mathbb{N}^{l}$. Throughout the paper, we use bold letters for vectors (small) and matrices (capital), e.g., {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{x}}=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}. For a given vector $\bm{x}$, we denote the $j$-th coordinate of $\bm{x}$ by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}x_{j}}\in\mathbb{R}, and {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{x}_{\setminus j}}\in\mathbb{R}^{n-1} as the vector obtained by $\bm{x}$ after removing $x_{j}$. Furthermore, let $\mathbb{R}[\bm{x}]=\mathbb{R}[x_{1},\ldots,x_{n}]$ be the ring of real $n$-variate polynomials. We define {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{Q}[\bm{x}]} similarly. We denote the set of all $n$-variate polynomials of degree less than or equal to $2d$ by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\mathbb{R}[\bm{x}]_{n,2d}}.

For a polynomial $p\in\mathbb{R}[\bm{x}]$, we denote its support by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}A\left(p\right)}\subset\mathbb{N}^{n}; we just write $A$ if it is unambiguous and use the convention {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}t}=\#A. Thus, $p$ is of the form {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}p(\bm{x})}=\sum_{\bm{\alpha}\in A}b_{\bm{\alpha}}\bm{x}^{\bm{\alpha}} with {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}b_{\bm{\alpha}}}\in\mathbb{R}\setminus\{0\} and {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{x}^{\bm{\alpha}}}=x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{n}}. We say, a polynomial is sparse, if $t\ll\binom{n+2d}{2d}=\dim\left(\mathbb{R}[\bm{x}]_{n,2d}\right)$. The support of $p$ can be expressed as an $n\times t$ matrix, which we denote by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{A}}, such that the $j$-th column of $\bm{A}$ is $\bm{\alpha(j)}$. Hence, $p$ is uniquely described by the pair $(\bm{A},\bm{b})$, written p={\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{poly}(\bm{A},\bm{b})}. If $\bm{0}\in A$, then $p(\bm{0})$ is called the constant term.

Let us denote by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{New}(p)}:=\operatorname{conv}\left(\{\bm{\alpha}\in\mathbb{N}^{n}:b_{\bm{\alpha}}\neq 0\}\right) the * Newton polytope* of $p$ and {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{Vert}\left(p\right)}:=\left\{\bm{\alpha}\in A\left(p\right):\bm{\alpha}\text{ is vertex of }\operatorname{New}\left(p\right)\right\} be its vertices. We define the exponents of monomial squares in the support of $p$ as {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{MoSq}\left(p\right)}:=\left\{\bm{\alpha}\in A\left(p\right):\bm{\alpha}\in(2\mathbb{N})^{n},b_{\bm{\alpha}}>0\right\}. Its complement in the support is {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{NoSq}\left(p\right)}=A\left(p\right)\setminus\operatorname{MoSq}\left(p\right). We indicate the elements of the support which are in the interior of $\operatorname{New}(p)$ by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\operatorname{int}\left(p\right)}=A\left(p\right)\setminus\partial\operatorname{New}(p).

2.1. SONC Polynomials

We now introduce the fundamental properties of SONC polynomials, as far as they are required in this article. SONC are composed of circuit polynomials, which were first introduced in [IdW16]:

Definition 2.1**.**

A circuit polynomial $p=\operatorname{poly}(\bm{A},\bm{b})\in\mathbb{R}[\bm{x}]$ is of the form $p(\bm{x})=\sum_{j=1}^{r}b_{\bm{\alpha}(j)}\bm{x}^{\bm{\alpha}(j)}+b_{\bm{\beta}}\bm{x}^{\bm{\beta}}$, with 0\leq{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}r}<n, coefficients {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}b_{\bm{\alpha}(j)}}\in\mathbb{R}_{>0}, {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}b_{\bm{\beta}}}\in\mathbb{R}, exponents {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{\alpha}(j)}\in(2\mathbb{Z})^{n}, {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\bm{\beta}}\in\mathbb{Z}^{n}, such that the following condition holds: there exist unique, positive barycentric coordinates $\lambda_{j}$ relative to the $\bm{\alpha}(j)$ with $j=1,\ldots,r$ satisfying

[TABLE]

For every circuit polynomial $p$ we define the corresponding circuit number as

[TABLE]

∎

By Condition (2.1), $\bm{A}(p)$ forms a minimal affine dependent set, which is called a circuit, see e.g., [Oxl11]. More specifically, Condition (2.1) implies that $\operatorname{New}(p)$ is a simplex with even vertices $\bm{\alpha}(1),\ldots,\bm{\alpha}(r)$ and that the exponent $\bm{\beta}$ is in the strict interior of $\operatorname{New}(p)$ if $\dim(\operatorname{New}(p))\geq 1$. Therefore, we call $p_{\bm{\beta}}\bm{x}^{\bm{\beta}}$ the inner term of $p$.

Since the circuit number alone determines whether they are nonnegative, circuit polynomials are proper building blocks for nonnegativity certificates

Theorem 2.2** ([IdW16], Theorem 3.8).**

Let $p$ be a circuit polynomial as in definition 2.1. Then $p$ is nonnegative if and only if:

1. (1)

$p$* is a sum of monomial squares, or* 2. (2)

the coefficient $b_{\bm{\beta}}$ of the inner term of $p$ satisfies $|b_{\bm{\beta}}|\leq\Theta_{p}$.

The set of sums of nonnegative circuit polynomials (SONC) is a convex cone. For further details about SONC see [dW15, IdW16, DIdW17].

Let us consider $p=\operatorname{poly}(\bm{A},\bm{b})=\sum_{\bm{\alpha}\in A}b_{\bm{\alpha}}\bm{x}^{\bm{\alpha}}$. To algorithmically determine a lower bound of $p$ via SONC, we take an approach, similar to the one described in [SdW18, § 3]. As initial relaxation, each monomial non-square is equipped with a negative sign. This allows us, to restrict ourselves to the positive orthant; see e.g., [IdW16, Section 3.1] for further details. Next, we compute a covering $\operatorname{Cov}$, which is a sequence of sets $\left(\operatorname{Cov}^{\bm{\beta}}\right)_{\bm{\beta}\in\operatorname{NoSq}\left(p\right)}\subseteq A$ such that $\operatorname{NoSq}\left(p\right)\subseteq\bigcup_{\bm{\beta}\in\operatorname{NoSq}\left(p\right)}\operatorname{Cov}^{\bm{\beta}}$ and each $\operatorname{Cov}^{\bm{\beta}}$ is the support of a nonnegative circuit polynomial $p_{\bm{\beta}}$ with interior point $\bm{\beta}$. To obtain a covering, we write each non-square as a minimal convex combination of monomial squares, by solving the following LP for each $\bm{\beta}\in\operatorname{NoSq}\left(p\right)$, as explained in [SdW18, Algorithm 3.3].

[TABLE]

If $\left\{\bm{\alpha}:\lambda_{\bm{\alpha}}>0\right\}$ is not minimal, then we can reduce it by applying the following lemma. The lemma is folklore, for a constructive proof see e.g. [SdW18, Lemma 3.1].

Lemma 2.3**.**

For every non-extremal point $\bm{v}\in A\left(p\right)\setminus\operatorname{Vert}\left(p\right)$, we can efficiently compute affinely independent $\bm{v}_{0},...,\bm{v}_{m}\in\operatorname{Vert}\left(p\right)$ with $m\leq n$ such that and $\bm{v}\in\operatorname{conv}\left(\{\bm{v}_{0},...,\bm{v}_{m}\}\right)$.

So for each $\bm{\beta}$, we obtain a vector $\bm{\lambda}^{\bm{\beta}}$ of barycentric coordinates, relative to the simplex $\operatorname{Cov}^{\bm{\beta}}=\left\{\bm{\alpha}\in\operatorname{MoSq}\left(p\right):\lambda_{\bm{\alpha}}^{\bm{\beta}}>0\right\}$. We denote these computations by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\left(\bm{\lambda},\operatorname{Cov}\right):=\texttt{cover}(p)}.

Then, we solve the following geometric program (GP):

[TABLE]

For an overview about GPs, see [BKVH07, BV04]. If $p_{\operatorname{SONC}}$ is attained at $\bm{X^{\star}}$, then one has $p_{\bm{\beta}}=\sum_{\bm{\alpha}\in\operatorname{Cov}^{\bm{\beta}}}X^{\star}_{\bm{\beta},\bm{\alpha}}\cdot\bm{x}^{\bm{\alpha}}+b_{\beta}\bm{x}^{\bm{\beta}}\geq 0$ by theorem 2.2, and $p+p_{\operatorname{SONC}}-b_{\bm{0}}=\sum_{\bm{\beta}\in\operatorname{NoSq}\left(p\right)}p_{\bm{\beta}}\geq 0$. Hence, $b_{\bm{0}}-p_{\operatorname{SONC}}$ is a lower bound of $p$ on $\mathbb{R}^{n}$.

These computations correspond to [SdW18, § 3.3.2], where each cover contains just a single non-square. The approach dispayed here simplifies the computations, while keeping the quality of the results.

2.2. SAGE Polynomials

Let {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}e}:=\exp\,(1). The relative entropy function is defined for $\bm{\nu},\bm{c}\in\mathbb{R}_{+}^{t}$ by {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}D(\bm{\nu},\bm{c})}:=\sum_{j=1}^{t}\nu_{j}\log\frac{\nu_{j}}{c_{j}}. A signomial $p$ is a weighted sum of exponentials composed with linear functionals of a variable $\bm{x}\in\mathbb{R}^{n}$: given $t\in\mathbb{N}$, $c_{1},\dots,c_{t}\in\mathbb{Q}$ and $\bm{\alpha}(1),\dots,\bm{\alpha}(t)\in\mathbb{N}^{n}$, we write $p(\bm{x})=\sum_{j=1}^{t}c_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})$. Note that for general signomials, one considers $c_{1},\dots,c_{t}\in\mathbb{R}$ and $\bm{\alpha}(1),\dots,\bm{\alpha}(t)\in\mathbb{R}^{n}$. However, for certification purposes, we restrict the coefficients to the set of rationals and the exponents to tuples of nonnegative integers. A globally nonnegative signomial with at most one negative coefficient is called an AM/GM exponential or arithmetic-geometric-mean-exponential (AGE). Certifying the nonnegativity of an AGE is done by verifying an arithmetic-geometric-mean inequality. This is recalled in the following result, stated in [CS16, Lemma 2.2].

Lemma 2.4**.**

Let $p(\bm{x})=\sum_{j=1}^{t}c_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})+\beta\exp\,(\bm{\alpha}(0)\cdot\bm{x})$, with $c_{1},\dots,c_{t}\in\mathbb{Q}_{>0}$, $\beta\in\mathbb{Q}$ and $\bm{\alpha}(0),\bm{\alpha}(1),\dots,\bm{\alpha}(t)\in\mathbb{N}^{n}$. Then $p(\bm{x})\geq 0$ for all $\bm{x}\in\mathbb{R}^{n}$ if and only if there exists $\bm{\nu}\in\mathbb{R}_{+}^{t}$ such that $D(\bm{\nu},e\bm{c})\leq\beta$ and $\sum_{j=1}^{t}\bm{\alpha}(j)\nu_{j}=(\bm{1}\cdot\bm{\nu})\,\bm{\alpha}(0)$.

Given $\bm{\alpha}(0),\bm{\alpha}(1),\dots,\bm{\alpha}(t)\in\mathbb{N}^{n}$, the set of AGE signomials is a convex cone, denoted by $\mathcal{C}_{\text{AGE}}$, and defined as follows:

[TABLE]

The set of sums of AGE (SAGE) polynomials is also a convex cone, denoted by $\mathcal{C}_{\text{SAGE}}$. By [CS16, Proposition 2.4], one has the following characterisation.

Theorem 2.5**.**

A signomial $f=\sum_{i=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})$ lies in $\mathcal{C}_{\text{SAGE}}$ if and only if there is $\bm{c}^{(1)},\dots,\bm{c}^{(t)}$, $\bm{\nu}^{(1)},\dots,\bm{\nu}^{(t)}\in\mathbb{R}^{t}$ satisfying the following conditions:

[TABLE]

One way to obtain lower bounds of a signomial $f$ is to solve the following REP:

[TABLE]

The constraints of equation SAGE correspond to equation SAGE-feas, after replacing $\bm{b}$ by the vector of coefficients of $f-\lambda$.

3. Exact Optimization via SONC/SAGE

In this section, we present two algorithms for converting a numerical solution for SONC and SAGE into a lower bound in exact arithmetic. For a polynomial $p=\sum_{\bm{\alpha}\in A}b_{\bm{\alpha}}\bm{x}^{\bm{\alpha}}$, we assume $\bm{0}\in A$, so there exists a constant term $p(\bm{0})\neq 0$. Furthermore, we require that every non-square monomial lies in the interior of $\operatorname{New}{(p)}$ or on a face of $\operatorname{New}{(p)}$ including the origin. We denote numerical solutions {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\widetilde{\bm{X}}},{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\widetilde{\bm{c}}},{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\widetilde{\bm{\nu}}} with a tilde, intermediate rational solutions {\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\widehat{\bm{X}}},{\color[rgb]{0.2,0.2,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.75}\widehat{\bm{c}}} with a hat, and our final rational solution with regular letters.

3.1. Symbolic Post-Processing for SONC

We focus on certifying exactly lower bounds of a given polynomial via SONC decompositions. We rely on the numerical procedure from section 2.1, which starts to compute a covering of this polynomial. Under these assumptions, we design an algorithm, called optsonc, to convert numerical lower bounds, corresponding to SONC decompositions obtained via GP, into exact lower bounds.

Algorithm 3.1**.**

optsonc

1:$p=\sum_{\bm{\alpha}\in A}b_{\bm{\alpha}}\bm{x}^{\bm{\alpha}}\in\mathbb{Q}[\bm{x}]$, rounding precision $\widehat{\delta}\in\mathbb{Q}_{>0}$, precision parameter $\widetilde{\delta}\in\mathbb{Q}_{>0}$ for the GP solver.

2:Matrix $\bm{X}$ of rational numbers, coefficients of the decomposition, certified lower bound $C\in\mathbb{Q}$ of $p$ on $\mathbb{R}^{n}$.

3:$(\bm{\lambda},\operatorname{Cov})\leftarrow\texttt{cover}(p)$

4:$\widetilde{\bm{X}}\leftarrow\texttt{{GP}}(p,\widetilde{\delta},\bm{\lambda},\operatorname{Cov})$ $\triangleright$ Solve equation SONC with accuracy $\widetilde{\delta}$

5:$\widehat{\bm{X}}\leftarrow\texttt{round}\left(\widetilde{\bm{X}},\widehat{\delta}\right)$ $\triangleright$ rounding step

6:for $\bm{\alpha}\in\operatorname{MoSq}\left(p\right)$ and $\bm{\beta}\in\operatorname{NoSq}\left(p\right)$ do

7:     $X_{\bm{\beta},\bm{\alpha}}\leftarrow b_{\bm{\alpha}}\cdot\widehat{X}_{\bm{\beta},\bm{\alpha}}/\sum_{\bm{\beta^{\prime}}\in\operatorname{NoSq}\left(p\right)}\widehat{X}_{\bm{\beta^{\prime}},\bm{\alpha}}$ $\triangleright$ projection step

8:end for

9:for $\bm{\beta}\in\operatorname{NoSq}\left(p\right)$ do

10:     $\operatorname{coeff}=\lambda^{\bm{\beta}}_{\bm{0}}\cdot\left(-b_{\bm{\beta}}\cdot\prod_{\bm{\alpha}\in\operatorname{Cov}^{\bm{\beta}}}\left(\frac{\lambda_{\bm{\alpha}}^{\bm{\beta}}}{X_{\bm{\beta},\bm{\alpha}}}\right)^{\lambda_{\bm{\alpha}}^{\bm{\beta}}}\right)^{\frac{1}{\lambda_{\bm{0}}^{\bm{\beta}}}}$

11:     $X_{\bm{\beta},\bm{0}}\leftarrow\texttt{round-up}\left(\operatorname{coeff},\widehat{\delta}\right)$

$\triangleright$ adjust constant term

12:end for

13:$C\leftarrow b_{\bm{0}}-\sum_{\bm{\beta}\in\operatorname{NoSq}\left(p\right)}X_{\bm{\beta},\bm{0}}$

14:return $\bm{X}$, $C$

In 4, the function GP calls a GP solver to compute a $\widetilde{\delta}$-approximation $\widetilde{\bm{X}}$ of equation SONC. This approximation is then rounded in 5 to a rational point $\widehat{\bm{X}}$ with a prescribed maximal relative error of $\widehat{\delta}$. The projection step from 7 scales the entries of $\widehat{\bm{X}}$, yielding $\sum_{\bm{\beta}}X_{\bm{\beta},\bm{\alpha}}=b_{\bm{\alpha}}$, for all $\bm{\beta}\in\operatorname{MoSq}\left(p\right)$, to satisfy the first set of equality constraints of equation SONC. In 11, we round the coefficient up, with relative error $\widehat{\delta}$, so that we have

[TABLE]

As in section 2.1, each $p_{\bm{\beta}}:=\sum_{\bm{\alpha}\in\operatorname{Cov}^{\bm{\beta}}}X_{\bm{\beta},\bm{\alpha}}\cdot\bm{x}^{\bm{\alpha}}+b_{\beta}\bm{x}^{\bm{\beta}}$ is a nonnegative circuit polynomial. Hence, $C$ is a lower bound for $p$. Our assumption that every circuit polynomial contains a constant term, is necessary to ensure that $\lambda_{\bm{0}}^{\bm{\beta}}\neq 0$ for all $\bm{\beta}\in\operatorname{NoSq}\left(p\right)$ in our computations above.

3.2. Symbolic Post-Processing for SAGE

Similarly to algorithm 3.1, our algorithm optsage takes a given polynomial as input, obtains a numerical lower bound related to a SAGE decomposition computed via REP, and applies a post-processing to find a certified lower bound .

Algorithm 3.2**.**

optsage

1:$g=\sum_{i=1}^{t}b_{i}\bm{x}^{\bm{\alpha(i)}}\in\mathbb{Q}[\bm{x}]$, rounding precision $\widehat{\delta}\in\mathbb{Q}_{>0}$, precision parameter $\widetilde{\delta}\in\mathbb{Q}_{>0}$ for the REP solver.

2:Matrices $\bm{c},\bm{\nu}$ of rational numbers, coefficients of the decomposition, certified lower bound $C\in\mathbb{Q}$ of $g$ on $\mathbb{R}^{n}$.

3:$f\leftarrow g(\exp\bm{x}-\exp(-\bm{x}))$

4:Build the $(n+1)\times t$ matrix $\bm{Q}$ with columns $(\bm{\alpha}(1),1),\dots,(\bm{\alpha}(t),1)$

5:$\widetilde{\bm{c}},\widetilde{\bm{\nu}}\leftarrow\texttt{REP}(f,\widetilde{\delta})$ $\triangleright$ Solve equation SAGE with accuracy $\widetilde{\delta}$

6:$\widehat{\bm{c}}\leftarrow\texttt{round}\left(\widetilde{\bm{c}},\widehat{\delta}\right)\,,\quad\widehat{\bm{\nu}}\leftarrow\texttt{round}\left(\widetilde{\bm{\nu}},\widehat{\delta}\right)$ $\triangleright$ rounding step

7:for $j\in\{1,\ldots,t\}$ do

8:     LP $\leftarrow\left\{\bm{Q}\cdot\bm{\nu}^{(j)}=\bm{0},\bm{\nu}^{(j)}_{\setminus j}\geq\bm{0},\|\bm{\nu}^{(j)}-\widetilde{\bm{\nu}}^{(j)}\|_{\infty}\leq\widehat{\delta},\nu_{1}^{(j)}\geq\widehat{\delta}\right\}$

9:     $\bm{\nu^{(j)}}\leftarrow$ some element from LP $\triangleright$ projection step

10:     $\bm{c}_{\setminus j}^{(j)}\leftarrow\widehat{\bm{c}}_{\setminus j}^{(j)}\quad,c^{(j)}_{j}\leftarrow b_{j}-\bm{1}\cdot\bm{c}_{\setminus j}^{(j)}$

11:end for

12:for $j\in\{1,\ldots,t\}$ do

13:     $\operatorname{power}\leftarrow 1-\log\nu_{1}^{(j)}-\frac{1}{\nu_{1}^{(j)}}\left(c_{j}^{(j)}-\sum_{i>1,i\neq j}\nu_{i}^{(j)}\log\frac{\nu_{i}^{(j)}}{ec_{i}^{(j)}}\right)$

14:     $c_{1}^{(j)}\leftarrow\texttt{round-up}\left(\exp\left(\operatorname{power}\right),\widehat{\delta}\right)$ $\triangleright$ adjust constant term

15:end for

16:$C\leftarrow b_{1}-\sum_{j=1}^{t}c_{1}^{(j)}$

17:return $\bm{c}$, $\bm{\nu}$, $C$

Given a polynomial $g(\bm{y})=\sum_{j=1}^{t}b_{j}\bm{y}^{\bm{\alpha}(j)}$, one could apply the change of variables $y_{i}:=\exp x_{i}$ when $\bm{y}\in\mathbb{R}^{n}_{>0}$. Since this transformation is only valid on the nonnegative orthant, one workaround used in optsage is to define the signomial $f(\bm{x})=g(\exp\bm{x}-\exp(-\bm{x}))$ from 3, in a such a way that a lower bound of $f$ yields a lower bound of $g$. The REP function in 5 calls an REP solver to compute a $\widetilde{\delta}$-approximation $(\widetilde{\bm{\nu}},\widetilde{\bm{c}})$ of equation SAGE. This approximation is then rounded to a rational point $(\widehat{\bm{\nu}},\widehat{\bm{c}})$ with a prescribed maximal relative error of $\widehat{\delta}$. The projection steps in 9 and 10 ensure that $(\bm{\nu},\bm{c})$ satisfies exactly the linear equality constraints of  equation SAGE, i.e., $\bm{Q}\bm{\nu}^{(j)}=\bm{0}$ and $\sum_{j=1}^{t}\bm{c}^{(j)}=b$. The first projection step boils down to exactly solve an LP with the constraint that $\nu_{1}^{(j)}>0$, for all $j=1,\dots,t$, to ensure that further computation in 14 are well-defined. Note that this projection could be done while relying on the pseudo-inverse of $\bm{Q}$, but one obtains better practical results via this procedure. To ensure that the relative entropy inequality constraints of equation SAGE are satisfied, the last step of optsage aims at finding $c_{j}^{(1)}$ such that $c_{j}^{(j)}\geq D\left(\bm{\nu}_{\setminus j}^{(j)},e\bm{c}_{\setminus j}^{(j)}\right)=\sum_{i>1,i\neq j}\nu_{i}^{(j)}\log\frac{\nu_{i}^{(j)}}{ec_{i}^{(j)}}+\nu_{1}^{(j)}\log\frac{\nu_{1}^{(j)}}{ec_{1}^{(j)}}$. Thus, one relies on the round-up procedure in 14 to compute $c_{1}^{(j)}\geq\exp\left(1-\log\nu_{1}^{(j)}-\frac{1}{\nu_{1}^{(j)}}\left(c_{j}^{(j)}-\sum_{i>1,i\neq j}\nu_{i}^{(j)}\log\frac{\nu_{i}^{(j)}}{ec_{i}^{(j)}}\right)\right)$. Eventually, one has $\sum_{j=1}^{t}c_{i}^{(j)}=b_{i}$, for all $i>1$ and $\sum_{j=1}^{t}c_{1}^{(j)}=b_{1}-C$, which certifies that $f-C\geq 0$ on $\mathbb{R}^{n}$. We refer to section A.2 for an example of exact SAGE decomposition obtained with optsage.

4. Deciding Nonnegativity via SAGE

We denote by $\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ the interior of the cone $\mathcal{C}_{\text{SAGE}}$ of SAGE signomials. A signomial $f=\sum_{j=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})$ lies in $\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ if and only there is $\bm{c}^{(1)},\dots,\bm{c}^{(t)}$, $\bm{\nu}^{(1)},\dots,\bm{\nu}^{(t)}\in\mathbb{R}^{t}$ such that

[TABLE]

In the case where $f\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$, we denote by $\tau_{\text{SAGE}}(f)$ an upper bound of the bit size of the distance of $f$ to the boundary of $\mathcal{C}_{\text{SAGE}}$.

Without the assumptions from section 3, we state and analyze a decision algorithm to certify nonnegativity of signomials belonging to the interior $\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ of the SAGE cone. The resulting hybrid numeric-symbolic algorithm, called intsage, computes exact rational SAGE decompositions of such signomials.

Lemma 4.1**.**

Let $f=\sum_{j=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ of degree $d$ with $\tau=\tau(f)$. Then, there exists $N\in\mathbb{N}$ such that for $\varepsilon:=2^{-N}$, $f-\varepsilon\sum_{j=1}^{t}\exp\,(\bm{\alpha}(j)\cdot\bm{x})\in\mathcal{C}_{\text{SAGE}}$, with $N=\tau_{\text{SAGE}}(f)$.

Proof.

Since $f\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$, there are $\bm{c}^{(1)},\dots,\bm{c}^{(t)}$, $\bm{\nu}^{(1)},\dots,\bm{\nu}^{(t)}$ $\in\mathbb{R}^{t}$ such that $\sum_{j=1}^{t}\bm{c}^{(j)}=\bm{b}$, $\sum_{i=1}^{t}\bm{\alpha}(i)\bm{\nu}_{i}^{(j)}=\bm{0}$, $-\bm{1}\cdot\bm{\nu}_{\setminus j}^{(j)}=\nu_{j}^{(j)}$, $\bm{c}_{\setminus j}^{(j)},\bm{\nu}_{\setminus j}^{(j)}>\bm{0}$ and $D\left(\bm{\nu}_{\setminus j}^{(j)},e\bm{c}_{\setminus j}^{(j)}\right)<c_{j}^{(j)}$, for all $j=1,\dots,t$. Therefore, there exists $N\in\mathbb{N}$ such that for $\varepsilon:=2^{-N}$, one has $D\left(\bm{\nu}_{\setminus j}^{(j)},e\bm{c}_{\setminus j}^{(j)}\right)+\varepsilon<c_{j}^{(j)}$, for all $j=1,\dots,t$. For all $i,j=1,\dots,t$, let us define $\mathring{\bm{b}}$ by $\mathring{b}_{i}:=b_{i}-\varepsilon$, as well as $\mathring{\bm{c}}$ by $\mathring{c}_{i}^{(j)}:=c_{i}^{(j)}$ for $i\neq j$ and $\mathring{c}_{j}^{(j)}:=c_{j}^{(j)}-\varepsilon$. Note that $\mathring{\bm{b}}$ is the coefficient vector of $f-\varepsilon\sum_{j=1}^{t}\exp\,(\bm{\alpha}(j)\cdot\bm{x})$. Then $\mathring{\bm{c}}^{(1)},\dots,\mathring{\bm{c}}^{(t)}$, $\bm{\nu}^{(1)},\dots,\bm{\nu}^{(t)}$ satisfy equation INTSAGE-feas after replacing $\bm{b}$ by $\mathring{\bm{b}}$, yielding the first claim.

The second claim follows from the definition of the interior of the SAGE cone. ∎

4.1. Algorithm intsage

We present our algorithm intsage computing exact rational SAGE decompositions for signomials in $\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$.

Algorithm 4.2**.**

intsage

1:$f=\sum_{j=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$, rounding precision $\widehat{\delta}\in\mathbb{Q}_{>0}$, precision parameter $\widetilde{\delta}\in\mathbb{Q}_{>0}$ for the REP solver.

2:Matrices $\bm{c},\bm{\nu}$ of rational numbers.

3:Build the $(n+1)\times t$ matrix $\bm{Q}$ with columns $(\bm{\alpha}(1),1),\dots,(\bm{\alpha}(t),1)$

4:$\bm{Q}^{+}\leftarrow\texttt{pseudoinv}(\bm{Q})$

5:ok $\leftarrow$ false

6:while not ok do

7:     $(\widetilde{\bm{c}},\widetilde{\bm{\nu}})\leftarrow\texttt{REP}(f,\widetilde{\delta})$

8:     $\widehat{\bm{c}}\leftarrow\texttt{round}\left(\widetilde{\bm{c}},\widehat{\delta}\right)\,,\quad\widehat{\bm{\nu}}\leftarrow\texttt{round}\left(\widetilde{\bm{\nu}},\widehat{\delta}\right)$ $\triangleright$ rounding step

9:     for $j\in\{1,\dots,t\}$ do $\triangleright$ projection step

10:         $\bm{\nu}^{(j)}\leftarrow(\bm{I}-\bm{Q}^{+}\bm{Q})\,\widehat{\bm{\nu}}^{(j)}$

11:         $\bm{c}_{\setminus j}^{(j)}\leftarrow\widehat{\bm{c}}_{\setminus j}^{(j)}\,,\quad c^{(j)}_{j}\leftarrow b_{j}-\bm{1}\cdot\bm{c}_{\setminus j}^{(j)}$

12:     end for

13:     if for all $j\in\{1,\dots,t\}$, $\bm{\nu}^{(j)}_{\setminus j},\bm{c}^{(j)}_{\setminus j}\geq\bm{0}$, $c_{j}^{(j)}\geq D\left(\bm{\nu}^{(j)}_{\setminus j},e\bm{c}_{\setminus j}^{(j)}\right)$,  then ok $\leftarrow$ true $\triangleright$ verification step

14:     else$\ \widetilde{\delta}\leftarrow\widetilde{\delta}/2$, $\widehat{\delta}\leftarrow\widehat{\delta}/2$

15:     end if

16:end while

17:return $\bm{c}$, $\bm{\nu}$

The routine pseudoinv in 4 computes the pseudo-inverse of $\bm{Q}$, i.e., a matrix $\bm{Q}^{+}$ such that $\bm{Q}\bm{Q}^{+}\bm{Q}=\bm{Q}$. Next, we enter in the loop starting from 6. The REP function calls an REP solver to compute a $\widetilde{\delta}$-approximation $(\widetilde{\bm{\nu}},\widetilde{\bm{c}})$ of equation INTSAGE-feas. The projection steps ensure that $(\bm{\nu},\bm{c})$ satisfies exactly the linear equality constraints of  equation SAGE-feas, i.e., $\bm{Q}\bm{\nu}^{(j)}=\bm{Q}(\bm{I}-\bm{Q}^{+}\bm{Q})\bm{\nu}^{(j)}=\bm{Q}-\bm{Q}\bm{Q}^{+}\bm{Q}=\bm{0}$ and $\sum_{j=1}^{t}\bm{c}^{(j)}=b$. If the inequality constraints are not verified in 13, the rounding-projection procedure is performed again with more accuracy.

4.2. Arithmetic Complexity

Before analyzing the arithmetic complexity of intsage, we first establish lower bounds for the nonnegative components of the solutions related to SAGE decompositions of polynomials in $\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$.

Lemma 4.3**.**

Let $f=\sum_{j=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ of degree $d$ with $\tau=\tau(f)$. Let $\varepsilon$ be as in lemma 4.1.

1. (1)

There exists a solution of $(\bm{\nu},\bm{c})$ of equation INTSAGE-feas and $\delta\in\mathbb{Q}_{>0}$ such that $\delta\leq 1$, $(\bm{\nu},{\bm{c}})$ satisfies, $D\left({\bm{\nu}}_{\setminus j}^{(j)},e\left({\bm{c}}_{\setminus j}^{(j)}+\delta\bm{1}\right)\right)+\frac{\varepsilon}{2}\leq{c}_{j}^{(j)}$, $\sum_{j=1}^{t}\bm{c}^{(j)}=\bm{b}$, for all $i,j=1,\dots,t$. 2. (2)

There exists a solution $(\bm{\nu},\bm{c})$ of equation INTSAGE-feas and $\delta\in\mathbb{Q}_{>0}$ such that $(\bm{\nu},{\bm{c}})$ satisfies $D\left({\bm{\nu}}_{\setminus j}^{(j)}+\delta\bm{1},e{\bm{c}}_{\setminus j}^{(j)}\right)+\frac{\varepsilon}{2}\leq{c}_{j}^{(j)}$, for all $j=1,\dots,t$. 3. (3)

There exists a solution $(\bm{\nu},\bm{c})$ of equation INTSAGE-feas and $\delta\in\mathbb{Q}_{>0}$ such that $(\bm{\nu},{\bm{c}})$ satisfies $D\left((1+\delta){\bm{\nu}}_{\setminus j}^{(j)},e{\bm{c}}_{\setminus j}^{(j)}\right)+\frac{\varepsilon}{2}\leq{c}_{j}^{(j)}$, for all $j=1,\dots,t$.

*In each case, $\tau(\delta)\in\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+\log t)}=\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+n\log d)}$. *

A proof for this lemma is provided in section A.1.

Theorem 4.4**.**

*Let $f=\sum_{j=1}^{t}b_{j}\exp\,(\bm{\alpha}(j)\cdot\bm{x})\in\operatorname{int}\left(\mathcal{C}_{\text{SAGE}}\right)$ of degree $d$ and $\tau=\tau(f)$.

There exist $\widehat{\delta}$ and $\widetilde{\delta}$ of bit size less than $\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+t\log t)}$, such that $\texttt{intsos}(f,\widehat{\delta},\widetilde{\delta})$ terminates and outputs a rational SAGE decomposition of $f$ within*

[TABLE]

*arithmetic operations. *

Proof.

We first show that the loop of Algorithm intsage terminates with $\widehat{\delta}$ and $\widetilde{\delta}$ of bit size bounded by $\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+t\log t)}$. Let $\varepsilon$ be as in lemma 4.1. When running the procedure REP, one solves equation SAGE-feas at precision $\widetilde{\delta}$, thus one finds an approximate solution $(\widetilde{\bm{\nu}},\widetilde{\bm{c}})$ such that $\|\sum_{j=1}^{t}\widetilde{\bm{c}}^{(j)}-\bm{b}\|_{\infty}\leq\widetilde{\delta}$, $D\left({\widetilde{\bm{\nu}}}_{\setminus j}^{(j)},e{\widetilde{\bm{c}}}_{\setminus j}^{(j)}\right)+\varepsilon\leq\widetilde{c}_{j}^{(j)}+\widetilde{\delta}$, and $\|\bm{Q}\widetilde{\bm{\nu}}^{(j)}\|_{\infty}\leq\widetilde{\delta}$, for all $j=1,\dots,t$. After the rounding and projection steps, one obtains $\bm{\nu}^{(j)}=(\bm{I}-\bm{Q}^{+}\bm{Q})\widehat{\bm{\nu}}^{(j)}$ and $\|\widehat{\bm{\nu}}^{(j)}-\widetilde{\bm{\nu}}^{(j)}\|_{\infty}\leq\widehat{\delta}$. Since $\bm{\alpha}_{j}(i)\leq d$, for all $i,j=1,\dots,t$, the bit size of the entries of the matrix $\bm{Q}$ is upper bounded by $\tau(d)$. Thus, the pseudo-inverse $\bm{Q}^{+}$ has rational entries of bit size bounded by $\operatorname{\mathcal{O}}{(t\log t+t\log d)}=\operatorname{\mathcal{O}}{(t\log t)}$, since the bit size is the same as for the determinant length, see [BPR98, Corollary 8.13]. This implies that the bit size of the difference between the entries of $\widetilde{\bm{\nu}}$ and $\bm{\nu}$ is upper bounded by \operatorname{\mathcal{O}}{\Big{(}t\log t+\tau(\widetilde{\delta})+\tau(\widehat{\delta})\Big{)}}. Similarly, the bit size of the difference between the entries of $\widetilde{\bm{c}}$ and $\bm{c}$ is upper bounded by \operatorname{\mathcal{O}}{\Big{(}\tau+\tau(\widehat{\delta})\Big{)}}. By lemma 4.3, one can perform any absolute or relative perturbation of $\widetilde{\bm{\nu}}$ and $\widetilde{\bm{c}}$, and still ensure that the resulting $(\bm{\nu},\bm{c})$ satisfies $D\left(\bm{\nu}^{(j)}_{\setminus j},e\bm{c}_{\setminus j}^{(j)}\right)\leq c_{j}^{(j)}$, if the perturbation is small enough with bit size at most $\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+\log t)}$. This implies that one must choose $\widetilde{\delta}$ and $\widehat{\delta}$ small enough, with an upper bound of $\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+t\log t)}$ on their bit sizes. The same reasoning applies to ensure that $\bm{\nu}_{\setminus j}^{(j)},\bm{c}_{\setminus j}^{(j)}\geq\bm{0}$.

Now, we give an upper bound on the number of arithmetic operations. For convex optimization problems having barrier complexity equal to $N$, the standard interior-point methods compute a $\widetilde{\delta}$-accurate solution in $\operatorname{\mathcal{O}}{(\tau(\widetilde{\delta})\sqrt{N}\log(N))}$ iterations; see e.g. [Ren01, Section 2.4]. For $(x,y,z)\in\mathbb{R}^{3}$, the standard barrier complexity of a single relative entropy constraint “$x\log(x/y)\leq z,x,y\geq 0$” is equal to 4. In addition, the barrier complexity of a set of constraints is upper bounded by the sum of the complexities of the individual constraints. Therefore, the relative entropy formulation given in equation SAGE-feas has a barrier complexity of $N\leq 4t^{2}$. At each iteration of the interior-point method, one needs to solve an LP involving $2t^{2}$ variables, which can be done within $\operatorname{\mathcal{O}}{(t^{6})}$ arithmetic operations. This yields the upper bound of $\operatorname{\mathcal{O}}{((\tau_{\text{SAGE}}(f)+\tau+t\log t)\cdot t^{7}\log t)}$ on the total number of arithmetic operations required while calling intsage. All other arithmetic operations performed by the algorithm have a negligible cost with respect to the REP procedure. ∎

5. Experimental Comparisons

We discuss the actual bit sizes and physical running time of optsonc and optsage procedures, given by algorithm 3.1 and algorithm 3.2. We describe the setup of our experiment and on which instances the algorithms were tested. Afterwards, we discuss our findings from running the algorithms on a large set of examples.

5.1. Experimental Setup

We give an overview about the experimental setup.

Software The entire experiment was steered by our Python 3.7 based software POEM 0.2.0.0(a) (Effective Methods in Polynomial Optimization), [SdW19], which we develop since July 2017. POEM is open source, under GNU public license, and available at:

https://www3.math.tu-berlin.de/combi/RAAGConOpt/poem.html

For our experiment, POEM calls a range of further software and solvers for computing the certificates. For the numerical solutions of SONC and SAGE, we use CVXPY 1.0.12 [DB16], to create the convex optimization problems, together with the solver ECOS 2.0.7 [DCB13]. The symbolic computations were done in SymPy 1.3 [JvMG12].

Investigated Data We carried out our experiments on 2020 randomly generated polynomials. The possible numbers of variables are $n=2,3,4,8,10$; the degree takes values $d=6,8,10,18,20,26,28$ and the number of terms can be $t=6,9,12,20,24,30,50$. For each combinations we create instances, where the number of negative terms is one of a few fixed ratios of $t$. In particular, the size of equation SAGE grows quadratically in $t$. We created the database using POEM, and it is available in full at the homepage cited above. Our instances are a subset of those from [SdW18]. In that paper, we also describe their creation in more detail. The overall running time for all our instances was 6780.0 seconds.

Hardware and System We used an Intel Core i7-8550U CPU with 1.8 GHz, 4 cores, 8 threads and 16 GB of RAM under Ubuntu 18.04 for our computations.

Stopping Criteria For the accuracy of the solver and the precision of the rounding in Python we used a tolerance of $\varepsilon=2^{-23}$. The restriction $t\leq 50$ was chosen, since otherwise we already encounter problem in the numerical solution of equation SAGE. The bound $d<30$ was chosen, because for large degree we had a significant increase in the memory required to perform the rounding. Both thresholds were obtained experimentally.

5.2. Evaluation of the Experiment

In this section we present and evaluate the results of our experiment and highlight our most important findings, when investigating the computational data. We focus on the results given by the procedure optsage (algorithm 3.2) via SAGE decompositions and in the end give a comparison to optsonc.

Running time decreases with growing number of variables. The formulation of equation SAGE shows that the size of the problem only depends on the number of terms $t$, but in the SAGE decomposition, the number of summands is the number of monomial non-squares. Most significantly, for more variables, our generating algorithm simply results in a smaller number of these terms. Additionally, for $n\geq 8$ and $d\leq 10$, most exponents lie on faces of the Newton polytope. This leads to a simpler combinatorial structure, which we believe to result in lower bit sizes and thus in faster solving faster the exact LP from 8 of optsonc. Next, we have more equality constraints in this LP, which could also improve the running time. Lastly, the exponential upper bound is just the worst case, which does not seem to actually happen among our examples.

For some selected parameters, we exhibit that behavior in Table 1.

Dependency of bit size and running time of degree and terms

To illustrate how bit size and running time of optsage vary for different degrees and numbers of terms, we restrict ourselves to at most 4 variables. Our numbers from the previous point show, that in these cases bit size and time are similar for fixed $(d,t)$, hence we may aggregate those instances. The results are shown in Table 2. We can see that running time and bit size roughly have a linear dependency. On the one hand, their growth is quadratic in the number of terms, which matches with the growth of the problem size in equation SAGE. On the other hand, bit size and running time are basically unaffected by the degree. This shows that the bound, given in the worst case analysis, usually is not met.

Quality of the rounding-projection Our experiments verify that in the majority of cases the symbolical lower bound does not diverge far from the numerical bound. The detailed distribution is shown in figure 1. Most notably, in 30 instances, the exact lower bound is even better than the numerical bound. In 81.9% of the instances, the exact bound differs by at most 0.001 from the numerical value. Only in 256 instances the difference lies above 1. Thus, in the clear majority of examples, the lower bound in exact arithmetic does not differ much from the numerical bound. Also, among the instances with large difference, it can also be that the numerical solution actually lies far away from an exact solution. So it is unclear, whether a large difference is due to bad behavior of the numerical solution, or a large error in the rounding algorithm.

Rounding time versus solving time In nearly every case the rounding procedure takes longer than the numerical solving. Only in 8 instances, the rounding took less time. The ratio of the rounding time to the total time ranges from 21.6% to 96.8%, with an average of 88.6%. However, one can implement the rounding procedure much closer to the hardware level, instead of working in Python. Thus, we expect that these ratios can be significantly improved.

Comparison between SONC and SAGE In their qualitative behavior, optsonc and optsage are similar. However, optsonc runs faster and has smaller certificates than optsage, as shown in Table 3. But one should note that optsonc only computes some lower bound (not necessarily the optimal SONC-bound), whereas optsage computes the best bound, that can be obtained via this approach. Still it shows, that for very large instances, SONC is the method of choice, when other approaches fail due to the problem size.

Comparison with SOS For polynomials lying in the interior of the SOS cone from [MED21, Table 2], we performed preliminary experiments with optsage and optsage, which are currently unable to provide nonnegativity certificates. For benchmarks from our database with $n\geq 8$ and $d\geq 10$, RealCertify often fails to provide SOS certificates. We plan to provide detailed experimental comparisons with SOS methods in the future.

6. Conclusion and Outlook

We make two main contributions in this paper. First, we present an algorithm to decide whether a given multivariate polynomial over the rationals lies in the interior of the SAGE cone. If that is the case, then the algorithm also computes a certificate in exact arithmetic. Additionally, we analyze the arithmetic complexity of the algorithm, which is singly exponential in the degree and the number of terms, and doubly exponential in the number of variables. Second, we use our numerical methods to obtain lower bounds via SAGE and apply a single iteration of the rounding-projection method, to obtain an exactly certified lower bound. This method, we run on a large number of test cases. Based on these experiments, we draw the following conclusions.

1. (1)

In the majority of cases, the exact solution lies close to the numerical solution, with a difference of at most 0.001. For few instances, the exact lower bound is even better than the numerical one. 2. (2)

The running time and the bit size grow quadratically in the number of terms, which corresponds to the growth of the problem size. 3. (3)

For the investigated parameters, increasing the degree or the number of variables does not increase the running time or the bit size. This also corresponds to the fact, that the size of the REP is independent of the degree and the number of variables. 4. (4)

For very large instances, SONC should be the first choice, to obtain a certified bound, since it runs significantly faster than the other methods.

For future work, the most interesting development would be to have an REP-solver with arbitrary precision, so that we can actually implement intsage and compare it to similar approaches. Furthermore, in a significant amount of instances, we encountered computational problem, when calling optsonc or optsage. So we would like to increase the robustness of our implementation. Another issue, we have left out so far, is the presence of exponents of monomial non-squares, which lie on a face of the Newton polytope, that does not include the origin. These result in values $\lambda_{1,j}=0$ (SONC) or $\nu_{1}^{(j)}=0$ (SAGE), so our computation is undefined. However, these problems can be circumvented and we plan to do so in a future version of the software. Next, we plan to extend our framework to constrained problems and provide more detailed experimental comparisons with SOS-based approaches from [PP08, MED21], as well as with methods based on critical points and cylindrical algebraic decomposition. Given a polynomial in the interior of the SAGE cone, our decision algorithm intsage is linear with respect to the distance of this polynomial to the border of the cone. In order to improve this bound, one could refine the bit size analysis for this distance. A further theoretical aim would be to analyze the boolean running time of intsage, which requires to prove bit complexity estimates for relative entropy optimization problems.

Appendix A Appendix

A.1. Proof of lemma 4.3.

We start with the first claim. By lemma 4.1, there exist $\widetilde{\bm{\nu}},\widetilde{\bm{c}}$ and $\varepsilon\in\mathbb{Q}_{>0}$, with $\tau(\varepsilon)\in\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f))}$, such that $(\widetilde{\bm{\nu}},\widetilde{\bm{c}})$ satisfies $\sum_{j=1}^{t}\bm{c}^{(j)}=\bm{b}$, $\widetilde{\bm{c}}_{\setminus j}^{(j)},\widetilde{\bm{\nu}}_{\setminus j}^{(j)}>\bm{0}$, $D\left(\widetilde{\bm{\nu}}_{\setminus j}^{(j)},e\widetilde{\bm{c}}_{\setminus j}^{(j)}\right)+\varepsilon\leq\widetilde{c}_{j}^{(j)}$, for all $j=1,\dots,t$. Let us define $\delta:=\frac{\varepsilon}{2(t-1)}$ and $\bm{c}$ such that $c_{i}^{(j)}:=\widetilde{c}_{i}^{(j)}+\frac{\varepsilon}{2}$ for all $i\neq j=1,\dots,t$, and $c_{j}^{(j)}:=\widetilde{c}_{j}^{(j)}-t\delta$, for $j=1,\dots,t$. Thus, one has $\sum_{j=1}^{t}\bm{c}^{(j)}=\bm{b}$. For all $j=1,\dots,t$, one has $\bm{\nu}^{(j)}_{\setminus j}>\bm{0}$. Combining this together with the fact that the $\log$ function is increasing yields

[TABLE]

Using that for all $d\geq 2$,

[TABLE]

one has $\tau(t)\in\operatorname{\mathcal{O}}{(n\log_{2}(d))}$. Since $\tau(\delta)\leq\tau(\varepsilon)+\tau(t)$, we obtain the first claim.

To prove the second claim, we rely on the following three auxiliary inequalities.

For $\nu,\delta$ such that $0<\nu<\delta\leq\frac{1}{2}$, one has

[TABLE]

since $0<\nu+\delta<1$ and the function $x\mapsto x\log x$ is negative on $(0,1)$.

For $\nu,\delta$ such that $0<\delta\leq\nu$, one has

[TABLE]

since $\nu+\delta\leq 2\nu$ and $\log(1+\delta/\nu)\leq\delta/\nu$.

For each $\nu,\delta,c>0$ such that $\nu\geq 1$, one has

[TABLE]

Indeed, if $\nu\leq ec$, the left hand side is less than 0. Otherwise, $\nu\geq 1$ implies that $\log(\frac{\nu}{ec})\leq\nu\log(\frac{\nu}{ec})$.

Now, by the first claim, there exist $\varepsilon_{\bm{c}},\varepsilon\in\mathbb{Q}_{>0}$, with $\tau(\varepsilon_{\bm{c}}),\tau(\varepsilon)\in\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\log t)}$, and $\widetilde{\bm{\nu}},\widetilde{\bm{c}}$ satisfying $\widetilde{\bm{\nu}}_{\setminus j}^{(j)}>\bm{0}$, $\widetilde{\bm{c}}_{\setminus j}^{(j)}>\varepsilon_{\bm{c}}\bm{1}$, and $D\left(\widetilde{\bm{\nu}}_{\setminus j}^{(j)},e\widetilde{\bm{c}}_{\setminus j}^{(j)}\right)+\varepsilon\leq\widetilde{c}_{j}^{(j)}$, for all $j=1,\dots,t$. For each $\delta\in\mathbb{Q}_{>0}$ with $\delta\leq\frac{1}{2}$, and all $j=1,\dots,t$, one has:

[TABLE]

We give an upper bound of each summand of the right hande side, depending on the value of $\widetilde{\nu}_{i}^{j}$, for all $i,j=1,\dots,t$ and $i\neq j$. Note that $-\log(ec_{i}^{j})\leq\log(e\varepsilon_{\bm{c}})$, for all $i,j=1,\dots,t$ and $i\neq j$.

If $0<\widetilde{\nu}_{i}^{j}\leq\delta\leq\frac{1}{2}$, one has by equation A.1

[TABLE]

If $\delta\leq\widetilde{\nu}_{i}^{j}\leq 1$, one has

[TABLE]

where we use the fact that $\delta\log\widetilde{\nu}_{i}^{(j)}\leq 0$ and bound the last term of the right hand side via equation A.2.

If $\widetilde{\nu}_{i}^{j}\geq 1$, we write the first inequality as in the former case and obtain

[TABLE]

where we rely on equation A.3 to bound the second term and the fact that $\delta+\widetilde{\nu}_{i}^{(j)}\leq 2\widetilde{\nu}_{i}^{(j)}$ together with $\log\left(1+\frac{\delta}{\widetilde{\nu}_{i}^{(j)}}\right)\leq\frac{\delta}{\widetilde{\nu}_{i}^{(j)}}$ to bound the last term.

In the worst case, we obtain

[TABLE]

using $\widetilde{c}_{j}^{(j)}=b_{j}-\bm{1}\cdot\widetilde{c}_{\setminus j}^{(j)}\leq b_{j}$, thus $\widetilde{c}_{j}^{(j)}\leq\max\{0,b_{j}\}\leq|b_{j}|$.

To ensure that $D\left(\widetilde{\bm{\nu}}_{\setminus j}^{(j)}+\delta\bm{1},e\widetilde{\bm{c}}_{\setminus j}^{(j)}\right)+\frac{\varepsilon}{2}\leq\widetilde{c}_{j}^{(j)}$, it is sufficient to have $-(1+4\delta)\varepsilon+4\delta(|b_{j}|+t|\log(e\varepsilon_{\bm{c}})|)\leq 0$, which is guaranteed by selecting the largest positive rational $\delta$ such that $\delta\leq\frac{\varepsilon}{4(|b_{j}|+t|\log(e\varepsilon_{\bm{c}})|-\varepsilon)}$. Since $\tau(|b_{j}|)\leq\tau$, $\tau(t)\in\operatorname{\mathcal{O}}{(d\log_{2}n)}$, $|\log(e\varepsilon_{\bm{c}})|)\leq\tau(\varepsilon_{\bm{c}})$, and $\tau(\varepsilon_{\bm{c}}),\tau(\varepsilon)\in\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\log t)}$, one can select $\delta$ with bit size at most $\operatorname{\mathcal{O}}{(\tau_{\text{SAGE}}(f)+\tau+\log t)}$.

The proof of the third claim is very similar and we omit it for the sake of conciseness.

A.2. An Example of Exact SAGE Decomposition

Let

[TABLE]

Our optimization algorithm optsage returns $(\bm{\nu},\bm{c})$ corresponding to the following exact rational SAGE decomposition: $f(\bm{x})=\sum_{j=1}^{12}f_{j}(\bm{x})$, where $f_{j}$ is the polynomial with coefficient vector $\bm{c}^{(j)}$ for $j\in\{1,\dots,12\}$, $f_{j}=0$ for $j\in\{1,3,4,7,11,12\}$, and

[TABLE]

[TABLE]