Initial values for the parameters to be optimized over: z=(x, lambda, mu).
dimx
a vector of dimension for x.
dimlam
a vector of dimension for lambda.
grobj
gradient of the objective function (to be minimized), see details.
arggrobj
a list of additional arguments of the objective gradient.
heobj
Hessian of the objective function, see details.
argheobj
a list of additional arguments of the objective Hessian.
constr
constraint function (g^i(x)<=0), see details.
argconstr
a list of additional arguments of the constraint function.
grconstr
gradient of the constraint function, see details.
arggrconstr
a list of additional arguments of the constraint gradient.
heconstr
Hessian of the constraint function, see details.
argheconstr
a list of additional arguments of the constraint Hessian.
dimmu
a vector of dimension for mu.
joint
joint function (h(x)<=0), see details.
argjoint
a list of additional arguments of the joint function.
grjoint
gradient of the joint function, see details.
arggrjoint
a list of additional arguments of the joint gradient.
hejoint
Hessian of the joint function, see details.
arghejoint
a list of additional arguments of the joint Hessian.
method
a character string specifying the method
"PR" or "AS".
control
a list with control parameters.
...
further arguments to be passed to the optimization routine.
NOT to the functions H and jacH.
silent
a logical to get some traces. Default to FALSE.
Details
GNE.ceq solves the GNE problem via a constrained equation reformulation of the KKT system.
This approach consists in solving the extended Karush-Kuhn-Tucker
(KKT) system denoted by H(z)=0, for z in Ω where eqnz is formed by the players strategy
x, the Lagrange multiplier lambda and the slate variable w.
The root problem H(z)=0 is solved by an iterative scheme z_{n+1} = z_n + d_n,
where the direction d_n is computed in two different ways. Let J(x)=Jac H(x).
There are two possible methods either "PR" for potential reduction algorithm
or "AS" for affine scaled trust reduction algorithm.
(a) potential reduction algorithm:
The direction solves the system
H(z_n) + J(z_n) d = sigma_n a^T H(z_n) / ||a||_2^2 a.
(b) bound-constrained trust region algorithm:
The direction solves the system
min_p ||J(z_n)^T p + H(z_n)||^2 ,
for p such that ||p|| <= Delta_n||.
... are further arguments to be passed to the optimization routine,
that is global, xscalm, silent.
A globalization scheme can be choosed using the global argument.
Available schemes are
(1) Line search:
if global is set to "qline" or "gline", a line search
is used with the merit function being half of the L2 norm of Phi, respectively with a
quadratic or a geometric implementation.
(3) Trust-region:
if global is set to "pwldog", the Powell dogleg method
is used.
(2) None:
if global is set to "none", no globalization is done.
The default value of global is "gline" when method="PR" and
"pwldog" when method="AS".
The xscalm is a scaling parameter to used, either "fixed" (default)
or "auto", for which scaling factors are calculated from the euclidean norms of the
columns of the jacobian matrix.
The silent argument is a logical to report or not the optimization process, default
to FALSE.
The control argument is a list that can supply any of the following components:
xtol
The relative steplength tolerance.
When the relative steplength of all scaled x values is smaller than this value
convergence is declared. The default value is 1e-8.
ftol
The function value tolerance.
Convergence is declared when the largest absolute function value is smaller than ftol.
The default value is 1e-8.
btol
The backtracking tolerance.
The default value is 1e-2.
maxit
The maximum number of major iterations. The default value is 100 if a
global strategy has been specified.
trace
Non-negative integer. A value of 1 will give a detailed report of the
progress of the iteration, default 0.
sigma, delta, zeta
Parameters initialized to 1/2,
1, length(init)/2, respectively, when method="PR".
Parameters initialized to 0.99995, 1, 0.1, 1e10,
1/2, 1/4, 2, 3/4, when method="AS".
Value
GNE.ceq returns a list with components:
par
The best set of parameters found.
value
The value of the merit function.
counts
A two-element integer vector giving the number of calls to
H and jacH respectively.
iter
The outer iteration number.
code
The values returned are
1
Function criterion is near zero.
Convergence of function values has been achieved.
2
x-values within tolerance. This means that the relative distance between two
consecutive x-values is smaller than xtol.
3
No better point found.
This means that the algorithm has stalled and cannot find an acceptable new point.
This may or may not indicate acceptably small function values.
4
Iteration limit maxit exceeded.
5
Jacobian is too ill-conditioned.
6
Jacobian is singular.
100
an error in the execution.
message
a string describing the termination code.
fvec
a vector with function values.
Author(s)
Christophe Dutang
References
J.E. Dennis and J.J. Moree (1977),
Quasi-Newton methods, Motivation and Theory,
SIAM review.
Monteiro, R. and Pang, J.-S. (1999),
A Potential Reduction Newton Method for Constrained equations,
SIAM Journal on Optimization 9(3), 729-754.
S. Bellavia, M. Macconi and B. Morini (2003),
An affine scaling trust-region approach to bound-constrained nonlinear systems,
Applied Numerical Mathematics 44, 257-280
A. Dreves, F. Facchinei, C. Kanzow and S. Sagratella (2011),
On the solutions of the KKT conditions of generalized Nash equilibrium problems,
SIAM Journal on Optimization 21(3), 1082-1108.
See Also
See GNE.fpeq, GNE.minpb and GNE.nseq
for other approaches; funCER and
jacCER for template functions of H and Jac H.
Examples
#-------------------------------------------------------------------------------
# (1) Example 5 of von Facchinei et al. (2007)
#-------------------------------------------------------------------------------
dimx <- c(1, 1)
#Gr_x_j O_i(x)
grobj <- function(x, i, j)
{
if(i == 1)
res <- c(2*(x[1]-1), 0)
if(i == 2)
res <- c(0, 2*(x[2]-1/2))
res[j]
}
#Gr_x_k Gr_x_j O_i(x)
heobj <- function(x, i, j, k)
2 * (i == j && j == k)
dimlam <- c(1, 1)
#constraint function g_i(x)
g <- function(x, i)
sum(x[1:2]) - 1
#Gr_x_j g_i(x)
grg <- function(x, i, j)
1
#Gr_x_k Gr_x_j g_i(x)
heg <- function(x, i, j, k)
0
x0 <- rep(0, sum(dimx))
z0 <- c(x0, 2, 2, max(10, 5-g(x0, 1) ), max(10, 5-g(x0, 2) ) )
#true value is (3/4, 1/4, 1/2, 1/2)
GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj,
constr=g, grconstr=grg, heconstr=heg, method="PR",
control=list(trace=0, maxit=10))
GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj,
constr=g, grconstr=grg, heconstr=heg, method="AS", global="pwldog",
xscalm="auto", control=list(trace=0, maxit=100))
#-------------------------------------------------------------------------------
# (2) Duopoly game of Krawczyk and Stanislav Uryasev (2000)
#-------------------------------------------------------------------------------
#constants
myarg <- list(d= 20, lambda= 4, rho= 1)
dimx <- c(1, 1)
#Gr_x_j O_i(x)
grobj <- function(x, i, j, arg)
{
res <- -arg$rho * x[i]
if(i == j)
res <- res + arg$d - arg$lambda - arg$rho*(x[1]+x[2])
-res
}
#Gr_x_k Gr_x_j O_i(x)
heobj <- function(x, i, j, k, arg)
arg$rho * (i == j) + arg$rho * (j == k)
dimlam <- c(1, 1)
#constraint function g_i(x)
g <- function(x, i)
-x[i]
#Gr_x_j g_i(x)
grg <- function(x, i, j)
-1*(i == j)
#Gr_x_k Gr_x_j g_i(x)
heg <- function(x, i, j, k)
0
#true value is (16/3, 16/3, 0, 0)
x0 <- rep(0, sum(dimx))
z0 <- c(x0, 2, 2, max(10, 5-g(x0, 1) ), max(10, 5-g(x0, 2) ) )
GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj, arggrobj=myarg,
argheobj=myarg, constr=g, grconstr=grg, heconstr=heg,
method="PR", control=list(trace=0, maxit=10))
GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj, arggrobj=myarg,
argheobj=myarg, constr=g, grconstr=grg, heconstr=heg,
method="AS", global="pwldog", xscalm="auto", control=list(trace=0, maxit=100))
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(GNE)
Loading required package: alabama
Loading required package: numDeriv
Loading required package: nleqslv
Loading required package: BB
Loading required package: SQUAREM
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GNE/GNE.ceq.Rd_%03d_medium.png", width=480, height=480)
> ### Name: GNE.ceq
> ### Title: Constrained equation reformulation of the GNE problem.
> ### Aliases: GNE.ceq
> ### Keywords: nonlinear optimize
>
> ### ** Examples
>
>
>
> #-------------------------------------------------------------------------------
> # (1) Example 5 of von Facchinei et al. (2007)
> #-------------------------------------------------------------------------------
>
> dimx <- c(1, 1)
> #Gr_x_j O_i(x)
> grobj <- function(x, i, j)
+ {
+ if(i == 1)
+ res <- c(2*(x[1]-1), 0)
+ if(i == 2)
+ res <- c(0, 2*(x[2]-1/2))
+ res[j]
+ }
> #Gr_x_k Gr_x_j O_i(x)
> heobj <- function(x, i, j, k)
+ 2 * (i == j && j == k)
>
> dimlam <- c(1, 1)
> #constraint function g_i(x)
> g <- function(x, i)
+ sum(x[1:2]) - 1
> #Gr_x_j g_i(x)
> grg <- function(x, i, j)
+ 1
> #Gr_x_k Gr_x_j g_i(x)
> heg <- function(x, i, j, k)
+ 0
>
>
> x0 <- rep(0, sum(dimx))
> z0 <- c(x0, 2, 2, max(10, 5-g(x0, 1) ), max(10, 5-g(x0, 2) ) )
>
> #true value is (3/4, 1/4, 1/2, 1/2)
> GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj,
+ constr=g, grconstr=grg, heconstr=heg, method="PR",
+ control=list(trace=0, maxit=10))
GNE: 0.7499995 0.2499995 0.5000009 0.5000009 1.800985e-06 1.800985e-06
with optimal norm 1.800857e-06
after 9 iterations with exit code 1 .
Output message: Function criterion near zero
Function/grad/hessian calls: 9 9
Optimal (vector) value: 0 -1.110223e-16 9.003626e-07 9.003626e-07 9.004943e-07 9.004943e-07
>
>
> GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj,
+ constr=g, grconstr=grg, heconstr=heg, method="AS", global="pwldog",
+ xscalm="auto", control=list(trace=0, maxit=100))
GNE: 0.5185801 0.4814199 0.9628398 0.03716016 6.344505e-20 1.82146e-17
with optimal norm 1.321855e-16
after 9 iterations with exit code 1 .
Output message: Function criterion near zero
Function/grad/hessian calls: 9 9
Optimal (vector) value: -1.110223e-16 -6.938894e-17 6.344505e-20 1.82146e-17 6.108742e-20 6.768573e-19
>
>
> #-------------------------------------------------------------------------------
> # (2) Duopoly game of Krawczyk and Stanislav Uryasev (2000)
> #-------------------------------------------------------------------------------
>
>
> #constants
> myarg <- list(d= 20, lambda= 4, rho= 1)
>
> dimx <- c(1, 1)
> #Gr_x_j O_i(x)
> grobj <- function(x, i, j, arg)
+ {
+ res <- -arg$rho * x[i]
+ if(i == j)
+ res <- res + arg$d - arg$lambda - arg$rho*(x[1]+x[2])
+ -res
+ }
> #Gr_x_k Gr_x_j O_i(x)
> heobj <- function(x, i, j, k, arg)
+ arg$rho * (i == j) + arg$rho * (j == k)
>
>
> dimlam <- c(1, 1)
> #constraint function g_i(x)
> g <- function(x, i)
+ -x[i]
> #Gr_x_j g_i(x)
> grg <- function(x, i, j)
+ -1*(i == j)
> #Gr_x_k Gr_x_j g_i(x)
> heg <- function(x, i, j, k)
+ 0
>
> #true value is (16/3, 16/3, 0, 0)
>
> x0 <- rep(0, sum(dimx))
> z0 <- c(x0, 2, 2, max(10, 5-g(x0, 1) ), max(10, 5-g(x0, 2) ) )
>
>
> GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj, arggrobj=myarg,
+ argheobj=myarg, constr=g, grconstr=grg, heconstr=heg,
+ method="PR", control=list(trace=0, maxit=10))
GNE: 5.333334 5.333334 1.960473e-06 1.960473e-06 5.333344 5.333344
with optimal norm 2.091e-05
after 7 iterations with exit code 3 .
Output message: No better point found (algorithm has stalled)
Function/grad/hessian calls: 7 7
Optimal (vector) value: -3.094584e-15 -3.094584e-15 1.045412e-05 1.045412e-05 1.045588e-05 1.045588e-05
>
> GNE.ceq(z0, dimx, dimlam, grobj=grobj, heobj=heobj, arggrobj=myarg,
+ argheobj=myarg, constr=g, grconstr=grg, heconstr=heg,
+ method="AS", global="pwldog", xscalm="auto", control=list(trace=0, maxit=100))
GNE: 5.333333 5.333333 5.047138e-21 3.625697e-20 5.333333 5.333333
with optimal norm 2.512177e-15
after 8 iterations with exit code 1 .
Output message: Function criterion near zero
Function/grad/hessian calls: 8 8
Optimal (vector) value: -1.776362e-15 -1.776393e-15 0 0 2.691807e-20 1.933705e-19
>
>
>
>
>
>
>
>
> dev.off()
null device
1
>