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

Question

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.

Answer 1

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

Answer 2

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"]
   ]
]