Fix #492: Optimization and Stability Enhancements in Logit´s path follower

src/solvers/logit/efglogit.cc

@@ -245,26 +245,47 @@ void EquationSystem::GetJacobian(const Vector<double> &p_point, Matrix<double> &
 class TracingCallbackFunction {
 public:
   TracingCallbackFunction(const Game &p_game, MixedBehaviorObserverFunctionType p_observer)
-    : m_game(p_game), m_observer(p_observer)
+    : m_game(p_game), m_observer(p_observer), m_best_lambda(-INFINITY), m_best_regret(INFINITY),
+    m_bestProfile_lambda(p_game),
+    m_bestProfile_regret(p_game)
   {
   }
   ~TracingCallbackFunction() = default;
 
   void AppendPoint(const Vector<double> &p_point);
   const std::list<LogitQREMixedBehaviorProfile> &GetProfiles() const { return m_profiles; }
+  double GetMinRegret() const { return m_best_regret; }
+  double GetMaxLambda() const { return m_best_lambda; }
+  LogitQREMixedBehaviorProfile GetBestProfileLambda() const { return m_bestProfile_lambda; }
+  LogitQREMixedBehaviorProfile GetBestProfileRegret() const { return m_bestProfile_regret; }
 
 private:
   Game m_game;
   MixedBehaviorObserverFunctionType m_observer;
   std::list<LogitQREMixedBehaviorProfile> m_profiles;
+  double m_best_lambda;
+  double m_best_regret;
+  LogitQREMixedBehaviorProfile m_bestProfile_lambda;
+  LogitQREMixedBehaviorProfile m_bestProfile_regret;
 };
 
 void TracingCallbackFunction::AppendPoint(const Vector<double> &p_point)
 {
   const MixedBehaviorProfile<double> profile(PointToProfile(m_game, p_point));
   m_profiles.emplace_back(profile, p_point.back(), 1.0);
   m_observer(m_profiles.back());
+  double lambda = m_profiles.back().GetLambda();
+  if (lambda > m_best_lambda) {
+    m_best_lambda = lambda;
+    m_bestProfile_lambda = m_profiles.back();
+  }
+  double regret = m_profiles.back().GetProfile().GetAgentMaxRegret();
+  if(regret < m_best_regret) {
+    m_best_regret = regret;
+    m_bestProfile_regret = m_profiles.back();
+  }
 }
+
 
 class EstimatorCallbackFunction {
 public:
@@ -338,7 +359,20 @@ LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret,
         return RegretTerminationFunction(game, p_point, p_regret);
       },
       [&callback](const Vector<double> &p_point) -> void { callback.AppendPoint(p_point); });
-  return callback.GetProfiles();
+
+      //Checking Nash Equilibrium has been reached
+      double regret = callback.GetMinRegret();
+      if(regret > p_regret) {               
+        std::cerr << "WARNING: Logit solver stopped at lambda = " << callback.GetBestProfileRegret().GetLambda()
+                  << " with regret = " << regret
+                  << " (Goal was regret = " << p_regret << ")" << std::endl;
+        std::cout << "Best profile found = "; 
+        p_observer(callback.GetBestProfileRegret());
+        std::cout << std::endl;   
+        return std::list<LogitQREMixedBehaviorProfile>();
+      }      
+
+      return callback.GetProfiles();
 }
 
 std::list<LogitQREMixedBehaviorProfile>
@@ -370,7 +404,18 @@ LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start,
           return x.back() - lam;
         });
     ret.push_back(callback.GetProfiles().back());
-  }
+    }
+  //Checking Nash Equilibrium has been reached
+  //In this case, we check that the lambda value is at least as large as the target lambda,
+  double lambda = callback.GetMaxLambda();
+  if(lambda < p_targetLambda.back())  {
+      std::cerr << "WARNING: Logit solver stopped at lambda = " << lambda
+                << " with regret = " << callback.GetBestProfileLambda().GetProfile().GetAgentMaxRegret()
+                << " (Goal was lambda = " << p_targetLambda.back() << ")" << std::endl;
+      std::cout << "Best profile found = "; p_observer(callback.GetBestProfileLambda());std::cout << std::endl;
+      return std::list<LogitQREMixedBehaviorProfile>();
+  }            
+
   return ret;
 }
 

src/solvers/logit/nfglogit.cc

@@ -284,18 +284,28 @@ void EquationSystem::GetJacobian(const Vector<double> &p_point, Matrix<double> &
 class TracingCallbackFunction {
 public:
   TracingCallbackFunction(const Game &p_game, const MixedStrategyObserverFunctionType &p_observer)
-    : m_game(p_game), m_observer(p_observer)
+    : m_game(p_game), m_observer(p_observer), m_best_lambda(-INFINITY), m_best_regret(INFINITY),
+      m_bestProfile_lambda(p_game),
+      m_bestProfile_regret(p_game)
   {
   }
   ~TracingCallbackFunction() = default;
 
   void AppendPoint(const Vector<double> &p_point);
   const std::list<LogitQREMixedStrategyProfile> &GetProfiles() const { return m_profiles; }
+  double GetMinRegret() const { return m_best_regret; }
+  double GetMaxLambda() const { return m_best_lambda; }
+  LogitQREMixedStrategyProfile GetBestProfileLambda() const { return m_bestProfile_lambda; }
+  LogitQREMixedStrategyProfile GetBestProfileRegret() const { return m_bestProfile_regret; }
 
 private:
   Game m_game;
   MixedStrategyObserverFunctionType m_observer;
   std::list<LogitQREMixedStrategyProfile> m_profiles;
+  double m_best_lambda;
+  double m_best_regret;
+  LogitQREMixedStrategyProfile m_bestProfile_lambda;
+  LogitQREMixedStrategyProfile m_bestProfile_regret;
 };
 
 void TracingCallbackFunction::AppendPoint(const Vector<double> &p_point)
@@ -304,6 +314,16 @@ void TracingCallbackFunction::AppendPoint(const Vector<double> &p_point)
   PointToProfile(profile, p_point);
   m_profiles.emplace_back(profile, p_point.back(), 1.0);
   m_observer(m_profiles.back());
+  double lambda = m_profiles.back().GetLambda();
+  if (lambda > m_best_lambda) {
+    m_best_lambda = lambda;
+    m_bestProfile_lambda = m_profiles.back();
+  }
+  double regret = m_profiles.back().GetProfile().GetMaxRegret();
+  if(regret < m_best_regret) {
+    m_best_regret = regret;
+    m_bestProfile_regret = m_profiles.back();
+  }
 }
 
 class EstimatorCallbackFunction {
@@ -372,10 +392,24 @@ LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret,
         system.GetJacobian(p_point, p_jac);
       },
       x, p_omega,
-      [p_start, p_regret](const Vector<double> &p_point) {
+      [p_start, p_regret](const Vector<double> &p_point) { //Termination
         return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret);
       },
       [&callback](const Vector<double> &p_point) -> void { callback.AppendPoint(p_point); });
+
+      //Checking Nash Equilibrium has been reached
+      double regret = callback.GetMinRegret();
+      if(regret > p_regret)  { 
+        std::cerr << "WARNING: Logit solver stopped at lambda = " << callback.GetBestProfileRegret().GetLambda()
+                  << " with max regret = " << regret
+                  << " (Goal was " << p_regret << ")" << std::endl;
+        std::cout << "Best profile found = "; 
+        p_observer(callback.GetBestProfileRegret());
+        std::cout << std::endl;
+        return std::list<LogitQREMixedStrategyProfile>();
+      }
+
+
   return callback.GetProfiles();
 }
 
@@ -409,6 +443,17 @@ LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start,
         });
     ret.push_back(callback.GetProfiles().back());
   }
+  //Checking Nash Equilibrium has been reached
+  double lambda = callback.GetMaxLambda();
+  if(lambda < p_targetLambda.back())  {
+    std::cerr << "WARNING: Logit solver stopped at lambda = " << lambda
+              << " with regret = " << callback.GetBestProfileLambda().GetProfile().GetMaxRegret()
+              << " (Goal was lambda = " << p_targetLambda.back() << ")" << std::endl;
+    std::cout << "Best profile found = "; 
+    p_observer(callback.GetBestProfileLambda());
+    std::cout << std::endl;
+    return std::list<LogitQREMixedStrategyProfile>();
+  }
   return ret;
 }
 
@@ -455,3 +500,4 @@ LogitStrategyEstimate(const MixedStrategyProfile<double> &p_frequencies, double
 }
 
 } // end namespace Gambit
+

src/solvers/logit/path.cc

@@ -110,6 +110,7 @@ void NewtonStep(Matrix<double> &q, Matrix<double> &b, Vector<double> &u, Vector<
   d = std::sqrt(d);
 }
 
+
 } // end anonymous namespace
 
 //----------------------------------------------------------------------------
@@ -125,6 +126,8 @@ void NewtonStep(Matrix<double> &q, Matrix<double> &b, Vector<double> &u, Vector<
 // 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(
     std::function<void(const Vector<double> &, Vector<double> &)> p_function,
     std::function<void(const Vector<double> &, Matrix<double> &)> p_jacobian, Vector<double> &x,
@@ -137,14 +140,15 @@ void PathTracer::TracePath(
   const double c_eta = 0.1;      // perturbation to avoid cancellation
                                  // in calculating contraction rate
   double h = m_hStart;           // initial stepsize
-  const double c_hmin = 1.0e-8;  // minimal stepsize
-  const int c_maxIter = 100;     // maximum iterations in corrector
+  const double c_hmin = 1.0e-11;  // minimal stepsize
+  const int c_maxIter = 400;     // maximum iterations in corrector
 
   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
+  const double c_pert = 0.0001; // 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
-
+  const double min_pert_countdown = 0.05; // Minimum amount of perturbation to apply.
+  double min_lambda = -1e-6;;     // Minimum value for lambda when in previous iteration it was positive. 
   Vector<double> u(x.size());
   // t is current tangent at x; newT is tangent at u, which is the next point.
   Vector<double> t(x.size()), newT(x.size());
@@ -179,11 +183,16 @@ void PathTracer::TracePath(
       double dist;
 
       p_function(u, y);
-      y[1] += pert;
+      if (pert != 0.0) {
+        for (size_t i = 1; i <= y.size(); i++) {
+            // Symmetry breaking, perturbing all directions with an altenating sign
+            y[i] += pert * (i % 2 == 0 ? 1.0 : -1.0); 
+        }
+      }
       NewtonStep(q, b, u, y, dist);
 
       if (dist >= c_maxDist) {
-        accept = false;
+        accept = false;   //H(u,y) is too far from zero; reject PC step and reduce stepsize
         break;
       }
 
@@ -204,6 +213,8 @@ void PathTracer::TracePath(
       disto = dist;
       iter++;
       if (iter > c_maxIter) {
+        std::cerr << "Warning: corrector failed to converge after " << c_maxIter
+             << " iterations; rejecting predictor step." << std::endl;
         return;
       }
     }
@@ -213,20 +224,26 @@ void PathTracer::TracePath(
     const 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);
+        pert_countdown = std::max(fabs(10.0 * h), min_pert_countdown);
       }
-      accept = false;
     }
 
+
+    // If lambda was positive in the previous iteration, and we are now in the region of negative lambda, 
+    // we are likely heading towards the wrong direction, so we reject this step and reduce the stepsize.     
+    if (u.back() < min_lambda && x.back() >= 0.0) 
+      accept = false; 
+
     if (!accept) {
       h /= m_maxDecel; // PC not accepted; change stepsize and retry
       if (fabs(h) <= c_hmin) {
-        return;
+      return;
       }
       continue;
     }
@@ -271,3 +288,5 @@ void PathTracer::TracePath(
 }
 
 } // end namespace Gambit
+
+

src/solvers/logit/path.h

@@ -83,7 +83,7 @@ class PathTracer {
             CriterionBracketFunctionType p_criterionBracker = NullCriterionBracketFunction) const;
 
 private:
-  double m_maxDecel{1.1}, m_hStart{0.03};
+  double m_maxDecel{1.1}, m_hStart{0.1};
 };
 
 } // end namespace Gambit

src/tools/logit/logit.cc

@@ -85,6 +85,7 @@ bool ReadProfile(std::istream &p_stream, MixedStrategyProfile<double> &p_profile
 template <class T>
 void PrintProfile(std::ostream &p_stream, int p_decimals, const T &p_profile, bool p_nash = false)
 {
+
   if (p_nash) {
     p_stream << "NE";
   }
@@ -113,7 +114,7 @@ int main(int argc, char *argv[])
   double maxregret = 1.0e-8;
   std::string mleFile;
   double maxDecel = 1.1;
-  double hStart = 0.03;
+  double hStart = 0.1;
   std::list<double> targetLambda;
   bool fullGraph = true;
   int decimals = 6;
@@ -195,6 +196,19 @@ int main(int argc, char *argv[])
           "Computing equilibria of games with imperfect recall is not supported.");
     }
 
+    if (game->IsTree()) {
+      bool has_decisions = false;
+      for (auto player : game->GetPlayers()) {
+        if (player->GetInfosets().size() > 0 && !has_decisions) {
+          has_decisions = true;
+        }
+      }        
+      if (!has_decisions) {
+        std::cout << "ERROR: The game does not contain decision nodes for the players (only chance). No equilibria to calculate." << std::endl;
+        return 0;
+      }
+    }
+
     if (!mleFile.empty() && (!game->IsTree() || useStrategic)) {
       MixedStrategyProfile<double> frequencies(game->NewMixedStrategyProfile(0.0));
       std::ifstream mleData(mleFile.c_str());
@@ -227,7 +241,8 @@ int main(int argc, char *argv[])
       }
       else {
         auto result = LogitStrategySolve(start, maxregret, 1.0, hStart, maxDecel, printer);
-        PrintProfile(std::cout, decimals, result.back(), true);
+        if(!result.empty()) 
+          PrintProfile(std::cout, decimals, result.back(), true);
       }
     }
     else {
@@ -246,7 +261,8 @@ int main(int argc, char *argv[])
       }
       else {
         auto result = LogitBehaviorSolve(start, maxregret, 1.0, hStart, maxDecel, printer);
-        PrintProfile(std::cout, decimals, result.back(), true);
+        if(!result.empty()) 
+          PrintProfile(std::cout, decimals, result.back(), true);
       }
     }
     return 0;