Programming a MaxMin Linear Optimization

Question

I want to program a function for a two-player game. Basically it's like this: . Each player has an array of options, and the result of the game is based on both players choices. So heres my most promising approach (i've also played around with LinearProgramming, which also works up to a point). c is like Length[x]. p[i,j] is a reference to my database, and gives a value between 0 and 1. It describes the probability of a win for player x if he chooses strategy i and player y chooses strategy j.

p := {{0.5, 0.7, 0.3}, {0.3, 0.5, 0.6}, {0.7, 0.4, 0.5}}

stratmin[c_, x__] := 
  NMinimize[{Sum[x[[i]]*y[j]*p[i, j], {i, c}, {j, c}], 
   Table[y[j] >= 0, {j, c}], Sum[y[j], {j, c}] == 1}, 
   Table[y[j], {j, c}]]

strat[c_] := 
 NMaximize[{stratmin[c, Table[x[i], {i, c}]], 
  Table[x[i] >= 0, {i, c}], Sum[x[i], {i, c}] == 1}, 
  Table[x[i], {i, c}]]

strat[3]

But it doesn't work! As you could have guessed, as I am asking a question here. NMinimize doesn't seem to use the x[i] as a parameter that NMaximize uses. If I program this with LinearProgramming, since I've got to "nest" variables in the same way, I've got the same problem. The solution btw. should look smth. like strat=0.5; x[1] ca. 0.5, x[2] ca. 0.4, x[3] ca. 0.6.

EDIT:

I looked into your program, and tried to understand and research all parts; had to read through this and this to understand all of the notation; but still dont really get how exactly the solution worked. I've reduced the problem to:

testmin[x_] := First[NMinimize[{-x^2 + y, y >= 0}, y]]

testmax := NMaximize[{testmin[x]}, x] 

testmax 

I don't know whether that is too much too ask, but would you do the minimal necessary correction for that example, too, please? Hopefully that helps me understand the necessary step. Also, if you know of any other good resources or tutorials, I would be thankful.

Answer 1

It takes some careful coding to make sure the right values are explicitly numeric at the time they need to be (in the inner optimization). Can be done as below. And there may be better ways, I'm no expert.

stratmin[p_ /; MatrixQ[p, Element[#, Reals] &], 
  xlist_List /; VectorQ[xlist, Element[#, Reals] &]] := Module[
   {y, c = Length[p], yvars, min, yvals},
   yvars = Array[y, c]; {min, yvals} = 
    NMinimize[{xlist.p.yvars, Map[# >= 0 &, yvars], 
      Total[yvars] == 1}, yvars];
   {min, yvars /. yvals}
   ] /; Length[p] === Length[p[[1]]] === Length[xlist]


strat[p_ /; MatrixQ[p, Element[#, Reals] &]] := Module[
   {x, c = Length[p], xvars, max, xvals},
   xvars = Array[x, c];
   obj[xl_ /; VectorQ[xl, Element[#, Reals] &]] := 
    First[stratmin[p, xl]];
   {max, xvals} = 
    NMaximize[{obj[xvars], Map[# >= 0 &, xvars], Total[xvars] == 1}, 
     xvars];
   {max, xvars /. xvals, stratmin[p, xvars /. xvals][[2]]}
   ] /; Length[p] === Length[p[[1]]]

Example:

SeedRandom[1234567];
c = 4;
p = RandomReal[1, {c, c}];

strat[p]

(* Out[203]= {0.457170468845, {0.823855940951, 0., 0.176144059049, 
  0.}, {1., 0., 0., 0.}} *)

This takes a couple of seconds. At c=8 it is already 14 seconds. So it might not scale as well as you want. Again, there may be better ways to go about this. One thing I note is that I have it coded in a way that manages to miss the fact that each level is a linear program (a message at c=16 makes that clear). This is almost certainly contributing to slowness and may well also cause suboptimal results. So this should just be viewed as a starting point for further refinement.

--- edit ---

Actually I do not see how to recast the outer level as an LP. And NMaximize seems to agree that this is not what it is. For any choice of x values the inner optimization is an LP (and is done a bit faster using FindMinimum instead of NMinimize). The fact that the optimal y values are always vertices of the c-simplex is just due to this LP nature of the inner problems.

That said, if I am seeing this correctly, you can just iterate over the c vertices of y vectors, and maximize for x.

--- end edit ---

--- edit #2 ---

I was not seeing that last part correctly. Player x has to use a probability vector in such a way that player y cannot make the expected value lower than the maximin. Knowing that a vertex vector will be chosen is not enough for x to resort to the idea put forward above.

--- end edit #2 ---