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$
  • 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$ – Henrik Schumacher Feb 5 '18 at 11:24
  • 1
    $\begingroup$ Have also a look at the option Method -> "InteriorPoint". This should speed up things considerably. $\endgroup$ – Henrik Schumacher Feb 5 '18 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$ – Héctor Feb 5 '18 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$ – Henrik Schumacher Feb 5 '18 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$ – Héctor Feb 5 '18 at 13:34
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$
  • $\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$ – Héctor Feb 10 '18 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@(N@Aineq.vars - 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@(N@Aineq.vars - 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$
  • $\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$ – Héctor Feb 6 '18 at 9:40
  • 1
    $\begingroup$ @Héctor, there was an error in the code, but I fixed it in the meantime. $\endgroup$ – Henrik Schumacher Feb 6 '18 at 10:21
  • 1
    $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Feb 6 '18 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$ – Henrik Schumacher Feb 6 '18 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$ – Henrik Schumacher Feb 6 '18 at 11:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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