Improve numerical stability and robustness of enumpoly

ChangeLog

@@ -8,6 +8,8 @@
 - Added `-R` option to `gambit-liap`, `gambit-simpdiv`, `gambit-ipa`, and `gambit-gnm` to reproducibly seed
   random starting point generation. (#448)
 - Added `Game.relabel_players`, which simultaneously reassigns the labels of the game's players.  (#1058)
+- Added `max_rectangles` to `enumpoly_solve` (and `-r` to `gambit-enumpoly`), bounding the number
+  of cells examined when searching for equilibria on a single support. (#1055)
 
 ### Fixed
 - `Game.add_player`, on a game in strategic form, no longer leaves the new player's sole
@@ -17,6 +19,7 @@
 - `enumpoly_solve` correctly handles failed perturbations at singular Jacobians (#591)
 - `enumpoly_solve` iterates sequences of actions in a deterministic order, pre-empting the
   possibility of equilibria being found in orders that differ across runs. (#1016)
+- Improved the numerical stability of `enumpoly_solve` handling less-well-conditioned systems.  (#1055)
 
 ### Changed
 - All long-running calculations are now done in a thread in the GUI process, rather than calling the

Makefile.am

@@ -352,6 +352,8 @@ enumpoly_SOURCES = \
 	${nashsupport_SOURCES} \
 	src/solvers/enumpoly/indexproduct.h \
 	src/solvers/enumpoly/rectangle.h \
+	src/solvers/enumpoly/lufactor.cc \
+	src/solvers/enumpoly/lufactor.h \
 	src/solvers/enumpoly/poly.cc \
 	src/solvers/enumpoly/poly.h \
 	src/solvers/enumpoly/polysystem.h \

doc/algorithms.rst

@@ -81,6 +81,17 @@ Subsequent analysis of unreached information sets can yield alternative
 profiles that specify different choices at unreached information sets
 while still satisfying the Nash equilibrium conditions.
 
+The subdivision process for a given support is guaranteed to terminate, because
+subdivision stops once a cell shrinks below a minimum size.  However, this bound
+is on the *depth* of subdivision, not on the *number* of cells examined to reach
+that depth.  As protection against pathological cases, the number of cells
+examined per support is capped (via ``max_rectangles`` in
+:py:func:`pygambit.nash.enumpoly_solve`, or :option:`gambit-enumpoly -r` for
+:program:`gambit-enumpoly`); if this limit is reached before the search of a
+support is complete, that support may be reported with incomplete equilibria,
+and (in verbose mode) is identified with the label "budget-exceeded" -- see
+:ref:`gambit-enumpoly <gambit-enumpoly>` for details.
+
 .. _lp:
 
 lp

doc/tools.enumpoly.rst

@@ -23,6 +23,15 @@ share the listed support.  However, the converse is not true:
 not all supports labeled as "singular" will necessarily be the
 support of some set of equilibria.
 
+Independently of degeneracy, the search on a single support can require
+examining an impractically large number of subdivided cells before it
+can either confirm or rule out an equilibrium.  To bound the time spent
+on any one support, the number of cells examined is capped (see
+:option:`gambit-enumpoly -r`).  When this limit is reached, the search for that support
+is abandoned, any equilibria already confirmed on it are still reported,
+and (in verbose mode) the support is identified with the label
+"budget-exceeded,".
+
 .. program:: gambit-enumpoly
 
 .. cmdoption:: -d
@@ -55,6 +64,14 @@ support of some set of equilibria.
    By default, the program will search all support profiles.
    This switch instructs the program to terminate when EQA equilibria have been found.
 
+.. cmdoption:: -r RECTANGLES
+
+   .. versionadded:: 17.0.0
+
+   Specify the maximum number of cells to examine when searching for equilibria
+   on a single support, before giving up on that support (default 20000).  See
+   the :ref:`algorithm description <enumpoly>` for why this is necessary.
+
 .. cmdoption:: -q
 
    Suppresses printing of the banner at program launch.

src/core/matrix.h

@@ -41,10 +41,6 @@ class SingularMatrixException final : public std::runtime_error {
 ///
 /// @tparam T Scalar element type. Must support arithmetic operations,
 ///           comparison with zero, and abs(T).
-///
-/// @note Inverse and Determinant use legacy Gaussian elimination algorithms.
-///       They are known not to be numerically optimal for near-singular matrices.
-///       Current behaviour is temporarily preserved for historical compatibility.
 template <class T> class Matrix {
   RectArray<T> m_data;
 
@@ -205,8 +201,6 @@ template <class T> class Matrix {
   /// @name Other operations
   ///@{
   Matrix Transpose() const;
-  Matrix Inverse() const;
-  T Determinant() const;
   ///@}
 };
 
@@ -412,158 +406,6 @@ template <class T> Matrix<T> Matrix<T>::Transpose() const
   return tmp;
 }
 
-// ----------------------------------------------------------------------------
-// Implementation of additional operations
-// ----------------------------------------------------------------------------
-
-template <class T> Matrix<T> Matrix<T>::Inverse() const
-{
-  if (!IsSquare()) {
-    throw DimensionException();
-  }
-  const int rmin = MinRow();
-  const int rmax = MaxRow();
-  const int cmin = MinCol();
-  const int cmax = MaxCol();
-  using Gambit::abs;
-
-  Matrix copy(*this);
-  Matrix inv(rmin, rmax, cmin, cmax);
-
-  inv = T{0};
-
-  // initialize inverse matrix and prescale row vectors
-  for (int i = rmin; i <= rmax; ++i) {
-    auto copy_row = copy.m_data.GetRowView(i);
-    auto inv_row = inv.m_data.GetRowView(i);
-
-    T max = maximize_function(copy_row, [](const T &v) { return abs(v); });
-    if (max == T{0}) {
-      throw SingularMatrixException();
-    }
-    const T scale = T{1} / max;
-    for (auto &v : copy_row) {
-      v *= scale;
-    }
-    inv_row[i] = scale;
-  }
-
-  for (int i = cmin; i <= cmax; ++i) {
-    // find pivot row
-    auto col_i = copy.m_data.GetColumnView(i);
-    T max = abs(col_i[i]);
-    int row = i;
-    for (int j = i + 1; j <= rmax; ++j) {
-      const T v = abs(col_i[j]);
-      if (v > max) {
-        max = v;
-        row = j;
-      }
-    }
-
-    if (max <= T{0}) {
-      throw SingularMatrixException();
-    }
-
-    copy.SwitchRows(i, row);
-    inv.SwitchRows(i, row);
-    // scale pivot row
-    T factor = T{1} / copy(i, i);
-    auto copy_row = copy.m_data.GetRowView(i);
-    auto inv_row = inv.m_data.GetRowView(i);
-    auto copy_it = copy_row.begin();
-    auto inv_it = inv_row.begin();
-    for (; copy_it != copy_row.end(); ++copy_it, ++inv_it) {
-      *copy_it *= factor;
-      *inv_it *= factor;
-    }
-
-    // reduce other rows
-    auto pivot_copy_row = copy.m_data.GetRowView(i);
-    auto pivot_inv_row = inv.m_data.GetRowView(i);
-
-    for (int j = rmin; j <= rmax; ++j) {
-      if (j == i) {
-        continue;
-      }
-      auto row_copy = copy.m_data.GetRowView(j);
-      auto row_inv = inv.m_data.GetRowView(j);
-      const T mult = row_copy[i];
-      auto pivot_copy_it = pivot_copy_row.begin();
-      auto pivot_inv_it = pivot_inv_row.begin();
-      auto row_copy_it = row_copy.begin();
-      auto row_inv_it = row_inv.begin();
-
-      for (; pivot_copy_it != pivot_copy_row.end();
-           ++pivot_copy_it, ++pivot_inv_it, ++row_copy_it, ++row_inv_it) {
-        *row_copy_it -= (*pivot_copy_it) * mult;
-        *row_inv_it -= (*pivot_inv_it) * mult;
-      }
-    }
-  }
-
-  return inv;
-}
-
-template <class T> T Matrix<T>::Determinant() const
-{
-  if (!IsSquare()) {
-    throw DimensionException();
-  }
-  const int rmin = MinRow();
-  const int rmax = MaxRow();
-  using Gambit::abs;
-
-  Matrix M(*this);
-
-  for (int row = rmin; row <= rmax; ++row) {
-    // Experience (as of 3/22/99) suggests that, in the interest of
-    // numerical stability, it might be best to do Gaussian
-    // elimination with respect to the row (of those feasible)
-    // whose entry has the largest absolute value.
-    int swap_row = row;
-    T max = abs(M(row, row));
-    for (int i = row + 1; i <= rmax; ++i) {
-      const T v = abs(M(i, row));
-      if (v > max) {
-        max = v;
-        swap_row = i;
-      }
-    }
-
-    if (swap_row != row) {
-      M.SwitchRows(row, swap_row);
-      auto pivot_row = M.m_data.GetRowView(row);
-      for (auto &v : pivot_row) {
-        v = -v;
-      }
-    }
-
-    if (M(row, row) == T{0}) {
-      return T{0};
-    }
-    // now do row operations to clear the row'th column
-    // below the diagonal
-    auto pivot_row = M.m_data.GetRowView(row);
-    for (int row1 = row + 1; row1 <= rmax; ++row1) {
-      auto elim_row = M.m_data.GetRowView(row1);
-      const T factor = -elim_row[row] / pivot_row[row];
-      auto pivot_it = pivot_row.begin();
-      auto elim_it = elim_row.begin();
-      for (; pivot_it != pivot_row.end(); ++pivot_it, ++elim_it) {
-        *elim_it += (*pivot_it) * factor;
-      }
-    }
-  }
-
-  // finally we multiply the diagonal elements
-  T det = T{1};
-  for (int row = rmin; row <= rmax; ++row) {
-    det *= M(row, row);
-  }
-  return det;
-}
-
 // ----------------------------------------------------------------------------
 // Implementation of operators
 // ----------------------------------------------------------------------------

src/gui/nashspec.cc

@@ -65,15 +65,15 @@ EnumPolyNashSpec::MakeSolver(NashRepresentation p_representation) const
     return [spec](const Game &p_game, const ProfileFoundCallback &p_callback,
                   const CancelToken &p_cancel) {
       Nash::EnumPolyBehaviorSolve(
-          p_game, spec.stopAfter, spec.maxRegret,
+          p_game, spec.stopAfter, spec.maxRegret, Nash::kDefaultEnumPolyMaxRectangles,
           [&p_callback](const MixedBehaviorProfile<double> &p) { p_callback(ComputedProfile(p)); },
           Nash::NullEnumPolyEventCallback<BehaviorSupportProfile>, p_cancel);
     };
   }
   return [spec](const Game &p_game, const ProfileFoundCallback &p_callback,
                 const CancelToken &p_cancel) {
     Nash::EnumPolyStrategySolve(
-        p_game, spec.stopAfter, spec.maxRegret,
+        p_game, spec.stopAfter, spec.maxRegret, Nash::kDefaultEnumPolyMaxRectangles,
         [&p_callback](const MixedStrategyProfile<double> &p) { p_callback(ComputedProfile(p)); },
         Nash::NullEnumPolyEventCallback<StrategySupportProfile>, p_cancel);
   };

src/pygambit/gambit.pxd

@@ -594,10 +594,10 @@ cdef extern from "solvers/nashsupport/nashsupport.h":
 
 cdef extern from "solvers/enumpoly/enumpoly.h":
     stdlist[c_MixedStrategyProfile[double]] EnumPolyStrategySolve(
-            c_Game, int, float
+            c_Game, int, float, size_t
     ) except +RuntimeError
     stdlist[c_MixedBehaviorProfile[double]] EnumPolyBehaviorSolve(
-            c_Game, int, float
+            c_Game, int, float, size_t
     ) except +RuntimeError
 
 cdef extern from "solvers/logit/logit.h":

src/pygambit/nash.pxi

@@ -180,16 +180,18 @@ def _enumpoly_strategy_solve(
         game: Game,
         stop_after: int,
         maxregret: float,
+        max_rectangles: int,
 ) -> list[MixedStrategyProfileDouble]:
-    return _convert_mspd(EnumPolyStrategySolve(game.game, stop_after, maxregret))
+    return _convert_mspd(EnumPolyStrategySolve(game.game, stop_after, maxregret, max_rectangles))
 
 
 def _enumpoly_behavior_solve(
         game: Game,
         stop_after: int,
         maxregret: float,
+        max_rectangles: int,
 ) -> list[MixedBehaviorProfileDouble]:
-    return _convert_mbpd(EnumPolyBehaviorSolve(game.game, stop_after, maxregret))
+    return _convert_mbpd(EnumPolyBehaviorSolve(game.game, stop_after, maxregret, max_rectangles))
 
 
 def _logit_strategy_solve(

src/pygambit/nash.py

@@ -662,7 +662,8 @@ def enumpoly_solve(
         game: libgbt.Game,
         use_strategic: bool = False,
         stop_after: int | None = None,
-        maxregret: float = 1.0e-4,
+        maxregret: float = 1.0e-8,
+        max_rectangles: int = 20_000,
         phcpack_path: pathlib.Path | str | None = None
 ) -> NashComputationResult:
     """:ref:`Compute Nash equilibria by enumerating all support profiles
@@ -682,11 +683,18 @@ def enumpoly_solve(
         Maximum number of equilibria to compute.  If not specified, examines
         all support profiles of the game.
 
-    maxregret : float, default 1e-4
+    maxregret : float, default 1e-8
         The acceptance criterion for approximate Nash equilibrium; the maximum
         regret of any player must be no more than `maxregret` times the
         difference of the maximum and minimum payoffs of the game
 
+    max_rectangles : int, default 20000
+        The maximum number of rectangles to examine when searching for equilibria on a
+        single support, before giving up on that support.  See the
+        :ref:`algorithm description <enumpoly>` for why this is necessary.
+
+        .. versionadded:: 17.0.0
+
     phcpack_path : str or pathlib.Path, optional
         If specified, use PHCpack to solve the systems of equations.
         This argument specifies the path to the PHCpack executable.
@@ -712,6 +720,11 @@ def enumpoly_solve(
         raise ValueError(
             f"enumpoly_solve(): maxregret must be a positive number; got {maxregret}"
         )
+    if max_rectangles <= 0:
+        raise ValueError(
+            f"enumpoly_solve(): "
+            f"max_rectangles argument must be a positive number; got {max_rectangles}"
+        )
     if phcpack_path is not None:
         if game.is_tree and not use_strategic:
             raise ValueError(
@@ -730,15 +743,16 @@ def enumpoly_solve(
         )
 
     if not game.is_tree or use_strategic:
-        equilibria = libgbt._enumpoly_strategy_solve(game, stop_after, maxregret)
+        equilibria = libgbt._enumpoly_strategy_solve(game, stop_after, maxregret, max_rectangles)
     else:
-        equilibria = libgbt._enumpoly_behavior_solve(game, stop_after, maxregret)
+        equilibria = libgbt._enumpoly_behavior_solve(game, stop_after, maxregret, max_rectangles)
     return NashComputationResult(
         game=game,
         method="enumpoly",
         rational=False,
         use_strategic=not game.is_tree or use_strategic,
-        parameters={"stop_after": stop_after, "maxregret": maxregret},
+        parameters={"stop_after": stop_after, "maxregret": maxregret,
+                    "max_rectangles": max_rectangles},
         equilibria=equilibria,
     )