EHN: Add minmax solver (#579)

* EHN: Add minmax solver

* MAINT: Player: Use minmax in is_dominated and dominated_actions

Co-authored-by: mmcky <[email protected]>

Author: oyamad

docs/source/optimize.rst

@@ -5,6 +5,7 @@ Optimize
    :maxdepth: 2
 
    optimize/linprog_simplex
+   optimize/minmax
    optimize/nelder_mead
    optimize/pivoting
    optimize/root_finding

docs/source/optimize/minmax.rst

@@ -0,0 +1,7 @@
+minmax
+======
+
+.. automodule:: quantecon.optimize.minmax
+    :members:
+    :undoc-members:
+    :show-inheritance:

quantecon/game_theory/normal_form_game.py

@@ -434,11 +434,10 @@ def is_dominated(self, action, tol=None, method=None):
             default to the value of the `tol` attribute.
 
         method : str, optional(default=None)
-            If None, `lemke_howson` from `quantecon.game_theory` is used
-            to solve for a Nash equilibrium of an auxiliary zero-sum
-            game. If `method` is set to `'simplex'`, `'interior-point'`,
-            or `'revised simplex'`, then `scipy.optimize.linprog` is
-            used with the method as specified by `method`.
+            If None, `minmax` from `quantecon.optimize` is used. If
+            `method` is set to `'simplex'`, `'interior-point'`, or
+            `'revised simplex'`, then `scipy.optimize.linprog` is used
+            with the method as specified by `method`.
 
         Returns
         -------
@@ -465,10 +464,9 @@ def is_dominated(self, action, tol=None, method=None):
             D.shape = (D.shape[0], np.prod(D.shape[1:]))
 
         if method is None:
-            from .lemke_howson import lemke_howson
-            g_zero_sum = NormalFormGame([Player(D), Player(-D.T)])
-            NE = lemke_howson(g_zero_sum)
-            return NE[0] @ D @ NE[1] > tol
+            from ..optimize.minmax import minmax
+            v, _, _ = minmax(D)
+            return v > tol
         elif method in ['simplex', 'interior-point', 'revised simplex']:
             from scipy.optimize import linprog
             m, n = D.shape
@@ -506,11 +504,10 @@ def dominated_actions(self, tol=None, method=None):
             default to the value of the `tol` attribute.
 
         method : str, optional(default=None)
-            If None, `lemke_howson` from `quantecon.game_theory` is used
-            to solve for a Nash equilibrium of an auxiliary zero-sum
-            game. If `method` is set to `'simplex'`, `'interior-point'`,
-            or `'revised simplex'`, then `scipy.optimize.linprog` is
-            used with the method as specified by `method`.
+            If None, `minmax` from `quantecon.optimize` is used. If
+            `method` is set to `'simplex'`, `'interior-point'`, or
+            `'revised simplex'`, then `scipy.optimize.linprog` is used
+            with the method as specified by `method`.
 
         Returns
         -------

quantecon/optimize/__init__.py

@@ -5,6 +5,7 @@
 from .linprog_simplex import (
     linprog_simplex, solve_tableau, get_solution, PivOptions
 )
+from .minmax import minmax
 from .scalar_maximization import brent_max
 from .nelder_mead import nelder_mead
 from .root_finding import newton, newton_halley, newton_secant, bisect, brentq

quantecon/optimize/minmax.py

@@ -0,0 +1,113 @@
+"""
+Contain a minmax problem solver routine.
+
+"""
+import numpy as np
+from numba import jit
+from .linprog_simplex import _set_criterion_row, solve_tableau, PivOptions
+from .pivoting import _pivoting
+
+
+@jit(nopython=True, cache=True)
+def minmax(A, max_iter=10**6, piv_options=PivOptions()):
+    r"""
+    Given an m x n matrix `A`, return the value :math:`v^*` of the
+    minmax problem:
+
+    .. math::
+
+        v^* = \max_{x \in \Delta_m} \min_{y \in \Delta_n} x^T A y
+            = \min_{y \in \Delta_n}\max_{x \in \Delta_m} x^T A y
+
+    and the optimal solutions :math:`x^* \in \Delta_m` and
+    :math:`y^* \in \Delta_n`: :math:`v^* = x^{*T} A y^*`, where
+    :math:`\Delta_k = \{z \in \mathbb{R}^k_+ \mid z_1 + \cdots + z_k =
+    1\}`, :math:`k = m, n`.
+
+    This routine is jit-compiled by Numba, using
+    `optimize.linprog_simplex` routines.
+
+    Parameters
+    ----------
+    A : ndarray(float, ndim=2)
+        ndarray of shape (m, n).
+
+    max_iter : int, optional(default=10**6)
+        Maximum number of iteration in the linear programming solver.
+
+    piv_options : PivOptions, optional
+        PivOptions namedtuple to set tolerance values used in the linear
+        programming solver.
+
+    Returns
+    -------
+    v : float
+        Value :math:`v^*` of the minmax problem.
+
+    x : ndarray(float, ndim=1)
+        Optimal solution :math:`x^*`, of shape (,m).
+
+    y : ndarray(float, ndim=1)
+        Optimal solution :math:`y^*`, of shape (,n).
+
+    """
+    m, n = A.shape
+
+    min_ = A.min()
+    const = 0.
+    if min_ <= 0:
+        const = min_ * (-1) + 1
+
+    tableau = np.zeros((m+2, n+1+m+1))
+
+    for i in range(m):
+        for j in range(n):
+            tableau[i, j] = A[i, j] + const
+        tableau[i, n] = -1
+        tableau[i, n+1+i] = 1
+
+    tableau[-2, :n] = 1
+    tableau[-2, -1] = 1
+
+    # Phase 1
+    pivcol = 0
+
+    pivrow = 0
+    max_ = tableau[0, pivcol]
+    for i in range(1, m):
+        if tableau[i, pivcol] > max_:
+            pivrow = i
+            max_ = tableau[i, pivcol]
+
+    _pivoting(tableau, n, pivrow)
+    _pivoting(tableau, pivcol, m)
+
+    basis = np.arange(n+1, n+1+m+1)
+    basis[pivrow] = n
+    basis[-1] = 0
+
+    # Modify the criterion row for Phase 2
+    c = np.zeros(n+1)
+    c[-1] = -1
+    _set_criterion_row(c, basis, tableau)
+
+    # Phase 2
+    solve_tableau(tableau, basis, max_iter-2, skip_aux=False,
+                  piv_options=piv_options)
+
+    # Obtain solution
+    x = np.empty(m)
+    y = np.zeros(n)
+
+    for i in range(m+1):
+        if basis[i] < n:
+            y[basis[i]] = tableau[i, -1]
+
+    for j in range(m):
+        x[j] = tableau[-1, n+1+j]
+        if x[j] != 0:
+            x[j] *= -1
+
+    v = tableau[-1, -1] - const
+
+    return v, x, y

quantecon/optimize/tests/test_minmax.py

@@ -0,0 +1,44 @@
+"""
+Tests for minmax
+
+"""
+import numpy as np
+from numpy.testing import assert_, assert_allclose
+
+from quantecon.optimize import minmax
+from quantecon.game_theory import NormalFormGame, Player, lemke_howson
+
+
+class TestMinmax:
+    def test_RPS(self):
+        A = np.array(
+            [[0, 1, -1],
+             [-1, 0, 1],
+             [1, -1, 0],
+             [-1, -1, -1]]
+        )
+        v_expected = 0
+        x_expected = [1/3, 1/3, 1/3, 0]
+        y_expected = [1/3, 1/3, 1/3]
+
+        v, x, y = minmax(A)
+
+        assert_allclose(v, v_expected)
+        assert_allclose(x, x_expected)
+        assert_allclose(y, y_expected)
+
+    def test_random_matrix(self):
+        seed = 12345
+        rng = np.random.default_rng(seed)
+        size = (10, 15)
+        A = rng.normal(size=size)
+        v, x, y = minmax(A)
+
+        for z in [x, y]:
+            assert_((z >= 0).all())
+            assert_allclose(z.sum(), 1)
+
+        g = NormalFormGame((Player(A), Player(-A.T)))
+        NE = lemke_howson(g)
+        assert_allclose(v, NE[0] @ A @ NE[1])
+        assert_(g.is_nash((x, y)))