Lemke on trees stops after finding 1 eq; stop_after/max_depth not used (#722)

Lemke's algorjthm, used for `lcp_solve` using the sequence form, finds only one equilibrium.  Previously this method was advertised as having the possibility of recursing to find more.  At least as implemented, this is not mathematically well-founded, and actually the implementation would not find another equilibrium after the first (unless possibly by pure luck).

This removes the `stop_after` and `max_depth` options from `lcp_solve` on sequence forms to clarify the operation of the method.

---------

Co-authored-by: Ted Turocy <ted.turocy@gmail.com>

Author: rahulsavani

ChangeLog

@@ -5,10 +5,11 @@
 ### Fixed
 - Sequence-form based equilibrium-finding methods returned incorrect output on games with
   outcomes at non-terminal nodes.  (#654)
+
 ### Added
 - Implement `IsAbsentMinded()` on information sets (C++) and `Infoset.is_absent_minded` (Python)
   to detect if an information is absent-minded.
-- Tests for EFG Nash solvers -- `enumpoly_solve`, `lp_solve`, `lcp_solve` -- in behavior stratgegies
+- Tests for EFG Nash solvers -- `enumpoly_solve`, `lp_solve`, `lcp_solve` -- in behavior strategies
 - In `pygambit`, `Node` objects now have a read-only property `own_prior_action` and `Infoset` objects
   have a read-only property `own_prior_actions` to retrieve the last action or the set of last actions
   taken by the player before reaching the node or information set, respectively. (#582)
@@ -27,6 +28,8 @@
   distinguish "agent" regret and liap values from their strategy-based analogs.  Methods which compute
   using the agent-form - specifically `enumpure_solve` and `liap_solve`, now clarify this by being named
   differently in `pygambit`. (#617)
+- For clarity, the `stop_after` and `max_depth` arguments to `lcp_solve` are no longer permitted when solving using
+  the sequence form.  These actually had no effect in previous versions. (#671)
 - In the graphical interface, removed option to configure information set link drawing; information sets
   are always drawn and indicators are always drawn if an information set spans multiple levels.
 - In `pygambit`, indexing the children of a node by a string inteprets the string as an action label,

doc/tools.lcp.rst

@@ -6,13 +6,20 @@
 
 :program:`gambit-lcp` reads a two-player game on standard input and
 computes Nash equilibria by finding solutions to a linear
-complementarity problem. For extensive games, the program uses the
+complementarity problem.
+
+For extensive games, the program uses the
 sequence form representation of the extensive game, as defined by
 Koller, Megiddo, and von Stengel [KolMegSte94]_, and applies the
-algorithm developed by Lemke. For strategic games, the program uses
-the method of Lemke and Howson [LemHow64]_.  There exist strategic
-games for which some equilibria cannot be located by this method; see
-Shapley [Sha74]_.
+algorithm developed by Lemke. In that case, the method will find
+one Nash equilibrium.
+
+For strategic games, the program uses the method of Lemke and Howson
+[LemHow64]_. In this case, the method will find all "accessible"
+equilibria, i.e., those that can be found as concatenations of Lemke-Howson
+paths that start at the artificial equilibrium.
+There exist strategic-form games for which some equilibria cannot be found
+by this method, i.e., some equilibria are inaccessible; see Shapley [Sha74]_.
 
 In a two-player strategic game, the set of Nash equilibria can be expressed
 as the union of convex sets. This program will find extreme points
@@ -53,9 +60,11 @@ game.
 
 .. cmdoption:: -e EQA
 
-   By default, the program will find all equilibria accessible from
-   the origin of the polytopes.  This switch instructs the program
-   to terminate when EQA equilibria have been found.
+   By default, when working with the reduced strategic game, the program
+   will find all equilibria accessible from the origin of the polytopes.
+   This switch instructs the program to terminate when EQA equilibria
+   have been found. This has no effect when using the extensive representation
+   of a game, in which case the method always only returns one equilibrium.
 
 .. cmdoption:: -h
 

src/pygambit/nash.pxi

@@ -74,15 +74,15 @@ def _enummixed_strategy_solve_rational(game: Game) -> list[MixedStrategyProfileR
 
 
 def _lcp_behavior_solve_double(
-        game: Game, stop_after: int, max_depth: int
+        game: Game
 ) -> list[MixedBehaviorProfileDouble]:
-    return _convert_mbpd(LcpBehaviorSolve[double](game.game, stop_after, max_depth))
+    return _convert_mbpd(LcpBehaviorSolve[double](game.game))
 
 
 def _lcp_behavior_solve_rational(
-        game: Game, stop_after: int, max_depth: int
+        game: Game
 ) -> list[MixedBehaviorProfileRational]:
-    return _convert_mbpr(LcpBehaviorSolve[c_Rational](game.game, stop_after, max_depth))
+    return _convert_mbpr(LcpBehaviorSolve[c_Rational](game.game))
 
 
 def _lcp_strategy_solve_double(

src/pygambit/nash.py

@@ -210,12 +210,13 @@ def lcp_solve(
         representation even if the game's native representation is extensive.
 
     stop_after : int, optional
-        Maximum number of equilibria to compute.  If not specified, computes all
-        accessible equilibria.
+        Maximum number of equilibria to compute when using the strategic representation.
+        If not specified, computes all accessible equilibria.
 
     max_depth : int, optional
-        Maximum depth of recursion.  If specified, will limit the recursive search,
-        but may result in some accessible equilibria not being found.
+        Maximum depth of recursion when using the strategic representation.
+        If specified, will limit the recursive search, but may result in some accessible
+        equilibria not being found.
 
     Returns
     -------
@@ -226,24 +227,32 @@ def lcp_solve(
     ------
     RuntimeError
         If game has more than two players.
+
+    ValueError
+        If stop_after or max_depth are supplied for use on the tree representation.
     """
-    if stop_after is None:
-        stop_after = 0
-    elif stop_after < 0:
+    if game.is_tree and not use_strategic:
+        if stop_after is not None:
+            raise ValueError(
+                "lcp_solve(): stop_after can only be used on the strategic representation"
+            )
+        if max_depth is not None:
+            raise ValueError(
+                "lcp_solve(): max_depth can only be used on the strategic representation"
+            )
+    if stop_after is not None and stop_after < 0:
         raise ValueError(
             f"lcp_solve(): stop_after argument must be a non-negative number; got {stop_after}"
         )
-    if max_depth is None:
-        max_depth = 0
     if not game.is_tree or use_strategic:
         if rational:
             equilibria = libgbt._lcp_strategy_solve_rational(game, stop_after or 0, max_depth or 0)
         else:
             equilibria = libgbt._lcp_strategy_solve_double(game, stop_after or 0, max_depth or 0)
     elif rational:
-        equilibria = libgbt._lcp_behavior_solve_rational(game, stop_after or 0, max_depth or 0)
+        equilibria = libgbt._lcp_behavior_solve_rational(game)
     else:
-        equilibria = libgbt._lcp_behavior_solve_double(game, stop_after or 0, max_depth or 0)
+        equilibria = libgbt._lcp_behavior_solve_double(game)
     return NashComputationResult(
         game=game,
         method="lcp",

src/solvers/lcp/efglcp.cc

@@ -29,9 +29,8 @@ namespace Gambit::Nash {
 
 template <class T> class NashLcpBehaviorSolver {
 public:
-  NashLcpBehaviorSolver(int p_stopAfter, int p_maxDepth,
-                        BehaviorCallbackType<T> p_onEquilibrium = NullBehaviorCallback<T>)
-    : m_onEquilibrium(p_onEquilibrium), m_stopAfter(p_stopAfter), m_maxDepth(p_maxDepth)
+  NashLcpBehaviorSolver(BehaviorCallbackType<T> p_onEquilibrium = NullBehaviorCallback<T>)
+    : m_onEquilibrium(p_onEquilibrium)
   {
   }
   ~NashLcpBehaviorSolver() = default;
@@ -40,7 +39,6 @@ template <class T> class NashLcpBehaviorSolver {
 
 private:
   BehaviorCallbackType<T> m_onEquilibrium;
-  int m_stopAfter, m_maxDepth;
 
   class Solution;
 
@@ -145,98 +143,23 @@ std::list<MixedBehaviorProfile<T>> NashLcpBehaviorSolver<T>::Solve(const Game &p
   solution.eps = tab.Epsilon();
 
   try {
-    if (m_stopAfter != 1) {
-      try {
-        AllLemke(p_game, solution.ns1 + solution.ns2 + 1, tab, 0, A, solution);
-      }
-      catch (EquilibriumLimitReached &) {
-        // Handle this silently; equilibria are recorded as found so no action needed
-      }
-    }
-    else {
-      tab.Pivot(solution.ns1 + solution.ns2 + 1, 0);
-      tab.SF_LCPPath(solution.ns1 + solution.ns2 + 1);
-      solution.AddBFS(tab);
-      Vector<T> sol(tab.MinRow(), tab.MaxRow());
-      tab.BasisVector(sol);
-      MixedBehaviorProfile<T> profile(p_game);
-      GetProfile(tab, profile, sol, p_game->GetRoot(), 1, 1, solution);
-      profile.UndefinedToCentroid();
-      solution.m_equilibria.push_back(profile);
-      this->m_onEquilibrium(profile, "NE");
-    }
+    tab.Pivot(solution.ns1 + solution.ns2 + 1, 0);
+    tab.SF_LCPPath(solution.ns1 + solution.ns2 + 1);
+    solution.AddBFS(tab);
+    Vector<T> sol(tab.MinRow(), tab.MaxRow());
+    tab.BasisVector(sol);
+    MixedBehaviorProfile<T> profile(p_game);
+    GetProfile(tab, profile, sol, p_game->GetRoot(), 1, 1, solution);
+    profile.UndefinedToCentroid();
+    solution.m_equilibria.push_back(profile);
+    this->m_onEquilibrium(profile, "NE");
   }
   catch (std::runtime_error &e) {
     std::cerr << "Error: " << e.what() << std::endl;
   }
   return solution.m_equilibria;
 }
 
-//
-// All_Lemke finds all accessible Nash equilibria by recursively
-// calling itself.  List maintains the list of basic variables
-// for the equilibria that have already been found.
-// From each new accessible equilibrium, it follows
-// all possible paths, adding any new equilibria to the List.
-//
-template <class T>
-void NashLcpBehaviorSolver<T>::AllLemke(const Game &p_game, int j, linalg::LemkeTableau<T> &B,
-                                        int depth, Matrix<T> &A, Solution &p_solution) const
-{
-  if (m_maxDepth != 0 && depth > m_maxDepth) {
-    return;
-  }
-
-  Vector<T> sol(B.MinRow(), B.MaxRow());
-  MixedBehaviorProfile<T> profile(p_game);
-
-  bool newsol = false;
-  for (int i = B.MinRow(); i <= B.MaxRow() && !newsol; i++) {
-    if (i == j) {
-      continue;
-    }
-
-    linalg::LemkeTableau<T> BCopy(B);
-    // Perturb tableau by a small number
-    A(i, 0) = static_cast<T>(-1) / static_cast<T>(1000);
-    BCopy.Refactor();
-
-    int missing;
-    if (depth == 0) {
-      BCopy.Pivot(j, 0);
-      missing = -j;
-    }
-    else {
-      missing = BCopy.SF_PivotIn(0);
-    }
-
-    newsol = false;
-
-    if (BCopy.SF_LCPPath(-missing) == 1) {
-      newsol = p_solution.AddBFS(BCopy);
-      BCopy.BasisVector(sol);
-      GetProfile(BCopy, profile, sol, p_game->GetRoot(), 1, 1, p_solution);
-      profile.UndefinedToCentroid();
-      if (newsol) {
-        m_onEquilibrium(profile, "NE");
-        p_solution.m_equilibria.push_back(profile);
-        if (m_stopAfter > 0 && p_solution.EquilibriumCount() >= m_stopAfter) {
-          throw EquilibriumLimitReached();
-        }
-      }
-    }
-    else {
-      // Dead end
-    }
-
-    A(i, 0) = static_cast<T>(-1);
-    if (newsol) {
-      BCopy.Refactor();
-      AllLemke(p_game, i, BCopy, depth + 1, A, p_solution);
-    }
-  }
-}
-
 template <class T>
 void NashLcpBehaviorSolver<T>::FillTableau(Matrix<T> &A, const GameNode &n, T prob, int s1, int s2,
                                            T payoff1, T payoff2, Solution &p_solution) const
@@ -342,16 +265,15 @@ void NashLcpBehaviorSolver<T>::GetProfile(const linalg::LemkeTableau<T> &tab,
 }
 
 template <class T>
-std::list<MixedBehaviorProfile<T>> LcpBehaviorSolve(const Game &p_game, int p_stopAfter,
-                                                    int p_maxDepth,
+std::list<MixedBehaviorProfile<T>> LcpBehaviorSolve(const Game &p_game,
                                                     BehaviorCallbackType<T> p_onEquilibrium)
 {
-  return NashLcpBehaviorSolver<T>(p_stopAfter, p_maxDepth, p_onEquilibrium).Solve(p_game);
+  return NashLcpBehaviorSolver<T>(p_onEquilibrium).Solve(p_game);
 }
 
-template std::list<MixedBehaviorProfile<double>> LcpBehaviorSolve(const Game &, int, int,
+template std::list<MixedBehaviorProfile<double>> LcpBehaviorSolve(const Game &,
                                                                   BehaviorCallbackType<double>);
 template std::list<MixedBehaviorProfile<Rational>>
-LcpBehaviorSolve(const Game &, int, int, BehaviorCallbackType<Rational>);
+LcpBehaviorSolve(const Game &, BehaviorCallbackType<Rational>);
 
 } // end namespace Gambit::Nash

src/solvers/lcp/lcp.h

@@ -34,7 +34,7 @@ LcpStrategySolve(const Game &p_game, int p_stopAfter, int p_maxDepth,
 
 template <class T>
 std::list<MixedBehaviorProfile<T>>
-LcpBehaviorSolve(const Game &p_game, int p_stopAfter, int p_maxDepth,
+LcpBehaviorSolve(const Game &p_game,
                  BehaviorCallbackType<T> p_onEquilibrium = NullBehaviorCallback<T>);
 
 } // end namespace Gambit::Nash

src/tools/lcp/lcp.cc

@@ -149,14 +149,14 @@ int main(int argc, char *argv[])
       if (useFloat) {
         auto renderer =
             MakeMixedBehaviorProfileRenderer<double>(std::cout, numDecimals, printDetail);
-        LcpBehaviorSolve<double>(game, stopAfter, maxDepth,
+        LcpBehaviorSolve<double>(game,
                                  [&](const MixedBehaviorProfile<double> &p,
                                      const std::string &label) { renderer->Render(p, label); });
       }
       else {
         auto renderer =
             MakeMixedBehaviorProfileRenderer<Rational>(std::cout, numDecimals, printDetail);
-        LcpBehaviorSolve<Rational>(game, stopAfter, maxDepth,
+        LcpBehaviorSolve<Rational>(game,
                                    [&](const MixedBehaviorProfile<Rational> &p,
                                        const std::string &label) { renderer->Render(p, label); });
       }