From 142845689c934465f881b6c77e02e47fe2ae6446 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sat, 5 Sep 2026 14:48:13 +0200 Subject: [PATCH 1/5] Doc --- doc/algorithms.rst | 19 +++++++++++++++++++ doc/gui.nash.rst | 1 + 2 files changed, 20 insertions(+) diff --git a/doc/algorithms.rst b/doc/algorithms.rst index 4238a1555..333a8e6f5 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -19,6 +19,7 @@ Algorithm Description :ref:`simpdiv` Compute equilibria via simplicial subdivision :py:func:`pygambit.nash.simpdiv_solve` :ref:`gambit-simpdiv ` :ref:`ipa` Compute equilibria using iterated polymatrix approximation :py:func:`pygambit.nash.ipa_solve` :ref:`gambit-ipa ` :ref:`gnm` Compute equilibria using a global Newton method :py:func:`pygambit.nash.gnm_solve` :ref:`gambit-gnm ` +:ref:`hp` Compute a specific Nash equilibrium using a homotopy path-following method :py:func:`pygambit.nash.hp_solve` ================ =========================================================================== ======================================== ========================================== .. _enumpure: @@ -234,3 +235,21 @@ The algorithm takes as a parameter a mixed strategy profile. This profile is interpreted as defining a ray in the space of games. The profile must have the property that, for each player, the most frequently played strategy must be unique. + +.. _hp: + +hp +--- +Computes the Nash equilibrium selected by the tracing procedure +of Harsanyi and Selten using a homotopy path-following method. The algorithm +was first described by P. Jean-Jacques Herings and Ronald J.A.P. Peeters +in :cite:p:`HerPee01`. + +The algorithm takes as a parameter a mixed strategy profile, which acts as +the subjective prior beliefs of the players. +The profile must have the property that, for each player, +there must only exist one best response. + +For generic games, the algorithm converges to the unique Nash equilibrium selected by +the tracing procedure of Harsanyi and Selten. For non-generic games, the algorithm may +converge to a Nash equilibrium that is not selected by the tracing procedure. diff --git a/doc/gui.nash.rst b/doc/gui.nash.rst index 2df8238e4..5a07164af 100644 --- a/doc/gui.nash.rst +++ b/doc/gui.nash.rst @@ -87,6 +87,7 @@ Method Parameters used by the graphical interface ``ipa`` One random perturbation. ``gnm`` One random perturbation; ending lambda ``-10``; 100 steps per support cell; local Newton refinement every 3 steps, with at most 10 iterations. +``hp`` No method-specific parameters. ================ ============================================================================ For extensive games, there is an option of whether to use the From 7f786b1131ee69c534759597a2dd1f11f7e5f717 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sat, 5 Sep 2026 18:14:57 +0200 Subject: [PATCH 2/5] Provisional termination --- src/solvers/hp/hp.cc | 58 ++++++++++++++-- src/solvers/path/path.cc | 142 ++++++++++++++++++++++++++++++++++----- src/solvers/path/path.h | 16 +++++ 3 files changed, 195 insertions(+), 21 deletions(-) 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..916d55c81 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; } @@ -279,18 +307,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) { + 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 From 3c0c44ec9fc74871f395a5971377e64340b5b735 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Sat, 5 Sep 2026 18:40:49 +0200 Subject: [PATCH 3/5] Remove doc Removed the section on the homotopy path-following method for computing Nash equilibrium, including its description and references. --- doc/algorithms.rst | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/doc/algorithms.rst b/doc/algorithms.rst index 333a8e6f5..4238a1555 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -19,7 +19,6 @@ Algorithm Description :ref:`simpdiv` Compute equilibria via simplicial subdivision :py:func:`pygambit.nash.simpdiv_solve` :ref:`gambit-simpdiv ` :ref:`ipa` Compute equilibria using iterated polymatrix approximation :py:func:`pygambit.nash.ipa_solve` :ref:`gambit-ipa ` :ref:`gnm` Compute equilibria using a global Newton method :py:func:`pygambit.nash.gnm_solve` :ref:`gambit-gnm ` -:ref:`hp` Compute a specific Nash equilibrium using a homotopy path-following method :py:func:`pygambit.nash.hp_solve` ================ =========================================================================== ======================================== ========================================== .. _enumpure: @@ -235,21 +234,3 @@ The algorithm takes as a parameter a mixed strategy profile. This profile is interpreted as defining a ray in the space of games. The profile must have the property that, for each player, the most frequently played strategy must be unique. - -.. _hp: - -hp ---- -Computes the Nash equilibrium selected by the tracing procedure -of Harsanyi and Selten using a homotopy path-following method. The algorithm -was first described by P. Jean-Jacques Herings and Ronald J.A.P. Peeters -in :cite:p:`HerPee01`. - -The algorithm takes as a parameter a mixed strategy profile, which acts as -the subjective prior beliefs of the players. -The profile must have the property that, for each player, -there must only exist one best response. - -For generic games, the algorithm converges to the unique Nash equilibrium selected by -the tracing procedure of Harsanyi and Selten. For non-generic games, the algorithm may -converge to a Nash equilibrium that is not selected by the tracing procedure. From 1f68874dada401ae1a8f176162f99fdd6b3eb294 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Sat, 5 Sep 2026 18:41:19 +0200 Subject: [PATCH 4/5] Remove doc2 Removed method-specific parameters description for 'hp'. --- doc/gui.nash.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/doc/gui.nash.rst b/doc/gui.nash.rst index 5a07164af..2df8238e4 100644 --- a/doc/gui.nash.rst +++ b/doc/gui.nash.rst @@ -87,7 +87,6 @@ Method Parameters used by the graphical interface ``ipa`` One random perturbation. ``gnm`` One random perturbation; ending lambda ``-10``; 100 steps per support cell; local Newton refinement every 3 steps, with at most 10 iterations. -``hp`` No method-specific parameters. ================ ============================================================================ For extensive games, there is an option of whether to use the From b072705bc5f8d579733215ac57e43d30e0f27d8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sun, 6 Sep 2026 15:28:03 +0200 Subject: [PATCH 5/5] Random game solution --- src/solvers/path/path.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/solvers/path/path.cc b/src/solvers/path/path.cc index 916d55c81..8003b10b9 100644 --- a/src/solvers/path/path.cc +++ b/src/solvers/path/path.cc @@ -294,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 @@ -392,7 +393,7 @@ PolishResult PolishPoint(std::function &, Vector eq_tol) { + 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}; }