Solve two-player game - minimize the l-infinity norm of a matrix-vector product

We can do some transformations to your problem to show that it's easily solvable via a linear program:

Given a matrix $M$ with non-negative real entries and a vector $v$ you wish to solve the problem: $$ begin{align} min_v quad & lVert Mv rVert_infty s.t. quad & v_ige0 & sum_i v_i = 1 end{align} $$ Now, note that $lVert Mv rVert_infty = max_i |(Mv)_i|$. With this in hand we get to the problem: $$ begin{align} min_v quad & max_i |(Mv)_i| s.t. quad & v_ige0 & sum_i v_i = 1 end{align} $$ We can treat an absolute value $|x|$ in a min-LP by replacing $|x|$ with a variable $y$ and adding the constraints $xle y$ and $-xle y$.

We can replace a maximum function in a min-LP by replacing $max_i (x)_i$ with a variable $y$ and adding the constraints $x_1le y, x_2le y, ldots, x_n le y$.

We can therefore rewrite the problem as $$ begin{align} min_v quad & y s.t. quad & v_ige0 & sum_i v_i = 1 & (M_{1,*} v) le y & -(M_{1,*} v) le y & (M_{i,*} v) le y & -(M_{i,*} v) le y end{align} $$ Where $M_{i,*}$ is the $i$-th row of the matrix $M$.

Since this is a convex problem, you can solve it using cvxpy like so:

import cvxpy as cp
import numpy as np

M = np.random.rand(10,10)
v = cp.Variable(10)

objective = cp.Minimize(cp.norm(M*v, 'inf'))
constraints = [sum(v)==1, v>=0]
problem = cp.Problem(objective, constraints)
objval = problem.solve()

print("Objective value = ", objval)
print("v values = ", v.value)

Notice that CVXPY has automagically performed all of the transformations we used above.

Now, efficiency. We can judge this by multiple metrics.

1. Your time. It's hopefully obvious that CVXPY and similar tools offer an extremely efficient and flexible way of solving this kind of problem. If you need to add constraints later you can do so quickly.
2. Compute time. LP solvers are often highly optimized. You should expect them to work quickly even with large numbers of variables and constraints.

Let's look at this second point by timing the above:

import cvxpy as cp
import numpy as np
import timeit

M = np.random.rand(1000,1000)
v = cp.Variable(1000)

objective = cp.Minimize(cp.norm(M*v, 'inf'))
constraints = [sum(v)==1, v>=0]
problem = cp.Problem(objective, constraints)
timeit.timeit(lambda: problem.solve(), number=4)

This gives:

     Size | Time
  10x10   | 0.39s
 100x100  | 3.37s
1000x1000 | 345s

A lot of this is Python overhead. If we instead use Julia, we get much better timing:

using Convex
using ECOS

M = rand(10,10);
v = Variable(10);
problem=minimize(norm_inf(M*v), [v>=0, sum(v)==1])
@time solve!(problem, ECOS.Optimizer)

Timing results:

     Size | Time
  10x10   | 0.0033s
 100x100  | 0.07s
1000x1000 | 96s

Much better! Note that we're also using the ECOS solver. Other options, especially commercial ones, might be much faster.

I'm skeptical that other approaches would improve much on the times for smaller (10x10, 100x100) problems, or that you'd be able to make meaningful use of those improvements (outside of some HPC context).

Dynamic programming, as another answer suggested, might also be tricky to implement here. DP alone is slow because the game tree expands exponentially for each additional level of recursion. You make DP fast by memoizing states, but that's impractical if your states are continuous (your problem) or don't overlap (chess, Go).

EDIT:

Brian Borchers comments:

Note that since M has nonnegative entries and v≥0, you don't actually need to handle the absolute values

I'd avoided making use of this information initially in order to provide a fully general answer, but if we do leverage it in Julia:

using Convex
using ECOS

N = 1000
M = rand(N,N);
v = Variable(N);
problem=minimize(maximum(M*v), [v>=0, sum(v)==1])
@time solve!(problem, ECOS.Optimizer)

With this simplification of the constraints, the 1000x1000 problem takes only 19s!