# Polynomials as a sum of squares

Is it possible find a decomposition of a polynomial as a sum of squares in Sage if such representation is possible? For example, if I want to prove that $x^6 - x^5 + x^4 - x^3 + x^2 - x + 2/5>0$ for all $x\in\mathbb{R}$ then Sage would return for example

$$\left (x^2\left (x - \frac{1}{2}\right)\right )^2+\left (\frac{\sqrt{3}x}{2}\left (x - \frac{2}{3}\right)\right )^2+\left (\sqrt{\frac{2}{3}}\left(x - \frac{3}{4}\right )\right )^2+\sqrt{\frac{1}{40}}^2$$

edit retag close merge delete

1

You want to prove that a certain polynomial is positive or to decompose it as a sum of squares ? For the first problem you may compute the number of real roots (which is algorithmically simple). For the second one is there an algorithm to do it ?

( 2013-04-24 15:09:09 +0100 )edit

Sorrry. I meant decompose the polynomial as sum of squares.

( 2013-04-24 15:14:21 +0100 )edit

Sort by » oldest newest most voted

There is nothing "out of the box" within Sage but it is realtively straightforward to design an algorithm. First of all, given a polynomial you need to compute its roots which can be done as follows

sage: R.<x> = PolynomialRing(QQ,'x')
sage: P = x^6 - x^5 + x^4 - x^3 + x^2 - x + 2/5
sage: r = P.roots(CDF)
sage: print r
[(-0.54474977598 - 0.816253990699*I, 1),
(-0.54474977598 + 0.816253990699*I, 1),
(0.355140182797 - 0.844015166249*I, 1),
(0.355140182797 + 0.844015166249*I, 1),
(0.689609593183 - 0.140734077916*I, 1),
(0.689609593183 + 0.140734077916*I, 1)]


Ideally, you would use QQbar (algebraic numbers: exact computation) instead of CDF (complex double field: approximate computation) but I found that it was just impossible to do the computation within QQbar... As you see in the result above, the roots are all complex (which proves that your polynomial is positive as its dominant term is positive). Moreover, the list is well ordered by conjugates and you can rewrite your polynomial as a product of polynomials of degree 2 as follows

sage: (r0,m0),(r1,m1),(r2,m2),(r3,m3),(r4,m4),(r5,m5) = r
sage: P1 = x^2 - (r0 + r1) * x + (r0*r1)
sage: P2 = x^2 - (r2 + r3) * x + (r2*r3)
sage: P3 = x^2 - (r4 + r5) * x + (r4*r5)
sage: print P1 * P2 * P3
x^6 - x^5 + x^4 - x^3 + x^2 - x + 0.4


The above form is not that useful, but each of the polynomial P1, P2, P3 may be rewritten as a sum of two squares :

sage: a1 = (r0+r1)/2
sage: b1 = (r0*r1 - (r0+r1)**2/4).real().sqrt()
sage: a2 = (r2+r3)/2
sage: b2 = (r2*r3 - (r2+r3)**2/4).real().sqrt()
sage: a3 = (r4+r5)/2
sage: b3 = (r4*r5 - (r4+r5)**2/4).real().sqrt()
sage: ((x-a1)**2 + b1**2) * ((x-a2)**2 + b2**2) * ((x-a3)**2 + b3**2)
x^6 - x^5 + x^4 - x^3 + x^2 - x + 0.4


The form $((x-a_1)^2 + b_1^2) ((x-a_2)^2 + b_2^2) ((x-a_3)^2 + b_3^2)$ is not exactly the one you asked for but you can obtain it using the identities $(a^2 + b^2)(c^2 + d^2) = (ac - bd)^2 + (ad + bc)^2$ several times.

more

Hi,

here is a possible solution, which is exact (the coefficients are assumed to be algebraic real numbers). The algorithm comes from the example given in those slides, page 12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]

def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True

def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)


To use it:

sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field


I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.

more