Scaled Diagonal Gradient-Type Method with Extra Update for Large-Scale Unconstrained Optimization

Abstract

We present a new gradient method that uses scaling and extra updating within the diagonal updating for solving unconstrained optimization problem. The new method is in the frame of Barzilai and Borwein (BB) method, except that the Hessian matrix is approximated by a diagonal matrix rather than the multiple of identity matrix in the BB method. The main idea is to design a new diagonal updating scheme that incorporates scaling to instantly reduce the large eigenvalues of diagonal approximation and otherwise employs extra updates to increase small eigenvalues. These approaches give us a rapid control in the eigenvalues of the updating matrix and thus improve stepwise convergence. We show that our method is globally convergent. The effectiveness of the method is evaluated by means of numerical comparison with the BB method and its variant.

1. Introduction

In this paper, we consider the unconstrained optimization problem min𝑓(𝑥),𝑥∈𝑅^𝑛, (1) where 𝑓(𝑥) is a continuously differentiable function from 𝑅^𝑛 to 𝑅. Given a starting point 𝑥₀, using notations 𝑔_𝑘=𝑔(𝑥_𝑘)=∇𝑓(𝑥_𝑘) and 𝐵_𝑘 as an approximation to the Hessian 𝐺_𝑘=[∇²𝑓(𝑥_𝑘)], the quasi-Newton-based methods for solving (1) are defined by the iteration 𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵_𝑘⁻¹𝑔_𝑘,𝑘=0,1,2,…, (2) where the stepsize 𝛼_𝑘 is determined through an appropriate selection. The updating matrix 𝐵_𝑘 is usually required to satisfy the quasi-Newton equation 𝐵_𝑘𝑠_𝑘−1=𝑦_𝑘−1, (3) where 𝑠_𝑘−1=𝑥_𝑘−𝑥_𝑘−1 and 𝑦_𝑘−1=𝑔_𝑘−𝑔_𝑘−1. One of the widely used quasi-Newton method to solve general nonlinear minimization is the BFGS method, which uses the following updating formula: 𝐵_𝑘+1=𝐵_𝑘−𝐵_𝑘𝑠_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘+𝑦_𝑘𝑦^𝑇_𝑘𝑠^𝑇_𝑘𝑦_𝑘. (4) On the numerical aspect, this method supersedes most of the optimization methods; however, it needs 𝑂(𝑛²) storage which makes it unsuitable for large-scale problems.

On the other hand, an ingenious stepsizes selection for gradient method was proposed by Barzilai and Borwein [1] in which the updating scheme is defined by 𝑥_𝑘+1=𝑥_𝑘−𝐷_𝑘⁻¹𝑔_𝑘, (5) where 𝐷_𝑘=(1/𝛼_𝑘)𝐼 and 𝛼_𝑘=𝑠^𝑇_𝑘−1𝑠_𝑘−1/𝑠^𝑇_𝑘−1𝑦_𝑘−1.

Since that, the study of new effective methods in the frame of BB-like gradient methods becomes an interesting research topic for a wide range of mathematical programming; for example, see [2–10]. However, it is well known that BB method cannot guarantee a descent in the objective function at each iteration and the extent of the nonmonotonicity depends in some way on the size of the condition number of objective function [11]. Therefore, the performance of BB method is greatly influenced by the condition of the problem (particularly, condition number of the Hessian matrix). Some new fixed stepsizes gradient-type methods of BB kind are proposed by [12–16] to overcome these difficulties. In contrast with the BB approach in which the stepsize is computed by means of a simple approximation of the Hessian in the form of scalar multiple of identity, these proposed methods consider approximation of the Hessian and its inverse in diagonal matrix form based on the weak secant equation and quasi-cauchy relation, respectively (for more details see [15, 16]). Though these diagonal updating methods are efficient, their performance can be greatly affected by solving ill-conditioned problems. Thus, there is room for improve on the quality of the diagonal updates formulation. Since methods as described in [15, 16] have useful theoretical and numerical properties, it is desirable to derive a new and more efficient updating frame for general functions. Therefore our aim is to improve the quality of diagonal updating when it is poor in approximating Hessian.

This paper is organized as follows. In the next section, we describe our motivation and propose our new-gradient type method. The global convergence of the method under mild assumption will be established in Section 3. Numerical evidence of the vast improvements due to the new approach is given in Section 4. Finally, conclusion is made in the last section.

2. Scaling and Extra Updating

Assume that 𝐵_𝑘 is positive definite, and let {𝑦_𝑘} and {𝑠_𝑘} be two sequences of 𝑛-vectors such that 𝑦^𝑇_𝑘𝑠_𝑘>0 for all 𝑘. Because it is usually difficult to satisfy the quasi-Newton equation (3) with a nonsingular 𝐵_𝑘+1 of the diagonal form, one can consider satisfying it in some directions. If we project the quasi-Newton equation (3) (also called the secant equation), in a direction 𝜐 such that 𝑦^𝑇_𝑘𝜐≠0, then it gives 𝑠^𝑇_𝑘𝐵_𝑘+1𝜐=𝑦^𝑇_𝑘𝜐. (6) If 𝜐=𝑠_𝑘 is chosen, it leads to the so-called weak-secant relation, 𝑠^𝑇_𝑘𝐵_𝑘+1𝑠_𝑘=𝑦^𝑇_𝑘𝑠_𝑘. (7)

Under this weak-secant equation, [15, 16] employ variational technique to derive updating matrix that approximates the Hessian matrix diagonally. The resulting update is derived to be the solution of the following variational problem: 1min2𝐵_𝑘+1−𝐵_𝑘²_𝐹s.t.𝑠^𝑇_𝑘𝐵_𝑘+1𝑠_𝑘=𝑠^𝑇_𝑘𝑦_𝑘,𝐵_𝑘+1isdiagonal (8) and gives the corresponding solution 𝐵_𝑘+1 as follows: 𝐵_𝑘+1=𝐵_𝑘+𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘, (9) where 𝐸_𝑘=diag(𝑠²_𝑘,1,𝑠²_𝑘,2,…,𝑠²_𝑘,𝑛), 𝑠_𝑘,𝑖 is the 𝑖th component of the vector 𝑠_𝑘, and tr denotes the trace operator.

Note that when 𝑠^𝑇_𝑘𝑦_𝑘<𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘, the resulting 𝐵_𝑘+1 is not necessarily positive definite and it is not appropriate for use within a quasi-Newton-based algorithm. Thus, it is desirable to propose a technique to measure the quality of 𝐵_𝑘 in terms of its Rayleigh quotient and try to find a way to improve “poor” quality 𝐵_𝑘 before calculating 𝐵_𝑘+1. For this purpose, it will be useful to propose, at first quality a criterion to distinguish between poor, and acceptable quality of 𝐵_𝑘.

Let us begin by considering the curvature of an objective function, 𝑓 in direction 𝑠_𝑘, which is represented by 𝑠^𝑇_𝑘𝐺_𝑘𝑠_𝑘=𝑠^𝑇_𝑘𝑦_𝑘, (10) where 𝐺_𝑘=∫¹₀∇²𝑓(𝑥_𝑘+𝑡𝑠_𝑘)𝑑𝑡 is the average Hessian matrix along 𝑠_𝑘. Since it is not practical to compute the eigenvalue of the Hessian matrix in each iteration, we can estimate its relative size on the basis of the scalar 𝜌_𝑘=𝑠^𝑇_𝑘𝐺_𝑘𝑠_𝑘/𝑠^𝑇_𝑘𝑠_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘/𝑠^𝑇_𝑘𝑠_𝑘=𝑠^𝑇_𝑘𝐺_𝑘𝑠_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘=𝑠^𝑇_𝑘𝑦_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘. (11) If 𝜌_𝑘>1, it implies that the eigenvalues of 𝐵_𝑘 approximated by its Rayleigh are relatively small compared to those of the local Hessian matrix at 𝑥_𝑘. In this condition, we find that the strategy of extra update [17] seems to be useful for improving the quality of 𝐵_𝑘 by rapidly increasing its eigenvalues up to those of the actual Hessian relatively. This is done by updating 𝐵_𝑘 twice to obtain 𝐵_𝑘+1,2: 𝐵_𝑘+1,1=𝐵_𝑘+𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘,𝐵_𝑘+1,2=𝐵_𝑘+1,1+𝑠^𝑇_𝑘−1𝑦_𝑘−1−𝑠^𝑇_𝑘−1𝐵_𝑘+1,1𝑠_𝑘−1tr𝐸²_𝑘−1𝐸_𝑘−1, (12)(13) and use it to obtain, finally, the updated 𝐵_𝑘+1: 𝐵_𝑘+1=𝐵_𝑘+1,2+𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘+1,2𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘. (14)

On the other hand, when 𝜌_𝑘<1, it implies that the eigenvalue of 𝐵_𝑘 represented by its Rayleigh is relatively large and we have 𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘<0. In this case, we should suggest a useful strategy to encounter this drawback. As we reviewed before, the updating scheme may generate nonpositive definite 𝐵_𝑘+1 when 𝐵_𝑘 has large eigenvalues relative to those possible values of 𝐺_𝑘, that is, when 𝑠^𝑇_𝑘𝑦_𝑘<𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘. On the contrary, this argument disappears when the eigenvalues of 𝐵_𝑘 are small (i.e., 𝑠^𝑇_𝑘𝑦_𝑘>𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘). This suggests that the scaling should be made to scale down 𝐵_𝑘, that is, choosing 𝜌_𝑘<1 only when 𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘<0 and take 𝜌_𝑘=1, whenever 𝑠^𝑇_𝑘𝑦_𝑘>𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘. Combining these two arguments, we choose the scaling parameter 𝜌_𝑘 such that 𝛾_𝑘=min𝜌_𝑘,1. (15) This scaling resembles the Al-Baali [18] scaling that is applied within the Broyden family. Because the value of 𝛾_𝑘 is always <1, then by incorporating the scaling to 𝐵_𝑘, it decreases the large eigenvalues of 𝐵_𝑘 constantly, and consequently we can keep positive definiteness of 𝐵_𝑘+1 (since 𝑠^𝑇_𝑘𝑦_𝑘>0), which is an important property in descent method. In this case, the following updating: 𝐵_𝑘+1=𝛾_𝑘𝐵_𝑘+𝑠^𝑇_𝑘𝑦_𝑘−𝛾_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘, (16) will be used. To this end, we have the following general updating scheme for 𝐵_𝑘+1: 𝐵_𝑘+1=𝛾_𝑘𝐵_𝑘+𝑠^𝑇_𝑘𝑦_𝑘−𝛾_𝑘𝑠^𝑇_𝑘𝐵_𝑘𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘;if𝜌_𝑘𝐵≤1,_𝑘+1,2+𝑠^𝑇_𝑘𝑦_𝑘−𝑠^𝑇_𝑘𝐵_𝑘+1,2𝑠_𝑘tr𝐸²_𝑘𝐸_𝑘;if𝜌_𝑘>1, (17) where 𝐵_𝑘+1,2 and 𝛾_𝑘 are given by (13) and (15), respectively.

An advantage of using (17) is that the positive definiteness of 𝐵_𝑘+1 can be guaranteed in all iterations. This property is not exhibited in the other diagonal updating formula such as those in [15, 16]. Note that there is no extra storage required to impose our strategy and the cost of computing is also not increased significantly throughout the entire iteration. Now we can state the steps of our new diagonal-gradient method algorithm with the safeguarding strategy for monotonicity as follows.

2.1. ESDG Algorithm

Step 1. Choose an initial point 𝑥₀∈𝑅^𝑛 and a positive definite matrix 𝐵₀=𝐼. Let 𝜃∈(1,2). Set 𝑘∶=0.

Step 2. Compute 𝑔_𝑘. If ‖𝑔_𝑘‖≤𝜖, stop.

Step 3. If 𝑘=0, set 𝑥₁=𝑥₀−𝑔₀/‖𝑔₀‖.

Step 4. Compute 𝑑_𝑘=−𝐵_𝑘⁻¹𝑔_𝑘, and calculate 𝛼_𝑘>0 such that the following condition holds: 𝑓(𝑥_𝑘+1)≤𝑓_𝑘^max+𝜎𝛼_𝑘𝑔^𝑇_𝑘𝑑_𝑘 where 𝑓_𝑘^max=max{𝑓(𝑥_𝑘),𝑓(𝑥_𝑘−1)} and 𝜎∈(0,1) is a given constant.

Step 5. If 𝑘≥1, let 𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵_𝑘⁻¹𝑔_𝑘 and compute 𝜌_𝑘 and 𝛾_𝑘 by (11) and (15), respectively. If 𝜌_𝑘<𝜃 then update 𝐵_𝑘+1 by (16).

Step 6. If 𝜌_𝑘≥𝜃 then compute 𝐵_𝑘+1,1 and 𝐵_𝑘+1,2 by (12), (13), respectively, and then update as defined 𝐵_𝑘+1 (14).

Step 7. Set 𝑘∶=𝑘+1, and return to Step 2.

In Step 4, we employ the nonmonotone line search of [19, 20] to ensure the convergence of the algorithm. However, some other line search strategies may also be used.

3. Convergence Analysis

This section is devoted to study the convergence behavior of ESDG method. We will establish the convergence of the ESDG algorithm when applied to the minimization of a strictly convex function. To begin, we give the convergence result, which is due to Grippo et al. [21] for the step generated by the nonmonotone line search algorithm. Here and elsewhere, ‖⋅‖ denotes the Euclidean norm.

Theorem 1. Assume that 𝑓 is a strictly convex function and its gradient 𝑔 satisfies the Lipschitz condition. Suppose that the nonmonotone line search algorithm is employed in a case that the steplength, 𝛼_𝑘, satisfies 𝑓𝑥_𝑘+1≤𝑓_𝑘^max+𝜎𝛼_𝑘𝑔^𝑇_𝑘𝑑_𝑘, (18) where 𝑓_𝑘^max=max{𝑓(𝑥_𝑘),𝑓(𝑥_𝑘−1),…,𝑓(𝑥_𝑘−𝑚)}, with 𝑚≤𝑘 and 𝜎∈(0,1), and the search direction 𝑑_𝑘 is chosen to obey the following conditions. There exist positive constants 𝑐₁ and 𝑐₂ such that −𝑔^𝑇_𝑘𝑑_𝑘≥𝑐₁𝑔_𝑘²,𝑑_𝑘≤𝑐₂𝑔_𝑘, (19) for all sufficiently large 𝑘. Then the iterates 𝑥_𝑘 generated by the nonmonotone line search algorithm have the property that liminf_𝑘→∞𝑔_𝑘=0. (20)

To prove that the ESDG algorithm is globally convergent, it is sufficient to show that the sequence {‖𝐵_𝑘‖} generated by (17) is bounded both above and below, for all finite 𝑘 so that its associated search direction satisfies condition (19). Since 𝐵_𝑘 is diagonal, it is enough to show that each element of 𝐵_𝑘, say 𝐵_𝑘^(𝑖), 𝑖=1,…,𝑛, is bounded above and below by some positive constants. The following theorem gives the boundedness of {‖𝐵_𝑘‖}.

Theorem 2. Assume that 𝑓 is strictly convex function where there exist positive constants 𝑚 and 𝑀 such that 𝑚‖𝑧‖²≤𝑧^𝑇∇²𝑓(𝑥)𝑧≤𝑀‖𝑧‖², (21) for all 𝑥,𝑧∈𝑅^𝑛. Let {‖𝐵_𝑘‖} be a sequence generated by the ESDG method. Then ‖𝐵_𝑘‖ is bounded above and below for all finite 𝑘, by some positive constants.

Proof. Let 𝐵_𝑘^(𝑖) be the 𝑖th element of 𝐵_𝑘. Suppose that 𝐵₀ is chosen such that 𝜔₁≤𝐵₀^(𝑖)≤𝜔₂,𝑖=1,…,𝑛 where 𝜔₁,𝜔₂ are some positive constants. It follows from (17) and the definition of 𝛾 in (15) that we have 𝐵₁=𝜌₀𝐵₀,if𝜌₀𝐵≤1,_1,2+𝑠^𝑇₀𝑦₀−𝑠^𝑇₀𝐵_1,2𝑠₀tr𝐸²₀𝐸₀,if𝜌₀>1, (22) where 𝐵_1,1=𝐵₀+𝑠^𝑇₀𝑦₀−𝑠^𝑇₀𝐵₀𝑠₀tr𝐸²₀𝐸₀,𝐵_1,2=𝐵_1,1+𝑠^𝑇₀𝑦₀−𝑠^𝑇₀𝐵_1,1𝑠₀tr𝐸²₀𝐸₀. (23) Moreover, by (21) and (11), we obtain 𝑚𝑠_𝑘²≤𝑠^𝑇_𝑘𝑦_𝑘≤𝑀𝑠_𝑘²,∀𝑘. (24)Case??1. When 𝜌₀≤1: by (24), one can obtain 𝑚𝜔₂≤𝜌₀=𝑠^𝑇₀𝑦₀𝑠^𝑇₀𝐵₀𝑠₀≤𝑀𝜔₁. (25) Thus, it implies that 𝑚𝜔₁/𝜔₂≤𝐵₁^(𝑖)=𝜌₀𝐵₀^(𝑖)≤𝑀𝜔₂/𝜔₁.
Case??2. When 𝜌₀>1: from (3), we have 𝐵^(𝑖)_1,1=𝐵₀^(𝑖)+𝑠^𝑇₀𝑦₀−𝑠^𝑇₀𝐵₀𝑠₀tr𝐸²₀𝑠²_0,𝑖. (26) Because 𝜌₀>1 also implies that 𝑠^𝑇₀𝑦₀−𝑠^𝑇₀𝐵₀𝑠₀>0, using this fact and (24) give 𝐵₀^(𝑖)≤𝐵^(𝑖)_1,1≤𝐵₀^(𝑖)+𝑀−𝜔₁𝑠₀²tr𝐸²₀𝑠²_0,𝑖. (27) Let 𝑠_0,𝑀 be the largest component in magnitude of 𝑠₀, that is, 𝑠²_0,𝑖≤𝑠²_0,𝑀,forall𝑖. Then it follows that ‖𝑠₀‖²≤𝑛𝑠²_0,𝑀, and (27) becomes 𝜔₁≤𝐵^(𝑖)_1,1≤𝜔₂+𝑛𝑀−𝜔₁tr𝐸²₀𝑠⁴_0,𝑀≤𝜔₂+𝑛𝑀−𝜔₁. (28) Using (28) and the same argument as previously mentioned, we can also show that 𝜔₁≤𝐵^(𝑖)_1,2≤𝜔₂+𝑛𝑀−𝜔₁+𝑛𝑀−𝜔₂+𝑛𝑀−𝜔₁. (29) Hence in both cases, 𝐵₁^(𝑖) is bounded above and below, by some positive constants. Since the upper and lower bounds for 𝐵₁^(𝑖) are independent of 𝑘, respectively, we can proceed by using induction to show that 𝐵_𝑘^(𝑖) is bounded, for all finite 𝑘.

4. Numerical Results

In this section we present the results of numerical investigation for ESDG method on different test problems. We also compare the performance of our new method with that of the BB method and that of MDGRAD method which is implemented using SMDQN of [22] with a same nonmonotone strategy as the ESDG method. Our experiments are performed on a set of 20 nonlinear unconstrained problems with dimensions ranging from 10 to 10⁴ (Table 1).

Table 1 

Test problem and its dimension.

These test problems are taken from [23, 24]. The codes are developed with Matlab 7.0. All runs are performed on a PC with Core Duo CPU. For each test, the termination condition is ‖𝑔(𝑥_𝑘)‖≤10⁻⁴. The maximum number of iterations is set to 1000.

Figures 1 and 2 show the efficiency of ESDG method when compared to MDGRAD and BB methods. Note that ESDG method increases the efficiency of Hessian approximation devoid of increasing the number of storages. Figure 2 also shows the implementation of the ESDG method with BB and MDGRAD methods using the CPU time as a measure. This figure shows that ESDG method is again faster than MDGRAD method in most problems and requires reasonable time to solve large-scale problems when compares to the BB method. Finally, we can conclude that our experimental comparisons indicate that our extension is very beneficial to the performance.

Figure 1 

Performance profile based on iterations.

Figure 2 

Performance profile based on CPU time per iteration.

5. Conclusion

We have presented a new diagonal gradient method for unconstrained optimization. Numerical study of the proposed method when compared with BB and MDGRAD methods is also performed. Based on our numerical experiments, we can conclude that ESDG method is significantly preferable compared to the BB and MDGRAD methods. Particularly, the ESDG method is proven to be a good option for large-scale problems when high-memory locations are required. In view of the remarkable performance of ESDG method, globally converged and with only 𝑂(𝑛) storage, we can expect that our proposed method would be useful for unconstrained large-scale optimization problems.