Exact Optimization via Sums of Nonnegative Circuits and Sums of AM/GM Exponentials
aa r X i v : . [ c s . S C ] F e b EXACT OPTIMIZATION VIA SUMS OF NONNEGATIVE CIRCUITSAND SUMS OF AM/GM EXPONENTIALS
VICTOR MAGRON, HENNING SEIDLER, AND TIMO DE WOLFF
Abstract.
We provide two hybrid numeric-symbolic optimization algorithms, com-puting exact sums of nonnegative circuits (SONC) and sums of arithmetic-geometric-exponentials (SAGE) decompositions. Moreover, we provide a hybrid numeric-symbolicdecision algorithm for polynomials lying in the interior of the SAGE cone. Each frame-work, inspired by previous contributions of Parrilo and Peyrl, is a rounding-projectionprocedure.For a polynomial lying in the interior of the SAGE cone, we prove that the decisionalgorithm terminates within a number of arithmetic operations, which is polynomial inthe degree and number of terms of the input, and singly exponential in the number ofvariables. We also provide experimental comparisons regarding the implementation ofthe two optimization algorithms. Introduction
In this paper, we focus on certifying the output of polynomial optimization problemsin a rigorous way. Finding the minimal value of a given polynomial in n variables underpolynomial constraints is known to be NP-hard in general [Lau09]. The related problemof 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 AlgebraicDecomposition algorithm [Col75], which runs in time doubly exponential in n and poly-nomial in the maximal total degree d of the input functions. Further improved algorithms[GV88, BPR98, BGHP05], relying on critical point methods, allow to decide nonnegativ-ity in singly exponential time in n . More generally, the complement of the problem lies inthe existential theory of the reals, which can be solved in polynomial space and single ex-ponential time [Ren88]. Safe validation of optimization problem results is mandatory forguaranted evaluation of mathematical functions [CHJL11], certified roundoff error bounds[MCD17] or computer assisted proofs [MAGW15, HAB + f is to decompose f as a sum of squares (SOS) ofpolynomials [Par00, Las01], which provides a certificate that f is nonnegative over the Date : February 10, 2019.2010
Mathematics Subject Classification.
Primary: 14P10, 68W30, 90C25; Secondary: 14Q20, 68R01
ACM Subject Classification:
Mathematical software performance.
Key words and phrases. nonnegative circuit polynomial, arithmetic-geometric-mean exponential, con-vex optimization, geometric programming, relative entropy programming, exact certificate, rounding-projection procedure, hybrid numeric-symbolic algorithm, real algebraic geometry. reals. An SOS decomposition can be computed by solving a semidefinite program (SDP) ofsize (cid:0) n + dn (cid:1) . In the constrained case, certificates can be provided by the prominent moment-SOS hierarchy , also called Lasserre’s hierarchy [Las01, Las10]. Each relaxation is solvedwith a semidefinite programming solver, implemented in finite-precision arithmetic, whoseoutput is an approximate certificate. A drawback of these methods is that the size (cid:0) n + dn (cid:1) 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/symmetrypattern arising in the input polynomials. A sparse version of Lasserre’s hierarchy hasbeen developed in [WKKM06, Las06] when the objective function can be written as asum of polynomials, each of them involving a small number of variables. The frame-work from [RTAL13] allows to take into account the symmetries of the polynomial op-timization 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 twoterms. The first term is an SOS polynomial of degree fixed in advance, while the secondone belongs to the set of Krivine-Stengle representations, that is, is a finite combinationof positive scalar weights and cross-products of the polynomials defining the set of con-straints. This allows to handle larger instances than with the standard SOS hierarchy. Thesparse variant of the BSOS hierarchy [WLT18] can handle even bigger problems, underthe 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 comparisonto (cid:0) n + dn (cid:1) , alternative relaxations based on geometric programming (GP) [DPZ67], and,more generally, relative entropy programming (REP) , potentially allow to obtain lowerbounds in a more efficient way than SDP relaxations. Both GP and REP are (equivalentto) convex optimization problems over the exponential cone. These alternative relax-ations also provide the possibility to obtain answers when the SDP relaxations cannot beimplemented 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 mono-mial squares, at the exception of at most one term, whose exponent is a strict convexcombination of the other exponents. In [IdW16], the authors derive a necessary and suf-ficient condition to prove that a given circuit polynomial is nonnegative. When the inputpolynomial has a more general support, a first attempt is given in [GM12, GM13] tocompute lower bounds while relying on GP. This approach is generalized in [DIdW16] tocompute SONC certificates when the set of constraints is defined as a finite conjunctionof polynomial inequalities. In [DIdW17] the authors provide a bounded degree hierar-chy, which can be computed via relative entropy programming. In the recent contribu-tion [SdW18a], the second and the third author develop an algorithm computing SONCcertificates for sparse unconstrained polynomials with arbitrary support, together witha software library [SdW18b], called POEM (Effective Methods in Polynomial Optimiza-tion) . Although this framework yields a very efficient way to obtain a lower bound for agiven polynomial, a drawback is that it currently remains unclear whether the number of
XACT OPTIMIZATION VIA SONC AND SAGE 3 circuits involved in a SONC relaxation is exponential in the number of terms of this poly-nomial. 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, whichis globally nonnegative with at most one negative coefficient. The framework from [CS16]derives a hierarchy of convex relaxations providing a sequence of increasing lower boundsfor the optimal value of signomial programs (SP). For an input polynomial belonging tothe SAGE cone, one can compute a SAGE decomposition by solving an REP, involvinglinear 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 interms of extreme rays. Namely, the extreme rays of the SAGE cone are supported oneither a single coordinate or a set of coordinates inducing a simplicial circuit (a circuitwith ℓ elements containing ℓ − extreme points). Hence, both cones contain the samepolynomials.However, these alternative schemes share the same certification issues than the onesbased on SDP relaxations. GP/REP/SP solvers rely on interior-point algorithms, im-plemented in finite-precision. Thus, they output only approximate certificates. In theunconstrained case, the rounding-projection procedure by Peryl and Parrilo [PP08] allowsto compute a weighted rational SOS decompositon of a polynomial f of degree d = 2 k ,belonging to the interior of the SOS cone. In the “rounding” step, the algorithm computesan approximate Gram matrix of f , i.e., a matrix e G such that f ≃ v Tk e Gv k , where v k is thevector of all monomials of degree at most k , then rounds e G in the space of rational matri-ces. In the “projection” step, the algorithm performs an orthogonal projection to obtain amatrix G , such that f = v Tk Gv k . With sufficient precision digits, there always exists anSDP 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 performingan exact LDL T procedure [GL96, § 4.1]. Another framework [MD18a, MSED18], pro-vides a hybrid numeric-symbolic framework, which computes exact SOS decompositionsunder the same assumptions. This is based on a “perturbation-compensation” algorithm.In the “perturbation” step, one considers an arbitrary small perturbation of the inputpolynomial, and computes an approximate SOS decomposition with an SDP solver. Inthe “compensation” step, one relies on the perturbation terms to recover an exact SOSdecomposition. By comparison with the rounding-projection procedure, this algorithmperturbates the input and provides an approximate LDL T decomposition of the approx-imate Gram matrix, instead of a projection. It is shown in [MSED18] that both proce-dures 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 therounding-projection algorithm is often larger than the one obtained with the perturbation-compensation algorithm. The perturbation-compensation algorithm is inspired from priorwork [MDS18], focusing on weighted SOS decompositions for nonnegative univariate poly-nomials, where the algorithm from [CHJL11] is formalized and analyzed. In the uncon-strained case, the framework from [MSED18] also extends the perturbation-compensation XACT OPTIMIZATION VIA SONC AND SAGE 4 and rounding-projection algorithms to compute exact Polya and Hilbert-Artin’s repre-sentations, respectively for positive definite forms and nonnegative polynomials, yieldingdecompositions into SOS of rational functions, under the assumption that the numera-tor belongs to the interior of the SOS cone. In the constrained case, further algorithmsare proposed to derive exact rational Putinar’s representations for positive polynomialsover 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 cer-tification frameworks, especially for large-size problems, which are currently out of reachwhen relying on SOS-based methods.
Contributions . In this paper, we provide a hybrid numeric-symbolic framework, in asimilar spirit as [PP08], to certify exactly lower bounds obtained after computing SON-C/SAGE decompositions with GP/REP relaxations. The resulting rounding-projectionalgorithms allow to handle unconstrained polynomial problems with such exact ratio-nal decompositions. Our first contributions, given in Section 3, are two procedures, called optsonc and optsage , providing exact rational SONC and SAGE decompositions, respec-tively. These two algorithms allow to certify exactly lower bounds of unconstrained poly-nomials. Our framework is inspired from [PP08], first by rounding the output of a givenGP/REP relaxation into rational numbers, next by performing an appropriate scaling ofthese numbers to obtain a solution satisfying exactly the (in)-equality constraints of therelaxation. We present another rounding-projection procedure called intsage in Section 4to handle the case of polynomials belonging to the interior of the SAGE cone. When theinput is an n -variate polynomial of degree d with t monomial terms, and integer coefficientsof maximum bit size τ , we prove that Algorithm intsage outputs SAGE decompositionswithin O ( τ · (4 d + 6) n +3 t log t ) arithmetic operations. This is in contrast with the deci-sion algorithm intsos from [MSED18], which certifies nonnegativity of polynomials lyingin the interior of the SOS cone in boolean time O ( τ · (4 d + 2) n +6 ) . The two optimiza-tion 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 ProgramPGMO (EPICS project) and EDF, Thales, Orange et Criteo. Timo de Wolff and HenningSeidler are supported by the DFG grant WO 2206/1-1. The authors would like to speciallyacknowledge the help of Riley Murray for providing insights about the barrier complexityof relative entropy programming. 2.
Preliminaries
Let Z be the set of integers and let R , R ≥ and R > be the set of real, nonnegativereal and positive real numbers, respectively. With Q being the set of rational numbers,one defines similarly Q > , Q ≥ . The bit size of i ∈ Z is denoted by τ ( i ) := ⌊ log ( | i | ) ⌋ + 1 with τ (0) := 1 . Given i ∈ Z and j ∈ Z \{ } with gcd ( i, j ) = 1 , we define τ ( i/j ) :=max { τ ( i ) , τ ( j ) } . For two mappings g, h : N l → R , we use the notation g ( v ) ∈ O ( h ( v )) tostate the existence of i ∈ N such that g ( v ) ≤ ih ( v ) , for all v ∈ N l . Throughout the paper,we use bold letters for vectors (small) and matrices (capital), e.g., x = ( x , . . . , x n ) ∈ R n . XACT OPTIMIZATION VIA SONC AND SAGE 5
For a given vector x , we denote the j -th coordinate of x by x j ∈ R , and x \ j ∈ R n − asthe vector obtained by x after removing x j . Furthermore, let R [ x ] = R [ x , . . . , x n ] bethe ring of real n -variate polynomials. We define Q [ x ] similarly. We denote the set of all n -variate polynomials of degree less than or equal to d by R [ x ] n, d .For a polynomial p ∈ R [ x ] , we denote its support by A ( p ) ⊂ N n ; we just write A if it isunambiguous and use the convention t = A . Thus, p is of the form p ( x ) = P α ∈ A b α x α with b α ∈ R \ { } and x α = x α · · · x α n n . We say, a polynomial is sparse , if t ≪ (cid:0) n +2 d d (cid:1) =dim ( R [ x ] n, d ) . The support of p can be expressed as an n × t matrix, which we denoteby A , such that the j -th column of A is α ( j ) . Hence, p is uniquely described by the pair ( A , b ) , written p = poly( A , b ) . If ∈ A , then p ( ) is called the constant term .Let us denote by New( p ) := conv ( { α ∈ N n : b α = 0 } ) the Newton polytope of p and Vert ( p ) := { α ∈ A ( p ) : α is vertex of New ( p ) } be its vertices. We define the exponentsof monomial squares in the support of p as MoSq ( p ) := { α ∈ A ( p ) : α ∈ (2 N ) n , b α > } .Its complement in the support is NoSq ( p ) = A ( p ) \ MoSq ( p ) . We indicate the elementsof the support which are in the interior of New( p ) by int ( p ) = A ( p ) \ ∂ New( p ) .2.1. SONC Polynomials.
We now introduce the fundamental properties of SONC poly-nomials, as far as they are required in this article. SONC are composed of circuit polyno-mials , which were first introduced in [IdW16]:
Definition 2.1. A circuit polynomial p = poly( A , b ) ∈ R [ x ] is of the form p ( x ) = P rj =1 b α ( j ) x α ( j ) + b β x β , with ≤ r < n , coefficients b α ( j ) ∈ R > , b β ∈ R , exponents α ( j ) ∈ (2 Z ) n , β ∈ Z n , such that the following condition holds: there exist unique,positive barycentric coordinates λ j relative to the α ( j ) with j = 1 , . . . , r satisfying β = r X j =1 λ j α ( j ) with λ j > and r X j =1 λ j = 1 . (2.1)For every circuit polynomial p we define the corresponding circuit number as Θ p = r Y j =1 (cid:18) b α ( j ) λ j (cid:19) λ j . (2.2) By Condition (2.1), A ( p ) forms a minimal affine dependent set, which is called a circuit ,see e.g., [Oxl11]. More specifically, Condition (2.1) implies that New( p ) is a simplex witheven vertices α (1) , . . . , α ( r ) and that the exponent β is in the strict interior of New( p ) if dim(New( p )) ≥ . Therefore, we call p β x β the inner term of p .Since the circuit number alone determines whether they are nonnegative, circuit poly-nomials 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) p is a sum of monomial squares, or(2) the coefficient b β of the inner term of p satisfies | b β | ≤ Θ p . XACT OPTIMIZATION VIA SONC AND SAGE 6
The set of sums of nonnegative circuit polynomials (SONC) is a convex cone. Forfurther details about SONC see [dW15, IdW16, DIdW17].Let us consider p = poly( A , b ) = P α ∈ A b α x α . To algorithmically determine a lowerbound of p via SONC, we take an approach, similar to the one described in [SdW18a,§ 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, Section3.1] for further details. Next, we compute a covering Cov , which is a sequence of sets (cid:0)
Cov β (cid:1) β ∈ NoSq( p ) ⊆ A such that NoSq ( p ) ⊆ S β ∈ NoSq( p ) Cov β and each Cov β is the supportof a nonnegative circuit polynomial p β with interior point β . To obtain a covering, wewrite each non-square as a minimal convex combination of monomial squares, by solvingthe following LP for each β ∈ NoSq ( p ) , as explained in [SdW18a, Algorithm 3.3]. X α ∈ MoSq( p ) λ βα · α = β X α ∈ MoSq( p ) λ βα = 1 λ βα ≥ for all α ∈ MoSq ( p ) (LP)If { α : λ α > } is not minimal, then we can reduce it by applying the following lemma.The lemma is folklore, for a constructive proof see e.g. [SdW18a, Lemma 3.1]. Lemma 2.3.
For every non-extremal point v ∈ A ( p ) \ Vert ( p ) , we can efficientlycompute affinely independent v , ..., v m ∈ Vert ( p ) with m ≤ n such that and v ∈ conv ( { v , ..., v m } ) . So for each β , we obtain a vector λ β of barycentric coordinates, relative to the simplex Cov β = (cid:8) α ∈ MoSq ( p ) : λ βα > (cid:9) . We denote these computations by ( λ , Cov) := cover ( p ) .Then, we solve the following geometric program (GP): p SONC = min X X β ∈ NoSq( p ) X β , s.t. X β ∈ NoSq( p ) X β , α ≤ b α , α ∈ MoSq ( p ) , α = , Y α ∈ Cov β (cid:18) X β , α λ βα (cid:19) λ βα = − b β , β ∈ NoSq ( p ) ,X β , α ≥ , α ∈ MoSq ( p ) , β ∈ NoSq ( p ) . (SONC)For an overview about GPs, see [BKVH07, BV04]. If p SONC is attained at X ⋆ , thenone has p β = P α ∈ Cov β X ⋆ β , α · x α + b β x β ≥ by Theorem 2.2, and p + p SONC − b = P β ∈ NoSq( p ) p β ≥ . Hence, b − p SONC is a lower bound of p on R n .These computations correspond to [SdW18a, § 3.3.2], where each cover contains just asingle non-square. The approach dispayed here simplifies the computations, while keepingthe quality of the results. XACT OPTIMIZATION VIA SONC AND SAGE 7
SAGE Polynomials.
Let e := exp (1) . The relative entropy function is defined for ν , c ∈ R t + by D ( ν , c ) := P tj =1 ν j log ν j c j . A signomial p is a weighted sum of exponentialscomposed with linear functionals of a variable x ∈ R n : given t ∈ N , c , . . . , c t ∈ Q and α (1) , . . . , α ( t ) ∈ N n , we write p ( x ) = P tj =1 c j exp ( α ( j ) · x ) . Note that for general sig-nomials, one considers c , . . . , c t ∈ R and α (1) , . . . , α ( t ) ∈ R n . However, for certificationpurposes, we restrict the coefficients to the set of rationals and the exponents to tuples ofnonnegative integers. A globally nonnegative signomial with at most one negative coeffi-cient 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-meaninequality. This is recalled in the following result, stated in [CS16, Lemma 2.2]. Lemma 2.4.
Let p ( x ) = P tj =1 c j exp ( α ( j ) · x ) + β exp ( α (0) · x ) , with c , . . . , c t ∈ Q > , β ∈ Q and α (0) , α (1) , . . . , α ( t ) ∈ N n . Then p ( x ) ≥ for all x ∈ R n if and only if thereexists ν ∈ R t + such that D ( ν , e c ) ≤ β and P tj =1 α ( j ) ν j = ( · ν ) α (0) . Given α (0) , α (1) , . . . , α ( t ) ∈ N n , the set of AGE signomials is a convex cone, denotedby C AGE , and defined as follows: C AGE := ( c , β ) ∈ R t + × R : There exists ν ∈ R t + with D ( ν , e c ) ≤ β , t X j =1 α ( j ) ν j = ( · ν ) α (0) . The set of sums of AGE (SAGE) polynomials is also a convex cone, denoted by C SAGE .By [CS16, Proposition 2.4], one has the following characterisation.
Theorem 2.5.
A signomial f = P ti =1 b j exp ( α ( j ) · x ) lies in C SAGE if and only if thereis c (1) , . . . , c ( t ) , ν (1) , . . . , ν ( t ) ∈ R t satisfying the following conditions: t X j =1 c ( j ) = b , t X i =1 α ( i ) ν ( j ) i = , − · ν ( j ) \ j = ν ( j ) j , c ( j ) \ j , ν ( j ) \ j ≥ , D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) ≤ c ( j ) j , j = 1 , . . . , t . (SAGE-feas)One way to obtain lower bounds of a signomial f is to solve the following REP: f SAGE = sup { C ∈ R : f − C ∈ C SAGE } . (SAGE)The constraints of (SAGE) correspond to (SAGE-feas), after replacing b by the vector ofcoefficients of f − λ . 3. Exact Optimization via SONC/SAGE
In this section, we present two algorithms for converting a numerical solution for SONCand SAGE into a lower bound in exact arithmetic. For a polynomial p = P α ∈ A b α x α ,we assume ∈ A , so there exists a constant term p ( ) = 0 . Furthermore, we requirethat every non-square monomial lies in the interior of New ( p ) or on a face of New ( p ) including the origin. We denote numerical solutions f X , e c , e ν with a tilde, intermediaterational solutions c X , b c with a hat, and our final rational solution with regular letters. XACT OPTIMIZATION VIA SONC AND SAGE 8
Symbolic Post-Processing for SONC.
We focus on certifying exactly lowerbounds of a given polynomial via SONC decompositions. We rely on the numerical proce-dure from Section 2.1, which starts to compute a covering of this polynomial. Under theseassumptions, 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
Require: p = P α ∈ A b α x α ∈ Q [ x ] , rounding precision b δ ∈ Q > , precision parameter e δ ∈ Q > for the GP solver. Ensure:
Matrix X of rational numbers, coefficients of the decomposition, certified lowerbound C ∈ Q of p on R n . ( λ , Cov) ← cover ( p ) f X ← GP ( p, e δ, λ , Cov) ⊲ Solve (SONC) with accuracy e δ c X ← round (cid:16) f X , b δ (cid:17) ⊲ rounding step for α ∈ MoSq ( p ) and β ∈ NoSq ( p ) do X β , α ← b α · b X β , α / P β ′ ∈ NoSq( p ) b X β ′ , α ⊲ projection step end for for β ∈ NoSq ( p ) do coeff = λ β · (cid:18) − b β · Q α ∈ Cov β (cid:16) λ βα X β , α (cid:17) λ βα (cid:19) λ β X β , ← round-up (cid:16) coeff , b δ (cid:17) ⊲ adjust constant term end for C ← b − P β ∈ NoSq( p ) X β , return X , C In line 2, the function GP calls a GP solver to compute a e δ -approximation f X of (SONC).This approximation is then rounded in line 3 to a rational point c X with a prescribedmaximal relative error of b δ . The projection step from line 5 scales the entries of c X ,yielding P β X β , α = b α , for all β ∈ MoSq ( p ) , to satisfy the first set of equality constraintsof (SONC). In line 9, we round the coefficient up , with relative error b δ , so that we have X β , ≥ λ β · − b β · Y α ∈ Cov β (cid:18) λ βα X β , α (cid:19) λ βα λ β . As in Section 2.1, each p β := P α ∈ Cov β X β , α · x α + b β x β is a nonnegative circuit poly-nomial. Hence, C is a lower bound for p . Our assumption that every circuit polynomialcontains a constant term, is necessary to ensure that λ β = 0 for all β ∈ NoSq ( p ) in ourcomputations 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
XACT OPTIMIZATION VIA SONC AND SAGE 9
SAGE decomposition computed via REP, and applies a post-processing to find a certifiedlower bound .
Algorithm 3.2. optsage
Require: g = P ti =1 b i x α ( i ) ∈ Q [ x ] , rounding precision b δ ∈ Q > , precision parameter e δ ∈ Q > for the REP solver. Ensure:
Matrices c , ν of rational numbers, coefficients of the decomposition, certifiedlower bound C ∈ Q of g on R n . f ← g (exp x − exp( − x )) Build the ( n + 1) × t matrix Q with columns ( α (1) , , . . . , ( α ( t ) , e c , e ν ← REP ( f, e δ ) ⊲ Solve (SAGE) with accuracy e δ b c ← round (cid:16)e c , b δ (cid:17) , b ν ← round (cid:16)e ν , b δ (cid:17) ⊲ rounding step for j ∈ { , . . . , t } do LP ← n Q · ν ( j ) = , ν ( j ) \ j ≥ , k ν ( j ) − e ν ( j ) k ∞ ≤ b δ, ν ( j )1 ≥ b δ o ν ( j ) ← some element from LP ⊲ projection step c ( j ) \ j ← b c ( j ) \ j , c ( j ) j ← b j − · c ( j ) \ j end for for j ∈ { , . . . , t } do power ← − log ν ( j )1 − ν ( j )1 (cid:18) c ( j ) j − P i> ,i = j ν ( j ) i log ν ( j ) i ec ( j ) i (cid:19) c ( j )1 ← round-up (cid:16) exp (power) , b δ (cid:17) ⊲ adjust constant term end for C ← b − P tj =1 c ( j )1 return c , ν , C Given a polynomial g ( y ) = P tj =1 b j y α ( j ) , one could apply the change of variables y i :=exp x i when y ∈ R n> . Since this transformation is only valid on the nonnegative orthant,one workaround used in optsage is to define the signomial f ( x ) = g (exp x − exp( − x )) from line 1, in a such a way that a lower bound of f yields a lower bound of g .The REP function in line 3 calls an REP solver to compute a e δ -approximation ( e ν , e c ) of (SAGE). This approximation is then rounded to a rational point ( b ν , b c ) with a pre-scribed maximal relative error of b δ . The projection steps in line 7 and line 8 ensurethat ( ν , c ) satisfies exactly the linear equality constraints of (SAGE), i.e., Qν ( j ) = and P tj =1 c ( j ) = b . The first projection step boils down to exactly solve an LP withthe constraint that ν ( j )1 > , for all j = 1 , . . . , t , to ensure that further computa-tion in line 12 are well-defined. Note that this projection could be done while rely-ing on the pseudo-inverse of Q , but one obtains better practical results via this pro-cedure. To ensure that the relative entropy inequality constraints of (SAGE) are sat-isfied, the last step of optsage aims at finding c (1) j such that c ( j ) j ≥ D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) = P i> ,i = j ν ( j ) i log ν ( j ) i ec ( j ) i + ν ( j )1 log ν ( j )1 ec ( j )1 . Thus, one relies on the round-up procedure in line 12 XACT OPTIMIZATION VIA SONC AND SAGE 10 to compute c ( j )1 ≥ exp (cid:18) − log ν ( j )1 − ν ( j )1 (cid:18) c ( j ) j − P i> ,i = j ν ( j ) i log ν ( j ) i ec ( j ) i (cid:19)(cid:19) . Eventually, onehas P tj =1 c ( j ) i = b i , for all i > and P tj =1 c ( j )1 = b − C , which certifies that f − C ≥ on R n . We refer to appendix A.2 for an example of exact SAGE decomposition obtainedwith optsage . 4. Deciding Nonnegativity via SAGE
We denote by int ( C SAGE ) the interior of the cone C SAGE of SAGE signomials. A sig-nomial f = P tj =1 b j exp ( α ( j ) · x ) lies in int ( C SAGE ) if and only there is c (1) , . . . , c ( t ) , ν (1) , . . . , ν ( t ) ∈ R t such that t X j =1 c ( j ) = b , t X i =1 α ( i ) ν ( j ) i = , − · ν ( j ) \ j = ν ( j ) j , c ( j ) \ j , ν ( j ) \ j > , D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) < c ( j ) j , j = 1 , . . . , t . (INTSAGE-feas)Without the assumptions from Section 3, we state and analyze a decision algorithmto certify nonnegativity of signomials belonging to the interior int ( C SAGE ) of the SAGEcone. The resulting hybrid numeric-symbolic algorithm, called intsage , computes exactrational SAGE decompositions of such signomials.For complexity analysis purpose, we recall the following bound for the roots of univariatepolynomials with integer coefficients: Lemma 4.1. [Mig92, Theorem 4.2 (ii)]
Let f ∈ Z [ x ] of degree d , with coefficient bit sizebounded from above by τ . If f ( e x ) = 0 and e x = 0 , then τ +1 ≤ | e x | ≤ τ + 1 . Lemma 4.2.
Let f = P tj =1 b j exp ( α ( j ) · x ) ∈ int ( C SAGE ) of degree d with τ = τ ( f ) .Then, there exists N ∈ N such that for ε := 2 − N , f − ε P tj =1 exp ( α ( j ) · x ) ∈ C SAGE , with N ≤ τ ( ε ) ∈ O ( τ · (4 d + 6) n +3 ) .Proof. Since f ∈ int ( C SAGE ) , there are c (1) , . . . , c ( t ) , ν (1) , . . . , ν ( t ) ∈ R t such that P tj =1 c ( j ) = b , P ti =1 α ( i ) ν ( j ) i = , − · ν ( j ) \ j = ν ( j ) j , c ( j ) \ j , ν ( j ) \ j > and D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) < c ( j ) j ,for all j = 1 , . . . , t . Therefore, there exists N ∈ N such that for ε := 2 − N , one has D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) + ε < c ( j ) j , for all j = 1 , . . . , t . For all i, j = 1 , . . . , t , let us define ˚ b by ˚ b i := b i − ε , as well as ˚ c by ˚ c ( j ) i := c ( j ) i for i = j and ˚ c ( j ) j := c ( j ) j − ε . Note that ˚ b is the coefficient vector of f − ε P tj =1 exp ( α ( j ) · x ) . Then ˚ c (1) , . . . , ˚ c ( t ) , ν (1) , . . . , ν ( t ) satisfy (INTSAGE-feas) after replacing b by ˚ b , yielding the first claim.For the second claim, we start to perform the change of variable y i := exp x i , for all i = 1 , . . . , n , define g ( y , z ) := P tj =1 b j y α ( j ) − z P tj =1 y α ( j ) , for all y ∈ R n> . It is enoughto select ε = 2 − N such that ε ≤ inf y ∈ R n> P tj =1 b j y α ( j ) P tj =1 y α ( j ) . Let us consider the algebraic set V defined by: V := (cid:26) ( y , z ) ∈ R n +1 : g ( y , z ) = ∂g∂y = · · · = ∂g∂y n = 0 (cid:27) . XACT OPTIMIZATION VIA SONC AND SAGE 11
Using [MSED18, Proposition A.1], there exists a polynomial in Z [ z ] of degree less than ( d + 1) n +1 with coefficients of bit size less than τ · (4 d + 6) n +3 such that its set of realroots contains V . By Lemma 4.1, it is enough to take N ≤ τ · (4 d + 6) n +3 , yielding thedesired result. (cid:3) Algorithm intsage . We present our algorithm intsage computing exact rationalSAGE decompositions for signomials in int ( C SAGE ) . Algorithm 4.3. intsage
Require: f = P tj =1 b j exp ( α ( j ) · x ) ∈ int ( C SAGE ) , rounding precision b δ ∈ Q > , precisionparameter e δ ∈ Q > for the REP solver. Ensure:
Matrices c , ν of rational numbers. Build the ( n + 1) × t matrix Q with columns ( α (1) , , . . . , ( α ( t ) , Q + ← pseudoinv ( Q ) ok ← false while not ok do ( e c , e ν ) ← REP ( f, e δ ) b c ← round (cid:16)e c , b δ (cid:17) , b ν ← round (cid:16)e ν , b δ (cid:17) ⊲ rounding step for j ∈ { , . . . , t } do ⊲ projection step ν ( j ) ← ( I − Q + Q ) b ν ( j ) c ( j ) \ j ← b c ( j ) \ j , c ( j ) j ← b j − · c ( j ) \ j end for if for all j ∈ { , . . . , t } , ν ( j ) \ j , c ( j ) \ j ≥ , c ( j ) j ≥ D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) , then ok ← true ⊲ verification step else e δ ← e δ/ , b δ ← b δ/ end if end while return c , ν The routine pseudoinv in line 2 computes the pseudo-inverse of Q , i.e., a matrix Q + such that QQ + Q = Q . Next, we enter in the loop starting from line 4. The REP function calls an REP solver to compute a e δ -approximation ( e ν , e c ) of (INTSAGE-feas).The projection steps ensure that ( ν , c ) satisfies exactly the linear equality constraints of(SAGE-feas), i.e., Qν ( j ) = Q ( I − Q + Q ) ν ( j ) = Q − QQ + Q = and P tj =1 c ( j ) = b . Ifthe inequality constraints are not verified in line 11, the rounding-projection procedure isperformed 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 toSAGE decompositions of polynomials in int ( C SAGE ) . Lemma 4.4.
Let f = P tj =1 b j exp ( α ( j ) · x ) ∈ int ( C SAGE ) of degree d with τ = τ ( f ) . Let ε be as in Lemma 4.2. XACT OPTIMIZATION VIA SONC AND SAGE 12 (1) There exists a solution of ( ν , c ) of (INTSAGE-feas) and δ ∈ Q > such that δ ≤ , ( ν , c ) satisfies, D (cid:16) ν ( j ) \ j , e (cid:16) c ( j ) \ j + δ (cid:17)(cid:17) + ε ≤ c ( j ) j , P tj =1 c ( j ) = b , for all i, j =1 , . . . , t .(2) There exists a solution ( ν , c ) of (INTSAGE-feas) and δ ∈ Q > such that ( ν , c ) satisfies D (cid:16) ν ( j ) \ j + δ , e c ( j ) \ j (cid:17) + ε ≤ c ( j ) j , for all j = 1 , . . . , t .(3) There exists a solution ( ν , c ) of (INTSAGE-feas) and δ ∈ Q > such that ( ν , c ) satisfies D (cid:16) (1 + δ ) ν ( j ) \ j , e c ( j ) \ j (cid:17) + ε ≤ c ( j ) j , for all j = 1 , . . . , t .In each case, τ ( δ ) ∈ O ( τ · (4 d + 6) n +3 ) . A proof for this lemma is provided in appendix A.1.
Theorem 4.5.
Let f = P tj =1 b j exp ( α ( j ) · x ) ∈ int ( C SAGE ) of degree d and τ = τ ( f ) .There exist b δ and e δ of bit size less than O ( τ · (4 d + 6) n +3 ) , such that intsos ( f, b δ, e δ ) ter-minates and outputs a rational SAGE decomposition of f within O ( τ · (4 d + 6) n +3 t log t ) arithmetic operations.Proof. We first show that the loop of Algorithm intsage terminates with b δ and e δ of bitsize bounded by O ( τ · (4 d + 6) n +3 ) . Let ε be as in Lemma 4.2. When running the pro-cedure REP , one solves (SAGE-feas) at precision e δ , thus one finds an approximate solution ( e ν , e c ) such that k P tj =1 e c ( j ) − b k ∞ ≤ e δ , D (cid:16)e ν ( j ) \ j , e e c ( j ) \ j (cid:17) + ε ≤ e c ( j ) j + e δ , and k Q e ν ( j ) k ∞ ≤ e δ , forall j = 1 , . . . , t . After the rounding and projection steps, one obtains ν ( j ) = ( I − Q + Q ) b ν ( j ) and k b ν ( j ) − e ν ( j ) k ∞ ≤ b δ . Since α j ( i ) ≤ d , for all i, j = 1 , . . . , t , the bit size of the entries ofthe matrix Q is upper bounded by τ ( d ) . Thus, the pseudo-inverse Q + has rational entriesof bit size bounded by O ( t log t + t log d ) = O ( t log t ) , since the bit size is the same as forthe determinant length, see [BPR98, Corollary 8.13]. This implies that the bit size of thedifference between the entries of e ν and ν is upper bounded by O (cid:16) t log t + τ ( e δ ) + τ ( b δ ) (cid:17) .Similarly, the bit size of the difference between the entries of e c and c is upper boundedby O (cid:16) τ + τ ( b δ ) (cid:17) . By Lemma 4.4, one can perform any absolute or relative perturbationof e ν and e c , and still ensure that the resulting ( ν , c ) satisfies D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) ≤ c ( j ) j , if theperturbation is small enough with bit size at most O ( τ · (4 d + 6) n +3 ) . This implies thatone must choose e δ and b δ small enough, with the same upper bound on their bit sizes. Thesame reasoning applies to ensure that ν ( j ) \ j , c ( j ) \ j ≥ .Now, we give an upper bound on the number of arithmetic operations. For convexoptimization problems having barrier complexity equal to N , the standard interior-pointmethods compute a e δ -accurate solution in O ( τ ( e δ ) √ N log( N )) iterations; see e.g. [Ren01,Section 2.4]. For ( x, y, z ) ∈ R , the standard barrier complexity of a single relative entropyconstraint “ x log( x/y ) ≤ z, x, y ≥ ” is equal to 4. In addition, the barrier complexityof a set of constraints is upper bounded by the sum of the complexities of the individualconstraints. Therefore, the relative entropy formulation given in (SAGE-feas) has a barriercomplexity of N ≤ t . At each iteration of the interior-point method, one needs to solve XACT OPTIMIZATION VIA SONC AND SAGE 13 an LP involving t variables, which can be done within O ( t ) arithmetic operations. Thisyields the upper bound of O ( τ · (4 d + 6) n +3 t log t ) on the total number of arithmeticoperations required while calling intsage . All other arithmetic operations performed bythe algorithm have a negligible cost with respect to the REP procedure. (cid:3) 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 ourexperiment and on which instances the algorithms were tested. Afterwards, we discussour 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
POEM
POEM is open source, under GNU public license, and available at:
For our experiment,
POEM calls a range of further software and solvers for computingthe certificates. For the numerical solutions of SONC and SAGE, we use
CVXPY
ECOS
SymPy
Investigated Data
We carried out our experiments on 2020 randomly generated poly-nomials. The possible numbers of variables are n = 2 , , , , ; the degree takes values d = 6 , , , , , , and the number of terms can be t = 6 , , , , , , . Foreach combinations we create instances, where the number of negative terms is one of afew fixed ratios of t . In particular, the size of (SAGE) grows quadratically in t . We cre-ated the database using POEM , and it is available in full at the homepage cited above.Our instances are a subset of those from [SdW18a]. In that paper, we also describe theircreation 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, 4cores, 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 roundingin
Python we used a tolerance of ε = 2 − . The restriction t ≤ was chosen, sinceotherwise we already encounter problem in the numerical solution of (SAGE). The bound d < was chosen, because for large degree we had a significant increase in the memoryrequired to perform the rounding. Both thresholds were obtained experimentally.5.2. Evaluation of the Experiment.
In this section we present and evaluate the resultsof our experiment and highlight our most important findings, when investigating the com-putational 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(SAGE) shows that the size of the problem only depends on the number of terms t , butin the SAGE decomposition, the number of summands is the number of monomial non-squares. Most significantly, for more variables, our generating algorithm simply results XACT OPTIMIZATION VIA SONC AND SAGE 14 d = 8 , t = 20 n bit size time2 27723 5.823 23572 4.994 22965 4.838 5678 1.1210 1749 0.35 d = 10 , t = 24 n bit size time2 38872 8.263 33938 7.264 31278 6.348 7042 1.4210 3778 0.91 d = 10 , t = 30 n bit size time2 61198 12.413 57833 11.814 53596 10.828 13343 2.5510 6974 1.41 Table 1.
Dependency of the average bit size and the average running timeof optsage , with the number of variables, for fixed values of degree d andnumber of terms t ; For n = 8 we observe a drastic drop both in runningtime and bit size. t \ d × × × × Table 2.
Bit size (upper part) and running time (lower part) of optsage independency of the degree d and the number of terms t for up to 4 variables;A “ × ” indicates, that we do not have instances with these parameters inour data set.in a smaller number of these terms. Additionally, for n ≥ and d ≤ , most exponentslie 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 LPfrom line 6 of optsonc . Next, we have more equality constraints in this LP, which couldalso improve the running time. Lastly, the exponential upper bound is just the worstcase, which does not seem to actually happen among our examples. For some selectedparameters, we exhibit that behavior in Table 1. Dependency of bit size and running time of degree and terms
To illustratehow 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 XACT OPTIMIZATION VIA SONC AND SAGE 15 e - . . ,
30 107 1098 420 109 256 nu m b e r o f i n s t a n ce s Figure 1.
Number of instances where the difference of numerical lowerbound and exact lower bound lies in the given interval; note that the exactbound sometimes is better.those instances. The results are shown in Table 2. We can see that running time and bitsize roughly have a linear dependency. On the one hand, their growth is quadratic in thenumber of terms, which matches with the growth of the problem size in (SAGE). On theother hand, bit size and running time are basically unaffected by the degree. This showsthat the bound, given in the worst case analysis, usually is not met.
Quality of the rounding-projection
Our experiments verify that in the majority ofcases the symbolical lower bound does not diverge far from the numerical bound. Thedetailed distribution is shown in Figure 1. Most notably, in 30 instances, the exact lowerbound is even better than the numerical bound. In 81.9% of the instances, the exactbound differs by at most 0.001 from the numerical value. Only in 256 instances thedifference lies above 1. Thus, in the clear majority of examples, the lower bound in exactarithmetic does not differ much from the numerical bound. Also, among the instanceswith large difference, it can also be that the numerical solution actually lies far away froman exact solution. So it is unclear, whether a large difference is due to bad behavior ofthe numerical solution, or a large error in the rounding algorithm.
Rounding time versus solving time
In nearly every case the rounding proceduretakes 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 anaverage of 88.6%. However, one can implement the rounding procedure much closer tothe hardware level, instead of working in Python. Thus, we expect that these ratios canbe 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
XACT OPTIMIZATION VIA SONC AND SAGE 16 t bit size SONC bit size SAGE time SONC time SAGE6 432 1005 0.06 0.269 806 2696 0.19 0.6612 1261 5568 0.37 1.2920 2592 19203 0.64 4.0024 3826 32543 0.97 6.6630 5029 53160 1.34 10.5850 10622 167971 3.95 32.78 Table 3.
Comparison of running time and bit size of the certificates be-tween optsonc and optsage ; optsonc runs faster and has significantlysmaller certificates than optsage .best bound, that can be obtained via this approach. Still it shows, that for very largeinstances, SONC is the method of choice, when other approaches fail due to the problemsize. Comparison with SOS
For polynomials lying in the interior of the SOS conefrom [MSED18, Table 2], we performed preliminary experiments with optsage and optsage , which are currently unable to provide nonnegativity certificates. For bench-marks from our database with n ≥ and d ≥ , RealCertify often fails to provide SOScertificates. We plan to provide detailed experimental comparisons with SOS methods inthe future. 6.
Conclusion and Outlook
We make two main contributions in this paper. First, we present an algorithm todecide whether a given multivariate polynomial over the rationals lies in the interior ofthe SAGE cone. If that is the case, then the algorithm also computes a certificate in exactarithmetic. Additionally, we analyze the arithmetic complexity of the algorithm, whichis polynomial in the degree and the number of terms, and exponential in the number ofvariables. Second, we use our numerical methods to obtain lower bounds via SAGE andapply a single iteration of the rounding-projection method, to obtain an exactly certifiedlower bound. This method, we run on a large number of test cases. Based on theseexperiments, we draw the following conclusions.(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 iseven better than the numerical one.(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) For the investigated parameters, increasing the degree or the number of variablesdoes not increase the running time or the bit size. This also corresponds to thefact, that the size of the REP is independent of the degree and the number ofvariables.
XACT OPTIMIZATION VIA SONC AND SAGE 17 (4) For very large instances, SONC should be the first choice, to obtain a certifiedbound, since it runs significantly faster than the other methods.For future work, the most interesting development would be to have an REP-solver witharbitrary precision, so that we can actually implement intsage and compare it to similarapproaches. Furthermore, in a significant amount of instances, we encountered compu-tational problem, when calling optsonc or optsage . So we would like to increase therobustness of our implementation. Another issue, we have left out so far, is the presenceof exponents of monomial non-squares, which lie on a face of the Newton polytope, thatdoes not include the origin. These result in values λ ,j = 0 (SONC) or ν ( j )1 = 0 (SAGE),so our computation is undefined. However, these problems can be circumvented and weplan to do so in a future version of the software. Next, we plan to extend our frame-work to constrained problems and provide more detailed experimental comparisons withSOS-based approaches from [PP08, MSED18], as well as with methods based on criticalpoints and cylindrical algebraic decomposition. Given a polynomial in the interior of theSAGE cone, our decision algorithm intsage is linear with respect to the distance of thispolynomial to the border of the cone. In order to improve this bound, one could refinethe bit size analysis for this distance. A further theoretical aim would be to analyze theboolean running time of intsage , which requires to prove bit complexity estimates forrelative entropy optimization problems. References [BCSS12] L. Blum, F. Cucker, M. Shub, and S. Smale.
Complexity and real computation . SpringerScience & Business Media, 2012.[BGHP05] B. Bank, M. Giusti, J. Heintz, and L.-M. Pardo. Generalized polar varieties: geometry andalgorithms.
Journal of complexity , 21(4):377–412, 2005.[BKVH07] S. Boyd, S. Kim, L. Vandenberghe, and A. Hassibi. A tutorial on geometric programming.
Optim. Eng. , 8(1):67–127, 2007.[BPR98] S. Basu, R. Pollack, and M.-F. Roy. A new algorithm to find a point in every cell defined bya family of polynomials. In
Quantifier elimination and cylindrical algebraic decomposition .Springer-Verlag, 1998.[BV04] S. Boyd and L. Vandenberghe.
Convex optimization . Cambridge University Press, Cambridge,2004.[CHJL11] S. Chevillard, J. Harrison, M. Joldes, and C. Lauter. Efficient and accurate computation ofupper bounds of approximation errors.
Theoretical Computer Science , 412(16):1523 – 1543,2011.[Col75] G. E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decompos-tion. In
ATFL 2nd GI Conf. Kaiserslautern , pages 134–183, 1975.[CS16] V. Chandrasekaran and P. Shah. Relative Entropy Relaxations for Signomial Optimization.
SIAM J. Optim. , 26(2):1147–1173, 2016.[DB16] S. Diamond and S. Boyd. CVXPY: A Python-embedded modeling language for convex op-timization.
Journal of Machine Learning Research , 2016.[DCB13] A. Domahidi, E. Chu, and S. Boyd. ECOS: An SOCP solver for embedded systems. In
European Control Conference (ECC) , pages 3071–3076, 2013.[DIdW16] M. Dressler, S. Iliman, and T. de Wolff. An approach to constrained polynomial optimizationvia nonnegative circuit polynomials and geometric programming, 2016. To appear in theJournal of Symbolic Computation (MEGA 2017 special issue); see also arXiv:1602.06180 . XACT OPTIMIZATION VIA SONC AND SAGE 18 [DIdW17] M. Dressler, S. Iliman, and T. de Wolff. A Positivstellensatz for Sums of Nonnegative CircuitPolynomials.
SIAM J. Appl. Algebra Geom. , 1(1):536–555, 2017.[DPZ67] R. Duffin, E. Peterson, and C. Zener.
Geometric programming: Theory and application . JohnWiley & Sons, Inc., New York-London-Sydney, 1967.[dW15] T. de Wolff. Amoebas, nonnegative polynomials and sums of squares supported on circuits.
Oberwolfach Rep. , (23):1308–1311, 2015.[GL96] G. H. Golub and C. F. V. Loan.
Matrix Computations (3rd Ed.) . Johns Hopkins UniversityPress, Baltimore, MD, USA, 1996.[GM12] M. Ghasemi and M. Marshall. Lower bounds for polynomials using geometric programming.
SIAM J. Optim. , 22(2):460–473, 2012.[GM13] M. Ghasemi and M. Marshall. Lower bounds for a polynomial on a basic closed semialgebraicset using geometric programming, 2013. Preprint, arxiv:1311.3726 .[GV88] D. Grigoriev and N. Vorobjov. Solving systems of polynomials inequalities in subexponentialtime.
Journal of Symbolic Computation , 5:37–64, 1988.[HAB +
17] T. Hales, M. Adams, G. Bauer, D. T. Dat, J. Harrison, H. L. Truong, C. Kaliszyk, V. Ma-gron, S. Mclaughlin, N. T. Thang, N. Q. Truong, T. Nipkow, S. Obua, J. Pleso, J. Rute,A. Solovyev, T. T. H. An, T. N. Trung, T. T. Diep, J. Urban, V. K. Ky, and R. Zumkeller.A Formal Proof of the Kepler Conjecture.
Forum of Mathematics, Pi , 5, 2017.[IdW16] S. Iliman and T. de Wolff. Amoebas, nonnegative polynomials and sums of squares supportedon circuits.
Res. Math. Sci. , 3:3:9, 2016.[JvMG12] D. Joyner, O. Čertík, A. Meurer, and B. E. Granger. Open source computer algebra systems:Sympy.
ACM Commun. Comput. Algebra , 45(3/4):225–234, January 2012.[Las01] J.-B. Lasserre. Global optimization with polynomials and the problem of moments.
SIAMJournal on Optimization , 11(3):796–817, 2001.[Las06] J.-B. Lasserre. Convergent SDP-Relaxations in Polynomial Optimization with Sparsity.
SIAM Journal on Optimization , 17(3):822–843, 2006.[Las10] J.-B. Lasserre.
Moments, positive polynomials and their applications , volume 1 of
ImperialCollege Press Optimization Series . Imperial College Press, London, 2010.[Lau09] M. Laurent. Sums of squares, moment matrices and optimization over polynomials. In
Emerg-ing applications of algebraic geometry , volume 149 of
IMA Vol. Math. Appl. , pages 157–270.Springer, New York, 2009.[LTY17] J.-B. Lasserre, K.-C. Toh, and S. Yang. A bounded degree SOS hierarchy for polynomialoptimization.
EURO Journal on Computational Optimization , 5(1-2):87–117, 2017.[MAGW15] V. Magron, X. Allamigeon, S. Gaubert, and B. Werner. Formal proofs for Nonlinear Opti-mization.
Journal of Formalized Reasoning , 8(1):1–24, 2015.[MCD17] V. Magron, G. Constantinides, and A. Donaldson. Certified Roundoff Error Bounds UsingSemidefinite Programming.
ACM Trans. Math. Softw. , 43(4):1–34, 2017.[MCW18] R. Murray, V. Chandrasekaran, and A. Wierman. Newton polytopes and relative entropyoptimization. arXiv preprint arXiv:1810.01614 , 2018.[MD18a] V. Magron and M. S. E. Din. On Exact Polya and Putinar’s Representations. In
ISSAC’18:Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Compu-tation . ACM, New York, NY, USA, 2018.[MD18b] V. Magron and M. S. E. Din. RealCertify: a Maple package for certifying non-negativity. In
ISSAC’18: Proceedings of the 2018 ACM International Symposium on Symbolic and Alge-braic Computation . ACM, New York, NY, USA, 2018.[MDS18] V. Magron, M. S. E. Din, and M. Schweighofer. Algorithms for weighted sum of squaresdecomposition of non-negative univariate polynomials.
Journal of Symbolic Computation ,2018.[Mig92] M. Mignotte.
Mathematics for Computer Algebra . Springer-Verlag New York, Inc., New York,NY, USA, 1992.
XACT OPTIMIZATION VIA SONC AND SAGE 19 [MSED18] V. Magron and M. Safey El Din. On Exact Polya, Hilbert-Artin and Putinar’s Representa-tions. arXiv preprint arXiv:1811.10062 , 2018.[Oxl11] J. Oxley.
Matroid theory , volume 21 of
Oxford Graduate Texts in Mathematics . Oxford Uni-versity Press, Oxford, second edition, 2011.[Par00] P. A. Parrilo.
Structured Semidefinite Programs and Semialgebraic Geometry Methods inRobustness and Optimization . PhD thesis, California Inst. Tech., 2000.[PP08] H. Peyrl and P. Parrilo. Computing sum of squares decompositions with rational coefficients.
Theoretical Computer Science , 409(2):269–281, 2008.[Ren88] J. Renegar. A faster PSPACE algorithm for deciding the existential theory of the reals. In
Foundations of Computer Science, 1988., 29th Annual Symposium on , pages 291–295. IEEE,1988.[Ren01] J. Renegar.
A Mathematical View of Interior-Point Methods in Convex Optimization . SIAM,2001.[RTAL13] C. Riener, T. Theobald, L. J. Andrén, and J.-B. Lasserre. Exploiting Symmetries in SDP-Relaxations for Polynomial Optimization.
Mathematics of Operations Research , 38(1):122–141, 2013, http://dx.doi.org/10.1287/moor.1120.0558.[SdW18a] H. Seidler and T. de Wolff. An experimental comparison of sonc and sos certificates forunconstrained optimization. arXiv preprint arXiv:1808.08431 , 2018.[SdW18b] H. Seidler and T. de Wolff. POEM: Effective methods in polynomial optimization, version0.1.1.0(a). , April 2018.[WKKM06] H. Waki, S. Kim, M. Kojima, and M. Muramatsu. Sums of Squares and Semidefinite Pro-gramming Relaxations for Polynomial Optimization Problems with Structured Sparsity.
SIAM Journal on Optimization , 17(1):218–242, 2006.[WLT18] T. Weisser, J.-B. Lasserre, and K.-C. Toh. Sparse-BSOS: a bounded degree SOS hierarchy forlarge scale polynomial optimization with sparsity.
Mathematical Programming Computation ,10(1):1–32, 2018.
XACT OPTIMIZATION VIA SONC AND SAGE 20
Appendix A. Appendix
A.1.
Proof of Lemma 4.4.
We start with the first claim. By Lemma 4.2, there exist e ν , e c and ε ∈ Q > , with τ ( ε ) ∈ O ( τ · (4 d + 6) n +3 ) , such that ( e ν , e c ) satisfies P tj =1 c ( j ) = b , e c ( j ) \ j , e ν ( j ) \ j > , D (cid:16)e ν ( j ) \ j , e e c ( j ) \ j (cid:17) + ε ≤ e c ( j ) j , for all j = 1 , . . . , t . Let us define δ := ε t − and c such that c ( j ) i := e c ( j ) i + ε for all i = j = 1 , . . . , t , and c ( j ) j := e c ( j ) j − tδ , for j = 1 , . . . , t . Thus,one has P tj =1 c ( j ) = b . For all j = 1 , . . . , t , one has ν ( j ) \ j > . Combining this togetherwith the fact that the log function is increasing yields D (cid:16) ν ( j ) \ j , e (cid:16) c ( j ) \ j + δ (cid:17)(cid:17) + ε ≤ D (cid:16) ν ( j ) \ j , e c ( j ) \ j (cid:17) ≤ e c ( j ) j − ε c ( j ) j . Using that for all d ≥ , t ≤ (cid:18) n + dn (cid:19) = ( n + d ) · · · ( d + 1) n ! = (cid:18) dn (cid:19) (cid:18) dn − (cid:19) · · · (1 + d ) ≤ d n − (1 + d ) ≤ d n , one has τ ( t ) ∈ O ( n log ( d )) . Since τ ( δ ) ≤ τ ( ε ) + τ ( t ) , we obtain the first claim.To prove the second claim, we rely on the following three auxiliary inequalities.For ν, δ such that < ν < δ ≤ , one has(A.1) ( ν + δ ) log( ν + δ ) ≤ , since < ν + δ < and the function x x log x is negative on (0 , .For ν, δ such that < δ ≤ ν , one has(A.2) ( ν + δ ) log(1 + δ/ν ) ≤ δ , since ν + δ ≤ ν and log(1 + δ/ν ) ≤ δ/ν .For each ν, δ, c > such that ν ≥ , one has(A.3) δ log (cid:16) νec (cid:17) ≤ δ max n , ν log (cid:16) νec (cid:17)o . Indeed, if ν ≤ ec , the left hand side is less than 0. Otherwise, ν ≥ implies that log( νec ) ≤ ν log( νec ) .Now, by the first claim, there exist ε c , ε ∈ Q > , with τ ( ε c ) , τ ( ε ) ∈ O ( τ · (4 d + 6) n +3 ) ,and e ν , e c satisfying e ν ( j ) \ j > , e c ( j ) \ j > ε c , and D (cid:16)e ν ( j ) \ j , e e c ( j ) \ j (cid:17) + ε ≤ e c ( j ) j , for all j = 1 , . . . , t .For each δ ∈ Q > with δ ≤ , and all j = 1 , . . . , t , one has: D (cid:16)e ν ( j ) \ j + δ , e e c ( j ) \ j (cid:17) = X i = j ( e ν ( j ) i + δ ) log e ν ( j ) i + δe e c ( j ) i . We give an upper bound of each summand of the right hande side, depending on the valueof e ν ji , for all i, j = 1 , . . . , t and i = j . Note that − log( ec ji ) ≤ log( eε c ) , for all i, j = 1 , . . . , t and i = j . XACT OPTIMIZATION VIA SONC AND SAGE 21 If < e ν ji ≤ δ ≤ , one has by (A.1) ( e ν ( j ) i + δ ) log e ν ( j ) i + δe e c ( j ) i ≤ ( e ν ( j ) i + δ ) log( e ν ( j ) i + δ ) − ( e ν ( j ) i + δ ) log( e e c ( j ) i ) ≤ δ log( eε c ) . If δ ≤ e ν ji ≤ , one has ( e ν ( j ) i + δ ) log e ν ( j ) i + δe e c ( j ) i ≤ e ν ( j ) i log e ν ( j ) i e e c ( j ) i + δ log e ν ( j ) i e e c ( j ) i + ( e ν ( j ) i + δ ) log δ e ν ( j ) i ! ≤ e ν ( j ) i log e ν ( j ) i e e c ( j ) i + 2 δ log( eε c ) + 2 δ , where we use the fact that δ log e ν ( j ) i ≤ and bound the last term of the right hand sidevia (A.2).If e ν ji ≥ , we write the first inequality as in the former case and obtain ( e ν ( j ) i + δ ) log e ν ( j ) i + δe e c ( j ) i ≤ e ν ( j ) i log e ν ( j ) i e e c ( j ) i + δ max ( , e ν ( j ) i log e ν ( j ) i e e c ( j ) i ) + 2 δ , where we rely on (A.3) to bound the second term and the fact that δ + e ν ( j ) i ≤ e ν ( j ) i togetherwith log (cid:18) δ e ν ( j ) i (cid:19) ≤ δ e ν ( j ) i to bound the last term.In the worst case, we obtain D (cid:16)e ν ( j ) \ j + δ , e e c ( j ) \ j (cid:17) ≤ (1 + 2 δ ) D (cid:16)e ν ( j ) \ j , e e c ( j ) \ j (cid:17) + 2 δ ( t − | log( eε c ) | + 1) ≤ e c ( j ) j + 2 δ e c ( j ) j − (1 + 2 δ ) ε + 2 δt | log( eε c ) |≤ e c ( j ) j − (1 + 2 δ ) ε + 2 δ ( | b j | + t | log( eε c ) | ) , using e c ( j ) j = b j − · e c ( j ) \ j ≤ b j , thus e c ( j ) j ≤ max { , b j } ≤ | b j | .To ensure that D (cid:16)e ν ( j ) \ j + δ , e e c ( j ) \ j (cid:17) + ε ≤ e c ( j ) j , it is sufficient to have − (1 + 4 δ ) ε +4 δ ( | b j | + t | log( eε c ) | ) ≤ , which is guaranteed by selecting the largest positive ra-tional δ such that δ ≤ ε | b j | + t | log( eε c ) |− ε ) . Since τ ( | b j | ) ≤ τ , τ ( t ) ∈ O ( d log n ) , | log( eε c ) | ) ≤ τ ( ε c ) , and τ ( ε c ) , τ ( ε ) ∈ O ( τ · (4 d + 6) n +3 ) , one can select δ with bit sizeat most O ( τ · (4 d + 6) n +3 ) .The proof of the third claim is very similar and we omit it for the sake of conciseness. XACT OPTIMIZATION VIA SONC AND SAGE 22
A.2.
An Example of Exact SAGE Decomposition.
Let f ( x ) = 277 − x + 159 x x + 275 x − x x x + 23 x x x + 338 x x + 166 x x x − x x x − x x x + 74 x x x + 268 x x . Our optimization algorithm optsage returns ( ν , c ) corresponding to the following exactrational SAGE decomposition: f ( x ) = P j =1 f j ( x ) , where f j is the polynomial withcoefficient vector c ( j ) for j ∈ { , . . . , } , f j = 0 for j ∈ { , , , , , } , and ν (2) = (cid:18) , − , , , , , , , , , , (cid:19) f = 70505161 − x + 51364821929347737990176335433100237218699413 x x + 128403578802267056788513490375489389305583 x +8306281050730062422200197261503152705103 x x x + 16345503627618597428155608650885 x x x + 15189310926109008012048313292258664398599 x x +68872628669902760563068799030848037634611567 x x x + 584984347065145672102911615960676281785795 x x x + 351808814017845417152094682478713 x x x +91006698046824055380305977850685287 x x x + 7052856072901897195672038441901303429467107 x x ν (5) = (cid:18) , , , , − , , , , , , , (cid:19) f = 48956532 + 39758344941359351362516496563648929 x + 4282065321967986445380570176335433100237218699413 x x + 5340068537046034265265107923003915114444664 x − x x x + 624998409184656940119485631121730177 x x x + 250559824383122027902128012048313292258664398599 x x +3177316622322934865488003068799030848037634611567 x x x + 798626491949623771695220582323192135256357159 x x x + 489457358510061917516251456284047436139 x x x +6200073611296304448855380305977850685287 x x x + 98753640749442855613080672038441901303429467107 x x ν (6) = (cid:18) , , , , , − , , , , , , (cid:19) f = 116124 + 1757518123308140626189133159801243971 x + 20993411202657969962132460176335433100237218699413 x x + 230789345356723546855107923003915114444664 x +349712347144442759232022200197261503152705103 x x x − x x x + 1903199631733604622132012048313292258664398599 x x +1881403227048869602158083068799030848037634611567 x x x + 27412383525960721370432102911615960676281785795 x x x + 6164401897310612801805717364894159571 x x x +23295511841744620296055380305977850685287 x x x + 239384185705575737209071672038441901303429467107 x x ν (8) = (cid:18) , , , , , , , − , , , , (cid:19) f = 165536389 + 801946315133135394046189133159801243971 x + 936097366115717531112852176335433100237218699413 x x + 96283800124071875502013490375489389305583 x +16851627907815928010784022200197261503152705103 x x x + 261920871973041930119485631121730177 x x x + 62219051633182216810752012048313292258664398599 x x − x x x + 1345600783114821896667952102911615960676281785795 x x x + 2332419503709903114405417152094682478713 x x x +184738012166616458823055380305977850685287 x x x + 140920372909036228167589650672038441901303429467107 x x XACT OPTIMIZATION VIA SONC AND SAGE 23 ν (9) = (cid:18) , , , , , , , , − , , , (cid:19) f = 1811076 + 32665647808648050562063044386600414657 x + 945568751423587451936415176335433100237218699413 x x + 19555286702192671237513490375489389305583 x +10859834179375768514928022200197261503152705103 x x x + 261624143096996295119485631121730177 x x x + 87681114365023463949264612048313292258664398599 x x +93303623560361276969456003068799030848037634611567 x x x − x x x + 19290963496573463552016251456284047436139 x x x +133955133103661412452055380305977850685287 x x x + 17532007073614874704137330672038441901303429467107 x x ν (10) = (cid:18) , , , , , , , , , − , , (cid:19) f = 1572744 + 54245379022417729126189133159801243971 x + 880139855950527034906380176335433100237218699413 x x + 57106936799639827353013490375489389305583 x +3316914557864893313328022200197261503152705103 x x x + 228162076394260360119485631121730177 x x x + 4869784075970617279361612048313292258664398599 x x +15878203082552055693731141022933010282679211537189 x x x + 5915050338968076718077720582323192135256357159 x x x − x x x +61625442506078228280055380305977850685287 x x x + 21315777567587124730738350672038441901303429467107 x x Victor Magron, CNRS LAAS, 7 avenue du Colonel Roche, F-31031 Toulouse Cédex 4,France
E-mail address : [email protected] Henning Seidler, Technische Universität Berlin, Institut für Mathematik, Straße des17. Juni 136, 10623 Berlin, Germany
E-mail address : [email protected] Timo de Wolff, Technische Universität Berlin, Institut für Mathematik, Straße des17. Juni 136, 10623 Berlin, Germany
E-mail address ::