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/9] 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 4238a15559..333a8e6f54 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 2df8238e4b..5a07164af0 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/9] 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 25b06325c2..460dcc6973 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 0a6fcba301..916d55c81e 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 d05b1c9adf..43cf9aff12 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/9] 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 333a8e6f54..4238a15559 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/9] 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 5a07164af0..2df8238e4b 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/9] 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 916d55c81e..8003b10b9d 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}; } From 31df22b0f88f8620fb8a6b47af10cd6a6139b57d Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 7 Sep 2026 15:44:49 +0100 Subject: [PATCH 6/9] Added maxRegret parameter for HP --- src/gui/dlnash.cc | 6 ++++-- src/gui/nashspec.cc | 2 +- src/gui/nashspec.h | 1 + src/pygambit/cli/hp.py | 11 ++++++++++- src/pygambit/gambit.pxd | 2 +- src/pygambit/nash.h | 5 +++-- src/pygambit/nash.pxi | 3 ++- src/pygambit/nash.py | 12 ++++++++++-- src/solvers/hp/hp.cc | 18 ++++++++++-------- src/solvers/hp/hp.h | 2 +- 10 files changed, 43 insertions(+), 19 deletions(-) diff --git a/src/gui/dlnash.cc b/src/gui/dlnash.cc index 702b3b63d2..5572e7b1f8 100644 --- a/src/gui/dlnash.cc +++ b/src/gui/dlnash.cc @@ -175,7 +175,8 @@ wxString ExternalCommand(const NashComputationSpec &p_spec) method.localNewtonMaxIterations); } else if constexpr (std::is_same_v) { - return prefix + wxString::Format("hp -d 10 -n %d", method.priors); + return prefix + + wxString::Format("hp -d 10 -n %d -m %.17g", method.priors, method.maxRegret); } else if constexpr (std::is_same_v) { return prefix + wxString::Format("ipa -d 10 -n %d", method.perturbations); @@ -274,7 +275,8 @@ wxString ParameterDescription(const NashMethodSpec &p_method) method.localNewtonMaxIterations); } else if constexpr (std::is_same_v) { - return wxString::Format(" (%d random priors)", method.priors); + return wxString::Format(" (%d random priors; maximum regret %.4g)", method.priors, + method.maxRegret); } else if constexpr (std::is_same_v) { return wxString::Format(" (%d perturbation)", method.perturbations); diff --git a/src/gui/nashspec.cc b/src/gui/nashspec.cc index 77e5d17f31..54be81996f 100644 --- a/src/gui/nashspec.cc +++ b/src/gui/nashspec.cc @@ -104,7 +104,7 @@ std::optional HPNashSpec::MakeSolver(NashRepresentation) const for (const auto &prior : NewRandomStrategyProfiles(p_game, spec.priors)) { p_cancel.Check(); Nash::HPStrategySolve( - prior, + prior, spec.maxRegret, [&p_callback](const MixedStrategyProfile &p) { p_callback(ComputedProfile(p)); }, Nash::NullHPEventCallback, p_cancel); } diff --git a/src/gui/nashspec.h b/src/gui/nashspec.h index de07b4d933..d138cf9f8b 100644 --- a/src/gui/nashspec.h +++ b/src/gui/nashspec.h @@ -102,6 +102,7 @@ struct HPNashSpec { // path-tracing can surface more), so this defaults high, like LiapNashSpec's // startingPoints, rather than to 1. int priors{10}; + double maxRegret{1.0e-8}; std::optional MakeSolver(NashRepresentation) const; }; diff --git a/src/pygambit/cli/hp.py b/src/pygambit/cli/hp.py index 0daa664b61..41acc7890c 100644 --- a/src/pygambit/cli/hp.py +++ b/src/pygambit/cli/hp.py @@ -80,6 +80,14 @@ default=None, help="file containing prior distributions (mutually exclusive with -n)", ) +@click.option( + "-m", + "maxregret", + default=1.0e-8, + show_default=True, + type=float, + help="maximum regret acceptable as a proportion of the range of payoffs in the game", +) @click.option("-q", "--quiet", is_flag=True, help="quiet mode (suppresses banner)") @click.option( "-V", @@ -95,6 +103,7 @@ def main( n_priors: int | None, seed: int | None, start_file: str | None, + maxregret: float, quiet: bool, verbose: bool, ) -> None: @@ -109,7 +118,7 @@ def render_event(event: gbt.HPStepEvent) -> None: prior = prior.as_float() if verbose: click.echo(render_profile_csv(prior, "prior", decimals)) - result = gbt.nash.hp_solve(prior, event_callback=render_event) + result = gbt.nash.hp_solve(prior, maxregret=maxregret, event_callback=render_event) for eq in result.equilibria: click.echo(render_profile_csv(eq, "NE", decimals)) diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index a29b0f08e5..50a2cf0cae 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -701,5 +701,5 @@ cdef extern from "nash.h": LogitEventCallbackType[c_LogitQREMixedStrategyProfile] ) except + stdlist[c_MixedStrategyProfile[double]] HPStrategySolveWrapper( - c_MixedStrategyProfile[double], HPEventCallbackType + c_MixedStrategyProfile[double], double, HPEventCallbackType ) except +RuntimeError diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 75dbfb2e20..815a27d786 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -143,8 +143,9 @@ LogitStrategyEstimateWrapper(std::shared_ptr> p_fre } std::list> -HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, +HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, double p_maxRegret, Nash::HPEventCallbackType p_onEvent = Nash::NullHPEventCallback) { - return Nash::HPStrategySolve(p_prior, Nash::NullStrategyCallback, p_onEvent); + return Nash::HPStrategySolve(p_prior, p_maxRegret, Nash::NullStrategyCallback, + p_onEvent); } diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index e319846a3e..595a422235 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -914,8 +914,9 @@ def _logit_behavior_branch(game: Game, def _hp_strategy_solve( prior: MixedStrategyProfileDouble, + maxregret: float, event_callback: object = None, ) -> list[MixedStrategyProfileDouble]: return _convert_mspd(HPStrategySolveWrapper( - deref(prior.profile), MakeHPEventCallback(event_callback) + deref(prior.profile), maxregret, MakeHPEventCallback(event_callback) )) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 32aa3aa07c..8789079022 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -1113,6 +1113,7 @@ def logit_solve( def hp_solve( prior: libgbt.MixedStrategyProfileDouble, + maxregret: float = 1.0e-8, event_callback: Callable[[libgbt.HPStepEvent], None] | None = None, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -1125,6 +1126,11 @@ def hp_solve( prior : MixedStrategyProfileDouble The prior distribution over strategies. + 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 + event_callback : Callable[[HPStepEvent], None], optional If specified, called with each point traced along the homotopy path, and the homotopy parameter ``t`` at which it was reached, on the way @@ -1137,12 +1143,14 @@ def hp_solve( res : NashComputationResult The result represented as a ``NashComputationResult`` object. """ - equilibria = libgbt._hp_strategy_solve(prior, event_callback) + if maxregret <= 0.0: + raise ValueError("hp_solve(): maxregret argument must be positive") + equilibria = libgbt._hp_strategy_solve(prior, maxregret, event_callback) return NashComputationResult( game=prior.game, method="hp", rational=False, use_strategic=True, equilibria=equilibria, - parameters={}, + parameters={"maxregret": maxregret}, ) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 460dcc6973..658a873e65 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -27,7 +27,7 @@ namespace Gambit::Nash { std::list> -HPStrategySolve(const MixedStrategyProfile &p_prior, +HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret, StrategyCallbackType p_onEquilibrium, HPEventCallbackType p_onEvent, const CancelToken &p_cancel) { @@ -40,16 +40,16 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, const PathTracer::TraceDirection direction = PathTracer::TraceDirection::Positive; const size_t tracking_index = 1; // Track the first variable (t) for orientation const double t_target = 1.0; - const double tol = 1e-8; double last_t = 0.0; bool has_crossed = false; - auto termination_condition = [t_target, &last_t, &has_crossed, tol, + auto termination_condition = [t_target, &last_t, &has_crossed, p_maxRegret, &system](const Vector &point) { const double t = point[1]; - // Path tracer reaches tol - if (system.ExtractEquilibrium(point).GetMaxRegret() <= tol && t >= t_target - tol) { + // Path tracer reaches maximum acceptable regret + if (system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret && + t >= t_target - p_maxRegret) { return true; } @@ -57,7 +57,8 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, if (!has_crossed) { has_crossed = true; } - else if (t > last_t + tol) { // Criterion function is not working; polish will do the job + else if (t > + last_t + p_maxRegret) { // Criterion function is not working; polish will do the job return true; } } @@ -75,8 +76,9 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, return point[1] - t_target; }; - auto polishing_termination_condition = [tol, &system](const Vector &point) -> bool { - return system.ExtractEquilibrium(point).GetMaxRegret() <= tol; + auto polishing_termination_condition = [p_maxRegret, + &system](const Vector &point) -> bool { + return system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret; }; const auto tracing_result = tracer.TracePath( diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index 4bc5a7280f..b600a1e020 100644 --- a/src/solvers/hp/hp.h +++ b/src/solvers/hp/hp.h @@ -45,7 +45,7 @@ inline void NullHPEventCallback(const HPEvent &) {} /// @brief Compute a Nash equilibrium of a game using the homotopy method of /// Herings and Peeters (2001) std::list> -HPStrategySolve(const MixedStrategyProfile &p_prior, +HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret = 1.0e-8, StrategyCallbackType p_onEquilibrium = NullStrategyCallback, HPEventCallbackType p_onEvent = NullHPEventCallback, const CancelToken &p_cancel = CancelToken()); From f5a66f30043882449bbe6638c3b59763dcb4bd46 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 7 Sep 2026 15:49:21 +0100 Subject: [PATCH 7/9] Added maxRegret parameter for HP to manpage --- doc/tools.hp.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/tools.hp.rst b/doc/tools.hp.rst index b2eac367c4..86770fa328 100644 --- a/doc/tools.hp.rst +++ b/doc/tools.hp.rst @@ -24,6 +24,11 @@ different equilibria being found. Prints a help message listing the available options. +.. cmdoption:: -m + + Specify the maximum regret criterion for acceptance as an approximate Nash equilibrium + (default is 1e-8). See :ref:`pygambit-nash-maxregret` for interpretation and guidance. + .. cmdoption:: -n Randomly generate the specified number of prior distributions. From 8abbb1673df9c7c8aadaa5887d1da99f4dd586f9 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 7 Sep 2026 16:28:16 +0100 Subject: [PATCH 8/9] Return maxregret docs to main docs from tutorial --- doc/algorithms.rst | 36 ++++++++++++++++++++++ doc/tutorials/03_stripped_down_poker.ipynb | 24 +-------------- 2 files changed, 37 insertions(+), 23 deletions(-) diff --git a/doc/algorithms.rst b/doc/algorithms.rst index 317306b029..db71bc87d2 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -22,6 +22,42 @@ Algorithm Description :ref:`hp` Compute a specific Nash equilibrium using a homotopy path-following method :py:func:`pygambit.nash.hp_solve` :ref:`gambit-hp ` ================ =========================================================================== ======================================== ========================================== +.. _pygambit-nash-maxregret: + +Acceptance criteria for approximate Nash equilibria +---------------------------------------------------- + +Some methods for computing Nash equilibria operate using floating-point arithmetic and/or +generate candidate equilibrium profiles using methods which involve some form of successive +approximation. The outputs of these methods therefore are in general +:math:`\varepsilon`-equilibria, for some positive :math:`\varepsilon`: a strategy profile at +which no player can gain more than :math:`\varepsilon` in expected payoff by unilaterally +deviating. Every Nash equilibrium is an :math:`\varepsilon`-equilibrium with +:math:`\varepsilon = 0`. + +To provide a uniform interface across methods, where relevant Gambit provides a parameter +`maxregret`, which specifies the acceptance criterion for labeling the output of the +algorithm as an equilibrium. This parameter is interpreted *proportionally* to the range of +payoffs in the game: any profile returned as an equilibrium is guaranteed to be an +:math:`\varepsilon`-equilibrium, for :math:`\varepsilon` no more than `maxregret` times the +difference of the game's maximum and minimum payoffs. For example, with the default +`maxregret` of :math:`10^{-8}` in a game whose payoffs range over 4 units, any equilibrium +returned is guaranteed to have a regret, measured directly in the game's own payoffs, of no +more than :math:`4 \times 10^{-8}`. + +Expressing `maxregret` scaled by the game's payoffs in this way standardises the behavior of +methods across games: for instance, doubling all the payoffs in a game does not change the +`maxregret` value needed to obtain equilibria of comparable quality. + +Methods differ in the guarantees they offer once a `maxregret` criterion is specified. +Globally-convergent methods, such as :ref:`logit` and :ref:`gnm`, are guaranteed to eventually +satisfy any `maxregret` criterion, though a tighter criterion generally requires more +computation. Other methods, such as :ref:`liap`, are not globally convergent, and may fail to +find any equilibrium satisfying the criterion from a given starting point. + +See the :doc:`stripped-down poker tutorial ` for a worked +example comparing `maxregret` across several methods. + .. _enumpure: enumpure diff --git a/doc/tutorials/03_stripped_down_poker.ipynb b/doc/tutorials/03_stripped_down_poker.ipynb index fab307ac26..cc890eec5a 100644 --- a/doc/tutorials/03_stripped_down_poker.ipynb +++ b/doc/tutorials/03_stripped_down_poker.ipynb @@ -675,29 +675,7 @@ "cell_type": "markdown", "id": "b2867dca", "metadata": {}, - "source": [ - "Acceptance criteria for Nash equilibria\n", - "---------------------------------------\n", - "\n", - "Some methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations.\n", - "The outputs of these methods therefore are in general $\\varepsilon$-equilibria, for some positive $\\varepsilon$.\n", - "\n", - "$\\varepsilon$-equilibria (from [Wikipedia](https://en.wikipedia.org/wiki/Epsilon-equilibrium)):\n", - "\n", - "> In game theory, an epsilon-equilibrium, or near-Nash equilibrium, is a strategy profile that approximately satisfies the condition of Nash equilibrium. In a Nash equilibrium, no player has an incentive to change his behavior. In an approximate Nash equilibrium, this requirement is weakened to allow the possibility that a player may have a small incentive to do something different.\n", - "\n", - "> Given a game and a real non-negative parameter $\\varepsilon$, a strategy profile is said to be an $\\varepsilon$-equilibrium if it is not possible for any player to gain more than $\\varepsilon$ in expected payoff by unilaterally deviating from his strategy. Every Nash Equilibrium is an $\\varepsilon$-equilibrium where $\\varepsilon = 0$.\n", - "\n", - "\n", - "To provide a uniform interface across methods, where relevant Gambit provides a parameter\n", - "`maxregret`, which specifies the acceptance criterion for labeling the output of the\n", - "algorithm as an equilibrium.\n", - "This parameter is interpreted *proportionally* to the range of payoffs in the game.\n", - "Any profile returned as an equilibrium is guaranteed to be an $\\varepsilon$-equilibrium, for $\\varepsilon$ no more than `maxregret`\n", - "times the difference of the game's maximum and minimum payoffs.\n", - "\n", - "As an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" - ] + "source": "Acceptance criteria for Nash equilibria\n---------------------------------------\n\nSome methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations, and so return only approximate equilibria.\nGambit expresses how close an approximation must be, before it is reported as an equilibrium, via a `maxregret` parameter, interpreted *proportionally* to the range of payoffs in the game.\nSee [Acceptance criteria for approximate Nash equilibria](../algorithms.html#pygambit-nash-maxregret) in the algorithms documentation for the full explanation of `maxregret` and the guarantees it provides.\n\nAs an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" }, { "cell_type": "code", From 585cfb4f88e5d32550fc3cce9b04994a3a788fda Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 7 Sep 2026 16:44:26 +0100 Subject: [PATCH 9/9] Adjustment to wire maxregret through consistently --- src/solvers/hp/hp.cc | 3 +++ src/solvers/path/path.cc | 14 ++------------ 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 658a873e65..95e640b875 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -108,6 +108,9 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret, return {}; } const MixedStrategyProfile equilibrium = system.ExtractEquilibrium(x); + if (equilibrium.GetMaxRegret() > p_maxRegret) { + return {}; + } p_onEquilibrium(equilibrium); equilibria.push_back(equilibrium); return equilibria; diff --git a/src/solvers/path/path.cc b/src/solvers/path/path.cc index 8003b10b9d..7db29c9bbd 100644 --- a/src/solvers/path/path.cc +++ b/src/solvers/path/path.cc @@ -331,8 +331,6 @@ PolishResult PolishPoint(std::function &, 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) @@ -340,8 +338,6 @@ PolishResult PolishPoint(std::function &, Vector 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; @@ -391,14 +387,8 @@ PolishResult PolishPoint(std::function &, Vector 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}; - } - } - + // The loop only exits here once p_terminate(x) holds for the actual, current x; + // there is nothing further to validate against a separate tolerance. return {x, true, "Polishing terminated successfully.", steps}; }