3
$\begingroup$

I need a (more or less consistent) way of computing a vector in Mathematica given a set of inequalities that define non-empty convex polytope. Because the question has a strong applied background (Game Theory), let me describe my problem intuitively by means of a simple example.

First, the preliminaries. Consider a function $v:2^N\longrightarrow \mathbb{R}$, where $N=\{1,2,3\}$ and $2^N$ denotes the Power Set of $N$. Then, consider the following set of inequalities: $ x_i \geqslant v(\{i\}) \forall i\in N$; $x_1+x_2+x_3=v(\{1,2,3\})$; $x_1+x_2\geqslant v(\{1,2\})$; $x_1+x_3\geqslant v(\{1,3\})$ and $x_2+x_3\geqslant v(\{2,3\})$. Let us suppose that the set of vectors $(x_1,x_2,x_3)$ that solve the previous set of inequalities is non-empty. If so, this set is known to be a convex polytope whose vertices can be derived from the previous set of inequalities.

To see so, consider now the following example: $v(\{1,2,3\})=1$, $v(\{1,2\})=0.75$, $v(\{1,3\})=0.50$, $v(\{2,3\})=0.25$ and $v(\{i\})=0\forall i \in \{1,2,3\}$. The following code solves for the corresponding set of inequalities, thus providing the set of vertices of the polytope given by $v$.

Reduce[{x + y + z == 1, x + y >= 0.75, x + z >= 0.5, y + z >= 0.25, x >= 0, y >= 0, z >= 0}, {x, y, z}, Reals]

From there, it can be readily verified that the vertices of my polytope (a polygon, in this particular case) are given by $(0.75,0.25,0),(0.75,0,0.25), (0.25,0.50,0.25)\text{ and }(0.50,0.50,0)$.

Now, the problem. Consider all orders over the set $N$. Now, consider the lexicographic maximum of each of these $|N|!$ orders with respect to the set of vertices. Let $\pi_i\in\Pi$ denote each order and $\lambda_{\pi_{i}}$ denote the lexicographic maximum for each order $\pi_i$. In other words, in our particular example: (i) for order $\pi_1=\{1,2,3\}$, $\lambda_{\pi_{1}}=(0.75,0.25,0)$; (ii) for orders $\pi_2=\{1,3,2\}$ and $\pi_6=\{3,1,2\}$, $\lambda_{\pi_{2}}=\lambda_{\pi_{6}}=(0.75,0,0.25)$; (iii) for order $\pi_3=\{2,1,3\}$, $\lambda_{\pi_{3}}=(0.50,0.50,0)$ and (iv) for orders $\pi_4=\{2,3,1\}$ and $\pi_5=\{3,2,1\}$, $\lambda_{\pi_{4}}=\lambda_{\pi_{5}}=(0,25,0.50,0.25)$.

Averaging all lexicographic maximums over all possible orders, we get $ \mathfrak{A}=\frac{1}{|N|!}\cdot\sum_{\pi_{i}\in\Pi} \lambda_{\pi_{i}}$. In our case, $\mathfrak{A}=\left(\frac{13}{24},\frac{7}{24},\frac{4}{24}\right)$. Then, I have two questions:

  1. Can anyone provide code in Mathematica that computes this vector $\mathfrak{A}$ for my particular example?
  2. Can anyone provide code in Mathematica that computes this vector $\mathfrak{A}$ given a set $N$ such that $|N|=4$ and a generic function $v:N\longrightarrow \mathbb{R}$?

Thank you all very much in advance for your time.

$\endgroup$
5
  • 1
    $\begingroup$ If $\nu$ is affine-linear then LinearProgramming might help you (see also this introductory site. It finds the minimum or (maximum) of a linear function restricted to a convex polytope, defined by affine-linear equalities. Using an arbitrary linear function, you can employ it to find feasible points, i.e., points within this polytope. $\endgroup$ Feb 5, 2018 at 11:24
  • 1
    $\begingroup$ Have also a look at the option Method -> "InteriorPoint". This should speed up things considerably. $\endgroup$ Feb 5, 2018 at 11:26
  • $\begingroup$ Thank you for your comments. I'm a total noob in Mathematica and I was not aware of such a resource. I'll dig into it later to see if I can get what I need from there. $\endgroup$
    – EoDmnFOr3q
    Feb 5, 2018 at 11:27
  • 1
    $\begingroup$ You're welcome! Let me know if you need further advice (write a comment that includes @HenrikSchumacher to ping me). $\endgroup$ Feb 5, 2018 at 11:28
  • $\begingroup$ @HenrikSchumacher: I've spent these two hours trying stuff and reading in the introductory site and I certainly don't see how "LinearProgramming" could help me. Could you please provide any further hint? $\endgroup$
    – EoDmnFOr3q
    Feb 5, 2018 at 13:34

2 Answers 2

2
$\begingroup$

If you have the defining inequalities then the idea is the same regardless of dimension. Take all orders, find successive maxima for a given order, then average. The code below has some redundancy, so is not as efficient as it might be.

opt[order_, constraints_, vars_] := Module[
  {cons = constraints, max, vals, var},
  Do[var = vars[[order[[j]]]]; {max, vals} = 
    Maximize[{var, cons}, vars]; 
   cons = Join[cons, {var == (var /. vals)}];
   , {j, Length[vars]}];
  vars /. vals
  ]

The example:

vars = Array[x, 3];
constraints = {x + y + z == 1, x + y >= 0.75, x + z >= 0.5, 
    y + z >= 0.25, x >= 0, y >= 0, z >= 0} /. 
   Thread[{x, y, z} -> vars];
orders = Permutations[Range[3]];

Mean[Table[opt[order, constraints, vars], {order, orders}]]

(* Out[654]= {0.541666666667, 0.291666666667, 0.166666666667} *)
$\endgroup$
1
  • $\begingroup$ I have computed manually the vector I should get for a couple of examples with $|N|=4$ and your code delivers the desired vector. Hence, I have accepted your answer. Thank you very much for your time. $\endgroup$
    – EoDmnFOr3q
    Feb 10, 2018 at 11:22
4
$\begingroup$

You may use LinearProgramming in order to find a point within a convex polytope defined by affine-linear equalities and inequalities. To this end, the conditions have to be encoded into a matrix m and a vector b such that the conditions are equivalent to m.vars >= b and var>=0; both inequalities are meant elementwise. This is often called standard form (beware that different authors use different conventions) and how every linear program can be transformed into standard form can be found, e.g., here.

The following code will perform that for you; we start with your example problem:

cond = {x + y + z == 1, x + y >= 0.75, x + z >= 0.5, y + z >= 0.25, 
   x >= 0, y >= 0, z >= 0};
vars = {x, y, z};

The following transforms these conditions into standard form:

cond1 = cond /. {
Equal[a_, b_] :> Equal[a - b, 0],
LessEqual[a_, b_] :> GreaterEqual[b - a, 0],
GreaterEqual[a_, b_] :> GreaterEqual[a - b, 0]
};
eq = Cases[cond1, _Equal];
ineq = Cases[cond1, _GreaterEqual];
c = ConstantArray[0., 2 Length[vars]];
With[{ToPack = Developer`ToPackedArray},
  Aeq = ToPack@N@D[eq[[All, 1]], {vars, 1}];
  beq = ToPack@N@(Aeq.vars - eq[[All, 1]]);
  Aineq = ToPack@N@D[ineq[[All, 1]], {vars, 1}];
  bineq = ToPack@N@([email protected] - ineq[[All, 1]]);
  ];
m = ArrayFlatten[{
    {Aeq, -Aeq},
    {-Aeq, Aeq},
    {Aineq, -Aineq}
    }];
b = Join[beq, -beq, bineq];

Now we are ready to fead it to LinearProgramming; it will look for a minimizer of the linear function Dot[c,vars] on the described polytope. Hence (c being the zero vector), it will just return any point in the polytope. With standard options, this will be always a corner of the polytope.

sol = With[{sol0 = LinearProgramming[c, m, b]},
 sol0[[1 ;; Length[vars]]] - sol0[[Length[vars] + 1 ;; 2 Length[vars]]]
 ]

{0.75, 0., 0.25}

Checking if the solution is feasible:

Show[
 DiscretizeRegion@ImplicitRegion[cond, Evaluate@vars],
 Graphics3D[Sphere[sol, 0.01]]
 ]

enter image description here

Aha, seems to be okay.

The number of variables and conditions is quite irrelevant; LinearProgramming is built for medium to large scale application.

Edit

If you need all corners of the polytope then the following should work: The idea is that every corner has to satisfy all linear equality and a suitable number if (linear independent) inequalities is active. A brute force method might be

p = DeleteDuplicates@Quiet@DeleteCases[
     Table[
      LinearSolve[Join[Aeq, Aineq[[idx]]], Join[beq, bineq[[idx]]]],
      {idx, 
       Subsets[Range[Length[Aineq]], {Length[vars] - Length[Aeq]}]}],
     _LinearSolve];
pts = Pick[p,
   Total[
    Subtract[1, 
     UnitStep[
      Subtract[Transpose[Aineq.Transpose[p]], 
       ConstantArray[bineq, Length[p]]]]],
    {2}
    ],
   0
   ];

Here, p is a list of all points that have suitably many active inequalities and from these, we sieve out the feasible ones (this is the list pts). Quick check:

Show[
 DiscretizeRegion@ImplicitRegion[cond, Evaluate@vars],
 Graphics3D[Sphere[pts, 0.01]]
 ]

enter image description here

Further edit

I think the following would solve your entire problem. At least for nn=3 and ν as given above, it returns your result. However for nn=4, I was not able to find a function ν that does not produce empty polytopes; but that's certainly something you can fix yourself.

nn = 4;
bigN = Range[nn];
vars = Array[x, nn];
powerset = Subsets[bigN, {1, nn}];
ν = AssociationThread[powerset,RandomReal[{0,1},Length[powerset]]];
(* ν = AssociationThread[powerset, {0., 0., 0., 0.75, 0.5, 0.25, 1.}]; *)
cond = Append[
   Table[Total[vars[[subset]]] >= ν[subset], {subset, 
     Most@powerset}],
   Total[vars] == ν[bigN]
   ];

cond1 = cond /. {
    Equal[a_, b_] :> Equal[a - b, 0],
    LessEqual[a_, b_] :> GreaterEqual[b - a, 0],
    GreaterEqual[a_, b_] :> GreaterEqual[a - b, 0]
    };
eq = Cases[cond1, _Equal];
ineq = Cases[cond1, _GreaterEqual];
With[{ToPack = Developer`ToPackedArray},
  Aeq = ToPack@N@D[eq[[All, 1]], {vars, 1}];
  beq = ToPack@N@(Aeq.vars - eq[[All, 1]]);
  Aineq = ToPack@N@D[ineq[[All, 1]], {vars, 1}];
  bineq = ToPack@N@([email protected] - ineq[[All, 1]]);
  ];

p = DeleteDuplicates@
   Quiet@DeleteCases[
     Table[LinearSolve[Join[Aeq, Aineq[[idx]]], 
       Join[beq, bineq[[idx]]]], {idx, 
       Subsets[Range[
         Length[Aineq]], {Length[vars] - 
          Length[Aeq]}]}], _LinearSolve];
pts = Pick[p, 
   Total[Subtract[1, 
     UnitStep[
      Subtract[Transpose[Aineq.Transpose[p]], 
       ConstantArray[bineq, Length[p]]]]], {2}], 0];

If[Length[pts] == 0,
 0,
 1/nn! Sum[
   pts[[Ordering[pts[[All, π]], -1][[1]]]], {π, 
    Permutations[Range[nn]]}]
 ]
$\endgroup$
12
  • $\begingroup$ Thank you very much for your Edit. Seeing your extended answer, let me be more specific. I do not need all vertices of the polytope. Instead, I need all the vertices that are lexicographically maximum for at least one order. Your original answer finds one of these lexicographical maximum vertices; the extended answer solves the Vertex Enumeration Problem. What I need is is right in between. Is my problem better explained now? I am sorry for any misunderstanding my previous comment may have generated. $\endgroup$
    – EoDmnFOr3q
    Feb 6, 2018 at 9:40
  • 1
    $\begingroup$ @Héctor, there was an error in the code, but I fixed it in the meantime. $\endgroup$ Feb 6, 2018 at 10:21
  • 1
    $\begingroup$ You're welcome! $\endgroup$ Feb 6, 2018 at 10:39
  • 1
    $\begingroup$ Sure. RandomReal in produces a list of ... erm... random reals. Replace that with a list of your choice. Then either remove all occurences of N@ or apply Rationalize to the final result. $\endgroup$ Feb 6, 2018 at 11:13
  • 1
    $\begingroup$ You just have to repace RandomReal[{0,1},Length[powerset]] by a list of numbers of the same length (the length is Length[powerset]). nn=4 tells the system that you want dimension 4. $\endgroup$ Feb 6, 2018 at 11:47

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.