Ask Your Question
2

Polynomials as a sum of squares

asked 2013-04-24 07:53:09 -0500

Jaakko Seppälä gravatar image

updated 2013-04-24 09:03:03 -0500

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 flag offensive close merge delete

Comments

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 ?

vdelecroix gravatar imagevdelecroix ( 2013-04-24 08:09:09 -0500 )edit

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

Jaakko Seppälä gravatar imageJaakko Seppälä ( 2013-04-24 08:14:21 -0500 )edit

2 answers

Sort by » oldest newest most voted
2

answered 2013-04-24 11:16:00 -0500

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.

edit flag offensive delete link more
1

answered 2013-04-24 11:27:44 -0500

tmonteil gravatar image

updated 2013-04-24 11:59:03 -0500

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)
    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.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2013-04-24 07:53:09 -0500

Seen: 410 times

Last updated: Apr 24 '13