A new spectral method with inertial technique for solving system of nonlinear monotone equations and applications

1.   Introduction

Spectral gradient methods are among the widely known first-order methods for unconstrained optimization problems $\min_{m \in \mathbb{R}^n} f(m),$ where $f:\mathbb{R}^n\to \mathbb{R}$ is a smooth nonlinear function that is bounded below. These methods generate a sequence of approximations using the iterative formula:

where $\Lambda_k$ is called the step size and is obtained using some line search procedures, and $t_k$ is called a search direction, defined as

where $g_k: = \nabla f(m_k),$ and the coefficient $\overline{\gamma}_k$ is a scalar known as the spectral parameter, which differentiates any two spectral gradient methods (see Barzilai and Borwein (BB) ^[1]). The classical forms of the parameter $\overline{\gamma}_k$ given in ^[1] are

Among the advantages of the spectral gradient method is its simplicity in implementation and low storage requirement; thus, it is suitable for large-scale problems. As a result, researchers have extended the spectral gradient method to solve systems of nonlinear equations (see, ^[2,3,4]). A system of nonlinear equations involves finding a vector $m$ in a nonempty, closed and convex set $C \subset \mathbb{R}^n$ such that

where $\Gamma:\mathbb{R}^n\to \mathbb{R}^n$ is continuous. When the mapping $\Gamma$ is monotone, i.e.,

Problem (1.5) becomes a system of nonlinear monotone equations. Many problems arising from science and engineering can be translated into the form of problem (1.5). For example, applications of system of nonlinear equations (1.5) have appeared in different fields ^{[5,6,7,8,9,10]}. In recent years, algorithms for solving systems of nonlinear monotone equations have been used in signal and image recovery, (see ^{[11,12,13,14,15,16,17]}).

Motivated by the hyperplane projection strategy of Solodov and Svaiter ^[18], spectral gradient-like methods for solving (1.5) have gained more attention. Zhang and Zhou ^[19] developed a spectral gradient projection method for unconstrained nonlinear equations based on the modified version of the $\overline{\gamma}_k^{long}$ (1.3). Later on, Yu et al.^[20] extended the work in ^[19] to solve convex constrained nonlinear equations. In addition, Yu et al. ^[21] proposed another spectral method for solving a system of nonlinear equations. The search direction in their work uses a convex combination of the modified $\overline{\gamma}_k^{long}$ (1.3) and $\overline{\gamma}_k^{short}$ (1.4). Their work revealed that combining the modified BB parameters gives better numerical performance than deploying them separately. This means that the efficiency of the convex combination of the $\overline{\gamma}_k^{long}$ (1.3) and $\overline{\gamma}_k^{short}$ (1.4) can be further explored.

Nowadays, there is a growing interest in incorporating the inertial technique into algorithms for solving systems of nonlinear equations (see, for example, ^[22,23]). The inertial step defined as $i_{k+1}: = m_{k+1}+\alpha_k(m_{k+1}-m_{k}), \; \alpha_k\in(0, 1),$ was proposed by Polyak ^[24] in order to speed up the performance of iterative algorithms. It can be observed that the inertial technique is using two previous iterates to compute the current iterate. This has also been shown to accelerate the iteration process of algorithms for solving nonlinear problems such as the proximal point method ^[25,26] and auxiliary problem principle ^[27].

Motivated by the contributions of the above-mentioned literature, this work seeks to explore the effect of the inertial technique on the convex combination of the $\overline{\gamma}_k^{long}$ and $\overline{\gamma}_k^{short}$ based on some modifications. This idea can be viewed as the modification of the work of Yu et al. ^[21]. Numerical experiments conducted in Section 3 reveal some level of improvement in numerical performance. Some of the notable contributions of this work include the following:

● A new spectral method for solving a system of nonlinear monotone equations based on the inertial technique is proposed.

● This work generalizes some existing algorithms in the literature.

● The global convergence of the proposed method is discussed under standard conditions.

● To depict the efficiency of the new method, a numerical experiment on a collection of test problems in comparison with some existing methods is presented.

● Subsequently, the proposed algorithm is applied to problems arising from robotic motion control.

In the remaining part of this work, the next section gives some definitions, details of the algorithms and global convergence. Section three gives some numerical experiments, while the application part is in Section four. In the final section, some concluding remarks are given.

2.   Proposed algorithm and convergence analysis

We begin this section by recalling the definition of the projection operator as follows:

Definition 2.1. Suppose $C\subset \mathbb{R}^n$ is a convex, nonempty and closed set. Then, any point $m\in \mathbb{R}^n$ can be projected onto the set $C$ using

The relation (2.1) satisfies the following useful property:

Definition 2.2. Any vector-valued map that satisfies

is said to be Lipschitz continuous.

In what follows, we present the proposed algorithm and subsequently give some remarks.

Remark 2.3. The choice of the $\theta_k$ defined by (2.10) is prompted by the work of Awwal et al. ^[28]. Now, since $\max \{\|\Gamma(m_{k-1})\|, \|\Gamma(m_k)\|\}\geq \|\Gamma(m_k)\|,$ then by the Cauchy-Schwarz inequality we have

and by the fact that $\mu\in(0, 1),$ we get

Remark 2.4. We note that if for all $k$, $i_k = m_k$ (2.8) and (2.9) and $\theta_k\in[0, 1],$ then the search direction (2.6) reduces to that of Yu et al. ^[21]. Moreover, if for all $k$, $i_k = m_k$ in (2.8) and $\theta_k = 0,$ then the search direction (2.6) reduces to that of Yu et al.^[20]. Furthermore, the new search direction (2.6) also reduces to DAIS1 and DAIS2 of Awwal et al. in ^[22] if $\theta_k = 1$ and $\theta_k = 0,$ respectively. Hence, our proposed work is regarded as an extension of the algorithms proposed in ^[21,22].

To show the global convergence of the proposed algorithms, we assumed the following:

($A_{1}$) The mapping $\Gamma$ is monotone.

($A_{2}$) The mapping $\Gamma$ is Lipschitz continuous.

($A_{3}$) The solution set of problem (1.5) is nonempty.

Lemma 2.5. Suppose that the assumptions ($A_{1}$)–($A_{3}$) hold, and that $\{m_k\}$ and $\{i_k\}$ are produced by Algorithm 1. Then,

(ⅰ) $\lim\limits_{k \to \infty }\|m_{k}-\widehat{m}\|$ exists.

(ⅱ) $\{m_k\}$, $\{i_k\}$ and $\|\Gamma(m_k)\|$ are bounded.

(ⅲ) The search direction is bounded, i.e.,

(ⅳ) The search direction satisfies

(ⅴ)

Proof. Let $\widehat{m}$ be a solution of problem (1.5), and by the assumption ($A_{1}$), we have

Now, since $0 < \eta < 2,$ from (2.2), (2.5) and (2.17) we have

This means that $\|m_k-\widehat{m}\|\leq \|m_{k-1}-\widehat{m}\|\leq \cdots\leq \|m_0-\widehat{m}\|,$ and thus $\lim\limits_{k\to\infty}\|m_k-\widehat{m}\|$ exists.

(ⅱ) Since $\lim\limits_{k\to\infty}\|m_k-\widehat{m}\|$ exists, $\{m_k\}$ is bounded. Combined with the fact that $0 < \alpha_k < 1, \; \forall k,$ this gives the boundedness of $\{i_k\}.$

In addition, since the mapping $\Gamma$ is Lipschitz continuous and $\{m_k\}$ is bounded, we can find a positive constant $c_3 > 0$ such that

(ⅲ) To prove that the search direction defined by (2.6) is bounded, we need to show that the parameter $\gamma_k$ (2.7) is bounded.

From assumption ($A_{1}$), the mapping $\Gamma$ is monotone. This gives us $\langle \Gamma(i_{k+1}) - \Gamma(i_k), \; i_{k+1}-i_k \rangle \geq 0,$ and thus

Moreover, using the assumption that the mapping $\Gamma$ is Lipschitz continuous, together with the Cauchy Schwarz inequality, we have

Thus, (2.20) and (2.21) imply

and therefore

On the other hand,

Using Lipschitz continuity, we have

Combining (2.24) and (2.25) gives

Using (2.22) and (2.26), we have

Therefore, setting $M = \frac{1}{r} + \frac{L+r}{r^2}$ yields

Since $\forall k,$ $\theta_k\leq 1$ (see Remark 1),

Combining this with (2.19) gives (2.14) with $c_1 = Mc_3.$

(ⅳ) By the definition of $\theta_k,$ we have three possibilities, and thus the parameter $\gamma_k$ may take any of the three different following forms: $\gamma_k = \overline{\beta}_k,$ $\gamma_k = \widehat{\beta}_k$ and $\gamma_k = (1-\theta_k) \overline{\beta}_k +\theta_k\widehat{\beta}_k$ for $\theta_k = 0,$ $\theta_k = 1$ and $\theta_k\in(0, 1),$ respectively. Therefore, we divide this proof into three cases:

Case Ⅰ: If $\theta = 0, \; \forall k,$ then the search direction (2.6) reduces to $t_k = -\overline{\beta}_k\Gamma_k(m_k),$ and therefore, using (2.23) gives

Case Ⅱ: If $\theta = 1, \; \forall k,$ then the search direction (2.6) becomes $t_k = -\widehat{\beta}_k\Gamma_k(m_k),$ and thus, using (2.27) yields

Case Ⅲ: If $0 < \theta < 1, \; \forall k,$ then we can find some constant $c_4 > 0$ such that $\theta_k > c_4$ and $(1-\theta_k) > 0.$ Therefore, from (2.7) and (2.27), we have

Thus, from the search direction (2.6) and (2.32), it holds that

Hence, from the three cases above, we see that (2.15) holds.

(ⅴ) Using the boundedness of $\{m_k\}$, (2.14) and the definition of $p_k$ in (2.3), $\{p_k\}$ is bounded. Also, using the assumption that $\Gamma$ is Lipschitz continuous, we get

Since $\min\left\{1, \|\Gamma(m_k + \Lambda_k t_k)\|^{\frac{1}{c}} \right\} \leq 1,$ squaring from both sides of (2.4) yields

Furthermore, since $0 < \eta < 2,$ from (2.18) we obtain

This together with (2.35) gives

Recall that $\|\Gamma(p_k)\|$ is bounded by $n_1$ (see (2.34)), and we obtain

Since $\lim\limits_{k\to \infty}\|m_{k}-\widehat{m}\|$ exists, taking the limit as $k \to \infty$ on both sides of (2.38) gives

which implies (2.16).

Lemma 2.6. Suppose that the Assumption ($A_{2}$) holds. Let the sequences $\{m_k\}$ and $\{p_k\}$ be generated by Algorithm 1. Then,

Proof. From (2.4), if $\Lambda_k\neq \kappa$, then $\tilde{\Lambda}_{k} = \Lambda_k\varsigma^{-1}$ violates (2.4), that is,

We know that $\min \{ 1, \|\Gamma(m_k+\varsigma^jt_k)\|^\frac{1}{c}\} \leq 1.$ Thus, from (2.15) and Assumption ($A_{2}$), we get

Therefore,

Substituting $\tilde{\Lambda}_{k} = \Lambda_k \varsigma^{-1}$ in (2.40) and solving for $\Lambda_k,$ we get

Thus, we have

Theorem 2.7. If the Assumptions $(A_{1}-A_{3})$ hold, and the sequence $\{m_k\}$ is produced by Algorithm 1, then

Proof. We show the proof by contradiction. Suppose (2.42) does not hold, and then there exists $s > 0$ such that $\forall k\geq0,$

From Eqs (2.15) and (2.43), we get $\forall k\geq0,$

Multiplying $\|t_k\|$ on both sides of (2.39), and using (2.14) and (2.43), we obtain

Taking the limit as $k\to \infty$ on both sides gives

which contradicts (2.16). Hence, $\liminf\limits_{k\to \infty}\|\Gamma(m_k)\| = 0.$

3.   Numerical experiments

In this section, we present the numerical experiments performed by solving a set of test problems taken from the literature. To depict the efficiency of the iSDFM algorithm, we perform numerical comparison with two other existing methods. The first one is the DAIS1 algorithm proposed in ^[22], which is an inertial-based algorithm for solving a system of nonlinear equations. The second algorithm is the MSGPALG proposed by Yu et al. in ^[21] based on the convex combination of the modified BB long and short parameters. As noted in Remark 2.4, both DAIS1 and MSGPALG can be viewed as special cases of the proposed iSDFM algorithm. We test the performances of all these three algorithms on seven test problems with eight different initial points (see Table 1) and five dimensions (1000, 5000, 10000, 50000, and 100000), thus making the total number of the test problems 280. These three algorithms are coded on MATLAB R2019b which runs on a PC of corei3-4005U processor with 4 GB RAM and 1.70 GHz CPU. The choice of parameters for MSGPALG and DAIS1 are maintained as reported in their respective references [21,22]. In the iSDFM algorithm, we choose $\varsigma = 0.47,$ $\alpha_k = 1/(k+1)^2,$ $\eta = 1.79,$ $\mu = 0.5,$ $\sigma = 0.01,$ $r = 0.001,$ $c = 2$ and $\kappa = 1$. The stopping criterion is set to be $\|\Gamma(m_k)\| < 10^{-6}.$

We consider the following test problems, where $\Gamma (m) = (g_1(m), g_2(m), \ldots, g_n(m))^T$:

Problem 1. ^[29]:

Problem 2. ^[28] Modified logarithmic function:

Problem 3. ^[30] Nonsmooth function:

Problem 4. ^[31] Strictly convex function:

Problem 5. ^[20] Nonsmooth function:

Problem 6.^[29]:

Problem 7. ^[32]:

Based on these settings, the results of the experiments are tabulated in Tables 2–8 with ITER, FVAL and TIME denoting the number of iterations, number of function evaluations and CPU time, respectively. Based on these metrics, it can be observed that the performances of the three algorithms varies in terms of the ITER, FVAL and TIME. However, taking the whole results of the experiment into consideration, it can be seen that the proposed iSDFM algorithm outperformed the DAIS1 and the MSGPALG algorithms in most instances. By outperforming we mean the new iSDFM recorded the least ITER, FVAL and TIME in most cases of the experiment.

With the help of the Dolan and Moré performance profile ^[33], we present the information in Tables 2–8 graphically for better and easier visualization of each algorithm's performance. These graphs are plotted in Figures 1–3. It can be clearly observed from Figure 1 that the iSDFM algorithm solved about 72% of the problems with the least ITER, as compared to the DAIS1 and MSGPALG with around 30% and 12%, respectively. Figure 2 shows that the iSDFM outperformed the other two methods by solving around 70% of the problems with the least FVAL. In terms of TIME, iSDFM algorithm competes favorably with the MSGPALG algorithm. In general, our proposed iSDFM algorithm shows better efficiency as compared to the DAIS1 and MSGPALG algorithms. This might not be unconnected with taking the convex combination of the BB-like parameters incorporated with the inertial technique.

4.   Application in motion control

In robotics, manipulators and effectors are the aspects in which some parts of the robots interact with other objects by performing different tasks such as picking from one point and placing on another. For stability and accuracy in the robots' movement, the characteristics of motor dynamics, which are contemporarily used as actuators in $n-$link and $1-$link robot systems, need to be considered ^[34,35]. This is a tracking control problem of a nonlinear system, and the motor dynamics are required to satisfy the condition that the actual output of the system can track the desired trajectory with least possible error ^[36].

Some of the developed methods for tracking control problems of nonlinear systems include proportional-integral-derivative (PID) control ^[37,38], feedback linearization ^[39,40] and optimal output tracking control by using approximation approach ^[41].

We present an application of our proposed algorithm in motion control of two planar robotic manipulators. Consider the following model:

where $f: \mathbb{R}^n \to \mathbb{R}$ is a smooth and convex function. In the case of motion control, the function in (4.1) has the form $f(v_k): = \frac{1}{2}\|\upsilon_k-y_{uk}\|_{2}^2.$ Subsequently, we minimize

at each computational time interval $\tau_k\in [0, \tau_f].$

As described in ^[42], the discrete-time kinematics equation of a two-joint planar robot manipulator at the position level is given as

The kinematics map $\psi(\cdot)$ is given as

where $l_1$ and $l_2$ are the lengths of the rod links, and $\theta_k \in \mathbb{R}^2$ is the joint angle vector.

From (4.2), the term $\upsilon_k$ is controlled to track a Lissajous curve:

For our algorithm to fit (4.1), we present its slight modification as follows:

Remark 4.1. Assume that the solution to problem (4.1) exists and that the level set $\{m\in \mathbb{R}^n: f(m)\leq f(m_0)\}$ is bounded. Using Assumption $A_{2}$ and Theorem 2 of ^[43], we conclude that Algorithm 2 converges, that is, $\liminf_{k\to \infty}\|\Gamma(m_k)\| = 0$ holds.

In this experiment, we chose initial joint angle vector $\theta_0 = [0, \frac{\pi}{3}]^T,$ $l_1 = l_2 = 1,$ $\varsigma = 0.2$ and $\sigma = 0.08.$ The task duration $[0, 20]$ is subdivided into 200 equal parts. The performance of Algorithm 2 in the motion control problem is shown in the Figure 4.

Figure 4(a) represents the robot trajectories synthesized by Algorithm 2. Figure 4(b) shows the end effector and the desired path. Figure 4(c, d) shows the performance errors on x and y axes, respectively. These figures indicate that Algorithm 2 efficiently performed the assignment with an error as low as $10^{-5}$ on both the x-axis and y-axis.

5.   Conclusions

In conclusion, we have proposed an inertial-based spectral method for solving a system of nonlinear equations. The inertial-step introduced is believed to have enhanced the performance of the new method. We have discussed the global convergence of the proposed algorithm under the monotonicity and Lipschitz continuity assumptions. We have also presented some numerical experiments which depicted the efficiency of the proposed algorithm. The proposed algorithm is reported to have better numerical performance than the methods in ^[21,22]. Subsequently, we have demonstrated the applicability of the new method in motion control problems arising from robotics. However, the numerical results showed that the proposed method won about 72% and 70% of the experiments in terms of ITER and FVAL. Therefore, we recommend further research in order to improve its performance to the optimum.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and suggestions. The sixth author was funded by Chiang Mai University and the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation (grant number B05F640183).

Conflict of interest

The authors declare that they have no conflicts of interest.