Search for squares

Question

Could anyone please help me find $a\left(x,y,z\right)$, $b_{i}\left(x,y,x\right),$$c_{i}\left(x,y,z\right)$ and $d\left(x,y,z\right)$ that are sums of squares of polynomials or rational functions such that

$$ \begin{array}a a\left(x,y,z\right)+b_{1}\left(x,y,z\right)\left(x+y\right)+b_{2}\left(x,y,z\right)\left(x+z\right)+b_{3}\left(x,y,z\right)\left(y+z\right)\\ +c_{1}\left(x,y,z\right)\left(4-x^{2}\right)+c_{2}\left(x,y,z\right)\left(4-y^{2}\right)+c_{3}\left(x,y,z\right)\left(4-z^{2}\right)\\ +d\left(x,y,z\right)\left(4+xyz-x^{2}-y^{2}-z^{2}\right) \end{array} $$

equals $15(2+x+y)^{2}(2+x+z)^{2}(2+y+z)^{2}-32(3+x+y+z)^{3}$ by Mathematica programming? or could anyone use the Mathematica codes on this page : Expressing a polynomial as a sum of squares to find $a\left(x,y,z\right)$, $b_{i}\left(x,y,x\right),$$c_{i}\left(x,y,z\right)$ and $d\left(x,y,z\right)$ for this polynomial?

Here are some Mathemtica codes that I wrote and you may use for the search by Mathematica programming:

  p = 15 (2 + x + y)^2 (2 + x + z)^2 (2 + y + z)^2 - 32 (3 + x + y + z)^3
 q11 = x + y
 q12 = x + z
 q13 = y + z
 q21 = 4 - x^2
 q22 = 4 - y^2
 q23 = 4 - z^2
 q3 = -4 + x^2 + y^2 - x y z + z^2

Thanks a lot.

Answer 1

The polynomial

$$ p(x,y,z) = 15 (x+y+2)^2 (x+z+2)^2 (y+z+2)^2-32 (x+y+z+3)^3 $$

has a Newton polytope hull which is

NP = {{4, 2, 0}, {4, 0, 2}, {4, 0, 0}, {0, 4, 2}, {0, 4, 0}, {0, 2, 4}, {0, 0, 4}, {0, 0, 0}, {2, 4, 0}, {2, 0, 4}}

has all the powers even. Follows a picture for the Newton polytope. In black, the hull.

The monomials that can generate this set of points (red and black) are

$$ Z = \left\{1,z,z^2,y,y z,y z^2,y^2,y^2 z,x,x z,x z^2,x y,x y z,x y^2,x^2,x^2 z,x^2 y\right\} $$

now

$$ \cases{ a = Z^T\cdot A\cdot Z\\ b_1 = Z^T\cdot B_1\cdot Z\\ \vdots\\ c_3 = Z^T\cdot C_3\cdot Z\\ d = Z^T\cdot D\cdot Z } $$

are feasible candidates with $\{A,B_1,\cdots,C_3,D\}$ positive definite matrices. The next step is to determine the conditions such that all the monomials for $p_0(x,y,z)$

$$ \begin{array}a p_0(x,y,z) = a\left(x,y,z\right)+b_{1}\left(x,y,z\right)\left(x+y\right)+b_{2}\left(x,y,z\right)\left(x+z\right)+b_{3}\left(x,y,z\right)\left(y+z\right)\\ +c_{1}\left(x,y,z\right)\left(4-x^{2}\right)+c_{2}\left(x,y,z\right)\left(4-y^{2}\right)+c_{3}\left(x,y,z\right)\left(4-z^{2}\right)\\ +d\left(x,y,z\right)\left(4+xyz-x^{2}-y^{2}-z^{2}\right) \end{array} $$

are contained in the former Newton polytope. Those linear conditions $R(A,B_1,\cdots,C_3,D)$ are obtained as

$$ p(x,y,z) - p_0(x,y,z) = 0, \ \ \forall \{x,y,z\} $$

so the problem is reduced to:

Determine $\{A,B_1,\cdots,C_3,D\}$ positive definite, subjected to $R(A,B_1,\cdots,C_3,D)$

Resuming, we need some procedures to handle the symbolic, and a positiveness solver like CXV.

NOTE

The Newton polytope can be extracted with the script.

ExtractElements[f_, vars_] := Module[{rf, rf0, nrf, ef, cf},
rf = CoefficientRules[f, vars];
nrf = Length[rf];
If[nrf == 1, rf0 = rf[[1]], rf0 = rf];
ef = Map[First, rf0];
cf = Map[Last, rf0];
Return[{ef, cf}]]

vars = {x, y, z};
pol = 15 (2 + x + y)^2 (2 + x + z)^2 (2 + y + z)^2 - 32 (3 + x + y + z)^3;

{elems, rels} = ExtractElements[pol, vars]

and the convex hull can be obtained using (found in a repository)

ConvexDepenentQ[corners_, cand_] := Module[{w, ws}, w = Array[ws, Length@corners];
1 == Length@FindInstance[w.corners == cand && Total[w] == 1 && And @@ Table[w[[i]] >= 0, {i, Length@w}], w]];

ConvexReduce[data_] := Module[{corners, ncorners, test}, corners = data;
Do[ncorners = Delete[corners, Position[corners, data[[i]]]];
test = ConvexDepenentQ[ncorners, data[[i]]];
If[test, corners = ncorners];, {i, Length@data}];corners];

convexHull[data_] := Module[{corners, rd}, corners = {};
Do[corners = Join[corners, Select[data, Min[data[[;; , i]]] == #[[i]] ||  Max[data[[;; , i]]] == #[[i]] &]];, {i, Length@data[[1]]}];
corners = DeleteDuplicates@corners;
rd = Delete[data, First@Position[data, #] & /@ corners];
Do[If[ConvexDepenentQ[corners, rd[[i]]], , AppendTo[corners, rd[[i]]]], {i, Length@rd}];
ConvexReduce@DeleteDuplicates@corners];


CH = convexHull[elems]

The linear restrictions can be obtained with ExtractElements in the output rels