Algorithms for Weighted Sums of Squares Decomposition of Non-negative Univariate Polynomials
AAlgorithms for Weighted Sums of Squares Decomposition ofNon-negative Univariate Polynomials
Victor Magron Mohab Safey El Din Markus Schweighofer November 11, 2018
Abstract
It is well-known that every non-negative univariate real polynomial can be written as the sumof two polynomial squares with real coefficients. When one allows a weighted sum of finitely many squares instead of a sum of two squares, then one can choose all coefficients in the representation tolie in the field generated by the coefficients of the polynomial. In particular, this allows an effectivetreatment of polynomials with rational coefficients.In this article, we describe, analyze and compare both from the theoretical and practical pointsof view, two algorithms computing such a weighted sums of squares decomposition for univariatepolynomials with rational coefficients.The first algorithm, due to the third author relies on real root isolation, quadratic approximationsof positive polynomials and square-free decomposition but its complexity was not analyzed. Weprovide bit complexity estimates, both on the runtime and the output size of this algorithm. Theyare exponential in the degree of the input univariate polynomial and linear in the maximum bitsizeof its complexity. This analysis is obtained using quantifier elimination and root isolation bounds.The second algorithm, due to Chevillard, Harrison, Joldes and Lauter, relies on complex rootisolation and square-free decomposition and has been introduced for certifying positiveness of poly-nomials in the context of computer arithmetics. Again, its complexity was not analyzed. We providebit complexity estimates, both on the runtime and the output size of this algorithm, which are poly-nomial in the degree of the input polynomial and linear in the maximum bitsize of its complexity.This analysis is obtained using Vieta’s formula and root isolation bounds.Finally, we report on our implementations of both algorithms and compare them in practice onseveral application benchmarks. While the second algorithm is, as expected from the complexityresult, more efficient on most of examples, we exhibit families of non-negative polynomials for whichthe first algorithm is better.
Keywords: non-negative univariate polynomials, Nichtnegativstellens¨atze, sums of squares decompo-sition, root isolation, real algebraic geometry.
Given a subfield K of R and a non-negative univariate polynomial f ∈ K [ X ], we consider the problem ofproving the existence and computing the weighted sums of squares decompositions of f with coefficientsalso lying in K . CNRS Verimag; 700 av Centrale 38401 Saint-Martin d’H`eres, France Sorbonne Universit´es, UPMC Univ. Paris 06, CNRS, Inria Paris Center, LIP6, Equipe PolSys, F-75005, Paris, France Fachbereich Mathematik und Statistik, Universit¨at Konstanz, 78457 Konstanz, Germany a r X i v : . [ c s . S C ] J un eyond the theoretical interest of this question, finding certificates of non-negative polynomials is manda-tory in many application fields. Among them, one can mention the stability proofs of critical controlsystems often relying on Lyapunov functions [30], the certified evaluation of mathematical functions inthe context of computer arithmetics (see for instance [4]), the formal verification of real inequalities [11]within proof assistants such as Coq [6] or
Hol-light [12] ; in these situations the univariate case isalready an important one. In particular, formal proofs of polynomial non-negativity can be handled withsums of squares certificates. These certificates are obtained with tools available outside of the proof assis-tants and eventually verified inside. Because of the limited computing power available inside such proofassistants, this is crucial to devise algorithms that produce certificates, whose checking is computationallyreasonably simple. In particular, we would like to ensure that such algorithms output sums of squarescertificates of moderate bitsize and ultimately with a computational complexity being polynomial withrespect to the input.
Related Works.
Decompositing non-negative univariate polynomials into sums of squares has a longstory ; very early quantitative aspects like the number of needed squares have been studied. For thecase K = Q , Landau shows in [18] that for every non-negative polynomial in Q [ X ], there exists adecomposition involving a weighted sum of (at most) eight polynomial squares in Q [ X ]. In [28], Pourchetimproves this result by showing the existence of a decomposition involving only a weighted sum of (atmost) five squares. This is done using approximation and valuation theory ; extracting an algorithm fromthese tools is not the subject of study of this paper.More recently, the use of semidefinite programming for computing sums-of-squares certificates of non-negativity for polynomials has become very popular since [19, 26]. Given a polynomial f of degree n ,this method consists in finding a real symmetric matrix G with non-negative eigenvalues (a semidefinitepositive matrix), such that f ( x ) = v ( x ) T Gv ( x ), where v is the vectors of monomials of degree less than n/
2. Hence, this leads to solve a so-called Linear Matrix Inequality and one can rely on semidefiniteprogramming (SDP) to find the coefficients of G . This task can be delegated to an SDP solver (e.g.SeDuMi, SDPA, SDPT3 among others). An important technical issue arises from the fact that suchSDP solvers are most of the time implemented with floating-point double precision. More accuratesolvers are available (e.g. SDPA-GMP [25]). However, note that these solvers always compute numericalapproximations of the algebraic solution to the semidefinite program under consideration. Hence, theyare not sufficient to provide algebraic certificates of posivity with rational coefficients. Hence, a processfor making exact and with rational coefficients the computed numerical approximations of sums of squarescertificates is needed. This issue has been tackled in [27, 17]. The certification scheme described in [21]allows to obtain lower bounds of non-negative polynomials over compact sets. However, despite theirefficiency, these methods do not provide any guarantee to output a rational solution to a Linear MatrixInequality when it exists (and especially when it is far from the computed numerical solution).A more systematic treatment of this problem has been brought by the symbolic computation community.Linear Matrix Inequalities can be solved as a decision problem over the reals with polynomial constraintsusing the Cylindrical Algebraic Decomposition algorithm [5] or more efficient critical point methods (seee.g. [1] for complexity estimates and [16, 9] for practical algorithms). But using such general algorithmsis an overkill and dedicated algorithms have been designed for computing exactly algebraic solutionsto Linear Matrix Inequalities [13, 14]. Computing rational solutions can also be considered thanks toconvexity properties [31]. In particular, the algorithm in [10] can be used to compute sums of squarescertificates with rational coefficients fpr a non-negative univariate polynomial of degree n with coefficientsof bit size bounded by τ using τ O (1) O ( n ) bit operations at most (see [10, Theorem 1]).For the case where K is an arbitrary subfield of R , Schweighofer gives in [32] a new proof of the existenceof a decomposition involving a sum of (at most) n polynomial squares in K [ X ]. This existence proofcomes together with a recursive algorithm to compute such decompositions. At each recursive step, the2lgorithm performs real root isolation and quadratic approximations of positive polynomials. Later on, asecond algorithm is derived in [4, Section 5.2], where the authors show the existence of a decompositioninvolving a sum of (at most) n +3 polynomial squares in K [ X ]. This algorithm is based on approximatingcomplex roots of perturbed positive polynomials.These both latter algorithms were not analyzed despite the fact that they were implemented and used.An outcome of this paper is a bit complexity analysis for both of them, showing that they have bettercomplexities than the algorithm in [10], the second algorithm being polynomial in n and τ . Notation for complexity estimates.
For complexity estimates, we use the bit complexity model. Foran integer b ∈ Z \{ } , we denote by τ ( b ) := log ( | b | ) + 1 the bitsize of b , with the convention τ (0) := 1.We write a given polynomial f ∈ Z [ X ] of degree n ∈ N as f = (cid:80) ni =0 b i X i , with b , . . . , b n ∈ Z . Inthis case, we define (cid:107) f (cid:107) ∞ := max ≤ i ≤ n | b i | and, using a slight abuse of notation, we denote τ ( (cid:107) f (cid:107) ∞ )by τ ( f ). Observe that when f has degree n , the bit size necessary to encode f is bounded by nτ ( f )(when storing each coefficients of f ). The derivative of f is f (cid:48) = (cid:80) ni =1 ib i X i − . For a rational number q = bc , with b ∈ Z , c ∈ Z \{ } and gcd( b, c ) = 1, we denote max { τ ( b ) , τ ( c ) } by τ ( q ). For two mappings g, h : N l → R , the expression “ g ( v ) = O ( h ( v ))” means that there exists a integer b ∈ N such that forall v ∈ N l , g ( v ) ≤ bh ( v ). The expression “ g ( v ) = ∼ O ( h ( v ))” means that there exists a integer c ∈ N suchthat for all v ∈ N l , g ( v ) ≤ h ( v ) log ( h ( v )) c . Contributions.
We present and analyze two algorithms, denoted by univsos1 and univsos2 , allowingto decompose a non-negative univariate polynomial f of degree n into sums of squares with coefficientslying in any subfield K of R . To the best of our knowledge, there was no prior complexity estimate for theoutput of such certification algorithms based on sums of squares in the univariate case. We summarizeour contributions as follows: • We describe in Section 3 the first algorithm, called univsos1 in the sequel. It was originally givenin [32, Chapter 2] ; Section 3 can be seen as a partial English translation of this text written inGerman since some proofs have been significantly simplified. In the same section, we analyze its bitcomplexity. When the input is a polynomial of degree n with integer coefficients of maximum bitsize τ , we prove that Algorithm univsos1 uses ∼ O (( n ) n τ ) boolean operations and returns polynomialsof bitsize bounded by O (( n ) n τ ). This is not restrictive: when f ∈ Q [ X ], one can multiply itby the least common multiple of the denominators of its coefficients and apply our statement forpolynomials in Z [ X ]. • We describe in Section 4 the second algorithm univsos2 , initially given in [4, Section 5.2]. Wealso analyze its bit complexity. When the input is a univariate polynomial of degree n with integercoefficients of maximum bitsize τ , we prove that Algorithm univsos2 returns a sums of squaredecompositions of n + 3 polynomials with coefficients of bitsize bounded by O ( n + n τ ) using ∼ O ( n + n τ ) boolean operations. • Both algorithms are implemented within the univsos tool. The first release of univsos is a Maplelibrary, is freely available and is integrated in the RAGlib (Real Algebraic Library) Maple package .The scalability of the library is evaluated in Section 5 on several non-negative polynomials issuedfrom the existing literature. Despite the significant difference of theoretical complexity betweenthe two algorithms, numerical benchmarks indicate that both may yield competitive performancew.r.t. specific sub-classes of problems. https://github.com/magronv/univsos Preliminaries
We first recall the proof of the following classical result for non-negative real-valued univariate polynomials(see e.g. [29, Section 8.1]).
Theorem 2.1.
Let f ∈ R [ X ] be a non-negative univariate polynomial, i.e. f ( x ) ≥ , for all x ∈ R .Then, f can be written as the sum of two polynomial squares in R [ X ] .Proof. Without loss of generality, one can assume that f is monic, i.e. the leading coefficient (nonzerocoefficient of highest degree) of f is 1. Then we decompose f as follows in C [ x ]: f = (cid:89) j ( X − a j ) r j (cid:89) k (( X − ( b k + ic k ))( X − ( b k − ic k ))) s k , with a j , b k , c k ∈ R , r j , s k ∈ N > , a j standing for the real roots of f and ( b k ± ic k ) standing for the complexconjugate roots of f . Since f is non-negative, all real roots must have even multiplicity r j , yielding theexistence of polynomials g, q, r ∈ R [ X ] satisfying the following: g = (cid:89) j ( X − a j ) r j , q + ir = (cid:89) k ( X − ( b k + ic k )) s k q − ir = (cid:89) k ( X − ( b k − ic k )) s k . Then, one has f = g ( q + ir )( q − ir ) = g ( q + r ) = ( gq ) + ( gr ) , which proves the claim.Let K be a field and g ∈ K [ X ]. One says that g is a square-free polynomial, when there is no primeelement p ∈ K [ X ] such that p divides g . Let now f ∈ K [ X ] − { } . A decomposition of f of the form f = ag g . . . g nn with a ∈ K and normalized pairwise coprime square-free polynomials g , g , . . . , g n iscalled a square-free decomposition of f in K [ X ].We recall the following useful classical bounds. Lemma 2.2. [2, Corollary 10.12] If p ∈ Z [ X ] and q ∈ Z [ X ] divides p in Z [ X ] , then one has τ ( q ) ≤ deg q + τ ( p ) + log (deg p + 1) . The algorithm of Yun [35] (also described in [7, Algorithm 14.21]) allows to compute a square-freedecomposition of polynomials with coefficients in a field of characteristic 0.
Lemma 2.3. [7, § f ∈ Z [ X ] with degree at most n with coefficient bitsize upper bounded by τ .Then the square-free decomposition of f using Yun’s Algorithm [35] can be computed using an expectednumber of ∼ O ( n τ ) boolean operations. Lemma 2.4. [24, § K be a field of characteristic 0 and L a field extension of K . Thesquare-free decomposition in L [ X ] of any polynomial f ∈ K [ X ] − { } is the same as the square-freedecomposition of f in K [ X ] . Any polynomial f ∈ K [ X ] − { } which is a square-free polynomial in K [ X ] is also a square-free polynomial in L [ X ] . The following lemma allows to obtain upper bounds on the magnitudes of the roots of a univariatepolynomial.
Lemma 2.5. (Cauchy Bound [3])
Let K be an ordered field. Let a , . . . , a n ∈ K with a n (cid:54) = 0 . Let x ∈ K such that (cid:80) ni =0 a i x i = 0 . Then, one has: | x | ≤ max (cid:26) , | a || a n | + · · · + | a n − || a n | (cid:27) . Lemma 2.6. [24] Let f ∈ Z [ X ] of degree n , with coefficient bitsize upper bounded by τ . If f ( x ) = 0 and x (cid:54) = 0 , then τ +1 ≤ | x | ≤ τ + 1 . The real (resp. complex) roots of a polynomial can be approximated using root isolation techniques. Tocompute the real roots one can use algorithms based on Uspensky’s method relying on Descartes’s ruleof sign, see e.g. [2, Chap. 10] for a general description of real root isolation algorithms.
Lemma 2.7. [22, Theorem 5] Let f ∈ Z [ X ] with degree at most n with coefficient bitsize upper boundedby τ . Isolating intervals (resp. disks) of radius less than − κ for all real (resp. complex) roots of f canbe computed in ∼ O ( n + n τ + nκ ) boolean operations. Vieta’s formulas provide relations between the coefficients of a polynomial and signed sums and productsof the complex roots of this polynomial:
Lemma 2.8. (Vieta’s formulas [8])
Let K be an ordered field. Given a polynomial f = (cid:80) i =0 n a i X i ∈ K [ X ] with a n (cid:54) = 0 with (not necessarily distinct) complex roots z , . . . , z n , one has for all j = 1 , . . . , n : (cid:88) ≤ i < ···
Let K be an ordered field. Let g = aX + bX + c ∈ K [ X ] with a, b, c ∈ K and a (cid:54) = 0 . Then, g can be rewritten as g = a (cid:0) X + b a (cid:1) + (cid:16) c − b a (cid:17) . Moreover, when g is non-negative over K , one has a > and c − b a ≥ .Proof. The decomposition of g is straightforward. Assume that g is non-negative over K . Remark that c − b a = g (cid:0) − b a (cid:1) ; hence since we assume that g is non-negative over K we deduce that c − b a ≥ a > a <
0. Then, this impliesthat for all x ∈ K , one has (cid:0) x + b a (cid:1) ≤ − a (cid:16) c − b a (cid:17) . Thus, there exists C ∈ K such that x ≤ C , foreach x ∈ K . This implies in particular for x = 2 that 4 ≤ C and for x = C that C ≤ C , thus C ≤ ≤ C ≤
1, yielding a contradiction.Let f ∈ K [ X ] be a square-free polynomial which is non-negative over R . Then, f is positive over R ,otherwise f would have at least one real root, implying that f would be neither a square-free polynomialin R [ X ] nor a square-free polynomial in K [ X ], according to Lemma 2.4. We want to find a polynomial g ∈ K [ X ] which fulfills the following conditions:(i) deg g ≤ g is non-negative over R ,(iii) f − g is non-negative over R , 5iv) f − g has a root t ∈ K .Assume that Property (i) holds. Then the existence of a sum of squares decomposition for g is ensuredfrom Property (ii). Property (iii) implies that h = f − g has only non-negative values over R . The aimof Property (iv) is to ensure the existence of a root t ∈ K of h , which is stronger than the existence of areal root. Note that the case where the degree of h = f − g is less than the degree of f occurs only whendeg f = 2. In this latter case, we can rely on Lemma 3.1.Now, we investigate the properties of a polynomial g ∈ K [ X ], which fulfills conditions (i)-(iii) and (iv)with t ∈ K . Using Property (i) and Taylor Decomposition, we obtain g = g ( t ) + g (cid:48) ( t )( X − t ) + c ( X − t ) for some c ∈ K . By Property (iv), one has g ( t ) = f ( t ). In addition, Property (iii) yields f ( x ) − g ( x ) ≥ f ( t ) − g ( t ), for all x ∈ K , which implies that ( f − g ) (cid:48) ( t ) = 0 and g (cid:48) ( t ) = f (cid:48) ( t ). By Property (ii),the quadratic polynomial g ( X + t ) = f ( t ) + f (cid:48) ( t ) X + cX has at most one root. This implies that thediscriminant of g , namely f (cid:48) ( t ) − cf ( t ) cannot be positive, thus one has c ≥ f (cid:48) ( t ) f ( t ) (since f ( t ) > g satisfying (i)-(iii) and (iv) for some t ∈ K , one necessarily has g = f t,c with f (cid:48) ( t ) f ( t ) ≤ c ∈ K and f t,c = f ( t ) + f (cid:48) ( t )( X − t ) + c ( X − t ) .In this case, one also has that the polynomial g = f t,c (cid:48) , with c (cid:48) = f (cid:48) ( t ) f ( t ) , fulfills (i)-(iii) and (iv). Indeed,(i) and (iv) trivially hold. Let us prove that (ii) holds: when deg f t,c (cid:48) = 0, then g = f ( t ) ≥ f t,c (cid:48) = 2, then g has a single root − f (cid:48) ( t )2 c and the minimum of g is g (cid:16) − f (cid:48) ( t )2 c (cid:17) = 0. The inequalities f t,c (cid:48) ≤ f t,c ≤ f over R yield (iii).Therefore, given f ∈ K [ X ] with f positive over R , we are looking for t ∈ K such that the inequality f ≥ f t holds over R , with f t := f ( t ) + f (cid:48) ( t )( X − t ) + f (cid:48) ( t ) f ( t ) ( X − t ) ∈ K [ X ] . The main problem is to ensure that t lies in K . If we choose t to be a global minimizer of f , then f t would be the constant polynomial min { f ( x ) | x ∈ R } . The idea is then to find t in the neighborhood ofa global minimizer of f . The following lemma shows that the inequality f t ≤ f can be always satisfiedfor t in some neighborhood of a local minimizer of f . Lemma 3.2.
Let f ∈ R [ X ] . Let a be a local minimizer of f and suppose that f ( a ) > . For all t ∈ R with f ( t ) (cid:54) = 0 , let us define the polynomial f t : f t := f ( t ) + f (cid:48) ( t )( X − t ) + f (cid:48) ( t ) f ( t ) ( X − t ) ∈ R [ X ] . Then, there exists a neighborhood U ⊂ R of a such that the inequality f t ( x ) ≤ f ( x ) holds for all ( t, x ) ∈ U × U .Proof. Set n := deg f . It is easy to see that we can suppose without loss of generality that a is the originand that f (0) = 1. Because of the Taylor formula f = n (cid:88) k =0 f ( k ) ( t ) k ! ( X − t ) k , we have f − f t = n (cid:88) k =2 f ( k ) ( t ) k ! ( X − t ) k − f (cid:48) ( t ) f ( t ) ( X − t ) = ( X − t ) (cid:32) n (cid:88) k =2 f ( k ) ( t ) k ! ( X − t ) k − − f (cid:48) ( t ) f ( t ) (cid:33) t ∈ R with f ( t ) (cid:54) = 0. Let h be the bivariate polynomial defined as follows: h := f ( T ) (cid:32) n (cid:88) k =2 f ( k ) ( T ) k ! ( X − T ) k − (cid:33) − f (cid:48) ( T ) ∈ R [ T, X ] . Let us prove that ( a, a ) is a local minimizer of h .Since f (0) = 1, there exists c (cid:54) = 0, α ∈ N and g ∈ R [ X ] such that f − cX α + X α +1 g . Therefore,lim x → f ( x ) − cx α = 1. Since f − R , one concludes that c > α is even. Let usconsider the lowest homogeneous part H of h , that is the sum of all monomials of lowest degree involvedin h . The lowest homogeneous part of f (cid:48) ( T ) is c ( α − T α − with degree 2 α − (cid:80) nk =2 f ( k ) ( T ) k ! ( X − T ) k − is c (cid:80) nk =2 (cid:0) αk (cid:1) T α − k ( X − T ) k − with degree α −
2. Then H = c n (cid:88) k =2 (cid:18) αk (cid:19) T α − k ( X − T ) k − , and thus( X − T ) H = c (( T + ( X − T )) α − T n − nT α − ( X − T )) = c ( X α − αT α − X + ( α − T α ) . Since lim (cid:107) ( x,t ) (cid:107)→ h ( x,t ) H ( x,t ) = 1, it is enough to prove that H is positive except at the origin in order toshow that ( a, a ) = (0 ,
0) is a local minimizer of h . Let us consider ( t, x ) ∈ R \ { } and show that H ( t, x ) >
0. If t = x , we have H ( t, x ) = H ( x, x ) = (cid:0) α (cid:1) x α − >
0. If t (cid:54) = x , then it is enough to show that( x − t ) H ( t, x ) = c ( x α − αt α − x + ( α − t α ) >
0. This is clear if t = 0 since c > α is even. Nowsuppose that t (cid:54) = 0 and define ξ := xt (cid:54) = 1. Then one has t − α ( x − t ) H ( t, x ) = c ( ξ α − αξ + α − > H follows from the fact that the univariate polynomial r := X α − αX + α − r (cid:48) = αX α − − α . The positivity of H implies that ( a, a ) is a local minimizer of h .Let us define q ( X, T ) := ( X − T ) h . Combining the fact that ( a, a ) is a local minimizer of the twopolynomials h , ( X − T ) and the fact that h ( a, a ) = f ( a ) f (cid:48)(cid:48) ( a ) − f (cid:48) ( a ) = 0, we conclude that ( a, a ) isalso a local minimizer of q . Since f ( x ) − f t ( x ) = f ( t ) q ( x, t ), this yields the existence of a neighborhood O ⊂ R of ( a, a ) such that the inequality f − f t ≥ x, t ) ∈ O . Since there exists someneighborhood U ⊂ R of a , such that the rectangle U × U is included in O , this proves the initial claim.Lemma 3.2 states the existence of a neighborhood U of a local minimizer of f such that the inequality f t ( x ) ≤ f ( x ) holds for all ( x, t ) ∈ U × U . Now, we show that with such a neighborhood U of the smallestglobal minimizer of f , the inequality f t ( x ) ≤ f ( x ) holds for all t ∈ U and for all x ∈ R . Proposition 3.3.
Let f ∈ R [ X ] with deg f > . Assume that f is positive over R . Then, there existsa smallest global minimizer a of f and a positive (cid:15) ∈ R such that for all t ∈ R with a − (cid:15) < t < a , thequadratic polynomial f t , defined by f t := f ( t ) + f (cid:48) ( t )( X − T ) + f (cid:48) ( t ) f ( t ) ( X − T ) ∈ R [ X ] , satisfies f t ≤ f over R .Proof. The existence of a is straightforward. First, we handle the case when deg f = 2. Using TaylorDecomposition of f at t , one obtains f = f ( t ) + f (cid:48) ( t )( X − T ) + f (cid:48)(cid:48) ( t )2 ( X − T ) . Since f has no real root,7he discriminant of f is negative, namely f (cid:48) ( t ) − f ( t ) f (cid:48)(cid:48) ( t )2 <
0. It implies that f (cid:48) ( t ) f ( t ) < f (cid:48)(cid:48) ( t ) , ensuringthat the inequality f t ≤ f holds over R .In the sequel, we assume that deg f >
2. We can find a neighborhood U as in Lemma 3.2 and withoutloss of generality, let us suppose that U = [ a − (cid:15) , a + (cid:15) ] for some positive (cid:15) , so that f (cid:48) has no real rootin [ a − (cid:15) , a ). Then, the inequality f t ( x ) ≤ f ( x ) holds for all x, t ∈ U . Next, we write f − f t = (cid:80) ni =0 a it x i ,with a it ∈ R and n = deg f > U → R : t (cid:55)→ C t := max (cid:26) , | a t || a nt | , . . . , | a ( n − t || a nt | (cid:27) . Note that the Cauchy bound (Lemma 2.5) implies that for all t ∈ U , all real roots of f − f t lie in[ − C t , C t ]. In addition, the closed interval domain U is compact, implying that the range values of thefunction U → R : t (cid:55)→ C t are bounded. Let C ∈ R with C ≥ C t for all t ∈ U . Then, for all t ∈ U ,all real roots of f − f t lie in the interval [ − C, C ] and we can assume without loss of generality that − C < a − (cid:15) < a < a + (cid:15) < C . Let us define M := min { f ( x ) | x ∈ [ − C, a − (cid:15) ] } . By definition, a is theglobal minimizer of f , ensuring that f ( a ) < M . For all t ∈ [ a − (cid:15) , a ), the quadratic polynomial f t hasone real root N t := − f ( t ) f (cid:48) ( t ) + t . When t ∈ [ a − (cid:15) , a ) converges to a , then f (cid:48) ( t ) < − f ( t ) converges towards − f ( a ) <
0. Thus, the corresponding limit of N t is + ∞ . In addition, f t ( − C )tends to f a ( − C ) = f ( a ) < M . Therefore, there exists some (cid:15) ∈ (0 , (cid:15) ] such that for all t ∈ ( a − (cid:15), a ), onehas N t ∈ [ C, ∞ ) and f t ( − C ) < M . For all t ∈ ( a − (cid:15), a ), we partition R into five intervals and prove thatthe inequality f t ≤ f holds on each interval: • The inequality f t ≤ f holds over ( −∞ , − C ]: it comes from the fact that f t ( − C ) < M ≤ f ( − C )and the fact f − f t has no real root in ( −∞ , − C ]. • The inequality f t ≤ f holds over ( − C, a − (cid:15) ]: f t is monotonically decreasing over ( −∞ , N t ]. Sinceone has − C < a − (cid:15) < C ≤ N t , then f t is monotonically decreasing over ( − C, a − (cid:15) ]. This impliesthat for all x ∈ ( − C, a − (cid:15) ], one has f t ( x ) ≤ f t ( − C ) < M ≤ f ( x ). • The inequality f t ≤ f holds over [ a − (cid:15) , a ): it follows from the fact that [ a − (cid:15) , a ) ⊆ U . • The inequality f t ≤ f holds over [ a, C ): f t is monotonically decreasing over ( −∞ , N t ]. Since onehas a < C ≤ N t , then f t is monotonically decreasing over [ a, C ). Since a is a global minimizer of f and a ∈ U , one has f t ( x ) ≤ f t ( a ) ≤ f ( a ) ≤ f ( x ) for all x ∈ [ a, C ) . • The inequality f t ≤ f holds over [ C, ∞ ]: the claim is implied by the fact that f t ( N t ) = 0 < f ( N t ), N t ∈ [ C, ∞ ] together with the fact that f − f t has no real root in [ C, ∞ ]. Proposition 3.4.
Let K be a subfield of R and f ∈ K [ X ] with deg f = n ≥ . Then f is non-negativeon R if and only if f is a weighted sum of n polynomial squares in K [ X ] , i.e. there exist a , . . . , a n ∈ K ≥ and g , . . . , g n ∈ K [ X ] such that f = (cid:80) ni =0 a i g i .Proof. The if part is straightforward. For the other direction, assume that f is non-negative on R and n is even. The proof is by induction over n . The base case n = 2 follows from Lemma 3.1. For theinduction case, let us consider n ≥ f is not a square-free polynomial, we show that f is a weighted sum of n − f = gh , for some polynomials g, h ∈ K [ X ] with deg g ≤ deg f −
2. This gives g ( x ) = f ( x ) h ( x ) ≥ x ∈ R such that h ( x ) (cid:54) = 0. Since h has a finite number of real roots, g is non-negative on R . Using8 nput: non-negative polynomial f ∈ K [ X ] of degree n ≥
2, with K a subfield of R Output: pair of lists of polynomials ( h list , q list ) with coefficients in K h list := [ ], q list := [ ]. while deg f > do ( g, h ) := sqrfree ( f ) (cid:46) f = gh if deg h > then h list := h list ∪ { h } , q list := q list ∪ { } , f := g else f t := parab ( f ) ( g, h ) := sqrfree ( f − f t ) h list := h list ∪ { h } , q list := q list ∪ { f t } , f := g end done h list := h list ∪ { } , q list := q list ∪ { f } return h list , q list Figure 1: univsos1 : algorithm to compute SOS decompositions of non-negative univariate polynomials.the induction hypothesis, g is a weighted sum of n − f is also a weightedsum of n − f is a square-free polynomial, then f has no real root, which implies by Lemma 2.4 that f is neithera square-free polynomial in K [ X ] nor in R [ X ]. Thus, f is positive on R . Using Proposition 3.3, thereexists some t ∈ K ( K is dense in R ) and a quadratic polynomial f t ∈ K [ X ] such that the inequalities0 ≤ f t ( x ) ≤ f ( x ) holds for all x ∈ R and f t ( t ) = f ( t ). The polynomial f − f t has degree n , takes onlynon-negative values. In addition ( f − f t )( t ) = 0, thus f − f t is not a square-free polynomial. Hence, we arein the above case, implying that f − f t is a weighted sum of n − f t is a weighted sum of 2 polynomial squares, implying that f is a weighted sum of n polynomial squares,as requested. univsos1 The global minimizer a is a real root of f (cid:48) ∈ K [ X ]. Therefore, by using root isolation techniques [2,Chap. 10], one can isolate all the real roots of f (cid:48) in distinct intervals with bounds in K . Such techniquesrely on applying successive bisections, so that one can arbitrarily reduce the width of every interval andsort them w.r.t. their lower bounds. Eventually, we apply this procedure to find a sequence of elementsin K converging from below to the smallest global minimizer of f in order to find a suitable t . We denoteby parab ( f ) the corresponding procedure which returns the polynomial f t := f (cid:48) ( t ) f ( t ) ( X − t ) + f (cid:48) ( t )( X − t ) + f ( t ) such that t ∈ K and f ≥ f t over R .Algorithm univsos1 , depicted in Figure 1, takes as input a polynomial f ∈ K [ X ] of even degree n ≥ parab performing root isolation(see Step ). The second one is denoted by sqrfree and performs square-free decomposition: for a givenpolynomial f ∈ K [ X ], sqrfree ( f ) returns two polynomials g and h in K [ X ] such that f = gh . When f is square-free, the procedure returns g = f and h = 1 (in this case deg h = 0). As in the proof ofProposition 3.4, this square-free decomposition procedure is performed either on the input polynomial f (Step ) or on the non-negative polynomial ( f − f t ) (Step ). The output of Algorithm univsos1 is apair of lists of polynomials in K [ X ], allowing to retrieve an SOS decomposition of f . By Proposition 3.4the length of all output lists, denoted by r , is bounded by n/
2. If we note h r , . . . , h the polynomials9elonging to h list and q r , . . . , q the positive quadratic polynomials belonging to q list , one obtainsthe following Horner-like decomposition: f = h r (cid:0) h r − ( h r − ( . . . ) + q r − ) + q r − (cid:1) + q r . Thus, each positivequadratic polynomial q i being a weighted SOS polynomial, this yields a valid weighted SOS decompositionfor f . Example 1.
Let us consider the polynomial f := X + X − X − X + X + 2 ∈ Q [ X ].We describe the different steps performed by Algorithm univsos • The polynomial f is square-free and the algorithm starts by providing the value t = − f . With f ( t ) = and f (cid:48) ( t ) = − , one obtains f − = ( − X + ) . • Next, after obtaining the square-free decomposition f ( X ) − f − = ( X + 1) g , the same procedureis applied on g . One obtains the value t = 1 as an approximation of the smallest minimizer of g and g = ( − X + ) . • Eventually, one obtains the square-free decomposition g ( X ) − g = ( X − h with h = ( X − ).Overall, Algorithm univsos h list = [1 , X +1 , , X − ,
0] and q list = [ ( − X + ) , , ( − X + ) , , ( X − )], yielding the following weighted SOS decomposition: f := (cid:32) ( X + 1) (cid:18) ( X − (cid:16)
116 ( X − (cid:17) + 502920237293 ( − X + 88411167640 ) (cid:19)(cid:33) + 7201397 ( − X + 271360 ) . In the sequel, we analyze the complexity of Algorithm univsos1 in the particular case K = Q . Weprovide bounds on the bitsize of related SOS decompositions as well as bounds on the arithmetic costrequired for computation and verification. Lemma 3.5.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n and τ be an upper boundon the bitsize of the coefficients of f . When applying Algorithm univsos1 to f , the sub-procedure parab outputs a polynomial f t such that τ ( t ) = O ( n τ ) .Proof. Let us consider the set S ⊆ Q defined by: S := { t ∈ Q | ∀ x ∈ R , f ( t ) + f (cid:48) ( t ) f ( t )( x − t ) + f (cid:48) ( t ) ( x − t ) ≤ f ( t ) f ( x ) } . The polynomial involved in S has degree 2 n , with maximum bitsize of the coefficients upper bounded by2 τ . Thanks to the complexity analysis of the quantifier elimination procedure described in [2, § S can be described by polynomials with maximum bitsize coefficients upper bounded by O ( n τ ).Since t is a rational root of one of these polynomials, the rational zero theorem [34] implies that τ ( t ) = O ( n τ ). Lemma 3.6.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n and τ be an upper boundon the bitsize of the coefficients of f . When applying Algorithm univsos1 to f , the sub-procedure parab outputs a polynomial f t such that τ ( f t ) = O ( n τ ) . Moreover, there exist polynomials ˆ f , ˆ f t , g ∈ Z [ X ] suchthat ˆ f − ˆ f t = ( X − t ) g and τ ( g ) = O ( n τ ) . roof. One can write f t = M ( t ) X + M ( t ) X + M ( t ) with M ( t ) := f (cid:48) ( t ) f ( t ) ,M ( t ) := 2 f (cid:48) ( t )(2 f ( t ) − tf (cid:48) ( t ))4 f ( t ) ,M ( t ) := (2 f ( t ) − tf (cid:48) ( t )) f ( t ) , and (cid:107) f t (cid:107) ∞ ≤ max { M ( t ) , | M ( t ) | , M ( t ) } . One has 0 ≤ M ( t ) = f t (0) ≤ f (0) ≤ (cid:107) f (cid:107) ∞ .In addition, 0 ≤ M ( t ) + M ( t ) + M ( t ) = f t (1) ≤ f (1) ≤ ( n + 1) (cid:107) f (cid:107) ∞ and 0 ≤ M ( t ) − M ( t ) + M ( t ) = f t ( − ≤ f ( − ≤ ( n + 1) (cid:107) f (cid:107) ∞ . Thus, one has M ( t ) + | M ( t ) | + M ( t ) ≤ ( n + 1) (cid:107) f (cid:107) ∞ , wich impliesthat (cid:107) f t (cid:107) ∞ ≤ ( n + 1) (cid:107) f (cid:107) ∞ .Now let us note t = t t , with t ∈ Z , t ∈ Z \{ } , t and t being coprime. Let us define the polynomialsˆ f ( X ) := t n f ( t ) f ( X ) and ˆ f t ( X ) := t n f ( t ) f t ( X ). By writing f ( X ) = (cid:80) ni =0 a i X i , one has t n f ( t ) = (cid:80) ni =0 a i t i t n − i ≤ (cid:107) f (cid:107) ∞ | t | i | t | n − i and τ ( ˆ f ) ≤ τ + τ ( t n ). By Lemma 3.5, one has τ ( ˆ f ) = O ( n τ ).The polynomials ˆ f ( X ) , ˆ f t ( X ) are polynomials in Z [ X ] and since (cid:107) ˆ f t (cid:107) ∞ ≤ ( n + 1) (cid:107) ˆ f (cid:107) ∞ , the triangularinequality (cid:107) ˆ f − ˆ f t (cid:107) ∞ ≤ (cid:107) ˆ f (cid:107) ∞ + (cid:107) ˆ f t (cid:107) ∞ implies that τ ( ˆ f − ˆ f t ) ≤ log ( n + 2) + τ ( ˆ f ). In addition, τ ( f t ) ≤ τ ( ˆ f t ) + τ ( t n f ( t )) = O ( n τ ).As in the proof of Proposition 3.4, one has ( ˆ f − ˆ f t )( t ) = 0 which allows to write the square-free decom-position of the polynomial ˆ f − ˆ f t ∈ Z [ X ] as ˆ f − ˆ f t = ( X − t ) g , with g ∈ Z [ X ]. By Lemma 2.2, one has τ ( g ) ≤ n − τ ( ˆ f − ˆ f t ) + log ( n + 1) ≤ n − ( n + 2) + τ ( ˆ f ) = O ( n τ ), which concludes theproof. Theorem 3.7.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n = 2 k and τ be an upperbound on the bitsize of the coefficients of f . Then the maximum bitsize of the coefficients involved in theSOS decomposition of f obtained with Algorithm univsos1 is upper bounded by O (( k !) τ ) = O (( n ) n τ ) .Proof. With k = n/ f , Algorithm univsos1 generates, in the worstcase scenario, two sequences of polynomials f k , . . . , f ∈ Z [ X ] , q k , . . . , q ∈ Z [ X ] as well as rationalnumbers t k , . . . , t ∈ Q such that f k = f , t i = t i t i , with t i ∈ Z , t i ∈ Z \{ } and t ii f i ( t i ) f i − q i = ( X − t i ) f i − , i = 2 , . . . , k . (1)From Lemma 3.6, for all i = 2 , . . . , k , one has τ ( f i − ) = O ( i τ ( f i )). This yields τ ( f ) = O (cid:0) ( k !) τ ( f ) (cid:1) .Using Stirling’s formula, we obtain k ! ≤ √ πk (cid:0) ke (cid:1) k and ( k !) ≤ √ π k (cid:0) ke (cid:1) k , where e denotesthe Euler number. Since k ≤ e k for each integer k ≥ <
3, one has ( k !) ∈ O ( k k ), yielding τ ( f i ) = O (( n ) n τ ), for all i = 1 , . . . , k . Similarly, we obtain τ ( q i ) = O (( n ) n τ ), for all i = 1 , . . . , k .Finally, using Lemma 3.5, one has τ ( t i ) = O ( i τ ( f i )), yielding the desired result. Theorem 3.8.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n = 2 k and τ be an upperbound on the bitsize of the coefficients of f . Then, on input f , Algorithm univsos1 runs in boolean time ∼ O ( k · ( k !) τ ) = ∼ O (cid:18)(cid:16) n (cid:17) n τ (cid:19) . roof. For i = 2 , . . . , k we obtain each polynomial f i − as in the proof of Theorem 3.7 by comput-ing the square-free decomposition of the polynomial t ii f i ( t i ) f i − q i . It follows by Lemma 2.3 thatthe polynomial f i − can be computed using an expected number of ∼ O ( i · i τ ( f i )) boolean opera-tions. The number of boolean operations to compute all polynomials f , . . . , f k − is thus bounded by ∼ O (cid:0) k · k τ ( lf ) + ( k − ( k − k τ ( lf ) + · · · + ( k !) τ ( lf ) (cid:1) .For each i = 2 , . . . , k , the bitsize of the rational number t i is upper bounded by O ( i τ ( f i )). There-fore, t i can be computed by approximating the roots of f (cid:48) i with isolating intervals of radius lessthan 2 − i τ ( f i ) . By Lemma 2.7, the corresponding computation cost is ∼ O ( i τ ( f i )) boolean opera-tions. The number of boolean operations to compute all rational numbers t , . . . , t k is bounded by ∼ O (cid:0) k · k τ ( lf ) + ( k − ( k − k τ ( lf ) + · · · + ( k !) τ ( lf ) (cid:1) .In addition, one has k · k + ( k − ( k − k + · · · + ( k !) = ( k !) (cid:80) ki =1 1( i !) ≤ k · ( k !) . UsingStirling’s formula, we obtain k · ( k !) ≤ √ π k (cid:0) ke (cid:1) k . Since k / ≤ e k for each integer k ≥
1, weobtain the announced complexity.For a given polynomial f of degree 2 k , one can check the correctness of the SOS decomposition obtainedwith Algorithm univsos1 by evaluating this SOS polynomial at 2 k + 1 distinct points and compare theresults with the ones obtained while evaluating f at the same points. Theorem 3.9.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n = 2 k and τ be an upperbound on the bitsize of the coefficients of f . Then one can check the correctness of the SOS decompositionof f obtained with Algorithm univsos1 within ∼ O ( k · ( k !) τ ) = ∼ O (( n n τ ) boolean operations.Proof. From [7, Corollary 8.27], the cost of polynomial multiplication in Z [ X ] of degree less than n = 2 k with coefficients of bitsize upper bounded by B is bounded by ∼ O ( k · B ). By Theorem 3.7, the maximalbitsize of the coefficients of the SOS decomposition of f obtained with Algorithm univsos1 is upperbounded by B = O (( k !) τ ). Let us consider 2 k + 1 distinct integers, with maximal bitsize upper boundedby log n . Therefore, from [7, Corollary 10.8], the cost of the evaluation of this decomposition at the2 k + 1 points can be performed using at most ∼ O ( k · ( k !) τ ) boolean operations, the desired result. Remark 3.10.
Let f k = f ∈ Z [ X ]. Under the strong assumption that all polynomials f k , . . . , f involvedin Algorithm univsos1 have at least one integer global minimizer, then Algorithm univsos1 has a poly-nomial complexity. Indeed, in this case, q i = f i ( t i ), τ ( t i ) = O ( τ ( f i )) and τ ( f i − ) = O (2( i −
1) + τ ( f i )),for all i = 2 , . . . , k . Hence, the maximal bitsize of the coefficients involved in the SOS decomposition of f is upper bounded by O ( k + τ ) and this decomposition can be computed using an expected numberof ∼ O ( k + k τ ) boolean operations. Here, we recall the algorithm given in [4, Section 5.2]. The description of this algorithm, denotedby univsos2 , is given in Figure 2. 12 nput: non-negative polynomial f ∈ K [ X ] of degree n ≥
2, with K a subfield of R , ε ∈ K such that0 < ε < f n , precision δ ∈ N for complex root isolation Output: list c list of numbers in K and list s list of polynomials in K [ X ] ( p, h ) := sqrfree ( f ) (cid:46) f = p h n (cid:48) := deg p , k := n (cid:48) / p ε := p − ε (cid:80) ki =0 X i while has real roots ( p ε ) do ε := ε , p ε := p − ε (cid:80) ki =0 X i done ε := ε ( s , s ) := sum two squares ( p ε , δ ) (cid:96) := f n , u := p ε − (cid:96)s − (cid:96)s , u − := 0, u k +1 := 0 (cid:46) u = (cid:80) k − i =0 u i X i while ε < min ≤ i ≤ k (cid:8) | u i +1 | − u i + | u i − | (cid:9) do δ := 2 δ , ( s , s ) := sum two squares ( p ε , δ ), u := p ε − (cid:96)s − (cid:96)s done c list := [ (cid:96), (cid:96) ], s list := [ h s , h s ] for i = 0 to k − do c list := c list ∪ {| u i +1 |} , s list := s list ∪ { h ( X i +1 + sgn ( u i +1 )2 X i ) } c list := c list ∪ { ε − | u i +1 | + u i − | u i − |} , s list := s list ∪ { h X i } done return c list ∪ { ε + u n − | u n − |} , s list ∪ { h X k } Figure 2: univsos2 : algorithm to compute SOS decompositions of non-negative univariate polynomials. univsos2
Given a subfield K of R and a non-negative polynomial f = (cid:80) ni =0 f i X i ∈ K [ X ] of degree n = 2 k , one firstobtains the square-free decomposition of f , yielding f = p h with p > R (see Step ).Then the ideais to find a positive number ε > K such that the perturbed polynomial p ε ( X ) := p ( X ) − ε (cid:80) ki =0 X i is also positive on R (see [4, Section 5.2.2] for more details). This number is computed thanks to the loopgoing from Step to Step and relies on the auxiliary procedure has real roots which checks whetherthe polynomial p ε has real roots using root isolation techniques. As mentioned in [4, Section 5.2.2], thenumber ε is divided by 2 again to allow a margin of safety (Step ).Note that one can always ensure that the leading coefficient (cid:96) := p n of p is the same as the leadingcoefficients f n of the input polynomial f .We obtain an approximate rational sums of squares decomposition of the polynomial p ε with the auxiliaryprocedure sum two squares (Step ), relying on an arbitrary precision complex root finder. RecallingTheorem 2.1, this implies that the polynomial p can be approximated as close as desired by the weightedsum of two polynomial squares in Q [ X ], that is (cid:96)s + (cid:96)s .Thus there exists a remainder polynomial u := p ε − (cid:96)s − (cid:96)s with coefficients of arbitrary small magnitude (as mentioned in [4, Section 5.2.3]). The magnitudeof the coefficients converges to 0 as the precision δ of the complex root finder goes to infinity. The precisionis increased thanks to the loop going from Step to Step until a condition between the coefficients of u and ε becomes true, ensuring that ε (cid:80) ki =0 X i + u ( X ) also admits a weighted SOS decomposition. Formore details, see [4, Section 5.2.4].The reason why Algorithm univsos2 terminates is the following: at first, one can always find a sufficiently13mall perturbation ε such that the perturbed polynomial p ε remains positive. Next, one can always findsufficiently precise approximations of the complex roots of p ε ensuring that the error between the initialpolynomial p and the approximate SOS decomposition is compensated thanks to the perturbation term.The output of Algorithm univsos2 are a list of numbers in K and a list of polynomials in K [ X ], allowingto retrieve a weighted SOS decomposition of f . The size r of both lists is equal to 2 k + 3 = n (cid:48) + 3 ≤ n + 3.If we note c r , . . . , c the numbers belonging to c list and s r , . . . , s the polynomials belonging to s list ,one obtains the following SOS decomposition: f = c r s r + · · · + c s . Example 2.
Let us consider the same polynomial f := X + X − X − X + X + 2 ∈ Q [ X ] asin Example 1. We describe the different steps performed by Algorithm univsos • The polynomial f is square-free so we obtain p = f (Step ). After performing the loop fromStep to Step , Algorithm univsos ε = at Step as well as the polynomial p ε := p − (1 + X + X + X ) which has no real root. • Next, after increasing three times the precision in the loop going from Step to Step , the resultof the approximate root computation yields s = X − X and s = 7 X − X − .Applying Algorithm univsos
2, we obtain the following two lists of size 6 + 3 = 9: c list = (cid:20) , , , , , , , , (cid:21) , s list = (cid:20) X − X, X − X − , , X, X , X , X + 12 , X ( X −
12 ) , X ( X + 12 ) (cid:21) , yielding the following weighted SOS decomposition: f = 132 (cid:18) X − X (cid:19) + 132 (cid:18) X − X − (cid:19) + 91315360 + 73192160 X + 71152 X + 132 X + 797680 (cid:18) X + 12 (cid:19) + 1576 X (cid:18) X − (cid:19) . First, we need the following auxiliary result:
Lemma 4.1.
Let p ∈ Z [ X ] be a positive polynomial over R , with deg p = n = 2 k and τ be an upperbound on the bitsize of the coefficients of p . Then, one has inf x ∈ R p ( x ) > ( n τ ) − n +2 − n log n − nτ . Proof.
Denoting by τ (cid:48) the maximum bit size of the coefficients of p (cid:48) and instantiating α = inf x ∈ R p ( x )with a global minimizer of p , Q with p and A with p (cid:48) in the third item of [23, Lemma 3.2], one obtains.inf x ∈ R p ( x ) > ( n τ ) − n +2 − nτ (cid:48) Now, remark that τ (cid:48) ≤ log n + τ . Using this inequality in the one above allows to conclude.14 emma 4.2. Let p ∈ Z [ X ] be a positive polynomial over R , with deg p = n = 2 k and let τ be an upperbound on the bitsize of the coefficients of p . Then there exists a positive integer N such that for ε := N ,the polynomial p ε := p − ε (cid:80) ki =0 X i is positive over R and N = τ ( ε ) ≤ O ( n log n + nτ ) .Proof. Let us first consider the polynomial r := p − p n (cid:80) ki =0 X i . Using [2, Corollary 10.4], the absolutevalue of each real root of the polynomial r is bounded by n τ ( r ) ≤ n τ . By defining R := 2 n τ , it followsthat the polynomial r is positive for all | x | > R . In addition, for all positive integer N and ε = N , onehas ε ≤ ≤ p n and p ε = p − ε (cid:80) ki =0 X i ≥ p − p n (cid:80) ki =0 X i = r , which implies that the polynomial p ε is also positive for all | x | > R . Since R = 2 n τ >
1, one has 1 + R · · · + R n < nR n . Let us choose thesmallest positive integer N such that nR n ≤ N inf | x |≤ R p . This implies that ε < inf | x |≤ R p R ··· + R n , ensuringthat the polynomial p ε is also positive for all | x | ≤ R . Applying Lemma 4.1, we obtain the followingupper bound: 2 N ≤ nR n ( n τ ) n − n log n + nτ = n n +1 n nτ ( n τ ) n − n log n + nτ . The announced estimate follows straightforwardly.In the sequel, we denote by z , . . . , z n the (not necessarily distinct) complex roots of the polynomial p ε .Assuming that we approximate each complex root with a relative precision of δ , we write ˆ z , . . . , ˆ z n forthe approximate complex root values satisfying ˆ z i = z i (1 + e i ), with | e i | ≤ − δ , for all i = 1 , . . . , n . Theorem 4.3.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n and τ be an upper boundon the bitsize of the coefficients of f . Then the maximal bitsize of the coefficients involved in the weightedSOS decomposition of f obtained with Algorithm univsos2 is upper bounded by O ( n + n τ ) .Proof. Let p be the square-free part of the polynomial f (see Step of Algorithm univsos2 ). Then byusing Lemma 2.2, one has τ ( p ) ≤ n + τ + log ( n + 1) = O ( n + τ ).Let ε = N as in Lemma 4.2 so that the polynomial p ε = p − ε (cid:80) ki =0 X i is positive over R . By Lemma 4.2,one has N = C ( n + nτ ) for some C >
1. Let us write p ε = (cid:80) ni =0 a i X i with a n = (cid:96) and prove that aprecision of δ := N + log (5 n (cid:107) p (cid:107) ∞ ) = C ( n + nτ ) + log (5 n (cid:107) p (cid:107) ∞ ) is enough to ensure that the coefficientsof u satisfy ε ≥ | u i +1 | − u i + | u i − | , for all i = 0 , . . . , k . First, note that e := 2 − δ < n ( n +1) holds. Byusing Vieta’s formulas provided in Lemma 2.8, one has for all j = 1 , . . . , n : (cid:88) ≤ i < ···
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n = 2 k and τ be an upperbound on the bitsize of the coefficients of f . Then, on input f , Algorithm univsos2 runs in boolean time ∼ O ( n + n τ ) . Proof.
By Lemma 2.3, the square-free decomposition of f can be computed using an expected numberof ∼ O ( n τ ) boolean operations. Checking that the polynomial p ε has no real root can be performed usingan expected number of ∼ O ( n · τ ( ε )) = ∼ O ( n τ ) boolean operations while relying on Sylvester-HabichtSequences [20, Corollary 5.2].As seen in the proof of Theorem 4.3, the complex roots of p ε must be approximated with isolating intervals(resp. disks) of radius less than 2 − τ ( p ε ) − δ . Thus, by Lemma 2.7, all real (resp. complex) roots of p ε canbe computed in ∼ O ( n + n τ ( p ε ) + n ( δ + τ ( p ε )) = ∼ O ( n + n τ ) boolean operations.As in the proof of Theorem 4.3, one can select | u n − j | = | a n − j || − (1 + 2 − δ ) j | , for all j = 1 , . . . , n .This implies that the computation of each coefficient of u can be performed with at most ∼ O ( n · τ ( δ )) = ∼ O ( n + n τ ) boolean operations. Eventually, we obtain a bound of ∼ O ( n + n τ ) for the computation ofall coefficients of u , which yields the desired result.We state now the complexity result for checking the SOS certificates output by Algorithm univsos2 . Asfor the output of Algorithm univsos1 , this is done through evaluation of the output at n + 1 distinctvalues where n is the degree of the output. Theorem 4.5.
Let f ∈ Z [ X ] be a positive polynomial over R , with deg f = n = 2 k and τ be an upperbound on the bitsize of the coefficients of f . Then one can check the correctness of the weighted SOSdecomposition of f obtained with Algorithm univsos2 using ∼ O ( n + n τ ) bit operations.Proof. From [7, Corollary 8.27], the cost of polynomial multiplication in Z [ X ] of degree less than n with coefficients of bitsize upper bounded by l is bounded by ∼ O ( n · l ). By Theorem 4.3, the maximalcoefficient bitsize of the SOS decomposition of f obtained with Algorithm univsos2 is upper boundedby l = O ( n + n τ ). Therefore, from [7, Corollary 10.8], the cost of the evaluation of this decompositionat n points can be performed using at most ∼ O ( n · ( n + nτ )) boolean operations as claimed.16 Practical experiments
Now, we present experimental results obtained by applying Algorithm univsos1 and Algorithm univsos2 ,respectively presented before in Sections 3 and 4. Both algorithms have been implemented in a tool, called univsos , written in Maple version 16. The interested reader can find more details about installation andbenchmark execution on the dedicated webpage. This tool is integrated to the RAGlib Maple package .We obtained all results on an Intel Core i7-5600U CPU (2.60 GHz) with 16Gb of RAM. SOS decomposition(resp. verification) times are provided after averaging over five (resp. thousand) runs.As mentioned in [4, Section 6], the SOS decomposition performed by Algorithm univsos2 has beenimplemented using the PARI/GP software tool and is freely available. To ensure fair comparison, wehave rewritten this algorithm in Maple. To compute approximate complex roots of univariate polynomials,we rely on the PARI/GP procedure polroots through an interface with our Maple library. We also triedto use the Maple procedure fsolve but the polroots routine from Pari/GPyields significantly betterperformance for the polynomials involved in our examples.The nine polynomial benchmarks presented in Table 1 allow to approximate some given mathematicalfunctions, considered in [4, Section 6]. Computation and verification of SOS certificates are a mandatorystep required to validate the supremum norm of the difference between such functions and their respectiveapproximation polynomials on given closed intervals. This boils down to certify two inequalities of theform ∀ x ∈ [ b, c ] , p ( x ) ≥
0, with p ∈ Q [ X ], b, c ∈ Q and deg p = n . As recalled in [4, Section 5.2.5], thislatter problem can be addressed by computing a weighted SOS decomposition of the polynomial q ( Y ) :=(1 + Y ) n p (cid:16) b + cY Y (cid:17) , with either Algorithm univsos1 or Algorithm univsos2 . For each benchmark,we indicate in Table 1 the degree n and the bitsize τ of the input polynomial, the bitsize τ of theweighted SOS decomposition provided by Algorithm univsos1 as well as the corresponding computation(resp. verification) time t (resp. t (cid:48) ). Similarly, we display τ , t , t (cid:48) for Algorithm univsos2 . The tableresults show that for all other eight benchmarks, Algorithm univsos2 yields better certification andverification performance, together with more concise SOS certificates. This observation confirms whatwe could expect after comparing the theoretical complexity results from Sections 3 and 4.Table 1: Comparison results of output size and performance between Algorithm univsos1 and Algo-rithm univsos2 for non-negative polynomial benchmarks from [4].Id n τ (bits) univsos1 univsos2 τ (bits) t (ms) t (cid:48) (ms) τ (bits) t (ms) t (cid:48) (ms) https://github.com/magronv/univsos http://pari.math.u-bordeaux.fr https://hal.archives-ouvertes.fr/ensl-00445343v2 n = 2 k with 10 ≤ n ≤ P n := 1 + X + · · · + X n . The rootsof this polynomial are the n + 1-th roots of unity, thus yielding the following SOS decomposition withreal coefficients: P n := (cid:81) kj =1 (( X − cos θ j ) + sin θ j ), with θ j := jπn +1 , for each j = 1 , . . . , k . By contrastwith the benchmarks from Table 1, Table 2 shows that Algorithm univsos1 produces output certificatesof much smaller size compared to Algorithm univsos2 , with a bitsize ratio lying between 6 and 38 forvalues of n between 10 and 200. This is due to the fact that Algorithm univsos1 outputs a value of tequal to 0 at each step. The execution performance of Algorithm univsos1 are also much better in thiscase. The lack of efficiency of Algorithm univsos2 is due to the computational bottleneck occurring inorder to obtain accurate approximation of the relatively close roots cos θ j ± i sin θ j , j = 1 , . . . , k . For n ≥ univsos2 did not succeed after two hours of computation, as meant bythe symbol − in the corresponding line.Table 2: Comparison results of output size and performance between Algorithm univsos1 and Algo-rithm univsos2 for non-negative power sums of increasing degrees. n univsos1 univsos2 τ (bits) t (ms) t (cid:48) (ms) τ (bits) t (ms) t (cid:48) (ms)10 84 7 0.03 567 264 0.0320 195 10 0.05 1 598 485 0.0640 467 26 0.09 6 034 2 622 0.1860 754 45 0.14 12 326 6 320 0.3280 1 083 105 0.18 21 230 12 153 0.47100 1 411 109 0.26 31 823 19 466 0.69200 3 211 444 0.48 120 831 171 217 2.08300 5 149 1 218 0.74 − − −
400 7 203 2 402 0.95500 9 251 4 292 1.191000 20 483 30 738 2.56Further experiments are summarized in Table 3 for modified Wilkinson polynomials W n of increasingdegrees n = 2 k with 10 ≤ n ≤
600 and W n := 1 + (cid:81) kj =1 ( X − j ) . The complex roots j ± i , j = 1 , . . . , k of W n are relatively close, which yields again a significant lack of performance of Algorithm univsos2 . Asobserved in the case of power sums, timeout behaviors occur for n ≥
60. In addition, the bitsize of theSOS decompositions returned by Algorithm univsos1 are much smaller. This is a consequence of thefact that in this case, a = 1 is the global minimizer of W n . Hence, the algorithm always terminates at thefirst iteration by returning the trivial quadratic approximation f t = f a = 1 together with the square-freedecomposition of W n − f t = (cid:81) kj =1 ( X − j ) .Finally, we consider experimentation performed on modified Mignotte polynomials defined by M n,m := X n +2(101 X − m and N n := ( X n +2(101 X − )( X n +2((101+ ) X − ), for even natural integers n and m ≤
2. The corresponding results are displayed in Table 4 for M n,m with m = 2 and 10 ≤ n ≤ m = n − ≤ n ≤
100 as well as for N n with 10 ≤ n ≤ univsos2 can only handle small size instances,due to limited scalablity of the polroots procedure. In this case a = is the unique global minimizerof M n, . Thus, Algorithm univsos1 always outputs weighted SOS decompositions of polynomials M n, within a single iteration by first computing the quadratic polynomial f t = f a = 2(101 X − and thetrivial square-free decomposition W n − f t = X n . In the absence of such minimizers, Algorithm univsos1 can only handle instances of polynomials M n,n − and N n with moderate degree (less than 100).18able 3: Comparison results of output size and performance between Algorithm univsos1 and Algo-rithm univsos2 for modified Wilkinson polynomials of increasing degrees. n τ (bits) univsos1 univsos2 τ (bits) t (ms) t (cid:48) (ms) τ (bits) t (ms) t (cid:48) (ms)10 140 47 17 0.01 2 373 751 0.0320 737 198 31 0.01 12 652 3 569 0.0840 3 692 939 35 0.01 65 404 47 022 0.1760 9 313 2 344 101 0.01 − − −
80 17 833 4 480 216 0.01100 29 443 7 384 441 0.01200 137 420 34 389 3 249 0.01300 335 245 83 859 11 440 0.01400 628 968 157 303 34 707 0.02500 1 022 771 255 767 73 522 0.02600 1 519 908 380 065 149 700 0.04Table 4: Comparison results of output size and performance between Algorithm univsos1 and Algo-rithm univsos2 for modified Mignotte polynomials of increasing degrees.Id n τ (bits) univsos1 univsos2 τ (bits) t (ms) t (cid:48) (ms) τ (bits) t (ms) t (cid:48) (ms) M n,
10 27 23 2 0.01 4 958 1 659 0.0410 − − − M n,n −
10 288 25 010 21 0.03 6 079 2 347 0.0420 1 364 182 544 138 0.04 26 186 10 922 0.0640 5 936 1 365 585 1 189 0.13 − − −
60 13 746 4 502 551 4 966 0.33100 39 065 20 384 472 38 716 1.66 N n
10 212 25 567 27 0.04 − − −
20 189 336 87 0.0540 5 027 377 1 704 0.1760 16 551 235 8 075 0.84100 147 717 572 155 458 11.1
We presented and analyzed two different algorithms univsos1 and univsos2 to compute weighted sumsof squares (SOS) decompositions of non-negative univariate polynomials. When the input polynomialhas rational coefficients, one feature shared by both algorithms is their ability to provide non-negativitycertificates whose coefficients are also rational. Our study shows that the complexity analysis of Algo-rithm univsos1 yields an upper bound which is exponential w.r.t. the input degree, while the complexityof Algorithm univsos2 is polynomial. However, comparison benchmarks emphasize the need for bothalgorithms to handle various classes of non-negative polynomials, e.g. in the presence of rational globalminimizers or when root isolation can be performed efficiently.19 first direction of further research is a variant of Algorithm univsos2 where one would compute approx-imate SOS decompositions of perturbed positive polynomials by using semidefinite programming (SDP)instead of root isolation. Preliminary experiments yield very promising results when the bitsize of thepolynomials is small, e.g. for power sums of degree up to 1000. However, the performance decrease whenthe bitsize becomes larger, either for polynomial benchmarks from [4] or modified Wilkinson polynomials.At the moment, we are not able to provide any SOS decomposition for all such benchmarks. Our SDP-based algorithm relies on the high-precision solver SDPA-GMP [25] but it is still challenging to obtainprecise values of eigenvalues/vectors of SDP output matrices. Another advantage of this technique isits ability to perform global polynomial optimization. A topic of interest would be to obtain the samefeature with the two current algorithms. We also plan to develop extensions to the non-polynomial case.
References [1] S. Basu, R. Pollack, and M.-F. Roy. On the combinatorial and algebraic complexity of quantifierelimination.
Journal of the ACM (JACM) , 43(6):1002–1045, 1996.[2] S. Basu, R. Pollack, and M.-F. Roy.
Algorithms in Real Algebraic Geometry (Algorithms and Com-putation in Mathematics) . Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006.[3] A. L. B. Cauchy. Calcul des indices des fonctions.
Journal de l’Ecole Polytechnique , 15(25):176 –229, 1832.[4] S. Chevillard, J. Harrison, M. Joldes, and C. Lauter. Efficient and accurate computation of upperbounds of approximation errors.
Theoretical Computer Science , 412(16):1523 – 1543, 2011.[5] G. E Collins. Quantifier elimination for real closed fields by cylindrical algebraic decompostion.In
Automata Theory and Formal Languages 2nd GI Conference Kaiserslautern, May 20–23, 1975 ,pages 134–183. Springer, 1975.[6] The Coq Proof Assistant, 2016. http://coq.inria.fr/ .[7] J. Gathen and J. Gerhard.
Modern Computer Algebra . Cambridge University Press, New York, NY,USA, 1999.[8] A. Girard.
Invention nouvelle en l’alg´ebre . Blauew, 1629.[9] A. Greuet and M. Safey El Din. Probabilistic algorithm for polynomial optimization over a realalgebraic set.
SIAM Journal on Optimization , 24(3):1313–1343, 2014.[10] Q. Guo, M. Safey El Din, and L. Zhi. Computing Rational Solutions of Linear Matrix Inequalities.In
Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation , ISSAC’13, pages 197–204, New York, NY, USA, 2013. ACM.[11] T. Hales, M. Adams, G. Bauer, D. T. Dat, J. Harrison, H. L. Truong, C. Kaliszyk, V. Magron,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 ofthe Kepler Conjecture, 2015. Submitted.[12] J. Harrison. HOL Light: A Tutorial Introduction. In Mandayam K. Srivas and Albert John Camilleri,editors,
FMCAD , volume 1166 of
Lecture Notes in Computer Science , pages 265–269. Springer, 1996.[13] D. Henrion, S. Naldi, and M. Safey El Din. Exact Algorithms for Linear Matrix Inequalities.
SIAMJournal on Optimization , 26(4):2512–2539, 2016.2014] D. Henrion, S/ Naldi, and M. Safey El Din. Spectra-a maple library for solving linear matrixinequalities in exact arithmetic. arXiv preprint arXiv:1611.01947 , 2016.[15] N.J. Higham.
Accuracy and Stability of Numerical Algorithms: Second Edition . SIAM, 2002.[16] H. Hong and M. Safey El Din. Variant quantifier elimination.
Journal of Symbolic Computation ,47(7):883–901, 2012.[17] E. L. Kaltofen, B. Li, Z. Yang, and L. Zhi. Exact certification in global polynomial optimization viasums-of-squares of rational functions with rational coefficients.
Journal of Symbolic Computation ,47(1):1 – 15, 2012.[18] E. Landau. ber die darstellung definiter funktionen durch quadrate.
Mathematische Annalen , 62:272–285, 1906.[19] J.-B. Lasserre. Global optimization with polynomials and the problem of moments.
SIAM Journalon Optimization , 11(3):796–817, 2001.[20] Thomas Lickteig and Marie-Franoise Roy. Sylvesterhabicht sequences and fast cauchy index compu-tation.
Journal of Symbolic Computation , 31(3):315 – 341, 2001.[21] V. Magron, X. Allamigeon, S. Gaubert, and B. Werner. Formal proofs for Nonlinear Optimization.
Journal of Formalized Reasoning , 8(1):1–24, 2015.[22] K. Mehlhorn, M. Sagraloff, and P. Wang. From Approximate Factorization to Root Isolation withApplication to Cylindrical Algebraic Decomposition.
J. Symb. Comput. , 66:34–69, January 2015.[23] S. Melczer and B. Salvy. Symbolic-Numeric Tools for Analytic Combinatorics in Several Variables.In
Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation ,ISSAC ’16, pages 333–340, New York, NY, USA, 2016. ACM.[24] M. Mignotte.
Mathematics for Computer Algebra . Springer-Verlag New York, Inc., New York, NY,USA, 1992.[25] M. Nakata. A numerical evaluation of highly accurate multiple-precision arithmetic version ofsemidefinite programming solver: SDPA-GMP, -QD and -DD. In
Computer-Aided Control SystemDesign (CACSD), 2010 IEEE International Symposium on , pages 29–34, Sept 2010.[26] P. Parrilo.
Structured semidefinite programs and semialgebraic geometry methods in robustness andoptimization . PhD thesis, California Institute of Technology, 2000.[27] H. Peyrl and P. Parrilo. Computing sum of squares decompositions with rational coefficients.
Theor.Comput. Sci. , 409(2):269–281, 2008.[28] Y. Pourchet. Sur la reprsentation en somme de carrs des polynmes une indtermine sur un corps denombres algbriques.
Acta Arithmetica , 19(1):89–104, 1971.[29] A. Prestel and C. Delzell.
Positive Polynomials: From Hilberts 17th Problem to Real Algebra . SpringerMonographs in Mathematics. Springer Berlin Heidelberg, 2001.[30] A. Rantzer and P. A. Parrilo. On convexity in stabilization of nonlinear systems. In
Proceedings ofthe 39th IEEE Conference on Decision and Control , volume 3, pages 2942–2945 vol.3, 2000.[31] M. Safey El Din and L. Zhi. Computing rational points in convex semialgebraic sets and sum ofsquares decompositions.
SIAM Journal on Optimization , 20(6):2876–2889, 2010.2132] M. Schweighofer. Algorithmische Beweise f¨ur Nichtnegativ- und Positivstellens¨atze. Master’s thesis,Diplomarbeit an der Universit¨at Passau, 1999.[33] A. Strzebonski and E. Tsigaridas. Univariate Real Root Isolation in an Extension Field. In
Pro-ceedings of the 36th International Symposium on Symbolic and Algebraic Computation , ISSAC ’11,pages 321–328, New York, NY, USA, 2011. ACM.[34] W. Swokowski.
Fundamentals of College Algebra . PWS-Kent Pub. Co., 1989, pages 216 - 221.[35] D. Y.Y. Yun. On Square-free Decomposition Algorithms. In