An efficient DY-type spectral conjugate gradient method for system of nonlinear monotone equations with application in signal recovery

1.   Introduction

Consider a nonempty closed convex set $\Omega \subset \mathbb{R}^n$ and a mapping $F:\mathbb{R}^n\to \mathbb{R}^n$ which is continuous and monotone. In this work, our interest is on finding the solution $x \in \Omega$ such that

System of nonlinear equation (1.1) arises in many practical applications such as in controlling the motion of a planar robot manipulator^[1], economic equilibrium problems ^[2], power flow equations ^[3] and chemical equilibrium systems ^[4]. Additionally, in mathematics, subproblems in the generalized proximal algorithms with Bregman distance ^[5] and monotone variational inequality problems by using fixed point map or normal map ^[6,7] can all be transformed into finding the solution of (1.1). Recently, algorithms for solving system of nonlinear monotone equations are proved to be efficient in signal and image recovery ^[8].

Due to the existence of nonlinear equations in various fields and their wide range of applications, a lot of methods have been proposed to find their solution. Some of the early and popular iterative methods are Newton method, Quasi-Newton method, Levenberg-Marquardt method and their variants ^{[9,10,11,12,13]}. These early methods are characterised by their advantage of rapid convergence from good initial points. However, they require solving linear systems using a Jacobian matrix or its approximation at every iteration. This problem affects their suitability to solve large-scale systems of nonlinear equations.

On the other hand, conjugate gradient methods, spectral gradient methods, and spectral conjugate gradient methods are class of methods for solving large scale unconstrained optimization problems. Among the advantages of these methods are their simplicity in implementation and low storage requirements. These advantages stimulate researchers to extend these methods in order to solve nonlinear equations. For example, the projection technique proposed by Solodov and Svaiter ^[14] motivated many researchers to extend conjugate gradient methods from solving unconstrained optimization problem to solve system of nonlinear equations. Inspired by the work of Solodov and Svaiter in ^[14], Wang et al. ^[15] proposed a projection method for solving system of nonlinear monotone equations. In their algorithm, a linear system of equations is solved approximately, at each iteration, to obtain a trial point and then a line search strategy is performed along the search direction determined by the current point and the trial point with the aim of getting a predictor-corrector point. Hence, the algorithm computes its next iterate by projection. They proved the global convergence of the proposed method and presented some numerical experiments in order to show the performance of the algorithm. In ^[16], Cheng combined the classical PRP method with hyperplane projection method to solve nonlinear monotone equations. Xiao and Zhou ^[17] extended the well-known CG_Descent for unconstrained minimization problems to solve large scale convex constraint nonlinear monotone equations. They achieved this by combining the CG_Descent with the projection technique in ^[14]. They proved the global convergence of this method and showed the numerical performance. Liu and Li ^[18] modified the work in ^[17] and proposed another extension of the CG_descent method to solve nonlinear monotone equations. Based on the popular Dai-Yuan (DY) conjugate gradient parameter ^[19], Liu ^[20] proposed a spectral DY-type method for solving nonlinear monotone equations. The method can be viewed as a combination of the DY conjugate gradient method, spectral gradient method and the projection technique. They showed that the method converges globally and presented some numerical experiments. Later on, Liu and Li ^[21] developed another DY-type algorithm, which is a multivariate spectral method for solving (1.1). In their work, the direction uses a combination of the multivariate spectral gradient method and DY conjugate gradient parameter. The numerical experiments of the method is reported and the global convergence is also proved. However, restriction is imposed on the lipschitz constant $L < 1 - r$ with $r \in (0, 1)$ before proving the global convergence. Motivated by this work, Liu and Feng ^[22] proposed another spectral conjugate gradient method for solving (1.1). Their work improved the computational effect of the DY conjugate gradient method and under some assumptions both the global and linear convergence of the method is proved. Most recently, a lot of algorithms have been developed for solving (1.1). Some of these algorithms can be found in ^{[23,24,25,26,27,28,29,30,31,32]}.

In this paper, motivated by the work of Liu and Feng ^[22] on the modification of DY conjugate gradient method, we propose an efficient spectral conjugate gradient algorithm for solving systems of nonlinear monotone equations with convex constraint. The search direction in our proposed approach uses a convex combination of the DY parameter and a modified CD parameter. Specifically, this paper gives the following contributions:

● We propose an efficient spectral conjugate gradient algorithm for solving systems of nonlinear monotone equations with convex constraint by taking the convex combination of the DY parameter and a modified CD parameter.

● This algorithm can be viewed as an extension of the work proposed by Yu et al in ^[33].

● The global convergence of the proposed algorithm is proved under some suitable assumptions.

● The proposed algorithm is applied to recover a distorted signal.

The organization of the paper is as follows: In the next section, we introduce the details of the algorithm, some important definitions and prove global convergence. In the third section, we provide numerical experiments of the proposed algorithm and compare it performance with an existing one. Finally, we apply the algorithm in signal recovery, and give conclusion in the last section.

2.   Algorithm: Motivation and convergence result

In this section, projection map, its properties, and some important assumptions needed for the convergence analysis are introduced.

Definition 2.1. Let $\Omega\subset \mathbb{R}^n$ be a nonempty, closed and convex set. The projection of any $x\in \mathbb{R}^n$ onto $\Omega$ is

The projection map has the following property:

Spectral and conjugate gradient algorithms generate sequence of iterates using the following formula:

where $\alpha_k$ is called the step length, and the direction $d_k$ defined in spectral and conjugate gradient method respectively as:

and

Different spectral and conjugate gradient directions are developed using different choice of the parameters $\nu_k$ and $\beta_k$ respectively. To ensure the global convergence of these methods, the direction $d_k$ needs to satisfy the sufficient decent property. That is:

where $\tau > 0.$

One of the well-known parameter proposed in this direction is the DY conjugate gradient parameter in ^[19] defined as

such that the direction in (2.4) becomes:

where $Y_{k-1} = F_k - F_{k-1}.$ Unfortunately, (2.7) does not satisfy the decency property (2.5). As a result, Liu and Feng ^[22] modified the work in ^[19] and proposed a spectral conjugate gradient algorithm for solving (1.1). The direction in ^[22] satisfies (2.5) and thus, the global convergence is proved successfully under some appropriate assumptions.

Motivated by the work in ^[22], and due to the limited number of DY-type conjugate gradient methods in literature, we proposed a new spectral DY-type conjugate gradient algorithm for solving (1.1). Interestingly, the proposed direction satisfies the sufficient descent property (2.5), and uses a convex combination of the DY parameter and a modified CD parameter as follows:

where

$\beta_k^{DY}$ is defined as (2.6), $\tilde{\beta_k} = \frac{\|F_k\|^2}{\max\{ -F_{k}^Td_{k-1}\;, \gamma \|d_{k-1}\;\|\}}$ and $\theta_k \in (0, 1).$ Substituting the value of $\beta_k^{DY}$ and $\tilde{\beta_k}$ in (2.8) we get

Throughout this work, the following assumptions are made.

● ($Assumption_{1}$) The mapping $F$ is monotone, that is,

● ($Assumption_{2}$) The mapping $F$ is Lipschitz continuous, that is there exists $L > 0$ such that

● ($Assumption_{3}$) The solution set of (1.1), denoted by $\Omega$, is nonempty.

We state the steps of our proposed algorithm as follows:

Algorithm 2.2. Step 0. Choose initial point $x_0\in \Omega$, $\theta_k \in(0, 1), \kappa \in (0, 1], \beta \in(0, 1)\; \mu > 1, \; \sigma, \; \gamma > 0$, $\delta \in (0, 2),$ and $Tol > 0$. Set $k: = 0$.

Step 1. If $\|F_k\|\leq Tol$, stop, otherwise proceed with Step 2.

Step 2. Compute $d_k = -\|F_k\|, \; k = 0$ and

Step 3. Compute $\Lambda_k = \max\{\kappa \beta^{i}\; : i = 0, 1, 2, \ldots\}$ such that

Step 4. Set $z_k = x_k +\Lambda_k d_k$. If $\|F(z_k)\| = 0,$ stop. Else compute

Step 5. Let $k = k+1$ and go to Step 1.

Remark 2.3. It is worth noting that when $Y_{k-1}^Td_{k-1} \leq \mu \|F_k\| \|d_{k-1}\|,$ our proposed search direction reduces to that of Yu et al. proposed in ^[33]. Thus, as a contribution, this work can be viewed as an extension of the work of Yu et al. ^[33].

We now state and proof the following Lemmas and Theorem for the convergence.

Lemma 2.4. The parameter $\nu_k$ given by (2.9) is well defined, and $\forall k \geq 0, \; d_k$ satisfies

Proof. Since F is monotone, then

which yields

Again, by Lipschitz continuity, we have

From (2.14) and (2.15) we get

Now, to show (2.13), for $k = 0, \; F_k^Td_k = -\|F_k\|^2,$ thus $\tau = 1$ and the result holds. When $k \neq 0$, If $Y_{k-1}^Td_{k-1} \leq \mu \|F_k\| \|d_{k-1}\|,$ then from (2.11),

using (2.16), we have

and (2.13) holds by taking $\tau = \frac{1}{L+r}.$

On the other hand, if $Y_{k-1}^Td_{k-1} > \mu \|F_k\| \|d_{k-1}\|,$ multiplying (2.11) by $F_k^T$ we obtain

The last inequality follows from (2.16). By letting $\tau = (\frac{1}{L+r}-\frac{1}{\mu}) > 0,$ the required result holds.

Lemma 2.5. Let $\{d_k\}$ be given by (2.11), then there are some constants $p_1 > 0, \; m_1 > 0$ and $m_2 > 0$ for which

Proof. If $Y_{k-1}^Td_{k-1} \leq \mu\|F_k\| \|d_{k-1}\|,$

Using (2.16), we have

where $p_1 = \dfrac{1}{r}.$ However, if $Y_{k-1}^Td_{k-1} > \mu \|F_k\| \|d_{k-1}\|,$ then

where $m_1 = p_1 + \frac{1}{\mu}, \; p_1 = \frac{1}{r}$ and $m_2 = \frac{1}{\gamma}.$

Lemma 2.6. Suppose ($Assumption_{1}$) - ($Assumption_{3}$) hold, then the sequences $\{x_k\}$ and $\{z_k\}$ generated by Algorithm 2.2 are bounded. Also,

and

Proof. Let $\tilde{x}$ be a solution of problem (1.1), using monotonicity from Assumption 1 we have

Using the above Eq (2.22) and $x_{k+1}$ we obtain

Showing that $\|x_k-\tilde{x}\|\leq \|x_0-\tilde{x}\|$ for all $k$ and hence $\{x_k\}$ is bounded and $\lim\limits_{k\to\infty}\|x_k-\tilde{x}\|$ exists. Since $\{x_k\}$ {is} bounded, and $F$ is Lipschitz continuous,

Using this and (2.18), we have

where $n_1 = p_1p$ and $n_2 = m_1p + m_2 p^2$ and taking $M = \min \{n_1, n_2\},$ we have that the direction $d_k$ is bounded. That is

To prove that $\{z_k\}$ is bounded, we know that

and since we have proved that $d_k$ is bounded. This implies $\{z_k\}$ is also bounded. Again, by Lipschitz continuity,

Now from our line search (2.12), let $\min\{1, \; \|F(x_k + \kappa\beta^{i} d_k)\|^{\frac{1}{c}}\} = \|F(x_k + \kappa\beta^{i} d_k)\|^{\frac{1}{c}},$ squaring from both sides of (2.12) we get

Also, since $0 < \delta < 2,$ then from (2.23) we have

This together with (2.28) gives

Since $\lim\limits_{k\to \infty}\|x_{k}-\tilde{x}\|$ exists and that (2.27) holds, taking the limit as $k \to \infty$ on both sides of (2.30) we have

but $\|F(z_k)\| \neq 0,$ therefore,

Note that if the $\min\{1, \|F(x_k + \Lambda_k d_k)\|^{\frac{1}{c}} \} = 1,$ then, (2.30) becomes

and thus (2.32) holds.

Using this and the definition of $z_k,$ we obtain

From the definition of projection operation, we get

Lemma 2.7. Suppose ($Assumption_{2}$) holds, and the sequences $\{x_k\}$ and $\{z_k\}$ are generated by Algorithm 2.2. Then

Proof. From (2.12), if $\Lambda_k\neq \kappa$, then $\widehat{\Lambda}_{k} = \Lambda_k\beta^{-1}$ does not satisfy (2.12), that is,

Now let the $\min \{ 1, \|F(x_k+ \widehat{\Lambda}_{k} d_k)\|^\frac{1}{c}\} = \|F(x_k+\widehat{\Lambda}_{k} d_k)\|^\frac{1}{c}.$ Using (2.13) and ($Assumption_{2}$), we have

Therefore,

substituting $\widehat{\Lambda}_{k} = \Lambda_k \beta^{-1}$ and solving for $\Lambda_k$ we get

On the other hand, if $\min \{ 1, \|F(x_k+\widehat{\Lambda}_{k} d_k)\|^\frac{1}{c}\} = 1,$ then (2.38) reduces to

Combining (2.38) and (2.39), we get

Theorem 2.8. Suppose that ($Assumption_{1}$-$Assumption_{3}$) hold and let the sequence $\{x_k\}$ be generated by Algorithm 2.2, then

Proof. We prove by contradiction. Suppose (2.41) is not satisfied, then there {exists} $\alpha > 0$ such that $\forall k\geq0,$

From Eqs (2.13) and (2.42), we obtain $\forall k\geq0,$

We multiply $\|d_k\|$ on both sides of (2.36), and from (2.26) and (2.42), we get

Taking limit as $k\to \infty$ on both sides, we obtain

which contradicts Eq (2.32). Therefore,

3.   Numerical experiments

In this section, we give the numerical experiments in order to depict the advantages and the performance of our proposed algorithm (MDY) in comparison with the projected Dai-Yuan derivative-free algorithm (PDY) by Liu and Feng ^[22]. All codes are written on Matlab R2019b and run on a PC of corei3-4005U processor, 4 GB RAM and 1.70 GHZ CPU.

In MDY, the parameters are choosen as follows: $r = 0.001, \; \theta_k = 1/(k+1), \; \mu = 1.9, \; \gamma = 0.9, \; \sigma = 0.02, \; c = 2, \; \kappa = 1, \; \beta = 0.70$ and $\delta = 1.1$. The parameters in the PDY algorithm are maintained as exactly as they are reported in ^[22]. Based on this setting, we consider nine test problems with eight different initial points and tested them on five different dimensions, $n = 1000, n = 5000, n = 10000, n = 50000\; \text{ and}\; n = 100000.$ We used $\|F_k\| < 10^{-6}$ as stopping criteria and denoted failure by "-" whenever the number of iterations exceeds 1000 and the stopping criterion is not satisfied. The test problems are listed below, where the function $F$ is taken as $F(x) = (f_1(x), f_2(x), \ldots, f_n(x))^T$.

Problem 1 ^[26].

Problem 2 ^[34] Modified Logarithmic Function.

Problem 3 ^[35] Nonsmooth Function.

Problem 4 ^[36]

Problem 5 ^[37] Strictly Convex Function.

Problem 6

Problem 7 ^[38] Tridiagonal Exponential Function

Problem 8 ^[22]

Problem 9 ^[26]

The results of our experiments are shown in Tables 1–9 based on the number of iterations denoted as (ITER), number of function evaluations (FVAL), CPU time (TIME), and the norm of the function (NORM) when the solution was obtained. Looking at the reported results, it can be observed that the proposed MDY {algorithm} outperformed the PDY algorithm in most of the problems by having the least ITER, FVAL and TIME.

In addition, to further visualize the comparison of the MDY algorithm with the PDY algorithm graphically, we adopt the well- known Dolan and Morè performance profile ^[39] as reported in Figures 1–3. From Figures 1 and 2, it can be seen that the MDY algorithm outperformed the PDY algorithm significantly. The MDY algorithm solves about 93% of the problems considered with least ITER and 99% of the problems with least FVAL as opposed to the PDY method with about 28% and 3% problems for the ITER and FVAL respectively. Moreover, in terms of TIME, Figure 3 indicated that the MDY algorithm still performs much better than the PDY algorithm by solving around 88% of the problems in lesser time. From these figures, we can conclude that the numerical performance of the MDY algorithm has great advantage when compared with the existing PDY algorithm.

4.   Application in compressive sensing

The problem of sparse signal reconstruction has attracted the attention of many researchers in the field of signal processing, machine learning and computer vision. This problem involves solving minimization of an objective function containing quadratic $\ell_2$ error term and a sparse $\ell_1$ regularization term as follows

where $x\in R^n$, $h\in R^m$, $A\in R^{m\times n}$ ($m < < n$) is a linear operator, $\rho \geq 0$, $\|x\|_2$ is the Euclidean norm of $x$ and $\|x\|_{1} = \sum_{i = 1}^n |x_i|$ is the $\ell_1-$norm of $x$.

A lot of methods have been developed for solving (4.1) some of which can be found in ^{[40,41,42,43,44,45]}. Figueiredo et al. ^[43] consider reformulating (4.1) into a quadratic problem by expressing $x\in\mathbb{R}^n$ into two parts as

where $t_i = (x_i)_+$, $y_i = (-x_i)_+$ for all $i = 1, 2, ..., n$, and $(.)_+ = \max\{0, .\}$. Also, we have $\|x\|_1 = e_{n}^T t + e_{n}^T y$, where $e_n = (1, 1, ..., 1)^T\in R^n$. From this reformulation, we can write (4.1) as

from ^[43], Eq (4.2) can be written as

where $z = \begin{pmatrix} t \\ y \\ \end{pmatrix}$, $c = \omega e_{2n}+\begin{pmatrix} -a \\ a \\ \end{pmatrix}$, $a = A^{T}h$, $E = \begin{pmatrix} A^T A & -A^T A \\ -A^T A & A^T A\\ \end{pmatrix}$.

It is not difficult to see that $E$ is a positive semi-definite showing that problem (4.3) is quadratic programming problem.

Xiao et al ^[17] further translated (4.3) into a linear variable problem, equivalently, a linear complementary problem and the variable $z$ solves the linear complementary problem provided that it solves the nonlinear equation:

where $A$ is a vector-valued function. In ^[8,46], it is proved that the function $A(z)$ is continuous and monotone. Thus, problem (4.1) is equivalent to problem (1.1). Therefore, the algorithm we proposed in this work to solve (1.1) can efficiently solve (4.1) as well.

As an application, we consider applying our proposed algorithm in reconstructing a sparse signal of length n from k observations using mean squared error (MSE) as a metric for assessing quality reconstruction. The (MSE) is defined as

where $s$ represents the original signal and $\tilde{s}$ the restored signal. We choose $n = 2^{12}, \; k = 2^{10}$ to be the size of the signal and the original signal contains $2^7$ randomly nonzero elements. The measurement y contains noise, $y = As+\omega$, where $\omega$ is the Gaussian noise distributed as $N(0, 10^{-4})$ and $A$ is the Gaussian matrix generated by command $randn(m, n),$ in Matlab.

We compared the performance of our proposed algorithm (MDY) with SGCS proposed in ^[8]. The parameters in SGCS are maintained as they are in ^[8], while in MDY we choose $r = 0.001, \; \theta = 1/(k+1)^2, \; \mu = 1.1, \; \gamma = 0.1, \; \sigma = 0.01, \; \Lambda = 1,$ and $\beta = 0.65$. Each code is run with same initial point and continuation technique on parameter $\mu.$ We only focused on the convergence behaviour of each method to obtain a solution with similar accuracy. We initialized the experiments by $x_0 = A^Ty$ and terminated the iteration when the relative change in the objective function satisfies

The performance of both MDY and SGCS are shown in Figures 4 and 5. Figure 4 shows that both the MDY and the SGCS methods recover the signal. However, looking at the reported metrics of each method, it can be observed that the MDY method is more efficient since it has a lesser MSE, and its recovery has fewer number of iterations and CPU time. To show the performance of both methods graphically, we plotted four graphs (see Figure 5) demonstrating the convergence behaviour of the MDY method and SGCS method based on the MSE, objective function, number of iterations and CPU time. From Figure 5, it can be observed that the proposed MDY method has faster convergence rate compared to the SGCS method. This shows that the MDY method can be a good alternative solver for signal recovery problems.

5.   Conclusions

In the work, a spectral conjugate gradient algorithm is proposed. The search direction uses a convex combination of the well known DY conjugate gradient parameter and a modified conjugate descent parameter. The search direction is suffiently descent, and global convergence of the proposed algorithm is proved under some assumptions. Numerical experiments are reported to show the efficiency of the algorithm in comparison with the PDY algorithm proposed in ^[22]. In addition, an application of the proposed algorithm is shown in signal recovery and the result is compared with SGCS algorithm proposed in ^[8]. Based on the results obtained, it can be observed that the proposed algorithm has a better performance than the PDY and SGCS {algorithms} in numerical and signal recovery {experiments} respectively. Future work include applying the new proposed algorithm to solve 2D robotic motion control as presented in ^[47].

Acknowledgments

The authors acknowledge the financial support provided by the Center of Excellence in Theoretical and Computational Science (TaCS-CoE), KMUTT and Center under Computational and Applied Science for Smart Innovation research Cluster (CLASSIC), Faculty of Science, KMUTT. {This research was funded by Thailand Science Research and Innovation Fund, and King Mongkut's University of Technology North Bangkok with Contract no. KMUTNB-BasicR-64-22.} The first author was supported by the Petchra Pra Jom Klao Ph.D. research Scholarship from King Mongkut's University of Technology Thonburi (KMUTT) Thailand (Grant No. 19/2562). Also, Aliyu Muhammed Awwal would like to thank the Postdoctoral Fellowship from King Mongkut's University of Technology Thonburi (KMUTT), Thailand. Moreover, this research project is supported by Thailand Science Research and Innovation (TSRI) Basic Research Fund: Fiscal year 2021 under project number 64A306000005.

Conflict of interest

The authors declare that they have no conflict of interest.