数理解析研究所講究録
第 829 巻 1993 年 152-167

152

Optimization-Based Iterative Methods for Solving
Nonlinear Complementarity Problems
TAJI Kouichi (田地宏一)
FUKUSHIMA Masao (福島雅夫)
IBARAKI Toshihide (茨木俊秀)
Department of Applied Mathematics and Physics,
Faculty of Engineering, Kyoto University, Kyoto 606, JAPANIntroduction

1

We consider the nonlinear complementarity problem, which is to find a vector $x\in R^{n}$such that$x\geq 0,$ $F(x)\geq 0$

and $x^{T}F(x)=0$ ,(1)where $F(x)=(F_{1}(x), F_{2}(x),$is a given continuously differentiable mappingfrominto itself and denotes transposition. This problem has been used to formulate
and study various equilibrium problems including the traffic equilibrium problem, thespatial economic equilibrium problem and Nash equilibrium problem [4, 12, 16].
The nonlinear complementarity problem (1) is a special case of the variational inequality
problem, which is to find a vector $x^{*}\in S$ such that$\ldots,$$R^{n}$$F_{n}(x))^{T}$$T$$(x-x^{*})^{T}F(x^{*})\geq 0$for 瓠 1 $x\in S$,(2)

is a nonempty closed convex set in. Problem (1) corresponds to the case
$S=R_{+}^{n}$ , the nonnegative orthant of. The variational inequality problem has also beenused to formulate and study various equilibrium problems arising in economics, operations
research, transportation and regional sciences.
To solve the nonlinear complementarity problem (1) and the variational inequality problem (2), various iterative algorithms, such as fixed point algorithms, projection methods,
nonlinear Jacobi method, successive over-relaxation methods and Newton method, havebeen proposed [6, 7, 13]. Among these, fixed point algorithms originally have been used
where

$S$

$R^{n}$$R^{n}$

153
in the proof that the problem has a solution. Fixed point algorithms are useless for practical computations because their convergence are extremely slow [14]. The other methods
are generalizations of methods for systems of nonlinear equations and their convergence
results have been obtained $[7, 13]$ . But, in general, these methods do not have globally
convergent property.
Assuming the monotonicity of mapping , Fukushima [5] has recently proposed a differentiable optimization formulation for variational inequality problem (2) and proposeda decent algorithm to solve (2). Based on this optimization formulation, Taji et al. [15]proposed a modification of Newton method for solving the variational inequality problem(2), and proved that, under the strong monotonicity assumption, the method is globallyand quadratically convergent.
In this paper we apply the methods of Fukushima [5] and Taji et al. [15] to the nonlinear complementarity problem. We show that those specialized methods can take fulladvantage of the special structure of problem (1), thereby yielding new globally convergent algorithms for solving monotone complementarity problems. We remark that theconstraint set of nonlinear complementarity problem iswhich is clearly not compact, while the descent method by Fukushima [5] has assumed that the constraint set iscompact. In this paper we show that the compactness assumption can be removed forthe nonlinear complementarity problem. We also present some computational results todemonstrate that those methods are practically efficient.
$F$$R_{+}^{n}$2

Equivalent optimization problem

We first review the results obtained by Fukushima [5] for the general variational inequality problem (2). Let be an $nxn$ symmetric positive definite matrix. The projection
under the G-norm of $x\in R^{n}$ onto the closed convex set , denoted Proj $s,c(x)$ , is definedas the unique solution of the following mathematical program:
$G$

$S$

minimize$||y-x||_{G}$subject to $y\in S$,, which is defined by 11
wheredenotes the G-norm in
Using this notation, we define function : $R^{n}arrow R$ by$||\cdot||_{G}$$R^{n}$$x||_{G}=(x^{T}Gx)^{\frac{1}{2}}$.

$f$$f(x)= \max\{-(y-x)^{T}F(x)-\frac{1}{2}(y-x)^{T}G(y-x)|y\in S\}$

(3)$=-(H(x)-x)^{T}F(x)- \frac{1}{2}(H(x)-x)^{T}G(H(x)-x)$ ,(4)

where the mapping $H:R^{n}arrow R^{n}$ is defined by$H(x)=Proj_{S,G}(x-G^{-1}F(x))$ .

(5)

154
It can be shown [5] that the function is continuously differentiable whenever so is the
mapping $F$ , and its gradient is given by$f$$\nabla f(x)=F(x)-[\nabla F(x)-G](H(x)-x)$ .

(6)

The function has the property that $f(x)\geq 0$ for $aUx\in S$ and $f(x)=0$ wheneveris a solution to the variational inequality problem (2). Hence problem (2) is equivalent
to the following optimization problem:
$f$

$x$

minimize $f(x)$ subject to $x\in S$.

(7)

In addition, when $\nabla F(x)$ is positive definite for all $x\in S$ , it can be shown that, if$x\in S$ is a stationary point of problem (7), i.e.,$(y-x)^{T}\nabla f(x)\geq 0$

for all $y\in S$,(8)then is a global optimal solution of problem (7), and hence it solves the nonlinear
complementarity problem (1).
Let us now specialize the above results to the complementarity problem. For generalvariational inequality problems, it may be expensive to evaluate the function unlessis tractable. In the case of the complementarity problem (1), the set turns out to be, where are. So, if we in particular let be a diagonal matrix$H$
positive constants, then the mapping takes the explicit form
$x$

$S$

$f$

$S$

$G$

$R_{+}^{n}$$D=diag(\delta_{1}, \ldots, \delta_{n})$$\delta_{:}$(9)$H(x)= \max(0, x-D^{-1}F(x))$

where maximum operator is taken component-wise, i.e.
$H_{i}(x)= \max(0, x_{i}-\delta_{i}^{-1}F_{i}(x))$

Hence the function$f$,$i=1,$$\ldots,$

$n$

.

(10)

and its gradient can be represented as$f(x)= \frac{1}{2}F(x)^{T}D^{-1}F(x)-\frac{1}{2}\max(0, D^{-1}F(x)-x)^{T}D\max(0, D^{-1}F(x)-x)$$= \sum_{=1}^{n}\frac{1}{2\delta_{1}}\{F_{i}(x)^{2}-(\max(0, F_{1}(x)-\delta_{i}x_{i}))^{2}\}$

,

(11)

and
(12)

$\nabla f(x)=\nabla F(x)D^{-1}F(x)+(I-\nabla F(x)D^{-1})\max(0, F(x)-Dx)$

respectively, and the optimization problem (7) becomesminimize $f(x)$ subject to $x\geq 0$ .
In this special case, it is therefore straight forward to evaluate the function
gradient. Furthermore, the optimization problem (13) has a simple constraint.

(13)
$f$

and its

155
Descent methods

3

In this section, we specialize the methods of Fukushima [5] and Taji et al. [15], whichwere originally proposed for variational inequality problems, to the nonlinear complementarity problem (1). Throughout this section, we let $D$ be a positive definite diagonal
matrix and the function be defined by (11). We so suppose that the mapping iswith modulus $\mu>0$ i.e.,strongly monotone on
$al$

$f$

$F$

$R_{+}^{n}$$(x-y)^{T}(F(x)-F(y))\geq\mu||x-y||^{2}$for all$x,$ $y\geq 0$

.

(14)

The first method uses the vector$d=$$H(x)-x$(15)$= \max(0, x-D^{-1}F(x))-x$

as a search direction at . When the mapping is strongly monotone with modulus ,it is shown [5] that the vector given by (15) satisfies the descent condition
$F$

$x$

$\mu$

$d$

$d^{T}\nabla f(x)\leq-\mu||d||^{2}$

Thus the direction can be used to determine the next iterate by using the following, whereis the smallest nonnegative integer
Armijo-type line search rule: Let$m$ such that$d$$\alpha$$:=\beta^{\dot{m}}$$f(x)-f(x+\beta^{m}d)\geq\sigma\beta^{m}||d||^{2}$$\hat{m}$,where $0<\beta<1$ and $0<\sigma$ .

Algorithm 1 a:choose$k$$x^{0}\geq 0,$$\beta\in(0,1)$,$\sigma>0$

and a positive diagonal matrix $D$ ;$:=0$while convergence criterion is not satisfied do$d^{k}$$:= \max(0, x^{k}-D^{-1}F(x^{k}))-x^{k}$

$m$

$:=0$

;

while $f(x^{k})-f(x^{k}+\beta^{m}d^{k})<\sigma\beta^{m}||d^{k}||^{2}$ do$m$$:=m+1$endwhile$x^{k+1}$$k$$:=x^{k}+\beta^{m}d^{k}$$:=k+1$endwhile;

156
In line search procedure of this algorithm, we examine only the points shorter than unit
step size. But we expect to decrease the value of the function when the longer step size
is chosen. Therefore we propose the algorithm in which we modify Algorithm la so thatthe longer step size can be selected.
$f$Algorithm 1 :$b$choose$k$$x^{0}\geq 0$,$\beta_{1}>1$,$\beta_{2}\in(0,1)$,$\sigma>0$

and a positive diagonal matrix

$D$

;

$:=0$

while convergence criterion is not satisfied do$:= \max(0, x^{k}-D^{-1}F(x^{k}))-x^{k}$ ;$:= \max\{t|x^{k}+td^{k}\geq 0, t\geq 0\}$ ;$d^{k}$$\hat{t}$

$m$

$:=0$

if $f(x^{k})-f(x^{k}+d^{k})\geq\sigma||d^{k}||^{2}$ then
while
and
and

$f(x^{k})-f(x^{k}+\beta_{1}^{m}d^{k})\geq\sigma\beta_{1}^{m}||d^{k}||^{2}$$\beta_{1}^{m}\leq\hat{t}$$f(x^{k}+\beta_{1}^{m+1}d^{k})\leq f(x^{k}+\beta_{1}^{m}d^{k})$

$m$

do

$:=m+1$endwhile$x^{k+1}$$:=x^{k}+\beta_{1}^{m}d^{k}$

else

while
$m$

$f(x^{k})-f(x^{k}+\beta_{2}^{m}d^{k})<\sigma\beta_{2}^{m}||d^{k}||^{2}$

do$:=m+1$endwhile$x^{k+1}$$:=x^{k}+\beta_{2}^{m}d^{k}$

endif
$k$

$:=k+1$endwhile
Note that, since an evaluation of at a given point is equivalent to evaluating thevector $\max(O, x-D^{-1}F(x))$ , the vector $H(x^{k})= \max(0, x^{k}-D^{-1}F(x^{k}))$ has already beenfound at the previous iteration as a by-product of evaluating $f(x^{k}+\beta^{m}d^{k})$ . Therefore
at the beginning of each iteration.
one need not compute the search direction
In the descent method proposed by Fukushima, to prove that the method is globally
convergent it is necessary not only the st ong monotonicity of mapping but the compactness of constraint set. However, in the case of nonlinear complementarity problem, we
can prove that Algorithms la and lb convergence globally if only is strongly monotone.
$f$

$x$

$d^{k}$$r$$F$

Proposition 3.1 If$\lim$$x\geq 0,$$||x||arrow\infty$$F$

is strongly monotone on$f(x)=+\infty$.

$R_{+}^{n}$, then157
Proof. For convenience, we define$f_{i}(x)=F_{i}(x)^{2}-( \max(0, F_{i}(x)-\delta_{i}x_{i}))^{2}$

hence

$f$

is written as$f(x)= \sum_{=1}^{n}\frac{1}{2\delta_{i}}f_{i}(x)$

,

(16)

. We first show that$f_{i}(x)\geq 0$

for all $x\geq 0$ . If

, then $f_{i}(x)=F_{1}(x)^{2}\geq 0$ . So we consider the case $F_{i}(x)-\delta_{i}x_{i}>0$ . Since; hold, we see, from (16),and$F_{i}(x)-\delta_{i}x_{i}\leq 0$$6_{i}>0,$$x_{i}\geq 0$$f_{t}(x)$$F_{i}(x)>\delta_{*}x$$=$$F_{i}(x)^{2}-(F_{1}(x)-\delta_{i}x_{i})^{2}$$=$$\delta_{i}x_{i}(2F_{1}(x)-\delta_{i}x_{i})$$\geq$$(\delta_{i}x_{i})^{2}$$\geq$$0$.

Let
be a sequence such that $x^{k}\geq 0$ and Il. Taking a subsequence, if
for $i\in I$ andisnecessary, there exists a set $I\subset\{1,2, \ldots, n\}$ such that
. Without loss of generality,
itself has a such set . From,bounded for
where
if $i\in I$ andif. From (14) and thewe define a sequence
definition of , we have$x^{k}||arrow\infty$$\{x^{k}\}$$x_{i}^{k}arrow+\infty$$I$$\{x^{k}\}$$i\not\in I$$y_{i}^{k}=0$$\{y^{k}\}$$\{x_{i}^{k}\}$$y_{i^{k}}=x_{i}^{k}$$\{x^{k}\}$$i\not\in I$$y^{k}$$\sum_{i\in I}(F_{1}(x^{k})-F_{1}(y^{k}))x^{\dot{k}}\geq\mu\sum_{:\in I}x^{k^{2}}$By Cauchy’s inequality$||F(x^{k})-F(y^{k})$$||||$$x^{k}-y^{k}||\geq(F(x^{k})-F(y^{k}))^{T}(x^{k}-y^{k})$,we have$( \sum_{i\in I}(F_{i}(x^{k})-F_{1}(y^{k}))^{2})^{\frac{1}{2}}(\sum_{i\in I}x^{k^{2}})^{\frac{1}{2}}\geq\sum_{i\in I}(F_{1}(x^{k})-F_{i}(y^{k}))x^{k}$,hence we have(17)$\sum_{:\in I}(F_{i}(x^{k})-F_{i}(y^{k}))^{2}\geq\mu^{2}\sum_{i\in I}x_{i}^{k^{2}}$Since, from the definition ofis bounded, $\{F_{i}(y^{k})\}are$ also bounded for all , andfor all $i\in I,$ (17) implies$y^{k},$$x^{k};arrow+\infty$$\sum_{i\in I}F_{1}(x^{k})^{2}arrow\infty$.

$\{y^{k}\}$

$i$

158
As shown at the beginning of the proof, $f_{1}(x)=F_{i}(x)^{2}\geq 0$ if
if $F_{1}(x)-\delta_{i}x_{i}>0$ . Therefore we have$F_{i}(x)-\delta_{i}x_{i}\leq 0$, and$f(x^{k})\geq(\delta_{i}x_{i})^{2}$$f(x^{k})$Since$= \sum_{=1}^{n}\frac{1}{2\delta_{i}}f_{1}(x)$$\geq$$\sum_{i\in I}\frac{1}{2\delta_{1}}f_{i}(x)$$\geq$$\sum_{:\in I}\frac{1}{2\delta_{*}}\min(F_{1}(x)^{2}, (\delta_{i}x.)^{2})$$x_{1}^{k}arrow+\infty$

for all $i\in I$ and.

$\sum_{i\in I}F_{1}(x^{k})^{2}arrow\infty$, it follows that $f(x^{k})arrow+\infty$ .

$\square$is continuously differentiable and strongly
Theorem 3.1 Suppose that the mapping
is Lipschitz continuous on any. Suppose also thatmonotone with modulus on, the sequence
generated. Then, for any starting point
bounded subset ofby Algorithm la or Algorithm $lb$ converges to the unique solution of problem (1) if the.
positive constant is chosen to be suff ciently small such that$F$$\nabla F$$R_{+}^{n}$$\mu$$x^{0}\geq 0$$R_{+}^{n}$$\{x^{k}\}$$\sigma<\mu$$\sigma$

Proof. By proposition 3.1, the level set $S=\{x|f(x)\leq f(x^{0})\}$ is shown to be bounded.
is Lipschitz continuous on . Since is continuously differentiable it is easy to
Henceshow that is also Lipschitz continuous on . Under these conditions, it is not difficult
is Lipschitz continuous on , i.e., there exists a constant $L>0$ suchto show that
that
$F$

$S$

$\nabla F$

$S$

$F$

$S$

$\nabla f$for all$||\nabla f(x)-\nabla f(y)||\leq L||x-y||$$x,$$y\in S$.

The remainder of this proof is the same as the proof of [5, theorem4.2].
The second method is a modification of Newton method, which incorporates the line
search strategy. The original Newton method for solving the nonlinear complementarity
is determined assuch that $x^{0}\geq 0$ andproblem (1) generates a sequence, where is a solution to the following linearized complementarity problem:$\square$$x^{k+1}$$\{x^{k}\}$$x^{k+1}$$:=\overline{x}$$\overline{x}$$x\geq 0,$ $F(x^{k})+\nabla F(x^{k})^{T}(x-x^{k})\geq 0$

and $x^{T}(F(x^{k})+\nabla F(x^{k})^{T}(x-x^{k}))=0$ . (18)

It is shown [13] that, under suitable as sumptions, the sequence generated by (18) quadratof the problem (1), provided that the starting point
ically converges to a solution
is chosen sufficiently close to . Taji et al. [15] have shown that, when the mapping
is strongly monotone with modulus , the vector
obtained by solving thelinearized complementarity problem (18) satisfies the inequality$x^{0}$$x^{*}$$x^{*}$

$F$

$\mu$

$d^{T}\nabla f(x^{k})<-$($\mu-\frac{1}{2}$

II

$D||$

) $||d||^{2}$$d:=\overline{x}-x^{k}$

159
Therefore, is actually a feasible descent direction of at , if the matrix $D$ is chosen
to satisfy $||D||= \max_{i}(\delta_{1})<2\mu$ . Using this result, we can construct a modified Newton
method for solving the nonlinear complementarity problem (1).
$d$

$f$

$x^{k}$

Algorithm 2:choose$k$$x^{0}\geq 0$,$\beta\in(0,1)$,$\gamma\in(0,1)$,$\sigma\in(0,1)$

and a positive diagonal matrix

$D$

;

$:=0$

while convergence criterion is not satisfied dofind hi such that$x^{T}(F(x^{k})+\nabla F(x^{k})^{T}(x-x^{k}))=0$ ;$x\geq 0,$ $F(x^{k})+\nabla F(x^{k})^{T}(x-x^{k})\geq 0$ and$d^{k}$$:=\overline{x}-x^{k}$

if $f(x^{k}+d^{k})\leq\gamma f(x^{k})$ then$\alpha_{k}$

$:=1$

else
$m$

$:=0$

while $f(x^{k})-f(x^{k}+\beta^{m}d^{k})<-\sigma\beta^{m}d^{k^{T}}\nabla f(x^{k})$ do$m$$:=m+1$endwhile$\alpha_{k}$$:=\beta^{m}$endif$x^{k+1}$$k$$:=x^{k}+\alpha_{k}d^{k}$;$:=k+1$endwhile
By applying the results proved in [15] for the general variational inequality problem, we
see that, when the mapping is strongly monotone with modulus , this algorithm isconvergent to the solution of (1) for any $x^{0}\geq 0$ if the matrix $D$ is chosen such that$||D||= \max(\delta_{i})<2\mu$ . From [15], we also see that the rate of convergence is quadratic ifof problem (1)is Lipschitz continuous on a neighborhood of the unique solution
and the strict complementarity condition holds at $x^{*}i.e.,$ $x^{*}=0$ implies $F_{i}(x^{*})>0$ for.
all $i=1,2,$$F$$\nabla F$$x^{*}$$\ldots,$

4

$\mu$

$n$

Computational results

In this section, we report some numerical results for Algorithms la, lb and 2 discussed
in the previous section. All computer programs were coded in FORTRAN and the runs
were made in double precision on a personal computer called Fujitsu FMR-70.

160
Throughout the $co$ mputational experiments, the parameters used in algorithms wereset as $\beta=\beta_{1}=\beta_{2}=0.5,$ $\gamma=0.5$ and $\sigma=0.0001$ . The positive diagonal matrix was. Therefore thechosen to be the identity matrix multiplied by a positive parameter
merit function (11) can be written more simply as$D$$\delta>0$$f(x)= \frac{1}{2\delta}\sum_{i=1}^{n}\{F_{1}(x)^{2}-(\max(0^{\cdot}, F_{1}(x)-\delta x_{i}))^{2}\}$

.

(19)

The search direction of Algorithms la and lb can also be written as$d^{k}$$:= \max(0,$ $x^{k}- \frac{1}{\delta}F(x^{k}))-x^{k}$.

The convergence criterion was$| \min(x;, F_{i}(x))|\leq 10^{-5}$for all $i=1,2,$$\ldots,$

$n$

.

(20)

For comparison purposes, we also tested two popular methods for solving the nonlinear
complementarity problem, the projection method [3] and the Newton method without
such thatandline search [9]. The projection method generates a sequenceby
is determined from$x^{0}\geq 0$$\{x^{k}\}$$x^{k+1}$$x^{k}$$x^{k+1}$$:= \max(0,$ $x^{k}- \frac{1}{\delta}F(x^{k}))$ ,(21)

for all . Note that this method may be considered a fixed step-size variant of Algorithms
la and lb. When the mapping $F$ is strongly monotone and Lipschitz continuous withconstants and , respectively, this method is globally convergent if 6 is chosen large
enough to satisfy $\delta>L^{2}/2\mu$ (see [13, Corollary 2.11.]).
The mappings used in our numerical experiments are of the form
$k$

$L$

$\mu$

$F(x)=Ix+\rho(N-N^{T})x+\phi(x)+c$ ,(22)

where is the $nxn$ identity matrix, $N$ is an $nxn$ matrix such that each row contains
is a nonlinear monotone mapping with components
only one nonzero element, and, where ; are positive constants. Elements of matrix $N$ and vector as
are randomly generated such that $-5\leq N_{*j}\leq 5,$ $-25\leq c_{t}\leq 25$
well as coefficientsand 0.001 $\leq p;\leq 0.006$ . The results are shown in Tables $1\sim 4$ . All starting points
is the total number of evaluating thewere chosen to be $(0,0, \ldots, 0)$ . In the tables,
merit function and all CPU times are in seconds and exclude input/output times. The
parameter is used to change the degree of asymmetry of , namely
deviates fromsymmetry as becomes large. Since the matrix $I+\rho(N-N^{T})$ is positive definite for any
are monotonically increasing for $x;\geq 0$ , the mapping defined by (22) isand.
strongly monotone on$I$$\phi(x)$$\phi_{i}(x_{i})=p_{i}x_{i}^{4}$

$c$

$p$

$p_{i}$$\# f$$f$$F$$\rho$$F$$\rho$$\rho$$F$$\phi_{1}(x_{i})$$R_{+}^{n}$

161
4.1

Comparison of Algorithms la, lb and the projection method

First we have compared Algorithms la, lb and the projection method (21) by using a
10-dimensional example, in which mapping is given by$F$

The results for this problem are shown in Table 1.
In general, the projection method is guaranteed to converge, only if the parameter 6
is chosen sufficiently large. In fact, Table 1 shows that when 6 is large, the projection
method is always convergent, but as becomes small, the behavior of the method becomes
unstable and eventually it fails to converge.
Table 1 also shows that Algorithms la and lb are always convergent even if 6 is chosen
small, since line search determines an adequate step size at each iteration. Note that,in Algorithm lb, the number of iterations is almost constant. This is because we may
choose a larger step size when the magnitude of vector
is small, i.e. 6 is large. This is
in contrast with Algorithm la, in which step size is bounded by 1 so that the number of
iterations increases as becomes large.
Algorithms la and lb spend more CPU times per iteration than the projection method,
because the former algorithms require overheads of evaluating the merit function . Moreover, when is between 0.1 and 20, Algorithm lb spends more CPU time than Algorithmla, though the number of iterations are almost equal. This is because Algorithm lb
attempts to find a larger step size at each iteration. But, when 6 becomes large, Algorithm lb tends to spend less CPU time than not only Algorithm la but also the projection
method, because the number of iterations of Algorithm lb does not increase so drastically.
$\delta$$d^{k}$$\delta$$f$$\delta$

4.2

Comparison of Algorithm 2 and Newton method

Next we have compared Algorithm 2 and the pure Newton method (18) without line
search. For each of the problem sizes $n=30,50$ and 90, we randomly generated five test
problems. The parameters and were set as $\rho=1$ and $\delta=10$ . The starting point waschosen to be $x=0$ . In solving the linearized subproblem at each iteration of Algorithm
2 and Newton method, we used Lemke’s complementarity pivoting method [10] codedby Fukushima [8]. The results are given in Table 2. All numbers shown in Table 2 are$\delta$$\rho$

162
averages of the results for five problems for each case and#Lemke is the total number ofpivotings in Lemke’s method.
Table 2 shows that the number of iterations of Newton method is consistently larger
than that of Algorithm 2 as far as the test problems used in the experiments are concerned.
Therefore, since it is time consuming to solve alinear subproblem at each iteration, Algorithm 2 required less CPU time than the pure Newton method in spite of the overheads
in line search. Finally we note that, the pure Newton method (18) is not guaranteed tobe globally convergent, although it actually converged for all test problems reported inTable 2.

Comparison of Algorithms la and 2

4.3

Finally we have compared Algorithms la and 2. For each of the problem sizes $n=$
$30,50$ and 90, we randomly generated five test problems. To see how these algorithms
behave for various degrees of asymunetry of the mapping , we have tested several valuesof between 0.1 and 2.0. The starting point was chosen to be $x=0$ . The results are
given in Table 3. All numbers shown in Table 3 are averages of the results for five test
problems for each case.
Table 3 shows that when the mapping is close to symmetry, Algorithm la converges
very fast, and when the mapping becomes asymmetric, the number of iterations and CPU
time of Algorithm la increase rapidly. On the other hand, in Algorithm 2, while the total
number of pivotings of Lemke’s method increases in proportion to problem size , thenumber of iterations stays constant even when problem size and the degree of asymmetryof are varied. Hence, when the degree of asymmetry of is relatively small, that iswhen is smaller than 1.0, Algorithm la requires less CPU time than Algorithm 2.
Note that, since the mapping used in our computational experience is sparse, complexity of each iteration in Algorithm la is small. On the other hand, the code [8] of
Lemke’s method used in Algorithm 2 to solve a hnear subproblem does not make use ofsparsity, so that it requires a significant amount of CPU time at each iteration for large
problems. If a method that can make use of sparsity is available to solve a linear subproblem, CPU time of Algorithm 2 may decrease. The projected Gauss-Seidel method[2, pp. 397] for solving linear complementarity problem is one of such methods. In Table
4, results of Algorithm 2 using the projected Gauss-Seidel method in place of Lemke’s
method are given. Table 4 shows that, if the mapping is almost symmetric, Algorithm
2 converges very fast. But Algorithm 2 fails to converge when the degree of asymmetry increased, because the projected Gauss-Seidel method could not to find a solution tolinear subproblem.
Figure 1 illustrates how Algorithms la and 2 converged for two typical test problems$F$$\rho$

$F$

$n$

$F$

$F$

$\rho$

$F$

$F$

163
with $n=30$ and 50. In the figure, the vertical axis represents the accuracy of a generated
point to the solution, which is evaluated by
ACC $= \max\{|\min(x_{i}, F_{1}(x))||i=1,2, \ldots, n\}$ .

(24)

Figure 1 indicates that Algorithm 2 is quadratically convergent near the solution. Figure
1 also indicates that Algorithm la is linearly convergent though it has not been proved
theoretically.

$\ovalbox{\tt\small REJECT}\#X\Re$[1] B. H. Ahn, “Solution of Nonsymmetric Linear Complementarity Problems by Iterative Methods,” Journal of optimization Theory and Applications 33 (1981) 175-185.
[2] R. W. Cottle, J. S. Pang and R. E. Stone, The Linear Complementarity Problem(Academic Press, San Diego, 1992)[3] S. Dafermos, (haffic Equilibrium and Variational Inequalities,” Transportation Science 14 (1980) 42-54.
[4] M. Florian, “Mathematical Programming Applications in National, Regional and
Urban Planning,” in: M. Iri and K. Tanabe eds., Mathematical Programming: Recent
Developments and Appli cations (KTK Scientific Publishers, Tokyo, 1989) pp.283-307.
[5] M. Fukushima, “Equivalent Differentiable optimization Problems and Descent Methods for Asymmetric Variational Inequality Problems,” Mathematical Programming 53(1992) 99-110.

[6] C. B. Garcia and W. I. Zangwill, Pathways to Solutions, Fixed Points, and Equilibria(Prentice-Hall, Eaglewood Cliffs N.J. 1981).
[7] P. T. Harker and J. S. Pang, (Finite-Dimensional Vaniational Inequality and Nonlinear Complementarity Problems: A Survey of Theory, Algorithms and Applications,”
Mathematical Programming 48 (1990) 161-220.
[8] T. Ibaraki and M. Fukushima, FORTRAN 77 optimization Programming (inJapanese) (Iwanami-Shoten, Tokyo, 1991).

[9] N. H. Josephy, “Newton’s Method for Generalized Equations,” Technical Report No.
1965, Mathematics Research Center, University of Wisconsin (Madison, WI, 1979).

164[10] C. E. Lemke, “Bimatrix Equilibrium Points and Mathematical Programming,” Management Science 11 (1965) 681-689.
[11] O. L. Mangasarian and M. V. Solodov, “Nonlinear Complementarity as Unconstrained and Constrained Minimization,” Technical Report No. 1074, Computer Sciences Department, University of Wisconsin (Madison, WI, 1992).
$\cdot$[12] L. Mathiesen, “An Algorithm based on a Sequence of Linear Complementarity Problems Applied to a Walrasian Equilibrium Model: An Example,” Mathematical Programming 37 (1987) 1-18.
[13] J. S. Pang and D. Chan, “Iterative Methods for Variational and Complementarity
Problems,” Mathematical Programming 24 (1982) 284-313.
[14] P. K. Subramanian, “Fixed Point Methods for The Complementarity Problem,”Technical Report No. 2857, Mathematics Research Center, University of Wisconsin (Madison, WI, 1985).
[15] K. Taji, M. Fukushima and T. Ibaraki, (A Globally Convergent Newton Method
for Solving Strongly Monotone Variational Inequalities,” to appear in Mathematical
Programming.
[16] R. L. Tobin, (A Variable Dimension Solution Approach for the General Spatial Price
Equilibrium Problem,” Mathematical Programming 40 (1988) 33-51.

165

Table l:Results for Algorithm la, lb and the projection method $(n=10, \rho=1)$

Table 2:Results for Algorithm 2 and Newton method $(\rho=1)$

166

167
ACC