# How to write a positive function as sum of positive terms only

A multivariate polynomial function can be written in a certain form (x^m (1-x)^n y^k + ... ) if the function is Positive in the whole domain and contains a finite set of minimum.

Lets take the following example:

fun=a + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d


where 0<= a,b,c,d <= 1.

Now the function fun is positive in the whole region spanned by a,b,c,d which can be confirmed through FindMinimum i.e.,

 FindMinimum[{fun, 0<=a<= 1, 0<=b<=1, 0<=c<=1, 0<=d<=1 },{a,b,c,d}]
{4.07048*10^-7, {a -> 8.09007*10^-8, b -> 0.00543937, c -> 0.99994, d -> 0.500736}}


In such scenario, it is guaranteed (at least one way) that fun can be organised as Sum of terms where each term looks like a beta integrand i.e. in the form

fun=$$\Sigma k_{i} a^{m1} (1-a)^{n1} b^{m2} (1-b)^{n2} c^{m3} (1-c)^{n3} d^{m4} (1-d)^{n4}$$

with $$m1,n1,...,m4,n4 >= 0$$ and $$k_{i} >=0$$ are real numbers . In this particular case

fun=a (1 - b)^2 + b (1 - c) (1 - a)^2 + 4 a b (1 - c) +  4 a b c (1 - d)


or also

fun=a (1 - b)^2 + b (1 - c)  + a^2 b (1 - c)  + 2 a b (1 - c) +   4 a b c (1 - d)


both are in desired form, where each term is in the form x^m(1-x)^n and importantly they are separated with +.

Any strategy how to achieve this?

Or may any functionality already available that I am not aware of! There is one functionality Factor with Modulus. However that does not work on multivariate functions.

NOTE: If the function has minimum for infinitely many points, then it is probably NOT possible to have beta-integrand representation. For example, for the simple case (a-b)^2 which has minimum(0) along the line a=b, it is NOT possible to have a representation in terms of beta-integrand.

• Aren't you missing a coefficient in front of $\text{fun}=\sum a^{m1}\dots$? Sep 15, 2020 at 11:49
• @anderstood, yes and edited. This is the important part :) . These $k_{i}$ should be either 0 or positive real numbers. Sep 15, 2020 at 11:57

## 2 Answers

One approach, maybe not the most efficient one though.

• use CoefficientArray to extract coefficients of your polynomial
• define a polynomial with unknown coefficients in the form you are looking for (beta integrand)
• use CoefficientArray on this second polynomial
• identify both; the system is linear, underdetermined, so I used FindInstance just to take one solution, but of course that's one among many.

There we go:

fun = a + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d;
coefs = Normal@CoefficientArrays[fun, {a, b, c, d}];

obj = Sum[
If[i + j + k + l + m + n + o + p > 4, 0, 1]*
alpha[i, j, k, l, m, n, o, p]*a^i (1 - a)^j*b^k*(1 - b)^l*
c^m*(1 - c)^n*d^o*(1 - d)^p, {i, 0, 2}, {j, 0, 2}, {k, 0, 2}, {l,
0, 2}, {m, 0, 2}, {n, 0, 2}, {o, 0, 2}, {p, 0, 2}];
vars = Table[
If[i + j + k + l + m + n + o + p <= 4, alpha[i, j, k, l, m, n, o, p],
Unevaluated[Sequence[]]], {i, 0, 2}, {j, 0, 2}, {k, 0, 2}, {l, 0,
2}, {m, 0, 2}, {n, 0, 2}, {o, 0, 2}, {p, 0, 2}] // Flatten;

coefs2 = Normal@CoefficientArrays[obj, {a, b, c, d}];
eqs = Thread[(Flatten@coefs - Flatten@coefs2) == 0] // DeleteDuplicates;
eqs = DeleteCases[eqs, True];
ineqs = Thread[vars >= 0];
sol = FindInstance[eqs~Join~ineqs, vars];

fun2 = obj /. sol // First
(* ((1 - a)^2 b (1 - c) + 3 a b (1 - c) + a b^2 (1 - d) +
a (1 - b) (1 - c) (1 - d) + a c (1 - d) + 2 a b c (1 - d) +
a (1 - b)^2 d + a b (1 - c) d *)

fun - fun2 // Simplify
(* 0 *)


You might instead want to transform the linear system into its matrix form with ArrayCoefficients and use LinearSolve, instead of FindInstance:

linearsystem = Normal@CoefficientArrays[eqs, vars];
sol = LinearSolve[#2, -#1] & @@ linearsystem;
fun2 /. vars -> sol
(* (1 - a)^2 b (1 - c) + 3 a b (1 - c) + a b^2 (1 - d) +
a (1 - b) (1 - c) (1 - d) + a c (1 - d) + 2 a b c (1 - d) +
a (1 - b)^2 d + a b (1 - c) d *)

fun - fun2 /. vars -> sol // Simplify
(* 0 *)

• A question: in principle your sol will have all the instances for which the system is solved for this particular form? One observation: if the expression (to begin with ) has a constant term say 1 + a b - ... this leads to two Lists of unequal length for coefs and coefs2. Sep 15, 2020 at 16:22
• @Boogeyman For the question, FindInstance finds a single instance that satisfies the equalities and inequalities, not all of them. Regarding the observation, could you give a complete example of fun that fails? I see no problem with fun = 1 + a*b + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d for example (i.e. your fun plus one). Maybe you tried with a higher order polynomial, in which case you need to adjust the 4. Sep 15, 2020 at 20:49
• For example, with fun = 1 + a*b + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d^2 (note the square at the end), you get an working solution by changing 4 to 5 (which is the higher order term ($abcd^2$). It would deserve automating. Sep 15, 2020 at 21:01
• yes indeed It was O(5) term. However for the following O(4) term it does not work... (it could not find any solution) fun=1 + 2 a b + a^2 b^2 - 2 c - 2 a c - 2 a b c - 2 a^2 b c + c^2 + 2 a c^2 + a^2 c^2 Sep 15, 2020 at 21:14
• , Now I know why this is not working. Its not the problem of your code. The problem is when I have an expression which does not have a finite number of minimum, for those cases even though the function is positive, it can NOT be written as beta integrand. Lets take the following simple positive function (a-b)^2 which has infinite set of minimum at 0 for a=b. I will edit the question. The previous example is also similar case. Sep 16, 2020 at 12:47

I included this code as a friendly version to the excellent code from @anderstood. This script handles the cases with number of variables from 1 to 4. Extension to more variables can be easily done by extending the Switch span.

Clear[ExtractFactors]
ExtractFactors[f_, vars_] := Module[{rf, ef},
rf = CoefficientRules[f, vars];
ef = Map[First, rf];
Return[ef]
]

pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d;
pol = 1 - a^2 b^2 c d^2;
pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d;
pol = 1 - a b c d;
pol = 1 - a^2 b^2 c^2 d^2;
pol = 1 - a b c;
pol = 1 - a - a^2;
pol = 1 - a b;
pol = a b + a c + c b - a b c;
pol = 1 - a b c d;
pol = 1 - a + a^2;
pol = 1 - a b c d e;
pol = 1 - a b + a^2 b^2;
pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d^2;

vars = Variables[pol];
monom = ExtractFactors[pol, vars];
numvars = Length[vars];
mexp = Table[Max[Transpose[monom][[All ;; k]]], {k, 1, numvars}];
table = Switch[Length[mexp],
1, amax = mexp[[1]]; Table[If[i + j > amax || i + j < 2, 0, 1]*a^i (1 - a)^j, {i, 0, amax}, {j, 0, amax}],
2, amax = mexp[[1]]; bmax = mexp[[2]]; Table[If[i + j > amax || k + l > bmax || i + j + k + l < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}],
3, amax = mexp[[1]]; bmax = mexp[[2]]; cmax = mexp[[3]];Table[If[i + j > amax || k + l > bmax || m + n > cmax || i + j + k + l + m + n < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l*c^m*(1 - c)^n, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}, {m, 0, cmax}, {n, 0, cmax}],
4, amax = mexp[[1]]; bmax = mexp[[2]]; cmax = mexp[[3]]; dmax = mexp[[4]]; Table[If[i + j > amax || k + l > bmax || m + n > cmax || o + p > dmax || i + j + k + l + m + n + o + p < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l*c^m*(1 - c)^n*d^o*(1 - d)^p, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}, {m, 0, cmax}, {n, 0, cmax}, {o, 0, dmax}, {p, 0, dmax}],
_, Print["Num. of variables should be less than 5"]];
If[Length[mexp] < 5,
monomials = Union[Flatten[table]];
np = Length[monomials];
A = Table[Subscript[\[Alpha], i], {i, 1, np}];
restrsA = Thread[A >= 0];
dif = A.monomials - pol;
coefs = CoefficientRules[dif, vars];
rels = Map[Last, coefs];
equs = Thread[rels == 0];
sol = Quiet@FindInstance[Join[equs, restrsA], A];
If[Length[sol] > 0,
polfound = A.monomials /. sol[[1]];
Print[polfound];
Print[Simplify[polfound - pol /. sol[[1]]]], Print["No match found"]
]
]

• thanks for the implementation. I have also adopted @anderstood method and made it a little nicer. However the cases with more than 5 variables are extremely slow (takes two days for 6-variables and degree-7 polynomial!!). Unfortunately FindInstance does not support Parallelize. However it is not clear, if I do not put explicit number of instances, whether it looks for all possible instances and then prints the first one!(that would be time-consuming). Btw, FindInstance also supports restricting the domain like FindInstance[equs,A,NonNegativeRationals] ... Saves one line ;) . Sep 24, 2020 at 9:59