Adv Syst Sci Appl 2024; 02:19–31

Published online at https://ijassa.ipu.ru.

Piecewise Levenberg–Marquardt Method
for Generalized Nash Equilibrium Problems

Alexey Izmailov1*, Evgeniy Uskov2, Yan Zhibai1

1Lomonosov Moscow State University, Moscow, Russia
2 Derzhavin Tambov State University, Tambov, Russia

Abstract: We consider the constrained piecewise Levenberg–Marquardt method globalized by
linesearch, and apply it to “min” reformulations of the optimality systems for generalized Nash
equilibrium problems. Numerical comparison of the performance of this method with some
relevant existing alternatives is provided.

Keywords: generalized Nash equilibrium problem, piecewise smooth equation, constrained
equation, Levenberg–Marquardt method, LP-Newton method, globalization of
convergence

1. INTRODUCTION

This paper aims at application of the globalized constrained piecewise Levenberg–Marquardt
method developed in [24] to “min” reformulations of the optimality systems for generalized
Nash equilibrium problems (GNEP), and comparison with the existing alternatives.

In Section 2, we recall the statement of a GNEP and the related first-order optimality
system, and we discuss some possible equivalent reformulations of the latter as constrained
systems of equations. Section 3 discusses the Levenberg–Marquardt method for piecewise
smooth equations, its linesearch-based globalization, and the related global convergence and
superlinear rate of convergence theory. Finally, in Section 4, we apply this method to the
piecewise smooth reformulation of the first-order optimality system for GNEP, and compare
this algorithm with the usual Levenberg–Marquardt method for a smooth constrained
reformulation, and with the corresponding two versions of the LP-Newton method [7]
globalized according to [13].

Some words about our notation and blanket arrangements are in order. Let all the norms
be Euclidian, unless something different is explicitly specified. For a given y ∈ Rm and an
index set I ⊂ {1, . . . , m}, we denote by yI the subvector of y with the components yi, i ∈ I .
For a sequence {uk} ⊂ Rp convergent to some u∗ ∈ Rp, the rate of convergence is referred to
as superlinear with Q-order ν > 1 if there exists c > 0 such that ∥uk+1 − u∗∥ ≤ c∥uk − u∗∥ν
for all k large enough.

∗Corresponding author: izmaf@cs.msu.ru



20 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

2. GENERALIZED NASH EQUILIBRIUM PROBLEM
AND RELATED CONSTRAINED EQUATIONS

In a GNEP, N players are involved, and each player indexed by ν ∈ {1, . . . , N} controls the
variable xν ∈ Rnν , and aims at minimizing a smooth objective function fν : Rn → R subject
to constraints given by a smooth mapping gν : Rn → Rmν , where n =

∑N
ν=1 nν . In particular,

both fν and gν may depend not only on xν , but also on the variables of rival players, denoted
by x−ν ∈ Rn−nν , thus forming the entire vector of variables x = (xν , x−ν) ∈ Rn. With this
notation, the optimization problem of the ν-th player is written as follows:

minimize xν fν(x
ν , x−ν) subject to gν(xν , x−ν) ≤ 0. (2.1)

The specified problem setting goes back to [30]; recent surveys can be found in [10, 15],
including various applications. Among other relevant references on theory and numerical
methods for GNEP are [5, 8, 9, 11, 19–22, 27–29]. GNEPs form a difficult problem class, in
particular, because their solutions are naturally nonisolated.

Considering x−ν as a parameter, we introduce the Lagrangian Lν : Rnν × Rn−nν ×
Rmν → R of (2.1) as

Lν(x
ν , x−ν , λν) = fν(x

ν , x−ν) + ⟨λν , gν(xν , x−ν)⟩.

Then the Karush–Kuhn–Tucker (KKT) system characterizing stationary points and Lagrange
multipliers of the ν-th player’s optimization problem (2.1) has the form:

∂Lν

∂xν
(xν , x−ν , λν) = 0, λν ≥ 0, gν(xν , x−ν) ≤ 0, ⟨λν , gν(xν , x−ν)⟩ = 0.

Concatenating these systems over all players yields the KKT-type system of the GNEP:

L(x, λ) = 0, λ ≥ 0, g(x) ≤ 0, ⟨λν , gν(x)⟩ = 0, ν = 1, . . . , N, (2.2)

where we define λ = (λ1, . . . , λN) ∈ Rm, g : Rn → Rm, g(x) = (g1(x), . . . , gN(x)), with
m =

∑N
ν=1mν , and

L(x, λ) =
(
∂L1

∂x1
(x1, x−1, λ1), . . . ,

∂LN

∂xN
(xN , x−N , λN)

)
.

Furthermore, the system (2.2) can be equivalently reformulated as a constrained equation

Φ(u) = 0, u ∈ P, (2.3)

with some mapping Φ : Rp → Rq and a nonempty closed convex set P ⊂ Rp. This can be
done in many different ways, and we restrict ourselves to the following two adopted, for
example, in [13]. In both, p = q = n+ 2m, u = (x, λ, y), with y = (y1, . . . , yN) ∈ Rm

being a slack variable, and the constraint set is

P = Rn × Rm
+ × Rm

+ . (2.4)

The difference is in how Φ is defined. One possibility is to adopt a piecewise smooth “min”
reformulation with

Φ(u) = (L(x, λ), g(x) + y, min{λ, y}) , (2.5)
where min is taken componentwise, while another reformulation is smooth, with

Φ(u) = (L(x, λ), g(x) + y, λ ◦ y) , (2.6)

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 21

employing the Hadamard product λ ◦ y = (λ1y1, . . . , λmym). Observe that the set P in (2.4)
is polyhedral, and moreover, given by simple bounds, and this is the main reason to use slack
variables in these reformulations. Note also that for Φ defined in (2.5), one could actually take
P = Rn × Rm × Rm, and (2.3) would still be equivalent to (2.2). However, the nonnegativity
constraints on λ and y are needed to ensure the so-called P -property and condition (3.13) that
will be essential for the results presented in Section 3 below.

Reformulation employing (2.5) will be of primary interest for us, and to that end, we
next recall some terminology related to piecewise smoothness. A mapping Φ is refereed to as
piecewise smooth if it is continuous, and there exists a finite collection of smooth selection
mappings Φ1, . . . , Φs : Rp → Rq such that

Φ(u) ∈ {Φ1(u), . . . , Φs(u)} ∀u ∈ Rp.

Considering s = 1 recovers the case of a smooth mapping Φ. Needless to say, a collection
of smooth selection mappings corresponding to a given piecewise smooth mapping is not
uniquely defined, and we will assume that some such collection is fixed, i.e., Φ is defined by
a given collection of smooth selection mappings. Obviously, the mapping Φ defined according
to (2.5) is piecewise smooth, with a natural collection of smooth selection mappings

Φj(u) =
(
L(x, λ), g(x) + y, λI(j), y{1, ...,m}\I(j)

)
, (2.7)

where a one-to-one mapping j 7→ I(j) from {1, . . . , 2m} to the set of all different subsets of
{1, . . . , m} is supposed to be fixed.

3. PIECEWISE LEVENBERG–MARQUARDT METHOD

Getting back to a general piecewise smooth mapping Φ, for each u ∈ Rp we define the set

A(u) = {j ∈ {1, . . . , s} | Φ(u) = Φj(u)} (3.8)

of indices of smooth selection mappings active at u. Let G : Rp → Rq×p be any mapping
satisfying

G(u) ∈ {(Φj)′(u) | j ∈ A(u)} ∀u ∈ Rp. (3.9)

For a current iterate uk ∈ P , the (constrained) piecewise Levenberg–Marquardt (LM)
method generates the next iterate as uk + vk, where vk is a solution of the problem

minimize
1

2
∥Φ(uk) +G(uk)v∥2 + 1

2
σ(uk)∥v∥2 subject to uk + v ∈ P. (3.10)

where σ : P → R+ defines the values of the regularization parameter. If σ(uk) > 0, and if
P is polyhedral, then (3.10) is a quadratic programming problem with a strongly convex
objective function, and in particular this subproblem is uniquely solvable.

It follows from (3.8) and (3.9) that the subproblem (3.10) can be written in the form

minimize
1

2
∥Φj(uk) + (Φj)′(uk)v∥2 + 1

2
σ(uk)∥v∥2 subject to uk + v ∈ P,

with some j ∈ A(uk) (that may be different for different k, of course). This means that the
iteration of this method coincides with the iteration of the constrained LM method applied to
the smooth constrained equation

Φj(u) = 0, u ∈ P. (3.11)

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



22 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

In the smooth case, the constrained LM method was proposed in [26], and sharp results
on its local superlinear convergence were obtained in [1]. The piecewise smooth case was
studied in [6, 14], and an overview of these developments can be found in [17]. See also the
very recent improvements of the local convergence analysis of the piecewise LM method
in [25].

We proceed with an exposition of the globalized piecewise LM method and the
related convergence theory developed in [24]. The algorithm below can be regarded as a
generalization of the proposal in [18], where the smooth case was considered, and θ = 2 was
taken in the algorithm.
Algorithm 3.1:
Choose the parameters θ > 0, ε ∈ (0, 1) and κ ∈ (0, 1). Choose u0 ∈ P , and set k = 0.

1. If Φ(uk) = 0, stop.
2. Set σ(uk) = ∥Φ(uk)∥θ, and compute vk as the solution of (3.10). If vk = 0, stop.
3. Set α = 1. If for j ∈ A(uk) such that G(uk) = (Φj)′(uk) (see (3.9)) the inequality

∥Φ(uk + αvk)∥2 ≤ ∥Φ(uk)∥2 − εσ(uk)α∥vk∥2 (3.12)

holds, set αk = α. Otherwise replace α by κα and check again the inequality in (3.12),
etc., until (3.12) is satisfied, and then set αk = α.

4. Set uk+1 = uk + αkv
k, increase k by 1 and go to Step 1.

The key ingredient of the global convergence analysis is the following assumption [24,
(1.5)] that appeared before in [13, (4.8)], [16, (32)], and in [3, (0.6)]:

∥Φ(u)∥ ≤ ∥Φj(u)∥ ∀ j ∈ {1, . . . , s}, ∀u ∈ P. (3.13)

Theorem 3.1:
Let Φ : Rp → Rq be a piecewise smooth mapping with continuously differentiable smooth
selection mappings Φ1, . . . , Φs : Rp → Rq. Let P ⊂ Rp be a nonempty closed convex set,
and assume that (3.13) holds. Let G : Rp → Rq×p be a fixed mapping satisfying (3.9).

Then Algorithm 3.1 uniquely defines the iterates u0, u1, . . ., and either terminates at a
point uk ∈ P satisfying

⟨((Φj)′(uk))⊤Φ(uk), u− uk⟩ ≥ 0 ∀u ∈ P, (3.14)

for at least one j ∈ A(uk), or generates an infinite sequence {uk}, and any accumulation
point ū of this sequence belongs to P and satisfies

⟨((Φj)′(ū))⊤Φ(ū), u− ū⟩ ≥ 0 ∀u ∈ P, (3.15)

for at least one j ∈ A(ū).
We complete this section with a result on asymptotic superlinear convergence rate, and

in order to do this, we need the following two assumptions. The P -property at ū, introduced
in [14, p. 434], consists of saying that Uj ⊂ U near ū for all j ∈ A(ū), where U stands for the
solution set of (2.3), while Uj stands for the solution set of (3.11). Evidently, (3.13) implies
the P -property at any solution.

Another assumption we need is the constrained local Lipschitzian error bound for each
active selection, that is,

∀ j ∈ A(ū) dist(u, Uj) = O(∥Φj(u)∥) as u ∈ P tends to ū. (3.16)

Theorem 3.2:
Let Φ : Rp → Rq be a given mapping, P ⊂ Rp a closed convex set, and ū ∈ U . Assume that Φ

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 23

is piecewise smooth and the derivatives of its smooth selection mappings Φ1, . . . , Φs : Rp →
Rq are Lipschitz-continuous near ū. Let the P -property at ū and condition (3.16) as u ∈ P
tends to ū be satisfied. Let G : Rp → Rq×p be a fixed mapping satisfying (3.9).

Then, if Algorithm 3.1 run with θ ∈ (0, 2] generates an iterate close enough to ū, it either
terminates with uk ∈ U , or generated an infinite sequence {uk} convergent to some u∗ ∈ U ,
and the rate of convergence is superlinear with Q-order min{θ + 1, 2}.

One can easily see that the assumption (3.13), and hence, the P -property at any solution
always hold for the GNEP KKT-type system reformulation using (2.4) and (2.5). As for the
assumption (3.16), some natural sufficient conditions for it were derived in [14, Theorems 4,
5], and we adapt them to our problem setting next.

Proposition 3.1:
Let fν : Rn → R and gν : Rn → Rmν , ν ∈ {1, . . . , N}, be twice differentiable near x̄ ∈ Rn,
with their second derivatives being Lipschitz-continuous near x̄. Assume that ū = (x̄, λ̄) with
some λ̄ ∈ Rm is a solution of (2.2), where the notation introduces in Section 2 for GNEP is
employed. Let Φ be defined in (2.5), with its smooth selection mappings Φj , j ∈ {1, . . . 2m},
defined according to (2.7), and let P be defined in (2.4). Furthermore, define the index set

A = {i ∈ {1, . . . , m} | gi(x̄) = 0},

and its partitions

A+ = {i ∈ A | ∃ ν ∈ {1, . . . , N} : λ̄ν
i > 0}, A0 = A \ A+,

and

Aν
+ = {i ∈ A | λ̄ν

i > 0}, Aν
0 = A \ Aν

+,

for every ν ∈ {1, . . . , N}.
If for any I ⊂ A0 and any Iν ⊂ Aν

0 , ν ∈ {1, . . . , N}, the matrices

∂L
∂x

(x, λ)

(
∂gA1

+∪I1

∂x1
(x)

)⊤

0 . . . 0

0

(
∂gA2

+∪I2

∂x2
(x)

)⊤

. . . 0

...
... . . . ...

0 . . . . . .

(
∂gAN

+∪IN

∂xN
(x)

)⊤

g′A+∪I(x) 0


have the same rank for all (x, λ) ∈ Rn × Rm near (x̄, λ̄), then the piecewise constrained
error bound condition (3.16) holds.

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



24 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

Corollary 3.1:
Under the assumptions of Proposition 3.1, if the matrix

∂L
∂x

(x̄, λ̄)

(
∂gA1

+

∂x1
(x̄)

)⊤

0 . . . 0

0

(
∂gA2

+

∂x2
(x̄)

)⊤

. . . 0

...
... . . . ...

0 . . . . . .

(
∂gAN

+

∂xN
(x̄)

)⊤

g′A(x̄) 0


has full row rank, then the piecewise constrained error bound condition (3.16) holds.

The sufficient condition for the piecewise constrained error bound in Proposition 3.1 is
automatically satisfied if fν is quadratic and gν is affine, for all ν ∈ {1, . . . , N}. Observe
also that any strict complementarity-like conditions are neither involved in Proposition 3.1,
nor in its Corollary 3.1.

4. NUMERICAL RESULTS

The test set used for the numerical comparisons presented below is a collection of GNEPs
from [5], later also used in [4, 13]. Information about these test problems and references to
their sources and detailed descriptions are provided in [5]. For each test problem, 20 random
starting points were used (the same for all algorithms involved, with random components of
x0 distributed uniformly within (0.1, 20), which guarantees that all the functions appearing
in all test problems are well-defined at x0. Furthermore, λ0 and y0 were chosen according
to [13, Section 5].

Throughout the rest of this section, P is defined according to (2.4). We compared the
performance of Algorithm 3.1 applied to (2.3) with Φ defined in (2.5) (abbreviated below as
PWLM) with the following alternatives:

• Globalized LP-Newton method from [13, Algorithm 1] supplied with all the perfor-
mance-improving modification proposed in [13, Section 5], and applied to the same
instances of (2.3) as PWLM (abbreviated as PWLPN).

• The LM method globalized according to Algorithm 3.1, but applied to (2.3) with smooth
Φ defined in (2.6) (abbreviated as LM-Had).

• Globalized LP-Newton method, the same as PWLPN, but applied to (2.3) with smooth
Φ defined in (2.6) (abbreviated as LPN-Had).

In case of a smooth Φ from (2.6), (3.9) implies that G(u) at any u = (x, λ, y) ∈ Rp equals
Φ′(u). For Φ from (2.5), the values of a mapping G satisfying (3.9) were computed according
to the following rule: for every i ∈ {1, . . . , m}, if λi > yi, then the (n+m+ i)-th row of
G(u) is Gn+m+i(u) = (0, ei), with ei being the i-th element of the standard basis in Rm;
otherwise Gn+m+i(u) = (ei, 0); the other rows of G(u) are the gradients of the corresponding
(smooth) components of Φ.

The experiments were performed in Matlab, with its built-in solver quadprog used
for quadratic programming subproblems of PWLM and LM-Had, and linprog for

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 25

100 101

0

0.2

0.4

0.6

0.8

1

PWLM 1e-5
PWLM 1e-8
PWLM 1e-10

(a) For “min” reformulation

100 101

0

0.2

0.4

0.6

0.8

1

LM-Had 1e-5
LM-Had 1e-8
LM-Had 1e-10

(b) For Hadamard product-based reformulation

Fig. 4.1. Comparison by average iteration counts (the same as evaluations of G) for different values of σ̄.

100 101

0

0.2

0.4

0.6

0.8

1

PWLM 1e-10
PWLPN
LM-Had 1e-10
LPN-Had

(a) By average iteration counts

100 101

0

0.2

0.4

0.6

0.8

1

PWLM 1e-10
PWLPN
LM-Had 1e-10
LPN-Had

(b) By average evaluations of Φ

Fig. 4.2. Comparison of different algorithms.

linear programming subproblems of PWLPN and LPN-Had, the former run with option
’OptimalityTolerance’, 1e-15, and the latter with options ’OptimalityTolerance’, 1e-8,
’ConstraintTolerance’, 1e-8.

The parameters in Algorithm 3.1: were chosen as follows: θ = 2, ε = 0.001, κ = 0.5. The
parameters of the globalized LP-Newton methods were chosen according to [13, Section 5].
In order to avoid large values of the regularization parameter far from solutions, the rule for
it in Step 2 of Algorithm 3.1 was replaced by σ(uk) = min{σ̄, ∥Φ(uk)∥θ}, where several
different values of σ̄ > 0 were tied. Runs were declared successful when terminated because
of

∥Φ(uk)∥∞ ≤ 10−6,

within 1000 iterations, where for all the algorithms involved, Φ defined in (2.5) was used
in this stopping test. Otherwise, failure was declared, as well as in the cases when the
backtracking procedure in Step 3 of Algorithm 3.1 was producing the trial stepsize parameter

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



26 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

100 101

0

0.2

0.4

0.6

0.8

1

PWLM 1e-10
PWLPN
LM-Had 1e-10
LPN-Had

(a) All test problems

100 101

0

0.2

0.4

0.6

0.8

1

PWLM 1e-10
PWLPN
LM-Had 1e-10
LPN-Had

(b) Large test problems

Fig. 4.3. Comparison by time.

value α ≤ 10−13, and also when ∥((Φj)′(uk))⊤Φ(uk)∥ ≤ 10−12 was encountered for j ∈
A(uk) such that G(uk) = (Φj)′(uk), and in particular, the stationarity condition (3.14) was
approximately satisfied. For the LP-Newton-based algorithms, similar rules were adopted, but
the approximate stationarity test was replaced by |∆(uk)| ≤ 10−12, where the quantity ∆(uk)
estimating from above the directional derivative of the merit function is defined according
to [13, (2.3)], where γ(uk) is taken the optimal value of the LP-Newton subproblem.

It turned out that in order to achieve reasonable robustness of the LM-based algorithms, σ̄
has to be taken quite small: larger values lead to long series of short steps far from solutions,
and this often ends up with a failure. This happens for LM-Had when run for problems
A7, Heu, Tr1b, Tr1c, and for PWLM when run for A10e, Tr1a. For A10c and A10d, runs
of PWLM are generally successful for larger σ̄ as well, but taking smaller values of this
parameter reduces the iteration count significantly.

Our numerical results are presented in the form of performance profiles originally
proposed in [2], and later adapted in [23] for the case when multiple starting points are used
for every test problem. Specifically, for each algorithm a, we present the plot the function
πa : [1, +∞) → [0, 1] constructed as follows. Let ka

τ be the average of some measure of
efficiency of algorithm a on problem τ , where the average is taken over successful runs, and
let saτ ∈ [0, 1] stand for the portion of successful runs of algorithm a on problem τ . Let rτ be
the best (say, smallest) value of ka

τ among all the algorithms. Then for each t ∈ [1, +∞), we
set

πa(t) =
1

T

∑
τ∈Ra(t)

saτ ,

where T the overall number of problems in the test set, and Ra(t) is the subset of problems
for which the performance of algorithm a is no more than t times worse than that of the best
algorithm:

Ra(t) = {τ ∈ {1, . . . , T} | ka
τ ≤ trτ}.

In particular, πa(t) for large t is the average portion of successful runs (averaged over all
test problems), while πa(1) is the quantity computed in a similar way but using the portions
of successful runs only on those problems for which the average performance of the given
algorithm over successful runs is the best among all algorithms being tested. In case of a
single run for each test problem, this agrees with the original proposal in [2], and in particular,

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 27

Table 4.1. Average elapsed times (seconds) and percentage of tails of last full steps out of the average iteration
counts, over successful runs

Test PWLM PWLPN LM-Had LPN-Had
Time Tail Time Tail Time Tail Time Tail

A1 0.080 100.0 0.068 100.0 0.016 76.3 0.057 100.0

A2 NaN NaN 0.044 100.0 NaN NaN 0.067 100.0

A3 0.003 100.0 0.043 100.0 0.006 100.0 0.066 100.0

A4 0.010 82.8 0.064 98.9 0.011 78.4 0.055 97.1

A5 0.005 100.0 0.042 100.0 0.009 100.0 0.053 100.0

A6 0.014 94.6 0.070 99.7 0.013 90.6 0.067 99.3

A7 0.022 100.0 0.092 100.0 0.018 86.5 0.077 100.0

A8 0.003 100.0 0.025 100.0 0.006 100.0 0.033 100.0

A9a 0.064 75.7 2.508 18.0 0.055 97.7 0.277 70.8

A9b 1.349 35.1 0.612 59.9 0.173 100.0 0.737 77.3

A10a 0.031 49.0 0.101 99.2 0.018 82.3 0.091 98.1

A10b 0.248 75.6 0.466 99.7 0.309 87.3 0.520 95.8

A10c 8.627 10.2 1.612 84.3 NaN NaN 4.097 81.0

A10d 6.777 36.3 9.110 55.1 2.800 92.3 11.134 60.7

A10e 25.550 37.2 42.698 42.6 NaN NaN 74.579 34.7

A11 0.004 100.0 0.021 100.0 NaN NaN 0.029 100.0

A12 0.003 100.0 0.022 100.0 0.003 100.0 0.025 100.0

A13 0.004 100.0 0.027 100.0 0.004 100.0 0.026 100.0

A14 0.010 100.0 0.040 98.9 0.012 81.8 0.049 100.0

A15 0.003 100.0 0.034 100.0 0.004 100.0 0.037 100.0

A16a 0.008 100.0 0.030 100.0 0.004 100.0 0.033 97.9

A16b 0.005 100.0 0.030 100.0 0.004 100.0 0.030 100.0

A16c 0.006 100.0 0.031 100.0 0.005 100.0 0.039 84.0

A16d 0.007 100.0 0.032 100.0 0.007 100.0 0.037 100.0

A17 0.004 100.0 0.024 100.0 0.005 100.0 0.036 100.0

A18 0.007 100.0 0.050 100.0 0.011 94.5 0.053 100.0

Harker 0.003 100.0 0.023 100.0 0.006 100.0 0.036 100.0

Heu 0.012 100.0 0.053 100.0 0.011 100.0 0.057 100.0

Lob 0.068 37.6 NaN NaN NaN NaN NaN NaN

NTF1 0.004 100.0 0.022 100.0 0.003 100.0 0.031 100.0

NTF2 0.007 90.5 0.036 99.4 0.006 94.8 0.034 100.0

Spam 1380.005 72.7 545.909 95.9 NaN NaN 166.393 83.7

Tr1a 0.059 71.1 0.111 90.4 0.050 94.4 0.071 96.7

Tr1b NaN NaN 1.617 66.8 1.188 43.2 1.636 46.1

Tr1c 3.878 26.2 2.245 71.4 2.231 39.7 1.365 62.6

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



28 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

Table 4.2. Failures

Test PWLM PWLPN LM-Had LPN-Had

A1 0 0 0 1 SS

A2 20 SS 0 20 (13 SF, 2 IL, 5 SS) 0

A3 0 12 SP 1 SF 9 (7 SP, 2 SS)

A4 4 (1 SF, 3 SS) 0 4 SF 0

A5 0 0 0 0

A6 2 SS 0 1 SF 0

A7 14 SS 8 SP 0 0

A8 0 0 5 SF 0

A9a 8 (2 SF, 6 SS) 12 IL 0 0

A9b 1 IL 0 0 0

A10a 0 0 3 SF 0

A10b 0 0 2 (1 SF, 1 SS) 0

A10c 0 0 20 (2 SF, 18 SS) 0

A10d 0 5 (4 IL, 1 SP) 19 (8 SF, 11 SS) 2 IL

A10e 2 (1 SF, 1 SS) 9 IL 20 (13 SF, 7 SS) 0

A11 0 0 20 SF 0

A12 0 0 5 SF 0

A13 5 SF 0 0 0

A14 0 0 0 0

A15 0 0 2 SF 0

A16a 1 SF 0 2 SF 0

A16b 2 SF 0 2 SF 0

A16c 1 SF 0 6 SF 0

A16d 0 0 2 SF 0

A17 0 0 2 SF 0

A18 0 0 1 SF 0

Harker 0 0 0 0

Heu 15 SS 0 6 SF 0

Lob 0 20 IL 20 SS 20 IL

NTF1 0 0 4 SF 0

NTF2 0 0 4 SF 0

Spam 0 15 SF 20 SF 17 (14 SF, 1 SP, 2 SS)

Tr1a 16 (9 SF, 7 SS) 2 SP 1 SF 0

Tr1b 20 (14 SF, 1 IL, 5 SS) 1 SP 12 (9 SF, 3 SS) 0

Tr1c 18 (5 SF, 3 IL, 10 SS) 4 (3 IL, 1 SP) 11 (5 SF, 1 IL, 5 SS) 0

πa(t) for large t is the portion of test problems the runs of algorithm a on which are successful,
while πa(1) is the portion of problems on which the runs of algorithm a are successful and
its efficiency is the best.

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 29

Figure 4.1 presents the comparison for three different values of σ̄. Figure 4.1a
demonstrates the advantages of running PWLM with σ̄ = 10−10, by iteration count, while
robustness of all compared variants is nearly the same. In Figure 4.1b, the comparison by
iteration count is similar, while robustness for this choice of σ̄ is only slightly lower than for
the value σ̄ = 10−8. That is why σ̄ = 10−10 was adopted in further experiments.

At this point we mention that we have also experimented with different values of the
exponent θ in the regularization parameter employed by PWLM. Specifically, we tried θ =
1.5 and θ = 1, along with the basic choice θ = 2, but much difference in performance was
detected neither by iteration counts not by evaluations of Φ. The variant with θ = 1 seriously
outperforms the alternatives by the elapsed time, and especially the variant with θ = 2. This
phenomenon is explained by the difference in solution times for PWLM subproblems formed
with different θ, which shows up mostly for relatively small problems. For large problems (as
selected below), the difference in behavior is negligible, by all the criteria, including times,
and we kept θ = 2 for the rest of the experiments.

In Figure 4.2a, we compare the specified versions of PWLM and LM-Had with each
other, and with PWLPN and LPN-Had, with average counts of iterations and evaluations of
Φ adopted as efficiency measures. PWLM demonstrates robustness somehow lower than that
of LP-Newton-based methods (thus providing an extra evidence of the use of the latter), but
definitely higher than that of LM-Had. Moreover, the performance of PWLPN and LPN-Had
is quite similar, while PWLM by far outperforms all the alternatives by efficiency.

Furthermore, in Figure 4.3, we provide the results of comparison by average elapsed
time as a measure of efficiency. Figure 4.3a is for the entire set of test problems, while
in Figure 4.3b, only selected large problems were left for comparison, namely, those with
more that 100 variables and/or constraints: A9b, A10b, A10c, A10d, A10e, Spam, Tr1b,
Tr1c. PWLM outperforms both PWLPN and LPN-Had by efficiency on the full test set,
although it is somewhat less robust. The picture with robustness remains the same for large
problems (and robustness of LM-Had appears very low in this case), but now the performance
of PWLPN and LPN-Had is much improved by efficiency as well, apparently because
large linear programming subproblems are solved faster than the quadratic programming
subproblems. The detailed information on average times is collected in Table 4.1, where
the identifiers of large problems are boldfaced. In addition, in that table, we also report on
the percentage of tails of last full (with αk = 1) steps of the algorithms in question before
successful termination, out of the average iteration counts. The issue of ultimate acceptance
of the full step is crucial for superlinear convergence guaranties, like those established for
Algorithm 3.1 in Theorem 3.2. In our experiments, the only cases where the last step was
not full were encountered only for LM-Had, in 5 runs out of 16 successful for NTF2, and for
LPN-Had, in 1 run out of 20 successful for A16c.

The largest problem in the test set is Spam, in which N = 101, n = 2020, m = 4040. It
is successfully solved by PWLM, and in a modest number of iterations, while all the other
algorithms involved typically fail on it. Many failures are caused by inability of quadprog or
linprog to solve the corresponding subproblems; for example, this is the case for PWLPN,
LM-Had, and LPN-Had, when run for Spam. LM-Had also systematically fails for the
specified reason for problems A2, A10d, A10e, A11, Tr1c. Typical failures of other kind
are caused by the stepzize parameter becoming too small.

Problem A2 is also troublesome for PWLM in a different way, with all runs ending
up with failures because of a too small stepsize generated, but that was actually caused
by convergence to stationary points for an active smooth selection being used, that are not
solutions. Observe that such scenario is of course not ruled out by Theorem 3.1.

Information about failures and their reasons is summarized in Table 4.2 making use of the
following abbreviations for failures’ reasons:

• SF: solving subproblem failed.
• IL: failure by iteration limit.
• SP: stationary point that is not a solution.

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)



30 A. F. IZMAILOV, E. I. USKOV, YAN ZHIBAI

• SS: stepsize too small.

At the end, we briefly mention some possibilities to further improve the performance
of PWLM for GNEP. Tuning σ̄ for each problem sometimes gives some positive effect.
Moreover, alternative rules for σ(·) might be tried, such as σ(u) = σ̄∥Φ(u)∥θ/(1 + ∥Φ(u)∥θ),
or σ(u) = σ̄∥(Φ′(u))⊤Φ(u)∥θ as proposed in [12], perhaps involving some extra parameters,
or some combinations of these rules. The cases when a short step is generated can be handled
in a more sophisticated way rather than just declaring failure, like trying some safeguarding
steps in such cases. Finally, the use of more advanced QP-solvers for subproblems should
certainly be helpful, e.g., those applying some pre-processing/scaling. All that will be the
subject of future research.

ACKNOWLEDGEMENTS

The authors thank Axel Dreves and Markus Herrich for their assistance with a Matlab
implementation of the GNEP test collection and of the LP-Newton method, used in Section 4
of this work. This work was funded by the Russian Science Foundation Grant 24-21-00015
(https://rscf.ru/en/project/24-21-00015/).

REFERENCES

1. Behling, R. & Fischer, A. (2012) A unified local convergence analysis of inexact
constrained Levenberg–Marquardt methods, Optim. Lett., 6, 927–940.

2. Dolan, E.D. & More, J.J. (2002) Benchmarking optimization software with
performance profiles, Math. Program., 91, 201–213.

3. Dorovskikh, D.I., Izmailov, A.F. & Uskov, E.I. (2024) Globalizing convergence of
piecewise Newton methods, Russian Universities Reports. Mathematics, 29, 149–163.

4. Dreves, A., Facchinei, F., Fischer, A. & Herrich, M. (2014) A new error bound result
for Generalized Nash Equilibrium Problems and its algorithmic application, Comput.
Optim. Appl., 59, 63–84.

5. Dreves, A., Facchinei, F., Kanzow, C. & Sagratella, S. (2011) On the solution of the
KKT conditions of generalized Nash equilibrium problems, SIAM J. Optim., 21, 1082–
1108.

6. Facchinei, F., Fischer, A. & Herrich, M. (2013) A family of Newton methods for
nonsmooth constrained systems with nonisolated solutions, Math. Methods Oper. Res.,
77, 433–443.

7. Facchinei, F., Fischer, A. & Herrich, M. (2014) An LP-Newton method: Nonsmooth
equations, KKT systems, and nonisolated solutions, Math. Program., 146, 1–36.

8. Facchinei, F., Fischer, A. & Piccialli, V. (2007) On generalized Nash games and
variational inequalities, Oper. Res. Lett., 35, 159–164.

9. Facchinei, F., Fischer, A. & Piccialli, V. (2009) Generalized Nash equilibrium
problems and Newton methods, Math. Program., 117, 163–194.

10. Facchinei, F. & Kanzow, C. (2010) Generalized Nash equilibrium problems, Ann. Oper.
Res., 175, 177–211.

11. Facchinei, F. & Kanzow, C. (2010) Penalty methods for the solution of generalized
Nash equilibrium problems, SIAM J. Optim., 20, 2228–2253.

12. Fischer, A. (2002) Local behavior of an iterative framework for generalized equations
with nonisolated solutions, Math. Program., 94, 91–124.

13. Fischer, A., Herrich, M., Izmailov, A.F. & Solodov, M.V. (2016) A globally convergent
LP-Newton method, SIAM J. Optim., 26, 2012–2033.

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)

https://rscf.ru/en/project/24-21-00015/


PIECEWISE LEVENBERG–MARQUARDT METHOD FOR GNEPS 31

14. Fischer, A., Herrich, M., Izmailov, A.F. & Solodov, M.V. (2016) Convergence
conditions for Newton-type methods applied to complementarity systems with
nonisolated solutions, Comput. Optim. Appl., 63, 425–459.

15. Fischer, A., Herrich, M. & Schönefeld, K. (2014) Generalized Nash equilibrium
problems – recent advances and challenges, Pesquisa Operacional, 34, 521–558.

16. Fischer, A., Izmailov, A.F. & Jelitte, M. (2021) Newton-type methods near critical
solutions of piecewise smooth nonlinear equations, Comput. Optim. Appl., 80, 587–
615.

17. Fischer, A., Izmailov, A.F. & Solodov, M.V. (2023) The Levenberg–Marquardt
method: an overview of modern convergence theories and more, Comput. Optim. Appl.,
2024.

18. Fischer, A. & Shukla, P.K. (2008) A Levenberg–Marquardt algorithm for
unconstrained multicriteria optimization, Oper. Res. Lett., 36, 643–646.

19. Fukushima, M. & Pang, J.-S. (2005) Quasi-variational inequalities, generalized Nash
equilibria, and multi-leader-follower games, Comput. Manag. Science, 2, 21–56.

20. Harker, P.T. (1991) Generalized Nash games and quasi-variational inequalities, Eur. J.
Oper. Res., 54, 81–94.

21. von Heusinger, A., Kanzow, C. & Fukushima, M. (2012) Newton’s method for
computing a normalized equilibrium in the generalized Nash game through fixed point
formulation, Math. Program., 132, 99–123.

22. Izmailov, A.F. & Solodov, M.V. (2014) On error bounds and Newton-type methods for
generalized Nash equilibrium problems, Comput. Optim. Appl., 59, 201–218.

23. Izmailov, A.F., Solodov, M.V. & Uskov, E.I. (2015) Combining stabilized SQP with
the augmented Lagrangian algorithm, Comput. Optim. Appl., 62, 405–429.

24. Izmailov, A.F., Uskov, E.I. & Yan Zhibai (2024) Globalization of convergence of the
constrained piecewise Levenberg–Marquardt method, Submitted.

25. Izmailov, A.F., Uskov, E.I. & Yan Zhibai (2024) The piecewise Levenberg–Marquardt
method, Adv. Syst. Sci. Appl., 24, 29–39.

26. Kanzow, C., Yamashita, N. & Fukushima, M. (2004) Levenberg–Marquardt methods
with strong local convergence properties for solving nonlinear equations with convex
constraints, J. Comput. Appl. Math., 172, 375–397.

27. Kubotam, K. & Fukushima, M. (2010) Gap function approach to the generalized Nash
equilibrium problem, J. Optim. Theory Appl., 144, 511–531.

28. Kulkarni, A.A. & Shanbhag, U.V. (2012) Revisiting generalized Nash games and
variational inequalities, J. Optim. Theory Appl., 154, 175–186.

29. Nabetani, K., Tseng, P. & Fukushima, M. (2011) Parametrized variational inequality
approaches to generalized Nash equilibrium problems with shared constraints,
Comput. Optim. Appl., 48, 423–452.

30. Rosen, J.B. (1965) Existence and uniqueness of equilibrium points for concave N -
person games, Econometrica, 33, 52–534.

Copyright © 2024 ASSA. Adv Syst Sci Appl (2024)


	Introduction
	Generalized Nash equilibrium problem  and related constrained equations
	Piecewise Levenberg–Marquardt method
	Numerical results