*Article* **A Family of Hybrid Stochastic Conjugate Gradient Algorithms for Local and Global Minimization Problems**

**Khalid Abdulaziz Alnowibet 1, Salem Mahdi 2, Ahmad M. Alshamrani 1, Karam M. Sallam <sup>3</sup> and Ali Wagdy Mohamed 4,\***


**Abstract:** This paper contains two main parts, Part I and Part II, which discuss the local and global minimization problems, respectively. In Part I, a fresh conjugate gradient (CG) technique is suggested and then combined with a line-search technique to obtain a globally convergent algorithm. The finite difference approximations approach is used to compute the approximate values of the first derivative of the function *f* . The convergence analysis of the suggested method is established. The comparisons between the performance of the new CG method and the performance of four other CG methods demonstrate that the proposed CG method is promising and competitive for finding a local optimum point. In Part II, three formulas are designed by which a group of solutions are generated. This set of random formulas is hybridized with the globally convergent CG algorithm to obtain a hybrid stochastic conjugate gradient algorithm denoted by HSSZH. The HSSZH algorithm finds the approximate value of the global solution of a global optimization problem. Five combined stochastic conjugate gradient algorithms are constructed. The performance profiles are used to assess and compare the rendition of the family of hybrid stochastic conjugate gradient algorithms. The comparison results between our proposed HSSZH algorithm and four other hybrid stochastic conjugate gradient techniques demonstrate that the suggested HSSZH method is competitive with, and in all cases superior to, the four algorithms in terms of the efficiency, reliability and effectiveness to find the approximate solution of the global optimization problem that contains a non-convex function.

**Keywords:** global optimization; unconstrained minimization; numerical approximations of gradients; meta-heuristics; stochastic parameters; conjugate gradient methods; efficient algorithm; performance profiles; comparisons; testing

**MSC:** 90C26

## **1. Introduction**

The major goal of this paper is to find the local and global minima of a convex and non-convex function. The local and global minimization problems are defined as follows.

**Definition 1.** *A local minimum xlo* <sup>∈</sup> <sup>S</sup> *of the function <sup>f</sup> , <sup>f</sup>* : <sup>S</sup> <sup>→</sup> <sup>R</sup> *is an input element with <sup>f</sup>*(*xlo*) <sup>≤</sup> *<sup>f</sup>*(*x*) *for all <sup>x</sup> neighboring xlo. If* <sup>S</sup> <sup>⊆</sup> <sup>R</sup>*n, it is formulated by*

$$\forall \mathbf{x}\_{lo} \; \exists \varepsilon > 0 : f(\mathbf{x}\_{lo}) \le f(\mathbf{x}) \; \forall \mathbf{x} \in \mathbb{B}, \left\| \mathbf{x} - \mathbf{x}\_{lo} \right\| \le \varepsilon. \tag{1}$$

**Citation:** Alnowibet, K.A.; Mahdi, S.; Alshamrani, A.M.; Sallam, K.M.; Mohamed, A.W. A Family of Hybrid Stochastic Conjugate Gradient Algorithms for Local and Global Minimization Problems. *Mathematics* **2022**, *10*, 3595. https://doi.org/ 10.3390/math10193595

Academic Editors: Humberto Rocha and Ana Maria Rocha

Received: 23 August 2022 Accepted: 24 September 2022 Published: 1 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Definition 2.** *The point xgl* <sup>∈</sup> <sup>S</sup> *is called the global minimizer of the function <sup>f</sup> ; <sup>f</sup>* : <sup>S</sup> <sup>→</sup> <sup>R</sup> *such that f*(*xgl*) <sup>≤</sup> *<sup>f</sup>*(*x*) <sup>∀</sup>*<sup>x</sup>* <sup>∈</sup> <sup>S</sup>*. When S* <sup>⊆</sup> <sup>R</sup>*n, then the problem can be formulated by*

$$\min\_{\mathbf{x}\in\mathbb{B}} f(\mathbf{x}) : \mathbb{B} \to \mathbb{R}, \tag{2}$$

In both problems (formulae) <sup>S</sup> <sup>⊆</sup> <sup>R</sup>*<sup>n</sup>* is the range in which we find the global minimizer of *f*(*x*). *f*(*x*) is continuously differentiable.

Global optimization (GO) attempts to find the approximate solution of the objective function are shown in Problem (2).

However, this task can be difficult since the knowledge about *f* is usually only local. On the other hand, the fastest algorithms (LO) prefer to find a local point since these algorithms are not capable of finding the global solution at each run.

The bottom line is that the core difference between the GO methods and the LO algorithms is as follows: the GO methods focus on solving Problem (2) over the given set, while the task of the LO methods is to solve (1). Consequently, solving Problem (1) is relatively simple by using deterministic (classical) local optimization methods. On the contrary, finding the global optimum of Problem (2) is an NP-hard problem.

Challenging problems arise in different application fields, for example, technical sciences, industrial engineering, economics, networks, chemical engineering, etc. See [1–11].

Recently, many optimization algorithms have been proposed to deal with these problems. The thoughts of those suggested methods rely on the standard of meta-heuristic strategies (random search).

There are different classifications for meta-heuristic methods [12].

Mohamed et al. [7] presented a brief description of these classifications.

In random algorithms, the minimization technique relies partly on probability.

In contrast, in the deterministic algorithms, a guessing scale is not utilized. Hence, deterministic techniques need an exhaustive examination over the research domain of function *f* to find the approximate solution to Problem (2) at each run. Otherwise, they fail in this task.

Therefore, finding the approximate solution to Problem (2) by using random techniques can be proved by the asymptotic convergence probability. See [13–15].

There are many deterministic methods that have been proposed for dealing with the local optimization problems. See, for example, Refs. [16–20].

The most popular deterministic method is the CG method [18]. CG methods are exceedingly utilized to find the local minimizer of Problem (1) [21].

However, the CG algorithms have a numerical weakness, so their subsequent actions might be low if a little step is created away from the local point. Hence, for solving this issue, a line-search technique is combined with the CG technique to create a globally convergent algorithm [22,23].

Therefore, many conjugant gradient line-search methods are suggested; see, for example, refs. [18,24–28].

The CG method is an efficient and inexpensive technique to deal with Problem (1).

The CG method is an iterative algorithm. Therefore, the candidate solutions are generated by the following recursive formula.

$$\mathbf{x}\_{k+1} = \mathbf{x}\_k + \mathbf{a}\_k \mathbf{d}\_{k'} \tag{3}$$

where the step size *α<sup>k</sup>* > 0, and the directions *dk* are created by the following formula:

$$d\_{k+1} = -\mathbb{g}\_{k+1} + \beta\_k d\_{k'} d\_{\mathfrak{o}} = -\mathbb{g}\_{\mathfrak{o}}.\tag{4}$$

where *gk* denotes the gradient vector of the function *f* at the point *xk* .

Several versions of the CG methods are suggested. The core difference between those CG algorithms relies on choosing the parameter *β<sup>k</sup>* [18,27–29]. The main features of the CG method are as follows: it has low memory requirements, it is strongly local, and it has global convergence properties [30].

Many authors presented several studies to analyze the CG method; see, for example, Refs. [31,32].

In 1964, the authors of [33] applied the CG methods to nonlinear problems, and they proposed the following parameter.

$$
\boldsymbol{\beta}\_k^{FR} = \frac{||\mathbf{g}\_{k+1}||^2}{||\mathbf{g}\_k||^2}.\tag{5}
$$

The authors of [34,35] established the global convergence of the scheme defined in (5); they used an exact line search and an inexact line search respectively.

However, the author of [36] showed that there are some cases that have some strays; these jamming occurrences happen when the search directions *dk* are almost orthogonal to the gradient vector *gk* [18].

The authors of [37,38] presented a modification of the parameter *βFR <sup>k</sup>* for treating the noise event denoted in [36]. Hence, they proposed the following parameter.

$$\boldsymbol{\beta}\_{k}^{PRP} = \frac{\mathbf{y}\_{k}^{T}\mathbf{g}\_{k+1}}{||\mathbf{g}\_{k}||^{2}},\tag{6}$$

where *yk* <sup>=</sup> *gk***+<sup>1</sup>** <sup>−</sup> *gk* . When a noise occurs *gk***+<sup>1</sup>** <sup>≈</sup> *gk* , *<sup>β</sup>PRP <sup>k</sup>* ≈ 0, and *dk***+<sup>1</sup>** ≈ −*gk***+1**, i.e., when jamming happens, the search direction *dk* is no longer perpendicular to the gradient vector *gk* , but it is aligned with the vector −*gk* . This built-in restart advantage of the *βPRP <sup>k</sup>* parameter usually has better quick convergence when compared to the parameter *βFR <sup>k</sup>* [18].

The authors of [39] proposed an approach closely related to *βPRP <sup>k</sup>* , and it is defined as follows.

$$\boldsymbol{\beta}\_{k}^{HS} = \frac{\boldsymbol{y}\_{k}^{T}\mathcal{G}\boldsymbol{k} + \mathbf{1}}{\boldsymbol{d}\_{k}^{T}\boldsymbol{y}\_{k}}.\tag{7}$$

in the case that step-size *α<sup>k</sup>* is found by an exact line search algorithm. Hence, by (4) and the orthogonality situation *g<sup>T</sup> <sup>k</sup>***+1***yk* = 0, the following can be obtained:

$$d\_k^T \mathbf{y}\_k = (\mathbf{g}\_{k+1} - \mathbf{g}\_k)^T d\_k = -d\_k^T \mathbf{g}\_k = ||\mathbf{g}\_k||^2. \tag{8}$$

Therefore, *βHS <sup>k</sup>* <sup>=</sup> *<sup>β</sup>PRP <sup>k</sup>* when the step size *α<sup>k</sup>* is calculated by an exact line search method. Other fundamentals formulas of the parameter *β<sup>k</sup>* which contain one term are listed as follows.

$$\boldsymbol{\beta}\_{k}^{LS} = \frac{\mathbf{g}\_{k+1}^{T} \mathbf{y}\_{k}}{-\mathbf{d}\_{k}^{T} \mathbf{g}\_{k}}.\tag{9}$$

Formula (9) was proposed by [40].

$$\boldsymbol{\beta}\_{k}^{DY} = \frac{||\mathbf{g}\_{k+1}||^2}{\mathbf{y}\_{k}^T \mathbf{d}\_{k}}.\tag{10}$$

Formula (10) was proposed by Dai and Yuan [41]. It is noteworthy that when the *f* is quadratic and step size *α<sup>k</sup>* is selected to reduce *f* along *dk* , the options of the parameter *β<sup>k</sup>* mentioned above are alike for the generic nonlinear function.

Different alternatives have fully different convergence possessions [18].

Many version of the parameter *β<sup>k</sup>* have been proposed in two- and three terms; see, for example, Refs. [32,42–50].

For example, in the following two approaches, we present some modifications to obtain a new CG method. See Section 2.

$$\boldsymbol{\beta}\_{k}^{HZ} = \frac{(\mathbf{y}\_{k}^{T}\mathbf{g}\_{k})(\mathbf{d}\_{k-1}^{T}\mathbf{y}\_{k}) - 2||\mathbf{y}\_{k}||^{2}(\mathbf{d}\_{k-1}^{T}\mathbf{g}\_{k})}{(\mathbf{d}\_{k-1}^{T}\mathbf{y}\_{k})^{2}}.\tag{11}$$

Formula (11) was proposed by [30].

$$\boldsymbol{\beta}\_{k}^{MHZ} = \frac{(\boldsymbol{y}\_{k}^{T}\mathbf{g}\_{k})(\boldsymbol{d}\_{k-1}^{T}\mathbf{y}\_{k}) - 2||\boldsymbol{y}\_{k}||^{2}(\boldsymbol{d}\_{k-1}^{T}\mathbf{g}\_{k})}{\max\{\sigma||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{d}\_{k}||^{2}, (\boldsymbol{d}\_{k-1}^{T}\mathbf{y}\_{k})^{2}\}},\tag{12}$$

where *σ* > 0.5 is a constant. Formula (12) was proposed by [49]. The denominator (*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )<sup>2</sup> in the *<sup>β</sup>HZ <sup>k</sup>* is modified to max{*σyk* 2*dk* 2,(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2} in the *<sup>β</sup>MHZ <sup>k</sup>* . This procedure may help the *dk* stay in a trusted area automatically beneath each iteration [49]. Furthermore, in a situation *<sup>σ</sup>*||*yk* ||2||*dk* ||<sup>2</sup> <sup>&</sup>lt; (*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2, *<sup>β</sup>MHZ <sup>k</sup>* decreases to *<sup>β</sup>HZ <sup>k</sup>* with *α<sup>k</sup>* calculated to satisfy the inexact line search. Moreover, *βHZ <sup>k</sup>* decreases to *<sup>β</sup>HS <sup>k</sup>* under the exact line search.

Consequently, by using a line search method, the CG method can satisfy the following descent condition:

$$\|\mathbf{g}\_k^T \mathbf{d}\_k \le -\mathbb{C} \|\mathbf{g}\_k\|^2,\tag{13}$$

where *C* > 0 is a constant.

The sufficient descent condition (13) has a core task in the convergence analysis of the algorithms. See [17,30–32,35,41,49,51,52].

However, the CG method has a numerical obstacle; its sub-sequential phases might be low if a little step is created away from the intended point [49].

Recently, the authors of [48,49] proved that the CG algorithm includes powerful convergence features if it satisfies the trust-region feature that is determined by

$$||\mathbf{d}\_{k}|| \, \leq \, \mathbb{C}\_{v} ||\mathbf{g}\_{k}|| \, \vert \, \tag{14}$$

where *Cv* > 0 is a constant. It is shown, therefore, that the trust-region property can enable the search direction *dk* to be bounded in the trust radius [49]. Numerous researchers proposed many CG algorithms that give perfect results and powerful convergence properties. See [30,48,49,51].

The selection of the right step size *α<sup>k</sup>* can help the CG algorithms to achieve global convergence.

The exact line search is defined as follows:

$$f(\mathbf{x}\_k + a\_k d\_k) = \min\_{a \ge 0} \theta(a) = f(\mathbf{x}\_k + a d\_k). \tag{15}$$

It is clear that in big-scale problems, the exact line search cannot be used.

Therefore, there are many techniques to achieve this task. Formula (15), for example, the weak Wolfe–Powell algorithm (WWP), is a popular technique, and it is exceedingly utilized. The WWP technique is designed to find the step size *α<sup>k</sup>* to satisfy the following inequalities:

$$f(\mathbf{x}\_k + \mathbf{a}\_k \mathbf{d}\_k) \le f(\mathbf{x}\_k) + \delta \mathbf{a}\_k \mathbf{g}\_k^T \mathbf{d}\_{k'} \tag{16}$$

and

$$\text{sg}(\mathbf{x}\_k + \mathbf{a}\_k \mathbf{d}\_k)^T \mathbf{d}\_k \ge \sigma \mathbf{g}\_k^T \mathbf{d}\_{k'} \tag{17}$$

where *δ* ∈ (0, 0.5) and *σ* ∈ (*δ*, 1) are constants.

Inequality (16) is named the Armijo condition, and the WWP line search decreases to strong Wolfe–Powell (SWP) by substituting Inequality (17) with the following inequality:

$$|\mathbf{g}(\mathbf{x}\_k + \boldsymbol{\alpha}\_k \mathbf{d}\_k)^T \mathbf{d}\_k| \le -\sigma \mathbf{g}\_k^T \mathbf{d}\_{k'} \tag{18}$$

Generally, under the WWP line search, it is assumed that the gradient *g***(***x***)** is Lipschitz continuous in the convergence analysis. Therefore, the following inequality is satisfied:

$$||\mathbf{g(x)} - \mathbf{g(y)}|| \le L||\mathbf{x} - \mathbf{y}||,\tag{19}$$

with *<sup>L</sup>* is a constant <sup>∀</sup> *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> <sup>R</sup>*n*.

In fact, the CG technique with the line search methods has proven notability in solving the local optimization problem [18,27,28]. However, in trying to solve Problem (2), the CG method fails to achieve this task per run because it is trapped to a local point. To prevent sticking in a local point, random parameters are used [53].

We can summarize the essence of the above discussions as follows.

Recently, there have been many and many proposed approaches presented to improve the performance of deterministic methods, such as CG methods, gradient descent methods, Newton methods, etc. Those new approaches are designed to deal with the local optimization problems. See, for example, Refs. [16–20].

On the other hand, a plentiful number of stochastic approaches are suggested to deal with the global optimization problems. See, for example, Refs. [1,2,4,5,7,54].

Therefore, to gain the features of both deterministic and stochastic methods, many studies presented several ideas and suggestions to combine deterministic and stochastic techniques to obtain a new technique that is efficient and effective in solving Problem (2). Numerical outcomes demonstrated that the interbreed between classical and stochastic techniques has been hugely successful. See [55–59].

This work focuses on solving the local and global minimization problems. So, the first part of this study trades with Problem (1) by suggesting a new modified CG method, while the second part of this paper presents a new random approach that includes three formulae by which the candidate solutions are generated randomly.

Therefore, the new proposed stochastic approach is combined with the new modified CG method that is proposed in the first part of this paper to obtain a new hybrid stochastic conjugate gradient algorithm that solves Problem (2). The new hybrid stochastic conjugate gradient algorithm has four formulae by which the candidate solutions are created. One of the four formulae is a purely deterministic formula, the second one is a mixture of deterministic and stochastic parameters, and the other two formulas contain parameters generated randomly. The bottom line is that we can claim that the main merit that makes the new hybrid algorithm capable of finding the approximate solution to the global minimum of a non-convex function comes from the hybridization of random and non-random parameters.

Consequently, the contribution of this paper is divided into two parts.

Part I presents the following contributions.


Part II presents the following contributions.


Consequently, the remainder of the study is arranged as follows.

Part I contains the following sections: Section 2 presents a new modified CG- SHZ technique with its convergence analysis.

In Section 3, the approximate value of the gradient vector is calculated by using the numerical differentiation. Section 4 presents the numerical investigations of the local minimization problem. Part II contains the following sections: Section 5 presents a random approach for unconstrained global optimization. Section 6 presents the hybridization of the conjugate gradient method with stochastic parameters. The numerical experiments of Problem (2) are presented in Section 7. Some concluding remarks are given in Section 8.

## *Part I: Local Minimization Problem*

In this part, a new modified CG technique is presented, the convergence analysis of this technique is designed, the numerical differentiation approach is utilized to calculate the approximate values of the first derivative, the five algorithms are designed to solve Problem (1), and their numerical experiments are analyzed by using the performance profiles.

#### **2. Suggested CG Method**

Recently, the authors of [49] suggested a new MHZ-CG method, relying on the study which was proposed by the authors of [30]. The MHZ method contains the sufficient descent and the trust-region features independent of a line search technique. The parameter of the MHZ is defined by (12).

Therefore, the story in this section begins with the authors of [30] who proposed a new CG-HZ method, where the parameter of the HZ method is defined by (11). The parameter *βHZ <sup>k</sup>* can ensure that *dk* satisfies the following inequality:

$$d\_k^T \mathbf{g}\_k \le -\frac{7}{8} ||\mathbf{g}\_k||^2,\tag{20}$$

where (20) is proved by [30]. If the step size *α<sup>k</sup>* is calculated by the true line search, then *βHZ <sup>k</sup>* decreases to the *<sup>β</sup>HS <sup>k</sup>* that was proposed by [39] because *<sup>d</sup><sup>T</sup> <sup>k</sup> gk* = 0 is true [49].

Hence, for obtaining the global convergence for a general function, Hager and Zhang [30] dynamically adjusted the down limitation of *βHZ <sup>k</sup>* by

$$d\_k = -\mathbf{g}\_k + \beta\_k^{lIZ^+} d\_{k-1}, \\ d\_\mathfrak{o} = -\mathbf{g}\_{\mathfrak{o}'} \tag{21}$$

*βHZ*<sup>+</sup> *<sup>k</sup>* <sup>=</sup> max{*βHZ*,*rk*}, *rk* <sup>=</sup> <sup>−</sup><sup>1</sup> ||*dk−***1**|| min{*r*,||*gk−***<sup>1</sup>** ||} , where *<sup>r</sup>* <sup>&</sup>gt; 0 is a constant.

Many researchers have suggested several modifications and refinements to improve the performance of the CG-HZ algorithm. The latest version of the CG-HZ method was offered by [49]. Yuan et al. [49] presented some modifications to the HZ-CG method, and the result was obtaining the new CG-MHZ algorithm.

The CG-MHZ algorithm contains a sufficient condition and the trust-region feature. The research direction of the MHZ-CG technique is designed as follows:

$$d\_k = -\mathbf{g}\_k + \boldsymbol{\beta}\_k^{MHZ} d\_{k-1}, \\ d\_\mathbf{o} = -\mathbf{g}\_{\mathbf{o}'} \tag{22}$$

where the *βMHZ* is defined by (12).

In this paper, the MHZ method is extended and modified to obtain a new proposed method called the SHZ method such that the SHZ method has a sufficient condition and the trust-region feature. This method is defined as follows:

$$d\_k = -\mathbf{g}\_k + \boldsymbol{\beta}\_k^{SHZ} d\_{k-1}, \\ d\_\mathbf{o} = -\mathbf{g}\_{\mathbf{o}\_{\mathbf{o}'}} \tag{23}$$

$$\boldsymbol{\beta}\_{k}^{SHZ} = \frac{(\boldsymbol{y}\_{k}^{T}\boldsymbol{y}\_{k})(\boldsymbol{d}\_{k-1}^{T}\boldsymbol{y}\_{k}) - 2||\boldsymbol{y}\_{k}||^{2}(\boldsymbol{d}\_{k-1}^{T}\boldsymbol{g}\_{k})}{\max\{\vartheta||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{d}\_{k}||^{2}, (\boldsymbol{d}\_{k-1}^{T}\boldsymbol{y}\_{k})^{2}\}},\tag{24}$$

where the *ϑ* = max{*ρ*, *Rk*}, the *ρ* and *Rk* are defined as follows. The parameter *ρ* is changed randomly at each iteration and its values are taken from the range [0.8, 2) and *Rk* = *f x*. The values of *f* and *x* are calculated by

$$
\triangle f = |f\_0 - f\_{tr}|\_\prime \tag{25}
$$

where *Itr* is the number of iterations, and after the *Itr* number of iterations, *fItr* and *f* are computed. Then, we set *f*<sup>0</sup> = *fItr*, while *x* is defined by

$$
\triangle \mathbf{x} = ||\mathbf{x}\_{k+1} - \mathbf{x}\_k||\_\prime \text{ for } k = 0, 1, \dots, Itr. \tag{26}
$$

Hence, when *ϑ* = *σ*, *βSHZ <sup>k</sup>* inevitably reduces to one of the following methods {*βMHZ <sup>k</sup>* , *<sup>β</sup>HZ <sup>k</sup>* , *<sup>β</sup>HS <sup>k</sup>* } as follows.

If *<sup>ϑ</sup>* <sup>=</sup> *<sup>σ</sup>* and *<sup>δ</sup>yk* 2*dk* <sup>2</sup> <sup>&</sup>gt; (*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2, the *<sup>β</sup>SHZ <sup>k</sup>* reduces to the *<sup>β</sup>MHZ <sup>k</sup>* . Otherwise, *βSHZ <sup>k</sup>* reduces to *<sup>β</sup>HZ <sup>k</sup>* or to *<sup>β</sup>HS <sup>k</sup>* under the exact line search [49]. This procedure gives the advantages of the MHZ, HZ and HS methods to the proposed SHZ method. In other words, the SHZ algorithm gains the characteristics of the three MHZ, HZ and HS algorithms. This is why the SHZ algorithm is superior to the four other MHZ, HZ, HS and FR methods.

**Note:** The authors of [49] imposed that the *σ* > 0.5 is a constant, while the parameter *ϑ* is modified dynamically at each iteration.

## *Convergence Analysis of Algorithm 1*

In this section, we present the features of Algorithm 1. We also present the convergence analysis of this algorithm, and we show that the search direction *dk* that is defined by Formula (23) satisfies the sufficient descent condition and the trust-region merit, which are defined by Formulae (13) and (14), respectively.

## **Algorithm 1** A conjugate gradient method (CG-SHZ).

**Input:** *<sup>f</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>, *<sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*1, *<sup>γ</sup>* <sup>∈</sup> (0, 1) ,*<sup>k</sup>* <sup>=</sup> 0, a starting point *xk* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* and *<sup>ε</sup>* <sup>&</sup>gt; 0. **Output:** *x<sup>∗</sup>* = *xloc* the local minimizer of *f* , *f*(*x*∗), the value of *f* at *x<sup>∗</sup>* 1: Set *d***<sup>0</sup>** = −*g***<sup>0</sup>** and *k* := 0. 2: **while** *gk* > *ε*. **do** 3: compute *α<sup>k</sup>* to satisfy (16) and (17). 4: Calculate a new point *xk*<sup>+</sup><sup>1</sup> = *xk* + *αkdk* . 5: compute *f <sup>k</sup>* = *f*(*xk***<sup>+</sup><sup>1</sup>** ), *gk* = *g***(***xk***<sup>+</sup><sup>1</sup>** ) 6: Set *k* = *k* + 1. 7: calculate the search direction *dk* by (23). 8: **end while** 9: **return** *xac* the local minimizer and its function value *fac*

Two sensible hypotheses are assumed as follows.

**Hypothesis 1.** *We suppose that Problems (1) and (2) contain an objective function f*(*x*) *with the following characteristics: continuity and differentiability properties.*

**Hypothesis 2.** *In some neighborhood* ℵ *of the level set*

$$\ell = \{ \mathbf{x} \in \mathbb{R}^n : f(\mathbf{x}) \le f(\mathbf{x\_0}) \},$$

*the gradient vector g***(***x***)** *is Lipschitz continuous. This means that there is a fixed real number L* < ∞ *such that*

$$\|\|\mathbf{g}(\mathbf{x}) - \mathbf{g}(y)\|\| \le L \|\mathbf{x} - \mathbf{y}\|\_{\prime}$$

*for all x*, *y* ∈ ℵ*.*

**Lemma 1.** *Suppose that the sequence* {*xk* } *is obtained by Algorithm 1. If <sup>d</sup><sup>T</sup> <sup>k</sup> yk* = 0*, then*

$$\|\mathbf{g}\_k^T \mathbf{d}\_k \le -c \|\mathbf{g}\_k\|^2,\tag{27}$$

*and*

$$||d\_k|| \le r\_v ||\mathcal{g}\_k||\_\prime \tag{28}$$

*where <sup>c</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>7</sup> <sup>9</sup>*<sup>ϑ</sup>* <sup>&</sup>gt; <sup>0</sup>*, <sup>ϑ</sup>* <sup>=</sup> max{*ρ*, *Rk*}*, <sup>ρ</sup> is taken randomly from* <sup>∈</sup> [ <sup>8</sup> <sup>10</sup> , 2) *at each iteration of Algorithm 1,* <sup>0</sup> <sup>≤</sup> *Rk* <sup>&</sup>lt; <sup>∞</sup>*, and rv* = (<sup>1</sup> <sup>+</sup> <sup>3</sup> *<sup>ϑ</sup>* ) *is the trust-region radius.*

**Proof.** If *<sup>k</sup>* <sup>=</sup> 0, *<sup>d</sup>***<sup>0</sup>** <sup>=</sup> <sup>−</sup>*g***<sup>0</sup>** , then *<sup>g</sup><sup>T</sup>* **<sup>0</sup>** *<sup>d</sup>***<sup>0</sup>** <sup>=</sup> −||*g***<sup>0</sup>** ||<sup>2</sup> and ||*d***<sup>0</sup>** || <sup>=</sup> *g***<sup>0</sup>** , which indicates (27) and (28) by picking *c* ∈ (0, 1] and *rv* ∈ [1, ∞) .

Merging (23) with (24), the result is obtaining the following:

$$\mathbf{g}\_{k}^{T}\mathbf{d}\_{k} = \frac{(\mathbf{y}\_{k}^{T}\mathbf{g}\_{k})(\mathbf{d}\_{k-1}^{T}\mathbf{y}\_{k})(\mathbf{g}\_{k}^{T}\mathbf{d}\_{k-1}) - 2||\mathbf{y}\_{k}||^{2}(\mathbf{g}\_{k}^{T}\mathbf{d}\_{k-1})^{2}}{\max\{\vartheta ||\mathbf{y}\_{k}||^{2}||\mathbf{d}\_{k-1}||^{2}, (\mathbf{d}\_{k-1}^{T}\mathbf{y}\_{k})^{2}\}} - \|\mathbf{g}\_{k}\|^{2}.\tag{29}$$

The following inequality *<sup>u</sup><sup>T</sup> <sup>v</sup>* <sup>≤</sup> <sup>1</sup> <sup>2</sup> (||*u*||<sup>2</sup> <sup>+</sup> *v*2) is applied to the first term of the numerator of Inequality (29), where *<sup>u</sup>* <sup>=</sup> *dk−***1***g<sup>T</sup> <sup>k</sup> yk* , *v* = *yk gk <sup>T</sup>dk−***1**, and it is clear that *<sup>u</sup><sup>T</sup> <sup>v</sup>* <sup>≤</sup> <sup>7</sup> <sup>9</sup> (||*u*||<sup>2</sup> <sup>+</sup> *v*2) is right.

Therefore, the following inequality obtains

*gT <sup>k</sup> dk* <sup>=</sup> (*y<sup>T</sup> <sup>k</sup> gk* )(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )(*gk <sup>T</sup>dk−***1**) <sup>−</sup> <sup>2</sup>||*yk* ||2(*gk <sup>T</sup>dk−***1**)<sup>2</sup> max{*ϑyk* 2*dk−***1**2,(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2} − *gk* <sup>2</sup> <sup>≤</sup> −*gk* <sup>2</sup> <sup>+</sup> 7 <sup>9</sup> ||*yk* ||2*gk* 2||*dk−***1**||<sup>2</sup> <sup>+</sup> <sup>7</sup> <sup>9</sup> ||*yk* ||2(*gk <sup>T</sup>dk−***1**)<sup>2</sup> <sup>−</sup> <sup>2</sup>||*yk* ||2(*gk <sup>T</sup>dk−***1**)<sup>2</sup> max{*ϑyk* 2*dk−***1**2,(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2} <sup>=</sup> −*gk* <sup>2</sup> <sup>+</sup> 7 <sup>9</sup> ||*yk* ||2*gk* 2||*dk−***1**||<sup>2</sup> <sup>−</sup> <sup>11</sup> <sup>9</sup> ||*yk* ||2(*gk <sup>T</sup>dk−***1**)<sup>2</sup> max{*ϑyk* 2*dk−***1**2,(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2} <sup>≤</sup> −*gk* <sup>2</sup> <sup>+</sup> 7 <sup>9</sup> ||*yk* ||2*gk* 2||*dk−***1**||<sup>2</sup> max{*ϑyk* 2*dk−***1**2,(*d<sup>T</sup> <sup>k</sup>−***<sup>1</sup>***yk* )2} <sup>≤</sup> ( <sup>7</sup> <sup>9</sup>*<sup>ϑ</sup>* <sup>−</sup> <sup>1</sup>)*gk* 2,

such that

$$\max \left\{ \theta \left\| \left| y\_k \right| \right\|^2 \left\| d\_{k-1} \right\| \right\|^2, \left( d\_{k-1}^T y\_k \right)^2 \right\} \ge \theta \left\| y\_k \right\|^2 \left\| d\_{k-1} \right\|^2,\tag{30}$$

where *<sup>ϑ</sup>* <sup>=</sup> max{*ρ*, *Rk*}. Since *<sup>ϑ</sup>* <sup>≥</sup> <sup>8</sup> <sup>10</sup> and *<sup>c</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>7</sup> <sup>9</sup>*<sup>ϑ</sup>* > 0, (27) is true. By using (30), it is obvious that

$$||\boldsymbol{d}\_{k}|| = \left|| -\mathbf{g}\_{k} + \frac{(\boldsymbol{y}\_{k}^{T}\mathbf{g}\_{k})(\boldsymbol{d}\_{k-1}^{T}\mathbf{y}\_{k}) - 2||\boldsymbol{y}\_{k}||^{2}(\boldsymbol{d}\_{k-1}^{T}\mathbf{g}\_{k})}{\max\{\vartheta ||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{d}\_{k}||^{2}/(\boldsymbol{d}\_{k-1}^{T}\mathbf{y}\_{k})^{2}\}}\boldsymbol{d}\_{k-1} \right|| \leq$$

$$||-\mathbf{g}\_{k}|| + \frac{||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{g}\_{k}|| ||\boldsymbol{d}\_{k-1}||^{2} + 2||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{g}\_{k}|| ||\boldsymbol{d}\_{k-1}||^{2}}{\vartheta ||\boldsymbol{y}\_{k}||^{2}||\boldsymbol{d}\_{k-1}||^{2}} = \left(1 + \frac{3}{\vartheta}\right) ||\boldsymbol{g}\_{k}||$$

Consequently, (28) is met, where *rv* <sup>∈</sup> [<sup>1</sup> <sup>+</sup> <sup>3</sup> *<sup>ϑ</sup>* , ∞). The proof is complete.

**Corollary 1.** *According to Formula (28) of Lemma 1, the following formula is met.*

$$\sum\_{k=0}^{\infty} \frac{||\mathbf{g}\_k||^4}{||\mathbf{d}\_k||^2} = \infty. \tag{31}$$

**Proof.** Since *dk* ≤ *rv gk* 2, where 1 <sup>&</sup>lt; *rv* <sup>&</sup>lt; <sup>∞</sup>, then *dk* <sup>2</sup> <sup>≤</sup> *<sup>r</sup>*<sup>2</sup> *<sup>v</sup> gk* 4, therefore, *dk* <sup>2</sup> *gk* <sup>4</sup> ≤ *r*2 *<sup>v</sup>* , hence *gk* <sup>4</sup> *dk* <sup>2</sup> <sup>≥</sup> <sup>1</sup> *r*2 *v* . Now, the final expression is summed as *k* → ∞. The result is obtaining the following inequality: ∞ ∑ *k*=0 *gk* <sup>4</sup> *dk* <sup>2</sup> <sup>≥</sup> ∞ ∑ *k*=0 1 *r*2 *v* <sup>=</sup> <sup>1</sup> *r*2 *v* ∞ ∑ *k*=0 1 = ∞. Therefore, (31) is met.

Under the assumptions, we give a helpful lemma that was basically proved by Zoutendijk [60] and Wolfe [61,62].

**Lemma 2.** *Assume that the x*<sup>0</sup> *is the initial point by which Assumption 1 is satisfied. Regarding any algorithm of Formula (23), dk is a descent direction, and α<sup>k</sup> satisfies the standard Wolfe conditions (16) and (17). Hence, the following inequality is met:*

$$\sum\_{k=0}^{\infty} \frac{(\mathbf{g}\_k^T \mathbf{d}\_k)^2}{||\mathbf{d}\_k||^2} < \infty \tag{32}$$

**Proof.** It tracks Formula (17), such that

$$d\_k^T \mathbf{y}\_k = d\_k^T (\mathbf{g}\_{k+1} - \mathbf{g}\_k) \ge (\sigma - 1)\mathbf{g}\_k^T d\_k. \tag{33}$$

On the other hand, the Lipschitz condition (19) implies

$$(\mathbf{g}\_{k+1} - \mathbf{g}\_k)^T \mathbf{d}\_k \le \boldsymbol{\alpha}\_k \boldsymbol{L} \|\boldsymbol{d}\_k\|^2. \tag{34}$$

The above two inequalities give

$$\alpha\_k \ge \frac{\sigma - 1}{L} \cdot \frac{\mathbf{g}\_k^{\top} \mathbf{d}\_k}{||\mathbf{d}\_k||^2},\tag{35}$$

which with (16) implies that

$$f\_k - f\_{k+1} \ge c \frac{(\mathbf{g}\_k \, ^T \mathbf{d}\_k)^2}{||\mathbf{d}\_k||^2},\tag{36}$$

where *c* = *<sup>δ</sup>*(1−*σ*) *<sup>L</sup>* . By summing (36) and with the observation that *f* is limited below, we see that (32) holds, which concludes the proof.

**Theorem 1.** *Suppose that Hypotheses 1 and 2 hold, and by utilizing the outcome of Corollary 1, the sequence* {*gk* } *that is generated by Algorithm 1 satisfies the following:*

$$\lim\_{k \to \infty} \inf \|\mathbf{g}\_k\| = 0,\tag{37}$$

**Proof.** By contradiction, suppose that (37) is not true; then, for some > 0, the following inequality is true:

$$\|\mathbf{g}\_k\| \ge \epsilon. \tag{38}$$

Hence, with inequality (38) and (27), we obtain

$$\left\|\mathbf{g}\_{k}\right\|^{T}\mathbf{d}\_{k}\leq -c\left\|\mathbf{g}\_{k}\right\|^{2}\leq -\epsilon^{2}.\tag{39}$$

Then, we have

$$\frac{\mathcal{g}\_k{}^T d\_k}{||d\_k||} \le \frac{-\epsilon^2}{||d\_k||};$$

$$\frac{\mathcal{g}\_k{}^T d\_k}{||d\_k||} \ge \frac{\epsilon^4}{||d\_k||^2};$$

and by summing the final expression, we obtain

$$\sum\_{k=0}^{\infty} \frac{(\mathbf{g}\_k^T \mathbf{d}\_k)^2}{||\mathbf{d}\_k||^2} \ge \sum\_{k=0}^{\infty} \frac{\epsilon^4}{||\mathbf{d}\_k||^2} = \infty. \tag{40}$$

Therefore, the above leads to a contradiction with (32). So, (37) is met.

**Note 1:** The search direction *dk* that is defined by Formula (23) satisfies the sufficient descent condition which is defined by Formula (13).

**Note 2:** Lemma 1 guarantees that Algorithm 1 has a sufficient descent property and the trust-region feature automatically.

**Note 3:** Theorem 1 confirms that the series {*gk* } that is obtained by Algorithm 1 approaches to 0 as long as *k* → ∞.

In the next section, the numerical differentiation approach is discussed by which the first derivative is estimated and the step size *α<sup>k</sup>* is computed.

#### **3. Numerical Differentiation**

We now turn our attention to the numerical approximation to compute the approximate value of the gradient vector. In precept, it can be possible to find an analytic form for the first derivative for any continuous and differentiable function. However, in some cases, the analytic form is very complicated. The numerical approximation of the derivative may be sufficient for some purposes.

In this paper, the values of the *α<sup>k</sup>* , *gk* and the direction *dk* are computed by using the numerical differentiation method. Moreover, we have another step size and research directions that are generated randomly.

Several suggested methods have given fair outcomes for computing the gradient vector values numerically. See [63–67].

The common approaches by which the first derivative is computed are the finite difference approximation methods. Therefore, the first derivative *f* (*x*) can be estimated by the following numerical differentiation formula:

$$D\_f f(\mathbf{x}\_i) = \frac{f(\mathbf{x}\_{i+1}) - f(\mathbf{x}\_i)}{\mathbf{x}\_{i+1} - \mathbf{x}\_i} = \frac{f(\mathbf{x}\_i + h) - f(\mathbf{x}\_i)}{h},\tag{41}$$

where *h* is limited and little, but it is not necessarily infinitesimally small.

Reasonably, if the value of the *h* is small, the approximated value of the first derivative may improve. The forward difference and the central difference are the familiar and common methods used in many studies; see for example, [68–72].

The Taylor series can be used to derive these formulas. Thus, 3, 4 and 5 points can be utilized to derive these formulas, but it will be more costly than utilizing 2 points. The central difference method is known to include aspects of both accuracy and precision [73] but it needs 2*n* function evaluations against the forward-difference approximation approach, which needs *n* function evaluations for each iteration. So, in this study, the forwarddifference approximation approach is used, because it is a cheap method and it has sensible precision [66,68].

The advantage of the finite difference approximation approaches relies on choosing the fit values of the *h*.

Error approximation of the first derivative is discussed in the next section.

Therefore, the discussion of the error analysis guides us to define an appropriate finitedifference interval for the forward-difference approximation that balances the truncation error that grows from the error in the Taylor formula, and the magnitude error that is obtained from noise during computing the function values [66].

## *3.1. Error Analysis*

Formula (41) contains the forward-difference approximation form that is used to estimate the first derivative of the function *f* . Its errors are proportional to some power of the values of *h*. Therefore, it appears that the errors go on to reduce if *h* is reduced. However, it is a part of the problem since it is assumed only the truncation error yielded by truncating the high-order terms in the Taylor series expansion and does not take into account the round-off error induced by quantization. The round-off error is beside the truncation error; all of them are discussed in this section as follows.

Regarding this goal, suppose that the function values *f*(*x*), *f*(*x* **+** *h*), are quantized to *θ*<sup>1</sup> = *f*(*x* **+** *h*) + <sup>1</sup> , *θ*<sup>0</sup> = *f*(*x*) + <sup>0</sup> , with the sizes of the round-off errors <sup>1</sup> and <sup>0</sup> all being smaller than some positive number *ε*, that is |*<sup>j</sup>* | ≤ *ε*; with *j* = 0, 1.

Hence, the total error of the forward difference approximation defined by (41) is derived by

$$D\_f f(\mathbf{x}) = \frac{\theta\_1 - \theta\_0}{h} = \frac{f(\mathbf{x} + h) + \mathfrak{e}\_1 - f(\mathbf{x}) - \mathfrak{e}\_0}{h} = f'(\mathbf{x}) + \frac{\mathfrak{e}\_1 - \mathfrak{e}\_0}{h} + \frac{T\_f}{2}h. \tag{42}$$

Hence,

$$\left|D\_f f(\mathbf{x}) - f'(\mathbf{x})\right| \le \left|\frac{\boldsymbol{\varepsilon}\_1 - \boldsymbol{\varepsilon}\_0}{h}\right| + \left|\frac{T\_f}{2}\right| h \le \frac{2\varepsilon}{h} + \frac{|T\_f|}{2}h,\tag{43}$$

with *Tf* = *f* (*x*). Therefore, the upper bound of the error is illustrated by the right-hand side of Formula (43). The maximum limited of error contains two expressions; the first comes from the rounding error and in inverse proportion to step-size *h*, whilst the second comes from the truncation error and in direct proportion to *h*. These two parts can be |*Tf* |

formulated as a function *φ*(*h*) with respect to *h* as follows *φ*(*h*) = <sup>2</sup>*<sup>ε</sup> <sup>h</sup>* + <sup>2</sup> *h*. Now, if we find the minimizer *h*∗ of the function *φ*(*h*), then the value *φ*(*h*∗) is the upper bound of the total error. Hence *<sup>d</sup>φ*(*h*) *dh* <sup>=</sup> <sup>−</sup>2*<sup>ε</sup> <sup>h</sup>*<sup>2</sup> + |*Tf* | <sup>2</sup> = 0, then

$$h^\* = 2\sqrt{\frac{\varepsilon}{|T\_{f}|}} = 2\sqrt{\frac{\varepsilon}{|f''(\mathbf{x})|}}.\tag{44}$$

Therefore, it can be concluded that as we create small values of *h*, the round-off error might grow, whilst the truncation error reduces. It is called the "step-size dilemma".

Consequently, there have to be some optimal values of the *h*∗ for the forward difference approximation formula, as derived analytically in (44). However, Formula (44) is only of theoretical value and cannot be used practically to determine *h*∗ because we do not have any information about the second derivative and, therefore, we cannot estimate the values of *Tf* .

Therefore, there are many approaches which have been presented to deal with the step-size dilemma.

Recently, Shi et al. [66] proposed a bisection search for finding a finite-difference interval for a finite-difference method. Their approach was presented to balance the truncation error that grows from the error in the Taylor formula and the measurement error obtained from noise in the function evaluation. According to their numerical experience, the finitedifference interval *<sup>h</sup>*<sup>∗</sup> are bounded between the following ranges [<sup>2</sup> <sup>×</sup> <sup>10</sup>−4, 6.32 <sup>×</sup> <sup>10</sup>−1], [2.72 <sup>×</sup> <sup>10</sup>−4, 8.26 <sup>×</sup> 100] and [8.44 <sup>×</sup> <sup>10</sup>−3, 3.94 <sup>×</sup> 100] by using the forward and central differences to estimate the values of the first derivative of the *f* .

Additionally, the authors of [68] gave a study of the theoretical and practical comparison of the approximate values of the gradient vector in derivative-free optimization. These authors analyzed some approaches for approximating gradients of noisy functions utilizing only function values; those techniques include a finite difference.

The values of the finite difference interval are as follows 10−<sup>8</sup> <sup>≤</sup> *<sup>h</sup>*<sup>∗</sup> <sup>≤</sup> 1.

According to the earlier investigations, the core of the difference between all approaches is to determine the step size *h*. Hence, the value of the step size is ranged between this range *<sup>h</sup>*<sup>∗</sup> <sup>∈</sup> [1, 12 <sup>×</sup> <sup>10</sup>−10].

In this paper, the *h* is designed in a way that makes its values generated randomly. Additionally, the values of the *h* are connected to the function values per iteration to cover this domain, thus the feature here is that the value of *h* is modified per iteration randomly.

Therefore, a fresh approach to define the *h*∗ is presented in the following section.

## *3.2. Selecting a Step-Size h*

The forward difference approach is a cheap method compared to the different techniques.

The forward difference approach has shown promising results for minimizing noisy black-box functions [66].

Depending on the hypotheses which are listed in Section 2, let *x***<sup>0</sup>** be any starting point, thus function *f* satisfies the following *f*<sup>0</sup> ≥ *f*<sup>1</sup> ≥ ... ≥ *f <sup>k</sup>* , for *k* = 0, 1, 2, .... The numerical outcomes that are given in the past papers denote that the values of step-size *h* belong to the following range [10−10, <sup>≤</sup> <sup>1</sup>].

Therefore, the next Algorithm 2 is created to generate the values of the *h*∗ randomly from the intervals [0.1, 10−8].

## **Algorithm 2** Algorithm for calculating the values of *h*∗.

**Step 1:** At each iteration *k*, we generate a set random values between 10−2, and 10−7, and this set of random values is denoted by *L* = {*l*<sup>1</sup> , *l*2 ,..., *l*<sup>10</sup> }. **Step 2:** The minimum and maximum of the set *L* are extracted, respectively, as follows

*M* = min{*l i* :*i*=1,2,...,10}, *N* = max{*l i* :*i*=1,2,...,10} and set *Mf* <sup>=</sup> *<sup>M</sup>*−<sup>1</sup> . **Step 3:** The function value *f* is calculated at each *k*; *f <sup>k</sup>* = *f*(*xk* ).

Now we determine two cases according to the function values of the | *f <sup>k</sup>* | as follows. **Case 1:** If | *f <sup>k</sup>* | ∈ [10−1, <sup>∞</sup>), the value of the *<sup>h</sup>* is determined by

$$\mathcal{M}\_k = \begin{cases} \sqrt{\frac{N\_\epsilon}{M\_f}} & \text{if } |f\_k| > M\_{f'}\\ \sqrt{\frac{M\_\epsilon}{|f\_f|}} & \text{otherwise.} \end{cases} \tag{45}$$

**Case 2:** If | *f <sup>k</sup>* | ∈ [0, 10−1), the value of the *<sup>h</sup>* is determined by a random way from the range [10−4, 10−8].

**Example:** In this example, we show how the above algorithm is run.

Let us suppose that the point *x*<sup>0</sup> has four different values as starting points with four different values of *<sup>f</sup>* , for example, *<sup>f</sup>*<sup>0</sup> <sup>=</sup> *<sup>f</sup>*(*x*<sup>0</sup> ) = {1010, 106, 103, 10−1} and suppose we generate the set *<sup>L</sup>* as random values between 10−1, and 10−<sup>7</sup> such that *<sup>L</sup>* <sup>=</sup> {1.50 <sup>×</sup> <sup>10</sup><sup>−</sup>4, 5.10 <sup>×</sup> <sup>10</sup><sup>−</sup>6, 1.01 <sup>×</sup> <sup>10</sup><sup>−</sup>6, 1.40 <sup>×</sup> <sup>10</sup><sup>−</sup>2, 1.78 <sup>×</sup> <sup>10</sup><sup>−</sup>7, 1.92 <sup>×</sup> <sup>10</sup><sup>−</sup>5, 1.09 <sup>×</sup> <sup>10</sup><sup>−</sup>3, 2.77 <sup>×</sup> <sup>10</sup>−4, 2.99 <sup>×</sup> <sup>10</sup>−04, 5.15 <sup>×</sup> <sup>10</sup>−4}, *<sup>M</sup>* <sup>=</sup> 1.78 <sup>×</sup> <sup>10</sup>−7; hence, *Mf* <sup>=</sup> 5.618 <sup>×</sup> 106, since *<sup>f</sup>*<sup>0</sup> <sup>=</sup> 1010 <sup>&</sup>gt; *Mf* <sup>=</sup> 5.618 <sup>×</sup> <sup>10</sup>6, then we set *<sup>F</sup>*<sup>0</sup> <sup>=</sup> *Mf* <sup>=</sup> 5.618 <sup>×</sup> 106 and *<sup>h</sup>*<sup>1</sup> <sup>=</sup> <sup>2</sup> 6 *<sup>M</sup> Mf* = 2 71.78×10−<sup>7</sup> 5.618×<sup>106</sup> <sup>=</sup> 3.56 <sup>×</sup> <sup>10</sup>−7. If *<sup>f</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup>6, *<sup>f</sup>*<sup>0</sup> <sup>=</sup> <sup>106</sup> <sup>&</sup>lt; *Mf* <sup>=</sup> 5.618 <sup>×</sup> 106, and then *h*<sup>1</sup> = 2 6 *<sup>M</sup> F*0 = 2 75.618×<sup>106</sup> <sup>106</sup> <sup>=</sup> 8.438 <sup>×</sup> <sup>10</sup>−7, and *<sup>f</sup>*<sup>0</sup> <sup>=</sup> {10<sup>3</sup> <sup>&</sup>lt; *Mf* 5.618 <sup>×</sup> 106, we set *F*<sup>0</sup> = 103, then *h*<sup>1</sup> = 2 6 *<sup>M</sup> F*0 = 2 75.618×<sup>106</sup> <sup>103</sup> <sup>=</sup> 2.6683 <sup>×</sup> <sup>10</sup>−5. Finally, if *f*<sup>0</sup> = 10−1, then *h*<sup>1</sup> = 2 75.618×<sup>106</sup> <sup>10</sup>−<sup>1</sup> <sup>=</sup> 2.67 <sup>×</sup> <sup>10</sup>−3. The above example shows how Case 1 is implemented by using Formula (45).

Regarding Case 2 when 0 ≤ | *f <sup>k</sup>* | < 0.1, the value of the *hk* is taken randomly from the range [10−4, 10−8].

## *3.3. Estimating Gradient Vector*

The forward finite difference (DFF) is utilized to compute the approximate value of the gradient vector of function *<sup>f</sup>* at *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* by

$$[DFF]\_i = \frac{f(\mathbf{x} + h\mathbf{e}\_i) - f(\mathbf{x})}{h}, \text{for } i = 1, 2, \dots, n. \tag{46}$$

where *<sup>h</sup>* <sup>&</sup>gt; 0 is the finite difference interval defined in Section 3.2, and *ei* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is the *<sup>i</sup> th* column of the identity matrix.

Therefore, *g***(***x***)** ≈ *DFF*(*x*), is the approximate value of the gradient vector of function *f* at point *x*.

Therefore, the step size *ϕ<sup>k</sup>* is defined in the following.

The function *f*(*x*) is estimated by utilizing Taylor's expansion up to the linear term around the point *xk* , for each iteration *k*. Then we have

$$f(\mathbf{x}\_k + \mathbf{p}) \approx f(\mathbf{x}\_k) + \mathbf{g}(\mathbf{x}\_k)^T \boldsymbol{p}.$$

We define the quadratic model of *f*(*x*) at *xk* as

$$m\_k(p) = \frac{1}{2} \left( f(\mathbf{x}\_k) + \mathbf{g}(\mathbf{x}\_k)^T p \right)^2 = \frac{1}{2} f(\mathbf{x}\_k)^2 + f(\mathbf{x}\_k) \mathbf{g}(\mathbf{x}\_k)^T p + \frac{1}{2} p^T \mathbf{g}(\mathbf{x}\_k) \mathbf{g}(\mathbf{x}\_k)^T p.$$

Set *p* = −*ϕg***(***xk* **)** where *ϕ* is the step size along the −*g***(***xk* **)**. The optimal value of the *ϕ* is picked by solving the following subproblem: min *<sup>ϕ</sup>*∈<sup>R</sup> *mk*(*ϕ*) = <sup>1</sup> <sup>2</sup> *<sup>f</sup>*(*xk* )<sup>2</sup> <sup>−</sup> *<sup>ϕ</sup>f*(*xk* )*g***(***xi***)***Tg***(***xk* **)** <sup>+</sup> 1 <sup>2</sup> *<sup>ϕ</sup>*2(*g***(***xk* **)** *<sup>T</sup>g***(***xk* **)**)2. This gives

$$\varphi\_k = \frac{f(\mathbf{x}\_k)}{||\!|g(\mathbf{x}\_k)||^2}. \tag{47}$$

Therefore,

$$\left\|\mathbf{g(x\_k)}\right\|^2 = \frac{f(x\_k)}{\varphi\_k}, \varphi\_k \neq 0,\tag{48}$$

where *g***(***xk* **)** ≈ *DFF*(*xk* ).

## *3.4. Convergence Analysis of DFF*

The condition which is usually utilized in the convergence analysis of first-order methods with inexact gradient (DFF) vectors is defined by

$$||DFF(\mathbf{x}) - \mathbf{g}(\mathbf{x})|| \le C||\mathbf{g}(\mathbf{x})||\_{\prime} \tag{49}$$

for some 0 ≤ *C* < 1. This condition is introduced by [74,75] and it is called a norm condition. This condition denotes that the *g***(***x***)** ≈ *DFF*(*x*) is a descent direction for the function *f* [68].

However, condition (49) cannot be applied, unless we know *g***(***x***)**; therefore, this condition might be hard or impossible to verify.

There are many authors who have attempted to deal with this issue; see, for example, Refs. [68,76–79]. Byrd et al. [76] suggested a practical approach to estimate *g*(*xk* ), and they utilized it to guarantee some approximation of (49). Cartis and Scheinberg [77] and Paquette and Scheinberg [79] replaced condition (49) by

$$\|DFF(\mathbf{x}) - \mathbf{g(x)}\| \le k a\_k \|\|\mathbf{g(x)}\|\|,\tag{50}$$

where *k* > 0, and convergence rate analysis were derived for a line search method that has access to deterministic function values in [77] and stochastic function values (with additional assumptions) in [79]. Berahas et al. [68] established conditions under which (49) holds. For the forward finite differences method (DFF), they set *h*∗ = 2 7 *<sup>M</sup><sup>ε</sup> L* .

Therefore, we present the following

**Theorem 2.** *Under Assumptions 1 and 2 of Section 2, let DFF***(***x***)** *denote the forward finite difference approximation to the gradient <sup>g</sup>***(***x***)***. Then, for all <sup>x</sup>* <sup>∈</sup> <sup>R</sup>*n, the following inequality is true:*

$$\left| \left| DFF(\mathbf{x}\_{k}) \right| \right|\_{\boldsymbol{\nu}} - \left| \| g(\mathbf{x}\_{k}) \| \_{\boldsymbol{\nu}} \right| \leq \left| f(\mathbf{x}\_{k})\_{h\_{i}} - f(\mathbf{x}\_{k}) \right| + \frac{f(\mathbf{x}\_{k})}{q\_{k}}, \boldsymbol{\varphi}\_{k} \neq \mathbf{0}, \tag{51}$$

*where the value of the ϕ<sup>k</sup> is estimated by (47). We know that X*<sup>∞</sup> *and X are the norm infinity and the 2-norm, respectively, and they are defined by*

$$\|\mathbf{X}\|\_{\infty} = \max\_{1 \le i \le n} |\mathbf{x}\_i|\_{\prime} \tag{52}$$

$$\|\mathbf{X}\| = \sqrt{\sum\_{i} \mathbf{x}\_{i}^{\;\;\prime}},\tag{53}$$

*and then*

$$\left\|\mathbf{X}\right\|\_{\alpha} = \max\_{1 \le i \le n} |\mathbf{x}\_i| \le \sqrt{\sum\_i \mathbf{x}\_i^2}. \tag{54}$$

According to (46) which defines the gradient approximation by forward differences, the vector of [*DFF*(*xk* )]*<sup>i</sup>* is described by [*DFF*(*xk* )]*<sup>i</sup>* = <sup>1</sup> *<sup>h</sup>* [ *f*(*xk* + *ei h*) − *f*(*xk* )]*<sup>i</sup>* , wherer *i* = 1, 2, . . . , *n*, then

$$\|\mathsf{DFF}(\mathbf{x}\_{k})\|\_{\infty} = \max\_{1 \le i \le n} \left| \left[ \frac{f(\mathbf{x}\_{k} + \boldsymbol{\varepsilon}\_{i} \boldsymbol{h}) - f(\mathbf{x}\_{k})}{h} \right]\_{i} \right| = \frac{1}{h} \max\_{1 \le i \le n} |[f(\mathbf{x}\_{k} + \boldsymbol{\varepsilon}\_{i} \boldsymbol{h}) - f(\mathbf{x}\_{k})]\_{i}|.$$

and therefore, the next inequality is true

$$\|DFF(\mathbf{x}\_{k})\|\_{\infty} = \frac{1}{h} \max\_{1 \le i \le n} |[f(\mathbf{x}\_{k} + \boldsymbol{\varepsilon}\_{i} \boldsymbol{h}) - f(\mathbf{x}\_{k})]\_{i}| \le |f(\mathbf{x}\_{k})\_{h\_{i}} - f(\mathbf{x}\_{k})|.\tag{55}$$

By using (48), (51), (54) and (55), we obtain # # #*DFF*(*xk* )<sup>∞</sup> − *g*(*xk* )<sup>∞</sup> # # # ≤ *DFF*(*xk* )<sup>∞</sup> <sup>+</sup> *g*(*xk* )<sup>∞</sup> ≤ |*f*(*xk* )*hi* <sup>−</sup> *<sup>f</sup>*(*xk* )<sup>|</sup> <sup>+</sup> *g*(*xk* )<sup>2</sup> <sup>=</sup> <sup>|</sup>*f*(*xk* )*hi* <sup>−</sup> *<sup>f</sup>*(*xk* )<sup>|</sup> <sup>+</sup> *<sup>f</sup>*(*xk* ) *ϕk* , *ϕ<sup>k</sup>* = 0. Therefore, the theorem holds.

#### **4. Numerical Experiments of Part I**

All experiments were run on a PC with Intel(R) Core(TM) i5-3230M CPU@2.60GHz 2.60 GHz with RAM 4.00 GB of memory on a Windows 10 operating system. The five methods were coded by utilizing MATLAB version 8.5.0.197613 (R2015a) and the machine epsilon was about 10−16.

The model optimization test problems are categorized into two types. The first type is the test problems that contain a convex function, while the second type include a nonconvex function. Both kinds of test problems are listed in Tables 1–8 such that the second type of the test problem is referred to by ∗. Columns 1–4 of Table 1 give the data of the test problems as follows: the abbreviation of the function *f* is given on Column 1, the number of variables *n* is listed on Column 2, the exact function value *f*(*x∗*) at the global point *x<sup>∗</sup>* is presented on Column 3, and the exact value of the norm of the gradient *g***(***x∗***)** vector is given by Column 4, where the mark "−" denotes that the value of the norm of the gradient *g***(***x∗*) for the convex function satisfies the stopping criterion *g***(***x∗***)** <sup>&</sup>lt; <sup>10</sup>−6. Columns 5–8 are as Columns 1–4.

The data in Table 1 are taken from [56].

The numerical results for the local minimizers of all test problems are listed in Tables 2–8. Columns 1–2 and 8–9 contain the abbreviation of the function *f* and the number of the variables *n*, respectively. Columns 3–7 contain the abbreviation of each algorithm of the five algorithm SHZ, MHZ, HZ, HS and FR, which present the number of worst iterations, number of worst function evaluations, number of best iterations, number of best function evaluations, average of time (CPU), average of the number of iterations and average of the number of function evaluations, respectively. Columns 10–14 are similar to Columns 3–7.

**Note 1:** It is worth noting that the full name for each test function is mentioned in Appendix A according to the reference in which the test problem is.

**Note 2:** F denotes that the algorithm has failed to find the local minimizer of the function *f* according to the stopping criteria of Algorithm 1 which are listed in Section 4.1 below.


**Table 1.** List of both kinds of test problems.

**Table 2.** The number of worst iterations.



**Table 3.** The number of worst function evaluations.

**Table 4.** The number of best iterations.



**Table 5.** The number of best function evaluations.

#### **Table 6.** The average of time.



**Table 7.** The average of number of iterations.

**Table 8.** The average of number of function evaluations.


The stopping criteria of Algorithm 1 are as follows.

## *4.1. Stopping Criteria of Algorithm 1*

Since this section focuses in finding a local minimizer of all test problems, the stopping criteria of Algorithm 1 can be defined as follows.

According to the discussions of the convergence analysis which are mentioned in the previous sections, the stopping criterion of Algorithm 1 is, if *g***(***xk* **)** ≤ *ε*<sup>1</sup> is satisfied, Algorithm <sup>1</sup> stops, where *<sup>ε</sup>*<sup>1</sup> <sup>∈</sup> [10−6, 10−8]. However, the exact value of the gradient vector is unknown since the value of the gradient vector is estimated by Formula (46); therefore, this condition is replaced by *DFFk* ≤ *<sup>ε</sup>*<sup>2</sup> or FEs <sup>=</sup> *<sup>n</sup>*104, i.e., if one of them

is met, Algorithm <sup>1</sup> stops, where *<sup>ε</sup>*<sup>2</sup> <sup>∈</sup> [10−7, 10−9], FEs denotes the maximum function evaluations and *n* is the number variables of the *f* .

In the following section, the performance profile is presented as an easy tool to compare the performance of our proposed method versus other methods in finding local minimizers of convex or non-convex functions regarding the worst and best numbers of iterations and function evaluations, the average of CPU time and the average of iterations and function evaluations, respectively.

## *4.2. Performance Profiles*

The performance profile is the best tool for testing the performance of the proposed algorithms [80–84].

In this paper, the five algorithms' performance evaluation standards are as follows: the worst and best numbers of iterations and function emulations, and the average of the CPU time, iterations and function emulations. They are abbreviated as itr.w, itr.be, FEs.w, FEs.be, time.a, itr.a and EFs.a, respectively. In the remainder of the paper, the set Fit will be used to denote the seven criteria; Fit = {itr.w, itr.be, FEs.w, FEs.be, time.a, itr.a, EFs.a}.

Therefore, the numerical outcomes are presented in the form of performance profiles, as depicted in [82]. The most important characteristic of the performance profiles is that they can be shown in one figure by plotting for the different solvers a cumulative distribution function *ρ<sup>s</sup>* (*τ*).

The performance ratio is defined by first setting *rp*,*<sup>s</sup>* <sup>=</sup> *tp*,*<sup>s</sup>* min{*tp*,*s*:*s*∈*S*} , where *<sup>p</sup>* <sup>∈</sup> *<sup>P</sup>*, *<sup>P</sup>* is a set of test problems, *S* is the set of solvers, and *tp*,*<sup>s</sup>* is the value obtained by solver *s* on test problem *p*.

Then, define *ρ<sup>s</sup>* (*τ*) = <sup>1</sup> |*P*| size{*p* ∈ *P* : *rp*,*<sup>s</sup>* ≤ *τ*}, where |*P*| is the number of test problems. The value of *ρ<sup>s</sup>* (1) is the probability that the solver will win over the remaining ones,

i.e., it will yield a value lower than the values of the remaining ones.

In the following, the performance profiles are utilized to evaluate the performance of the five methods: SHZ, MHZ, HZ, SH and FR.

Therefore, in this paper, the term *tp*,*<sup>s</sup>* indicates one element of the set Fit, |*P*| = 46 is the number of test problems. We have 46 unconstrained test problems, 14 of which include non-convex functions. The group of solvers *S* = {*SHZ*, *MHZ*, *HZ*, *SH*, *FR*} finds the local minimizers of the 46 test problems; therefore, the values of the *Fit* are taken from the results of the 46 test problems as follows.

Each solver *s* of the set *S* is run 51 times for each of the 46 problems; at each run, every element of the set Fit has owned its value. So, they are analyzed in the following.

$$r\_{p,s} = \begin{cases} \frac{\text{fit}\_{p,s}}{\min\{\text{fit}\_{p,s} \colon s \in \mathbb{S}\}} & \text{if the s pass to solve the p}\_{\prime} \\ \infty & \text{otherwise,} \end{cases} \tag{56}$$

where fit*p*,*s* is an element of the Fit for the test problem *p* by using the solver *s*.

**Note:** Formula (56) means that if the final result, obtained by a solver *s* ∈ *S*, satisfies Inequality (57), then the first branch of (56) is computed. Otherwise, we set *rp*,*<sup>s</sup>* = ∞.

$$\|DFF\_k\| \le \varepsilon\_{2^\prime} \tag{57}$$

where *<sup>ε</sup>*<sup>2</sup> <sup>∈</sup> [10−5, 10−9].

Therefore, the performance profile of solver s is defined as follows:

$$\delta(r\_{p,s}, \tau) = \begin{cases} 1 & \text{if } r\_{p,s} \le \tau\_{\prime} \\ 0 & \text{otherwise} \end{cases} \tag{58}$$

Therefore, the performance profile for solver s is then given by the following function:

$$\rho\_\*(\tau) = \frac{1}{|P|} \left\{ \sum\_{p \in P} \delta(r\_{p,s\_\tau} \tau) \right\}, \ \tau \ge 1. \tag{59}$$

As we mentioned above, |*P*| = 46 and *τ* ∈ [1, 60].

By definition of Fit*p*,*s*, *ρ<sup>s</sup>* (1) denotes the fraction of test problems for which solver *s* performs the best. In general, *ρ<sup>s</sup>* (*τ*) can be explained as the probability for solver *s* ∈ *S* that the performance ratio *rp*,*<sup>s</sup>* is within a factor *τ* of the best possible ratio. Additionally, the essential characteristic of performance profiles is that they present data on the proportional performance of numerous solvers [82,83].

The numerical outcomes of the five methods are analyzed by using the performance profiles as follows. Figures 1–4 show the performance profiles of the set solvers *S*, for each element of the set Fit, respectively.

The performance profile depicted on the left of Figure 1 (in the term itr.w) compares the five techniques for a set of the 46 test problems.

The SHZ method has the best performance for the 46 test problems; this means that our suggested approach is capable of finding a local minimizer to the 46 test problems as fast as, or faster than, the other four approaches.

For instance, if *τ* = 1, the SHZ technique is capable of finding the local minimizer for 65% of problems versus the 33%, 20%, 20% and 13% of a set of test problems solved by the MHZ, HS, FR and HZ methods, respectively.

In general, the term itr.w, *τ* = 60 displays that all test problems are solved by SHZ against 96% of test problems solved by the MHZ, HZ and FR methods respectively, while 93% of test problems are solved by the HS method. At *τ* ≥ 400, all test problems are solved by the MHZ, HZ and FR methods respectively, while 98% of test problems are solved by the HS.

The right graph of Figure 1 shows that the method SHZ is capable of finding the local minimum of all test problems regarding term FEs.w.

The rest of Figures 2–4 show that the SHZ algorithm is superior to the four algorithms regarding the rest of the terms of the set Fit.

Therefore, the SHZ technique includes the characteristics of efficiency, reliability and effectiveness in solving Problem (1) compared to the other four methods.

**Note:** The power of the SHZ technique comes from the fact that the SHZ method gains the features of the four methods MHZ, HZ and HS, as we mentioned in Section 2.

**Figure 1.** *Cont*.

**Figure 1.** Plotting the results of the terms itr.w and FEs.w for 5 algorithms.

**Figure 2.** Plotting the results of the terms itr.be and FEs.be for 5 algorithms.

Part II: Global Minimization Problem

It is worth mentioning that the final results of Part I for the second set of test problems contain some global minimizers at some runs for some non-convex functions. This means that the pure CG technique could not find the global minimizer of the second type of test problems for each run because it is a local method.

Therefore, to make this method capable of solving Problem (2) per run, the random technique is proposed and it is added to the CG approach to gain a new PS-CG hybrid technique that solves Problem (2). In many studies, the numerical outcomes indicated that the interbreed between a classical method and a random technique is very successful in overcoming the weakness of these methods. See [55–59].

**Figure 4.** Plotting the results of the terms itr.a and FEs.a for 5 algorithms.

Consequently, this part of the paper seeks to solve Problem (2).

Therefore, each method of the five CG methods mentioned in Part I is hybridized with the stochastic technique to obtain five algorithms to try to solve Problem (2).

In the next section, a stochastic technique is presented.

#### **5. Random Technique**

In this section, a new random parameter "SP" is presented. This stochastic technique contains three different formulas by which three different points are generated. This set of formulas is combined with the CG method to obtain a new algorithm that solves Problem (2).

*Random Parameters (SP Technique)*

**Step 1:** The first point is computed as follows, generate *Vk* ∼ [−1, 1] *<sup>n</sup>* is as a random vector, set *<sup>γ</sup><sup>k</sup>* <sup>=</sup> <sup>10</sup>*ψ<sup>k</sup>* , *<sup>ψ</sup><sup>k</sup>* <sup>∈</sup> [0.01, 1), where the interval [0.01, 1) is divided into Itr of fractions and at every iteration *k*, the parameter *ψ<sup>k</sup>* takes one value of the Itr and then computes *<sup>λ</sup><sup>k</sup>* <sup>=</sup> (1+*γ<sup>k</sup>* ) |*Vi* | *γk SVi* as a research direction with the step lengths, where *i* = 1, 2, ... , *n*, *n* is the of number variables, Itr is the number of iterations, and *SVi* denotes the signs of the *V* and is defined by

$$SV\_i = \begin{cases} -1 & \text{if } V\_i < 0, \\ 1 & \text{otherwise.} \end{cases} \tag{60}$$

Thus, a point is calculated as follows:

$$\mathbf{x}\_1 = \mathbf{x}\_{\text{ac}} + \boldsymbol{\lambda}\_{\text{k}'} \tag{61}$$

where *xac* is the best point obtained yet, and then we compute *f*<sup>1</sup> = *f*(*x*<sup>1</sup> ).

**Step 2:** The second point is defined by

$$
\mathfrak{x}\_2 = \mathfrak{x}\_{ac} + \eta\_k \mathfrak{B}\_{k'} \tag{62}
$$

where *Bk* = *ϕkdk* , *ϕ<sup>k</sup>* is defined by (47), *η<sup>k</sup>* ∈ (0, 2) is a random number, and the *dk* is defined by (23). Then, we compute *f*<sup>2</sup> = *f*(*x***<sup>2</sup>** ).

**Step 3:** This point is defined by

$$\mathbf{x}\_{\mathfrak{z}} = \mathbf{X}\_{w} + \frac{1}{2} \mathbf{D} \mathbf{x}\_{\mathfrak{z}} \tag{63}$$

where *Dx* <sup>=</sup> (1+*μ<sup>k</sup>* ) |*Vi* | −1 *<sup>μ</sup>k*+0.1 *<sup>S</sup>Vi* , *<sup>μ</sup><sup>k</sup>* <sup>=</sup> <sup>|</sup> *fac*<sup>|</sup> 2, *fac* is the function value at the point *xac* that has been accepted, and *Xw* is a stochastic variable picked from the feasible range of the objective function. This means that for *Xw* ∼ [*a*, *b*] *<sup>n</sup>*, *a* and *b* are the lower and upper bounds of the feasible range, respectively, and the random vector *V* with its signs *SVi* is defined by the first step.

Therefore, we calculate *f*<sup>3</sup> = *f*(*x***<sup>3</sup>** ).

For finding the global minimizer of a non-convex function, the above stochastic technique is used since Algorithm 1 is not capable of finding the global solution at each run. In other words, in some runs, Algorithm 1 fails to find the global solution to this function due to it sticking to a local point.

In the following example, we show how the SP algorithm is run.

**Example:** This example shows how the three steps of the SP algorithm are implemented. We use the first test problem of the list of the test problems that are listed in Appendix A. *R*<sup>2</sup> (*x*) = 100(*x*<sup>2</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> )<sup>2</sup> + (*x*<sup>1</sup> <sup>−</sup> <sup>1</sup>)2, to facilitate an explanation of the mechanism of using the Sp algorithm (Formulas (61)–(63)), we use the following easy information about the function *R*<sup>2</sup> (*x*), *n* = 2 is the number of the variables, *xac* = [2; −1], or *xac* = [2; 1], where *xac* represents the best solution has been accepted so far or the starting point; hence, the function values at the two points are *<sup>R</sup>*<sup>2</sup> (*xac*) = <sup>100</sup>(22 <sup>+</sup> <sup>1</sup>)<sup>2</sup> + (<sup>2</sup> <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> <sup>2500</sup> <sup>+</sup> <sup>1</sup> <sup>=</sup> <sup>2501</sup> and *<sup>R</sup>*<sup>2</sup> (*xac*) = <sup>100</sup>(2<sup>2</sup> <sup>−</sup> <sup>1</sup>)<sup>2</sup> + (<sup>2</sup> <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> <sup>900</sup> <sup>+</sup> <sup>1</sup> <sup>=</sup> 901.

Supposing *Itr* = 5 is the number of iterations, the interval [0.01; 1) is divided into five fractions with step size <sup>1</sup>−0.01 <sup>5</sup> = 0.198, and thus the set of this fractions is *A* = {0.01, 0.208, 0.406, 0.604, 0.802}, let *k* be 3 which means the algorithm is at the third iteration. Then, *<sup>ψ</sup>*<sup>3</sup> <sup>=</sup> 0.406, *<sup>γ</sup>*<sup>3</sup> <sup>=</sup> <sup>10</sup>*ψ*<sup>3</sup> <sup>=</sup> <sup>10</sup>0.406 <sup>=</sup> 2.5468. Let *<sup>V</sup>*<sup>3</sup> be [−0.5; 1], then *λ***<sup>3</sup>** = (1+2.5468)|*−***0.5**<sup>|</sup> 2.5468 × −1; (1+2.5468)|**1**<sup>|</sup> 2.5468 × 1 = <sup>−</sup> 1.8833 2.5468 ; 3.5468 2.5468 = <sup>−</sup> 0.73948; 1.3926 .

Therefore, the new solution is computed by Formula (61) as follows.

*x*<sup>1</sup> = *xac* + *λ*<sup>3</sup> = [2; −1]+[−0.73948; 1.3926]=[1.2605; 0.3926] or *x*<sup>1</sup> = *xac* + *λ*<sup>3</sup> = [2; 1]+[−0.73948; 1.3926]=[1.2605; 2.3926].

The function values at both points are as follows.

*<sup>R</sup>*<sup>2</sup> (*x*<sup>1</sup> ) = <sup>100</sup>(1.2605<sup>2</sup> <sup>−</sup> 0.3926)<sup>2</sup> + (1.2605 <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> 143.1 <sup>+</sup> 0.06786 <sup>=</sup> 143.17 or *<sup>R</sup>*<sup>2</sup> (*x*<sup>1</sup> ) = <sup>100</sup>(1.2605<sup>2</sup> <sup>−</sup> 2.3926)<sup>2</sup> + (1.2605 <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> 64.6 <sup>+</sup> 0.06786 <sup>=</sup> 64.668.

Therefore, *R*<sup>2</sup> (*x*<sup>1</sup> ) < *R*<sup>2</sup> (*xac*); this means the solution that is generated by Formula (61) reduces the function value.

In the following, we explain how the candidate solution is generated by Formula (62).

Let *<sup>M</sup>* <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup>−6. By using Formula (45), we obtain *<sup>h</sup>*<sup>3</sup> <sup>=</sup> 4.381 <sup>×</sup> <sup>10</sup>−<sup>5</sup> as the step size *h* (a random interval) to the difference approximations method, and then we have *xh*<sup>1</sup> = [*xac*(1) + *<sup>h</sup>*<sup>3</sup> ; *xac*(2)] = [<sup>2</sup> <sup>+</sup> 4.381 <sup>×</sup> <sup>10</sup>−5; <sup>−</sup>1], *xh*<sup>2</sup> = [*xac*(1); *xac*(2) + *<sup>h</sup>*<sup>3</sup> ] = [2; <sup>−</sup><sup>1</sup> <sup>+</sup> 4.381 <sup>×</sup> <sup>10</sup>−5].

Therefore, the values of the function at the three points *xac*, *xh*<sup>1</sup> and *xh*<sup>2</sup> are listed in the following.

*<sup>R</sup>*<sup>2</sup> (*xac*) = 2501, *<sup>R</sup>*<sup>2</sup> (*xh*<sup>1</sup> ) = 2501.175 and *<sup>R</sup>*<sup>2</sup> (*xh*<sup>2</sup> ) = 2500.956.

We compute the approximate value of the gradient vector by Formula (46) as follows:

$$DFF(\mathbf{x}\_{\rm ac}) = \left[\frac{2501.175 - 2501}{4.381 \times 10^{-5}}; \frac{2500.956 - 2501}{4.381 \times 10^{-5}}\right] = \left[3994.522; -1004.34\right].$$

*ϕ*<sup>3</sup> = <sup>2501</sup> *DFF*<sup>2</sup> <sup>=</sup> 0.0002, where *<sup>ϕ</sup>*<sup>3</sup> is defined by (47).

We consider *d*<sup>3</sup> = −*g*(*xac*) ≈ − 3994.522; 1004.34 because we do not have information about the value of the *d*<sup>2</sup> in this illustration example.

Now, we apply Formula (62), as follows *B*<sup>3</sup> = *ϕ*<sup>3</sup> *d*<sup>2</sup> = [−0.799; 0.201], we take *η*<sup>3</sup> = 0.971 as a random number from the range (0, 2), then *x*<sup>2</sup> = [2; −1] + 0.971 × [−0.799; 0.201]=[1.2242; −0.80483], the function value at the point *x*<sup>2</sup> is *R*<sup>2</sup> (*x*<sup>2</sup> ) = 530.66. We note that the *R*<sup>2</sup> (*x*<sup>2</sup> ) = 530.66 < *R*<sup>2</sup> (*xac*) = 2501, i.e., the function value is reduced by the point *x*<sup>2</sup> .

In the following, we explain how the candidate solution is generated by Formula (63). *μ*<sup>3</sup> = | *fac*| <sup>2</sup> <sup>=</sup> 25012 <sup>=</sup> 6,255,001, *Dx* <sup>=</sup> (1+6,255,001)|−0.5<sup>|</sup> −1 6,255,001+0.1 × −1; (1+6,255,001)|1<sup>|</sup> −1 6,255,001+0.1 × 1 = <sup>−</sup> <sup>2501</sup>−<sup>1</sup> 6,255,001.1 ; 6,255,002−<sup>1</sup> 6,255,001.1 = − 0.0004; 0.999 . *Xw* = [−3.095; 8.701] is as a random vector picked from the range [−5, 10] 2, and then *<sup>x</sup>*<sup>3</sup> = [−3.095; 8.701] + <sup>1</sup> <sup>2</sup> [−0.0004; 0.999] = [−3.095; 8.701]+[−0.0002; 0.4995]=[−3.0952; 9.2005].

We compute the function value at the point *<sup>x</sup>*<sup>3</sup> ; *<sup>R</sup>*<sup>2</sup> (*x*<sup>3</sup> ) = <sup>100</sup>((−3.0952)<sup>2</sup> <sup>−</sup> 9.2005)<sup>2</sup> <sup>+</sup> (−3.0952 <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> 14.422 <sup>+</sup> 16.771 <sup>=</sup> 31.193.

We note that the *R*<sup>2</sup> (*x*<sup>3</sup> ) = 31.193 < *R*<sup>2</sup> (*xac*) = 2501. Therefore, the point *x*<sup>3</sup> minimizes the function value.

According to the above example that illustrates the mechanism of Formulas (61)–(63), we deduce the following results.

**Remark 1.** *Formulas (3), (61) and (62) are the main formulas which are used in the new hybrid proposed algorithm that is described in Section 6. However, Formula (63) is used when* Δ*f* = 0 *that is defined by Formula (25); in this case, Algorithm 3 reaches a critical point, thus if this point is the approximate value of the global minimizer point of the f , then Algorithm 3 stops according to the condition in Line 4 or Line 1 of Algorithm 3. Otherwise, the candidate solution is generated by Formula (63); see Section 6. Consequently, in this example, at iteration k* = 3*, the result which is obtained by Formula (63) cannot be taken into account due to the* Δ*f* = 0*.*

**Remark 2.** *All Formulas (61)–(63) minimize the function value from any starting point.*

## **6. Hybridization of the CG Method with Stochastic Parameters**

When a stochastic method as a global optimization algorithm is combined with a globally convergent method (deterministic method), the result is a global optimization algorithm [55,56].

Therefore, the SP technique is hybridized with each of the five conjugate gradient methods SHZ, MHZ, HZ, HS and FR to obtain five techniques.

Our proposed algorithm is called a hybrid stochastic CG method abbreviated by HSSHZ that solves Problem (2). However, Algorithm 3 represents five alternative algorithms when the SHZ method is hybridized with the PS technique, then we obtain a new algorithm abbreviated by HSSHZ. When we combine any method of MHZ, HZ, HS or FR, we obtain four other abbreviations of algorithms as follows: HSMHZ, HSHZ, HSHS and HSFR, respectively.

In general, the outputs of this paper are five algorithms that solve Problem (2), where the best one is the HSSHZ algorithm as illustrated by the numerical experiments section of Part II.

In the following, Algorithm 1 is combined with SP technique to obtain Algorithm 3.

The SP method permits conducting an exhaustive wipe of the search range to guarantee that the global minimizer point is visited at least once per run.

## **Algorithm 3** Hybrid stochastic CG method.

**Input:** *<sup>f</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>, *<sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*1, *fac* <sup>=</sup> *fcg* gained by Algorithm <sup>1</sup> and *<sup>ε</sup>* <sup>&</sup>gt; 0. **Output:** *xgl* = *xac* the global minimizer of *f* , *f*(*xgl*), the value of *f* at *xgl* . 1: **while** <sup>|</sup> *fac* <sup>−</sup> *<sup>f</sup>* ∗| <sup>&</sup>gt; *<sup>ε</sup>* or FEs<sup>&</sup>lt; *<sup>n</sup>*104 **do** 2: *fcg* is a function value *f* gained by Algorithm 1. 3: *fac* = min{ *fcg*, *f*<sup>1</sup> , *f*<sup>2</sup> } and *xac* the best point gives the *fac*. 4: **if** | *fac* − *f* ∗| ≤ *ε* **then** 5: Stop. 6: **end if** 7: **if** *f* == 0 **then** 8: calculate the *x***<sup>3</sup>** and the *f*<sup>3</sup> = *f*(*x***<sup>3</sup>** ) by Formula (63). 9: **if** *f*<sup>3</sup> < *fac* **then** 10: the *x***<sup>3</sup>** is accepted, compute the *xac* → *x***<sup>3</sup>** , *fac* → *f*<sup>3</sup> , and go to Line 1. 11: **else** 12: generate another point *x***<sup>3</sup>** by Formula (63). 13: **end if** 14: **else** 15: go to Line 1. 16: **end if** 17: **end while** 18: **return** *xac* the best point and its function value *fac*

## *A Mechanism Running Algorithm 3*

As we mentioned above, Algorithm 3 is a combination of two methods; the first is a CG method of the five techniques CG = {*SHZ*, *MHZ*, *HZ*, *SH*, *FR*} that are discussed in Part I, and the second is a random method is depicted by Section 5. The point *xcg* is obtained by Algorithm 1 and it will be an input to Algorithm 3.

Algorithm 3 begins with Line 1 that is the stopping standard of the algorithm. Therefore, Algorithm 3 ends if one of the following standards is satisfied: The first standard is <sup>|</sup> *fac* <sup>−</sup> *<sup>f</sup>* ∗| ≤ *<sup>ε</sup>*, and the second standard is FEs<sup>≥</sup> *<sup>n</sup>*104, where *fac* the best value of the function *f* is gained, the *f* <sup>∗</sup> is the true solution, *ε* = 10−6, FEs is the number of function evaluations, and FEs = *n*104 is a stopping standard indicated by [85,86].

In Line 3, the best value of *f* is selected from the three values of the function *fcg*, *f*<sup>1</sup> and *f*<sup>2</sup> , and indicated by *fac*, the three values of the function *f* are calculated by Algorithms (1), (61) and (62), respectively, and *xac* indicates this.

In Line 4, if | *fac* − *f* ∗| ≤ *ε* is fulfilled, the algorithm ends. The standard that is listed in Line 7 gives the algorithm an opportunity to flee from the local points. Consequently, if *f* = 0, then the algorithm has reached a crucial point. Therefore, if the norm of the gradient vector is 0 or ≈0, this point is either a local point or the global point. According to the above actions, the hybrid algorithm has been granted sequential opportunities to escape out of a snare (a local point). Thus, the procedures in Lines 8–12 are eligible for helping the

algorithm to flee this snare, especially since the second stopping standard guarantees that most of the research domain is scanned.

The numerical outcomes of the five methods are given in the next section.

#### **7. Numerical Experiments of Part II**

The numerical results for the second test problems (non-convex functions) are presented, and these results are obtained by Algorithm 3.

The performance profiles tool that is described in Part I is used here for assessing the achievement of Algorithm 3 that contains five alternatives of algorithms as we mentioned above in Section 6.

The numerical results of the second type of the test problems are listed in Tables 9–15. Columns 1–2 and 8–9 contain the abbreviation of the function *f* and the number of the unknowns *n*, respectively. Columns 3–7 contain the abbreviation of each algorithm of the five algorithm HSSHZ, HSMHZ, HSHZ, HSHS and HSFR, which present the number of worst iterations, number of worst function evaluations, number of best iterations, number of best function evaluations, average of time (CPU), average of number of iterations and average of number of function evaluations, respectively. Columns 10–14 are similar to Columns 3–7.

**Note:** F denotes that the algorithm has failed to find the local minimizer of the function *f* according to the stopping criteria of Algorithm 3 which are listed in Section 6.

**Table 9.** The number of worst iterations.


**Table 10.** The number of worst function evaluations.


**Table 11.** The number of best iterations.



**Table 12.** The number of best function evaluations.



**Table 14.** The average of number of iterations.




The performance profiles for the five algorithms are analyzed as follows.

Figures 5–8 show the performance profiles of the five set solvers *S* regarding the set standard Fit that is mentioned in Section 4.2.

The performance profiles which are drawn on the left of Figure 5 (in the term itr.w) compares 5 methods for the 14 test problems.

The HSSHZ technique has a good achievement (for the term itr.w) for all test problems, which indicates that the HSSHZ technique is capable of solving Problem (2) as fast as or faster than the four techniques.

For instance, if *τ* = 1, the HSSHZ algorithm solves 71% of the 14 test problems against 14%, 14%, 7% and 0%, of the 14 test problems solved by the HSMHZ, HSHZ, HSFR and HSHS algorithms, respectively.

In general, for the term itr.w, *τ* ≥ 60 exhibits that the second type of the test problems are solved by HSSHZ, while 64%, 71%, 43% and 50% of test problems are solved by the HSMHZ, HSHZ, HSHS and HSFR algorithms respectively.

Figures 5–8 demonstrate that the performance of the HSSHZ technique is better than the performance of the four techniques regarding the seven standards listed in the set Fit, respectively.

Therefore, the HSSHZ technique includes the characteristics of efficiency, reliability and effectiveness in finding the global minimizer of the non-convex function *f* compared to the other four methods.

It is worth observing that the power of the HSSHZ algorithm comes from the fact that the SHZ method gains the features of the four methods, MHZ, HZ, HS and FR, as mentioned in Section 2.

**Note 1:** In Algorithm 3, a run is considered successful if Inequality (64) is met.

$$|f\_{ac} - f^\*| \le 10^{-5},\tag{64}$$

where *f* ∗ is the exact global solution that is listed in Columns 3 and 7 of Table 1, respectively, and the *fac* is the final result obtained by Algorithm 3.

**Note 2:** Formula (56) means if the final result *fac* , obtained by Algorithm 3 satisfies Inequality (64), then the first branch of (56) is computed; otherwise, we set *rp*,*<sup>s</sup>* = ∞.

**Figure 5.** Plotting the results of the terms itr.w and FEs.w for 5 algorithms.

**Figure 6.** Plotting the results of the terms itr.be and FEs.be for 5 algorithms.

**Figure 7.** Plotting the results of the term time.a "CPU" for 5 algorithms.

**Figure 8.** Plotting the results of the terms itr.a and FEs.a for 5 algorithms.

## **8. Conclusions and Future Work**

A new modified CG algorithm is suggested, named SHZ. The SHZ finds the local minimizers of unconstrained optimization problems. The modernized formulae of the SHZ algorithm are more complicated than previous approaches; nevertheless, the numerical experiments of the SHZ are very strong. The convergence analysis of the SHZ algorithm is designed. We also analyzed the gradient approximation *g***(***x***)** ≈ DFF constructed by finite differences (the forward differences method). This method includes a new approach for selecting the fit value of the *h* according to the value of the objective function and it is updated dynamically at each iteration. The numerical results demonstrate that the performance of the SHZ method is positively competitive with the other four conjugate gradient methods based on performance profiles.

Comparing the final results of the gradient vector that were obtained by the method DFF to the exact values of the gradient vector demonstrates that the fresh technique succeeded in picking the right value of *h*. The proposed random approach recreates a critical role to make the SHZ method capable of finding the global minimizers of unconstrained optimization test problems, especially when the objective function is non-convex.

It can be worth observing that the power of the HSSHZ algorithm comes from the fact that the SHZ method gains the characteristics of the four methods, MHZ, HZ, HS and FR.

The suggested approach can be improved and modified to deal with constrained, multi-objective optimization problems, and it will be used for image restorations.

**Author Contributions:** Conceptualization, K.A.A.; Data curation, A.M.A.; Formal analysis, A.M.A. and K.M.S.; Funding acquisition, K.A.A.; Investigation, A.W.M.; Methodology, S.M.; Project administration, K.A.A. and K.M.S.; Software, S.M.; Supervision, A.W.M.; Validation, A.M.A.; Writing original draft, S.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research is funded by Researchers Supporting Program at King Saud University (Project# RSP-2021/305).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors present their appreciation to King Saud University for funding the publication of this research through Researchers Supporting Program (RSP-2021/305), King Saud University, Riyadh, Saudi Arabia.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A. List of Test Problems**

1 *Rn*: Rosenbrock functions [57,87,88]

$$\min\_{\mathbf{x}} \left\{ \sum\_{i=1}^{n-1} \left[ 100(\mathbf{x}\_i^2 - \mathbf{x}\_{i+1})^2 + (\mathbf{x}\_i - 1)^2 \right] \right\}.$$

Range of starting points −5 < *xi* < 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 1, 1, . . . , 1 .

2 *Zn*: Zakharov functions [57,80,87,88]

$$\min\_{\mathbf{x}} \left\{ \sum\_{i=1}^{n} \mathbf{x}\_{i}^{2} + \left( \sum\_{i=1}^{n} 0.5i\mathbf{x}\_{i} \right)^{2} + \left( \sum 0.5i\mathbf{x}\_{i} \right)^{4} \right\}.$$

Range of starting points −5 < *xi* < 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 0, 0, . . . , 0 .


$$\min\_{x} \left\{ \sum\_{i=1}^{n} x\_i^2 \right\}.$$

Range of starting points −10 ≤ *xi* ≤ 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 0, 0, . . . , 0 . 5 Tr: Trid function [80]

$$\min \sum\_{i=1}^{n} (\mathbf{x}\_i - \mathbf{1})^2 - \sum\_{i=2}^{n} \mathbf{x}\_i \mathbf{x}\_{i-1}.$$

Range of starting points <sup>−</sup>*n*<sup>2</sup> <sup>&</sup>lt; *xi* <sup>&</sup>lt; *<sup>n</sup>*2, *<sup>i</sup>* <sup>=</sup> 1, 2, . . . , *<sup>n</sup>*. Global minima : *f*(*x*∗) = *<sup>n</sup>*(*n*+4)(*n*−1) <sup>6</sup> . at *x*<sup>∗</sup> = *i*(*n* + 1 − *i*) 6: Sum Squares function [90]

$$\min\_{x} \left\{ \sum\_{i=1}^{n} i x\_i^2 \right\}.$$

Range of starting points −100 < *xi* < 100, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 0, 0, . . . , 0 .

7 CV: Colville function [57,80,91]

$$\begin{split} & \min\_{\mathbf{x}} \left\{ 100 \left( \mathbf{x}\_1^2 - \mathbf{x}\_2 \right)^2 + \left( \mathbf{x}\_1 - 1 \right)^2 + \left( \mathbf{x}\_3 - 1 \right)^2 + 90 \left( \mathbf{x}\_3^2 - \mathbf{x}\_4 \right)^2 \right\} \\ & + 10.1 \left( \left( \mathbf{x}\_2 - 1 \right)^2 + \left( \mathbf{x}\_4 - 1 \right)^2 \right) + 19.8 \left( \mathbf{x}\_2 - 1 \right) \left( \mathbf{x}\_4 - 1 \right)^2 \right\}. \end{split}$$

Range of starting points −10 < *xi* < 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = (1, 1, 1, 1). 8 BR: Branin function [57,92,93] min*<sup>x</sup>* (*x*<sup>2</sup> <sup>−</sup> 5.1 <sup>4</sup>*π*<sup>2</sup> *<sup>x</sup>*<sup>2</sup> <sup>1</sup> + <sup>5</sup> *<sup>π</sup> <sup>x</sup>*<sup>1</sup> <sup>−</sup> <sup>6</sup>)<sup>2</sup> <sup>+</sup> <sup>10</sup>(<sup>1</sup> <sup>−</sup> <sup>1</sup> <sup>8</sup>*<sup>π</sup> cos*(*x*1)) + <sup>10</sup> . Range of starting points −5 < *xi* < 15, *i* = 1, 2. Only one global minima: *f*(*x*∗) = 0.397887. at *x*<sup>∗</sup> = {(−*π*, 12.275), (9.42478, 2.475),(*π*, 2.275)}. 9 DJ: De Joung function [57,87,88] min*<sup>x</sup> x*2 <sup>1</sup> + *<sup>x</sup>*<sup>2</sup> <sup>2</sup> + *<sup>x</sup>*<sup>2</sup> 3 . Range of starting points −5 < *xi* < 15, *i* = 1, 2, 3. Number of local minima: no local minima. Global minima: *f*(*x*∗) = 0 at *x*∗ = (0, 0, 0). 10 BO: Booth function [89] min*<sup>x</sup>* (*x*<sup>1</sup> <sup>+</sup> <sup>2</sup>*x*<sup>2</sup> <sup>−</sup> <sup>7</sup>)<sup>2</sup> + (2*x*<sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> <sup>−</sup> <sup>5</sup>)<sup>2</sup> . Range of starting points −10 < *xi* < 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 1, 3 . 11 Ma: Matyas function [90] min*<sup>x</sup>* 0.26 *x*2 <sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> 2 − 0.48*x*<sup>1</sup> *x*<sup>2</sup> . Range of starting points −10 < *xi* < 10, *i* = 1, 2, . . . , *n*. Global minima: *f*(*x*∗) = 0 at *x*∗ = 0, 0 . 12 Sm∗: Shekel functions [57,80,87,88,92–94] min*<sup>x</sup>* − *m* ∑ *j*=1 <sup>4</sup> ∑ *i*=1 *xi* <sup>−</sup> *Aij*<sup>2</sup> <sup>+</sup> *cj* −<sup>1</sup> . where *c* = 0.1[1, 2, 2, 4, 4, 6, 3, 7, 5, 5], *A* = ⎡ ⎢ ⎢ ⎣ 4.0 1.0 8.0 6.0 3.0 2.0 5.0 8.0 6.0 7.0 4.0 1.0 8.0 6.0 7.0 9.0 3.0 1.0 2.0 3.0 4.0 1.0 8.0 6.0 3.0 2.0 5.0 8.0 6.0 7.0 4.0 1.0 8.0 6.0 7.0 9.0 3.0 1.0 2.0 3.0 ⎤ ⎥ ⎥ ⎦ Range of starting points 0 < *xi* < 10, *i* = 1, . . . , *n*. Number of local minima: *m* local minima. Global minima: *f*(*x*∗)*n*,*<sup>m</sup>* = ⎧ ⎨ ⎩ −10.1532 , when *m* = 5, −10.4029 , when *m* = 7, −10.5364 , when *m* = 10. Global minima for three functions at *x*∗ = 4, 4, 4, 4 . 13 GP∗: Goldstein and Price function [57,80,87,88,92,94] *<sup>u</sup>*(*x*) = <sup>1</sup> + (*x*<sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> <sup>+</sup> <sup>1</sup>)2(<sup>19</sup> <sup>−</sup> <sup>14</sup>*x*<sup>1</sup> <sup>+</sup> <sup>3</sup>*x*<sup>2</sup> <sup>1</sup> <sup>−</sup> <sup>14</sup>*x*<sup>2</sup> <sup>+</sup> <sup>6</sup>*x*1*x*<sup>2</sup> <sup>+</sup> <sup>3</sup>*x*<sup>2</sup> 2) *<sup>v</sup>*(*x*) = <sup>30</sup> + (2*x*<sup>1</sup> <sup>−</sup> <sup>3</sup>*x*2)2(<sup>18</sup> <sup>−</sup> <sup>32</sup>*x*<sup>1</sup> <sup>+</sup> <sup>12</sup>*x*<sup>2</sup> <sup>1</sup> <sup>+</sup> <sup>48</sup>*x*<sup>2</sup> <sup>−</sup> <sup>36</sup>*x*1*x*<sup>2</sup> <sup>+</sup> <sup>27</sup>*x*<sup>2</sup>

$$\min\_{\mathbf{x}} \left\{ v(\mathbf{x}) u(\mathbf{x}) \right\}.$$

2).

Range of starting points −2 < *xi* < 2, *i* = 1, 2. Number of local minima: 4 local minima. Global minima: *f*(*x*∗) = 3 at *x*<sup>∗</sup> = (0, −1).

14 Ras∗: Rastrigin function [93]

$$\min\_{\mathbf{x}} \left\{ \mathbf{x}\_1^2 + \mathbf{x}\_2^2 - \cos(18\mathbf{x}\_1) - \cos(18\mathbf{x}\_2) \right\}.$$

Range of starting points −1 < *xi* < 1, *i* = 1, 2. Number of local minima: many local minima. Global minima: *f*(*x*∗) = −2 at *x*<sup>∗</sup> = (0, 0).

15 Bh1∗: Bohachevsky function [80]

$$\min\_{\mathbf{x}} \left\{ \mathbf{x}\_1^2 + 2\mathbf{x}\_2^2 - 0.3\cos(3\pi\mathbf{x}\_1) - 0.4\cos(4\pi\mathbf{x}\_2) + 0.7 \right\}.$$

Range of starting points −100 < *xi* < 100, *i* = 1, 2. Number of local minima: many local minima. Global minima: *f*(*x*∗) = 0 at *x*∗ = 0, 0 .

16 SH∗: Shubert function in [57,80,87,88,92]

$$\min\_{\mathbf{x}} \left\{ \left( \sum\_{i=1}^{5} i \cos \left( (i+1)\mathbf{x}\_1 + i \right) \right) \left( \sum\_{i=1}^{5} i \cos \left( (i+1)\mathbf{x}\_2 + i \right) \right) \right\}.$$

Range of starting points −5.12 < *xi* < 5.12, *i* = 1, 2. Number of local minima: 760 local minima. Global minima: *f*(*x*∗) = −186.7309 at 18 point different of *x*∗.

17 P8∗ Ref. [92]

$$\min\_{\pi} \left\{ \frac{\pi}{n} \left( k\_1 \sin(\pi y\_1)^2 + \sum\_{i=1}^{n-1} \left( y\_i - k\_2 \right)^2 \left[ 1 + k\_1 \sin(\pi y\_{i+1})^2 \right] + (y\_n - k\_2)^2 \right) \right\},$$

with *yi* = 1 + <sup>1</sup> 4 *xi* + 1 , *k*<sup>1</sup> = 10 and *k*<sup>2</sup> = 1. Range of starting points −10 ≤ *xi* ≤ 10, *i* = 1, 2, 3. Number of local minima: 5<sup>3</sup> local minima. Global minima: *f*(*x*∗) = 0 at *x*<sup>∗</sup> = (−1, −1, −1).

18 P16∗ Ref. [92]

$$\min\_{\mathbf{x}} k\_3 \left\{ \sin^2(\pi k\_4 \mathbf{x}\_1) + \sum\_{i=1}^{n-1} \left( \mathbf{x}\_i - k \mathbf{y} \right)^2 \left[ 1 + k\_6 \sin^2(\pi k\_4 \mathbf{x}\_{i+1}) \right] + (\mathbf{x}\_n - k \mathbf{y})^2 \left[ 1 + k\_6 \sin^2(\pi k\_7 \mathbf{x}\_n) \right] \right\},$$

$$\min\_{\mathbf{x}} k\_3 \left\{ \mathbf{x}\_i - \mathbf{x}\_i \cos \alpha\_i \mathbf{x}\_{i+1} + \mathbf{x}\_i \cos \alpha\_i \mathbf{x}\_{i+1} \right\},$$

where *k*<sup>3</sup> = 0.1,*k*<sup>4</sup> = 3, *k*<sup>5</sup> = 1, *k*<sup>6</sup> = 1,*k*<sup>7</sup> = 2. Range of starting points −5 ≤ *xi* ≤ 5, *i* = 1, .., *n*. Number of local minima: 15<sup>5</sup> local minima. Global minima: *f*(*x*∗) = 0 at *x*∗ = (1, 1, 1, 1, 1). 19 CB∗: Camel back in [80] and camel function in [93]

$$\min\_{\mathbf{x}} \left\{ 4\mathbf{x}\_1^2 - 2.1\mathbf{x}\_1^4 + \frac{1}{3}\mathbf{x}\_1^6 + \mathbf{x}\_1\mathbf{x}\_2 - 4\mathbf{x}\_2^2 + 4\mathbf{x}\_2^4 \right\}$$

Range of starting points −5 < *xi* < 5, *i* = 1, 2. Number of local minima: many local minima. Global minima: *f*(*x*∗) = −1.0316285 at *x*<sup>∗</sup> = {(0.089842, −0.71266), (−0.089842, 0.71266)}.

20 H3∗: Hartmann function [57,80,87,88,92–94]

$$\min\_{\mathbf{x}} \left\{ -\sum\_{i=1}^{4} c\_i \varepsilon x p \left( -\sum\_{j=1}^{3} a\_{ij} \left( x\_j - p\_{ij} \right)^2 \right) \right\}.$$

2 .

Range of starting points −1 < *xj* < 1, *j* = 1, 2, 3. Number of local minima: 4 local minima. Global minima: *f*(*x*∗) = −3.86278 at *x*<sup>∗</sup> = 0.114614, 0.555649, 0.852547 . 21 H6∗: Hartmann function [57,80,87,88,92–94]

$$\min\_{\mathbf{x}} \left\{ -\sum\_{i=1}^{4} c\_i \mathbf{x} \mathbf{x} p \left( -\sum\_{j=1}^{6} a\_{ij} \left( \mathbf{x}\_j - p\_{ij} \right)^2 \right) \right\}.$$

Range of starting points −1 < *xj* < 1, *j* = 1, 2, . . . , *n*. Number of local minima: 4 local minima.

Global minima: *f*(*x*∗) = −3.32237 at *x*<sup>∗</sup> = 0.201690, 0.150011, 0.476874, 0.275332,0.311652, 0.657300 .

22 HM∗: hump Function [57]

$$\min\_{\mathbf{x}} \left\{ 1.0316285 + 4\mathbf{x}\_1^2 - 2.1\mathbf{x}\_1^4 + \frac{1}{3}\mathbf{x}\_1^6 + \mathbf{x}\_1\mathbf{x}\_2 - 4\mathbf{x}\_2^2 + 4\mathbf{x}\_2^4 \right\}.$$

Range of starting points −5 < *xi* < 5, *i* = 1, 2. Number of local minima: 3 local minima.

Global minima: *f*(*x*∗) = 0 at *x*<sup>∗</sup> = {(0.0898, −0.7126),(−0.0898, 0.7126)}.

23 Le∗: Levy function [95]

$$\min\_{\pi} \left\{ \sin^2(\pi w\_1) + \sum\_{i=1}^{n-1} (w\_i - 1)^2 \left( 1 + 10 \sin^2(\pi w\_i + 1) \right) + (w\_n - 1)^2 \left[ 1 + \sin^2(2\pi w\_n) \right] \right\},$$

where *wi* = <sup>1</sup> + *xi*−<sup>1</sup> <sup>4</sup> , for *i* = 1, . . . , *n*. Range of starting points −10 < *xi* < 10, *i* = 1, 2, . . . , *n*. Number of local minima: many local minima. Global minima: *f*(*x*∗) = 0 at *x*∗ = 1, 1, . . . , 1 .

#### **References**

