1 Introduction

The standard LCP is a classical problem of the optimization theory. Given a matrix \(M\in \mathbb {R}^{n\times n}\) and a vector \(\mathbf {q}\in \mathbb {R}^n\) our task is to find a vector pair \((\mathbf {x},\mathbf {s})\in \mathbb {R}^n\times \mathbb {R}^n\) satisfying

$$\begin{aligned} -M\mathbf {x}+ \mathbf {s}= \mathbf {q}, \qquad \mathbf {x}\mathbf {s}= \mathbf {0}, \qquad \mathbf {x}, \mathbf {s}\ge 0, \end{aligned}$$
(1)

where \(\mathbf {x}\mathbf {s}\) stands for the Hadamard-product (or componentwise product) of two vectors.

There are different versions of LCPs, such as horizontal-LCP, mixed LCP, geometric LCP. Anitescu et al. [4] proved that the above versions of LCP can be transformed into the standard LCP. The primal-dual pair of linear or convex quadratic optimization problems can be reformulated as an LCP. Furthermore, several important practical problems lead to LCP. For a comprehensive treatment of LCPs theory and applications, we refer the reader to the monographs of Cottle et al. [8] and Kojima et al. [35].

Cottle et al. [9] introduced the concept of sufficient matrix. A matrix \(M\in \mathbb {R}^{n\times n}\) is called column sufficient if for all \(\mathbf {x}\in \mathbb {R}^n\) we have

\(x_i(M\mathbf {x})_i\le 0\;\) for all \(i\in \mathcal {I}\; \Longrightarrow \; x_i(M\mathbf {x})_i= 0\;\) for all \(i\in \mathcal {I},\)

where \(\mathcal {I} = \{1,2,\ldots , n\}\) is an index set. Moreover, M is row sufficient matrix if \(M^T\) is column sufficient. A matrix is called sufficient if it is both row and column sufficient.

The sufficient matrix class can be defined in a different way through the \(\mathcal {P}_*\) matrices. The concept of \(\mathcal {P}_*\)-matrices was introduced by Kojima et al. [35]. Let \(\kappa \ge 0\) be a nonnegative number. We call the matrix \(M\in \mathbb {R}^{n\times n}\) a \(\mathcal {P}_*(\kappa )\)-matrix if

$$\begin{aligned} (1+4\kappa )\sum _{i\in \mathcal {I}_+(\mathbf {x})}x_i(M\mathbf {x})_i + \sum _{i\in \mathcal {I}_-(\mathbf {x})}x_i(M\mathbf {x})_i\ge 0 \end{aligned}$$
(2)

holds for any \(\mathbf {x}\in \mathbb {R}^n\), where \(\mathcal {I}_+(\mathbf {x})=\{i\in \mathcal {I}\;|\; x_i(M\mathbf {x})_i>0\}\) and \(\mathcal {I}_-(\mathbf {x})=\{i\in \mathcal {I}\;|\;x_i(M\mathbf {x})_i<0\}\).

We may define the set of \(\mathcal {P_*}\)-matrices as the union of sets of \(\mathcal {P}_*(\kappa )\)-matrices in the following way:

$$\begin{aligned} \mathcal {P_*}:=\bigcup \limits _{\kappa \ge 0}\mathcal {P}_*(\kappa ). \end{aligned}$$

Kojima et al. [35] proved that a \(\mathcal {P_*}\)-matrix is also column sufficient. Subsequently Guu and Cottle [28] proved that a \(\mathcal {P_*}\)-matrix is also row sufficient. Väliaho [42] proved the reverse inclusion, therefore the class of sufficient matrices and that of \(\mathcal {P_*}\)-matrices are the same. Considering the definition of \(\mathcal {P}_*\)-matrix class, a matrix is sufficient if and only if it is a \(\mathcal {P}_*(\kappa )\)-matrix for some \(\kappa \ge 0\). The smallest \(\kappa \) with this property is called the handicap of the matrix [43].

If the matrix of the LCP has the \(P_*(\kappa )\) property the LCP is called \(P_*(\kappa )\)-LCP. There are several approaches to solve \(\mathcal {P}_*(\kappa )\)-LCPs. Among them IPAs received much more attention than pivot algorithms [10,11,12, 22, 26, 27]. Since the primal-dual linear optimization (LO) problem pair can be formulated as a special class of LCPs, it is—still—an interesting question whether any IPA could be extended to solve efficiently \(P_*(\kappa )\)-LCPs. First Kojima et al. [35] followed this way and generalized a primal–dual path-following IPA to \(P_*(\kappa )\)-LCP. Their algorithm has \(O((1+\kappa )\sqrt{n}L)\) iteration complexity, where L is the input size of the problem with integer data. This complexity is still the best one for solving \(P_*(\kappa )\)-LCP. Since then a lot of known IPAs introduced for solving LO problems were extended to \(P_*(\kappa )\)-LCPs [32].

In 2002 Darvay offered a new technique for finding search directions of IPAs for LO problems [13,14,15]. Main idea of Darvay’s approach is very simple: perform a one-to-one transformation on the nonlinear equation of the central path problem and apply Newton’s method to the transformed system. In [13, 15] the transformation is achieved by the square root function. Since then, this technique has been extended to different types of optimization problems, for LCP [1, 3, 5, 6, 34, 39, 44, 46], for convex quadratic programming (CQP) [2], for semidefinite programing (SDP) [7, 45], for second-order cone programming (SOCP) [46] and symmetric optimization (SO) [33, 47].

Let us define the transformation that Darvay introduced and analysed for the central path of the LO problem as a continuously differentiable monotone increasing function

$$\begin{aligned} \varphi : (\xi , \infty )\rightarrow \mathbb {R}, \end{aligned}$$

where \(0 \le \xi < 1\) is a fixed value, which depends on the function \(\varphi \). In the view of [13, 18], we may say that the first interesting transformation of the nonlinear equation of the central path is

$$\begin{aligned} \varphi (t) = t \end{aligned}$$
(3)

that has been used widely in the area of LO (for details see Roos et al. [40]). First monotone function analysed by Darvay [13] was

$$\begin{aligned} \varphi (t) = \sqrt{t}. \end{aligned}$$
(4)

Recently, Darvay et al. [18] used

$$\begin{aligned} \varphi (t) = t - \sqrt{t}. \end{aligned}$$
(5)

Naturally, for functions used in [13, 18], the transformed central path equation and hence the right-hand side of the Newton system becomes more complicated, with considerable effects on the computations and the complexity analysis of the algorithm.

Kheirfam [33] and Pirhaji et al. [39] investigated the function \(\varphi (t) = \frac{\sqrt{t}}{2(1+\sqrt{t})}\) for \(\mathcal {P}_*(\kappa )\)-LCPs and circular cone optimization, respectively.

Due to the fact that function \(\varphi \) used by Darvay and others need to be a continuously differentiable monotone increasing function, we refer to this transformation technique as algebraic equivalent transformation of the central path equation.

Motivated by the mentioned works we propose a full-Newton step IPA for solving the \(\mathcal {P}_*(\kappa )\)-LCP based on the direction introduced by Darvay et al. [18]. We analyze the algorithm and obtain the complexity bound, which coincides with the best known result for \(\mathcal {P}_*(\kappa )\)-LCP, in this way showing that the direction presented in [18] can be used in IPAs for \(\mathcal {P}_*(\kappa )\)-LCPs. We mention that our approach is significantly different from [12], where criss-cross algorithms are analyzed. There are also important differences between our algorithm and the one presented in [32], which focuses on general LCPs, because we apply the AET technique, which is not considered in [32].

The outline of the rest of the paper is as follows. In Sect. 2 we recall some basic concepts and describe Darvay’s method for LCPs. In Sect. 3, we present a full-Newton step feasible IPA for \(\mathcal {P}_* (\kappa )-LCP\)s using the direction induced by the function given in (5) and state two theorems about the convergence and complexity of the algorithm. In Sect. 4 we provide some numerical results to illustrate performance of the proposed IPA and compare it, with other variants. Finally, we give some conclusions in Sect. 5.

The notations used in this paper are presented in “Appendix A.1”. Moreover, some important well-known results used in this paper are summerized in “Appendix A.2”.

In the next section we consider the AET technique presented by Darvay [13, 15] in order to obtain new search directions.

2 Newton system for the transformed central path equation

Let \(0 \le \xi < 1\) and \(\varphi :(\xi , \infty )\rightarrow \mathbb {R}\) be a continuously differentiable and invertible function. For a vector \(\mathbf {v}\in \mathbb {R}^n\) we define \(\varphi (\mathbf {v})\) in the following way:

$$\begin{aligned} \varphi (\mathbf {v})\in \mathbb {R}^n,\quad [\varphi (\mathbf {v})]_i := \varphi (v_i). \end{aligned}$$

We use an AET on the centrality condition of (26) getting the following equivalent system

$$\begin{aligned} -M\mathbf {x}+ \mathbf {s}= \mathbf {q}, \qquad \varphi \left( \frac{\mathbf {x}\mathbf {s}}{\mu }\right) = \varphi (\mathbf {e}), \qquad \mathbf {x}, \mathbf {s}\ge 0, \end{aligned}$$
(6)

where \(\mu >0\). Furthermore, we apply Newton’s method in order to approximate the solutions of (6). A simple calculus yields the following system

$$\begin{aligned} \begin{aligned} -M\varDelta \mathbf {x}+ \varDelta \mathbf {s}&= \mathbf {0}\\ \mathbf {s}\varDelta \mathbf {x}+ \mathbf {x}\varDelta \mathbf {s}&= \mu \;\frac{\varphi (\mathbf {e}) - \varphi \left( \frac{\mathbf {x}\mathbf {s}}{\mu }\right) }{\varphi '\left( \frac{\mathbf {x}\mathbf {s}}{\mu }\right) }. \end{aligned} \end{aligned}$$
(7)

From Lemma 3 the above transformed Newton system has a unique solution if M is a \(\mathcal {P}_{*}(\kappa )\)-matrix.

This type of search direction was introduced by Zsolt Darvay for linear optimization. In [15] he used the function \(\varphi (t) = \sqrt{t}\) and later in [18] the function \(\varphi (t) = t - \sqrt{t}\). There are several generalization of Darvay’s method with the square root function [34, 46, 47].

In this paper we use the function \(\varphi (t) = t - \sqrt{t}\) to obtain a new search direction for short step IPA that solves \(\mathcal {P}_{*}(\kappa )\)-LCPs. We mention that this direction has been applied for the first time for LCPs in [17] and for horizontal LCPs in [19].

2.1 Scaling the Newton system

Scaling is a usual technique in the area of IPAs (see [40]). The aim of it is to simplify the analysis of IPAs. We use the following notations

$$\begin{aligned} \mathbf {v}:= \sqrt{\frac{\mathbf {x}\mathbf {s}}{\mu }},\quad \mathbf {d_x}:= \frac{\mathbf {v}\varDelta \mathbf {x}}{\mathbf {x}},\quad \mathbf {d_s}:= \frac{\mathbf {v}\varDelta \mathbf {s}}{ \mathbf {s}},\quad \mathbf {d}:= \sqrt{\frac{\mathbf {x}}{\mathbf {s}}}. \end{aligned}$$
(8)

Using (8) a simple computation yields the following scaled Newton system

$$\begin{aligned} \begin{aligned} -DMD \mathbf {d_x}+ \mathbf {d_s}&= \mathbf {0}, \\ \mathbf {d_x}+ \mathbf {d_s}&= \mathbf {p_v}, \end{aligned} \end{aligned}$$
(9)

where

$$\begin{aligned} \mathbf {p_v}= \frac{\varphi (\mathbf {e})-\varphi (\mathbf {v}^2)}{\mathbf {v}\varphi '(\mathbf {v}^2)}\quad \text { and }\quad D = \mathrm{diag}(\mathbf {d}). \end{aligned}$$

For more details on obtaining the vector \(\mathbf {p_v}\) we refer to [13].

2.2 Centrality measure and some bounds

We use the following centrality measure

$$\begin{aligned} \delta (\mathbf {x},\mathbf {s};\mu ) = \frac{||\mathbf {p_v}||}{2}. \end{aligned}$$

This measure meets our expectations: it is nonnegative and it is zero if and only if the \((\mathbf {x},\mathbf {s})\) vector pair lies on the central path \(\mathcal {C}\). In this paper we use the function given by (5) which was investigated for the first time in [16]. This function is not invertible on \(\mathbb {R}_+\), but it is invertible on the interval \(\displaystyle \left[ \frac{1}{4},\infty \right) \). Thus, we can apply \(\varphi \) for \(\mathbf {v}^2 = \frac{\mathbf {x}\mathbf {s}}{\mu }\) if all coordinates of this vector are greater than or equal to \(\displaystyle \frac{1}{4}\). For the function (5) we get

$$\begin{aligned} \mathbf {p_v}= \frac{2(\mathbf {v}-\mathbf {v}^2)}{2\mathbf {v}- \mathbf {e}} \end{aligned}$$

and

$$\begin{aligned} \delta (\mathbf {x},\mathbf {s};\mu ) = \frac{||\mathbf {p_v}||}{2} = \left| \left| \frac{\mathbf {v}-\mathbf {v}^2}{2\mathbf {v}- \mathbf {e}}\right| \right| . \end{aligned}$$
(10)

Using Lemma 3 for the system (7) we can easily prove the following inequalities.

Corollary 1

Let \(M\in \mathcal {P}_*(\kappa )\), \((\mathbf {x},\mathbf {s})\in \mathcal {F^+}\) and \(\mu >0\) be given. Moreover, let \(\delta := \delta (\mathbf {x}, \mathbf {s}; \mu )\).

  1. (i)

    System (7) has a unique solution \((\varDelta \mathbf {x}, \varDelta \mathbf {s})\) that satisfies the following bounds

    $$\begin{aligned}&||\varDelta \mathbf {x}\varDelta \mathbf {s}||_{\infty }\le \left( \frac{1}{4} + \kappa \right) 4\,\mu \,\delta ^2, \quad ||\varDelta \mathbf {x}\varDelta \mathbf {s}||_1\le \left( \frac{1}{2} + \kappa \right) 4\,\mu \,\delta ^2,\\&||\varDelta \mathbf {x}\varDelta \mathbf {s}||_2\le \sqrt{\left( \frac{1}{4} + \kappa \right) \left( \frac{1}{2} + \kappa \right) }4\,\mu \,\delta ^2. \end{aligned}$$
  2. (ii)

    System (9) has a unique \((\mathbf {d_x}, \mathbf {d_s})\) solution for which

    $$\begin{aligned}&||\mathbf {d_x}\mathbf {d_s}||_{\infty }\le \left( 1 + 4\,\kappa \right) \delta ^2, \quad ||\mathbf {d_x}\mathbf {d_s}||_1\le \left( 2 + 4\,\kappa \right) \delta ^2,\\&||\mathbf {d_x}\mathbf {d_s}||_2\le \sqrt{\left( 1 + 4\,\kappa \right) \left( 2 + 4\,\kappa \right) }\delta ^2 \end{aligned}$$

    hold.

3 The IPA and its analysis

In this section, we present a full Newton-step IPA for the \(\mathcal {P}_*(\kappa )\)-LCP. It is based on a new search direction coming from the system (7) using the function \(\varphi (t) = t - \sqrt{t}\).

figure a

3.1 Analysis of the short step IPA

The next theorem states that the Algorithm 1 is well-defined. This means that the full-Newton step is feasible in each iteration and the initial conditions for the new \((\mathbf {x}, \mathbf {s})\in \mathcal {F}^+\) and \(\mu \) hold, i. e. \(\delta := \delta (\mathbf {x},\mathbf {s};\mu )<\tau \) and \(\mathbf {v}^2>\frac{1}{4} \mathbf{e}\).

Theorem 1

Let \((\mathbf {x},\mathbf {s})\in \mathcal {F}^+\) be given such that \(\delta (\mathbf {x}, \mathbf {s}; \mu ) < \tau \) and \(\mathbf {v}^2 > \frac{1}{4} \, \mathbf {e}\) where \(\tau =\frac{1}{4\kappa + 2}\) and \(\mathbf {v}^2 = \frac{\mathbf {x}\mathbf {s}}{\mu }\). The Newton directions \((\varDelta \mathbf {x},\varDelta \mathbf {s})\) are the unique solutions of the system (7). Then \((\mathbf {x}^+,\mathbf {s}^+)\in \mathcal {F}^+\) such that \(\left( \mathbf {v}^+\right) ^2>\frac{1}{4}\,\mathbf {e}\) and \(\delta (\mathbf {x}^+,\mathbf {s}^+,\mu ^+) < \tau \), where \(\mathbf {x}^+ := \mathbf {x}+ \varDelta \mathbf {x}, \mathbf {s}^+:=\mathbf {s}+ \varDelta \mathbf {s}\), \(\mu ^+ = (1-\theta )\mu \) and \(\mathbf {v}^+ := \sqrt{\frac{\mathbf {x}^+ \mathbf {s}^+}{\mu ^+}}\).

Proof

Let us prove that \(\mathbf {x}+ \alpha \varDelta \mathbf {x}\) and \(\mathbf {s}+ \alpha \varDelta \mathbf {s}\) are strictly feasible solutions for any \(\alpha \in [0,1]\). Then,

$$\begin{aligned} \mathbf {x}(\alpha )\mathbf {s}(\alpha )&:= (\mathbf {x}+ \alpha \varDelta \mathbf {x})(\mathbf {s}+ \alpha \varDelta \mathbf {s}) = \mathbf {x}\mathbf {s}+ \alpha (\mathbf {s}\varDelta \mathbf {x}+ \mathbf {x}\varDelta \mathbf {s}) + \alpha ^2\varDelta \mathbf {x}\varDelta \mathbf {s}\\&=\mu \mathbf {v}^2 + \mu \mathbf {v}\alpha (\mathbf {d_x}+ \mathbf {d_s}) + \alpha ^2 \mu \mathbf {d_x}\mathbf {d_s}\\&= \mu ((1-\alpha )\mathbf {v}^2 + \alpha (\mathbf {v}^2 + \mathbf {v}\mathbf {p_v}) + \alpha ^2\mathbf {d_x}\mathbf {d_s})\\&= \mu ((1-\alpha )\mathbf {v}^2 + \alpha (\mathbf {v}^2 + \mathbf {v}\mathbf {p_v}+ \alpha \mathbf {d_x}\mathbf {d_s})). \end{aligned}$$

Now \(\mathbf {x}(\alpha )\mathbf {s}(\alpha )>0\) if \(\mathbf {v}^2 + \mathbf {v}\mathbf {p_v}+ \alpha \mathbf {d_x}\mathbf {d_s}>0\). Bounding from below and using \(\mathbf {v}>\frac{1}{2}\mathbf {e}\) and Corollary 1 we get

$$\begin{aligned} \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}+ \alpha \mathbf {d_x}\mathbf {d_s}&\ge \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}- \alpha ||\mathbf {d_x}\mathbf {d_s}||_{\infty } \mathbf {e}\ge \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}- \alpha (1 + 4 \kappa ) \delta ^2 \mathbf {e}\\&\ge \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}- \alpha \frac{1+4\kappa }{(2+4\kappa )^2} \mathbf {e}\\&> \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}- \alpha \mathbf {e}\ge \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}- \mathbf {e}\\&= \mathbf {v}^2 + \mathbf {v}\frac{2(\mathbf {v}- \mathbf {v}^2)}{2\mathbf {v}- \mathbf {e}} - \mathbf {e}= \frac{(\mathbf {v}- \mathbf {e})^2}{2\mathbf {v}- \mathbf {e}} \ge 0. \end{aligned}$$

The functions \(\alpha \rightarrow \mathbf {x}(\alpha )\) and \(\alpha \rightarrow \mathbf {s}(\alpha )\) are continuous, \(\mathbf {x}(0)> 0,\; \mathbf {s}(0) > 0\) and \(\mathbf {x}(\alpha )\, \mathbf {s}(\alpha ) > 0\) for each \(\alpha \in [0,1]\), thus we proved that \((\mathbf {x}^+,\mathbf {s}^+)\in \mathcal {F}^+\).

Now we need to show that \((\mathbf {v}^+)^2>\frac{1}{4}\mathbf {e}\). Actually we prove a stronger statement, namely that

$$\begin{aligned} v^+_{\min } {:= \min _i \; v^+_i} \ge \frac{\sqrt{3}}{2\eta }, \end{aligned}$$
(11)

where \(\eta = \sqrt{1-\theta }\) and \(\displaystyle \theta =\frac{1}{(8 + 9\kappa )\sqrt{n}}\).

We use a weaker condition for \(\delta \)

$$\begin{aligned} \delta ^2\le \tau ^2\le \frac{1}{(4\kappa + 2)^2} = \frac{1}{16\kappa ^2 + 16\kappa + 4}\le \frac{1}{16\kappa + 4}=\frac{1}{4 (1 + 4\kappa )}. \end{aligned}$$
(12)

We have

$$\begin{aligned} (1-\theta )(\mathbf {v}^+)^2=\frac{\mathbf {x}^+\mathbf {s}^+}{\mu } = \mathbf {v}^2 + \mathbf {v}\mathbf {p_v}+ \mathbf {d_x}\mathbf {d_s}= \frac{\mathbf {v}^2}{2\mathbf {v}-\mathbf {e}} + \mathbf {d_x}\mathbf {d_s}\ge \mathbf {e}+ \mathbf {d_x}\mathbf {d_s}, \end{aligned}$$
(13)

since \(\displaystyle \frac{\mathbf {v}^2}{2\mathbf {v}-\mathbf {e}}\ge \mathbf {e}.\) Using inequalities (12) and (13) we get

$$\begin{aligned} (\mathbf {v}^+_{\min })^2\ge & {} \frac{1}{1-\theta }(1 - ||\mathbf {d_x}\mathbf {d_s}||_{\infty })\ge \frac{1}{1-\theta }(1 - (1 + 4\kappa )\delta ^2)\\\ge & {} \frac{1}{1-\theta }\left( 1 - \frac{1}{4}\right) =\frac{3}{4(1-\theta )}. \end{aligned}$$

Taking square root we obtain (11), thus the inequality \(\displaystyle \mathbf {v}^+>\frac{1}{2}\,\mathbf {e}\) follows. Now we need to show that \(\delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+)<\tau \). We have

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+)&= \left| \left| \frac{\mathbf {v}^+ - (\mathbf {v}^+)^2}{2\mathbf {v}^+ - \mathbf {e}}\right| \right| = \left| \left| \frac{\mathbf {v}^+}{(2\mathbf {v}^+ - \mathbf {e})(\mathbf {e}+ \mathbf {v}^+)} \, (\mathbf {e}- (\mathbf {v}^+)^2)\right| \right| \\&= \left| \left| \frac{\mathbf {v}^+}{(2\mathbf {v}^+ - \mathbf {e})(\mathbf {e}+ \mathbf {v}^+)}\right| \right| \; \left| \left| (\mathbf {e}- (\mathbf {v}^+)^2)\right| \right| . \end{aligned}$$

Based on the first term of the previous expression, we introduce the function

$$\begin{aligned} f(t):= \frac{t}{(2t - 1)(t + 1)}. \end{aligned}$$

Using f we can write

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+) = ||f(\mathbf {v}^+)||\cdot ||(\mathbf {e}- (\mathbf {v}^+)^2)|| \le ||f(\mathbf {v}^+)||_{\infty }||\mathbf {e}- (\mathbf {v}^+)^2|| \end{aligned}$$
(14)

The function f is continuous, monotonically decreasing and positive on the interval \(\left( \frac{1}{2},\infty \right) \). Combining the properties of the function f with the inequality (11) we obtain

$$\begin{aligned} 0 < [f(\mathbf {v}^+)]_i =[f(v^+_i)]\le f(v^+_{\min })\le f\left( \frac{\sqrt{3}}{2\eta }\right) ,\quad \text {for all indices }\, i\in \mathcal {I}. \end{aligned}$$
(15)

Using the inequalities (14) and (15) we get

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+)\le f\left( \frac{\sqrt{3}}{2 \eta }\right) \, ||\mathbf {e}- (\mathbf {v}^+)^2|| = \frac{\sqrt{3} \, \eta }{(\sqrt{3} - \eta )(\sqrt{3} + 2\eta )} \; ||\mathbf {e}- (\mathbf {v}^+)^2||. \end{aligned}$$
(16)

Bounding the last term yields

$$\begin{aligned} ||\mathbf {e}- (\mathbf {v}^+)^2||= & {} \left| \left| \mathbf {e}-\frac{1}{\eta ^2}\frac{\mathbf {x}^+\mathbf {s}^+}{\mu }\right| \right| = \left| \left| -\frac{1}{\eta ^2}\left( \frac{\mathbf {x}^+\mathbf {s}^+}{\mu }-\eta ^2 \mathbf {e}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( \mathbf {v}^2+\mathbf {v}\mathbf {p_v}+\mathbf {d_x}\mathbf {d_s}-\eta ^2\mathbf {e}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( \mathbf {v}^2+\mathbf {v}\frac{2(\mathbf {v}- \mathbf {v}^2)}{2\mathbf {v}- \mathbf {e}}+\mathbf {d_x}\mathbf {d_s}-\eta ^2\mathbf {e}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( \frac{\mathbf {v}^2}{2\mathbf {v}-\mathbf {e}}+\mathbf {d_x}\mathbf {d_s}-\eta ^2\mathbf {e}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( \mathbf {e}+ \frac{(\mathbf {e}-\mathbf {v})^2}{2\mathbf {v}-\mathbf {e}}+\mathbf {d_x}\mathbf {d_s}-\eta ^2\mathbf {e}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( (1-\eta ^2)\mathbf {e}+ \frac{(\mathbf {v}-\mathbf {v}^2)^2}{(2\mathbf {v}-\mathbf {e})^2}\cdot \frac{2\mathbf {v}-\mathbf {e}}{\mathbf {v}^2}+\mathbf {d_x}\mathbf {d_s}\right) \right| \right| \\= & {} \left| \left| -\frac{1}{\eta ^2}\left( \theta \mathbf {e}+ \frac{\mathbf {p_v}^2}{4}\cdot \frac{2\mathbf {v}-\mathbf {e}}{\mathbf {v}^2}+\mathbf {d_x}\mathbf {d_s}\right) \right| \right| \\\le & {} \frac{1}{\eta ^2}\Bigg (||\theta \mathbf {e}|| + \Bigg |\Bigg |\frac{\mathbf {p_v}^2}{4}\cdot \frac{2\mathbf {v}- \mathbf {e}}{\mathbf {v}^2}\Bigg |\Bigg | + ||\mathbf {d_x}\mathbf {d_s}||\Bigg )\\\le & {} \frac{1}{\eta ^2}(\theta \sqrt{n} + \delta ^2 + \delta ^2\sqrt{2(1+4\kappa )(1+2\kappa )}). \end{aligned}$$

In the last inequality we used Corollary 1, the inequalities \(\mathbf {v}^2 \ge 2\mathbf {v}- \mathbf {e}\ge \mathbf {0}\) and the definition of \(\delta \) given in (10). Using (16) and the previous estimation we are ready to bound \(\delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+)\) as follows

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+) \le g(\eta ) \, (\theta \, \sqrt{n} + \left( 1 + c(\kappa )\right) \, \delta ^2), \end{aligned}$$
(17)

where

$$\begin{aligned} g(\eta )= \frac{\sqrt{3}}{-2\eta ^3 + \sqrt{3} \, \eta ^2 + 3 \, \eta } \qquad \hbox {and} \qquad c(\kappa ) = \sqrt{\left( 1 + 4 \,\kappa \right) \left( 2 + 4 \, \kappa \right) }. \end{aligned}$$

Let us bound from above \(\theta \) as

$$\begin{aligned} \displaystyle \theta =\frac{1}{(9\kappa + 8)\sqrt{n}}\le \frac{1}{8} < \frac{1}{4}. \end{aligned}$$

Therefore,

$$\begin{aligned} \displaystyle \eta = \sqrt{1-\theta } \ge \sqrt{\frac{3}{4}} = \frac{\sqrt{3}}{2}.\ \end{aligned}$$

Thus, using that the function g is decreasing we get

$$\begin{aligned} \displaystyle g(\eta ) \le g\left( \frac{\sqrt{3}}{2}\right) = \frac{2}{3}. \end{aligned}$$

Using the inequality (17) we obtain the following bound

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+) \le \frac{2}{3} \, (\theta \, \sqrt{n} + \left( 1 + c(\kappa )\right) \, \tau ^2). \end{aligned}$$
(18)

Substituting the values of the parameters \(\theta \) and \(\tau \) and using the fact that \(c(\kappa )\le \frac{3}{2} + 4\kappa \) we have

$$\begin{aligned} \delta (\mathbf {x}^+,\mathbf {s}^+;\mu ^+)&\le \frac{2}{3}\left( \frac{1}{8 + 9\kappa } + \frac{1 + c(\kappa )}{(2 + 4\kappa )^2}\right) \le \frac{2}{3}\left( \frac{1}{8 + 9\kappa } + \frac{\frac{5}{2} + 4\kappa }{(2 + 4\kappa )^2}\right) \\&=\frac{1}{3}\left( \frac{2}{8 + 9\kappa } + \frac{5 + 8\kappa }{(2 + 4\kappa )^2}\right) =\frac{104\kappa ^2 + 141\kappa + 48}{3(8 + 9\kappa )(2 + 4\kappa )^2}\\&\le \frac{108\kappa ^2 + 150\kappa + 48}{3(8 + 9\kappa )(2 + 4\kappa )^2} = \frac{3(8 + 9\kappa )(2 + 4\kappa )}{3(8 + 9\kappa )(2 + 4\kappa )^2}=\frac{1}{2 + 4\kappa }=\tau . \end{aligned}$$

So far we proved that the Algorithm 1 is well defined. The iterations are repeatable because of the inheritance of initial conditions. \(\square \)

3.2 Complexity analysis of the short step IPA

In this subsection we prove a theorem which deals with the complexity of the Algorithm 1. We show that this algorithm enjoys the best theoretical polynomial complexity.

First we examine the duality gap. If the Newton-step took us exactly to the \(\mu \)-center then \(\displaystyle \frac{(\mathbf {x}^+)^T\mathbf {s}^+}{n} = \mu \,\) would hold. The next lemma examines the relationship between the terms \(\,\displaystyle \frac{(\mathbf {x}^+)^T\mathbf {s}^+}{n}\,\) and \(\,\mu \).

Lemma 1

Let \((\mathbf {x},\mathbf {s})\in \mathcal {F}^+\) be given such that \(\delta (\mathbf {x}, \mathbf {s}; \mu ) < \tau \) and \(\mathbf {v}^2 = \frac{\mathbf {x}\mathbf {s}}{\mu }> \frac{1}{4} \, \mathbf {e}\), where \(\tau =\frac{1}{4\kappa + 2}\). The solution of the Newton system (7) is denoted by \((\varDelta \mathbf {x},\varDelta \mathbf {s})\). Then,

$$\begin{aligned} \frac{(\mathbf {x}^+)^T\mathbf {s}^+}{n}\le \mu \left( 1 + \frac{2\delta ^2}{n}\right) , \end{aligned}$$

where \(\mathbf {x}^+ = \mathbf {x}+ \varDelta \mathbf {x}\), \(\mathbf {s}^+=\mathbf {s}+ \varDelta \mathbf {s}\) and \(\delta =\delta (\mathbf {x}, \mathbf {s}; \mu )\).

Proof

Let us compute \(\mathbf {x}^+\mathbf {s}^+\). We have

$$\begin{aligned} \mathbf {x}^+\mathbf {s}^+&=\mu (\mathbf {v}^2 + \mathbf {v}\mathbf {p_v}+ \mathbf {d_x}\mathbf {d_s})=\mu \Bigg (\mathbf {e}+ \frac{\mathbf {p_v}^2}{4}\cdot \frac{2\mathbf {v}-\mathbf {e}}{\mathbf {v}^2} + \mathbf {d_x}\mathbf {d_s}\Bigg ) \nonumber \\&\le \mu \Bigg (\mathbf {e}+ \frac{\mathbf {p_v}^2}{4} + \mathbf {d_x}\mathbf {d_s}\Bigg ). \end{aligned}$$
(19)

The last inequality follows from \(\mathbf {v}\ge \frac{1}{2}\,\mathbf {e}\), which implies \(\displaystyle \frac{2\mathbf {v}-\mathbf {e}}{\mathbf {v}^2}\in [0,1]^n\). From (19) we get

$$\begin{aligned} (\mathbf {x}^+)^T\mathbf {s}^+&=\mathbf {e}^T(\mathbf {x}^+\mathbf {s}^+)\le \mu \left( \mathbf {d_x}^T\mathbf {d_s}+ \frac{||\mathbf {p_v}||^2}{4} + n\right) =\mu \left( \mathbf {d_x}^T\mathbf {d_s}+ \delta ^2 + n\right) . \end{aligned}$$
(20)

We use the following simple estimation:

$$\begin{aligned} \mathbf {d_x}^T\mathbf {d_s}&= \frac{||\mathbf {d_x}+ \mathbf {d_s}||^2 - ||\mathbf {d_x}- \mathbf {d_s}||^2}{4} \le \frac{||\mathbf {d_x}+ \mathbf {d_s}||^2}{4} = \frac{||\mathbf {p_v}||^2}{4} = \delta ^2. \end{aligned}$$
(21)

From (20) and (21) we obtain the statement of the lemma. \(\square \)

Based on the previous lemma it is easy to derive the complexity result of the algorithm.

Theorem 2

Let \(\theta = \frac{1}{(9\kappa + 8)\sqrt{n}}\) and \(\tau = \frac{1}{4\kappa + 2}\). If \(M\in \mathcal {P}_*(\kappa )\), \((\mathbf {x}^0, \mathbf {s}^0)\in \mathcal {F}^+\), \(\mu ^0>0\), \(\;\delta (\mathbf {x}^0,\mathbf {s}^0;\mu ^0)<\tau \) and \(\frac{\mathbf {x}^0\mathbf {s}^0}{\mu ^0}>\frac{1}{4} \mathbf {e}\) then the Algorithm 1 provides an \(\varepsilon \)-solution of the LCP after at most

$$\begin{aligned} \left\lceil 9(\kappa + 1)\,\sqrt{n}\,\log \left( \frac{2 n \mu ^0}{\varepsilon }\right) \right\rceil \end{aligned}$$

iterations.

Proof

Let us assume that we are in the \((k+1)^{\text {th}}\) iteration of our algorithm. At the beginning of the iteration we have \(\mathbf {x}=\mathbf {x}_{k},\,\mathbf {s}=\mathbf {s}_{k}\), \(\mu = \mu _{k}\) and \(\delta _k = \delta (\mathbf {x}_k,\mathbf {s}_k;\mu _k)\). Using Lemma 1 we have

$$\begin{aligned} (\mathbf {x}_{k+1})^T\mathbf {s}_{k+1}\le & {} \mu _{k}\,(2\delta _k^2 + n) \le \mu _{k}\,(2\tau ^2 + n)\nonumber \\= & {} \mu _{k}\,\left( \frac{2}{(4\kappa + 2)^2} + n\right) \le \mu _{k}\,\left( \frac{1}{2} + n\right) . \end{aligned}$$
(22)

Using the equality \(\mu _{k}=(1-\theta )^{k}\mu _0\) and the stopping criteria of our IPA we conclude that the algorithm will surely stop if the inequality

$$\begin{aligned} (1-\theta )^{k}\,\mu _0\,\left( \frac{1}{2} + n\right) \le \varepsilon \end{aligned}$$
(23)

holds. Taking logarithm of (23) we get

$$\begin{aligned} k\,\log (1-\theta ) + \log \left( \mu _0\,\left( \frac{1}{2} + n\right) \right) \le \log \varepsilon . \end{aligned}$$
(24)

Using the inequality \(\log (1-\theta )\le -\theta \) and the fact that \(\theta =\frac{1}{(9\kappa + 8)\sqrt{n}}\) we get that the algorithm stops if

$$\begin{aligned} k\ge \frac{1}{\theta }\,\log \left( \frac{\mu ^0}{\varepsilon }\,\left( \frac{1}{2} + n\right) \right) = (9\kappa + 8)\,\sqrt{n}\,\log \left( \frac{\mu ^0}{\varepsilon }\,\left( \frac{1}{2} + n\right) \right) . \end{aligned}$$
(25)

Now using simple computation we have

$$\begin{aligned} (9\kappa + 8)\,\sqrt{n}\,\log \left( \frac{\mu ^0}{\varepsilon }\,\left( \frac{1}{2} + n\right) \right)\le & {} 9(\kappa + 1)\,\sqrt{n}\,\log \left( \frac{\mu ^0}{\varepsilon }\,\left( \frac{1}{2} + n\right) \right) \\\le & {} 9(\kappa + 1)\,\sqrt{n}\,\log \left( \frac{2n\mu ^0}{\varepsilon }\right) , \end{aligned}$$

which proves the theorem. \(\square \)

Corollary 2

Let \(\theta = \frac{1}{(9\kappa + 8)\sqrt{n}}\) and \(\tau = \frac{1}{4\kappa + 2}\). If \(M\in \mathcal {P}_*(\kappa )\), \((\mathbf {x}^0, \mathbf {s}^0)\in \mathcal {F}^+\), \(\mu ^0>0\), \(\;\delta (\mathbf {x}^0,\mathbf {s}^0;\mu ^0)<\tau \) and \(\frac{\mathbf {x}^0\mathbf {s}^0}{\mu ^0}>\frac{1}{4} \mathbf {e}\) then the Algorithm 1 provides an \(\varepsilon \)-solution of the LCP after at most

$$\begin{aligned} O\left( (1 + \kappa )\,\sqrt{n}\,\log \left( \frac{n\,\mu ^0}{\varepsilon }\right) \right) \end{aligned}$$

iterations.

Proof

Using \(n \ge 2\) and \(\mu ^0 \ge \varepsilon \) we obtain \(\frac{2n\mu ^0}{\varepsilon } \le \frac{n^2(\mu ^0)^2}{\varepsilon }\). Thereofore, we have

$$\begin{aligned} 9(\kappa + 1)\,\sqrt{n}\,\log \left( \frac{2n\mu ^0}{\varepsilon }\right) \le 18(\kappa + 1)\,\sqrt{n}\,\log \left( \frac{n\,\mu ^0}{\varepsilon }\right) , \end{aligned}$$

which implies the corollary. \(\square \)

Using the \(\varepsilon \)-solution provided by Corollary 2 we can apply the rounding procedure presented in Illés et al. [30] to obtain a maximally complementary solution in strongly polynomial time.

4 Numerical results

We implemented our algorithm in the C++ programming language using the object-oriented code developed by Darvay and Takó [20]. The main difference between our implementation and the analysed theoretical version of the algorithm is that in the implemented version of the PD IPA, the barrier update parameter, \(\theta \) is a given constant between 0 and 1, while in the theoretical version of the algorithm \(\theta = \frac{1}{(9 \, \kappa + 8) \, \sqrt{n}}\). The reason of this modification lies in the fact that usually the handicap is not exactly known or might be very (exponentially) large. In the latter case the barrier update parameter \(\theta \) is much smaller than it is necessary for a reliable implementation.

However, the price of this deviation from the theoretical version of the algorithm is that there is no guarantee that the full step can be performed. In each iteration we calculated the maximum step length \(\alpha _{\max }\) using the classical ratio test, which guarantees that the new iterates \(\mathbf {x}^+ = \mathbf {x}+ \alpha \, \varDelta \mathbf {x}\), \(\mathbf {s}^+=\mathbf {s}+ \alpha \, \varDelta \mathbf {s}\) remain nonnegative. To ensure the strict feasibility of the new iterates we used a factor \(\rho =0.95\) to shorten the step length, thus the used step length is \(\alpha ^+ = \rho \, \alpha _{\max }\).

Another problem caused by choosing arbitrary \(\theta \in (0, 1)\) is related to the selection of the new value of the \(\mu \). Since the \(\delta (\mathbf {x}^+, \mathbf {s}^+; \mu ) < \tau \) may not hold for any \(0< \tau < 1\) the selection of the target \(\mu \) value need to be modified comparing the theoretical version of the PD IPA. Namely, the new target value of \(\mu \) is computed as \(\mu = (1 - \theta ) \, \frac{(\mathbf {x}^+)^T \mathbf {s}^+}{n}\), where \(\theta \in \{ 0.999, 0.1 \}\) in our computations presented in this paper.

We compared three primal–dual algorithms based on the AET technique. The only difference between these methods were in the selection of the function \(\varphi \). We considered the identity map, the square root function and the our one, given in (5).

We tested the above mentioned algorithms on two different types of test problems. The first set of problems used sufficient matrices with positive handicap provided by Illés and Morapitiye [29]. For these problems we set the initial interior point as \(\mathbf {x}^0 = \mathbf {s}^0 = \mathbf {e}\), because the righthand side vector has been computed as \(\mathbf {q}:= - M \, \mathbf {e}+ \mathbf {e}\).

The second test set problem is based on the P-matrix proposed by Zsolt Csizmadia and presented in [23]: \(A_n=\left( a_{ij}\right) _{1\le i, j \le n}\), where

$$\begin{aligned} a_{ij} = \left\{ \begin{array}{rl} 0~~~~~ &{} \text {if i<j,}\\ 1~~~~~ &{} \text {if i=j,}\\ -1~~~~~ &{} \text {if i>j.} \end{array} \right. \end{aligned}$$

The handicap \(\kappa (A_n) \ge 2^{2n-8}-\frac{1}{4}\) as de Klerk and E.-Nagy [23] proved. The righthand side vectors of sufficient LCPs using matrix \(A_n\) has been computed as \(\mathbf {q}:= - A_n \, \mathbf {e}+ \mathbf {e}\), thus \(\mathbf {x}^0 = \mathbf {s}^0 = \mathbf {e}\) are ideal initial interior points for starting analized IPAs.

Based on the definition of \(P_*(\kappa )\)-matrices, (2), for a given vector \({\mathbf {y}}\in \mathbb {R}^n\) and \(P_*(\kappa )\)-matrix, M with the property that \({\mathbf {y}}^T M \, {\mathbf {y}}< 0\), corresponds to a \(\kappa ({\mathbf {y}}) > 0\) parameter such that (2) holds with equality for the given vector \({\mathbf {y}}\). Interestingly enough, \(\kappa (\varDelta \mathbf {x})\) may take large values when the handicap of the given \(P_*(\kappa )\)-matrix M, is large, like for the Csizmadia-matrix, \(A_n\).

Before presenting and analysing our computational results let us shortly discuss from complexity point of view the three implemented algorithms that differ only in the AET function.

Observe that independently from the used AET functions for transforming the central path equation (for \(\varphi (t) = t\) see [35]; \(\varphi (t)=\sqrt{t}\) see [5], and for \(\varphi (t)=t - \sqrt{t}\) see Corollary 2 of this paper), the best complexity results for primal-dual IPAs for sufficient LCPs are just the same, namely \(O((1+\kappa )\sqrt{n}\log \frac{n}{\varepsilon })\). Let us discuss computational performance of the PD IPAs differing only in the used \(\varphi \) functions on those earlier described two sufficient LCP test sets.

4.1 Computational results on LCPs with positive handicap

The sufficient matrices generated by Illés and Morapitiye [29] have positive handicap, but not very large. Furthermore, these matrices were generated using same rules and similar starting matrices therefore does not differ significantly from each other. This observation has been also confirmed by the behavior of the different versions of the PD IPAs on the first test set of problems.

We tested PD IPAs on sufficient LCPs for sizes \(n=10, 20, 50, 100, 200, 500\). For each of these sizes the test set contains 10 problems. Due to the fact that for a given size the iteration number and the CPU time was very similar, instead of giving results for separate problems, we present a table containing the necessary average iteration number for different sizes (Table 1).

Table 1 The average number of iterations for LCPs based on sufficient matrices presented in [29], which have positive handicap
Full size table
Table 2 The number of iterations for LCPs with Csizmadia-matrices
Full size table

4.2 Computational results on LCPs with Csizmadia-matrix

The handicap of Csizmadia matrix \(\kappa (A_n) \ge 2^{2n-8}-\frac{1}{4}\) (see [23]). The worst case complexity result with the default value of \(\theta \) (see Algorithm 1) ensures that the IPA has polynomial iteration complexity depending on \(\kappa \). Since we use much larger values for the update parameter \(\theta \) in our computations, no theoretical guarantee follows for polynomial number of iteration of the IPA from Theorem 2.

On the other hand, the implemented variants of theoretical IPAs for different AETs with several safeguard modifications, as we expected, performed very well (see Tables 2 and 3). From the obtained iteration numbers no one could see any effect of the exponentially large handicap on the performance of the implemented IPAs. However, if we analyse each iteration, we could observe an interesting behaviour of the implemented algorithm. Namely, there is a strict relation between the \(\kappa (\varDelta \mathbf {x})\) and the corresponding step length, \(\alpha \).

For instance, at \(A_{100}, \, \theta = 0.999\), the \(\kappa (\varDelta \mathbf {x}_k)\) was continuously decreasing from a very large number \(1.88905 \cdot 10^{34}\) to 0, and in the meantime the step length \(\alpha \) was continuously increasing from a very small number \(7.49944 \cdot 10^{-18}\) to 1, that is a full Newton-step (see [21]). Naturally, the gap decreases more when the step length is larger.

We tested the effect of the barrier parameter update \(\theta \) by using different values of it. Results with \(\theta = 0.1\) (see Table 3) show that the selection of AET function could have minor effect on the iteration number of the PD IPA depending on the value of \(\theta \).

However, in case of \(\theta = 0.1\) the influence of the handicap of the given matrix on the increase of the iteration numbers can be observed, as well.

Table 3 The number of iterations for LCPs with Csizmadia-matrices
Full size table

5 Conclusions

In this article we have developed a new short-step path-following algorithm for solving LCPs. Our approach is a generalization of the results presented in [15]. We dealt with \(\mathcal {P}_*(\kappa )-LCP\)s and derived an algorithm with the currently best known complexity. We provided some numerical results which prove the practical efficiency of the method as well.

Darvay and his coauthors together with their followers demonstrated that this approach can be used widely. Moreover, some other interesting questions can be raised, which need further investigation. Can any Darvay-type AET be used to derive search directions for IPAs for different optimization problem classes? Is it possible to characterize those one-to-one transformations that are applicable to obtain new search directions for IPAs of different problem classes? Can we rank different (Darvay-type) search directions as one is better than the other in the sense of generating significantly less iterations in a practical implementation? We think that there are possibilities to modify the complexity analysis of Darvay’s IPAs for LO [13, 18], with careful applications of ideas presented in the book of Roos et al. [40]. As a further research we would like to extend the AET approach to the case of general LCPs analysed by Illés, Nagy and Terlaky [32]. Another question of interest is the following: is there a short-step IPA for \(\mathcal {P}_*(\kappa )\)-LCPs (or for a well defined subclass of it), where the centrality parameter \(\tau \) doesn’t depend on \(\kappa \)? We suspect that there isn’t. Understanding this phenomenon requires further research as well.