A Limited-Memory BFGS Algorithm Based on a Trust-Region Quadratic Model for Large-Scale Nonlinear Equations.

In this paper, a trust-region algorithm is proposed for large-scale nonlinear equations, where the limited-memory BFGS (L-M-BFGS) update matrix is used in the trust-region subproblem to improve the effectiveness of the algorithm for large-scale problems. The global convergence of the presented method is established under suitable conditions. The numerical results of the test problems show that the method is competitive with the norm method.

Introduction

Consider where b:ℜ → ℜ is a system of nonlinear equations and n denotes the large-scale dimension. Problems involving non-linear equations can be found in numerous applications in engineering, such as nonlinear fitting, function approximating, parameter estimating and leaning models [1-5]. Large-scale nonlinear equations are believed to be difficult to solve because the relations with x are complex and because the dimension is high. Currently, many effective algorithms exist for addressing (1), such as the trust-region method [6-11], Levenberg-Marquardt method [12-14], quasi-Newton method [15-19], and Gauss-Newton method [20-25]. Let θ be the norm function defined by , and suppose that b(x) has a zero point. Then, (1) is equivalent to the following global optimization problem:

In this paper, we will focus on the solution of (2) via the trust-region (TR) methods. The traditional TR methods, at each iterative point x , obtain the trial step d using the following TR subproblem model: where △ > 0 is the so-called TR radius and ‖⋅‖ denotes the Euclidean norm of vectors or their induced matrix. The TR method is very useful when the exact Jacobian or Hessian computation is inexpensive or possible. However, this case is rare in practice. It is not difficult to see that the Jacobian matrix ∇b(x) at every iteration must be computed, which obviously increases the workload and CPU time. Moreover, the global and superlinear convergence of the TR methods often requires that the nondegenerate assumption about ∇b(x*) holds, where x* is a solution of (1) and nondegeneracy means nonsingularity. The nondegeneracy of ∇b(x*) seems to be too stringent of a requirement for the purpose of ensuring the global and superlinear convergence. To overcome the above drawbacks, Yuan et al. [9] presented a new TR subproblem model defined by where the TR radius △ = c ‖b(x )‖, c ∈ (0, 1), p is a nonnegative integer and B is generated by the following BFGS formula: where s = x − x , y = b(x ) − b(x ), x is the next iteration and B 0 is an initial symmetric, positive definite matrix. They established the global and superlinear convergence without nondegeneracy. The numerical results show that the given algorithm is competitive with the normal TR method up to a dimension of 600. The subproblem (4) can be also rewritten as follows Yuan et al. [8] consider the case with the symmetric Jacobian matrix ∇b(x ) and propose the following TR subproblem model: where is defined by the special BFGS update: where δ = b(x + y ) − b(x ). About the above TR methods, the TR subproblems models will be repeatedly computed if the calculated trial step d does not satisfy the given conditions, which obviously increase the workload. To avoid this drawback, Yuan et al. [10] give the following TR model: where In this model, the next point is x = x + d if the trial step d satisfies the successful iteration; otherwise, the next point is x = x + α d , where α > 0 is the so-called steplength designed by where δ ∈ (0, 1) and the technique was proposed by Yuan and Lu [16]. However, the above model can not ensure that the number of searching α > 0 is small. Many TR model methods (see [7, 12] etc.) have been developed for solving nonlinear equations. However, these TR methods require the matrix information (the Jacobian matrix or the BFGS update matrix), which will increase the computational complexity. This observation motivate us find another approach to generate the matrix information requiring minimal memory.

Based on the problem (4), we will use the limited-memory BFGS (L-M-BFGS) method instead of the BFGS update because the former often provides a fast rate of linear convergence and requires minimal memory. The L-BFGS method is an adaptation of the standard BFGS method [26]. The implementation is identical to that of the normal BFGS method, the difference between the L-M-BFGS method and the BFGS method is that the inverse Hessian approximation is not explicitly formed, but defined by a small number of BFGS updates. At every iteration x , the L-M-BFGS method stores a small number (say, m) of correction pairs {s , y }(i = k − 1, …, k − m) to obtain H , instead of storing the matrices H , where with H being the inverse of B . The L-M-BFGS update formula is defined as where and . It has been proved that the L-M-BFGS method is very suitable for large-scale systems of nonlinear equations [17]. Therefore, combining the L-M-BFGS update (8) and (4), we design the L-M-BFGS TR subproblem model as follows where △ = c ‖b(x )‖.

This paper is organized as follows. In the next section, the algorithm is discussed. In Section 3, the global convergence is established. Finally, the numerical results are reported in Section 4.

The L-M-BFGS Trust-Region Method

The algorithm is given as follows.

Algorithm 1.

Initial: Given the constants ρ, c ∈ (0, 1), p = 0, ϵ > 0, x 0 ∈ ℜ, H 0 ∈ ℜ × ℜ is symmetric and positive definite. Set k: = 0;

Step 1: If ‖b ‖ < ϵ, stop;

Step 2: Solve the TR subproblem (9) with △ = △ to obtain d ;

Step 3: Calculate the following ratio of the actual reduction over the predicted reduction: If then we let p = p + 1, go to Step 2. Otherwise, go to Step 4;

Step 4: Set x = x + d , y = b − b . Use (8) to generate the updated matrix H with positive definiteness.

Step 5: Let k ≔ k + 1, p = 0, and go to Step 1.

Remark 1. (i) In Algorithm 1, the “Step 2-Step 3-Step 2” procedure is called the inner cycle.

(ii) It is well known that the positive definiteness of the update matrix H is very important when analyzing the convergence of the algorithm. Byrd et al. [27] showed that the limited-memory BFGS matrix has this property if the curvature (s ) y > 0 is satisfied. Similarly, Powell [28] proposed that y should be as follows: where B is an approximation of ∇b(x ) and Therefore, Step 4 of Algorithm 1 is reasonable.

To compare with the normal TR method, another algorithm is given and we call it Algorithm 2.

Algorithm 2.(BFGS Trust-Region Method)

Step 2: Solve the TR subproblem (4) with △ = △ to obtain d ;

Step 4: Set x = x + d , y = b − b . Use (5) to generate the updated matrix with positive definiteness.

Remark 2. It is easy to see that the difficult step of these two algorithms is the Step 2. The following Dogleg method [29] is used to solve the TR subproblem models (4) and (9) to obtain d .

At the kth iteration:

set ϖ = H b(x );

If ‖ϖ ‖ ≤ △, let d = −ϖ ; Otherwise, set

If ς ≥ 0 and ς ≤ 1 hold, set d = ςϱ; If ς ≥ 1 and ς ≤ 2 hold, set d = ϱ + (ς − 1)(τ − ϱ).

Convergence Analysis

This section will provide some convergence results under the following assumptions.

Assumption A (i) Let the level set Ω be bounded, and let b(x) be twice continuously differentiable on an open convex set Ω1 containing Ω.

(ii) The relation holds.

Assumption A(i) implies that there exists M 1 > 0 such that Similar to Moré [30], Yuan et al. [9], and Yuan and Sun [11], we can obtain the following lemma.

Lemma 0.1 Let the predicted reduction be and let d be the solution of (9). Then, we have

Proof. For any α ∈ [0, 1], by the property of the solution of (9) named d , we obtain Thus, the inequality holds. This completes the proof.

Lemma 0.2 Let d be the solution of (9), and let the actual reduction Ared (d ) be Suppose that Assumption A holds and that the sequence {x , d } is generated by Algorithm 1. Then, the inner cycle of Algorithm 1 will stop in a finite number of steps.

Proof. We will establish this lemma by contradiction. Suppose that the inner cycle of Algorithm 1 infinitely circles at x , namely, p → ∞, and c → 0 hold. It is easy to obtain ‖b ‖ ≥ ϵ, or the algorithm stops. Thus, we can conclude that ‖d ‖ ≤ △ = c ‖b ‖ → 0 is true.

By the definitions of d , Ared (d ), Pred (d ), and Assumption A, we obtain where the second equality follows Taylor’s formula. It follows from Lemma 0.1 and 7) that Thus, if p is sufficiently large, we have This relation contradicts The proof is complete.

The above lemma shows that the given algorithm is well-defined.

Theorem 0.1 Let Assumption A hold, and let {x } be generated by Algorithm 1. Then, {x } ⊂ Ω, and {θ(x )} converges. In particular, Algorithm 1 either stops in a finite number of steps or generates an infinite sequence {x } such that

Proof. By the definition of Algorithm 1, we obtain Combining this with Lemma 0.1 generates Thus, {x } ⊂ Ω holds. Considering θ(x ) ≥ 0, we then conclude that {θ(x )} converges.

If Algorithm 1 does not stop after a finite number of steps, again by (9) and θ(x ) ≥ 0, we easily conclude that which shows that b(x ) → 0, k → ∞. The proof is complete.

Theorem 0.1 proves that the sequence {x } of Algorithm 1 converges such that ‖b(x )‖ → 0 without the assumption that ∇b(x*) is nondegenerate, where x* is a cluster point of {x }.

Numerical Results

This section will report numerical results obtained using Algorithm 1. The test functions with initial points x 0 are listed as follows:

Function 1. Exponential function 2 Initial guess:

Function 2. Trigonometric function Initial guess:

Function 3. Singular function Initial guess: x 0 = (1, 1, ⋯, 1).

Function 4. Logarithmic function Initial guess: x 0 = (1, 1, ⋯, 1).

Function 5. Broyden Tridiagonal function [[31], pp. 471–472] Initial guess: x 0 = (−1, −1, ⋯, −1).

Function 6. Trigexp function Initial guess: x 0 = (0, 0, ⋯, 0).

Function 7. Strictly convex function 1 [[32], p. 29]

b(x) is the gradient of Initial guess:

Function 8. Variable dimensional function Initial guess:

Function 9. Discrete boundary value problem [33]. Initial guess: x 0 = (h(h − 1), h(2h − 1), ⋯, h(nh − 1)).

Function 10. The discretized two-point boundary value problem that is similar to the problem in [34] when A is the n × n tridiagonal matrix given by and F(x) = (F 1(x), F 2(x), …, F (x)), with F (x) = sin x − 1, i = 1, 2, …, n, and x = (50, 0, 50, 0, ⋯).

All codes for the experiments were written in MATLAB r2009a and run on a PC with an E5507@2.27 GHz CPU and 2.99 GB of memory using the Windows XP operating system. The parameters were set as ρ = 0.0001, ϵ = 10−5, c = 0.1, γ = 0.7, where B 0 and H 0 are unit matrices and m = 6. In the inner iterations of Algorithm 1 and Algorithm 2, the trial step will be accepted when p > 5. We also stop the program if the number of iterations is larger than one thousand or when θ(x) < 10−5. The columns of the tables have the following meaning:

the dimension of the problem.

the number of norm function evaluations.

the total number of iterations.

the normal value of θ(x) when the program stops.

fails to find the final values of θ(x).

From Table 1, it is clear that the new algorithm performs quite well in terms of solving these problems and that the dimension does not have an obvious influence on the performance of the presented method. Algorithm 2 can also successfully solve some of these problems. However, Algorithm 2 fails to solve problems 3 and 6.

Table 1

Numerical results using Algorithms 1 and 2.

Algorithm 1				Algorithm 2
Functions	Dim	NI/NG	GN	NI/NG	GN
1	800	0/1	8.348973e-006	0/1	8.348973e-006
	1000	0/1	6.676674e-006	0/1	6.676674e-006
	2000	0/1	3.335834e-006	0/1	3.335834e-00
2	800	7/18	1.766154e-006	8/19	4.764396e-006
	1000	7/18	1.578116e-006	8/19	4.228277e-006
	2000	6/17	8.507453e-006	8/19	2.609172e-006
3	800	82/108	5.997218e-006	1000/1716	3.453004e+018
	1000	68/89	8.159424e-006	1000/1636	3.372369e+019
	2000	95/152	8.892174e-006	1000/1701	1.757490e+022
4	800	5/6	6.440159e-008	5/6	6.440159e-008
	1000	5/6	7.954854e-008	5/6	7.954854e-008
	2000	5/6	1.553487e-007	5/6	1.553487e-007
5	800	61/67	9.546109e-006	47/83	7.073607e-006
	1000	58/64	9.424256e-006	35/66	5.881295e-006
	2000	61/67	9.412218e-006	44/75	8.483566e-006
6	800	62/73	6.368574e-008	6/22	NaN
	1000	59/75	1.493209e-007	6/22	NaN
	2000	71/92	2.616014e-007	6/22	NaN
7	800	6/7	1.665516e-009	5/6	1.718325e-006
	1000	6/7	2.070362e-009	5/6	2.139904e-006
	2000	6/7	4.094910e-009	5/6	4.247951e-006
8	800	1/2	0.000000e+000	1/2	0.000000e+000
	1000	1/2	0.000000e+000	1/2	0.000000e+000
	2000	1/2	0.000000e+000	1/2	0.000000e+000
9	800	2/3	3.679023e-006	4/10	1.336824e-006
	1000	2/3	2.358640e-006	4/10	9.011723e-007
	2000	2/3	5.917002e-007	2/3	3.637130e-006
10	800	15/21	1.212419e-007	32/40	6.985962e-006
	1000	15/21	1.391409e-007	34/42	5.778079e-006
	2000	15/21	1.824225e-007	33/44	5.520880e-006

To clearly demonstrate these algorithms’ efficiencies, a tool developed by Dolan and Moré [35] is used to analyze them. In this tool, the so-called (cumulative) distribution function ϱ is defined by where p is a problem, s is its solver, n is the number of problems, n is the number of existing solvers, P is the test set, κ denotes the spending time or the iteration number, and π is the performance ratio given by where S is the set of solvers. According to the above formulas and the rules given in [35], the performance profile plot of Table 1 is given in the following two figures.

Figs 1 and 2 show that the performance of these methods is a function of NI and NG, respectively, in Table 1. It is clear that the given method performs better because it has a higher probability of being the optimal solver. Fig 1 shows that Algorithm 1 is superior to Algorithm 2 for κ ≤ 14, approximately. From Fig 2, it is clear that the given method can successfully solve 100% of the test problems at κ ≈ 5 and that Algorithm 2 successfully solves the test problems at κ ≈ 17. Overall, the numerical results are interesting.

Fig 1

Performance profiles of these methods(NI).

Fig 2

Performance profiles of these methods(NG).

Supporting Information:Legends of Figs 1 and 2

Dolan and Moré [35] developed a new tool to analyze the efficiency of algorithms. They introduced the notion of a performance profile for evaluating and comparing the performance of the set of solvers S on a test set P. Assuming that n solvers and n problems exist, for each problem p and solver s, they defined κ as the computing time (the number of function evaluations or some other metric) required to solve problem p using solver s.

Requiring a baseline for comparison, they compared the performance on problem p using solver s with the best performance by any solver on this problem, giving the the performance ratio Suppose that a parameter π ≥ π for all p, s is chosen, and π = π if and only if solver s does not solve problem p.

The performance of solver s on any given problem might be of interest. More importantly, one would like to obtain an overall assessment of the performance of the solver. With this motivation, these authors defined Specifically, ϱ(κ) is the probability for solver s ∈ S that a performance ratio π is within a factor κ ∈ ℜ of the best possible ratio. Then, function ϱ is the (cumulative) distribution function for the performance ratio. The performance profile ϱ : ℜ ↦ [0, 1] for a solver is a nondecreasing, piecewise constant function and is continuous from the right at each breakpoint. The value of ϱ(1) is the probability that the solver would outperform the remaining solvers.

According to the above rules, we know that one solver whose performance profile plot is on the top right outperforms the remaining solvers. Similar legends can be found in [10, 36–38].

In this paper, a TR model combining with the L-M-BFGS technique is presented for nonlinear equations and the global convergence is established. Fact, this work is an extension of the paper [9] and the difference between these two papers is the update matrix: this paper chooses the L-M-BFGS formula and the paper [9] chooses the normal BFGS formula.

Since the L-M-BFGS technique requires minimal storage and is suitable for certain classes of large scale optimization problems, so we extend it to large scale nonlinear equations. Numerical results turn out the given algorithm is competitive to similar method for large-scale nonlinear equations.

It is well known that the conjugate gradient methods are effective techniques for large-scale optimization problems since they do not need the matrix information, which motivates us to use the conjugate gradient technique for large scale nonlinear equations. We think the numerical performance is also very interesting.