The Hybrid BFGS-CG Method in Solving Unconstrained Optimization Problems

Abstract

In solving large scale problems, the quasi-Newton method is known as the most efficient method in solving unconstrained optimization problems. Hence, a new hybrid method, known as the BFGS-CG method, has been created based on these properties, combining the search direction between conjugate gradient methods and quasi-Newton methods. In comparison to standard BFGS methods and conjugate gradient methods, the BFGS-CG method shows significant improvement in the total number of iterations and CPU time required to solve large scale unconstrained optimization problems. We also prove that the hybrid method is globally convergent.

1. Introduction

The unconstrained optimization problem only requires the objective function as min_{𝑥∈𝑅𝑛}𝑓(𝑥), (1) where 𝑅^𝑛 is an 𝑛-dimensional Euclidean space and 𝑓∶𝑅^𝑛→𝑅 is continuously differentiable. The iterative methods are used to solve (1). On the 𝑖th iteration, an approximation point 𝑥_𝑖 and the (𝑖+1)th iteration are given by 𝑥_𝑖+1=𝑥_𝑖+𝛼_𝑖𝑑_𝑖, (2) where 𝑑_𝑖 denotes the search direction and 𝛼_𝑖 denotes the step size. The search direction must satisfy the relation 𝑔^𝑇_𝑖𝑑_𝑖<0, which guarantees that 𝑑_𝑖 is a descent direction of 𝑓(𝑥) at 𝑥_𝑖. The different choices of 𝑑_𝑖 and 𝛼_𝑖 yield the different convergence properties. Generally the first order condition ∇𝑓(𝑥_∗)=0 is used to check for local convergence to stationary point 𝑥_∗. There are many ways to calculate the search direction depending on the method used, such as the steepest descent method, conjugate gradient (CG) method, Newton-Raphson method, and quasi-Newton method.

The different choices of the step size ensure that the sequence of iterates 𝑥_𝑖 defined by (2) is globally convergent with some rates of convergence. There are two ways to determine the values of the step size, the exact line search, and the inexact line search. For the exact line search, 𝛼_𝑖 is calculated by using the formula 𝛼_𝑖=argmin_𝑎>0(𝑓(𝑥_𝑖+𝛼_𝑖𝑑_𝑖)). However, it is difficult and often impossible to find the value of step size in practical computation using the exact line search. Hence, the inexact line search is proposed by previous researchers like Armijo [1], Wolfe [2, 3], and Goldstein [4] to overcome the problem. Recently Shi proposed a new inexact line search rule similar to the Armijo line search and analysed the global converge [5]. Shi also claimed that among several well-known inexact line search procedures published by previous researchers, the Armijo line search rule is one of the most useful and the easiest to be implemented in computational calculations. The Armijo line search rule can be described as follows: ()(),𝛼Given𝑠>0,𝛽∈0,1,𝜎∈0,1_𝑖=max𝑠,𝑠𝛽,𝑠𝛽²𝑓𝑥,…suchthat_𝑖−𝑓𝑥_𝑖+𝛼_𝑖𝑑_𝑖≥−𝜎𝛼_𝑖𝑔^𝑇_𝑖𝑑_𝑖, (3)𝑖=0,1,2,…. Then, the sequence of {𝑥_𝑖}^∞_𝑖=0 is converged to the optimal point, 𝑥^∗, which minimises 𝑓 [6]. Hence, we will use the Armijo line search in this research associated with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method and the new hybrid method.

This paper is organised as follows. In Section 2, we elaborate the step size and search direction that are used in this research. Here, the BFGS method and CG method also will be presented. Then, the new hybrid method and convergence analysis will be discussed in Section 3. An explanation about the numerical results is provided in Section 4 and the paper ends with a short conclusion in Section 5.

2. The Search Direction

The different methods in solving unconstrained optimization problems depend on the calculation of search direction, 𝑑_𝑖 in (2). In this paper, we will focus on the CG method and quasi-Newton methods. The CG method is useful for finding the minimum value of functions or unconstrained optimization problems, which are introduced by [7]. The search direction of the CG method is 𝑑_𝑖=−𝑔_𝑖,𝑖=0,−𝑔_𝑖+𝛽_𝑖𝑑_𝑖−1,𝑖≥2, (4) where 𝑔_𝑖=∇𝑓(𝑥_𝑖) and 𝛽_𝑖 is known as the CG coefficient. There are many ways to calculate 𝛽_𝑖 and some well-known formulas are 𝛽^FR_𝑖=𝑔^𝑇_𝑖𝑔_𝑖𝑔_𝑖−1²,𝛽^PR_𝑖=𝑔^𝑇_𝑖𝑔_𝑖−𝑔_𝑖−1𝑔_𝑖−1²,𝛽^HS_𝑖=𝑔^𝑇_𝑖𝑔_𝑖−𝑔_𝑖−1𝑔_𝑖−𝑔_𝑖−1^𝑇𝑑_𝑖−1, (5) where 𝑔_𝑖 and 𝑔_𝑖−1 are gradients of 𝑓(𝑥) at points 𝑥_𝑖 and 𝑥_𝑖−1, respectively, while ‖⋅‖ is a norm of vectors and 𝑑_𝑖−1 is a search direction for the previous iteration. The above corresponding coefficients are known as Fletcher-Reeves (CG-FR) [7], Polak-Ribière (CG-PR) [8–11], and Hestenes-Stiefel (CG-HS) [12].

In quasi-Newton methods, the search direction is the solution of linear system 𝑑_𝑖=−𝐻_𝑖𝑔_𝑖, (6) where 𝐻_𝑖 is an approximation of Hessian. Initial matrix 𝐻₀ is chosen by the identity matrix, which subsequently updates by an update formula. There are a few update formulas that are widely used like Davidon-Fletcher-Powell (DFP), BFGS, and Broyden family formula. This research uses a BFGS formula in a classical algorithm and the new hybrid method. The update formula for BFGS is 𝐻_𝑖+1=𝐻_𝑖−𝐻_𝑖𝑠_𝑖𝑠^𝑇_𝑖𝐻_𝑖𝑠^𝑇_𝑖𝐻_𝑖𝑠_𝑖+𝑦_𝑖𝑦^𝑇_𝑖𝑠^𝑇_𝑖𝑦_𝑖, (7) with 𝑠_𝑖=𝑥_𝑖−𝑥_𝑖−1 and 𝑦_𝑖=𝑔_𝑖−𝑔_𝑖−1. The approximation that the Hessian must fulfil is 𝐻_𝑖+1𝑠_𝑖=𝑦_𝑖. (8) This condition is required to hold for the updated matrix 𝐻_𝑖+1. Note that it is only possible to fulfil the secant equation if 𝑠^𝑇_𝑖𝑦_𝑖>0, (9) which is known as the curvature condition.

3. The New Hybrid Method

The modification of the quasi-Newton method based on a hybrid method has already been introduced by previous researchers. One of the studies is a hybridization of the quasi-Newton and Gauss-Seidel methods, aimed at solving the system of linear equations in [13]. Luo et al. [14] suggest the new hybrid method, which can solve the system of nonlinear equations by combining the quasi-Newton method with chaos optimization. Han and Neumann [6] combine the quasi-Newton methods and Cauchy descent method to solve unconstrained optimization problems, which is recognised as the quasi-Newton-SD method.

Hence, the modification of the quasi-Newton method by previous researchers spawned the new idea of hybridizing the classical method to yield the new hybrid method. Hence, this study proposes a new hybrid search direction that combines the concept of search direction of the quasi-Newton and CG methods. It yields a new search direction of the hybrid method which is known as the BFGS-CG method. The search direction for the BFGS-CG method is 𝑑_𝑖=−𝐻_𝑖𝑔_𝑖,𝑖=0,−𝐻_𝑖𝑔_𝑖+𝜂−𝑔_𝑖+𝛽_𝑖𝑑_𝑖−1,𝑖≥1, (10) where 𝜂>0 and 𝛽_𝑖=(𝑔^𝑇_𝑖𝑔_𝑖−1/𝑔^𝑇_𝑖𝑑_𝑖−1).

Hence, the complete algorithms for the BFGS method, CG-HS, CG-PR, and CG-FR methods, and the BFGS-CG method will be arranged in Algorithms 1, 2, and 3, respectively.

Algorithm 1 (BFGS method). States the following.

Step??0. Given a starting point 𝑥₀ and 𝐻₀=𝐼_𝑛, choose values for 𝑠, 𝛽, and, 𝜎 and set 𝑖=1.

Step??1. Terminate if ‖𝑔(𝑥_𝑖+1)‖<10⁻⁶ or 𝑖≥10000.

Step??2. Calculate the search direction by (6).

Step??3. Calculate the step size 𝛼_𝑖 by (3).

Step??4. Compute the difference between 𝑠_𝑖=𝑥_𝑖−𝑥_𝑖−1 and 𝑦_𝑖=𝑔_𝑖−𝑔_𝑖−1.

Step??5. Update 𝐻_𝑖−1 by (7) to obtain 𝐻_𝑖.

Step??6. Set 𝑖=𝑖+1 and go to Step 1.

Algorithm 2 (CG-HS, CG-PR, and CG-FR). States the following.

Step??0. Given a starting point 𝑥₀, choose values for 𝑠, 𝛽, and 𝜎 and set 𝑖=1.

Step??1. Terminate if ‖𝑔(𝑥_𝑘+1)‖<10⁻⁶ or 𝑖≥10000.

Step??2. Calculate the search direction by (4) with respect to the coefficient of CG.

Step??3. Calculate the step size 𝛼_𝑖 by (3).

Step??4. Compute the difference between 𝑠_𝑖=𝑥_𝑖−𝑥_𝑖−1 and 𝑦_𝑖=𝑔_𝑖−𝑔_𝑖−1.

Step??5. Set 𝑖=𝑖+1 and go to Step 1.

Algorithm 3 (BFGS-CG method). States the following.

Step??0. Given a starting point 𝑥₀ and 𝐻₀=𝐼_𝑛, choose values for 𝑠, 𝛽, and 𝜎 and set 𝑖=1.

Step??1. Terminate if ‖𝑔(𝑥_𝑖+1)‖<10⁻⁶ or 𝑖≥10000.

Step??2. Calculate the search direction by (10).

Step??3. Calculate the step size 𝛼_𝑖 by (3).

Step??4. Compute the difference between 𝑠_𝑖=𝑥_𝑖−𝑥_𝑖−1 and 𝑦_𝑖=𝑔_𝑖−𝑔_𝑖−1.

Step??5. Update 𝐻_𝑖−1 by (7) to obtain 𝐻_𝑖.

Step??6. Set 𝑖=𝑖+1 and go to Step 1.

Based on Algorithms 1, 2, and 3 we assume that every search direction 𝑑_𝑖 satisfied the descent condition 𝑔^𝑇_𝑖𝑑_𝑖<0, (11) for all 𝑖≥0. If there exists a constant 𝑐₁>0 such that 𝑔^𝑇_𝑖𝑑_𝑖≤𝑐₁𝑔_𝑖² (12) for all 𝑖≥0, then the search directions satisfy the sufficient descent condition which can be proved in Theorem 6. Hence, we need to make a few assumptions based on the objective function.

Assumption 4. Consider the following.H1:the objective function 𝑓 is twice continuously differentiable.H2:the level set 𝐿 is convex. Moreover, positive constants 𝑐₁ and 𝑐₂ exist, satisfying 𝑐₁‖𝑧‖²≤𝑧^𝑇𝐹(𝑥)𝑧≤𝑐₂‖𝑧‖², (13)?for all 𝑧∈𝑅^𝑛 and 𝑥∈𝐿, where 𝐹(𝑥) is the Hessian matrix for 𝑓.H3:the Hessian matrix is Lipschitz continuous at the point 𝑥^∗; that is, there exists the positive constant 𝑐₃ satisfying 𝑔(𝑥)−𝑔𝑥^∗≤𝑐₃𝑥−𝑥^∗ (14)?for all 𝑥 in a neighbourhood of 𝑥^∗.

Theorem 5 (see [15, 16]). Let {𝐵_𝑖} be generated by the BFGS formula (8), where 𝐵₁ is symmetric and positive definite, and 𝑦^𝑇_𝑖𝑠_𝑖>0 for all 𝑖. Furthermore, assume that {𝑠_𝑖} and {𝑦_𝑖} are such that 𝑦_𝑖−𝐺_∗𝑠_𝑖𝑠_𝑖≤𝜀_𝑖 (15) for some symmetric and positive definite matrix 𝐺(𝑥^∗) and for some sequence {𝜀_𝑖} with the property ∑^∞_𝑖=1𝜀_𝑖<∞. Then lim_𝑖→∞𝐵_𝑖−𝐺_∗𝑑_𝑖𝑑_𝑖=0 (16) and the sequence {‖𝐵_𝑖‖}, {‖𝐵_𝑖⁻¹‖} are bound.

Theorem 6. Suppose that Assumption 4 and Theorem 5 hold. Then condition (12) holds for all 𝑖≥0.

Proof. From (12), we see that 𝑔^𝑇_𝑖𝑑_𝑖=−𝑔^𝑇_𝑖𝐵_𝑖⁻¹𝑔_𝑖+𝜂𝑔^𝑇_𝑖−𝑔_𝑖+𝑔^𝑇_𝑖𝑔_𝑖−1𝑔^𝑇_𝑖𝑑_𝑖−1𝑑_𝑖−1=−𝑔^𝑇_𝑖𝐵_𝑖⁻¹𝑔_𝑖+𝜂−𝑔^𝑇_𝑖𝑔_𝑖+𝑔^𝑇_𝑖𝑔_𝑖−1𝑔^𝑇_𝑖𝑑_𝑖−1𝑔^𝑇_𝑖𝑑_𝑖−1=−𝑔^𝑇_𝑖𝐵_𝑖⁻¹𝑔_𝑖+𝜂−𝑔^𝑇_𝑖𝑔_𝑖+𝑔^𝑇_𝑖𝑔_𝑖−1. (17) Based on Powell [17], 𝑔^𝑇_𝑖𝑔_𝑖−1≥𝜀‖𝑔_𝑖‖² with 𝜀=(0,1], and 𝑔^𝑇_𝑖𝑑_𝑖=−𝑔^𝑇_𝑖𝐵_𝑖⁻¹𝑔_𝑖+𝜂−𝑔_𝑖²+𝜀𝑔_𝑖²≤−𝜆_𝑖𝑔_𝑖²+−𝜂+𝜂𝜀𝑔_𝑖²≤𝑐₁𝑔_𝑖², (18) where 𝑐₁=−(𝜆_𝑖+𝜂−𝜂𝜀) which is bound away from zero. Hence, 𝑔^𝑇_𝑖𝑑_𝑖≤𝑐₁‖𝑔_𝑖‖² holds. The proof is completed.

Lemma 7. Under Assumption 4, positive constants 𝜛₁ and 𝜛₂ exist, such that for any 𝑥_𝑖 and any 𝑑_𝑖 with 𝑔^𝑇_𝑖𝑑_𝑖<0, the step size 𝑎_𝑖 produced by Algorithm 2 will satisfy either 𝑓𝑥_𝑖+𝛼_𝑖𝑑_𝑖−𝑓_𝑖≤−𝜛₁𝑔^𝑇_𝑖𝑑_𝑖²𝑑_𝑖² (19) or 𝑓𝑥_𝑖+𝛼_𝑖𝑑_𝑖−𝑓_𝑖≤𝜛₁𝑔^𝑇_𝑖𝑑_𝑖. (20)

Proof. Suppose that 𝑎_𝑖<1, which means that (3) failed for a step size 𝑎^′≤𝑎_𝑖/𝜏: 𝑓𝑥_𝑖+𝛼^′_𝑖𝑑_𝑖−𝑓𝑥_𝑖≤𝜛𝑎^′𝑔^𝑇_𝑖𝑑_𝑖. (21) Then, by using the mean value theorem, we obtain 𝑓𝑥_𝑖+1−𝑓𝑥_𝑖=𝑔^𝑇𝑥_𝑖+1−𝑥_𝑖, (22) where 𝑔=∇𝑓(𝑥), for some 𝑥∈(𝑥_𝑖,𝑥_𝑖+1). Now, by the Cauchy-Schwartz inequality, we get 𝑔^𝑇𝑥_𝑖+1−𝑥_𝑖=𝑔^𝑇𝑥_𝑖+1−𝑥_𝑖+𝑔−𝑔_𝑖^𝑇𝑥_𝑖+1−𝑥_𝑖=𝑔^𝑇𝑥_𝑖+1−𝑥_𝑖+𝑔−𝑔_𝑖𝑥_𝑖+1−𝑥_𝑖≤𝑔^𝑇𝑥_𝑖+1−𝑥_𝑖+𝜇𝑥_𝑖+1−𝑥_𝑖²≤𝑔^𝑇𝑎^′𝑑_𝑖+𝜇𝑎^′𝑑²≤𝑔^𝑇𝑎^′𝑑_𝑖+𝜇𝑎^′‖𝑑‖². (23) Thus, from H3 (𝜛−1)𝑎^′𝑔^𝑇_𝑖𝑑_𝑖<𝑎^′𝑔−𝑔_𝑖^𝑇𝑑_𝑖≤𝑀𝑎^′𝑑_𝑖², (24) which implies that 𝑎_𝑖≥𝜏𝑎^′>𝜏(1−𝜛)−𝑔^𝑇_𝑖𝑑_𝑖𝑀𝑎^′𝑑_𝑖². (25) Substituting this into (21), we have 𝑓𝑥_𝑖+𝛼^′_𝑖𝑑_𝑖−𝑓𝑥_𝑖≤𝑐₂−𝑔^𝑇_𝑖𝑑_𝑖𝑎^′𝑑_𝑖², (26) where 𝑐₂=𝜏(1−𝜛)/𝑀, which gives (19).

Theorem 8 (global convergence). Suppose that Assumption 4 and Theorem 5 hold. Then lim_𝑖→∞𝑔_𝑖²=0. (27)

Proof. Combining descent property (12) and Lemma 7 gives ^∞_𝑖=0𝑔_𝑖⁴𝑑_𝑖²<∞. (28) Hence, from Theorem 6, we can define that ‖𝑑_𝑖‖≤−𝑐‖𝑔_𝑖‖. Then, (28) will be simplified as ∑^∞_𝑖=0‖𝑔_𝑖‖²<∞. Therefore, the proof is completed.

4. Numerical Result

In this section, we use the test problem considered by Andrei [18], Michalewicz [19], and Moré et al. [20] in Table 1 to analyse the improvement of the BFGS-CG method compared with the BFGS method and CG method. Each of the test problems is tested with dimensions varying from 2 to 1,000 variables. This represents a total of 159 test problems. As suggested by [20], for each of the test problems, the initial point 𝑥₀ will further subtract from the minimum point. In doing so, this leads us to test the global convergence properties and the robustness of our method. For the Armijo line search, we use 𝑠=1, 𝛽=0.5, and 𝜎=0.1. The stopping criteria we use are ‖𝑔_𝑖‖≤10⁻⁶ and the number of iterations exceeds its limit, which is set to be 10,000. In our implementation, the numerical tests were performed on an Acer Aspire with a Windows 7 operating system and using Matlab 2012 languages.

Table 1 

A list of problem functions.

The performance results will be shown in Figures 1 and 2, respectively, using the performance profile introduced by Dolan and Moré [21]. The performance profile seeks to find how well the solvers perform relative to the other solvers on a set of problems. In general, 𝑃(𝜏) is the fraction of problems with performance ratio 𝜏; thus, a solver with high values of 𝑃(𝜏) or one that is located at the top right of the figure is preferable.

Figure 1 

Performance profile in a log10 scale based on iteration.

Figure 2 

Performance profile in a log10 scale based on CPU time.

Figures 1 and 2 show that the BFGS-CG method has the best performance since it can solve 99% of the test problems compared with the BFGS (84%), CG-HS (65%), CG-PR (80%), and CG-FR (75%) methods. Moreover, we can also say that the BFGS-CG is the fastest solver on approximately 68% of the test problems for iteration and 52% of CPU time.

5. Conclusion

We have presented a new hybrid method for solving unconstrained optimization problems. The numerical results for a broad class of test problems show that the BFGS-CG method is efficient and robust in solving the unconstrained optimization problem. We also note that, as the size and complexity of the problem increase, greater improvements could be realised by our BFGS-CG method. Our future research will be to try the BFGS-CG method with coefficients of CG like Fletcher-Reeves, Hestenes-Stiefel, and Polak-Ribiére.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This research was supported by Fundamental Research Grant Scheme (FRGS Vote no. 59256).