Use perturbations in gambit-logit to avoid infinite loops

The previous behaviour of the numerical continuation used to trace the logit QRE correspondence
attempted to deal with simple bifurcations by accepting the reversal of the orientation of the curve
and continuing "straight through".  However, in some cases, depending on where the tracing landed,
this could lead to an infinite loop (where the tracing repeatedly jumped back and forth between
two branches), or reversing the tracing entirely (and ending up back at the starting point).

Instead, the continuation now uses a perturbation approach and attempts to follow a path which has the same orientation.

It is still possible for the tracing to get stuck if it is too close to a bifurcation point.
This results in early termination of the tracing but not an infinite loop, and in the event this occurs
a limiting QRE is not marked as an epsilon-equilibrium.

Closes #3.

ChangeLog

@@ -10,6 +10,10 @@
   problem where the method could return non-Nash output (#406)
 - `gambit-enumpoly` could get stuck in an infinite loop, and/or fail to report some equilibria,
   due to floating-point rounding/tolerance issues; this has been fixed on known cases (#198)
+- `gambit-logit` now uses perturbations to attempt to resolve correspondences that have
+  bifurcations, and instead tries always to follow a curve that has the same orientation.
+  This should eliminate cases in which tracing gets stuck in a loop or reverses itself
+  when encountering bifurcations (#3)
 
 ### Added
 - MixedStrategyProfile and MixedBehaviorProfile objects in pygambit can now be iterated in

src/pygambit/nash.pxi

@@ -33,7 +33,6 @@ def _convert_mspd(
 ) -> typing.List[MixedStrategyProfileDouble]:
     ret = []
     for i in range(inlist.Length()):
-        print(i)
         p = MixedStrategyProfileDouble()
         p.profile = copyitem_list_mspd(inlist, i+1)
         ret.append(p)

src/solvers/logit/efglogit.cc

@@ -316,7 +316,7 @@ namespace {
 bool RegretTerminationFunction(const Game &p_game, const Vector<double> &p_point, double p_regret)
 {
   if (p_point.back() < 0.0) {
-    return false;
+    return true;
   }
   MixedBehaviorProfile<double> profile(p_game);
   for (int i = 1; i < p_point.Length(); i++) {
@@ -350,6 +350,9 @@ AgentQREPathTracer::TraceAgentPath(const LogitQREMixedBehaviorProfile &p_start,
         return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret);
       },
       func);
+  if (!m_fullGraph && func.GetProfiles().back().GetProfile().GetMaxRegret() >= p_regret) {
+    return {};
+  }
   return func.GetProfiles();
 }
 

src/solvers/logit/efglogit.h

@@ -82,7 +82,9 @@ inline List<MixedBehaviorProfile<double>> LogitBehaviorSolve(const Game &p_game,
   auto result =
       tracer.TraceAgentPath(LogitQREMixedBehaviorProfile(p_game), ostream, p_epsilon, 1.0);
   auto ret = List<MixedBehaviorProfile<double>>();
-  ret.push_back(result[1].GetProfile());
+  if (!result.empty()) {
+    ret.push_back(result.back().GetProfile());
+  }
   return ret;
 }
 

src/solvers/logit/nfglogit.cc

@@ -214,7 +214,7 @@ namespace {
 bool RegretTerminationFunction(const Game &p_game, const Vector<double> &p_point, double p_regret)
 {
   if (p_point.back() < 0.0) {
-    return false;
+    return true;
   }
   MixedStrategyProfile<double> profile(p_game->NewMixedStrategyProfile(0.0));
   for (int i = 1; i < p_point.Length(); i++) {
@@ -247,6 +247,9 @@ StrategicQREPathTracer::TraceStrategicPath(const LogitQREMixedStrategyProfile &p
         return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret);
       },
       func);
+  if (!m_fullGraph && func.GetProfiles().back().GetProfile().GetMaxRegret() >= p_regret) {
+    return {};
+  }
   return func.GetProfiles();
 }
 

src/solvers/logit/nfglogit.h

@@ -103,16 +103,17 @@ class StrategicQREEstimator : public StrategicQREPathTracer {
   class CallbackFunction;
 };
 
-inline List<MixedStrategyProfile<double>> LogitStrategySolve(const Game &p_game,
-                                                             double p_regret = 0.0001)
+inline List<MixedStrategyProfile<double>> LogitStrategySolve(const Game &p_game, double p_regret)
 {
   StrategicQREPathTracer tracer;
   tracer.SetFullGraph(false);
   std::ostringstream ostream;
   auto result =
       tracer.TraceStrategicPath(LogitQREMixedStrategyProfile(p_game), ostream, p_regret, 1.0);
   auto ret = List<MixedStrategyProfile<double>>();
-  ret.push_back(result[1].GetProfile());
+  if (!result.empty()) {
+    ret.push_back(result.back().GetProfile());
+  }
   return ret;
 }
 

src/solvers/logit/path.cc

@@ -109,6 +109,15 @@ void NewtonStep(Matrix<double> &q, Matrix<double> &b, Vector<double> &u, Vector<
 //             PathTracer: Implementation of path-following engine
 //----------------------------------------------------------------------------
 
+// To handle possible (simple) bifurcations in the graph, TracePath detects a
+// change in orientation of the curve, and, if one is found, implements a
+// perturbation on the first equation.  This perturbation is maintained only
+// long enough to traverse past the apparent bifurcation.
+// Preliminary experience suggests this works fairly well (when applied to QRE).
+// It is possible for the path-following to land sufficiently close to the
+// bifurcation point that the tracing gets stuck there as it is not possible
+// to find a small enough step size to avoid stepping over the bifurcation
+// point.
 void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, double &p_omega,
                            TerminationFunctionType p_terminate, const CallbackFunction &p_callback,
                            const CriterionFunction &p_criterion) const
@@ -122,7 +131,10 @@ void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, do
   const double c_hmin = 1.0e-8;  // minimal stepsize
   const int c_maxIter = 100;     // maximum iterations in corrector
 
-  bool newton = false; // using Newton steplength (for zero-finding)
+  bool newton = false;             // using Newton steplength (for zero-finding)
+  const double c_pert = 0.0000001; // The size of perturbation to apply to avoid bifurcation traps
+  double pert = 0.0;               // The current version of the perturbation being applied
+  double pert_countdown = 0.0;     // How much longer (in arclength) to apply perturbation
 
   Vector<double> u(x.Length()), restart(x.Length());
   // t is current tangent at x; newT is tangent at u, which is the next point.
@@ -162,6 +174,7 @@ void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, do
       double dist;
 
       p_system.GetValue(u, y);
+      y[1] += pert;
       NewtonStep(q, b, u, y, dist);
 
       if (dist >= c_maxDist) {
@@ -195,6 +208,21 @@ void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, do
       }
     }
 
+    // Obtain the tangent at the next step
+    q.GetRow(q.NumRows(), newT);
+    double omega_flip = (t * newT < 0.0) ? -1.0 : 1.0;
+
+    if (omega_flip == -1.0) {
+      // The orientation of the curve has changed, indicating a bifurcation.
+      // Switch on perturbation and attempt to continue following the branch that
+      // is oriented in the same direction as we were originally following
+      if (pert_countdown == 0.0) {
+        pert = c_pert;
+        pert_countdown = abs(2 * h);
+      }
+      accept = false;
+    }
+
     if (!accept) {
       h /= m_maxDecel; // PC not accepted; change stepsize and retry
       if (fabs(h) <= c_hmin) {
@@ -205,24 +233,17 @@ void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, do
         }
         return;
       }
-
       continue;
     }
 
     // Determine new stepsize
-    if (decel > m_maxDecel) {
-      decel = m_maxDecel;
-    }
-
-    // Obtain the tangent at the next step
-    q.GetRow(q.NumRows(), newT);
+    decel = std::min(decel, m_maxDecel);
 
     // If we are at a bifurcation point, the orientation of the tangent
     // will flip.  This will confuse many criterion functions, especially
     // those which are using derivatives to maximize or minimize an objective.
     // This ensures the criterion function is called with both the old and
     // new tangent oriented in the same sense.
-    double omega_flip = (t * newT < 0.0) ? -1.0 : 1.0;
     if (!newton && p_criterion(x, t) * p_criterion(u, newT * omega_flip) < 0.0) {
       newton = true;
       restart = u;
@@ -239,16 +260,18 @@ void PathTracer::TracePath(const EquationSystem &p_system, Vector<double> &x, do
 
     // PC step was successful; update and iterate
     x = u;
+    t = newT;
     p_callback(x, false);
 
-    if (t * newT < 0.0) {
-      // Bifurcation detected; for now, just "jump over" and continue,
-      // taking into account the change in orientation of the curve.
-      // Someday, we need to do more here!
-
-      p_omega = -p_omega;
+    if (pert_countdown > 0.0) {
+      // If we are currently perturbing in the neighborhood of a bifurcation, check to see
+      // whether we think we are likely past it, and switch off if we are.
+      pert_countdown -= abs(h);
+      if (pert_countdown < 0.0) {
+        pert = 0.0;
+        pert_countdown = 0.0;
+      }
     }
-    t = newT;
   }
 
   // Cleanup after termination