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

Question

I have a matrix $M$ with non-negative real entries, and I would like to minimize the objective function $$\Phi(v) = \|Mv\|_\infty,$$ where $v$ is constrained to be a probability vector, i.e., $v_1+\dots+v_n=1$ and $v_i\ge 0$. Is there an efficient algorithm for this?

Motivation. This comes up in solving a two-player game, where the first player chooses $v$ (which represents a probability distribution on the columns of $M$), then the second player sees $v$ and chooses a row $i$, and the first player loses $(Mv)_i$ and the second player gains $(Mv)_i$. I want to find an efficient algorithm to find the optimal strategy for the first player.

Best I can do. I can see how to solve this with linear programming. We introduce a variable $\ell$ to represent the final loss of the first player, then we add the linear inequalities $M_i v \le \ell$ (and the constraints on $v$) and we minimize $\ell$. Is there a more efficient algorithm than this?

Answer 1

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_i\ge0 \\ & \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_i\ge0 \\ & \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 $x\le y$ and $-x\le 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_1\le y, x_2\le y, \ldots, x_n \le y$.

We can therefore rewrite the problem as $$ \begin{align} \min_v \quad & y \\ s.t. \quad & v_i\ge0 \\ & \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!

Answer 2

You can and should solve this problem without linear programming and apply the Bellman equation instead.

Actually, the minmax theorem -- handled numerically via LP -- is only required to solve the problem where both players simultaneously choose an action.

In contrast, your game consists of a two-step process, and the mathematical model should incorporate this structure. This can be realized by a Markov decision process that is optimized via the Bellman equation. Basically, you there solve two "max" problems instead of one "minmax" problem, which is way easier from both the mathematical and the computational perspective.