Optimization benchmark with the GNE package

As usual, the package is loaded via the function. In the following, we assume that the line below has been called

library(GNE)

Introduction

If we have parametrized action space Xi(x−i) = {yi, gi(yi, x−i) ≤ 0}, we denote the GNEP by GNEP(N, θi, gi).

We denote by X(x) the action set X(x) = X1(x−1) × … × XN(x−N). For standard NE, this set does not depend on x.

The following example seems very basic, but in fact it has particular features, one of them is to have four solutions, i.e. four GNEs. Let N = 2. The objective functions are defined as θ1(x) = (x1 − 2)2(x2 − 4)4textandθ2(x) = (x2 − 3)2(x1)4, for $x\in \R^2$, while the constraint functions are given by g1(x) = x1 + x2 − 1 ≤ 0 and g2(x) = 2x1 + x2 − 2 ≤ 0. Objective functions can be rewritten as θi(x) = (xi − ci)2(x−idi)4, with c = (2, 3) and d = (4, 0). First-order derivatives are ∇jθi(x) = 2(xi − ci)(x−idi)4δij + 4(xi − ci)2(x−idi)3(1 − δij), and ∇jg1(x) = 1 and ∇jg2(x) = 2δj1 + δj2. Second-order derivatives are and ∇k∇jg1(x) = ∇k∇jg2(x) = 0.

GNEP as a nonsmooth equation

Notation and definitions

From , assuming differentiability and a constraint qualification hold, the first-order necessary conditions of player i’s subproblem state there exists a Lagrangian multiplier $\lambda^i \in {\R}^{m_i}$ such that Regrouping the N subproblems, we get the following system.

Using complementarity function ϕ(a, b) (e.g. min (a, b)), we get the following nonsmooth equation $$ \Phi(z) = \left( \begin{matrix} \tilde L(x, \lambda) \\ \phi_.(-G(x), \lambda) \\ \end{matrix} \right) = 0 , $$ where ϕ. is the component-wise version of the function ϕ and L̃ is the Lagrangian function of the extended system. The generalized Jacobian is given in Appendix .

A classic example

Returning to our example, we define the Φ as $$ \Phi(x) = \left( \begin{matrix} 2(x_1-2) (x_2-4)^4 + \lambda_1 \\ 2(x_2-3) (x_1)^4 + \lambda_2 \\ \phi(\lambda_1, 1- x_1-x_2) \\ \phi(\lambda_2, 2- 2x_1-x_2) \\ \end{matrix} \right) , $$ where ϕ denotes a complementarity function. In , we use

myarg <- list(C=c(2, 3), D=c(4,0))
dimx <- c(1, 1)
#Gr_x_j O_i(x)
grobj <- function(x, i, j, arg)
{
    dij <- 1*(i == j)
    other <- ifelse(i == 1, 2, 1)
    res <- 2*(x[i] - arg$C[i])*(x[other] - arg$D[i])^4*dij 
    res + 4*(x[i] - arg$C[i])^2*(x[other] - arg$D[i])^3*(1-dij) 
}

dimlam <- c(1, 1)
#g_i(x)
g <- function(x, i)
    ifelse(i == 1, sum(x[1:2]) - 1, 2*x[1]+x[2]-2)
#Gr_x_j g_i(x)
grg <- function(x, i, j)
    ifelse(i == 1, 1, 1 + 1*(i == j))

Note that the triple dot arguments … is used to pass arguments to the complementarity function.

Elements of the generalized Jacobian of Φ have the following form $$ \partial \Phi(x) = \left\{ \left( \begin{matrix} 2(x_2-4)^4 & 8(x_1-2) (x_2-4)^3 & 1 & 0 \\ 8(x_2-3) (x_1)^3 & 2(x_1)^4 & 0 & 1 \\ -\phi_b'(\lambda_1, 1- x_1-x_2) & -\phi_b'(\lambda_1, 1- x_1-x_2) & \phi_a'(\lambda_1, 1- x_1-x_2) & 0 \\ -2\phi_b'(\lambda_2, 2- 2x_1-x_2) & - \phi_b'(\lambda_2, 2- 2x_1-x_2) & 0 & \phi_a'(\lambda_2, 2- 2x_1-x_2) \\ \end{matrix} \right) \right\} , $$ where ϕa′ and ϕb′ denote elements of the generalized gradient of the complementarity function. The corresponding code is

#Gr_x_k Gr_x_j O_i(x)
heobj <- function(x, i, j, k, arg)
{
    dij <- 1*(i == j)
    dik <- 1*(i == k)
    other <- ifelse(i == 1, 2, 1)
    res <- 2*(x[other] - arg$D[i])^4*dij*dik 
    res <- res + 8*(x[i] - arg$C[i])*(x[other] - arg$D[i])^3*dij*(1-dik)
    res <- res + 8*(x[i] - arg$C[i])*(x[other] - arg$D[i])^3*(1-dij)*dik
    res + 12*(x[i] - arg$C[i])^2*(x[other] - arg$D[i])^2*(1-dij)*(1-dik)
}
#Gr_x_k Gr_x_j g_i(x)
heg <- function(x, i, j, k) 0

Usage example

Therefore, to compute a generalized Nash equilibrium, we use

set.seed(1234)
z0 <- rexp(sum(dimx)+sum(dimlam))
GNE.nseq(z0, dimx, dimlam, grobj=grobj, myarg, heobj=heobj, myarg, 
    constr=g, grconstr=grg, heconstr=heg, 
    compl=phiFB, gcompla=GrAphiFB, gcomplb=GrBphiFB, method="Newton", 
    control=list(trace=0))
## GNE: 2 -1.999999 3.054699e-17 79.99999 
## with optimal norm 5.086687e-07 
## after  25 iterations with exit code 1 .
## Output message: Function criterion near zero 
## Function/grad/hessian calls: 28 25 
## Optimal (vector) value: 3.054699e-17 0 0 5.086687e-07

Recalling that the true GNEs are

##   x1 x2 lam1 lam2
## 1  2 -2    0  160
## 2 -2  3    8    0
## 3  0  1  324    0
## 4  1  0  512    6

Localization of the GNEs

On figure , we draw contour plots of the function $\frac{1}{2} || \Phi(z) ||^2$ with respect to x1 and x2, given λ1 and λ2. The second figure just plots the initial points and the 6 GNEs.

Benchmark of the complementarity functions and the computation methods

Using the following function, we compare all the different methods with different initial points and different complementarity functions. We consider the following complementarity functions.

Firstly, we define a function calling the benchmark function for the five complementarity functions under consideration.

wholebench <- function(z0)
{
  #min function
  resMin <- bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiMin), 
                           argjac=list(gphia= GrAphiMin, gphib= GrBphiMin), echo=FALSE)
  
  #FB function
  resFB <- bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiFB), 
                          argjac=list(gphia= GrAphiFB, gphib= GrBphiFB), echo=FALSE)
  
  #Mangasarian function
  resMan <- bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiMan, f=function(t) t^3), 
                        argjac=list(gphia= GrAphiMan, gphib= GrBphiMan, fprime=function(t) 3*t^2),
                        echo=FALSE, control=list(maxit=200))
  
  #LT function
  resLT <- bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiLT, q=4), 
                          argjac=list(gphia= GrAphiLT, gphib= GrBphiLT, q=4))
  
  #KK function
  resKK <- bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiKK, lambda=3/2), 
                          argjac=list(gphia= GrAphiKK, gphib= GrBphiKK, lambda=3/2))
  
  list(resMin=resMin, resFB=resFB, resMan=resMan, resLT=resLT, resKK=resKK)
}

Then the following call give us a list of result tables.

initialpt <- cbind(c(4, -4), c(-4, 4), c(3, 0), c(0, 3), c(-1, -1), c(0, 0))
mytablelist <- list()
for(i in 1: NCOL(initialpt))
{
    z0 <- c(initialpt[, i], 1, 1)
    mybench <- wholebench(z0)

    cat("z0", z0, "\n") 

    mytable12 <- data.frame(method=mybench[[1]]$compres[, 1], 
    round( 
        cbind(mybench[[1]]$compres[,c(-1, -4)], mybench[[2]]$compres[,c(-1, -4)])
        , 3) )

    mytable35 <- data.frame(method=mybench[[1]]$compres[, 1], 
    round( 
        cbind(mybench[[3]]$compres[,c(-1, -4)], mybench[[5]]$compres[,c(-1, -4)])
        , 3) )

    mytablelist <- c(mytablelist, z0=list(z0), MINFB=list(mytable12), MANKK=list(mytable35))
}

Note that one result table given by the function reports the computation results for 10 methods given an initial point and a complementarity function. Below an example

z0 <- c(-4, 4, 1, 1)
bench.GNE.nseq(z0, F, JacF, argPhi=list(phi=phiMin), 
               argjac=list(gphia= GrAphiMin, gphib= GrBphiMin), echo=FALSE)$compres

The following subsections report the computation for 4 complementarity functions, the Luo-Tseng being discarded due to non convergence. We also remove the final estimates zn when the method has not converged, ||Φ(zn)||2 ≠ 0. Tables are put in appendix, except the first one.

Initial point z0 = (4, −4, 1, 1)

We work on the initial point z0 = (4, −4, 1, 1), close the GNE (2, −2, 0, 160). Clearly, we observe the Mangasarian complementarity function ϕMan does not converge except in the pure Newton method, for which the sequence converges to (−2, 3, 8, 0) quite far from the initial point. So the sequence converged by a chance! For ϕMin function, when it converges, the GNEs found are (2, −2, 0, 160) or (1, 0, 512, 6). ϕFB and ϕKK associated sequences converge mostly to (2, −2, 0, 160). In terms of function/Jacobian calls, ϕFB is significantly better when used with the Newton scheme.

Initial point z0 = (−4, 4, 1, 1)

We work on the initial point z0 = (−4, 4, 1, 1), close the GNE (−2, 3, 8, 0). Again, we observe the Mangasarian complementarity function ϕMan does not converge. All other sequences converge the closest GNE (−2, 3, 8, 0). ϕMin sequence with Newton scheme is particularly good, then comes ϕFB and finally ϕKK.

Initial point z0 = (3, 0, 1, 1)

We work on the initial point z0 = (3, 0, 1, 1) close to the GNE (1, 0, 512, 6). As always, the sequence converges by chance with the pure Newton method to a GNE (−2, 3, 8, 0). Otherwise the other sequences, namely , and converges to the expected GNE. As the previous subsection, Broyden updates of the Jacobian is less performant than the true Jacobian (i.e. Newton scheme). The convergence speed order is preserved.

Initial point z0 = (0, 3, 1, 1)

We work on the initial point z0 = (0, 3, 1, 1) close to the GNE (0, 1, 324, 0). As always, the sequence converges by chance with the pure Newton method to a GNE (−2, 3, 8, 0). Others sequences have difficulty to converge the closest GNE. Local methods (i.e. pure) find the GNE (0, 1, 324, 0), while global version converges to (1, 0, 512, 6). It is logical any method will have difficulty to choose between these two GNEs, because they are close.

Initial point z0 = (−1, −1, 1, 1)

We work on the initial point z0 = (−1, −1, 1, 1) equidistant to the GNEs (0, 1, 324, 0) and (1, 0, 512, 6). Despite being closer to these GNEs, the pure Newton version of the sequence converges unconditionally to the GNE (−2, 3, 8, 0). All other sequences converges to the GNE (0, 1, 324, 0) except for the Broyden version of the sequence, converging to the farthest GNEs. In terms of function calls, the Newton line search version of the sequence is the best, followed by the Newton trust region version of the sequence.

Initial point z0 = (0, 0, 1, 1)

We work on the initial point z0 = (0, 0, 1, 1) equidistant to the GNEs (0, 1, 324, 0) and (1, 0, 512, 6). Both the and the sequences do not converge. The sequence diverges because the Jacobian at the initial point is exactly singular. Indeed, we have

z0 <- c(0, 0, 1, 1)
jacSSR(z0, dimx, dimlam, heobj=heobj, myarg, constr=g, grconstr=grg, 
    heconstr=heg, gcompla=GrAphiMin, gcomplb=GrBphiMin)
##      [,1] [,2] [,3] [,4]
## [1,]  512 1024    1    0
## [2,]    0    0    0    2
## [3,]   -1   -1    1    0
## [4,]    0    0    0    1

For the and sequences, we do not have this problem.

jacSSR(z0, dimx, dimlam, heobj=heobj, myarg, constr=g, grconstr=grg, 
    heconstr=heg, gcompla=GrAphiFB, gcomplb=GrBphiFB)
##             [,1]         [,2]       [,3]       [,4]
## [1,] 512.0000000 1024.0000000  1.0000000  0.0000000
## [2,]   0.0000000    0.0000000  0.0000000  2.0000000
## [3,]   0.2928932    0.2928932 -0.2928932  0.0000000
## [4,]   0.1055728    0.2111456  0.0000000 -0.5527864
jacSSR(z0, dimx, dimlam, heobj=heobj, myarg, constr=g, grconstr=grg, 
    heconstr=heg, gcompla=GrAphiKK, gcomplb=GrBphiKK, argcompl=3/2)
##             [,1]         [,2]       [,3]       [,4]
## [1,] 512.0000000 1024.0000000  1.0000000  0.0000000
## [2,]   0.0000000    0.0000000  0.0000000  2.0000000
## [3,]   0.2679492    0.2679492 -0.2679492  0.0000000
## [4,]   0.1101776    0.2203553  0.0000000 -0.4881421

So the sequence converge to a GNE, either (0, 1, 324, 0) or (−2, 3, 8, 0). Again the sequence converges faster.

Conclusions

In conclusion to this analysis with respect to initial point, the computation method and the complementarity function, we observe the strong difference in terms of convergence, firstly and in terms of convergence speed. Clearly the choice of the complementarity function is crucial, the Luo-Tseng and the Mangasarian are particularly inadequate in our example. Regarding the remaining three complementarity functions (the minimum, the Fisher-Burmeister and the Kanzow-Kleinmichel functions) generally converge irrespectively of the computation method. However, the sequences are particularly efficient and most of the time the Newton trust region method is the best in terms of function/Jacobian calls.

Special case of shared constraints with common multipliers

Let $h: \R^{n} \mapsto \R^{m_l}$ be a constraint function shared by all players. The total constraint function and the Lagrange multiplier for the ith player is $$ \tilde g^i(x) = \left( \begin{matrix} g^i(x) \\ h(x) \end{matrix} \right) \text{ and } \tilde \lambda^i = \left( \begin{matrix} \lambda^i \\ \mu \end{matrix} \right), $$ where $\mu \in \R^l$. This could fall within the previous framework, if we have not required the bottom part of λ̃i to be common among all players. The Lagrangian function of the ith player is given by $$ L ^i(x, \lambda^i, \mu) = O_i(x) + \sum_{k=1}^{m_i} g^i_k(x) \lambda_k^i + \sum_{p=1}^l h_p(x) \mu_p. $$

The generalized Jacobian is given in Appendix .

Constrained-equation reformulation of the KKT system

This subsection aims to present methods specific to solve constrained (nonlinear) equations, first proposed by in the GNEP context. The root function $H: \R^n \times \R^{2m} \mapsto \R^n \times \R^{2m}$ is defined as $$ H(x, \lambda, w) = \left( \begin{matrix} \tilde L(x, \lambda) \\ g(x) + w \\ \lambda \circ w \end{matrix} \right) , $$ where the dimensions n, m correspond to the GNEP notation (λ = (λ1, …, λN)) and (a, σ̄) is given by $((0_{n}, \II_{m}), 1)$. The potential function is given by $$ p\left(u \right) = \zeta \log\left( ||x||_2^2 + ||\lambda ||_2^2+ ||w||_2^2 \right) - \sum_{k=1}^{m} \log (\lambda_{k}) - \sum_{k=1}^{m} \log (w_{k}), $$ where $u=(x, \lambda, w) \in \R^n \times \R_{+}^{m} \times \R_{+}^{m}$ and ζ > m. The Jacobian is given in Appendix .

When there is a constraint function h shared by all players, the root function is given by $$ \widetilde H(x, \tilde \lambda, \tilde w) = \left( \begin{matrix} \bar L(x, \tilde \lambda) \\ \tilde g(x) + \tilde w \\ \tilde \lambda \circ \tilde w \end{matrix} \right) , \text{ with } \tilde \lambda = \left( \begin{matrix} \lambda^1 \\ \vdots \\ \lambda^N \\ \mu \end{matrix} \right) , \tilde w = \left( \begin{matrix} w^1 \\ \vdots \\ w^N \\ y \end{matrix} \right) \text{ and } \tilde g(x) = \left( \begin{matrix} g^1(x) \\ \vdots \\ g^N(x) \\ h(x) \end{matrix} \right) . $$ The Jacobian is given in Appendix .

A classic example

Using the classic example presented above, we get

Therefore, to compute a generalized Nash equilibrium, we use

z0 <- 1+rexp(sum(dimx)+2*sum(dimlam))
GNE.ceq(z0, dimx, dimlam, grobj=grobj, myarg, heobj=heobj, myarg, 
    constr=g, grconstr=grg, heconstr=heg, 
    method="PR", control=list(trace=0))
## GNE: 1.741725 -0.6156581 235.2884 32.61971 0.0002134447 0.001281449 
## with optimal norm 1.787033 
## after  100 iterations with exit code 4 .
## Output message: Iteration limit exceeded 
## Function/grad/hessian calls: 743 100 
## Optimal (vector) value: 0.8399171 -1.308634 0.1262801 0.8690728 0.05022106 0.0418005

GNEP as a fixed point equation or a minimization problem

We present another reformulation of the GNEP, which was originally introduced in the context of standard Nash equilibrium problem. The fixed-point reformulation arise from two different problem: either using the Nikaido-Isoda (NI) function or the quasi-varational inequaltiy (QVI) problem. We present both here. We also present a reformulation of the GNEP through a minimization problem. The gap minimization reformulation is closed linked to the fixed-equation reformulation.

NI reformulation

We define the Nikaido-Isoda function as the function ψ from $\R^{2n}$ to $\R$ by This function represents the unilateral player improvement of the objective function between actions x and y. Let V̂ be the gap function $$ \hat V(x) = \underset{ y \in X(x) }{\sup}~ \psi(x,y). $$ Theorem 3.2 of shows the relation between GNEPs and the Nikaido-Isoda function. If objective functions θi are continuous, then x⋆ solves the GNEP if and only if x⋆ is a minimimum of V̂ such that where the set $X(x) = \{y \in \R^n, \forall i, g^i(y_i, x_{-i}) \leq 0 \}$ and V̂ defined in (). Furthermore, the function V̂ is such that ∀x ∈ X(x), V̂(x) ≥ 0. There is no particular algorithm able to solve this problem for a general constrained set X(x). But a simplification will occur in a special case: the jointly convex case.

QVI reformulation

Assuming the differentiability of objective functions, the GNEP can be reformulated as a QVI problem
and a constrained set $X(x) = \{y \in \R^n, \forall i,~ g^i(y_i, x_{-i}) \leq 0 \}$. The following theorem states the equivalence between the GNEP and the QVI, see Theorem 3.3 of .

propose to refomulate the QVI problem as a minimization of a (regularized) gap function. The regularized gap function of the QVI () is $$ V_{QVI}(x) = \underset{y\in X(x)}{\sup}~\psi_{\alpha VI}(x, y), $$ where ψαVI is given by for a regularization parameter α > 0. Note that the minimisation problem appearing in the definition of VQVI is a quadratic problem. The theorem of given below shows the equivalence a minimizer of VQVI and the GNEP.

For each x ∈ X(x), the regularized gap function VQVI is non-negative VQVI(x) ≥ 0. If objective functions are continuous, then x⋆ solves the GNEP if and only if x⋆ is a minimum of VQVI such that

The jointly convex case

In this subsection, we present reformulations for a subclass of GNEP called jointly convex case. Firstly, the jointly convex setting requires that the constraint function is common to all players g1 = … = gN = g. Then, we assume, there exists a closed convex subset $X \subset \R^n$ such that for all player i, $$ \{y_i \in \R^{n_i}, g(y_i, x_{-i}) \leq 0 \} = \{y_i \in \R^{n_i}, (y_i, x_{-i}) \in X \} . $$ In our context parametrized context, the jointly convex setting requires that the constraint function is common to all players g1 = … = gN = g and is convex.

We consider the following example based on the previous example. Let N = 2. The objective functions are defined as θ1(x) = (x1 − 2)2(x2 − 4)4 and θ2(x) = (x2 − 3)2(x1)4, for $x\in \R^2$, while the constraint function g(x) = (g1(x), g2(x)) is given by g1(x) = x1 + x2 − 1 ≤ 0 and g2(x) = 2x1 + x2 − 2 ≤ 0. Objective functions can be rewritten as θi(x) = (xi − ci)2(x−i − di)4, with c = (2, 3) and d = (4, 0). First-order and second-order derivatives are given in the introduction. ∇jg1(x) = 1 and ∇jg2(x) = 2δj1 + δj2.

#O_i(x)
obj <- function(x, i, arg)
  (x[i] - arg$C[i])^2*(x[-i] - arg$D[i])^4
#g(x)
gtot <- function(x)
  sum(x[1:2]) - 1
#Gr_x_j g(x)
jacgtot <- function(x)
    cbind(1, 1)

z0 <- rexp(sum(dimx))

GNE.fpeq(z0, dimx, obj, myarg, grobj, myarg, heobj, myarg, gtot, NULL, 
         jacgtot, NULL, silent=TRUE, control.outer=list(maxit=10), 
         problem="NIR", merit="NI")
## GNE: 1.91041 -0.9104103 
## with optimal norm 1.372768e-07 
## after  iterations with exit code 1 .
## Output message: 
## Outer Function/grad/hessian calls: 5 3 
## Inner Function/grad/hessian calls: 2604 388
GNE.fpeq(z0, dimx, obj, myarg, grobj, myarg, heobj, myarg, gtot, NULL, 
         jacgtot, NULL, silent=TRUE, control.outer=list(maxit=10), 
         problem="VIR", merit="VI")
## GNE: -134.7119 135.7119 
## with optimal norm 7.205928e+22 
## after  iterations with exit code 6 .
## Output message: 
## Outer Function/grad/hessian calls: 19 10 
## Inner Function/grad/hessian calls: 454 148

NIF formulation for the jointly convex case

In the jointly convex case, the gap function becomes $$ V_{\alpha NI}(x)= \underset{ y \in X }{\max}~ \psi_{\alpha NI}(x,y). $$ Since y ↦ ψαNI(x, y) is strictly concave as long as objective functions θi are player-convex, the supremum is replaced by the maximum. Using two regularization parameters 0 < α < β, the constrained minimization problem can be further simplified to the unconstrained problem see .

Furthermore, a generalized equilibrium also solves a fixed-point equation, see Property 3.4 of . Assuming θi and g are C functions and g is convex and θi player-convex. x⋆ is a normalized equilibrium if and only if x⋆ is a fixed-point of the function where X is defined in () and ψαNI called the regularized Nikaido-Isoda function is defined as for a regularization parameter α > 0.

QVI formulation for the jointly convex case

The regularized gap function also simplifies and becomes $$ V_{\alpha VI}(x) = \underset{y\in X}{\sup}~\psi_{\alpha VI}(x, y), $$ where ψαVI is in (). Constrained equation () simplifies to a nonlinear equation VαVI(x⋆) = 0 and x⋆ ∈ X. Using two regularization parameters 0 < α < β, x⋆ is the global minimum of the unconstrained minimization problem

Furthermore, the VI reformulation leads to a fixed-point problem as shown in the following proposition. Assuming that θi and g are C functions, g is convex and θi player-convex, then x⋆ solves the VI (VαVI(x⋆) = 0 and x⋆ ∈ X) if and only if x⋆ is a fixed point of the function where X is defined in () and ψαVI is defined in ().

List of examples

We consider a two-player game defined by O1(x) = (x1 − 1)2 and O2(x) = (x2 − 1/2)2, with a shared constraint function g(x) = x1 + x2 − 1 ≤ 0. Solutions are given by (α, 1 − α) with α ∈ [1/2, 1] with Lagrange multipliers given by λ1 = 2 − 2α and λ2 = 2α − 1. But there is a unique normalized equilibrium for which λ1 = λ2 = 1/2. The nonsmooth reformulation of the KKT system uses the following terms ∇1O1(x) = 2(x1 − 1), ∇2O2(x) = 2(x2 − 1/2),  and ∇1g(x) = ∇2g(x) = 1. and ∇i2Oi(x) = 2, ∇j∇kOi(x) = 0,  and ∇j∇kg(x) = 0.

We consider a two-player game defined by Oi(x) = −(d − λ − ρ(x1 + x2))xi, with gi(x) = −xi ≤ 0, where d = 20, λ = 4, ρ = 1. Derivatives are given by ∇jOi(x) = −(−ρxi + (d − λ − ρ(x1 + x2))δij) and ∇jgi(x) = −δij, and ∇k∇jOi(x) = −(−ρδik − ρδij) and ∇k∇jgi(x) = 0. There is a unique solution given by x⋆ = (d − λ)/(3ρ).

We consider a two-player game defined by Oi(x) = −(d1 − d2(x1 + x2 + x3) − c1i − c2ixi)xi, and $$ g(x) = \left( \begin{matrix} \sum\limits_{l=1}^3 u_{l1} e_l x_l - K_1 \\ \sum\limits_{l=1}^3 u_{l2} e_l x_l - K_2 \end{matrix} \right). $$ Derivatives are given by $$ \nabla_j O_i(x) = - ( - d_2 - c_{2i} \delta_{ij})x_i - (d_1 - d_2 (x_1+x_2+x_3) - c_{1i} - c_{2i} x_i)\delta_{ij} \text{ and } \nabla_j g(x) = \left( \begin{matrix} u_{j1} e_j \\ u_{j2} e_j \end{matrix} \right), $$ and $$ \nabla_k \nabla_j O_i(x) = -( - d_2\delta_{ik} - d_2\delta_{ij} - 2 c_{2i} \delta_{ij}\delta_{ik}) \text{ and } \nabla_k \nabla_j g(x) = \left( \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} \right). $$

Tables for the nonsmooth reformulation

Appendix for the nonsmooth reformulation

The generalized Jacobian of the complementarity formulation has the following form $$ J(z) = \left( \begin{array}{c|c} \begin{matrix} \Jac_{x_1} L_1(x, \lambda^1) & \dots & \Jac_{x_N} L_1(x, \lambda^1) \\ \vdots & & \vdots \\ \Jac_{x_1} L_N(x, \lambda^N) & \dots & \Jac_{x_N} L_N(x, \lambda^N) \\ \end{matrix} & \begin{matrix} \Jac_{x_1} g^1(x)^T & & 0 \\ & \ddots & \\ 0 & & \Jac_{x_N} g^N(x)^T \\ \end{matrix} \\ \hline \begin{matrix} - D_1^a(x, \lambda^1) \Jac_{x_1} g^1(x) & \dots & - D_1^a(x, \lambda^1) \Jac_{x_N} g^1(x) \\ \vdots & & \vdots \\ - D_N^a(x, \lambda^N) \Jac_{x_1} g^N(x) & \dots & - D_N^a(x, \lambda^N) \Jac_{x_N} g^N(x) \\ \end{matrix} & \begin{matrix} D^b_{1}(x, \lambda^1) & & 0 \\ & \ddots & \\ 0 & & D^b_{N}(x, \lambda^N) \\ \end{matrix} \end{array} \right) . $$ The diagonal matrices Dia and Dib are given by $$ D_i^a(x, \lambda^i) = \diag[ a^i(x, \lambda^i) ] \text{ and } D_i^b(x, \lambda^i) = \diag[ b^i(x, \lambda^i) ], $$ with $a^i(x, \lambda^i), b^i(x, \lambda^i) \in \R^{m_i}$ defined as $$ ( a^i_j(x, \lambda_j^i) , b^i_j(x, \lambda_j^i) ) = \left\{ \begin{array}{cl} %\frac{ \partial \phi }{ \partial a} \left( \phi'_a( -g^i_j(x), \lambda_j^i ), \phi'_b ( -g^i_j(x), \lambda_j^i ) \right) & \text{if } ( -g^i_j(x), \lambda_j^i ) \neq (0, 0), \\ (\xi_{ij} , \zeta_{ij} ) & \text{if } ( -g^i_j(x), \lambda_j^i ) = (0, 0), \end{array} \right. $$ where ϕ′a (resp. ϕ′b) denotes the derivative of ϕ with respect to the first (second) argument a (b) and (ξij, ζij) ∈ B̄(pϕ, cϕ), the closed ball at pϕ of radius cϕ.

The generalized Jacobian of the complementarity formulation has the following form J(z)= $$ \left( \begin{array}{c|c|c} \begin{matrix} \Jac_{x_1} L_1(x, \lambda^1, \mu) & \dots & \Jac_{x_N} L_1(x, \lambda^1, \mu) \\ \vdots & & \vdots \\ \Jac_{x_1} L_N(x, \lambda^N, \mu) & \dots & \Jac_{x_N} L_N(x, \lambda^N, \mu) \\ \end{matrix} & \begin{matrix} \Jac_{x_1} g^1(x)^T & & 0 \\ & \ddots & \\ 0 & & \Jac_{x_N} g^N(x)^T \\ \end{matrix} & \begin{matrix} \Jac_{x_1} h(x)^T \\ \vdots \\ \Jac_{x_N} h(x)^T \\ \end{matrix} \\ \hline \begin{matrix} - D_1^a(x, \lambda^1) \Jac_{x_1} g^1(x) & \dots & - D_1^a(x, \lambda^1) \Jac_{x_N} g^1(x) \\ \vdots & & \vdots \\ - D_N^a(x, \lambda^N) \Jac_{x_1} g^N(x) & \dots & - D_N^a(x, \lambda^N) \Jac_{x_N} g^N(x) \\ \end{matrix} & \begin{matrix} D^b_{1}(x, \lambda^1) & & 0 \\ & \ddots & \\ 0 & & D^b_{N}(x, \lambda^N) \\ \end{matrix} & \begin{matrix} 0 \\ \vdots \\ 0 \\ \end{matrix} \\ \hline \begin{matrix} - D_h^a(x, \mu) \Jac_{x_1} h(x) & \dots & - D_h^a(x, \mu) \Jac_{x_N} h(x) \\ \end{matrix} & \begin{matrix} 0 & \dots & 0 \\ \end{matrix} & \begin{matrix} D^b_{h}(x, \mu) \\ \end{matrix} \end{array} \right) . $$ The diagonal matrices Da and Db are given by $$ D_h^a(x, \mu) = \diag[ \tilde a(x, \mu) ] \text{ and } D_h^b(x, \mu) = \diag[ \tilde b(x, \mu) ], $$ with $\tilde a(x, \mu), \tilde b(x, \mu) \in \R^{l}$ defined as $$ ( \tilde a_j(x, \mu) , \tilde b_j(x, \mu) ) = \left\{ \begin{array}{cl} %\frac{ \partial \phi }{ \partial a} \left( \phi'_a( -h_j(x), \mu_j ), \phi'_b ( -h_j(x), \mu_j ) \right) & \text{if } ( -h_j(x), \mu_j ) \neq (0, 0), \\ (\tilde \xi_{j} ,\tilde \zeta_{j} ) & \text{if } ( -h_j(x), \mu_j ) = (0, 0), \end{array} \right. $$ where (ξ̃j, ζ̃j) ∈ B̄(pϕ, cϕ).

For the line-search, the gradient ∇p is given by $$ \nabla p(x, \lambda, w)= \left( \begin{matrix} \frac{2 \zeta }{ ||x||_2^2 + ||\lambda ||_2^2+ ||w||_2^2 } x \\ \frac{2 \zeta }{ ||x||_2^2 + ||\lambda ||_2^2+ ||w||_2^2 } \lambda - \lambda^{-1} \\ \frac{2 \zeta }{ ||x||_2^2 + ||\lambda ||_2^2+ ||w||_2^2 } w - w^{-1} \\ \end{matrix} \right), $$ where λ and w have positive components and terms λ−1 and w−1 correspond to the component-wise inverse vector. Compared to the semismooth reformulation, the root function H is now C. The Jacobian is given by $$ \Jac H(x, \lambda, w) = \left( \begin{matrix} \Jac_x \tilde L(x, \lambda) & \diag\left[ \left( \nabla_{x_i} g^i(x) \right)_i \right] & 0 \\ \Jac_x g(x) & 0 & I \\ 0 & \diag[w] & \diag[\lambda] \end{matrix} \right). $$ As reported in , the computation of the direction dk = (dx, k, dλ, k, dw, k) can be simplified due to the special structure of the above Jacobian matrix. The system reduces to a linear system of n equations to find dx, k and the 2m components dλ, k, dw, k are simple linear algebra. Using the classic chain rule, the gradient of the merit function is given by $$ \nabla \psi(x, \lambda, w) = \Jac H(x, \lambda, w)^T \nabla p(H(x, \lambda, w)). $$ Again the computation of this gradient can be simplified due to the sparse structure of $\Jac H$.

The Jacobian is given by $$ \Jac \widetilde H(x, \tilde \lambda, \tilde w) = \left( \begin{matrix} \Jac_x \bar L (x, \tilde \lambda) & \Jac_{\tilde \lambda} \bar L (x, \tilde \lambda) & 0 \\ \Jac_x \tilde g(x) & 0 & I \\ 0 & \diag[\tilde w] & \diag[\tilde \lambda] \end{matrix} \right), $$ where $$ \Jac_{\tilde \lambda} \bar L (x, \tilde \lambda) = \left( \begin{matrix} \nabla_{x_1} g^1(x) & 0 & & \nabla_{x_1} h(x) \\ 0 & \ddots & 0 & \vdots \\ & 0 & \nabla_{x_N} g^N(x) & \nabla_{x_N} h(x) \\ \end{matrix} \right), $$ and $$ \Jac_x \tilde g(x) = \left( \begin{matrix} \Jac_x \tilde g^1(x) \\ \dots \\ \Jac_x \tilde g^N(x) \\ \Jac_x \tilde h(x) \end{matrix} \right). $$