Implementing Lemke-Howson in Python

Daisuke Oyama
Faculty of Economics, University of Tokyo



In [1]:

    
import numpy as np



In [2]:

    
np.set_printoptions(precision=5)  # Reduce the number of digits printed



In [3]:

    
A = np.array([[3, 3],
              [2, 5],
              [0 ,6]])
B_T = np.array([[3, 2, 3],
                [2, 6, 1]])

m, n = A.shape  # Numbers of actions of the players

To be consistent with the 0-based indexing in Python, we call the players 0 and 1.

Complementary pivoting

Build the tableau for each player:



In [4]:

    
# Player 0
tableau0 = np.empty((n, m+n+1))
tableau0[:, :m] = B_T
tableau0[:, m:m+n] = np.identity(n)
tableau0[:, -1] = 1



In [5]:

    
# One-line commamd
# tableau0 = np.hstack((B_T, np.identity(n), np.ones((n, 1))))



In [6]:

    
tableau0









    Out[6]:





array([[ 3.,  2.,  3.,  1.,  0.,  1.],
       [ 2.,  6.,  1.,  0.,  1.,  1.]])



In [7]:

    
# Player 1
tableau1 = np.empty((m, n+m+1))
tableau1[:, :m] = np.identity(m)
tableau1[:, m:m+n] = A
tableau1[:, -1] = 1



In [8]:

    
# One-line command
# tableau1 = np.hstack((np.identity(m), A, np.ones((m, 1))))



In [9]:

    
tableau1









    Out[9]:





array([[ 1.,  0.,  0.,  3.,  3.,  1.],
       [ 0.,  1.,  0.,  2.,  5.,  1.],
       [ 0.,  0.,  1.,  0.,  6.,  1.]])

Denote the player 0's variables by $x_0, x_1, x_2$ and $s_3, s_4$, and the player 1's variables by $r_0, r_1, r_2$ and $y_3, y_4$.

The initial basic variables are $s_3, s_4$ and $r_0, r_1, r_2$.



In [10]:

    
basic_vars0 = np.arange(m, m+n)
basic_vars0









    Out[10]:





array([3, 4])



In [11]:

    
basic_vars1 = np.arange(0, m)
basic_vars1









    Out[11]:





array([0, 1, 2])

Let the initial pivot index be 1, so that $x_1$ is to enter the basis:



In [12]:

    
init_pivot = 1



In [13]:

    
# Current pivot
pivot = init_pivot

Step 1

Determine the basic variable to leave by the minimum ratio test:



In [14]:

    
ratios = tableau0[:, -1] / tableau0[:, pivot]
ratios









    Out[14]:





array([ 0.5    ,  0.16667])



In [15]:

    
row_min = ratios.argmin()
row_min









    Out[15]:





1



In [16]:

    
basic_vars0[row_min]









    Out[16]:





4

$s_4$ is the basic variable that leaves the basis.

Update the tableau:



In [17]:

    
tableau0[row_min, :] /= tableau0[row_min, pivot]



In [18]:

    
tableau0









    Out[18]:





array([[ 3.     ,  2.     ,  3.     ,  1.     ,  0.     ,  1.     ],
       [ 0.33333,  1.     ,  0.16667,  0.     ,  0.16667,  0.16667]])



In [19]:

    
for i in range(tableau0.shape[0]):
    if i != row_min:
        tableau0[i, :] -= tableau0[row_min, :] * tableau0[i, pivot]



In [20]:

    
# Another approach by a NumPy trick
# ind = np.ones(tableau0.shape[0], dtype=bool)
# ind[row_min] = False
# tableau0[ind, :] -= tableau0[row_min, :] * tableau0[ind, :][:, [pivot]]



In [21]:

    
tableau0









    Out[21]:





array([[ 2.33333,  0.     ,  2.66667,  1.     , -0.33333,  0.66667],
       [ 0.33333,  1.     ,  0.16667,  0.     ,  0.16667,  0.16667]])

Update the basic variables and the pivot for the next step:



In [22]:

    
basic_vars0[row_min], pivot = pivot, basic_vars0[row_min]



In [23]:

    
basic_vars0









    Out[23]:





array([3, 1])



In [24]:

    
basic_vars0[row_min]









    Out[24]:





1



In [25]:

    
pivot









    Out[25]:





4

That is, $x_1$ has become a basic variable, while $s_4$ becomes a nonbasic variable (i.e., $s_4 = 0$).

If the new pivot is equal to the initial pivot, we are done.



In [26]:

    
pivot == init_pivot









    Out[26]:





False

But this is False, so we continue.

In the next step, the variable $y_4$ which is complementary to $s_4$ (i.e., $y_4 s_4 = 0$) becomes a basic variable.

Step 2

Repeat the same exercise as above for tableau1.



In [27]:

    
tableau1









    Out[27]:





array([[ 1.,  0.,  0.,  3.,  3.,  1.],
       [ 0.,  1.,  0.,  2.,  5.,  1.],
       [ 0.,  0.,  1.,  0.,  6.,  1.]])



In [28]:

    
ratios = tableau1[:, -1] / tableau1[:, pivot]
row_min = ratios.argmin()
row_min









    Out[28]:





2



In [29]:

    
tableau1[row_min, :] /= tableau1[row_min, pivot]



In [30]:

    
tableau1









    Out[30]:





array([[ 1.     ,  0.     ,  0.     ,  3.     ,  3.     ,  1.     ],
       [ 0.     ,  1.     ,  0.     ,  2.     ,  5.     ,  1.     ],
       [ 0.     ,  0.     ,  0.16667,  0.     ,  1.     ,  0.16667]])



In [31]:

    
ind = np.ones(tableau1.shape[0], dtype=bool)
ind[row_min] = False
tableau1[ind, :] -= tableau1[row_min, :] * tableau1[ind, :][:, [pivot]]



In [32]:

    
tableau1









    Out[32]:





array([[ 1.     ,  0.     , -0.5    ,  3.     ,  0.     ,  0.5    ],
       [ 0.     ,  1.     , -0.83333,  2.     ,  0.     ,  0.16667],
       [ 0.     ,  0.     ,  0.16667,  0.     ,  1.     ,  0.16667]])



In [33]:

    
basic_vars1[row_min], pivot = pivot, basic_vars1[row_min]



In [34]:

    
pivot









    Out[34]:





2



In [35]:

    
basic_vars1









    Out[35]:





array([0, 1, 4])



In [36]:

    
basic_vars1[row_min]









    Out[36]:





4

That is, $y_4$ has become a basic variable, while $r_2$ becomes a nonbasic variable.

If the new pivot is equal to the initial pivot, we are done.



In [37]:

    
pivot == init_pivot









    Out[37]:





False

But this is False, so we continue.

In the next step, the variable $x_2$ which is complementary to $r_2$ becomes a basic variable.

Step 3



In [38]:

    
tableau0









    Out[38]:





array([[ 2.33333,  0.     ,  2.66667,  1.     , -0.33333,  0.66667],
       [ 0.33333,  1.     ,  0.16667,  0.     ,  0.16667,  0.16667]])



In [39]:

    
ratios = tableau0[:, -1] / tableau0[:, pivot]
row_min = ratios.argmin()
row_min









    Out[39]:





0



In [40]:

    
tableau0[row_min, :] /= tableau0[row_min, pivot]
ind = np.ones(tableau0.shape[0], dtype=bool)
ind[row_min] = False
tableau0[ind, :] -= tableau0[row_min, :] * tableau0[ind, :][:, [pivot]]



In [41]:

    
tableau0









    Out[41]:





array([[ 0.875 ,  0.    ,  1.    ,  0.375 , -0.125 ,  0.25  ],
       [ 0.1875,  1.    ,  0.    , -0.0625,  0.1875,  0.125 ]])



In [42]:

    
basic_vars0[row_min], pivot = pivot, basic_vars0[row_min]



In [43]:

    
pivot









    Out[43]:





3



In [44]:

    
pivot == init_pivot









    Out[44]:





False

Step 4



In [45]:

    
tableau1









    Out[45]:





array([[ 1.     ,  0.     , -0.5    ,  3.     ,  0.     ,  0.5    ],
       [ 0.     ,  1.     , -0.83333,  2.     ,  0.     ,  0.16667],
       [ 0.     ,  0.     ,  0.16667,  0.     ,  1.     ,  0.16667]])



In [46]:

    
ratios = tableau1[:, -1] / tableau1[:, pivot]
row_min = ratios.argmin()
row_min









    



/Applications/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in true_divide
  if __name__ == '__main__':






    Out[46]:





1

Note on the warning:

tableau1[:, pivot] has a zero entry, so we get a "divide by zero" warning.



In [47]:

    
tableau1[:, pivot]









    Out[47]:





array([ 3.,  2.,  0.])

We can just ignore it, but we can also suppress it by np.errstate with a with clause:



In [48]:

    
with np.errstate(divide='ignore'):
    ratios = tableau1[:, -1] / tableau1[:, pivot]



In [49]:

    
ratios









    Out[49]:





array([ 0.16667,  0.08333,      inf])



In [50]:

    
tableau1[row_min, :] /= tableau1[row_min, pivot]
ind = np.ones(tableau1.shape[0], dtype=bool)
ind[row_min] = False
tableau1[ind, :] -= tableau1[row_min, :] * tableau1[ind, :][:, [pivot]]



In [51]:

    
tableau1









    Out[51]:





array([[ 1.     , -1.5    ,  0.75   ,  0.     ,  0.     ,  0.25   ],
       [ 0.     ,  0.5    , -0.41667,  1.     ,  0.     ,  0.08333],
       [ 0.     ,  0.     ,  0.16667,  0.     ,  1.     ,  0.16667]])



In [52]:

    
basic_vars1[row_min], pivot = pivot, basic_vars1[row_min]



In [53]:

    
pivot









    Out[53]:





1



In [54]:

    
pivot == init_pivot









    Out[54]:





True

Now we have complete labeling, so we are done.

Obtaining the Nash equilibrium

The basic variables are:



In [55]:

    
basic_vars0









    Out[55]:





array([2, 1])

$x_2$ and $x_1$, and



In [56]:

    
basic_vars1









    Out[56]:





array([0, 3, 4])

$r_0$, $y_3$, and $y_4$.

The indices of the basic variables corresponding to $x$:



In [57]:

    
basic_vars0[basic_vars0 < m]









    Out[57]:





array([2, 1])

The indices of the basic variables corresponding to $y$:



In [58]:

    
basic_vars1[basic_vars1 >= m]









    Out[58]:





array([3, 4])

The values of the basic variables are stored in the last columns of the tableaux.

The values of $x_2$ and $x_1$ are:



In [59]:

    
tableau0[basic_vars0 < m, -1]









    Out[59]:





array([ 0.25 ,  0.125])

The values of $y_3$ and $y_4$ are:



In [60]:

    
tableau1[basic_vars1 >= m, -1]









    Out[60]:





array([ 0.08333,  0.16667])

We need to normalize these values so that $x$ and $y$ are probability distributions.



In [61]:

    
x = np.zeros(m)
x[basic_vars0[basic_vars0 < m]] = tableau0[basic_vars0 < m, -1]
x /= x.sum()



In [62]:

    
x









    Out[62]:





array([ 0.     ,  0.33333,  0.66667])



In [63]:

    
y = np.zeros(n)
y[basic_vars1[basic_vars1 >= m] - m] = tableau1[basic_vars1 >= m, -1]
y /= y.sum()



In [64]:

    
y









    Out[64]:





array([ 0.33333,  0.66667])

The Nash equilibrium we have found is:



In [65]:

    
(x, y)









    Out[65]:





(array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667]))

Wrapping the procedure in functions



In [66]:

    
def min_ratio_test(tableau, pivot):
    ind_nonpositive = tableau[:, pivot] <= 0
    with np.errstate(divide='ignore', invalid='ignore'):
        ratios = tableau[:, -1] / tableau[:, pivot]
    ratios[ind_nonpositive] = np.inf
    row_min = ratios.argmin()
    return row_min



In [67]:

    
def pivoting(tableau, pivot, pivot_row):
    """
    Perform a pivoting step.
    
    Modify `tableau` in place (and return its view).
    
    """
    # Row indices except pivot_row
    ind = np.ones(tableau.shape[0], dtype=bool)
    ind[pivot_row] = False
    
    # Store the values in the pivot column, except for row_min
    # Made 2-dim by np.newaxis
    multipliers = tableau[ind, pivot, np.newaxis]
    
    # Update the tableau
    tableau[pivot_row, :] /= tableau[pivot_row, pivot]
    tableau[ind, :] -= tableau[pivot_row, :] * multipliers
    
    return tableau



In [68]:

    
def lemke_howson_tbl(tableau0, tableau1, basic_vars0, basic_vars1, init_pivot):
    m, n = tableau1.shape[0], tableau0.shape[0]
    tableaux = (tableau0, tableau1)
    basic_vars = (basic_vars0, basic_vars1)
    init_player = int((basic_vars[0]==init_pivot).any())
    players = [init_player, 1 - init_player]

    pivot = init_pivot

    while True:
        for i in players:
            # Determine the leaving variable
            row_min = min_ratio_test(tableaux[i], pivot)
            
            # Pivoting step: modify tableau in place
            pivoting(tableaux[i], pivot, row_min)
            
            # Update the basic variables and the pivot
            basic_vars[i][row_min], pivot = pivot, basic_vars[i][row_min]
            
            if pivot == init_pivot:
                break
        else:
            continue
        break

    out_dtype = np.result_type(*tableaux)
    out = np.zeros(m+n, dtype=out_dtype)
    for i, (start, num) in enumerate(zip((0, m), (m, n))):
        ind = basic_vars[i] < start + num if i == 0 else start <= basic_vars[i]
        out[basic_vars[i][ind]] = tableaux[i][ind, -1]
        
    return out

Note: There is no nested break in Python; see e.g., break two for loops.



In [69]:

    
def normalize(unnormalized, m, n):
    normalized = np.empty(m+n)
    for (start, num) in zip((0, m), (m, n)):
        s = unnormalized[start:start+num].sum()
        if s != 0:
            normalized[start:start+num] = unnormalized[start:start+num] / s
        else:
            normalized[start:start+num] = 0
    return normalized[:m], normalized[m:]



In [70]:

    
def lemke_howson(A, B_T, init_pivot=0, return_tableaux=False):
    m, n = A.shape
    tableaux = (np.hstack((B_T, np.identity(n), np.ones((n, 1)))),
                np.hstack((np.identity(m), A, np.ones((m, 1)))))
    basic_vars = (np.arange(m, m+n), np.arange(0, m))
    unnormalized = lemke_howson_tbl(*tableaux, *basic_vars, init_pivot)
    normalized = normalize(unnormalized, m, n)
    if return_tableaux:
        return normalized, tableaux, basic_vars
    return normalized



In [71]:

    
init_pivot = 1
x, y = lemke_howson(A, B_T, init_pivot)
print("Nash equilibrium found\n", (x, y))









    



Nash equilibrium found
 (array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667]))



In [72]:

    
lemke_howson(A, B_T, init_pivot, return_tableaux=True)









    Out[72]:





((array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667])),
 (array([[ 0.875 ,  0.    ,  1.    ,  0.375 , -0.125 ,  0.25  ],
         [ 0.1875,  1.    ,  0.    , -0.0625,  0.1875,  0.125 ]]),
  array([[ 1.     , -1.5    ,  0.75   ,  0.     ,  0.     ,  0.25   ],
         [ 0.     ,  0.5    , -0.41667,  1.     ,  0.     ,  0.08333],
         [ 0.     ,  0.     ,  0.16667,  0.     ,  1.     ,  0.16667]])),
 (array([2, 1]), array([0, 3, 4])))

Enumerating all equilibria that are reached by Lemke-Howson paths



In [73]:

    
def lemke_howson_all(A, B_T):
    m, n = A.shape
    k = 0
    NEs = []
    basic_vars_list = []
    player = (m <= n)
    init_pivot = k
    actions, tableaux, basic_vars = \
        lemke_howson(A, B_T, init_pivot, return_tableaux=True)
    NEs.append(actions)
    basic_vars_list.append(np.sort(basic_vars[player]))

    for a in range(m+n):
        if a == k:
            continue
        init_pivot = a
        actions, tableaux, basic_vars = \
            lemke_howson(A, B_T, init_pivot, return_tableaux=True)
        basic_vars_sorted = np.sort(basic_vars[player])
        for arr in basic_vars_list:
            if np.array_equal(basic_vars_sorted, arr):
                break
        else:
            NEs.append(actions)
            basic_vars_list.append(basic_vars_sorted)

            unnormalized = \
                lemke_howson_tbl(*tableaux, *basic_vars, init_pivot=k)
            NEs.append(normalize(unnormalized, m, n))
            basic_vars_list.append(np.sort(basic_vars[player]))
            
    return NEs



In [74]:

    
lemke_howson_all(A, B_T)









    Out[74]:





[(array([ 1.,  0.,  0.]), array([ 1.,  0.])),
 (array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667])),
 (array([ 0.8,  0.2,  0. ]), array([ 0.66667,  0.33333]))]

There are games in which some of the Nash equilibria cannot be reached by Lemke-Howson path from the origin, as in the game in (3.7) in von Stengel (2007).



In [75]:

    
payoff_matrix_3_7 = np.array([[3, 3, 0], [4, 0, 1], [0, 4, 5]])
lemke_howson_all(payoff_matrix_3_7, payoff_matrix_3_7)









    Out[75]:





[(array([ 0.,  0.,  1.]), array([ 0.,  0.,  1.]))]

Integer pivoting



In [76]:

    
def pivoting_int(tableau, pivot, pivot_row, prev_pivot_el):
    """
    Perform a pivoting step with integer input data.
    
    Modify `tableau` in place (and return its view).
    
    """
    # Row indices except pivot_row
    ind = np.ones(tableau.shape[0], dtype=bool)
    ind[pivot_row] = False
    
    # Store the values in the pivot column, except for row_min
    # Made 2-dim by np.newaxis
    multipliers = tableau[ind, pivot, np.newaxis]
    
    # Update the tableau
    tableau[ind, :] *= tableau[pivot_row, pivot]
    tableau[ind, :] -= tableau[pivot_row, :] * multipliers
    tableau[ind, :] //= prev_pivot_el  # Floor division: return int
    
    return tableau



In [77]:

    
def lemke_howson_tbl_int(tableau0, tableau1, basic_vars0, basic_vars1, init_pivot):
    m, n = tableau1.shape[0], tableau0.shape[0]
    tableaux = (tableau0, tableau1)
    basic_vars = (basic_vars0, basic_vars1)
    init_player = int((basic_vars[0]==init_pivot).any())
    players = [init_player, 1 - init_player]

    pivot = init_pivot
    prev_pivot_els = np.ones(2, dtype=np.int_)

    while True:
        for i in players:
            # Determine the leaving variable
            row_min = min_ratio_test(tableaux[i], pivot)
            
            # Pivoting step: modify tableau in place
            pivoting_int(tableaux[i], pivot, row_min, prev_pivot_els[i])
            
            # Backup the pivot element
            prev_pivot_els[i] = tableaux[i][row_min, pivot]
            
            # Update the basic variables and the pivot
            basic_vars[i][row_min], pivot = pivot, basic_vars[i][row_min]
            
            if pivot == init_pivot:
                break
        else:
            continue
        break

    out = np.zeros(m+n, dtype=np.int_)
    for i, (start, num) in enumerate(zip((0, m), (m, n))):
        ind = basic_vars[i] < start + num if i == 0 else start <= basic_vars[i]
        out[basic_vars[i][ind]] = tableaux[i][ind, -1]
        
    return out

Let us use the Rational class in SymPy to represent mixed actions in rational numbers.



In [78]:

    
import sympy



In [79]:

    
def normalize_rational(unnormalized, m, n):
    """
    Normalize the integer array `unnormalized` with ratioal numbers.
    
    """
    normalized = np.empty(m+n, np.object_)
    for (start, num) in zip((0, m), (m, n)):
        s = unnormalized[start:start+num].sum()
        if s != 0:
            for k in range(start, start+num):
                normalized[k] = sympy.Rational(sympy.S(unnormalized[k]), sympy.S(s))
        else:
            normalized[start:start+num] = sympy.Rational(0)
    return normalized[:m], normalized[m:]



In [80]:

    
def lemke_howson_int(A, B_T, init_pivot=0, rational=False, return_tableaux=False):
    m, n = A.shape
    tableaux = (np.hstack((B_T, np.identity(n, dtype=np.int_), np.ones((n, 1), dtype=np.int_))),
                np.hstack((np.identity(m, dtype=np.int_), A, np.ones((m, 1), dtype=np.int_))))
    basic_vars = (np.arange(m, m+n), np.arange(0, m))
    unnormalized = lemke_howson_tbl_int(*tableaux, *basic_vars, init_pivot)
    if rational:
        normalized = normalize_rational(unnormalized, m, n)
    else:
        normalized = normalize(unnormalized, m, n)
    if return_tableaux:
        return normalized, tableaux, basic_vars
    return normalized



In [81]:

    
def lemke_howson_all_int(A, B_T, rational=False):
    m, n = A.shape
    k = 0
    NEs = []
    basic_vars_list = []
    player = (m <= n)
    init_pivot = k
    actions, tableaux, basic_vars = \
        lemke_howson_int(A, B_T, init_pivot, rational=rational, return_tableaux=True)
    NEs.append(actions)
    basic_vars_list.append(np.sort(basic_vars[player]))

    for a in range(m+n):
        if a == k:
            continue
        init_pivot = a
        actions, tableaux, basic_vars = \
            lemke_howson_int(A, B_T, init_pivot, rational=rational, return_tableaux=True)
        basic_vars_sorted = np.sort(basic_vars[player])
        for arr in basic_vars_list:
            if np.array_equal(basic_vars_sorted, arr):
                break
        else:
            NEs.append(actions)
            basic_vars_list.append(basic_vars_sorted)

            unnormalized = \
                lemke_howson_tbl_int(*tableaux, *basic_vars, init_pivot=k)
            if rational:
                normalized = normalize_rational(unnormalized, m, n)
            else:
                normalized = normalize(unnormalized, m, n)
            NEs.append(normalized)
            basic_vars_list.append(np.sort(basic_vars[player]))
            
    return NEs



In [82]:

    
lemke_howson_int(A, B_T, init_pivot=1)









    Out[82]:





(array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667]))



In [83]:

    
lemke_howson_int(A, B_T, init_pivot=1, return_tableaux=True)









    Out[83]:





((array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667])),
 (array([[14,  0, 16,  6, -2,  4],
         [ 3, 16,  0, -1,  3,  2]]), array([[ 12, -18,   9,   0,   0,   3],
         [  0,   6,  -5,  12,   0,   1],
         [  0,   0,   2,   0,  12,   2]])),
 (array([2, 1]), array([0, 3, 4])))



In [84]:

    
lemke_howson_int(A, B_T, init_pivot=1, rational=True)









    Out[84]:





(array([0, 1/3, 2/3], dtype=object), array([1/3, 2/3], dtype=object))



In [85]:

    
lemke_howson_int(A, B_T, init_pivot=0, rational=True)









    Out[85]:





(array([1, 0, 0], dtype=object), array([1, 0], dtype=object))



In [86]:

    
lemke_howson_all_int(A, B_T)









    Out[86]:





[(array([ 1.,  0.,  0.]), array([ 1.,  0.])),
 (array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667])),
 (array([ 0.2,  0.8,  0. ]), array([ 0.66667,  0.33333]))]



In [87]:

    
lemke_howson_all_int(A, B_T, rational=True)









    Out[87]:





[(array([1, 0, 0], dtype=object), array([1, 0], dtype=object)),
 (array([0, 1/3, 2/3], dtype=object), array([1/3, 2/3], dtype=object)),
 (array([1/5, 4/5, 0], dtype=object), array([2/3, 1/3], dtype=object))]

Lexico-minimum ratio test



In [88]:

    
def min_ratio_test_no_tie_breaking(tableau, pivot, test_col, argmins, num_argmins):
    idx = 0
    i = argmins[idx]
    if tableau[i, pivot] > 0:
        min_ratio = tableau[i, test_col] / tableau[i, pivot]
    else:
        min_ratio = np.inf
    for k in range(1, num_argmins):
        i = argmins[k]
        if tableau[i, pivot] <= 0:
            continue
        ratio = tableau[i, test_col] / tableau[i, pivot]
        if ratio > min_ratio:
            continue
        elif ratio < min_ratio:
            min_ratio = ratio
            idx = 0
        elif ratio == min_ratio:
            idx += 1
        argmins[idx] = k
    return idx + 1



In [89]:

    
def lex_min_ratio_test(tableau, pivot, slack_start):
    num_rows = tableau.shape[0]
    argmins = np.arange(num_rows)
    num_argmins = num_rows
    
    num_argmins = min_ratio_test_no_tie_breaking(tableau, pivot, -1, argmins, num_argmins)
    if num_argmins == 1:
        return argmins[0]
    
    for j in range(slack_start, slack_start+num_rows):
        if j == pivot:
            continue
        num_argmins = min_ratio_test_no_tie_breaking(tableau, pivot, j, argmins, num_argmins)
        if num_argmins == 1:
            break
    return argmins[0]

Caveat: Because of rounding errors, one should not rely on equality between floating point numbers. For example:



In [90]:

    
2/3









    Out[90]:





0.6666666666666666



In [91]:

    
1 - 1/3









    Out[91]:





0.6666666666666667



In [92]:

    
1 - 1/3 == 2/3









    Out[92]:





False

Note: In comparing $\frac{a}{b}$ and $\frac{a'}{b'}$, one may instead compare $a b'$ and $a' b$: when these are integers, the latter involves only integers. (This, of course, does not apply only for lexico-minimum test.)



In [93]:

    
def lemke_howson_tbl_int_lex_min(tableau0, tableau1, basic_vars0, basic_vars1, init_pivot):
    m, n = tableau1.shape[0], tableau0.shape[0]
    tableaux = (tableau0, tableau1)
    basic_vars = (basic_vars0, basic_vars1)
    init_player = int((basic_vars[0]==init_pivot).any())
    players = [init_player, 1 - init_player]

    pivot = init_pivot
    prev_pivot_els = np.ones(2, dtype=np.int_)
    slack_starts = (m, 0)

    while True:
        for i in players:
            # Determine the leaving variable
            row_min = lex_min_ratio_test(tableaux[i], pivot, slack_starts[i])
            
            # Pivoting step: modify tableau in place
            pivoting_int(tableaux[i], pivot, row_min, prev_pivot_els[i])
            
            # Backup the pivot element
            prev_pivot_els[i] = tableaux[i][row_min, pivot]
            
            # Update the basic variables and the pivot
            basic_vars[i][row_min], pivot = pivot, basic_vars[i][row_min]
            
            if pivot == init_pivot:
                break
        else:
            continue
        break

    out = np.zeros(m+n, dtype=np.int_)
    for i, (start, num) in enumerate(zip((0, m), (m, n))):
        ind = basic_vars[i] < start + num if i == 0 else start <= basic_vars[i]
        out[basic_vars[i][ind]] = tableaux[i][ind, -1]
        
    return out



In [94]:

    
def lemke_howson_int_lex_min(A, B_T, init_pivot=0, rational=False, return_tableaux=False):
    m, n = A.shape
    tableaux = (np.hstack((B_T, np.identity(n, dtype=np.int_), np.ones((n, 1), dtype=np.int_))),
                np.hstack((np.identity(m, dtype=np.int_), A, np.ones((m, 1), dtype=np.int_))))
    basic_vars = (np.arange(m, m+n), np.arange(0, m))
    unnormalized = lemke_howson_tbl_int_lex_min(*tableaux, *basic_vars, init_pivot)
    if rational:
        normalized = normalize_rational(unnormalized, m, n)
    else:
        normalized = normalize(unnormalized, m, n)
    if return_tableaux:
        return normalized, tableaux, basic_vars
    return normalized



In [95]:

    
def lemke_howson_all_int_lex_min(A, B_T, rational=False):
    m, n = A.shape
    k = 0
    NEs = []
    basic_vars_list = []
    player = (m <= n)
    init_pivot = k
    actions, tableaux, basic_vars = \
        lemke_howson_int_lex_min(A, B_T, init_pivot, rational=rational, return_tableaux=True)
    NEs.append(actions)
    basic_vars_list.append(np.sort(basic_vars[player]))

    for a in range(m+n):
        if a == k:
            continue
        init_pivot = a
        actions, tableaux, basic_vars = \
            lemke_howson_int_lex_min(A, B_T, init_pivot, rational=rational, return_tableaux=True)
        basic_vars_sorted = np.sort(basic_vars[player])
        for arr in basic_vars_list:
            if np.array_equal(basic_vars_sorted, arr):
                break
        else:
            NEs.append(actions)
            basic_vars_list.append(basic_vars_sorted)

            unnormalized = \
                lemke_howson_tbl_int_lex_min(*tableaux, *basic_vars, init_pivot=k)
            if rational:
                normalized = normalize_rational(unnormalized, m, n)
            else:
                normalized = normalize(unnormalized, m, n)
            NEs.append(normalized)
            basic_vars_list.append(np.sort(basic_vars[player]))
            
    return NEs

Consider the following degenerate game:



In [96]:

    
C = np.array([[3, 3],
              [2, 5],
              [0 ,6]])
D_T = np.array([[3, 2, 3],
                [3, 6, 1]])

lemke_howson_all fails to work properly:



In [97]:

    
lemke_howson_all(C, D_T)









    Out[97]:





[(array([ 1.,  0.,  0.]), array([ 1.,  0.])),
 (array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667])),
 (array([ 0.,  0.,  0.]), array([ 0.,  0.]))]

With lexico-minimum test:



In [98]:

    
lemke_howson_all_int_lex_min(C, D_T)









    Out[98]:





[(array([ 0.     ,  0.33333,  0.66667]), array([ 0.33333,  0.66667]))]

Due to the paricular fixed way of introducing the perturbations $(\varepsilon^1, \ldots, \varepsilon^k)$, the output does depend on the ordering of the actions:



In [99]:

    
# Change the order of the actions of player 1
E = np.array([[3, 3],
              [5, 2],
              [6 ,0]])
F_T = np.array([[3, 6, 1],
                [3, 2, 3]])



In [100]:

    
lemke_howson_all_int_lex_min(E, F_T)









    Out[100]:





[(array([ 1.,  0.,  0.]), array([ 0.,  1.])),
 (array([ 0.     ,  0.33333,  0.66667]), array([ 0.66667,  0.33333])),
 (array([ 1.,  0.,  0.]), array([ 0.33333,  0.66667]))]

Essentially, the exercise corresponds to considering $$ \begin{pmatrix} \dfrac{3}{1+\varepsilon_1} & \dfrac{2}{1+\varepsilon_1} & \dfrac{3}{1+\varepsilon_1} \\ \dfrac{3}{1+\varepsilon_2} & \dfrac{6}{1+\varepsilon_2} & \dfrac{1}{1+\varepsilon_2} \end{pmatrix}, $$ where $(\varepsilon_1, \varepsilon_2) =(\varepsilon^1, \varepsilon^2)$ or $(\varepsilon_1, \varepsilon_2) =(\varepsilon^2, \varepsilon^1)$. In fact:



In [101]:

    
eps = 0.01
D_T_eps = np.array([[3 / (1 + eps), 2 / (1 + eps), 3 / (1 + eps)],
                    [3 / (1 + eps**2), 6 / (1 + eps**2), 1 / (1 + eps**2)]])
F_T_eps = np.array([[3 / (1 + eps**2), 2 / (1 + eps**2), 3 / (1 + eps**2)],
                    [3 / (1 + eps), 6 / (1 + eps), 1 / (1 + eps)]])



In [102]:

    
lemke_howson_all(C, D_T_eps)









    Out[102]:





[(array([ 0.     ,  0.32897,  0.67103]), array([ 0.33333,  0.66667]))]



In [103]:

    
lemke_howson_all(C, F_T_eps)









    Out[103]:





[(array([ 1.,  0.,  0.]), array([ 1.,  0.])),
 (array([ 0.     ,  0.33773,  0.66227]), array([ 0.33333,  0.66667])),
 (array([ 0.99259,  0.00741,  0.     ]), array([ 0.66667,  0.33333]))]



In [ ]: