1.
Introduction
In this study, we focus on the nonnegative cone horizontal linear complementarity problem (HLCP): Given the matrices A,B∈Rn×n and a source term q∈Rn, our goal is to find a pair of vectors, denoted as x,y, that meet the following criteria:
where Rn+={x∈Rn∣x⩾0} is the so-called nonnegative cone. Here, ⟨⋅,⋅⟩ is the Euclidean inner product and ‖⋅‖ denotes the Euclidean norm. Obviously, if either A=In or B=In, the HLCP simplifies to the linear complementarity problem (LCP).
The HLCP has many applications in noncooperative games, optimization problems, economic equilibrium problems, and traffic assignment problems [1]. The finite-difference discretization of hydrodynamic lubrication also gives rise to the HLCP [2]. Many scholars have studied this problem. Various methods can be employed to address this issue, including interior point approaches [3], neural network techniques [4], and verification techniques [5]. In [6], Mezzadri and Galligani introduced a class of projected splitting methods.
Recently, the modulus-based matrix-splitting approach, a widely used technique, has been applied to many kinds of linear/nonlinear complementarity problems. For example, Mezzadri and Galligani in 2019 presented the modulus-based matrix-splitting method to solve the HLCP [7]. Furthermore, the convergence of the methods proposed in [7] has been enlarged by Zheng and Vong [8]. Zhang et al. [9] have introduced a two-step parallel iterative approach for solving large sparse HLCP problems.
The first-order methods have been the subject of considerable attention in recent years, largely due to their minimal storage requirements. One of the most widely known first-order methods is the spectral gradient method, which was first introduced by Barzilai and Borwein in reference [10] and subsequently extended by Raydan and La Cruz in references [11,12,13]. In their study, Zhang and Zhou [14] employed the spectral gradient projection method for unconstrained monotone equations, as proposed by Solodov and Svaiter [15]. This approach is based on the inexact Newton method. Furthermore, Yu et al. demonstrated the efficacy of this approach in solving constrained monotone equations, as evidenced by the favorable outcomes reported in reference [16]. In a recent study, Ibrahim et al. proposed the projection method with inertial step [17] for nonlinear equations, Li et al. studied the modified spectral gradient projection-based algorithm [18] for large-scale constrained nonlinear equations, and Ibrahim et al. proposed the two-step inertial derivative-free projection method [19] for solving nonlinear equations.
Jiang [20] developed a derivative-free descent technique for the nonlinear complementarity problem (NCP) when the nonlinear mapping is directionally differentiable and strongly monotone, utilizing an equivalent system of nonlinear equations derived from the squared Fischer-Burmeister (FB) function. Mangasarian et al. [21] introduced another derivative-free descent method for the strongly monotone NCP by minimizing the implicit Lagrangian function, establishing its global convergence. Ma et al. [22] proposed a smooth Broyden-like method for the NCP, which uses a smooth approximation of the FB function and a derivative-free line search rule, demonstrating global convergence under suitable conditions. This smooth Broyden-like method has also been shown to achieve global and superlinear convergence under appropriate conditions.
Yu et al. [23] converted the absolute value equation into a monotone system and resolved it using a multivariate spectral gradient method. Drawing inspiration from this work, we first reframe the HLCP (1.1) into a fixed-point equation based on the modulus defined within the nonnegative cone. By transforming the implicit fixed-point equation into a monotone system, we introduce a new class of modified multivariate spectral gradient projection methods for its solution. We outline the conditions under which the new equation remains monotone and continuous. It is shown that the proposed iterative method converges to the solution of the HLCP in equation (1.1) under the specified assumptions. Additionally, numerical examples illustrate that the modulus-based modified multivariate spectral gradient projection method is both feasible and efficient.
The proposed algorithm is contrasted with the modulus-based matrix-splitting iterative method, which generally requires matrix inversion, a step not needed here. Since matrix inversion is computationally intensive, the algorithm described in this paper holds a distinct advantage. Compared to conventional gradient descent methods, there is potential for improving the line search technique and the distribution strategy of the spectral gradient in existing algorithms. This paper's algorithm has been refined by incorporating a new line search technique [24] and a more efficient allocation of the spectral gradient, leading to better performance in reducing the number of iterations and CPU time.
The paper is structured as follows: In Section 2, we develop a monotone system of equations that is equivalent to the HLCP. In Section 3, we introduce a modified multivariate spectral gradient projection method grounded in the modulus approach. Section 4 focuses on the convergence analysis of the proposed algorithm. Numerical experiments and their outcomes are detailed and analyzed in Section 5. Finally, concluding observations are provided in Section 6.
2.
The establishment of equivalent equations
This section introduces some notations and auxiliary results.
Theorem 2.1. Let matrices A,B∈Rn×n and vector q∈Rn. Then,
− if (x,y) is a solution to the HLCP in (1.1), then z=12(x−y) fulfills
− if z satisfies (2.1), then
is a solution of the HLCP, where z∈Rn. |z| denotes the absolute value of z.
Proof. With regard to the first statement, since (x,y) is a solution of the HLCP in (1.1), it follows that both variables are nonnegative and can be expressed as Eq (2.2) with z∈Rn. If we put (2.2) into (1.1), we obtain (x,y) as a solution to the complementarity problem if and only if
Rearranging the above equation can easily be written as (2.1).
About the second statement, we commence with (2.1), which can be rewritten as Eq (2.3) through the simple rearrangement of terms. By defining x and y as specified in Eq (2.2), we derive
where x, y is nonnegative. It can be readily confirmed that xi>0 and yi=0 when zi>0, whereas xi=0 and yi>0 when zi<0. If zi=0, then trivially xi=yi=0. Here, zi denotes the i-th component of z, with a similar notation used for xi and yi. Consequently, the complementarity condition is met, and (x,y) constitutes a solution to the HLCP as stated in Eq (1.1). This completes the proof. □
Let x=|z|+z, y=|z|−z, where z∈Rn. According to Theorem 2.1, the HLCP in (1.1) is equivalent to the fixed-point equation
If B−A is invertible, we have
We set
where z∈Rn.
Definition 2.1. A mapping F:Rn→Rn is defined as monotone if, for ∀α,β∈Rn,
is satisfied.
Definition 2.2. A mapping F:Rn→Rn is considered Lipschitz continuous if, for ∀α,β∈Rn, there exists a constant L>0 such that
Assumption 2.1. In Eq (2.4), (B−A) is invertible and (B−A)−1A is semi-positive definite.
Theorem 2.2. F(z) is monotone as long as (B−A)−1A is semi-positive definite.
Proof. For any α,β∈Rn, the dot product αTβ=∑ni=1αiβi. Specifically, we have defined the index set
and
On one side,
and expanding it out, we notice that ∑i∈I2⟨2βi,αi−βi⟩⩽0, and ∑i∈I3⟨2αi,βi−αi⟩⩽0. Hence,
On the flip side,
Obviously, F(z) will be monotone, as long as (B−A)−1A is semi-positive definite. This completes the proof. □
Theorem 2.3. The function F(z) is Lipschitz continuous, meaning there exists a nonnegative constant L such that the following inequality is satisfied:
Proof. For ∀α,β∈Rn, by the Lipschitz condition and (2.4), we have
Since ‖(B−A)−1(B+A)‖⩾0, thus, set L=‖(B−A)−1(B+A)‖+1, and it can be demonstrated that the Lipschitz condition is satisfied. This completes the proof. □
According to Theorems 2.2 and 2.3, under the condition of Assumption 2.1, HLCP (1.1) can be equivalently reformulated to a monotone system:
3.
Algorithm
For monotone system (2.5), the matrix-splitting method was commonly used in the past. In each iteration, the fixed-point equation needs to be inverted to obtain a new iteration point. As we know, it takes a certain amount of time to find the inverse. In order to simplify the time cost, we introduce a new multivariate spectral gradient method to solve (2.5), called the modulus-based modified multivariate spectral gradient projection method (M-MMSGP). The M-MMSGP algorithm employs a novel line search methodology [24] that differs from existing gradient algorithms. It utilizes the multivariate spectral gradient vector as the descent direction and assigns the multivariate spectral gradient as a whole, which significantly reduces the number of iterations of the algorithm and CPU time.
Algorithm 3.1. M-MMSGP algorithm for (2.5)
Step 0: We are given ε>0, z∈Rn. Let σ, β∈(0,1), α0=1, r,αmin∈[0,1]. Set k:=0.
Step 1: If ‖F(zk)‖⩽ε, then halt and consider zk as an approximate solution.
Step 2: Calculate search direction
Step 3: Find vk=zk+θkdk, where θk£=βmk, and mk is the smallest nonnegative integer such that
where γk£=‖F(zk)‖1+‖F(zk)‖.
Step 4: Compute the new iterate by
Step 5: Update the spectral vector by
where sk=zk+1−zk and yk=F(zk+1)−F(zk)+rsk.
Step 6: Increase k by 1 and return to Step 1.
Remark 3.1. If Algorithm 3.1 terminates at finite iteration k and F(zk)=0, then zk is the solution to F(z)=0. If k approaches infinity and F(zk)≠0, Algorithm 3.1 produces an infinite sequence {zk}. For simplicity, we can assume that {zk} is an infinite sequence.
Remark 3.2. (1) Based on the monotonic behavior of F(z), we can derive
(2) From the Lipschitz continuity of F(z), we can deduce that
(3) From (3.4)–(3.6), we can deduce that
Here, L=‖(B−A)−1(B+A)‖+1 is given by Theorem 2.3.
The subsequent lemma illustrates that Algorithm 3.1 is well-defined.
Lemma 3.1. According to the propositions outlined above, there exists a nonnegative integer mk ensuring that (3.2) is satisfied for all k⩾0.
Proof. Assume there is a k0⩾0 for which (3.2) fails to hold for any nonnegative integer m, i.e.,
By letting m→∞ and utilizing the continuity of F(z), we obtain
Based on Step 1 of Algorithm 3.1 and (3.1), we find that
for any k⩾0. Consequently, it follows that
and
where k⩾0. This result conflicts with (3.8), thus completing the proof. □
4.
Convergence analysis
In this part, we demonstrate that every convergent subsequence {zkj} of the sequence {zk} produced by Algorithm 3.1 approaches a solution of system (2.5). Due to the equivalence, it also converges to a solution of the HLCP (1.1).
Lemma 4.1. [15] Let F(z) be monotone, and z,v∈Rn such that F(v)T(z−v)>0. Let z∗ be a solution of F(z)=0 and
Then we have
Theorem 4.1. Given Assumption 2.1, assume the sequence {zk} is generated by the M-MMSGP method with ε=0. In this case, either the sequence {zk} is finite, and the final iterate {zk} is a solution to F(z)=0, or the sequence {zk} is infinite, and
Proof. By Assumption 2.1 and Theorem 2.2, F(z) in (2.4) is monotone. By (3.2) of Algorithm 3.1, we have
By Lemma 4.1, one has
where z∗ solves F(z)=0. Hence, {‖zk−z∗‖} is a decreasing sequence and ‖zk−z∗‖⩽‖z0−z∗‖ holds for all k⩾1 which implies that {zk}⊂{z∈Rn:‖z−z∗‖⩽‖z0−z∗‖}. Hence, {zk} is a bounded sequence.
Given the continuity of F(z), the sequence {‖F(zk)‖} is bounded.
By (4.2), {dk} is also bounded. From Step 3 of Algorithm 3.1, we have vk=zk+θkdk. Since {zk} and {dk} are bounded, θk∈(0,1), so {vk} is bounded. Given the continuity of F(v), the sequence {‖F(vk)‖} is bounded. There exists a positive constant M such that ‖F(vk)‖⩽M for all k⩾1.
Assuming without loss of generality that the sequence {zk} is infinite, from (4.1) we obtain
By (3.3) we have
Combining the above inequality and (4.3) yields limk→∞θk‖dk‖=0. This concludes the proof. □
The subsequent theorem illustrates the global convergence of the proposed M-MMSGP method.
Theorem 4.2. Under the condition of Assumption 2.1, suppose that sequence {zk} is generated by the proposed M-MMSGP method with ε=0. Then, {zk} converges to some point z∗ satisfying F(z∗)=0.
Proof. By Theorem 4.1, we get
(1) If lim infk→∞‖dk‖=0, then from (3.7), it follows that lim infk→∞‖F(zk)‖=0. Consequently, there exists a subsequence {zkj}⊂{zk} such that limkj→∞‖F(zkj)‖=0. Given the continuity of F(z), this implies limkj→∞‖zkj−z∗‖=0, where z∗ is a point satisfying F(z∗)=0.
(2) If lim infk→∞‖dk‖>0, then according to Theorem 4.1, limk→∞θk=0. From Eq (4.4), it follows that limk→∞‖F(zk)‖>0.
Based on Algorithm 3.1, we have
Given that the sequences {zk} and {dk} are bounded, they must have at least one cluster point. Therefore, there exist subsequences {zkj}⊂{zk} and {dkj}⊂{dk} that converge to this cluster point.
In the above equation, let k→∞, and we can get
where ˆz and ˆd are limit points of subsequences {zkj}⊂{zk} and {dkj}⊂ {dk}, respectively. On the other hand, let k→∞ in (3.9), and we have
This contradicts (4.5). Therefore, inequality lim infk→∞‖dk‖>0 does not hold. This completes the proof. □
5.
Numerical results
In this section, we present some numerical results to demonstrate the efficiency of Algorithm 3.1. We choose a modified gradient projection algorithm (MGP) [25] used to solve (2.5), a modified multivariate spectral gradient algorithm (MMSGP) [23] used to solve (2.5), a modulus-based Jacobi algorithm (MJ) [6] used to solve (2.1), and a modified spectral gradient projection-based algorithm (PGP) [18] used to solve (2.5) as the comparison. All methods were implemented in MATLAB R2018a and executed on a personal computer with an Intel Core i7 processor operating at 1.80 GHz and 8GB of RAM.
All experiment results include three aspects: the elapsed CPU time in seconds (CPU), the norm of absolute residual vectors (RES), and the number of iteration steps (IT), respectively. RES is defined as RES:=‖Ax(k)−By(k)−q‖.
In the following experiments, when the prescribed iteration number ITmax = 600 is exceeded or the residual vector satisfies RES⩽10−6, all runs are terminated. We will consider the problems with six dimensions, i.e., n=100, 400, 900, 1600, 2500, and 3600.
The primary numerical outcomes are presented in Tables 1–5, as well as Figures 1–4, to facilitate easy comparison. In these tables and figures, the algorithm parameters are set as follows:
(1) For the MGP algorithm, set λ0=0.4,δ=1.01,α=0.4δ;
(2) For the MMSGP algorithm, set α0=1,β=0.2,τ=0.001,σ=0.01,r=0.001;
(3) For the MJ algorithm, set Ω=0.5In;
(4) For the PGP algorithm, set β=1,ρ=0.8,σ=10−5,l=10−4,u=1030;
(5) For the M-MMSGP algorithm, set λmin=0.1ones(n,1),β=0.618,σ=0.01,r=0.001.
Example 5.1. [7] To solve the HLCP, consider the tridiagonal matrix
where m is a given positive integer. We define the matrices A,B∈Rn×n, with n=m2, as A=ˆA+μI and B=ˆB+νI with μ,v as the real parameters and
Let z∗=(1,−1,1,−1,…)T, x∗=|z∗|+z∗, and y∗=|z∗|−z∗, and then q=Ax∗−By∗.
Example 5.1 can be derived from the discretization form of a two-dimensional boundary problem, which is specifically described as
where z(u,v),w(u,v), and q(u,v) are three two-dimensional maps. Therefore, under appropriate boundary conditions, the boundary problem is discretized by the five-point difference discretization method, and its discretization scheme is in the form of an HLCP.
Example 5.2. [10] To solve the HLCP, consider the tridiagonal matrix
where m is a given positive integer. We define the matrices A,B∈Rn×n, with n=m2, as A=ˆA+μI and B=ˆB+νI with μ,v as the real parameters and
Let z∗=(1,−1,1,−1,…)T, x∗=|z∗|+z∗, and y∗=|z∗|−z∗, and then q=Ax∗−By∗.
Example 5.3. To solve the HLCP, where A∈Rn×n and B∈Rn×n, consider the block-tridiagonal matrices
and
with the tridiagonal matrix
Here, n=m2. Let z∗=(1,−1,1,−1,…)T, x∗=|z∗|+z∗, and y∗=|z∗|−z∗, and then q=Ax∗−By∗. This example is adapted from the literature [10].
Example 5.4. To solve the HLCP, where A∈Rn×n and B∈Rn×n, consider the block-tridiagonal matrices
and
with the tridiagonal matrix
Here, n=m2. Let z∗=(1,−1,1,−1,…)T, x∗=|z∗|+z∗, and y∗=|z∗|−z∗, and then q=Ax∗−By∗. This example is adapted from the literature [10].
Example 5.5. To solve the HLCP, where A∈Rn×n and A is a random matrix with eigenvalues {1,1,3,3,3,...}, consider
Here, n=m2. Let z∗=(1,−1,1,−1,…)T, x∗=|z∗|+z∗, and y∗=|z∗|−z∗, and then q=Ax∗−By∗.
For Examples 5.1–5.5, the numeric results are list in Tables 1–5, respectively.
Table 1 contains the results of the M-MMSGP method, MJ method, MGP method, and PGP method for Example 5.1 under n=100, 400, 900, 1600, 2500, and 3600, respectively. Here, we take μ=2, and ν=3. Among the four algorithms, it can be seen that the M-MMSGP algorithm has obvious advantages in both the number of iterations and the CPU time. The M-MMSGP algorithm has good numerical performance.
Tables 2 and 3 contain the results of the M-MMSGP method, MJ method, MGP method, and PGP method for Examples 5.2 and 5.3 under n=100, 400, 900, 1600, 2500, and 3600, respectively. The M-MMSGP algorithm has excellent numerical performance. The numerical results show that the M-MMSGP method excels in both CPU time and iteration count. Among the four algorithms, the MGP algorithm has the most iteration steps but consumes less CPU time. The MJ algorithm has the highest CPU time but fewer iteration steps. The MJ algorithm requires inversion in each iteration, while the other three algorithms do not require inversion, which requires more CPU time. Therefore, the MJ algorithm has the highest CPU time. The selection of gradient descent direction in the MGP algorithm is not effective, while the other three algorithms have better descent directions, so the MGP algorithm has more iteration steps. The PGP algorithm always spends more CPU time than the M-MMSGP algorithm although the number of iteration steps is less than the M-MMSGP algorithm in certain dimensions. The M-MMSGP algorithm only needs to compare and allocate spectral gradients, which provides better CPU time. Especially for the M-MMSGP algorithm, the larger the matrix dimension, the faster the convergence speed.
Table 4 contains the results of the M-MMSGP method, MJ method, MGP method, and MMSGP method for Example 5.4 under n=400, 900, 1600, 2500, and 3600, respectively. The M-MMSGP algorithm has excellent numerical performance. The numerical results indicate that the M-MMSGP method has the optimal CPU time. The iteration times of the M-MMMGP algorithm are only inferior to the MJ algorithm in certain dimensions. Overall, the M-MMSGP algorithm also has better numerical performance in terms of iteration times. The reason why the MJ algorithm takes more time than the M-MMSGP algorithm is that the MJ algorithm requires inversion in each iteration, which requires a certain amount of CPU time. The M-MMSGP algorithm only needs to compare and allocate spectral gradients, which provides better CPU time.
The experimental results demonstrate that despite having a lower number of iterations, the MMSGP algorithm consumes the most CPU time. Specifically, when n=900, the CPU time surpasses 100 seconds. The reason why MMSGP has fewer iterations is that the algorithm needs to assign values to each element in the gradient direction in order to find the optimal gradient descent direction. However, due to the allocation of each element, more CPU time is required when the dimensionality is high. The M-MMSGP algorithm only requires an overall comparison of the distribution spectral gradient, which occupies less CPU time. Especially for the M-MMSGP algorithm, the larger the matrix dimension, the faster the convergence speed.
From Table 5, it can be seen that the M-MMSGP method is sensitive to solving the HLCP. The numerical calculation results include the M-MMSGP algorithm, MJ algorithm, and MGP algorithm, which have good performance. The M-MMSGP method has the optimal CPU time and iteration times. The experimental data indicates that for the M-MMSGP algorithm, an increase in dimensionality does not significantly affect the number of iterations, suggesting robust numerical performance. Compared to the other two algorithms, the MJ algorithm outperforms the MGP algorithm, while the M-MMSGP algorithm surpasses both.
Figure 1 shows the relationship between the number of iterations and RES for Example 5.1. It can be seen that the RES descent speed of the M-MMSGP algorithm is faster than that of the MJ algorithm and MGP algorithm. This also verifies the effectiveness of the M-MMSGP algorithm. Figure 2 shows the relationship between the number of iterations and RES for Example 5.3. It can be seen that the RES drop speed of the M-MMSGP algorithm is the fastest among the three algorithms. Figure 3 shows the relationship between the number of iterations and RES for Example 5.4. It can be seen that the RES drop speed of the M-MMSGP algorithm is the fastest among the three algorithms. The relationship between the number of iterations and RES in Example 5.5 is shown in Figure 4. The above proves that the M-MMSGP algorithm has good numerical performance.
The aforementioned numerical results highlight that the M-MMSGP algorithm demonstrates greater efficiency in terms of CPU time and iteration count under specific conditions. Therefore, our proposed method could be well-suited for solving the HLCP.
The M-MMSGP algorithm employs a novel linear search method, which is distinct from existing gradient algorithms. It utilizes multivariate spectral gradient vectors as the descent direction and allocates the multivariate spectral gradients as a whole, resulting in a significant increase in the number of iterations and CPU time required by the algorithm.
6.
Conclusions
The focus of this study is the HLCP. First, the HLCP is reconstructed into a fixed-point equation based on the modulus defined in the nonnegative cone. By rewriting the fixed-point equation as a monotone system, we put forward a new class of modified multivariate spectral gradient projection method for solving it. It is shown that the proposed iterative method converges to the HLCP solution, assuming the specified conditions are met. Furthermore, the viability and efficacy of the modulus-based modified multivariate spectral gradient projection method are illustrated through numerical examples. The algorithm presented in this paper has been developed to address the horizontal complementarity problem within the theoretical framework of this study. Further efforts could be directed toward solving this problem under more simplified assumptions.
Author contributions
Ting Lin: Conceptualization, Methodology, Validation, Formal analysis, Resources, Writing–review and editing, Supervision, Project administration; Hong Zhang: Data curation, Visualization; Chaofan Xie: Conceptualization, Software, Validation, Data curation, Writing–original draft preparation, Visualization. All authors have read and agreed to the published version of the manuscript.
Use of Generative-AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
This research was funded by the National Natural Science Foundation of China under Grant 62071123 and 61601125, Natural Science Foundation of Fujian Province of China (number: 2023J011117), the Fujian Province Education Hall Youth Project (number: JAT220258), and the Fujian Natural Science Foundation Project (number: 2019J01887).
Conflict of interest
The authors declare no conflicts of interest.