Learning Approximate Solutions to Multiparametric Generalized Nash Equilibrium Problems

Abstract: We propose a learning-based approach for approximating solution mappings of multiparametric generalized Nash equilibrium problems (GNEPs) with coupling in both objectives and constraints. Rather than solving a standard regression problem on a training dataset of GNEP solutions, which are expensive and possibly difficult to collect, we use the Nikaido-Isoda (NI) gap function as a training loss, which requires only best-response data. To avoid bilevel optimization, a value-function surrogate approximates each agent's optimal best-response cost and is substituted into the NI loss, yielding a single-level learning problem. Learning approximate solutions to standard multiparametric programming problems is a special case of the approach. We also establish new sufficient conditions for the existence of a continuous parametric variational GNE selection under a strong variational stability assumption that generalizes strong monotonicity. The trained neural network delivers approximate GNE solutions with speedups of several orders of magnitude over online solvers. Numerical experiments on different problem classes confirm the effectiveness of the approach. A Python library is available at https://github.com/bemporad/mpfit.

• The paper introduces a two-stage machine learning framework that uses best-response value data, Nikaido–Isoda gap losses, and constraint penalties to approximate GNE solution mappings without solving exact equilibria for every parameter.
• The method achieves microsecond-scale inference with low best-response errors across linear-quadratic, QCQP, and nonlinear games, including non-monotone settings, while revealing accuracy and feasibility trade-offs controlled by penalty weights.
• The paper proves a sufficient continuity condition for variational GNE selections under strong variational stability and shows that the approach can approximate standard multiparametric programs with substantial speedups over exact explicit solvers.

The paper develops a machine learning framework for approximating solution mappings of multiparametric generalized Nash equilibrium problems (GNEPs) with coupling in both objective functions and shared constraints (2605.28757). The central contribution is a training methodology that avoids both the expensive collection of exact GNE solutions and the bilevel optimization structure inherent in equilibrium formulations, replacing them with a single-level learning problem built on the Nikaido–Isoda (NI) gap function and learned value-function surrogates.

Problem setting and motivation

The authors consider $N$ agents, each minimizing a cost $J_i(x_i, x_{-i}, p)$ subject to shared inequality and equality constraints $g(x,p)\le 0$, $h(x,p)=0$, where $p\in\mathcal{P}$ is a parameter vector. A generalized Nash equilibrium $x^\star(p)$ satisfies the fixed-point condition that each agent plays its best response given the others. Computing such equilibria online is prohibitive in game-theoretic MPC and receding-horizon settings, particularly when structural assumptions such as monotonicity—required by provably efficient distributed methods—do not hold. Exact explicit solution maps are known only for restricted classes such as linear-quadratic (LQ) GNEPs, and they suffer from exponential growth in the number of polyhedral critical regions, limiting them to small-scale problems. The stated goal is therefore to learn an approximate map $\hat{x}(p)\approx x^\star(p)$ without imposing convexity or monotonicity on the underlying game.

Learning via the Nikaido–Isoda gap function

Rather than performing regression on a dataset of GNEs—which requires solving a full GNEP for every training parameter and is ambiguous when equilibria are non-unique—the method exploits the fact that GNEs are exactly the zero points of the NI gap function over the shared feasible set. Training data consist only of tuples $(x_k, p_k, \bar{J}_k)$, where $\bar{J}_{k,i}$ is the optimal value of agent $i$'s best-response problem at $J_i(x_i, x_{-i}, p)$0, $J_i(x_i, x_{-i}, p)$1; best responses are substantially cheaper to compute than equilibria.

A direct NI-based training objective would be bilevel, because evaluating $J_i(x_i, x_{-i}, p)$2 requires solving an optimization problem. The key algorithmic device is a two-stage decomposition:

• Value-function learning: a surrogate model $J_i(x_i, x_{-i}, p)$3 is trained by standard regression on the best-response value data, with $J_i(x_i, x_{-i}, p)$4 and $J_i(x_i, x_{-i}, p)$5 regularization.
• GNE model learning: the surrogate $J_i(x_i, x_{-i}, p)$6 is substituted into the NI terms $J_i(x_i, x_{-i}, p)$7, yielding a single-level unconstrained problem in $J_i(x_i, x_{-i}, p)$8 combining the NI-based loss, a smooth log-sum-exponential penalty $J_i(x_i, x_{-i}, p)$9 approximating the worst constraint violation (with $g(x,p)\le 0$0 recovering the exact max), and regularization.

Three variants of the NI loss are proposed: the raw sum of $g(x,p)\le 0$1, the sum of positive parts $g(x,p)\le 0$2, and a smoothed version of the latter. The positive-part variants are important: they prevent the optimizer from exploiting imperfect value-function surrogates to make the predicted solution "super-optimal" (driving $g(x,p)\le 0$3 negative by underestimating $g(x,p)\le 0$4). When no GNE exists for a given $g(x,p)\le 0$5, the learned model still returns a decision minimizing the trade-off between aggregate equilibrium gap and constraint violation, tunable via $g(x,p)\le 0$6. A post-hoc projection step onto the shared constraint set is available when strict feasibility is required; for box constraints this reduces to clipping, or can be built into the model via an output saturation.

Justification of continuous models and a new continuity result

The use of continuous feedforward networks is justified in two steps. For the value-function surrogates, convexity of $g(x,p)\le 0$7, $g(x,p)\le 0$8, $g(x,p)\le 0$9 in $h(x,p)=0$0 yields a convex, continuous best-response value function; under merely partial convexity (each $h(x,p)=0$1 convex in $h(x,p)=0$2 alone), continuity of $h(x,p)=0$3 follows from known results on generalized derivatives of parametric value functions, provided uniqueness of the best response, constraint qualifications, and strong Slater's condition hold.

For the solution model, the paper establishes a new sufficient condition for the existence of a continuous parametric selection of variational GNEs (v-GNEs) of the associated variational inequality $h(x,p)=0$4. The main assumptions are: (i) strong variational stability (SVS) of the pseudo-gradient $h(x,p)=0$5 with respect to $h(x,p)=0$6, uniformly in $h(x,p)=0$7—a condition strictly weaker than strong monotonicity and satisfied by some non-monotone games; (ii) a Lipschitz-type bound on the distance of $h(x,p)=0$8 from $h(x,p)=0$9 as $p\in\mathcal{P}$0, which is shown to hold under Slater's condition and Lipschitz-continuous parametric constraint functions; and (iii) continuity of the set-valued projection $p\in\mathcal{P}$1. Under these assumptions, a continuity argument based on the SVS inequality and limit-point analysis shows that any limit point of $p\in\mathcal{P}$2 as $p\in\mathcal{P}$3 equals $p\in\mathcal{P}$4, hence the v-GNE selection is continuous. Notably, SVS also implies uniqueness of the variational solution for each $p\in\mathcal{P}$5. Outside these assumptions the learning procedure can still be run, but the learned map should then be interpreted as a best approximate equilibrium over an "almost" feasible set rather than an approximation to an exact solution.

The single-agent special case

For $p\in\mathcal{P}$6 the framework reduces to learning approximate solutions of standard multiparametric programs. The NI gap collapses to the suboptimality gap $p\in\mathcal{P}$7, and with the linear NI loss the value-function surrogate becomes unnecessary: the training problem reduces to jointly minimizing $p\in\mathcal{P}$8 over all samples under the feasibility penalty and regularization. Remarkably, dataset preparation then has zero optimization cost—only parameter samples $p\in\mathcal{P}$9 are required. This connects the approach to approximate explicit MPC, where exact multiparametric QP solutions require enumerating critical regions.

Numerical results

All experiments use the accompanying mpfit Python library with JAX, Adam followed by L-BFGS from 32 restarts, and report mean-squared best-response error (MSE$x^\star(p)$0) and constraint violations on held-out test sets.

Linear-quadratic GNEPs. On a 2-agent strongly monotone LQ-GNEP with 5 shared constraints, the smooth-max NI loss with $x^\star(p)$1 achieves average constraint violation $x^\star(p)$2 and MSE$x^\star(p)$3 of $x^\star(p)$4 in about 37 s of training; the raw-sum loss with small $x^\star(p)$5 exhibits exactly the predicted super-optimality pathology, with violations up to $x^\star(p)$6. Inference takes approximately $x^\star(p)$7s per parameter value—several orders of magnitude faster than solving the GNEP online. Scaling to random LQ-GNEPs with $x^\star(p)$8 agents, up to 96 shared constraints, and up to 4 parameters yields average violations in the range $x^\star(p)$9–$\hat{x}(p)\approx x^\star(p)$0 and MSE$\hat{x}(p)\approx x^\star(p)$1 up to $\hat{x}(p)\approx x^\star(p)$2, with accuracy degrading as the number of parameters grows.

Non-monotone LQ-GNEP. A 2-agent game with non-monotone pseudo-gradient (symmetric part eigenvalues $\hat{x}(p)\approx x^\star(p)$3 and $\hat{x}(p)\approx x^\star(p)$4) satisfying the SVS-based continuity conditions is learned with average best-response error $\hat{x}(p)\approx x^\star(p)$5 and zero constraint violation, illustrating that the method is not confined to monotone games.

QCQP and nonlinear GNEPs. A 3-agent quadratically constrained GNEP achieves mean violation $\hat{x}(p)\approx x^\star(p)$6 and best-response error $\hat{x}(p)\approx x^\star(p)$7. A nonlinear internet-switching game with an analytical best response is learned for $\hat{x}(p)\approx x^\star(p)$8 to $\hat{x}(p)\approx x^\star(p)$9 agents with zero constraint violation in all cases and squared best-response errors from $(x_k, p_k, \bar{J}_k)$0 ($(x_k, p_k, \bar{J}_k)$1) to $(x_k, p_k, \bar{J}_k)$2 ($(x_k, p_k, \bar{J}_k)$3), with training times growing to about 348 s. A second nonconvex nonlinear GNEP with tanh-coupled objectives achieves an average best-response error of $(x_k, p_k, \bar{J}_k)$4.

Multiparametric programming. For a 10-variable, 6-parameter mpQP with 50 constraints, the learned model attains an average relative suboptimality of $(x_k, p_k, \bar{J}_k)$5 versus the exact solution (which has 1627 critical regions), with average constraint violation $(x_k, p_k, \bar{J}_k)$6; inference takes $(x_k, p_k, \bar{J}_k)$7s versus $(x_k, p_k, \bar{J}_k)$8 ms for the exact multiparametric solver. A multiparametric QCQP with 20 additional convex quadratic constraints—admitting no exact explicit representation—is approximated with relative error $(x_k, p_k, \bar{J}_k)$9. These results indicate that the learned model trades a small suboptimality and occasional constraint violation for roughly a three-order-of-magnitude inference speedup.

Limitations and open questions

The paper concedes several points. The theoretical guarantees for the solution model require the SVS, projection-continuity, and constraint-regularity assumptions; outside these, the learned map is only a heuristic aggregate-gap minimizer, and the authors do not provide quantitative approximation-error bounds relating training loss to equilibrium accuracy. Constraint satisfaction is enforced only softly during training; worst-case violations remain non-negligible in the harder instances (up to $\bar{J}_{k,i}$0 for the mpQCQP and $\bar{J}_{k,i}$1 in the multi-agent LQ case), and the projection remedy is itself an optimization problem, albeit a cheaper one. The experiments also show systematic degradation with the number of parameters, and the approach assumes a full-information setting in which all costs and constraints are known to a central trainer. Open questions raised by the paper include integrating the projection step into training via differentiable optimization layers, using the learned models to warm-start GNE solvers, and exploiting sensitivity properties of variational inequalities to target specific v-GNE selections.

Conclusion

The paper offers a practical, single-level learning formulation for multiparametric GNEPs built on the NI gap function with learned value-function surrogates, requiring only best-response data rather than equilibrium solutions, and extending naturally to standard multiparametric programming. It is accompanied by a new sufficient condition for continuous parametric v-GNE selections under strong variational stability, generalizing strong monotonicity. The numerical evidence—microsecond-scale inference with small best-response errors and controllable constraint violations across LQ, QCQP, and nonlinear games—supports the approach as a viable alternative to exact multiparametric solution maps and online equilibrium solvers, with accuracy–feasibility trade-offs governed by the penalty weight $\bar{J}_{k,i}$2.