∎

11institutetext:  Xin-long Luo, Corresponding author 22institutetext: School
of Artificial Intelligence, Beijing University of Posts and
Telecommunications, P. O. Box 101, Xitucheng Road No. 10, Haidian District,
100876, Beijing China, 22email<EMAIL_ADDRESS>33institutetext: Sen
Zhang 44institutetext: School of Artificial Intelligence,
Beijing University of Posts and Telecommunications, P. O. Box 101,
Xitucheng Road No. 10, Haidian District, 100876, Beijing China
44email<EMAIL_ADDRESS>55institutetext: Hang Xiao 66institutetext:
School of Artificial Intelligence,
Beijing University of Posts and Telecommunications, P. O. Box 101,
Xitucheng Road No. 10, Haidian District, 100876, Beijing China
66email<EMAIL_ADDRESS>
# Residual regularization path-following methods for linear complementarity
problems

Xin-long Luo $\ast$ Sen Zhang  Hang Xiao

(Received: date / Accepted: date)

###### Abstract

In this article, we consider the residual regularization path-following method
with the trust-region updating strategy for the linear complementarity
problem. This time-stepping selection based on the trust-region updating
strategy overcomes the shortcoming of the line search method, which consumes
the unnecessary trial steps in the transient-state phase. In order to improve
the robustness of the path-following method, we use the residual
regularization parameter to replace the traditional complementarity
regularization parameter. Moreover, we prove the global convergence of the new
method under the standard assumptions without the traditional assumption
condition of the priority to feasibility over complementarity. Numerical
results show that the new method is robust and efficient for the linear
complementarity problem, especially for the dense cases. And it is more robust
and faster than some state-of-the-art solvers such as the built-in subroutines
PATH and MILES of the GAMS v28.2 (2019) environment. The computational time of
the new method is about 1/3 to 1/10 of that of PATH for the dense linear
complementarity problem.

###### Keywords:

Continuation Newton method residual regularization trust-region updating
strategy complementarity path-following method

###### MSC:

90C33 65K05 65L05 65L20

††journal: Journal of XXXX

##  1 Introduction

In this article, we are mainly concerned with the linear complementarity
problem as follows:

$\displaystyle y=Mx+q,\;x_{i}y_{i}=0,\;i=1,\,2,\,\ldots,\,n,\;x\geq 0,\;y\geq
0,$ (1)

where $q\in\Re^{n}$ is a vector and $M$ is an $n\times n$ positive semi-
definite matrix. For the linear complementarity problem (1), there are many
practical applications such as the equilibrium of forces Erleben2007  and the
economic equilibrium problem Mathiesen1985 ; Rutherford1995 . And the
solutions of many problems such as the linear programming and the convex
quadratic programming can be obtained by solving it CPS2009 ; Wright1997 .
Furthermore, there are many efficient methods to solve it such as the Lemke’s
complementary pivoting algorithm (MILES) CPS2009 ; LH1964 ; Rutherford1995 ,
the path-following methods Fisher1992 ; Wright1994 ; Wright1997 ; Zhang1994
and their mixture method (PATH) Dirkse1994 ; DF1995 .

In this paper, we consider another path-following method based on the Newton
flow with nonnegative constraints, which is the variant of the primal-dual
path-following method. In order to improve its robustness and efficiency, we
use the residual regularization technique to avoid the singularity of the
Jacobian matrix, and adopt the trust-region updating strategy to adjust the
time step adaptively. Firstly, we construct the regularization Newton flow
with nonnegative constraints for the linear complementarity problem (1) based
on the primal-dual path-following method. Then, we use the implicit Euler
method and the linear approximation of the quadratic function to obtain the
regularization path-following method for following the the trajectory of the
Newton flow. Finally, we adopt the trust-region updating strategy to adjust
the time step adaptively. The advantage of this time-stepping selection
compared with the line search method is that it overcomes the shortcoming of
the line search method, which consumes the unnecessary trial steps in the
transient-state phase. The other advantage is that it improves the robustness
of the path-following method and it does not require the traditional
assumption condition of the priority to feasibility over complementarity of
the path-following method Wright1994 ; Zhang1994  when we prove its global
convergence.

The rest of this article is organized as follows. In the next section, we
consider the regularization path-following method and the adaptive trust-
region updating strategy for the linear complementarity problem. In section 3,
we prove the global convergence of the new method under the standard
assumptions without the traditional assumption condition of the priority to
feasibility over complementarity. In section 4, we compare the new method with
two state-of-the-art solvers, i.e. PATH Dirkse1994 ; DF1995 ; PATH  and MILES
(a Mixed Inequality and nonLinear Equation Solver) Mathiesen1985 ;
Rutherford1995 ; Rutherford2022  for sparse problems and dense problems, test
matrices of which come from the linear programming subset of NETLIB NETLIB .
The new method is coded with the MATLAB language and executed in MATLAB
(R2020a) environment MATLAB . PATH and MILES are executed in the GAMS v28.2
(2019) environment GAMS . Numerical results show that the new method is robust
and efficient for solving the linear complementarity problem. It is more
robust and faster than PATH and MILES for the dense problems. Finally, some
discussions are given in section 5. $\|\cdot\|$ denotes the Euclidean vector
norm or its induced matrix norm throughout the paper.

##  2 Regulation path-following methods

###  2.1 The continuous Newton flow

For convenience, we rewrite the linear complementarity problem (1) as the
following nonlinear system of equations with nonnegativity constraints:

$\displaystyle F(z)=\begin{bmatrix}y-(Mx+q)\\\
XYe\end{bmatrix}=0,\;(x,\,y)\geq 0,\;z=(x,\,y),$ (2)

where $X=diag(x),\;Y=diag(y)$ and all components of vector $e$ equal one. It
is not difficult to know that the Jacobian matrix $J(z)$ of $F(z)$ has the
following form:

$\displaystyle J(z)=\begin{bmatrix}-M&I\\\ Y&X\end{bmatrix}.$ (3)

From the second block $XYe=0$ of equation (2), we know that $x_{i}=0$ or
$y_{i}=0\,(i=1:n)$. Thus, the Jacobian matrix $J(z)$ of equation (3) may be
singular, which leads to numerical difficulties near the solution of the
nonlinear system (2) for the Newton method or its variants. In order to
overcome this difficulty, we consider its perturbed system AG2003 ; Tanabe1988
as follows:

$\displaystyle F_{\mu}(z)=F(z)-\begin{bmatrix}0\\\ \mu
e\end{bmatrix}=0,\;z=(x,\,y)>0,\;\mu>0.$ (4)

It is not difficult to verify that the Jacobian matrix $J(z)$ defined by
equation (3) is nonsingular when $M$ is a positive semi-definite matrix and
$(x,\;y)>0$ (see Lemma 5.9.8, p. 469 in CPS2009 ). Thus, by using the implicit
theorem and the inequality $(a-b)(c-d)\leq|ac-bd|\;(a>0,\,b>0,\,c>0,\,d>0)$,
the perturbed system (4) has a unique solution when $M$ is a positive definite
matrix (see Theorem 5.9.13, p. 471 in CPS2009 ) for its detailed proof). The
solution $z(\mu)$ of the perturbed system (4) defines the central path, and
$z(\mu)$ approximates the solution $z^{\ast}$ of the nonlinear system (2) when
$\mu$ tends to zero Wright1997 ; CPS2009 .

If the damped Newton method is applied to the perturbation system (4)
Kelley2003 ; NW1999 , we have

$\displaystyle z_{k+1}=z_{k}-\alpha_{k}J(z_{k})^{-1}F_{\mu}(z_{k}),$ (5)

where $J(z_{k})$ is the Jacobian matrix of $F_{\mu}(z)$. We regard
$z_{k+1}=z(t_{k}+\alpha_{k}),\;z_{k}=z(t_{k})$ and let $\alpha_{k}\to 0$, then
we obtain the continuous Newton flow with nonnegativity constraints AS2015 ;
Davidenko1953 ; LY2021 ; PHP1975 ; Tanabe1988  of the perturbed system (4) as
follows:

$\displaystyle\frac{dz(t)}{dt}=-J(z)^{-1}F_{\mu}(z),\hskip
5.69054ptz=(x,\,y)>0.$ (6)

Actually, if we apply an iteration with the explicit Euler method SGT2003  for
the continuous Newton flow (6), we obtain the damped Newton method (5).

Since the Jacobian matrix $J(z)=F^{\prime}_{\mu}(z)$ may be singular, we
reformulate the continuous Newton flow (6) as the following general formula
for the linear complementarity problem (2) :

$\displaystyle-M\frac{dx(t)}{dt}+\frac{dy(t)}{dt}=-r_{q}(x,\,y),$ (7)
$\displaystyle
Y\frac{dx(t)}{dt}+X\frac{dy(t)}{dt}=-(XYe-\mu(t)e),\;(x,\,y)>0,$ (8)

where the residual $r_{q}(x,\,y)=y-(Mx+q)$. The continuous Newton flow (7)-(8)
has some nice properties. We state one of them as the following property 1
Branin1972 ; LXL2022 ; LY2021 ; LX2021 ; Tanabe1979 .

###### Property 1

Assume that $(x(t),\,y(t))$ is the solution of the continuous Newton flow
(7)-(8), then $r_{q}(x(t),\,y(t))$ converges to zero and
$x_{i}(t)y_{i}(t)\,(i=1,\,2,\,\ldots,\,n)$ converge to zero when
$0\leq\mu(t)\leq\sigma\min_{1\leq i\leq
n}\\{x_{i}(t)y_{i}(t)\\}\;(0<\sigma<1)$ and $t\to\infty$. That is to say, for
every limit point $(x^{\ast},\,y^{\ast})$ of $(x(t),\,y(t))$, it is also a
solution of the linear complementarity problem (2). Furthermore, $x(t)$ and
$y(t)$ keep positive values when the initial point
$\left(x^{0}_{i},\;y^{0}_{i}\right)>0\,(i=1,\,2,\,\ldots,\,n)$.

Proof. Assume that $z(t)$ is the continuous solution of the continuous Newton
flow (7)-(8), then we have

$\displaystyle\frac{d}{dt}r_{q}(x,\,y)=-M\frac{dx}{dt}+\frac{dy}{dt}=-r_{q}(x,\,y),$
(9)
$\displaystyle\frac{d}{dt}(XYe)=X\frac{dy}{dt}+Y\frac{dx}{dt}=-(XYe-\mu(t)e).$
(10)

Consequently, from equations (9)-(10) and $0\leq\mu(t)\leq\sigma\min_{1\leq
i\leq n}\\{x_{i}(t)y_{i}(t)\\}$, we obtain

$\displaystyle r_{q}(x(t),\,y(t))=r_{q}\left(x^{0},\,y^{0}\right)\exp(-t),$
(11) $\displaystyle-XYe\leq\frac{d}{dt}(XYe)\leq-(1-\sigma)XYe.$ (12)

From equations (11)-(12), we know that $r_{q}(x(t),\,y(t))$ converges to zero
with the linear rate of convergence when $t$ tends to infinity. Furthermore,
from equation (12) and the Gronwall inequality Gronwall1919 , we have

$\displaystyle x_{i}^{0}y_{i}^{0}\exp(-t)\leq x_{i}(t)y_{i}(t)\leq
x_{i}^{0}y_{i}^{0}\exp(-(1-\sigma)t),\;i=1,\,2,\,\ldots,\,n.$ (13)

Consequently, from equation (13), we know that $x_{i}(t)y_{i}(t)\geq
0\;(i=1,\,2,\,\ldots,\,n)$ and
$\lim_{t\to\infty}x_{i}(t)y_{i}(t)=0\;(i=1,\,2,\,\ldots,\,n)$.

From equation (13), we know that $x_{i}(t)y_{i}(t)>0\;(t\geq 0)$ when
$\left(x^{0},\;y^{0}\right)>0$. If we have $x_{i}(T)<0$ or $y_{i}(T)<0$ for a
finite value $T>0$, there exists $\bar{t}\in(0,\;T)$ such that
$x_{i}(\bar{t})=0$ or $y_{i}(\bar{t})=0$. Consequently, we have
$x_{i}(\bar{t})y_{i}(\bar{t})=0$, which contradicts
$x_{i}(\bar{t})y_{i}(\bar{t})>0$. Thus, we have $(x(t),\,y(t))>0$ for all
$t>0$. Therefore, if the solution $(x(t),\,y(t))$ of the continuous Newton
flow (7)-(8) belongs to a compact set, there exists a limit point
$(x^{\ast},\,y^{\ast})$ when $t$ tends to infinity, and this limit point
$(x^{\ast},\,y^{\ast})$ is also a solution of the linear complementarity
problem (2). ∎

###### Remark 1

The inverse $J(z)^{-1}$ of the Jacobian matrix $J(z)$ can be regarded as the
pre-conditioner of $F_{\mu}(z)$ such that the every element $z_{i}(t)$ of
$z(t)$ has roughly the same rate of convergence and it mitigates the stiffness
of the ODE (6) LXL2022 ; LX2021 . This property is very useful since it makes
us adopt the explicit ODE method to follow the trajectory of the Newton flow
(6) efficiently.

###  2.2 The residual regularization path-following method

From property 1, we know that the continuous Newton flow (7)-(8) has the
global convergence. However, when the Jacobian matrix $J(z)$ is singular or
nearly singular, the ODE (7)-(8) is the system of differential-algebraic
equations AP1998 ; BCP1996 ; HW1996  and its trajectory can not be efficiently
solved by the general ODE method BJ1998 ; JT1995  such as the backward
differentiation formulas (the built-in subroutine ode15s.m of the MATLAB
environment MATLAB ; SGT2003 ). Thus, we need to construct the special method
to solve this problem. Furthermore, we expect that the new method has the
global convergence as the homotopy continuation methods AG2003 ; OR2000  and
the fast rate of convergence as the traditional optimization methods. In order
to achieve these two aims, we consider the continuation Newton method and the
trust-region updating strategy for problem (7)-(8).

We apply the implicit Euler method AP1998 ; BCP1996  to the continuous Newton
flow (7)-(8), then we obtain

$\displaystyle-M\frac{x^{k+1}-x^{k}}{\Delta t_{k}}+\frac{y^{k+1}-y^{k}}{\Delta
t_{k}}=-r_{q}\left(x^{k+1},\,y^{k+1}\right),$ (14) $\displaystyle
Y^{k+1}\frac{x^{k+1}-x^{k}}{\Delta t_{k}}+X^{k+1}\frac{y^{k+1}-y^{k}}{\Delta
t_{k}}=-\left(X^{k+1}Y^{k+1}e-\mu(t_{k+1})e\right).$ (15)

Since equation (15) is a nonlinear system, it is not directly solved, and we
seek for its explicit approximation formula. We replace $Y^{k+1}$ and
$X^{k+1}$ with $Y^{k}$ and $X^{k}$ in the left-hand side of equation (15),
respectively. And we substitute $X^{k+1}Y^{k+1}e$ with its linear
approximation $X^{k}Y^{k}e+\frac{\Delta t_{k}}{1+\Delta t_{k}}(Y^{k}\Delta
x^{k}+X^{k}\Delta y^{k})$ in the right-hand side of equation (15). We set
$\mu(t_{k+1})=\sigma_{k}\mu_{k}$. Then, we obtain the continuation Newton
method (one of path-following methods) as follows:

$\displaystyle-M\Delta x^{k}+\Delta y^{k}=-r_{q}^{k},$ (16) $\displaystyle
Y^{k}\Delta x^{k}+X^{k}\Delta y^{k}=-r_{c}^{k},$ (17) $\displaystyle
x^{k+1}=x^{k}+\frac{\Delta t_{k}}{1+\Delta t_{k}}\Delta
x^{k},\;y^{k+1}=y^{k}+\frac{\Delta t_{k}}{1+\Delta t_{k}}\Delta y^{k},$ (18)

where $r_{q}^{k}=r_{q}(x^{k},\,y^{k})$,
$r_{c}^{k}=X^{k}y^{k}-\sigma_{k}\mu_{k}e$,
$0<\sigma_{min}\leq\sigma_{k}\leq\sigma_{max}<1$ and the perturbation
parameter $\mu_{k}$ is selected as follows:

$\displaystyle\mu_{k}=\frac{(x^{k})^{T}y^{k}+\|r_{q}^{k}\|}{2n}.$ (19)

The selection of $\mu_{k}$ in equation (19) is slightly different to the
traditional selection $\mu_{k}=\left(x^{k}\right)^{T}y^{k}/n$ Wright1994 ;
Wright1997 ; Zhang1994 . According to our numerical experiments, the selection
of $\mu_{k}$ in equation (19) can improve the robustness of the path-following
method, in comparison to the traditional selection selection
$\mu_{k}=\left(x^{k}\right)^{T}y^{k}/n$.

###### Remark 2

If we set $\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$ in equation (18), we
obtain the damped Newton method (the primal-dual path-following method) (5).
However, from the view of the ODE method, they are different. The damped
Newton method (5) is derived from the explicit Euler method applied to the
continuous Newton flow (7)-(8). Its time step $\alpha_{k}$ is restricted by
the numerical stability HW1996 ; SGT2003 . That is to say, for the linear test
equation $dx/dt=-\lambda x\;(\lambda>0)$, its time step $\alpha_{k}$ is
restricted by the stable region $|1-\lambda\alpha_{k}|\leq 1$. Therefore, the
large step $\alpha_{k}$ can not be adopted in the steady-state phase. The
continuation Newton method (16)-(18) is derived from the implicit Euler method
applied to the continuous Newton flow (7)-(8) and the linear approximation of
$F_{\mu(t_{k+1})}(z_{k+1})$, and its time step $\Delta t_{k}$ is not
restricted by the numerical stability for the linear test equation. Therefore,
the large time step $\Delta t_{k}$ can be adopted in the steady-state phase,
and the continuation Newton method (16)-(18) mimics the Newton method.
Consequently, it has the fast rate of convergence near the solution $z^{\ast}$
of the nonlinear system (2). The most of all, the substitution $\alpha_{k}$
with $\Delta t_{k}/(\Delta t_{k}+1)$ is favourable to adopt the trust-region
updating strategy for adaptively adjusting the time step $\Delta t_{k}$ such
that the continuation Newton method (16)-(18) accurately follows the
trajectory of the continuation Newton flow (7)-(8) in the transient-state
phase and achieves the fast rate of convergence in the steady-state phase.

When the diagonal matrix $X^{k}$ is positive definite, from equations
(16)-(17), $\Delta x^{k}$ and $\Delta y^{k}$ can be solved by the following
two subsystems:

$\displaystyle\left(M+(X^{k})^{-1}Y^{k}\right)\Delta
x^{k}=r_{q}^{k}-(X^{k})^{-1}r_{c}^{k},$ (20) $\displaystyle\Delta
y^{k}=M\Delta x^{k}-r_{q}^{k}.$ (21)

When matrix $M$ is positive semi-definite and $(x^{k},\,y^{k})>0$, the left
hand-side matrix is positive definite. Thus, the system (20) can be solved by
the partial pivoting Gaussian elimination method (see pp. 125-130, GV2013 ).

###  2.3 The time-stepping control and the initial point selection

Another issue is how to adaptively adjust the time step $\Delta t_{k}$ at
every iteration. A popular and efficient time-stepping control is based on the
trust-region updating strategy CGT2000 ; Deuflhard2004 ; Higham1999 ; Luo2010
; Luo2012 ; LXLZ2021 ; LX2021 ; LX2021B ; LX2022 ; LY2021 ; LLS2022 ; LXL2022
; Yuan2015 . Its main idea can be described as follows. Firstly, we construct
a merit function reflecting the feasibility $r_{q}(x,\,y)=0$ and the
complementarity $x_{i}y_{i}=0\;(i=1,\,2,\,\ldots,\,n)$ as

$\displaystyle\phi(x,\,y)=x^{T}y+\|r_{q}(x,\,y)\|,$ (22)

where $r_{q}(x,\,y)=y-(Mx+q)$ and $(x,\,y)\geq 0$.

Then, we define the linear approximation model $m_{k}$ of
$\phi(x^{k}+\alpha_{k}\Delta x^{k},\,y^{k}+\alpha_{k}\Delta y^{k})$ as

$\displaystyle m_{k}(\alpha)$
$\displaystyle=\phi(x^{k},\,y^{k})+\alpha\left((y^{k})^{T}\Delta
x^{k}+(x^{k})^{T}\Delta y^{k}-\|r_{q}^{k}\|\right),$ (23)

where $\alpha=\Delta t/(1+\Delta t)$ and $(\Delta x^{k},\,\Delta y^{k})$ is
the Newton step and solved by equations (16)-(17).

Finally, we adaptively adjust the time step $\Delta t_{k}$ according to the
difference between $m_{k}(\alpha_{k})$ and $\phi(x^{k}+\alpha_{k}\Delta
x^{k},\,y^{k}+\alpha_{k}\Delta y^{k})$. Namely, the time step $\Delta t_{k+1}$
will be enlarged when $m_{k}(\alpha_{k})$ approximates
$\phi(x^{k}+\alpha_{k}\Delta x^{k},\,y^{k}+\alpha_{k}\Delta y^{k})$ well, and
$\Delta t_{k+1}$ will be reduced when $m_{k}(\alpha_{k})$ approximates
$\phi(x^{k}+\alpha_{k}\Delta x^{k},\,y^{k}+\alpha_{k}\Delta y^{k})$ badly.

We define the predicted reduction $Pred_{k}$ and the actual reduction
$Ared_{k}$ as follows:

$\displaystyle
Pred_{k}=\phi(x^{k},\,y^{k})-m_{k}(\alpha_{k})=\alpha_{k}\left(\|r_{q}^{k}\|-(y^{k})^{T}\Delta
x^{k}-(x^{k})^{T}\Delta y^{k}\right),$ (24) $\displaystyle
Ared_{k}=\phi(x^{k},\,y^{k})-\phi(x^{k}+\alpha_{k}\Delta
x^{k},\,y^{k}+\alpha_{k}\Delta y^{k})$ $\displaystyle\hskip
5.69054pt=\alpha_{k}\left(\|r_{q}^{k}\|-(y^{k})^{T}\Delta
x^{k}-(x^{k})^{T}\Delta y^{k}\right)-\alpha_{k}^{2}(\Delta x^{k})^{T}\Delta
y^{k}.$ (25)

Then, we enlarge or reduce the time step $\Delta t_{k+1}$ at every iteration
according to the following ratio:

$\displaystyle\rho_{k}=\frac{Ared_{k}}{Pred_{k}}=\frac{\|r_{q}^{k}\|-(y^{k})^{T}\Delta
x^{k}-(x^{k})^{T}\Delta y^{k}-\alpha_{k}(\Delta x^{k})^{T}\Delta
y^{k}}{\|r_{q}^{k}\|-(y^{k})^{T}\Delta x^{k}-(x^{k})^{T}\Delta y^{k}},$ (26)

where $\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$. A particular adjustment
strategy is given as follows:

$\displaystyle\Delta t_{k+1}=\begin{cases}2\Delta
t_{k},\;\text{if}\;\rho_{k}\geq\eta_{2}\;\text{and}\;(x^{k+1},\,y^{k+1})>0,\\\
\Delta t_{k},\;\text{else
if}\;\eta_{1}\leq\rho_{k}<\eta_{2}\;\text{and}\;(x^{k+1},\,y^{k+1})>0,\\\
\frac{1}{2}\Delta t_{k},\;\text{others},\end{cases}$ (27)

where the constants $\eta_{1},\;\eta_{2}$ are selected as
$\eta_{1}=0.25,\;\eta_{2}=0.75$, respectively. We set

$\displaystyle(x^{k+1},\,y^{k+1})=(x^{k},\,y^{k})+\frac{\Delta t_{k}}{1+\Delta
t_{k}}(\Delta x^{k},\,\Delta y^{k}).$ (28)

When $\rho_{k}\geq\eta_{a}$ and $(x^{k+1},\,y^{k+1})>0$, we accept the trial
step, otherwise we discard the trial step and set

$\displaystyle(x^{k+1},\,y^{k+1})=(x^{k},\,y^{k}),$ (29)

where $\eta_{a}$ is a small positive number such as $\eta_{a}=1.0\times
10^{-6}$.

###### Remark 3

This new time-stepping control based on the trust-region updating strategy has
some advantages compared to the traditional line search strategy. If we use
the line search strategy and the damped Newton method (5) to follow the
trajectory $z(t)$ of the continuous Newton flow (6), in order to achieve the
fast rate of convergence in the steady-state phase, the time step $\alpha_{k}$
of the damped Newton method is tried from 1 and reduced by the half with many
times at every iteration. Since the linear model
$F_{\sigma_{k}\mu_{k}}(z_{k})+J(z_{k})\Delta z_{k}$ may not approximate
$F_{\sigma_{k}\mu_{k}}(z_{k}+\Delta z_{k})$ well in the transient-state phase,
the time step $\alpha_{k}$ will be small. Consequently, the line search
strategy consumes the unnecessary trial steps in the transient-state phase.
However, the selection of the time step $\Delta t_{k}$ based on the trust-
region updating strategy (26)-(27) can overcome this shortcoming.

In order to ensure that the algorithm works well for the general linear
complementarity problem, the initial point selection is also a key point. We
select the starting point $(x^{0},\,y^{0})$ as follows:

$\displaystyle
x^{0}=10*e,\;v^{0}=Mx^{0}+q,\;y^{0}_{i}=\begin{cases}v_{i}^{0},\;\text{if}\;v_{i}^{0}>0,\\\
10^{-3},\;\text{otherwise}.\end{cases}$ (30)

In order to improve the stability of the algorithm, we add a small
regularization item $\upsilon I$ to matrix $M$ CPS2009 ; Gana1982 ; Hansen1994
; Venkateswaran1993  when $\mu_{k}$ is not small, where $\mu_{k}$ is defined
by equation (19). Specifically, we adopt the following strategy as the
regularizer $M_{\upsilon}$ of matrix $M$:

$\displaystyle M_{\upsilon}=\begin{cases}M+\upsilon
I,\;\text{if}\;\mu_{k}\geq\upsilon,\\\ M,\;\text{otherwise},\end{cases}$ (31)

where $\upsilon=10^{-3}$.

According to the above discussions, we give the detailed descriptions of the
path-following method and the trust-region updating strategy for the monotone
linear complementarity problem (1) in Algorithm 1.

Algorithm 1 Regularization path-following methods with the trust-region
updating strategy for linear complementarity problems (The RPFMTr method)

0: matrix $M\in\Re^{n\times n}$ and vector $q\in\Re^{n}$ for the problem:
$y=Mx+q,\;x_{i}y_{i}=0\;(i=1:n),\;x\geq 0,\;y\geq 0$.

0: the linear complementarity solution: $(solx,\,soly)$.

1: Initialize parameters: $\eta_{a}=10^{-6}$, $\eta_{1}=0.25$,
$\eta_{2}=0.75$, $\epsilon=10^{-6}$, $\Delta t_{0}=10^{-2}$, bigMfac = 10,
$\upsilon=10^{-3}$, $\sigma_{0}=0.5$, maxit = 600.

2: Initialize $x^{0}$ = bigMfac*ones(n, 1); $y^{0}=Mx^{0}+q$;
$y^{0}(y^{0}<0)=10^{-3}$.

3: Regularize matrix: $M_{\upsilon}=M+\upsilon I$.

4: Set flag_success_trialstep = 1, $\text{itc}=0,\;k=0$.

5: while (itc $<$ maxit) do

6: if (flag_success_trialstep == 1) then

7: Set itc = itc + 1.

8: Compute $r_{q}^{k}=y^{k}-(M_{\upsilon}x^{k}+q)$;
$\mu_{k}=(\|r_{q}^{k}\|+(x^{k})^{T}y^{k})/(2n)$.

9: Compute
$\text{Resk}=\max\\{\|x^{k}.y^{k}\|_{\infty},\;\|y^{k}-(Mx^{k}+q)\|_{\infty}\\}$.

10: if ($\text{Resk}<\epsilon$) then

11: break;

12: end if

13: Set $\sigma_{k}=\min\\{\sigma_{k},\;\mu_{k}\\}$.

14: Compute $r_{c}^{k}=x^{k}.y^{k}-\sigma_{k}\mu_{k}*\text{ones}(n,1)$.

15: By solving the linear system (20)-(21), we obtain $\Delta x^{k}$ and
$\Delta y^{k}$.

16: end if

17: Set $(x^{k+1},\,y^{k+1})=(x^{k},\,y^{k})+\frac{\Delta t_{k}}{1+\Delta
t_{k}}(\Delta x^{k},\,\Delta y^{k})$.

18: Compute the ratio $\rho_{k}$ from equation (26) and adjust $\Delta
t_{k+1}$ according to the formula (27).

19: if (($\rho_{k}\geq\eta_{a}$) && ($x^{k+1},\,y^{k+1})>0$)) then

20: Accept the trial point $(x^{k+1},\,y^{k+1})$; Set flag_success_trialstep =
1.

21: if ($\|x^{k+1}-x^{k}\|_{\infty}>0.1$) then

22: Set $\sigma_{k+1}=0.5$.

23: else

24: Set $\sigma_{k+1}=0.1$.

25: end if

26: if ($\mu_{k}<\upsilon$) then

27: Set $M_{\upsilon}=M$.

28: end if

29: else

30: Set $(x^{k+1},\,y^{k+1})=(x^{k},\,y^{k})$; flag_success_trialstep = 0.

31: end if

32: Set $k\leftarrow k+1$.

33: end while

34: Return $(solx,\,soly)=(x^{k},\,y^{k})$.

##  3 Convergence analysis

In this section, we analyze the global convergence of Algorithm 1. Without
loss of generality, we assume $M_{\upsilon}=M$ in the following analysis.
Namely, we do not discriminate between $M_{\upsilon}$ and $M$. Similarly to
the analysis techniques of the references Wright1994 ; Xu1991 ; Zhang1994 , we
firstly construct an auxiliary sequence $(u^{k},\,v^{k})$ as follows:

$\displaystyle u^{k+1}=u^{k}+\alpha_{k}(x^{k}+\Delta x^{k}-u^{k}),$ (32)
$\displaystyle v^{k+1}=v^{k}+\alpha_{k}(y^{k}+\Delta y^{k}-v^{k}),$ (33)

where $(\Delta x^{k},\,\Delta y^{k})$ are solved by equations (16)-(17) and
$(u^{0},\,v^{0})\leq(x^{0},\,y^{0})$ satisfies the feasibility
$r_{q}(u^{0},\,v^{0})=0$. Then, it is not difficult to verify

$\displaystyle r_{q}(u^{k},\,v^{k})=v^{k}-(Mu^{k}+q)=0,$ (34) $\displaystyle
x^{k+1}-u^{k+1}=(1-\alpha_{k})(x^{k}-u^{k})\geq 0,$ (35) $\displaystyle
y^{k+1}-v^{k+1}=(1-\alpha_{k})(y^{k}-v^{k})\geq 0,$ (36)

where $\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$.

Meanwhile, in order to obtain the global convergence, we need to enforce the
condition specified below. Firstly, we select a constant $\gamma\in(0,\,1)$
that satisfies

$\displaystyle 0<\gamma\leq\frac{\mu_{0}}{\min_{1\leq i\leq
n}\\{x^{0}_{i}y^{0}_{i}\\}}=\frac{(x^{0})^{T}y^{0}+\|r_{q}^{0}\|}{2n\,\min_{1\leq
i\leq n}\\{x^{0}_{i}y^{0}_{i}\\}}.$ (37)

Then, we select $\alpha_{k}$ such that

$\displaystyle
x^{k}_{i}(\alpha)y^{k}_{i}(\alpha)\geq\gamma\mu_{k}(\alpha),\;i=1,\,2,\,\ldots,n$
(38)

holds for all $\alpha\in(0,\;\alpha_{k}]\subset(0,\,1]$, where
$x^{k}(\alpha)$, $y^{k}(\alpha)$, $r_{q}^{k}(\alpha)$ and $\mu_{k}(\alpha)$
are defined by

$\displaystyle x^{k}(\alpha)=x^{k}+\alpha\Delta
x^{k},\;y^{k}(\alpha)=y^{k}+\alpha\Delta y^{k},$ (39) $\displaystyle
r_{q}^{k}(\alpha)=r_{q}(x^{k}(\alpha),y^{k}(\alpha)),\;\mu_{k}(\alpha)=\frac{(x^{k}(\alpha))^{T}y^{k}(\alpha)+\|r_{q}^{k}(\alpha)\|}{2n}.$
(40)

Condition (38) is to prevent iterations from prematurely getting too close to
the boundary of the positive quadrant and its restriction on $\gamma$ is very
mild. In the practice, it can be selected to be very small.

In order to establish the main global convergence of Algorithm 1, similarly to
the results Wright1994 ; Zhang1994 , we prove the following several technique
lemmas.

###### Lemma 1

When $(x^{k}(\alpha),\,y^{k}(\alpha))$ satisfies the condition (38), where
$(\Delta x^{k},\,\Delta y^{k})$ is the solution of equation (16)-(17), we have

$\displaystyle(x^{k}(\alpha),\,y^{k}(\alpha))\geq 0.$ (41)

Furthermore, when $x^{k}_{i}y^{k}_{i}\geq\gamma\mu_{k}\;(i=1,\,2,\ldots,\,n)$
and $\alpha$ satisfies

$\displaystyle
0\leq\alpha\leq\min\left\\{1,\;\frac{1-\gamma/2}{1+\gamma/2}\frac{\sigma_{k}\mu_{k}}{\left\|\Delta
X^{k}\Delta y^{k}\right\|_{\infty}}\right\\},$ (42)

the inequality (38) holds.

Proof. By summing two sides of the condition (38), we obtain

$\displaystyle(x^{k}(\alpha))^{T}y^{k}(\alpha)\geq
n\gamma\mu_{k}(\alpha)=\frac{1}{2}\gamma\left((x^{k}(\alpha))^{T}y^{k}(\alpha)+(1-\alpha)\|r_{q}^{k}\|\right)$
$\displaystyle\hskip
5.69054pt\geq\frac{1}{2}\gamma(x^{k}(\alpha))^{T}y^{k}(\alpha).$ (43)

Consequently, from equation (43) and $0<\gamma<1$, we obtain
$(x^{k}(\alpha))^{T}y^{k}(\alpha)\geq 0$. By substituting it into the
condition (38), we have $x^{k}_{i}(\alpha)y^{k}_{i}(\alpha)\geq
0\,(i=1,\,2,\,\ldots,\,n)$. From equation (17), we obtain

$\displaystyle\alpha\sigma_{k}\mu_{k}=\alpha x^{k}_{i}y^{k}_{i}+\alpha
x^{k}_{i}\Delta y^{k}_{i}+\alpha y^{k}_{i}\Delta x^{k}_{i}$
$\displaystyle\hskip
5.69054pt=(\alpha-2)x^{k}_{i}y^{k}_{i}+y^{k}_{i}x^{k}_{i}(\alpha)+x^{k}_{i}y^{k}_{i}(\alpha),\;i=1,\,2,\ldots,\,n.$
(44)

If $(x^{k}_{i}(\alpha),\,y^{k}_{i}(\alpha))<0$ or
$x^{k}_{i}(\alpha)=0,\,y^{k}_{i}(\alpha)<0$ or
$x^{k}_{i}(\alpha)<0,\,y^{k}_{i}(\alpha)=0$, by combining it with
$(x^{k}_{i},\,y^{k}_{i})\geq 0$, it contradicts equation (44). Therefore, we
obtain $(x^{k}_{i}(\alpha),\,y^{k}_{i}(\alpha))\geq
0\,(i=1,\,2,\,\ldots,\,n)$, which proves equation (41).

From equations (16)-(17), we have

$\displaystyle
x^{k}_{i}(\alpha)y^{k}_{i}(\alpha)-\gamma\mu_{k}(\alpha)=x_{i}^{k}y_{i}^{k}+\alpha(x^{k}_{i}\Delta
y^{k}_{i}+y^{k}_{i}\Delta x^{k}_{i})+\alpha^{2}\Delta x^{k}_{i}\Delta
y^{k}_{i}$ $\displaystyle\hskip
5.69054pt-\gamma\mu_{k}-\frac{1}{2n}\gamma\alpha\left((x^{k})^{T}\Delta
y^{k}+(y^{k})^{T}\Delta x^{k}-\|r_{q}^{k}\|+\alpha(\Delta x^{k})^{T}\Delta
y^{k}\right)$
$\displaystyle=x_{i}^{k}y_{i}^{k}+\alpha(\sigma_{k}\mu_{k}-x^{k}_{i}y^{k}_{i})+\alpha^{2}\Delta
x^{k}_{i}\Delta y^{k}_{i}$ $\displaystyle\hskip
5.69054pt-\gamma\mu_{k}-\frac{1}{2n}\gamma\alpha\left(n\sigma_{k}\mu_{k}-(x^{k})^{T}y^{k}-\|r_{q}^{k}\|+\alpha(\Delta
x^{k})^{T}\Delta y^{k}\right)$
$\displaystyle=(1-\alpha)\left(x^{k}_{i}y^{k}_{i}-\gamma\mu_{k}\right)+\alpha\left(\sigma_{k}\mu_{k}\left(1-\frac{\gamma}{2}\right)+\alpha\left(\Delta
x^{k}_{i}\Delta y^{k}_{i}-\frac{\gamma}{2n}(\Delta x^{k})^{T}\Delta
y^{k}\right)\right)$
$\displaystyle\geq(1-\alpha)\left(x^{k}_{i}y^{k}_{i}-\gamma\mu_{k}\right)+\alpha\left(\sigma_{k}\mu_{k}\left(1-\frac{\gamma}{2}\right)-\alpha\left(1+\frac{\gamma}{2}\right)\left\|\Delta
X^{k}\Delta y^{k}\right\|_{\infty}\right),$ (45)

where $\Delta X^{k}=\text{diag}(\Delta x^{k})$. By substituting
$x^{k}_{i}y^{k}_{i}\geq\gamma\mu_{k}$ and $0<\alpha\leq 1$ into equation (45),
we have

$\displaystyle
x^{k}_{i}(\alpha)y^{k}_{i}(\alpha)-\gamma\mu_{k}(\alpha)\geq\alpha\left(\sigma_{k}\mu_{k}\left(1-\frac{\gamma}{2}\right)-\alpha\left(1+\frac{\gamma}{2}\right)\left\|\Delta
X^{k}\Delta y^{k}\right\|_{\infty}\right).$ (46)

Therefore, when $\alpha$ satisfies equation (42), from equation (46), we
obtain

$\displaystyle x^{k}_{i}(\alpha)y^{k}_{i}(\alpha)-\gamma\mu_{k}(\alpha)\geq
0,\;i=1,\,2,\,\ldots,\,n,$

which gives the inequality (38). $\quad\square$

###### Lemma 2

Assume that $\\{(x^{k},\,y^{k})\\}$ is generated by Algorithm 1. Then, for any
feasible solution $(\bar{x},\,\bar{y})$ of the linear complementarity problem
(2), there exists a positive constant $C_{1}$ such that

$\displaystyle(y^{k})^{T}(x^{k}-u^{k})+(x^{k})^{T}(y^{k}-v^{k})\leq C_{1}$
(47)

holds for all $k=0,\,1,\,\ldots$. Furthermore, if $(\bar{x},\,\bar{y})$ is
strictly feasible, $\\{(x^{k},\,y^{k})\\}$ is bounded.

Proof. From equation (34), we know that $(u^{k},\,v^{k})$ satisfies the
feasibility. By combining it with the semi-definite positivity of $M$, we have

$\displaystyle(\bar{x}-u^{k})^{T}(\bar{y}-v^{k})=(\bar{x}-u^{k})^{T}M(\bar{x}-u^{k})\geq
0.$ (48)

Consequently, from equation (48), we have

$\displaystyle 0$ $\displaystyle\leq(\bar{x}-u^{k})^{T}(\bar{y}-v^{k})$
$\displaystyle=(\bar{x}-x^{k}+(x^{k}-u^{k}))^{T}(\bar{y}-y^{k}+(y^{k}-v^{k}))$
$\displaystyle=\bar{x}^{T}(y^{k}-v^{k})-(x^{k})^{T}(y^{k}-v^{k})+\bar{x}^{T}\bar{y}+(x^{k})^{T}y^{k}-\bar{x}^{T}y^{k}-\bar{y}^{T}x^{k}$
$\displaystyle\hskip
5.69054pt+(x^{k}-u^{k})^{T}(y^{k}-v^{k})+\bar{y}^{T}(x^{k}-u^{k})-(y^{k})^{T}(x^{k}-u^{k}).$
(49)

From equations (35)-(36), we have

$\displaystyle 0\leq x^{k}-u^{k}\leq x^{0}-u^{0},\;0\leq y^{k}-v^{k}\leq
y^{0}-v^{0}.$ (50)

By substituting equation (50) into equation (49), from $(x^{k},\,y^{k})\geq
0,\;(\bar{x},\,\bar{y})\geq 0$ and
$\phi(x^{k},\,y^{k})\leq\phi(x^{0},\,y^{0})$, we obtain

$\displaystyle(y^{k})^{T}(x^{k}-u^{k})+(x^{k})^{T}(y^{k}-v^{k})\leq\bar{x}^{T}(y^{k}-v^{k})+\bar{y}^{T}(x^{k}-u^{k})$
$\displaystyle\hskip
5.69054pt+(x^{k})^{T}y^{k}+(x^{k}-u^{k})^{T}(y^{k}-v^{k})+\bar{x}^{T}\bar{y}$
$\displaystyle\hskip
5.69054pt\leq\bar{x}^{T}(y^{0}-v^{0})+\bar{y}^{T}(x^{0}-u^{0})+\phi(x^{k},\,y^{k})+(x^{0}-u^{0})^{T}(y^{0}-v^{0})+\bar{x}^{T}\bar{y}$
$\displaystyle\hskip
5.69054pt\leq\bar{x}^{T}(y^{0}-v^{0})+\bar{y}^{T}(x^{0}-u^{0})+\phi(x^{0},\,y^{0})+(x^{0}-u^{0})^{T}(y^{0}-v^{0})+\bar{x}^{T}\bar{y}$
$\displaystyle\hskip
5.69054pt=\bar{x}^{T}(y^{0}-v^{0})+\bar{y}^{T}(x^{0}-u^{0})+(x^{0}-u^{0})^{T}(y^{0}-v^{0})+2n\mu_{0}+\bar{x}^{T}\bar{y}.$
(51)

We set

$\displaystyle
C_{1}=(y^{0}-v^{0})^{T}(x^{0}-u^{0})+\bar{y}^{T}(x^{0}-u^{0})+\bar{x}^{T}(y^{0}-v^{0})+\bar{x}^{T}\bar{y}+2n\mu_{0}.$

Then, from equation (51), we obtain the inequality (47).

When $(\bar{x},\,\bar{y})$ is strictly feasible, from inequalities (49)-(50),
we have

$\displaystyle\bar{y}^{T}x^{k}+\bar{x}^{T}y^{k}\leq\bar{x}^{T}(y^{k}-v^{k})+\bar{y}^{T}(x^{k}-u^{k})+\bar{x}^{T}\bar{y}+(x^{k})^{T}y^{k}+(x^{k}-u^{k})^{T}(y^{k}-v^{k})$
$\displaystyle\hskip
5.69054pt\leq\bar{x}^{T}(y^{0}-v^{0})+\bar{y}^{T}(x^{0}-u^{0})+\bar{x}^{T}\bar{y}+\phi(x^{k},\,y^{k})+(x^{0}-u^{0})^{T}(y^{0}-v^{0})$
$\displaystyle\hskip
5.69054pt\leq\bar{x}^{T}(y^{0}-v^{0})+\bar{y}^{T}(x^{0}-u^{0})+\bar{x}^{T}\bar{y}+2n\mu_{0}+(x^{0}-u^{0})^{T}(y^{0}-v^{0})\triangleq
C_{2},$ (52)

where we use the property
$\phi(x^{k},\,y^{k})\leq\phi(x^{0},\,y^{0})=2n\mu_{0}$ in the third
inequality. Consequently, from equation (52) and $(\bar{x},\,\bar{y})>0$, we
have

$\displaystyle 0\leq x^{k}_{i}\leq\frac{C_{2}}{\min_{1\leq i\leq
n}\\{\bar{y}_{i}\\}},\;0\leq y^{k}_{i}\leq\frac{C_{2}}{\min_{1\leq i\leq
n}\\{\bar{x}_{i}\\}},\;i=1,\,2,\,\ldots,\,n,$

which prove the boundedness of $(x^{k},\,y^{k})$. $\quad\square$

###### Lemma 3

Assume that $\\{(x^{k},\,y^{k})\\}$ is generated by Algorithm 1 and satisfies
the condition (38). Then, when $\mu_{k}\geq\epsilon>0$, there exists a
positive constant $\omega^{\ast}$ such that

$\displaystyle\|D^{k}\Delta x^{k}\|^{2}+\|(D^{k})^{-1}\Delta
y^{k}\|^{2}\leq(\omega^{\ast})^{2},$ (53)

holds for all $k=0,\,1,\,\ldots$, where
$D^{k}=\text{diag}\left((y^{k}./x^{k})^{1/2}\right)$,
$(D^{k})^{-1}=\text{diag}\left((x^{k}./y^{k})^{1/2}\right)$ and $(\Delta
x^{k},\,\Delta y^{k})$ is the solution of equations (16)-(17).

Proof. We denote

$\displaystyle\omega_{k}=\left(\|D^{k}\Delta x^{k}\|^{2}+\|(D^{k})^{-1}\Delta
y^{k}\|^{2}\right)^{1/2}.$ (54)

Then, from equation (54), we have

$\displaystyle\|D^{k}\Delta x^{k}\|_{\infty}\leq\|D^{k}\Delta
x^{k}\|\leq\omega_{k},\;\|(D^{k})^{-1}\Delta
y^{k}\|_{\infty}\leq\|(D^{k})^{-1}\Delta y^{k}\|\leq\omega_{k}.$ (55)

From equation (34), we know that $(u^{k},\,v^{k})$ is feasible. By combining
it with equation (16), we obtain

$\displaystyle y^{k}+\Delta y^{k}-v^{k}=M(x^{k}+\Delta x^{k}-u^{k}).$ (56)

Therefore, from equation (56) and the semi-definite positivity of $M$, we have

$\displaystyle(y^{k}+\Delta y^{k}-v^{k})^{T}(x^{k}+\Delta
x^{k}-u^{k})=(x^{k}+\Delta x^{k}-u^{k})^{T}M(x^{k}+\Delta x^{k}-u^{k})\geq 0.$
(57)

From equations (35)-(36), we know

$\displaystyle x^{0}-u^{0}\geq x^{k}-u^{k},\;y^{0}-v^{0}\geq y^{k}-v^{k}.$

By substituting them into equation (57), we obtain

$\displaystyle(x^{0}-u^{0})^{T}(y^{0}-v^{0})+(y^{k}-v^{k})^{T}\Delta
x^{k}+(x^{k}-u^{k})^{T}\Delta y^{k}+(\Delta x^{k})^{T}\Delta y^{k}$
$\displaystyle\geq(x^{k}-u^{k})^{T}(y^{k}-v^{k})+(y^{k}-v^{k})^{T}\Delta
x^{k}+(x^{k}-u^{k})^{T}\Delta y^{k}+(\Delta x^{k})^{T}\Delta y^{k}\geq 0.$
(58)

By using the inequality $|x^{T}y|\leq\|x\|\|y\|\leq\|x\|_{1}\|y\|$, from
equation (55), we have

$\displaystyle\left|(y^{k}-v^{k})^{T}\Delta
x^{k}\right|=\left|\left((D^{k})^{-1}(y^{k}-v^{k})\right)^{T}(D^{k}\Delta
x^{k})\right|$ $\displaystyle\hskip
5.69054pt\leq\left\|(D^{k})^{-1}(y^{k}-v^{k})\right\|_{1}\|D^{k}\Delta
x^{k}\|\leq\left\|(D^{k})^{-1}(y^{k}-v^{k})\right\|_{1}\omega_{k}.$ (59)

From the condition (38) and $y^{k}-v^{k}\geq 0$, we have

$\displaystyle({x^{k}_{i}}/{y^{k}_{i}})^{1/2}(y^{k}_{i}-v^{k}_{i})=x^{k}_{i}(y^{k}_{i}-v^{k}_{i})/(x^{k}_{i}y^{k}_{i})^{1/2}\leq
x^{k}_{i}(y^{k}_{i}-v^{k}_{i})/(\gamma\mu_{k})^{1/2}.$ (60)

By substituting equation (60) into equation (59), we obtain

$\displaystyle\left|(y^{k}-v^{k})^{T}\Delta
x^{k}\right|\leq\left(\sum_{i=1}^{n}({x^{k}_{i}}/{y^{k}_{i}})^{1/2}(y^{k}_{i}-v^{k}_{i})\right)\omega_{k}$
$\displaystyle\hskip
5.69054pt\leq(x^{k})^{T}(y^{k}-v^{k})\omega_{k}/(\gamma\mu_{k})^{1/2}.$ (61)

Similarly to the proof of equation (61), we obtain

$\displaystyle\left|(x^{k}-u^{k})^{T}\Delta
y^{k}\right|\leq(y^{k})^{T}(x^{k}-u^{k})\omega_{k}/(\gamma\mu_{k})^{1/2}.$
(62)

By substituting inequalities (47) and (61)-(62) into inequality (58), we
obtain

$\displaystyle(\Delta x^{k})^{T}\Delta
y^{k}\geq-(x^{0}-u^{0})^{T}(y^{0}-v^{0})-\left|(y^{k}-v^{k})^{T}\Delta
x^{k}\right|-\left|(x^{k}-u^{k})^{T}\Delta y^{k}\right|$ $\displaystyle\hskip
5.69054pt\geq-(x^{0}-u^{0})^{T}(y^{0}-v^{0})-((y^{k})^{T}(x^{k}-u^{k})+(x^{k})(y^{k}-v^{k}))\omega_{k}/(\gamma\mu_{k})^{1/2}$
$\displaystyle\hskip
5.69054pt\geq-(x^{0}-u^{0})^{T}(y^{0}-v^{0})-\left(C_{1}/(\gamma\epsilon)^{1/2}\right)\omega_{k}.$
(63)

From equation (17) and the condition (38), we have

$\displaystyle\|D^{k}\Delta x^{k}\|^{2}+\|(D^{k})^{-1}\Delta
y^{k}\|^{2}+2(\Delta x^{k})^{T}\Delta y^{k}=\|D^{k}\Delta
x^{k}+(D^{k})^{-1}\Delta y^{k}\|^{2}$ $\displaystyle\hskip
5.69054pt=(\sigma_{k}\mu_{k})^{2}\sum_{i=1}^{n}\frac{1}{x^{k}_{i}y^{k}_{i}}+(x^{k})^{T}y^{k}-2\sigma_{k}\mu_{k}n$
$\displaystyle\hskip 5.69054pt\leq
n\sigma_{k}^{2}\mu_{k}/\gamma+\phi(x^{k},\,y^{k})-2\sigma_{k}\mu_{k}n\leq
n\mu_{0}\left(\sigma_{max}^{2}/\gamma+2-2\sigma_{min}\right).$ (64)

By substituting inequality (63) into inequality (64), we obtain

$\displaystyle
q(\omega_{k})\triangleq\omega_{k}^{2}-2\left(C_{1}/(\gamma\epsilon)^{1/2}\right)\omega_{k}-\zeta\leq
0,$ (65)

where the positive constant $\zeta$ is defined by

$\displaystyle\zeta=2(x^{0}-u^{0})^{T}(y^{0}-v^{0})+n\mu_{0}(\sigma_{max}^{2}/\gamma+2-2\sigma_{min}).$
(66)

The quadratic function $q(\omega)$ is convex and has a unique positive root at

$\displaystyle\omega^{\ast}=C_{1}/(\gamma\epsilon)^{1/2}+\sqrt{C_{1}^{2}/(\gamma\epsilon)+\zeta}.$

This implies

$\displaystyle\omega_{k}\leq\omega^{\ast}.$ (67)

Consequently, we obtain inequality (53). $\quad\square$

In order to prove the global convergence of Algorithm 1, we need to estimate
the positive lower bound of $\Delta t_{k}\,(k=0,\,1,\,\ldots)$.

###### Lemma 4

Assume that $\\{(x^{k},\,y^{k})\\}$ is generated by Algorithm 1 and satisfies
the condition (38). Then, when $\mu_{k}\geq\epsilon>0$, there exists a
positive constant $\delta_{\Delta t}$ such that

$\displaystyle\Delta t_{k}\geq\frac{1}{2}\delta_{\Delta t}>0$ (68)

holds for all $k=0,\,1,\,\ldots$.

Proof. From inequality (53), we have

$\displaystyle|\Delta x^{k}_{i}\Delta
y^{k}_{i}|=\left|(y^{k}_{i}/x^{k}_{i})^{1/2}\Delta
x^{k}_{i}\right|\left|(x^{k}_{i}/y^{k}_{i})^{1/2}\Delta
y^{k}_{i}\right|\leq\|D^{k}\Delta x^{k}\|\|(D^{k})^{-1}\Delta y^{k}\|$
$\displaystyle\hskip
5.69054pt\leq\frac{1}{2}\omega_{k}^{2}\leq\frac{1}{2}(\omega^{\ast})^{2},\;i=1,\,2,\,\ldots,n,$

which gives

$\displaystyle\|\Delta X^{k}\Delta
y^{k}\|_{\infty}\leq\frac{1}{2}(\omega^{\ast})^{2},\;k=0,\,1,\,2,\,\ldots.$
(69)

Consequently, we have

$\displaystyle\frac{\sigma_{k}\mu_{k}}{\|\Delta X^{k}\Delta
y^{k}\|_{\infty}}\geq\frac{2\sigma_{min}\,\epsilon}{(\omega^{\ast})^{2}},\;k=0,\,1,\,2,\,\ldots.$
(70)

We denote

$\displaystyle\alpha^{\ast}_{\mu}=\min\left\\{1,\;\frac{1-\gamma/2}{1+\gamma/2}\,\frac{(2\sigma_{min}\,\epsilon)}{(\omega^{\ast})^{2}}\right\\}.$
(71)

Then, when $\Delta t_{k}$ satisfies

$\displaystyle\Delta
t_{k}\leq\frac{\alpha^{\ast}_{\mu}}{1-\alpha^{\ast}_{\mu}},$ (72)

from equations (70)-(71) and (42), we know that the condition (38) holds,
where $\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$.

From equation (53) and the Cauchy-Schwartz inequality
$|x^{T}y|\leq\|x\|\|y\|$, we have

$\displaystyle|(\Delta x^{k})^{T}\Delta y^{k}|=\left|(D^{k}\Delta
x^{k})^{T}((D^{k})^{-1}\Delta y^{k})\right|\leq\|D^{k}\Delta
x^{k}\|\|(D^{k})^{-1}\Delta y^{k}\|$ $\displaystyle\hskip
5.69054pt\leq\frac{1}{2}\left(\|D^{k}\Delta x^{k}\|^{2}+\|(D^{k})^{-1}\Delta
y^{k}\|^{2}\right)\leq\frac{1}{2}(\omega^{\ast})^{2}.$ (73)

From equation (17), we have

$\displaystyle(x^{k})^{T}\Delta y^{k}+(y^{k})^{T}\Delta
x^{k}=n\sigma_{k}\mu_{k}-(x^{k})^{T}y^{k}.$ (74)

By substituting equations (73)-(74) into equation (26), we obtain

$\displaystyle\rho_{k}=\frac{Ared_{k}}{Pred_{k}}=1-\frac{\alpha_{k}(\Delta
x^{k})^{T}\Delta y^{k}}{\|r_{q}^{k}\|+(x^{k})^{T}y^{k}-n\sigma_{k}\mu_{k}}$
$\displaystyle\hskip 5.69054pt=1-\frac{\alpha_{k}(\Delta x^{k})^{T}\Delta
y^{k}}{(2-\sigma_{k})n\mu_{k}}\geq
1-\frac{(\omega^{\ast})^{2}}{2(2-\sigma_{max})n\epsilon}\alpha_{k}.$ (75)

We denote

$\displaystyle\alpha_{\rho}^{\ast}=\min\left\\{1,\;\frac{2(2-\sigma_{max})(1-\eta_{2})n\epsilon}{(\omega^{\ast})^{2}}\right\\}.$
(76)

Then, when $\Delta t_{k}$ satisfies

$\displaystyle\Delta
t_{k}\leq\frac{\alpha_{\rho}^{\ast}}{1-\alpha_{\rho}^{\ast}},$ (77)

from equations (75)-(76), we know $\rho_{k}\geq\eta_{2}$, where
$\alpha_{k}=\Delta t_{k}/(1+\Delta t_{k})$.

We denote

$\displaystyle\delta_{\Delta
t}\triangleq\min\left\\{\frac{\sigma_{\mu}^{\ast}}{1-\sigma_{\mu}^{\ast}},\;\frac{\alpha_{\rho}^{\ast}}{1-\alpha_{\rho}^{\ast}},\;\Delta
t_{0}\right\\}.$ (78)

Assume that $K$ is the first index such that $\Delta t_{K}\leq\delta_{\Delta
t}$ holds. Then, from equations (72)-(78), we know that the condition (38)
holds and $\rho_{K}\geq\eta_{2}$. Consequently, according to the time-stepping
adjustment scheme (27), $\Delta t_{K+1}$ will be enlarged. Therefore, $\Delta
t_{k}\geq(1/2)\delta_{\Delta t}$ holds for all $k=0,\,1,\,2,\,\dots$. ∎

By using the estimation results of Lemma 4, we prove that $\\{\mu_{k}\\}$
converges to zero when $k$ tends to infinity as follows.

###### Theorem 3.1

Assume that $\\{(x^{k},\,y^{k})\\}$ is generated by Algorithm 1 and satisfies
the condition (38). Then, we have

$\displaystyle\lim_{k\to\infty}\mu_{k}=0.$ (79)

Proof. We prove the result (79) by contradiction. Assume that there exists a
positive constant $\epsilon$ such that

$\displaystyle\mu_{k}\geq\epsilon>0$ (80)

holds for all $k=0,\,1,\,\ldots$. Then, according to Algorithm 1 and Lemma 4,
we know that there exists an infinite subsequence
$\\{(x^{k_{l}},\,y^{k_{l}})\\}$ such that

$\displaystyle\frac{\phi(x^{k_{l}},\,y^{k_{l}})-\phi(x^{k_{l}}+\alpha_{k_{l}}\Delta
x^{k_{l}},\,y^{k_{l}}+\alpha_{k_{l}}\Delta
y^{k_{l}})}{\phi(x^{k_{l}},\,y^{k_{l}})-m_{k_{l}}(\alpha_{k_{l}})}\geq\eta_{1}$
(81)

holds for all $l=0,\,1,\,2,\,\ldots$, where $\alpha_{k_{l}}=\Delta
t_{k_{l}}/(1+\Delta t_{k_{l}})$. Otherwise, all steps are rejected after a
given iteration index, then the time step will keep decreasing, which
contradicts equation (68).

From equations (68), (75) and (81), we have

$\displaystyle\phi(x^{k_{l}},\,y^{k_{l}})-\phi(x^{k_{l}}+\alpha_{k_{l}}\Delta
x^{k_{l}},\,y^{k_{l}}+\alpha_{k_{l}}\Delta y^{k_{l}})$ $\displaystyle\hskip
5.69054pt\geq\frac{\eta_{1}\Delta t_{k_{l}}}{1+\Delta
t_{k_{l}}}(2-\sigma_{k_{l}})n\mu_{k_{l}}\geq\frac{\eta_{1}\delta_{\Delta
t}/2}{1+\delta_{\Delta t}/2}(2-\sigma_{max})n\mu_{k_{l}}.$ (82)

Therefore, from equation (82) and
$\phi(x^{k+1},\,y^{k+1})\leq\phi(x^{k},\,y^{k})$, we have

$\displaystyle\phi(x^{0},\,y^{0})\geq\phi(x^{0},\,y^{0})-\lim_{k\to\infty}\phi(x^{k},\,y^{k})=\sum_{k=0}^{\infty}\left(\phi(x^{k},\,y^{k})-\phi(x^{k+1},\,y^{k+1})\right)$
$\displaystyle\hskip
5.69054pt\geq\sum_{l=0}^{\infty}\left(\phi(x^{k_{l}},\,y^{k_{l}})-\phi(x^{k_{l}+1},\,y^{k_{l}+1})\right)\geq\frac{\eta_{1}\delta_{\Delta
t}/2}{1+\delta_{\Delta t}/2}\sum_{l=0}^{\infty}(2-\sigma_{max})n\mu_{k_{l}}.$
(83)

Consequently, from equation (83), we obtain

$\displaystyle\lim_{l\to\infty}\mu_{k_{l}}=0,$

which contradicts the assumption $\mu_{k}\geq\epsilon>0$ for all
$k=0,\,1,\ldots$. Therefore, we have

$\displaystyle\lim_{k\to\infty}\inf\mu_{k}=0.$ (84)

Since $\mu_{k}=\phi(x^{k},\,y^{k})/(2n)$ and $\\{\phi(x^{k},\,y^{k})\\}$ is
monotonically decreasing, we know that $\\{\mu_{k}\\}$ is monotonically
decreasing. By combining it with equation (84), we know that the result (79)
is true. ∎

###### Remark 4

In the analysis framework of the global convergence, in comparison to that of
the known path-following algorithm, we do not need to select $\alpha_{k}$ such
that the condition of the priority to feasibility over complementarity

$\displaystyle\|r_{q}^{k}(\alpha)\|\leq(x^{k}(\alpha))^{T}y^{k}(\alpha)\frac{\|r_{q}^{0}\|}{(x^{0})^{T}y^{0}}$
(85)

holds for all $\alpha\in(0,\;\alpha_{k}]\subset(0,\;1]$.

###### Theorem 3.2

Assume that $\\{(x^{k},\,y^{k})\\}$ is the sequence of iterations generated by
Algorithm 1 and satisfies the condition (38). If the linear complementarity
problem has a strictly feasible solution $(\bar{x},\,\bar{y})$,
$\\{(x^{k},\,y^{k})\\}$ exists a limit point $(x^{\ast},\,y^{\ast})$ and this
limit point satisfies the linear complementarity equation (1).

Proof. Since the linear complementarity problem has a strictly feasible
solution, from Lemma 2, we know that $\\{(x^{k},\,y^{k})\\}$ is bounded.
Consequently, the sequence $\\{(x^{k},\,y^{k})\\}$ exists a convergence
subsequence $\\{(x^{k_{l}},\,y^{k_{l}})\\}$ and we denote its limit point as
$(x^{\ast},\,y^{\ast})$. By combining it with Theorem 3.1, we obtain

$\displaystyle\lim_{l\to\infty}\mu_{k_{l}}=\frac{1}{2n}\left(\lim_{l\to\infty}(x^{k_{l}})^{T}y^{k_{l}}+\lim_{l\to\infty}\|y^{k_{l}}-(Mx^{k_{l}}+q)\|\right)$
$\displaystyle\hskip
5.69054pt=\frac{1}{2n}\left((x^{\ast})^{T}y^{\ast}+\|y^{\ast}-(Mx^{\ast}+q)\|\right)=0.$
(86)

Since $(x^{k_{l}},\,y^{k_{l}})\geq 0$, we have $(x^{\ast},\,y^{\ast})\geq 0$.
By substituting it into equation (86), we conclude

$\displaystyle
y^{\ast}-(Mx^{\ast}+q)=0,\;(x^{\ast})^{T}y^{\ast}=0,\;(x^{\ast},\,y^{\ast})\geq
0.$

It implies that $(x^{\ast},\,y^{\ast})$ is the solution of the linear
complementarity problem (1). $\quad\square$

##  4 Numerical experiments

In this section, some numerical experiments are conducted to test the
performance of RPFMTr (Algorithm 1) for the linear complementarity problems
(LCPs), in comparison to two state-of-the-art commercial solvers, i.e. PATH
DF1995 ; PATH ; FM2000 ; FM2022  and MILES (a Mixed Inequality and nonLinear
Equation Solver) Mathiesen1985 ; Rutherford1995 ; Rutherford2022 . RPFMTr
(Algorithm 1) is coded with the MATLAB language and executed in MATLAB
(R2020a) environment MATLAB . PATH and MILES are executed in the GAMS v28.2
(2019) environment GAMS . All numerical experiments are performed by a HP
notebook with the Intel quadcore CPU and 8Gb memory. MILES solves the LCP
based on the Lemke’s almost complementary pivoting algorithm CPS2009 ; LH1964
; Rutherford1995 . PATH is a solver based on the path-following procedure
Dirkse1994 ; DF1995 , the Fisher’s non-smooth regularization technique
Fisher1992  and the Lemke’s pivoting algorithm CPS2009 ; LH1964 . PATH and
MILES are two robust and efficient solvers for the complementarity problems
and used in many general modelling systems such as AMPL (A Mathematical
Programming Language) AMPL  and GAMS (General Algebraic Modelling System) GAMS
. Therefore, we select these two solvers as the basis for comparison.

The construction approach of test problems is described as follows. Firstly,
we generate a test matrix $M$ with the following structure:

$\displaystyle M=\begin{bmatrix}0&-A^{T}\\\ A&0\end{bmatrix},$ (87)

where the test matrix $A$ comes from the linear programming subset of NETLIB
NETLIB . It is easy to verify that the test matrix $M$ defined by equation
(87) is semi-definite positive. Then, we generate two complementarity vectors
$x$ and $y$ as follows:

$\displaystyle
x=[1\;0\;1\;0\;\cdots\;1\;0]^{T},\;y=[0\;1\;0\;1\;\cdots\;0\;1]^{T}.$ (88)

Finally, by using these two vectors $x,\,y$ defined by equation (88) and the
test matrix $M$ defined by equation (87), we generate the following test
vector $q$:

$\displaystyle q=y-Mx.$ (89)

The scales of test problems vary from dimension 78 to 40216. And the
termination conditions of Algorithm 1 (RPFMTr), PATH and MILES for finding a
solution of the nonlinear system (2) are

$\displaystyle\|y-(Mx+q)\|_{\infty}\leq 10^{-6},\;\|Xy\|_{\infty}\leq
10^{-6},\;(x,\,y)\geq 0,$

where $X=\text{diag}(x)$.

According to the manual of GAMS FM2022 , an LCP will be generated a list of
all variables appearing in the equations found in the model statement, and the
number of equations equals the number of variables. In other words, for an LCP
solved in the GAMS environment, its matrix $M$ should not include the row with
all zeros. In order to avoid this error, we add $\epsilon=1.0\times 10^{-6}$
to the first element of the row or column with all zeros in the test sub-
matrix $A$ of matrix $M$ defined by equation (87).

Most of matrices $A$ coming from the linear programming subset of NETLIB
NETLIB  are sparse. In order to test the performance of RPFMTr, PATH and MILES
for the dense LCP further, we add a small random perturbation to the elements
of the sub-matrix $A$ in equation (87) and generate the dense test matrix
$\overline{M}$ to replace the test matrix $M$ in equation (87) as follows:

$\displaystyle\overline{M}=\begin{bmatrix}0&-A_{\epsilon}^{T}\\\
A_{\epsilon}&0\end{bmatrix},\;A_{\epsilon}=A+\text{rand}(m_{A},\,n_{A})*\epsilon,\;[m_{A},\,n_{A}]=\text{sizes}(A),$
(90)

where $\epsilon=10^{-3}$ and matrix $A$ comes from the linear programming
subset of NETLIB NETLIB .

Numerical results of the sparse test problems are arranged in Tables 1-3 and
Figures 1-2. And numerical results of the dense test problems are arranged in
Tables 4-6 and Figures 3-4. “major” in the fourth column of Tables 1-6
represents the number of the linear mixed complementarity problems solved by a
pivotal Lemke’s method of PATH FM2022 . “minor” in the fourth column of Tables
1-6 represents the the number of pivots performed per major iteration of PATH
FM2022 . “grad” in the fourth column of Tables 1-6 represents the cumulative
number of Jacobian evaluations used in PATH FM2022 . “major” in the sixth
column of Tables 1-6 represents the number of Newton iterations of MILES
Rutherford2022 . “pivots” in the sixth column of Tables 1-6 represents the
number of Lemke pivots of MILES Rutherford2022 . “refactor” in the sixth
column of Tables 1-6 represents the number of re-factorizations in the LUSOL
solver GMSW1991 , which is called by MILES for solving the linear systems of
equations.

From Tables 1-6 and Figures 1-4, we find that RPFMTr and PATH work well for
the sparse LCPs and the dense LCPs. However, MILES is not robust for solving
the sparse LCPs or the dense LCPs. For the sparse LCPs, PATH performs better
than RPFMTr and MILES from Tables 1-3 and Figures 1-2. For the dense LCPs,
from Tables 4-6, we find that PATH and MILES fail to solve $5$ problems and
$57$ problems of 73 test problems, respectively. RPFMTr solves all the sparse
test LCPs and the dense test LCPs. Furthermore, from Tables 4-6 and Figures
3-4, we find that the computational time of RPFMTr is about $1/3$ to $1/10$ of
that of PATH for the dense LCPs. Therefore, RPFMTr is a robust and efficient
solver for the LCPs, especially for the dense LCPs.

Table 1: Numerical results of small-scale sparse LCPs (no. 1-27).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 1 lp_25fv47
---
(n = 2697)
| 49
---
(4.13)
| 8.84e-07
---
| 8+949+23
---
(0.30)
| 5.82e-08
---
| 1+1684+24
---
(1.34)
| 7.48e-12
---
| Exam. 2 lp_adlittle
---
(n = 194)
| 45
---
(0.03)
| 9.61e-07
---
| 4+80+9
---
(0.02)
| 4.44e-12
---
| 1+1+3
---
(0.02)
| 4.65e-14
---
| Exam. 3 lp_afiro
---
(n = 78)
| 41
---
(0.01)
| 6.20e-07
---
| 3+5+7
---
(0.02)
| 1.01e-10
---
| 1+1+3
---
(0.00)
| 3.91e-14
---
| Exam. 4 lp_agg
---
(n = 1103)
| 44
---
(0.19)
| 9.85e-07
---
| 11+413+16
---
(0.08)
| 2.87e-09
---
| 1+1+6
---
(0.05)
| 3.74e-12
---
| Exam. 5 lp_agg2
---
(n = 1274)
| 46
---
(0.84)
| 8.24e-07
---
| 10+303+13
---
(0.03)
| 6.96e-10
---
| 1+959+6
---
(0.06)
| 5.08e-12
---
| Exam. 6 lp_agg3
---
(n = 1274)
| 47
---
(0.50)
| 6.64e-07
---
| 9+259+12
---
(0.05)
| 2.41e-07
---
| 1+955+6
---
(0.06)
| 9.37e-12
---
| Exam. 7 lp_bandm
---
(n = 777)
| 51
---
(0.27)
| 8.58e-07
---
| 11+613+14
---
(0.05)
| 3.31e-07
---
| 1+1+8
---
(0.06)
| 1.19e-12
---
| Exam. 8 lp_beaconfd
---
(n = 468)
| 50
---
(0.17)
| 6.06e-07
---
| 8+124+17
---
(0.03)
| 3.76e-09
---
| 1+1+6
---
(0.03)
| 5.46e-12
---
| Exam. 9 lp_blend
---
(n = 188)
| 46
---
(0.03)
| 6.12e-07
---
| 9+143+12
---
(0.02)
| 6.50e-08
---
| 1+1+4
---
(0.02)
| 1.96e-12
---
| Exam. 10 lp_bnl1
---
(n = 2229)
| 50
---
(1.24)
| 8.60e-07
---
| 12+1727+17
---
(0.50)
| 5.22e-11
---
| 1+1+10
---
(0.27)
| 1.62e-12
---
| Exam. 11 lp_bore3d
---
(n = 567)
| 48
---
(0.13)
| 8.63e-07
---
| 11+423+16
---
(0.03)
| 6.99e-07
---
| 1+1+6
---
(0.03)
| 1.31e-12
---
| Exam. 12 lp_brandy
---
(n = 523)
| 49
---
(0.31)
| 7.42e-07
---
| 10+384+13
---
(0.06)
| 2.84e-10
---
| 1+1+6
---
(0.03)
| 1.31e-12
---
| Exam. 13 lp_capri
---
(n = 753)
| 46
---
(0.16)
| 9.07e-07
---
| 11+573+15
---
(0.03)
| 7.18e-07
---
| 1+1+6
---
(0.05)
| 2.60e-12
---
| Exam. 14 lp_czprob
---
(n = 4491)
| 51
---
(3.37)
| 8.74e-07
---
| 4+25+14
---
(0.05)
| 1.90e-10
---
| 1+1+13
---
(1.19)
| 1.82e-12
---
| Exam. 15 lp_degen2
---
(n = 1201)
| 42
---
(0.60)
| 8.53e-07
---
| 4+447+13
---
(0.11)
| 1.71e-07
---
| 1+1+10
---
(0.14)
| 2.22e-16
---
| Exam. 16 lp_degen3
---
(n = 4107)
| 47
---
(5.85)
| 6.85e-07
---
| 4+185+26
---
(0.42)
| 3.00e-07
---
| 1+1+51
---
(6.13)
| 0.00e-00
---
| Exam. 17 lp_e226
---
(n = 695)
| 54
---
(0.28)
| 9.05e-07
---
| 9+412+17
---
(0.03)
| 2.94e-07
---
| 1+1+7
---
(0.05)
| 5.86e-12
---
| Exam. 18 lp_etamacro
---
(n = 1216)
| 49
---
(0.73)
| 6.31e-07
---
| 6+440+15
---
(0.05)
| 1.87e-12
---
| 1+58+6
---
(0.05)
| 6.82e-12
---
| Exam. 19 lp_fffff800
---
(n = 1552)
| 64
---
(1.27)
| 6.05e-07
---
| 15+1111+20
---
(0.17)
| 1.39e-07
---
| 100+72624+770
---
(13.80)
| 1.78e+01
---
(failed)
| Exam. 20 lp_finnis
---
(n = 1561)
| 47
---
(0.26)
| 6.06e-07
---
| 8+313+16
---
(0.03)
| 5.24e-09
---
| 1+1+9
---
(0.14)
| 4.55e-13
---
| Exam. 21 lp_fit1d
---
(n = 1073)
| 65
---
(1.16)
| 6.45e-07
---
| 8+671+12
---
(0.24)
| 2.34e-11
---
| 1+1+3
---
(0.03)
| 1.75e-10
---
| Exam. 22 lp_ganges
---
(n = 3015)
| 43
---
(0.49)
| 8.88e-07
---
| 4+98+21
---
(0.05)
| 4.95e-07
---
| 1+1+11
---
(0.47)
| 2.40e-14
---
| Exam. 23 lp_gfrd_pnc
---
(n = 1776)
| 61
---
(0.28)
| 7.85e-07
---
| 18+1250+21
---
(0.17)
| 3.80e-11
---
| 1+716+13
---
(0.25)
| 3.64e-12
---
| Exam. 24 lp_grow15
---
(n = 945)
| 33
---
(0.24)
| 7.23e-07
---
| 7+703+9
---
(0.13)
| 8.66e-07
---
| 1+927+6
---
(0.08)
| 3.33e-14
---
| Exam. 25 lp_grow22
---
(n = 1386)
| 33
---
(0.42)
| 7.27e-07
---
| 8+1054+10
---
(0.23)
| 3.86e-11
---
| 1+1403+8
---
(0.13)
| 1.20e-14
---
| Exam. 26 lp_lotfi
---
(n = 519)
| 52
---
(0.20)
| 6.05e-07
---
| 13+263+17
---
(0.03)
| 8.42e-07
---
| 1+1+4
---
(0.05)
| 9.10e-13
---
| Exam. 27 lp_maros
---
(n = 2812)
| 68
---
(2.80)
| 6.18e-07
---
| 13+3222+14
---
(0.50)
| 2.16e-09
---
| 1+776+20
---
(1.13)
| 2.95e-09
---

Table 2: Numerical results of small-scale sparse LCPs (no. 28-53).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 28 lp_modszk1
---
(n = 2307)
| 42
---
(0.46)
| 8.61e-07
---
| 5+370+14
---
(0.05)
| 7.91e-11
---
| 1+1+9
---
(0.22)
| 8.70e-14
---
| Exam. 29 lp_perold
---
(n = 2131)
| 67
---
(2.62)
| 9.15e-07
---
| 12+1498+23
---
(0.42)
| 1.72e-08
---
| 100+18168+1156
---
(35.88)
| 6.64e+01
---
(failed)
| Exam. 30 lp_pilot_ja
---
(n = 3207)
| 78
---
(7.43)
| 7.94e-07
---
| 14+2468+23
---
(1.81)
| 3.29e-07
---
| 100+171048+1900
---
(140.47)
| 3.92e+04
---
(failed)
| Exam. 31 lp_qap8
---
(n = 2544)
| 43
---
(4.78)
| 8.80e-07
---
| 9+2413+11
---
(1.31)
| 1.83e-11
---
| 6+50530+477
---
(27.30)
| 2.13e-14
---
| Exam. 32 lp_lp_recipe
---
(n = 295)
| 53
---
(0.04)
| 7.81e-07
---
| 8+103+11
---
(0.00)
| 7.33e-12
---
| 1+1+4
---
(0.02)
| 5.12e-13
---
| Exam. 33 lp_sc50a
---
(n = 128)
| 40
---
(0.01)
| 8.30e-07
---
| 3+59+6
---
(0.05)
| 9.21e-07
---
| 1+1+3
---
(0.03)
| 1.33e-15
---
| Exam. 34 lp_scagr7
---
(n = 314)
| 42
---
(0.03)
| 6.28e-07
---
| 5+110+8
---
(0.02)
| 1.41e-10
---
| 1+1+3
---
(0.00)
| 1.07e-14
---
| Exam. 35 lp_scagr25
---
(n = 1142)
| 42
---
(0.08)
| 6.30e-07
---
| 4+205+13
---
(0.02)
| 4.24e-07
---
| 1+1+6
---
(0.05)
| 1.07e-14
---
| Exam. 36 lp_scfxm1
---
(n = 930)
| 52
---
(0.51)
| 9.05e-07
---
| 12+721+15
---
(0.11)
| 7.23e-10
---
| 1+1+7
---
(0.06)
| 2.62e-12
---
| Exam. 37 lp_scfxm2
---
(n = 1860)
| 52
---
(0.61)
| 8.84e-07
---
| 14+1357+18
---
(0.14)
| 3.36e-09
---
| 1+1+11
---
(0.24)
| 1.00e-11
---
| Exam. 38 lp_scfxm3
---
(n = 2790)
| 53
---
(1.14)
| 6.17e-07
---
| 14+2110+17
---
(0.27)
| 4.96e-07
---
| 1+200+14
---
(0.70)
| 7.38e-12
---
| Exam. 39 lp_scorpion
---
(n = 854)
| 40
---
(0.08)
| 8.40e-07
---
| 5+111+15
---
(0.03)
| 1.35e-07
---
| 1+1+5
---
(0.03)
| 6.22e-15
---
| Exam. 40 lp_shell
---
(n = 2313)
| 45
---
(0.40)
| 6.34e-07
---
| 9+2448+11
---
(0.38)
| 5.86e-10
---
| 1+1+9
---
(0.24)
| 0.00e-00
---
| Exam. 41 lp_ship04l
---
(n = 2568)
| 50
---
(1.06)
| 7.03e-07
---
| 11+3359+14
---
(0.77)
| 4.05e-07
---
| 1+1+8
---
(0.23)
| 1.81e-14
---
| Exam. 42 lp_ship04s
---
(n = 1908)
| 49
---
(0.57)
| 9.60e-07
---
| 12+2376+15
---
(0.24)
| 4.04e-07
---
| 1+1+6
---
(0.09)
| 1.81e-14
---
| Exam. 43 lp_ship08s
---
(n = 3245)
| 47
---
(0.91)
| 5.31e-07
---
| 12+3328+15
---
(0.44)
| 9.32e-07
---
| 1+1+9
---
(0.41)
| 5.68e-14
---
| Exam. 44 lp_ship12s
---
(n = 4020)
| 48
---
(0.91)
| 2.68e-07
---
| 13+3446+16
---
(0.52)
| 3.14e-07
---
| 1+1+9
---
(0.53)
| 4.17e-14
---
| Exam. 45 lp_sierra
---
(n = 3962)
| 63
---
(1.35)
| 6.79e-07
---
| 12+2135+23
---
(0.53)
| 7.90e-07
---
| 1+1+11
---
(0.72)
| 4.37e-11
---
| Exam. 46 lp_standgub
---
(n = 1744)
| 60
---
(0.70)
| 9.29e-07
---
| 4+99+12
---
(0.02)
| 2.88e-12
---
| 1+1+7
---
(0.09)
| 1.16e-13
---
| Exam. 47 lp_tuff
---
(n = 961)
| 61
---
(0.69)
| 3.85e-07
---
| 9+828+13
---
(0.13)
| 4.81e-07
---
| 1+1+10
---
(0.11)
| 3.82e-11
---
| Exam. 48 lp_wood1p
---
(n = 2839)
| 53
---
(6.00)
| 7.48e-07
---
| 7+1285+12
---
(19.56)
| 9.65e-10
---
| 1+3237+19
---
(1.36)
| 5.09e-09
---
| Exam. 49 lpi_box1
---
(n = 492)
| 39
---
(0.03)
| 6.17e-07
---
| 7+290+9
---
(0.02)
| 7.28e-07
---
| 1+279+3
---
(0.02)
| 1.11e-16
---
| Exam. 50 lpi_cplex2
---
(n = 602)
| 41
---
(0.07)
| 7.33e-07
---
| 9+422+11
---
(0.03)
| 1.45e-09
---
| 1+1+4
---
(0.03)
| 4.44e-15
---
| Exam. 51 lpi_ex72a
---
(n = 412)
| 40
---
(0.03)
| 6.16e-07
---
| 7+164+9
---
(0.00)
| 1.24e-12
---
| 1+1+3
---
(0.02)
| 0.00e-00
---
| Exam. 52 lpi_ex73a
---
(n = 404)
| 40
---
(0.03)
| 7.56e-07
---
| 6+125+8
---
(0.00)
| 6.75e-10
---
| 1+1+3
---
(0.02)
| 0.00e-00
---
| Exam. 53 lpi_mondou2
---
(n = 916)
| 38
---
(0.09)
| 6.85e-07
---
| 6+555+8
---
(0.08)
| 5.02e-07
---
| 1+1+5
---
(0.03)
| 0.00e-00
---

Table 3: Numerical results of large-scale sparse LCPs (no. 54-78).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 54 lp_80bau3b
---
(n = 14323)
| 49
---
(7.56)
| 6.82e-07
---
| 9+858+23
---
(0.30)
| 1.52e-10
---
| 41+101114+1107
---
(1015.16)
| 6.19e+01
---
(failed)
| Exam. 55 lp_bnl2
---
(n = 6810)
| 47
---
(4.83)
| 7.01e-07
---
| 12+4597+17
---
(1.47)
| 2.72e-11
---
| 1+1+27
---
(6.84)
| 1.24e-12
---
| Exam. 56 lp_cre_a
---
(n = 10764)
| 53
---
(5.14)
| 8.39e-07
---
| 12+5755+17
---
(2.22)
| 4.09e-07
---
| 1+1+42
---
(22.98)
| 1.48e-12
---
| Exam. 57 lp_cre_c
---
(n = 9479)
| 50
---
(4.50)
| 8.58e-07
---
| 11+3604+14
---
(1.34)
| 1.77e-09
---
| 1+490+40
---
(18.56)
| 1.59e-12
---
| Exam. 58 lp_cycle
---
(n = 5274)
| 59
---
(5.21)
| 4.59e-07
---
| 17+4532+24
---
(2.58)
| 8.67e-09
---
| 100+197608+2800
---
(461.45)
| 7.15e+03
---
(failed)
| Exam. 59 lp_d6cube
---
(n = 6599)
| 56
---
(48.11)
| 2.56e-07
---
| 5+941+19
---
(0.53)
| 7.37e-07
---
| 1+6085+39
---
(10.75)
| 2.73e-12
---
| Exam. 60 lp_fit2d
---
(n = 10549)
| 69
---
(68.27)
| 6.72e-07
---
| 7+2114+10
---
(21.23)
| 1.53e-09
---
| 1+1+3
---
(0.72)
| 3.49e-10
---
| Exam. 61 lp_greenbea
---
(n = 7990)
| 50
---
(13.68)
| 7.29e-07
---
| 12+3899+28
---
(2.91)
| 1.01e-09
---
| 100+74070+1754
---
(513.08)
| 9.90e+01
---
(failed)
| Exam. 62 lp_greenbeb
---
(n = 7990)
| 50
---
(14.48)
| 7.29e-07
---
| 12+3899+28
---
(3.19)
| 1.01e-09
---
| 100+74070+1754
---
(514.13)
| 9.90e+01
---
(failed)
| Exam. 63 lp_ken_07
---
(n = 6028)
| 42
---
(1.09)
| 6.24e-07
---
| 8+4004+11
---
(0.98)
| 9.71e-09
---
| 1+1+21
---
(3.31)
| 0.00e-00
---
| Exam. 64 lp_pds_02
---
(n = 10669)
| 38
---
(2.66)
| 9.65e-07
---
| 5+3264+26
---
(3.27)
| 7.09e-09
---
| 1+1+47
---
(24.70)
| 0.00e-00
---
| Exam. 65 lp_pilot87
---
(n = 8710)
| 59
---
(59.97)
| 7.96e-07
---
| 10+3002+21
---
(2.36)
| 6.47e-07
---
| 100+121000+1500
---
(498.08)
| 1.00e+03
---
(failed)
| Exam. 66 lp_qap12
---
(n = 12048)
| 45
---
(223.00)
| 6.25e-07
---
| 7+9138+9
---
(205.47)
| 8.54e-07
---
| 12+145940+1427
---
(1064.42)
| 1.07e+01
---
(failed)
| Exam. 67 lp_ship08l
---
(n = 5141)
| 48
---
(1.87)
| 5.86e-07
---
| 10+1676+21
---
(0.80)
| 5.38e-07
---
| 1+1+12
---
(1.31)
| 5.68e-14
---
| Exam. 68 lp_ship12l
---
(n = 6684)
| 47
---
(1.88)
| 8.51e-07
---
| 14+8950+17
---
(4.24)
| 1.22e-10
---
| 1+1+16
---
(3.03)
| 4.17e-14
---
| Exam. 69 lp_truss
---
(n = 9806)
| 44
---
(4.72)
| 6.22e-07
---
| 6+5607+9
---
(6.92)
| 4.99e-08
---
| 98+299152+2379
---
(1011.18)
| 1.67e+00
---
(failed)
| Exam. 70 lp_woodw
---
(n = 9516)
| 86
---
(201.33)
| 6.31e-07
---
| 10+4715+19
---
(9.89)
| 1.78e-09
---
| 100+232954+1904
---
(735.08)
| 2.00e+05
---
(failed)
| Exam. 71 lpi_bgindy
---
(n = 13551)
| 44
---
(50.77)
| 8.49e-07
---
| 4+614+25
---
(1.28)
| 1.20e-07
---
| 1+1+39
---
(34.14)
| 9.06e-08
---
| Exam. 72 lpi_gran
---
(n = 5183)
| 81
---
(7.50)
| 9.94e-07
---
| 15+1988+28
---
(1.27)
| 8.14e-07
---
| 100+70036+700
---
(74.33)
| 5.68e+03
---
(failed)
| Exam. 73 lpi_greenbea
---
(n = 7989)
| 50
---
(11.91)
| 9.94e-07
---
| 11+3173+29
---
(2.94)
| 4.68e-07
---
| 100+303396+2731
---
(921.31)
| 6.72e+01
---
(failed)
| Exam. 74 lp_ken_11
---
(n = 36043)
| 43
---
(16.37)
| 6.95e-07
---
| 9+23833+12
---
(172.78)
| 2.65e-08
---
| 1+1+142
---
(1081.75)
| 0.00e-00
---
| Exam. 75 lp_osa_07
---
(n = 26185)
| 45
---
(451.93)
| 6.55e-07
---
| 4+954+17
---
(2.64)
| 1.77e-11
---
| 1+1+23
---
(78.39)
| 3.87e-09
---
| Exam. 76 lp_pds_06
---
(n = 39232)
| 38
---
(115.16)
| 9.02e-07
---
| 5+6809+37
---
(43.47)
| 9.00e-11
---
| 1+1+240
---
(1885.67)
| 0.00e-00
---
| Exam. 77 lp_qap15
---
(n = 28605)
| 46
---
(824.18)
| 6.34e-07
---
| 5+20058+8
---
(18004.75)
| 2.11e+01
---
(failed)
| 3+39920+342
---
(1410.55)
| 1.20e+01
---
(failed)
| Exam. 78 lp_stocfor3
---
(n = 40216)
| 54
---
(11.08)
| 8.52e-07
---
| 21+1623+50
---
(5.69)
| 2.78e-12
---
| 2+4776+135
---
(1169.09)
| 5.14e+01
---
(failed)

Figure 1: Iterations of RPFMTr, PATH and MILES for sparse LCPs. Figure 2: CPU
time of RPFMTr, PATH and MILES for sparse LCPs.

Table 4: Numerical results of small-scale dense LCPs (no. 1-27).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 1 lp_25fv47
---
(n = 2697)
| 56
---
(16.83)
| 5.85e-07
---
| 12+3457+22
---
(200.42)
| 2.78e-09
---
| 100+61000+600
---
(132.30)
| 3.56e+02
---
(failed)
| Exam. 2 lp_adlittle
---
(n = 194)
| 43
---
(0.26)
| 7.90e-07
---
| 8+203+10
---
(0.03)
| 1.33e-11
---
| 1+327+3
---
(0.02)
| 1.38e-12
---
| Exam. 3 lp_afiro
---
(n = 78)
| 40
---
(0.03)
| 7.49e-07
---
| 9+73+11
---
(0.03)
| 2.58e-11
---
| 1+77+2
---
(0.08)
| 2.22e-14
---
| Exam. 4 lp_agg
---
(n = 1103)
| 53
---
(2.78)
| 7.48e-07
---
| 12+1284+17
---
(5.20)
| 5.24e-08
---
| 100+26000+600
---
(47.64)
| 1.43e+00
---
(failed)
| Exam. 5 lp_agg2
---
(n = 1274)
| 54
---
(2.80)
| 7.30e-07
---
| 10+1071+13
---
(4.61)
| 5.44e-07
---
| 1+132+8
---
(2.30)
| 3.58e-09
---
| Exam. 6 lp_agg3
---
(n = 1274)
| 44
---
(1.64)
| 7.03e-07
---
| 10+1053+13
---
(4.58)
| 1.51e-07
---
| 1+1335+8
---
(2.25)
| 3.42e-11
---
| Exam. 7 lp_bandm
---
(n = 777)
| 51
---
(0.78)
| 9.17e-07
---
| 14+1075+18
---
(2.80)
| 8.96e-11
---
| 100+23867+901
---
(37.73)
| 2.74e+02
---
(failed)
| Exam. 8 lp_beaconfd
---
(n = 468)
| 51
---
(0.25)
| 8.57e-07
---
| 14+646+17
---
(0.45)
| 9.99e-07
---
| 100+25300+500
---
(4.72)
| 1.76e+02
---
(failed)
| Exam. 9 lp_blend
---
(n = 188)
| 44
---
(0.03)
| 5.39e-07
---
| 13+258+15
---
(0.03)
| 4.73e-08
---
| 1+547+4
---
(0.03)
| 1.83e-09
---
| Exam. 10 lp_bnl1
---
(n = 2229)
| 52
---
(9.88)
| 6.71e-07
---
| 13+3083+17
---
(56.92)
| 1.65e-09
---
| 100+26600+800
---
(121.78)
| 7.79e+01
---
(failed)
| Exam. 11 lp_bore3d
---
(n = 567)
| 50
---
(0.40)
| 9.26e-07
---
| 13+745+18
---
(0.81)
| 8.51e-08
---
| 2+1525+12
---
(0.38)
| 1.05e-11
---
| Exam. 12 lp_brandy
---
(n = 523)
| 53
---
(0.34)
| 5.13e-07
---
| 13+697+16
---
(0.53)
| 2.73e-07
---
| 100+71600+600
---
(11.34)
| 2.72e+01
---
(failed)
| Exam. 13 lp_capri
---
(n = 753)
| 52
---
(0.75)
| 5.41e-07
---
| 10+761+19
---
(1.45)
| 2.37e-09
---
| 100+79200+800
---
(37.36)
| 2.10e+01
---
(failed)
| Exam. 14 lp_czprob
---
(n = 4491)
| 44
---
(42.87)
| 8.31e-07
---
| 9+3102+16
---
(650.25)
| 3.36e-07
---
| 100+6200+500
---
(120.50)
| 5.33e+00
---
(failed)
| Exam. 15 lp_degen2
---
(n = 1201)
| 49
---
(2.32)
| 7.63e-07
---
| 12+1582+14
---
(7.38)
| 2.83e-07
---
| 3+566+16
---
(2.25)
| 2.47e-09
---
| Exam. 16 lp_degen3
---
(n = 4107)
| 45
---
(35.17)
| 6.85e-07
---
| 9+4003+20
---
(388.39)
| 9.91e-07
---
| 100+66800+600
---
(129.42)
| 9.50e+01
---
(failed)
| Exam. 17 lp_e226
---
(n = 695)
| 55
---
(0.65)
| 6.15e-07
---
| 14+953+17
---
(1.14)
| 1.45e-10
---
| 100+169800+1100
---
(34.97)
| 3.32e+01
---
(failed)
| Exam. 18 lp_etamacro
---
(n = 1216)
| 51
---
(2.46)
| 8.41e-07
---
| 7+640+18
---
(6.08)
| 5.85e-07
---
| 2+120+9
---
(1.17)
| 1.36e-08
---
| Exam. 19 lp_fffff800
---
(n = 1552)
| 70
---
(5.76)
| 6.93e-07
---
| 14+2172+21
---
(15.77)
| 2.09e-08
---
| 100+99000+700
---
(58.69)
| 1.09e+05
---
(failed)
| Exam. 20 lp_finnis
---
(n = 1561)
| 46
---
(3.87)
| 9.59e-07
---
| 10+1583+18
---
(14.41)
| 1.92e-07
---
| 100+24400+500
---
(28.81)
| 5.84e+01
---
(failed)
| Exam. 21 lp_fit1d
---
(n = 1073)
| 54
---
(2.09)
| 8.83e-07
---
| 8+681+13
---
(0.53)
| 3.46e-10
---
| 1+152+4
---
(0.06)
| 5.64e-09
---
| Exam. 22 lp_ganges
---
(n = 3015)
| 40
---
(14.89)
| 5.23e-07
---
| 8+1910+28
---
(111.67)
| 6.15e-08
---
| 82+104468+1476
---
(1003.45)
| 9.53e+00
---
(failed)
| Exam. 23 lp_gfrd_pnc
---
(n = 1776)
| 64
---
(7.15)
| 6.02e-07
---
| 14+1694+20
---
(40.06)
| 2.69e-09
---
| 100+6800+300
---
(14.97)
| 2.16e+03
---
(failed)
| Exam. 24 lp_grow15
---
(n = 945)
| 42
---
(1.17)
| 7.34e-07
---
| 8+856+10
---
(2.91)
| 9.04e-13
---
| 1+1+11
---
(1.02)
| 1.44e-09
---
| Exam. 25 lp_grow22
---
(n = 1386)
| 41
---
(2.65)
| 9.94e-07
---
| 9+1535+11
---
(12.94)
| 4.57e-08
---
| 100+2600+400
---
(18.05)
| 1.68e+00
---
(failed)
| Exam. 26 lp_lotfi
---
(n = 519)
| 46
---
(0.29)
| 6.47e-07
---
| 11+408+18
---
(0.48)
| 3.33e-11
---
| 100+47400+500
---
(13.56)
| 4.44e-02
---
(failed)
| Exam. 27 lp_maros
---
(n = 2812)
| 69
---
(21.44)
| 8.93e-07
---
| 15+4178+17
---
(193.91)
| 1.04e-08
---
| 100+14400+300
---
(35.28)
| 2.35e+04
---
(failed)

Table 5: Numerical results of small-scale dense LCPs (no. 28-53).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 28 lp_modszk1
---
(n = 2307)
| 46
---
(9.37)
| 8.17e-07
---
| 13+3281+16
---
(89.69)
| 4.94e-07
---
| 100+34000+500
---
(33.89)
| 9.03e+00
---
(failed)
| Exam. 29 lp_perold
---
(n = 2131)
| 60
---
(10.21)
| 7.57e-07
---
| 16+2738+24
---
(66.41)
| 1.33e-09
---
| 100+4200+400
---
(16.06)
| 2.88e+04
---
(failed)
| Exam. 30 lp_pilot_ja
---
(n = 3207)
| 76
---
(32.29)
| 9.82e-07
---
| 18+4063+29
---
(272.11)
| 7.42e-07
---
| 100+71396+779
---
(233.45)
| 5.77e+06
---
(failed)
| Exam. 31 lp_qap8
---
(n = 2544)
| 42
---
(10.69)
| 8.04e-07
---
| 12+3423+14
---
(116.20)
| 1.81e-07
---
| 100+23000+400
---
(29.77)
| 1.18e+01
---
(failed)
| Exam. 32 lp_lp_recipe
---
(n = 295)
| 53
---
(0.09)
| 7.61e-07
---
| 10+240+15
---
(0.16)
| 6.28e-08
---
| 100+5144+595
---
(4.11)
| 1.91e+01
---
(failed)
| Exam. 33 lp_sc50a
---
(n = 128)
| 42
---
(0.02)
| 7.79e-07
---
| 5+92+7
---
(0.02)
| 2.25e-08
---
| 1+135+2
---
(0.02)
| 7.11e-14
---
| Exam. 34 lp_scagr7
---
(n = 314)
| 41
---
(0.07)
| 8.19e-07
---
| 12+391+14
---
(0.16)
| 1.24e-08
---
| 2+452+8
---
(0.11)
| 1.77e-13
---
| Exam. 35 lp_scagr25
---
(n = 1142)
| 41
---
(1.72)
| 8.73e-07
---
| 9+648+19
---
(3.58)
| 1.65e-07
---
| 100+53500+700
---
(33.66)
| 1.09e+01
---
(failed)
| Exam. 36 lp_scfxm1
---
(n = 930)
| 54
---
(1.37)
| 5.12e-07
---
| 12+1260+16
---
(2.67)
| 1.22e-07
---
| 100+102400+1462
---
(75.61)
| 2.23e+02
---
(failed)
| Exam. 37 lp_scfxm2
---
(n = 1860)
| 58
---
(7.08)
| 8.13e-07
---
| 15+2685+19
---
(29.86)
| 3.10e-10
---
| 100+104383+902
---
(94.95)
| 2.35e+02
---
(failed)
| Exam. 38 lp_scfxm3
---
(n = 2790)
| 61
---
(18.68)
| 9.23e-07
---
| 13+3864+17
---
(135.97)
| 2.04e-07
---
| 100+27200+500
---
(46.27)
| 3.89e+02
---
(failed)
| Exam. 39 lp_scorpion
---
(n = 854)
| 40
---
(0.79)
| 7.15e-07
---
| 13+1094+15
---
(2.64)
| 2.19-10
---
| 100+84572+800
---
(31.31)
| 1.38e+00
---
(failed)
| Exam. 40 lp_shell
---
(n = 2313)
| 50
---
(10.06)
| 6.00e-07
---
| 13+3371+15
---
(85.95)
| 4.39e-08
---
| 100+56406+497
---
(45.75)
| 1.17e+02
---
(failed)
| Exam. 41 lp_ship04l
---
(n = 2568)
| 42
---
(10.56)
| 5.64e-07
---
| 12+4243+15
---
(96.83)
| 1.66e-09
---
| 100+31600+500
---
(61.70)
| 1.31e+02
---
(failed)
| Exam. 42 lp_ship04s
---
(n = 1908)
| 42
---
(5.50)
| 7.15e-07
---
| 12+3170+15
---
(32.59)
| 4.05e-07
---
| 100+20000+500
---
(32.59)
| 5.69e+01
---
(failed)
| Exam. 43 lp_ship08s
---
(n = 3245)
| 43
---
(18.97)
| 9.32e-07
---
| 12+5907+15
---
(279.95)
| 6.91e-08
---
| 100+70700+700
---
(162.72)
| 2.38e+02
---
(failed)
| Exam. 44 lp_ship12s
---
(n = 4020)
| 42
---
(30.77)
| 6.94e-07
---
| 9+3805+21
---
(333.53)
| 2.36e-10
---
| 77+205665+1472
---
(1002.61)
| 1.97e+02
---
(failed)
| Exam. 45 lp_sierra
---
(n = 3962)
| 47
---
(33.55)
| 8.12e-07
---
| 12+4069+16
---
(319.80)
| 7.15e-09
---
| 100+5000+300
---
(35.28)
| 9.99e+04
---
(failed)
| Exam. 46 lp_standgub
---
(n = 1744)
| 53
---
(5.65)
| 6.12e-07
---
| 12+2064+16
---
(19.86)
| 3.20e-08
---
| 100+94960+792
---
(146.17)
| 2.44e+03
---
(failed)
| Exam. 47 lp_tuff
---
(n = 961)
| 61
---
(1.72)
| 5.46e-07
---
| 15+1404+18
---
(3.91)
| 1.50e-10
---
| 100+24300+939
---
(62.39)
| 7.78e+03
---
(failed)
| Exam. 48 lp_wood1p
---
(n = 2839)
| 64
---
(20.82)
| 6.11e-07
---
| 58+3749+64
---
(473.97)
| 1.45e-09
---
| 100+2700+300
---
(10.44)
| 7.27e+05
---
(failed)
| Exam. 49 lpi_box1
---
(n = 492)
| 30
---
(0.20)
| 6.72e-07
---
| 13+657+15
---
(0.56)
| 7.62e-08
---
| 2+2115+24
---
(0.66)
| 2.67e-15
---
| Exam. 50 lpi_cplex2
---
(n = 602)
| 43
---
(0.45)
| 5.29e-07
---
| 12+627+14
---
(0.95)
| 6.53e-07
---
| 100+46800+798
---
(25.13)
| 1.48e+02
---
(failed)
| Exam. 51 lpi_ex72a
---
(n = 412)
| 34
---
(0.12)
| 5.52e-07
---
| 13+570+15
---
(0.41)
| 5.36e-09
---
| 1+867+6
---
(0.14)
| 3.13e-14
---
| Exam. 52 lpi_ex73a
---
(n = 404)
| 33
---
(0.12)
| 9.00e-07
---
| 13+552+15
---
(0.38)
| 7.35e-09
---
| 100+6739+600
---
(9.94)
| 1.50e-03
---
(failed)
| Exam. 53 lpi_mondou2
---
(n = 916)
| 38
---
(0.94)
| 5.86e-07
---
| 11+773+13
---
(2.36)
| 6.92e-10
---
| 2+1151+9
---
(0.63)
| 1.55e-14
---

Table 6: Numerical results of large-scale dense LCPs (no. 54-73).

Problems | RPFMTr | PATH | MILES
---|---|---|---
| (n*n)
---
| steps
---
(time)
| Terr
---
| major+minor+grad
---
(time)
| Terr
---
| major+pivots+refactor
---
(time)
| Terr
---
| Exam. 54 lp_80bau3b
---
(n = 14323)
| 45
---
(1385.62)
| 5.38e-07
---
| 5+8283+8
---
(18049.72)
| 3.92e+01
---
(failed)
| 2+8245+54
---
(1432.94)
| 1.32e+03
---
(failed)
| Exam. 55 lp_bnl2
---
(n = 6810)
| 55
---
(153.40)
| 8.21e-07
---
| 11+8765+18
---
(4549.92)
| 3.93e-08
---
| 43+27624+298
---
(1009.19)
| 9.22e+01
---
(failed)
| Exam. 56 lp_cre_a
---
(n = 10764)
| 61
---
(583.43)
| 6.36e-07
---
| 4+3835+7
---
(18099.73)
| 7.38e+01
---
(failed)
| 138+49799+558
---
(1003.36)
| 1.51e+03
---
(failed)
| Exam. 57 lp_cre_c
---
(n = 9479)
| 55
---
(373.29)
| 8.33e-07
---
| 7+8042+10
---
(18008.97)
| 2.68e+00
---
(failed)
| 140+28754+599
---
(1004.88)
| 3.55e+03
---
(failed)
| Exam. 58 lp_cycle
---
(n = 5274)
| 62
---
(94.96)
| 5.76e-07
---
| 12+7861+17
---
(2322.98)
| 2.77e-07
---
| 100+15300+400
---
(148.55)
| 2.07e+03
---
(failed)
| Exam. 59 lp_d6cube
---
(n = 6599)
| 54
---
(141.97)
| 6.46e-07
---
| 10+6621+16
---
(725.33)
| 5.15e-07
---
| 80+66400+800
---
(1009.78)
| 2.19e+03
---
(failed)
| Exam. 60 lp_fit2d
---
(n = 10549)
| 59
---
(550.98)
| 7.06e-07
---
| 8+5337+13
---
(185.28)
| 2.82e-09
---
| 2+147+4
---
(1.11)
| 3.78e-10
---
| Exam. 61 lp_greenbea
---
(n = 7990)
| 50
---
(223.83)
| 8.66e-07
---
| 12+11911+16
---
(7877.61)
| 5.44e-09
---
| 100+15200+400
---
(373.86)
| 1.14e+02
---
(failed)
| Exam. 62 lp_greenbeb
---
(n = 7990)
| 50
---
(225.27)
| 9.43e-07
---
| 12+12033+16
---
(10248.31)
| 2.46e-08
---
| 100+17600+500
---
(575.45)
| 1.14e+02
---
(failed)
| Exam. 63 lp_ken_07
---
(n = 6028)
| 42
---
(96.95)
| 8.68e-07
---
| 8+4721+28
---
(2922.75)
| 1.13e-07
---
| 56+6944+448
---
(1008.24)
| 5.01e+01
---
(failed)
| Exam. 64 lp_pds_02
---
(n = 10669)
| 48
---
(437.80)
| 7.48e-07
---
| 8+8486+26
---
(15408.33)
| 1.82e-08
---
| 100+9000+300
---
(504.81)
| 1.86e+01
---
(failed)
| Exam. 65 lp_pilot87
---
(n = 8710)
| 58
---
(307.63)
| 5.54e-07
---
| 11+8828+23
---
(7239.50)
| 8.07e-10
---
| 56+48048+560
---
1014.17)
| 9.99e+02
---
(failed)
| Exam. 66 lp_qap12
---
(n = 12048)
| 48
---
(722.82)
| 8.45e-07
---
| 3+5789+6
---
(18016.27)
| 2.73e+01
---
(failed)
| 100+23400+400
---
(771.77)
| 2.14e+01
---
(failed)
| Exam. 67 lp_ship08l
---
(n = 5141)
| 47
---
(66.37)
| 9.12e-07
---
| 12+9419+15
---
(1537.80)
| 1.03-09
---
| 100+80800+900
---
(464.31)
| 4.34e+02
---
(failed)
| Exam. 68 lp_ship12l
---
(n = 6684)
| 47
---
(124.64)
| 7.67e-07
---
| 12+12138+15
---
(7256.50)
| 1.51-07
---
| 29+11948+145
---
(1001.00)
| 7.53e+02
---
(failed)
| Exam. 69 lp_truss
---
(n = 9806)
| 45
---
(339.85)
| 6.61e-07
---
| 10+6890+12
---
(10715.00)
| 2.30e-08
---
| 23+8685+201
---
(1022.70)
| 8.80e+00
---
(failed)
| Exam. 70 lp_woodw
---
(n = 9516)
| 71
---
(494.92)
| 7.61e-07
---
| 14+7490+23
---
(18080.28)
| 1.35e-05
---
| 100+15186+299
---
(340.14)
| 2.06e+02
---
(failed)
| Exam. 71 lpi_bgindy
---
(n = 13551)
| 53
---
(1297.80)
| 7.46e-07
---
| 1+663+6
---
(18181.77)
| 3.12e+02
---
(failed)
| 9+11664+81
---
(1025.18)
| 2.08e+02
---
(failed)
| Exam. 72 lpi_gran
---
(n = 5183)
| 60
---
(81.64)
| 5.32e-07
---
| 14+5542+23
---
(1450.42)
| 1.29e-08
---
| 500+145000+2500
---
(719.47)
| 2.357e+03
---
(failed)
| Exam. 73 lpi_greenbea
---
(n = 7989)
| 50
---
(211.05)
| 5.37e-07
---
| 12+12005+16
---
(12504.91)
| 7.46e-09
---
| 47+15546+227
---
(1006.56)
| 1.13e+02
---
(failed)

Figure 3: Iterations of RPFMTr, PATH and MILES for dense LCPs. Figure 4: CPU
time of RPFMTr, PATH and MILES for dense LCPs.

##  5 Conclusions

In this paper, we give the residual regularization path-following method with
the trust-region updating strategy (RPFMTr) for the linear complementarity
problem. The new residual regularization parameter improves the robustness of
the path-following method, in comparison to the traditional complementarity
regularization parameter. Meanwhile, we prove the global convergence of the
new method under the standard assumptions without the traditional assumption
condition of the priority to feasibility over complementarity. Numerical
results show that RPFMTr is a robust and efficient solver for the linear
complementarity problem, especially for the dense case. Furthermore, it is
more robust and faster than some state-of-the-art solvers such as PATH DF1995
; PATH ; FM2000 ; FM2022  and MILES Mathiesen1985 ; Rutherford1995 ;
Rutherford2022  (the built-in subroutines of the GAMS v28.2 (2019) environment
GAMS ). The computational time of RPFMTr is about 1/3 to 1/10 of that of PATH
for the dense linear complementarity problem. Therefore, RPFMTr is an
alternating solver for the linear complementarity problem and worth exploring
further for the nonlinear complementarity problem.

## Acknowledgments

This work was supported in part by Grant 61876199 from National Natural
Science Foundation of China, Grant YBWL2011085 from Huawei Technologies Co.,
Ltd., and Grant YJCB2011003HI from the Innovation Research Program of Huawei
Technologies Co., Ltd..

Conflicts of interest / Competing interests: Not applicable.

Availability of data and material (data transparency): If it is requested, we
will provide the test data.

Code availability (software application or custom code): If it is requested,
we will provide the code.

## References

  * (1) AMPL, A Mathematical Programming Language, http://www.ampl.com, 2022.
  * (2) E. L. Allgower and K. Georg, _Introduction to Numerical Continuation Methods_ , SIAM, Philadelphia, 2003.
  * (3) U. M. Ascher and L. R. Petzold, _Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations_ , SIAM, Philadelphia, 1998.
  * (4) O. Axelsson and S. Sysala, _Continuation Newton methods_ , Comput. Math. Appl. 70 (2015), 2621-2637.
  * (5) F. H. Branin, _Widely convergent method for finding multiple solutions of simultaneous nonlinear equations_ , IBM J. Res. Dev. 16 (1972), 504-521.
  * (6) J. C. Butcher and Z. Jackiewicz, _Construction of high order diagonally implicit multistage integration methods for ordinary differential equations_ , Appl. Numer. Math. 27 (1998), 1-12.
  * (7) K. E. Brenan, S. L. Campbell, and L. R. Petzold, _Numerical solution of initial-value problems in differential-algebraic equations_ , SIAM, Philadelphia, 1996\.
  * (8) R. W. Cottle, J.-S. Pang and R. E. Stone, _The Linear Complementarity Problem_ , SIAM, Philadelphia, 2009.
  * (9) A. R. Conn, N. Gould, and Ph. L. Toint, _Trust-Region Methods_ , SIAM, Philadelphia, 2000.
  * (10) M. T. Chu and M. M. Lin, _Dynamical system characterization of the central path and its variants- a vevisit_ , SIAM J. Appl. Dyn. Syst. 10 (2011), 887-905.
  * (11) D. F. Davidenko, _On a new method of numerical solution of systems of nonlinear equations_ (in Russian), Dokl. Akad. Nauk SSSR 88 (1953), 601–602.
  * (12) P. Deuflhard, _Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms_ , Springer, Berlin, 2004.
  * (13) P. Deuflhard, H. J. Pesch, and P. Rentrop, _A modified continuation method for the numerical solution of nonlinear two-point boundary value problems by shooting techniques_ , Numer. Math. 26 (1975), 327-343.
  * (14) S. P. Dirkse, _Robust Solution of Mixed Complementarity Problems_ , PhD thesis, Computer Sciences Department, University of Wisconsin, Madison, Wisconsin, 1994.
  * (15) S. Dirkse and M. C. Ferris, _The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems_ , Optim. Methods Softw. 5 (1995), 123-156.
  * (16) E. J. Doedel, _Lecture notes in numerical analysis of nonlinear equations_ , in _Numerical Continuation Methods for Dynamical Systems_ , B. Krauskopf, H.M. Osinga, and J. Galán-Vioque, eds., Springer, Berlin, 2007, pp. 1-50.
  * (17) K. Erleben, _Velocity-based shock propagation for multibody dynamics animation_ , ACM Trans. Graph. 26 (2):12 (2007), https://doi.org/10.1145/1243980.1243986.1.
  * (18) M. C. Ferris and T. S. Munson, _Complementarity problems in GAMS and the PATH solver_ , J. Econ. Dyn. Control 24 (2000), 165-188.
  * (19) M. C. Ferris and T. S. Munson, _The manual of PATH_ , GAMS Corporation, https://www.gams.com/latest/docs/S_PATH.html, 2022.
  * (20) A. Fischer, _A special Newton-type optimization method_ , Optimization 24 (1992), 269-284.
  * (21) GAMS v28.2, GAMS Corporation, https://www.gams.com/, 2019.
  * (22) A. Gana, _Studies in the Complementarity Problem_ , Ph.D. dissertation, Department of Industrial and Operations Engineering, University of Michigan, Ann Arbor, MI, 1982.
  * (23) P. Gill, W. Murray, M. A. Saunders, M. H. Wright, _Mamtaining LU factors of a general sparse matrix_ , Linear Algebra Appl. 88-89 (1991), 239-270.
  * (24) G. H. Golub and C. F. Van Loan, _Matrix Computation_ , 4th ed., The John Hopkins University Press, Baltimore, 2013.
  * (25) T. H. Gronwall, _Note on the derivatives with respect to a parameter of the solutions of a system of differential equations_ , Ann. of Math. 20 (1919), pp. 292-296, JFM 47.0399.02 (https://zbmath.org/?format=complete&q=an:47.0399.02), JSTOR 1967124 (https://www.jstor.org/stable/1967124), MR 1502565 (https://www.ams.org/mathscinet-getitem?mr=1502565).
  * (26) E. Hairer and G. Wanner, _Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems_ , 2nd ed., Springer, Berlin, 1996.
  * (27) D. J. Higham, _Trust region algorithms and timestep selection_ , SIAM J. Numer. Anal. 37 (1999), 194-210.
  * (28) P. C. Hansen, _Regularization Tools: A MATLAB package for analysis and solution of discrete ill posed problems_. Numer. Algorithms 6 (1994), 23-29.
  * (29) Z. Jackiewicz, S. Tracogna, _A general class of two-step Runge-Kutta methods for ordinary differential equations_ , SIAM J. Numer. Anal. 32 (1995), 1390-1427.
  * (30) C. T. Kelley, Solving Nonlinear Equations with Newton’s Method, SIAM, Philadelphia, 2003.
  * (31) C. T. Kelley, _Numerical methods for nonlinear equations_ , Acta Numer. 27 (2018), 207-287.
  * (32) C. E. Lemke and J. T. Howson, Jr. _Equilibrium points of bimatrix games_ , Journal of the Society for Industrial and Applied Mathematics 12 (1964), 413-423.
  * (33) X.-L. Luo, _A second-order pseudo-transient method for steady-state problems_ , Appl. Math. Comput. 216 (2010), 1752-1762.
  * (34) X.-L. Luo, _A dynamical method of DAEs for the smallest eigenvalue problem_ , J. Comput. Sci. 3 (2012), 113-119.
  * (35) X.-L. Luo, H. Xiao, J.-H. Lv and S. Zhang, _Explicit pseudo-transient continuation and the trust-region updating strategy for unconstrained optimization_ , Appl. Numer. Math. 165 (2021), 290-302, available at http://doi.org/10.1016/j.apnum.2021.02.019.
  * (36) X.-L. Luo, J.-H. Lv and G. Sun, _Continuation method with the trusty time-stepping scheme for linearly constrained optimization with noisy data_ , Optim. Eng. 23 (2022), 329-360, available at http://doi.org/10.1007/s11081-020-09590-z.
  * (37) X.-L. Luo, H. Xiao and J.-H. Lv, _Continuation Newton methods with the residual trust-region time-stepping scheme for nonlinear equations_ , Numer. Algorithms 89 (2022), 223-247, available at http://doi.org/10.1007/s11075-021-01112-x.
  * (38) X.-L. Luo and H. Xiao, _Generalized continuation Newton methods and the trust-region updating strategy for the underdetermined system_ , J. Sci. Comput. 88 (2021), article 56, published online at http://doi.org/10.1007/s10915-021-01566-0, pp. 1-22, July 13, 2021.
  * (39) X.-L. Luo and Y.-Y. Yao, _Primal-dual path-following methods and the trust-region updating strategy for linear programming with noisy data_ , J. Comp. Math. 40 (2022), 760-780, published online at http://doi.org/10.4208/jcm.2101-m2020-0173, or available at http://arxiv.org/abs/2006.07568.
  * (40) X.-L. Luo and H. Xiao, _Continuation Newton methods with deflation techniques and quasi-genetic evolution for global optimization problems_ , arXiv preprint available at http://arxiv.org/abs/2107.13864, or Research Square preprint available at https://doi.org/10.21203/rs.3.rs-1102775/v1, July 30, 2021. Software available at https://teacher.bupt.edu.cn/luoxinlong/zh_CN/zzcg/41406/list/index.htm.
  * (41) X.-L. Luo and H. Xiao, _The regularization continuation method with an adaptive time step control for linearly constrained optimization problems_ , Appl. Numer. Math. 181 (2022), 255-276, available at https://doi.org/10.1016/j.apnum.2022.06.008.
  * (42) MATLAB v9.8.0 (R2020a), The MathWorks Inc., http://www.mathworks.com, 2020.
  * (43) L. Mathiesen, _Computation of economic equilibria by a sequence of linear complementarity problems_ , Math. Program. 23 (1985), 144-162.
  * (44) T. F. Rutherford, _The manual of MILES_ , GAMS Corporation, https://www.gams.com/latest/docs/S_MILES.html, 2022.
  * (45) The NETLIB collection of linear programming problems, available at http://www.netlib.org.
  * (46) J. Nocedal and S. J. Wright, _Numerical Optimization_ , Springer, Berlin, 1999.
  * (47) J. M. Ortega and W. C. Rheinboldt, _Iteration Solution of Nonlinear Equations in Several Variables_ , SIAM, Philadelphia, 2000.
  * (48) M. C. Ferris et al., _PATH solver 5.0.00_ , software available at https://pages.cs.wisc.edu/~ferris/path.html, 2019.
  * (49) M. J. D. Powell, _Convergence properties of a class of minimization algorithms_ , in _Nonlinear Programming 2_ , O. L. Mangasarian, R. R. Meyer, and S. M. Robinson, eds., Academic Press, New York, 1975, pp. 1-27.
  * (50) T. F. Rutherford, _Extension of GAMS for complementarity problems arising in applied economic analysis_ , J. Econ. Dyn. Control 19 (1995), 1299-1324.
  * (51) L. F. Shampine, I. Gladwell, and S. Thompson, _Solving ODEs with MATLAB_ , Cambridge University Press, Cambridge, 2003.
  * (52) K. Tanabe, _Continuous Newton-Raphson method for solving an underdetermined system of nonlinear equations_ , Nonlinear Anal. 3 (1979), 495-503.
  * (53) K. Tanabe, _Centered Newton method for mathematical programming_ , in _System Modeling and Optimization: Proceedings of the 13th IFIP conference_ , vol. 113 of Lecture Notes in Control and Information Systems, Springer, Berlin, 1988, pp. 197-206.
  * (54) V. Venkateswaran, _An algorithm for the linear complementarity problem with a $P_{0}$-matrix_, SIAM J. Matrix Anal. Appl. 14 (1993), 967-977.
  * (55) S. J. Wright, _An infeasible-interior-point algorithm for linear complementarity problems_ , Math. Program. 67 (1994), 29-51.
  * (56) S. J. Wright, _Primal-dual Interior Point Methods_ , SIAM, Philadelphia, 1997.
  * (57) Q. Xu, _On the New Linear Programming Algorithms-New Sequential Penalty Function Method and Two Point Barrier Function Method (in Chinese)_ , Ph.D. thesis, Institute of Nuclear Technology, Tsinghua University, Beijing, China, 1991.
  * (58) Y. X. Yuan, _Recent advances in trust region algorithms_ , Math. Program. 151 (2015), 249-281.
  * (59) Y. Zhang, _On the convergence of a class of infeasible interior-point methods for the horizontal linear complementarity problem_ , SIAM J. Optim. 4 (1994), 208-227.