Nash equilibria in bimatrix games
¶

Alfred Galichon (NYU & Sciences Po) and Antoine Jacquet (Sciences Po)
¶

'math+econ+code' masterclass series
¶

With python code examples
¶

© 2018–2025 by Alfred Galichon and Antoine Jacquet. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.

If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/

References¶

  • Mangasarian and Stone (1964). "Two-person nonzero-sum games and quadratic programming." Journal of Mathematical Analysis and Applications.
  • Lemke and Howson (1964). "Equilibrium points of bimatrix games." SIAM Journal on Applied Mathematics.
  • Bruno Codenotti. Computational Aspects of Game Theory. Bertinoro Spring School 2011, Lecture 11: The Lemke–Howson Algorithm

http://wwwold.iit.cnr.it/staff/bruno.codenotti/lecture11p.pdf

Learning objectives¶

  • KKT conditions for Nash equilibrium
  • Mangasarian–Stone quadratic programming formulation
  • Lemke–Howson algorithm
In [1]:
#!pip install mec --upgrade
import numpy as np
import gurobipy as grb
import mec

Nash equilibrium in a two-player bimatrix game¶

Consider a two-player game where if player 1 plays $i \in \mathcal I$ and player 2 plays $j \in \mathcal J$, the payoff to player 1 is $A_{ij} > 0$ and the payoff to player 2 is $B_{ij} > 0$.

Remark. Assuming positive payoffs is without loss of generality.

Let's begin by defining a Python class for bimatrix games.

In [2]:
class Bimatrix_game:
    def __init__(self, A_i_j, B_i_j):
        if A_i_j.shape != B_i_j.shape:
            raise ValueError("A_i_j and B_i_j must be of the same size.")
        self.A_i_j, self.B_i_j = A_i_j, B_i_j
        self.nbi,self.nbj = A_i_j.shape

Now recall the definition of a Nash equilibrium: it is a pair of vectors $p = (p_i)_{i \in \mathcal I}, q = (q_j)_{j \in \mathcal J}$ which simultaneously solve

\begin{equation} \max_{p \geq 0} \big\{ p^\top A q \;\big|\; {\textstyle\sum_i} p_i = 1 \big\}, \qquad\qquad \max_{q \geq 0} \big\{ p^\top B q \;\big|\; {\textstyle\sum_j} q_j = 1 \big\}. \end{equation}

By writing the KKT conditions for these two linear programs we obtain necessary and sufficient conditions for $p$ and $q$ to be a Nash equilibrium: there must exist two real numbers $\alpha$ and $\beta$ such that

\begin{align} p_i &\geq 0 \quad (\forall i) \\ q_j &\geq 0 \quad (\forall j) \\ \textstyle\sum_i p_i &= 1 \\ \textstyle\sum_j q_j &= 1 \\ \alpha - (A q)_i &\geq 0 \quad (\forall i) \\ \beta - (B^\top p)_j &\geq 0 \quad (\forall j) \\ \textstyle\sum_i p_i \big( \alpha - (A q)_i \big) &= 0 \\ \textstyle\sum_j q_j \big( \beta - (B^\top p)_j \big) &= 0. \end{align}

Question. What is the interpretation of $\alpha$ and $\beta$?

Mangasarian and Stone formulation¶

Now let's look at the following quadratic program (QP) formulation proposed by Mangasarian and Stone:

\begin{align} \min_{p \geq 0, q \geq 0, \alpha \geq 0, \beta \geq 0} &\left\{ \alpha + \beta - p^\top (A + B) q \right\} \\ \text{s.t.} ~ & \alpha 1_{\mathcal I} - Aq \geq 0 \\ & \beta 1_{\mathcal J} - B^\top p \geq 0 \\ & -1 + 1_{\mathcal I}^\top p = 0 \\ & -1 + 1_{\mathcal J}^\top q = 0. \end{align}

Claim. The solutions to the Mangasarian–Stone quadratic program are exactly the Nash equilibria of the bimatrix game $(A, B)$.

Proof. Because of the two equality constraints, the objective can be rewritten as $p^\top \big( \alpha 1_{\mathcal I} - Aq \big) + \big(\beta 1_{\mathcal J} - B^\top p\big)^\top q$, which is always non-negative under the constraints of the program.

Furthermore, any Nash equilibrium satisfies the quadratic program's constraints, and attains an objective value of 0. Hence Nash equilibria are solutions to this quadratic program.

Conversely, since we know that Nash equilibria exist, the value of the program must be 0. As a consequence, any solution of this program must satisfy the system derived above, so must be a Nash equilibrium.

Zero-sum games¶

Assume $B = c - A$ with $c$ a constant.

In [3]:
penalty_data = np.array([[53.21, 71.35, 93.80], 
                         [90.26, 42.81, 86.12], 
                         [96.88, 100.0, 75.43]])

penalty_zero_sum = Bimatrix_game(A_i_j = penalty_data,
                                 B_i_j = 100 - penalty_data)
print('A_i_j =\n', penalty_zero_sum.A_i_j)

A_i_j =
 [[ 53.21  71.35  93.8 ]
 [ 90.26  42.81  86.12]
 [ 96.88 100.    75.43]]

Here only $A_{ij}$ is used. $B_{ij}$ is ignored, and the game is solved as a zero-sum game.

Recall that in a zero-sum game, we can solve for the equilibrium using a LP from one of the two players' perspective.
For instance, from the perspective of player 1:

\begin{align} \min_{x_i \geq 0} ~& \sum_i x_i \\ \text{s.t.} ~& \sum_i A_{ij} x_i \geq 1 \quad [y_j \geq 0]. \end{align}

We then recover the strategies using $p_i = \dfrac{x_i}{\sum_i x_i}$, $q_j = \dfrac{y_j}{\sum_i x_i}$.

Solution using Gurobi¶

We shall solve the LP problem using Gurobi.

In [4]:
def Bimatrix_game_zero_sum_solve(self):
    model=grb.Model()
    model.Params.OutputFlag = 0
    x = model.addMVar(shape=self.nbi, name="x")
    model.setObjective(np.ones(self.nbi) @ x, grb.GRB.MINIMIZE)
    model.addConstr(self.A_i_j.T @ x >= np.ones(self.nbj))
    model.optimize() 
    U = 1 / model.ObjVal
    p_i = U * np.array(model.getAttr('x'))
    q_j = U * np.array(model.getAttr('pi'))
    sol_dict = {'val': U, 'p_i': p_i, 'q_j': q_j}
    return(sol_dict)

Bimatrix_game.zero_sum_solve = Bimatrix_game_zero_sum_solve
In [5]:
penalty_zero_sum.zero_sum_solve()

Restricted license - for non-production use only - expires 2025-11-24
Out[5]:
{'val': 82.62606457928055,
 'p_i': array([0.3033884 , 0.15180727, 0.54480433]),
 'q_j': array([0.21908541, 0.10161508, 0.67929951])}

Nonzero-sum games¶

Consider a version of the game with another $B$, which is not zero-sum.
Here, imagine for instance that the goalkeeper gets some bonus points if he dives where the player shoots, whether or not he catches the ball.

In [6]:
penalty_nonzero_sum = Bimatrix_game(A_i_j = penalty_data,
                                    B_i_j = np.array([[150, 100, 100],
                                                      [100, 150, 100],
                                                      [100, 100, 150]])
                                            - penalty_data)
print('A_i_j =\n', penalty_nonzero_sum.A_i_j)
print('B_i_j =\n', penalty_nonzero_sum.B_i_j)

A_i_j =
 [[ 53.21  71.35  93.8 ]
 [ 90.26  42.81  86.12]
 [ 96.88 100.    75.43]]
B_i_j =
 [[ 96.79  28.65   6.2 ]
 [  9.74 107.19  13.88]
 [  3.12   0.    74.57]]

Solution using the Mangasarian–Stone formulation¶

Recall the Mangasarian–Stone formulation:

\begin{align} \min_{p \geq 0, q \geq 0, \alpha \geq 0, \beta \geq 0} &\left\{ \alpha + \beta - p^\top (A + B) q \right\} \\ \text{s.t.} ~ & \alpha 1_{\mathcal I} - Aq \geq 0 \\ & \beta 1_{\mathcal J} - B^\top p \geq 0 \\ & -1 + 1_{\mathcal I}^\top p = 0 \\ & -1 + 1_{\mathcal J}^\top q = 0. \end{align}
In [7]:
def Bimatrix_game_mangasarian_stone_solve(self, verbose=0):
    model=grb.Model()
    model.Params.OutputFlag = 0
    model.params.NonConvex = 2
    p_i = model.addMVar(shape = self.nbi)
    q_j = model.addMVar(shape = self.nbj)
    α = model.addMVar(shape = 1)
    β = model.addMVar(shape = 1)
    model.setObjective(α + β - p_i @ (self.A_i_j + self.B_i_j) @ q_j, sense = grb.GRB.MINIMIZE)
    model.addConstr(α * np.ones((self.nbi,1)) - self.A_i_j @ q_j >= 0)
    model.addConstr(β * np.ones((self.nbj,1)) - self.B_i_j.T @ p_i >= 0)
    model.addConstr(p_i.sum() == 1)
    model.addConstr(q_j.sum() == 1)
    model.optimize() 
    sol = np.array(model.getAttr('x'))
    if verbose > 0: print('p_i =', sol[:self.nbi], '\nq_j =', sol[self.nbi:(self.nbi+self.nbj)])
    return {'p_i': sol[:self.nbi], 'q_j': sol[self.nbi:(self.nbi+self.nbj)],
            'val1': sol[-2], 'val2': sol[-1]}

Bimatrix_game.mangasarian_stone_solve = Bimatrix_game_mangasarian_stone_solve
In [8]:
penalty_zero_sum.mangasarian_stone_solve()
Out[8]:
{'p_i': array([0.3033884 , 0.15180727, 0.54480433]),
 'q_j': array([0.21908541, 0.10161508, 0.67929951]),
 'val1': 82.62606457928055,
 'val2': 17.373935420719445}
In [9]:
penalty_nonzero_sum.mangasarian_stone_solve()
Out[9]:
{'p_i': array([0.3374337 , 0.24917907, 0.41338723]),
 'q_j': array([0.21908541, 0.10161508, 0.67929951]),
 'val1': 82.62606457928055,
 'val2': 36.37698012002112}

Let's verify that this is indeed a Nash equilibrium:

In [10]:
def Bimatrix_game_is_NashEq(self, p_i, q_j, tol=1e-5):
    for i in range(self.nbi):
        if np.eye(self.nbi)[i] @ self.A_i_j @ q_j > p_i @ self.A_i_j @ q_j + tol:
            print('Pure strategy i =', i, 'beats p_i.')
            return False
    for j in range(self.nbj):
        if p_i @ self.B_i_j @ np.eye(self.nbj)[j] > p_i @ self.B_i_j @ q_j + tol:
            print('Pure strategy j =', j, 'beats q_j.')
            return False
    return True

Bimatrix_game.is_NashEq = Bimatrix_game_is_NashEq
In [11]:
sol = penalty_nonzero_sum.mangasarian_stone_solve()
p_i, q_j = sol['p_i'], sol['q_j']

penalty_nonzero_sum.is_NashEq(p_i, q_j)
Out[11]:
True

The Lemke–Howson algorithm¶

The Lemke–Howson algorithm was developed in 1964 to solve bimatrix games. Like the simplex, it is a path-following method. As we will see, it relies on the notion of complementarity.

Abridged system formulation¶

Recall our Nash equilibrium system, which we write here in matrix form:

\begin{align} p &\geq 0 \\ q &\geq 0 \\ 1_{\mathcal I}^\top p &= 1 \\ 1_{\mathcal J}^\top q &= 1 \\ \alpha 1_{\mathcal I} - A q &\geq 0 \\ \beta 1_{\mathcal J} - B^\top p &\geq 0 \\ p^\top \big( \alpha 1_{\mathcal I} - A q \big) &= 0 \\ q^\top \big( \beta 1_{\mathcal J} - B^\top p \big) &= 0. \end{align}

We know that any solution has $\alpha > 0$ and $\beta > 0$, so we can use the change of variables $x = p / \beta$ and $y = q / \alpha$ to obtain

\begin{align} x &\geq 0 \\ y &\geq 0 \\ 1_{\mathcal I}^\top x &= 1/\beta \\ 1_{\mathcal J}^\top y &= 1/\alpha \\ 1_{\mathcal I} - A y &\geq 0 \\ 1_{\mathcal J} - B^\top x &\geq 0 \\ x^\top \big( 1_{\mathcal I} - A y \big) &= 0 \\ y^\top \big( 1_{\mathcal J} - B^\top x \big) &= 0. \end{align}

The variables $x$ and $y$ are simply scaled versions of the players' mixed strategies: if we find a solution $(x,y)$ to this system, then we recover a Nash equilibrium with

$p = \frac{x}{\sum_i x_i}$ and $q = \frac{y}{\sum_j y_j}$.

(Note that this is the same normalization we saw in the zero-sum case!)

Since we can always adjust $\alpha$ and $\beta$ to verify $1/\beta = 1_{\mathcal I}^\top x$ and $1/\alpha = 1_{\mathcal J}^\top y$, we can actually drop these equations. Our system boils down to the abridged representation

\begin{align} x &\geq 0 \\ y &\geq 0 \\ 1_{\mathcal I} - A y &\geq 0 \\ 1_{\mathcal J} - B^\top x &\geq 0 \\ x^\top \big( 1_{\mathcal I} - A y \big) &= 0 \\ y^\top \big( 1_{\mathcal J} - B^\top x \big) &= 0. \end{align}

Observe that our abridged system has $x = 0$, $y = 0$ as a trivial solution. But this solution does not correspond to a Nash equilibrium! This is because the change of variable we performed above is not possible in this case.

This 'fake' solution, however, is good news: it will be our entry point to a path towards an actual solution.

Remark. The system above can be interpreted as a Linear Complementarity Problem (LCP):

$0 \leq x \perp 1_{\mathcal I} - A y \geq 0$

$0 \leq y \perp 1_{\mathcal J} - B^\top x \geq 0$.

This LCP has the trivial solution $(x,y) = (0,0)$. To avoid this, one may instead consider that $A > 0$ is the loss matrix of player 1, in which case the problem becomes

$0 \leq x \perp -1_I + A y \geq 0$
$0 \leq y \perp 1_J - B^\top x \geq 0$

which does not have a trivial solution anymore.

For more info, see the math+econ+code notebook on LCPs.

Lemke–Howson on an example¶

Let's take a simple example to illustrate:

\begin{equation} A = \begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}, \quad B = \begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}. \end{equation}

This is a battle-of-the-sexes game, which has 2 equilibria in pure strategies, and 1 in mixed strategies.

In [12]:
battle_game = Bimatrix_game(A_i_j = np.array([[3, 1], [1, 2]]),
                            B_i_j = np.array([[2, 1], [1, 3]]))

Step 1: Initialize.¶

To start with, let's introduce the slack variables $s$ and $t$ defined by $s_i = 1 - (A y)_i$ and $t_j = 1 - (B^\top x)_j$, i.e.

\begin{align} s_1 &= 1 & & &- 3 y_1 &- y_2 \\ s_2 &= 1 & & &- y_1 &- 2 y_2 \\ t_1 &= 1 &- 2 x_1 &- x_2 & & \\ t_2 &= 1 &- x_1 &- 3 x_2 & & \end{align}

With these notations, our problem becomes the following: find vectors $x, y, s, t \geq 0$ which satisfy the equalities above, as well as the complementarity conditions

\begin{align} x_i s_i &= 0 \quad \text{for all $i$} \\ y_j t_j &= 0 \quad \text{for all $j$}. \end{align}

As in the case of linear programming, we reason in terms of basic and non-basic variables.
Our initial solution is $(x=0, y=0)$, so our initial basic variables are $s_1, s_2, t_1, t_2$, and the non-basic variables are $x_1, x_2, y_1, y_2$.

Let's code this using the Dictionary method we built for LPs.

In [13]:
import mec.lp
from mec.lp import Dictionary
from sympy import *

battle_dict = Dictionary(slack_var_names_i = ['s_1', 's_2', 't_1', 't_2'],
                         decision_var_names_j = ['x_1', 'x_2', 'y_1', 'y_2'],
                         A_i_j = np.block([[np.zeros((battle_game.nbi, battle_game.nbi)),
                                            battle_game.A_i_j],
                                           [battle_game.B_i_j.T,
                                            np.zeros((battle_game.nbj, battle_game.nbj))]]),
                         b_i = np.array([1,1,1,1]))

battle_dict.display()

-------------------------- 
Objective and constraints:
s_1 = -3.0*y_1 - 1.0*y_2 + 1
s_2 = -1.0*y_1 - 2.0*y_2 + 1
t_1 = -2.0*x_1 - 1.0*x_2 + 1
t_2 = -1.0*x_1 - 3.0*x_2 + 1

Let's also list the variables which are complements:

In [14]:
def Dictionary_make_complements(self, verbose=0):
    comp_vars = {Symbol(name): Symbol(self.decision_var_names_j[i])
                  for (i,name) in enumerate(self.slack_var_names_i) }
    comp_vars.update( {Symbol(name): Symbol(self.slack_var_names_i[j])
                        for (j,name) in enumerate(self.decision_var_names_j) } )
    self.complements = comp_vars
    return

Dictionary.make_complements = Dictionary_make_complements

battle_dict.make_complements()
battle_dict.complements
Out[14]:
{s_1: x_1,
 s_2: x_2,
 t_1: y_1,
 t_2: y_2,
 x_1: s_1,
 x_2: s_2,
 y_1: t_1,
 y_2: t_2}

Step 2: Departing variable.¶

Now to start off the algorithm, we pick any non-basic variable to enter the basis. Let's choose $y_1$ for instance.

Since a variable entered the basis, another one needs to leave it. To determine the departing variable, we use the same method as in the simplex algorithm. When $y_1$ enters, we increase its value until we hit one of the constraints $s_i \geq 0$ or $t_j \geq 0$. Since $y_1$ only appears in the expressions of the $s$ variables, it will necessarily be one of those which leaves the basis:

\begin{align} s_1 &= 1 - 3 y_1 - y_2 \\ s_2 &= 1 - y_1 - 2 y_2. \end{align}

Here it is the constraint $s_1 \geq 0$ that we hit first, for $y_1 = 1/3$. Hence $s_1$ is our departing variable.

This is the minimum-ratio rule, and the good news is that we already coded this step in the simplex algorithm! We can simply reuse the function determine_departing we built for linear programs.

In [15]:
entering_var = Symbol('y_1')
departing_var = battle_dict.determine_departing(entering_var)
departing_var
Out[15]:
$\displaystyle s_{1}$

In general, as in the simplex algorithm, if $y_{j^\star}$ is the entering variable, then the departing variable $s_{i^\star}$ is determined with:

$i^\star = \arg\min_{\, i} \left\{ \frac {1} {A_{i \, j^\star}} \right\}$.

Now, you may remark that in making $y_1$ enter the basis, we violated one complementarity condition, namely $y_1 t_1 = 0$.

If $t_1$ had left the basis this wouldn't be a problem, but this is not the case (in fact it cannot be the case at the first step, since we saw that it was an $s$ variable which had to leave).

We say that the basis is now almost-complementary, because it violates exactly one complementarity condition.
This is a property we will maintain throughout the algorithm. And if we manage to reestablish full complementarity, then it means we found a solution!

First we update our tableau using the pivot method we built for LPs.

In [16]:
battle_dict.pivot(entering_var, departing_var, verbose=1)

battle_dict.display()

Entering = y_1; departing = s_1
-------------------------- 
Objective and constraints:
s_2 = 0.33*s_1 - 1.67*y_2 + 0.67
t_1 = -2.0*x_1 - 1.0*x_2 + 1
t_2 = -1.0*x_1 - 3.0*x_2 + 1
y_1 = -0.33*s_1 - 0.33*y_2 + 0.33

Step 3: Entering variable.¶

Recall that in the simplex algorithm, the entering variable was determined by updating the objective with the new non-basic variables, and taking one of the non-basic variables which had a positive coefficient in this updated objective.

Here there is no objective to optimize, but instead we will choose the entering variable using our complementarity conditions.

Here, since $s_1$ left the basis, the complementarity conditions allows $x_1$ to enter.

In [17]:
entering_var = battle_dict.complements[departing_var]
entering_var
Out[17]:
$\displaystyle x_{1}$

Now that we have an entering variable, we can go back to finding a departing variable using the minimum-ratio rule, etc.

Step 4: Stopping condition.¶

The algorithm stops when we have reestablished a complementary basis. (It is possible to show that this always happens.)

In [18]:
def Dictionary_is_basis_complementary(self, verbose=0):
    for var in self.base.keys():
        comp_var = self.complements[var]
        
        if comp_var in self.base.keys():
            if verbose > 0:
                print("Basis contains " + str(var) + " and " + str(comp_var))
            return False
    
    print("Complementary basis found!")
    print(self.base.keys())
    return True

Dictionary.is_basis_complementary = Dictionary_is_basis_complementary
In [19]:
battle_dict.is_basis_complementary(verbose=1)

Basis contains t_1 and y_1
Out[19]:
False

Now we can iterate these steps.

In [20]:
while not battle_dict.is_basis_complementary(verbose=1):
    departing_var = battle_dict.determine_departing(entering_var)
    battle_dict.pivot(entering_var, departing_var)
    entering_var = battle_dict.complements[departing_var]

Basis contains t_1 and y_1
Complementary basis found!
dict_keys([s_2, t_2, y_1, x_1])

Finally, we recover our strategies by normalizing the $x$ and $y$ we found:

In [21]:
z_sol, _ = battle_dict.solution()

x_sol = z_sol[:battle_game.nbi]
y_sol = z_sol[battle_game.nbi:(battle_game.nbi+battle_game.nbj)]

print("p =", x_sol/x_sol.sum())
print("q =", y_sol/y_sol.sum())

p = [1. 0.]
q = [1. 0.]

Let's compare with Mangasarian–Stone:

In [22]:
battle_game.mangasarian_stone_solve()
Out[22]:
{'p_i': array([0.66666667, 0.33333333]),
 'q_j': array([0.33333333, 0.66666667]),
 'val1': 1.6666666666666667,
 'val2': 1.6666666666666667}

Full algorithm¶

Now we can code the full algorithm and test it on the penalty game:

In [23]:
def Bimatrix_game_lemke_howson_solve(self, verbose = 0):
    dictionary = Dictionary(slack_var_names_i = ['s_'+str(i) for i in range(1,self.nbi+1)]
                                                + ['t_'+str(j) for j in range(1,self.nbj+1)],
                            decision_var_names_j= ['x_'+str(i) for i in range(1,self.nbi+1)]
                                                  + ['y_'+str(j) for j in range(1,self.nbj+1)],
                            A_i_j = np.block([[np.zeros((self.nbi, self.nbi)), self.A_i_j],
                                            [self.B_i_j.T, np.zeros((self.nbj, self.nbj))]]),
                            b_i = np.ones(self.nbi+self.nbj))
    dictionary.make_complements()
    entering_var = dictionary.nonbasic[0]
    departing_var = dictionary.determine_departing(entering_var)
    dictionary.pivot(entering_var, departing_var)
    entering_var = dictionary.complements[departing_var]
    
    while not dictionary.is_basis_complementary(verbose):
        departing_var = dictionary.determine_departing(entering_var)
        dictionary.pivot(entering_var, departing_var, verbose = 2)
        entering_var = dictionary.complements[departing_var]
    
    z_sol, _ = dictionary.solution()
    x_sol = z_sol[:self.nbi]
    y_sol = z_sol[self.nbi:(self.nbi+self.nbj)]
    p, q = x_sol/x_sol.sum(), y_sol/y_sol.sum()
    
    return p, q

Bimatrix_game.lemke_howson_solve = Bimatrix_game_lemke_howson_solve
In [24]:
penalty_nonzero_sum.lemke_howson_solve(verbose=1)

Basis contains s_1 and x_1
Entering = y_1; departing = s_3
y_1 = -0.01*s_3 - 1.03*y_2 - 0.78*y_3 + 0.01
Basis contains s_1 and x_1
Entering = x_3; departing = t_3
x_3 = -0.01*t_3 - 0.18*x_2 + 0.01
Basis contains s_1 and x_1
Entering = y_3; departing = s_2
y_3 = -0.06*s_2 + 0.06*s_3 + 3.18*y_2
Basis contains s_1 and x_1
Entering = x_2; departing = t_2
x_2 = 0.01 - 0.e-2*t_2
Basis contains s_1 and x_1
Entering = y_2; departing = s_1
y_2 = -0.e-2*s_1 + 0.02*s_2 - 0.01*s_3
Complementary basis found!
dict_keys([x_1, y_1, x_3, y_3, x_2, y_2])
Out[24]:
(array([0.3374337 , 0.24917907, 0.41338723]),
 array([0.21908541, 0.10161508, 0.67929951]))

Compare with the Mangasarian–Stone method:

In [25]:
penalty_nonzero_sum.mangasarian_stone_solve()
Out[25]:
{'p_i': array([0.3374337 , 0.24917907, 0.41338723]),
 'q_j': array([0.21908541, 0.10161508, 0.67929951]),
 'val1': 82.62606457928055,
 'val2': 36.37698012002112}

We also compare with the method we built to solve LCPs:

In [26]:
from mec.gt import LCP

def Bimatrix_game_to_LCP(self):
    M = 1 + self.A_i_j.max()
    M_i_j = np.block([[np.zeros((self.nbi, self.nbi)), M-self.A_i_j],
                      [self.B_i_j.T, np.zeros((self.nbj, self.nbj))]])
    q_i = np.concatenate((-np.ones(self.nbi), np.ones(self.nbj)))
    return LCP(M_i_j, q_i)

Bimatrix_game.to_LCP = Bimatrix_game_to_LCP
In [27]:
penalty_nonzero_sum.to_LCP().lemke_solve(verbose=1)

==========
Solution not found: Ray termination.

Since the Lemke–Howson algorithm is certain to find a solution, it is more robust than Lemke's algorithm to find solutions to bimatrix games.