Source code for quantecon.optimize.minmax

"""
Contain a minmax problem solver routine.

"""
import numpy as np
from numba import jit
from .linprog_simplex import solve_tableau, PivOptions
from .pivoting import _pivoting


[docs]@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 tableau[-1, n] = -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 # 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