A lower bound on the iterative complexity of the Harker and Pang globalization technique of the Newton-min algorithm for solving the linear complementarity problem

Abstract

The plain Newton-min algorithm for solving the linear complementarity problem (LCP) “$0⩽x\perp (Mx+q)⩾0$” can be viewed as an instance of the plain semismooth Newton method on the equational version “$min(x,Mx+q)=0$” of the problem. This algorithm converges for any q when M is an $\mathbf{M}$-matrix, but not when it is a $\mathbf{P}$-matrix. When convergence occurs, it is often very fast (in at most n iterations for an $\mathbf{M}$-matrix, where n is the number of variables, but often much faster in practice). In 1990, Harker and Pang proposed to improve the convergence ability of this algorithm by introducing a stepsize along the Newton-min direction that results in a jump over at least one of the encountered kinks of the min-function, in order to avoid its points of nondifferentiability. This paper shows that, for the Fathi problem (an LCP with a positive definite symmetric matrix M, hence a $\mathbf{P}$-matrix), an algorithmic scheme, including the algorithm of Harker and Pang, may require n iterations to converge, depending on the starting point.

Introduction

Let $n⩾1$ be an integer, $M\in {\mathbb{R}}^{n\times n}$ be a real matrix, $q\in {\mathbb{R}}^n$ be a real vector, and $[1:n]:=\{1,\dots ,n\}$ be the set of the first n positive integers. The linear complementarity problem (LCP) consists in searching a vector $x\in {\mathbb{R}}^n$ such that

$$\begin{array}{c}\hfill 0⩽x\perp (Mx+q)⩾0.\end{array}$$

This means that the sought x must satisfy $x⩾0$, $Mx+q⩾0$ (vectorial inequalities must be understood componentwise), and $x^{\mathsf{T}}(Mx+q)=0$ (the exponent “${}^{\mathsf{T}}$” is used to denote matrix transposition). The problem has a combinatorial aspect, which lies in this last equation, since, by the nonnegativity of x and $Mx+q$, it amounts to the set of n complementarity conditions $x_i{(Mx+q)}_i=0$ for all indices $i\in [1:n]$. The term complementarity comes from the fact that, for all $i\in [1:n]$, either $x_i$ or ${(Mx+q)}_i$ must vanish; these conditions may be realized in $2^n$ different ways. Actually, the problem of determining whether a particular instance of the LCP has a solution is strongly NP-complete (Chung 1989), and NP-complete for a ${\mathbf{P}}_{\mathbf{0}}$-matrix (i.e., when M has nonegative principal minors) (Kojima et al. 1991a).

Let $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ be the min-function associated with the LCP (1.1), which is the function that takes at $x\in {\mathbb{R}}^n$ the value

$$\begin{array}{c}\hfill F(x)=min(x,Mx+q).\end{array}$$

The Newton-min algorithm can be viewed as an instance of the semismooth Newton method (Qi and Sun 1993; Kojima and Shindo 1986) to solve the equational equivalent form of (1.1) (Kostreva 1976; Mangasarian 1977; Cottle et al. 2009) that reads $F(x)=0$. To write compactly the algorithm, it is useful to introduce, for $I\subseteq [1:n]$ and its complement $A:=[1:n]\backslash I$, the point $x^{(I)}$ defined by

$$\begin{array}{c}\hfill x_A^{(I)}=0\text{and}{(M{x}^{(I)}+q)}_I=0.\end{array}$$

This point is well defined when M is nondegenerate, meaning that the principal minors of M do not vanish. The plain Newton-min algorithm computes the next iterate by

$$\begin{array}{c}\hfill \widehat{x}:=x^{(\mathcal{S}(x))},\end{array}$$

where the index selector$\mathcal{S}:{\mathbb{R}}^n⊸[1:n]$ is the multifunction defined at $x\in {\mathbb{R}}^n$ by

$$\begin{array}{c}\hfill \mathcal{S}(x):=\{i\in [1:n]:x_i>{(Mx+q)}_i\}.\end{array}$$

In some versions of the algorithm, $\mathcal{S}(x)$ also contains some or all the indices in $\{i\in [1:n]:x_i={(Mx+q)}_i\}$. See paragraph 7 of the introduction of Ben Gharbia and Gilbert (2012) for more details on the origin of this algorithm and a discussion on the contributions from (Chandrasekaran 1970; Kojima and Shindo 1986; Harker and Pang 1990; Fischer and Kanzow 1996; Bergounioux et al. 1999, 2000; Hintermüller et al. 2003; Kanzow 2004). When the current iterate $x\in {\mathbb{R}}^n$ is not on a kink of F, like in this paper, the Newton-min algorithm is identical to the Newton method to find a zero of F, which is then well defined.

Even though the Newton-min algorithm uses no globalization technique, like line searches or trust regions (Bonnans et al. 2006; Conn et al. 2000), it may converge globally, i.e., from any starting point. This is due to the very particular piecewise linearity of F. For example, global convergence occurs, whatever q is, when M is an $\mathbf{M}$-matrix (Aganagić 1984), which is a $\mathbf{P}$-matrix (i.e., with positive principal minors) with nonpositive off-diagonal elements. It also occurs when M is close enough to an $\mathbf{M}$-matrix (Hintermüller et al. 2003). However, this global convergence property does not extend up to the larger class of $\mathbf{P}$-matrices (Ben Gharbia and Gilbert 2012, 2013; Curtis et al. 2015; Ben Gharbia and Gilbert 2019). This is unfortunate for the Newton-min algorithm, since $\mathbf{P}$-matrices are exactly those ensuring the existence and uniqueness of the solution to the LCP, whatever q is (Samelson et al. 1958; Cottle et al. 2009).

A natural idea to enlarge the class of matrices, for which the global convergence of the Newton-min algorithm can be guaranteed, is to introduce a line search on the associated least-square merit function, which is the function $\Theta :{\mathbb{R}}^n\to \mathbb{R}$ defined at $x\in {\mathbb{R}}^n$ by

$$\begin{array}{c}\hfill \Theta (x)=\frac{1}{2}{\Vert F(x)\Vert }^2,\end{array}$$

where $\Vert ·\Vert$ denotes the Euclidean norm. This least-square function is natural, since it has been used, often with success, for globalizing the Newton method when the function F is smooth (Dennis and Schnabel 1983; Conn et al. 2000; Bonnans et al. 2006; Izmailov and Solodov 2014). In the presence of nonsmoothness of F, like here, this technique is more difficult to implement since the Newton-min direction $d:=\widehat{x}-x$ may not be a descent direction of $\Theta$ at a kink of F (Dussault et al. 2019) (this fact was already observed during the Ph.D. thesis of Ben Gharbia (2012, example 5.8). To overcome this difficulty, Harker and Pang (1990, p. 275) proposed a method named the Modified Damped-Newton Algorithm, which consists in taking for the next iterate the point

$$\begin{array}{c}\hfill x^{+}:=x+({\check{\alpha }}_1+{\epsilon }_{\mathrm{H}\mathrm{P}})d,\end{array}$$

where ${\check{\alpha }}_1>0$ is a stepsize so that $x+{\check{\alpha }}_{1}d$ is on the first kink of F encountered along d from x, and ${\epsilon }_{\mathrm{H}\mathrm{P}}>0$ is a number such that the new iterate $x^{+}$ is not on a kink of F and ensures a sufficient decrease in the least-square merit function $\Theta$. Consequently, this algorithm avoids the points of nondifferentiability of F, generates descent directions of $\Theta$, and forces $\Theta$ to decrease sufficiently at each iteration. To the best of our knowledge, the only convergence result for any line search algorithm using a semismooth Newton direction uses the assumption that ${lim\; inf}_k{\alpha }_k>0$ (Pang 1990), which is a very weak result since this strong assumption relates to the algorithm products rather than the problem’s data.

In this research field, sparing of theoretical results, this paper provides the value n as a worse case lower bound on the number of iterations of the Harker and Pang algorithm when the extra stepsize ${\epsilon }_{\mathrm{H}\mathrm{P}}$ is taken sufficiently small, which is allowed by the description of the method given in Harker and Pang (1990, p. 275). This lower bound is obtained on the Fathi problem for a set of starting points, including the one of Fathi (1979), which is zero. To extend the applicability of this result, we describe an algorithmic scheme, for which this worse case lower bound is valid, a scheme that includes the Harker and Pang algorithm for sufficiently small positive ${\epsilon }_{\mathrm{H}\mathrm{P}}$. In this scheme, the iterates avoid the kinks of F and the stepsizes are chosen arbitrarily between the first two break-stepsizes ${\check{\alpha }}_1$ and ${\check{\alpha }}_2$ (to be defined). Now, on many practical problems, an algorithm using the Newton-min direction and a stepsize that is not forced to be in $({\check{\alpha }}_1,{\check{\alpha }}_2)$ usually finds a solution in much fewer iterations than n; in the experiments of Dussault et al. (2019), it is not uncommon to encounter LCPs having up to ${10}^5$ variables that are solved in fewer than 10 iterations. Nevertheless, the Fathi problem remains a difficult instance of LCP for this family of methods, independently of the chosen stepsizes. To illustrate this, we show in the numerical experiment section that, surprisingly, doing exact line searches hardly modifies the iteration counter. Finally, this worse case lower bound and the numerical experiments of Sect. 5 suggest us that it is unlikely that the improvement in the Newton-min algorithm can lie in a better determination of the stepsizes. This observation paves the way for the proposals made in Dussault et al. (2019).

To conclude Introduction, let us mention that there are a large number of contributions related to the complexity of algorithms for solving the LCP. Most of them are related to interior point methods, and it is out of the scope of this paper to review them (they can be found by looking at those citing one of the first accounts on the subject, which is Kojima et al. 1991a, b). Other approaches are sometimes qualified as noninterior path-following/continuation methods and are based on the smoothing of equational versions of the LCP: the function $(a,b)\in {\mathbb{R}}^2\mapsto a+b-{[{(a-b)}^2+4{\mu }^2]}^{1/2}$ is considered in Chen and Harker (1993), Kanzow (1996), and the smooth Fischer–Burmeister function $(a,b)\in {\mathbb{R}}^2\mapsto a+b-{[a^2+b^2+2{\mu }^2]}^{1/2}$ is used in Kanzow (1996). The complexity of these approaches has been studied in Burke and Xu (1999), Hotta et al. (2000), Burke and Xu (2002), for instance.

The paper is structured as follows: The algorithmic scheme, for which the lower bound on the iterative complexity is obtained, is presented in Sect. 2. The Fathi problem and two properties of its matrix are given in Sect. 3. The iterative complexity result is proved in Sect. 4. Finally, some numerical experiments are reported in Sect. 5 and the paper ends with Conclusion Sect. 6.

Notation For the $n\times n$ matrix M and index sets I and $J\subseteq [1:n]$, we denote by $M_{\mathit{IJ}}$ the submatrix of M formed of its elements with row indices in I and column indices in J. We also define $M_{I:}:=M_{I[1:n]}$ and $M_{\mathit{II}}^{-1}:={(M_{\mathit{II}})}^{-1}$.