The Modified HZ Conjugate Gradient Algorithm for Large-Scale Nonsmooth Optimization.

In this paper, the Hager and Zhang (HZ) conjugate gradient (CG) method and the modified HZ (MHZ) CG method are presented for large-scale nonsmooth convex minimization. Under some mild conditions, convergent results of the proposed methods are established. Numerical results show that the presented methods can be better efficiency for large-scale nonsmooth problems, and several problems are tested (with the maximum dimensions to 100,000 variables).

Introduction

Consider the following optimization problem: If the constrained set satisfies Φ = ℜ, then Eq (1) is called unconstrained optimization; if Φ = {x ∣ l ≤ x ≤ u, x ∈ ℜ}, where n is the number of variables and the vectors l and u represent the lower and upper bounds on the variables, then Eq (1) is called box-constrained optimization; and if Φ = {x ∣ h(x) = 0, g(x) ≤ 0, x ∈ ℜ, i = 1, ⋯, r, j = 1, ⋯, k}, where r and k are positive integers, then Eq (1) is called normal constrained optimization. If the objective function f: ℜ → ℜ is continuously differentiable, then Eq (1) is called smooth optimization; if f: ℜ → ℜ is a nondifferentiable function, then Eq (1) is called nonsmooth optimization. For a given objective function f and constrained set Φ, Eq (1) will be called the corresponding optimization problem.

Optimization problems are encountered in many fields, including engineering, management, finance, medicine, and biology, and similarly, optimization models can be used in many fields (see [1-13]). At present, there are many efficient methods available for solving optimization problems [14-24]. However, many challenging optimization problems exist, for example, large-scale problems and nonsmooth problems. The workload increases greatly as the dimension of the problem increases, causing the required CPU time to become very long. For certain large-scale problems, a computer may fail to solve them. To address this issue, more efficient algorithms should be designed. It is well known that spectral gradient approaches, conjugate gradient (CG) techniques and limited-memory quasi-Newton methods can cope with large-scale optimization problems. CG methods in particular have been widely used in many practical large-scale optimization problems because of their simplicity and low memory requirements.

Nonsmooth problems are believed to be very difficult to solve, even when they are unconstrained. The direct application of smooth gradient-based methods to nonsmooth problems may lead to a failure in optimality conditions, convergence, or gradient approximation [25]. Haarala et al. (see, e.g., [26, 27]) presented limited-memory bundle methods for large-scale nonsmooth unconstrained and constrained minimization and demonstrated their application to test problems of up to one thousand variables in dimension. Karmitsa et al. [28] tested and compared various methods of both types as well as several methods that may be regarded as hybrids of these two approaches and/or others; the dimensions of the tested nonsmooth problems ranged from 20 to 4000, and the most effective method for large- and extra-large-scale problems was found to be that of [27]. Therefore, special tools for solving nonsmooth optimization problems are needed.

This paper is organized as follows. In the next section, the Hager and Zhang (HZ) CG method is presented, and a modified HZ (MHZ) CG formula is proposed. In Section 3, the application of the HZ and MHZ methods to large-scale nonsmooth problems is discussed, global convergence is established, and numerical experiments on nonsmooth problems are reported. In the last section, conclusions are presented. Throughout the paper, ‖⋅‖ denotes the Euclidean norm.

The HZ CG formula [29, 30] and a modification thereof

For convenience, we rewrite Eq (1) as the following special case: where f: ℜ → ℜ is continuously differentiable; Eq (2) is called an unconstrained optimization problem. CG methods are a class of effective line search methods for solving Eq (2), especially when the dimension n is large. The iterative formula for CG methods is defined as follows: where x is the current point in the iteration, α > 0 is the step length, and d is the search direction; the latter is determined as where g = ∇f(x) is the gradient of f(x) at point x and β ∈ ℜ is a scalar. The parameter β is chosen such that when the method is applied to minimize a strongly convex quadratic function, the directions d and d are conjugate with respective to the Hessian of the quadratic function. Let where from [29] with y = g − g; we call this formula for [29] the HZ formula. If , then Eq (4) with satisfies This method exhibits global convergence for the strongly convex function f. To obtain a similar result for a general nonlinear function, Hager and Zhang [30] proposed the formula where η > 0 is a constant; in their experiments, they set η = 0.01. This new parameter also has the property expressed in Eq (6).

Based on the formula for , if we let , we can obtain the more general formula In this paper, we set where c ∈ (0, 1) is a scalar. It is easy to deduce that ; moreover, if , then Eq (7) is identical to the HZ formula. The modified formula has the following features:

(i) The new formula can overcome the shortcomings of the CG parameter β. If , it is not difficult to find that

(ii) Another property of this formula is that it can ensure that the new direction d in Eq (4) with belongs to a trust region without the need for any line search technique. By Eq (4), for , we obtain By combining this result with Step 1 of Algorithm 2.2, we find that for all k.

(iii) If T ≠ 0, then the new direction d in Eq (4) when possesses the following sufficient descent property: which holds for all k. Now let us analyze this result. If k = 0, we have d1 = −g1 and satisfying Eq (10). For k ≥ 1, because T ≠ 0, multiplying Eq (4) by yields Let and ; then, by the inequality , we have

Nonsmooth Convex Problems and Their Results

Consider the unconstrained convex optimization problem where f: ℜ → ℜ is a possibly nonsmooth convex function. For the special case that f is continuously differentiable, this optimization problem has been well studied for several decades. However, nonsmooth optimization problems of the form of Eq (12) also arise in many applications, such as image restoration [31] and optimal control [32]. The Moreau-Yosida regularization of f generates where λ is a positive parameter; then, Eq (12) is equivalent to the following problem: Let p(x) = argmin θ(z) and define a function Thus, p(x) is well defined and unique because the function θ(z) is strongly convex. Therefore, F(x) can be expressed as follows: F(x) possesses many known features (see, e.g., [33-35]). The generalized Jacobian of F(x) and its property of BD-regularity are demonstrated in [36, 37]. Here, we list some additional findings regarding the function F(x), as follows:

(i) x is an optimal solution to Eq (12) ⇔ ∇F(x) = 0, i.e., p(x) = x.

(ii) where is the gradient of F.

(iii) The set of generalized Jacobian matrices ∂g(x) = {V ∈ ℜ: V = lim∇g(x), x ∈ D} is nonempty and compact, where

Furthermore, every V ∈ ∂g(x) is a symmetric positive semidefinite matrix for each x ∈ ℜ because g is a gradient mapping of the convex function F.

(iv) There exist two constants, μ1 > 0 and μ2 > 0, and a neighborhood Ω of x that satisfy

Two algorithms for nonsmooth problems

Based on the above discussion, we present two algorithms for application to nonsmooth problems of the form Eq (14); afterward, we analyze the solution of Eq (12). In the following, unless otherwise noted, g = g(x) is defined as in Eq (16).

Algorithm 1 Algorithm 4.1.

Require:

An initial point x0 ∈ ℜ, λ > 0, σ, η ∈ (0, 1), ρ ∈ (0, 1/2], ϵ ∈ [0, 1).

Specify g0 by solving the subproblem Eq (13);

d0 ← −g0, k ← 0.

while ‖g‖ > ϵ do

Determine the step size α = max{ρ|j = 0, 1, 2, ⋯} satisfying the following Armijo line search condition

x = x + α d;

Compute g by solving the subproblem Eq (13);

if ‖g‖ ≤ ϵ then

break.

else

Compute d as

end if

x ← x, d ← d, k ← k + 1.

end while

Algorithm 4.2. Eq (18) of Algorithm 4.1 is replaced with the following: Compute d as and let k: = k + 1, go back to Step 2.

The following assumptions are needed to ensure the global convergence of Algorithms 4.1 and 4.2.

Assumption 4.1. (i) F is bounded from below and the sequence {V} is bounded, namely, there exists a constant M > 0 such that

(ii) g is BD-regular at x, namely, item (iv) in Section 3 above holds.

By Assumption 4.1, it is not difficult to deduce that there exists a constant M* > 1 such that

Lemma 1 The sequence {x} is generated by Algorithm 4.1 (or Algorithm 4.2). Let Assumption 4.1 hold; then, for sufficiently large k, there exists a constant α* > 0 that satisfies

Proof. Suppose that α satisfies the Armijo line search condition Eq (17). The proof is complete if α = 1 holds. Otherwise, let ; then, we have By performing a Taylor expansion, we obtain where and the last inequality follows from Eq (20). By combining Eqs (21) and (23), we obtain Thus, we find that Let , and the proof is complete.

Now, let us prove the global convergence of Algorithm 4.1.

Theorem 1 Suppose that the conditions in Lemma 1 hold. Then, we have moreover, any accumulation point of x is an optimal solution of Eq (12).

Proof. Suppose that Eq (25) is not true. Then, there must exist two constants, ϵ0 > 0 and k* > 0, that satisfy By combining , Eqs (17), (22) and (26), we obtain Because F(x) is bounded from below for all k, it follows from Eq (27) that This contradicts Therefore, Eq (25) holds. Let x* be an accumulation point of {x} and, without loss of generality, suppose that there exists a subsequence {x} that satisfies From Eq (28), we find that ‖g(x*)‖ = ‖∇F(x*)‖ = 0. Then, from property (i) of F(x) as given in Section 3, x* is an optimal solution of Eq (12). The proof is complete.

In a manner similar to Theorem 4.1 in [38], we can establish the linear convergence rate of Algorithm 4.1 (or Algorithm 4.2). Here, we simply state this property, as follows, but omit the proof.

Theorem 2 Let Assumptions 4.1 (i) and (ii) hold, and let x* be the unique solution of Eq (14). Then, there exist two constants b > 0 and r ∈ (0, 1) that satisfy namely, the sequence {x} generated by Algorithm 4.1 (or Algorithm 4.2) linearly converges to x*.

Numerical results for nonsmooth problems

In this section, we present several numerical experiments using Algorithms 4.1 and 4.2 and a modified Polak-Ribière-Polyak conjugate gradient algorithm (called MPRP) [18]. It is well known that the CG method is very effective for large-scale smooth problems. We will show that these two algorithms are also applicable to large-scale nonsmooth problems. The nonsmooth academic test problems that are listed, along with their initial points, in Table 1 are described in [27], where “Problem” is the name of the test problem, “x0” is the initial point, and the corresponding numbers of variables are also given, “f” is the optimal value of the test problem. Problems 1-5 are convex functions, and the others are nonconvex functions. The detailed characteristics of these problems can be found in [27]. Because we wished to test the three considered methods for application to large-scale nonsmooth problems, problem dimensions of 5000, 10000, 50000, and 100000 were chosen. In our experiments, we found that problem 2 required a considerably amount of time to solve; therefore, we set its largest dimension to 50000.

Table 1

Test problems and their initial points and optimal value.

No.	Problem	x0	fops
1	Generalization of MAXQ	(1, 2, ⋯, n/2, −(n/2 + 1), ⋯, −n)	0
2	Generalization of MXHILB	(1, 1, ⋯, 1)	0
3	Chained LQ	(−0.5, −0.5, ⋯, −0.5)	-2(n-1)
4	Chained CB3 I	(2, 2, ⋯, 2)	2(n − 1)
5	Chained CB3 II	(2, 2, ⋯, 2)	2(n − 1)
6	Number of active faces	(1, 0, ⋯, 1, 0)	0
7	Nonsmooth generalization of Brown function 2	(−1, −1, ⋯, −1)	0
8	Chained Mifflin 2	(−1.5, 2, ⋯, −1.5, 2)	varies
9	Chained Crescent I	(1, 0, ⋯, 1, 0)	0
10	Chained Crescent II	(1, 0, ⋯, 1, 0)	0

Both algorithms were implemented using Fortran PowerStation 4.0 with double-precision arithmetic, and all experiments were run on a PC with a Core 2 Duo E7500 CPU @2.93 GHz with 2.00 GB of memory and the Windows XP operating system. The following parameters were chosen for Algorithms 4.1 and 4.2 and MPRP: σ = 0.8, ρ = 0.5, c = 0.01, ϵ = 1E − 15, and η = 0.01. We stopped the algorithms when the condition ‖g(x)‖ ≤ 1E − 5 or ∣F(x) − F(x)∣ ≤ 1E − 8 or |f(x) − f| ≤ 1E − 4 was satisfied. If the number of iterations exceeded ten thousand, the program would also terminate. Because a line search cannot always ensure the descent condition , an uphill search direction may arise in real numerical experiments, which may cause the line search rule to fail. To avoid this condition, the step size α was accepted only if the search number in the line search was greater than five.

In the experiments, the subproblem Eq (13) was solved using the PRP CG method (called the sub-algorithm), and its numbers of iterations and function evaluations were added to those of Algorithm 4.1, Algorithm 4.2 or MPRP (called the main algorithm). In the sub-algorithm, if ‖∂f(x)‖ ≤ 1E − 4 or f(x) − f(x) + ‖∂f(x)‖2 − ‖∂f(x)‖2 ≤ 1E − 3 (see [39]) holds, where ∂f(x) is the subgradient of f(x) at point x, then the algorithm terminates. The sub-algorithm will also terminate when the iteration number exceeds ten. For the line search, the Armijo line search technique was used and the step length was accepted if the search number was greater than five. The columns in Tables 2, 3 and 4 have the following meanings:

Table 2

Test results for Algorithm 4.1.

No.	Dim	NI/NF	‖g(x)‖	Time	ffinal
1	5000	122 / 2304	4.963444E-05	3.109375E-01	2.977173E-08
1	50000	142 / 2720	5.295883E-04	3.451688E+00	3.177435E-08
1	100000	148 / 2844	8.533346E-04	8.389563E+00	2.559965E-08
2	5000	48 / 902	0	5.279200E+02	0
2	10000	52 / 978	0	2.290124E+03	0
2	50000	63 / 1194	3.689089E-10	1.577677E+04	9.789607E-07
3	3000	11 / 55	6.505226E-16	7.031086E-04	-4.241179E+03
4	3000	5 / 50	0	2.943751E-02	5.998031E+03
5	3000	14 / 65	0	4.543751E-02	5.998000E+03
6	5000	64 / 1240	7.081679E-16	5.470469E-01	1.881710E-06
6	10000	69 / 1345	7.428122E-16	1.187422E+00	2.725454E-06
6	50000	82 / 1618	6.935890E-16	7.172547E+00	5.888909E-06
6	100000	87 / 1723	7.275201E-16	1.548309E+01	8.529442E-06
7	5000	12 / 59	2.710511E-16	2.024688E-01	2.327843E-06
7	10000	13 / 61	5.421024E-16	4.373125E-01	4.656153E-06
7	50000	23 / 82	6.776293E-16	2.733234E+00	1.164132E-05
7	100000	41 / 119	3.388159E-16	8.077172E+00	1.164146E-05
8	3000	12 / 58	6.505227E-16	1.404690E-02	-2.120705E+03
9	5000	16 / 124	5.750809E-16	1.084375E-01	7.202892E-07
9	10000	16 / 125	2.875404E-16	3.114688E-01	7.198567E-07
9	50000	26 / 146	3.594263E-16	3.202235E+00	1.798780E-06
9	100000	27 / 148	7.188527E-16	9.281172E+00	3.597344E-06
10	5000	13 / 61	8.470350E-16	3.274688E-01	4.365231E-06
10	10000	14 / 64	4.235176E-16	5.623125E-01	4.365406E-06
10	50000	42 / 121	5.293999E-16	3.483234E+00	1.091389E-05
10	100000	51 / 140	2.647004E-16	9.937172E+00	1.091395E-05

Table 3

Test results for Algorithm 4.2.

No.	Dim	NI/NF	‖g(x)‖	Time	ffinal
1	5000	186 / 3878	4.403006E-06	4.539688E-01	2.641011E-09
1	50000	242 / 5051	8.447711E-05	6.235391E+00	5.068475E-09
1	100000	259 / 5377	3.400454E-04	1.389234E+01	1.020121E-08
2	5000	98 / 2010	5.769694E-06	1.189811E+03	3.089375E-04
2	10000	107 / 2199	3.611318E-06	5.205405E+03	1.859985E-04
2	50000	129 / 2679	3.710004E-06	1.349417E+04	9.817318E-05
3	3000	11 / 55	6.505226E-16	1.590626E-02	-4.241179E+03
4	3000	7 / 89	0	2.552654E-02	5.998031E+03
5	3000	16 / 104	0	7.656254E-02	5.998000E+03
6	5000	64 / 1240	7.081679E-16	5.306719E-01	1.881710E-06
6	10000	69 / 1345	7.428122E-16	1.157047E+00	2.725454E-06
6	50000	82 / 1618	6.935890E-16	7.391297E+00	5.888909E-06
6	100000	87 / 1723	7.275201E-16	1.539180E+01	8.529442E-06
7	5000	12 / 42	4.656623E-16	2.339532E-01	2.327843E-06
7	10000	13 / 44	9.313248E-16	4.377500E-01	4.656153E-06
7	50000	23 / 66	2.910396E-16	2.734688E+00	1.164132E-05
7	100000	41 / 102	5.820813E-16	7.891281E+00	1.164146E-05
8	3000	12 / 58	6.505227E-16	3.190626E-02	-2.120705E+03
9	5000	16 / 107	9.879814E-16	2.969532E-01	7.202892E-07
9	10000	16 / 108	4.939907E-16	5.157501E-01	7.198567E-07
9	50000	26 / 129	6.174896E-16	3.202688E+00	1.798780E-06
9	100000	27 / 132	3.087449E-16	9.000281E+00	3.597344E-06
10	5000	13 / 45	3.637988E-16	3.119532E-01	4.365231E-06
10	10000	14 / 47	7.275977E-16	5.627500E-01	4.365406E-06
10	50000	42 / 104	9.095022E-16	3.468688E+00	1.091389E-05
10	100000	51 / 123	4.547519E-16	9.625281E+00	1.091395E-05

Table 4

Test results for MPRP.

No.	Dim	NI/NF	‖g(x)‖	Time	ffinal
1	5000	250 / 5197	1.146977E-04	6.209688E-01	6.879798E-08
1	50000	286 / 5991	1.099941E-03	8.980437E+00	6.599447E-08
1	100000	297 / 6222	2.127262E-03	2.067906E+01	6.381689E-08
2	5000	98 / 2025	5.373446E-15	1.105497E+03	9.428024E-09
2	10000	107 / 2214	3.363302E-15	4.828827E+03	5.676222E-09
2	50000	129 / 2693	1.382084E-14	1.683294E+04	5.992015E-09
3	3000	11 / 55	6.505226E-16	3.221878E-02	2.793039E-06
4	3000	7 / 89	0	4.721878E-02	7.714992E+03
5	3000	16 / 104	0	7.821879E-02	7.714994E+03
6	5000	64 / 1240	7.081679E-16	6.085625E-01	1.881710E-06
6	10000	69 / 1345	7.428122E-16	1.265188E+00	2.725454E-06
6	50000	82 / 1618	6.935890E-16	7.234313E+00	5.888909E-06
6	100000	87 / 1723	7.275201E-16	1.573519E+01	8.529442E-06
7	5000	12 / 59	2.710511E-16	1.105805E+03	2.327843E-06
7	10000	13 / 61	5.421024E-16	4.829446E+03	4.656153E-06
7	50000	23 / 82	6.776293E-16	1.037413E+01	1.164132E-05
7	100000	41 / 119	3.388159E-16	1.103525E+01	1.164146E-05
8	3000	12 / 58	6.505227E-16	3.259378E-02	-2.120705E+03
9	5000	16 / 124	5.750809E-16	1.105867E+03	7.202892E-07
9	10000	16 / 125	2.875404E-16	4.829565E+03	7.198567E-07
9	50000	26 / 146	3.594263E-16	1.115213E+01	1.798780E-06
9	100000	27 / 148	7.188527E-16	1.286225E+01	3.597344E-06
10	5000	13 / 61	8.470350E-16	1.105892E+03	4.365231E-06
10	10000	14 / 64	4.235176E-16	4.829612E+03	4.365406E-06
10	50000	42 / 121	5.293999E-16	1.158712E+01	1.091389E-05
10	100000	51 / 140	2.647004E-16	1.381325E+01	1.091395E-05

Dim: the dimensions of problem.

NI: the total number of iterations.

NF: the number of function evaluations.

‖g(x)‖: the norm of g(x) at the final iteration.

Time: the CPU time in seconds.

f: the value of f(x) at the final iteration.

From the above three tables, it is not difficult to see that Algorithm 4.1 is superior to Algorithm 4.2 and MPRP in terms of NI, NF, and ‖g(x)‖. Algorithm 4.1 yields smaller values of NI and NF when the program terminates; moreover, the value of ‖g(x)‖ is smaller than that for Algorithm 4.2 in most cases. However, we also note that our algorithms can efficient solve the 3000 dimensional(with maximum dimensional) case for problems 3, 4, 5 and 8, if we increase the dimensional, these algorithms fail to converge to good minima and become stuck at local. To directly illustrate the performances of these two methods, we used the tool developed by Dolan and Moré [40] to analyze their efficiencies in terms of the number of iterations, number of function evaluations, and CPU time. In the following, Figs 1, 2 and 3 represent the results presented in Tables 2, 3 and 4 in terms of NI, NF, and Time, respectively.

Fig 1

Performance profiles of these methods (NI).

Fig 2

Performance profiles of these methods (NF).

Fig 3

Performance profiles of these methods (Time).

From Figs 1 and 2, we can conclude that Algorithm 4.1 performs better than Algorithm 4.2 and MPRP do in terms of the numbers of iterations and function evaluations. Moreover, the excellence of Algorithm 4.1 can obviously be attributed to the superior theoretical properties of the MHZ method compared with the usual HZ method. However, Fig 3 indicates that Algorithm 4.2 is superior to Algorithm 4.1 and MPRP in terms of CPU time. Overall, all methods are very effective for application to large-scale nonsmooth optimization problems.

Conclusion

(i) In this paper, we focus on the HZ CG method and study the application of this method to solve nonsmooth optimization problems. Several results are presented that prove the efficiency of this method for application to large-scale problems of nonsmooth unconstrained optimization.

(ii) Motivated by the HZ formula, we also present a modified HZ CG formula. The modified HZ formula not only possesses the sufficient descent property of the HZ formula but also belongs to a trust region and has the non-negative scale

(iii) We report the results of applying three methods to solve large-scale nonsmooth convex minimization problems. Global convergence is achieved, and numerical experiments verify that both methods can be successfully used to solve large-scale nonsmooth problems.

(iv) Although the HZ and MHZ methods offer several key achievements for large-scale nonsmooth optimization, we believe that there are at least five issues that could be addressed to gain further improvements. The first is the scale c in the modified HZ CG algorithm, which could be adjusted. The second is the application of other CG methods for this type of optimization areas; perhaps there exists a more suitable CG method for this purpose. Regarding the third issue, it is well known that limited-memory quasi-Newton methods are effective techniques for solving certain classes of large-scale optimization problems because they require minimal storage; this inspires us to combine limited-memory quasi-Newton methods with the HZ CG technique to solve large-scale nonsmooth optimization. In future, we will also use the HZ CG method to investigate large-scale nonsmooth optimization with constraints; this is the fourth issue that we believe must be addressed. The last issue is the most important one, namely, the consideration of other optimality conditions and convergence conditions in nonsmooth problems should be paid. All of these issues will be addressed in our future work.