FIX: support_enumeration: Use _numba_linalg_solve (#311)

* support_enumeration: Remove fallback for Numba < 0.28

* support_enumeration: Add a test

"LinAlgError: Matrix is singular to machine precision.” raised

* FIX: support_enumeration: Use `_numba_linalg_solve`

Remove `is_singular` by svd

* util: Add `_numba_linalg_solve`

For use in a jitted function in nopython mode
* Call directly Numba internal `numba_xgesv`
* Return nonzero int if input matrix is singular, allowing alternative to try-except np.linalg.LinAlgError

* support_enumeration: Remove `any()`

Allow `cache=True`, close #285

Author: oyamad

docs/source/util.rst

@@ -7,5 +7,6 @@ Utilities
    util/array
    util/common_messages
    util/notebooks
+   util/numba
    util/random
    util/timing

docs/source/util/numba.rst

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

quantecon/game_theory/support_enumeration.py

@@ -11,19 +11,9 @@
 Tardos, and V. Vazirani eds., Algorithmic Game Theory, 2007.
 
 """
-from distutils.version import LooseVersion
 import numpy as np
-import numba
 from numba import jit
-
-
-least_numba_version = LooseVersion('0.28')
-is_numba_required_installed = True
-if LooseVersion(numba.__version__) < least_numba_version:
-    is_numba_required_installed = False
-nopython = is_numba_required_installed
-
-EPS = np.finfo(float).eps
+from ..util.numba import _numba_linalg_solve
 
 
 def support_enumeration(g):
@@ -46,11 +36,6 @@ def support_enumeration(g):
     list(tuple(ndarray(float, ndim=1)))
         List containing tuples of Nash equilibrium mixed actions.
 
-    Notes
-    -----
-    This routine is jit-complied if Numba version 0.28 or above is
-    installed.
-
     """
     return list(support_enumeration_gen(g))
 
@@ -80,7 +65,7 @@ def support_enumeration_gen(g):
                                     g.players[1].payoff_array)
 
 
-@jit(nopython=nopython)  # cache=True raises _pickle.PicklingError
+@jit(nopython=True)  # cache=True raises _pickle.PicklingError
 def _support_enumeration_gen(payoff_matrix0, payoff_matrix1):
     """
     Main body of `support_enumeration_gen`.
@@ -105,32 +90,28 @@ def _support_enumeration_gen(payoff_matrix0, payoff_matrix1):
 
     for k in range(1, n_min+1):
         supps = (np.arange(k), np.empty(k, np.int_))
-        actions = (np.empty(k), np.empty(k))
+        actions = (np.empty(k+1), np.empty(k+1))
         A = np.empty((k+1, k+1))
-        A[:-1, -1] = -1
-        A[-1, :-1] = 1
-        A[-1, -1] = 0
-        b = np.zeros(k+1)
-        b[-1] = 1
+
         while supps[0][-1] < nums_actions[0]:
             supps[1][:] = np.arange(k)
             while supps[1][-1] < nums_actions[1]:
                 if _indiff_mixed_action(payoff_matrix0, supps[0], supps[1],
-                                        A, b, actions[1]):
+                                        A, actions[1]):
                     if _indiff_mixed_action(payoff_matrix1, supps[1], supps[0],
-                                            A, b, actions[0]):
+                                            A, actions[0]):
                         out = (np.zeros(nums_actions[0]),
                                np.zeros(nums_actions[1]))
                         for p, (supp, action) in enumerate(zip(supps,
                                                                actions)):
-                            out[p][supp] = action
+                            out[p][supp] = action[:-1]
                         yield out
                 _next_k_array(supps[1])
             _next_k_array(supps[0])
 
 
-@jit(nopython=nopython)
-def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out):
+@jit(nopython=True, cache=True)
+def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, out):
     """
     Given a player's payoff matrix `payoff_matrix`, an array `own_supp`
     of this player's actions, and an array `opp_supp` of the opponent's
@@ -139,8 +120,7 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out):
     among the actions in `own_supp`, if any such exists. Return `True`
     if such a mixed action exists and actions in `own_supp` are indeed
     best responses to it, in which case the outcome is stored in `out`;
-    `False` otherwise. Arrays `A` and `b` are used in intermediate
-    steps.
+    `False` otherwise. Array `A` is used in intermediate steps.
 
     Parameters
     ----------
@@ -154,17 +134,11 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out):
         Array containing the opponent's action indices, of length k.
 
     A : ndarray(float, ndim=2)
-        Array used in intermediate steps, of shape (k+1, k+1). The
-        following values must be assigned in advance: `A[:-1, -1] = -1`,
-        `A[-1, :-1] = 1`, and `A[-1, -1] = 0`.
-
-    b : ndarray(float, ndim=1)
-        Array used in intermediate steps, of shape (k+1,). The following
-        values must be assigned in advance `b[:-1] = 0` and `b[-1] = 1`.
+        Array used in intermediate steps, of shape (k+1, k+1).
 
     out : ndarray(float, ndim=1)
-        Array of length k to store the k nonzero values of the desired
-        mixed action.
+        Array of length k+1 to store the k nonzero values of the desired
+        mixed action in `out[:-1]` (and the payoff value in `out[-1]`.)
 
     Returns
     -------
@@ -175,15 +149,22 @@ def _indiff_mixed_action(payoff_matrix, own_supp, opp_supp, A, b, out):
     m = payoff_matrix.shape[0]
     k = len(own_supp)
 
-    A[:-1, :-1] = payoff_matrix[own_supp, :][:, opp_supp]
-    if _is_singular(A):
-        return False
-
-    sol = np.linalg.solve(A, b)
-    if (sol[:-1] <= 0).any():
+    for i in range(k):
+        for j in range(k):
+            A[j, i] = payoff_matrix[own_supp[i], opp_supp[j]]  # transpose
+    A[:-1, -1] = 1
+    A[-1, :-1] = -1
+    A[-1, -1] = 0
+    out[:-1] = 0
+    out[-1] = 1
+
+    r = _numba_linalg_solve(A, out)
+    if r != 0:  # A: singular
         return False
-    out[:] = sol[:-1]
-    val = sol[-1]
+    for i in range(k):
+        if out[i] <= 0:
+            return False
+    val = out[-1]
 
     if k == m:
         return True
@@ -280,39 +261,3 @@ def _next_k_array(a):
         pos += 1
 
     return a
-
-
-if is_numba_required_installed:
-    @jit(nopython=True, cache=True)
-    def _is_singular(a):
-        s = numba.targets.linalg._compute_singular_values(a)
-        if s[-1] <= s[0] * EPS:
-            return True
-        else:
-            return False
-else:
-    def _is_singular(a):
-        s = np.linalg.svd(a, compute_uv=False)
-        if s[-1] <= s[0] * EPS:
-            return True
-        else:
-            return False
-
-_is_singular_docstr = \
-"""
-Determine whether matrix `a` is numerically singular, by checking
-its singular values.
-
-Parameters
-----------
-a : ndarray(float, ndim=2)
-    2-dimensional array of floats.
-
-Returns
--------
-bool
-    Whether `a` is numerically singular.
-
-"""
-
-_is_singular.__doc__ = _is_singular_docstr

quantecon/game_theory/tests/test_support_enumeration.py

@@ -4,10 +4,33 @@
 Tests for support_enumeration.py
 
 """
+import numpy as np
 from numpy.testing import assert_allclose
+from nose.tools import eq_
+from quantecon.util import check_random_state
 from quantecon.game_theory import Player, NormalFormGame, support_enumeration
 
 
+def random_skew_sym(n, m=None, random_state=None):
+    """
+    Generate a random skew symmetric zero-sum NormalFormGame of the form
+    O    B
+    -B.T O
+    where B is an n x m matrix.
+
+    """
+    if m is None:
+        m = n
+    random_state = check_random_state(random_state)
+    B = random_state.random_sample((n, m))
+    A = np.empty((n+m, n+m))
+    A[:n, :n] = 0
+    A[n:, n:] = 0
+    A[:n, n:] = B
+    A[n:, :n] = -B.T
+    return NormalFormGame([Player(A) for i in range(2)])
+
+
 class TestSupportEnumeration():
     def setUp(self):
         self.game_dicts = []
@@ -35,10 +58,18 @@ def setUp(self):
     def test_support_enumeration(self):
         for d in self.game_dicts:
             NEs_computed = support_enumeration(d['g'])
+            eq_(len(NEs_computed), len(d['NEs']))
             for actions_computed, actions in zip(NEs_computed, d['NEs']):
                 for action_computed, action in zip(actions_computed, actions):
                     assert_allclose(action_computed, action)
 
+    def test_no_error_skew_sym(self):
+        # Test no LinAlgError is raised.
+        n, m = 3, 2
+        seed = 7028
+        g = random_skew_sym(n, m, random_state=seed)
+        NEs = support_enumeration(g)
+
 
 if __name__ == '__main__':
     import sys

quantecon/util/numba.py

@@ -0,0 +1,71 @@
+"""
+Utilities to support Numba jitted functions
+
+"""
+import numpy as np
+from numba import generated_jit, types
+from numba.targets.linalg import _LAPACK
+
+
+# BLAS kinds as letters
+_blas_kinds = {
+    types.float32: 's',
+    types.float64: 'd',
+    types.complex64: 'c',
+    types.complex128: 'z',
+}
+
+
+@generated_jit(nopython=True, cache=True)
+def _numba_linalg_solve(a, b):
+    """
+    Solve the linear equation ax = b directly calling a Numba internal
+    function. The data in `a` and `b` are interpreted in Fortran order,
+    and dtype of `a` and `b` must be the same, one of {float32, float64,
+    complex64, complex128}. `a` and `b` are modified in place, and the
+    solution is stored in `b`. *No error check is made for the inputs.*
+
+    Parameters
+    ----------
+    a : ndarray(ndim=2)
+        2-dimensional ndarray of shape (n, n).
+
+    b : ndarray(ndim=1 or 2)
+        1-dimensional ndarray of shape (n,) or 2-dimensional ndarray of
+        shape (n, nrhs).
+
+    Returns
+    -------
+    r : scalar(int)
+        r = 0 if successful.
+
+    Notes
+    -----
+    From github.com/numba/numba/blob/master/numba/targets/linalg.py
+
+    """
+    numba_xgesv = _LAPACK().numba_xgesv(a.dtype)
+    kind = ord(_blas_kinds[a.dtype])
+
+    def _numba_linalg_solve_impl(a, b):  # pragma: no cover
+        n = a.shape[-1]
+        if b.ndim == 1:
+            nrhs = 1
+        else:  # b.ndim == 2
+            nrhs = b.shape[-1]
+        F_INT_nptype = np.int32
+        ipiv = np.empty(n, dtype=F_INT_nptype)
+
+        r = numba_xgesv(
+            kind,         # kind
+            n,            # n
+            nrhs,         # nhrs
+            a.ctypes,     # a
+            n,            # lda
+            ipiv.ctypes,  # ipiv
+            b.ctypes,     # b
+            n             # ldb
+        )
+        return r
+
+    return _numba_linalg_solve_impl

quantecon/util/tests/test_numba.py

@@ -0,0 +1,59 @@
+"""
+Tests for Numba support utilities
+
+"""
+import numpy as np
+from numpy.testing import assert_array_equal
+from numba import jit
+from nose.tools import eq_, ok_
+from quantecon.util.numba import _numba_linalg_solve
+
+
+@jit(nopython=True)
+def numba_linalg_solve_orig(a, b):
+    return np.linalg.solve(a, b)
+
+
+class TestNumbaLinalgSolve:
+    def setUp(self):
+        self.dtypes = [np.float32, np.float64]
+        self.a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
+        self.b_1dim = np.array([2, 4, -1])
+        self.b_2dim = np.array([[2, 3], [4, 1], [-1, 0]])
+        self.a_singular = np.array([[0, 1, 2], [3, 4, 5], [3, 5, 7]])
+
+    def test_b_1dim(self):
+        for dtype in self.dtypes:
+            a = np.asfortranarray(self.a, dtype=dtype)
+            b = np.asfortranarray(self.b_1dim, dtype=dtype)
+            sol_orig = numba_linalg_solve_orig(a, b)
+            r = _numba_linalg_solve(a, b)
+            eq_(r, 0)
+            assert_array_equal(b, sol_orig)
+
+    def test_b_2dim(self):
+        for dtype in self.dtypes:
+            a = np.asfortranarray(self.a, dtype=dtype)
+            b = np.asfortranarray(self.b_2dim, dtype=dtype)
+            sol_orig = numba_linalg_solve_orig(a, b)
+            r = _numba_linalg_solve(a, b)
+            eq_(r, 0)
+            assert_array_equal(b, sol_orig)
+
+    def test_singular_a(self):
+        for b in [self.b_1dim, self.b_2dim]:
+            for dtype in self.dtypes:
+                a = np.asfortranarray(self.a_singular, dtype=dtype)
+                b = np.asfortranarray(b, dtype=dtype)
+                r = _numba_linalg_solve(a, b)
+                ok_(r != 0)
+
+
+if __name__ == '__main__':
+    import sys
+    import nose
+
+    argv = sys.argv[:]
+    argv.append('--verbose')
+    argv.append('--nocapture')
+    nose.main(argv=argv, defaultTest=__file__)