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.Sun, 25 Oct 2020 20:20:39 +0100Positive polynomial as a sum of squareshttps://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/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/10062/polynomials-as-a-sum-of-squares/
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]$.Fri, 26 Jan 2018 20:53:40 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/Comment by dan_fulea for <p>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 <a href="https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/">https://ask.sagemath.org/question/100...</a></p>
<pre><code>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?
</code></pre>
<p>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]$.</p>
https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40802#post-id-40802Which 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"?Sat, 27 Jan 2018 13:37:14 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40802#post-id-40802Comment by mathhobbyist for <p>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 <a href="https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/">https://ask.sagemath.org/question/100...</a></p>
<pre><code>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?
</code></pre>
<p>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]$.</p>
https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40815#post-id-40815I 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.Sat, 27 Jan 2018 21:31:26 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40815#post-id-40815Comment by dan_fulea for <p>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 <a href="https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/">https://ask.sagemath.org/question/100...</a></p>
<pre><code>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?
</code></pre>
<p>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]$.</p>
https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40828#post-id-40828None 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.Sun, 28 Jan 2018 18:32:39 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?comment=40828#post-id-40828Answer by dan_fulea for <p>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 <a href="https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/">https://ask.sagemath.org/question/100...</a></p>
<pre><code>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?
</code></pre>
<p>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]$.</p>
https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?answer=40793#post-id-40793The 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$.
Sat, 27 Jan 2018 00:55:53 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?answer=40793#post-id-40793Answer by Max Alekseyev for <p>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 <a href="https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/">https://ask.sagemath.org/question/100...</a></p>
<pre><code>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?
</code></pre>
<p>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]$.</p>
https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?answer=54028#post-id-54028If 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
Sun, 25 Oct 2020 20:20:39 +0100https://ask.sagemath.org/question/40790/positive-polynomial-as-a-sum-of-squares/?answer=54028#post-id-54028