diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 25b06325c..460dcc697 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -39,11 +39,47 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, const PathTracer tracer; const PathTracer::TraceDirection direction = PathTracer::TraceDirection::Positive; const size_t tracking_index = 1; // Track the first variable (t) for orientation - auto termination_condition = [](const Vector &point) { return point[1] >= 1.5; }; - auto criterion_function = [](const Vector &point, - const Vector &tangent) -> double { return point[1] - 1.0; }; + const double t_target = 1.0; + const double tol = 1e-8; + double last_t = 0.0; + bool has_crossed = false; - tracer.TracePath( + auto termination_condition = [t_target, &last_t, &has_crossed, tol, + &system](const Vector &point) { + const double t = point[1]; + + // Path tracer reaches tol + if (system.ExtractEquilibrium(point).GetMaxRegret() <= tol && t >= t_target - tol) { + return true; + } + + if (t > t_target) { + if (!has_crossed) { + has_crossed = true; + } + else if (t > last_t + tol) { // Criterion function is not working; polish will do the job + return true; + } + } + else { + // Criterion function might take t back to being less than t_target + has_crossed = false; + } + + last_t = t; + return false; + }; + + auto criterion_function = [t_target](const Vector &point, + const Vector &tangent) -> double { + return point[1] - t_target; + }; + + auto polishing_termination_condition = [tol, &system](const Vector &point) -> bool { + return system.ExtractEquilibrium(point).GetMaxRegret() <= tol; + }; + + const auto tracing_result = tracer.TracePath( [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, [&system](const Vector &point, Matrix &jac) { system.GetJacobian(point, jac); @@ -55,6 +91,20 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, }, criterion_function, NullCriterionBracketFunction, p_cancel); + const PolishResult polishing_result = PolishPoint( + [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, + [&system](const Vector &point, Matrix &jac) { + system.GetJacobian(point, jac); + }, + x, t_target, 1, polishing_termination_condition, 100, + [&system, &p_onEvent](const Vector &point) { + const MixedStrategyProfile profile = system.ExtractEquilibrium(point); + p_onEvent(HPStepEvent{.profile = profile, .t = point[1]}); + }); + + if (!polishing_result.status) { + return {}; + } const MixedStrategyProfile equilibrium = system.ExtractEquilibrium(x); p_onEquilibrium(equilibrium); equilibria.push_back(equilibrium); diff --git a/src/solvers/path/path.cc b/src/solvers/path/path.cc index 0a6fcba30..8003b10b9 100644 --- a/src/solvers/path/path.cc +++ b/src/solvers/path/path.cc @@ -132,14 +132,15 @@ TracePathResult PathTracer::TracePath( CallbackFunctionType p_callback, CriterionFunctionType p_criterion, CriterionBracketFunctionType p_criterionBracket, const CancelToken &p_cancel) const { - const double c_tol = 1.0e-4; // tolerance for corrector iteration - const double c_maxDist = 0.4; // maximal distance to curve - const double c_maxContr = 0.6; // maximal contraction rate in corrector - 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_tol = 1.0e-4; // tolerance for corrector iteration + const double c_maxDist = 0.4; // maximal distance to curve + const double c_maxContr = 0.6; // maximal contraction rate in corrector + 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_newtonTol = 1.0e-8; // tolerance for Newton convergence 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 @@ -147,6 +148,9 @@ TracePathResult PathTracer::TracePath( double pert_countdown = 0.0; // How much longer (in arclength) to apply perturbation const double c_orientTol = 1.0e-8; // tolerance for detecting change in orientation + const double b_tol = 1.0e-10; // Tolerance for perturbing the b matrix in case of singularity + const double b_pert = 1.0e-8; // Perturbation of the b matrix in case of singularity + Vector u(x.size()); // t is current tangent at x; newT is tangent at u, which is the next point. Vector t(x.size()), newT(x.size()); @@ -158,11 +162,22 @@ TracePathResult PathTracer::TracePath( QRDecomp(b, q); q.GetRow(q.NumRows(), t); p_callback(x); + int steps = 0; + + auto stepsizeBelowMinimum = [&]() -> TracePathResult { + if (newton && std::abs(p_criterion(x, t)) < c_newtonTol) { + return {x, true, + "Path following terminated successfully at point satisfying criterion function.", + steps}; + } + return {x, false, "Stepsize fell below minimum threshold.", steps}; + }; + bool first_step = true; double omega = (p_direction == TraceDirection::Positive) ? 1.0 : -1.0; if (p_trackingIndex > x.size() || p_trackingIndex < 1) { - return {x, false, "Tracking index exceeds dimension of point vector."}; + return {x, false, "Tracking index exceeds dimension of point vector.", steps}; } while (!p_terminate(x)) { @@ -171,16 +186,17 @@ TracePathResult PathTracer::TracePath( bool accept = true; if (std::abs(h) <= c_hmin) { - return {x, false, "Stepsize fell below minimum threshold."}; + return stepsizeBelowMinimum(); } if (first_step) { if (std::abs(t[p_trackingIndex]) <= c_orientTol) { - return {x, false, "Initial tangent vector is orthogonal to path-following direction."}; + return {x, false, "Initial tangent vector is orthogonal to path-following direction.", + steps}; } // Ensure that the tangent is oriented in the same direction as // the path-following direction. - else if (t[p_trackingIndex] < -c_orientTol) { + if (t[p_trackingIndex] < -c_orientTol) { omega *= -1.0; } first_step = false; @@ -195,6 +211,18 @@ TracePathResult PathTracer::TracePath( p_jacobian(u, b); QRDecomp(b, q); + // Perturb the b matrix if it is singular or nearly singular + for (size_t i = 1; i < b.NumRows(); i++) { + if (std::abs(b(i, i)) < b_tol) { + if (b(i, i) < 0) { + b(i, i) -= b_pert; + } + else { + b(i, i) += b_pert; + } + } + } + int iter = 1; double disto = 0.0; while (true) { @@ -226,7 +254,7 @@ TracePathResult PathTracer::TracePath( disto = dist; iter++; if (iter > c_maxIter) { - return {x, false, "Maximum iterations exceeded."}; + return {x, false, "Maximum iterations exceeded.", steps}; } } @@ -240,7 +268,7 @@ TracePathResult PathTracer::TracePath( // 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::abs(2 * h); } accept = false; } @@ -248,7 +276,7 @@ TracePathResult PathTracer::TracePath( if (!accept) { h /= m_maxDecel; // PC not accepted; change stepsize and retry if (std::abs(h) <= c_hmin) { - return {x, false, "Stepsize fell below minimum threshold."}; + return stepsizeBelowMinimum(); } continue; } @@ -266,9 +294,10 @@ TracePathResult PathTracer::TracePath( p_criterionBracket(x, u); } - if (newton) { + const double diff = p_criterion(u, newT) - p_criterion(x, t); + if (newton && std::abs(diff) > c_newtonTol) { // Newton-type steplength adaptation, secant method - h *= -p_criterion(u, newT) / (p_criterion(u, newT) - p_criterion(x, t)); + h *= -p_criterion(u, newT) / diff; } else { // Standard steplength adaptation @@ -279,18 +308,98 @@ TracePathResult PathTracer::TracePath( x = u; t = newT; p_callback(x); + steps++; 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); + pert_countdown -= std::abs(h); if (pert_countdown < 0.0) { pert = 0.0; pert_countdown = 0.0; } } } - return {x, true, "Path tracing terminated successfully."}; + return {x, true, "Path tracing terminated successfully.", steps}; +} + +PolishResult PolishPoint(std::function &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &x, double fixed_value, size_t fixed_index, + TerminationFunctionType p_terminate, int max_iter, + CallbackFunctionType p_callback) +{ + x[fixed_index] = fixed_value; + + const Vector original_x = x; + + const size_t N = x.size() - 1; + Vector y(N); // Equations results + Matrix jac_full(N + 1, N); // Full Jacobian matrix (N+1 unknowns, N equations) + Matrix jac_square(N, N); // Jacobian matrix with fixed_index row removed + Matrix Q(N, N); // Orthogonal matrix from QR decomposition + Vector x_reduced(N); // Reduced x vector with fixed_index removed + + double const eq_tol = 1e-2; + + int steps = 0; + double dist = 0.0; + + while (!p_terminate(x)) { + if (steps >= max_iter) { + return {x, false, "Polishing exceeded maximum iterations.", steps}; + } + + p_function(x, y); + p_jacobian(x, jac_full); + + size_t row_index = 1; + for (size_t i = 1; i <= N + 1; ++i) { // Newton step expects the transposed Jacobian + if (i != fixed_index) { + for (size_t j = 1; j <= N; ++j) { + jac_square(row_index, j) = jac_full(i, j); + } + row_index++; + } + } + + // Reduced x vector removing fixed_index + size_t temp_idx = 1; + for (size_t i = 1; i <= N + 1; ++i) { + if (i != fixed_index) { + x_reduced[temp_idx++] = x[i]; + } + } + + QRDecomp(jac_square, Q); + + // Solve jac_square * x_reduced = -y + NewtonStep(Q, jac_square, x_reduced, y, dist); + + // Update x, keeping fixed_index constant + temp_idx = 1; + for (size_t i = 1; i <= N + 1; ++i) { + if (i != fixed_index) { + x[i] = x_reduced[temp_idx++]; + } + } + + steps++; + + if (p_callback) { + p_callback(x); + } + } + + // Checking that the profile satisfies the system of equations + for (size_t i = 1; i <= N; ++i) { + if (std::abs(y[i]) > eq_tol || std::isnan(y[i]) || std::isinf(y[i])) { + x = original_x; + return {x, false, "Polishing converged to an invalid mathematical state. Reverted.", steps}; + } + } + + return {x, true, "Polishing terminated successfully.", steps}; } } // end namespace Gambit diff --git a/src/solvers/path/path.h b/src/solvers/path/path.h index d05b1c9ad..43cf9aff1 100644 --- a/src/solvers/path/path.h +++ b/src/solvers/path/path.h @@ -63,8 +63,15 @@ struct TracePathResult { Vector final_point; bool status; // true if path tracing terminated successfully, false if it terminated due to error std::string message; // error message if status is false + int steps; // Step at which the tracing terminated }; +struct PolishResult { + Vector final_point; + bool status; // true if polishing terminated successfully, false if it terminated due to error + std::string message; // error message if status is false + int steps; // Step at which the polishing terminated +}; // // This class implements a generic path-following algorithm for smooth curves. // It is based on the ideas and codes presented in Allgower and Georg's @@ -96,6 +103,15 @@ class PathTracer { double m_maxDecel{1.1}, m_hStart{0.03}; }; +// This function reduces the regret of a point that is close to an equilibrium that has been found +// by the path-following algorithm. Fixing the value of a component of the point, it uses a +// Newton-type method to find a nearby point with lower regret. +PolishResult PolishPoint(std::function &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &p_x, double fixed_value, size_t fixed_index, + TerminationFunctionType p_terminate, int max_iter = 100, + CallbackFunctionType p_callback = NullCallbackFunction); + } // end namespace Gambit #endif // GAMBIT_SOLVERS_LOGIT_PATH_H