A Generalized Alternating Linearization Bundle Method for Structured Convex Optimization with Inexact First-Order Oracles

: In this paper, we consider a class of structured optimization problems whose objective function is the summation of two convex functions: f and h , which are not necessarily differentiable. We focus particularly on the case where the function f is general and its exact ﬁrst-order information (function value and subgradient) may be difﬁcult to obtain, while the function h is relatively simple. We propose a generalized alternating linearization bundle method for solving this class of problems, which can handle inexact ﬁrst-order information of on-demand accuracy. The inexact information can be very general, which covers various oracles, such as inexact, partially inexact and asymptotically exact oracles, and so forth. At each iteration, the algorithm solves two interrelated subproblems: one aims to ﬁnd the proximal point of the polyhedron model of f plus the linearization of h ; the other aims to ﬁnd the proximal point of the linearization of f plus h . We establish global convergence of the algorithm under different types of inexactness. Finally, some preliminary numerical results on a set of two-stage stochastic linear programming problems show that our method is very encouraging.


Introduction
In this paper, we consider the following structured convex optimization problem where f : dom h → R and h : R n → (−∞, ∞] are closed proper convex functions, but not necessarily differentiable, and dom h := {x : h(x) < ∞} is the effective domain of h. Problems of this type frequently arise in practice, such as compressed sensing [1], image reconstruction [2], machine learning [3], optimal control [4] and power system [5][6][7][8], and so forth. The following are three interesting examples.

Example 1.
( 1 minimization in compressed sensing). The signal recovery problems in compressed sensing [1] usually take the following form min x∈R n where A ∈ R m×n , b ∈ R m , λ > 0 and x 1 := ∑ n i=1 |x i |, which aims to get a sparse solution x of the linear system Ax = b. Note that by defining f (x) = 1 2 Ax − b 2 and h(x) = λ x 1 , (2) is of the form of (1).

Example 2. (Regularized risk minimization)
. At the core of many machine learning problems is to minimize a regularized risk function [9,10]: where R emp (x) := 1 is the empirical risk, {(u i , v i ), i = 1, · · · , m} is a training set, and l is a convex loss function measuring the gap between v and the predicted values generated by using x. In general, R emp (x) is a nondifferentiable and computationally expensive convex function, whereas the regularization term Ω(x) is a simple convex function, say Ω(x) = 1 2 x 2 2 or Ω(x) = x 1 . By defining f (x) = R emp (x) and h(x) = λΩ(x), (3) is also of the form of (1).

Example 3. (Unconstrained transformation of a constrained problem). Given a constrained problem
where f is a convex function and X ⊆ R n is a convex set. By introducing the indicator function δ X of X, that is, δ X (x) equals 0 on X and infinity elsewhere, problem (4) can be written equivalently as min x∈R n f (x) + δ X (x). (5) Clearly, by setting h(x) = δ X (x), (5) becomes the form of (1). We note that such transformation could be very efficient in practice if the set X has some special structure [11,12].
The design of methods for solving problems of the form (1) has attracted the attention of many researchers. We mention here four classes of these methods-operator splitting methods [13][14][15], alternating direction methods of multipliers [5,[16][17][18][19], alternating linearization methods [20,21], and alternating linearization bundle method [22]. They all fall within the well-known class of first-order black-box methods, that is, it is assumed that there is an oracle that can return the (approximate) function value and one arbitrary (approximate) subgradient at any given point. Regarding the above methods, we are concerned about the following three points:

•
Smoothness of one or both functions in the objective has been assumed for many of the methods. • Except for Reference [22], they all require the exact computation of the function values and (sub)gradients.

•
The alternating linearization methods [20,21] essentially assume that both functions f and h are "simple" in the sense that minimizing the function plus a separable convex quadratic function is easy.
However, for some practical problems, the functions may be nondifferentiable (nonsmooth), not easy to handle, and computationally expensive in the sense that the exact first-order information may be impossible to calculate, or be computable but difficult to obtain. For example, if f has the form where U is an infinite set and φ u (x) : R n → R is convex for any u ∈ U, then it is often difficult to calculate the exact function value f (x). But for any tolerance > 0, we may usually find a lower approximation f Then we can take a subgradient of φ u at x as an approximate subgradient of f at x. Another example is two-stage stochastic programming (see, e.g., References [23,24]), in which the function value is generated after solving a series of linear programs (details will be given in the section of numerical experiments), therefore its accuracy is determined by the tolerance of the linear programming solver. Some other practical examples, such as Lagrangian relaxation, chance-constrained programs and convex composite functions, can be found in Reference [25].
Based on the above observation, in this paper, we focus particularly on the case where the function f is general, possibly nonsmooth and its exact function values and subgradients may be difficult to obtain, whereas the function h is assumed to be relatively simple. Our main goal is to provide an efficient method, namely, the improved alternating linearization bundle method, for such kind of structured convex optimization problems. The basic tool we used here to handle nonsmoothness and inexactness is the bundle technique, since in the nonsmooth optimization community, bundle methods [26][27][28][29] are regarded as the most robust and reliable methods, whose variants have been well studied for handling inexact oracles [23,25,[30][31][32][33].
Roughly speaking, our method generalizes the alternating linearization bundle method of Kiwiel [22] from exact and inexact oracles to various oracles, including exact, inexact, partially inexact, asymptotically exact and partially asymptotically exact oracles. These oracles are covered by the so-called on-demand accuracy oracles proposed by de Oliveira and Sagastizábal [23], whose accuracy is controlled by two parameters: a descent target and an error bound. More precisely, the accuracy is bounded by an error bound whenever the function estimation reaches a certain descent target. The most advantage of oracles with on-demand accuracy is that the function and subgradient estimations can be rough without an accuracy control for some "non-critical" iterates, thus the computational effort can be saved.
At each iteration, the proposed algorithm alternately solves two interrelated subproblems: one is to find the proximal point of the polyhedron model of f plus the linearization of h; the other is to find the proximal point of the linearization of f plus h. We establish global convergence of the algorithm under different types of inexactness. Finally, some preliminary numerical results on a set of two-stage stochastic linear programming problems show that our method is very encouraging. This paper is organized as follows. In Section 2, we recall the condition of the inexact frist-order oracles. In Section 3, we present an improved alternating linearization bundle method for structured convex optimization with inexact first-order oracles and show some properties of the algorithm. In Section 4, we establish global convergence of the algorithm under different types of inexactness. In Section 5, we provide some numerical experiments on two-stage stochastic linear programming problems. The notations are standard. The Euclidean inner product in R n is denoted by x, y := x T y, and the associated norm by · .

Preliminaries
For a given constant ≥ 0, the -subdifferential of function f at x is defined by (see Reference [34]) being the usual subdifferential in convex analysis [35]. Each element g ∈ ∂ f (x) is called a subgradient. For simplicity, we use the following notations: Aiming at the special structure of problem (1), we present a slight variant of the oracles with on-demand accuracy proposed in Reference [23] as follows: for a given x ∈ R n , a descent target γ x and an error bound ε x ≥ 0, the approximate values f x , g x and F x satisfy the following condition From the relations in (6), we see that although the error η(γ x ) is unknown, it has to be restricted within the error bound ε x whenever the descent target F x ≤ γ x is reached. This ensures that the exact and inexact function values satisfy The advantages of oracle (6) are that: (1) if the descent target is not reached, the calculation of oracle information can be rough without an accuracy control, which can potentially reduce the computational cost; (2) by properly choosing the parameters γ x and ε x , the oracle (6) covers various existing oracles: • Exact (Ex) [12,21]: set γ x = +∞ and ε x = 0; • Partially Inexact (PI) [24]: set γ x < +∞ and ε x = 0; • Inexact (IE) [11,25,32,36,37]: set γ x = +∞ and ε x ≡ ε > 0 (possibly unknown); • Asymptotically Exact (AE) [38,39]: set γ x = +∞ and ε x → 0 along the iterative process; • Partially Asymptotically Exact (PAE) [23]: set γ x < +∞ and ε x → 0.

The Generalized Alternating Linearization Bundle Method
In this section, we present our generalized alternating linearization bundle method with inexact first-order oracles for solving (1).
Let k be the current iteration index, x j , j ∈ J k ⊆ {1, · · · , k} be given points generated by previous iterations, and the corresponding approximate values f x j /g x j be produced by the oracle (6). For notational convenience, we denote The approximate linearizations of f at x j are given by From the second relation in (6), we have which implies that f j is a lower approximation to f . Next, it is natural to define the polyhedral inexact cutting-plane model of f byf which is obviously a lower polyhedral model for f , that is,f k (·) ≤ f (·). Letx k (called stability center) be the "best" point obtained so far, which satisfies thatx (7), we have By applying the bundle idea to the "complex" function f , and keeping the simple function h unchanged, similar to traditional proximal bundle methods (see, e.g., Reference [28]), we may solve the following subproblem to obtain a new iterate x k+1 : where t k > 0 is a proximal parameter. However, subproblem (10) is generally not easy to solve, so by making use of the alternating linearization idea of Kiwiel [22], we solve two easier subproblems instead of (10). These two subproblems are interrelated: one is to find the proximal point of the polyhedron modelf plus the linearization of h, aiming at generating an aggregate linear model of f for use in the second subproblem; the other is to find the proximal point of the aggregate linear model of f plus h, aiming at obtaining a new trial point. Now, we are ready to present the details of our algorithm, which generalizes the work of Kiwiel [22]. We note that the choice of the model functionf k in the algorithm may be different from the form of (8), since the subgradient aggregation strategy [40] is used to compress the bundle. The algorithm generates three sequences of iterates as follows: {y k }, the sequence of intermediate points, at which the aggregate linear models of f are generated; {x k }, the sequence of trial points; {x k }, the sequence of stability centers.
We make some comments about Algorithm 1 as follows.
Call the oracle (6) at x 1 to compute the approximate values f 1 x and g 1 x . Choose an initial error bound ε 1 x ≥ 0 and a descent target . Let i 1 t := 0, l := 1, k(l) := 1 and k := 1.
Step 5 (Noise attenuation). If v k < − k , set t k := 10t k , i k t := k, and go back to Step 2.

Remark 1. (i)
Theoretically speaking, the model functionf k can be the simplest form max{f k−1 , f k }, but in order to keep numerical stability, it may additionally consist of some active linearizations.
(iii) Iff k is a polyhedral function, then subproblem (11) is equivalent to a convex quadratic programming and thus can be solved efficiently. In addition, if h is simple, subproblem (13) can also be solved easily, or even has a closed-form solution (say h(x) = 1 2 x 2 ). (iv) The role of Step 5 is to reduce the impact of inexactness. The algorithm loops between steps 2-5 by increasing the step size t k until v k ≥ −ε k .
(v) The stability center, descent target and error bound keep unchanged in the loop between Steps 2 and 5 (vi) In order to establish global convergence of the algorithm, the descent target and error bound at Step 6 should be suitably updated. Some detailed rules are presented in the next section.
The following lemma summarizes some fundamental properties of Algorithm 1, whose proof is a slight modification of that in [22], Lemma 2.2.  (12) and (14) satisfy The linearizationsf k ,h k ,F k satisfy the following inequalities (ii) The aggregate subgradient p k defined in (15) and the above linearizationF k can be expressed as follows (iii) The predicted descent v k and the aggregate linearization error k of (15) satisfy v k = t k p k 2 + k and k = F k (iv) The aggregate linearizationF k is also expressed (v) Denote the optimality measure by and Proof. (i) From the optimality condition of subproblem (11), we obtain which implies p k f ∈ ∂f k (y k+1 ). In addition, the fact thatf k (y k+1 ) =f k (y k+1 ) yieldsf k ≤f k . Similarly, by the optimality condition of (14), we have (ii) By (14), we obtain Utilizing the linearity ofF k (·), (12) and (19), we derivē (iii) We obtain directly v k = k + t k p k 2 from (15). Combining (15) and (ii), we have (v) Using the Cauchy-Schwarz inequality in the definition (22) gives (vi) By (iii), it is easy to get (25). Next, by (18), (20) and (9), we conclude that, if F k x .
Relation (26) follows from the facts that v k ≥ k and v k ≥ t k p k 2 /2. Relation (27) follows x shows that p k < ( 2ε k(l) x t k ) 1/2 , and therefore (28) holds. (17) shows that p k f is a subgradient of the model functionf k at y k+1 and that p k h is a subgradient of h at x k+1 . V k defined by (22) can be viewed as an optimality measure of the iterates, which will be proved to converge to zero in the next section. Relation (24) is also a test for optimality, in thatx k is an approximate optimal solution to problem (1) whenever V k is sufficiently small.

Global Convergence
This section aims to establish the global convergence of Algorithm 1 for various oracles. These oracles are controlled by two parameters: the error bound ε x and the descent target γ x . In Table 1, we present the choices of these two parameters for different type of instances described in Section 2, including Exact (Ex), Partially Inexact (PI), Inexact (IE), Asymptotically Exact (AE) and Partially Asymptotically Exact (PAE) oracles, where the constants are selected as θ, κ ∈ (0, 1), and κ ∈ (0, κ).
The following lemma is crucial to guarantee the global convergence of Algorithm 1.

Lemma 2.
The descent target is always reached at the stability centers, that is, F k x ≤ γ kx for all k ≥ 1.
Proof. For instances Ex, IE and AE, since γ k x = +∞, the claim holds immediately. For instances PI and PAE, we have γ k+1 Step 0 it follows that x . In addition, for k ≥ 2, since θ ∈ (0, 1), once the descent test (16) is satisfied at iteration k − 1, we have The following lemma shows that an (approximate) optimal solution can be obtained whenever the algorithm terminates finitely or loops infinitely between Steps 2 and 5.

Lemma 3.
If either Algorithm 1 terminates at the kth iteration at Step 4, or loops between Steps 2 and 5 infinitely, then (i)x k is an optimal solution to problem (1) for instances Ex and PI.
x , for instances AE and PAE.
Proof. Firstly, suppose that Algorithm 1 terminates at Step 4 with iteration k. Then from (23), we have V k = 0. This together with (24) shows that Thus, from (7), we can conclude that: x for instances AE and PAE. Secondly, suppose that Algorithm 1 loops between Steps 2 and 5 infinitely. Then from Lemma 2 and the condition at Step 5, it follows that (28) holds and t k ↑ ∞. Thus, we obtain V k → 0. This along with (24) implies (29), and therefore the claims hold by repeating the corresponding lines in first case.
From the above lemma, in what follows, we may assume that Algorithm 1 neither terminates finitely nor loops infinitely between Steps 2 and 5. In addition, as in Reference [22], we assume that the model subgradients p k f ∈ ∂f k (y k+1 ) in (17) satisfy that {p k f } is bounded if {y k } is bounded. Algorithm 1 must take only one of the following two cases: (i) the algorithm generates finitely many descent steps; (ii) the algorithm generates infinitely many descent steps.
We first consider case (i), in which two subcases may occur: t ∞ := lim k t k = ∞ and t ∞ < ∞. The first subcase of t ∞ = ∞ is analyzed in the following lemma.

Lemma 4.
Suppose that Algorithm 1 generates finitely many descent steps, that is, there exists an indexk such that only null steps occur for all k ≥k, and that t ∞ = ∞. Denote K := {k ≥k : t k+1 > t k }, then V k → 0 as k ∈ K, k → ∞.
Proof. For the last time t k increases before Step 5 for k ∈ K, one has which along with t k → ∞ shows the lemma.
The following lemma analyzes the second subcase of t ∞ < ∞.

Lemma 5.
Suppose that there existsk such thatx k =x¯k and t min ≤ t k+1 ≤ t k for all k ≥k. If the descent criterion (16) fails for all k ≥k, then V k → 0.
Proof. In view of the facts that t min ≤ t k+1 ≤ t k andx k =x¯k for all k ≥k, we know that only null steps occur and t k does not increase at Step 5. By Taylor's expansion, Cauchy-Schwarz inequality, and the properties of subproblems (11) and (13), we can conclude that v k → 0, so the conclusion holds from (27). For more details, one can refer to [11], Lemma 3.2.
By combining Lemmas 4 and 5, we have the following lemma.

Lemma 6.
Suppose that there existsk such that only null steps occur for all k ≥k.
Now, we can present the main convergence result for the case where the algorithm generates finitely many descent steps. Theorem 1. Suppose that Algorithm 1 generates finitely many descent steps, and thatx¯k is the last stability center. Then,x¯k is an optimal solution to problem (1) for instances Ex and PI; an ε-optimal solution for IE; and an ε¯k (l) x -optimal solution for AE and PAE.
Proof. Under the stated assumption, we know thatx k ≡x¯k and f k x ≡ f¯k x for all k ≥k. This together with (24) and Lemma 6 shows that Hence, similar to the proof of Lemma 3, we obtain the results of the theorem.
Next, we consider the second case where the algorithm generates infinitely many descent steps.

Lemma 7.
Suppose that Algorithm 1 generates infinitely many descent steps, and that F ∞ x := lim k F k x > −∞.
Proof. From the descent test condition (16), we may first prove that v k K → 0, and therefore k , t k p k 2 K → 0 from (26) and p k K → 0 from the fact that t k ≥ t min . It can be further proved that k , p k K → 0, so it follows lim k∈K V k = 0 from the definition of V k . Moreover, under the condition that {x k } is bounded, we have V k K → 0 by (v) of Lemma 1. For more details, one can refer to [11], Lemma 3.4.
Finally, we present the convergence results for the second case.

Theorem 2.
Suppose that Algorithm 1 generates infinitely many descent steps, F ∞ x > −∞, and that the index set K is defined in Lemma 7. Then Table 1; x for the remaining instances in Table 1; Proof. It is obvious that (i) For instance IE, it follows that ε k+1 x = ε and F k+1 x ≤ γ k+1 x = +∞ for all k ∈ K. Then from (7), This along with (30) shows part (i).
(ii) Next, the other four instances in Table 1 are considered separately. For instance Ex, we have ε k+1 For instance PI, we have ε k+1 This implies F(x k+1 ) = F k+1 x , and therefore For instance AE, we have ε k+1 This along with Lemma 7 (v k For instance PAE, we have ε k+1 Again from Lemma 7, we obtain Summarizing the above analysis and noticing (30), we complete the proof of part (ii).
(iii) From Lemma 7, we know that lim k∈K V k = 0. This together with (24) shows part (iii).

Numerical Experiments
In this section, we aim to test the numerical efficiency of the proposed algorithm. In the fields of production and transportation, finance and insurance, power industry, and telecommunications, decision makers usually need to solve problems with uncertain information. As an effective tool to solve such problems, stochastic programming (SP) has attracted more and more attention and research on its practical instances and theories; see, for example, References [41,42]. We consider a class of two-stage SP problems with fixed recourse, whose discretization of uncertainty into N scenarios has the form (see e.g., References [23,43]) where x is the first-stage decision variable, c ∈ R n 1 , A ∈ R m 1 ×n 1 , and b ∈ R m 1 . In addition, the recourse function is where corresponding to the ith scenario (h i , T i ), with probability p i > 0 for h i ∈ R m 2 and T i ∈ R m 2 ×n 1 .
Here π is the second-stage decision variable. Clearly, by introducing the indicator function δ X , problem (31) can be written as the form of (5), and then becomes the form of (1) by setting h(x) = δ X (x).
The above recourse function can be written as its dual form: where q ∈ R n 2 and W ∈ R m 2 ×n 2 . By solving these linear programming problems to return solutions with precision up to a given tolerance, one can establish an inexact oracle in the form (6), see Reference [23] for more detailed description. The instances of SP problems are downloaded from the link: http://pwp.gatech.edu/guanghuilan/computer-codes/.
Four instances are tested, namely, SSN(50), SSN(100), 20-term(50), 20-term(100), where the integers in the brackets mean the number of scenarios N. Here, the SSN instances come from the telecommunications and have been studied by Sen, Doverspike, and Cosares [44]. And the 20-term instances come from the motor freight carrier's problem and have been studied by Mak, Morton, and Wood [45]. The dimensions of these instances are listed in Table 2. The parameters are selected as: κ = 0.04, t min = 0.1 and t 1 = 1.1. The maximum bundle size is set to be 35. All the tests were performed in MATLAB (R2014a) on a PC with Intel(R) Core(TM) i7-4790 CPU 3.60GHz, 4GB RAM. The quadratic programming and linear programming subproblems were solved by the MOSEK 8 toolbox for MATLAB; see http://www.mosek.com.
We first compare our algorithm (denoted by GALBM) with the accelerated prox-level method (APL) in Reference [43], where the tolerances of the linear programming solver of MOSEK are set by default. The results are listed in Table 3, in which the number of iterations (NI), the consumed CPU time in seconds (Time), and the returned minimum values ( f * ) are compared. Note that we use the MATLAB commands tic and toc to measure the consumed CPU time. For each instance, we run 10 times and report the average CPU time. From Table 3, we see that, when similar solution quality is achieved, GALBM can significantly outperform than APL in terms of the number of iterations and CPU time.

Conclusions
In this paper, we have proposed a generalized alternating linearization bundle method for solving structured convex optimization with inexact first-order oracles. Our method can handle various inexact data by making use of the so-called on-demand accuracy oracles. At each iteration, two interrelated subproblems are solved alternately, aiming to reduce the computational cost. Global convergence of the algorithm is established under different types of inexactness. Numerical results show that the proposed algorithm is promising.
Author Contributions: C.T. mainly contributed to the algorithm design and convergence analysis; Y.L. and X.D. mainly contributed to the convergence analysis and numerical results; and B.H. mainly contributed to the numerical results. All authors have read and agree to the published version of the manuscript.
Funding: This work was supported by the National Natural Science Foundation (11761013) and Guangxi Natural Science Foundation (2018GXNSFFA281007) of China.