4
$\begingroup$

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.

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

2 Answers 2

3
+50
$\begingroup$

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 *)
$\endgroup$
11
  • $\begingroup$ 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. $\endgroup$
    – BabaYaga
    Sep 15, 2020 at 16:22
  • 1
    $\begingroup$ @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. $\endgroup$
    – anderstood
    Sep 15, 2020 at 20:49
  • 1
    $\begingroup$ 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. $\endgroup$
    – anderstood
    Sep 15, 2020 at 21:01
  • $\begingroup$ 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 $\endgroup$
    – BabaYaga
    Sep 15, 2020 at 21:14
  • 1
    $\begingroup$ , 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. $\endgroup$
    – BabaYaga
    Sep 16, 2020 at 12:47
1
$\begingroup$

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"]
   ]
]
$\endgroup$
1
  • 1
    $\begingroup$ 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 ;) . $\endgroup$
    – BabaYaga
    Sep 24, 2020 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.