GNE.nseq: Non smooth equation reformulation of the GNE problem.

Description

Non smooth equation reformulation via the extended KKT system of the GNE problem.

Usage

GNE.nseq(init, dimx, dimlam, grobj, arggrobj, heobj, argheobj, 
	constr, argconstr, grconstr, arggrconstr, heconstr, argheconstr,
	compl, gcompla, gcomplb, argcompl, 
	dimmu, joint, argjoint, grjoint, arggrjoint, hejoint, arghejoint, 
	method="default", control=list(), silent=TRUE, ...)

Value

GNE.nseq 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 phi and jacphi 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.

bench.GNE.nseq returns a list with components:

compres: a data.frame summarizing the different computations.
reslist: a list with the different results from GNE.nseq.

Arguments

init: 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.
compl: the complementarity function with (at least) two arguments: compl(a,b).
argcompl: list of possible additional arguments for compl.
gcompla: derivative of the complementarity function w.r.t. the first argument.
gcomplb: derivative of the complementarity function w.r.t. the second argument.
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 "Newton", "Broyden", "Levenberg-Marquardt" or "default" which is "Newton".
control: a list with control parameters.
...: further arguments to be passed to the optimization routine. NOT to the functions phi and jacphi.
silent: a logical to get some traces. Default to FALSE.

Author

Christophe Dutang

Details

Functions in argument must respect the following template:

constr must have arguments the current iterate z, the player number i and optionnally additional arguments given in a list.
grobj, grconstr must have arguments the current iterate z, the player number i, the derivative index j and optionnally additional arguments given in a list.
heobj, heconstr must have arguments the current iterate z, the player number i, the derivative indexes j, k and optionnally additional arguments given in a list.
compl, gcompla, gcomplb must have two arguments a, b and optionnally additional arguments given in a list.
joint must have arguments the current iterate z and optionnally additional arguments given in a list.
grjoint must have arguments the current iterate z, the derivative index j and optionnally additional arguments given in a list.
hejoint must have arguments the current iterate z, the derivative indexes j, k and optionnally additional arguments given in a list.

GNE.nseq solves the GNE problem via a non smooth reformulation of the KKT system. bench.GNE.nseq carries out a benchmark of the computation methods (Newton and Broyden direction with all possible global schemes) for a given initial point. bench.GNE.nseq.LM carries out a benchmark of the Levenberg-Marquardt computation method.

This approach consists in solving the extended Karush-Kuhn-Tucker (KKT) system denoted by $\Phi(z)=0$, where $z$ is formed by the players strategy $x$ and the Lagrange multiplier $\lambda$. The root problem $\Phi(z)=0$ is solved by an iterative scheme $z_{n+1} = z_n + d_n$, where the direction $d_n$ is computed in three different ways. Let $J(x)=Jac\Phi(x)$.

(a) Newton:: The direction solves the system $J(z_n) d = - \Phi(z_n) $, generally called the Newton equation.
(b) Broyden:: It is a quasi-Newton method aiming to solve an approximate version of the Newton equation $d = -\Phi(z_n) W_n$ where $W_n$ is computed by an iterative scheme. In the current implementation, $W_n$ is updated by the Broyden method.
(c) Levenberg-Marquardt:: The direction solves the system $$ \left[ J(z_n)^T J(z_n) + \lambda_n^\delta I \right] d = - J(z_n)^T\Phi(x_n) $$ where $I$ denotes the identity matrix, $\delta$ is a parameter in [1,2] and $\lambda_n = ||\Phi(z_n)|| $ if LM.param="merit", $||J(z_n)^T \Phi(z_n)|| $ if LM.param="jacmerit", the minimum of both preceding quantities if LM.param="min", or an adatpive parameter according to Fan(2003) if LM.param="adaptive".

In addition to the computation method, a globalization scheme can be choosed using the global argument, via the ... 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.
(2) Trust region:: if global is set to "dbldog" or "pwldog", a trust region is used respectively with a double dogleg or a Powell (simple) dogleg implementation. This global scheme is not available for the Levenberg-Marquardt direction.
(3) None:: if global is set to "none", no globalization is done.

The default value of global is "gline". Note that in the special case of the Levenberg-Marquardt direction with adaptive parameter, the global scheme must be "none".

In the GNEP context, details on the methods can be found in Facchinei, Fischer & Piccialli (2009), "Newton" corresponds to method 1 and "Levenberg-Marquardt" to method 3. In a general nonlinear equation framework, see Dennis & Moree (1977), Dennis & Schnabel (1996) or Nocedal & Wright (2006),

The implementation relies heavily on the nleqslv function of the package of the same name. So full details on the control parameters are to be found in the help page of this function. We briefly recall here the main parameters. 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 $10^{-8}$.
ftol: The function value tolerance. Convergence is declared when the largest absolute function value is smaller than ftol. The default value is $10^{-8}$.
delta: A numeric delta in [1, 2], default to 2, for the Levenberg-Marquardt method only.
LM.param: A character string, default to "merit", for the Levenberg-Marquardt method only.
maxit: The maximum number of major iterations. The default value is 150 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.

... are further arguments to be passed to the optimization routine, that is global, xscalm, silent. See above for the globalization scheme. 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. See nleqslv for details. The silent argument is a logical to report or not the optimization process, default to FALSE.

References

J.E. Dennis and J.J. Moree (1977), Quasi-Newton methods, Motivation and Theory, SIAM review.

J.E. Dennis and R.B. Schnabel (1996), Numerical methods for unconstrained optimization and nonlinear equations, SIAM.

F. Facchinei, A. Fischer and V. Piccialli (2009), Generalized Nash equilibrium problems and Newton methods, Math. Program.

J.-Y. Fan (2003), A modified Levenberg-Marquardt algorithm for singular system of nonlinear equations, Journal of Computational Mathematics.

B. Hasselman (2011), nleqslv: Solve systems of non linear equations, R package.

A. von Heusinger and C. Kanzow (2009), Optimization reformulations of the generalized Nash equilibrium problem using Nikaido-Isoda-type functions, Comput Optim Appl .

J. Nocedal and S.J. Wright (2006), Numerical Optimization, Springer Science+Business Media

Examples

Run this code



#-------------------------------------------------------------------------------
# (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



#true value is (3/4, 1/4, 1/2, 1/2)

z0 <- rep(0, sum(dimx)+sum(dimlam))

funSSR(z0, dimx, dimlam, grobj=grobj, constr=g, grconstr=grg, compl=phiFB, echo=FALSE)

	
jacSSR(z0, dimx, dimlam, heobj=heobj, constr=g, grconstr=grg, 
	heconstr=heg, gcompla=GrAphiFB, gcomplb=GrBphiFB)


GNE.nseq(z0, dimx, dimlam, grobj=grobj, NULL, heobj=heobj, NULL, 
	constr=g, NULL, grconstr=grg, NULL, heconstr=heg, NULL, 
	compl=phiFB, gcompla=GrAphiFB, gcomplb=GrBphiFB, method="Newton", 
	control=list(trace=1))

GNE.nseq(z0, dimx, dimlam, grobj=grobj, NULL, heobj=heobj, NULL, 
	constr=g, NULL, grconstr=grg, NULL, heconstr=heg, NULL, 
	compl=phiFB, gcompla=GrAphiFB, gcomplb=GrBphiFB, method="Broyden", 
	control=list(trace=1))



#-------------------------------------------------------------------------------
# (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) 

z0 <- rep(0, sum(dimx)+sum(dimlam))

funSSR(z0, dimx, dimlam, grobj=grobj, myarg, constr=g, grconstr=grg, compl=phiFB, echo=FALSE)

jacSSR(z0, dimx, dimlam, heobj=heobj, myarg, constr=g, grconstr=grg, 
	heconstr=heg, gcompla=GrAphiFB, gcomplb=GrBphiFB)


GNE.nseq(z0, dimx, dimlam, grobj=grobj, myarg, heobj=heobj, myarg, 
	constr=g, NULL, grconstr=grg, NULL, heconstr=heg, NULL, 
	compl=phiFB, gcompla=GrAphiFB, gcomplb=GrBphiFB, method="Newton", 
	control=list(trace=1))

GNE.nseq(z0, dimx, dimlam, grobj=grobj, myarg, heobj=heobj, myarg, 
	constr=g, NULL, grconstr=grg, NULL, heconstr=heg, NULL, 
	compl=phiFB, gcompla=GrAphiFB, gcomplb=GrBphiFB, method="Broyden", 
	control=list(trace=1))

Run the code above in your browser using DataLab