First follow the lecture on NumPy.
Import NumPy:
In [1]:
import numpy as np
Consider the following 1-dimensional array:
In [2]:
n = 5
a = np.arange(n)
In [3]:
a
Out[3]:
Obtain the value in the last entry:
In [4]:
a[-1]
Out[4]:
Obtain the values in all entries except the last:
In [5]:
a[:-1]
Out[5]:
Assign -10 to the last entry:
In [6]:
a[:-1] = -10
In [7]:
a
Out[7]:
Next consider the following 2-dimensional array:
In [8]:
m, n = 3, 3
b = np.arange(m*n).reshape(m, n)
In [9]:
b
Out[9]:
Compute the matrix-vector mutiplication of b and [1/3, 1/3, 1/3]:
In [10]:
b @ [1/3, 1/3, 1/3]
Out[10]:
Let
In [11]:
I = [1, 2]
J = [0, 1]
We want to obtain the submatrix b_IJ of b with rows in I and columns in J.
Obtain the submatrix b_I of b with rows in I:
In [12]:
b_I = b[I, :]
In [13]:
b_I
Out[13]:
Obtain the submatrix b_IJ of b_I with columns in J:
In [14]:
b_IJ = b_I[:, J]
In [15]:
b_IJ
Out[15]:
Thus, we can obtain b_IJ from b by:
In [16]:
b_IJ = b[I, :][:, J]
In [17]:
b_IJ
Out[17]:
Assign 100 to the last column of b:
In [18]:
b[:, -1] = 100
In [19]:
b
Out[19]:
Assign -20 to all the entries except those in the last row and the last column:
In [20]:
b[:-1, :-1] = -20
In [21]:
b
Out[21]:
Consider the example from von Stengel (2007): $$ A = \begin{bmatrix} 3 & 3 \\ 2 & 5 \\ 0 & 6 \end{bmatrix}, \quad B = \begin{bmatrix} 3 & 2 \\ 2 & 6 \\ 3 & 1 \end{bmatrix}. $$
It will be convenient to denote the actions by indices, so let the action spaces of players 1 and 2 be $$ M = \{0, 1, 2\}, \quad N = \{0, 1\}. $$ (Note that in Python, indices start with 0.)
Create NumPy arrays for the payoff matrices:
In [22]:
A = np.array([[3, 3],
[2, 5],
[0 ,6]])
B = np.array([[3, 2],
[2, 6],
[3, 1]])
In [23]:
m, n = A.shape # Numbers of actions of players 1 and 2, respectively
M = np.arange(m)
N = np.arange(n)
It will be more convenient to work with the transpose of matrix B:
In [24]:
B_T = B.T
In [25]:
B_T
Out[25]:
For each $k = 1, \ldots, \min\{m, n\}$ and each support pair $(I, J)$ with $\lvert I\rvert = \lvert J\rvert = k$:
Solve the systems of linear equations: $$ C \begin{pmatrix} y_J \\ u \end{pmatrix} = e $$ and $$ D \begin{pmatrix} x_I \\ v \end{pmatrix} = e $$ with $$ C = \begin{pmatrix} A_{IJ} & -\mathbf{1} \\ \mathbf{1}' & 0 \end{pmatrix}, \quad D = \begin{pmatrix} B'_{JI} & -\mathbf{1} \\ \mathbf{1}' & 0 \end{pmatrix}, \quad e = \begin{pmatrix}\mathbf{0} \\ 1\end{pmatrix}, $$ where $A_{IJ}$ is the submatrix of $A$ given by rows $I$ and columns $J$, $B'_{JI}$ is the submatrix of $B'$ given by rows $J$ and columns $I$, and $\mathbf{0}$ and $\mathbf{1}$ are the $k$-dimensional vectors of zeros and ones, respectively.
Then check that $x_i > 0$ for all $i \in I$ and $y_j > 0$ for all $j \in J$, and that $u \geq (A_J y_J)_i$ for all $i \notin I$ and $v \geq (B'_I x_I)_j$ for all $j \notin J$, where $A_J$ is the submatrix of $A$ given by columns $J$ and $B'_I$ is the submatrix of $B'$ given by columns $I$.
As an example, let $I = \{0, 1\} \subset M$ and $J = \{0, 1\} \subset N$, and get the corresponding submatices of $A$ and $B'$.
In [26]:
k = 2
I = [0, 1]
J = [0, 1]
Submatrix $A_{IJ}$:
In [27]:
A_IJ = A[I, :][:, J]
In [28]:
A_IJ
Out[28]:
Submatrix $B'_{JI}$:
In [29]:
B_T_JI = B_T[J, :][:, I]
In [30]:
B_T_JI
Out[30]:
Construct the matrices $C$ and $D$ and the vector $e$,
and pass them the linear equation solver np.linalg.solve.
First construct the vector $e = (0, \ldots, 0, 1)$:
In [31]:
e = np.zeros(k+1) # Vector of zeros of length k+1
e[-1] = 1 # Assign 1 to the last entry of e
In [32]:
e
Out[32]:
Here's how to construct $C = \begin{pmatrix} A_{IJ} & -\mathbf{1} \\ \mathbf{1}' & 0 \end{pmatrix}$:
In [33]:
C = np.empty((k+1, k+1)) # Prepare an empty arry of shape (k+1, k+1)
C[:-1, :-1] = A_IJ # Assign A_IJ to the top left k x k submatrix
C[-1, :-1] = 1 # Assign 1 to the first k entries of the row at the bottom
C[:-1, -1] = -1 # Assign -1 to the first k entries of the right most column
C[-1, -1] = 0 # Assign 0 to the bottom right corner
In [34]:
C
Out[34]:
Now solve the equation $C \begin{pmatrix} y_J \\ u \end{pmatrix} = e$:
In [35]:
sol_1 = np.linalg.solve(C, e)
In [36]:
sol_1
Out[36]:
sol_1[:-1] corresponds to $y_J$, and sol_1[-1] to $u$:
In [37]:
y_J = sol_1[:-1]
y_J
Out[37]:
In [38]:
u = sol_1[-1]
u
Out[38]:
Check that $y_J \gg 0$:
In [39]:
(sol_1[:-1] > 0).all()
Out[39]:
One way to check that $u \geq (A_J y_J)_i$ for all $i \notin I$:
In [40]:
I_c = np.setdiff1d(M, I) # Complement of I in M
A[I_c, :][:, J] @ y_J # Vector given by (A_J y_J)_i for i in complement of I
Out[40]:
In [41]:
(u >= A[I_c, :][:, J] @ y_J).all()
Out[41]:
So this is a candidate for the mixed action for player 2 in an equilibrium:
In [42]:
y = np.zeros(n)
y[J] = y_J
y
Out[42]:
(y is the same as y_J because $J = N$,
but let's make the code general.)
Do the same exercise for player 2:
Construct $D = \begin{pmatrix} B'_{JI} & -\mathbf{1} \\ \mathbf{1}' & 0 \end{pmatrix}$:
In [43]:
D = np.empty((k+1, k+1)) # Prepare an empty arry of shape (k+1, k+1)
D[:-1, :-1] = B_T_JI # Assign A_IJ to the top left k x k submatrix
D[-1, :-1] = 1 # Assign 1 to the first k entries of the row at the bottom
D[:-1, -1] = -1 # Assign -1 to the first k entries of the right most column
D[-1, -1] = 0 # Assign 0 to the bottom right corner
In [44]:
D
Out[44]:
Solve $D \begin{pmatrix} x_I \\ v \end{pmatrix} = e$:
In [45]:
sol_2 = np.linalg.solve(D, e)
sol_2
Out[45]:
Check that $x_I \gg 0$:
In [46]:
(sol_2[:-1] > 0).all()
Out[46]:
Check that $v \geq (B'_I x_I)_j$ for all $j \notin J$:
In [47]:
J_c = np.setdiff1d(N, J)
(sol_2[-1] >= B[J_c, :][:, I] @ sol_2[:-1]).all()
Out[47]:
So this is an equilibrium mixed action for player 1:
In [48]:
x = np.zeros(m)
x[I] = sol_2[:-1]
x
Out[48]:
We have found a Nash equilibrium:
In [49]:
(x, y)
Out[49]:
Let's do the same exercise for $I = \{0, 2\}$ and $J = \{0, 1\}$.
In [50]:
k = 2
I = [0, 2]
J = [0, 1]
To get the matrix C, let's use the previous C,
just replacing the submatrix in C[:-1, :-1]:
In [51]:
C[:-1, :-1] = A[I, :][:, J]
In [52]:
C
Out[52]:
In [53]:
sol_1 = np.linalg.solve(C, e)
In [54]:
sol_1
Out[54]:
Check the positivity:
In [55]:
(sol_1[:-1] > 0).all()
Out[55]:
Check the optimality:
In [56]:
I_c = np.setdiff1d(M, I)
(sol_1[-1] >= A[I_c, :][:, J] @ sol_1[:-1]).all()
Out[56]:
This means that $v \neq \max_i (A_J y_J)_i$.
So this implies that there is no Nash equilibrium for the support pair $(I, J)$. But let's just check also the conditions for player 2.
In [57]:
D[:-1, :-1] = B_T[J, :][:, I]
In [58]:
D
Out[58]:
In [59]:
sol_2 = np.linalg.solve(D, e)
sol_2
Out[59]:
In [60]:
(sol_2[:-1] > 0).all()
Out[60]:
This means that the solution cannot be a probability distribution.
Let's write a function that wraps the above procedure.
The following is just one possible implementation design:
payoff_matrix;own_supp for the support of the player in consideration;opp_supp for the support of the opponent player;out that stores the candidate mixed action.opp_supp
against which the actions in own_supp are best responses,
then return True; otherwise return False.
In the former case, the mixed action is stored in out.out must be of length equal to the number of the opponent's actions.
In [61]:
def indiff_mixed_action(payoff_matrix, own_supp, opp_supp, out):
# (number of own actions, number of opponent's actions)
nums_actions = payoff_matrix.shape
# Support size
k = len(own_supp)
# Matrix in the left hand side of the linear equation
a = np.empty((k+1, k+1))
a[:-1, :-1] = payoff_matrix[own_supp, :][:, opp_supp]
a[-1, :-1] = 1
a[:-1, -1] = -1
a[-1, -1] = 0
# Vector in the right hand side of the linear equation
b = np.zeros(k+1)
b[-1] = 1
try:
sol = np.linalg.solve(a, b)
except np.linalg.LinAlgError:
return False
# Return False immediately if any of the "probabilities" is not positive
if (sol[:-1] <= 0).any():
return False
own_supp_c = np.setdiff1d(np.arange(nums_actions[0]), own_supp)
# Return False immediately if the solution mixed action is not optimal
if (sol[-1] < payoff_matrix[own_supp_c, :][:, opp_supp] @ sol[:-1]).any():
return False
out.fill(0)
out[opp_supp] = sol[:-1]
return True
In [62]:
I = [0, 1]
J = [0, 1]
In [63]:
out = np.empty(n)
indiff_mixed_action(A, I, J, out)
Out[63]:
In [64]:
out
Out[64]:
In [65]:
out = np.empty(m)
indiff_mixed_action(B_T, J, I, out)
Out[65]:
In [66]:
out
Out[66]:
In [67]:
I = [0, 2]
J = [0, 1]
In [68]:
out = np.empty(n)
indiff_mixed_action(A, I, J, out)
Out[68]:
In [69]:
out = np.empty(m)
indiff_mixed_action(B_T, J, I, out)
Out[69]:
Note on the try-except block:
There are cases where the coeffiecient matrix in the linear equation is singular,
in which case np.linalg.solve(a, b) raises an error np.linalg.LinAlgError.
For example:
In [70]:
np.linalg.solve(np.zeros((2, 2)), np.ones(2))
If this happens, the program stops.
But singularity just implies that there is no Nash equilibrium for the given support pair,
so we want to just return False without stopping the program.
A try-except statement can be used for this purpose;
see this section
in the Python tutorial.
In the final step,
we need a loop interating over all the equal-size support pairs
in which indiff_mixed_action is applied to each pair.
One approach is to use itertools.combinations:
In [71]:
import itertools
In [72]:
m, n = 3, 3
k = 2
for I in itertools.combinations(range(m), k):
for J in itertools.combinations(range(n), k):
print((I, J))
In [73]:
def support_enumeration(A, B_T):
"""
Given a nondegenerate m x n bimatrix game (A, B_T), return a list of
all Nash equilibria computed by the support enumeration algorithm.
Parameters
----------
A : ndarray(float, ndim=2)
Payoff matrix for player 1, of shape (m, n).
B_T : ndarray(float, ndim=2)
Payoff matrix for player 2, of shape (n, m).
Returns
-------
NEs : list(tuple(ndarray(float, ndim=1)))
List containing tuples of Nash equilibrium mixed actions.
"""
m, n = A.shape
NEs = []
for k in range(1, min(m, n)+1):
for I in itertools.combinations(range(m), k):
for J in itertools.combinations(range(n), k):
y = np.empty(n)
if indiff_mixed_action(A, list(I), list(J), y):
x = np.empty(m)
if indiff_mixed_action(B_T, list(J), list(I), x):
NEs.append((x, y))
return NEs
In [74]:
support_enumeration(A, B_T)
Out[74]:
""" are called "docstrings",
which offer a description about how the function works.
It is a good practice to always write docstrings, which often times helps yourself in the future;
it is very likely that you will forget what you have written in the past.
(I should provide docstrings for the function indiff_mixed_action, too.)
Refer to the docstring guideline from QuantEcon.py.indiff_mixed_action are appended to the list NEs,
which starts with an empty list [].itertools.combinations yields tuples,
so that they have to be converted to lists by list(I) and list(J).Consider the $6 \times 6$ game from von Stengel (1997) "New Lower Bounds for the Number of Equilibria in Bimatrix Games," page 12:
In [75]:
A = np.array(
[[ 9504, -660, 19976, -20526, 1776, -8976],
[ -111771, 31680, -130944, 168124, -8514, 52764],
[ 397584, -113850, 451176, -586476, 29216, -178761],
[ 171204, -45936, 208626, -263076, 14124, -84436],
[ 1303104, -453420, 1227336, -1718376, 72336, -461736],
[ 737154, -227040, 774576, -1039236, 48081, -300036]]
)
B_T = np.array(
[[ 72336, -461736, 1227336, -1718376, 1303104, -453420],
[ 48081, -300036, 774576, -1039236, 737154, -227040],
[ 29216, -178761, 451176, -586476, 397584, -113850],
[ 14124, -84436, 208626, -263076, 171204, -45936],
[ 1776, -8976, 19976, -20526, 9504, -660],
[ -8514, 52764, -130944, 168124, -111771, 31680]]
)
Count the number of Nash equilibria of this game.
In [76]:
len(support_enumeration(A, B_T))
Out[76]:
Background:
The $n \times n$ game where the payoff matrices are given by the identy matrix has $2^n - 1$ equilibria. It had been conjectured by Quint and Shubik (1997) that this is the maximum number of equilibria of any nondegenerate $n \times n$ game. The game above was presented as a counter-example to this conjecture.
In [ ]: