# Implementing Lemke-Howson in Python

**Daisuke Oyama**  
*Faculty of Economics, University of Tokyo*


```python
import numpy as np
```


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


```python
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:


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


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


```python
tableau0
```




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




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


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


```python
tableau1
```




    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$.


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




    array([3, 4])




```python
basic_vars1 = np.arange(0, m)
basic_vars1
```




    array([0, 1, 2])



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


```python
init_pivot = 1
```


```python
# Current pivot
pivot = init_pivot
```

### Step 1

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


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




    array([ 0.5    ,  0.16667])




```python
row_min = ratios.argmin()
row_min
```




    1




```python
basic_vars0[row_min]
```




    4



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

Update the tableau:


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


```python
tableau0
```




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




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


```python
# 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]]
```


```python
tableau0
```




    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:


```python
basic_vars0[row_min], pivot = pivot, basic_vars0[row_min]
```


```python
basic_vars0
```




    array([3, 1])




```python
basic_vars0[row_min]
```




    1




```python
pivot
```




    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.


```python
pivot == init_pivot
```




    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`.


```python
tableau1
```




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




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




    2




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


```python
tableau1
```




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




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


```python
tableau1
```




    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]])




```python
basic_vars1[row_min], pivot = pivot, basic_vars1[row_min]
```


```python
pivot
```




    2




```python
basic_vars1
```




    array([0, 1, 4])




```python
basic_vars1[row_min]
```




    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.


```python
pivot == init_pivot
```




    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


```python
tableau0
```




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




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




    0




```python
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]]
```


```python
tableau0
```




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




```python
basic_vars0[row_min], pivot = pivot, basic_vars0[row_min]
```


```python
pivot
```




    3




```python
pivot == init_pivot
```




    False



### Step 4


```python
tableau1
```




    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]])




```python
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__':





    1



Note on the warning:

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


```python
tableau1[:, pivot]
```




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



We can just ignore it,
but we can also suppress it by
[`np.errstate`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.errstate.html)
with a `with` clause:


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


```python
ratios
```




    array([ 0.16667,  0.08333,      inf])




```python
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]]
```


```python
tableau1
```




    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]])




```python
basic_vars1[row_min], pivot = pivot, basic_vars1[row_min]
```


```python
pivot
```




    1




```python
pivot == init_pivot
```




    True



Now we have complete labeling, so we are done.

### Obtaining the Nash equilibrium

The basic variables are:


```python
basic_vars0
```




    array([2, 1])



$x_2$ and $x_1$, and


```python
basic_vars1
```




    array([0, 3, 4])



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

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


```python
basic_vars0[basic_vars0 < m]
```




    array([2, 1])



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


```python
basic_vars1[basic_vars1 >= m]
```




    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:


```python
tableau0[basic_vars0 < m, -1]
```




    array([ 0.25 ,  0.125])



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


```python
tableau1[basic_vars1 >= m, -1]
```




    array([ 0.08333,  0.16667])



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


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


```python
x
```




    array([ 0.     ,  0.33333,  0.66667])




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


```python
y
```




    array([ 0.33333,  0.66667])



The Nash equilibrium we have found is:


```python
(x, y)
```




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



## Wrapping the procedure in functions


```python
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
```


```python
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
```


```python
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](http://stackoverflow.com/questions/9038160/break-two-for-loops).


```python
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:]
```


```python
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
```


```python
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]))



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




    ((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


```python
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
```


```python
lemke_howson_all(A, B_T)
```




    [(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).


```python
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)
```




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



## Integer pivoting


```python
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
```


```python
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](http://docs.sympy.org/latest/modules/core.html#rational) class
in [SymPy](http://www.sympy.org)
to represent mixed actions in rational numbers.


```python
import sympy
```


```python
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:]
```


```python
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
```


```python
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
```


```python
lemke_howson_int(A, B_T, init_pivot=1)
```




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




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




    ((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])))




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




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




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




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




```python
lemke_howson_all_int(A, B_T)
```




    [(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]))]




```python
lemke_howson_all_int(A, B_T, rational=True)
```




    [(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


```python
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
```


```python
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:


```python
2/3
```




    0.6666666666666666




```python
1 - 1/3
```




    0.6666666666666667




```python
1 - 1/3 == 2/3
```




    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.)


```python
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
```


```python
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
```


```python
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:


```python
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:


```python
lemke_howson_all(C, D_T)
```




    [(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:


```python
lemke_howson_all_int_lex_min(C, D_T)
```




    [(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:


```python
# 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]])
```


```python
lemke_howson_all_int_lex_min(E, F_T)
```




    [(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:


```python
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)]])
```


```python
lemke_howson_all(C, D_T_eps)
```




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




```python
lemke_howson_all(C, F_T_eps)
```




    [(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]))]




```python

```

86f7df33145c8abdf768767844bda6ffd6776da7
ipynb
Jupyter Notebook
lemke_howson/lemke_howson_py.ipynb
oyamad/theory16
09d0d948b801f1e7d5005ddc1b5582235d5abe0c
MIT
2016-10-26T14:54:49.000Z
2017-02-19T19:15:25.000Z
lemke_howson/lemke_howson_py.ipynb
oyamad/theory16
09d0d948b801f1e7d5005ddc1b5582235d5abe0c
MIT
2016-09-30T05:42:44.000Z
2016-11-26T01:01:13.000Z
lemke_howson/lemke_howson_py.ipynb
oyamad/theory16
09d0d948b801f1e7d5005ddc1b5582235d5abe0c
MIT
2016-09-30T02:20:20.000Z
2020-08-13T20:26:33.000Z
Qwen/Qwen-72B
1. YES
2. YES