Unconstrained Optimization Methods: Conjugate Gradient Methods and Trust-Region Methods

Abstract

Here, we consider two important classes of unconstrained optimization methods: conjugate gradient methods and trust region methods. These two classes of methods are very interesting; it seems that they are never out of date. First, we consider conjugate gradient methods. We also illustrate the practical behavior of some conjugate gradient methods. Then, we study trust region methods. Considering these two classes of methods, we analyze some recent results.

1. Introduction

Remind to the unconstrained optimization problem which we can present as

$$\underset{x\in R^n}{min}f\left(x\right),$$ (1) E1

where $f:{\mathbb{R}}^n\to \mathbb{R}$ is a smooth function.

Here, we consider two classes of unconstrained optimization methods: conjugate gradient methods and trust region methods. Both of them are made with the aim to solve the unconstrained optimization problem (1).

In this chapter, at first, we consider the conjugate gradient methods. Then, we study trust region methods. Also, we try to give some of the most recent results in these areas.

2. Conjugate gradient method (shortly CG)

The conjugate gradient method is the method between the steepest descent method and the Newton method.

The conjugate gradient method in fact deflects the direction of the steepest descent method by adding to it a positive multiple of the direction used in the last step.

The restarting and the preconditioning are very important to improve the conjugate gradient method [47].

Some of well-known CG methods are [12, 19, 20, 23, 24, 31, 39, 40, 49]:

$${\beta }_k^{\mathit{HS}}=\frac{g_{k+1}^T{y}_k}{d_k^T{y}_k}$$

$${\beta }_k^{\mathit{FR}}=\frac{\parallel g_{k+1}{\parallel }^2}{\parallel g_k{\parallel }^2}$$

$${\beta }_k^{\mathit{PRP}}=\frac{g_{k+1}^T{y}_k}{\parallel g_k{\parallel }^2}$$

$${\beta }_k^{\mathit{CD}}=\frac{\parallel g_{k+1}{\parallel }^2}{-d_k^T{g}_k}$$

$${\beta }_k^{\mathit{LS}}=\frac{g_{k+1}^T{y}_k}{-d_k^T{g}_k}$$

$${\beta }_k^{\mathit{DY}}=\frac{\parallel g_{k+1}{\parallel }^2}{d_k^T{y}_k}$$

$${\beta }_k^N={\left(y_k-2d_k\frac{\parallel y_k{\parallel }^2}{d_k^T{y}_k}\right)}^T\frac{g_{k+1}}{d_k^T{y}_k}$$

$${\beta }_k^{\mathit{WYL}}=\frac{g_k^T\left(g_k-\frac{\parallel g_k\parallel }{\parallel g_{k-1}\parallel }{g}_{k-1}\right)}{\parallel g_{k-1}{\parallel }^2}$$

Consider positive definite quadratic function

$$f\left(x\right)=\frac{1}{2}{x}^T\mathit{Gx}+b^{T}x+c,$$ (2) E2

where G is an $n\times n$ symmetric positive definite matrix, $b\in {\mathbb{R}}^n$ , and $c$ is a real number.

Theorem 1.2.1. [47] (Property theorem of conjugate gradient method) For positive definite quadratic function (2), FR conjugate gradient method with the exact line search terminates after $m\le n$ steps, and the following properties hold for all $i$ , $0\le i\le m$ :

$$d_i^T{\mathit{Gd}}_j=0,j=0,1,\dots ,i-1;$$

$$g_i^T{g}_j=0,j=0,1,\dots ,i-1;$$

$$d_i^T{g}_i=-g_i^T{g}_i;$$

$$\left[g_0,g_1,\dots ,g_i\right]=\left[g_0,{\mathit{Gg}}_0,\dots ,G^i{g}_0\right];$$

$$\left[d_0,d_1,\dots ,d_i\right]=\left[g_0,{\mathit{Gg}}_0,\dots ,G^i{g}_0\right],$$

where $m$ is the number of distinct eigenvalues of $G$ .

Now, we give the algorithm of conjugate gradient method.

Algorithm 1.2.1. (CG method).

Assumptions: $\epsilon <0$ and $x_0\in {\mathbb{R}}^n$ . Let $k=0$ , $t_0=0$ , $d_{-1}=0$ , $d_0=-g_0$ , ${\beta }_{-1}=0$ , and ${\beta }_0=0$ .

Step 1. If $\parallel g_k\parallel \le \epsilon$ , then STOP.

Step 2. Calculate the step-size $t_k$ by a line search.

Step 3. Calculate ${\beta }_k$ by any of the conjugate gradient method.

Step 4. Calculate $d_k=-g_k+{\beta }_{k-1}{d}_{k-1}$ .

Step 5. Set $x_{k+1}=x_k+t_k{d}_k$ .

Step 6. Set $k=k+1$ and go to Step 1.

2.1 Convergence of conjugate gradient methods

Theorem 1.2.2. [47] (Global convergence of FR conjugate gradient method) Suppose that $f:{\mathbb{R}}^n\to \mathbb{R}$ is continuously differentiable on a bounded level set

$$L=\left\{x\in {\mathbb{R}}^n|f\left(x\right)\le f\left(x_0\right)\right\},$$

and let FR method be implemented by the exact line search. Then, the produced sequence $\left\{x_k\right\}$ has at least one accumulation point, which is a stationary point, i.e.:

1. When $\left\{x_k\right\}$ is a finite sequence, then the final point $x^{\ast }$ is a stationary point of $f$ .

2. When $\left\{x_k\right\}$ is an infinite sequence, then it has a limit point, and it is a stationary point.

In [35], a comparison of two methods, the steepest descent method and the conjugate gradient method which are used for solving systems of linear equations, is illustrated. The aim of the research is to analyze, which method is faster in solving these equations and how many iterations are needed by each method for solving.

The system of linear equations in the general form is considered:

$$\mathit{Ax}=B,$$ (3) E3

where matrix $A$ is symmetric and positive definite.

The conclusion is that the $\mathit{SD}$ method is a faster method than the $\mathit{CG}$ , because it solves equations in less amount of time.

By the other side, the authors find that the $\mathit{CG}$ method is slower but more productive than the $\mathit{SD}$ , because it converges after less iterations.

So, we can see that one method can be used when we want to find solution very fast and another can converge to maximum in less number of iterations.

Again, we consider the problem (1), where $f:{\mathbb{R}}^n\to \mathbb{R}$ is a smooth function and its gradient is available.

A hybrid conjugate gradient method is a certain combination of different conjugate gradient methods; it is made to improve the behavior of these methods and to avoid the jamming phenomenon.

An excellent survey of hybrid conjugate gradient methods is given in [5].

Three-term conjugate gradient methods were studied in the past (e.g., see [8, 32, 34], etc.); but, from recent papers about $\mathit{CG}$ methods, we can conclude that maybe the mainstream is made by three-term and even four-term conjugate gradient methods. An interesting paper about a five-term hybrid conjugate gradient method is [1]. Also, from recent papers we can conclude that different modifications of the existing $\mathit{CG}$ methods are made, as well as different hybridizations of $\mathit{CG}$ and $\mathit{\text{BFGS}}$ methods.

Consider unconstrained optimization problem (1), where $f:{\mathbb{R}}^n\to \mathbb{R}$ is a continuously differentiable function, bounded from below. Starting from an initial point $x_0\in {\mathbb{R}}^n$ , the three-term conjugate gradient method with line search generates a sequence $\left\{x_k\right\}$ , given by the next iterative scheme:

$$x_{k+1}=x_k+t_k{d}_k,$$ (4) E4

where $t_k$ is a step-size which is obtained from the line search, and

$$d_0=-g_0,d_{k+1}=-g_{k+1}+{\delta }_k{s}_k+{\eta }_k{y}_k.$$

In the last relation, ${\delta }_k$ and ${\eta }_k$ are the conjugate gradient parameters, $s_k=x_{k+1}-x_k$ , $g_k=\nabla f\left(x_k\right)$ , and $y_k=g_{k+1}-g_k$ . We can see that the search direction $d_{k+1}$ is computed as a linear combination of $-g_{k+1}$ , $s_k$ , and $y_k$ .

In [6], the author suggests another way to get three-term conjugate gradient algorithms by minimization of the one-parameter quadratic model of the function $f$ . The idea is to consider the quadratic approximation of the function $f$ in the current point and to determine the search direction by minimization of this quadratic model. It is assumed that the symmetrical approximation of the Hessian matrix $B_{k+1}$ satisfies the general quasi-Newton equation which depends on a positive parameter:

$$B_{k+1}{s}_k={\omega }^{-1}{y}_k,\omega \ne 0.$$ (5) E5

In this paper the quadratic approximation of the function $f$ is considered:

$${\mathrm{\Phi }}_{k+1}\left(d\right)=f_{k+1}+g_{k+1}^{T}d+\frac{1}{2}{d}^T{B}_{k+1}d.$$

The direction $d_{k+1}$ is computed as

$$d_{k+1}=-g_{k+1}+{\beta }_k{s}_k,$$ (6) E6

where the scalar ${\beta }_k$ is determined as the solution of the following minimizing problem:

$$\underset{{\beta }_k\in \mathbb{R}}{min}{\mathrm{\Phi }}_{k+1}\left(d_{k+1}\right).$$ (7) E7

From (6) and (7), the author obtains

$${\beta }_k=\frac{g_{k+1}^T{B}_{k+1}{s}_k-g_{k+1}^T{s}_k}{s_k^T{B}_{k+1}{s}_k}.$$ (8) E8

Using (5), from (7), the next expression for ${\beta }_k$ is obtained:

$${\beta }_k=\frac{g_{k+1}^T{y}_k-\omega g_{k+1}^T{s}_k}{y_k^T{s}_k}.$$ (9) E9

Using the idea of Perry [36], the author obtains

$$d_{k+1}=-g_{k+1}+\frac{y_k^T{g}_{k+1}-\omega s_k^T{g}_{k+1}}{y_k^T{s}_k}{s}_k-\frac{s_k^T{g}_{k+1}}{y_k^T{s}_k}{y}_k.$$

In fact, in this approach the author gets a family of three-term conjugate gradient algorithms depending of a positive parameter $\omega$ .

Next, in [52], the $\mathit{WYL}$ conjugate gradient ( $\mathit{CG}$ ) formula, with ${\beta }_k^{\mathit{WYL}}\ge 0$ , is further studied. A three-term $\mathit{WYL}\mathit{CG}$ algorithm is presented, which has the sufficiently descent property without any conditions. The global convergence and the linear convergence are proven; moreover, the $n$ -step quadratic convergence with a restart strategy is established if the initial step length is appropriately chosen.

The first three-term Hestenes-Stiefel ( $\mathit{HS}$ ) method ( $\mathit{\text{TTHS}}$ method) can be found in [55].

Baluch et al. [7] describe a modified three-term Hestenes-Stiefel ( $\mathit{HS}$ ) method. Although the earliest conjugate gradient method $\mathit{HS}$ achieves global convergence using an exact line search, this is not guaranteed in the case of an inexact line search. In addition, the $\mathit{HS}$ method does not usually satisfy the descent property. The modified three-term conjugate gradient method from [7] possesses a sufficient descent property regardless of the type of line search and guarantees global convergence using the inexact Wolfe-Powell line search [50, 51]. The authors also prove the global convergence of this method. The search direction, which is considered in [7], has the next form:

$$d_k=\left\{\begin{array}{l}-g_k,\text{if}k=0,\\ -g_k+{\beta }_k^{\mathit{BZA}}{d}_{k-1}-{\theta }_k^{\mathit{BZA}}{y}_{k-1},\text{if}k\ge 1,\end{array}\right.$$

where ${\beta }_k^{\mathit{BZA}}=\frac{g_k^T\left(g_k-g_{k-1}\right)}{d_{k-1}^T{y}_{k-1}+\mu \mid g_k^T{d}_{k-1}\mid }$ , ${\theta }_k^{\mathit{BZA}}=\frac{g_k^T{d}_{k-1}}{d_{k-1}^T{y}_{k-1}+\mu \mid g_k^T{d}_{k-1}\mid },$ $\mu >1.$

In [13], an accelerated three-term conjugate gradient method is proposed, in which the search direction satisfies the sufficient descent condition as well as extended Dai-Liao conjugacy condition:

$$d_k^T{y}_{k-1}=-{\mathit{tg}}_k^T{s}_{k-1},t\ge 0.$$

This method seems different from the existent methods.

Next, Li-Fushikuma quasi-Newton equation is

$${\nabla }^{2}f\left(x_k\right)s_{k-1}=z_{k-1},$$ (10) E10

where

$$z_{k-1}=y_{k-1}+C\parallel g_{k-1}{\parallel }^r{s}_{k-1}+max\left\{-\frac{s_{k-1}^T{y}_{k-1}}{\parallel s_{k-1}{\parallel }^2},0\right\}{s}_{k-1},$$

where $C$ and $r$ are two given positive constants. Based on (10), Zhou and Zhang [56] propose a modified version of $\mathit{DL}$ method, called $\mathit{ZZ}$ method in [13].

In [30], some new conjugate gradient methods are extended, and then some three-term conjugate gradient methods are constructed. Namely, the authors remind to [41, 42], with its conjugate gradient parameters, respectively:

$${\beta }_k^{\mathit{\text{RMIL}}}=\frac{g_k^T{y}_{k-1}}{\parallel d_{k-1}{\parallel }^2},$$ (11) E11

$${\beta }_k^{\mathit{\text{MRMIL}}}=\frac{g_k^T\left(g_k-g_{k-1}-d_{k-1}\right)}{\parallel d_{k-1}{\parallel }^2},$$ (12) E12

wherefrom it is obvious that ${\beta }_k^{\mathit{\text{MRMIL}}}={\beta }_k^{\mathit{\text{RMIL}}}$ for the exact line search. Let us say that these methods, presented in [41, 42], are $\mathit{\text{RMIL}}$ and $\mathit{\text{MRMIL}}$ methods.

The three-term $\mathit{\text{RMIL}}$ and $\mathit{\text{MRMIL}}$ methods are introduced in [30].

The search direction $d_k$ can be expressed as

$$d_0=-g_0,d_k=-g_k+{\beta }_k{d}_{k-1}+{\theta }_k{y}_{k-1},$$

where ${\beta }_k$ is given by (11) or (12), and

$${\theta }_k=-\frac{g_k^T{d}_{k-1}}{\parallel d_{k-1}{\parallel }^2}.$$

An important property of the proposed methods is that the search direction always satisfies the sufficient descent condition without any line search, that is, the next relation always holds

$$g_k^T{d}_k\le -\parallel g_k{\parallel }^2.$$

Under the standard Wolfe line search and the classical assumptions, the global convergence properties of the proposed methods are proven.

Having in view the conjugate gradient parameter suggested in [49], in [45] the next two conjugate gradient parameters are presented:

$${\beta }_k^{\mathit{MHS}}=\frac{\parallel g_k{\parallel }^2-\frac{\parallel g_k\parallel }{\parallel g_{k-1}\parallel }{g}_k^T{g}_{k-1}}{d_{k-1}^T\left(g_k-g_{k-1}\right)},$$ (13) E13

$${\beta }_k^{\mathit{MLS}}=\frac{\parallel g_k{\parallel }^2-\frac{\parallel g_k\parallel }{\parallel g_{k-1}\parallel }{g}_k^T{g}_{k-1}}{-d_{k-1}^T{g}_{k-1}}.$$ (14) E14

Motivated by [49], as well as by [45], in [1], a new hybrid nonlinear CG method is proposed; it combines the features of five different CG methods, with the aim of combining the positive features of different non-hybrid methods. The proposed method generates descent directions independently of the line search. Under some assumptions on the objective function, the global convergence is proven under the standard Wolfe line search. Conjugate gradient parameter, proposed in [1], is

$${\beta }_k^{\mathit{hAO}}=\frac{\parallel g_k{\parallel }^2-max\left\{0,\frac{\parallel g_k\parallel }{\parallel g_{k-1}\parallel }{g}_k^T{g}_{k-1}\right\}}{max\left\{\parallel g_{k-1}{\parallel }^2,d_{k-1}^T\left(g_k-g_{k-1}\right),-d_{k-1}^T{g}_{k-1}\right\}}.$$ (15) E15

Let’s note that the proposed method is hybrid of $\mathit{FR}$ , $\mathit{DY}$ , $\mathit{WYL}$ , $\mathit{MHS}$ , and $\mathit{MLS}$ .

The behaviors of the methods $\mathit{BZA}$ , $\mathit{\text{TTRMIL}}$ , $\mathit{\text{MRMIL}}$ , $\mathit{MHS}$ , $\mathit{MLS}$ , and $\mathit{hAO}$ are illustrated by the next tables.

The test criterion is CPU time.

The tests are performed on the computer Workstation Intel Celeron CPU 1,9 GHz.

The experiments are made on the test functions from [3].

Each problem is tested for a number of variables $n=1000$ and $n=5000$ .

The average CPU time values are given in the last rows of these tables (Tables 1, 2, 3, 4).

function	BZA	TTRMIL	MRMIL	MHS	MLS	hAO
Ext.Pen.	21.793340	20.966534	16.036903	19.812127	21.933741	20.326930
Pert.Quad.	21.855740	22.542144	15.506499	20.904134	22.230142	18.954121
Raydan1	6.801644	7.066845	6.349241	7.098045	7.066845	7.332047
Raydan2	0.608404	0.592804	0.577204	0.592804	0.608404	0.639604
Diag.1	0.608404	0.608404	0.577204	0.608404	0.514803	0.577204
Diag.2	5.163633	5.600436	4.695630	4.758031	5.662836	4.851631
Diag.3	5.616036	5.694037	5.241634	5.756437	5.584836	5.506835
Gen.Tridiag.-1	3.042019	2.932819	2.683217	2.948419	2.792418	2.808018
Hager	2.917219	2.932819	2.620817	3.042019	2.917219	2.886019
Ext.Tridiag.-1	2.886019	2.932819	2.761218	2.932819	2.730018	2.917219
Ext.ThreeExp.	2.979619	2.964019	2.605217	2.886019	3.042019	2.714417
Diag.4	2.901619	2.870418	2.574016	2.792418	2.948419	2.652017
Diag.5	2.792418	2.917219	2.574016	2.901619	3.026419	2.901619
Ext.Himm.	2.761218	2.714417	2.667617	2.964019	2.995219	2.854818
Ext.PSC1	2.932819	2.745618	2.714417	2.511616	3.026419	2.792418
FullHess.FH2	2.870418	2.948419	2.886019	2.839218	3.010819	2.948419
Ext.Bl.Diag.BD1	2.979619	2.886019	2.948419	2.886019	2.901619	2.542816
Quad.QF1	2.854818	2.870418	3.057620	2.964019	2.964019	2.886019
Ext.Quad.Pen.QP1	2.948419	2.808018	2.605217	2.964019	2.823618	2.542816
Quad.QF2	2.839218	2.620817	2.886019	2.979619	2.901619	2.683217
Ext.EP1	2.730018	2.402415	2.932819	2.698817	2.792418	2.652017
Ext.Tridiag.-2	2.683217	2.605217	2.839218	2.870418	2.886019	2.542816
Tridia	2.683217	2.511616	2.964019	2.823618	2.823618	2.511616
Arwhead	2.917219	2.995219	2.745618	2.823618	2.745618	2.012413
Dqdrtic	2.761218	2.995219	2.901619	2.823618	2.730018	2.589617
Quartc(Cute)	2.886019	2.776818	2.886019	2.776818	2.870418	2.839218
Dixon3dq(Cute)	2.808018	2.948419	2.948419	2.839218	2.917219	2.605217

Table 1.

n = 1000.

function	BZA	TTRMIL	MRMIL	MHS	MLS	hAO
Biggsb1(Cute)	2.792418	2.870418	2.870418	2.917219	2.979619	2.901619
Gen.quart.	2.917219	2.932819	2.464816	2.948419	2.808018	2.620817
Diag.7	2.574016	2.589617	2.870418	2.620817	3.026419	2.698817
Diag.8	2.730018	2.979619	2.839218	2.964019	2.792418	2.979619
Full Hess.FH3	2.948419	2.574016	2.698817	3.026419	2.636417	2.745618
Himmelbg	2.854818	3.010819	2.901619	2.854818	2.995219	2.730018
Ext.Pow.	2.901619	2.854818	2.761218	2.808018	2.870418	2.995219
Ext.Maratos	2.854818	2.948419	2.870418	2.995219	2.870418	2.917219
Ext.Cliff	2.964019	3.042019	2.854818	2.932819	2.886019	2.854818
Pert.quad.diag.	2.714417	3.104420	2.683217	2.964019	2.667617	2.901619
Ext.Wood	2.995219	2.932819	2.948419	2.948419	2.964019	2.948419
Ext.Trigon.	2.792418	2.995219	2.839218	3.010819	2.995219	2.745618
Ext.Rosenbr.	2.964019	2.839218	2.948419	2.932819	2.995219	2.776818
Average	3.915625	3.928105	3.533423	3.868045	3.973345	3.722184

Table 2.

n = 1000.

function	BZA	TTRMIL	MRMIL	MHS	MLS	hAO
Ext.Pen.	46.160696	46.831500	48.656712	66.284825	65.863622	63.695208
Pert.Quad.	48.375910	45.801894	52.307135	66.612427	66.113224	65.551620
Raydan1	12.994883	12.105678	13.759288	16.972909	16.598506	16.754507
Raydan2	1.170008	1.029607	1.076407	1.154407	1.092007	1.107607
Diag.1	8.845257	0.904806	1.076407	1.123207	1.170008	1.092007
Diag.2	8.658055	7.831250	7.924851	9.094858	10.358466	10.327266
Diag.3	8.361654	9.141659	8.673656	10.686068	10.358466	10.514467
Gen.Tridiag.-1	5.616036	5.382034	5.865638	6.021639	6.489642	6.364841
Hager	5.241634	4.851631	5.881238	6.286840	5.304034	6.021639
Ext.Tridiag.-1	5.007632	4.804831	5.740837	5.787637	6.224440	5.803237
Ext.ThreeExp.	4.882831	4.820431	5.522435	6.115239	6.333641	5.834437
Diag.4	4.929632	4.898431	5.179233	5.803237	6.177640	6.427241
Diag.5	5.694037	4.851631	5.538036	5.709637	5.896838	6.115239
Ext.Himm.	5.834437	5.116833	5.382034	6.099639	5.772037	6.411641
Ext.PSC1	5.023232	5.054432	5.163633	6.411641	6.115239	5.990438
FullHess.FH2	5.210433	4.929632	4.851631	6.068439	6.349241	6.349241
Ext.Bl.Diag.BD1	4.851631	5.007632	5.226033	6.364841	6.364841	5.569236
Quad.QF1	5.475635	5.662836	6.302440	6.177640	6.146439	6.286840
Ext.Quad.Pen.QP1	5.226033	5.163633	4.929632	6.130839	5.818837	5.943638
Quad.QF2	5.335234	4.836031	5.990438	6.084039	6.084039	6.084039
Ext.EP1	5.070032	5.038832	6.052839	6.115239	4.992032	6.177640
Ext.Tridiag.-2	4.851631	4.976432	4.851631	6.349241	5.990438	6.099639
Tridia	5.413235	4.820431	5.475635	5.569236	5.818837	6.021639
Arwhead	4.867231	4.882831	5.023232	6.099639	6.380441	6.177640
Dqdrtic	5.163633	4.945232	5.023232	5.428835	6.006038	5.850038
Quartc(Cute)	5.912438	5.350834	5.834437	5.787637	5.896838	6.193240
Dixon3dq(Cute)	5.428835	4.789231	5.163633	6.162039	5.616036	5.881238

Table 3.

n = 5000.

function	BZA	TTRMIL	MRMIL	MHS	MLS	hAO
Biggsb1(Cute)	5.148033	4.695630	5.413235	5.912438	6.052839	6.349241
Gen.quart.	5.288434	4.758031	5.023232	6.349241	6.052839	4.960832
Diag.7	5.163633	4.664430	5.054432	5.959238	6.193240	6.255640
Diag.8	5.787637	4.742430	4.898431	6.099639	5.600436	6.208840
Full Hess.FH3	5.444435	4.789231	5.569236	6.177640	6.162039	6.224440
Himmelbg	5.584836	6.130839	5.475635	5.475635	6.006038	5.912438
Ext.Pow.	5.569236	4.789231	4.773631	5.990438	5.772037	6.162039
Ext.Maratos	5.148033	5.740837	4.976432	6.021639	6.286840	6.130839
Ext.Cliff	5.943638	5.850038	4.976432	5.990438	5.304034	6.286840
Pert.quad.diag.	5.912438	6.427241	4.976432	6.318041	6.115239	6.068439
Ext.Wood	5.584836	5.647236	4.789231	6.255640	5.350834	6.021639
Ext.Trigon.	5.366434	5.709637	4.773631	6.115239	6.021639	5.787637
Ext.Rosenbr.	6.177640	5.319634	4.617630	6.333641	6.021639	6.021639
Average	7.79302995	7.327367	7.694749	9.287519525	9.206789	9.225899

Table 4.

n = 5000.

In [2], based on the numerical efficiency of Hestenes-Stiefel ( $\mathit{HS}$ ) method, a new modified $\mathit{HS}$ algorithm is proposed for unconstrained optimization. The new direction independent of the line search satisfies the sufficient descent condition. Motivated by theoretical and numerical features of three-term conjugate gradient ( $\mathit{CG}$ ) methods proposed by [33], similar to the approach in [10], the new direction is computed by minimizing the distance between the $\mathit{CG}$ direction and the direction of the three-term $\mathit{CG}$ methods proposed by [33]. Under some mild conditions, the global convergence of the new method for general functions is established when the standard Wolfe line search is used. In this paper the conjugate gradient parameter is given by

$${\beta }_k={\beta }_k^{\mathit{HS}}{\theta }_k,$$ (16) E16

where

$${\theta }_k=1-\frac{{\left(g_k^T{d}_{k-1}\right)}^2}{\parallel g_k{\parallel }^2\parallel d_{k-1}{\parallel }^2}.$$

But this new CG direction does not fulfill a descent condition, so further modification is made, namely, having in view [53], the authors [2] introduce

$${\overline{\beta }}_k={\beta }_k-\lambda {\left(\frac{\parallel y_{k-1}\parallel {\theta }_k}{d_{k-1}^T{y}_{k-1}}\right)}^2{g}_k^T{d}_{k-1},$$

where $\lambda >\frac{1}{4}$ is a parameter. Also, the global convergence is proven under standard conditions.

It is worth to mention the next papers about this theme, which can be interesting [4, 14, 15, 16, 17, 25, 26, 27].

3. Trust region methods

We remind that the basic idea of Newton method is to approximate the objective function $f\left(x\right)$ around $x_k$ by using a quadratic model:

$$q^{\left(k\right)}\left(s\right)=f\left(x_k\right)+g_k^{T}s+\frac{1}{2}{s}_k^T{G}_{k}s,$$

where $g_k=\nabla f\left(x_k\right)$ , $G_k={\nabla }^{2}f\left(x_k\right)$ , and also use the minimizer $s_k$ of $q^{\left(k\right)}\left(s\right)$ to set $x_{k+1}=x_k+s_k$ .

Also, remind that Newton method can only guarantee the local convergence, i.e., when $s$ is small enough and the method is convergent locally.

Further, Newton method cannot be used when Hessian is not positive definite.

There exists another class of methods, known as trust region methods. It does not use the line search to get the global convergence, as well as it avoids the difficulty which is the consequence of the nonpositive definite Hessian in the line search.

Furthermore, it produces greater reduction of the function $f$ than line search approaches.

Here, we define the region around the current iterate:

$${\mathrm{\Omega }}_k=\left\{x:\parallel x-x_k\parallel \le {\mathrm{\Delta }}_k\right\},$$

where ${\mathrm{\Delta }}_k$ is the radius of ${\mathrm{\Omega }}_k$ , inside which the model is trusted to be adequate to the objective function.

Our further intention is to choose a step which should be the approximate minimizer of the quadratic model in the trust region. In fact, $x_k+s_k$ should be the approximately best point on the sphere:

$$\left\{x_k+s|\parallel s\parallel \le {\mathrm{\Delta }}_k\right\},$$

with the center $x_k$ and the radius ${\mathrm{\Delta }}_k$ .

In the case that this step is not acceptable, we reduce the size of the step, and then we find a new minimizer.

This method has the rapid local convergence rate, and that’s the property of Newton method and quasi-Newton method, too, but the important characteristic of trust region method is also the global convergence.

Since the step is restricted by the trust region, this method is also called the restricted step method.

The model subproblem of the trust region method is

$$min{q}^{\left(k\right)}\left(s\right)=f\left(x_k\right)+g_k^{T}s+\frac{1}{2}{s}^T{B}_{k}s,$$ (17) E17

$$\mathrm{s}.\mathrm{t}.\parallel s\parallel \le {\mathrm{\Delta }}_k,$$ (18) E18

where ${\mathrm{\Delta }}_k$ is the trust region radius and $B_k$ is a symmetric approximation of the Hessian $G_k$ .

In the case that we use the standard $l_2$ norm $\parallel \cdot {\parallel }_2$ , $s_k$ is the minimizer of $q^{\left(k\right)}\left(s\right)$ in the ball of radius ${\mathrm{\Delta }}_k$ . Generally, different norms define the different shapes of the trust region.

Setting $B_k=G_k$ in (17)–(18), the method becomes a Newton-type trust region method.

The problem by itself is the choice of ${\mathrm{\Delta }}_k$ at each single iteration.

If the agreement between the model $q^{\left(k\right)}\left(s\right)$ and the objective function $f\left(x_k+s\right)$ is satisfactory enough, the value ${\mathrm{\Delta }}_k$ should be chosen as large as it is possible. The expression ${\mathit{\text{Ared}}}_k=f\left(x_k\right)-f\left(x_k+s_k\right)$ is called the actual reduction, and the expression ${\mathit{\text{Pred}}}_k=q^{\left(k\right)}\left(0\right)-q^{\left(k\right)}\left(s_k\right)$ is called the predicted reduction; here, we emphasize that

$$r_k=\frac{{\mathit{\text{Ared}}}_k}{{\mathit{\text{Pred}}}_k}$$

measures the agreement between the model function $q^{\left(k\right)}\left(s\right)$ and the objective function $f\left(x_k+s\right)$ .

If $r_k$ is close to 0 or it is negative, the trust region is going to shrink; otherwise, we do not change the trust region.

The conclusion is that $r_k$ is important in making the choice of new iterate $x_{k+1}$ as well as in updating the trust region radius ${\mathrm{\Delta }}_k$ . Now, we give the trust region algorithm.

Algorithm 1.3.1. (Trust region method).

Assumptions: $x_0$ , $\overline{\mathrm{\Delta }}$ , ${\mathrm{\Delta }}_0\in \left(0,\overline{\mathrm{\Delta }}\right)$ , $\epsilon \ge 0$ , $0<{\eta }_1\le {\eta }_2<1$ , and $0<{\gamma }_1<1<{\gamma }_2$ .

Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \epsilon$ , then STOP.

Step 2. Approximately solve the problem (17)–(18) for $s_k$ .

Step 3. Compute $f\left(x_k+s_k\right)$ and $r_k$ . Set

$$x_{k+1}=\left\{\begin{array}{l}{x}_k+s_k,\mathit{\text{if}}{r}_k\ge {\eta }_1,\\ x_k,\mathit{\text{otherwise}}.\end{array}\right.$$

Step 4. If $r_k<{\eta }_1$ , then ${\mathrm{\Delta }}_{k+1}\in \left(0,{\gamma }_1{\mathrm{\Delta }}_k\right)$ .

If $r_k\in \left[{\eta }_1,{\eta }_2\right)$ , then ${\mathrm{\Delta }}_{k+1}\in \left({\gamma }_1{\mathrm{\Delta }}_k,{\mathrm{\Delta }}_k\right)$ .

If $r_k\ge {\eta }_2$ and $\parallel s_k\parallel ={\mathrm{\Delta }}_k$ , then ${\mathrm{\Delta }}_{k+1}\in \left[{\mathrm{\Delta }}_k,min\left\{{\gamma }_2{\mathrm{\Delta }}_k,\overline{\mathrm{\Delta }}\right\}\right]$ .

Step 5. Generate $B_{k+1}$ , update $q^{\left(k\right)}$ , set $k=k+1$ , and go to Step 1.

In Algorithm 1.3.1, $\overline{\mathrm{\Delta }}$ is a bound for all ${\mathrm{\Delta }}_k$ . Those iterations with the property $r_k\ge {\eta }_2$ (and so those for which ${\mathrm{\Delta }}_{k+1}\ge {\mathrm{\Delta }}_k$ ) are called very successful iterations; the iterations with the property $r_k\ge {\eta }_1$ (and so those for which $x_{k+1}=x_k+s_k$ ) are called successful iterations; and the iterations with the property $r_k<{\eta }_1$ (and so those for which $x_{k+1}=x_k$ ) are called unsuccessful iterations. Generally, the iterations from the two first cases are called successful iterations.

Some choices of parameters are ${\eta }_1=0,01$ , ${\eta }_2=0,75$ , ${\gamma }_1=0,5$ , ${\gamma }_2=2$ , ${\mathrm{\Delta }}_0=1$ , and ${\mathrm{\Delta }}_0=\frac{1}{10}\parallel g_0\parallel$ . The algorithm is insensitive to change of these parameters.

Next, if $r_k<0,01$ , then ${\mathrm{\Delta }}_{k+1}$ can be chosen from $\left(\mathrm{0.01,0.5}\right)\parallel s_k\parallel$ on the basis of a polynomial interpolation.

In the case of quadratic interpolation, we set

$${\mathrm{\Delta }}_{k+1}=\lambda \parallel s_k\parallel ,$$

where

$$\lambda =\frac{-g_k^T{s}_k}{2\left(f\left(x_k+s_k\right)-f\left(x_k\right)-g_k^T{s}_k\right)}.$$

3.1 Convergence of trust region methods

Assumption 1.3.1 (Assumption $A_0$ ).

We assume that the approximations of Hessian $\left\{B_k\right\}$ are uniformly bounded in norm and the level set $L=\left\{x|f\left(x\right)\le f\left(x_0\right)\right\}$ is bounded, as well as $f:{\mathbb{R}}^n\to \mathbb{R}$ is continuously differentiable on $L$ . We allow the length of the approximate solution $s_k$ of the subproblem (17)–(18) to exceed the bound of the trust region, but we also assume that

$$\parallel s_k\parallel \le \tilde{\eta }{\mathrm{\Delta }}_k,$$

where $\tilde{\eta }$ is a positive constant.

In this kind of trust region way of thinking, generally we do not seek an accurate solution of the subproblem (17)–(18); we are satisfied by finding a nearly optimal solution of the subproblem (17)–(18).

Strong theoretical as well as numerical results can be obtained if the step $s_k$ , produced by Algorithm 1.3.1, satisfies

$$q_k\left(0\right)-q_k\left(s_k\right)\ge {\beta }_1\parallel g_k{\parallel }_{2}min\left\{{\mathrm{\Delta }}_k,\frac{\parallel g_k{\parallel }_2}{\parallel B_k{\parallel }_2}\right\},{\beta }_1\in \left(0,1\right).$$

Theorem 1.3.1 [47] Under Assumption $A_0$ , if Algorithm 3.1 has finitely many successful iterations, then it converges to the first-order stationary point.

Theorem 1.3.2 [47] Under Assumption $A_0$ , if Algorithm 3.1 has infinitely many successful iterations, then

$$\underset{k\to \infty }{liminf}\parallel g_k\parallel =0.$$

In [44], it is emphasized that trust region methods are very effective for optimization problems and a new adaptive trust region method is presented. This method combines a modified secant equation with the $\mathit{\text{BFGS}}$ update formula and an adaptive trust region radius, where the new trust region radius makes use of not only the function information but also the gradient information. Let ${\widehat{B}}_k$ be a positively definite matrix based on modified Cholesky factorization [43]. Under suitable conditions, in [44] the global convergence is proven; also, the local superlinear convergence of the proposed method is demonstrated. Motivated by the adaptive technique, the proposed method possesses the following nice properties:

1. The trust region radius uses not only the gradient value but also the function value.

2. Computing the matrix ${\widehat{B}}_k$ of the inverse and the value of $\parallel {{\widehat{B}}_k}^{-1}\parallel$ , at each iterative point $x_k$ , is not required.

3. The computational time is reduced.

A modified secant equation is introduced:

$$B_{k+1}{d}_k=q_k,$$ (19) E19

where $q_k=y_k+h_k{d}_k$ , $f_k=f\left(x_k\right)$ , and $h_k=\frac{{\left(g_{k+1}+g_k\right)}^T{d}_k+2\left(f_k-f_{k+1}\right)}{\parallel d_k{\parallel }^2}$ .

When $f$ is twice continuously differentiable and $B_{k+1}$ is generated by the $\mathit{\text{BFGS}}$ formula, where $B_0=I$ , this modified secant Eq. (19) possesses the following nice property:

$$f_k=f_{k+1}-g_{k+1}^T{d}_k+\frac{1}{2}{d}_k^T{B}_{k+1}{d}_k,$$

and this property holds for all $k$ .

Under classical assumptions, the global convergence of the method presented in [44] is also proven in this paper.

In [28], the hybridization of monotone and non-monotone approaches is made; a modified trust region ratio is used, in which more information is provided about the agreement between the exact and the approximate models. An adaptive trust region radius is used, as well as two accelerated Armijo-type line search strategies to avoid resolving the trust region subproblem whenever a trial step is rejected. It is shown that the proposed algorithm is globally and locally superlinearly convergent. In this paper trust region methods are denoted shortly by $\mathit{TR}$ ; it is emphasized that in $\mathit{TR}$ method, having in view that the iterative scheme is

$$x_0\in {\mathbb{R}}^n,x_{k+1}=x_k+s_k,k=0,1,\dots ,$$

and it often happens that $s_k$ is an approximate solution of the following quadratic subproblem:

$$\underset{s\in {\mathbb{R}}^n,\parallel s_k\parallel \le {\mathrm{\Delta }}_k}{min}{m}_k\left(s\right)=g_k^{T}s+\frac{1}{2}{s}_k^T{B}_{k}s.$$ (20) E20

Performance of the $\mathit{TR}$ methods is much influenced by the strategy of choosing the $\mathit{TR}$ radius at each iteration. To determine the radius ${\mathrm{\Delta }}_k$ , in the standard $\mathit{TR}$ method, the agreement between $f\left(x_k+s\right)$ and $m_k\left(s\right)$ is evaluated by the so-called $\mathit{TR}$ ratio ${\rho }_k$ :

$${\rho }_k=\frac{f\left(x_k\right)-f\left(x_k+s_k\right)}{m_k\left(0\right)-m_k\left(s_k\right)}.$$

When ${\rho }_k$ is negative or a small positive number near to zero, the quadratic model is a poor approximation of the objective function. In such situation, ${\mathrm{\Delta }}_k$ should be decreased and, consequently, the subproblem (20) should be solved again. However, when ${\rho }_k$ is close to $1$ , it is reasonable to use the quadratic model as an approximation of the objective function. So, the step $s_k$ should be accepted and ${\mathrm{\Delta }}_k$ can be increased. Here, the authors use the modified version of ${\rho }_k$ :

$${\overline{\rho }}_k=\frac{R_k-f\left(x_k+s_k\right)}{P_k-m_k\left(s_k\right)},$$

where $R_k={\eta }_k{f}_{l\left(k\right)}+\left(1-{\eta }_k\right)f_k$ , ${\eta }_k\in \left[{\eta }_{\mathit{min}},{\eta }_{\mathit{max}}\right]$ , ${\eta }_{\mathit{min}}\in \left[0,1\right)$ , and ${\eta }_{\mathit{max}}\in \left[{\eta }_{\mathit{min}},1\right]$ . Also,

$$f_{l\left(k\right)}=\underset{0\le j\le q\left(k\right)}{max}\left\{f_{k-j}\right\},f_i=f\left(x_i\right),q\left(0\right)=0,0\le q\left(k\right)\le min\left\{q\left(k-1\right)+1,N\right\},$$

where $N\in \mathbb{N}$ which is originally used by Toint [48].

Something more about trust region methods can be found in [9, 18, 21, 22, 54].

4. Conclusion

The conjugate gradient methods and trust region methods are very popular now.

Many scientists consider these methods.

Namely, different modifications of these methods are made, with the aim to improve them.

Next, the scientists try to make not only new methods but also whole new classes of methods. For the specific values of the parameters, individual methods are distinguished from these classes. It is always more desirable to make a class of methods instead of individual methods.

Hybrid conjugate gradient methods are made in many different ways; this class of conjugate gradient methods is always actual.

Further, one of the contemporary trends is to use BFGS update in constructions of new conjugate gradient methods (e.g., see [46]).

Finally, let us emphasize that contemporary papers often use the Picard-Mann-Ishikawa iterative processes and they make the connection of these kinds of processes with the unconstrained optimization (see [29, 37, 38]).