There exists a well-known efficient method for solving bimatrix games, namely, the Lemke-Howson method (LH). The proposed Hybrid method using a 2LP algorithm first obtains a point with a small Nash criterion value. Then it transforms this point to a starting point for the LH method. The numerical results of solving bimatrix games of the dimension up to 5000 are also presented.
X = {x = (x1, . . . , xm) ∈ E+m, xTem = 1},
Y = {y = (y1, . . . , yn) ∈ En+, yTen = 1}, and the players’ payoff functions f1(x, y) = xTAy, f2(x, y) = xTBy, defined on the compact subset Z = X × Y of the Euclidean space E = Em × En. Here and onward, for any l > 0, the subset El+ is the nonnegative orthant of the space El, the vector el has all its components equal to 1. Since the payoff functions fi(x, y) (i = 1, 2) are bilinear, the game Γ is convex, hence the set of its Nash points Z∗ is nonempty (but not necessarily convex). The definition of convex games and their properties can be found in [Osborne, 2004]. (1) (2)
For the game Γ, one can introduce the standard (unnormed) Nash function as follows
F (z) = max f1(x, y) − f1(z) + max f2(x, y) − f2(z), x∈X y∈Y z ∈ Z.
Firstly such a function was introduced by H. Mills [Mills, 1960]. The inequality F (z) > 0 is true for all z = (x, y) ∈ Z. Moreover, F (z) = 0 if, and only if z is a Nash equilibrium or point of the game Γ. Therefore, the set of global minimum points of the function F = F (z) on Z coincides with the set Z∗ of the Nash points of the game Γ, with the global minimum value 0. In other words, solving the game Γ is tantamount to finding the global minimum of the function F = F (z) over the subset Z, and the value of F (z) can be treated as a natural measure of deviation of the point z from the set Z∗. 3
The idea of consecutive solution of the LPs (see (4),(5) below) followed from the initial bilinear problem was firstly introduced by H. Konno [Konno, 1976]. For bimatrix games this idea was used by B. M. Mukhamediev [Mukhamediev, 1978] , A. V. Orlov and A. S. Strekalovsky [Orlov & Strekalovsky, 2005] and E. G. Golshtein et al. [Golshtein et al., 2013]. Even though the problem of minimization of the function F (z) with respect to z ∈ Z is extremely difficult because this function has a lot of local minima distinct from the global one, it is easy to reduce its (partial) minimization with respect to, e.g., the variable x (as part of the global variable z = (x, y)) while the other part (the variable y) is fixed, to a linear program. (The same is true if one decides to minimize F with respect to y with the first strategy x fixed.)
Indeed, let us fix the vector y (strategy of player 2) with values y = y′ = (y1′, . . . , yn′) ∈ Y and consider the following linear programming problem P1(y′) with m + n + 1 variables (x, v, α), where v = (v1, . . . , vn) is a vector of auxiliary variables, and with n + 1 equality constraints: −xT(A + B)y′ + α → min, xTB + v = α, xTem = 1, x ∈ E+m, v ∈ En+, α ∈ E1.
Denote as x∗ = x∗(y′) the first m components of an optimal solution of problem (4), then F (x∗, y′) 6 F (x′, y′) for every x′ ∈ X.
Similarly, having fixed the first strategy x with some values x = x′ ∈ X, one can introduce the linear program P2(x′) with n + m + 1 variables (y, u, β) (where u = (u1, . . . , um) is a vector of auxiliary variables) and m + 1 equality constraints: −(x′)T(A + B)y + β → min,
Ay + v = β, yTen = 1, y ∈ En+, u ∈ E+m, β ∈ E1.
Again, if y∗ = y∗(x′) is the vector of the first n components of an optimal plan of problem P2(x′), then clearly F (x′, y∗) 6 F (x′, y′) for an arbitrary vector y′ ∈ Y .
Now starting with an arbitrary initial strategy y0 ∈ Y of player 2, let us generate the sequence of iterations {zt} = {(xt, yt)}, t > 1, as follows: xt+1 is a solution of P1(yt), while yt+1 solves problem P2(xt+1), t = 0, 1, 2, . . .
It is evident the inequality F (zt+1) 6 F (zt) holds for all t > 1.
Let T be the maximum permitted number of iterations (pairs of tasks (4)-(5)), εf > 0 and εF > 0 be small numbers, S be the set of pure starting strategies. For any z ∈ S define t(z) = min{t ∈ {1, . . . , T }: F (zt) − F (zt−1) 6 εf }; We put t(z) = T at the failure of the inequality for all t. Denoting Q(z) = F (zt(z)), we have Q(S) = minz∈S Q(z) the lowest value that can be obtained. The game is approximately solved at the point z with accuracy εF if F (z) 6 εF . If Q(S) > εF we believe that the solution to the game not found. 4
The LH-algorithm replaces the solution of a bimatrix game by the search of a solution of the system of linear equations in a linear complementarity problem related to the game. By making the following change of variables x := x/α, y := y/β in the system (4)–(5), we obtain the system of linear equalities and inequalities (call it the LH-problem):
Ay + u = em, xTB + v = en, x, u ∈ E+m, y, v ∈ En+,
(6) (3) (4) (5) lacking the equalities xTem = 1, yTen = 1. Introduce the notation
H = (Om
BT
A ) On , z = (z1, . . . , zm+n) ≡ (x, y),
e ≡ em+n, w = (w1, . . . , wm+n) ≡ (u, v), o = (0, 0, . . . , 0) ∈ Em+n,
E = diag (e), , where Om and On are the square null matrices of dimension m and n, respectively; E is the unit matrix of dimension m + n, and thus yield the linear complementarity problem: Find a nonnegative (z, w ∈ E+m+n) solution of the linear system
Hz + Ew = e, (7) such that ziwi = 0 for each 1 6 i 6 m + n (complementarity conditions). Denote by ∆(z, w) the number of complementarity conditions broken at the point (z, w) ∈ E2+(m+n). We will say that a solution (z, w) of system (7) satisfies almost all complementarity conditions if ziwi ̸= 0 for only one index i, 1 6 i 6 m + n.
The initial basis of the LH-problem is composed of all slack variables w. The latter satisfies the complementarity conditions, however, z = o, which means that the equalities xTem = 1, yTen = 1 don’t hold. By having introduced into the basis one of the structural variables (say, z1) and violated therewith one of the complementarity conditions, the LH-method starts from the initial point (z1, w1) and makes further steps (iterations) of the simplex-method-type until the broken complementarity condition becomes true. The generated sequence of vertices (satisfying almost all complementarity conditions) forms the so-called Lemke path running the feasible bases of system (7). After several iterations of the LH-method, the structural variables x and y will boast various non-zero entries. One can make the inverse transformation of the variables back to their original form and continue the Lemke steps by pivoting the bases of the initial problem taking into account the constraints xTem = 1, yTen = 1. In this manner, we exclude the possibility of returning to the initial point z = o.
The algorithm stops if either it has run into a previously generated point (that is, a loop has been closed) or the current point satisfies all complementarity conditions, i.e., it solves the linear system. As a rule (but not always) the algorithm finds an exact solution of system (7), hence, it solves the game Γ. In order to come back to the original variables, one must normalize the just obtained solution of the system as follows: x := x/(xTem), y; = y/(yTen).
The LH-method (as well as its derivative algorithms LH1 and LH2) as steps of Hybrid algorithm is founded on the LH-procedure that generates a sequence of pivoted bases under the principal rule: the leaving variable is selected so that the other basic variables keep having nonnegative values, and the entering variable is either zk or wk (1 6 k 6 n + m), which is complementary to the leaving variable (that is, wk or zk, respectively).
Thus, the LH-procedure determines the Lemke path uniquely after the starting structural variable has been selected to introduce into the initial basis. The pitfalls of the LH-method are, on the one hand, a possibility of loops (returns to the initial point), and on the other hand, a forbiddingly large number of iterations. If no solution has been found for the selection of z = z1 as the initial point, we start the method again with the initial z = z2, and so on. The game is considered as unsolved if no solution is found after all the initial points zi, i = 1, . . . , n + m, have been tried. 5
We say that a point z = (x, y) ∈ Em+n approximately solves the bimatrix game Γ, whenever F (z) 6 εF and ∆(z, w) 6 1 (recall that ∆(z, w) = 0 at the exact solution of Γ). Next, we say that a point z = (x, y) is promising if ε<F (z) 6 F and 1 6 ∆(z, w) 6 ∆, where εF , F , and ∆ are some preliminary selected parameters.
Let z be a promising starting point for the game Γ. With the former, we start the LH-method’s version that finds the solution (Nash point) with a lower number of broken complementarity conditions. This version is named as the LH1-algorithm.
In contrast to the “classical” LH-method where the unit matrix determines the initial basis and the complementarity condition ∆(z, w) = 0 holds for z = o, w = e, in the LH1-algorithm, ∆(z, w) = k > 0, i. e., the initial basis comprises (in an arbitrary order) 2k basic variables zi, wi, i = i1, . . . , ik (called “offenders”), for which ziwi > 0. At the same time, there are 2k nonbasic variables (“candidates” for entering the basis) such that zj = wj = 0, j = j1, . . . , jk. In order to diminish the number of broken complementarity conditions, some of the candidate variables are to enter the basis to replace the offenders.
An initial basis of an LH-problem for a “hot start” is combined from the bases (without variables α and β) for problems (4) and (5), obtained by the 2LP-method, what we use as first step of Hybrid algorithm.
We try to introduce all the candidate variables one by one instead of the offenders. If we have managed to do that with success and all the complementarity conditions have been satisfied, then the LH-problem (6) has been solved. Otherwise, that is, when entering a candidate variable into the basis, the leaving basic variable doesn’t break the corresponding complementarity condition (i.e., it isn’t an offender), then we again use the LH-procedure and introduce the variable complementary to the variable that has just left the basis. We repeat these operations until either the offender has been expelled from the basis and ∆0 has dropped, or the number of those attempts has reached the upper limit of I(m + n), where the constant I has been preliminary selected. It is also possible that a loop closes before the upper limit has been achieved: in this case, we simply switch to the next candidate variable without changing the value of ∆.
In our test numerical experiments, in the majority of the considered examples, the proposed algorithm did manage to exclude from the basis one variable from each pair of offenders, that is, all the broken complementarity conditions were successfully repaired. If there was left broken complementary conditions we go to step 6 as third step of Hybrid algorithm. 6
If the LH1-algorithm has failed with excluding all offenders from the basis, we will try to do that with the aid of a so-called supplementary column.
First, we will forcefully replace one of the variables in each offending pair by a candidate variable from the list of candidate pairs with the aim of repairing the k broken complementarity conditions that the LH1-method has failed to remove. Let kx > 0 and ky > 0 be the numbers of the invalid complementarity conditions related to the variables x and y, respectively.
If kx > 0 then among all possible kx2 candidates for entering the basis we select the variable that maximizes the minimum negative value of the current basic variables among such negative basic values that must appear after the forceful change of the basis. Now we update kx by subtracting 1 and repeat the procedure until we reach kx = 0. The similar actions are undertaken towards ky > 0.
Having selected the next candidate variable, decompose the corresponding column with respect to the current basis. Among the coefficients of this decomposition, find those belonging to the positions related to the offenders. Detect the maximum (or large enough) absolute value among the just selected coefficients and introduce the candidate variable exactly in this position.
A supplementary column should be generated in such a way that as a result of having introduced it into the basis, all the negative values of the basic variables disappear. In order to do that, we compose the supplementary column as a sum of the columns that correspond to the variables having negative values in the current basic solution. This new column multiplied by the minimum negative weight (that is, having the maximum absolute value among the above mentioned negative values) is introduced into the basis and occupies the position corresponding to the same basic variable with the most negative value. It is easy to check that after that, all basic variables will acquire positive values, while the variable related to the supplementary column will get the value 1.
Next, we begin the LH-procedure starting from the variable complementary to the variable that has left the basis when introducing into it the supplementary column. We repeat the steps of the LH-procedure until either the supplementary column leaves the basis (that is, a Nash point has been obtained), or the number of iterations exceeds the upper limit. 7
Now that all the necessary parts (the 2LP-method, LH1- , and LH2-algorithms) have been described we are in a position to present our Hyb-Algorithm for solving bimatrix games as a solid alternative to the LH-method and even as a possible tool for finding a solution to the game when the LH-method fails.
1. 1.
⇑ ⇑ 0. ⇒ 1. ⇒ 2. ⇒ 3. ⇒ 4. ⇒ stop ⇓ ⇓ stop stop Step 0. Define the start point, for example, the pure strategy y1 = 1, yi = 0, i = 2, n.
Step 1. Going through the set of initial points (pure strategies) look for a prospective starting point (F (x, y) 6 εF ) by an 2LP algorithm, otherwise use a new start (pure strategy) point. If the set of starting points has exhausted, the game Γ is reported as unsolved.
Scheme of hybrid algorithm
Step 2. Collect the basis of LH problem by making use of the solution obtained by the 2LP-algorithm. Detect the broken complementarity conditions and list the offenders (the basic variables in broken complementary conditions) and candidates (the pairs of nonbasic variables related by the complementary conditions).
Step 3. (Key) Attempt to remove the broken complementary conditions by introducing one candidate from each pair as starting points into the basis in LH (LH1) method with \hot start". If in the end, all the offenders have been excluded from the basis, a Nash point has been obtained and the game solved. Otherwise, go to Step 4 to generate a supplementary column.
Step 4. It will get a equilibrium point by using the supplementary column if step 3 fails. In all pairs offenders will forced replace one of them by one of the candidates. In the obtained solution will be negative values. An supplementary column is collected as the sum of columns that correspond to variables with negative values in the basic solution. This column is used as starting point for LH (LH2) method with \hot start".
If the solution was not obtained through steps 1, 3 or 5, repeat from 1 step. 8
The tested algorithms (2LP, LH, and Hyb) were coded and implemented with the software MATLAB 7.14.0.739 R2012a. The linear subprograms were solved by the procedure of cplexlp from the IBM ILOG CPLEX 12.6.2, which is much more efficient than the procedure linprog of MATLAB. The personal computer IntelCore i52400 CPU (3.1 GHz), 4 Gb RAM was used for the computation in the majority of cases.
The algorithms were tested and their efficiency was compared mainly on the base of a family of bimatrix games containing 15 series of medium-size games with the parameters m, n ∈ {20, 40, 60, 80, 100}, m 6 n. For each series, R = 100 of problems were initially generated. The main parameters of the algorithm(s) (if not specified otherwise) are as follows: S = m, εF = 10−4, F = 0.01, I = ∆ = 5, I = 1.
With the purpose of saving the efforts and running time, the ordered set X0 = {e1m, . . . , emm} of pure strategies of player 1 was selected as the set of possible starting points. Here, eim is the standard unit vector of the space Em having all entries zero except entry i, which is 1.
The matrices A = (aij ) and B = (bij ) were generated in two stages. First, the entries a′ij and b′ij were generated independently by the MATLAB generator of pseudo-random numbers uniformly distributed in the interval (0, 1). After that, we applied the transformation aij = a′ij + φb′ij , bij = b′ij + φa′ij , 1 6 i 6 m, 1 6 j 6 n, where φ is the coefficient of the matrix mutual interdependency, 0 6 φ 6 1. The latter transformation implies that the entries of both matrices lie within the range of (0, 1 + φ).
The numerical experiments investigated the performance and compared the efficiency of the algorithms 2LP, LH, and Hyb when solving bimatrix games. We also examined whether it was possible, starting with the 2LP method, to obtain such a starting point for the LH-algorithm that would help find an exact solution of the game with a save on the computational efforts (iterations, running time, etc.).
The numerical and testing results are given in the three tables below, where, in Tables 1–2, for the reader’s convenience, the results are shown for 5 series of games, while the lines Sum provides the numbers obtained by summing up over all the series of games.
The games having solutions in pure strategies needn’t be solved by the LH and Hybrid methods; Since analyzing preliminary payoff matrixes we can simply establish whether the game has Nash equilibrium in pure strategies.
Table 1 accumulates the testing results obtained for the family of 1500 games treated by the LH-method for three values of the matrix mutual interdependency parameter (φ ∈ {0; 0.1; 0.2}). The notation of Table 1 means: m, n are the game’s dimension parameters; the parameter I determines the upper limit (equal to I (m + n)) of iterations allowed in the LH-algorithm;kPS is the number of games having solutions in pure strategies kLH is the number of games solved by the LH-algorithm; kNO is the number of games on which the LH-algorithm failed. To resume, kPS + kLH + kNO = 100; In addition, for each unsolved game, the total number of tested starting points was m and the total number of iterations for every starting point didn’t exceed I(m + n); the value of Ssum equals the total number of starting points tested, JLH is the total number of simplex-type iterations, Time is the running time (in seconds) that the LH-method spent to find the solution, Smax is the maximum number of the starting points tested when solving one (solved) game.
Analysis of the data from Table 1 allows one to make the following conclusions.
1. Generation of the independent matrices A and B gets 935 games solvable in the pure strategies. Therefore, when φ = 0 one needs to solve only 565 games. For φ = 0.1 we had to solve 1477 games, for φ = 0.2 we treated φ 0 0.1 0.2 0,5 5 5
Sum
935 1496 games. For φ 6 0.1, ALL the games from the family were solved by the LH-algorithm, while for φ = 0.2, the LH-method failed on 25 games of the total of 1500.
2. The value of the parameter I determining the maximum allowed number of iterations when solving a game seriously affects the efficiency of the LH-method. When selecting an appropriate value for this parameter I, one should take into account the value of φ as well as the dimension of the solved problem. For our test examples, the following values proved to be well-suited (but not necessarily optimal): I = 0.5 for φ = 0 and I = 5 for φ either 0.1 or 0.2.
Our test computations have proven that it makes sense to classify the initial points zi, i = 1, . . . , m + n, by trying to produce no more than I(n + m) iterations with each one, where 0.01 < I 6 5. As one can see from the numerical results given below, the efficiency of the algorithm vitally depends on the constant value of I, which is desirable to select particularly for each game. kPS kLH kNO
kPS k2LP kLH1 kLH2 kNO kLH2 = 65; for φ = 0.1 there happened k2LP = 429, kLH1 = 760, kLH2 = 288; for φ = 0.2 there appeared k2LP = 387, kLH1 = 717, kLH2 = 332.
3. The average number Se of the used starting points for one solved problem for φ = 0; 0.1; 0.2 was, respectively, 2.3; 5; 14. As φ grows the ratio SHyb/Ssum of the promising starting points to the total number of the (used) starting points drops as follows: (0.9; 0.8; 0.5). The total running time needed to solve all tested games prove to be 20, 138, and 469, respectively.
Table 3 provides the numerical results obtained by the LH- and Hyb-algorithms when solving several bimatrix games of a large size. Here, the following notation is used: the parameter I determines the upper bound (equal to I(m + n)) for the number of iterations allowed to produce by the LH-algorithm; m, n define the size of the bimatrix game; ϱ is the number of the games in the serie, while kNO denotes the number of unsolved games in the series in question. Ssum gives the total number of the used starting points, Jsum is the total number of the simplex-type iterations conducted, and Time denotes the aggregate running time (in seconds).
Having analyzed the data of Table 3, we come to the following conclusions.
1. Once again, the results have shown that a lucky choice of the parameter I greatly affects the velocity of the LH-algorithm. For instance, when solving the 400 × 800 game, for the value of I = 500 (that is, the upper bound of iterations equalling I(m + n) = 600000) the totality of iterations (substitutions of the basis) resulted in 4 × 600000 + 221 = 2400221, whereas for I = 0, 25 only 4 × 0, 25(400 + 800) + 221 (iterations on the 5-th start point) u = 1421 were made.
In both cases, the final result was achieved after the fifth starting point has been employed and altogether 221 iterations has been made, but the difference in the running time was enormous. However, the Hyb-algorithm made use of only three starting points, and the final solution was obtained by the LH1-procedure, after 1464 iterations performed for 17 seconds.
2. When applied to the two games listed in Table 3 with the dimensions 400 × 800 and 500 × 500, both with the value of φ = 0.1, the LH-algorithm failed for various values of the parameter I. At the same time, these games were solved successfully with the Hyb-method: the game of the size 400 × 800 was solved after having checked 208 starting points and conducted 381653 iterations for1382 seconds; the solution of other game of the Num kNO dimensions 500 × 500 was similarly finished after having tried 39 starting points and applied 74978 iterations for 179 seconds.
3. Both algorithms efficiently solve the games with independent matrices at about the same computational cost provided that the value of the parameter I has been appropriately selected. For instance, when solving the problem of the size 3000 × 3000 the Hyb-algorithm made use of 328 starting points and performed 485359 iterations for 17238 seconds. The LH-method had to check 6004 starting points and perform 358826 iterations for 19374 seconds, with the value of the parameter I = 0.01 (the upper bound for the number of iterations thus being I (m + n) = 60).
You can see on Figure 1. the results of solving two problems (with interdependent matrices) by LH method and the hybrid algorithm. LH method failed to solve these problems using all the initial points (400 and 500, respectively). The hybrid algorithm got the solution of the first problem on the 208-th start point using the additional column, and the second problem was solved on the 39th starting point on the local search step. All the above presented leads to the following conclusion: The proposed Hyb-algorithm has proven to be competitive as compared to the 2LP-method and LH-algorithm. Moreover, in some instances, the Hyb-method manages to solve the games on which both the 2LP-method and LH-algorithm fail.