Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Chen, Jinghong
Master of Science - MSc
University of British Columbia
We conduct a study on an alternative method to Derivative-Free Optimization (DFO). It is a method by applying numerical differentiation techniques to the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. We apply different gradient approximation strategies on the BFGS, particularly a complex-step method to the gradient approximation. We conduct an experiment to test the feasibility of the approximation-based-BFGS methods. The experiment also explores the efficiency and reliability of BFGS algorithm with different strategies and different settings. Results find that the approximation-based-BFGS method is feasible alternative to DFO. With appropriate parameter settings, approximation-based-BFGS methods achieve high level accuracy relative to the original BFGS method. Among different numerical differentiation techniques, the approximation-based-BFGS method with complex-step achieves outstanding and stable performance in accuracy.
Using Complex Step NumericalDifferentiation in theBroyden-Fletcher-Goldfard-ShannoAlgorithmbyJinghong ChenA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)May 2023© Jinghong Chen, 2023The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Using Complex Step Numerical Differentiation in the Broy-den-Fletcher-Goldfard-Shanno Algorithmsubmitted by Jinghong Chen in partial fulfilment of the requirements ofthe degree of Master of ScienceWarren Hare, I. K. Barber Faculty of ScienceSupervisorPeter J. Simpson, I. K. Barber Faculty of ScienceSupervisory Committee MemberMohamed Shehata, I. K. Barber Faculty of ScienceSupervisory Committee MemberYang Cao, School of EngineeringUniversity ExamineriiAbstractWe conduct a study on an alternative method to Derivative-Free Opti-mization (DFO). It is a method by applying numerical differentiation tech-niques to the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. Weapply different gradient approximation strategies on the BFGS, particularlya complex-step method to the gradient approximation. We conduct an ex-periment to test the feasibility of the approximation-based-BFGS methods.The experiment also explores the efficiency and reliability of BFGS algo-rithm with different strategies and different settings.Results find that the approximation-based-BFGS method is feasible al-ternative to DFO. With appropriate parameter settings, approximation-based-BFGS methods achieve high level accuracy relative to the originalBFGS method. Among different numerical differentiation techniques, theapproximation-based-BFGS method with complex-step achieves outstand-ing and stable performance in accuracy.iiiLay SummaryThe Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an out-standing method for optimization problems with smooth objective function.However, the BFGS algorithm, as well as other derivative-based methods,cannot be applied directly to those problems without derivative information.We explore an alternative approach to overcome this weakness by applyingnumerical differentiation to BFGS algorithm. The experiment verifies it isa feasible approach which inherits the speed and efficiency of the originalBFGS method. The experimental result also suggests complex step approxi-mation is a good choice among different numerical differentiation techniques.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . 4Chapter 2: Methodology . . . . . . . . . . . . . . . . . . . . . . 52.1 The BFGS Method . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Gradient approximation using forward difference . . . . . . . 142.3 Gradient approximation using central difference . . . . . . . . 142.4 Gradient approximation using complex step . . . . . . . . . . 152.5 Approximation-based-BFGS methods . . . . . . . . . . . . . . 16Chapter 3: Experimental set up . . . . . . . . . . . . . . . . . . 173.1 Test Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2 Experimental parameters . . . . . . . . . . . . . . . . . . . . 193.2.1 Stopping tolerance . . . . . . . . . . . . . . . . . . . . 193.2.2 Optimal tolerance . . . . . . . . . . . . . . . . . . . . 193.2.3 Sampling radius . . . . . . . . . . . . . . . . . . . . . 203.3 Performance measurements . . . . . . . . . . . . . . . . . . . 203.4 Performance profiles . . . . . . . . . . . . . . . . . . . . . . . 213.5 Data profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . 21Chapter 4: Experimental results and discussions . . . . . . . . 24vTABLE OF CONTENTS4.1 Result on ε = 10−3 and h = 10−3 . . . . . . . . . . . . . . . . 244.2 Result on ε = 10−3 and h = 10−6 . . . . . . . . . . . . . . . . 254.3 Result on ε = 10−3 and h = 10−9 . . . . . . . . . . . . . . . . 264.4 Result on ε = 10−6 and h = 10−3 . . . . . . . . . . . . . . . . 274.5 Result on ε = 10−6 and h = 10−6 . . . . . . . . . . . . . . . . 284.6 Result on ε = 10−6 and h = 10−9 . . . . . . . . . . . . . . . . 294.7 Result on ε = 10−9 and h = 10−3 . . . . . . . . . . . . . . . . 304.8 Result on ε = 10−9 and h = 10−6 . . . . . . . . . . . . . . . . 314.9 Result on ε = 10−9 and h = 10−9 . . . . . . . . . . . . . . . . 324.10 Comprehensive result analysis . . . . . . . . . . . . . . . . . . 354.10.1 The BFGS solver . . . . . . . . . . . . . . . . . . . . . . 354.10.2 The BFGS-AG-Oh solver . . . . . . . . . . . . . . . . . 364.10.3 The BFGS-AG-Oh2 solver . . . . . . . . . . . . . . . . . 374.10.4 The BFGS-AG-CS solver . . . . . . . . . . . . . . . . . 38Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 405.1 Further research directions . . . . . . . . . . . . . . . . . . . . 42Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43viList of TablesTable 2.1 Run-time difference on function evaluations . . . . . . 6Table 2.2 Run-time difference on function evaluations . . . . . . 6Table 2.3 Notation definitions . . . . . . . . . . . . . . . . . . . 8Table 2.4 Solver set . . . . . . . . . . . . . . . . . . . . . . . . . 16Table 3.1 Test problems . . . . . . . . . . . . . . . . . . . . . . . 18Table 4.1 Stopping tolerance affect on the BFGS solver . . . . . . 36Table 4.2 Impact on accuracy of different sampling radius h . . . 38viiList of FiguresFigure 3.1 An example of performance profiles . . . . . . . . . . 22Figure 3.2 An example of data profiles . . . . . . . . . . . . . . . 23Figure 4.1 Performance profiles and data profiles at ε = 10−3and h = 10−3 . . . . . . . . . . . . . . . . . . . . . . . 25Figure 4.2 Performance profiles and data profiles at ε = 10−3and h = 10−6 . . . . . . . . . . . . . . . . . . . . . . . 26Figure 4.3 Performance profiles and data profiles at ε = 10−3and h = 10−9 . . . . . . . . . . . . . . . . . . . . . . . 27Figure 4.4 Performance profiles and data profiles at ε = 10−6and h = 10−3 . . . . . . . . . . . . . . . . . . . . . . . 28Figure 4.5 Performance profiles and data profiles at ε = 10−6and h = 10−6 . . . . . . . . . . . . . . . . . . . . . . . 29Figure 4.6 Performance profiles and data profiles at ε = 10−6and h = 10−9 . . . . . . . . . . . . . . . . . . . . . . . 30Figure 4.7 Performance profiles and data profiles at ε = 10−9and h = 10−3 . . . . . . . . . . . . . . . . . . . . . . . 31Figure 4.8 Full graph of performance profiles and data profilesat ε = 10−9 and h = 10−3 . . . . . . . . . . . . . . . . 32Figure 4.9 Performance profiles and data profiles at ε = 10−9and h = 10−6 . . . . . . . . . . . . . . . . . . . . . . . 33Figure 4.10 Full graph of performance profiles and data profilesat ε = 10−9 and h = 10−6 . . . . . . . . . . . . . . . . 33Figure 4.11 Performance profiles and data profiles at ε = 10−9and h = 10−9 . . . . . . . . . . . . . . . . . . . . . . . 34Figure 4.12 Full graph of performance profiles and data profilesat ε = 10−9 and h = 10−9 . . . . . . . . . . . . . . . . 35Figure 4.13 Performance profiles on setting h to 10−12 and 10−15 39viiiChapter 1IntroductionOptimization refers to the process of finding the best possible solution fora problem given a set of constraints. It involves maximizing or minimizinga particular objective function subject to certain constraints. Optimiza-tion can be applied to various fields such as engineering [MA04], finance[Vas12], system design [HLS+22], and engineering management [LTY12],among many others. Optimization is also becoming increasingly importantin machine learning, where it is used to optimize the performance of modelsand algorithms, as well as to improve the efficiency and accuracy of trainingalgorithms [SCZZ19]. Researchers found optimization methods are powerfultools that can help solve many complex problems and improve efficiency.Derivative-free optimization (DFO) addresses methods of optimizationwithout examining derivatives or other higher-order information, such asgradients and Hessians. The research field of DFO attracts more and moreattention as there is an increasing need of solving optimization problemswith no derivative information provided in their objective function. DFOis essential in real-world problems where the objective function is noisy,black-box, or computationally expensive to evaluate, making it impossibleor impractical to compute derivatives. DFO methods are widely used invarious fields such as engineering [HNT13], healthcare [Har10], and science[WSS15], where the optimization of complex systems is essential. A typicalscenario is when the objective function is provided by a computer simu-lation. An example application is optimization on improving earthquakeretrofits [AH17, pp. 8-9]. A common engineering solution to protect build-ings during an earthquake is using dampers to join two adjacent buildingstogether. Mathematical models via computer simulations have been builtto compute the cost of installing dampers and to quantify the effective-ness of their protection. Those mathematical models are generated throughcomputer simulation, so we do not have direct access to its derivatives orgradients. This type of optimization problem urges the need for good DFOmethods to be explored.The study of DFO covers a large scope. Our research attention focuseson unconstrained optimization of continuous, single-objective problems with1deterministic objective function. That is,argminf(x) subject to x ∈ Rn, where f is smooth. (1.1)Here, argmin means the ‘argument of the minimum’, that is, the point atwhich the minimum occurs. A function is smooth if it is two times contin-uously differentiable.Like other DFO research, the objective function f of the problem maybe differentiable, but the derivative of f may not available or may be toocomputationally expensive to obtain.There are two major categories of DFO methods [CSV09, Chapter 1].The first category is direct search methods. This kind of method main-tains a set of sampling points and evaluates them by the objective function.By comparing function values, these methods keep track of a series of trialpoints in the search space and using the results to guide the search towardsthe optimal solution. The Nelder-Mead simplex method [NM65] is a widelyused direct-search method. It keeps track of the vertices of a simplex inn-dimensional space and manipulates those vertices to search for the mini-mizer. Directional direct-search (DDS) [KLT03], such as coordinate search[Fer52], is another subtype of methods in this category. DDS methods gen-erate a set of poll directions and control the step-size parameters to reach alocal minimizer. Recent research like the generalized pattern-search method(GPS) by [AD02] introduces flexible and adaptive patterns to search forthe optimal solution. The mesh adaptive direct search (MADS) method[AD06] improves GPS methods by incorporating a mesh to work with thecases when the objective function is non-smooth. All these research lead tosuccessful solutions of general problems without derivative information.Another major category of DFO methods is model-based methods. Thiskind of method builds models to predict new points to minimize the objec-tive function. The models can be linear models [Reg15], quadratic interpo-lation models [CT96], regression models [CSV08], and so on. Some methods([WRS08] and [VKPS17]) in this category apply a trust-region to the model,which helps build a reasonable neighborhood to maintain the model.Model-based methods normally have a faster convergence in smoothproblems than direct search methods, since they take some advantage inusing models to detect the shape of the objective functions. General speak-ing, methods with higher order models have faster convergence (in termsof iteration) on smooth functions. However, higher order models also re-quire heavier computation. For example, to build a determined quadraticmodel in a n-dimensional problem, it needs to solve (n+1)(n+2)2 coefficients2and therefore needs (n+1)(n+2)2 function evaluations. Model-based methodsnormally take relatively heavy computation at each iteration, especially forhigh dimensional problems.Researchers keep making progress in these two categories of DFO meth-ods these days. However, traditional derivative-based methods, like theBroyden-Fletcher-Goldfarb-Shanno (BFGS) method, still have an advantagein solving problems with smooth objective function, since using derivativeinformation is a highly efficient way to detect the shape of the function.Nevertheless, they cannot be applied directly to problems whose derivativeare not analytically available.We remark that there are two alternatives to DFO method mentionedin [LMW19]. The first one is automated differentiation, which generatesderivatives using algorithms that automatically apply differentiation tech-niques within a simulation. The disadvantage of this approach is that itneeds to have the simulation coded in a language that supports automateddifferentiation. The second one is numerical differentiation, which appliesnumerical differentiation, such as central differences, to a derivative-basedoptimization method. Considering the feasibility of this alternative and thehigh performance of modern computers in calculation accuracy, we extendour interest in this approach. However, it is also worth noting that thequality of the derivative estimate may affect the accuracy and the perfor-mance of the optimization method. We therefore carry out an experimentand analysis the result in this thesis to explore the following questions.1. Is applying numerical differentiation within the BFGS, a derivative-based optimization method, a feasible alternative approach to DFOmethods in solving problems with a smooth objective function?2. What is the difference among those approximation-based-BFGS meth-ods in terms of function evaluations and run-time performance?3. There are some parameter settings used within numerical differentia-tion, such as sampling radius and stopping tolerance. How does thechoices of these parameters affect the performance of different solvers?We select the BFGS algorithm as the derivative-based method in ourexperiment. The BFGS algorithm is a popular method of a group of algo-rithms that are an extension to Newton’s method and referred as Quasi-Newton methods. The BFGS uses an approximation of the Hessian usedin Newton methods and has some improvement on the calculation, whichmakes it exceed other Quasi-Newton methods in run-time performance. The31.1. THESIS OVERVIEWBFGS is widely used in many application and incorporated in many softwarelibraries. It is particularly popular for solving problems when the Hessianis hard to compute [GNS09, pp.411].For numerical differentiation, we apply finite difference methods, suchas forward difference and central difference, in our experiment. Also, we areinterested in applying complex-step derivative approximation in our exper-iment, by the inspiration of [HS22]. The original work of applying complexvariables to develop approximation of derivatives started by Lyness andMoler in [LM67]. Following their theory, Squire and Trapp gave a simpleexpression for estimating the gradient of real functions in [ST98]. Fromtheir research results, we see that the complex-step approximation showssome advantages on accuracy performance by avoiding catastrophic cancel-lation, which is when numerical errors are caused by a near-zero subtractioncombined with a near-zero division. We will discuss more details about themethodology in Chapter 2.We conduct our experiments and implement all methods in Matlab.All source codes and test files are organized in the following repository,https://github.com/timtimchen/BFGS_Optimization.1.1 Thesis overviewThe remainder of this thesis is organized as follows. In Chapter 2, wedescribe the methods and algorithms that are used in our experiment. InChapter 3, we explain the details of the experimental set up, including in-formation on test problems collection, parameter settings in algorithms andanalysis techniques. In Chapter 4, we illustrate and interpret the experi-mental results, and discuss the significance of our findings. In Chapter 5,we present the conclusion of our research and a preview of future possiblestudy.4Chapter 2MethodologyIn this chapter, we will discuss the methods that are used in our exper-iments. First, we present some mathematical concepts and theorems thatsupport the methods we used in our experiments.Definition 2.1. Complex variable and complex function are defined as,z = x+ iy and f(z) = f(x+ iy) = u(x, y) + iv(x, y).Here, i is the imaginary unit, that is i2 = −1. Variable x is the realpart of z, and y is the imaginary part of z. Function u(x, y) is the realpart of f(z), and v(x, y) is the imaginary part of f(z), and they are single-valued real functions [AF03, Chapter 1]. Both x and y can be a real-valuednumber, or n-dimensional real-valued column vector, which is suitable to theobjective function. In our experiments, the size of optimization problemsare various, so we do not specify the dimension n in our discussion, and allmethods and formulas work the same way for dimension of any sizes.From Definition 2.1, we can see a complex function takes more steps toevaluate, comparing to a real-valued function with the same form. Therefore,in our experiment, we expect that it will take more CPU-time to evaluatean objective function at a complex-valued input than at a real-valued input.We use a simple example to illustrate this thought.Example 2.2. Consider the function f(x) = x2. To evaluate f at x = 5, itonly takes one multiplication operation. On the other hand, to evaluate atx = 5 + 3i, it takes four multiplications and one subtraction floating-pointoperations to finish, since the actual function form is f(5+3i) = (5+3i)2 =(52 − 32) + (2)(5)(3)i.Although the real running time would be different based on the functionand the running platform, computing in complex value takes longer thancomputing in real value is expected. Note that some functions, such asf(x) = ex, have some advantage in complex value evaluation, but most of thearithmetic calculation of complex value still takes more time than the realvalue. We carried out a small test to verify that the above situation applies5to the running platform we experimented on. The test takes a number ofiterations to run. In each iteration, it generate a real number and a complexnumber randomly. Then it pass the numbers to the test functions and recordthe evaluation time of each test function. The details of test functions aregiven in Table 2.1. The accumulated run-time result of this test is given inTable 2.2.Table 2.1: Run-time difference on function evaluationsfunction dimension implementation in Matlabf(x) = x2 1 x′ ∗ xf(x) =∑ni=1 x2i 10 sum(x.ˆ2)f(x) = ex 1 exp(x)Rosenbrock[Ros60] 2 100 ∗ (x(2)− x(1)ˆ2)ˆ2 + (1− x(1))ˆ2Table 2.2: Run-time difference on function evaluationsfunction evaluating real value evaluating complex value(second in 107 runs) (second in 107 runs)f(x) = x2 6.257 23.304f(x) =∑ni=1 x2i 7.269 15.286f(x) = ex 7.194 7.224Rosenbrock 5.854 8.237Definition 2.3. Big-O notation O() is an asymptotic notation that de-scribes an upper limit of a function. If there exists a constant numberM > 0 such that|f(x)| ≤Mg(x) as x→ 0,then we say f(x) = O(g(x)).Note that the limit of the function can also be x → ∞ or x → x0,whereas x → 0 is the case which is applied to our methods in this paper.Big-O notation is often used to denote an upper bound on the growth rateof a function or its remainder. For example, when we use O(h) to denote theremainder a function f(x + h), it means its remainder R(h) has order of hand satisfies R(h) < Mh as h→ 0 for someM > 0. In the following sectionsof this chapter, we also use Big-O notation in a column vector form, whichmeans each component of the multi-valued function satisfies the conditionin Definition 2.3.6The following are some mathematical facts of Big-O notation which willbe applied to develop our methods. More details about the definition ofBig-O notation and its properties can be found in [dB81, Section 1.2].Fact 2.4.O(h) = −O(h)O(hn)hm= O(hn−m)O((ih)n) = inO(hn)Theorem 2.5. Taylor’s Theorem (in Big-O notation form). If f is a (n+1)-times continuously differentiable function, thenf(x+ h) = f(x) +n∑k=1f (k)(x)hkk!+O(hn+1) as h→ 0.Details can be found in [Ken13, Chapter 31]. The polynomial appearingin Taylor’s theorem is referred as nth-order Taylor expansion.Example 2.6. Consider the function f(x) = e−x. The 3rd-order Taylorpolynomial at x isPx(h) = e−x − he−x + 12h2e−x − 16h3e−x.It related to f byf(x+ h) = Px(h) +O(h4).Note that Taylor’s Theorem also works for complex functions if the func-tion is complex analytic [BN12, Chapter 5]. For our research purpose, we aremainly interested in getting a second-order Taylor expansion of f(x + ih),where x is a real number and ih is a small step on the imaginary unitdirection with radius h. That is,f(x+ ih) = f(x) + ihf ′(x) +12(ih)2f ′′(x) +O((ih)3)= f(x) + ihf ′(x) +12i2h2f ′′(x) + i3O(h3)= f(x) + ihf ′(x)− 12h2f ′′(x) + iO(h3).(2.1)We also collect other mathematical notations used throughout this chap-ter and present them in Table 2.3.72.1. THE BFGS METHODTable 2.3: Notation definitionsNotation Meaningf objective function, f : Rn 7→ R∇f the gradient of the objective function (n-dimensional vector)∇2f the Hessian of the objective function (n× n matrix)I n× n identity matrixIm(f) the imaginary part of f∥A∥F the Frobenius norm of A, ∥A∥F =√∑ni,j=1 a2ij2.1 The BFGS MethodThe BFGS is a Quasi-Newton line-search method, which is based onNewton’s method in optimization, but using an approximation matrix in-stead of the true Hessian [NW99, pp. 23]. The BFGS algorithm is widelyrecognized as the gold-standard in the category of derivative-based optimiza-tion methods due to its efficiency, robustness, and scalability [ABSM14].Our description of the BFGS is mainly based on Chapter 4 of [Kel99].The descent direction of BFGS in each iteration is generated from a quadraticmodel mk(x), as follows,mk(x) = fk +∇fTk (x− xk) +12(x− xk)TBk(x− xk). (2.2)Here fk is the function value of the current point xk and ∇fk is the gradient.The matrix Bk is an approximation of the Hessian ∇2fk at xk and will beupdated at each iteration. Note that this model is an approximation of thesecond-order Taylor polynomial of the objective function f at the currentpoint xk. The only difference is we are using an approximation Bk insteadof the true Hessian ∇2fk. So mk(x) is also a good approximation of f nearxk. To find the minimizer x∗ of mk(x), we may take the derivative of mk(x)and set it to be equal to zero (in column vector form), that is,∇fk +Bk(x− xk) = 0, (2.3)and it yields,x∗ = xk −B−1∇fk. (2.4)This result gives us a good search direction to select the next point xk+1.Therefore, the BFGS set the descent direction dk as,dk = −B−1k ∇fk. (2.5)82.1. THE BFGS METHODThen it calculates the iterate xk+1 as,xk+1 = xk + αkdk, (2.6)where αk is the step size, and it is chosen to satisfy the Wolfe conditions,provided in Definition 2.7.Definition 2.7. Given algorithmic parameters c1, c2 such that 0 < c1 <c2 < 1, the Wolfe conditions require,f(xk + αkdk) ≤ f(xk) + c1αk∇fTk dk, and∇f(xk + αkdk)Tdk ≥ c2∇fTk dk.(2.7)Now we illustrate how to get a matrix Bk as an approximation of ∇2fk.As it is a n × n matrix, if we compute a fresh Bk at every iteration, itneeds to solve n(n+1)2 unknowns (since Bk is symmetry) and it is a heavycomputation. So it is proposed to update it in a manner to get a betterapproximation of ∇2fk iteratively. We have the first-order Taylor expansionof ∇f centered at xk as follow,∇f = ∇fk +∇2f(xk)(x− xk) +O((x− xk)2). (2.8)So we are seeking a matrix Bk+1 such that∇fk+1 −∇fk = Bk+1(xk+1 − xk), (2.9)and Bk+1 is symmetric and positive-definite. To simplify the formula, welet s = xk+1 − xk, and y = ∇fk+1 −∇fk, then we have,y = Bk+1s. (2.10)We then apply the least change principle to Equation (2.10), that is mini-mizing ∥B − Bk∥, which can be examined using weighted Frobenius norm,that is,∥B −Bk∥2F =n∑i,j=1wij(aij − bij)2, (2.11)where W = [wij ] and it satisfies W =WT , Ws = y. According to the proofin [NW99, pp. 137-139], the unique solution of B∗ yields,Bk+1 =(I − ysTyT s)Bk(I − syTyT s)+yyTyT s. (2.12)92.1. THE BFGS METHODFormula (2.12) is actually used in the iterative process of the Davidon-Fletcher-Powell (DFP) method. The BFGS method notices that it actuallyneeds matrix B−1 in the iterative process, instead of directly using thematrix B. To avoid computing inverse matrix in the algorithm, the BFGSimproves the process by letting Hk = B−1k and applies inverse operation toFormula (2.12) [NW99, pp. 140-141]. That yields,Hk+1 =(I − syTsT y)Hk(I − ysTsT y)+ssTsT y. (2.13)Now we have shown all the components of the BFGS method, and wewill illustrate its iterative process in Algorithms 1. Theoretically, the opti-mization process should converge to a point with zero gradient. Therefore,in the main process of the BFGS, we introduce a predefined small positivenumber ε as the stopping tolerance and test if the norm of current gradient isless than it. We will discuss more details of this stopping criteria in Section3.2.1. Algorithms 2 and 3 are subroutines of the BFGS algorithm, whichare supposed to find a proper step size α for current iteration that satisfiesthe Wolfe conditions. They are referred to the algorithms in Page 60 and61 of [NW99], and have been modified to fit for multi-valued functions. Theparameters c1, c2 in these subroutines should satisfy 0 < c1 < c2 < 1 andbe kept as the same value throughout the process. The values of them, aswell as other parameters such as m1 and amax, are empirical settings. Theymight affect on the efficiency of the solver slightly, but as we keep them thesame of all solvers in our experiments, it will not affect the comparison andour experimental results.Note that we mention a subroutine Grad() in Line 6 of Algorithm 1, Lines2 and 10 of Algorithm 2, and Lines 2 and 9 of Algorithm 3. It is supposed toreturn the gradient of the current point. For the original BFGS algorithm,it is obtained by evaluating the true gradient of the objective function atthe current point. For other approximation-based-BFGS methods in ourexperiment, it returns an approximate gradient at the current point by usingsome form of gradient approximation techniques, which we will discuss inthe following sections.102.1. THE BFGS METHODAlgorithm 1: BFGS solverBFGS (x0, f, ε,m)Input : x0: initial pointf : objective functionε: stopping tolerancem: max number of iterationsOutput: x∗: estimated minimizer in the last iterationbegin1 let H0 = I;2 set c1 = 1e− 6, c2 = 0.9;3 set amax = 2;4 set m1 = 20;5 for k = 0 to m do6 set gk = Grad(f, xk);7 if ∥gk∥ < ε then8 break ; // stopping condition triggeredend9 set dk = −Hkgk;10 set αk = WolfeSearch(f, xk, dk,m1, amax, c1, c2);11 set xk+1 = xk + αkdk;12 set s = xk+1 − xk;13 set y = ∇fk+1 −∇fk;14 set Hk+1 =(I − syTsT y)Hk(I − ysTsT y)+ ssTsT y;end15 return x∗ = xk;end112.1. THE BFGS METHODAlgorithm 2: Wolfe line searchWolfeSearch (f, x0, d,m1, amax, c1, c2)Input : f : objective functionx0: current pointd: search directionm1: max iteration numberamax: a given max step sizec1: 1st Wolfe condition parameterc2: 2nd Wolfe condition parameterOutput: α∗: appropriate step sizebegin1 evaluate f0 = f(x0);2 set g0 = Grad(f, x0);3 let α0 = 0, α1 = 0.5αmax;4 for j = 1 to m1 do5 set xj = x0 + αjd;6 evaluate fj = f(xj);7 if fj > f0 + c1αjgT0 d or (fj ≥ fj−1 and j > 1) then8 set αj = Zoom(f, x0, d, αj−1, αj , c1, c2) ; // searchforward9 break;end10 set gj = Grad(f, xj);11 if |gTj d| ≤ −c2gT0 d then12 break;end13 if gTj d ≥ 0 then14 αj = Zoom(f, x0, d, αj , αj−1, c1, c2) ; // searchbackward15 break;endset αj+1 = αj + 0.5(αmax − αj);end16 return α∗ = αj ;end122.1. THE BFGS METHODAlgorithm 3: Wolfe ZoomZoom (f, x0, d, αlo, αhi, c1, c2)Input : f : objective functionx0: current pointd: search directionαlo: lower bound of search stepαhi: upper bound of search stepc1: 1st Wolfe condition parameterc2: 2nd Wolfe condition parameterOutput: α∗: appropriate step sizebegin1 evaluate f0 = f(x0);2 set g0 = Grad(f, x0);3 for j = 1 to m2 do4 set αj =αlo+αhi2 ;5 set xj = x0 + αjd, xlo = x0 + αlod;6 evaluate fj = f(xj), flo = f(xlo);7 if fj > f0 + c1αjgT0 d or fj ≥ flo then8 set αhi = αj ;else9 set gj = Grad(f, xj);10 if |gTj d| ≤ −c2gT0 d then11 break;end12 if gTj d(αhi − αlo) ≥ 0 then13 set αhi = αlo;end14 set αlo = αj ;endend15 return α∗ = αj ;end132.2. GRADIENT APPROXIMATION USING FORWARD DIFFERENCE2.2 Gradient approximation using forwarddifferenceThere are many methods for estimating derivatives. Forward differenceis the most basic one of them. By the Taylor’s theorem, an approximatederivative of single-variable function f(x) can be derived from a first-orderTaylor expansion, as follows,f(x+ h) = f(x) + hf ′(x) +O(h2),hf ′(x) = f(x+ h)− f(x) +O(h2),f ′(x) =f(x+ h)− f(x)h+O(h).If we choose the sampling radius h to be sufficiently small, then we havea small remainder term and get a good approximation of the derivative asthe formula shows above. For a multi-variable function, We can apply thisdivided difference formula to each partial derivative component to get agradient approximation,AG1 =f(x+he1)−f(x)hf(x+he2)−f(x)h...f(x+hen)−f(x)h = ∇f(x) +O(h), (2.14)where ei is the unit vector in ith coordinate unit direction. We then get anapproximated gradient AG1 using forward difference, and its gradient erroris of O(h).2.3 Gradient approximation using centraldifferenceAnother common way to estimate a derivative is central difference. Themethod can be derived from the difference of two second-order Taylor ex-pansions, as follows,f(x+ h) = f(x) + hf ′(x) +12h2f ′′(x) +O(h3),f(x− h) = f(x)− hf ′(x) + 12h2f ′′(x) +O(h3),142.4. GRADIENT APPROXIMATION USING COMPLEX STEP2hf ′(x) = f(x+ h)− f(x− h) +O(h3),f ′(x) =f(x+ h)− f(x− h)2h+O(h2).If we choose the sampling radius h to be sufficiently small, then we havea small remainder and get a good approximation of the derivative. For amulti-variable function, We can apply this formula to each partial derivativecomponent to get a gradient approximation,AG2 =f(x+he1)−f(x−he1)2hf(x+he2)−f(x−he2)2h...f(x+hen)−f(x−hen)2h = ∇f(x) +O(h2). (2.15)Here we get an approximated gradient AG2 using central difference, and itsgradient error is of O(h2).2.4 Gradient approximation using complex stepThe third gradient approximation method that we apply to our researchuses a complex step. The expression can be derived from a second-orderTaylor expansion, as follows,f(x+ ih) = f(x) + ihf ′(x)− 12h2f ′′(x) + iO(h3),Im(f(x+ ih)) = hf ′(x) +O(h3),Im(f(x+ ih))h= f ′(x) +O(h2),f ′(x) =Im(f(x+ ih))h+O(h2).If we choose the sampling radius h to be sufficiently small, then we geta good approximation of the derivative. We can also apply this formulato each partial derivative component of a multi-valued function to get agradient approximation,AG3 =Im(f(x+ihe1))hIm(f(x+ihe2))h...Im(f(x+ihen))h = ∇f(x) +O(h2). (2.16)152.5. APPROXIMATION-BASED-BFGS METHODSHere we get an approximated gradient AG3 using central difference, and itsgradient error is of O(h2).2.5 Approximation-based-BFGS methodsFrom the discussion Section 2.1, we can see that the BFGS is a derivative-based optimization method. If we replace the gradient calculations with gra-dient approximation, particularly for Equations (2.5), (2.7), and (2.13) inthe BFGS algorithm, then we can get an optimization solver without usingthe derivative. We apply three different gradient approximation strategiesthat have been indicated in Equations (2.14), (2.15), and (2.16) to formthree approximation-based-BFGS methods. We then have four optimiza-tion solvers: the BFGS with exact gradients, BFGS with forward-differenceapproximated gradient, BFGS with central-difference approximated gradi-ent, and BFGS with complex-step approximated gradient. We label theseas provided in Table 2.4. We numerically compare these four methods inthe following chapters.Table 2.4: Solver setSovlers name Source Code Gradient errororiginal BFGS BFGS.m O(0)BFGS with forward-difference BFGS AP Oh.m O(h)BFGS with central-difference BFGS AP Oh2.m O(h2)BFGS with complex-step BFGS AP CS.m O(h2)16Chapter 3Experimental set up3.1 Test SetWe build a test set with optimization problems retrieved from [MGH81].A total of 35 unconstrained optimization problems are provided in [MGH81],with a recommended initial point and the true optimal value of each prob-lem. We implement the 35 objective functions with their gradient functionsin Matlab, and verify the correctness of each function and its gradient bycomparing the function values at initial point and the minimizer (if it isprovided). Implementations of these function can be found in the projectrepository by the function name, e.g. file Rosenbrock.m is the Rosenbrockfunction and file D Rosenbrock.m is its gradient.We first apply the original BFGS solvers to each function to find outwhether it can be solved from the recommended starting point. The solversfail or find a wrong minimum on six of the test problems, these are Jennrichand Sampson, Meyer, Gulf research and development, Brown and Dennis,Penalty I, and Broyden Tridiagonal. We remove the Gulf research and de-velopment function as its objective function is nonsmooth. For other fivefunctions above, the BFGS solver can find the true minimum if we choose acloser starting point. However, this change the difficulties of the problems,so we do not use these five functions in our test set.Table 3.1 lists the 29 problems that are used in the final testing. Wenote that the authors only provide six significant digits for the true optimalvalue of each problem in [MGH81]. We may need a higher precision to testthe reliability of the solvers. Therefore we collect the best optimal valuefor those non-zero minimums by the original BFGS, and use those values as“true minimums”. The minimum value provided by [MGH81] and the bestquality minimum value found by BFGS appear in columns 3 and 4 of Table3.1.173.1. TEST SETTable 3.1: Test problemsFunction Name dim. minimum value used minimumfrom [MGH81] valueRosenbrock 2 0 0Freudenstein and Roth 2 48.9842 48.9842536792Powell badly scaled 2 0 0Brown badly scaled 2 0 0Beale 2 0 0Helical valley 3 0 0Bard 3 8.21487e-3 8.21487730658e-3Guassian 3 1.12793e-8 1.12793276962e-8Box three-dimensional 3 0 0Powell singular 4 0 0Wood 4 0 0Kowalik and Osborne 4 3.07505e-4 3.07505603849e-4Osborne 1 5 5.46489e-5 5.46489469748e-5Biggs EXP6 6 5.65565e-3 5.6556499255e-3Osborne 2 11 4.01377e-2 4.01377362935e-2Watson 6 2.28767e-3 2.28767005355e-3Extended Rosenbrock 8 0 0Extended Powell singular 12 0 0Penalty II 10 2.93660e-4 293660537457e-4Variable dimensioned 4 0 0Trigonnometric 5 0 0Brown almost-llinear 4 0 0Discrete boundary value 10 0 0Discrete integral equation 10 0 0Broyden banded 4 0 0Linear - full rank 4 6 6Linear - rank 1 4 15/7 15/7Linear - r1 with 0 cols/rows 4 62/17 62/17Chebyquad 10 6.50395e-3 6.50395480088e-3183.2. EXPERIMENTAL PARAMETERS3.2 Experimental parametersFor each problem in the test set, we apply all four solvers we discussedin Section 2.5 with different parameter settings, then record the optimalresults, running time, and total number of function evaluation for furtheranalysis. The experiment code can be found in the file experiment run.min the project repository.3.2.1 Stopping toleranceAs discussed in the Section 2.1, we will set a small positive number εto be the stopping tolerance to test if the solver finds a local minimum.Considering the dimension n of the problem will also affect the precisionof testing the tolerance, we normalize all test problems to the same levelprecision, we define the stopping condition as∥∇˜f(x¯)∥√n< ε, (3.1)where ∇˜f(x¯) is the approximation gradients of the current point. If thesolver runs over the maximum iteration setting and still does not meet thiscondition, we say the solver fails to solve the test problem.3.2.2 Optimal toleranceIf the algorithm terminates properly, then we calculate the optErrorbetween the returned function value and the true optimal value. We thencheck if the optError is within a predefined optimal tolerance optTol. TheoptError is defined byoptError =∣∣f(x¯)− f(x∗)∣∣, if minimum = 0∣∣∣∣f(x¯)−f(x∗)f(x∗) ∣∣∣∣, otherwise, (3.2)where f(x∗) is the true optimal value of the program, and f(x¯) is the functionvalue of last iteration.If optError < optTol, then we say the solver has solved the optimizationproblem, since the process has found a local minimum and the solution isclose to the true optimal value within a desired precision. We can see theoptimal tolerance is also related to the stopping tolerance we select, for ifthe stopping tolerance is large, the solver will stop too early to find a highprecision minimum function value. Therefore in our experiments, we choose193.3. PERFORMANCE MEASUREMENTSthe optimal tolerance as the same number as the stopping tolerance so thatthe solver can reach the same level of precision of the true optimal.3.2.3 Sampling radiusFor the three versions of the BFGS solver using gradient approximation,the gradient error is of order O(h) and O(h2). Therefore the sampling radiush is an important parameter for these three solvers. Since we will comparethe norm of the gradient to the stopping tolerance, if the gradient erroris larger than the stopping tolerance, then the algorithm is unlikely to con-verge. Technically, selecting a small sampling radius h will increase the preci-sion of the gradient approximation, but it could also encounter some numeri-cal errors, like round-off error, if the sampling radius is too small. To explorethe impact of sampling radius along with different stopping tolerance, weconduct 9 experiments on different combinations of stopping tolerance ε andsampling radius h, where (ε, h) ∈ {10−3, 10−6, 10−9} × {10−3, 10−6, 10−9}.3.3 Performance measurementsTo explore and compare the efficiencies of the four solvers, we collecttwo primary measurements, the number of fundamental evaluations and therunning time of each solver uses on each test problem.For the number of fundamental evaluations, we count all the functionevaluations and gradient evaluations used during the solution process. Eachfunction evaluation will be counted as one fundamental evaluation. Thegradient evaluations will be counted sightly differently.For the BFGS with forward-difference gradient approximation, each gra-dient approximation will be counted as n fundamental evaluations, wheren is the dimension of the problem. This follows from the need to evaluatef(x+ hei) for each i ∈ {1, 2, ...n} to construct one gradient approximation.Note that, although f(x) is also used during the gradient approximation, itis already counted through function evaluations. Thus, there is no need toinclude it again in the fundamental evaluation counter.For the BFGS with central-difference gradient approximation, each gra-dient approximation will be counted as 2n fundamental evaluations. Thisfollows from the need to evaluate f(x + hei) and f(x − hei)for each i ∈{1, 2, ...n} to construct one gradient approximate evaluation.For the BFGS with complex-step gradient approximation, each gradientapproximation will be counted as n fundamental evaluations, where n is the203.4. PERFORMANCE PROFILESdimension of the problem. This follows from the need to evaluate f(x+ihei)for each i ∈ {1, 2, ...n} to construct one gradient approximation.For the original BFGS solver, which calculates the real gradient, weacknowledge that some Mathematicians suggest each gradient evaluationshould be counted as one fundamental evaluation, while others may suggestas n fundamental evaluations. To keep a consistent standard with the otherthree versions, we decide to count it as n fundamental evaluation.For the running time, we collect the elapsed time that each solver workson a test problem. To minimize the uncertainty and influence of the runningplatform, we run the same solver 20 times on the same optimal problem,and then record the average the elapsed time of the fastest 5 runs.3.4 Performance profilesWith these measurements, we then create graphical comparisons of theseoptimization solvers, which are from the performance profiles benchmarkingmethod in [DM02]. That is, for a problem p from a problem set P and asolver s from a solver set S, we first define its performance ratio asrp,s ={tp,smin{tp,s:s∈S} , if convergence test passed,∞, if convergence test failed. (3.3)Then we define the performance profile of a solver s as followsρs(α) =|{p ∈ P, rp,s ≤ α}||P| , (3.4)where |P| denotes the cardinality of the set P. Note that the ρs(α) representsthe portion of problems solved by solver s within a performance ratio α ∈R+.Figure 3.1 shows an example of performance profiles. The x-axis indi-cates the run-time ration comparing the run-time of a solver to the timeof the fastest solver solving the same test problem. The y-axis indicatesthe portion of the test problems being solved by a solver s within a certainrun-time ratio.3.5 Data profilesAnother benchmarking method, data profiles, was proposed in [MW09]as an refinement to performance profiles, particularly for interest in per-formance in terms of the number of fundamental evaluations of a solver.213.5. DATA PROFILES2 4 6 8 10 12 1400.20.40.60.81sBFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2Figure 3.1: An example of performance profilesThe main concern of this method is the required number of fundamentalevaluations to satisfy the convergence test is likely to increase as dimensionof problem (i.e. the number of variables) grows. The data profile of anoptimization algorithm s is defined as,ds(k) =|{p ∈ P, tp,snp+1 ≤ k}||P| , (3.5)where tp,s is the number of fundamental evaluations of a convergence test,or is infinity for a divergence test, np is the dimension of the problem p ∈ P,and ds(k) is the portion of problems in the test set that can be solvedwith k(np + 1) fundamental evaluations. The value k(np + 1) is used as acomparison of the equivalent of k simplex gradient estimates [MW09]. Withthis a scaling factor (np + 1), the data profiles reflect the true evaluationcost of problems in a test set despite of its dimension. Therefore, We applyit to the fundamental evaluation measurement of our experiment results.Figure 3.2 shows an example of data profiles. The x-axis indicates nor-malized function evaluations, which is the total function evaluations of a223.5. DATA PROFILESsolver divided by n + 1, where n is the problem dimension. The y-axis in-dicates the portion of the test problems being solved by a solver s within acertain amount of function evaluations.0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2Figure 3.2: An example of data profilesWe generate graphic plots for each experiments with performance profilesand data profiles to demonstrate and compare performance of four differentsolvers. Results are discussed in the next chapter.23Chapter 4Experimental results anddiscussionsIn our experiments, we define that a solver solves a test problem if thealgorithm converges and the final objective function value is close to thetrue minimum. In particular, using the terminology from Section 3.2 in eachexperiment, optTol is set to ε, so we define a solver solves a test problem ifoptError < ε. With this criteria, we discuss all results of our experimentsin this chapter.4.1 Result on ε = 10−3 and h = 10−3Figure 4.1 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−3 and sam-pling radius h = 10−3. From the graph, we can see that the BFGS, theBFGS-AG-Oh2 and the BFGS-AP-CS solvers have a good reliability using thesesettings, that is, all these three solvers solve over 80% of the problems in thetest set. Conversely, the BFGS-AG-Oh solver only solved about 20% of testproblems using these settings. The BFGS-AG-Oh solver’s poor reliability isprobably due to it has a O(h) gradient error. That is, the sampling radiusis not small enough, so it often failed to meet the stopping criteria.Note that the BFGS solver can not solve every problem in the test set.As discussed in Section 3.1, all problems we chose are solvable by the BFGSsolver when ε is set to 10−9. However, in this experiment we found thatwhen the stopping tolerance is not set to a sufficient small number, thenthe BFGS solver occasionally stopped at a false minimum. In particular,in these cases the norm of the gradient met the stopping criterion in theBFGS method, but the final answer was not sufficiently accurate to pass thedefinition of solves. We will discuss this issue further in the last section ofthis chapter after we go through all the experiment results.For the efficiency, the BFGS solver has the best performance in run-timeand fundamental evaluation counts, since it calculates the true gradient.244.2. RESULT ON ε = 10−3 AND H = 10−62 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.1: Performance profiles and data profiles at ε = 10−3 and h = 10−3The BFGS-AG-CS solver reaches a very good efficiency comparing to the BFGSsolver, as it converged as quickly as the original BFGS algorithm and useda similar amount of fundamental evaluations to obtain the solutions. Itsrun-time performance is worse than the BFGS solver, since it takes longer toobtain function values at n complex valued points than to obtain a functionand gradient value at 1 real valued point. The BFGS-AG-Oh2 solver alsoworks with an efficiency in run-time, almost as good as the BFGS-AG-CSsolver, but from the view of fundamental evaluations, it is not as efficient asthe BFGS-AG-CS. The reason is it needs double the number of fundamentalevaluations in its approximate gradient.4.2 Result on ε = 10−3 and h = 10−6Figure 4.2 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−3 and sam-pling radius h = 10−6. The BFGS solver has the same reliability and per-formance as the previous experiment. This is expected, since the originalBFGS algorithm does not use the sampling radius so the change of it hasno effect on the BFGS solver.The reliability of the BFGS-AG-Oh solver increases significantly in thissettings, where it can solve over 70% of the test problems. Applying asmaller sampling radius h to the BFGS-AG-Oh solver increases the precisionof the gradient approximation, and its gradient error is now much smaller254.3. RESULT ON ε = 10−3 AND H = 10−92 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.2: Performance profiles and data profiles at ε = 10−3 and h = 10−6than the stopping tolerance in the case, which improves reliability of theBFGS-AG-Oh solver. The BFGS-AG-CS solver and the BFGS-AG-Oh2 solverkeep the same level of reliability as the BFGS solver.For the running time, the BFGS-AG-Oh solver performs about 20% fasterthan the other two approximation-based solvers on most problems that itsolved. The BFGS-AG-CS solver and the BFGS-AG-Oh2 solver have simi-lar performance on running time. For the fundamental evaluations, theBFGS-AG-CS solver performs almost as good as the BFGS solver, and it workssightly better than the BFGS-AG-Oh solver. The BFGS-AG-Oh2 solver takesabout 2 times of the number of fundamental evaluations when compared tothe other three solvers.4.3 Result on ε = 10−3 and h = 10−9Figure 4.3 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−3 and sam-pling radius h = 10−9. The the BFGS-AG-CS solver and the BFGS-AG-Oh2solver have the same reliability as they did in the previous experiment. TheBFGS-AG-Oh solver improves its reliability to solve about 80% of the prob-lems in the test set.For the running time performance, all solvers keep a similar trend asthey did in the previous experiment. The BFGS-AG-Oh solver runs a slightlybit faster than the other two solvers with gradient approximation. For the264.4. RESULT ON ε = 10−6 AND H = 10−32 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.3: Performance profiles and data profiles at ε = 10−3 and h = 10−9performance on the fundamental evaluations, all solvers keep almost thesame as the previous experiment.4.4 Result on ε = 10−6 and h = 10−3Figure 4.4 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−6 and sam-pling radius h = 10−3. The BFGS solver works more reliable in this setting,and it can solve over 95% of the problems in the test set. This is due tosetting the stopping tolerance to a smaller value, so the BFGS solver conductsmore iterations to search a more accurate minimum.The reliability of the BFGS-AG-CS solver and the BFGS-AG-Oh2 solverdrop a bit comparing to the previous experiments. The BFGS-AG-Oh2 solvercan solve about 70% of the test problems, and the BFGS-AG-CS solver cansolver about 60% of them. Both of these two solvers have a O(h2) gradienterror which is at the same level precision as the stopping tolerance setting.From the result we can see this setting will affect their reliability. TheBFGS-AG-Oh2 solver failed to converge on 6 test problems out of total 29problems in the test set, and the The BFGS-AG-CS solver failed to convergeon 9 test problems in this setting.The reliablity of the BFGS-AG-Oh solver drops tremendously, and it cansolve almost none of the test problems. The BFGS-AG-Oh solver has gradienterror of O(h), which is a worse level of precision than the stopping tolerance.274.5. RESULT ON ε = 10−6 AND H = 10−62 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.4: Performance profiles and data profiles at ε = 10−6 and h = 10−3As expected, this has a negative effect on reliability.For the running time performance, the BFGS solver remains the fastest.The BFGS-AG-CS solver and BFGS-AG-Oh2 solver take about 2 times to 10times of the running time on various test problems comparing to the fastestsolver. The BFGS-AG-Oh2 solver works slightly faster then the BFGS-AG-CSsolver in this setting. For the performance of fundamental evaluations, theBFGS solver has the best efficiency in this setting, the BFGS-AG-CS solvertakes a bit more evaluations in solving, and the BFGS-AG-Oh2 solver takesabout double evaluations than the BFGS-AG-CS solver in most those problemsit can solve.Note that the BFGS-AG-CS solver takes about half fundamental eval-uations of the BFGS-AG-Oh2 solver does, but still runs slower than theBFGS-AG-Oh2 solver. This is due to the computation time of function evalu-ation at a complex value is longer than the evaluation time at a real value,as we discussed after Definition 2.1.4.5 Result on ε = 10−6 and h = 10−6Figure 4.5 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−6 andsampling radius h = 10−6. Reliability of the BFGS-AG-CS solver and theBFGS-AG-Oh2 solver increase to over 90% of solved problems in the test set.The BFGS-AG-Oh solver can solve about 40% of the problems in the test set,284.6. RESULT ON ε = 10−6 AND H = 10−92 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.5: Performance profiles and data profiles at ε = 10−6 and h = 10−6which is still in a low reliability level. Among the 18 test problems it fails tosolved, the BFGS-AG-Oh solver fails to converge on 16 of them. This resultshows that if the gradient error is not significantly smaller than the stop-ping tolerance, it has a negative effect on the convergence of the BFGS-AG-Ohsolver.For the running time performance, the BFGS-AG-Oh2 solver works almostas fast as the BFGS-AG-CS solver, but they are still about 2 times to 10 timesof the running time comparing to the fastest solver BFGS. All solvers remaina similar trend on their performance of fundamental evaluations.4.6 Result on ε = 10−6 and h = 10−9Figure 4.6 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−6 and sam-pling radius h = 10−9. Comparing to the last experiment, the samplingradius h is set to 10−9 which is smaller than the stopping tolerance. Thishelps the BFGS-AG-Oh solver to elevate its reliability to a good level, andit can solve about 80% of the problems in the test set. Other solvers havesimilar reliability. Note that the BFGS-AG-CS solver solves more test prob-lems than it does in the last experiment, but the BFGS-AG-Oh2 solver actuallysolves one less test problems in this setting. It is probably due to the impactof catastrophic cancellation in the finite difference used in the BFGS-AG-Oh2solver. In [MSA03], the authors show the impact of catastrophic cancella-294.7. RESULT ON ε = 10−9 AND H = 10−32 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ationFigure 4.6: Performance profiles and data profiles at ε = 10−6 and h = 10−9tion errors on a case study. As the sampling radius (which is referred as“step size” in [MSA03]) decreases below a value of 10−5, the relative er-ror of the central difference derivative approximation increases. We will seehow this affects on the BFGS-AG-Oh2 solver more significantly in the laterexperiments.For the running time performance, the BFGS-AG-Oh solver works slightlyfaster than the BFGS-AG-CS solver and the BFGS-AG-Oh2 solver on most ofthe problems it can solve. For the performance of fundamental evaluations,the BFGS solver and the BFGS-AG-CS solver are the best out of the foursolvers.4.7 Result on ε = 10−9 and h = 10−3Figure 4.7 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−9 and sam-pling radius h = 10−3. The BFGS solver solved all the test problems. Thisis expected as our problem selection was based on BFGS successfully solvingthe problem with this stopping tolerance.Conversely, the other three solvers with gradient approximation havepoor reliability in this setting. The BFGS-AG-Oh solver can not solve any ofthe problems in the test test, and the other two solvers can only solve about30% of the problems in the test set. As the sampling radius h is set to 10−3,the gradient error of the BFGS-AG-Oh solver is much larger than the stopping304.8. RESULT ON ε = 10−9 AND H = 10−62 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time(partial graph)0 200 400 600 800 1000Number of normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (partial graph)Figure 4.7: Performance profiles and data profiles at ε = 10−9 and h = 10−3tolerance. For the other two solvers, theirO(h2) gradient error are also largerthan the stopping tolerance. This results in all these approximation-basedalgorithms failing to meet the convergent conditions in many test problems.Given the poor performance of the algorithms, run times are not particu-larly relevant. However, for the running time performance, the BFGS-AG-Oh2solver works slightly faster than the BFGS-AG-CS solver. For the performanceof fundamental evaluations, the efficiency of the BFGS-AG-Oh2 solver are closeto the BFGS-AG-CS solver.Note that the graphs in Figure 4.7 do not show the whole result in thisexperiment. We cut the far end of the profiles in this figure, in order toshow more detail trends of the majority test problems. Readers may findthe full result profiles figures of this experiment in Figure 4.8.4.8 Result on ε = 10−9 and h = 10−6Figure 4.9 shows the performance profiles of run-time and the data pro-files of fundamental evaluations using stopping tolerance ε = 10−9 and sam-pling radius h = 10−6. As the sampling radius decreases and the gradienterror of the BFGS-AG-CS solver and the BFGS-AG-Oh2 solver get smaller thanthe stopping tolerance, the reliability these two solvers increase to a verygood level. Both solvers can solve about 90% of the problems in the testset. The BFGS-AG-Oh solver still has gradient error worse than the stoppingtolerance, and can only solve about 10% of the test problems.314.9. RESULT ON ε = 10−9 AND H = 10−91000 2000 3000 4000 5000 600000.20.40.60.81s(a) Performance profiles of run-time(full graph)0 0.5 1 1.5 2Normalized function evaluations k 10500.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (full graph)Figure 4.8: Full graph of performance profiles and data profiles at ε = 10−9and h = 10−3For the running time performance, the BFGS-AG-Oh2 solver works as fastas the BFGS-AG-CS solver on about half of the test problems, but it getless efficient than the BFGS-AG-CS solver on the other problems. All solversremain a similar trend on their performance of fundamental evaluations.Note that in Figure 4.9 we cut the far end of the data profiles to showmore details in the near end. Readers may find the full graph of profiles ofthis experiment in Figure 4.10.4.9 Result on ε = 10−9 and h = 10−9Figure 4.11 shows the performance profiles of run-time and the dataprofiles of fundamental evaluations using stopping tolerance ε = 10−9 andsampling radius h = 10−9. Comparing to the last experiment, the samplingradius h is set to 10−9. The BFGS-AG-Oh solver can solve about 40% of theproblems in the test set, but is still not in a good level of reliability. Thecatastrophic cancellation impact we discussed in Section 4.6 will also affectthe forward difference approximation. The example in [MSA03] shows thereis a turning point at sampling radius h = 10−8 where the relative error ofthe derivative estimate becomes increasing as the sampling radius decreases.By this impact, the gradient error of the BFGS-AG-Oh solver is expected tobe larger than the stopping tolerance in this setting, and hence lower thereliability of it.324.9. RESULT ON ε = 10−9 AND H = 10−92 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 200 400 600 800 1000Normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (partial graph)Figure 4.9: Performance profiles and data profiles at ε = 10−9 and h = 10−62 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time0 2000 4000 6000 8000 10000Normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (full graph)Figure 4.10: Full graph of performance profiles and data profiles at ε = 10−9and h = 10−6334.9. RESULT ON ε = 10−9 AND H = 10−92 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles of run-time(partial graph)0 200 400 600 800 1000Normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (partial graph)Figure 4.11: Performance profiles and data profiles at ε = 10−9 and h = 10−9The BFGS-AG-CS solver elevates its reliability and solves over 95% of testproblems. Conversely, the BFGS-AG-Oh2 solver actually drops its reliabilitycomparing to the last experiment, and it can now solve about 70% prob-lems. The catastrophic cancellation affects the BFGS-AG-Oh2 solver moresignificantly in this setting than it does in Section 4.6. As the stoppingtolerance is set to a smaller value, the solver is expected to find a more ac-curate minimum, but if the gradient error is not sufficiently small comparedto the stopping tolerance, then the algorithm will fail to converge in moretest problems. Therefore, although the gradient error of the BFGS-AG-Oh2solver remains at the same level in this experiment and the one discussed inSection 4.6, it actually solves less problems in the test set.For the running time performance, the BFGS-AG-Oh solver works slightlyfaster than the other two approximation-based solvers on simple problemsin the test set. The BFGS-AG-CS solver exceeds the BFGS-AG-Oh2 solver onrunning time efficiency. For the performance of fundamental evaluations,the BFGS-AG-CS solver is almost as good as the BFGS solver.Note that in Figure 4.11 we cut the far end of the performance profilesand data profiles to show more details in the near end. Readers may findthe full graph of profiles of this experiment in the Figure 4.12.344.10. COMPREHENSIVE RESULT ANALYSIS5 10 15 20 2500.20.40.60.81s(a) Performance profiles of run-time(full graph)0 1000 2000 3000 4000 5000Normalized function evaluations k00.20.40.60.81d s(k)BFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Data profiles of fundamental evalu-ation (full graph)Figure 4.12: Full graph of performance profiles and data profiles at ε = 10−9and h = 10−94.10 Comprehensive result analysisBy the results of these nine experiments, we analyze the overall trendsin reliability and efficiency of these four solvers.4.10.1 The BFGS solverWe choose the BFGS solver as a comparing base to conduct these ex-periments, and selected test problems that were all solvable by the BFGSsolver, so it is not a surprise that the BFGS solver is quite reliable. However,from these experiments we can see the setting of the stopping tolerance isstill an important factor for this algorithm. The BFGS solver can not solveevery problem in the test set when the stopping tolerance is not set to asufficient small number. By checking all the test problems it failed to solvein our experiments, we found the reason is that the algorithm stopped ata false minimum where the norm of its gradient met the stopping criteria,which is the only condition to test the convergence of the algorithm. Ta-ble 4.1 shows solving percentage of the BFGS solver over different stoppingtolerance settings.From Table 4.1 we can see that the BFGS is a very efficient solver. Evenwhen the stopping tolerance is set to a small value, it takes only a fewmore fundamental evaluations for the algorithm to reach the more accurateminimum.354.10. COMPREHENSIVE RESULT ANALYSISFrom the results of performance profiles, we can also see that the BFGSsolver is the most efficient solver, as it is always the fastest for the runningtime on the test problems it solved. Other solvers take around 2 times to10 times with respect to the running time of the BFGS solver on most of theproblems in the test set. Even when a solver, e.g. the BFGS-AG-CS solver,has a very similar number of fundamental evaluation on some problem, it isstill slower than the BFGS solver in running time, as the BFGS solver takesadvantage of evaluating the true gradient by a single function call. This isexpected, as it is almost always advantageous to use true gradient valuesinstead of approximations when available.Table 4.1: Stopping tolerance affect on the BFGS solverStopping tolerance Solved test Median of fundamentalproblem percentage evaluations10−3 82.76% 102.8510−6 96.55% 112.7510−9 100% 1214.10.2 The BFGS-AG-Oh solverFor the BFGS-AG-Oh solver, the sampling radius h is also an importantparameter to the solver. As the BFGS-AG-Oh solver has a gradient error ofO(h), it performs poorly when the sampling radius is set to be larger thanthe stopping tolerance (see Figures 4.4, 4.7 and 4.9). When the samplingradius is set to be equal to the stopping tolerance, the reliability of theBFGS-AG-Oh solver is around 40% and it is still in a low level (see Figures4.1, 4.5 and 4.11). Although in these three experiments, the gradient errorof the BFGS-AG-Oh solver is at the same level as the stopping tolerance, it isstill too large to guarantee the algorithm converges for some test problems,since the gradient error will become a significant part when the true gradientis close to zero. When the sampling radius is set to be less than the stoppingtolerance, the reliability of the BFGS-AG-Oh solver improves to an acceptablelevel, where the solver can solve about 70% to 80% of the problems in thetest set (see Figures 4.2, 4.3 and 4.6). However, even in these experimentsettings, the reliability of the BFGS-AG-Oh solver is still worse than the othertwo approximation-based solvers. This is likely due to when being in thesame setting, the BFGS-AG-CS solver and the BFGS-AG-Oh2 solver alwayshave a smaller gradient error comparing to the BFGS-AG-Oh solver.364.10. COMPREHENSIVE RESULT ANALYSISWhen focusing on those settings that it works well in, we can see theBFGS-AG-Oh solver is actually pretty good at efficiency compared to theother two approximation-based solvers. For the fundamental evaluations, itusually catches up with the top solver on those problems it can solve. Forrunning time performance, it runs slightly faster than the other two on mostof the test problems it can solve. Comparing to the BFGS-AG-Oh2 solver, ithas less fundamental evaluations and both of them evaluate the objective atreal values. Comparing to the BFGS-AG-CS solver, the BFGS-AG-Oh solver hasslightly more fundamental evaluations but it evaluates the objective functionat real values which takes less computing time than the BFGS-AG-CS solverwhich evaluates the objective function at complex values.4.10.3 The BFGS-AG-Oh2 solverThe BFGS-AG-Oh2 solver works well in the majority of the experimentsettings. Often the reliability of it is quite close to the BFGS solver. Insome particular cases it even solved one more test problem than the BFGSsolver (see Figures 4.2, 4.2 and 4.3). However, the BFGS-AG-Oh2 solverworks poorly in experiments where it has a gradient error larger than, orequal to, the stopping tolerance (see Figures 4.4 and 4.7). The reliabilityof the BFGS-AG-Oh2 solver in these two experiments is far behind the BFGSsolver. For the running time performance, all three solvers are actually ata very similar efficiency level. In some settings, the BFGS-AG-Oh2 solver isfaster than the BFGS-AG-CS solver (see Figure 4.4), whereas in most of theother cases their running time performances are close to each other. Onedisadvantage of the BFGS-AG-Oh2 solver is that it takes more fundamentalevaluations than the other solvers in most of the test problems. This is dueto the central-difference gradient approximation it uses, where each gradientapproximation takes 2n function calls to evaluate. This is reflected in mostof the results, which we can see in the data profiles. The BFGS-AG-Oh2 solverusually takes about twice the amount of fundamental evaluations comparingto the top solver, which is normally the BFGS solver or the BFGS-AG-CS solver.One issue of the BFGS-AG-Oh2 solver is that it actually got worse whenthe sampling radius h was set to the tiny number 10−9. We see this trendby comparing the reliability of Figure 4.5 to Figure 4.6, and also the reli-ability of Figure 4.9 to Figure 4.11. In both comparisions, the portion ofsolved test problems drops when h changes from 10−6 to 10−9. This indi-cates the BFGS-AG-Oh2 solver has some difficulties in handle tiny samplingradius, which is likely caused by the numerical error of the catastrophic can-cellation. In [MSA03], Martins, Sturdza, and Alonso provides an example374.10. COMPREHENSIVE RESULT ANALYSISshowing how this catastrophic cancellation affects three different derivativeapproximation methods. In that example, the relative error of the centraldifference approximation increases as the sampling radius decreases below10−5.To verify if this impact is also applied to the BFGS-AG-Oh2 solver, we dotwo extra experiments on setting sampling radius h to 10−12 and 10−15, withthe stopping tolerance ε = 10−9. We show the performance profiles of thesetwo experiments in Figure 4.13, and put the reliability of the BFGS-AG-Oh2solver and the BFGS-AG-CS solver with the results from two other relatedexperiments (see Figures 4.9 and 4.11) in Table 4.2 for comparison. Theresults show the catastrophic cancellation does affect the BFGS-AG-Oh2 solverand amplifies its gradient error similarly to the examples in [MSA03]. Thecatastrophic cancellation does not affect the BFGS-AG-CS solver since thereis no subtraction needed in the complex-step gradient approximation.Table 4.2: Impact on accuracy of different sampling radius hExperimental settings Solved percentage of Solved percentage ofthe BFGS-AG-Oh2 the BFGS-AG-CSε = 10−9 and h = 10−6 89.66% 93.76%ε = 10−9 and h = 10−9 72.41% 96.55%ε = 10−9 and h = 10−12 58.62% 93.76%ε = 10−9 and h = 10−15 31.03% 100%4.10.4 The BFGS-AG-CS solverFor the BFGS-AG-CS solver, it also uses gradient approximation withO(h2) gradient error. Therefore, it has bad performance when its gradienterror is larger than, or equal to, the stopping tolerance (see Figures 4.4 and4.7). This is similar to the BFGS-AG-Oh2 solver. Except these two cases,the BFGS-AG-CS solver works well in all other experiments. The reliabilityof the BFGS-AG-CS solver is very close to the BFGS solver, and it can solvealmost all problems that solved by the BFGS, sometimes only missing one ortwo problems, in the same experiment setting. Also, from the extra resultsshown in Table 4.2, the BFGS-AG-CS solver is very stable when handling tinysampling radius. This is due to the complex-step derivative approximationhaving the advantage of being step-size insensitive and its relative estimateerror remains in a good level for small enough sampling radius (see [MSA03]).For the fundamental evaluations, the BFGS-AG-CS solver is one of the384.10. COMPREHENSIVE RESULT ANALYSIS2 4 6 8 10 12 1400.20.40.60.81s(a) Performance profiles at ε = 10−9and h = 10−122 4 6 8 10 12 1400.20.40.60.81sBFGSBFGS-AP-CSBFGS-AG-OhBFGS-AG-Oh2(b) Performance profiles at ε = 10−9and h = 10−15Figure 4.13: Performance profiles on setting h to 10−12 and 10−15top performing solvers, tying the BFGS solver on most of the problems inthe test set. For running time performance, the BFGS-AG-CS solver runsslower than the BFGS solver, but it is on the same level of the other twoapproximation-based solvers. Although it takes longer for a single functionevaluation at complex values, the BFGS-AG-CS solver actually catches upwith the other two in running time by taking the advantage of its goodefficiency on fundamental evaluations.39Chapter 5ConclusionIn this research, we have explored an alternative of DFO methods forproblems with a smooth objective function by applying gradient approxima-tion to the BFGS methods. We conducted an experiment on approximation-based-BFGS methods with three different techniques of gradient approxi-mation, comparing to the original BFGS method. For the experiment, weselected 29 problems with a smooth objective function from [MGH81] asthe test set, and applied different settings on the stopping tolerance andthe sampling radius on different techniques. Our analysis is based on theexperimental results of performance profiles and data profiles on differentsettings. We conclude by answering the three questions which we broughtforward in Chapter 1.1. Is applying numerical differentiation to the BFGS method a feasiblealternative approach to DFOmethods in solving problems with smoothobjective function?The overall result suggests that the approximation-based-BFGS methodis a feasible alternative to DFO methods. With appropriate parameter set-tings, approximation-based-BFGS methods can achieve the same level ac-curacy as the original BFGS method does, with adding a small amount offunction evaluations and machine run-time comparing to the BFGS method.However, we also notice that the setting of sampling radius affect the ac-curacy of the approximation-based-BFGS methods significantly, due to theexistence of numerical differentiation errors.We applied forward difference, central difference, and complex step gra-dient approximation to the BFGS methods. In general, the approximation-based-BFGS with forward difference gets lower accuracy than the other twoapproximation-based-BFGS solvers, due to it has relatively large numericaldifferentiation error in the same setting. The approximation-based-BFGSwith central difference achieved outstanding accuracy at setting sampling ra-dius h = 10−6, but gets worse gradually as the sampling radius gets smaller,affected by the catastrophic cancellation. The approximation-based-BFGS40with complex step gradient approximation obtains excellent and stable ac-curacy it most of the settings.2. What is the difference among those approximation-based-BFGS meth-ods in terms of function evaluations and run-time performance?By the result of data profiles, the approximation-based-BFGS with com-plex step only takes a few more function evaluations than the originalBFGS method. The approximation-based-BFGS with central differenceneeds about double the function evaluations of the original BFGS, due tothe way it calculates approximated gradient. For the run-time performance,the approximation-based-BFGS methods takes about 2 to 8 times of theCPU time of the original BFGS. The approximation-based-BFGS with for-ward difference runs slightly faster than the other two approximation-based-BFGS solvers in the problems it can solve. The approximation-based-BFGSwith complex step and the one with central difference are at about the samelevel of run-time performance. The complex-step-approximated-BFGS takesmore computational time in each function evaluation at complex variables,while the central-difference-approximated-BFGS takes more function evalu-ations at each iteration.3. There are some parameter settings used within numerical differentia-tion, such as sampling radius and stopping tolerance. How does thechoices of these parameters affect the performance of different solvers?General speaking, the settings of sampling radius and stopping toleranceare important to the approximation-based-BFGS methods. When choosingthese parameters, we should try to guarantee the numerical differentiationerrors of the methods are much smaller than the stopping tolerance. On theother hand, the stopping tolerance setting also controls the search processof the solvers and limits the precision of the minimizer it found. When werequire a higher precision solution of the optimization problems, we need asmaller stopping tolerance, and a smaller sampling radius correspondingly.The approximation-based-BFGS with complex step gradient approximationmay be a good choice in this case, as it performs well and is still stable whenthe sampling radius is less than 10−9.In summary, the approximation-based-BFGS methods will inherit thespeed and efficiency of the original BFGS method, and will reach as goodof an accuracy as the BFGS can do. We consider the approximation-based-BFGS method a good option in solving optimization problems with a smoothobjective function.415.1. FURTHER RESEARCH DIRECTIONS5.1 Further research directionsWe have achieved some good results in the study of the feasibility of theapproximation-based-BFGS methods, whereas there are still improvementsthat can be made in our research, and further study of applying approxima-tion techniques in other DFO methods can also be explored. We summarizedsome possible directions of further research as follows.First, as we have seen the parameter settings have an obvious impacton the approximation-based-BFGS method, there may be some refinementthat can be made to the parameter settings, especially for the samplingradius. One idea is to decrease the sampling radius gradually during thesearch process of the algorithm. In such way that there will be no need forpre-set sampling radius in the solver, and it might improve the efficiency ofthe search process.Another possible research direction is to study how the approximation-based-BFGS method works in real world optimization problems. The ul-timate goal is to test the effectiveness of the approximation-based-BFGSmethod as an alternative of DFO methods working on problems withoutderivative information. So it is indispensable to apply this method to realcases.Lastly, as shown in the experimental result, we found that the complex-step approach has some advantage in handling high precision gradient ap-proximation, due to its immunity from catastrophic cancellation. In futureresearch, we may apply it to other DFO methods which also require numer-ical differentiation in their process, such as the Implicit Filtering [Kel11], toseek improvement of those methods.42Bibliography[ABSM14] M. Al-Baali, E. Spedicato, and F. Maggioni. Broyden’s quasi-newton methods for a nonlinear system of equations and uncon-strained optimization: a review and open problems. Optimiza-tion Methods and Software, 29(5):937–954, 2014. → pages 8[AD02] C. Audet and J. E. Dennis. Analysis of generalized patternsearches. SIAM Journal on Optimization, 13(3):889–903, 2002.→ pages 2[AD06] C. Audet and J. E. Dennis. Mesh adaptive direct search al-gorithms for constrained optimization. SIAM Journal on Opti-mization, 17(1):188–217, 2006. → pages 2[AF03] M. J. Ablowitz and A. S. Fokas. Complex Variables: Introductionand Applications. Cambridge Texts in Applied Mathematics.Cambridge University Press, 2003. → pages 5[AH17] C. Audet and W. Hare. Derivative-Free and Blackbox Optimiza-tion. Springer, 2017. → pages 1[BN12] J. Bak and D. J. Newman. Complex Analysis. UndergraduateTexts in Mathematics. Springer, 2012. → pages 7[CSV08] A. R. Conn, K. Scheinberg, and L. N. Vicente. Geometry of sam-ple sets in derivative free optimization: Polynomial regressionand underdetermined interpolation. IMA journal of numericalanalysis, 28(4):721–748, 2008. → pages 2[CSV09] A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction toDerivative-Free Optimization. SIAM, 2009. → pages 2[CT96] A. R. Conn and P. L. Toint. An algorithm using quadratic inter-polation for unconstrained derivative free optimization. Springer,1996. → pages 243Bibliography[dB81] N. G. de Bruijn. Asymptotic Methods in Analysis, volume 4.Courier Corporation, 1981. → pages 7[DM02] E. D. Dolan and J. J. More´. Benchmarking optimization soft-ware with performance profiles. Mathematical Programming,91(1):201–213, 2002. → pages 21[Fer52] E. Fermi. Numerical solution of a minimum problem. Technicalreport, Los Alamos Scientific Lab., Los Alamos, NM, 1952. →pages 2[GNS09] I. Griva, S. G. Nash, and A. Sofer. Linear and Nonlinear Opti-mization, volume 108. Siam, 2009. → pages 4[Har10] W. Hare. Using derivative free optimization for constrained pa-rameter selection in a home and community care forecastingmodel. In International perspectives on operations research andhealth care, Proceedings of the 34th meeting of the EURO workinggroup on operational research applied to health sciences, pages61–73, 2010. → pages 1[HLS+22] E. Henrichs, V. Lesch, M. Straesser, S. Kounev, andC. Krupitzer. A literature review on optimization techniquesfor adaptation planning in adaptive systems: State of the artand research directions. Information and Software Technology,149:106940, 2022. → pages 1[HNT13] W. Hare, J. Nutini, and S. Tesfamariam. A survey of non-gradient optimization methods in structural engineering. Ad-vances in Engineering Software, 59:19–28, 2013. → pages 1[HS22] W. L. Hare and K. Srivastava. A numerical study of applyingcomplex-step gradient and hessian approximations in derivative-free optimization. Pacific Journal of Optimization, 2022. →pages 4[Kel99] C. T. Kelley. Iterative Methods for Optimization. Frontiers inApplied Mathematics. SIAM, Philadelphia, 1999. → pages 8[Kel11] C. T. Kelley. Implicit filtering. SIAM, 2011. → pages 42[Ken13] A. R. Kenneth. Elementary Analysis: The Theory of Calculus.Spinger, 2013. → pages 744Bibliography[KLT03] T. G. Kolda, R. M. Lewis, and V. Torczon. Optimization bydirect search: New perspectives on some classical and modernmethods. SIAM review, 45(3):385–482, 2003. → pages 2[LM67] J. N. Lyness and C. B. Moler. Numerical differentiation of ana-lytic functions. SIAM Journal on Numerical Analysis, 4(2):202–210, 1967. → pages 4[LMW19] J. Larson, M. Menickelly, and S. M. Wild. Derivative-free opti-mization methods. Acta Numerica, 28:287–404, 2019. → pages3[LTY12] M. Lin, J. Tsai, and C. Yu. A review of deterministic optimiza-tion methods in engineering and management. MathematicalProblems in Engineering, 2012, 2012. → pages 1[MA04] R. T. Marler and J. S. Arora. Survey of multi-objective optimiza-tion methods for engineering. Structural and multidisciplinaryoptimization, 26:369–395, 2004. → pages 1[MGH81] J. J. More´, B. S. Garbow, and K. E. Hillstrom. Testing uncon-strained optimization software. ACM Transactions on Mathe-matical Software, 7(1):17–41, 1981. → pages 17, 18, 40[MSA03] J. R. R. A. Martins, P. Sturdza, and J. J. Alonso. The complex-step derivative approximation. ACM Transactions on Mathe-matical Software, 29(3):245–262, 2003. → pages 29, 30, 32, 37,38[MW09] J. J. More´ and S. M. Wild. Benchmarking derivative-free opti-mization algorithms. SIAM Journal on Optimization, 20(1):172–191, 2009. → pages 21, 22[NM65] J. A. Nelder and R. Mead. A simplex method for function min-imization. The Computer Journal, 7(4):308–313, 1965. → pages2[NW99] J. Nocedal and S. J. Wright. Numerical Optimization. SpringerSeries in Operations Research and Financial Engineering.Springer, 1999. → pages 8, 9, 10[Reg15] R. G. Regis. The calculus of simplex gradients. OptimizationLetters, 9:845–865, 2015. → pages 245Bibliography[Ros60] H. Rosenbrock. An automatic method for finding the greatest orleast value of a function. The computer journal, 3(3):175–184,1960. → pages 6[SCZZ19] S. Sun, Z. Cao, H. Zhu, and J. Zhao. A survey of optimizationmethods from a machine learning perspective. IEEE transactionson cybernetics, 50(8):3668–3681, 2019. → pages 1[ST98] W. Squire and G. Trapp. Using complex variables to estimatederivatives of real functions. SIAM review, 40(1):110–112, 1998.→ pages 4[Vas12] P. M. Vasant. Meta-heuristics optimization algorithms in engi-neering, business, economics, and finance. IGI Global, 2012. →pages 1[VKPS17] A. Verde´rio, E. W. Karas, L. G. Pedroso, and K. Scheinberg.On the construction of quadratic models for derivative-free trust-region algorithms. EURO Journal on Computational Optimiza-tion, 5(4):501–527, 2017. → pages 2[WRS08] S. M. Wild, R. G. Regis, and C. A. Shoemaker. Orbit: Opti-mization by radial basis function interpolation in trust-regions.SIAM Journal on Scientific Computing, 30(6):3197–3219, 2008.→ pages 2[WSS15] S. M. Wild, J. Sarich, and N. Schunck. Derivative-free op-timization for parameter estimation in computational nuclearphysics. Journal of Physics G: Nuclear and Particle Physics,42(3):034031, 2015. → pages 146
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Hare, Warren
Using complex step numerical differentiation in the Broyden-Fletcher-Goldfard-Shanno algorithm
Text