Algorithms for Weighted Sums of Squares Decomposition of Non-negative Univariate Polynomials

TL;DR

This paper compares two algorithms for computing weighted sums of squares decompositions of non-negative univariate polynomials, analyzing their theoretical complexity, practical performance, and implementation results.

Contribution

It provides the first bit complexity analysis of both algorithms, highlighting their differences and practical efficiencies for polynomial decomposition.

Findings

The second algorithm generally outperforms the first in efficiency.

The first algorithm is better for certain families of polynomials.

Complexity estimates are exponential for the first and polynomial for the second algorithm.

Abstract

It is well-known that every non-negative univariate real polynomial can be written as the sum of 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 to lie in the field generated by the coefficients of the polynomial. In this article, we describe, analyze and compare both from the theoretical and practical points of view, two algorithms computing such a weighted sums of squares decomposition for univariate polynomials with rational coefficients. The first algorithm, due to the third author relies on real root isolation, quadratic approximations of positive polynomials and square-free decomposition but its complexity was not analyzed. We provide bit complexity estimates, both on runtime and output size of this algorithm. They are exponential in the degree of the input univariate polynomial and linear in the maximum bitsize of 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 root isolation and square-free decomposition and has been introduced for certifying positiveness of polynomials in the context of computer arithmetics. Again, its complexity was not analyzed. We provide bit complexity estimates, both on runtime and output size of this algorithm, which are polynomial 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. While the second algorithm is, as expected from the complexity result, more efficient on most of examples, we exhibit families of non-negative polynomials for which the first algorithm is better.

Tables

Table 1. Table 1: Comparison results of output size and performance between Algorithm  univsos1 and Algorithm  univsos2 for non-negative polynomial benchmarks from  [ 4 ] .

Id	$n$	$\tau$ (bits)	univsos1			univsos2
			${\tau }_1$ (bits)	$t_1$ (ms)	$t_1^{\prime }$ (ms)	${\tau }_2$ (bits)	$t_2$ (ms)	$t_2^{\prime }$ (ms)
# 1	13	22 682	3 403 218	2 723	0.40	51 992	824	0.14
# 3	32	269 958	11 613 480	13 109	1.18	580 335	2 640	0.68
# 4	22	47 019	1 009 507	4 063	1.45	106 797	1 776	0.31
# 5	34	117 307	8 205 372	102 207	20.1	265 330	5 204	0.60
# 6	17	26 438	525 858	1 513	0.74	59 926	1 029	0.21
# 7	43	67 399	62 680 827	217 424	48.1	152 277	11 190	0.87
# 8	22	27 581	546 056	1 979	0.77	63 630	1 860	0.38
# 9	20	30 414	992 076	964	0.44	68 664	1 605	0.25
# 10	25	42 749	3 146 982	1 100	0.38	98 926	2 753	0.39

Table 2. Table 2: Comparison results of output size and performance between Algorithm  univsos1 and Algorithm  univsos2 for non-negative power sums of increasing degrees.

$n$	univsos1			univsos2
	${\tau }_1$ (bits)	$t_1$ (ms)	$t_1^{\prime }$ (ms)	${\tau }_2$ (bits)	$t_2$ (ms)	$t_2^{\prime }$ (ms)
10	84	7	0.03	567	264	0.03
20	195	10	0.05	1 598	485	0.06
40	467	26	0.09	6 034	2 622	0.18
60	754	45	0.14	12 326	6 320	0.32
80	1 083	105	0.18	21 230	12 153	0.47
100	1 411	109	0.26	31 823	19 466	0.69
200	3 211	444	0.48	120 831	171 217	2.08
300	5 149	1 218	0.74	$-$	$-$	$-$
400	7 203	2 402	0.95
500	9 251	4 292	1.19
1000	20 483	30 738	2.56

Table 3. Table 3: Comparison results of output size and performance between Algorithm  univsos1 and Algorithm  univsos2 for modified Wilkinson polynomials of increasing degrees.

$n$	$\tau$ (bits)	univsos1			univsos2
		${\tau }_1$ (bits)	$t_1$ (ms)	$t_1^{\prime }$ (ms)	${\tau }_2$ (bits)	$t_2$ (ms)	$t_2^{\prime }$ (ms)
10	140	47	17	0.01	2 373	751	0.03
20	737	198	31	0.01	12 652	3 569	0.08
40	3 692	939	35	0.01	65 404	47 022	0.17
60	9 313	2 344	101	0.01	$-$	$-$	$-$
80	17 833	4 480	216	0.01
100	29 443	7 384	441	0.01
200	137 420	34 389	3 249	0.01
300	335 245	83 859	11 440	0.01
400	628 968	157 303	34 707	0.02
500	1 022 771	255 767	73 522	0.02
600	1 519 908	380 065	149 700	0.04

Table 4. Table 4: Comparison results of output size and performance between Algorithm  univsos1 and Algorithm  univsos2 for modified Mignotte polynomials of increasing degrees.

Id	$n$	$\tau$ (bits)	univsos1			univsos2
			${\tau }_1$ (bits)	$t_1$ (ms)	$t_1^{\prime }$ (ms)	${\tau }_2$ (bits)	$t_2$ (ms)	$t_2^{\prime }$ (ms)
$M_{n,2}$	10	27	23	2	0.01	4 958	1 659	0.04
	${10}^2$			3		$-$	$-$	$-$
	${10}^3$			85
	${10}^4$			3 041
$M_{n,n-2}$	10	288	25 010	21	0.03	6 079	2 347	0.04
	20	1 364	182 544	138	0.04	26 186	10 922	0.06
	40	5 936	1 365 585	1 189	0.13	$-$	$-$	$-$
	60	13 746	4 502 551	4 966	0.33
	100	39 065	20 384 472	38 716	1.66
$N_n$	10	212	25 567	27	0.04	$-$	$-$	$-$
	20		189 336	87	0.05
	40		5 027 377	1 704	0.17
	60		16 551 235	8 075	0.84
	100		147 717 572	155 458	11.1

Equations

$$f=\prod _j(X-a_j{)}^{r_j}\prod _k((X-(b_k+i{c}_k))(X-(b_k-i{c}_k)){)}^{s_k},$$

$$g^2=\prod _j(X-a_j{)}^{r_j},q+ir=\prod _k(X-(b_k+i{c}_k){)}^{s_k}q-ir=\prod _k(X-(b_k-i{c}_k){)}^{s_k}.$$

$$|x|\le \max \left\{1,\frac{|a_0|}{|a_n|}+\cdots +\frac{|a_{n-1}|}{|a_n|}\right\}.$$

$$\sum _{1\le i_1<\cdots <i_j\le n}{z}_{i_1}\dots z_{i_j}=(-1{)}^j\frac{a_{n-j}}{a_n}.$$

$$f_t:=f(t)+f^{\prime }(t)(X-t)+\frac{f^{\prime }(t{)}^2}{4f(t)}(X-t{)}^2\in K[X].$$

$$f_t:=f(t)+f^{\prime }(t)(X-t)+\frac{f^{\prime }(t{)}^2}{4f(t)}(X-t{)}^2\in \mathbb{R}[X].$$

$$f=\sum _{k=0}^n\frac{f^{(k)}(t)}{k!}(X-t{)}^k,$$

$$f-f_t=\sum _{k=2}^n\frac{f^{(k)}(t)}{k!}(X-t{)}^k-\frac{f^{\prime }(t{)}^2}{4f(t)}(X-t{)}^2=(X-t{)}^2\left(\sum _{k=2}^n\frac{f^{(k)}(t)}{k!}(X-t{)}^{k-2}-\frac{f^{\prime }(t{)}^2}{4f(t)}\right)$$

$$h:=f(T)\left(\sum _{k=2}^n\frac{f^{(k)}(T)}{k!}(X-T{)}^{k-2}\right)-\frac{1}{4}{f}^{\prime }(T{)}^2\in \mathbb{R}[T,X].$$

$$H=c\sum _{k=2}^n\left(\genfrac{}{}{0ex}{}{\alpha }{k}\right)T^{\alpha -k}(X-T{)}^{k-2},$$

$$(X-T{)}^{2}H=c((T+(X-T){)}^{\alpha }-T^n-n{T}^{\alpha -1}(X-T))=c(X^{\alpha }-\alpha T^{\alpha -1}X+(\alpha -1)T^{\alpha }).$$

$$f_t:=f(t)+f^{\prime }(t)(X-T)+\frac{f^{\prime }(t{)}^2}{4f(t)}(X-T{)}^2\in \mathbb{R}[X],$$

$$U\to \mathbb{R}:t\mapsto C_t:=\max \left\{1,\frac{|a_{0t}|}{|a_{nt}|},\dots ,\frac{|a_{(n-1)t}|}{|a_{nt}|}\right\}.$$

$$f:=\Biggl{(}(X+1)^{2}\biggl{(}(X-1)^{2}\Bigl{(}\frac{1}{16}(X-\frac{19108973}{17085096})^{2}\Bigr{)}+\frac{502920}{237293}(-\frac{1}{18}X+\frac{88411}{167640})^{2}\biggr{)}\Biggr{)}+\frac{720}{1397}(-\frac{19}{16}X+\frac{271}{360})^{2}.$$

$$S:=\{t\in \mathbb{Q}\mid \forall x\in \mathbb{R},f(t{)}^2+f^{\prime }(t)f(t)(x-t)+f^{\prime }(t{)}^2(x-t{)}^2\le 4f(t)f(x)\}.$$

$$M_2(t)$$

$$M_1(t)$$

$$M_0(t)$$

$$t_{i2}^{4i}{f}_i(t_i)f_i-q_i=(X-t_i{)}^2{f}_{i-1},i=2,\dots ,k.$$

$$\stackrel{\sim }{\mathcal{O}}(k^3\cdot (k!{)}^3\tau )=\stackrel{\sim }{\mathcal{O}}\left({\left(\frac{n}{2}\right)}^{\frac{3n}{2}}\tau \right).$$

$$\stackrel{\sim }{\mathcal{O}}(k\cdot (k!{)}^3\tau )=\stackrel{\sim }{\mathcal{O}}((\frac{n}{2}{)}^{\frac{3n}{2}}\tau )$$

$$\text{c\_list}=$$

$$\text{s\_list}=$$

$$\displaystyle f=\frac{1}{32}\biggl{(}X^{3}-\frac{69}{8}X\biggr{)}^{2}+\frac{1}{32}\biggl{(}7X^{2}-\frac{1}{4}X-\frac{63}{8}\biggr{)}^{2}+\frac{913}{15360}+\frac{731}{92160}X^{2}$$

$$\displaystyle+\frac{7}{1152}X^{4}+\frac{1}{32}X^{6}+\frac{79}{7680}\biggl{(}X+\frac{1}{2}\biggr{)}^{2}+\frac{1}{576}X^{2}\biggl{(}X-\frac{1}{2}\biggr{)}^{2}\,.$$

$$\underset{x\in \mathbb{R}}{\inf }p(x)>(n{2}^{\tau }{)}^{-n+2}{2}^{-n{\log }_{2}n-n\tau }.$$

$$\underset{x\in \mathbb{R}}{\inf }p(x)>(n{2}^{\tau }{)}^{-n+2}{2}^{-n{\tau }^{\prime }}$$

$$2^N\le n{R}^n(n{2}^{\tau }{)}^{n-2}{2}^{n{\log }_{2}n+n\tau }=n^{n+1}{2}^n{2}^{n\tau }(n{2}^{\tau }{)}^{n-2}{2}^{n{\log }_{2}n+n\tau }.$$

$$\sum _{1\le i_1<\cdots <i_j\le n}{z}_{i_1}\dots z_{i_j}=(-1{)}^j\frac{a_{n-j}}{\ell }.$$

$$u_{n-j}={\ell}\sum_{1\leq i_{1}<\dots<i_{j}\leq n}(z_{i_{1}}\dots z_{i_{j}}-\hat{z}_{i_{1}}\dots\hat{z}_{i_{j}})=\sum_{1\leq i_{1}<\dots<i_{j}\leq n}z_{i_{1}}\dots z_{i_{j}}\bigl{(}1-(1+e_{i_{1}})\dots(1+e_{i_{j}})\bigr{)}\,.$$

$$\frac{|u_{2i+1}|}{4}-u_{2i}+|u_{2i-1}|\leq{e\|p_{\varepsilon}\|_{\infty}}\bigl{(}\frac{2n}{4}+2n-1+2n-2\bigr{)}\leq 5ne\|p_{\varepsilon}\|_{\infty}\leq 5ne\|p\|_{\infty}\,.$$

$$\stackrel{\sim }{\mathcal{O}}(n^4+n^3\tau ).$$

Full text

In this article, we describe, analyze and compare both from the theoretical and practical points of view, two algorithms computing such a weighted sums of squares decomposition for univariate polynomials with rational coefficients.

The first algorithm, due to the third author relies on real root isolation, quadratic approximations of positive polynomials and square-free decomposition but its complexity was not analyzed. We provide bit complexity estimates, both on the runtime and the output size of this algorithm. They are exponential in the degree of the input univariate polynomial and linear in the maximum bitsize of 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 root isolation and square-free decomposition and has been introduced for certifying positiveness of polynomials in the context of computer arithmetics. Again, its complexity was not analyzed. We provide bit complexity estimates, both on the runtime and the output size of this algorithm, which are polynomial 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 on several application benchmarks. While the second algorithm is, as expected from the complexity result, more efficient on most of examples, we exhibit families of non-negative polynomials for which the first algorithm is better.

Keywords:

non-negative univariate polynomials, Nichtnegativstellensätze, sums of squares decomposition, root isolation, real algebraic geometry.

1 Introduction

Given a subfield $K$ of $\mathbb{R}$ and a non-negative univariate polynomial $f\in K[X]$, we consider the problem of proving the existence and computing the weighted sums of squares decompositions of $f$ with coefficients also lying in $K$.

Beyond the theoretical interest of this question, finding certificates of non-negative polynomials is mandatory in many application fields. Among them, one can mention the stability proofs of critical control systems often relying on Lyapunov functions [30], the certified evaluation of mathematical functions in the 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 is already an important one. In particular, formal proofs of polynomial non-negativity can be handled with sums of squares certificates. These certificates are obtained with tools available outside of the proof assistants and eventually verified inside. Because of the limited computing power available inside such proof assistants, this is crucial to devise algorithms that produce certificates, whose checking is computationally reasonably simple. In particular, we would like to ensure that such algorithms output sums of squares certificates of moderate bitsize and ultimately with a computational complexity being polynomial with respect to the input.

Related Works.

Decompositing non-negative univariate polynomials into sums of squares has a long story ; very early quantitative aspects like the number of needed squares have been studied. For the case $K=\mathbb{Q}$, Landau shows in [18] that for every non-negative polynomial in $\mathbb{Q}[X]$, there exists a decomposition involving a weighted sum of (at most) eight polynomial squares in $\mathbb{Q}[X]$. In [28], Pourchet improves this result by showing the existence of a decomposition involving only a weighted sum of (at most) five squares. This is done using approximation and valuation theory ; extracting an algorithm from these 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 semidefinite positive 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 semidefinite programming (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 such SDP solvers are most of the time implemented with floating-point double precision. More accurate solvers are available (e.g. SDPA-GMP [25]). However, note that these solvers always compute numerical approximations of the algebraic solution to the semidefinite program under consideration. Hence, they are not sufficient to provide algebraic certificates of posivity with rational coefficients. Hence, a process for making exact and with rational coefficients the computed numerical approximations of sums of squares certificates 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 their efficiency, these methods do not provide any guarantee to output a rational solution to a Linear Matrix Inequality 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 constraints using the Cylindrical Algebraic Decomposition algorithm [5] or more efficient critical point methods (see e.g. [1] for complexity estimates and [16, 9] for practical algorithms). But using such general algorithms is an overkill and dedicated algorithms have been designed for computing exactly algebraic solutions to Linear Matrix Inequalities [13, 14]. Computing rational solutions can also be considered thanks to convexity properties [31]. In particular, the algorithm in [10] can be used to compute sums of squares certificates with rational coefficients fpr a non-negative univariate polynomial of degree $n$ with coefficients of bit size bounded by $\tau$ using $\tau^{\operatorname{\mathcal{O}}{(1)}}2^{\operatorname{\mathcal{O}}{(n^{3})}}$ bit operations at most (see [10, Theorem 1]).

For the case where $K$ is an arbitrary subfield of $\mathbb{R}$, Schweighofer gives in [32] a new proof of the existence of a decomposition involving a sum of (at most) $n$ polynomial squares in $K[X]$. This existence proof comes together with a recursive algorithm to compute such decompositions. At each recursive step, the algorithm performs real root isolation and quadratic approximations of positive polynomials. Later on, a second algorithm is derived in [4, Section 5.2], where the authors show the existence of a decomposition involving a sum of (at most) $n+3$ polynomial squares in $K[X]$. This algorithm is based on approximating complex 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 better complexities than the algorithm in [10], the second algorithm being polynomial in $n$ and $\tau$.

Notation for complexity estimates.

For complexity estimates, we use the bit complexity model. For an integer $b\in\mathbb{Z}\backslash\{0\}$, we denote by $\tau(b):=\log_{2}(|b|)+1$ the bitsize of $b$, with the convention $\tau(0):=1$. We write a given polynomial $f\in\mathbb{Z}[X]$ of degree $n\in\mathbb{N}$ as $f=\sum_{i=0}^{n}b_{i}X^{i}$, with $b_{0},\dots,b_{n}\in\mathbb{Z}$. In this case, we define $\|f\|_{\infty}:=\max_{0\leq i\leq n}|b_{i}|$ and, using a slight abuse of notation, we denote $\tau(\|f\|_{\infty})$ by $\tau(f)$. Observe that when $f$ has degree $n$, the bit size necessary to encode $f$ is bounded by $n\tau(f)$ (when storing each coefficients of $f$). The derivative of $f$ is $f^{\prime}=\sum_{i=1}^{n}ib_{i}X^{i-1}$. For a rational number $q=\frac{b}{c}$, with $b\in\mathbb{Z},c\in\mathbb{Z}\backslash\{0\}$ and ${\rm gcd}(b,c)=1$, we denote $\max\{\tau(b),\tau(c)\}$ by $\tau(q)$. For two mappings $g,h:\mathbb{N}^{l}\to\mathbb{R}$, the expression “$g(v)=\operatorname{\mathcal{O}}{(h(v))}$” means that there exists a integer $b\in\mathbb{N}$ such that for all $v\in\mathbb{N}^{l}$, $g(v)\leq bh(v)$. The expression “$g(v)=\operatorname{\overset{\sim}{\mathcal{O}}}{(h(v))}$” means that there exists a integer $c\in\mathbb{N}$ such that for all $v\in\mathbb{N}^{l}$, $g(v)\leq h(v)\log_{2}(h(v))^{c}$.

Contributions.

We present and analyze two algorithms, denoted by univsos1 and univsos2, allowing to decompose a non-negative univariate polynomial $f$ of degree $n$ into sums of squares with coefficients lying in any subfield $K$ of $\mathbb{R}$. To the best of our knowledge, there was no prior complexity estimate for the output of such certification algorithms based on sums of squares in the univariate case. We summarize our contributions as follows:

• •

We describe in Section 3 the first algorithm, called univsos1 in the sequel. It was originally given in [32, Chapter 2] ; Section 3 can be seen as a partial English translation of this text written in German since some proofs have been significantly simplified. In the same section, we analyze its bit complexity. When the input is a polynomial of degree $n$ with integer coefficients of maximum bitsize $\tau$, we prove that Algorithm univsos1 uses $\operatorname{\overset{\sim}{\mathcal{O}}}{((\frac{n}{2})^{\frac{3n}{2}}\tau)}$ boolean operations and returns polynomials of bitsize bounded by $\operatorname{\mathcal{O}}{((\frac{n}{2})^{\frac{3n}{2}}\tau)}$. This is not restrictive: when $f\in\mathbb{Q}[X]$, one can multiply it by the least common multiple of the denominators of its coefficients and apply our statement for polynomials in $\mathbb{Z}[X]$.

• •

We describe in Section 4 the second algorithm univsos2, initially given in [4, Section 5.2]. We also analyze its bit complexity. When the input is a univariate polynomial of degree $n$ with integer coefficients of maximum bitsize $\tau$, we prove that Algorithm univsos2 returns a sums of square decompositions of $n+3$ polynomials with coefficients of bitsize bounded by $\operatorname{\mathcal{O}}{(n^{3}+n^{2}\tau)}$ using $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{4}+n^{3}\tau)}$ boolean operations.

• •

Both algorithms are implemented within the univsos tool. The first release of univsos is a Maple library, is freely available111https://github.com/magronv/univsos and is integrated in the RAGlib (Real Algebraic Library) Maple package222http://www-polsys.lip6.fr/~safey/RAGLib/. The scalability of the library is evaluated in Section 5 on several non-negative polynomials issued from the existing literature. Despite the significant difference of theoretical complexity between the two algorithms, numerical benchmarks indicate that both may yield competitive performance w.r.t. specific sub-classes of problems.

2 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\in\mathbb{R}[X]$ be a non-negative univariate polynomial, i.e. $f(x)\geq 0$, for all $x\in\mathbb{R}$. Then, $f$ can be written as the sum of two polynomial squares in $\mathbb{R}[X]$.

Proof.

Without loss of generality, one can assume that $f$ is monic, i.e. the leading coefficient (nonzero coefficient of highest degree) of $f$ is $1$. Then we decompose $f$ as follows in $\mathbb{C}[x]$:

[TABLE]

with $a_{j},b_{k},c_{k}\in\mathbb{R}$, $r_{j},s_{k}\in\mathbb{N}^{>0}$, $a_{j}$ standing for the real roots of $f$ and $(b_{k}\pm ic_{k})$ standing for the complex conjugate roots of $f$. Since $f$ is non-negative, all real roots must have even multiplicity $r_{j}$, yielding the existence of polynomials $g,q,r\in\mathbb{R}[X]$ satisfying the following:

[TABLE]

Then, one has $f=g^{2}(q+ir)(q-ir)=g^{2}(q^{2}+r^{2})=(gq)^{2}+(gr)^{2}$, which proves the claim. ∎

Let $K$ be a field and $g\in K[X]$. One says that $g$ is a square-free polynomial, when there is no prime element $p\in K[X]$ such that $p^{2}$ divides $g$. Let now $f\in K[X]-\{0\}$. A decomposition of $f$ of the form $f=ag_{1}^{1}g_{2}^{2}\dots g_{n}^{n}$ with $a\in K$ and normalized pairwise coprime square-free polynomials $g_{1},g_{2},\dots,g_{n}$ is called 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\in\mathbb{Z}[X]$ and $q\in\mathbb{Z}[X]$ divides $p$ in $\mathbb{Z}[X]$, then one has $\tau(q)\leq\deg q+\tau(p)+\log_{2}(\deg p+1)$.*

The algorithm of Yun [35] (also described in [7, Algorithm 14.21]) allows to compute a square-free decomposition of polynomials with coefficients in a field of characteristic 0.

Lemma 2.3**.**

[7*, § 11.2]**

Let $f\in\mathbb{Z}[X]$ with degree at most $n$ with coefficient bitsize upper bounded by $\tau$. Then the square-free decomposition of $f$ using Yun’s Algorithm [35] can be computed using an expected number of $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{2}\tau)}$ boolean operations.*

Lemma 2.4**.**

[24*, § 6.3.1]**

Let $K$ be a field of characteristic 0 and $L$ a field extension of $K$. The square-free decomposition in $L[X]$ of any polynomial $f\in K[X]-\{0\}$ is the same as the square-free decomposition of $f$ in $K[X]$. Any polynomial $f\in K[X]-\{0\}$ 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 univariate polynomial.

Lemma 2.5**.**

**(Cauchy Bound [3])

Let $K$ be an ordered field. Let $a_{0},\dots,a_{n}\in K$ with $a_{n}\neq 0$. Let $x\in K$ such that $\sum_{i=0}^{n}a_{i}x^{i}=0$. Then, one has:**

[TABLE]

For polynomial with integer coefficients, one has the following:

Lemma 2.6**.**

[24*]**

Let $f\in\mathbb{Z}[X]$ of degree $n$, with coefficient bitsize upper bounded by $\tau$. If $f(x)=0$ and $x\neq 0$, then $\frac{1}{2^{\tau}+1}\leq|x|\leq 2^{\tau}+1$.*

The real (resp. complex) roots of a polynomial can be approximated using root isolation techniques. To compute the real roots one can use algorithms based on Uspensky’s method relying on Descartes’s rule of sign, see e.g. [2, Chap. 10] for a general description of real root isolation algorithms.

Lemma 2.7**.**

[22*, Theorem 5]**

Let $f\in\mathbb{Z}[X]$ with degree at most $n$ with coefficient bitsize upper bounded by $\tau$. Isolating intervals (resp. disks) of radius less than $2^{-\kappa}$ for all real (resp. complex) roots of $f$ can be computed in $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{3}+n^{2}\tau+n\kappa)}$ boolean operations.*

Vieta’s formulas provide relations between the coefficients of a polynomial and signed sums and products of the complex roots of this polynomial:

Lemma 2.8**.**

**(Vieta’s formulas [8])

Let $K$ be an ordered field. Given a polynomial $f=\sum_{i=0^{n}}a_{i}X^{i}\in K[X]$ with $a_{n}\neq 0$ with (not necessarily distinct) complex roots $z_{1},\dots,z_{n}$, one has for all $j=1,\dots,n$:**

[TABLE]

3 Nichtnegativstellensätze with quadratic approximations

3.1 A proof of the existence of SOS decompositions

Lemma 3.1**.**

Let $K$ be an ordered field. Let $g=aX^{2}+bX+c\in K[X]$ with $a,b,c\in K$ and $a\neq 0$. Then, $g$ can be rewritten as $g=a\left(X+\frac{b}{2a}\right)^{2}+\left(c-\frac{b^{2}}{4a}\right)$. Moreover, when $g$ is non-negative over $K$, one has $a>0$ and $c-\frac{b^{2}}{4a}\geq 0$.

Proof.

The decomposition of $g$ is straightforward. Assume that $g$ is non-negative over $K$. Remark that $c-\frac{b^{2}}{4a}=g\left(-\frac{b}{2a}\right)$ ; hence since we assume that $g$ is non-negative over $K$ we deduce that $c-\frac{b^{2}}{4a}\geq 0$.

It remains to prove that $a>0$ which we do by contradiction, assuming that $a<0$. Then, this implies that for all $x\in K$, one has $\left(x+\frac{b}{2a}\right)^{2}\leq-\frac{1}{a}\left(c-\frac{b^{2}}{4a}\right)$. Thus, there exists $C\in K$ such that $x^{2}\leq C$, for each $x\in K$. This implies in particular for $x=2$ that $4\leq C$ and for $x=C$ that $C^{2}\leq C$, thus $C\leq 1$. Finally, one obtains $4\leq C\leq 1$, yielding a contradiction. ∎

Let $f\in K[X]$ be a square-free polynomial which is non-negative over $\mathbb{R}$. Then, $f$ is positive over $\mathbb{R}$, otherwise $f$ would have at least one real root, implying that $f$ would be neither a square-free polynomial in $\mathbb{R}[X]$ nor a square-free polynomial in $K[X]$, according to Lemma 2.4. We want to find a polynomial $g\in K[X]$ which fulfills the following conditions:

1. (i)

$\deg g\leq 2$, 2. (ii)

$g$ is non-negative over $\mathbb{R}$, 3. (iii)

$f-g$ is non-negative over $\mathbb{R}$, 4. (iv)

$f-g$ has a root $t\in K$.

Assume that Property (i) holds. Then the existence of a sum of squares decomposition for $g$ is ensured from Property (ii). Property (iii) implies that $h=f-g$ has only non-negative values over $\mathbb{R}$. The aim of Property (iv) is to ensure the existence of a root $t\in K$ of $h$, which is stronger than the existence of a real root. Note that the case where the degree of $h=f-g$ is less than the degree of $f$ occurs only when $\deg f=2$. In this latter case, we can rely on Lemma 3.1.

Now, we investigate the properties of a polynomial $g\in K[X]$, which fulfills conditions (i)-(iii) and (iv) with $t\in K$. Using Property (i) and Taylor Decomposition, we obtain $g=g(t)+g^{\prime}(t)(X-t)+c(X-t)^{2}$ for some $c\in K$. By Property (iv), one has $g(t)=f(t)$. In addition, Property (iii) yields $f(x)-g(x)\geq 0=f(t)-g(t)$, for all $x\in K$, which implies that $(f-g)^{\prime}(t)=0$ and $g^{\prime}(t)=f^{\prime}(t)$. By Property (ii), the quadratic polynomial $g(X+t)=f(t)+f^{\prime}(t)X+cX^{2}$ has at most one root. This implies that the discriminant of $g$, namely $f^{\prime}(t)^{2}-4cf(t)$ cannot be positive, thus one has $c\geq\frac{f^{\prime}(t)^{2}}{4f(t)}$ (since $f(t)>0$).

Finally, given a polynomial $g$ satisfying (i)-(iii) and (iv) for some $t\in K$, one necessarily has $g=f_{t,c}$ with $\frac{f^{\prime}(t)^{2}}{4f(t)}\leq c\in K$ and $f_{t,c}=f(t)+f^{\prime}(t)(X-t)+c(X-t)^{2}$.

In this case, one also has that the polynomial $g=f_{t,c^{\prime}}$, with $c^{\prime}=\frac{f^{\prime}(t)^{2}}{4f(t)}$, fulfills (i)-(iii) and (iv). Indeed, (i) and (iv) trivially hold. Let us prove that (ii) holds: when $\deg f_{t,c^{\prime}}=0$, then $g=f(t)\geq 0$ and when $\deg f_{t,c^{\prime}}=2$, then $g$ has a single root $\frac{-f^{\prime}(t)}{2c}$ and the minimum of $g$ is $g\left(\frac{-f^{\prime}(t)}{2c}\right)=0$. The inequalities $f_{t,c^{\prime}}\leq f_{t,c}\leq f$ over $\mathbb{R}$ yield (iii).

Therefore, given $f\in K[X]$ with $f$ positive over $\mathbb{R}$, we are looking for $t\in K$ such that the inequality $f\geq f_{t}$ holds over $\mathbb{R}$, with

[TABLE]

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)\mid x\in\mathbb{R}\}$. The idea is then to find $t$ in the neighborhood of a global minimizer of $f$. The following lemma shows that the inequality $f_{t}\leq f$ can be always satisfied for $t$ in some neighborhood of a local minimizer of $f$.

Lemma 3.2**.**

Let $f\in\mathbb{R}[X]$. Let $a$ be a local minimizer of $f$ and suppose that $f(a)>0$. For all $t\in\mathbb{R}$ with $f(t)\neq 0$, let us define the polynomial $f_{t}$:

[TABLE]

Then, there exists a neighborhood $U\subset\mathbb{R}$ of $a$ such that the inequality $f_{t}(x)\leq f(x)$ holds for all $(t,x)\in U\times U$.

Proof.

Set $n:=\deg f$. It is easy to see that we can suppose without loss of generality that $a$ is the origin and that $f(0)=1$. Because of the Taylor formula

[TABLE]

we have

[TABLE]

for all $t\in\mathbb{R}$ with $f(t)\neq 0$. Let $h$ be the bivariate polynomial defined as follows:

[TABLE]

Let us prove that $(a,a)$ is a local minimizer of $h$.

Since $f(0)=1$, there exists $c\neq 0$, $\alpha\in\mathbb{N}$ and $g\in\mathbb{R}[X]$ such that $f-1=cX^{\alpha}+X^{\alpha+1}g$. Therefore, $\lim_{x\to 0}\frac{f(x)-1}{cx^{\alpha}}=1$. Since $f-1$ is non-negative over $\mathbb{R}$, one concludes that $c>0$ and $\alpha$ is even. Let us consider the lowest homogeneous part $H$ of $h$, that is the sum of all monomials of lowest degree involved in $h$. The lowest homogeneous part of $f^{\prime}(T)^{2}$ is $c^{2}(\alpha-1)^{2}T^{2\alpha-2}$ with degree $2\alpha-2$ while the lowest homogeneous part of $\sum_{k=2}^{n}\frac{f^{(k)}(T)}{k!}(X-T)^{k-2}$ is $c\sum_{k=2}^{n}\binom{\alpha}{k}T^{\alpha-k}(X-T)^{k-2}$ with degree $\alpha-2$. Then

[TABLE]

and thus

[TABLE]

Since $\lim_{\|(x,t)\|\to 0}\frac{h(x,t)}{H(x,t)}=1$, it is enough to prove that $H$ is positive except at the origin in order to show that $(a,a)=(0,0)$ is a local minimizer of $h$. Let us consider $(t,x)\in\mathbb{R}^{2}\setminus\{0\}$ and show that $H(t,x)>0$. If $t=x$, we have $H(t,x)=H(x,x)=\binom{\alpha}{2}x^{\alpha-2}>0$. If $t\neq x$, then it is enough to show that $(x-t)^{2}H(t,x)=c(x^{\alpha}-\alpha t^{\alpha-1}x+(\alpha-1)t^{\alpha})>0$. This is clear if $t=0$ since $c>0$ and $\alpha$ is even. Now suppose that $t\neq 0$ and define $\xi:=\frac{x}{t}\neq 1$. Then one has $t^{-\alpha}(x-t)^{2}H(t,x)=c(\xi^{\alpha}-\alpha\xi+\alpha-1)>0$. The positivity of $H$ follows from the fact that the univariate polynomial $r:=X^{\alpha}-\alpha X+\alpha-1$ is positive except at $1$ since $r^{\prime}=\alpha X^{\alpha-1}-\alpha$. The positivity of $H$ implies that $(a,a)$ is a local minimizer of $h$.

Let us define $q(X,T):=(X-T)^{2}h$. Combining the fact that $(a,a)$ is a local minimizer of the two polynomials $h$, $(X-T)^{2}$ and the fact that $h(a,a)=f(a)f^{\prime\prime}(a)-\frac{1}{4}f^{\prime}(a)^{2}=0$, we conclude that $(a,a)$ is also a local minimizer of $q$. Since $f(x)-f_{t}(x)=f(t)\,q(x,t)$, this yields the existence of a neighborhood $O\subset\mathbb{R}^{2}$ of $(a,a)$ such that the inequality $f-f_{t}\geq 0$ holds for all $(x,t)\in O$. Since there exists some neighborhood $U\subset\mathbb{R}$ of $a$, such that the rectangle $U\times 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)\leq f(x)$ holds for all $(x,t)\in U\times U$ . Now, we show that with such a neighborhood $U$ of the smallest global minimizer of $f$, the inequality $f_{t}(x)\leq f(x)$ holds for all $t\in U$ and for all $x\in\mathbb{R}$.

Proposition 3.3**.**

Let $f\in\mathbb{R}[X]$ with $\deg f>0$. Assume that $f$ is positive over $\mathbb{R}$. Then, there exists a smallest global minimizer $a$ of $f$ and a positive $\epsilon\in\mathbb{R}$ such that for all $t\in\mathbb{R}$ with $a-\epsilon<t<a$, the quadratic polynomial $f_{t}$, defined by

[TABLE]

satisfies $f_{t}\leq f$ over $\mathbb{R}$.

Proof.

The existence of $a$ is straightforward. First, we handle the case when $\deg f=2$. Using Taylor Decomposition of $f$ at $t$, one obtains $f=f(t)+f^{\prime}(t)(X-T)+\frac{f^{\prime\prime}(t)}{2}(X-T)^{2}$. Since $f$ has no real root, the discriminant of $f$ is negative, namely $f^{\prime}(t)^{2}-4f(t)\frac{f^{\prime\prime}(t)}{2}<0$. It implies that $\frac{f^{\prime}(t)^{2}}{4f(t)}<\frac{f^{\prime\prime}(t)^{2}}{2}$, ensuring that the inequality $f_{t}\leq f$ holds over $\mathbb{R}$.

In the sequel, we assume that $\deg f>2$. We can find a neighborhood $U$ as in Lemma 3.2 and without loss of generality, let us suppose that $U=[a-\epsilon_{0},a+\epsilon_{0}]$ for some positive $\epsilon_{0}$, so that $f^{\prime}$ has no real root in $[a-\epsilon_{0},a)$. Then, the inequality $f_{t}(x)\leq f(x)$ holds for all $x,t\in U$. Next, we write $f-f_{t}=\sum_{i=0}^{n}a_{it}x^{i}$, with $a_{it}\in\mathbb{R}$ and $n=\deg f>2$ and define the following function:

[TABLE]

Note that the Cauchy bound (Lemma 2.5) implies that for all $t\in 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 the function $U\to\mathbb{R}:t\mapsto C_{t}$ are bounded. Let $C\in\mathbb{R}$ with $C\geq C_{t}$ for all $t\in U$. Then, for all $t\in 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-\epsilon_{0}<a<a+\epsilon_{0}<C$. Let us define $M:=\min\{f(x)\mid x\in[-C,a-\epsilon_{0}]\}$. By definition, $a$ is the global minimizer of $f$, ensuring that $f(a)<M$. For all $t\in[a-\epsilon_{0},a)$, the quadratic polynomial $f_{t}$ has one real root $N_{t}:=\frac{-2f(t)}{f^{\prime}(t)}+t$. When $t\in[a-\epsilon_{0},a)$ converges to $a$, then $f^{\prime}(t)<0$ converges towards 0 and $-2f(t)$ converges towards $-2f(a)<0$. Thus, the corresponding limit of $N_{t}$ is $+\infty$. In addition, $f_{t}(-C)$ tends to $f_{a}(-C)=f(a)<M$. Therefore, there exists some $\epsilon\in(0,\epsilon_{0}]$ such that for all $t\in(a-\epsilon,a)$, one has $N_{t}\in[C,\infty)$ and $f_{t}(-C)<M$. For all $t\in(a-\epsilon,a)$, we partition $\mathbb{R}$ into five intervals and prove that the inequality $f_{t}\leq f$ holds on each interval:

• •

The inequality $f_{t}\leq f$ holds over $(-\infty,-C]$: it comes from the fact that $f_{t}(-C)<M\leq f(-C)$ and the fact $f-f_{t}$ has no real root in $(-\infty,-C]$.

• •

The inequality $f_{t}\leq f$ holds over $(-C,a-\epsilon_{0}]$: $f_{t}$ is monotonically decreasing over $(-\infty,N_{t}]$. Since one has $-C<a-\epsilon_{0}<C\leq N_{t}$, then $f_{t}$ is monotonically decreasing over $(-C,a-\epsilon_{0}]$. This implies that for all $x\in(-C,a-\epsilon_{0}]$, one has $f_{t}(x)\leq f_{t}(-C)<M\leq f(x)$.

• •

The inequality $f_{t}\leq f$ holds over $[a-\epsilon_{0},a)$: it follows from the fact that $[a-\epsilon_{0},a)\subseteq U$.

• •

The inequality $f_{t}\leq f$ holds over $[a,C)$: $f_{t}$ is monotonically decreasing over $(-\infty,N_{t}]$. Since one has $a<C\leq N_{t}$, then $f_{t}$ is monotonically decreasing over $[a,C)$. Since $a$ is a global minimizer of $f$ and $a\in U$, one has $f_{t}(x)\leq f_{t}(a)\leq f(a)\leq f(x)$ for all $x\in[a,C)$ .

• •

The inequality $f_{t}\leq f$ holds over $[C,\infty]$: the claim is implied by the fact that $f_{t}(N_{t})=0<f(N_{t})$, $N_{t}\in[C,\infty]$ together with the fact that $f-f_{t}$ has no real root in $[C,\infty]$.

∎

Proposition 3.4**.**

Let $K$ be a subfield of $\mathbb{R}$ and $f\in K[X]$ with $\deg f=n\geq 1$. Then $f$ is non-negative on $\mathbb{R}$ if and only if $f$ is a weighted sum of $n$ polynomial squares in $K[X]$, i.e. there exist $a_{1},\dots,a_{n}\in K^{\geq 0}$ and $g_{1},\dots,g_{n}\in K[X]$ such that $f=\sum_{i=0}^{n}a_{i}g_{i}^{2}$.

Proof.

The if part is straightforward. For the other direction, assume that $f$ is non-negative on $\mathbb{R}$ and $n$ is even. The proof is by induction over $n$. The base case $n=2$ follows from Lemma 3.1. For the induction case, let us consider $n\geq 4$.

When $f$ is not a square-free polynomial, we show that $f$ is a weighted sum of $n-2$ polynomial squares. We can write $f=gh^{2}$, for some polynomials $g,h\in K[X]$ with $\deg g\leq\deg f-2$. This gives $g(x)=\frac{f(x)}{h(x)^{2}}\geq 0$ for all $x\in\mathbb{R}$ such that $h(x)\neq 0$. Since $h$ has a finite number of real roots, $g$ is non-negative on $\mathbb{R}$. Using the induction hypothesis, $g$ is a weighted sum of $n-2$ polynomial squares. Therefore, $f$ is also a weighted sum of $n-2$ polynomial squares.

When $f$ is a square-free polynomial, then $f$ has no real root, which implies by Lemma 2.4 that $f$ is neither a square-free polynomial in $K[X]$ nor in $\mathbb{R}[X]$. Thus, $f$ is positive on $\mathbb{R}$. Using Proposition 3.3, there exists some $t\in K$ ($K$ is dense in $\mathbb{R}$) and a quadratic polynomial $f_{t}\in K[X]$ such that the inequalities $0\leq f_{t}(x)\leq f(x)$ holds for all $x\in\mathbb{R}$ and $f_{t}(t)=f(t)$. The polynomial $f-f_{t}$ has degree $n$, takes only non-negative values. In addition $(f-f_{t})(t)=0$, thus $f-f_{t}$ is not a square-free polynomial. Hence, we are in the above case, implying that $f-f_{t}$ is a weighted sum of $n-2$ polynomial squares. From Lemma 3.1, $f_{t}$ is a weighted sum of 2 polynomial squares, implying that $f$ is a weighted sum of $n$ polynomial squares, as requested. ∎

3.2 Algorithm univsos1

The global minimizer $a$ is a real root of $f^{\prime}\in K[X]$. Therefore, by using root isolation techniques [2, Chap. 10], one can isolate all the real roots of $f^{\prime}$ in distinct intervals with bounds in $K$. Such techniques rely on applying successive bisections, so that one can arbitrarily reduce the width of every interval and sort them w.r.t. their lower bounds. Eventually, we apply this procedure to find a sequence of elements in $K$ converging from below to the smallest global minimizer of $f$ in order to find a suitable $t$. We denote by $\texttt{parab}(f)$ the corresponding procedure which returns the polynomial $f_{t}:=\frac{f^{\prime}(t)^{2}}{f(t)}(X-t)^{2}+f^{\prime}(t)(X-t)+f(t)$ such that $t\in K$ and $f\geq f_{t}$ over $\mathbb{R}$.

Algorithm univsos1, depicted in Figure 1, takes as input a polynomial $f\in K[X]$ of even degree $n\geq 2$. The steps performed by this algorithm correspond to what is described in the proof of Proposition 3.4 and relies on two auxiliary procedures. The first one is the procedure parab performing root isolation (see Step 8). The second one is denoted by sqrfree and performs square-free decomposition: for a given polynomial $f\in K[X]$, $\texttt{sqrfree}(f)$ returns two polynomials $g$ and $h$ in $K[X]$ such that $f=gh^{2}$. When $f$ is square-free, the procedure returns $g=f$ and $h=1$ (in this case $\deg h=0$). As in the proof of Proposition 3.4, this square-free decomposition procedure is performed either on the input polynomial $f$ (Step 5) or on the non-negative polynomial $(f-f_{t})$ (Step 9).

The output of Algorithm univsos1 is a pair of lists of polynomials in $K[X]$, allowing to retrieve an SOS decomposition of $f$. By Proposition 3.4 the length of all output lists, denoted by $r$, is bounded by $n/2$. If we note $h_{r},\dots,h_{1}$ the polynomials belonging to h_list and $q_{r},\dots,q_{1}$ the positive quadratic polynomials belonging to q_list, one obtains the following Horner-like decomposition: f=h_{r}^{2}\bigl{(}h_{r-1}^{2}(h_{r-2}^{2}(\dots)+q_{r-2})+q_{r-1}\bigr{)}+q_{r}. Thus, each positive quadratic polynomial $q_{i}$ being a weighted SOS polynomial, this yields a valid weighted SOS decomposition for $f$.

Example 1**.**

Let us consider the polynomial $f:=\frac{1}{16}X^{6}+X^{4}-\frac{1}{9}X^{3}-\frac{11}{10}X^{2}+\frac{2}{15}X+2\in\mathbb{Q}[X]$.

We describe the different steps performed by Algorithm $\texttt{univsos}1$:

• •

The polynomial $f$ is square-free and the algorithm starts by providing the value $t=-1$ as an approximation of the smallest minimizer of $f$. With $f(t)=\frac{1397}{720}$ and $f^{\prime}(t)=\frac{-19}{8}$, one obtains $f_{-1}=\frac{720}{1397}(-\frac{19}{16}X+\frac{271}{360})^{2}$.

• •

Next, after obtaining the square-free decomposition $f(X)-f_{-1}=(X+1)^{2}g$, the same procedure is applied on $g$. One obtains the value $t=1$ as an approximation of the smallest minimizer of $g$ and $g_{1}=\frac{502920}{237293}(-\frac{1}{18}X+\frac{88411}{167640})^{2}$.

• •

Eventually, one obtains the square-free decomposition $g(X)-g_{1}=(X-1)^{2}h$ with $h=\frac{1}{16}(X-\frac{19108973}{17085096})$.

Overall, Algorithm $\texttt{univsos}1$ provides the lists $\texttt{h\_list}=[1,X+1,1,X-1,0]$ and $\texttt{q\_list}=[\frac{720}{1397}(-\frac{19}{16}X+\frac{271}{360})^{2},0,\frac{502920}{237293}(-\frac{1}{18}X+\frac{88411}{167640})^{2},0,\frac{1}{16}(X-\frac{19108973}{17085096})]$, yielding the following weighted SOS decomposition:

[TABLE]

In the sequel, we analyze the complexity of Algorithm univsos1 in the particular case $K=\mathbb{Q}$. We provide bounds on the bitsize of related SOS decompositions as well as bounds on the arithmetic cost required for computation and verification.

3.3 Bit size of the output

Lemma 3.5**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. When applying Algorithm univsos1 to $f$, the sub-procedure parab outputs a polynomial $f_{t}$ such that $\tau(t)=\operatorname{\mathcal{O}}{(n^{2}\tau)}$.

Proof.

Let us consider the set $S\subseteq\mathbb{Q}$ defined by:

[TABLE]

The polynomial involved in $S$ has degree $2n$, with maximum bitsize of the coefficients upper bounded by $2\tau$. Thanks to the complexity analysis of the quantifier elimination procedure described in [2, §11.1.1] the set $S$ can be described by polynomials with maximum bitsize coefficients upper bounded by $\operatorname{\mathcal{O}}{(n^{2}\tau)}$. Since $t$ is a rational root of one of these polynomials, the rational zero theorem [34] implies that $\tau(t)=\operatorname{\mathcal{O}}{(n^{2}\tau)}$. ∎

Lemma 3.6**.**

*Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. When applying Algorithm univsos1 to $f$, the sub-procedure parab outputs a polynomial $f_{t}$ such that $\tau(f_{t})=\operatorname{\mathcal{O}}{(n^{3}\tau)}$. Moreover, there exist polynomials $\hat{f},\hat{f_{t}},g\in\mathbb{Z}[X]$ such that $\hat{f}-\hat{f_{t}}=(X-t)^{2}g$ and $\tau(g)=\operatorname{\mathcal{O}}{(n^{3}\tau)}$. *

Proof.

One can write $f_{t}=M_{2}(t)X^{2}+M_{1}(t)X+M_{0}(t)$ with

[TABLE]

and $\|f_{t}\|_{\infty}\leq\max\{M_{2}(t),|M_{1}(t)|,M_{0}(t)\}$. One has $0\leq M_{0}(t)=f_{t}(0)\leq f(0)\leq\|f\|_{\infty}$.

In addition, $0\leq M_{0}(t)+M_{1}(t)+M_{2}(t)=f_{t}(1)\leq f(1)\leq(n+1)\|f\|_{\infty}$ and $0\leq M_{0}(t)-M_{1}(t)+M_{2}(t)=f_{t}(-1)\leq f(-1)\leq(n+1)\|f\|_{\infty}$. Thus, one has $M_{0}(t)+|M_{1}(t)|+M_{2}(t)\leq(n+1)\|f\|_{\infty}$, wich implies that $\|f_{t}\|_{\infty}\leq(n+1)\|f\|_{\infty}$.

Now let us note $t=\frac{t_{1}}{t_{2}}$, with $t_{1}\in\mathbb{Z},t_{2}\in\mathbb{Z}\backslash\{0\}$, $t_{1}$ and $t_{2}$ being coprime. Let us define the polynomials $\hat{f}(X):=t_{2}^{2n}f(t)f(X)$ and $\hat{f_{t}}(X):=t_{2}^{2n}f(t)f_{t}(X)$. By writing $f(X)=\sum_{i=0}^{n}a_{i}X^{i}$, one has $t_{2}^{2n}f(t)=\sum_{i=0}^{n}a_{i}t_{1}^{i}t_{2}^{2n-i}\leq\|f\|_{\infty}|t_{1}|^{i}|t_{2}|^{2n-i}$ and $\tau(\hat{f})\leq\tau+\tau(t^{2n})$. By Lemma 3.5, one has $\tau(\hat{f})=\operatorname{\mathcal{O}}{(n^{3}\tau)}$.

The polynomials $\hat{f}(X),\hat{f_{t}}(X)$ are polynomials in $\mathbb{Z}[X]$ and since $\|\hat{f_{t}}\|_{\infty}\leq(n+1)\|\hat{f}\|_{\infty}$, the triangular inequality $\|\hat{f}-\hat{f_{t}}\|_{\infty}\leq\|\hat{f}\|_{\infty}+\|\hat{f_{t}}\|_{\infty}$ implies that $\tau(\hat{f}-\hat{f_{t}})\leq\log_{2}(n+2)+\tau(\hat{f})$. In addition, $\tau(f_{t})\leq\tau(\hat{f_{t}})+\tau(t_{2}^{2n}f(t))=\operatorname{\mathcal{O}}{(n^{3}\tau)}$.

As in the proof of Proposition 3.4, one has $(\hat{f}-\hat{f_{t}})(t)=0$ which allows to write the square-free decomposition of the polynomial $\hat{f}-\hat{f_{t}}\in\mathbb{Z}[X]$ as $\hat{f}-\hat{f_{t}}=(X-t)^{2}g$, with $g\in\mathbb{Z}[X]$. By Lemma 2.2, one has $\tau(g)\leq n-2+\tau(\hat{f}-\hat{f_{t}})+\log_{2}(n+1)\leq n-2+2\log_{2}(n+2)+\tau(\hat{f})=\operatorname{\mathcal{O}}{(n^{3}\tau)}$, which concludes the proof. ∎

Theorem 3.7**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then the maximum bitsize of the coefficients involved in the SOS decomposition of $f$ obtained with Algorithm univsos1 is upper bounded by $\operatorname{\mathcal{O}}{((k!)^{3}\tau)}=\operatorname{\mathcal{O}}{((\frac{n}{2})^{\frac{3n}{2}}\tau)}$.

Proof.

With $k=n/2$ and starting from the polynomial $f$, Algorithm univsos1 generates, in the worst case scenario, two sequences of polynomials $f_{k},\dots,f_{1}\in\mathbb{Z}[X],q_{k},\dots,q_{2}\in\mathbb{Z}[X]$ as well as rational numbers $t_{k},\dots,t_{2}\in\mathbb{Q}$ such that $f_{k}=f$, $t_{i}=\frac{t_{i1}}{t_{i2}}$, with $t_{i1}\in\mathbb{Z}$, $t_{i2}\in\mathbb{Z}\backslash\{0\}$ and

[TABLE]

From Lemma 3.6, for all $i=2,\dots,k$, one has $\tau(f_{i-1})=\operatorname{\mathcal{O}}{(i^{3}\tau(f_{i}))}$. This yields \tau(f_{1})=\operatorname{\mathcal{O}}{\bigl{(}(k!)^{3}\tau(f)\bigr{)}}.

Using Stirling’s formula, we obtain k!\leq 2\sqrt{2\pi k}\bigl{(}\frac{k}{e}\bigr{)}^{k} and (k!)^{3}\leq 1024\sqrt{2}\pi^{\frac{3}{2}}k^{\frac{3}{2}}\bigl{(}\frac{k}{e}\bigr{)}^{3k}, where $e$ denotes the Euler number. Since $k\leq e^{k}$ for each integer $k\geq 1$ and $\frac{3}{2}<3$, one has $(k!)^{3}\in\operatorname{\mathcal{O}}{(k^{3k})}$, yielding $\tau(f_{i})=\operatorname{\mathcal{O}}{((\frac{n}{2})^{\frac{3n}{2}}\tau)}$, for all $i=1,\dots,k$. Similarly, we obtain $\tau(q_{i})=\operatorname{\mathcal{O}}{((\frac{n}{2})^{\frac{3n}{2}}\tau)}$, for all $i=1,\dots,k$. Finally, using Lemma 3.5, one has $\tau(t_{i})=\operatorname{\mathcal{O}}{(i^{2}\tau(f_{i}))}$, yielding the desired result. ∎

3.4 Bit complexity analysis

Theorem 3.8**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then, on input $f$, Algorithm univsos1 runs in boolean time

[TABLE]

Proof.

For $i=2,\ldots,k$ we obtain each polynomial $f_{i-1}$ as in the proof of Theorem 3.7 by computing the square-free decomposition of the polynomial $t_{i2}^{4i}f_{i}(t_{i})f_{i}-q_{i}$. It follows by Lemma 2.3 that the polynomial $f_{i-1}$ can be computed using an expected number of $\operatorname{\overset{\sim}{\mathcal{O}}}{(i^{2}\cdot i^{3}\tau(f_{i}))}$ boolean operations. The number of boolean operations to compute all polynomials $f_{1},\dots,f_{k-1}$ is thus bounded by \operatorname{\overset{\sim}{\mathcal{O}}}{\bigl{(}k^{2}\cdot k^{3}\tau(lf)+(k-1)^{2}(k-1)^{3}k^{3}\tau(lf)+\dots+(k!)^{3}\tau(lf)\bigr{)}}.

For each $i=2,\dots,k$, the bitsize of the rational number $t_{i}$ is upper bounded by $\operatorname{\mathcal{O}}{(i^{2}\tau(f_{i}))}$. Therefore, $t_{i}$ can be computed by approximating the roots of $f_{i}^{\prime}$ with isolating intervals of radius less than $2^{-i^{2}\tau(f_{i})}$. By Lemma 2.7, the corresponding computation cost is $\operatorname{\overset{\sim}{\mathcal{O}}}{(i^{3}\tau(f_{i}))}$ boolean operations. The number of boolean operations to compute all rational numbers $t_{2},\dots,t_{k}$ is bounded by \operatorname{\overset{\sim}{\mathcal{O}}}{\bigl{(}k^{3}\cdot k^{3}\tau(lf)+(k-1)^{3}(k-1)^{3}k^{3}\tau(lf)+\dots+(k!)^{3}\tau(lf)\bigr{)}}.

In addition, one has $k^{3}\cdot k^{3}+(k-1)^{3}(k-1)^{3}k^{3}+\dots+(k!)^{3}=(k!)^{3}\sum_{i=1}^{k}\frac{1}{(i!)^{3}}\leq 2k^{3}\cdot(k!)^{3}$. Using Stirling’s formula, we obtain k^{3}\cdot(k!)^{3}\leq 1024\sqrt{2}\pi^{\frac{3}{2}}k^{\frac{9}{2}}\bigl{(}\frac{k}{e}\bigr{)}^{3k}. Since $k^{3/2}\leq e^{k}$ for each integer $k\geq 1$, we obtain the announced complexity. ∎

For a given polynomial $f$ of degree $2k$, one can check the correctness of the SOS decomposition obtained with Algorithm univsos1 by evaluating this SOS polynomial at $2k+1$ distinct points and compare the results with the ones obtained while evaluating $f$ at the same points.

Theorem 3.9**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then one can check the correctness of the SOS decomposition of $f$ obtained with Algorithm univsos1 within

[TABLE]

boolean operations.

Proof.

From [7, Corollary 8.27], the cost of polynomial multiplication in $\mathbb{Z}[X]$ of degree less than $n=2k$ with coefficients of bitsize upper bounded by $B$ is bounded by $\operatorname{\overset{\sim}{\mathcal{O}}}{(k\cdot B)}$. By Theorem 3.7, the maximal bitsize of the coefficients of the SOS decomposition of $f$ obtained with Algorithm univsos1 is upper bounded by $B=\operatorname{\mathcal{O}}{((k!)^{3}\tau)}$. Let us consider $2k+1$ distinct integers, with maximal bitsize upper bounded by $\log_{2}n$. Therefore, from [7, Corollary 10.8], the cost of the evaluation of this decomposition at the $2k+1$ points can be performed using at most $\operatorname{\overset{\sim}{\mathcal{O}}}{(k\cdot(k!)^{3}\tau)}$ boolean operations, the desired result. ∎

Remark 3.10**.**

Let $f_{k}=f\in\mathbb{Z}[X]$. Under the strong assumption that all polynomials $f_{k},\dots,f_{1}$ involved in Algorithm univsos1 have at least one integer global minimizer, then Algorithm univsos1 has a polynomial complexity. Indeed, in this case, $q_{i}=f_{i}(t_{i})$, $\tau(t_{i})=\operatorname{\mathcal{O}}{(\tau(f_{i}))}$ and $\tau(f_{i-1})=\operatorname{\mathcal{O}}{(2(i-1)+\tau(f_{i}))}$, for all $i=2,\dots,k$. Hence, the maximal bitsize of the coefficients involved in the SOS decomposition of $f$ is upper bounded by $\operatorname{\mathcal{O}}{(k^{2}+\tau)}$ and this decomposition can be computed using an expected number of $\operatorname{\overset{\sim}{\mathcal{O}}}{(k^{4}+k^{3}\tau)}$ boolean operations.

4 Nichtnegativstellensätze with perturbed polynomials

Here, we recall the algorithm given in [4, Section 5.2]. The description of this algorithm, denoted by univsos2, is given in Figure 2.

4.1 Algorithm univsos2

Given a subfield $K$ of $\mathbb{R}$ and a non-negative polynomial $f=\sum_{i=0}^{n}f_{i}\,X^{i}\in K[X]$ of degree $n=2k$, one first obtains the square-free decomposition of $f$, yielding $f=p\,h^{2}$ with $p>0$ on $\mathbb{R}$ (see Step 3).Then the idea is to find a positive number $\varepsilon>0$ in $K$ such that the perturbed polynomial $p_{\varepsilon}(X):=p(X)-\varepsilon\sum_{i=0}^{k}X^{2i}$ is also positive on $\mathbb{R}$ (see [4, Section 5.2.2] for more details). This number is computed thanks to the loop going from Step 6 to Step 8 and relies on the auxiliary procedure has_real_roots which checks whether the polynomial $p_{\varepsilon}$ has real roots using root isolation techniques. As mentioned in [4, Section 5.2.2], the number $\varepsilon$ is divided by 2 again to allow a margin of safety (Step 9).

Note that one can always ensure that the leading coefficient ${\ell}:=p_{n}$ of $p$ is the same as the leading coefficients $f_{n}$ of the input polynomial $f$.

We obtain an approximate rational sums of squares decomposition of the polynomial $p_{\varepsilon}$ with the auxiliary procedure sum_two_squares (Step 10), relying on an arbitrary precision complex root finder. Recalling Theorem 2.1, this implies that the polynomial $p$ can be approximated as close as desired by the weighted sum of two polynomial squares in $\mathbb{Q}[X]$, that is ${\ell}s_{1}^{2}+{\ell}s_{2}^{2}$.

Thus there exists a remainder polynomial $u:=p_{\varepsilon}-\\ {\ell}s_{1}^{2}-{\ell}s_{2}^{2}$ with coefficients of arbitrary small magnitude (as mentioned in [4, Section 5.2.3]). The magnitude of the coefficients converges to 0 as the precision $\delta$ of the complex root finder goes to infinity. The precision is increased thanks to the loop going from Step 12 to Step 14 until a condition between the coefficients of $u$ and $\varepsilon$ becomes true, ensuring that $\varepsilon\sum_{i=0}^{k}X^{2i}+u(X)$ also admits a weighted SOS decomposition. For more details, see [4, Section 5.2.4].

The reason why Algorithm univsos2 terminates is the following: at first, one can always find a sufficiently small perturbation $\varepsilon$ such that the perturbed polynomial $p_{\varepsilon}$ remains positive. Next, one can always find sufficiently precise approximations of the complex roots of $p_{\varepsilon}$ ensuring that the error between the initial polynomial $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]$, allowing to retrieve a weighted SOS decomposition of $f$. The size $r$ of both lists is equal to $2k+3=n^{\prime}+3\leq n+3$. If we note $c_{r},\dots,c_{1}$ the numbers belonging to c_list and $s_{r},\dots,s_{1}$ the polynomials belonging to s_list, one obtains the following SOS decomposition: $f=c_{r}s_{r}^{2}+\dots+c_{1}s_{1}^{2}$.

Example 2**.**

Let us consider the same polynomial $f:=\frac{1}{16}X^{6}+X^{4}-\frac{1}{9}X^{3}-\frac{11}{10}X^{2}+\frac{2}{15}X+2\in\mathbb{Q}[X]$ as in Example 1. We describe the different steps performed by Algorithm $\texttt{univsos}1$:

• •

The polynomial $f$ is square-free so we obtain $p=f$ (Step 3). After performing the loop from Step 6 to Step 8, Algorithm $\texttt{univsos}2$ provides the value $\varepsilon=\frac{1}{32}$ at Step 9 as well as the polynomial $p_{\varepsilon}:=p-\frac{1}{32}(1+X^{2}+X^{4}+X^{6})$ which has no real root.

• •

Next, after increasing three times the precision in the loop going from Step 12 to Step 14, the result of the approximate root computation yields $s_{1}=X^{3}-\frac{69}{8}X$ and $s_{2}=7X^{2}-\frac{1}{4}X-\frac{63}{8}$.

Applying Algorithm $\texttt{univsos}2$, we obtain the following two lists of size $6+3=9$:

[TABLE]

yielding the following weighted SOS decomposition:

[TABLE]

4.2 Bit size of the output

First, we need the following auxiliary result:

Lemma 4.1**.**

*Let $p\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg p=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $p$. Then, one has *

[TABLE]

Proof.

Denoting by $\tau^{\prime}$ the maximum bit size of the coefficients of $p^{\prime}$ and instantiating $\alpha=\inf_{x\in\mathbb{R}}p(x)$ with a global minimizer of $p$, $Q$ with $p$ and $A$ with $p^{\prime}$ in the third item of [23, Lemma 3.2], one obtains.

[TABLE]

Now, remark that $\tau^{\prime}\leq\log_{2}n+\tau$. Using this inequality in the one above allows to conclude. ∎

Lemma 4.2**.**

Let $p\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg p=n=2k$ and let $\tau$ be an upper bound on the bitsize of the coefficients of $p$. Then there exists a positive integer $N$ such that for $\varepsilon:=\frac{1}{2^{N}}$, the polynomial $p_{\varepsilon}:=p-\varepsilon\sum_{i=0}^{k}X^{2i}$ is positive over $\mathbb{R}$ and $N=\tau(\varepsilon)\leq\operatorname{\mathcal{O}}{(n\log_{2}n+n\tau)}$.

Proof.

Let us first consider the polynomial $r:=p-\frac{p_{n}}{2}\sum_{i=0}^{k}X^{2i}$. Using [2, Corollary 10.4], the absolute value of each real root of the polynomial $r$ is bounded by $n\,2^{\tau(r)}\leq 2n2^{\tau}$. By defining $R:=2n2^{\tau}$, it follows that the polynomial $r$ is positive for all $|x|>R$. In addition, for all positive integer $N$ and $\varepsilon=\frac{1}{2^{N}}$, one has $\varepsilon\leq\frac{1}{2}\leq\frac{p_{n}}{2}$ and $p_{\varepsilon}=p-\varepsilon\sum_{i=0}^{k}X^{2i}\geq p-\frac{p_{n}}{2}\sum_{i=0}^{k}X^{2i}=r$, which implies that the polynomial $p_{\varepsilon}$ is also positive for all $|x|>R$. Since $R=2n2^{\tau}>1$, one has $1+R^{2}\dots+R^{n}<nR^{n}$. Let us choose the smallest positive integer $N$ such that $nR^{n}\leq 2^{N}\inf_{|x|\leq R}p$. This implies that $\varepsilon<\frac{\inf_{|x|\leq R}p}{1+R^{2}\dots+R^{n}}$, ensuring that the polynomial $p_{\varepsilon}$ is also positive for all $|x|\leq R$. Applying Lemma 4.1, we obtain the following upper bound:

[TABLE]

The announced estimate follows straightforwardly. ∎

In the sequel, we denote by $z_{1},\dots,z_{n}$ the (not necessarily distinct) complex roots of the polynomial $p_{\varepsilon}$. Assuming that we approximate each complex root with a relative precision of $\delta$, we write $\hat{z_{1}},\dots,\hat{z_{n}}$ for the approximate complex root values satisfying $\hat{z_{i}}=z_{i}(1+e_{i})$, with $|e_{i}|\leq 2^{-\delta}$, for all $i=1,\dots,n$.

Theorem 4.3**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then the maximal bitsize of the coefficients involved in the weighted SOS decomposition of $f$ obtained with Algorithm univsos2 is upper bounded by $\operatorname{\mathcal{O}}{(n^{3}+n^{2}\tau)}$.

Proof.

Let $p$ be the square-free part of the polynomial $f$ (see Step 3 of Algorithm univsos2). Then by using Lemma 2.2, one has $\tau(p)\leq n+\tau+\log_{2}(n+1)=\operatorname{\mathcal{O}}{(n+\tau)}$.

Let $\varepsilon=\frac{1}{2^{N}}$ as in Lemma 4.2 so that the polynomial $p_{\varepsilon}=p-\varepsilon\sum_{i=0}^{k}X^{2i}$ is positive over $\mathbb{R}$. By Lemma 4.2, one has $N=C(n^{2}+n\tau)$ for some $C>1$. Let us write $p_{\varepsilon}=\sum_{i=0}^{n}a_{i}X^{i}$ with $a_{n}={\ell}$ and prove that a precision of $\delta:=N+\log_{2}(5n\|p\|_{\infty})=C(n^{2}+n\tau)+\log_{2}(5n\|p\|_{\infty})$ is enough to ensure that the coefficients of $u$ satisfy $\varepsilon\geq\frac{|u_{2i+1}|}{4}-u_{2i}+|u_{2i-1}|$, for all $i=0,\dots,k$. First, note that $e:=2^{-\delta}<\frac{1}{n(n+1)}$ holds. By using Vieta’s formulas provided in Lemma 2.8, one has for all $j=1,\dots,n$:

[TABLE]

Then one has for all $j=1,\dots,n$:

[TABLE]

Since $e<\frac{1}{n}$, one can apply [15, Lemma 3.3], which yields $\prod_{1\leq i_{1}<\dots<i_{j}\leq n}(1+e_{i_{j}})\leq 1+\theta_{j}$, with $|\theta_{j}|\leq\frac{je}{1-je}$. In addition, one has $(j+1)e-\frac{je}{1-je}=\frac{e(1-j(j+1)e)}{1-je}\geq 0$ since $e<\frac{1}{n(n+1)}<\frac{1}{j(j+1)}$, for all $j=1,\dots,n$. Hence, one has $|u_{n-j}|\leq|a_{n-j}|(j+1)e$, for all $j=1,\dots,n$.

This implies that for all $i=0,\dots,k$:

[TABLE]

Since $\delta=N+\log_{2}(5n\|p\|_{\infty})$, one has $5ne\|p\|_{\infty}=\varepsilon$. Thus, for all $i=0,\dots,k$, $\varepsilon\geq\frac{|u_{2i+1}|}{4}-u_{2i}+|u_{2i-1}|$ holds with $\delta=\operatorname{\mathcal{O}}{(n^{2}+n\tau+\log_{2}n+n+\tau)}=\operatorname{\mathcal{O}}{(n^{2}+n\tau)}$.

For each $j=1,\dots,n$, choosing $e_{j}=e=2^{-\delta}$ and $\hat{z}_{j}=z_{j}(1+2^{-\delta})$, yields $|u_{n-j}|=|a_{n-j}||1-(1+2^{-\delta})^{j}|$. Next, we bound the size of the weighted SOS decomposition. One has $\tau(\delta)=\operatorname{\mathcal{O}}{(n^{2}+n\tau)}$ and for all $i=1,\dots,n$, $\tau(a_{n-i})\leq\tau(\varepsilon)=\operatorname{\mathcal{O}}{(n^{2}+n\tau)}$. Therefore, for all $j=1,\dots,n$, $\tau(u_{n-j})\leq\operatorname{\mathcal{O}}{(n^{2}+n\tau+j(n^{2}+n\tau))}$ and the maximal bitsize of the coefficients of $u$ is bounded by $\operatorname{\mathcal{O}}{(n^{3}+n^{2}\tau)}$.

From Lemma 2.6, one has $|\hat{z}_{j}|=|z_{j}|(1+2^{-\delta})\geq\frac{1}{2^{\tau(p_{\varepsilon})}+1}(1+2^{-\delta})|$, so that it is enough to perform root isolation for the polynomial $p_{\varepsilon}$ with a precision bounded by $\operatorname{\mathcal{O}}{(\tau(p_{\varepsilon})+\delta)}=\operatorname{\mathcal{O}}{(n^{2}+n\tau)}$.

Finally, the weighted SOS decomposition of $f$ has coefficients of maximal bitsize bounded by $\operatorname{\mathcal{O}}{(n^{3}+n^{2}\tau)}$ as claimed. ∎

4.3 Bit complexity analysis

Theorem 4.4**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then, on input $f$, Algorithm univsos2 runs in boolean time

[TABLE]

Proof.

By Lemma 2.3, the square-free decomposition of $f$ can be computed using an expected number of $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{2}\tau)}$ boolean operations. Checking that the polynomial $p_{\varepsilon}$ has no real root can be performed using an expected number of $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{2}\cdot\tau(\varepsilon))}=\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{3}\tau)}$ boolean operations while relying on Sylvester-Habicht Sequences [20, Corollary 5.2].

As seen in the proof of Theorem 4.3, the complex roots of $p_{\varepsilon}$ must be approximated with isolating intervals (resp. disks) of radius less than $2^{-\tau(p_{\varepsilon})-\delta}$. Thus, by Lemma 2.7, all real (resp. complex) roots of $p_{\varepsilon}$ can be computed in $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{3}+n^{2}\tau(p_{\varepsilon})+n(\delta+\tau(p_{\varepsilon}))}=\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{4}+n^{3}\tau)}$ boolean operations.

As in the proof of Theorem 4.3, one can select $|u_{n-j}|=|a_{n-j}||1-(1+2^{-\delta})^{j}|$, for all $j=1,\dots,n$. This implies that the computation of each coefficient of $u$ can be performed with at most $\operatorname{\overset{\sim}{\mathcal{O}}}{(n\cdot\tau(\delta))}=\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{3}+n^{2}\tau)}$ boolean operations. Eventually, we obtain a bound of $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{4}+n^{3}\tau)}$ for the computation of all coefficients of $u$, which yields the desired result. ∎

We state now the complexity result for checking the SOS certificates output by Algorithm univsos2. As for the output of Algorithm univsos1, this is done through evaluation of the output at $n+1$ distinct values where $n$ is the degree of the output.

Theorem 4.5**.**

Let $f\in\mathbb{Z}[X]$ be a positive polynomial over $\mathbb{R}$, with $\deg f=n=2k$ and $\tau$ be an upper bound on the bitsize of the coefficients of $f$. Then one can check the correctness of the weighted SOS decomposition of $f$ obtained with Algorithm univsos2 using $\operatorname{\overset{\sim}{\mathcal{O}}}{(n^{4}+n^{3}\tau)}$ bit operations.

Proof.

From [7, Corollary 8.27], the cost of polynomial multiplication in $\mathbb{Z}[X]$ of degree less than $n$ with coefficients of bitsize upper bounded by $l$ is bounded by $\operatorname{\overset{\sim}{\mathcal{O}}}{(n\cdot l)}$. By Theorem 4.3, the maximal coefficient bitsize of the SOS decomposition of $f$ obtained with Algorithm univsos2 is upper bounded by $l=\operatorname{\mathcal{O}}{(n^{3}+n^{2}\tau)}$. Therefore, from [7, Corollary 10.8], the cost of the evaluation of this decomposition at $n$ points can be performed using at most $\operatorname{\overset{\sim}{\mathcal{O}}}{(n\cdot(n^{3}+n\tau))}$ boolean operations as claimed. ∎

5 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 and benchmark execution on the dedicated webpage.333https://github.com/magronv/univsos This tool is integrated to the RAGlib Maple package444http://www-polsys.lip6.fr/~safey/RAGLib/. 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 been implemented using the PARI/GP software tool555http://pari.math.u-bordeaux.fr and is freely available.666https://hal.archives-ouvertes.fr/ensl-00445343v2 To ensure fair comparison, we have 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 tried to use the Maple procedure fsolve but the polroots routine from Pari/GPyields significantly better performance for the polynomials involved in our examples.

The nine polynomial benchmarks presented in Table 1 allow to approximate some given mathematical functions, considered in [4, Section 6]. Computation and verification of SOS certificates are a mandatory step required to validate the supremum norm of the difference between such functions and their respective approximation polynomials on given closed intervals. This boils down to certify two inequalities of the form $\forall x\in[b,c],p(x)\geq 0$, with $p\in\mathbb{Q}[X]$, $b,c\in\mathbb{Q}$ and $\deg p=n$. As recalled in [4, Section 5.2.5], this latter problem can be addressed by computing a weighted SOS decomposition of the polynomial q(Y):=(1+Y^{2})^{n}\,p\Bigl{(}\frac{b+cY^{2}}{1+Y^{2}}\Bigr{)}, with either Algorithm univsos1 or Algorithm univsos2. For each benchmark, we indicate in Table 1 the degree $n$ and the bitsize $\tau$ of the input polynomial, the bitsize $\tau_{1}$ of the weighted SOS decomposition provided by Algorithm univsos1 as well as the corresponding computation (resp. verification) time $t_{1}$ (resp. $t_{1}^{\prime}$). Similarly, we display $\tau_{2},t_{2},t_{2}^{\prime}$ for Algorithm univsos2. The table results show that for all other eight benchmarks, Algorithm univsos2 yields better certification and verification performance, together with more concise SOS certificates. This observation confirms what we could expect after comparing the theoretical complexity results from Sections 3 and 4.

The comparison results available in Table 2 are obtained for power sums of increasing degrees. For a given natural integer $n=2k$ with $10\leq n\leq 500$, we consider the polynomial $P_{n}:=1+X+\dots+X^{n}$. The roots of this polynomial are the $n+1$-th roots of unity, thus yielding the following SOS decomposition with real coefficients: $P_{n}:=\prod_{j=1}^{k}((X-\cos\theta_{j})^{2}+\sin^{2}\theta_{j})$, with $\theta_{j}:=\frac{2j\pi}{n+1}$, for each $j=1,\dots,k$. By contrast with the benchmarks from Table 1, Table 2 shows that Algorithm univsos1 produces output certificates of much smaller size compared to Algorithm univsos2, with a bitsize ratio lying between $6$ and $38$ for values of $n$ between $10$ and $200$. This is due to the fact that Algorithm univsos1 outputs a value of t equal to 0 at each step. The execution performance of Algorithm univsos1 are also much better in this case. The lack of efficiency of Algorithm univsos2 is due to the computational bottleneck occurring in order to obtain accurate approximation of the relatively close roots $\cos\theta_{j}\pm i\sin\theta_{j}$, $j=1,\dots,k$. For $n\geq 300$, execution of Algorithm univsos2 did not succeed after two hours of computation, as meant by the symbol $-$ in the corresponding line.

Further experiments are summarized in Table 3 for modified Wilkinson polynomials $W_{n}$ of increasing degrees $n=2k$ with $10\leq n\leq 600$ and $W_{n}:=1+\prod_{j=1}^{k}(X-j)^{2}$. The complex roots $j\pm i$, $j=1,\dots,k$ of $W_{n}$ are relatively close, which yields again a significant lack of performance of Algorithm univsos2. As observed in the case of power sums, timeout behaviors occur for $n\geq 60$. In addition, the bitsize of the SOS decompositions returned by Algorithm univsos1 are much smaller. This is a consequence of the fact that in this case, $a=1$ is the global minimizer of $W_{n}$. Hence, the algorithm always terminates at the first iteration by returning the trivial quadratic approximation $f_{t}=f_{a}=1$ together with the square-free decomposition of $W_{n}-f_{t}=\prod_{j=1}^{k}(X-j)^{2}$.

Finally, we consider experimentation performed on modified Mignotte polynomials defined by $M_{n,m}:=X^{n}+2(101X-1)^{m}$ and $N_{n}:=(X^{n}+2(101X-1)^{2})(X^{n}+2((101+\frac{1}{101})X-1)^{2})$, for even natural integers $n$ and $m\leq 2$. The corresponding results are displayed in Table 4 for $M_{n,m}$ with $m=2$ and $10\leq n\leq 10000$, $m=n-2$ and $10\leq n\leq 100$ as well as for $N_{n}$ with $10\leq n\leq 100$. Note that similar benchmarks are used in [33] to anayze the efficiency of (real) root isolation techniques over polynomial with relatively close roots. As for modified Wilkinson polynomials, Algorithm univsos2 can only handle small size instances, due to limited scalablity of the polroots procedure. In this case $a=\frac{1}{100}$ is the unique global minimizer of $M_{n,2}$. Thus, Algorithm univsos1 always outputs weighted SOS decompositions of polynomials $M_{n,2}$ within a single iteration by first computing the quadratic polynomial $f_{t}=f_{a}=2(101X-1)^{2}$ and the trivial 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-2}$ and $N_{n}$ with moderate degree (less than 100).

6 Conclusion and perspectives

We presented and analyzed two different algorithms univsos1 and univsos2 to compute weighted sums of squares (SOS) decompositions of non-negative univariate polynomials. When the input polynomial has rational coefficients, one feature shared by both algorithms is their ability to provide non-negativity certificates whose coefficients are also rational. Our study shows that the complexity analysis of Algorithm univsos1 yields an upper bound which is exponential w.r.t. the input degree, while the complexity of Algorithm univsos2 is polynomial. However, comparison benchmarks emphasize the need for both algorithms to handle various classes of non-negative polynomials, e.g. in the presence of rational global minimizers or when root isolation can be performed efficiently.

A first direction of further research is a variant of Algorithm univsos2 where one would compute approximate 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 the polynomials is small, e.g. for power sums of degree up to 1000. However, the performance decrease when the 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 obtain precise values of eigenvalues/vectors of SDP output matrices. Another advantage of this technique is its ability to perform global polynomial optimization. A topic of interest would be to obtain the same feature with the two current algorithms. We also plan to develop extensions to the non-polynomial case.