Some Unconstrained Optimization Methods

Abstract

Although it is a very old theme, unconstrained optimization is an area which is always actual for many scientists. Today, the results of unconstrained optimization are applied in different branches of science, as well as generally in practice. Here, we present the line search techniques. Further, in this chapter we consider some unconstrained optimization methods. We try to present these methods but also to present some contemporary results in this area.

1. Introduction

Optimization is a very old subject of a great interest; we can search deep into a human history to find important examples of applying optimization in the usual life of a human being, for example, the need of finding the best way to produce food yielded finding the best piece of land for producing, as well as (later on, how the time was going) the best ways of treatment of the chosen land and the chosen seedlings to get the best results.

From the very beginning of manufacturing, the manufacturers were trying to find the ways to get maximum income with minimum expenses.

There are plenty of examples of optimization processes in pharmacology (for determination of the geometry of a molecule), in meteorology, in optimization of a trajectory of a deep-water vehicle, in optimization of power management (optimization of the production of electrical power plants), etc.

Optimization presents an important tool in decision theory and analysis of physical systems.

Optimization theory is a very developed area with its wide application in science, engineering, business management, military, and space technology.

Optimization can be defined as the process of finding the best solution to a problem in a certain sense and under certain conditions.

Along with the passage of time, optimization was evolving. Optimization became an independent area of mathematics in 1940, when Dantzig presented the so-called simplex algorithm for linear programming.

The development of nonlinear programming became great after presentation of conjugate gradient methods and quasi-Newton methods in the 1950s.

Today, there exist many modern optimization methods which are made to solve a variety of optimization problems. Now, they present the necessary tool for solving problems in diverse fields.

At the beginning, it is necessary to define an objective function, which, for example, could be a technical expense, profit or purity of materials, time, potential energy, etc.

The object function depends on certain characteristics of the system, which are known as variables. The goal is to find the values of those variables, for which the object function reaches its best value, which we call an extremum or an optimum.

It can happen that those variables are chosen in such a way that they satisfy certain conditions, i.e., restrictions.

The process of identifying the object function, variables, and restrictions for the given problem is called modeling.

The first and the most important step in an optimization process is the construction of the appropriate model, and this step can be the problem by itself. Namely, in the case that the model is too much simplified, it cannot be a faithful reflection of the practical problem. By the other side, if the constructed model is too complicated, then solving the problem is also too complicated.

After the construction of the appropriate model, it is necessary to apply the appropriate algorithm to solve the problem. It is no need to emphasize that there does not exist a universal algorithm for solving the set problem.

Sometimes, in the applications, the set of input parameters is bounded, i.e., the input parameters have values within the allowed space of input parameters $D_x$ ; we can write

$$x\in D_x.$$ (1) E1

Except (1), the next conditions can also be imposed:

$${\phi }_l\left(x_1,\dots ,x_n\right)={\phi }_{0l},l=1,\dots ,m_1<n,$$ (2) E2

$${\psi }_j\left(x_1,\dots ,x_n\right)\le {\psi }_{0j},j=1,\dots ,m_2.$$ (3) E3

Optimization task is to find the minimum (maximum) of the objective function $f\left(x\right)=f\left(x_1,\dots ,x_n\right)$ , under the conditions (1), (2), and (3).

If the object function is linear, and the functions ${\phi }_l\left(x_1,\dots ,x_n\right)l=1,\dots ,m_1$ and ${\psi }_j\left(x_1,\dots ,x_n\right)j=1,\dots ,m_2$ are linear, then it is about the linear programming problem, but if at least one of the mentioned functions is nonlinear, it is about the nonlinear programming problem.

Unconstrained optimization problem can be presented as

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

where $f\in R^n$ is a smooth function.

Problem (4) is, in fact, the unconstrained minimization problem. But, it is well known that the unconstrained minimization problem is equivalent to an unconstrained maximization problem, i.e.

$$minf\left(x\right)=-max\left(-f\left(x\right)\right),$$ (5) E5

as well as

$$maxf\left(x\right)=-min\left(-f\left(x\right)\right).$$ (6) E6

Definition 1.1.1 $x^{\ast }$ is called a global minimizer of $f$ if $f\left(x^{\ast }\right)\le f\left(x\right)$ for all $x\in R^n$ .

The ideal situation is finding a global minimizer of $f$ . Because of the fact that our knowledge of the function $f$ is usually only local, the global minimizer can be very difficult to find. We usually do not have the total knowledge about $f$ . In fact, most algorithms are able to find only a local minimizer, i.e., a point that achieves the smallest value of $f$ in its neighborhood.

So, we could be satisfied by finding the local minimizer of the function $f$ . We distinguish weak and strict (or strong) local minimizer.

Formal definitions of local weak and strict minimizer of the function $f$ are the next two definitions, respectively.

Definition 1.1.2 $x^{\ast }$ is called a weak local minimizer of $f$ if there exists a neighborhood $N$ of $x^{\ast }$ , such that $f\left(x^{\ast }\right)\le f\left(x\right)$ for all $x\in N$ .

Definition 1.1.3 $x^{\ast }$ is called a strict (strong) local minimizer of $f$ if there exists a neighborhood $N$ of $x^{\ast }$ , such that $f\left(x^{\ast }\right)<f\left(x\right)$ for all $x\in N$ .

Considering backward definitions 1.1.2 and 1.1.3, the procedure of finding local minimizer (weak or strict) does not seem such easy; it seems that we should examine all points from the neighborhood of $x^{\ast }$ , and it looks like a very difficult task.

Fortunately, if the object function $f$ satisfies some special conditions, we can solve this task in a much easier way.

For example, we can assume that the object function $f$ is smooth or, furthermore, twice continuously differentiable. Then, we concentrate to the gradient $\nabla f\left(x^{\ast }\right)$ as well as to the Hessian ${\nabla }^{2}f\left(x^{\ast }\right)$ .

All algorithms for unconstrained minimization require the user to start from a certain point, so-called the starting point, which we usually denote by $x_0$ . It is good to choose $x_0$ such that it is a reasonable estimation of the solution. But, to find such estimation, a little more knowledge about the considered set of data is needed, and the systematic investigation is needed also. So, it seems much simpler to use one of the algorithms to find $x_0$ or to take it arbitrarily.

There exist two important classes of iterative methods—line search methods and trust-region methods—made in the aim to solve the unconstrained optimization problem (4).

In this chapter, at first, we discuss different kinds of line search. Then, we consider some line search optimization methods in details, i.e., we study steepest descent method, Barzilai-Borwein gradient method, Newton method, and quasi-Newton method.

Also, we try to give some of the most recent results in these areas.

2. Line search

Now, let us consider the problem

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

where $f:R^n\to R$ is a continuously differentiable function, bounded from below.

There exists a great number of methods made in the aim to solve the problem (7).

The optimization methods based on line search utilize the next iterative scheme:

$$x_{k+1}=x_k+t_k{d}_k,$$ (8) E8

where $x_k$ is the current iterative point, $x_{k+1}$ is the next iterative point, $d_k$ is the search direction, and $t_k$ is the step size in the direction $d_k$ .

At first, we consider the monotone line search.

Now, we give the iterative scheme of this kind of search.

Algorithm 1.2.1. (Monotone line search).

Assumptions: $\mathrm{ϵ}>0$ , $x_0$ , $k≔0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP.

Step 2. Find the descent direction $d_k$ .

Step 3. Find the step size $t_k$ , such that $f\left(x_k+t_k{d}_k\right)<f\left(x_k\right)$ .

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

Step 5. Take $k≔k+1$ and go to Step 1.

Denote

$$\mathrm{\Phi }\left(t\right)=f\left(x_k+{\mathit{td}}_k\right).$$

Trying to solve the minimization problem, we are going to search for the step size $t=t_k$ , in the direction $d_k$ , such that the next relation holds:

$$\mathrm{\Phi }\left(t_k\right)<\mathrm{\Phi }\left(0\right).$$

That procedure is called the monotone line search.

We can search for the step size $t_k$ in such a way that the next relation holds:

$$f\left(x_k+t_k{d}_k\right)=\underset{t\ge 0}{min}f\left(x_k+t_k{d}_k\right),$$ (9) E9

i.e.

$$\mathrm{\Phi }\left(t_k\right)=\underset{t\ge 0}{min}\mathrm{\Phi }\left(t\right),$$ (10) E10

or we can use the next formula:

$$t_k=min\left\{t|g{\left(x_k+{\mathit{td}}_k\right)}^T{d}_k=0,t\ge 0\right\}.$$ (11) E11

In this case we are talking about the exact or the optimal line search, where the parameter $t_k$ , which is received as the solution of the one-dimensional problem (10), is the optimal step size.

By the other side, instead of using the relation (9), or the relation (11), we can be satisfied by searching for such $t_k$ , which is acceptable if the next relation suits us:

$$f\left(x_k\right)-f\left(x_k+t_k{d}_k\right)>{\delta }_k>0.$$

Then, we are talking about the inexact or the approximate or the acceptable line search, which is very much utilized in the practice.

There are several reasons to use the inexact instead of the exact line search. One of them is that the exact line search is expensive. Further, in the cases when the iteration is far from the solution, the exact line search is not efficient. Next, in the practice, the convergence rate of many optimization methods (such as Newton or quasi-Newton) does not depend on the exact line search.

First, we are going to mention so-called basic and, by the way, very well-known inexact line searches.

Algorithm 1.2.2. (Backtracking).

Assumptions: $x_k$ , the descent direction $d_k$ , $0<\delta <\frac{1}{2}$ , $\eta \in \left(0,1\right)$ .

Step 1. $t≔1$ .

Step 2. While $f\left(x_k+{\mathit{td}}_k\right)>f\left(x_k\right)+\delta {\mathit{tg}}_k^T{d}_k$ , $t≔t\cdot \eta$ .

Step 3. Set $t_k=t$ .

Now, we describe the Armijo rule.

Theorem 1.2.1. [1] Let $f\in C^1\left(R^n\right)$ and let $d_k$ be the descent direction. Then, there exists the nonnegative number $m_k$ , such that

$$f\left(x_k+{\eta }^{m_k}{d}_k\right)\le f\left(x_k\right)+c_1{\eta }^{m_k}{g}_k^T{d}_k,$$

where $c_1\in \left(0,1\right)$ and $\eta \in \left(0,1\right)$ .

Next, we describe the Goldstein rule [2].

The step size $t_k$ is chosen in such a way that

$$\begin{array}{l}f\left(x_k+{\mathit{td}}_k\right)\le f\left(x_k\right)+\delta {\mathit{tg}}_k^T{d}_k,\\ f\left(x_k+{\mathit{td}}_k\right)>f\left(x_k\right)+\left(1-\delta \right){\mathit{tg}}_k^T{d}_k,\end{array}$$

where $0<\delta <\frac{1}{2}$ .

Now, Wolfe line search rules follow [3], [4].

Standard Wolfe line search conditions are

$$f\left(x_k+t_k{d}_k\right)-f\left(x_k\right)\le \delta t_k{g}_k^T{d}_k,$$ (12) E12

$$g_{k+1}^T{d}_k\ge \sigma g_k^T{d}_k,$$ (13) E13

where $d_k$ is a descent direction and $0<\delta \le \sigma <1$ .

This efficient strategy means that we should accept a positive step length $t_k$ , if conditions (12)–(13) are satisfied.

Strong Wolfe line search conditions consist of (12) and the next, stronger version of (13):

$$\mid g_{k+1}^T{d}_k\mid \le -\sigma g_k^T{d}_k.$$ (14) E14

In the generalized Wolfe line search conditions, the absolute value in (14) is replaced by the inequalities:

$${\sigma }_1{g}_k^T{d}_k\le g_{k+1}^T{d}_k\le -{\sigma }_2{g}_k^T{d}_k,0<\delta \le {\sigma }_1<1,{\sigma }_2\ge 0.$$ (15) E15

By the other side, in the approximate Wolfe line search conditions, the inequalities (15) are changed into the next ones:

$$\sigma g_k^T{d}_k\le g_{k+1}^T{d}_k\le \left(2\delta -1\right)g_k^T{d}_k,0<\delta <\frac{1}{2},\delta <\sigma <1.$$ (16) E16

The next lemma is very important.

Lemma 1.2.1. [5] Let $f\in C\left(R^n\right)$ . Let $d_k$ be a descent direction at the point $x_k$ , and assume that the function $f$ is bounded from below along the direction $\left\{x_k+{\mathit{td}}_k|t>0\right\}$ . Then, if $0<\delta <\sigma <1$ , there exist the intervals inside which the step length satisfies standard Wolfe conditions and strong Wolfe conditions.

By the other side, the introduction of the non-monotone line search is motivated by the existence of the problems where the search direction does not have to be a descent direction. This can happen, for example, in stochastic optimization [6].

Next, some efficient quasi-Newton methods, for example, $\mathit{SR}1$ update, do not produce the descent direction in every iteration [5].

Further, some efficient methods like spectral are not monotone at all.

Some numerical results given in [7, 8, 9, 10, 11] show that non-monotone techniques are better than the monotone ones if the problem is to find the global optimal values of the object function.

Algorithms of the non-monotone line search do not insist on a descent of the object function in every step. But, even these algorithms require the reduction of the object function after a predetermined number of iterations.

The first non-monotone line search technique is presented in [12]. Namely, in [12], the problem is to find the step size which satisfies

$$f\left(x_k+t_k{d}_k\right)\le \underset{0\le j\le m\left(k\right)}{max}f\left(x_{k-j}\right)+\delta t_k{g}_k^T{d}_k,$$

where $m\left(0\right)=0$ , $0\le m\left(k\right)\le min\left\{m\left(k-1\right)+1,M\right\}$ , for $k\ge 1$ , $\delta \in \left(0,1\right)$ , where $M$ is a nonnegative integer.

This strategy is in fact the generalization of Armijo line search. In the same work, the authors suppose that the search directions satisfy the next conditions for some positive constants $b_1$ and $b_2$ :

$$\begin{array}{l}{g}_k^T{d}_k\le -b_1\parallel g_k{\parallel }^2,\\ \parallel d_k\parallel \le b_2\parallel g_k\parallel .\end{array}$$

The next non-monotone line search is described in [11].

Let $x_0$ be the starting point, and let

$$0\le {\eta }_{\mathit{min}}\le {\eta }_{\mathit{max}}\le 1,0<\delta <\sigma <1<\rho ,\mu >0.$$

Let $C_0=f\left(x_0\right)$ , $Q_0=1$ .

The step size has to satisfy the next conditions:

$$f\left(x_k+t_k{d}_k\right)\le C_k+\delta t_k{g}_k^T{d}_k,$$ (17) E17

$$g\left(x_k+t_k{d}_k\right)\ge \sigma g_k^T{d}_k.$$ (18) E18

The value ${\eta }_k$ is chosen from the interval $\left[{\eta }_{\mathit{min}},{\eta }_{\mathit{max}}\right]$ and then

$$Q_{k+1}={\eta }_k{Q}_k+1,C_{k+1}=\frac{{\eta }_k{Q}_k{C}_k+f\left(x_{k+1}\right)}{Q_{k+1}}.$$

Non-monotone rules which contain the sequence of nonnegative parameters $\left\{{\mathrm{ϵ}}_k\right\}$ are used firstly in [13], and they are successfully used in many other algorithms, for example, in [14]. The next property of the parameters ${\mathrm{ϵ}}_k$ is assumed:

$${\mathrm{ϵ}}_k>0,\sum _k{\mathrm{ϵ}}_k=\mathrm{ϵ}<\infty ,$$

and the corresponding rule is

$$f\left(x_k+t_k{d}_k\right)\le f\left(x_k\right)+c_1{t}_k{g}_k^T{d}_k+{\mathrm{ϵ}}_k.$$

Now, we give the non-monotone line search algorithm, shortly $\mathit{\text{NLSA}}$ , presented in [11].

Algorithm 1.2.3. ( $\mathit{\text{NLSA}}$ ).

Assumptions: $x_0$ , $0\le {\eta }_{\mathit{min}}\le {\eta }_{\mathit{max}}\le 1$ , $0<\delta <\sigma <1<\rho$ , $\mu >0$ .

Set $C_0=f\left(x_0\right)$ , $Q_0=1$ , $k=0$ .

Step 1. If $\parallel \nabla f\left(x_k\right)\parallel$ is sufficiently small, then STOP.

Step 2. Set $x_{k+1}=x_k+t_k{d}_k$ , where $t_k$ satisfies either the (non-monotone) Wolfe conditions (17) and (18) or the (non-monotone) Armijo conditions: $t_k={\overline{t}}_k{\rho }^{h_k}$ , where ${\overline{t}}_k>0$ is the trial step and $h_k$ is the largest integer such that (17) holds and $t_k\le \mu$ .

Step 3. Choose ${\eta }_k\in \left[{\eta }_{\mathit{min}},{\eta }_{\mathit{max}}\right]$ , and set

$$Q_{k+1}={\eta }_k{Q}_k+1,C_{k+1}=\left({\eta }_k{Q}_k{C}_k+f\left(x_{k+1}\right)\right)/Q_{k+1}.$$

Step 4. Set $k≔k+1$ and go to Step 1.

We can notice [11] that $C_{k+1}$ is a convex combination of $f\left(x_0\right),f\left(x_1\right),\dots ,f\left(x_k\right)$ . The parameter ${\eta }_k$ controls the degree of non-monotonicity.

If ${\eta }_k=0$ for all $k$ , then this non-monotone line search becomes monotone Wolfe or Armijo line search.

If ${\eta }_k=1$ for all $k$ , then $C_k=A_k$ , where

$$A_k=\frac{1}{k+1}\sum _{i=0}^{k}f\left(x_i\right).$$

Lemma 1.2.2. [11] If $\nabla f{\left(x_k\right)}^T{d}_k\le 0$ for each $k$ , then for the iterates generated by the non-monotone line search algorithm, we have $f_k\le C_k\le A_k$ for each $k$ . Moreover, if $\nabla f{\left(x_k\right)}^T{d}_k<0$ and $f\left(x\right)$ are bounded from below, then there exists $t_k$ satisfying either Wolfe or Armijo conditions of the line search update.

This study would be very incomplete unless we mention that there are many modifications of the abovementioned line searches. All these modifications are made to improve the previous results.

For example, in [15], the new inexact line search is described by the next way.

Let $\beta \in \left(0,1\right)$ , $\sigma \in \left(0,\frac{1}{2}\right)$ ; let $B_k$ be a symmetric positive definite matrix which approximates ${\nabla }^{2}f\left(x_k\right)$ and $s_k=-\frac{g_k^T{d}_k}{d_k^T{B}_k{d}_k}$ . The step size $t_k$ is the largest one in $\left\{s_k,s_k\beta ,s_k{\beta }^2,\dots \right\}$ such that

$$f\left(x_k+{\mathit{td}}_k\right)-f\left(x_k\right)\le \mathit{\sigma t}\left[g_k^T{d}_k+\frac{1}{2}{\mathit{td}}_k^T{B}_k{d}_k\right].$$

Further, in [16], a new inexact line search rule is presented. This rule is a modified version of the classical Armijo line search rule. We describe it now.

Let $g=\nabla f\left(x\right)$ be a Lipschitz continuous function and $L$ the Lipschitz constant. Let $L_k$ be an approximation of $L$ . Set

$${\beta }_k=-\frac{g_k^T{d}_k}{L_k\parallel d_k{\parallel }^2}.$$

Find a step size $t_k$ as the largest component in the set $\left\{{\beta }_k,{\beta }_k\rho ,{\beta }_k{\rho }^2\dots \right\}$ such that the inequality

$$f\left(x_k+t_k{d}_k\right)\le f\left(x_k\right)+\sigma t_k\left(g_k^T{d}_k-\frac{1}{2}{t}_k\mu L_k\parallel d_k{\parallel }^2\right)$$

holds, where $\sigma \in \left(0,1\right)$ , $\mu \in \left[0,\infty \right)$ , and $\rho \in \left(0,1\right)$ are given constants.

Next, in [17], a new, modified Wolfe line search is given in the next way.

Find $t_k>0$ such that

$$\begin{array}{l}f\left(x_k+t_k{d}_k\right)-f\left(x_k\right)\le min\left\{\delta t_k{g}_k^T{d}_k,-\gamma t_k^2\parallel d_k{\parallel }^2\right\},\\ g{\left(x_k+t_k{d}_k\right)}^T{d}_k\ge \sigma g_k^T{d}_k,\end{array}$$

where $\delta \in \left(0,1\right)$ , $\sigma \in \left(\delta ,1\right)$ , and $\gamma >0$ .

More recent results on this topic can be found, for example, in [18, 19, 20, 21, 22, 23].

2.1 Steepest descent ( $\mathit{SD}$ )

The classical steepest descent method which is designed by Cauchy [24] can be considered as one among the most important procedures for minimization of real-valued function defined on ${\mathbb{R}}^n$ .

Steepest descent is one of the simplest minimization methods for unconstrained optimization. Since it uses the negative gradient as its search direction, it is known also as the gradient method.

It has low computational cost and low matrix storage requirement, because it does not need the computations of the second derivatives to be solved to calculate the search direction [25].

Suppose that $f\left(x\right)$ is continuously differentiable in a certain neighborhood of a point $x_k$ and also suppose that $g_k\triangleq \nabla f\left(x_k\right)\ne 0$ .

Using Taylor expansion of the function $f$ near $x_k$ as well as Cauchy-Schwartz inequality, one can easily prove that the greatest fall of $f$ exists if and only if $d_k=-g_k$ , i.e., $-g_k$ is the steepest descent direction.

The iterative scheme of the $\mathit{SD}$ method is

$$x_{k+1}=x_k-t_k{g}_k.$$ (19) E19

The classical steepest descent method uses the exact line search.

Now, we give the algorithm of the steepest descent method which refers to the exact as well as to the inexact line search.

Algorithm 1.2.4. (Steepest descent method, i.e., $\mathit{SD}$ method).

Assumptions: $0<\mathrm{ϵ}\ll 1$ , $x_0\in {\mathbb{R}}^n$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \epsilon$ , then STOP, else set $d_k=-g_k$ .

Step 2. Find the step size $t_k$ , which is the solution of the problem

$$\underset{t\ge 0}{min}f\left(x_k+{\mathit{td}}_k\right),$$ (20) E20

else find the step size $t_k$ by any of the inexact line search methods.

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

Step 4. Set $k≔k+1$ and go to Step 1.

The classical and the oldest steepest descent step size $t_k$ , which was designed by Cauchy (in the case of the exact line search), is computed as [26]

$$t_k=\frac{g_k^T{g}_k}{g_k^T{\mathit{Gg}}_k},$$

where $g_k=\nabla f\left(x_k\right)$ and $G={\nabla }^{2}f\left(x_k\right)$ .

Theorem 1.2.2. [27] (Global convergence theorem of the $\mathit{SD}$ method) Let $f\in C^1$ . Then, each accumulation point of the iterative sequence $\left\{x_k\right\}$ , generated by Algorithm 1.2.4, is a stationary point.

Remark 1.2.1. The steepest descent method has at least the linear convergence rate.

More information about the convergence of the SD method can be found in [5, 27].

Although known as the first unconstrained optimization method, this method is still a theme considered by scientists.

Different modifications of this method are made, for example, see [25, 28, 29, 30, 31, 32].

In [28], the authors presented a new search direction from Cauchy’s method in the form of two parameters known as Zubai’ah-Mustafa-Rivaie-Ismail method, shortly, $\mathit{\text{ZMRI}}$ method:

$$d_k=-g_k-\parallel g_k\parallel g_{k-1}.$$ (21) E21

So, in [28], a new modification of $\mathit{SD}$ method is suggested using a new search direction, $d_k$ , given by (21). The numerical results are presented based on the number of iterations and CPU time. It is shown that this new method is efficient when it is compared to the classical $\mathit{SD}$ .

In [25], a new scaled search direction of $\mathit{SD}$ method is presented. The inspiration for this new method is the work of Andrei [33], in which the author presents and analyzes a new scaled conjugate gradient algorithm, based on an interpretation of the secant equation and on the inexact Wolfe line search conditions.

The method proposed in [25] is known as Rashidah-Rivaie-Mamat ( $\mathit{RRM}$ ) method, and it suggests the direction $d_k$ given by the next relation:

$$d_k=\left\{\begin{array}{l}-g_k,\text{if}k=0,\\ -{\theta }_k{g}_k-\parallel g_k\parallel g_{k-1},\end{array}\right.$$ (22) E22

where ${\theta }_k$ is a scaling parameter, ${\theta }_k=\frac{d_{k-1}^T{y}_{k-1}}{\parallel g_{k-1}{\parallel }^2}$ , $y_{k-1}=g_k-g_{k-1}$ .

Further, in [25], a comparison among $\mathit{RRM}$ , $\mathit{\text{ZMRI}}$ , and $\mathit{SD}$ methods is made; it is shown that $\mathit{RRM}$ method is better than $\mathit{\text{ZMRI}}$ and $\mathit{SD}$ methods.

It is interesting that the exact line search is used in [25].

In [34], the properties of steepest descent method from the literature are reviewed together with advantages and disadvantages of each step size procedure.

Namely, the step size procedures, which are compared in this paper, are:

1. $t_k=\frac{g_k^T{g}_k}{g_k^T{H}_k{g}_k}$ : Step size method by Cauchy [24], computed by exact line search ( $C$ step size).

2. Given $s>0,\beta ,\sigma \in \left(0,1\right),t_k=max\left\{s,\mathit{s\beta },s{\beta }^2,\dots \right\}$ such that

$$f\left(x_k+t_k{d}_k\right)\le f\left(x_k\right)+\sigma t_k{g}_k^T{d}_k-{\text{Armijo}}^{’}\mathrm{s}\text{line} \text{search}\left(\mathrm{A}\text{step} \text{size}\right).$$

3. Given $\beta ,\sigma \in \left(0,1\right)$ , ${\tilde{t}}_0=1$ , and $t_k=\beta {\tilde{t}}_k$ such that

$$f\left(x_k+t_k{d}_k\right)\le f\left(x_k\right)+\sigma t_k{g}_k^T{d}_k-\text{Backtracking} \text{line} \text{search}\left(\mathrm{B}\text{step} \text{size}\right).$$

4. $t_k=\frac{s_{k-1}^T{y}_{k-1}}{\parallel y_{k-1}{\parallel }^2}$ , ( $\mathit{BB}1$ ), $t_k=\frac{\parallel s_{k-1}{\parallel }^2}{s_{k-1}^T{y}_{k-1}}$ , ( $\mathit{BB}2$ ), $s_{k-1}=x_k-x_{k-1}$ $y_{k-1}=g_k-g_{k-1},$ : Barzilai and Borwein’s formula. The convergence is R-superlinear.

5. $t_k=\frac{t_{k-1}^2{g}_k^T{g}_k}{2(f\left(x_k+t_k{d}_k\right)-f\left(x_k\right)+t_{k-1}{g}_k^T{g}_k},$ : Elimination line search ( $\mathit{EL}$ step size), which estimates the step size without computation of the Hessian.

The comparison is based on time execution, number of total iteration, total percentage of function, gradient and Hessian evaluation, and the most decreased value of objective function obtained.

From the numerical results, the authors conclude that the $A$ method and $\mathit{BB}1$ method are the best methods among others.

Further, in [34], the general conclusions about the steepest descent method are given:

1. This method is sensitive to the initial point.

2. This method has a descent property, and it is a logical starting procedure for all gradient based methods.

3. $x_k$ approaches the minimizer slowly, in fact in a zigzag way.

In [35], in the aim to achieve fast convergence and the monotone property, a new step size for the steepest descent method is suggested.

In [36], for quadratic positive definite problems, an over-relaxation has been considered. Namely, Raydan and Svaiter [36] proved that the poor behavior of the steepest descent method is due to the optimal Cauchy choice of step size and not to the choice of the search direction. These results are extended in [29] to convex, well-conditioned functions. Further, in [29], it is shown that a simple modification of the step length by means of a random variable uniformly distributed in $\left(0,1\right]$ , for the strongly convex functions, represents an improvement of the classical gradient descent algorithm. Namely, in this paper, the idea is to modify the gradient descent method by introducing a relaxation of the following form:

$$x_{k+1}=x_k+{\theta }_k{t}_k{d}_k,$$ (23) E23

where ${\theta }_k$ is the relaxation parameter, a random variable uniformly distributed between $0$ and $1$ .

In the recent years, the steepest descent method has been applied in many branches of science; one can be inspired, for example, by [37, 38, 39, 40, 41, 42, 43].

2.2 Barzilai and Borwein gradient method

Remind to the fact that $\mathit{SD}$ method performs poorly, converges linearly, and is badly affected by the ill-conditioning.

Also, remind to the fact that this poor behavior of $\mathit{SD}$ method is due to the optimal choice of the step size and not to the choice of the steepest descent direction $-g_k$ .

Barzilai and Borwein presented [44] a two-point step size gradient method, which is well known as $\mathit{BB}$ method.

The step size is derived from a two-point approximation to the secant equation.

Consider the gradient iteration form:

$$x_{k+1}=x_k-t_k{g}_k.$$

It can be rewritten as $x_{k+1}=x_k-D_k{g}_k,$ where $D_k=t_{k}I$ .

To make the matrix $D_k$ having quasi-Newton property, the step size $t_k$ is computed in such a way that we get

$$min\parallel s_{k-1}-D_k{y}_{k-1}\parallel .$$

This yields that

$$t_k^{\mathit{BB}1}=\frac{s_{k-1}^T{y}_{k-1}}{y_{k-1}^T{y}_{k-1}},s_{k-1}=x_k-x_{k-1},y_{k-1}=g_k-g_{k-1}.$$ (24) E24

But, using symmetry, we may minimize $\parallel D_k^{-1}{s}_{k-1}-y_{k-1}\parallel$ , with respect to $t_k$ , and we get:

$$t_k^{\mathit{BB}2}=\frac{\parallel s_{k-1}{\parallel }^2}{s_{k-1}^T{y}_{k-1}},s_{k-1}=x_k-x_{k-1},y_{k-1}=g_k-g_{k-1}.$$ (25) E25

Now, we give the algorithm of $\mathit{BB}$ method.

Algorithm 1.2.5. (Barzilai-Borwein gradient method, i.e., $\mathit{BB}$ method).

Assumptions: $0<\mathrm{ϵ}\ll 1$ , $x_0\in {\mathbb{R}}^n$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP, else set $d_k=-g_k$ .

Step 2. If $k=0$ , then find the step size $t_0$ by the line search, else compute $t_k$ using the formula (24) or (25).

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

Step 4. Set $k≔k+1$ and go to Step 1.

Considering Algorithm 1.2.5, we can conclude that this method does not require any matrix computation or any line search.

The Barzilai-Borwein method is in fact the gradient method, which requires less computational work than $\mathit{SD}$ method, and it speeds up the convergence of the gradient method. Barzilai and Borwein proved that $\mathit{BB}$ algorithm is $R-$ superlinearly convergent for the quadratic case.

In the general non-quadratic case, a globalization strategy based on non-monotone line search is applied in this method.

In this general case, $t_k$ , computed by (24) or (25), may be unacceptably large or small. That is the reason why we assume that there exist the numbers $t^l$ and $t^r$ , such that

$$0<t^l\le t_k\le t^r,\text{for}\mathrm{all}k.$$

Using the iteration

$$x_{k+1}=x_k-\frac{1}{t_k}{g}_k=x_k-{\lambda }_k{g}_k,$$ (26) E26

with

$$t_k=\frac{s_{k-1}^T{y}_{k-1}}{s_{k-1}^T{s}_{k-1}},{\lambda }_k=\frac{1}{t_k},$$

$$s_k=-\frac{1}{t_k}{g}_k=-{\lambda }_k{g}_k,$$

we get

$$t_{k+1}=\frac{s_k^T{y}_k}{s_k^T{s}_k}=\frac{-{\lambda }_k{g}_k^T{y}_k}{{\lambda }_k^2{g}_k^T{g}_k}=-\frac{g_k^T{y}_k}{{\lambda }_k{g}_k^T{g}_k}.$$

Now, we give the algorithm of the Barzilai-Borwein method with non-monotone line search.

Algorithm 1.2.6. ( $\mathit{BB}$ method with non-monotone line search).

Assumptions: $0<\mathrm{ϵ}\ll 1$ , $x_0\in {\mathbb{R}}^n$ , $M\ge 0$ is an integer, $\rho \in \left(0,1\right)$ , $\delta >0$ , $0<{\sigma }_1<{\sigma }_2<1$ , $t^l$ , $t^r$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP.

Step 2. If $t_k\le t^l$ , or $t_k\ge t^r$ , then set $t_k=\delta$ .

Step 3. Set $\lambda =\frac{1}{t_k}$ .

Step 4. (non-monotone line search) If

$$f\left(x_k-\lambda g_k\right)\le \underset{0\le j\le min\left(k,M\right)}{max}f\left(x_{k-j}\right)-\mathit{\rho \lambda }{g}_k^T{g}_k,$$

then set

$${\lambda }_k=\lambda ,x_{k+1}=x_k-{\lambda }_k{g}_k,$$

and go to Step 6.

Step 5. Choose $\sigma \in \left[{\sigma }_1,{\sigma }_2\right]$ , set $\lambda =\mathit{\sigma \lambda }$ , and go to Step 4.

Step 6. Set $t_{k+1}=-\frac{g_k^T{y}_k}{{\lambda }_k{g}_k^T{g}_k}$ and $k≔k+1$ , and return to Step 1.

Obviously, the above algorithm is globally convergent.

Several authors paid attention to the Barzilai-Borwein method, and they proposed some variants of this method.

In [8], the globally convergent Barzilai-Borwein method is proposed by using non-monotone line search by Grippo et al. [12]. In the same paper, Raydan proves the global convergence of the non-monotone Barzilai-Borwein method.

Further, Grippo and Sciandrone [45] propose another type of the non-monotone Barzilai-Borwein method.

Dai [7] gives the basic analysis of the non-monotone line search strategy.

Moreover, in [46] numerical results are presented, using

$$t_k=\frac{s_{\nu \left(k\right)}^T{y}_{\nu \left(k\right)}}{s_{\nu \left(k\right)}^T{s}_{\nu \left(k\right)}}.$$ (27) E27

and

$$\nu \left(k\right)=M_c\cdot ⌞\frac{k-1}{M_c}⌟,$$

where for $r\in \mathbb{R}$ , $⌞r⌟$ denotes the largest integer $j$ such that $j\le r$ and Mc is a positive integer. The gradient method with (27) is called the cyclic Barzilai-Borwein method. Numerical results in [46] prove that their method performs better than the Barzilai-Borwein method.

Many researchers study the gradient method for minimizing a strictly convex quadratic function, namely,

$$minf\left(x\right)=\frac{1}{2}{x}^T\mathit{Ax}-b^{T}x,$$ (28) E28

where $A\in {\mathbb{R}}^{n\times n}$ is a symmetric positive definite matrix and $b\in {\mathbb{R}}^n$ is a given vector. For an application of the Barzilai-Borwein method to the problem (28), Raydan [47] establishes global convergence, and Dai and Liao [48] prove $\mathbb{R}$ -linear rate of convergence. Friedlander, Martinez, Molina, and Raydan [49] propose a new gradient method with retards, in which $t_k$ is defined by

$$t_k=\frac{g_{\nu \left(k\right)}^T{A}^{\rho \left(k\right)+1}{g}_{\nu \left(k\right)}}{g_{\nu \left(k\right)}^T{A}^{\rho \left(k\right)}{g}_{\nu \left(k\right)}},\nu \left(k\right)\in \left\{k,k-1,\dots ,max\left\{0,k-m\right\}\right\}$$ (29) E29

and $\rho \left(k\right)\in \left\{q_1,\dots ,q_m\right\}$ , where $m$ is a positive integer and $q_1,\dots ,q_m\ge -2$ are integers. In the same paper, they establish its global convergence for problem (28) and prove the $Q$ -superlinear rate of convergence in the special case.

In [50], the authors extend the Barzilai-Borwein method, and they give extended Barzilai-Borwein method, which they denote $\mathit{EBB}$ . They also establish global and $Q-$ superlinear convergence properties of the proposed method for minimizing a strictly convex quadratic function. Furthermore, they discuss an application of their method to general objective functions. In [50], a new step size is proposed by extending (29). Namely, in this paper, following Friedlander et al. [49], a new step size is proposed as follows:

$$\begin{array}{l}{t}_k=\sum _{i=1}^l{\varphi }_i\frac{g_{{\nu }_{i\left(k\right)}}^T{A}^{{\rho }_i\left(k\right)+1}{g}_{{\nu }_i\left(k\right)}}{g_{{\nu }_{i\left(k\right)}}^T{A}^{{\rho }_i\left(k\right)}{g}_{{\nu }_i\left(k\right)}},\\ {\varphi }_i\ge 0,\sum _{i=1}^n{\varphi }_i=1,\\ {\nu }_i\left(k\right)\in \left\{k,k-1,\dots ,max\left\{0,k-m\right\}\right\}\end{array}$$

and

$${\varphi }_i\left(k\right)\in \left\{q_1,\dots ,q_m\right\},$$

where $l$ and $m$ are positive integers and $q_1,\dots ,q_m$ are integers.

Also, an application of algorithm $\mathit{EBB}$ to general unconstrained minimization problems (4) is considered.

Following Raydan [8], the authors [50] further combine the non-monotone line search and algorithm $\mathit{EBB}$ to get the algorithm called $\mathit{\text{NEBB}}$ . They also prove the global convergence of the algorithm $\mathit{\text{NEBB}}$ , under some classical assumptions.

The Barzilai-Borwein method and its related methods are reviewed by Dai and Yuan [51] and Fletcher [52].

In [53], a new concept of the approximate optimal step size for gradient method is introduced and used to interpret the $\mathit{BB}$ method; an efficient gradient method with the approximate optimal step size for unconstrained optimization is presented. The next definition is introduced in [53].

Definition 1.2.1. Let $\mathrm{\Phi }\left(t\right)$ be an approximation model of $f\left(x_k-{\mathit{tg}}_k\right)$ . A positive constant $t^{\ast }$ is called approximate optimal step size associated to $\mathrm{\Phi }\left(t\right)$ for gradient method, if $t^{\ast }$ satisfies

$$t^{\ast }=\mathit{arg}\underset{t>0}{min}\mathrm{\Phi }\left(t\right).$$

The approximate optimal step size is different from the steepest descent step size, which will lead to the expensive computational cost. The approximate optimal step size is generally calculated easily, and it can be applied to unconstrained optimization.

Due to the effectiveness of $t_k^{\mathit{BB}1}$ and the fact that $t_k^{\mathit{BB}1}=\mathit{arg}{min}_{t>0}\mathrm{\Phi }\left(t\right),$ we can naturally ask if more suitable approximation models can be constructed to generate more efficient approximate optimal step-sizes.

This is the purpose of work [53]. Further, if the objective function $f\left(x\right)$ is not close to a quadratic function on the line segment between $x_{k-1}$ and $x_k$ , in this paper a conic model is developed to generate the approximate optimal step size if the conic model is suitable to be used. Otherwise, the authors consider two cases:

1. If $s_{k-1}^T{y}_{k-1}>0$ , the authors construct a new quadratic model, to derive the approximate optimal step size.

2. If $s_{k-1}^T{y}_{k-1}\le 0$ , they construct a new quadratic model or two other new approximation models to generate the approximate optimal step size for gradient method. They also analyze the convergence of the proposed method under some suitable conditions. Numerical results show the proposed method is better than the BB method.

In [54], derivative-free iterative scheme that uses the residual vector as search direction for solving large-scale systems of nonlinear monotone equations is presented.

The Barzilai-Borwein method is widely used; some interesting results can be found in [55, 56, 57].

2.3 Newton method

The basic idea of Newton method for unconstrained optimization is the iterative usage of the quadratic approximation $q^{\left(k\right)}$ to the objective function $f$ at the current iterate $x_k$ and then minimization of such approximation $q^{\left(k\right)}$ .

Let $f:{\mathbb{R}}^n\to \mathbb{R}$ be twice continuously differentiable, $x_k\in {\mathbb{R}}^n$ , and let the Hessian ${\nabla }^{2}f\left(x_k\right)$ be positive definite.

We model $f$ at the current point $x_k$ by the quadratic approximation $q^{\left(k\right)}$ :

$$f\left(x_k+s\right)\approx q^{\left(k\right)}\left(s\right)=f\left(x_k\right)+\nabla f{\left(x_k\right)}^{T}s+\frac{1}{2}{s}^T{\nabla }^{2}f\left(x_k\right)s,s=x-x_k.$$

Minimization of $q^{\left(k\right)}\left(s\right)$ gives the next iterative scheme:

$$x_{k+1}=x_k-{\left({\nabla }^{2}f\left(x_k\right)\right)}^{-1}\nabla f\left(x_k\right),$$

which is known as Newton formula.

Denote $G_k={\nabla }^{2}f\left(x_k\right)$ , $g_k=\nabla f\left(x_k\right)$ .

Then, we have a simpler form:

$$x_{k+1}=x_k-G_k^{-1}{g}_k.$$ (30) E30

A Newton direction is

$$s_k=x_{k+1}-x_k=-G_k^{-1}{g}_k.$$ (31) E31

We have supposed that $G_k$ is positive definite. So, the Newton direction is a descent direction. This we can conclude from

$$g_k^T{s}_k=-g_k^T{G}_k^{-1}{g}_k<0.$$

Now, we give the algorithm of the Newton method.

Algorithm 1.2.7. (Newton method).

Assumptions: $\mathrm{ϵ}>0$ , $x_0\in {\mathbb{R}}^n$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP.

Step 2. Solve $G_{k}s=-g_k$ for $s_k$ .

Step 3. Set $x_{k+1}=x_k+s_k$ .

Step 4. $k≔k+1$ , return to Step 1.

The next theorem shows the local convergence and the quadratic convergence rate of Newton method.

Theorem 1.2.3. [27] (Convergence theorem of Newton method) Let $f\in C^2$ and $x_k$ be close enough to the solution $x^{\ast }$ of the minimization problem with $g\left(x^{\ast }\right)=0$ . If the Hessian $G\left(x^{\ast }\right)$ is positively definite and $G\left(x\right)$ satisfies Lipschitz condition

$$\mid G_{\mathit{ij}}\left(x\right)-G_{\mathit{ij}}\left(y\right)\mid \le \beta \parallel x-y\parallel ,\mathit{for\hspace{0.17em}some}\beta ,\mathit{\text{for}}\mathit{all}i,j,$$

where $G_{\mathit{ij}}\left(x\right)$ is the $\left(i,j\right)$ element of $G\left(x\right)$ and then for all $k$ , Newton direction (31) is well-defined; the generated sequence $\left\{x_k\right\}$ converges to $x^{\ast }$ with a quadratic rate.

But, in spite of this quadratic rate, the Newton method is a local method: when the starting point is far away from the solution, there is a possibility that $G_k$ is not positive definite, as well as Newton direction is not a descent direction.

So, to guarantee the global convergence, we can use Newton method with line search. We can remind to the fact that only when the step size sequence $\left\{t_k\right\}$ tends to 1, Newton method is convergent with the quadratic rate.

Newton iteration with line search is as follows:

$$d_k=-G_k^{-1}{g}_k,$$ (32) E32

$$x_{k+1}=x_k+t_k{d}_k.$$ (33) E33

Now, we give the algorithm.

Algorithm 1.2.8. (Newton method with line search).

Assumptions: $\mathrm{ϵ}>0$ , $x_0\in {\mathbb{R}}^n$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP.

Step 2. Solve $G_{k}d=-g_k$ for $d_k$ .

Step 3. Line search step: find $t_k$ such that

$$f\left(x_k+t_k{d}_k\right)=\underset{t\ge 0}{min}f\left(x_k+{\mathit{td}}_k\right),$$

or find $t_k$ such that (inexact) Wolfe line search rules hold.

Step 4. Set $x_{k+1}=x_k+t_k{d}_k$ and $k=k+1$ , and go to Step 1.

The next theorems claim that Algorithm 1.2.8 with the exact line search, as well as Algorithm 1.2.8 with the inexact line search, are globally convergent.

Theorem 1.2.4. [27] Let $f:{\mathbb{R}}^n\to \mathbb{R}$ be twice continuously differentiable on open convex set $D\subset {\mathbb{R}}^n$ . Assume that for any $x_0\in D$ there exists a constant $m>0$ , such that $f\left(x\right)$ satisfies

$$u^T{\nabla }^{2}f\left(x\right)u\ge m\parallel u{\parallel }^2,\mathit{\text{for}}\mathit{all}u\in {\mathbb{R}}^n,x\in L\left(x_0\right),$$ (34) E34

where $L\left(x_0\right)=\left\{x|f\left(x\right)\le f\left(x_0\right)\right\}$ is the corresponding level set. Then, the sequence $\left\{x_k\right\}$ , generated by Algorithm 1.2.8, with the exact line search, satisfies:

1. When $\left\{x_k\right\}$ is a finite sequence, $g_k=0$ for some $k$ .

2. When $\left\{x_k\right\}$ is an infinite sequence, $\left\{x_k\right\}$ converges to the unique minimizer $x^{\ast }$ of $f$ .

Note that the next relation holds from the standard Wolfe line search:

$$f\left(x_k\right)-f\left(x_k+t_k{d}_k\right)\ge \overline{\eta }\parallel g_k{\parallel }^2{cos}^2\angle \left(d_k,-g_k\right),$$ (35) E35

where the constant $\overline{\eta }$ does not depend on $k$ .

Theorem 1.2.5. [27] Let $f:{\mathbb{R}}^n\to \mathbb{R}$ be twice continuously differentiable on open convex set $D\subset {\mathbb{R}}^n$ . Assume that for any $x_0\in D$ there exists a constant $m>0$ , such that $f\left(x\right)$ satisfies the relation (34) on the level set $L\left(x_0\right)$ . If the line search satisfies the relation (35), then the sequence $\left\{x_k\right\}$ , generated by Algorithm 1.2.8, with the inexact Wolfe line search, satisfies

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

and $\left\{x_k\right\}$ converges to the unique minimizer of $f\left(x\right)$ .

2.4 Modified Newton method

The main problem in Newton method could be the fact that the Hessian $G_k$ may be not positive definite. In that case, we are not sure that the objective function $f$ has its minimizers; furthermore, when $G_k$ is indefinite, the objective function $f$ is unbounded.

So, many modified schemes are made. Now, we describe the next two methods shortly.

In [58], Goldstein and Price use the steepest descent method when $G_k$ is not positive definite. Denoting the angle between $d_k$ and $-g_k$ by $\theta$ , as well as having in view the angle rule, $\theta \le \frac{\pi }{2}-\mu$ , where $\mu >0$ , they determine the direction $d_k$ as

$$d_k=\left\{\begin{array}{l}-G_k^{-1}{g}_k,\text{if}cos\theta \ge \eta ,\\ -g_k,\text{otherwise},\end{array}\right.$$

where $\eta >0$ is a given constant.

In [59], the authors present another modified Newton method. When $G_k$ is not positive definite, Hessian $G_k$ is changed into $G_k+{\nu }_{k}I$ , where ${\nu }_k>0$ is chosen in such a way that $G_k+{\nu }_{k}I$ is positive definite and well-conditioned. Otherwise, when $G_k$ is positive definite, ${\nu }_k=0$ .

To consider the other modified Newton methods, such as finite difference Newton method, negative curvature direction method, Gill-Murray stable Newton method, etc., one can see [27], for example.

2.5 Inexact Newton method

By the other side, because of the high cost of the exact Newton method, especially when the dimension $n$ is large, the inexact Newton method might be a good solution. This type of method means that we only approximately solve the Newton equation.

Consider solving the nonlinear equations:

$$F\left(x\right)=0,$$ (36) E36

where $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is assumed to have the next properties:

A₁ There exists $x^{\ast }$ such that $F\left(x^{\ast }\right)=0$ .

A₂ $F$ is continuously differentiable in the neighborhood of $x^{\ast }$ .

A₃ $F^{\prime }\left(x^{\ast }\right)$ is nonsingular.

Remind that the basic Newton step is obtained by solving

$$F^{\prime }\left(x_k\right)s_k=-F\left(x_k\right)$$

and setting

$$x_{k+1}=x_k+s_k.$$

The inexact Newton method means that we solve

$$F^{\prime }\left(x_k\right)s_k=-F\left(x_k\right)+r_k,$$ (37) E37

where

$$\parallel r_k\parallel \le {\eta }_k\parallel F\left(x_k\right)\parallel .$$ (38) E38

Set

$$x_{k+1}=x_k+s_k.$$ (39) E39

Here, $r_k$ denotes the residual, and the sequence $\left\{{\eta }_k\right\}$ , where $0<{\eta }_k<1$ , is the sequence which controls the inexactness.

Now, we give two theorems; the first of them claims the linear convergence, and the second claims the superlinear convergence of the inexact Newton method.

Theorem 1.2.6. [27] Let $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ satisfy the assumptions $A_1$ – $A_3$ . Let the sequence $\left\{{\eta }_k\right\}$ satisfies $0\le {\eta }_k\le \eta <t<1$ . Then, for some $\mathrm{ϵ}>0$ , if the starting point is sufficiently near $x^{\ast }$ , the sequence $\left\{x_k\right\}$ generated by inexact Newton’s method (37)–(39) converges to $x^{\ast }$ , and the convergence rate is linear, i.e.

$$\parallel x_{k+1}-x^{\ast }{\parallel }_{\ast }\le t\parallel x_k-x^{\ast }{\parallel }_{\ast },$$

where $\parallel y{\parallel }_{\ast }=\parallel F^{\prime }\left(x^{\ast }\right)y\parallel$ .

Theorem 1.2.7. [27] Let all assumptions of Theorem 1.2.6 hold. Assume that the sequence $\left\{x_k\right\}$ , generated by the inexact Newton method, converges to $x^{\ast }$ . then

$$\parallel r_k\parallel =o\left(\parallel F\left(x_k\right)\parallel \right),k\to \infty ,$$

if and only if $\left\{x_k\right\}$ converges to $x^{\ast }$ superlinearly.

The relation

$$x_{k+1}=x_k-\frac{f^{\prime }\left(x_k\right)}{f^{\prime }\left(x_k\right)-f^{\prime }\left(x_{k-1}\right)}\cdot \left(x_k-x_{k-1}\right),$$ (40) E40

presents the secant method.

In [60], a modification of the classical secant method for solving nonlinear, univariate, and unconstrained optimization problems based on the development of the cubic approximation is presented. The iteration formula including an approximation of the third derivative of $f\left(x\right)$ by using the Taylor series expansion is derived. The basic assumption on the objective function $f\left(x\right)$ is that $f\left(x\right)$ is a real-valued function of a single, real variable $x$ and that $f\left(x\right)$ has a minimum at $x^{\ast }$ . Furthermore, in this chapter it is noted that the secant method is the simplification of Newton method. But, the order of the secant method is lower than one of the Newton methods; it is $Q$ -superlinearly convergent, and its order is $p=\frac{\sqrt{5}+1}{2}\approx \mathrm{1,618}$ .

This modified secant method is constructed in [60], having in view, as it is emphasized, that it is possible to construct a cubic function which agrees with $f\left(x\right)$ up to the third derivatives. The third derivative of the objective function $f$ is approximated as

$$f^{‴}\left(x\right)=\frac{3\left\{\frac{2\left[f^{\prime }\left(x_k\right)-\frac{f\left(x_k\right)-f\left(x_{k-1}\right)}{x_k-x_{k-1}}\right]}{x_k-x_{k-1}}-f^{″}\left(x_k\right)\right\}}{x_{k-1}-x_k}.$$

In [61], the authors propose an inexact Newton-like conditional gradient method for solving constrained systems of nonlinear equations. The local convergence of the new method as well as results on its rate is established by using a general majorant condition.

2.6 Quasi-Newton method

Consider the Newton method.

For various practical problems, the computation of Hessian may be very expensive, or difficult, or Hessian can be unavailable analytically. So, the class of so-called quasi-Newton methods is formed, such that it uses only the objective function values and the gradients of the objective function and it is close to Newton method. Quasi-Newton method is such a class of methods which does not compute Hessian, but it generates a sequence of Hessian approximations and maintains a fast rate of convergence.

So, we would like to construct Hessian approximation $B_k$ in quasi-Newton method. Naturally, it is desirable that the sequence $\left\{B_k\right\}$ possesses positive definiteness, as well as its direction $d_k=-B_k^{-1}{g}_k$ should be a descent one.

Now, let  $f:{\mathbb{R}}^n\to \mathbb{R}$ be twice continuously differentiable function on an open set $D\subset {\mathbb{R}}^n$ . Consider the quadratic approximation of $f$ at $x_{k+1}$ :

$$f\left(x\right)\approx f\left(x_{k+1}\right)+g_{k+1}^T\left(x-x_{k+1}\right)+\frac{1}{2}{\left(x-x_{k+1}\right)}^T{G}_{k+1}\left(x-x_{k+1}\right).$$

Finding the derivatives, we get

$$g\left(x\right)\approx g_{k+1}+G_{k+1}\left(x-x_{k+1}\right).$$

Setting $x=x_k$ and using the standard notation: $s_k=x_{k+1}-x_k$ , $y_k=g_{k+1}-g_k$ , from the last relation, we get

$$G_{k+1}^{-1}{y}_k\approx s_k.$$ (41) E41

Relation (41) transforms into the next one if $f$ is the quadratic function:

$$G_{k+1}^{-1}{y}_k=s_k.$$ (42) E42

Let $H_k$ be the approximation of the inverse of Hessian. Then, we want $H_k$ to satisfy the relation (42). In this way, we come to the quasi-Newton condition or quasi-Newton equation:

$$H_{k+1}{y}_k=s_k.$$ (43) E43

Let $B_{k+1}=H_{k+1}^{-1}$ be the approximation of Hessian $G_{k+1}$ . Then

$$B_{k+1}{s}_k=y_k$$ (44) E44

is also the quasi-Newton equation.

If

$$s_k^T{y}_k>0,$$ (45) E45

then the matrix $B_{k+1}$ is positive definite. The condition (45) is known as the curvature condition.

Algorithm 1.2.9. (A general quasi-Newton method).

Assumptions: $0\le \mathrm{ϵ}<1$ , $x_0\in {\mathbb{R}}^n$ , $H_0\in {\mathbb{R}}^{n\times n}$ . Let $k=0$ .

Step 1. If $\parallel g_k\parallel \le \mathrm{ϵ}$ , then STOP.

Step 2. Compute $d_k=-H_k{g}_k$ .

Step 3. Find $t_k$ by line search and set $x_{k+1}=x_k+t_k{d}_k$ .

Step 4. Update $H_k$ into $H_{k+1}$ such that quasi-Newton equation (43) holds.

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

In Algorithm 1.2.9, usually we take $H_0=I$ , where $I$ is an identity matrix.

Sometimes, instead of $H_k$ , we use $B_k$ in Algorithm 1.2.9.

Then, Step 2 becomes

Step $2^{\ast }$ . Solve

$$B_{k}d=-g_k,\text{for}{d}_k.$$

By the other side, Step 4 becomes

Step $4^{\ast }$ . Update $B_k$ into $B_{k+1}$ in such a way that quasi-Newton equation (44) holds.

2.7 Symmetric rank-one ( $\mathit{SR}1$ ) update

Let $H_k$ be the inverse Hessian approximation of the $k$ th iteration. We are trying to update $H_k$ into $H_{k+1}$ , i.e.

$$H_{k+1}=H_k+E_k,$$

where $E_k$ is a matrix with a lower rank. If it is about a rank-one update, we get

$$H_{k+1}=H_k+{\mathit{uv}}^T,$$ (46) E46

where $u,v\in {\mathbb{R}}^n$ . Using quasi-Newton equation (43), we can get

$$H_{k+1}{y}_k=\left(H_k+{\mathit{uv}}^T\right)y_k=s_k,$$

wherefrom

$$\left(v^T{y}_k\right)u=s_k-H_k{y}_k.$$ (47) E47

Further, from (46) and (47), we have

$$H_{k+1}=H_k+\frac{1}{v^T{y}_k}\left(s_k-H_k{y}_k\right)v^T.$$

Having in view that the inverse Hessian approximation $H_k$ has to be the symmetric one, we use $v=s_k-H_k{y}_k$ , so we get the symmetric rank-one update (i.e., $\mathit{SR}1$ update):

$$H_{k+1}=H_k+\frac{\left(s_k-H_k{y}_k\right){\left(s_k-H_k{y}_k\right)}^T}{{\left(s_k-H_k{y}_k\right)}^T{y}_k}.$$ (48) E48

Theorem 1.2.8. [27] (Property theorem of $\mathit{SR}1$ update) Let $s_0$ , $s_1$ , and $s_{n-1}$ be linearly independent. Then, for quadratic function with a positive definite Hessian, $\mathit{SR}1$ method terminates at $n+1$ steps, i.e., $H_n=G^{-1}$ .

More information about SR1 update can be found.

2.8 Davidon-Fletcher-Powell ( $\mathit{DFP}$ ) update

There exists another type of update, which is a rank-two update. In fact, we get $H_{k+1}$ using two symmetric, rank-one matrices:

$$H_{k+1}=H_k+{\mathit{auu}}^T+{\mathit{bvv}}^T,$$ (49) E49

where $u,v\in {\mathbb{R}}^n$ and $a,b$ are scalars which have to be determined.

Using quasi-Newton equation (43), we can get

$$H_k{y}_k+{\mathit{auu}}^T{y}_k+{\mathit{bvv}}^T{y}_k=s_k.$$ (50) E50

The values of $u,v$ are not determined in a unique way, but the good choice is

$$u=s_k,v=H_k{y}_k.$$

Now, from (50), we get:

$$a=\frac{1}{s_k^T{y}_k},b=-\frac{1}{y_k^T{H}_k{y}_k}.$$

Hence, we get the formula

$$H_{k+1}=H_k+\frac{s_k{s}_k^T}{s_k^T{y}_k}-\frac{H_k{y}_k{y}_k^T{H}_k}{y_k^T{H}_k{y}_k},$$ (51) E51

which is $\mathit{DFP}$ update.

Theorem 1.2.9. [27] (Positive definiteness of $\mathit{DFP}$ update) $\mathit{DFP}$ update (51) retains positive definiteness if and only if $s_k^T{y}_k>0$ .

Theorem 1.2.10. [27] (Quadratic termination theorem of $\mathit{DFP}$ method) Let $f\left(x\right)$ be a quadratic function with positive definite Hessian $G$ . Then, if the exact line search is used, the sequence $\left\{s_j\right\}$ , generated from $\mathit{DFP}$ method, satisfies, for $i=0,1,\dots ,m$ , where $m\le n-1$ :

1. $H_{i+1}{y}_j=s_j,j=0,1,\dots ,i$ $\left(\mathit{hereditary\hspace{0.17em}property}\right)$ .

2. $s_i^T{\mathit{Gs}}_j=0,j=0,1,\dots ,i-1$ $\left(\mathit{conjugate\hspace{0.17em}direction\hspace{0.17em}property}\right)$ .

3. $\mathit{The\hspace{0.17em}method\hspace{0.17em}terminates}\mathit{at}m+1\le n\mathit{\text{steps}}.\mathit{\text{If}}m=n-1,\mathit{\text{then}}{H}_n=G^{-1}$ .

2.9 Broyden-Fletcher-Goldfarb-Shanno ( $\mathit{\text{BFGS}}$ ) update

$\mathit{\text{BFGS}}$ update is given by the formula

$$B_{k+1}^{\mathit{\text{BFGS}}}=B_k+\frac{y_k{y}_k^T}{y_k^T{s}_k}-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}.$$ (52) E52

The $\mathit{\text{BFGS}}$ update is also said to be a complement to $\mathit{DFP}$ update.

In [62], an adaptive scaled $\mathit{\text{BFGS}}$ method for unconstrained optimization is presented. In this paper, the author emphasizes that the $\mathit{\text{BFGS}}$ method is one of the most efficient quasi-Newton methods for solving small-size and medium-size unconstrained optimization problems. The third term in the standard $\mathit{\text{BFGS}}$ update formula is scaled in order to reduce the large eigenvalues of the approximation to the Hessian of the minimizing function. In fact, in [62], the general scaling $\mathit{\text{BFGS}}$ updating formula is considered:

$$B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+{\gamma }_k\frac{y_k{y}_k^T}{y_k^T{s}_k},$$ (53) E53

where ${\gamma }_k$ is a positive parameter. Obviously, using ${\gamma }_k=1$ for all $k=0,1,\dots$ , we get the standard $\mathit{\text{BFGS}}$ formula. By the way, there exist several procedures created to select the scaling parameter ${\gamma }_k$ , for example, see [62, 63, 64, 65, 66, 67, 68, 69]. The approach for determining the scaling parameters of the terms of the $\mathit{\text{BFGS}}$ update in [62] is to minimize the Byrd and Nocedal measure function.

Namely, in [70], the next function was introduced:

$$\phi \left(A\right)=\mathit{tr}\left(A\right)-ln\left(\mathit{det}\left(A\right)\right),$$ (54) E54

which is defined on positive definite matrices.

This function is a measure of matrices involving all the eigenvalues of $A$ , not only the smallest one and the largest one, as it is traditionally used in the analysis of the quasi-Newton method based on the condition number of matrices.

Observe that function $\phi$ works simultaneously with the trace and the determinant, thus simplifying the analysis of the quasi-Newton methods. Fletcher [71] proves that this function is strictly convex on the set of symmetric and positive definite matrices, and it is minimized by $A=I$ . Besides, this function becomes unbounded when $A$ becomes singular or infinite, and therefore it works as a barrier function that keeps $A$ positive definite. It is worth saying that the $\mathit{\text{BFGS}}$ update tends to generate updates with large eigenvalues.

Further, in [62], a double-parameter scaling $\mathit{\text{BFGS}}$ update is considered, in which the first two terms on the right-hand side of the $\mathit{\text{BFGS}}$ update (52) are scaled with a positive parameter, while the third one is scaled with another positive parameter:

$$B_{k+1}={\delta }_k\left[B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}\right]+{\gamma }_k\frac{y_k{y}_k^T}{y_k^T{s}_k},$$ (55) E55

where ${\delta }_k$ and ${\gamma }_k$ are the two positive parameters that have to be determined.

In [62], the next proposition is proved.

Proposition 1.2.1. If the step size $t_k$ is determined by the standard Wolfe line search (12) and (13), $B_k$ is positive definite and ${\gamma }_k>0$ , and then $B_{k+1}$ , given by (55), is also positive definite.

From (55), it can be seen that $\phi \left(B_{k+1}\right)$ depends on the scaling parameters ${\delta }_k$ and ${\gamma }_k$ . In [62], these scaling parameters are determined as solution of the minimizing problem:

$$\underset{{\delta }_k>0,{\gamma }_k>0}{min}\phi \left(B_{k+1}\right).$$ (56) E56

Further, the next values of the scaling parameters ${\delta }_k$ and ${\gamma }_k$ are reached:

$${\delta }_k=\frac{n-1}{\mathit{tr}\left(B_k\right)-\frac{\parallel B_k{s}_k{\parallel }^2}{s_k^T{B}_k{s}_k}}$$ (57) E57

$${\gamma }_k=\frac{y_k^T{s}_k}{\parallel y_k{\parallel }^2}.$$ (58) E58

Consider the relation

$$x_{k+1}=x_k+t_k{d}_k,$$ (59) E59

where $d_k$ is the $\mathit{\text{BFGS}}$ search direction obtained as solution of the linear algebraic system

$$B_k{d}_k=-g_k,$$

where the matrix $B_k$ is the $\mathit{\text{BFGS}}$ approximation to the Hessian ${\nabla }^{2}f\left(x_k\right)$ , being updated by the classical formula (52).

The next theorems are also given in [62].

Theorem 1.2.11. If the step size in (59) is determined by the Wolfe search conditions (12)–(13), then the scaling parameters given by (57) and (58) are the unique global solutions of the problem (56).

Theorem 1.2.12. Let ${\delta }_k$ be computed by (57). Then, for any $k=0,1,\dots$ , ${\delta }_k$ is positive and close to $1$ .

Next, in [72], using chain rule, a modified secant equation is given, to get a more accurate approximation of the second curvature of the objective function. Then, based on this modified secant equation, a new $\mathit{\text{BFGS}}$ method is presented. The proposed method makes use of both gradient and function values, and it utilizes information from two most recent steps, while the usual secant relation uses only the latest step information. Under appropriate conditions, it is shown that the proposed method is globally convergent without convexity assumption on the objective function.

Some interesting applications of Newton, modified Newton, inexact Newton, and quasi-Newton methods can be found, for example, in [73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83], etc.

A very interesting paper is [84].

An interesting application of $\mathit{\text{BFGS}}$ method can be found in [85].

3. Conclusion

Today, the modifications of the line search techniques are very actual and all in the aim to create new, better optimization methods.

Further, following recent trends in unconstrained optimization, we can notice that almost all optimization methods, which are considered in this chapter, are still actual.

They are applied in the other areas of Mathematics, as well as in practice. Also, different modifications of these methods are made, in the aim to improve them.

Let us emphasize that $\mathit{\text{BFGS}}$ update is very popular now.