A new three-term conjugate gradient algorithm with modified gradient-differences for solving unconstrained optimization problems

1.   Introduction

Since unconstrained optimization problems often arise from scientific computing and mining of big data ^[1,2,3,4], it is valuable to develop efficient numerical algorithms to solve these problems. However, it seems that there is no any algorithm available in the literature which is in commanding position when it is used to solve all the unconstrained optimization problems, compared with other similar algorithms ^{[5,6,7,8,9,10]}. For this reason, many researchers have been studying new numerical methods to solve the unconstrained optimization problems ^[1,11].

Mathematically, a unconstrained optimization problem is written as

where $f: R^{n}\to R$ is continuously differentiable such that its gradient function $g$: $R^{n}\to R^{n}$ is available. By $g_{k}$ we denote the gradient vector of $g$ at $x_{k}$.

Owing to smaller capacity of computation and storage, conjugate gradient methods (CG) are usually used to solve problem (1.1). By CG, the iterative format to generate a sequence of approximate optimal solutions is

where $x_{0}$ is an arbitrarily chosen initial solution, $d_{k}$ is a search direction to efficiently seek for an optimal solution of problem (1.1), and $\alpha_{k} > 0$ is a step size found by line search along $d_{k}$. In general, the search directions in the classical CG methods are given by

where $\beta_{k}$ is the so called conjugate parameter, often being computed by the following classical methods ^[2,12]:

In (1.4), $y_{k-1} = g_{k}-g_{k-1}$ when $k\ge 1$. When $\alpha_k$ in (1.2) is the exact step size and problem (1.1) is a strictly convex quadratic minimization problem, the values of all $\beta_{k}$ in (1.4) are the same. However, for a generic nonlinear objective function, it is often difficult to find the exact step size. Thus, an inexact line search with lower computational cost is generally adopted. For instance, using the strong Wolfe inexact line search, Riahi and Qattan^[13] established global convergence theory of the Fletcher-Reeves CG method and proved its property of local linear convergence. Unfortunately, in most cases, when the Armijo-type line search is used to find the step size $\alpha_k$, it is often difficult to establish global convergence of the classical CG methods, where the search direction $d_k$ is not necessarily descent. For this reason, many variants of CG methods have been proposed to overcome the above difficulty. For instance, using a modified Armijo-type line search, an improved spectral conjugate algorithm was developed in ^[6] and its global convergence was proved. Numerical tests also showed the advantages of this algorithm.

As remarkable extensions of the classical CG methods, three-term CG methods have been attracting extensive research interest ^{[8,9,10,14,15,16]}. The first three-term CG method was proposed in ^[14], which chooses the search directions by

where $\beta_{k} = \beta_k^{HS}$ (or $\beta_k^{FR}$, $\beta_k^{DY}$ etc.), $d_{t}(t \leq k-1)$ was a restart direction, and

By numerical tests, it was shown ^[14] that in the third term of (1.5), the automatical restarts of using the gradient information in (1.6) may improve convergence of the algorithm.

Nazareth ^[15] presented another method of choosing the directions, given by

where $d_{-1} = d_{0} = 0$. It was proved that without requirement of the exact line search, the developed algorithm based on (1.7) can maintain finite termination as applied to solve convex quadratic minimization problems.

In ^[10,17], two three-term conjugate gradient methods were given by

and

respectively. Independent of any line search, it was proved that the directions in (1.8) and (1.9) are sufficiently descent. Since (1.8) and (1.9) can reduce to the standard PRP and HS conjugate gradient methods in (1.4) under the exact line search, respectively, they are regarded as two modified versions of the standard CG methods. It is noteworthy that the search directions in the standard PRP and HS conjugate gradient methods are not necessarily descent in general.

In ^[8,9,16], Andrei suggested three descent three-term CG methods, which computed the search directions by the following different formats:

and

All of them satisfy the conjugacy condition, and except for (1.10), the search directions in (1.11) and (1.12) are descent when the Wolfe line search are used. Numerical experiments indicated that the CG method in ^[8] outperforms the other six algorithms available in the literature.

Recently, Liu et al. ^[18] constructed two three-term CG methods, specified by:

and

respectively. A remarkable property of these search directions is that they were proved to be sufficiently descent under any line search. However, it is unclear whether these directions in (1.13) and (1.14) satisfy any conjugacy condition or not.

Motivated by a need to further improve numerical efficiency of algorithms, we intend to develop a novel three-term CG algorithm such that the search directions may simultaneously possess descent and conjugate properties. Then, using the standard Wolfe inexact line search, we attempt to prove its global and local convergence under appropriate assumptions, and test its numerical performance by solving benchmark test problems.

The remainder of this article is organized as follows. The new three-term CG algorithm is first developed in Section 2. In Section 3, global convergence of this algorithm is proved. Section 4 is devoted to testing of its numerical performance. Conclusions are drawn in the last section.

2.   Development of a new algorithm

In this section, we state ideas to develop a new algorithm, and then present its framework of computer procedures.

Combining the ideas in ^[6,18], we construct the search direction by

where $w_k = \max\{ \left|s_k^T\bar{y}_{k}\right|, s_{k}^{T}y_k\}$, $d_{0} = -g_{0}$, and

is defined as done in ^[6]. Clearly, $\bar{y}_{k}$ can be regarded as a modified difference of gradients. For this reason, we call the proposed CG method in this paper, where the search directions are defined by (2.1), a three-term CG method with modified gradient-differences. In essence, this three-term CG method is an extension of the two-term spectral conjugate gradient method in ^[6].

We first prove the following property of the search directions in (2.1).

Proposition 1. Let $d_{k}$ be given by (2.1). Then, for any $k\geq0$, the following inequality holds:

Proof. By definition, when $k = 0$, we have $g_0^Td_0 = -\|g_0\|^2$. When $k > 0$, we have $g_{k}^{T}d_{k} = -\|g_{k}\|^{2}$ if $w_{k-1} = 0$; Otherwise, it is true that

Consequently, for any $k \geq0$,

i.e., $d_{k}$ is always sufficiently descent.

Remark 1. As pointed out in ^{[19,20,21,22]}, such a sufficiently descent condition like (2.4) plays a critical role in proving global convergence of CG methods.

Based on the above nice property of the search directions (2.1) in Propositions 1, we come to state a framework of computer procedures for solving unconstrained optimization problems (1.1).

Algorithm 1. (New three-term conjugate gradient algorithm (NTTCG))

Step 0. Take an initial (approximate) solution $x_{0}\in R^{n}$ and an initial search direction $d_0 = -g_0$. Choose the parameters $0 < \rho < \sigma < 1$ used in the line search. The tolerance error is $\varepsilon \in (0, 1)$. Set $k: = 0$.

Step 1. If $\|g_k\|_\infty \leq \varepsilon$, then the algorithm stops.

Step 2. Determine the step size $\alpha_k$ by the following standard Wolfe line search:

Step 3. Update the solution by $x_{k+1}: = x_{k}+\alpha_kd_k$. Compute $g_{k+1}$, and compute a search direction $d_{k+1}$ given in (2.1).

Step 4. Set $k: = k+1$. Return to Step 1.

Remark 2. By Proposition 1, we know that $d_{k+1}$ in Step 3 of Algorithm 1 is a sufficiently descent direction at $x_{k+1}$, which ensures that the conducted line search in Step 2 of Algorithm 1 stops in finitely many steps ^[23].

Remark 3. It follows from (2.4) that the inequalities:

hold for any $k\ge 0$. Thus, $\|d_k\| \geq \|g_k\|$ is true for any $k\ge 0$.

3.   Conjugacy properties and convergence analysis

In this section, we study conjugacy property of the search directions defined by (2.1), and establish global and local convergence theory of Algorithm 1.

3.1.   Convex cases

We first study the conjugacy property of the search directions generated by Algorithm 1 in the case that the objective function in problem (1.1) is convex quadratic.

Specifically, a problem of convex quadratic minimization is written as:

where $Q \in R^{n\times n}$ is a given positive definite matrix and $q \in R^n$ is a given vector. When Algorithm 1 is applied to solve problem (3.1), we can prove that the search directions in (2.1) have the following property.

Proposition 2. For problem (3.1), let $d_k$ be chosen by (2.1). Then, by the exact line search, $d_{k+1}$ and $d_{k}$ are conjugate with respect to $Q$ for any $k\geq0$.

Proof. With the exact line search, we have

and

Consequently,

Thus, for $k = 0$,

For $k > 0$, we have

In other words, $d_{k+1}$ and $d_{k}$ are conjugate with respect to $Q$ for any $k\geq0$.

Remark 4. The basic idea to derive the search directions (2.1) is to guarantee their sufficiently descent property by appropriately modifying the steepest descent direction $-g_{k+1}$ (see Proposition 1). It is well known that the conjugate directions in the classic conjugate direction method are not necessarily descent. Proposition 2 demonstrates that our search directions also satisfy the so-called conjugacy condition in the case that the objective function is convex quadratic, although it is not true that any two directions (for example, the two directions $d_{k+2}$ and $d_k$) are conjugate with respect to the matrix $Q$, as in the classic conjugate direction method. In one word, compared with majority of the existing three-term conjugate gradient methods in the literature, the search directions (2.1) have an advantage of simultaneously possessing descent and conjugate properties.

The following result further states global convergence of Algorithm 1 when it is implemented to solve a uniformly convex optimization problem, an extention of the convex quadratic minimization problem (3.1).

Theorem 1. Let $f: R^{n}\to R$ be twice continuously differentiable and uniformly convex on a level set $\Omega = \{x\in R^n | f(x)\leq f(x_0)\}$, i.e., there exists a positive constant $\mu$ such that for all $x, y \in \Omega$,

holds. Let $\{g_k\}$ be the gradient sequence generated by Algorithm 1. Then,

Proof. For the sake of contradiction, we suppose that there exists a constant $\varepsilon > 0$ such that $\|g_{k}\|\geq \varepsilon$ for all $k\in N$.

Since the step size satisfies (2.5), the sequence $\{f(x_k)\}$ generated by Algorithm 1 is decreasing, and all the iterative points $x_k$ are in the level set $\Omega$. Since $f: R^{n}\to R$ is twice continuously differentiable and uniformly convex on $\Omega$, it follows from Steps 2 and 3 in Algorithm 1 that the level set $\Omega$ is a bounded closed convex set, i.e., there exists a positive constant $B > 0$ such that

In addition, the gradient of $f$ is also Lipschitz continuous on $\Omega$, i.e., there exists a constant $L > 0$ such that for all $x, y \in \Omega$, the following inequality holds:

Furthermore, we can prove boundedness of the sequence $\{d_k\}$. Actually, from (3.9) and the uniformly convex property of $f$, we have $\|y_k\| \leq L\|s_k\|$ and $s_k^Ty_k \geq \mu \|s_k\|^2$. Thus, by the Cauchy-Schwarz inequality, we have

Take $M = 1+\frac{2L+1}{\mu }$. From (3.10) and $\|g_k\| \geq \varepsilon$, it follows that

Using (3.9) and the second inequality in (2.5), it yields

Consequently,

hence,

From (3.14) and the boundedness of $f$ on $\Omega$, we know that

which contradicts (3.11). The proof is completed.

We can also prove that Algorithm 1 is R-linearly convergent when the objective function is uniformly convex.

Theorem 2. Suppose that $f: R^{n}\to R$ is twice continuously differentiable and uniformly convex on the level set $\Omega$, and that the sequence $\{x_k\}$ generated by Algorithm 1 converges to the unique optimal solution $x^*$. Then for all $k > 0$, there exist constants $a > 0$ and $b \in (0, 1)$ such that

Proof. Since $f$ is twice continuously differentiable and uniformly convex on $\Omega$, it follows from (3.2)-(3.4) and (3.12) in ^[24] that there exist constants $\hat{\lambda} > \lambda > 0$, $\hat{\zeta} > \zeta > 0$ such that for all $x \in \Omega$, the following inequalities hold:

Thus, from the first inequality in (2.5), we have

where the second inequality follows from (3.12) and (2.2), the third inequality follows from (3.10), and the last inequality follow from (3.17). Consequently,

Taking $a = f(x_{0})-f(x^*)$ and $b = 1-\rho \dfrac{1-\sigma}{\hat{\lambda} LM^2}$, the desired result (3.16) has been proved.

3.2.   Non-convex cases

For non-convex minimization problems, we can prove that the search directions in (2.1) satisfy an approximate Dai-Liao conjugate condition.

Proposition 3. Suppose that $s_{k}^{T}y_k > \left|s_k^T\bar{y}_{k}\right|$. Then, $d_{k+1}$ in (2.1) satisfies the following approximate Dai-Liao conjugate condition:

Proof. When $s_{k}^{T}y_k > \left|s_k^T\bar{y}_{k}\right|$, it holds that

where $t_k = 1+\dfrac{\|y_{k}\|^2}{s_{k}^{T}y_k} > 0$. The result (3.20) has been proved.

Remark 5. Although the condition (3.20) in Proposition 3 does not always hold at any iteration, it does not affect global convergence of Algorithm 1. By a simple example, we can show that this condition is often satisfied (see Table 1). In addition, our numerical tests will also show advantages of the search directions given by (2.1).

Example 1. For the Rosenbrock problem:

Initial point $x_{0} = (-1.2, 1)$.

We implement Algorithm 1 to solve the Rosenbrock problem, and partly present the obtained values of $s_k^Ty_k$ and $\left|s_k^T\bar{y}_{k}\right|$ in Table 1.

In Table 1, it is easy to see that the directions at the 10th, 12th and 15th iterations, the inequality (3.20) holds.

Before stating global convergence of Algorithm 1 in the non-convex case, we first make the following mild assumptions.

Assumption 1. The level set $\Omega = \{x\in R^n | f(x)\leq f(x_0)\}$ is bounded, i.e., there exists a positive constant $B > 0$ such that (3.8) holds for all $x \in \Omega$.

Assumption 2. In some neighborhood $N$ of $\Omega$, $f$ is continuously differentiable and its gradient is Lipschitz continuous. That is to say, there exists a constant $L > 0$ such that (3.9) holds for all $x, y \in N$.

Theorem 3. Let $\{g_k\}$ be a gradient sequence generated by Algorithm 1. Suppose that there exists a constant $\tau > 0$ such that $s_k^Ty_k\geq\tau$ for any $k\geq1$. Under Assumptions 1 and 2, it is true that

Proof. From Assumptions 1 and 2, we have

where the first and second inequalities follow from the Cauchy-Schwarz inequality and $s_k^Ty_k\geq\tau$, and the last inequality follows from (3.8) and (3.9).

Similar to the proof of Theorem 1, we can also prove that (3.15) holds under Assumptions 1 and 2. Together with (3.23), it is concluded that (3.22) holds.

In order to present R-linear convergence of Algorithm 1 in the non-convex case, we need the following assumption:

Assumption 3. (1) $f: R^{n}\to R$ is twice continuously differentiable.

(2) The sequence $\{x_k\}$ generated by Algorithm 1 satisfies $x_k \rightarrow x^*$, where $\nabla f(x^*) = 0$ and $\nabla^2 f(x^*)$ is positive definite.

(3) There exists a positive constant $\tau$ such that $s_k^Ty_k > \tau$ holds for all the sufficiently large $k > 0$.

Theorem 4. Let $\{x_k\}$ be the sequence generated by Algorithm 1, which converges to a solution $x^*$ satisfying (2) in Assumption 3. Suppose that Assumptions 1 and 2 hold. Then, for the sufficiently large $k > 0$, there exist constants $a > 0$ and $b\in (0, 1)$ such that

Proof. From (3) in Assumption 3, we can know that $\|d_k\|$ is bounded for all $k > 0$. Moreover, from (4.1)–(4.3) in ^[18], it is clear that there exists a neighborhood of $x^*$, denoted by $U(x^*)$, such that (3.17) holds for all $x \in U(x^*)$. The rest of the proof is similar to that of Theorem 2, we omit it here.

4.   Numerical tests

In this section, by numerical experiments, we study effectiveness and robustness of Algorithm 1 when it is employed to solve unconstrained optimization problems.

Algorithm 1 (NTTCG) is tested through solution of the 75 benchmark test problems with variable dimensions from 1000 to 10000. These problems are from ^[25] or CUTE ^[26]. Its computer codes are written using the language of Fortran 77, and run on a personal computer with a 2.2GHZ CPU processor, 8GB memory and Windows 10 operation system.

To show advantages of our algorithm (NTTCG), we compare it with the other four similar algorithms, including TMRMIL in ^[18], ISCG in ^[6], CG_DESCENT in ^[7] and THREECG in ^[8]. For all the compared algorithms, the termination condition is $\varepsilon = 10^{-6}$ or the number of iterations exceeds 10,000. In Algorithm 1, $\rho = 0.0001, \; \; \sigma = 0.01$, and the parameters not mentioned here are consistent with the corresponding literature. We show the numerical performance differences among these five algorithms by the Dolan and Moré performance profiles ^[27]. Let $S$ be a set of all methods, $P$ be a set of test problems, $n_p$ be the size of the set $P$ and $t_{p, s}$ be the number of iterations or the CPU time needed to solve problem $p \in P$ by method $s \in S$. Then, the performance ratio is computed by $r_{p, s} = \dfrac{t_{p, s}}{\min\{t_{p, s}:s \in S\} }$, and the overall performance of Algorithm $s$ is given by $\rho_s(\tau) = \frac{1}{n_p}size\{p\in P:r_{p, s}\leq \tau\}$. In fact, $\rho_s(\tau)$ is the probability for Algorithm $s$ that a performance ratio $r_{p, s}$ is within a factor $\tau \in R$ of the best possible ratio. The function $\rho_s(\tau)$ is the distribution function for the performance ratio $r_{p, s}$. We report the numerical results in Figures 1 and 2. From Figures 1 and 2, we can known that our algorithm (NTTCG) performs the best among the five algorithms, either with respect to the number of iterations, or with respect to the elapsed CPU time.

We also underline good numerical results in Table 3. In Table 3, P, Ni and CPU stand for the number of problems in Table 2, the number of iterations and the consumed CPU time in centisecond (cs), respectively. From the underlined results in Table 3, we know that our algorithm (NTTCG) performs very well in some test problems.

5.   Conclusions

In this paper, we have developed a novel three-term CG algorithm (NTTCG) based on modified gradient-differences for solving unconstrained optimization problems. Global convergence has been proved for this algorithm.

By applying our method to solve the 750 benchmark test problems, the numerical results have demonstrated that NTTCG outperforms the compared four algorithms in the literature. Especially, compared with the existing methods, NTTCG can find the optimal solutions of the unconstrained optimization problems, using less number of iterations, or less CPU time consumed.

In future research, it is valuable to study the method of obtaining the iteration direction by minimizing a quadratic approximate model of the objective function or the conic model in a specific subspace spanned by $g_{k+1}$, $s_{k}$ and $y_{k}$.

It is also interesting to extend the algorithm proposed in this paper to solve nonlinear system of monotone equations since it has been shown in ^{[28,29,30,31,32,33,34]} that recovering sparse signals and restoring blurred images can be formulated as a system of equations, and the CG methods can solve nonlinear system of monotone equations efficiently. Furthermore, as done in ^[35], we can also modify our algorithm to solve symmetric system of nonlinear equations.

Acknowledgments

The authors are very grateful to the editors and referees for their constructive advice.

This research is supported by the National Social Science Foundation of China (Grant No. 21BGL122).

Availability of data and materials

The data and material used to support the findings in this research can be provided by the corresponding author upon request.

Conflict of interest

We declare that all the authors have no any conflicts of interest about this submission and publication of this article.