Diagonal Hessian Approximation for Limited Memory Quasi-Newton via Variational Principle

Abstract

This paper proposes some diagonal matrices that approximate the (inverse) Hessian by parts using the variational principle that is analogous to the one employed in constructing quasi-Newton updates. The way we derive our approximations is inspired by the least change secant updating approach, in which we let the diagonal approximation be the sum of two diagonal matrices where the first diagonal matrix carries information of the local Hessian, while the second diagonal matrix is chosen so as to induce positive definiteness of the diagonal approximation at a whole. Some numerical results are also presented to illustrate the effectiveness of our approximating matrices when incorporated within the L-BFGS algorithm.

1. Introduction

Our investigation begins by seeking effective ways to diagonally scale an identity matrix, which is often used to initialize the L-BFGS method. For this purpose, it is useful to state the weak quasi-Newton equation: 𝑠^𝑇_𝑘𝑀_𝑘+1𝑠_𝑘=𝑠^𝑇_𝑘𝑦_𝑘, (1) where the 𝑛-dimensional vector 𝑠_𝑘=𝑥_𝑘+1−𝑥_𝑘 denotes the step corresponding to two different points 𝑥_𝑘 and 𝑥_𝑘+1, 𝑦_𝑘=𝑔_𝑘+1−𝑔_𝑘 denotes the gradient change corresponding to the gradients 𝑔_𝑘 and 𝑔_𝑘+1 at the two points, 𝑀_𝑘+1 is a full 𝑛×𝑛 matrix that approximates the Hessian of 𝑓. Both 𝑠 and 𝑦 are used explicitly and 𝑂(𝑛²) storage is required for the matrix 𝑀_𝑘+1. If 𝑀 is further chosen to be a diagonal matrix, say 𝐷_𝑘+1, then we can establish the so-called quasi-Cauchy relation denwol 1993, [1] by which 𝑠^𝑇_𝑘𝐷_𝑘+1𝑠_𝑘=𝑠^𝑇_𝑘𝑦_𝑘. (2) Different from the standard quasi-Newton which requires 𝑂(𝑛²) of storage, here only 𝑂(𝑛) storage is required to store the diagonal update that satisfies the quasi-Cauchy relation. In addition, we suppose the matrix 𝐷_𝑘+1 to be positive definite by which it is able to define a metric.

The first example of diagonal updating that satisfies the quasi-Cauchy relation is the well-known Oren-Luenberger OREN1974 scaling matrix, given by 𝐷_𝑘+1=𝑠^𝑇_𝑘𝑦_𝑘𝑠^𝑇_𝑘𝑠_𝑘𝐼, (3) where 𝐼 is the identity matrix. Expression (3) would be obtained from the quasi-Cauchy relation with the further restriction where the diagonal matrix is a scalar multiple of the identity matrix. Therefore, scaling matrices that are derived from the quasi-Cauchy relation are a natural generalization of Oren-Luenberger scaling.

In general, a procedure to obtain diagonal updating formulae via quasi-Cauchy relation can be summarized as follow. Suppose that 𝐷_𝑘>0 is a positive definite diagonal matrix and 𝐷_𝑘+1 is the updated version of 𝐷, which is also diagonal. Then the updated 𝐷_𝑘+1 is required to satisfy the quasi-Cauchy relation. It is also essential to require the deviation between 𝐷_𝑘 and 𝐷_𝑘+1 being minimized under some variational principle, which in return will encourage numerical stability. Very often the Frobenius matrix norm is used to measure this deviation. As noted earlier, a diagonal matrix uses the same computer storage as a vector, resulting in potential use in limited memory algorithms.

This paper is structured as follows. In Section 2, we consider a briefly review on the diagonal updating via quasi-Cauchy relation and the preconditioning strategy. New diagonal initial approximation for L-BFGS method is introduced in this section as well. It follows that in Section 3, the convergence properties of the proposed L-BFGS methods are investigated. Numerical results are presented in Section 4, on a large set of unconstrained minimization problems, mainly from the collections of Moré et al. [2], CUTE [3], and Toint [4]. In Section 5, through the numerical results, we do the discussion regarding the new diagonal initial approximating which matrix is more effective at improving the performance of the L-BFGS update compared to the standard preconditioner in the literature. Finally, Section 6 ends the paper by presenting a summary and conclusion.

2. Diagonal Updating via Quasi-Cauchy Relation

The performance of the L-BFGS method is depending on a good approximations of the actual Hessian. In the basic implementation of the L-BFGS method, the correction pairs are as follows: 𝑠_𝑖,𝑦_𝑖for𝑘−𝑚+1≤𝑖≤𝑘, (4) to correct 𝐻⁰_𝑘. The choices of 𝐻⁰_𝑘 often influences the behaviour of the method. It is worth the investigation on the choices of 𝐻⁰_𝑘.

Throughout this section, when we mention a direct initial matrix, we mean a matrix that is a rough approximation to the Hessian; otherwise an initial matrix is an approximation to the inverse of Hessian.

Our approach is inspired by [5, 6] which employed a variational technique that is analogue to the one used to derive the Powell Symmetric Broyden (PSB) quasi-Newton update (see, e.g, Dennis and Schnabel [7]). The resulting update for approximating the Hessian matrix diagonally is derived as follows: 𝐷_𝑘=𝐷_𝑘−1+𝑠^𝑇_𝑘−1𝑦_𝑘−1−𝑠^𝑇_𝑘−1𝐷_𝑘−1𝑠_𝑘−1tr𝐸²_𝑘−1𝐸_𝑘−1, (5) where 𝐸_𝑘−1=diag(𝑠²_𝑘−1,1,𝑠²_𝑘−1,2,…,𝑠²_{𝑘−1,𝑛}), 𝑠_{𝑘−1,𝑖} denotes the 𝑖th component of the vector 𝑠_𝑘−1, and tr is the trace operator.

Note that when 𝑠^𝑇_𝑘−1𝑦_𝑘−1<𝑠^𝑇_𝑘−1𝐷_𝑘−1𝑠_𝑘−1, the resulting 𝐷_𝑘+1 is not necessarily positive definite. Thus, like their counterpart of PSB update in the quasi-Newton setting, the foregoing update may suffer from the loss of positive definiteness and it is not appropriate for use within a quasi-Newton-based algorithm.

In this study, our approach in finding an efficient diagonal Hessian approximation is done through letting the diagonal approximating matrix be a combination of two diagonal matrices. This gives us a freedom to incorporate curvature information into one of these diagonal matrices, while the property of hereditary positive definiteness is carried over to the second matrix.

To begin, suppose that the Hessian matrix of an objective function has positive diagonal elements. Let us divide the Hessian matrix into two parts: ∇²𝑓(𝑥)=∇²𝑓_𝐴(𝑥)+∇²𝑓_𝐵(𝑥), (6) where ∇²𝑓_𝐴(𝑥) is a diagonal matrix consisting the diagonal entries of the Hessian and ∇²𝑓_𝐵(𝑥) would resemble the actual Hessian except that its diagonal entries are all zero. Thus, we intend to form two diagonal approximating matrices to approximate each part of the Hessian, respectively; that is, ∇²𝑓(𝑥)≈𝐷⁰=Ψ₁+Ψ₂. (7) Since it is assume that the entries of the actual Hessian are all positive, an excellent choice would be to let a positive definite diagonal matrix, says Ψ₁, to approximate ∇²𝑓_𝐴(𝑥). Meanwhile, ∇²𝑓_𝐵(𝑥) is expected to be dense and expensive to compute and would be approximated by Ψ₂. To preserve positive-definiteness, we introduce additional term in the form of 𝜃𝐼, that is, a Levenberg-Marquat-like step to maintain positive definite of Ψ₂ in a way that Ψ₂ is expressed as Ψ₂=𝜃𝐼+Ψ₃. Thus, the technique for calculating 𝐷⁰_𝑘 subjected to the weak-quasi-Newton relation is as follows.

Theorem 1. Suppose that 𝑠_𝑘−1≠0; then the optimal solution of the following minimization problem: 1min2Ψ₃²_𝐹,𝑠.𝑡.𝑠^𝑇_𝑘−1Ψ₁+𝜃𝐼+Ψ₃𝑠_𝑘−1=𝑠^𝑇_𝑘−1𝑦_𝑘−1,Ψ₃𝑖𝑠𝑑𝑖𝑎𝑔𝑜𝑛𝑎𝑙, (8) is given by Ψ₃=𝑠^𝑇_𝑘−1𝑦_𝑘−1−𝑠^𝑇_𝑘−1Ψ₁𝑠_𝑘−1−𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1tr𝐸²_𝑘−1𝐸_𝑘−1. (9)

Proof. Since the objective function in (8) is convex and the feasible set is also convex, then (8) has a unique solution. Its Lagrangian function is given by 𝐿Ψ₃=1,𝜇2Ψ₃²_𝐹+𝜇𝑠^𝑇_𝑘−1Ψ₃𝑠_𝑘−1+𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1+𝑠^𝑇_𝑘−1Ψ₁𝑠_𝑘−1−𝑠^𝑇_𝑘−1𝑦_𝑘−1, (10) where 𝜇 is the Lagrange multiplier associated with the constraint. Differentiating (10) with respect to each of elements of Ψ₃ and setting the results to zero yields, Ψ₃=−𝜇𝐸_𝑘−1. (11) Pre- and postmultiplying (11) by 𝑠_𝑘−1 and invoking the constraint, 𝑠^𝑇_𝑘−1Ψ₃𝑠_𝑘−1=𝑠^𝑇_𝑘−1𝑦_𝑘−1−𝑠^𝑇_𝑘−1Ψ₁𝑠_𝑘−1−𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1=−𝜇𝑠^𝑇_𝑘−1𝐸_𝑘−1𝑠_𝑘−1. (12) We have 𝑠𝜇=^𝑇𝑘−1Ψ₁𝑠_𝑘−1+𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1−𝑠^𝑇_𝑘−1𝑦_𝑘−1𝑠^𝑇_𝑘−1𝐸_𝑘−1𝑠_𝑘−1. (13) Finally, by substituting (13) into (11) and using the fact that 𝑠^𝑇_𝑘−1𝐸_𝑘−1𝑠_𝑘−1=tr(𝐸²_𝑘−1) yields (9), Ψ₃=−𝜇𝐸_𝑘−1𝑠=−^𝑇𝑘−1Ψ₁𝑠_𝑘−1+𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1−𝑠^𝑇_𝑘−1𝑦_𝑘−1tr𝐸²_𝑘−1𝐸_𝑘−1. (14)

A direct result of Theorem 1 leads to the following diagonal preconditioning formulation: 𝐷⁰_𝑘=Ψ₁+𝑠𝜃𝐼+^𝑇𝑘−1𝑦_𝑘−1−𝑠^𝑇_𝑘−1Ψ₁𝑠_𝑘−1−𝜃𝑠^𝑇_𝑘−1𝑠_𝑘−1tr𝐸²_𝑘−1𝐸_𝑘−1. (15) An analogue of the above approach can be used to derive the initial inverse Hessian approximation 𝐻⁰, which is more useful for algorithmic purpose. Once again by letting the inverse Hessian of an objective be separated by parts; ∇²𝑓(𝑥)⁻¹=𝐻⁽¹⁾+𝐻⁽²⁾, (16) then the initial approximating matrix 𝐻⁰ can be expressed as 𝐻⁰=Δ₁+𝜃𝐼+Δ₃. (17) Since it is our intention to derive an initial approximation for the L-BFGS method, then an excellent choice of Δ_𝑘+1,𝑖 for approximating the diagonal entries of the inverse Hessian would be the diagonalized inverse BFGS formula of Gilbert and Lemaréchal [8] updated from a multiple of identity matrix; that is, Δ_𝑘+1,𝑖=Δ_𝑘,𝑖+1𝑦^𝑇_𝑘𝑠_𝑘+𝑦^𝑇_𝑘Δ_𝑘,𝑖𝑦_𝑘𝑦^𝑇_𝑘𝑠_𝑘²𝑠_𝑘,𝑖²−2Δ_𝑘,𝑖𝑠_𝑘,𝑖𝑦_𝑘,𝑖𝑦^𝑇_𝑘𝑠_𝑘, (18) where we choose Δ_𝑘,𝑖=𝑦^𝑇_𝑘𝑠_𝑘/𝑦^𝑇_𝑘𝑦_𝑘 as the Oren-Luenberger scalar. Thus, diagonal matrix Δ₁ can be obtained in (18). Note that, to safeguard very small or very large ‖𝑠_𝑘−1‖, we impose the additional following condition: if 𝜀₂≤‖𝑠_𝑘−1‖≤𝜀₁, for some small and large positive 𝜀₁ and 𝜀₂, we set 𝐵₀=𝐼. Here, we used 𝜀₁=10⁸ and 𝜀₂=10⁻⁸.

By interchanging the role of 𝑠 and 𝑦 in Theorem 1, one can obtain the formula for 𝐻⁰ at step 𝑘 as follows: 𝐻⁰_𝑘=Δ₁+𝑦𝜃𝐼+^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1−𝜃𝑦^𝑇_𝑘−1𝑦_𝑘−1tr𝑌²_𝑘−1𝑌_𝑘−1, (19) where Δ₁ is given by (18) and 𝑌_𝑘−1=diag(𝑦²_𝑘−1,1,…,𝑦²_{𝑘−1,𝑛}) with 𝑦_{𝑘−1,𝑖} is the 𝑖th component of the vector 𝑦_𝑘−1.

For purposes of numerical illustrations, the latter diagonal formula (19) is used to initialize the L-BFGS method, although the potential use of (15) should not be neglected. Furthermore, one can observe that Δ₃ involved in the solution of the variational problem is isolated from the solution and its value does not affect the quality of the solution. This allows us to choose the value of 𝜃 freely to ensure that 𝐻⁰_𝑘 are positive definite while satisfying the weak-secant equation.

Note that maintaining positive definiteness for 𝐻⁰_𝑘 is crucial for L-BFGS method to generate descent direction. For this purpose, the following lemma suggests a possible choice on 𝜃.

Lemma 2. Assume that 𝑦_𝑘≠0 for all 𝑘. Then 𝐻⁰_𝑘 is a positive definite, if 𝑦𝜃=min1,^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1. (20)

Proof. Note that to keep positive definiteness of 𝐻⁰_𝑘, we should choose 𝜃 such that 𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1−𝜃𝑦^𝑇_𝑘−1𝑦_𝑘−1≥0, (21) which also implies 𝑦𝜃≤^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1. (22)
Therefore by our choice on 𝜃 as (20), we have the following.
Case 1. If (𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1)/𝑦^𝑇_𝑘−1𝑦_𝑘−1≥1, we will let 𝜃=1, and thus, 𝐻⁰_𝑘=Δ₁𝑦+𝐼+^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1−𝑦^𝑇_𝑘−1𝑦_𝑘−1𝑡𝑟𝑌²_𝑘−1𝑌_𝑘−1. (23)
Since 𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1≥𝑦^𝑇_𝑘−1𝑦_𝑘−1, then 𝐻⁰_𝑘 would be positive definite.
Case 2. If (𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1)/𝑦^𝑇_𝑘−1𝑦_𝑘−1<1, 𝜃₁=(𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1)/𝑦^𝑇_𝑘−1𝑦_𝑘−1 is set, which leads to 𝐻⁰_𝑘=Δ₁+𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1𝐼. (24) It is clear that for both cases, 𝐻⁰_𝑘 would maintain positive definiteness.

Hence, 𝐻⁰_𝑘 can be expressed in the following form: 𝐻⁰_𝑘=Δ₁+𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1𝑦𝐼,if^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1Δ<1,₁𝑦+𝐼+^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1−𝑦^𝑇_𝑘−1𝑦_𝑘−1tr𝑌²_𝑘−1𝑌_𝑘−1,𝑦if^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1≥1. (25) Now, we can set up the basic algorithm of our L-BFGS methods using 𝐻⁰_𝑘 as (25).

2.1. LBFGS-USD Algorithm

Step 1. Choose an initial point 𝑥₀∈ℝ^𝑛 and a positive definite matrix 𝐻⁰_𝑘=𝐼. Let 𝑘=0.

Step 2. Compute 𝑑_𝑘=−𝐻_𝑘𝑔_𝑘; then 𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑑_𝑘 where 𝛼_𝑘 satisfies the Wolfe conditions given by 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝑓𝑥_𝑘+𝛿₁𝛼_𝑘𝑔^𝑇_𝑘𝑑_𝑘,𝑔𝑥_𝑘+𝛼_𝑘𝑑_𝑘^𝑇𝑑_𝑘≥𝛿₂𝑔^𝑇_𝑘𝑑_𝑘 (26) (we always try the step length 𝛼_𝑘=1).

Step 3. Let 𝑚=min{𝑘,𝑚−1}. Update 𝐻⁰_𝑘𝑚+1 times using the pairs {𝑦_𝑗,𝑠_𝑗}^𝑘_𝑗=𝑘−𝑚; that is, let 𝐻_𝑘+1=𝑉^𝑇_𝑘⋯𝑉^𝑇_𝑘−𝑚𝐻⁰_𝑘𝑉_𝑘−𝑚⋯𝑉_𝑘+𝜌_𝑘−𝑚𝑉^𝑇_𝑘⋯𝑉^𝑇_𝑘−𝑚₊₁𝑠_𝑘−𝑚𝑠^𝑇_𝑘−𝑚𝑉_𝑘−𝑚⁺¹⋯𝑉_𝑘+𝜌_𝑘−𝑚₊₁𝑉^𝑇_𝑘⋯𝑉^𝑇_𝑘−𝑚₊₂×𝑠_𝑘−𝑚₊₁𝑠^𝑇_𝑘−𝑚⁺¹𝑉_𝑘−𝑚₊₂⋯𝑉_𝑘⋮+𝜌_𝑘𝑠_𝑘𝑠^𝑇_𝑘. (27)

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

Remark 3. The LBFGS-USD algorithm is exactly the L-BFGS algorithm of Liu and Nocedal [9], except that 𝐻⁰_𝑘 is computed by (25).

3. Convergence Analysis

We shall also establish the convergence of the LBFGS-USD algorithm. The following standard assumptions are made on the objective function.(a) The objective function 𝑓 is twice continuously differentiable.(b) The level set £=𝑥∈ℝ^𝑛∣𝑓(𝑥)<𝑓(𝑥₀) is convex, and there exist positive constants 𝑚₁ and 𝑚₂ such that, 𝑚₁‖𝑧‖²≤𝑧^𝑇𝐺(𝑥)𝑧≤𝑚₂‖𝑧‖², (28) for all 𝑧∈ℝ^𝑛 and 𝑥∈£.

Theorem 4. Let 𝑥₀ be a starting point for which 𝑓 satisfies assumptions above, and assume that the matrices 𝐵⁰_𝑘 are chosen so that {‖𝐵⁰_𝑘‖}({‖(𝐻⁰_𝑘)⁻¹‖}) are bounded. Subsequently, for any positive definite 𝐵⁰ (most often, 𝐵⁰=𝐼 is chosen), LBFGS-USD algorithm generates a sequence {𝑥_𝑘} which converges to 𝑥^∗. Moreover, there is a constant 0≤𝑟≤1 such that 𝑓𝑥_𝑘−𝑓𝑥^∗≤𝑟^𝑘𝑓𝑥₀−𝑓𝑥^∗, (29) which implies that {𝑥_𝑘} converges ℝ-linearly.

Proof. See [9].

Theorem 4 suggests that as long as 𝐻⁰_𝑘 is chosen such that ‖𝐻⁰_𝑘‖ is bounded for all 𝑘, then the corresponding L-BFGS algorithm restarts by 𝐻⁰_𝑘 would generate {𝑥_𝑘} that converges globally and ℝ-linearly. For this purpose, we give the following result which ensures that ‖𝐻⁰_𝑘‖ where 𝐻⁰_𝑘 is given by (25) is upper and lower bounded by some constants.

Lemma 5. Let 𝑥₀ be a starting point for which 𝑓 satisfies assumptions (a)-(b). Consider the sequence {𝑥_𝑘} generated by L-BFGS algorithms subject to the diagonal initial approximation, 𝐻⁰_𝑘 given by (25). If 𝑦_𝑘≠0 for all 𝑘, then ‖𝐻⁰_𝑘‖ is upper and lower bounded for all 𝑘.

Proof. Let 𝐺_𝑘−1=∫¹₀𝐺(𝑥_𝑘−1+𝑡𝑠_𝑘−1)𝑑𝑡, then 𝑦_𝑘−1=𝐺_𝑘−1𝑠_𝑘−1, and assumption (a) implies that 𝑚₁≤𝑠^𝑇_𝑘−1𝐺_𝑘−1𝑠_𝑘−1𝑠_𝑘−1²≤𝑚₂. (30) Subsequently, (28) also gives 𝑠^𝑇_𝑘−1𝑦_𝑘−1𝑦_𝑘−1²=𝑠^𝑇_𝑘−1𝐺_𝑘−1𝑠_𝑘−1𝑠^𝑇_𝑘−1𝐺²_𝑘−1𝑠_𝑘−1≤𝑚₂𝑠_𝑘−1²𝑚²₁𝑠_𝑘−1², (31) and therefore, 𝑠^𝑇_𝑘−1𝑦_𝑘−1≤𝑚₂𝑚²₁𝑦_𝑘−1². (32)
First, we begin by showing that every component of Δ₁ in (18) is bounded, so that ‖Δ₁‖ is bounded. It is more convenient to show that each part of Δ₁, namely, (i)𝜎_𝑘−1=𝑦^𝑇_𝑘−1𝑠_𝑘−1/𝑦^𝑇_𝑘−1𝑦_𝑘−1,(ii)(1/𝑦^𝑇_𝑘−1𝑠_𝑘−1+𝑦^𝑇_𝑘−1𝜎_𝑘−1𝑦_𝑘−1/(𝑦^𝑇_𝑘−1𝑠_𝑘−1)²)(𝑠^𝑇_𝑘−1𝑒_𝑖)², and(iii)−(2𝜎_𝑘−1(𝑦^𝑇_𝑘−1𝑒_𝑖)(𝑠^𝑇_𝑘−1𝑒_𝑖))/𝑦^𝑇_𝑘−1𝑠_𝑘−1, is bounded, such that Δ₁ is bounded at whole. (i)Since 1𝑚₂≤𝑦^𝑇_𝑘−1𝐺_𝑘−1⁻¹𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1≤1𝑚₁, (33) then, we obtain 1𝑚₂≤𝜎_𝑘−1≤1𝑚₁. (34) (ii)Note that 𝑠^𝑇_𝑘−1𝑒_𝑖²=𝑠²_{𝑘−1,𝑖}, (35) where 𝑠𝑘_𝑘,𝑖 is the 𝑖th component of 𝑠_𝑘−1. Then, 1𝑦^𝑇_𝑘−1𝑠_𝑘−1+𝑦^𝑇_𝑘−1𝜎_𝑘−1𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝑒_𝑖²≥𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝑒_𝑖²≥𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1𝑦_𝑘−1²𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝑒_𝑖²≥𝜎_𝑘−1𝜀²₂𝑠_𝑘−1²≥𝜎_𝑘−1𝜀²₂𝜀²₁≥1𝑚₂𝜀²₂𝜀²₁. (36) On the other hand, 1𝑦^𝑇_𝑘−1𝑠_𝑘−1+𝑦^𝑇_𝑘−1𝜎_𝑘−1𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝑒_𝑖²≤∑^𝑛_𝑖=1𝑠_{𝑘−1,𝑖}²𝑦^𝑇_𝑘−1𝑠_𝑘−1+𝑦_𝑘−1²𝑠_𝑘−1²𝑚₁𝑦^𝑇_𝑘−1𝐺_𝑘−1⁻¹𝑦_𝑘−1²≤𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝐺_𝑘−1𝑠_𝑘−1+1𝑚₁𝑚²₂𝑠_𝑘−1²𝑦_𝑘−1²≤1𝑚₁+𝑚²₂𝑚₁𝑦^𝑇_𝑘−1𝐺⁻²_𝑘−1𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1≤1𝑚₁+𝑚²₂𝑚₁1𝑚²₁=1𝑚₁+𝑚²₂𝑚³₁. (37) Therefore, it gives, in overall, 1𝑚₂𝜀²₂𝜀²₁≤1𝑦^𝑇_𝑘−1𝑠_𝑘−1+𝑦^𝑇_𝑘−1𝜎_𝑘−1𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−1²𝑠^𝑇_𝑘−1𝑒_𝑖²≤1𝑚₁+𝑚²₂𝑚³₁. (38) (iii)In a similar way, we can establish −2𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑒_𝑖𝑠^𝑇_𝑘−1𝑒_𝑖𝑦^𝑇_𝑘−1𝑠_𝑘−12≥−𝑚₁𝑦^𝑇_𝑘−1𝑠_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−12≥−𝑚₁. (39) and subsequently, −2𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑒_𝑖𝑠^𝑇_𝑘−1𝑒_𝑖𝑦^𝑇_𝑘−1𝑠_𝑘−1≤2𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑒_𝑖𝑠^𝑇_𝑘−1𝑒_𝑖𝑦^𝑇_𝑘−1𝑠_𝑘−1≤21/𝑚₁∑^𝑛_𝑖=1𝑦_{𝑘−1,𝑖}𝑠_{𝑘−1,𝑖}𝑦^𝑇_𝑘−1𝑠_𝑘−1=2𝑚₁. (40) Then, we have −2𝑚₁≤−2𝜎_𝑘−1𝑦^𝑇_𝑘−1𝑒_𝑖𝑠^𝑇_𝑘−1𝑒_𝑖𝑦^𝑇_𝑘−1𝑠_𝑘−1≤2𝑚₁. (41)
Following that, the component of Δ₁, namely Δ_1,𝑖, is bounded as follow: 1𝑚₂+1𝑚₂𝜀²₂𝜀²₁−2𝑚₁≤Δ_1,𝑖≤4𝑚₁+𝑚²₂𝑚³₁, (42) Note that, by the definition of 𝐻⁰_𝑘, we have 𝐻⁰_𝑘,𝑖=Δ_1,𝑖+𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1𝑦if^𝑇𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1𝑦^𝑇_𝑘−1𝑦_𝑘−1<1. (43) Hence, 𝐻⁰_𝑘,𝑖 satisfies 1𝑚₂+1𝑚₂𝜀²₂𝜀²₁−2𝑚₁≤Δ_1,𝑖≤𝐻⁰_𝑘,𝑖≤Δ_1,𝑖4+1≤𝑚₁+𝑚²₂𝑚³₁+1. (44) Whereas, when (𝑦^𝑇_𝑘−1𝑠_𝑘−1−𝑦^𝑇_𝑘−1Δ₁𝑦_𝑘−1)/𝑦^𝑇_𝑘−1𝑦_𝑘−1≥1, we have 1𝑚₂+1𝑚₂𝜀²₂𝜀²₁−2𝑚₁≤Δ_1,𝑖≤𝐻⁰_𝑘,𝑖≤Δ_1,𝑖𝑦+1+^𝑇𝑘−1𝑠_𝑘−1tr𝑌²_𝑘−1𝑦²_{𝑘−1,𝑀}≤4𝑚₁+𝑚²₂𝑚³₁𝑚+1+²𝑦_𝑘−1²𝑚²₁∑^𝑛_𝑖=1𝑦⁴_{𝑘−1,𝑖}𝑦²_{𝑘−1,𝑀}𝑦⁴_{𝑘−1,𝑀}≤4𝑚₁+𝑚²₂𝑚³₁+1+𝑛𝑚₂𝑦²_{𝑘−1,𝑀}𝑚²₁𝑦²_{𝑘−1,𝑀}4=1+𝑚₁+𝑚²₂𝑚³₁+𝑛𝑚₂𝑚²₁,∀𝑖, (45) where 𝑦²_{𝑘−1,𝑀} is the largest component among all 𝑦²_{𝑘−1,𝑖}. Therefore, we can conclude that ‖(𝐻⁰_𝑘,𝑖)‖ is upper and lower bounded for all 𝑘.

4. Numerical Experiences

Our test used a large set of unconstrained minimization problem consisting of 50 test problems where the list of problems is given in Table 1. These test problems are selected from Moré et al. [2], CUTE [3], Toint [4], and various other test function collections such as in [10]. The subroutine of test problems is available at http://camo.ici.ro/forum/SCALCG/evalfg.for (accessed on Jan 2012). The method tested is as follows:(1)LBFGS-I: L-BFGS method with the initial matrix, 𝐻⁰_𝑘=𝐼,(2)LBFGS-𝑚I: L-BFGS method with the initial matrix, 𝐻⁰_𝑘=𝜃𝐼 where 𝜃_𝑘 is the Oren-Luenberger scaling at 𝑘,(3)LBFGS-USD: L-BFGS method with the initial matrix, 𝐻⁰_𝑘 is given by (19).

Table 1 

Selected test problems.

5. Discussion

In general, Figure 1 indicates that the new diagonal initial approximating matrix are substantially better, followed by both standard initializations of the L-BFGS method in terms of number of iterations, function/gradient calls, and CPU time, respectively. To better study the effect of our initial approximation, we include Tables 2, 3, and 4 that give ratio of function/gradient calls over iteration counts for the L-BFGS method with standard initial matrices and our diagonal initial approximating matrix. As conclusion, our diagonal initial approximation, LBFGS-USD, performs better in which the ratio is close to one and hence would likely to accept unit step length compared to the LBFGS-I (without any preconditioning) method.

Table 2 

Ratio of function/gradient calls over iteration counts for 𝑚=3.

Table 3 

Ratio of function/gradient calls over iteration counts for 𝑚=5.

Table 4 

Ratio of function/gradient calls over iteration counts for 𝑚=7.

(a)

(b)

(c)

Figure 1 

Performance profile of LBFGS-USD based on iterations, function/gradient call, and CPU time for all choices of 𝑚 and 𝑛.

Moreover for LBFGS-USD method, it requires, in general less iteration counts as much as 38% and 19% than LBFGS-I and LBFGS-𝑚I, respectively. Meanwhile, in terms of number of function/gradient, LBFGS-MTD needs 46% and 16% less function/gradient counts, respectively. Finally, LBFGS-MTD requires 33% and 6% less CPU time in second, respectively. In conclusion, the numerical results for a broad class of the test problems show that the LBFGS-USD algorithm is efficient and vast superior in solving small to large size problems.

6. Concluding Remarks

We proposed technique that exploits the presence of the Hessian in the diagonal matrix form. Under some standard assumption on the objective function, we observe that the convergence of the diagonal initial approximation of the L-BFGS scheme is ℝ-linear. Based on our numerical results, we believe that the following conclusions can be made on the diagonal approximation that is derived in this study. (i)Our diagonal approximating matrix is able to maintain positive definiteness in a very simple way and give 𝑂(𝑛) storage requirement. (ii)The numerical experiments show that our diagonal initial approximating matrix is generally effective for the L-BFGS method compared to the standard initial matrices in the literature.

Conflict of Interests

The author(s) declare(s) that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the Malaysian MOHE-FRGS Grant no. 01-11-09-722FR.