- Email: [email protected]
An approximate decomposition algorithm for convex minimization
Journal of Computational and Applied Mathematics 234 (2010) 658–666
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An approximate decomposition algorithm for convex minimizationI Yuan Lu a , Li-Ping Pang b,∗ , Xi-Jun Liang b , Zun-Quan Xia b a
School of Sciences, Shenyang University, Shenyang 110044, China
b
School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China
article
info
Article history: Received 11 September 2009 Received in revised form 3 January 2010 MSC: 65Kxx Keywords: Nonsmooth convex optimization VU-decomposition Approximate U-Lagrangian Proximal bundle method Smooth path
abstract For nonsmooth convex optimization, Robert Mifflin and Claudia Sagastizábal introduce a VU-space decomposition algorithm in Mifflin and Sagastizábal (2005) [11]. An attractive property of this algorithm is that if a primal–dual track exists, this algorithm uses a bundle subroutine. With the inclusion of a simple line search, it is proved to be globally and superlinearly convergent. However, a drawback is that it needs the exact subgradients of the objective function, which is expensive to compute. In this paper an approximate decomposition algorithm based on proximal bundle-type method is introduced that is capable to deal with approximate subgradients. It is shown that the sequence of iterates generated by the resulting algorithm converges to the optimal solutions of the problem. Numerical tests emphasize the theoretical findings. © 2010 Elsevier B.V. All rights reserved.
1. Introduction In the area of continuous optimization many real-life problems can be written as min f (x), x∈Rn
(1)
with a convex and not necessarily differentiable function f : Rn → R. This fact gives reason for studying and developing algorithms to solve it numerically. If f is not differentiable at a solution x¯ of (1), constructing fast algorithm to compute x¯ is a challenge. Essentially, this has to do with the intrinsic difficulty in using appropriate second-order objects. One line of research [1] suggests introducing second-order information about f by means of its Moreau–Yosida regularization, which is a smooth function. Implementable forms of this method can be obtained by using a bundle technique, alternating serious steps with a sequence of null steps [2]. The other line of research parts from the viewpoint that nonsmoothness in practical applications is usually structured; see [3,4]. Nonsmooth functions may behave smoothly and even have appropriate secondorder representations on certain manifolds along certain directions. For the second line of research, new conceptual schemes have been developed, which are based on the VU-theory introduced in [5]; see also [6,7]. The idea is to decompose Rn into two orthogonal subspace V and U at a point x¯ that the nonsmoothness of f is concentrated essentially on V , and the smoothness of f appears on U-subspace. More precisely, for a given g¯ ∈ ∂ f (¯x), where ∂ f (¯x) denotes the subdifferential of f at x¯ in the sense of convex analysis. Then Rn can be decomposed as direct sum of two orthogonal subspaces, i.e., Rn = U ⊕ V , where V = lin(∂ f (¯x) − g¯ ), and U = V ⊥ . They define the U-Lagrangian, an approximation of the original function, which can be used to create a second-order expansion of f along certain manifolds. The manifolds can be identified and called primal track if f has a strongly transversal primal–dual
I The research is supported by the National Natural Science Foundation of China under project No. 10771026.
∗
Corresponding author. E-mail address: [email protected] (L.-P. Pang).
0377-0427/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2010.01.003
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
659
gradient structure; see [8]. The conceptual VU-algorithm in [5] finds points on primal track, by generating minimizing steps in the V -subspace. Alternating with these corrector steps are U-Newton predictor steps in order to provide for superlinear convergence. Accordingly, there is a dual track which provides the U-gradients needed for the U-Newton steps. Therefore, the study of approximation of primal–dual track is important. A first important step for this approximation is given in [9], where it is shown that, for a point near a minimizer, its corresponding proximal point lands on the primal track. This result is fundamental for implementation of the minimization algorithm, because a sequence of null steps from a bundle mechanism can approximate proximal points with any desired accuracy [10]. Bundle mechanism also can approximate a dual track by finding the minimum norm vector in the convex hull of the subgradients collected for estimating the proximal point. Based on the analysis above, Robert Mifflin and Claudia Sagastizábal design a VU-algorithm in [11]. This algorithm brings the iterate to the primal track with the help of bundle subroutine. Then the U-Newton step is performed to gain superlinear decrease of the distance to solution. However, this algorithm is conceptual in the sense that it makes use of the exact subgradients of the objective function, which is not easily to compute. The purpose of this paper is to explore the possibility of utilizing approximate subgradients instead of the exact ones. We present an approximate decomposition algorithm frame by defining the approximate U-Lagrangian. By introducing the definitions of smooth path and dual path, we obtain an implemental approximate decomposition algorithm based on proximal bundle subroutine [12], which can deal with approximate subgradients. This modification makes the algorithm easier to implement. The paper is organized as follows. In the next section we present the approximate U-Lagrangian. Based on this definition, we design a conceptual Algorithm I. It can deal with the approximate subgradients and only need to compute the ε -optimal solutions in approximate V -step. In the first part of Section 3 we propose a smooth path. The second part of Section 3 is devoted to establish an implemental Algorithm II which uses approximate subgradients and substitutes the approximate V -step in Algorithm I with proximal bundle subroutine. Numerical testing of resulting Algorithm II is reported in the final section. 2. A conceptual approximate decomposition algorithm 2.1. Approximate U-Lagrangian and its properties Since the U-Lagrangian itself is defined through a minimization problem involving f , exact evaluation of it in general practically is impossible. In this section we will briefly outline the definition of an approximation of U-Lagrangian and some related properties as well. For more details, one can refer to [13]. To begin with, we recall the concepts of ε -minimizer, ε -minimum and ε - min ϕ , where ϕ : Rp → R1 . We say that xε-min is an ε -minimizer of ϕ if xε-min ∈ Arg(ε - min ϕ)
= {x ∈ Rp | ϕ(z ) ≥ ϕ(x) − ε, ∀z ∈ Rp }, and that ϕε-min is an ε -minimum of ϕ if there exists an ε -minimizer xε-min such that
(2)
ϕ(z ) ≥ ϕε-min − ε, ∀z ∈ Rp , (3) where ϕε-min = ϕ(xε-min ), and that ε - min ϕ is ε -minimization of ϕ if minimizing ϕ is in the sense of (3). Under these concepts the definition of an approximation U-Lagrangian of f associated with some ε can be given below. Definition 2.1. The approximate U-Lagrangian of f with respect to ε (and g¯V ), denoted by Lg¯ ,ε (u), is defined as follows, Lg¯ ,ε (u) = ε - min{f (¯x + u ⊕ v) − h¯gV , viV },
(4)
v∈V
and Wg¯ ,ε (u) is defined by
Wg¯ ,ε (u) = Arg ε - min{f (¯x + u ⊕ v) − h¯gV , viV } . v∈V
Let ε > 0 be given. Introducing ε to [5, Prop. 2.4], one can restate it as the following lemma. Lemma 2.1. Suppose g¯ ∈ ∂ f (¯x)
T
ri∂ε f (¯x). Then Wg¯ ,ε (u) is nonempty.
Theorem 2.1. Suppose g¯ ∈ ∂ f (¯x) and g¯ ∈ ∂ f (¯x) ri∂ε f (¯x). Then the following assertions are true: (i) the function Lg¯ ,ε defined in (4) is convex and finite everywhere; (ii) Wg¯ ,ε (u) = {wε,u | ∃g ε ∈ ∂ε f (¯x + u ⊕ wε,u ) : gVε = g¯V }; (iii) 0 ∈ Wg¯ ,ε (0) and Lg¯ ,ε (0) = f (¯x).
T
Theorem 2.2. Suppose Wg¯ ,ε (u) 6= ∅. Then one has that
∂ Lg¯ ,ε (u) = {u∗ε | u∗ε ⊕ g¯V ∈ ∂ε f (¯x + u ⊕ wε,u )}, ∂ε Lg¯ ,ε (u) = {u∗ε | u∗ε ⊕ g¯V ∈ ∂2ε f (¯x + u ⊕ wε,u )}.
(5)
660
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Corollary 2.1. If 0 ∈ Wg¯ ,ε (0), then
∂ Lg¯ ,ε (0) = (∂ε f (¯x) − g¯ ) ∩ U. Theorem 2.3. Suppose g¯ ∈ ri∂ f (¯x). Then the following conclusion holds,
kwε,u k ≤ o(kuk) + ε,
∀wε,u ∈ Wg¯ ,ε (u).
2.2. Approximate decomposition algorithm frame It is well known that Hessian matrix plays a very important role in constructing numerical methods. In order to give an approximate decomposition algorithm frame, we restate the definition of Hessian matrix in [5] as follows. Definition 2.2. Assume that f is finite and x¯ and g¯ ∈ ∂ f (¯x) are fixed. We say that f has at x¯ a U-Hessian HU f (¯x) associated with g¯ if Lg¯ (u) has a generalized Hessian at 0, setting HU f (¯x) = HLg¯ (0). Based on the definition of approximate U-Lagrangian, we investigate an approximate decomposition algorithm frame. Algorithm I. Begin of Algorithm Step 0 Initiation x0 , η > 0, 0 < τ < 1, ε0 > 0, g 0 ∈ ∂ε0 f (x0 ), k = 0. Step 1 Stop if kg k k ≤ η. Step 2 Approximate V -step. Compute an ε -optimal solution δv ∈ V satisfying
δv ∈ Arg{ε - min f (xk + 0 ⊕ δv)}. δv∈V
Set x0 = xk + 0 ⊕ δv . Step 3 U-step. Make a Newton step in x0 . Compute g k ∈ ∂εk f (x0 ): gVk = 0,
k gU ∈ ∂εk Lg k ,εk ((x0 − x¯ )U ).
Compute the solution δ u ∈ U satisfying k gU + HU f (¯x)δ u = 0.
Set xk+1 = x0 + δ u ⊕ 0 = xk + δ u ⊕ δv and εk+1 = τ εk . Step 4 Update-set. Compute g k+1 ∈ ∂εk+1 f (xk+1 ). Set k = k + 1 and go to Step 1. End of algorithm Theorem 2.4. Assume that g¯ = 0 ∈ ∂ f (¯x) and f has a positive definite U-Hessian at x¯ , a minimizer of f . Then the iterate points {xk } constructed by Algorithm I satisfies
kxk+1 − x¯ k = o(kxk − x¯ k). Remark 2.1. If {εk }∞ k=1 ≡ 0, this algorithm is the same as Algorithm 4.5 in [5]. However, it uses approximate subgradients of the objective function and only needs to compute ε -optimal solutions in V -step, which make the algorithm easier to implement. 3. Approximate decomposition algorithm Since Algorithm I in Section 2 relies on knowing the subspaces U and V and converges only locally, it needs significant modification. In [9], Robert Mifflin and Claudia Sagastizábal show that a proximal point sequence follows primal track near a minimizer. This opens the way for defining a VU-algorithm where V -steps are replaced by proximal steps. In addition, the proximal step can be estimated with a bundle technique which also can approximate the unknown U and V subspaces as a computational by-product. Therefore, they establish Algorithm 6 in [11] by combining the bundle subroutine with the VU-space decomposition method. However, this algorithm needs the exact subgradients of the objective function, which is expensive to compute. Therefore, the study of using approximate subgradients instead of the exact ones is deserve.
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
661
3.1. Smooth path to ε -minimizer Given a positive scalar parameter µ, the ε -proximal point function depending on f , is defined by
1
pµ,ε (x) := arg ε − min f (p) +
2
p∈Rn
µkp − xk2
,
x ∈ Rn ,
where k · k stands for the Euclidean norm. It has the property: gµ,ε (x) := µ(x − pµ,ε (x)) ∈ ∂ε f (pµ,ε (x)). Similarly to the definition of primal track, we give the definition of smooth path. Definition 3.1. For any ε > 0, µ = µ(x) > 0, we say that Θε (u) = x¯ + u ⊕ vε (u) is a smooth path leading to x¯ , an ε -minimizer of f , if for all u ∈ Rdim U small enough, it satisfies the following: (i) (ii) (iii) (iv)
Θε (uµ (x)) = x¯ + uµ (x) ⊕ vε (uµ (x)), where uµ (x) = (pµ,ε (x) − x¯ )U ; vε : Rdim U → Rdim V is a C 2 -function satisfying V¯ vε ∈ Wg¯ ,ε (u) for all g¯ ∈ ∂ f (¯x) ∩ ri∂ε f (¯x); the Jacobian J Θε (u) is a basis matrix for V (Θε (u))⊥ ; the particular approximate U-Lagrangian L0,ε (u) is a C 2 -function.
Accordingly, we have a dual path denoted by Γε (u) corresponding to smooth path. More precisely,
Γε (u) = arg min{kgε k2 : gε ∈ ∂ε f (Θε (u))}.
(6)
In fact, the smooth path and dual path are approximation of primal–dual track in [11]. If ε = 0, they are exactly primal–dual track. As is shown in [11], if 0 ∈ ri∂ f (¯x) and f has a pdg structure about x¯ satisfying strong transversality, then f has a primal–dual track. Hence, f has smooth path and dual path. The next lemma addresses that making an approximate V -step in Algorithm I essentially amounts to finding a corresponding smooth path point.
¯ := ∇ 2 L0,ε (0). Suppose 0 ∈ ∂ f (¯x) ∩ Lemma 3.1. Let Θε (uµ (x)) be a smooth path leading to x¯ , an ε -minimizer of f , and let H ri∂ε f (¯x). Then for all u sufficiently small Θε (uµ (x)) is the unique ε -minimizer of f on the affine set Θε (uµ (x)) + V (Θε (uµ (x))). Proof. Since J Θε (u) is a basis for V (Θε (u))⊥ , Theorem 3.4 in [9] with BU (u) = J Θε (u) gives the result.
3.2. The proximal bundle-type subroutine The study of approximate subgradients of convex function is deserve since in some cases a subgradient g (x) ∈ ∂ f (x) is expensive to compute. But if we know an already computed subgradient g (y) ∈ ∂ f (y), where y is near to x, then we have g (y) ∈ ∂ε f (x) because f (x) + g (y)T (z − x) = f (y) + g (y)T (z − y) + ε
≤ f (z ) + ε,
∀z ∈ Rn ,
where ε = f (x) − f (y) − g (y)T (x − y) ≥ 0. On the basis of the above observation, we attempt to explore the possibility of utilizing the approximate subgradients instead of the exact values. We approximate f from below by a piecewise linear convex function ϕ of the form:
ϕ(z ) := f (x) + max{gεjT (z − x) − αj − ε}, j∈B
z ∈ Rn ,
(7)
where B is some index set containing an index j such that yj = x, gεj ∈ ∂ε f (yj ) and αj = f (x)− f (yj )− gεjT (x − yj ) ≥ −ε, j ∈ B . Since gεj ∈ ∂ε f (yj ), j ∈ B , we have f (z ) ≥ f (yj ) + gεjT (z − yj ) − ε,
∀z ∈ Rn .
Adding 0 = αj − αj to this inequality gives f (z ) ≥ f (x) + gεjT (z − x) − αj − ε,
∀z ∈ Rn ,
which means that gεj is an (αj + ε)-subgradient of f at x, i.e., gεj ∈ ∂αj +ε f (x), j ∈ B . Since (7) becomes more and more crude if an approximation of f farther away from x, we add the proximal term 1 µkp − xk2 , µ > 0, to it. To approximate an ε -proximal point, we solve the first quadratic programming subproblem 2
(
min s.t.
1
µkp − xk2 2 r ≥ f (x) + gεjT (p − x) − αj − ε, r+
j ∈ B.
(8)
662
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Its corresponding dual problem is
min s.t.
2
X 1 X
λj gεjT + λj α j
2µ j∈B X j∈B λj = 1, λj ≥ 0, j ∈ B .
(9)
j∈B
ˆ = (λˆ 1 , . . . , λˆ |B | ) denote the optimal solution of (8) and (9), it is easily seen that Let (ˆr , pˆ ) and λ X 1 λˆ j gεj . rˆ = ϕ(ˆp) and pˆ = x − gˆε , where gˆε := µ j∈B
(10)
ˆ j = 0 for all j ∈ B such that r > f (x) + gεjT (ˆp − x) − αj − ε and In addition, λ ϕ(ˆp) = f (x) +
X
= f ( x) −
X
λˆ j (gεjT (ˆp − x) − αj − ε)
j∈B
λˆ j αj − ε −
j∈B
1
µ
kˆgε k2 .
(11)
The vector pˆ is an estimate of an approximate ε -proximal point. Hence, approximates a smooth path point when the latter j exists. To proceed further we let yj+ := pˆ and compute f (ˆp), gε+ ∈ ∂ε f (ˆp) and eˆ := f (ˆp) − ϕ(ˆp) = f (ˆp) − rˆ . An approximate dual path point, denoted by sˆ, is constructed by solving a second quadratic problem, which depends on a new index set
Bˆ := {j ∈ B : rˆ = f (x) + gεjT (ˆp − x) − αj − ε} ∪ {j+ }.
(12)
The second quadratic programming problem is
min
s.t.
1
kp − xk2 2 r ≥ gεjT (p − x),
r+
(13) j ∈ Bˆ .
It has a dual problem
min s.t.
2
j
λ g j ε
2
j∈Bˆ X λj = 1, λj ≥ 0, j ∈ Bˆ . 1 X
(14)
j∈Bˆ
¯ , satisfy Similar to (10), the respective solutions, denoted by (¯r , p¯ ) and λ r¯ = ϕ(¯p)
and p¯ − x = −ˆs,
where sˆ =
X
λ¯ j gεj .
(15)
j∈Bˆ
Given σ ∈ 0, if eˆ ≤
1 2
, the proximal bundle subprocedure is terminated and pˆ is declared to be a σ -approximation of pµ,ε (x)
σ kˆsk2 . µ
(16)
Otherwise, B above is replaced by Bˆ and new iterate data are computed by solving updated subproblems (8) and (13). This j update, appending (αj+ , gε+ ) to active data at (8), ensures convergence to an ε -minimizing point x¯ in case of nontermination. Remark 3.1. From the remarks above, the following results are true: (i) sˆ = arg min{ksk2 : s ∈ co{gεj : j ∈ Bˆ }}; (ii) since pµ,ε (x) is a smooth path point Θε (uµ (x)) approximated by pˆ and co{gεj : j ∈ Bˆ } approximates ∂ε f (Θε (uµ (x))), from (6) the corresponding Γε (uµ (x)) is estimated by sˆ; (iii) we can obtain the Uˆ by means of the following iteration. Let
Bˆ act := {j ∈ Bˆ : r¯ = gεjT (¯p − x)}. Then, from (15), r¯ = −gεjT sˆ, j ∈ Bˆ act , so
(gεj − gεl )T sˆ = 0,
(17)
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
663
for all such j and for a fixed l ∈ Bˆ act . Define a full column rank matrix Vˆ by choosing the largest number of indices j satisfying (17) such that the corresponding vectors gεj − gεl are linearly independent and by letting these vectors be the
columns of Vˆ . Then let Uˆ be a matrix whose columns form an orthonormal basis for the null-space of Vˆ T with Uˆ = I if Vˆ is vacuous. Theorem 3.1. Each iteration of the above proximal bundle subprocedure satisfies the following: (i) (ii) (iii) (iv) (v)
gεj ∈ ∂ε+ˆe f (ˆp), j ∈ Bˆ ; sˆ ∈ ∂ε+ˆe f (ˆp) and gˆε ∈ ∂ε+ˆe f (ˆp); µkˆp − pµ,ε (x)k2 ≤ 2ε + eˆ ; kˆsk ≤ kˆgε k, where gˆε = µ(x − pˆ ); for any parameter m ∈ (0, 1), (16) implies f (ˆp) − f (x) ≤ −
m 2µ
kˆgε k2 .
(18)
j
j
Proof. (i) Since gε+ ∈ ∂ε f (ˆp) and eˆ = f (ˆp) − ϕ(ˆp) ≥ 0, gε+ ∈ ∂ε+ˆe f (ˆp), so the result of item (i) holds for j = j+ . From the definition of pˆ , rˆ and Bˆ we have that for all j 6= j+ in Bˆ
ϕ(ˆp) = rˆ = f (x) + gεjT (ˆp − x) − αj − ε, so for all such j, eˆ = f (ˆp) − ϕ(ˆp) = f (ˆp) − f (x) − gεjT (ˆp − x) + αj + ε. In addition, f (z ) ≥ f (x) + gεjT (z − x) − αj − ε,
z ∈ Rn .
Adding 0 = eˆ − eˆ to this inequality gives f (z ) ≥ f (ˆp) + gεjT (z − pˆ ) − eˆ
≥ f (ˆp) + gεjT (z − pˆ ) − eˆ − ε,
∀z ∈ Rn ,
(19)
which means that gεj ∈ ∂ε+ˆe f (ˆp) for j+ 6= j ∈ Bˆ . ¯ j ≥ 0 and summing these inequalities, we have (ii) Multiplying each inequality in (19) by its corresponding multiplier λ
X j∈Bˆ
λ¯ j f (z ) ≥
X
λ¯ j f (ˆp) +
j∈Bˆ
X
λ¯ j gεjT (z − pˆ ) −
j∈Bˆ
λ¯ j eˆ −
j∈Bˆ
Using the definition of sˆ from (15) and the fact that f (z ) ≥ f (ˆp) + sˆT (z − pˆ ) − eˆ − ε,
X
P
j∈Bˆ
X
λ¯ j ε,
∀z ∈ Rn .
j∈Bˆ
λ¯ j = 1 gives
∀z ∈ Rn ,
ˆ j that solve dual problem (9) and define which means that sˆ ∈ ∂ε+ˆe f (ˆp). In a similar manner, this time using the multipliers λ ˆ j+ := 0, obtain the result. gˆε in (10), together with λ (iii) Since gµ,ε ∈ ∂ε f (pµ,ε (x)), we have
f (pµ,ε (x)) ≤ f (ˆp) − gµ,ε (x)T (ˆp − pµ,ε (x)) + ε. From (ii): gˆε ∈ ∂ε+ˆe f (ˆp), we get f (pµ,ε (x)) ≥ f (ˆp) + gˆεT (pµ,ε (x) − pˆ ) − ε − eˆ . Therefore, f (ˆp) + gˆεT (pµ,ε (x) − pˆ ) − ε − eˆ ≤ f (pµ,ε (x)) ≤ f (ˆp) − gµ,ε (x)T (ˆp − pµ,ε ) + ε, i.e., (ˆgε − gµ,ε (x))T (pµ,ε (x) − pˆ ) − 2ε − eˆ ≤ 0. Then, since the expression for gˆε from (10) written in the form gˆε = −µ(ˆp − x),
(20)
combined with gµ,ε (x) = µ(x − pµ,ε (x)) implies that gˆε − gµ,ε (x) = µ(pµ,ε (x) − pˆ ), we obtain item (iii). (iv) From (9), (10), (20) and the definition of Bˆ , we have µ(x − pˆ ) = gˆε is in the convex hull of {gεj , j ∈ Bˆ }. We obtain the result by virtue of the minimum norm property of sˆ.
664
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
(v) Since σ ≤
1 2
and m ∈ (0, 1), we have σ ≤ 1 −
ˆ j αj gives definition of eˆ , (22) and the nonnegativity of λ
m . 2
Thus if (16) holds then eˆ ≤ [(1 −
m 2
)/µ]kˆsk2 . Together with the
f (ˆp) − f (x) = eˆ + ϕ(ˆp) − f (x) X 1 = eˆ − λˆ j αj − kˆgε k2 − ε j∈B
µ
1
≤ eˆ − kˆgε k2 µ m 1 ≤ 1− µ kˆsk2 − kˆgε k2 . 2 µ Finally, combining this inequality with item (iv) gives (18).
3.3. Approximate decomposition algorithm and convergence analysis Substituting the approximate V -step in Algorithm I with proximal bundle-type subroutine, we present an approximate decomposition algorithm as follows. Afterwards a detailed convergence analysis is given. The main statement comprises the fact that each cluster point of the sequence of iterates generated by the algorithm is an optimal solution. Algorithm II. Begin of Algorithm Choose a starting point p0 ∈ Rn and positive parameters η, ε0 , µ, τ and m with τ < 1, m < 1. Step 0 Compute gε00 ∈ ∂ε0 f (p0 ). Let U0 be a matrix with orthonormal n-dimensional columns estimating an optimal U-basis.
Set s0 = gε00 and k := 0. Step 1 Stop if ksk k2 ≤ η. Step 2 Choose an nk × nk positive definite matrix Hk , where nk is the number of columns of Uk . Step 3 Compute an approximate U-Newton step by solving the linear system Hk 1uk = −UkT sk .
Set xk+1 := pk + Uk 1uk . Step 4 Choose µk+1 > µ, σk+1 ∈ 0, 21 , initialize B and run the bundle subprocedure with x = x0k+1 . Compute recursively, and set (e0k+1 , p0k+1 , s0k+1 , Uk0 +1 ) := (ˆe, pˆ , sˆ, Uˆ ). Step 5 If m f (p0k+1 ) − f (pk ) ≤ − ks0k+1 k2 , 2µk+1 0
then set
(xk+1 , ek+1 , pk+1 , sk+1 , Uk+1 , εk+1 ) := (x0k+1 , e0k+1 , p0k+1 , s0k+1 , Uk0 +1 , εk ). Otherwise, execute a line search xk+1 := arg min{f (pk ), f (p0k+1 )}; reinitialize B and rerun the bundle subroutine with x = xk+1 , to find new values for (ˆe, pˆ , sˆ, Uˆ ), then set (ek+1 , pk+1 , sk+1 , Uk+1 , εk+1 ) = (ˆe, pˆ , sˆ, Uˆ , τ εk ). Step 6 Replace k by k + 1 and go to Step 1. End of Algorithm ∞ Remark 3.2. In this algorithm, {εk }∞ k=0 ↓ 0. If {εk }k=0 ≡ 0, this algorithm is the same as Algorithm 6 in [11]. However, this algorithm uses approximate subgradients and proximal bundle-type subroutine which can deal with the approximate subgradients.
Theorem 3.2. One of the following two cases is true: (i) if the proximal bundle procedure in Algorithm II does not terminate, i.e., if (16) never hold, then the sequence of pˆ -values converges to pµ,ε (x) and pµ,ε (x) is a minimizer of f ; (ii) if the procedure terminates with sˆ = 0, then the corresponding pˆ equals pµ,ε (x) and is a minimizer of f . Proof. In [10, Prop. 4.3], if this procedure does not terminate then it generates an infinite sequence of eˆ -values converging to zero. Since (16) does not hold, the sequence of kˆsk-values also converges to 0. Thus, item (iii) in Theorem 3.1 implies that {ˆp} → pµ,ε (x). And Theorem 3.1(ii) gives f (z ) ≥ f (ˆp) + sˆT (z − pˆ ) − (ˆe + ε),
z ∈ Rn .
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
665
Table 1 Problem data. Problem
n
dim V
dim U
f¯
x¯
Starting point
F2d F3d-U3 F3d-U2 F3d-U1 F3d-U0
2 3 3 3 3
1 0 1 2 3
1 3 2 1 0
0 0 0 0 0
(0, 0) (0, 1, 10) (0, 0, 10) (0, 0, 0) (1, 0, 0)
x¯ + (0.9, 1.9) x¯ + (100, 33, −100) x¯ + (100, 33, −100) x¯ + (100, 33, −100) x¯ + (100, 33, −100)
By the continuity of f , this becomes f (z ) ≥ f (pµ,ε (x)) − ε,
z ∈ Rn .
The termination case with sˆ = 0 follows in a similar manner, since (16) implies eˆ = 0 in this case.
The next theorem establishes the convergence of Algorithm II, and the proof of which is similar to Theorem 9 in [11]. Theorem 3.3. Suppose that the sequence {µk } in Algorithm II is bounded above by µ ¯ . Then the following hold: (i) the sequence {f (pk )} is decreasing and either {f (pk )} → −∞ or {ksk k} and {ek } both converge to 0; (ii) if f is bounded from below, then any accumulation point of {pk } is a minimizer of f . Proof. In this paper, the inequalities of (20), (21), (22) in [11] become f (pk+1 ) − f (pk ) ≤ −
m 2µk+1
ksk+1 k2 ,
f (z ) ≥ f (pk ) + sTk (z − pk ) − εk − ek ,
(21)
∀z ∈ Rn ,
(22)
and ek ≤
σk ks k k2 . µk
(23)
Since ksk k 6= 0, (21) implies that {f (pk )} is decreasing. Suppose {f (pk )} 9 −∞. Then summing (21) over k and using the ≥ 2mµ¯ for all k implies that {ksk k} → 0. Then (23) with σk ≤ 12 and µk ≥ µ > 0 implies that {ek } → 0, which fact that 2m µk establishes (i). Now suppose f is bounded from below and p¯ is any accumulation point of {pk }. Then, because {ksk k} and {ek } converge to 0 by item (i), (22) together with the continuity of f implies that f (¯p) ≤ f (z ) for all z ∈ Rn and (ii) is proved. 4. An illustrative numerical example Now we report numerical result to illustrate Algorithm II. Our numerical experiment is carried out in Matlab 7.1 running on a PC Intel Pentium IV of 1.70 GHz CPU and 256 MB memory. All test examples are of the form f = max fi , i∈I
where I is finite and each fi is C 2 on Rn . For our runs we used the following examples:
• F2d: the objective function is given in [14], defined for x ∈ Rn by 1 2 2 F2d(x) := max (x1 + x2 ) − x2 , x2 ; 2
• F3d-U v : four functions of three variables, where v = 3, 2, 1, 0 denotes the corresponding dimension of the U-subspace. Given e := (0, 1, 1)T and four parameter vectors β v ∈ R4 , for x ∈ R3 1 2 F3d-U v(x) := max (x1 + x22 + 0.1x23 ) − eT x − β1v , x21 − 3x1 − β2v , x2 − β3v , x2 − β4v , 2
where β := (−5.5, 10, 11, 20), β 2 := (−5, 10, 0, 10), β 1 := (0, 10, 0, 0) and β 0 := (0.5, −2, 0, 0). 3
In Table 1, we show some relevant data for the problems, including the dimensions of V and U, the (known) optimal values and solutions and the starting points. We calculate an ε -subgradient at x by using the method in [12]: gε (x) = λg (x)+(1 −λ)g (x1 ), where g (x) is a subgradient at x and g (x1 ) is a subgradient at a point x1 such that 0 < α0 (x, x1 ) = f (x) − f (x1 ) − g (x1 )T (x − x1 ) ≤ ε . Here x1 ∈ Br (x) = {z ∈ Rn | kz − xk ≤ r } and λ ∈ [0, 1] are randomly chosen. The radius r is adjusted iteratively in the
666
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Table 2 Summary of the results. F2d
F3d-U3
F3d-U2
F3d-U1
F3d-U0
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
37 12
42 9
33 –
43 –
#f /g 34 Acc 9
20 9
5 16
24 10
15 13
20 12
following way: If we find the linearization error α0 (x, x1 ) > ε then r is reduced by a multiple smaller than one. On the other hand, if α0 (x, x1 ) is significantly smaller than ε , then r is increased by a multiple greater than one. When sk = sˆ in the P ¯ j H j )Uk , where H j = ∇ 2 fij (yj ), ij is an algorithm, then U-Hessian at x is computed in the following form: Hk = UkT ( j∈Bˆ λ
¯ j correspond to sˆ via (15). active index such that fij (yj ) = f (yj ), λ
The parameters have values η = 1.0 × 10−4 , m = 1.0 × 10−1 , τ = 1.0 × 10−1 and U0 equal to the n × n identity matrix. As for σk , µk , one can refer to [11]. Table 2 shows computational results of Algorithm II and Algorithm 6 in [11]. We use #f /g to denote the number of function and subgradient in Algorithm 6 and ε -subgradient in Algorithm II evaluation respectively. Acc stands for accuracy measure which is equal to the number of correct optimal objective value digits after the decimal point. For Algorithm II, F3d − U1(xk ) = −1.74 and F3d − U0(xk ) = −0.25, where k indicates the iteration in which the stopping criterion in Algorithm II is satisfied. This favorable results demonstrate that it is suitable to use approximate decomposition algorithm to solve (1) numerically. References [1] R.T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization 14 (1976) 877–898. [2] J.-B. Hirriart-Urruty, C. Lemaréchal, Convex Analysis and Minimization Algorithm, in: Grund. der Math. Wiss., two volumes, Springer-verlag, 1993, pp. 305–306. [3] S.J. Wright, Identifiable surfaces in constrained oprimization, SIAM Journal on Control and Optimization 31 (1993) 1063–1079. [4] A. Lewis, Active sets, nonsmoothness and sensitivity, SIAM Journal on Optimization 13 (2002) 702–725. [5] C. Lemaréchal, F. Oustry, C. Sagastizábal, The U-Lagrangian of a convex function, Transactions of the American Mathematical Society 352 (2000) 711–729. [6] R. Mifflin, C. Sagastizábal, VU-decomposition derivatives for convex max-functions, in: Tichatschke, Thra (Eds.), Ill-posed Variational Problems and Regularization Techniques, in: Lecture Notes in Economics and Mathematical Systems, vol. 477, 1999, pp. 167–186. [7] R. Mifflin, C. Sagastizábal, On VU-theory for functions with primal-dual gradient structure, SIAM Journal on Optimization 11 (2) (2000) 547–571. [8] R. Mifflin, C. Sagastizábal, VU-smoothness and proximal point results for some nonconvex functions, Optimization Methods and Software 19 (5) (2004) 463–478. [9] R. Mifflin, C. Sagastizábal, Proximal points are on the fast track, Journal of Convex Analysis 9 (2) (2002) 563–579. [10] R. Correa, C. Lemaréchal, Convergence of some algorithm for convex minimization, Mathematical Programming 62 (2) (1993) 261–275. [11] R. Mifflin, C. Sagastizábal, A VU-algorihtm for convex minimization, Mathematical Programming, Series B 104 (2005) 583–608. [12] M. Hintermuller, A proximal bundle method based on approximate subgradients, Computational Optimization and Applications 20 (2001) 245–266. [13] F. Shan, L.P. Pang, Z.Q. Xia, An approximate U-Lagrangian and algorithm to UV decomposition, Applied Mathematics and Computation 184 (2) (2007) 924–930. [14] C. Lemaréchal, C. Sagastizábal, Practical aspects of the Moreau–Yosida regulation: Theoretical preliminaries, SIAM Journal on Optimization 7 (2) (1997) 367–385.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An approximate decomposition algorithm for convex minimizationI Yuan Lu a , Li-Ping Pang b,∗ , Xi-Jun Liang b , Zun-Quan Xia b a
School of Sciences, Shenyang University, Shenyang 110044, China
b
School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China
article
info
Article history: Received 11 September 2009 Received in revised form 3 January 2010 MSC: 65Kxx Keywords: Nonsmooth convex optimization VU-decomposition Approximate U-Lagrangian Proximal bundle method Smooth path
abstract For nonsmooth convex optimization, Robert Mifflin and Claudia Sagastizábal introduce a VU-space decomposition algorithm in Mifflin and Sagastizábal (2005) [11]. An attractive property of this algorithm is that if a primal–dual track exists, this algorithm uses a bundle subroutine. With the inclusion of a simple line search, it is proved to be globally and superlinearly convergent. However, a drawback is that it needs the exact subgradients of the objective function, which is expensive to compute. In this paper an approximate decomposition algorithm based on proximal bundle-type method is introduced that is capable to deal with approximate subgradients. It is shown that the sequence of iterates generated by the resulting algorithm converges to the optimal solutions of the problem. Numerical tests emphasize the theoretical findings. © 2010 Elsevier B.V. All rights reserved.
1. Introduction In the area of continuous optimization many real-life problems can be written as min f (x), x∈Rn
(1)
with a convex and not necessarily differentiable function f : Rn → R. This fact gives reason for studying and developing algorithms to solve it numerically. If f is not differentiable at a solution x¯ of (1), constructing fast algorithm to compute x¯ is a challenge. Essentially, this has to do with the intrinsic difficulty in using appropriate second-order objects. One line of research [1] suggests introducing second-order information about f by means of its Moreau–Yosida regularization, which is a smooth function. Implementable forms of this method can be obtained by using a bundle technique, alternating serious steps with a sequence of null steps [2]. The other line of research parts from the viewpoint that nonsmoothness in practical applications is usually structured; see [3,4]. Nonsmooth functions may behave smoothly and even have appropriate secondorder representations on certain manifolds along certain directions. For the second line of research, new conceptual schemes have been developed, which are based on the VU-theory introduced in [5]; see also [6,7]. The idea is to decompose Rn into two orthogonal subspace V and U at a point x¯ that the nonsmoothness of f is concentrated essentially on V , and the smoothness of f appears on U-subspace. More precisely, for a given g¯ ∈ ∂ f (¯x), where ∂ f (¯x) denotes the subdifferential of f at x¯ in the sense of convex analysis. Then Rn can be decomposed as direct sum of two orthogonal subspaces, i.e., Rn = U ⊕ V , where V = lin(∂ f (¯x) − g¯ ), and U = V ⊥ . They define the U-Lagrangian, an approximation of the original function, which can be used to create a second-order expansion of f along certain manifolds. The manifolds can be identified and called primal track if f has a strongly transversal primal–dual
I The research is supported by the National Natural Science Foundation of China under project No. 10771026.
∗
Corresponding author. E-mail address: [email protected] (L.-P. Pang).
0377-0427/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2010.01.003
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
659
gradient structure; see [8]. The conceptual VU-algorithm in [5] finds points on primal track, by generating minimizing steps in the V -subspace. Alternating with these corrector steps are U-Newton predictor steps in order to provide for superlinear convergence. Accordingly, there is a dual track which provides the U-gradients needed for the U-Newton steps. Therefore, the study of approximation of primal–dual track is important. A first important step for this approximation is given in [9], where it is shown that, for a point near a minimizer, its corresponding proximal point lands on the primal track. This result is fundamental for implementation of the minimization algorithm, because a sequence of null steps from a bundle mechanism can approximate proximal points with any desired accuracy [10]. Bundle mechanism also can approximate a dual track by finding the minimum norm vector in the convex hull of the subgradients collected for estimating the proximal point. Based on the analysis above, Robert Mifflin and Claudia Sagastizábal design a VU-algorithm in [11]. This algorithm brings the iterate to the primal track with the help of bundle subroutine. Then the U-Newton step is performed to gain superlinear decrease of the distance to solution. However, this algorithm is conceptual in the sense that it makes use of the exact subgradients of the objective function, which is not easily to compute. The purpose of this paper is to explore the possibility of utilizing approximate subgradients instead of the exact ones. We present an approximate decomposition algorithm frame by defining the approximate U-Lagrangian. By introducing the definitions of smooth path and dual path, we obtain an implemental approximate decomposition algorithm based on proximal bundle subroutine [12], which can deal with approximate subgradients. This modification makes the algorithm easier to implement. The paper is organized as follows. In the next section we present the approximate U-Lagrangian. Based on this definition, we design a conceptual Algorithm I. It can deal with the approximate subgradients and only need to compute the ε -optimal solutions in approximate V -step. In the first part of Section 3 we propose a smooth path. The second part of Section 3 is devoted to establish an implemental Algorithm II which uses approximate subgradients and substitutes the approximate V -step in Algorithm I with proximal bundle subroutine. Numerical testing of resulting Algorithm II is reported in the final section. 2. A conceptual approximate decomposition algorithm 2.1. Approximate U-Lagrangian and its properties Since the U-Lagrangian itself is defined through a minimization problem involving f , exact evaluation of it in general practically is impossible. In this section we will briefly outline the definition of an approximation of U-Lagrangian and some related properties as well. For more details, one can refer to [13]. To begin with, we recall the concepts of ε -minimizer, ε -minimum and ε - min ϕ , where ϕ : Rp → R1 . We say that xε-min is an ε -minimizer of ϕ if xε-min ∈ Arg(ε - min ϕ)
= {x ∈ Rp | ϕ(z ) ≥ ϕ(x) − ε, ∀z ∈ Rp }, and that ϕε-min is an ε -minimum of ϕ if there exists an ε -minimizer xε-min such that
(2)
ϕ(z ) ≥ ϕε-min − ε, ∀z ∈ Rp , (3) where ϕε-min = ϕ(xε-min ), and that ε - min ϕ is ε -minimization of ϕ if minimizing ϕ is in the sense of (3). Under these concepts the definition of an approximation U-Lagrangian of f associated with some ε can be given below. Definition 2.1. The approximate U-Lagrangian of f with respect to ε (and g¯V ), denoted by Lg¯ ,ε (u), is defined as follows, Lg¯ ,ε (u) = ε - min{f (¯x + u ⊕ v) − h¯gV , viV },
(4)
v∈V
and Wg¯ ,ε (u) is defined by
Wg¯ ,ε (u) = Arg ε - min{f (¯x + u ⊕ v) − h¯gV , viV } . v∈V
Let ε > 0 be given. Introducing ε to [5, Prop. 2.4], one can restate it as the following lemma. Lemma 2.1. Suppose g¯ ∈ ∂ f (¯x)
T
ri∂ε f (¯x). Then Wg¯ ,ε (u) is nonempty.
Theorem 2.1. Suppose g¯ ∈ ∂ f (¯x) and g¯ ∈ ∂ f (¯x) ri∂ε f (¯x). Then the following assertions are true: (i) the function Lg¯ ,ε defined in (4) is convex and finite everywhere; (ii) Wg¯ ,ε (u) = {wε,u | ∃g ε ∈ ∂ε f (¯x + u ⊕ wε,u ) : gVε = g¯V }; (iii) 0 ∈ Wg¯ ,ε (0) and Lg¯ ,ε (0) = f (¯x).
T
Theorem 2.2. Suppose Wg¯ ,ε (u) 6= ∅. Then one has that
∂ Lg¯ ,ε (u) = {u∗ε | u∗ε ⊕ g¯V ∈ ∂ε f (¯x + u ⊕ wε,u )}, ∂ε Lg¯ ,ε (u) = {u∗ε | u∗ε ⊕ g¯V ∈ ∂2ε f (¯x + u ⊕ wε,u )}.
(5)
660
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Corollary 2.1. If 0 ∈ Wg¯ ,ε (0), then
∂ Lg¯ ,ε (0) = (∂ε f (¯x) − g¯ ) ∩ U. Theorem 2.3. Suppose g¯ ∈ ri∂ f (¯x). Then the following conclusion holds,
kwε,u k ≤ o(kuk) + ε,
∀wε,u ∈ Wg¯ ,ε (u).
2.2. Approximate decomposition algorithm frame It is well known that Hessian matrix plays a very important role in constructing numerical methods. In order to give an approximate decomposition algorithm frame, we restate the definition of Hessian matrix in [5] as follows. Definition 2.2. Assume that f is finite and x¯ and g¯ ∈ ∂ f (¯x) are fixed. We say that f has at x¯ a U-Hessian HU f (¯x) associated with g¯ if Lg¯ (u) has a generalized Hessian at 0, setting HU f (¯x) = HLg¯ (0). Based on the definition of approximate U-Lagrangian, we investigate an approximate decomposition algorithm frame. Algorithm I. Begin of Algorithm Step 0 Initiation x0 , η > 0, 0 < τ < 1, ε0 > 0, g 0 ∈ ∂ε0 f (x0 ), k = 0. Step 1 Stop if kg k k ≤ η. Step 2 Approximate V -step. Compute an ε -optimal solution δv ∈ V satisfying
δv ∈ Arg{ε - min f (xk + 0 ⊕ δv)}. δv∈V
Set x0 = xk + 0 ⊕ δv . Step 3 U-step. Make a Newton step in x0 . Compute g k ∈ ∂εk f (x0 ): gVk = 0,
k gU ∈ ∂εk Lg k ,εk ((x0 − x¯ )U ).
Compute the solution δ u ∈ U satisfying k gU + HU f (¯x)δ u = 0.
Set xk+1 = x0 + δ u ⊕ 0 = xk + δ u ⊕ δv and εk+1 = τ εk . Step 4 Update-set. Compute g k+1 ∈ ∂εk+1 f (xk+1 ). Set k = k + 1 and go to Step 1. End of algorithm Theorem 2.4. Assume that g¯ = 0 ∈ ∂ f (¯x) and f has a positive definite U-Hessian at x¯ , a minimizer of f . Then the iterate points {xk } constructed by Algorithm I satisfies
kxk+1 − x¯ k = o(kxk − x¯ k). Remark 2.1. If {εk }∞ k=1 ≡ 0, this algorithm is the same as Algorithm 4.5 in [5]. However, it uses approximate subgradients of the objective function and only needs to compute ε -optimal solutions in V -step, which make the algorithm easier to implement. 3. Approximate decomposition algorithm Since Algorithm I in Section 2 relies on knowing the subspaces U and V and converges only locally, it needs significant modification. In [9], Robert Mifflin and Claudia Sagastizábal show that a proximal point sequence follows primal track near a minimizer. This opens the way for defining a VU-algorithm where V -steps are replaced by proximal steps. In addition, the proximal step can be estimated with a bundle technique which also can approximate the unknown U and V subspaces as a computational by-product. Therefore, they establish Algorithm 6 in [11] by combining the bundle subroutine with the VU-space decomposition method. However, this algorithm needs the exact subgradients of the objective function, which is expensive to compute. Therefore, the study of using approximate subgradients instead of the exact ones is deserve.
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
661
3.1. Smooth path to ε -minimizer Given a positive scalar parameter µ, the ε -proximal point function depending on f , is defined by
1
pµ,ε (x) := arg ε − min f (p) +
2
p∈Rn
µkp − xk2
,
x ∈ Rn ,
where k · k stands for the Euclidean norm. It has the property: gµ,ε (x) := µ(x − pµ,ε (x)) ∈ ∂ε f (pµ,ε (x)). Similarly to the definition of primal track, we give the definition of smooth path. Definition 3.1. For any ε > 0, µ = µ(x) > 0, we say that Θε (u) = x¯ + u ⊕ vε (u) is a smooth path leading to x¯ , an ε -minimizer of f , if for all u ∈ Rdim U small enough, it satisfies the following: (i) (ii) (iii) (iv)
Θε (uµ (x)) = x¯ + uµ (x) ⊕ vε (uµ (x)), where uµ (x) = (pµ,ε (x) − x¯ )U ; vε : Rdim U → Rdim V is a C 2 -function satisfying V¯ vε ∈ Wg¯ ,ε (u) for all g¯ ∈ ∂ f (¯x) ∩ ri∂ε f (¯x); the Jacobian J Θε (u) is a basis matrix for V (Θε (u))⊥ ; the particular approximate U-Lagrangian L0,ε (u) is a C 2 -function.
Accordingly, we have a dual path denoted by Γε (u) corresponding to smooth path. More precisely,
Γε (u) = arg min{kgε k2 : gε ∈ ∂ε f (Θε (u))}.
(6)
In fact, the smooth path and dual path are approximation of primal–dual track in [11]. If ε = 0, they are exactly primal–dual track. As is shown in [11], if 0 ∈ ri∂ f (¯x) and f has a pdg structure about x¯ satisfying strong transversality, then f has a primal–dual track. Hence, f has smooth path and dual path. The next lemma addresses that making an approximate V -step in Algorithm I essentially amounts to finding a corresponding smooth path point.
¯ := ∇ 2 L0,ε (0). Suppose 0 ∈ ∂ f (¯x) ∩ Lemma 3.1. Let Θε (uµ (x)) be a smooth path leading to x¯ , an ε -minimizer of f , and let H ri∂ε f (¯x). Then for all u sufficiently small Θε (uµ (x)) is the unique ε -minimizer of f on the affine set Θε (uµ (x)) + V (Θε (uµ (x))). Proof. Since J Θε (u) is a basis for V (Θε (u))⊥ , Theorem 3.4 in [9] with BU (u) = J Θε (u) gives the result.
3.2. The proximal bundle-type subroutine The study of approximate subgradients of convex function is deserve since in some cases a subgradient g (x) ∈ ∂ f (x) is expensive to compute. But if we know an already computed subgradient g (y) ∈ ∂ f (y), where y is near to x, then we have g (y) ∈ ∂ε f (x) because f (x) + g (y)T (z − x) = f (y) + g (y)T (z − y) + ε
≤ f (z ) + ε,
∀z ∈ Rn ,
where ε = f (x) − f (y) − g (y)T (x − y) ≥ 0. On the basis of the above observation, we attempt to explore the possibility of utilizing the approximate subgradients instead of the exact values. We approximate f from below by a piecewise linear convex function ϕ of the form:
ϕ(z ) := f (x) + max{gεjT (z − x) − αj − ε}, j∈B
z ∈ Rn ,
(7)
where B is some index set containing an index j such that yj = x, gεj ∈ ∂ε f (yj ) and αj = f (x)− f (yj )− gεjT (x − yj ) ≥ −ε, j ∈ B . Since gεj ∈ ∂ε f (yj ), j ∈ B , we have f (z ) ≥ f (yj ) + gεjT (z − yj ) − ε,
∀z ∈ Rn .
Adding 0 = αj − αj to this inequality gives f (z ) ≥ f (x) + gεjT (z − x) − αj − ε,
∀z ∈ Rn ,
which means that gεj is an (αj + ε)-subgradient of f at x, i.e., gεj ∈ ∂αj +ε f (x), j ∈ B . Since (7) becomes more and more crude if an approximation of f farther away from x, we add the proximal term 1 µkp − xk2 , µ > 0, to it. To approximate an ε -proximal point, we solve the first quadratic programming subproblem 2
(
min s.t.
1
µkp − xk2 2 r ≥ f (x) + gεjT (p − x) − αj − ε, r+
j ∈ B.
(8)
662
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Its corresponding dual problem is
min s.t.
2
X 1 X
λj gεjT + λj α j
2µ j∈B X j∈B λj = 1, λj ≥ 0, j ∈ B .
(9)
j∈B
ˆ = (λˆ 1 , . . . , λˆ |B | ) denote the optimal solution of (8) and (9), it is easily seen that Let (ˆr , pˆ ) and λ X 1 λˆ j gεj . rˆ = ϕ(ˆp) and pˆ = x − gˆε , where gˆε := µ j∈B
(10)
ˆ j = 0 for all j ∈ B such that r > f (x) + gεjT (ˆp − x) − αj − ε and In addition, λ ϕ(ˆp) = f (x) +
X
= f ( x) −
X
λˆ j (gεjT (ˆp − x) − αj − ε)
j∈B
λˆ j αj − ε −
j∈B
1
µ
kˆgε k2 .
(11)
The vector pˆ is an estimate of an approximate ε -proximal point. Hence, approximates a smooth path point when the latter j exists. To proceed further we let yj+ := pˆ and compute f (ˆp), gε+ ∈ ∂ε f (ˆp) and eˆ := f (ˆp) − ϕ(ˆp) = f (ˆp) − rˆ . An approximate dual path point, denoted by sˆ, is constructed by solving a second quadratic problem, which depends on a new index set
Bˆ := {j ∈ B : rˆ = f (x) + gεjT (ˆp − x) − αj − ε} ∪ {j+ }.
(12)
The second quadratic programming problem is
min
s.t.
1
kp − xk2 2 r ≥ gεjT (p − x),
r+
(13) j ∈ Bˆ .
It has a dual problem
min s.t.
2
j
λ g j ε
2
j∈Bˆ X λj = 1, λj ≥ 0, j ∈ Bˆ . 1 X
(14)
j∈Bˆ
¯ , satisfy Similar to (10), the respective solutions, denoted by (¯r , p¯ ) and λ r¯ = ϕ(¯p)
and p¯ − x = −ˆs,
where sˆ =
X
λ¯ j gεj .
(15)
j∈Bˆ
Given σ ∈ 0, if eˆ ≤
1 2
, the proximal bundle subprocedure is terminated and pˆ is declared to be a σ -approximation of pµ,ε (x)
σ kˆsk2 . µ
(16)
Otherwise, B above is replaced by Bˆ and new iterate data are computed by solving updated subproblems (8) and (13). This j update, appending (αj+ , gε+ ) to active data at (8), ensures convergence to an ε -minimizing point x¯ in case of nontermination. Remark 3.1. From the remarks above, the following results are true: (i) sˆ = arg min{ksk2 : s ∈ co{gεj : j ∈ Bˆ }}; (ii) since pµ,ε (x) is a smooth path point Θε (uµ (x)) approximated by pˆ and co{gεj : j ∈ Bˆ } approximates ∂ε f (Θε (uµ (x))), from (6) the corresponding Γε (uµ (x)) is estimated by sˆ; (iii) we can obtain the Uˆ by means of the following iteration. Let
Bˆ act := {j ∈ Bˆ : r¯ = gεjT (¯p − x)}. Then, from (15), r¯ = −gεjT sˆ, j ∈ Bˆ act , so
(gεj − gεl )T sˆ = 0,
(17)
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
663
for all such j and for a fixed l ∈ Bˆ act . Define a full column rank matrix Vˆ by choosing the largest number of indices j satisfying (17) such that the corresponding vectors gεj − gεl are linearly independent and by letting these vectors be the
columns of Vˆ . Then let Uˆ be a matrix whose columns form an orthonormal basis for the null-space of Vˆ T with Uˆ = I if Vˆ is vacuous. Theorem 3.1. Each iteration of the above proximal bundle subprocedure satisfies the following: (i) (ii) (iii) (iv) (v)
gεj ∈ ∂ε+ˆe f (ˆp), j ∈ Bˆ ; sˆ ∈ ∂ε+ˆe f (ˆp) and gˆε ∈ ∂ε+ˆe f (ˆp); µkˆp − pµ,ε (x)k2 ≤ 2ε + eˆ ; kˆsk ≤ kˆgε k, where gˆε = µ(x − pˆ ); for any parameter m ∈ (0, 1), (16) implies f (ˆp) − f (x) ≤ −
m 2µ
kˆgε k2 .
(18)
j
j
Proof. (i) Since gε+ ∈ ∂ε f (ˆp) and eˆ = f (ˆp) − ϕ(ˆp) ≥ 0, gε+ ∈ ∂ε+ˆe f (ˆp), so the result of item (i) holds for j = j+ . From the definition of pˆ , rˆ and Bˆ we have that for all j 6= j+ in Bˆ
ϕ(ˆp) = rˆ = f (x) + gεjT (ˆp − x) − αj − ε, so for all such j, eˆ = f (ˆp) − ϕ(ˆp) = f (ˆp) − f (x) − gεjT (ˆp − x) + αj + ε. In addition, f (z ) ≥ f (x) + gεjT (z − x) − αj − ε,
z ∈ Rn .
Adding 0 = eˆ − eˆ to this inequality gives f (z ) ≥ f (ˆp) + gεjT (z − pˆ ) − eˆ
≥ f (ˆp) + gεjT (z − pˆ ) − eˆ − ε,
∀z ∈ Rn ,
(19)
which means that gεj ∈ ∂ε+ˆe f (ˆp) for j+ 6= j ∈ Bˆ . ¯ j ≥ 0 and summing these inequalities, we have (ii) Multiplying each inequality in (19) by its corresponding multiplier λ
X j∈Bˆ
λ¯ j f (z ) ≥
X
λ¯ j f (ˆp) +
j∈Bˆ
X
λ¯ j gεjT (z − pˆ ) −
j∈Bˆ
λ¯ j eˆ −
j∈Bˆ
Using the definition of sˆ from (15) and the fact that f (z ) ≥ f (ˆp) + sˆT (z − pˆ ) − eˆ − ε,
X
P
j∈Bˆ
X
λ¯ j ε,
∀z ∈ Rn .
j∈Bˆ
λ¯ j = 1 gives
∀z ∈ Rn ,
ˆ j that solve dual problem (9) and define which means that sˆ ∈ ∂ε+ˆe f (ˆp). In a similar manner, this time using the multipliers λ ˆ j+ := 0, obtain the result. gˆε in (10), together with λ (iii) Since gµ,ε ∈ ∂ε f (pµ,ε (x)), we have
f (pµ,ε (x)) ≤ f (ˆp) − gµ,ε (x)T (ˆp − pµ,ε (x)) + ε. From (ii): gˆε ∈ ∂ε+ˆe f (ˆp), we get f (pµ,ε (x)) ≥ f (ˆp) + gˆεT (pµ,ε (x) − pˆ ) − ε − eˆ . Therefore, f (ˆp) + gˆεT (pµ,ε (x) − pˆ ) − ε − eˆ ≤ f (pµ,ε (x)) ≤ f (ˆp) − gµ,ε (x)T (ˆp − pµ,ε ) + ε, i.e., (ˆgε − gµ,ε (x))T (pµ,ε (x) − pˆ ) − 2ε − eˆ ≤ 0. Then, since the expression for gˆε from (10) written in the form gˆε = −µ(ˆp − x),
(20)
combined with gµ,ε (x) = µ(x − pµ,ε (x)) implies that gˆε − gµ,ε (x) = µ(pµ,ε (x) − pˆ ), we obtain item (iii). (iv) From (9), (10), (20) and the definition of Bˆ , we have µ(x − pˆ ) = gˆε is in the convex hull of {gεj , j ∈ Bˆ }. We obtain the result by virtue of the minimum norm property of sˆ.
664
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
(v) Since σ ≤
1 2
and m ∈ (0, 1), we have σ ≤ 1 −
ˆ j αj gives definition of eˆ , (22) and the nonnegativity of λ
m . 2
Thus if (16) holds then eˆ ≤ [(1 −
m 2
)/µ]kˆsk2 . Together with the
f (ˆp) − f (x) = eˆ + ϕ(ˆp) − f (x) X 1 = eˆ − λˆ j αj − kˆgε k2 − ε j∈B
µ
1
≤ eˆ − kˆgε k2 µ m 1 ≤ 1− µ kˆsk2 − kˆgε k2 . 2 µ Finally, combining this inequality with item (iv) gives (18).
3.3. Approximate decomposition algorithm and convergence analysis Substituting the approximate V -step in Algorithm I with proximal bundle-type subroutine, we present an approximate decomposition algorithm as follows. Afterwards a detailed convergence analysis is given. The main statement comprises the fact that each cluster point of the sequence of iterates generated by the algorithm is an optimal solution. Algorithm II. Begin of Algorithm Choose a starting point p0 ∈ Rn and positive parameters η, ε0 , µ, τ and m with τ < 1, m < 1. Step 0 Compute gε00 ∈ ∂ε0 f (p0 ). Let U0 be a matrix with orthonormal n-dimensional columns estimating an optimal U-basis.
Set s0 = gε00 and k := 0. Step 1 Stop if ksk k2 ≤ η. Step 2 Choose an nk × nk positive definite matrix Hk , where nk is the number of columns of Uk . Step 3 Compute an approximate U-Newton step by solving the linear system Hk 1uk = −UkT sk .
Set xk+1 := pk + Uk 1uk . Step 4 Choose µk+1 > µ, σk+1 ∈ 0, 21 , initialize B and run the bundle subprocedure with x = x0k+1 . Compute recursively, and set (e0k+1 , p0k+1 , s0k+1 , Uk0 +1 ) := (ˆe, pˆ , sˆ, Uˆ ). Step 5 If m f (p0k+1 ) − f (pk ) ≤ − ks0k+1 k2 , 2µk+1 0
then set
(xk+1 , ek+1 , pk+1 , sk+1 , Uk+1 , εk+1 ) := (x0k+1 , e0k+1 , p0k+1 , s0k+1 , Uk0 +1 , εk ). Otherwise, execute a line search xk+1 := arg min{f (pk ), f (p0k+1 )}; reinitialize B and rerun the bundle subroutine with x = xk+1 , to find new values for (ˆe, pˆ , sˆ, Uˆ ), then set (ek+1 , pk+1 , sk+1 , Uk+1 , εk+1 ) = (ˆe, pˆ , sˆ, Uˆ , τ εk ). Step 6 Replace k by k + 1 and go to Step 1. End of Algorithm ∞ Remark 3.2. In this algorithm, {εk }∞ k=0 ↓ 0. If {εk }k=0 ≡ 0, this algorithm is the same as Algorithm 6 in [11]. However, this algorithm uses approximate subgradients and proximal bundle-type subroutine which can deal with the approximate subgradients.
Theorem 3.2. One of the following two cases is true: (i) if the proximal bundle procedure in Algorithm II does not terminate, i.e., if (16) never hold, then the sequence of pˆ -values converges to pµ,ε (x) and pµ,ε (x) is a minimizer of f ; (ii) if the procedure terminates with sˆ = 0, then the corresponding pˆ equals pµ,ε (x) and is a minimizer of f . Proof. In [10, Prop. 4.3], if this procedure does not terminate then it generates an infinite sequence of eˆ -values converging to zero. Since (16) does not hold, the sequence of kˆsk-values also converges to 0. Thus, item (iii) in Theorem 3.1 implies that {ˆp} → pµ,ε (x). And Theorem 3.1(ii) gives f (z ) ≥ f (ˆp) + sˆT (z − pˆ ) − (ˆe + ε),
z ∈ Rn .
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
665
Table 1 Problem data. Problem
n
dim V
dim U
f¯
x¯
Starting point
F2d F3d-U3 F3d-U2 F3d-U1 F3d-U0
2 3 3 3 3
1 0 1 2 3
1 3 2 1 0
0 0 0 0 0
(0, 0) (0, 1, 10) (0, 0, 10) (0, 0, 0) (1, 0, 0)
x¯ + (0.9, 1.9) x¯ + (100, 33, −100) x¯ + (100, 33, −100) x¯ + (100, 33, −100) x¯ + (100, 33, −100)
By the continuity of f , this becomes f (z ) ≥ f (pµ,ε (x)) − ε,
z ∈ Rn .
The termination case with sˆ = 0 follows in a similar manner, since (16) implies eˆ = 0 in this case.
The next theorem establishes the convergence of Algorithm II, and the proof of which is similar to Theorem 9 in [11]. Theorem 3.3. Suppose that the sequence {µk } in Algorithm II is bounded above by µ ¯ . Then the following hold: (i) the sequence {f (pk )} is decreasing and either {f (pk )} → −∞ or {ksk k} and {ek } both converge to 0; (ii) if f is bounded from below, then any accumulation point of {pk } is a minimizer of f . Proof. In this paper, the inequalities of (20), (21), (22) in [11] become f (pk+1 ) − f (pk ) ≤ −
m 2µk+1
ksk+1 k2 ,
f (z ) ≥ f (pk ) + sTk (z − pk ) − εk − ek ,
(21)
∀z ∈ Rn ,
(22)
and ek ≤
σk ks k k2 . µk
(23)
Since ksk k 6= 0, (21) implies that {f (pk )} is decreasing. Suppose {f (pk )} 9 −∞. Then summing (21) over k and using the ≥ 2mµ¯ for all k implies that {ksk k} → 0. Then (23) with σk ≤ 12 and µk ≥ µ > 0 implies that {ek } → 0, which fact that 2m µk establishes (i). Now suppose f is bounded from below and p¯ is any accumulation point of {pk }. Then, because {ksk k} and {ek } converge to 0 by item (i), (22) together with the continuity of f implies that f (¯p) ≤ f (z ) for all z ∈ Rn and (ii) is proved. 4. An illustrative numerical example Now we report numerical result to illustrate Algorithm II. Our numerical experiment is carried out in Matlab 7.1 running on a PC Intel Pentium IV of 1.70 GHz CPU and 256 MB memory. All test examples are of the form f = max fi , i∈I
where I is finite and each fi is C 2 on Rn . For our runs we used the following examples:
• F2d: the objective function is given in [14], defined for x ∈ Rn by 1 2 2 F2d(x) := max (x1 + x2 ) − x2 , x2 ; 2
• F3d-U v : four functions of three variables, where v = 3, 2, 1, 0 denotes the corresponding dimension of the U-subspace. Given e := (0, 1, 1)T and four parameter vectors β v ∈ R4 , for x ∈ R3 1 2 F3d-U v(x) := max (x1 + x22 + 0.1x23 ) − eT x − β1v , x21 − 3x1 − β2v , x2 − β3v , x2 − β4v , 2
where β := (−5.5, 10, 11, 20), β 2 := (−5, 10, 0, 10), β 1 := (0, 10, 0, 0) and β 0 := (0.5, −2, 0, 0). 3
In Table 1, we show some relevant data for the problems, including the dimensions of V and U, the (known) optimal values and solutions and the starting points. We calculate an ε -subgradient at x by using the method in [12]: gε (x) = λg (x)+(1 −λ)g (x1 ), where g (x) is a subgradient at x and g (x1 ) is a subgradient at a point x1 such that 0 < α0 (x, x1 ) = f (x) − f (x1 ) − g (x1 )T (x − x1 ) ≤ ε . Here x1 ∈ Br (x) = {z ∈ Rn | kz − xk ≤ r } and λ ∈ [0, 1] are randomly chosen. The radius r is adjusted iteratively in the
666
Y. Lu et al. / Journal of Computational and Applied Mathematics 234 (2010) 658–666
Table 2 Summary of the results. F2d
F3d-U3
F3d-U2
F3d-U1
F3d-U0
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
Algorithm II Algorithm 6
37 12
42 9
33 –
43 –
#f /g 34 Acc 9
20 9
5 16
24 10
15 13
20 12
following way: If we find the linearization error α0 (x, x1 ) > ε then r is reduced by a multiple smaller than one. On the other hand, if α0 (x, x1 ) is significantly smaller than ε , then r is increased by a multiple greater than one. When sk = sˆ in the P ¯ j H j )Uk , where H j = ∇ 2 fij (yj ), ij is an algorithm, then U-Hessian at x is computed in the following form: Hk = UkT ( j∈Bˆ λ
¯ j correspond to sˆ via (15). active index such that fij (yj ) = f (yj ), λ
The parameters have values η = 1.0 × 10−4 , m = 1.0 × 10−1 , τ = 1.0 × 10−1 and U0 equal to the n × n identity matrix. As for σk , µk , one can refer to [11]. Table 2 shows computational results of Algorithm II and Algorithm 6 in [11]. We use #f /g to denote the number of function and subgradient in Algorithm 6 and ε -subgradient in Algorithm II evaluation respectively. Acc stands for accuracy measure which is equal to the number of correct optimal objective value digits after the decimal point. For Algorithm II, F3d − U1(xk ) = −1.74 and F3d − U0(xk ) = −0.25, where k indicates the iteration in which the stopping criterion in Algorithm II is satisfied. This favorable results demonstrate that it is suitable to use approximate decomposition algorithm to solve (1) numerically. References [1] R.T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization 14 (1976) 877–898. [2] J.-B. Hirriart-Urruty, C. Lemaréchal, Convex Analysis and Minimization Algorithm, in: Grund. der Math. Wiss., two volumes, Springer-verlag, 1993, pp. 305–306. [3] S.J. Wright, Identifiable surfaces in constrained oprimization, SIAM Journal on Control and Optimization 31 (1993) 1063–1079. [4] A. Lewis, Active sets, nonsmoothness and sensitivity, SIAM Journal on Optimization 13 (2002) 702–725. [5] C. Lemaréchal, F. Oustry, C. Sagastizábal, The U-Lagrangian of a convex function, Transactions of the American Mathematical Society 352 (2000) 711–729. [6] R. Mifflin, C. Sagastizábal, VU-decomposition derivatives for convex max-functions, in: Tichatschke, Thra (Eds.), Ill-posed Variational Problems and Regularization Techniques, in: Lecture Notes in Economics and Mathematical Systems, vol. 477, 1999, pp. 167–186. [7] R. Mifflin, C. Sagastizábal, On VU-theory for functions with primal-dual gradient structure, SIAM Journal on Optimization 11 (2) (2000) 547–571. [8] R. Mifflin, C. Sagastizábal, VU-smoothness and proximal point results for some nonconvex functions, Optimization Methods and Software 19 (5) (2004) 463–478. [9] R. Mifflin, C. Sagastizábal, Proximal points are on the fast track, Journal of Convex Analysis 9 (2) (2002) 563–579. [10] R. Correa, C. Lemaréchal, Convergence of some algorithm for convex minimization, Mathematical Programming 62 (2) (1993) 261–275. [11] R. Mifflin, C. Sagastizábal, A VU-algorihtm for convex minimization, Mathematical Programming, Series B 104 (2005) 583–608. [12] M. Hintermuller, A proximal bundle method based on approximate subgradients, Computational Optimization and Applications 20 (2001) 245–266. [13] F. Shan, L.P. Pang, Z.Q. Xia, An approximate U-Lagrangian and algorithm to UV decomposition, Applied Mathematics and Computation 184 (2) (2007) 924–930. [14] C. Lemaréchal, C. Sagastizábal, Practical aspects of the Moreau–Yosida regulation: Theoretical preliminaries, SIAM Journal on Optimization 7 (2) (1997) 367–385.