Positive polynomial as a sum of squares

asked 2018-01-26 20:53:40 +0200

mathhobbyist
123 ●3 ●8 ●17

updated 2018-01-27 08:34:44 +0200

I found a question to prove the polynomial $58x^{10}-42x^9+11x^8+42x^7+53x^6-160x^5+118x^4+22x^3-56x^2-20x+74$ to be positive on $\mathbb{R}$. The solution was $(7x^5-4x^3+6x^2+2x-5)^2+(-3x^5+7x^4-3x^3+7)^2$. Is there a method in Sage to find such a representations? I tried the code modified from https://ask.sagemath.org/question/100...

R = AA[x]
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)
R = AA[x]
P = R(58*x^10-42*x^9+11*x^8+42*x^7+53*x^6-160*x^5+118*x^4+22*x^3-56*x^2-20*x+74)
A,B = sum_of_two_squares(P)
A
7.615773105863908?*x^5 - 2.75743509005418?*x^4 - 44.2356981337271?*x^3 + 11.9854478203666?*x^2 + 25.5052589522359?*x - 3.05770635146248?

It looks like the code does not always find the simplest form of solution. So is there an algorithm that tries to find a representation such that it first tries to represent the polynomial as a sum of two constants, then sum of two polynomials of degree at most 1, then sum of two polynomials of degree at most 2 and so on and in every case the one can see the minimal polynomial of the coefficients over $\mathbb Q[x]$.

edit retag flag offensive close merge delete

Comments

Which is the question exactly?

A decomposition $P=A^2+B^2$ over integers is needed?
Do we really need the above code?
An algorithm is needed? If yes, for what, for which problem specifically? The algorithm should work in general, or only in case there is a "simpl decomposition"?

dan_fulea ( 2018-01-27 13:37:14 +0200 )edit

I saw problems where was given polynomials with real coefficients and one variable and was asked to prove them to be positive. I was wondering a method to show them positive by completing them to a sum of squares. But I was not sure in what cases that gives a full proof as the coefficients are not always in closed form. I managed to do such proofs using Samuelson's inequality and then Sturm's theorem but I was wondering in which cases the sum of squares method provides a shorter proof.

mathhobbyist ( 2018-01-27 21:31:26 +0200 )edit

None of the three questions above is answered, moreover, some new actors come into play, making the situation more unclear. Let us go back to the question, as is stated. To see why, just take a look at the <-1> i have got because somebody considered my answer to be misleading, or unrelated, or offensive, or mathematically wrong, or possibly abusively using bold face words, or all together. Not my feeling, so this is the reason why i still politely insist to make things clear. The points again, starting a new try to make things clear:

"a sum of squares" means "a sum of two squares"?
completing (them to a sum of squares) makes no sense above, maybe "writing" is the intention.
"gives a full proof" for what above? Clearly $A^2+B^2+\dots Z^2\ge 0$ for real polynomials involved.

dan_fulea ( 2018-01-28 18:32:39 +0200 )edit

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2018-01-27 00:55:53 +0200

dan_fulea
5337 ●5 ●43 ●90

The question is not really precise at the points that matter for an algorithm. Please let me complete with some words...

So is there an algorithm that returns the representation as sum of two (?) polynomials with coefficients (which) are (in) some "simple" form - like (integers or) roots of some integer coefficient polynomial with degree not too large?

In this situation, i'd like to be more precise, thus fixing the frame:

We give a polynomial $P\in\mathbb Z[x]$ and know it has a simple representation as a sum of two polynomials $A,B\in\mathbb Z[x]$. (Any real polynomial in one variable takes non-negative values on $\mathbb R$ if and only if it is a sum of two squares of real polynomials.)

A first question is now if this representation is unique, else we cannot ask for "the" representation. Simple examples show that we may have more representations for more reasons, e.g. $65 = 1^2+8^2 = 7^2 + 4^2$, or $65x^2 = x^2+(8x)^2 = (7x)^2 + (4x)^2$, and it is hard to make sentences, that eliminate this "trivial case", but even then one may come with the not so trivial, but in an other sense trivial example $(5x)^2 + 5^2 = 25x^2 + 25 =(3x-4)^2+(4x+3)^2. So i think it is clear that we cannot assume uniqueness.

Then the relation $$ P = A^2 + B^2 $$ is a relation in $\mathbb Q[x]$, which transferred to the ring $R=\mathbb Q[i][x]$ becomes: $$ P = (A+iB)(A-iB) \ ,$$ so an algorithm would be to factor over this ring!

Sage can do this. From this decomposition and from each choice of an element $u+iv$ of norm one, $u,v\in \mathbb Q$, $u^2 + v^2 = 1$, we can then use the "contorsion" $P=(u+iv)(A+iB)\cdot (u-iv)(A-iB)$ to get some other solution $A'+iB'=(u+iv)(A+iB)$ instead of the $A+iB$.

In the given case we can proceed as follows:

K.<j> = QuadraticField( -1 )
R.<x> = PolynomialRing( K, sparse=False )
P = R( 58*x^10 - 42*x^9 + 11*x^8 + 42*x^7 + 53*x^6 - 160*x^5 + 118*x^4 + 22*x^3 - 56*x^2 - 20*x + 74 ) / 58
# P.factor()
print "In our case we have exactly two simple factors, so we take f to be the first one..." 
print "The the other one, g, is its Galois conjugate..." 
f = P.factor()[0][0]    # the first factor, taken without multiplicity
fcoeffs = f.coefficients()
gcoeffs = [ c.galois_conjugate() for c in fcoeffs ]
g = R( gcoeffs )

print "The factors are:"
print "f = %s" % f
print "g = %s" % g
print "\nWe associate A = (f+g)/2 and B = (f-g)/2/j and get"
A = (f+g)/2
B = (f-g)/2/j
print "A = %s" % A
print "B = %s" % B

print "Final check: Is P == A^2 + B^2  ? %s" % bool( P == A^2+B^2 )

I took the monic version of $P$. Results:

In our case we have exactly two simple factors, so we take f to be the first one...
The the other one, g, is its Galois conjugate...
The factors are:
f = x^5 + (-49/58*j - 21/58)*x^4 + (33/58*j - 19/58)*x^3 + (-9/29*j + 21/29)*x^2 + (-3/29*j + 7/29)*x - 17/29*j - 28/29
g = x^5 + (49/58*j - 21/58)*x^4 + (-33/58*j - 19/58)*x^3 + (9/29*j + 21/29)*x^2 + (3/29*j + 7/29)*x + 17/29*j - 28/29

We associate A = (f+g)/2 and B = (f-g)/2/j and get
A = x^5 - 21/58*x^4 - 19/58*x^3 + 21/29*x^2 + 7/29*x - 28/29
B = -49/58*x^4 + 33/58*x^3 - 9/29*x^2 - 3/29*x - 17/29
Final check: Is P == A^2 + B^2 ? True

To find the "other solution", we have to combine in an other way...

In this case we knew there is a "beautiful solution", but in general...

I will pick an other example, where the positivity is clear, but where we cannot expect to factorize over the Gaussian integers. $$ P = x^6 + 4x^4 + 4x + 2= (x^2)^3 +(2x^2-1)^2+(2x+1)^2\ . $$ Then

sage: P = R( x^6 + 4*x^4 + 4*x + 2 )
sage: P.factor()
x^6 + 4*x^4 + 4*x + 2

So we need a bigger ring. The following still works:

R.<x> = PolynomialRing( QQbar, sparse=False )
P = R( x^6 + 4*x^4 + 4*x + 2 )
P.factor()
f1, f2, f3, f4, f5, f6 = [ f[0] for f in P.factor() ]

Then we have the polynomials:

sage: f1
x - 0.7197390712152117? - 0.8347870198157105?*I
sage: f2
x - 0.7197390712152117? + 0.8347870198157105?*I
sage: f3
x + 0.10620918972850244? - 2.047873213713062?*I
sage: f4
x + 0.10620918972850244? + 2.047873213713062?*I
sage: f5
x + 0.6135298814867092? - 0.12275817772089546?*I
sage: f6
x + 0.6135298814867092? + 0.12275817772089546?*I

see how the complex conjugation "builds pairs", consider then:

sage: f = f1*f3*f5
sage: g = f2*f4*f6
sage: f
x^3 + (0.?e-18 - 3.005418411249669?*I)*x^2 + (-2.516269913339240? - 0.3080100693808800?*I)*x - 0.9256991333675843? + 1.069149715653755?*I
sage: g
x^3 + (0.?e-18 + 3.005418411249669?*I)*x^2 + (-2.516269913339240? + 0.3080100693808800?*I)*x - 0.9256991333675843? - 1.069149715653755?*I

and obtain in a similar way two polynomials $A,B$. But:

sage: A = f+g
sage: A

already takes a looong time. Of course, understanding such coefficients is possibly / in general as hard as understanding the Galois closure of the splitting field of $P$.

edit flag offensive delete link

add a comment

answered 2020-10-25 20:20:39 +0200

Max Alekseyev

5701 ●5 ●36 ●122

If one just wants to test whether a given polynomial $f(x)$ is positive for all real $x$, this can be done without representing $f(x)$ as the sum of squares, using the following function:

def is_pol_positive(f):
    return f.leading_coefficient() > 0 and f.number_of_real_roots() == 0

edit flag offensive delete link

add a comment

Your Answer

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

Add Answer

Positive polynomial as a sum of squares

Comments

2 Answers

Your Answer

Question Tools

Stats

Positive polynomial as a sum of squares edit

Comments

2 Answers

Your Answer

Question Tools

Stats

Positive polynomial as a sum of squares