Convergence analysis of a three-term extended RMIL CGP-based algorithm for constrained nonlinear equations and image denoising applications

1.   Introduction

Nonlinear equations are integral to various practical applications, such as image segmentation ^[1], financial forecasting ^[2], compressed sensing ^[3,4], machine learning ^[5], neural networks ^[6], Lyapunov matrix equations ^[7], and traffic assignment ^[8]. In this paper, we specifically address the study of such equations, focusing on their formulation and properties. A typical form of these nonlinear equations with convex constraints is defined as follows:

Here, the function $F:\Gamma \to \mathbb{R}^n$ is both monotonic and continuous, and $\Gamma$ represents a closed convex subset of $\mathbb{R}^n$. The monotonicity condition, in particular, implies that for any $a, b\in \Gamma$, the following inequations hold:

There are numerous iterative methods designed to solve nonlinear equations with convex constraints, which can be broadly categorized into two major classes. The first class involves utilizing the Jacobian matrix or its approximations, such as Newton's methods ^[9], quasi-Newton methods ^[10], and Levenberg-Marquardt methods ^[11]. These methods are well-known for their rapid local convergence, with the key advantage being their quadratic convergence rate near the solution. However, a significant drawback arises from the computation and storage requirements of the Jacobian matrix or its approximation. Specifically, the need to explicitly compute and store the Jacobian, or even to approximate it iteratively, can result in high computational costs and memory demands. In contrast, the second class involves utilizing first-order gradient-based information, such as spectral gradient methods and conjugate gradient methods (e.g., ^{[12,13,14,15,16,17]}). These methods are favored for their simpler structure and significantly lower storage requirements, making them well-suited for solving large-scale nonlinear equations. Unlike the first class, first-order gradient-based methods do not require the computation or storage of the full Jacobian, thereby drastically reducing the computational burden. However, a critical limitation of these methods is that they usually require the function to satisfy the Lipschitz continuity condition. This constraint has motivated the development of the new methods that operate under weaker theoretical assumptions.

Recently, some researchers have turned their attention to exploring novel methods that do not require Lipschitz continuity. This is driven by the limitations of traditional optimization algorithms that heavily rely on this condition for global convergence. For example, Jian et al. ^[18] proposed an inertial derivative-free projection method that ensures global convergence without requiring Lipschitz continuity. Their numerical experiments of this method demonstrate that it has promising results on test problems and compressed sensing problems. Similarly, Liu et al. ^[19] introduced a new three-term conjugate gradient projection (CGP) method with an adaptive line search strategy for solving large-scale nonlinear monotone equations with convex constraints, while guaranteeing global convergence without Lipschitz continuity. Liu et al. ^[20] also developed an efficient algorithm for solving convex constrained equations, while emphasizing its global convergence without the need for Lipschitz continuity. Li et al. ^[21] improved the three-term conjugate gradient algorithm for solving nonlinear equations with convex constraints. The authors established global convergence without Lipschitz continuity and demonstrated its effectiveness in benchmark problems and image denoising. In parallel, other researchers have focused on extending classical conjugate gradient methods for unconstrained optimization, aiming to enhance their applicability and performance in a wider range of problems. For example, Ibrahim et al. ^[22] introduced a derivative-free method that extends the three-term Polak-Ribière-Polyak (PRP) and the three-term Hestenes-Stiefel (HS) methods. Zheng et al. ^[23] enhanced the Fletcher-Reeves (FR) conjugate coefficient by incorporating a shrinkage multiplier, which leads to a derivative-free two-term search direction and its extended spectral version. Similarly, Wu et al. ^[24] extended the conjugate coefficient with a derivative-free projection technique to solve nonlinear monotone equations, proving global convergence without Lipschitz continuity.

Motivated by these methods ^[25,26,27] and the projection technique, this paper extends these methods and proposed a three-term extended Rivaie-Mustafa-Ismail-Leong (RMIL) conjugate gradient-based projection algorithm for solving the problem (1.1). In this algorithm, the global convergence is proven without requiring the Lipschitz continuity condition. The structure of this paper is as follows. In Section 2, we provide a detailed formulation of the search direction. Section 3 outlines the illustration of the proposed algorithm and establishes the global convergence of the proposed algorithm. Sections 4 and 5 present the numerical results and performance profiles on benchmark problems and image denoising problems. Finally, the conclusion is given in Section 6. For notational clarity, we use $\|\cdot\|$ to denote the Euclidean vector norm throughout the paper.

2.   The construction of search direction

Before presenting the construction of our search direction, we begin by reviewing several known different conjugate coefficients commonly used in unconstrained optimization $\min \{f(a) \mid a \in \mathbb{R}^n\}$. Among these, one of the most efficient and widely recognized methods is the PRP-type conjugate coefficient, which is defined as:

where $g_k: = \nabla f(a_k)$ represents the gradient, and $y_k = g_k - g_{k-1}$ is the difference between successive gradients. From a computational efficiency standpoint, the PRP method has been shown to outperform many other conjugate gradient methods. In an effort to further enhance the performance of this coefficient, Rivaie et al. ^[25] proposed a modification by replacing the denominator $\|g_{k-1}\|^2$ with $\|d_{k-1}\|^2$, resulting in the RMIL conjugate coefficient:

where $d_{k-1}$ denotes the previous search direction. The resulting coefficient is not only simple in form but also easy to implement in practice. Moreover, theoretical analysis has shown that the conjugate gradient method employing this coefficient enjoys global convergence with a linear rate. While this adjustment aims to improve convergence properties, it raises a concern regarding the nonnegativity of $\beta_k^\text{RMIL}$, which is essential for ensuring global convergence under certain conditions. To mitigate this issue, Dai et al. ^[26] introduced a safeguard mechanism and proposed the RMIL variant:

This modification guarantees nonnegativity. However, when the condition fails, it comes at the cost of potentially degenerating the search direction to the steepest descent. This can significantly slow down convergence in practice. In response to this limitation, Xia et al. ^[27] recently proposed a more robust conjugate coefficient that not only preserves nonnegativity but also avoids degeneration to solve the problem (1.1). For convenient writing, let $F_k$ and $w_{k-1}$ denote $F(a_k)$ and $F_k-F_{k-1}$, respectively. Then, their coefficient is defined as:

where $\mu_1 > 1$ is a positive parameter that enhances numerical stability. Utilizing this idea, we explore the efficient variant of conjugate gradient methods known as the three-term conjugate gradient method. In this approach, the search direction is given by:

where $\beta_k$ and $\theta_k$ are conjugate and scale coefficients, respectively. Based on this structure, we extend the conjugate coefficient $\beta_k$ defined in (2.1) to propose a novel three-term RMIL-type CGP algorithm designed to solve the problem (1.1). To ensure the sufficient descent property of the search direction and the trust-region feature of the search direction, we further define the scale coefficient $\theta_k$ as:

where $\mu_2 > 1$.

The following lemma shows that the proposed search direction possesses the sufficient descent property and trust-region feature without the need for any line search approaches.

Lemma 1. If the search direction sequence $\{d_k\}$ is generated by (2.1), (2.2) and (2.3), then the following conclusions hold:

Proof. For $k = 0$, it is easy to deduce that (2.4) holds. For $k \geq 1$, utilizing the Cauchy-Schwartz inequality, and applying the definitions of (2.1) and (2.3), we obtain the following relationships for $\beta_k$ and $\theta_k$:

Multiplying both sides of (2.2) by $F_k^\text{T}$ and substituting the bounds from (2.6) and (2.7), we obtain:

which proves that the search direction $d_k$ satisfies sufficient descent (2.4).

To establish the trust-region feature, we further apply the Cauchy-Schwartz inequality and (2.4) to derive:

which simplifies to: $\left\|d_k\right\|\geq \left(1-\frac{1}{2\mu_1}-\frac{1}{2\mu_2} \right)\left\|F_k\right\|$. Together with (2.6) and (2.7), we can further refine (2.2) as follows:

Thus, the search direction $d_k$ satisfies the trust-region feature (2.5). Therefore, the proof is complete. □

3.   Algorithmic framework and convergence analysis

In this section, we propose a novel three-term RMIL CGP algorithm (abbreviated as TTRMIL algorithm) for solving constrained nonlinear equations. Before presenting the algorithm, we first define the line search approach and the projection operator, which are essential to ensuring the stability of the algorithm. The line search approach is crucial for maintaining the algorithm's stability, and we employ the following condition to guide the search direction ^[28]:

where $\kappa > 0$, and $t_k : = \max\left\{ \xi\rho^i\; |\; i = 0, 1, 2, \cdots \right\}$ with $\xi > 0$ and $\rho\in (0, 1)$. This ensures that the algorithm maintains stability while avoiding excessively small step sizes that would slow down convergence. To ensure that each iterative point remains feasible in the problem's constraints, we define the projection operator ^[14,16] as follows:

This projection operator possesses a well-known non-expansive property, i.e.,

With these foundational concepts established, we now proceed to present the detailed steps of the TTRMIL algorithm, which is summarized in Algorithm 1.

Remark 1. The parameter $l$ in (3.4) serves as a relaxation factor, which is introduced to potentially enhance the convergence speed of the algorithm. In practical applications, choosing a relatively larger relaxation $l\in [1, 2)$ has been observed to yield better empirical performance in terms of convergence speed.

In the following, we analyze the global convergence of the TTRMIL algorithm. We assume that the sequences $\{a_k\}$, $\{b_k\}$ and $\{d_k\}$ generated by the TTRMIL algorithm are infinite; otherwise, the final iterations will yield the optimal solution to the problem described in (1.1). To facilitate the following analysis, we provide the following general assumptions:

(H1) The solution set $\Gamma_*$ of the problem described in (1.1) is nonempty. (H2) The function $F(\cdot)$ is continuous and monotone on $\mathbb{R}^n$.

Lemma 2. There exists a positive integer $i_k$ such that $t_k = \xi\rho^{i_k}$ satisfies the condition (3.1) for any $k \geq 0$.

Proof. The existence of such an integer $i_k$ follows from a standard backtracking line search approach. In particular, the structure of the condition (3.1) and the properties of the parameter sequence align closely with those analyzed in Lemma 2 of ^[28]. As the proof involves routine analysis and follows similar steps, we omit the detailed derivation here for brevity. □

Lemma 3. Suppose that Assumptions (H1) and (H2) are satisfied. For a solution $a_* \in \Gamma_*$, the sequence $\{\|a_k-a_*\|\}$ converges, and the sequence $\{a_k\}$ is bounded, where $\{a_k\}$ is generated by the TTRMIL algorithm.

Proof. We begin by applying the monotone property from (1.2) with $a = b_k$ and $b = a_*$. This yields the inequality $F(b_k)^\text{T}(b_k-a_*)\geq 0$, which can be rewritten as:

Next, using the projection operator defined in (3.2) and the non-expansive property defined in (3.3), we obtain the following inequality:

Substituting (3.5) into this inequality and applying the definition of $\varpi_k$, we derive:

which yields that $\left\| a_{k+1} - a_* \right\| \leq \left\| a_{k} - a_* \right\|$ with $l \in (0, 2)$. Therefore, the sequence $\{\|a_k-a_*\|\}$ is convergent. Since the norm sequence converges, it follows that the sequence $\{a_k\}$ is bounded. Therefore, the proof is complete. □

Lemma 4. Suppose that Assumptions (H1) and (H2) are satisfied. Then, the following relationship holds:

where the sequences $\{t_k\}$ and $\{d_k\}$ are generated by the TTRMIL algorithm.

Proof. We begin by recalling the definition of $b_k$ and multiplying both sides of (3.1) by $t_k$, which leads to the following inequality:

Substituting this result into (3.6), we get:

This can be rewritten as:

To further summarize this result over all iterations, we sum the inequalities for $k$ from 0 to $\infty$, yielding:

Thus, we conclude that $\lim\limits_{k \to \infty} \|a_k - b_k \| = 0$. From the definition of $b_k$, this implies $\lim\limits_{k \to \infty} t_k \|d_k\| = 0$. Therefore, the proof is complete. □

Theorem 1. Suppose that Assumptions (H1) and (H2) are satisfied. Then, the following conclusion holds:

where the sequence $\{a_k\}$ is generated by the TTRMIL algorithm.

Proof. We proceed by contradiction. Assume that (3.8) is not true, meaning there exists a constant $\delta_0 > 0$ such that $\|F(a_k)\| > \delta_0$ for all $k \geq 0$. From (2.5), we deduce that:

Additionally, since $F(\cdot)$ is continuous and the sequence $\{a_k\}$ is bounded, it follows that the sequence $\{F(a_k)\}$ is also bounded. Therefore, there exists a constant $\delta_1 > 0$ such that $\|F(a_k)\| < \delta_1$ for all $k \geq 0$. From (2.5), we also deduce that:

Together with (3.9) and (3.10), we conclude that the sequence $\{d_k\}$ is bounded. Therefore, we can apply the result from (3.7), which gives us:

Given that $\{a_k\}$ and $\{d_k\}$ are both bounded, there exists an infinite index set $\Omega$ such that

Now, observe that for $t_{k_n}$, the line search approach defined in (3.1) is violated:

Taking the limit as $n \to \infty$, we find:

In addition, together with (2.4), we have:

which leads to:

Thus, (3.11) is contradicts with (3.12). Therefore, the proof is complete. □

4.   Numerical experiments

In this section, we evaluate the numerical performance of the TTRMIL, self-adaptive three-term conjugate gradient method (SATTCGM) ^[29], modified spectral conjugate gradient (MSCG) ^[30], and modified RMIL (MRMIL) ^[27] algorithms in solving large-scale constrained nonlinear equations. All experiments are performed on a Lenovo PC with an Intel (R) Core(TM) i7-12700F processor (2.10 GHz), 16 GB of RAM, and Windows 11 operation system. The parameters for the TTRMIL algorithm are chosen as follows:

For the SATTCGP, MSCG, and MRMIL algorithms, the parameters remain consistent with those presented in their original literatures. The algorithms stopped when either of the following conditions is met: (ⅰ) $\|F_k\| \leq \delta$ or (ⅱ) NI $> 10,000$. Here, NI denotes the number of iterations. The benchmark problems are solved by eight different initial points, as detailed in Table 1, and the problem dimensions are selected from the set: [1000, 5000, 10,000, 50,000,100,000]. The benchmark problems include the following:

Problem 1 ^[30]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 2 ^[30]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 3 ^[30]:

and $\Gamma = [-1, +\infty)$.

Problem 4 ^[28]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 5 ^[28]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 6 ^[30]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 7 ^[28]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 8 ^[30]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 9. ^[31]:

and $\Gamma = \mathbb{R}^n_+, \bar{a} = \sum\limits_{i = 1}^n a_i^2, c = 10^{-5}$.

Problem 10 ^[31]:

and $\Gamma = \mathbb{R}^n_+$.

Problem 11 ^[32]:

and $m = 0.3$, $\Gamma = \mathbb{R}^n_+$.

The experimental results are compared and presented in Tables 2–12, where the following notations are used: "Init" refers to the initial points; "$n$" refers to the problem dimensions, multiplied by 1000; "CPUT" refers to the central processing unit (CPU) time in seconds; "NF" refers to the number of function evaluations; "NI" refers to the number of iterations; "NORM" refers to the norm of $F_k$ at the final iteration. To clearly illustrate the numerical performance, we employ the performance profiles developed by Dolan and Morè ^[33] for comparing the performance in terms of CPUT, NF and NI. The performance profiles used in our evaluation are defined as follows:

In this expression, $\Omega_p$ denotes the test set, and $|\Omega_p|$ denotes the total number of test problems. The set $\Pi$ includes all candidate solvers, while $\omega_{p, q}$ represents the performance indicator of solver $q \in \Pi$ on problem instance $\omega_p \in \Omega_p$. From Tables 2–12 and Figures 1–3, we can draw the following conclusions: (1) As shown in Tables 2–12, all benchmark problems are solved successfully by the TTRMIL, SATTCGM, MSCG, and MRMIL algorithms, demonstrating their effectiveness. (2) Figure 1 illustrates that the TTRMIL algorithm solves 37.73% of the benchmark problems in the least amount of CPU time, outperforming the SATTCGM, MSCG and MRMIL algorithms, which solve 28.64%, 6.14%, and 28.41%, respectively, in the same time. (3) Figure 2 highlights that the TTRMIL algorithm achieves a lower number of function evaluations for 50.23% of the benchmark problems, surpassing the SATTCGM, MSCG, and MRMIL algorithms, which solve 34.09%, 6.14%, and 27.27%, respectively, with fewer evaluations. (4) Figure 3 demonstrates that the TTRMIL algorithm requires fewer iterations to solve 70.68% of the benchmark problems, outperforming the SATTCGM, MSCG, and MRMIL algorithms, which solve 14.55%, 6.36%, and 35.68%, respectively, with fewer iterations. Moreover, to further highlight the robustness of the TTRMIL algorithm, we examine its performance on a particularly challenging class of problems, involving non-Lipschitz functions. In particular, Problem 11 serves as a representative case. As illustrated in Table 12, the TTRMIL algorithm solves 77.5% of the 40 cases with the least CPU time, and achieves the lowest number of function evaluations and iteration counts in 87.5% of the cases for each metric.

These results consistently highlight the efficiency and convergence strength of the TTRMIL algorithm across multiple metrics. Notably, the enhanced performance can be attributed to the computational contribution of the third term in its search direction, which effectively accelerates convergence and reduces overall resource consumption.

5.   Image denoising applications

Image denoising is a prominent inverse problem in compressive sensing, where various sources of noise complicate the process, including defective pixels in camera sensors, hardware storage errors, and interference during transmission through noisy channels. Among these, salt-and-pepper noise poses a particular challenge, as it can obscure vital image details and edges, which are crucial for applications in fields such as medical imaging, remote sensing, and object recognition. To tackle this issue, earlier studies (e.g., ^[34,35]) proposed a robust two-phase approach designed specifically to detect and remove salt-and-pepper noise. In the following, we provide a comprehensive explanation of this approach, which has proven effective in addressing these challenges.

Let an original image $x$ have dimensions $m \times n$, where $x_{i, j}$ denotes the grayscale value at pixel location $(i, j) \in \mathcal{I} = \{1, 2, \ldots, m\} \times \{1, 2, \ldots, n\}$. For the purposes of image processing, we define the neighborhood of each pixel $(i, j)$, denoted by $\mathcal{M}_{i, j} = \{(i, j-1), (i, j+1), (i-1, j), (i+1, j)\}$. When the image $x$ is corrupted by the salt-and-pepper noise, the resulting noisy image is represented by $o$. The two-phase approach begins with the application of an adaptive median filter to the noisy image $o$, resulting in an intermediate image $\tilde{o}$. In the first phase, noise candidates are identified based on the difference between $o$ and $\tilde{o}$. The set of noise candidates is defined as:

The second phase focuses on recovering the noisy pixels in the set $\mathcal{N}$. For each pixel $(i, j) \in \mathcal{N}$, if it is determined to be non-noisy, its value is retained as $u^*_{i, j} = o_{i, j}$. For noisy pixels, their values are recovered by considering the values of their non-noisy neighbors. Specifically, non-noisy neighboring pixels are preserved as $u^*_{m, n} = o_{m, n}$ for $(m, n) \in \mathcal{M}_{i, j} \setminus \mathcal{N}$, while the values of noisy neighbors, $(m, n) \in \mathcal{M}_{i, j} \cap \mathcal{N}$, are also recovered. To achieve the restoration, the following objective function is minimized:

Here, $\phi_{\alpha}$ is an edge-preserving potential function with a coefficient $\alpha$. According to Corolly 1 in ^[35], if $\phi_{\alpha}$ is convex and first order Lipschitz continuous, then $\nabla f(u)$ is also Lipschitz continuous, where $\nabla f(u)$ is the gradient of the function $f(u)$. Furthermore, Proposition 6 in the same work states that if $\phi_{\alpha}$ is convex, then $\nabla f(u)$ is monotone. To satisfy these conditions, we adopt the well-known Huber function as the form of $\phi_{\alpha}$, defined as:

It is straightforward to verify that the Huber function is both convex and Lipschitz continuous. As a result, the proposed algorithm is well-suited to solve the underlying optimization problem with both theoretical soundness and numerical stability.

We utilize several well-known test images: Boat (512 $\times$ 512), Clown (512 $\times$ 512), Couple (512 $\times$ 512), Crowd (512 $\times$ 512), Tank (512 $\times$ 512), Trucks (512 $\times$ 512), Woman (512 $\times$ 512), Zelda (512 $\times$ 512), which can be found in https://hlevkin.com/. Each of these images is corrupted by 30% salt-and-pepper noise. The experiments are repeated 10 times with different noise samples to ensure robustness. The noisy images, along with their recovered versions using TTRMIL, SATTCGM, and MSCG algorithms, are presented in Figures 4 and 5. From these figures, it is evident that all images affected by 30% salt-and-pepper are successfully recovered by TTRMIL, SATTCGM, and MSCG algorithms. In addition, the corresponding numerical results, shown in Table 13, include the average number of iterations ($\overline{NI}$), the average CPU time in seconds ($\overline{CPUT}$), the average peak signal-to-noise-ratio ($\overline{PSNR}$), and the average structural similarly index ($\overline{SSIM}$), respectively. The table highlights that, despite achieving similar average structural similarity indices, the TTRMIL algorithm outperforms the SATTCGM and MSCG algorithms by requiring fewer CPU time, and fewer iterations. These finding suggest that the TTRMIL algorithm is not only efficient but also highly competitive in the domain of image denoising.

6.   Conclusions

In this paper, we propose a three-term extended RMIL CGP-based algorithm designed to solve nonlinear equations with convex constraints. The key contribution of this work lies in the extension of the traditional two-term search direction to a more robust three-term direction. By carefully designing the scale coefficient, we ensure both the sufficient descent property and the trust-region feature, thereby obviating the need for computationally expensive line search techniques. Importantly, we establish the global convergence of the proposed algorithm under weaker assumptions without relying on the Lipschitz continuity condition commonly required in optimization literature. To validate the practical effectiveness of the proposed algorithm, we conducted comprehensive numerical experiments on eleven benchmark problems, including Lipschitz and non-Lipschitz functions. The results reveal that the proposed algorithm not only demonstrates superior performance but also consistently outperforms other comparable three-term methods in terms of computational efficiency. Notably, its robustness in handling large-scale problems and non-Lipschitz functions highlights its practical relevance and competitiveness in real-world applications. In the future work, we will integrate the inertial-relax scheme into the proposed algorithm for accelerating the convergent speed, and apply the proposed algorithm to specific application domains, such as machine learning, signal processing, or control systems.

Use of 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 Natural Science Foundation in Guangxi Province, PR China (grant number 2024GXNSFAA010478), the Special projects in key areas of ordinary universities in Guangdong Province, PR China (grant number 2023ZDZX4069), the Research Team Project of Guangzhou Huashang College (grant number 2021HSKT01), and the Guangzhou Huashang College Daoshi Project (grant number 2024HSDS15, 2024HSDS28)

Conflict of interest

The authors declare there is no conflicts of interest.