Global convergence in a hybrid conjugate gradient projection method for finding solutions of constrained nonlinear equations with applications

TL;DR

This paper introduces a new method for solving constrained nonlinear equations that uses less memory and performs better than existing methods on benchmark problems.

Contribution

A hybrid conjugate gradient projection method with global convergence and improved performance on constrained nonlinear equations.

Findings

The method achieves global convergence under reasonable assumptions.

It outperforms existing methods on 75.71% to 86.43% of benchmark problems in efficiency metrics.

Successfully applied to sparse signal restoration problems.

Abstract

In this paper, a hybrid conjugate gradient projection method for finding solutions of constrained nonlinear equations is proposed by integrating both hyperplane projection and hybrid techniques. The key features of this method are as follows: (1) It is characterized by a low storage requirement and relies solely on function values; (2) The designed search direction ensures the sufficient descent property without the need for line search approaches; (3) Under certain reasonable assumptions, the global convergence of the method is established; (4) Experimental results demonstrate that the proposed method outperforms the two existing methods about 75.71%, 85.36%, and 86.43% of benchmark problems in terms of CPU time, the number of function evaluations, and iterations. Furthermore, it is applied to successfully solve the sparse signal restoration problems.

Figures

Performance profiles on CPUT.

Performance profiles on NFE.

Performance profiles on NI.

From top to bottom: the original signal, the observed signal, and the recovered signals by HCGP, HSDY, and DCG methods.

The comparison results of HCGP, HSDY, and DCG methods.The x-axis represents the number of iterations (top right and bottom right) and the CPU time in seconds (top left and bottom left). The y-axis represents the MSE (top left and top right) and the objection value (bottom left and bottom right

Equations

Full text

Introduction

Numerical applications for nonlinear equations span various domains, including financial forecasting [1], compressed sensing problems [2], subproblems in generalized proximal methods with Bergman distances [3], machine learning [4], and robot arm control [5,6]. Motivated by these applications, the goal of this paper is to finding solutions of constrained nonlinear equations:

$$[eqn]$$

where $[eqn]$ is a closed convex and non-empty set, and $[eqn]$ is a continuous and monotone function. Throughout the paper, we denote $[eqn]$ as $[eqn]$ , and $[eqn]$ as the Euclidean norm on $[eqn]$ . To establish a foundation for the discussion, we now review some key concepts related to nonlinear equations:

(1) A function $[eqn]$ is said to be monotone if:

$$[eqn]$$

(2) A function $[eqn]$ satisfies the Lipschitz continuity condition with L > 0 if:

$$[eqn]$$

Since conjugate gradient methods [7,8] are known for their simplicity and low storage requirements, they are widely used to solve the problem (1). The sequence of iterates {xk} generated by these methods is updated as:

$$[eqn]$$

where $[eqn]$ is the step-size determined by a line search approach, and dk is the search direction. Various methods can be employed to generate the search direction dk, such as Newton methods [9–11], quasi-Newton methods [12,13], Levenberg-Marquardt methods [14,15], and their modifications. For example, in [11], Waziri et al. proposed a method that approximates the Jacobian matrix as a diagonal matrix, significantly reducing storage requirements, computational cost, and CPU time. In [12], Ai-Baali et al. explored the properties of quasi-Newton methods, originally introduced by Broyden as an alternative to Newton’s method and later extended to nonlinear equations. In [14], Ou and Lin presented a hybrid method for solving non-smooth equations with box constraints, combining the Levenberg-Marquardt-like techniques and non-monotone strategies. Their method demonstrates both global and local superlinear convergence with promising numerical results. Despite the advantages of these methods, they generally require solving linear equations using either the Jacobian matrix or its approximation at each iteration, which makes them less suitable for large-scale nonlinear problems. Consequently, many researchers have turned to conjugate gradient methods for generating the search direction, as these methods do not rely on such costly operations and can efficiently handle large-scale equations.

Recently, the hybridization of conjugate gradient methods [16,17] has emerged as an interesting and effective technique for designing search directions. Several researches have explored different combinations of conjugate parameters to improve performance across various applications. For example, in [18], Waziri et al. applied the convex combination technique to develop an improved search direction based on Fletcher-Reeves (FR) and Polak-Ribière-Polyak (PRP) conjugate parameters. Similarly, in [19], Sun et al. introduced a hybrid conjugate gradient-based projection method that combines Hestenes-Stiefel (HS) and Dai-Yuan (DY) conjugate parameters. In another approach, Yusuf et al. [20] proposed a hybrid method for solving large-scale nonlinear equations, which merges FR and PRP parameters. Additionally, in [21], Danmalam et al. developed a hybrid conjugate residual method for nonlinear monotone equations with convex constraints, combining FR and a new residual parameter to ensure descent direction without line search approaches. In [22], Koorapetse et al. proposed a hybrid conjugate gradient projection method for large-scale nonlinear monotone equations, ensuring descent property and global convergence. Furthermore, in [23], Li et al. presented a hybrid conjugate gradient-based projection method with modified PRP and DY parameters, and proved global convergence and superior efficiency in numerical experiments, particularly in sparse signal and image restoration. In [24], Li et al. also developed a derivative-free descent method, which incorporates spectral gradient methods, projection techniques, and monotone line searches. In a different approach, Abubakar et al. [25] proposed a hybrid three-term conjugate gradient method that ensures sufficient descent and trust region properties without line search approaches. Lastly, in [26], Yin et al. introduced a hybrid three-term conjugate gradient projection method with adaptive line search for large-scale nonlinear monotone equations with convex constraints. These diverse hybrid methods highlight the growing potential and versatility of conjugate gradient methods in optimization, particularly for large-scale and constrained problems.

Motivated by the works in [27–29], we adopt the hybridization techniques to modify a novel conjugate parameter, which is designed to ensure the global convergence of the method. Based on this parameter, we then propose a search direction that satisfies the sufficient descent property. By incorporating effective line search approaches and projection techniques, we propose a hybrid conjugate gradient projection method specifically tailored for finding solutions of constrained nonlinear equations.

Algorithm

To establish a coherent framework for formulating the search direction in constrained optimization, we first review the foundational works in unconstrained optimization, specifically addressing the problem of minimizing a function f(x), where $[eqn]$ . The conjugate gradient method, a widely-used iterative approach in unconstrained optimization, depends heavily on the choice of conjugate parameters that govern the search direction. Hager et al. [27] developed a novel conjugate parameter, which is defined as:

$$[eqn]$$

where $[eqn]$ is the gradient of the function at the current point xk, $[eqn]$ , and $[eqn]$ . Similarly, Yu et al. [28] developed a modified version of the PRP conjugate parameter:

$$[eqn]$$

where the first term promotes the classical PRP conjugacy, while the second term introduces a correction to account for the history of the gradients, providing a more adaptive search direction. Moreover, Li et al. [29] presented a modified Liu-Story (LS) conjugate parameter, which is expressed as:

$$[eqn]$$

Notably, all of these conjugate parameters share a common structural pattern, which suggests that their forms can be unified into a more general framework for optimization. Motivated by these findings, we extend these conjugate parameters to address constrained optimization problems, particularly those involving nonlinear equations as in problem (1). This insight leads us to propose a hybrid conjugate parameter to ensure the global convergence of the method:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . To ensure that the search direction satisfies the sufficient descent property, we define the following search direction:

$$[eqn]$$

In this section, we propose a hybrid conjugate gradient projection (HCGP) method to find solutions of constrained nonlinear equations. To ensure clarity and precision, we first discuss the line search strategy employed in our method, as well as the projection operator used to enforce the feasibility of iterates.

We adopt a line search strategy based on the approach outline in the work [30]. This approach determines the optimal step-size $[eqn]$ through the following inequality condition:

$$[eqn]$$

where $[eqn]$ is a small constant, and the step-size is selected according to the rule $[eqn]$ with parameters $[eqn]$ and $[eqn]$ . This approach ensures a gradual and controlled descent along the search direction while maintaining stability in the optimization process.

A critical aspect of constrained optimization is ensuring that each iterative point lies within the feasible region E. To address this, we define the projection operator PE that projects any point $[eqn]$ onto the feasible set. This operator minimizes the Euclidean distance between x and the feasible region E:

$$[eqn]$$

The projection operator exhibits a well-known non-expansive property, which ensures that the distance between projections of two points is no greater than the distance between the original points themselves:

$$[eqn]$$

With the line search strategy and projection operator, our method computes the trial point $[eqn]$ at each iteration. If the trial point wk does not lie within the feasible region E, we update the next iterative point xk + 1 using the following update rule:

$$[eqn]$$

where $[eqn]$ is a scaling factor.

The iterative process continues until the term $[eqn]$ falls below a given threshold δ, indicating convergence to a solution. The detailed steps of the method are summarized in Algorithm 1.

Algorithm 1. HCGP algorithm.

1: Initialization: Parameters: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ ,

$[eqn]$ , $[eqn]$ . Set $[eqn]$ .

2: while $[eqn]$ do

3:   Compute the search direction dk by (5).

4:   Determine the step-size $[eqn]$ by (6) and compute the trial

point $[eqn]$ .

5:   if $[eqn]$ and $[eqn]$ then

6:    Break.

7:   else

8:    Update the next iteration xk + 1 by (8).

9:   end if

10:   Set $[eqn]$ .

11: end while

Analysis of convergence properties

In order to demonstrate that the sequence {xk} generated by the HCGP method converges to its optimal solution of the problem (1), we must establish the following general assumptions:

(A1) The solution set $[eqn]$ of the problem (1) is non-empty;

(A2) The function e(x) is both monotone and Lipschitz continuous on $[eqn]$ .

Search direction analysis

With these assumptions in place, we proceed to analyze the properties of the search direction sequence {dk}, which is a key element in proving convergence. The following lemma characterizes the sufficient descent property of the search direction sequence.

Lemma 1. The search direction sequence {dk} satisfies the sufficient descent property for any $[eqn]$ :

$$[eqn]$$

Moreover, for $[eqn]$ , we have:

$$[eqn]$$

Proof: When k = 0, we have $[eqn]$ . When $[eqn]$ , multiplying both sides of (5) by ek, we have:

$$[eqn]$$

Thus, the sufficient descent property (9) holds. Additionally, by applying the Cauchy-Schwartz inequality and using the sufficient descent property, we have:

$$[eqn]$$

which implies that $[eqn]$ . □

Line search analysis

Lemma 2. Suppose that Assumptions (A1) and (A2) hold. Then the following conclusion are satisfied:

(1) The line search strategy (6) is well-defined.

(2) The step-size $[eqn]$ satisfies the following relation:

$$[eqn]$$

Proof: To prove the lemma, we address each part of the statement separately.

Part 1: Suppose, by contradiction, that there exists a positive integer $[eqn]$ such that the line search strategy (6) dose not hold for any integer i > 0. This means that for all i > 0, we obtain:

$$[eqn]$$

Given the continuity of $[eqn]$ , we can take the limit as $[eqn]$ , yielding:

$$[eqn]$$

Besides, from the sufficient descent property (6), we know that:

$$[eqn]$$

which contradicts (12). Thus, such a k0 cannot exist, and we conclude that the line search strategy (6) is well-defined.

Part 2: Assume that $[eqn]$ . In this case, we define $[eqn]$ , which does not satisfy the line search strategy (6), implying that:

$$[eqn]$$

Using the sufficient descent property (13) and combining with (9), we obtain the following inequality:

$$[eqn]$$

which is rewritten as:

$$[eqn]$$

□

Global convergence analysis

Lemma 3. Suppose that Assumptions (A1) and (A2) hold. Then, the sequence $[eqn]$ is convergent for any $[eqn]$ . Furthermore, the sequences {xk} and {wk} are bounded, and

$$[eqn]$$

Proof: The result follows directly from a similar result presented in Lemma 4 of [23]. Therefore, the proof is omitted here for brevity.

Lemma 4. Suppose that Assumptions (A1) and (A2) hold. If $[eqn]$ , there exists a positive constant N such that, for all $[eqn]$ ,

$$[eqn]$$

Proof: From (3), (8), we have the following relationship:

$$[eqn]$$

Utilizing this in conjunction with (4), we obtain:

$$[eqn]$$

Additionally, using Lemma 3, we know that the sequence {ek} is bounded due to the continuity of $[eqn]$ and the boundedness of {xk}. That is, there exists a constant r>0 such that $[eqn]$ holds. From this result, (5), and (16), it follows that

$$[eqn]$$

Thus, we have $[eqn]$ with $[eqn]$ . □

Theorem 1. Suppose that Assumptions (A1) and (A2) hold. Then, we have:

$$[eqn]$$

Proof: We proceed by assuming the contrary, i.e., (17) does not hold. This implies the existence of a constant $[eqn]$ such that $[eqn]$ holds for all k. Given this assumption, and using the result (10), we immediately deduce the following inequality:

$$[eqn]$$

Next, by multiplying both sides of (11) by $[eqn]$ , we obtain the following expression:

$$[eqn]$$

Substituting (15) and (18) into (19), we obtain:

$$[eqn]$$

This leads to a contradiction with (14), which shows that our assumption that $[eqn]$ for all k must be false. □

Numerical experiments

In this section, we compare the performance of the HCGP method with that of two similar methods: HSDY method [31] and DCG method [32]. All experiments are conducted on a Lenovo PC equipped with an Intel(R) Core(TM) i7-12700F processor (2.10GHz), 16GB of RAM, and Windows 11 operation system. For the HCGP method, the following parameters are used:

$$[eqn]$$

The parameters for the HSDY and DCG methods are chosen according to the settings specified in the original literature.

In this experiment, we compare the performance of the HCGP, HSDY, and DCG methods on eight benchmark problems. For all methods, the programs are terminated when one of the following criteria is met:

$$[eqn]$$

where “NI” is the number of iterations. The dimensions of the benchmark problems is chosen from the set: [1000, 5000, 10000, 50000, 100000]. The initial points for the benchmark problems are selected as follows: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ . The specific formulations of the benchmark problems are as follows:

Problem 1:

$$[eqn]$$ $$[eqn]$$

and $[eqn]$ .

Problem 2:

$$[eqn]$$

and $[eqn]$ . Clearly, this problem has a unique solution $[eqn]$ .

Problem 3:

$$[eqn]$$

and $[eqn]$ .

Problem 4;

$$[eqn]$$

and $[eqn]$ .

Problem 5:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

and $[eqn]$ .

Problem 6:

$$[eqn]$$

and $[eqn]$ .

Problem 7:

$$[eqn]$$

and $[eqn]$ .

Problem 8:

$$[eqn]$$

and $[eqn]$ .

The experiment results are summarized in [Tables 1–8](#pone.0335265.t001 pone.0335265.t002 pone.0335265.t003 pone.0335265.t004 pone.0335265.t005 pone.0335265.t006 pone.0335265.t007 pone.0335265.t008), where “CPUT” represents the CPU time in second, “NFE” represents the number of function evaluations, “NI” represents the number of iterations, “NORM” represents the norm value of the function at the approximative point, and “Init(n)” represents the initial point and the dimension (multiplied by 1000). To provide a visual representation of the computational differences, we utilize the performance profiles introduced by Dolan and Moré [33]. These profiles allow us to assess the relative performance of the methods, where higher curves indicate better numerical efficiency. In Fig 1, the HCGP method demonstrates the lowest CPU time in nearly 75.71% of the benchmark problems, outperforming the HSDY and DCG methods, which achieve this in approximately 8.21%, 16.79% of cases, respectively. Similarly, in Fig 2, the HCGP method consistently shows the fewest number of function evaluations, excelling in almost 85.36% of the benchmark problems. In contrast, the HSDY and DCG methods reach this efficiency in 18.21% and 13.21% of cases, respectively. Lastly, Fig 3 illustrates that the HCGP method also requires the fewest number of iterations in 86.43% of the benchmark problems, while the HSDY and DCG methods perform better in only in 17.86% and 23.57% of cases, respectively.

Performance profiles on CPUT.

Performance profiles on NFE.

Performance profiles on NI.

Table 1: Numerical results for Problem 1.

Table 2: Numerical results for Problem 2.

Table 3: Numerical results for Problem 3.

Table 4: Numerical results for Problem 4.

Table 5: Numerical results for Problem 5.

Table 6: Numerical results for Problem 6.

Table 7: Numerical results for Problem 7.

Table 8: Numerical results for Problem 8.

Application on sparse signal restoration

In statistical inference and signal processing, solving linear systems of equations that yield sparse solutions poses a significant challenge. These systems are frequently ill-conditioned, which implies that they are highly sensitive to small perturbations in the input data. This characteristic can lead to instability in the solution, which makes it challenging to compute the solution accurately. To address this, a common approach is to formulate an objective function that combines a sparse $[eqn]$ -norm regularization term and a quadratic $[eqn]$ -norm error term. This approach ensures that the solution remains sparse while minimizing the reconstruction error. The corresponding optimization problem can be expressed as follows:

$$[eqn]$$

where $[eqn]$ represents the observed signal, $[eqn]$ is the unknown signal to be estimated, and $[eqn]$ is a regularization parameter that controls the balance between the error and the sparsity of the solution. The matrix $[eqn]$ defines a linear transformation.

To simplify the optimization, we introduce the decomposition:

$$[eqn]$$

Here, $[eqn]$ and $[eqn]$ for all $[eqn]$ , where $[eqn]$ . By using this decomposition, the $[eqn]$ -norm of v can be rewritten as:

$$[eqn]$$

where $[eqn]$ is a vector of ones. This leads to a reformulation of the original problem (20) as:

$$[eqn]$$

which is a bound-constrained quadratic program. The next step is to express this problem in a form that is more conducive to numerical solution. By defining $[eqn]$ , we can rewrite the problem (21) as:

$$[eqn]$$

where

$$[eqn]$$

Since H is positive semi-definite, the problem (22) represents a convex quadratic program, which guarantees the existence of a unique optimal solution. Xiao et al. [34] showed that this quadratic program can be transformed into the following simpler form:

$$[eqn]$$

They further demonstrated that $[eqn]$ is a continuous and monotone function, which provides useful properties for solving the optimization problem. Consequently, solving the original problem (20) is equivalent to solving the simpler problem (23), which can be expressed in the form of the problem (1). Thus, the proposed method efficiently solves the sparse signal restoration problem.

In this experiment, we compare the performance of the HCGP method against the HSDY and DCG methods. The goal is to restore a sparse signal of length n from an observed signal of length m. The mean squared error (MSE), used to evaluate the quality of the restoration, is defined as:

$$[eqn]$$

where $[eqn]$ and v represent the original signal and the recovered signal, respectively. The dimensions of the signal is set to n = 6144, m = 1536, and k = 192 non-zero elements are randomly selected. The matrix A is randomly generated in MATLAB, and $[eqn]$ is the observed data where $[eqn]$ is the Gaussian noise. The value η is obtained by the same continuation technique for the tested methods, i.e., $[eqn]$ . The initial point is selected as $[eqn]$ and the numerical experiments terminate when the relative change between successive iterative falls below 10^−5^. The stopping criterion for all methods is defined as:

$$[eqn]$$

where fk denotes the function value f(v) at $[eqn]$ .

The sparse signal restored by all methods is shown in Fig 4, where it is clear that the HCGP, HSDY, and DCG methods are all capable of successfully restoring the original signal. In addition, Fig 5 illustrates the comparative performance of these methods in terms of convergence properties, including MSE, the CPU time (CPUT), the number of iterations (NI), and the values of the objective functions (ObjFun). To further assess their effectiveness, the experiment (as shown in Table 9) is repeated multiple times, allowing for a robust comparison of these methods. Specifically, in terms of NI and CPUT, the HCGP method outperforms the other two methods, requiring fewer iterations and less computational time to recover the signal.

From top to bottom: the original signal, the observed signal, and the recovered signals by HCGP, HSDY, and DCG methods.

The comparison results of HCGP, HSDY, and DCG methods.The x-axis represents the number of iterations (top right and bottom right) and the CPU time in seconds (top left and bottom left). The y-axis represents the MSE (top left and top right) and the objection value (bottom left and bottom right

Table 9: Ten experimental results of sparse signal restoration for HCGP, HSDY, and DCG methods.

Conclusions

In this paper, we propose an extension of a hybrid conjugate gradient projection method for finding solutions of constrained nonlinear equations, leveraging hyperplane projection and hybrid techniques. The proposed method is distinguished by its low storage requirement and its exclusive reliance on function values, making it both efficient and practical for large-scale problems. The search direction designed in this method ensures the sufficient descent property, effectively eliminating the need for line search methods. Under suitable assumptions, we prove the global convergence of the proposed method. The experimental results, including those for large-scale constrained nonlinear equations and sparse signal restoration, demonstrate that the proposed method is numerically competitive and efficient, which outperforms similar methods from the literature in terms of both convergence rate and computational efficiency.