NumPivots blows up (only for Poly2) in EnumMixed::Solve

I am using a stripped version of gambit that uses only the files necessary to run the EnumMixed solver. Precisely:

It removes the Game representations and populates the matrices A1, A2 from my own Matrix struct with entry type double
It removes the conversions from double (in the input file) to Number to Rational back to double. The arithmetic of the Tableau is all of type double anyway (at least when T is the same). I can't imagine that it makes any difference in regards to vertex enumeration algorithm.
It stops enumerating the vertices/extreme equilibria after the first one is found to be nash. I also don't bother with calculating the Cliques since for my purposes I just need a single NE.

I don't think the above considerations matter at all in regards to the bug. I've posted the modified Solve function below. All the other includes are untouched

   const int rows = row_payoff_matrix.rows;
    const int cols = row_payoff_matrix.cols;
    /*
    Accept surksit style matrices and use gambit protocol to solve
    */

    // Construct matrices A1, A2 // T = double
    Gambit::Matrix<double> A1(1, rows,
                              1, cols);
    Gambit::Matrix<double> A2(1, cols,
                              1, rows);

    double max = std::max(row_payoff_matrix.max(), col_payoff_matrix.max());
    double min = std::min(row_payoff_matrix.min(), col_payoff_matrix.min());

    if (min > 0)
    {
        min = 0;
    }
    min -= 1;
    if (max < 0)
    {
        max = 0;
    }

    const double fac = 1 / (max - min);

    for (size_t i = 1; i <= rows; i++)
    {
        for (size_t j = 1; j <= cols; j++)
        {
            A1(i, j) = fac * (row_payoff_matrix.get(i - 1, j - 1) - min);
            A2(j, i) = fac * (col_payoff_matrix.get(i - 1, j - 1) - min);
        }
    }

    // Construct vectors b1, b2 // Length rows, cols with default value -1
    Gambit::Vector<double> b1(1, rows);
    Gambit::Vector<double> b2(1, cols);
    b1 = (double)-1;
    b2 = (double)-1;

    // enumerate vertices of A1 x + b1 <= 0 and A2 x + b2 <= 0
    Gambit::linalg::VertexEnumerator<double> poly1(A1, b1); // Instantiantion calls method 'Enum()'
    Gambit::linalg::VertexEnumerator<double> poly2(A2, b2);

    const Gambit::List<BFS> &verts1(poly1.VertexList()); // List is computed in prior instantiation. This simply returns ref to that LIst member
    const Gambit::List<BFS> &verts2(poly2.VertexList());
    int m_v1 = verts1.Length();
    int m_v2 = verts2.Length();

    Gambit::Array<int> vert1id(m_v1);
    Gambit::Array<int> vert2id(m_v2);
    for (int i = 1; i <= vert1id.Length(); vert1id[i++] = 0)
        ;
    for (int i = 1; i <= vert2id.Length(); vert2id[i++] = 0)
        ; // Initialize vert1id, 2id lists with resp. lengths as 0

    int i = 0;
    int id1 = 0, id2 = 0;

    for (int i2 = 2; i2 <= m_v2; i2++)
    {
        const BFS &bfs1 = verts2[i2];
        i++;
        for (int i1 = 2; i1 <= m_v1; i1++)
        {
            const BFS &bfs2 = verts1[i1];

            // check if solution is nash
            // need only check complementarity, since it is feasible
            bool nash = true;
            for (size_t k = 1; nash && k <= rows; k++)
            {
                if (bfs1.count(k) && bfs2.count(-k))
                {
                    nash = EqZero(bfs1[k] * bfs2[-k]); // each value for
                }
            }
            for (size_t k = 1; nash && k <= cols; k++)
            {
                if (bfs2.count(k) && bfs1.count(-k))
                {
                    nash = EqZero(bfs2[k] * bfs1[-k]);
                }
            }

            if (nash)
            {
                double row_sum = 0;
                double col_sum = 0;
                for (size_t k = 1; k <= rows; k++)
                {
                    if (bfs1.count(k))
                    {
                        const double x = -bfs1[k];
                        row_strategy[k - 1] = x;
                        row_sum += x;
                    }
                    else
                    {
                        row_strategy[k - 1] = 0;
                    }
                }
                for (size_t k = 1; k <= cols; k++)
                {
                    if (bfs2.count(k))
                    {
                        const double x = -bfs2[k];
                        col_strategy[k - 1] = x;
                        col_sum += x;
                    }
                    else
                    {
                        col_strategy[k - 1] = 0;
                    }
                }
                for (size_t k = 1; k <= rows; k++)
                {
                    row_strategy[k - 1] /= row_sum;
                }
                for (size_t k = 1; k <= cols; k++)
                {
                    col_strategy[k - 1] /= col_sum;
                }
                return; // TODO continue calculating?
            }
        }
    }

The problem is that for larger matrices (9x9 in my case, I did not reproduce this error with 4x4), the instantiation of poly2 is MUCH slower than poly1. Recall that the constructor for these objects runs the Enum function, which does pretty much all of the Linear Programmin legwork.

Not only that, but in my use case the instantiation of poly2 will get slower and slower. I've done some debugging step-throughs and the problem seems to be that the number of recursive VertexEnumerator::Search() calls will explode (this can also be seen in the exploding npivots). Again this only affects poly2 and persists even if I swap the order of their initialization.

Additionally, this problem only exists in the context of using the solve function on Bimatrices that occur within my tree search. I have a test that will call the solve function on random bimatrices of varying sizes (n x m for 2 <= n, m <= 9) and I do not get this bug in the tests. There is nothing wrong with the bimatrices being solved in the Tree search context, and I've even modified the test so that it's bimatrices are more like the tree searches. I've noted that even with the slow down, the matrices are still being solved correctly (very low exploitability of the resulting strategies)

I understand that this is not a Gambit bug per se, but I was wondering if you might have any ideas as to what might be causing this. I don't understand Linear Programming enough to guess why the number of pivots is so large and constantly increasing, only for the second polytope!

The stripped down Gambit is this branch of my fork https://github.com/baskuit/gambit/tree/rewrite

It removes the conversions from double (in the input file) to Number to Rational back to double. The arithmetic of the Tableau is all of type double anyway (at least when T is the same). I can't imagine that it makes any difference in regards to vertex enumeration algorithm.

It is necessary to do the vertex enumeration in exact arithmetic; it is literally make or break. We are not solving a Linear Program here but a Linear Complementarity Problem, and without exact arithmetic one could readily get nonsense answers (i.e., strategies that are not Nash or even close to Nash reported as Nash) or no answers, i.e, the enumeration doesn't find any "Nash", real or fake.

For an alternative implementation of full vertex enumeration (withe the maximal bipartite clique enumeration for components too), you could take a look at: https://github.com/rahulsavani/bimatrix_solver

I am not totally sure why you are using a stripped down version, but the comments on the fork suggest it is for speed. In that case, and given that you only want a sample Nash, your best bet may be to use lrsnash directly:

http://cgm.cs.mcgill.ca/~avis/C/lrslib/USERGUIDE72.html#nash

(c.f. the -L switch of gambit-enummixed).

In any case, don't waste any time trying to find actual Nash equilibria with floating point arithmetic, especially if you are ultimately interesting in solving reasonably large games. (If you are interested in approximate equilibria then that's a totally different question, but that's not for a Gambit github issue.)

@rahulsavani

It is necessary to do the vertex enumeration in exact arithmetic; it is literally make or break.

I don't believe this is the case; It seems Gambit genuinely has a solver for double type. It uses fuzzy equality. This is also my guess as to why Gambit rarely fails to solve a matrix.

I am not totally sure why you are using a stripped down version, but the comments on the fork suggest it is for speed. In that case, and given that you only want a sample Nash, your best bet may be to use lrsnash directly:

Yes, I am implementing MatrixUCB for my tree search repo. For this algorithm, we only need an approximate NE. There may be value in a full vertex enumeration algorithm, since the Bimatrices here are not constant-sum, so it could be beneficial to search to chose a NE with maximal total payoff (my suspicion, anyway). Calculation of Cliques might also be useful for choosing a more 'regular' NE by convex combination of extreme equilibria.

I tried using lrsnash with low weight rational approximations, but even the 128 bit would frequently overflow. I have not gotten around to trying the GMP version, but I suspect its slower than EnumMixedSolver<double>.

Any suggestions for a fast approximate (Exploitability < .05 for 9x9 matrices with payoffs in [0, 1]) NE solver are very welcome! The other gambit tools didn't work in my testing and empirical strategies of joint adversarial bandits are slow and not very convergent.

I've exported the matrix in question in .nfg format below. In my fork of Gambit, this requires 177, 10227 pivots in poly1, poly2 respectively. In the .nfg form with Gambit's vanilla input functionality, it requires 12012, 5257 pivots.

NFG 1 R "" {"" ""} {9 9} {{ "" 4.0743378875130709 3.0743378875130705 } { "" 3.8243378875130705 3.3243378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 3.0743378875130705 4.0743378875130709 } { "" 3.5743378875130705 3.5743378875130705 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 3.5743378875130705 3.5743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.8243378875130705 3.3243378875130705 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 4.0743378875130709 3.0743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 4.0743378875130709 3.0743378875130705 } { "" 4.3477703358384350 4.3477703358384350 } { "" 3.5743378875130705 3.5743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 4.0743378875130709 3.0743378875130705 } { "" 3.0743378875130705 4.0743378875130709 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 4.3477703358384350 4.3477703358384350 } { "" 5.3477703358384350 4.3477703358384350 } }

I don't believe this is the case;

It is simply a fact that if you want exact equilibria you need exact arithmetic or things can go wrong. I didn't say things would always go wrong, they won't, but they can sometimes or even frequently go wrong. I know this both in theory and from concrete experiments.

It seems Gambit genuinely has a solver for double type. It uses fuzzy equality. This is also my guess as to why Gambit rarely fails to solve a matrix.

The vertex enumeration is done in exact arithmetic. I think you may be talking about EqZero being a fuzzy inequality. Sure, it is, but this is only the complementarity check (done on vertex pairs, after two vertex enumerations are done), not the vertex enumeration itself. Anyway, whatever, my point is just that doing things in exact arithmetic (in lrsnash everything is done in exact arithmetic, in the native gambit-enummixed not everything) is for good reason.

BTW, for large game, lrsnash should be much quicker than the native gambit-enummixed, since the latter enumerate both "full size" best response polytopes and then does a check on all |vertices in poly 1| * |vertices in poly 2| possible pairs. lrsnash enumerates just one "full size" polytope, solving the implied low dimensions problem for finding a corresponding vertex on the other side to complete a Nash if possible; details are in our paper: D. Avis, G. Rosenberg, R. Savani, and B. von Stengel (2010), Enumeration of Nash equilibria for twoplayer games. Economic Theory 42, 9-37. Note: lrsnash only finds extreme equilibria, the clique stuff is separate.

Since you only want a sample equilibrium, you could try the Lemke-Howson (LH) algorithm, but that will only work in exact arithmetic, but should in general be very fast (the bad examples are really pathological), and it will give you an exact equilibrium. If you have the time to try it, it's definitely worth benchmarking against. You said the other gambit tools didn't work, but LH will always find one exact Nash (possibly more if you vary the starting label, but you only need one) -- perhaps you mean it is too slow. But, I would not expect it to be slower than "lrsnash set to terminate when it finds the first equilibrium" (let's call it lrsnash-sample). I would make sure that LH and lrsnah-sample (both with exact arithmetic) are really definitely too slow for what you need before moving on. I am not sure, but perhaps your focus on doing floating point for speed was misplaced, but maybe not, maybe exact is simply too slow for your application, which is fine, but then look to other algorithms, not to make algorithms that rely crucially on exact arithmetic (e.g. lrsnash and LH) fast by a quick swap to floating. It is an interesting question how to do floating point implementations of LH / lrsnash, but this will need some period checks and backtracking, and is a research question in itself.

For algorithms designed specifically for approximation, you can try the Tsaknakis-Spirakis (TS) gradient-descent based algorithm , which is only guaranteed to find a 1/3-NE (additive eps, all payoffs in [0,1]) in polynomial time but normally does much much better. For a link to implementations and a reasonably comprehensive empirical study, which includes both the TS and LH algorithms see:

John Fearnley, Tobenna Peter Igwe, Rahul Savani: An Empirical Study of Finding Approximate Equilibria in Bimatrix Games. SEA 2015: 339-351.

This doesn't really belong here, so if you want to continue to talk please email me. Cheers.

gambitproject / gambit

NumPivots blows up (only for Poly2) in EnumMixed::Solve #334