Scarf's ordinal basis algorithm

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

'math+econ+code' masterclass series

With python code examples

© 2018–2024 by Alfred Galichon. Research assistance from Alessandro Facchini is gratefully acknowledged. 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 Antoine Jacquet, 'Scarf's ordinal basis algorithm', 'math+econ+code' masterclass series. https://www.math-econ-code.org/

References

• Scarf (1967). "The core of an N person game." Econometrica.
• Biró, Fleiner and Irving (2016). "Matching couples with Scarf's algorithm." Annals of Mathematics and Artificial Intelligence.
• Nguyen and Vohra (2018). "Near-feasible stable matchings with couples." American Economic Review.

#!pip install mec --upgrade

import numpy as np
from mec.lp import Tableau

Motivation: core of a cooperative game

Description of the problem

Consider the game from Scarf (1967), section 5. This is a cooperative game with three players $z \in \mathcal Z = \{1,2,3\}$. Players can remain in autarky and get a payoff of zero, or they can sign arrangements between two players.

(Scarf ignores the possible arrangements between all three players because...)

Players 1 and 2 can sign three different arrangements. In the first one their payoffs are $(12,6)$ (meaning that 1 gets 12 and 2 gets 6), in the second one, $(3,7)$, and in the third one, $(2,9)$. Players 1 and 3 can also sign three possible arrangements with respective payoffs $(9,3)$, $(5,8)$, and $(4,10)$. Lastly, three arrangements are also available to players 2 and 3, yielding respectively $(5,6)$, $(2,9)$ and $(8,4)$. Let $\mathcal A$ be the set of all arrangements, which also includes autarky as arrangements of a single player with payoff 0.

We build a matrix whose rows are the players and columns are the arrangements. Given player $z \in \mathcal Z$ and arrangement $a \in \mathcal A$, we denote $\Phi_{za}$ the payoff of a player $z$ in arrangement $a$. If $z$ is not part of arrangement $a$, the value of $\Phi_{za}$ is not relevant, but it will be useful to define it as a very large number $K_a$ which will differ across arrangements.

The matrix $\Phi$ is therefore defined as \begin{equation} \Phi = \begin{bmatrix} 0 & K_1 & K_2 & 12 & 3 & 2 & 9 & 5 & 4 & K_9 & K_{10} & K_{11} \\ K_0 & 0 & K_2 & 6 & 7 & 9 & K_6 & K_7 & K_8 & 5 & 2 & 8 \\ K_0 & K_1 & 0 & K_3 & K_4 & K_5 & 3 & 8 & 10 & 6 & 9 & 4 \\ \end{bmatrix} \end{equation} with $K_0 > K_1 > \dots > K_{11}$.

Next, we define $M$ the margining-out matrix (MOM), which is a $|\mathcal Z| \times |\mathcal A|$ matrix such that $M_{za}$ is 1 if $z$ appears in arrangement $a$, and 0 otherwise. In the present example, we have:

\begin{equation} M = \begin{bmatrix} 1 \quad & 0 \quad & 0 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \quad & 0 \quad & 0 \quad & 0 \\ 0 \quad & 1 \quad & 0 \quad & 1 \quad & 1 \quad & 1 \quad & 0 \quad & 0 \quad & 0 \quad & 1 \quad & 1 \quad & 1 \\ 0 \quad & 0 \quad & 1 \quad & 0 \quad & 0 \quad & 0 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \quad & 1 \\ \end{bmatrix}. \end{equation}

Finally, we define $q$ a vector of ones, the same size as the number of players.

Remark. We could reinterpret this with a large number of individuals divided into a finite number of types $z \in \mathcal Z$. Then, $q_z$ would be the mass of individuals of type $z$.

We set this up as follows:

class OrdinalBasis:
    def __init__(self, Φ_z_a, M_z_a, q_z=None, K=None, eps=1e-5):
        if not (Φ_z_a.shape == M_z_a.shape):
            raise ValueError('Φ_z_a and M_z_a must have the same size.')
        self.Φ_z_a, self.M_z_a = Φ_z_a, M_z_a
        if K is None:
            K = Φ_z_a.max()
        self.nbstep, self.K, self.eps = 0, K, eps
        self.nbz, self.nba = self.Φ_z_a.shape
        if q_z is None:
            q_z = np.ones(self.nbz)
        self.q_z = q_z

def create_scarf_ex_sec5():
    K = 100
    Φ_z_a = 1.0 * np.array([[0,K,K,12,3,2,9,5, 4,K,K,K],
                            [K,0,K, 6,7,9,K,K, K,5,2,8],
                            [K,K,0, K,K,K,3,8,10,6,9,4]])
    M_z_a = (Φ_z_a != K) * 1.
    return (Φ_z_a, M_z_a)

scarf_ex_sec5 = OrdinalBasis(*create_scarf_ex_sec5())

print(scarf_ex_sec5.Φ_z_a)
print(scarf_ex_sec5.M_z_a)

[[  0. 100. 100.  12.   3.   2.   9.   5.   4. 100. 100. 100.]
 [100.   0. 100.   6.   7.   9. 100. 100. 100.   5.   2.   8.]
 [100. 100.   0. 100. 100. 100.   3.   8.  10.   6.   9.   4.]]
[[1. 0. 0. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 1. 0. 1. 1. 1. 0. 0. 0. 1. 1. 1.]
 [0. 0. 1. 0. 0. 0. 1. 1. 1. 1. 1. 1.]]

def OrdinalBasis_remove_degeneracies(self):
    self.Φ_z_a += np.arange(self.nba,0,-1) * (self.Φ_z_a == self.K)
    self.q_z = self.q_z + np.arange(1,self.nbz+1)*self.eps

OrdinalBasis.remove_degeneracies = OrdinalBasis_remove_degeneracies

scarf_ex_sec5.remove_degeneracies()
print(scarf_ex_sec5.Φ_z_a)
print(scarf_ex_sec5.q_z)

[[  0. 111. 110.  12.   3.   2.   9.   5.   4. 103. 102. 101.]
 [112.   0. 110.   6.   7.   9. 106. 105. 104.   5.   2.   8.]
 [112. 111.   0. 109. 108. 107.   3.   8.  10.   6.   9.   4.]]
[1.00001 1.00002 1.00003]

Feasible and stable matchings

We wish to find a matching $\mu = (\mu_a)_{a \in \mathcal A}$ which is both feasible and stable. The matching $\mu$ is feasible if

\begin{equation} M \mu = q. \end{equation}

The matching $\mu$ is stable if there is no arrangement $a$ which would make all its members strictly better off than they currently are. If $\mu$ is integral (i.e. if its elements are either 0 or 1) then each player belongs to a single arrangement and his payoff is therefore

\begin{equation} p_z = \min_{a | \mu_a > 0} \Phi_{za}. \end{equation}

This is by construction, since we have set $\Phi_{za} = K_a$ large enough whenever $z$ is not part of $a$.
Hence the stability condition can be written

\begin{equation} \not\exists a, \forall z, \; p_z < \Phi_{za}, \qquad \text{or equivalently,} \qquad \forall a, \exists z, \; p_z \geq \Phi_{za}. \end{equation}

If instead we interpret $z$ as types, then a matching does not have to be integral, and two individuals of the same type can have different payoffs depending on which arrangement they end up assigned to. In this case, $p_z$ is the worst-case utility for an individual of type $z$. Still, an arrangement $a$ will block $\mu$ if and only if $p_z < \Phi_{za}$ for all $z$. Indeed, if the match can be blocked by individuals who are not the worse-off within their type (so with payoff more than $p_z$), then a fortiori the match can also be blocked by those who are the worse-off in their type.

Lastly, what if we really have $N$ players and $\mu$ is not integral? How do we understand the matching then?

Scarf's algorithm

Overview

Scarf's lemma states that a matching $\mu$ which is both feasible and stable always exists (although it may not be integral). Scarf's ordinal-basis algorithm finds such a matching.

We provide a rough description here — the operations and terms will be defined precisely below.

---

Scarf's algorithm

Step 0: Initialization.
Let $n = 0$. Initialize the $M$ basis to the autarky columns $\{0,1,2\}$ and the $\Phi$ basis to columns $\{1,2,a_0\}$, where $a_0 \geq 3$ is chosen to form an ordinal basis.

Step 1: Cardinal pivot.
Increment $n$. Column $a_{n-1}$ enters the $M$ basis. Determine the column $a'_n$ which departs from the $M$ basis via the cardinal pivot.

Step 2: Ordinal pivot.
Column $a'_n$ also departs from the $\Phi$ basis. Determine the column $a_n$ which enters the $\Phi$ basis via the ordinal pivot .

Step 3: Check if bases coincide.
If the $M$ basis coincides with the $\Phi$ basis, go to step 4; if not, go back to step 1.

Step 4: Find the solution.
Solve $M \mu = q$ using only the columns of $M$ corresponding to the $M$ basis.

---

Let's implement this algorithm.

Initialization

We start with the basis formed by the first three columns of $M$, $\{0,1,2\}$, corresponding to the three autarky arrangements.

def OrdinalBasis_init_basis_M(self):
    self.tableau_M = Tableau(A_i_j=self.M_z_a[:,self.nbz:self.nba], b_i=self.q_z)
    self.basis_M = list(self.tableau_M.k_b)   # k_b are the indices of columns associated with basic variables
    return 
    
OrdinalBasis.init_basis_M = OrdinalBasis_init_basis_M

scarf_ex_sec5.init_basis_M()
scarf_ex_sec5.tableau_M.display()
scarf_ex_sec5.basis_M

+-----+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+---------+
|     |   s_0 |   s_1 |   s_2 |   x_0 |   x_1 |   x_2 |   x_3 |   x_4 |   x_5 |   x_6 |   x_7 |   x_8 |     RHS |
+=====+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=========+
| s_0 |     1 |     0 |     0 |     1 |     1 |     1 |     1 |     1 |     1 |     0 |     0 |     0 | 1.00001 |
+-----+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+---------+
| s_1 |     0 |     1 |     0 |     1 |     1 |     1 |     0 |     0 |     0 |     1 |     1 |     1 | 1.00002 |
+-----+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+---------+
| s_2 |     0 |     0 |     1 |     0 |     0 |     0 |     1 |     1 |     1 |     1 |     1 |     1 | 1.00003 |
+-----+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+---------+

[0, 1, 2]

We verify that this basis yields a feasible solution to the equation $M \mu = q$.

def OrdinalBasis_μsol_a(self, basis=None):
    if basis is None:
        basis = self.basis_M
    B = self.M_z_a[:,list(basis)]
    μ_a = np.zeros(self.nba)
    μ_a[list(basis)] = np.linalg.solve(B,self.q_z)
    return μ_a
    
OrdinalBasis.μsol_a = OrdinalBasis_μsol_a

scarf_ex_sec5.μsol_a([0,1,2])

array([1.00001, 1.00002, 1.00003, 0.     , 0.     , 0.     , 0.     ,
       0.     , 0.     , 0.     , 0.     , 0.     ])

In this solution, each player is assigned to his autarky arrangement.
Next, we check if that matching is stable. For this we need to define $p_z$, the worst-case payoff of each player $z$ in the basis $\{0,1,2\}$.

def OrdinalBasis_p_z(self, basis=None):
    if basis is None:
        basis = self.basis_Φ
    return self.Φ_z_a[:,list(basis)].min(axis=1)

OrdinalBasis.p_z = OrdinalBasis_p_z
scarf_ex_sec5.p_z([0,1,2])

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

This makes sense – in autarky, every player gets zero. Now we would like to check if there are blocking arrangements, that is, if there are columns $a$ such that $p_z < \Phi_{za}$ for every $z$. If such columns don't exist, the basis $\{0,1,2\}$ will be called an ordinal basis.

def OrdinalBasis_is_ordinal_basis(self, basis):
    ans, which = False, None
    if len(basis)==self.nbz:
        blocking = (self.Φ_z_a[:,list(basis)].min(axis = 1)[:,None] < self.Φ_z_a).all(axis = 0)
        if blocking.any():
            which = np.where(blocking)
        else:
            ans = True
    return ans, which

OrdinalBasis.is_ordinal_basis = OrdinalBasis_is_ordinal_basis
scarf_ex_sec5.is_ordinal_basis([0,1,2])

(False, (array([ 3,  4,  5,  6,  7,  8,  9, 10, 11]),))

One sees that $\{0,1,2\}$ is not an ordinal basis — far from it in fact, since every other column corresponds to a blocking arrangement.

We pick up an arbitrary arrangement of the current basis, for example column 0, and we wonder if it is possible to exchange it against another column $a_0 \geq 3$ so that $\{1,2, a_0\}$ is an ordinal basis. In fact, we don't have much choice, as we shall have $p_z = 0$ for $z \neq 1$ irrespective of $a_0$. Therefore we need to choose $a_0$ such that

$p_1 = \min \{ \Phi_{11}, \Phi_{12}, \Phi_{1 a_0} \} \geq \Phi_{1 a}$ for all $a \notin \{1,2,a_0\}$.

But since $\min \{ \Phi_{11}, \Phi_{12}, \Phi_{1 a_0} \} = \Phi_{1 a_0}$, this yields

$\Phi_{1 a_0} \geq \Phi_{1 a}$ for all $a \notin \{1,2,a_0\}$, hence $a_0 = \arg\max_{a \notin \{0,1,2\}} ~ \Phi_{1a}$.

def OrdinalBasis_init_basis_Φ(self, a_departing):
    self.basis_Φ = list(range(self.nbz))
    self.basis_Φ.remove(a_departing)
    a_entering = self.nbz + self.Φ_z_a[a_departing,self.nbz:].argmax()
    self.basis_Φ.append(a_entering)
    self.entcol = a_entering
    return a_entering

OrdinalBasis.init_basis_Φ = OrdinalBasis_init_basis_Φ

We run:

entcol = scarf_ex_sec5.init_basis_Φ(0)
print('Column ' + str(entcol) + ' is entering the Φ basis.')
print('Φ basis:', scarf_ex_sec5.basis_Φ)
print('M basis:', scarf_ex_sec5.basis_M)

Column 9 is entering the Φ basis.
Φ basis: [1, 2, 9]
M basis: [0, 1, 2]

Cardinal pivot

Column 9 has entered the $\Phi$ basis. As a result, we shall want column 9 to enter the $M$ basis as well.
The standard pivot step on matrix $M$ will determine which column will depart from the $M$ basis.

def OrdinalBasis_determine_departing(self, a_entering):
    a_departing = self.tableau_M.determine_departing(a_entering)
    self.tableau_M.pivot(a_entering,a_departing)
    self.basis_M.remove(a_departing)
    self.basis_M.append(a_entering)
    return a_departing

OrdinalBasis.determine_departing = OrdinalBasis_determine_departing

depcol = scarf_ex_sec5.determine_departing(entcol)
print('Column ' + str(depcol) + ' is leaving the M basis.')
print('Φ basis:', scarf_ex_sec5.basis_Φ)
print('M basis:', scarf_ex_sec5.basis_M)

Column 1 is leaving the M basis.
Φ basis: [1, 2, 9]
M basis: [0, 2, 9]

Ordinal pivot

Column 1 has left the $M$ basis. As the two bases still differ, we need to remove column 1 from the $\Phi$ basis, and to determine which column enters instead to ensure that it remains an ordinal basis.

We begin by considering the submatrix of $\Phi$ consisting of the columns associated with the current $\Phi$ basis, $\{1,2,9\}$:

\begin{equation} \begin{bmatrix} K_1 & K_2 & K_9 \\ 0 & K_2 & 5 \\ K_1 & 0 & 6 \\ \end{bmatrix}. \end{equation}

The associated vector of payoffs is $p = (K_9, 0, 0)$.

Note that by construction, each of these columns has one, and only one, row minimizer (indeed, when we brought in column 9 it verified $\min(\Phi_{10}, \Phi_{11}, \Phi_{12}) \geq \Phi_{19}$). Hence, when we remove column 1, there will be exactly one column which will have two row minimizers.
Here it is column 9: its row minimizers are $K_9$, which was already the row minimizer for the previous basis; and $5$, which is the new row minimizer after removing column 1. We call $z^*$ the row associated with the old row minimizer: $z^* = 1$.

We now consider all columns $a$ such that $\Phi_{za} > \min (\Phi_{z2}, \Phi_{z9})$ for all $z \neq z^*$.
In our example, those are all the columns $a$ such that $\Phi_{1a} > 5$ and $\Phi_{3a} > 0$: the eligible columns are $A = \{3, 4, 5, 6, 7, 8, 11\}$.

Finally, we choose the entering column as that column in $A$ which maximizes $\Phi_{z^* a}$:

$a^* = \arg\max_{a \in A} \Phi_{z^* a}$.

For us, this is column 11 (with $\Phi_{1, 11} = K_{11}$).

Let's prove formally the following:

Proposition. The ordinal pivot maintains an ordinal basis with each column containing exactly one row minimizer.

Proof. Assume we start from an ordinal basis $B \subset \mathcal A$. The vector of payoffs $p$ is defined by $p_z = \min_{a \in B} \Phi_{za}$, and it verifies: $\forall a, \exists z, \; p_z \geq \Phi_{za}$. Furthermore, each $a \in B$ has a unique row minimizer.

Let $a_0 \in B$ be the column which is leaving the basis, and denote $z_0$ the row which was minimized at $a_0$. Now row $z_0$ is minimized at some other column $a_1 \in B - \{a_0\}$, which already minimized a row $z_1 \neq z_0$ beforehand.

To choose an entering column $a_2$, we restrict attention to candidates in $A = \{a \notin B : \forall z \neq z_1, \; \Phi_{za} > \min_{\tilde a \in B - \{a_0\}} \Phi_{z \tilde a} \}$.
We need to prove that such candidates exist. The column $a$ corresponding to $z_1$ in autarky is actually always a suitable candidate (by construction, the corresponding $K_a$ is among the $|\mathcal Z|-1$ largest numbers of any row $z \neq z_1$) – that is, except if it is already in the basis. But if that is the case, then $a_1$ must be this column (because $a_1$ achieves the minimum of row $z_1$ in the full matrix $\Phi$). Since $a_1$ is also a row minimizer for $z_0$ in $B - \{a_0\}$, it must be that $B - \{a_0\}$ corresponds to the first $|\mathcal Z| - 1$ columns of $\Phi$, with $a_1$ being number $|\mathcal Z| - 1$ (and the column where $z_0$ is in autarky as number $|\mathcal Z|$).
However, the only way this can occur is if $B = \{0, 1, 2, \dots, |\mathcal Z| - 2, a_0\}$

Note that the departing column $a_0$ cannot be a suitable candidate since for $z = z_0$, $\Phi_{z_0 a_0} < \Phi_{z_0 a_1}$.

Finally, we choose $a_2 = \arg\max_{a \in A} \Phi_{z_1 a}$. We necessarily have $\Phi_{z_1 a_2} \leq p_{z_1}$, otherwise $B$ would not be an ordinal basis. Since there are no degeneracies, we even have $\Phi_{z_1 a_2} < p_{z_1}$.
We define $B' = B - \{a_0\} + \{a_2\}$ and $p_z'$ the corresponding vector of payoffs. Then:
$p_{z_0}' = \Phi_{z_0 a_1} > p_{z_0}$,
$p_{z_1}' = \Phi_{z_1 a_2} < p_{z_1}$,
$p_z' = p_z$ for all $z \notin \{z_0, z_1\}$.
Thus each column still has a single row minimizer in $B'$.

Now we want to prove that for all $a$, there is some $z$ such that $p_z' \geq \Phi_{za}$. We have four cases.

(i) If $a \in B'$, then it minimizes one row among the columns of the basis, i.e. there is some $z$ such that $p_z' = \Phi_{za}$.
(ii) If $a = a_0$, then $p_{z_0}' > p_{z_0} = \Phi_{z_0 a_0}$.
(iii) If $a \in A$, then $p_{z_1}' = \Phi_{z_1 a_2} \geq \Phi_{z_1 a}$ by definition of $a_2$ as the maximizer.
(iv) Otherwise, $a$ is neither in $B$ nor in $A$. Hence by definition of $A$, there is some $z \neq z_1$ such that $\Phi_{za} \leq \min_{\tilde a \in B - \{a_0\}} \Phi_{z \tilde a}$. But since $z \neq z_1$, $p_z' = \min_{\tilde a \in B - \{a_0\}} \Phi_{z \tilde a}$ (only $z_1$ reaches the minimum at $a = a_2$ which is not in $B - \{a_0\}$) hence $p_z' \geq \Phi_{za}$.

We code this as follows:

def OrdinalBasis_determine_entering(self, a_departing):
    pbefore_z = self.p_z(self.basis_Φ)
    self.basis_Φ.remove(a_departing)
    pafter_z = self.p_z(self.basis_Φ)
    z0 = np.where(pbefore_z < pafter_z)[0][0]
    a0 = self.basis_Φ[np.argmin(self.Φ_z_a[z0,self.basis_Φ])]
    zstar = np.where(pbefore_z == self.Φ_z_a[:,a0])[0][0]
    eligible_columns = np.where(np.all(np.delete(self.Φ_z_a, zstar, axis=0) > np.delete(pafter_z, zstar)[:, None], axis=0))[0]
    a_entering = eligible_columns[np.argmax(self.Φ_z_a[zstar,eligible_columns])]
    self.basis_Φ.append(a_entering)
    return a_entering

OrdinalBasis.determine_entering = OrdinalBasis_determine_entering

We run:

entcol = scarf_ex_sec5.determine_entering(depcol)
print('Column ' + str(entcol) + ' is entering the Φ basis.')
print('Φ basis:', scarf_ex_sec5.basis_Φ)
print('M basis:', scarf_ex_sec5.basis_M)

Column 11 is entering the Φ basis.
Φ basis: [2, 9, 11]
M basis: [0, 2, 9]

Implementation of the algorithm

We now encode a full step of the algorithm.

def OrdinalBasis_step(self, a_entering, verbose=0):
    self.nbstep += 1
    a_departing = self.determine_departing(a_entering)

    if self.basis_M == self.basis_Φ:
        if verbose>0:
            print('Solution found in '+ str(self.nbstep) +' steps.\nBasis:', self.basis_Φ)
        return False
    else:
        new_entcol = self.determine_entering(a_departing)
        if verbose>1:
            print('-- Step ' + str(self.nbstep) + ' --')
            print('Column ' + str(a_entering) + ' enters the M basis.')
            print('Column ' + str(a_departing) + ' leaves the two bases.')
            print('Column ' + str(new_entcol) + ' enters the Φ basis.')
            print('M basis: ', self.basis_M)
            print('Φ basis: ', self.basis_Φ)
            print('p_z =', self.p_z(self.basis_Φ))
            
    return new_entcol

OrdinalBasis.step = OrdinalBasis_step

entcol = scarf_ex_sec5.step(entcol,verbose=2)

-- Step 1 --
Column 11 enters the M basis.
Column 9 leaves the two bases.
Column 6 enters the Φ basis.
M basis:  [0, 2, 11]
Φ basis:  [2, 11, 6]
p_z = [9. 8. 0.]

entcol = scarf_ex_sec5.step(entcol,verbose=2)

-- Step 2 --
Column 6 enters the M basis.
Column 2 leaves the two bases.
Column 7 enters the Φ basis.
M basis:  [0, 11, 6]
Φ basis:  [11, 6, 7]
p_z = [5. 8. 3.]

entcol = scarf_ex_sec5.step(entcol,verbose=2)

-- Step 3 --
Column 7 enters the M basis.
Column 6 leaves the two bases.
Column 3 enters the Φ basis.
M basis:  [0, 11, 7]
Φ basis:  [11, 7, 3]
p_z = [5. 6. 4.]

entcol = scarf_ex_sec5.step(entcol,verbose=2)

Solution found in 4 steps.
Basis: [11, 7, 3]

np.around(scarf_ex_sec5.μsol_a(),2)

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

To summarize:

def OrdinalBasis_solve(self, a_departing=0, verbose=0):
    Φ_z_a_store, q_z_store = self.Φ_z_a.copy(), self.q_z.copy()
    self.remove_degeneracies()
    self.init_basis_M()
    a_entering = self.init_basis_Φ(a_departing)
    if verbose>1:
        print('M basis: ', self.basis_M)
        print('Φ basis: ', self.basis_Φ)
    while a_entering:
        a_entering = self.step(a_entering, verbose)
    self.Φ_z_a, self.q_z = Φ_z_a_store, q_z_store
    return {'basis': self.basis_Φ, 'μ_a': self.μsol_a(), 'p_z': self.p_z()}

OrdinalBasis.solve = OrdinalBasis_solve

scarf_ex_sec5 = OrdinalBasis(*create_scarf_ex_sec5())
scarf_ex_sec5.solve(verbose=2)

M basis:  [0, 1, 2]
Φ basis:  [1, 2, 9]
-- Step 1 --
Column 9 enters the M basis.
Column 1 leaves the two bases.
Column 11 enters the Φ basis.
M basis:  [0, 2, 9]
Φ basis:  [2, 9, 11]
p_z = [101.   5.   0.]
-- Step 2 --
Column 11 enters the M basis.
Column 9 leaves the two bases.
Column 6 enters the Φ basis.
M basis:  [0, 2, 11]
Φ basis:  [2, 11, 6]
p_z = [9. 8. 0.]
-- Step 3 --
Column 6 enters the M basis.
Column 2 leaves the two bases.
Column 7 enters the Φ basis.
M basis:  [0, 11, 6]
Φ basis:  [11, 6, 7]
p_z = [5. 8. 3.]
-- Step 4 --
Column 7 enters the M basis.
Column 6 leaves the two bases.
Column 3 enters the Φ basis.
M basis:  [0, 11, 7]
Φ basis:  [11, 7, 3]
p_z = [5. 6. 4.]
Solution found in 5 steps.
Basis: [11, 7, 3]

{'basis': [11, 7, 3],
 'μ_a': array([0. , 0. , 0. , 0.5, 0. , 0. , 0. , 0.5, 0. , 0. , 0. , 0.5]),
 'p_z': array([5., 6., 4.])}

More detail

Scarf's lemma assumes that $M$ and $\Phi$ are two $n \times m$ matrices in standard form, meaning that $M$ can be written by block:

$M = \begin{pmatrix} \text I_n & N \end{pmatrix} $

where $\text I_n$ is the identity matrix of size $n$, $N$ is a $(m-n)\times n$ matrix, and

for all $i \in \left\{1, \dots, n\right\}$ and $j \in \left\{1, \dots, n\right\}$ and $j' \in \left\{n+1, \dots, m\right\}$, $\Phi_{ii} \leq \Phi_{ij'} \leq \Phi_{ij}$.

Let's create some other examples, still taken from Scarf (1967):

def create_scarf_ex_sec3():
    K = 100
    payoffs = [[6, 6, K], [2, 8, K], [12, K, 2], [3, K, 8], [K, 7, 5], [K, 2, 9]]
    nbz=len(payoffs[0])
    Φ_z_a = np.hstack([K*(np.ones((nbz,nbz)) - np.eye(nbz)),np.array(payoffs).T])
    M_z_a = (Φ_z_a != K) * 1.0
    return (Φ_z_a, M_z_a)

scarf_ex_sec3 = OrdinalBasis(*create_scarf_ex_sec3())
scarf_ex_sec3.solve(verbose=1)

Solution found in 4 steps.
Basis: [7, 6, 3]

{'basis': [7, 6, 3],
 'μ_a': array([0. , 0. , 0. , 0.5, 0. , 0. , 0.5, 0.5, 0. ]),
 'p_z': array([3., 6., 5.])}

LCP formulation

Idea: define $D_a =\max_{z \in Z} \left\{ p_z - \Phi_{za} \right\}$

We have the following mixed LCP:

$0 \leq \mu_a \perp D_a \geq 0$

$0 \leq \pi_{za} \perp D_a - p_z + \Phi_{za} \geq 0$

$p_z \in \mathbf R, \left( A\mu \right)_z = q_z$

$D_a \in \mathbf R,\sum_z \pi_{za} = 1$

import gurobipy as grb

def solve_lcp(self,verbose=0):
    m = grb.Model()
    if verbose==0:
        m.Params.OutputFlag = 0
    m.Params.NonConvex = 2
    mu_a = m.addMVar(self.nba)
    d_a = m.addMVar(self.nba)
    p_z = m.addMVar(self.nbz, lb = -grb.GRB.INFINITY)
    delta_z_a = m.addMVar((self.nbz,self.nba))
    pi_z_a = m.addMVar((self.nbz,self.nba))
    m.addConstr(self.M_z_a @ mu_a == self.q_z)
    m.addConstr(d_a[None,:] - p_z[:,None] + self.Φ_z_a == delta_z_a)
    m.addConstr(pi_z_a.sum(axis=0) == 1)
    m.setObjective( (d_a * mu_a).sum() + (delta_z_a*pi_z_a).sum(), sense = grb.GRB.MINIMIZE)
    m.optimize()
    print('pi_z_a:'), print(pi_z_a.x)
    print('d_a:'), print(d_a.x)
    return {'μ_a': mu_a.x, 'p_z': p_z.x}

OrdinalBasis.solve_lcp = solve_lcp

scarf_ex_sec3.solve_lcp(verbose=0)

Set parameter Username
Academic license - for non-commercial use only - expires 2025-12-03
pi_z_a:
[[1. 0. 0. 0. 1. 0. 1. 0. 0.]
 [0. 1. 0. 1. 0. 0. 0. 0. 1.]
 [0. 0. 1. 0. 0. 1. 0. 1. 0.]]
d_a:
[3. 6. 5. 0. 1. 3. 0. 0. 4.]

{'μ_a': array([0. , 0. , 0. , 0.5, 0. , 0. , 0.5, 0.5, 0. ]),
 'p_z': array([3., 6., 5.])}

We would like to simplify this formulation.

If we assume $\Phi_{za} > 0$ for all $za$ (this is without loss of generality), then $p_z = \min_a \{ \Phi_{za} + D_a \} > 0$. Hence we can restrict $p_z$ to positive values and relax the equality $\sum_a M_{za} \mu_a = q_z$, to rewrite

$0 \leq p_z \perp q_z - \sum_a M_{za} \mu_a \geq 0$.

Furthermore, in $D_a$ free, $\sum_z \pi_{za} = 1$, we already have $D_a \geq 0$ and we can relax the equality (the $\pi_{za}$ are indicators, all that matters is whether they are positive or 0). We write instead:

$0 \leq D_a \perp \sum_z \pi_{za} - 1 \geq 0$.

To summarize, we now have a proper LCP:

$0 \leq \mu_a \perp D_a \geq 0$

$0 \leq \pi_{za} \perp \Phi_{za} + D_a - p_z \geq 0$

$0 \leq p_z \perp q_z - \sum_a M_{za} \mu_a \geq 0$

$0 \leq D_a \perp \sum_z \pi_{za} - 1 \geq 0$.

Since the value of $\sum_z \pi_{za}$ does not matter (except being 0 or positive), we can introduce without loss the following constraint into the problem:

$\mu_a = \sum_z \pi_{za} - 1$.

Then the first and last lines of the LCP become redundant, so we eliminate the first, and we replace $\mu_a$ by its expression in the third line. We get:

$0 \leq \pi_{za} \perp \Phi_{za} - p_z + D_a \geq 0$

$0 \leq p_z \perp q_z + \sum_a M_{za} - \sum_a M_{za} \sum_x \pi_{xa} \geq 0$

$0 \leq D_a \perp \sum_z \pi_{za} - 1 \geq 0$.

In matrix form this is:

$0 \leq \pi \perp \Phi - M_Z^\top p + M_A^\top D \geq 0$

$0 \leq p \perp q + M 1_A - (1_Z^\top \otimes M) \pi \geq 0$

$0 \leq D \perp M_A \pi - 1_A \geq 0$

where $M_Z$ and $M_A$ are the usual margining-out matrices, such that $(M_Z \pi)_z = \sum_a \pi_{za}$ and $(M_A \pi)_a = \sum_z \pi_{za}$,

\begin{equation} M_Z = \text I_Z \otimes 1_A^\top, \qquad\qquad M_A = 1_Z^\top \otimes \text I_A. \end{equation}

def solve_lcp_six(self,verbose=0):
    m = grb.Model()
    M_z_a = self.M_z_a.copy()
    Φ_z_a = self.Φ_z_a.copy()
    print('M_z_a:'), print(M_z_a)
    print('Φ_z_a:'), print(Φ_z_a)
    Φ_z_a += 1
    if verbose==0:
        m.Params.OutputFlag = 0
    m.Params.NonConvex = 2
    π_z_a = m.addMVar(self.nbz*self.nba)
    p_z = m.addMVar(self.nbz)
    d_a = m.addMVar(self.nba)
    w_z_a = m.addMVar(self.nbz*self.nba)
    w_z = m.addMVar(self.nbz)
    w_a = m.addMVar(self.nba)
    M_Z = np.kron(np.eye(self.nbz), np.ones((1,self.nba)))
    M_A = np.kron(np.ones((1,self.nbz)), np.eye(self.nba))
    kronprod = np.kron(np.ones((1,self.nbz)), M_z_a)
    m.addConstr(w_z_a == Φ_z_a.flatten() - M_Z.T @ p_z + M_A.T @ d_a)
    m.addConstr(w_z == self.q_z + M_z_a @ np.ones(self.nba) - kronprod @ π_z_a)
    m.addConstr(w_a == M_A @ π_z_a - 1)
    m.setObjective(π_z_a @ w_z_a + p_z @ w_z + d_a @ w_a, sense = grb.GRB.MINIMIZE)
    m.optimize()
    print('π_z_a:'), print(π_z_a.x.round(3).reshape(self.nbz, self.nba))
    return {'π_z_a': π_z_a.x.round().reshape(self.nbz, self.nba), 'p_z': p_z.x-1, 'd_a': d_a.x, 
           'μ_a': M_A @ π_z_a.x - 1}

OrdinalBasis.solve_lcp_six = solve_lcp_six

scarf_ex_sec3.solve_lcp_six(verbose=0)

M_z_a:
[[1. 0. 0. 1. 1. 1. 1. 0. 0.]
 [0. 1. 0. 1. 1. 0. 0. 1. 1.]
 [0. 0. 1. 0. 0. 1. 1. 1. 1.]]
Φ_z_a:
[[  0. 100. 100.   6.   2.  12.   3. 100. 100.]
 [100.   0. 100.   6.   8. 100. 100.   7.   2.]
 [100. 100.   0. 100. 100.   2.   8.   5.   9.]]
π_z_a:
[[1.  0.  0.  0.  1.  0.  1.5 0.  0. ]
 [0.  1.  0.  1.5 0.  0.  0.  0.  1. ]
 [0.  0.  1.  0.  0.  1.  0.  1.5 0. ]]

{'π_z_a': array([[1., 0., 0., 0., 1., 0., 2., 0., 0.],
        [0., 1., 0., 2., 0., 0., 0., 0., 1.],
        [0., 0., 1., 0., 0., 1., 0., 2., 0.]]),
 'p_z': array([3., 6., 5.]),
 'd_a': array([3., 6., 5., 0., 1., 3., 0., 0., 4.]),
 'μ_a': array([0. , 0. , 0. , 0.5, 0. , 0. , 0.5, 0.5, 0. ])}

from mec.gt import LCP

def OrdinalBasis_lemke_solve(self,verbose=0): 
    Φ_z_a = self.Φ_z_a.copy() + 1 # +1 to ensure p_z > 0
    M_Z = np.kron(np.eye(self.nbz), np.ones((1,self.nba)))
    M_A = np.kron(np.ones((1,self.nbz)), np.eye(self.nba))
    kronprod = np.kron(np.ones((1,self.nbz)), self.M_z_a)
    π_names_z_a = ['π_'+str(z+1)+'_'+str(a+1) for z in range(self.nbz) for a in range(self.nba)]
    p_names_z = ['p_'+str(z+1) for z in range(self.nbz)]
    d_names_a = ['d_'+str(a+1) for a in range(self.nba)]
    z_names_i = π_names_z_a + p_names_z + d_names_a
    w_names_i = ['comp_'+z_names_i[i] for i in range(self.nbz*self.nba+self.nbz+self.nba)]
    lcp = LCP(M_i_j = np.block([[np.zeros((self.nbz*self.nba,self.nbz*self.nba)), -M_Z.T, M_A.T],
                                [-kronprod, np.zeros((self.nbz,self.nbz+self.nba))],
                                [M_A, np.zeros((self.nba,self.nbz+self.nba))]]),
              q_i = np.hstack([Φ_z_a.flatten(), self.q_z + self.M_z_a @ np.ones(self.nba), -np.ones(self.nba)]),
              z_names_i = z_names_i, w_names_i = w_names_i)
    sol = lcp.lemke_solve(verbose)
    if sol is not None:
        zsol, _ = sol
        π_z_a = zsol[:(self.nbz*self.nba)]
        p_z = zsol[-(self.nbz+self.nba):(-self.nba)]
        d_a = zsol[-self.nba:]
        return {'π_z_a': π_z_a.reshape(self.nbz, self.nba), 'p_z': p_z-1, 'd_a': d_a, 
                'μ_a': M_A @ π_z_a - 1}
    else:
        return None

OrdinalBasis.lemke_solve = OrdinalBasis_lemke_solve

scarf_ex_sec3.lemke_solve(verbose=2)

z_0 enters, comp_d_1 departs
==========
Solution not found: Ray termination.

LCP formulation 2

We use the duality:

$p_z = \min_a \{ \Phi_{za} + D_a \}$

$D_a = \max_z \{ p_z - \Phi_{za} \}$

Consider the right-hand parts of $\Phi$ and $M$ only (without the autarky columns). Consider the LCP:

$0 \leq \mu_a \perp D_a \geq 0$

$0 \leq \pi_{za} \perp \Phi_{za} + D_a - p_z \geq 0$

$0 \leq p_z \perp q_z - \sum_a M_{za} \mu_a \geq 0$

$0 \leq D_a \perp \sum_z \pi_{za} - 1 \geq 0$.

Use the following algorithm:

Step 0. Initially, set all variables to 0 except for $D_a$ which is set to a very large value $K$.

Step 1. Choose a $a$ and increase $\mu_a$ until one of the constraints $q_z - \sum_a M_{za} \mu_a$ hits.

Step 2. Take the corresponding $z$ and increase $p_z$ until one of the constraints $\Phi_{za} + D_a - p_z \geq 0$ hits.

Step 3. Take the corresponding $za$ and increase $\pi_{za}$ until one of the constraints $\sum_z \pi_{za} \leq 1$ hits.

Step 4. Take the corresponding $D_a$ and set it to 0.

Step 5. Now either this $a$ corresponds to a $\mu_a$ that we have already increased, or it does not. If it does, we stop and we have found a solution. If not, we take the corresponding $\mu_a$ and go back to step 1.