ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 24 Apr 2013 18:27:44 +0200Polynomials as a sum of squareshttps://ask.sagemath.org/question/10062/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$$
Wed, 24 Apr 2013 14:53:09 +0200https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/Comment by Jaakko Seppälä for <p>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 </p>
<p>$$\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$$</p>
https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?comment=17810#post-id-17810Sorrry. I meant decompose the polynomial as sum of squares.Wed, 24 Apr 2013 15:14:21 +0200https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?comment=17810#post-id-17810Comment by vdelecroix for <p>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 </p>
<p>$$\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$$</p>
https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?comment=17812#post-id-17812You 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 ?Wed, 24 Apr 2013 15:09:09 +0200https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?comment=17812#post-id-17812Answer by tmonteil for <p>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 </p>
<p>$$\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$$</p>
https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?answer=14838#post-id-14838Hi,
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](http://www.stanford.edu/class/ee364b/lectures/sos_slides.pdf), 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)
LC = P.leading_coefficient()
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)
LC = P.leading_coefficient()
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.
Wed, 24 Apr 2013 18:27:44 +0200https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?answer=14838#post-id-14838Answer by vdelecroix for <p>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 </p>
<p>$$\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$$</p>
https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?answer=14837#post-id-14837There 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.Wed, 24 Apr 2013 18:16:00 +0200https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/?answer=14837#post-id-14837