Noname manuscript No. (will be inserted by the editor) Optimal subgradient algorithms with application to large-scale linear inverse problems Masoud Ahookhosh Received: date / Accepted: date Abstract This study addresses some algorithms for solving structured unconstrained convex optimization problems using first-order information where the underlying function includes high-dimensional data. The primary aim is to develop an implementable algorithmic framework for solving problems with multiterm composite objective functions involving linear mappings using the optimal subgradient algorithm, OSGA, proposed by Neumaier in [49]. To this end, we propose some prox-functions for which the corresponding subproblem of OSGA is solved in a closed form. Considering various inverse problems arising in signal and image processing, machine learning, statistics, we report extensive numerical and comparisons with several state-of-the-art solvers proposing favourably performance of our algorithm. We also compare with the most widely used optimal first-order methods for some smooth and nonsmooth convex problems. Surprisingly, when some Nesterov-type optimal methods originally proposed for smooth problems are adapted for solving nonsmooth problems by simply passing a subgradient instead of the gradient, the results of these subgradient-based algorithms are competitive and totally interesting for solving nonsmooth problems. Finally, the OSGA software package is available. Keywords Structured convex optimization · Linear inverse problems · Nonsmooth optimization · Sparse optimization · First-order black-box oracle · Optimal complexity · Subgradient methods · High-dimensional data · Compressed sensing · Image restoration Mathematics Subject Classification (2000) 90C25 · 90C60 · 49M37 · 65K05 1 Introduction In many applications, e.g., those arising in signal and image processing, machine learning, compressed sensing, geophysics and statistics, key features cannot be studied by straightforward investigations, but must be indirectly inferred from some observable quantities. Thanks to this characteristic, they are usually referred as inverse problems. If a linear relevance between the features of interest and observed data can be prescribed, this leads to linear inverse problems. If y ∈ Rm be an indirect observation of an original object x ∈ Rn and A : Rn → Rm is a linear operator, then the linear inverse problem is defined using y = Ax + ν, (1) where ν ∈ Rm represents an additive noise about which little is known apart from which qualitative knowledge. Motivation & history. In practice, the systems (1) is typically underdetermined , rank-deficient or ill-conditioned . The primary difficulty with linear inverse problems is that the inverse object is extremely sensitive to y due to small or zero singular values of A meaning ill-conditioned systems behave like singular cases. Indeed, in the case that A−1 for square problems or pseudo-inverse A† = (A∗ A)−1 A∗ for full rank over-determined systems is exist, then analyzing the singular value decomposition shows that Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria E-mail: masoud.ahookhosh@univie.ac.at 2 Masoud Ahookhosh x̃ = A−1 y or x̃ = A† y be inaccurate and meaningless approximation for x, see [50]. Furthermore, when the vector δ is not known, one cannot solve (1) directly. From underdetermined or rank-deficient feature of inverse problems, we know that if a solution exists, then there exist infinitely many solutions. Hence some additional information is required to determine a satisfactory solution of (1). It is usually interested to determine a solution by minimizing the residual kAx − yk2 , leading to the well-known least-square problem min x∈Rn 1 kAx − yk22 , 2 (2) where k.k2 denotes l2 -norm. In view of ill-condition feature of the problem (1), the solution of (2) is usually improper. Hence Tikhonov in [53] proposed the following penalized minimization problem λ 1 kAx − yk22 + kxk22 , 2 2 min x∈Rn (3) where λ is the regularization parameter that controls the trade-off between the least square data fitting term and the regularization term. The reformulated problem is also smooth and convex, and selecting a suitable regularization parameter leads to a well-posed problem. In many applications, we seek the sparsest solution of (1) among all the solutions provided that y is acquired from a highly sparse observation. To do so, the following problem is proposed min kxk0 s.t. kAx − yk ≤ ǫ where kxk0 is the number of all nonzero elements of x and ǫ is a nonnegative constant. It is well-known that the objective function is nonconvex, however its convex relaxation is much more preferable. As investigated in [18], this is achievable by considering the problem min kxk1 s.t. kAx − yk ≤ ǫ. or its unconstrained version 1 kAx − yk22 + λkxk1 , (4) 2 which is referred as basis pursuit denoising or lasso. Due to non-differentiability of the l1 -norm, this problem is nonsmooth and its solving is more difficult than solving (3). In general, linear inverse problems can be modelled as a special case of the following affine composite minimization problem (5) minn f (Ax) + ϕ(Wx). min x∈Rn x∈R n Usually, f : R → R is a smooth convex loss function defined by a black-box oracle and ϕ : Rn → R is a regularizer or regularization function, which is usually nonsmooth and possibly nonconvex. This means that the general underlying regularizer can be smooth or nonsmooth and convex or nonconvex. In this paper, we consider structured convex optimization problems of the form (5) when regularizers can be smooth or nonsmooth convex terms, so nonconvex regularizers are not the matter of this study. Over the past few decades, because of emerging many applications of inverse problems in applied sciences and engineering, solving structured composite optimization problems of the form (5) attracts many interesting researches. To exploit different features of objective functions, a variety of techniques for solving this class of problems are introduced, where proposed optimization schemes employ particular structures of objective functions. It is also known that solving nonsmooth problems are usually more difficult than solving smooth problems, even for unconstrained cases. In the meantime, black-box convex optimization is theoretically and computationally regarded as a mature area in optimization and frequently employed to solve problems with large number of variables appearing in various applications. Since problems of the form (5) mostly include high-dimensional data, optimization schemes should avoid of exploiting costly operations and also high amount of memory. Therefore, the first-order black-box oracle is intensively used to construct appropriate schemes for handling large-scale problems, where just function and subgradient evaluations are needed. Some popular first-order references are gradient-based algorithms [25, 45, 46, 48], subgradient-based schemes [8, 12, 39, 42, 47], bundle-type procedures [35], smoothing techniques [9, 43] and proximal gradient methods [9, 11, 52]. Optimal subgradient algorithms with application to large-scale linear inverse problems 3 In the context of convex optimization, the analytical complexity refer to the number of calls of an oracle, for example iterations or function and subgradient evaluations, to reach an ε-solution of the optimum, for more information see [39, 42]. In 1983, Nemirovski & Yudin in [39] proved that the lower complexity bound of first-order methods for finding √ an ε-solution of the optimum for smooth convex functions with Lipschitz continuous gradient is O(1/ ε), while for Lipschitz continuous nonsmooth convex function it is O(1/ε2 ). Indeed, the interesting feature of these error bounds is that they are independent of the problem’s dimension. For the established class of problems, an algorithm called optimal if it can achieve these worst-case complexities. In a case of nonsmooth convex problems, the well-known subgradient decent and mirror descent methods are optimal, see [39]. The seminal optimal first-order method for smooth convex problems with Lipschitz continuous gradient is introduced by Nesterov [40] in 1983, and since then he has developed the idea of optimal methods for various classes of problems in [41, 42, 46, 48]. It is clear that the gap between the complexity of smooth and nonsmooth convex problems is considerably big. This fact motivates Nesterov to develop some schemes with better performance for nonsmooth convex problems so called smoothing techniques where the worst-case complexity of these methods is in order O(1/ε) and can be improved using the idea of optimal gradient methods, see [43, 44, 45]. The theoretical analysis and surprising computational results of optimal methods are totally interesting, especially for solving large-scale problems. These appealing features along with increasing the interest of solving large-scale problems attract many authors to develop the idea of the optimal firstorder methods, for examples Beck & Teboulle [8], Becker et al. [11], Devolder et al. [26], Gonzaga et al. [31, 32], Lan [35], Meng & Chen [37], Nemirovski & Nesterov [38] and Tseng [54]. In most of the studies mentioned, the first difficulty in implementation is the algorithms need the knowledge of global parameters such as the Lipschitz constant of objective function for nonsmooth problems and the Lipschitz constant for its gradient in smooth cases. Nesterov in [48] proposed an adaptive procedure to approximate the Lipschitz constant, however it still needs a suitable initial guess. Recently Lan in [35] based on employing bundle techniques avoids of requiring the global parameters. More recently, Neumaier in [49] introduced a novel optimal algorithm on the basis of incorporating a linear relaxation of the objective function and a prox-function into a fractional subproblem to construct a framework that can be employed for both smooth and nonsmooth convex optimization in the same time and to avoid of needing Lipschitz constants. It called OSGA, where it only needs the first-order blackbox oracle. In a case that no strong convexity is assumed, it shares the uniformity, freeness of global parameters and complexity properties of Lan’s approach, but has a far simpler structure and derivation. Thanks to achieving the optimal complexity of first-order methods, it can be regarded as a fully adaptive alternative of Nesterov-type optimal methods. Furthermore, low memory requirements of this algorithm make it appropriate for solving large-scale practical problems. Contribution. The main focus of this paper is to extend and develop the basic idea of the optimal subgradient algorithm for minimizing general composite functions to deal with structured convex optimization problems in applications. There exist many applications of linear inverse problems with a complicated objective function consisting of various linear operators and regularizers, which are not supported by algorithms and solvers using the proximal mapping in their structures. For example, consider the scaled lasso problem that is defined the same as (4) in which the regularization term kxk1 is substituted by kW xk1 for a general linear operator W . In this case, the following proximal mapping should be solved as an auxiliary subproblem 1 proxλ (y) = argmin kx − yk22 + kWxk1 . n 2 x∈R Indeed, this proximity problem cannot be solved explicitly using iterative shrinkage-thresholding except for the case that W is orthonormal, i.e., it cannot be solved using usual proximal-based algorithms. Thanks to accessible requirements of OSGA needing just function and subgradient evaluations, solving most of such complicated problems are tractable in reasonable costs. Furthermore, some appropriate prox-functions are proposed for unconstrained problems of the form (5) such that the corresponding subproblem of OSGA can be solved explicitly. We also adapt some Nesterov-type optimal methods originally proposed for smooth problems to solve nonsmooth problems by passing subgradient instead of gradient. The surprising results indicate that the adapted algorithms perform competitively or even better than optimal methods specifically proposed for solving nonsmooth problems. Finally, the paper is accompanied with a software release. 4 Masoud Ahookhosh Contents. The rest of the paper is organized as follows. In next section, we establish the basic idea of the optimal subgradient algorithm for the general affine composite functions. Section 3 devotes to implementing OSGA and comparing it with some popular first-order methods and state-of-the-art solvers for some practical problems in applications. The OSGA software package is described in Section 4. Finally, we conclude our results in Section 5. 2 Optimal subgradient framework In this section, we shall give a short sketch of the optimal subgradient framework introduced by Neumaier in [49] for composite functions in the presence of linear operators which is appropriate for solving considered applications of linear inverse problems. This paper considers composite functions of the form (5) or a more general form that will be discussed in the sequel, however we can still solve constrained problems by defining an indicator function. Let us to consider the following convex constrained problem min f (Ax), (6) x∈C where f : C → R is a convex function defined on a nonempty convex subset C of a finite-dimensional real vector space X with a bilinear pairing hh, zi, where h is in the dual space X ∗ which is formed by all linear functions on R. The indicator function of the set C is defined by  0 x ∈ C; IC (x) = (7) ∞ otherwise. Since C is convex set, it is straightforward to show IC is a convex function. Thus, by setting ϕ(x) = IC (x), we can reformulate the problem (6) in an equivalent unconstrained form (5), which is desired problem of the current paper. In some applications, one has no easy access to the matrices used in (5), instead it is assembled in a subroutine without ever forming explicit matrices so that it is preferred to employ linear operators in (5) instead of their matrix representations. We generally consider functions of the form Ψ (x) = n1 X fi (Ai x) + n2 X ϕj (Wj x), (8) j=1 i=1 where fi : Ui → R for i = 1, 2, · · · , n1 are convex smooth functions, ϕj : Vj → R for j = 1, 2, · · · , n2 are smooth or nonsmooth convex functions and Ai : X → Ui for i = 1, 2, · · · , n1 and Wj : X → Vj for j = 1, 2, · · · , n2 are linear operators. Here, Ψ is called the multi-term affine composite function. Since each affine term Ai x or Wj x and each function fi or ϕj for i = 1, · · · , n1 and j = 1, · · · , n2 is convex and the domain of Ψ defined by  ! n n1 2 \ \ \  dom ϕj (Wj x) dom fi (Ai x) dom Ψ = j=1 i=1 is a convex set, then the function Ψ is convex on its domain. Considering (8), we are generally interested in solving the following composite convex optimization problem Ψ ∗ := min Ψ (x) x∈X (9) Under considered properties of Ψ , the problem (9) has a global optimizer denoted by x∗ and Ψ ∗ = Ψ (x∗ ) is its global optimum. The functions of the form Ψ are frequently appeared in the recent interests of employing hybrid regularizations or mixed penalty functions for solving problems in the fields such as signal and image processing, machine learning and statistics. In the sequel, we will see that many well studied structured optimization problems are special cases of (8). It is also clear that (9) is a generalization of (5). For example, the following scaled elastic net problem Ψ (x) = λ2 1 kAx − bk22 + λ1 kW1 xk1 + kW2 xk22 2 2 (10) Optimal subgradient algorithms with application to large-scale linear inverse problems 5 can be considered in the form (8) by setting Ψ (x) = f1 (A1 x) + ϕ1 (W1 x) + ϕ2 (W2 x) where f1 (A1 x) = λ2 1 2 2 2 kA1 x − bk2 , ϕ1 (W1 x) = λ1 kW1 xk1 and ϕ2 (W2 x) = 2 kW2 xk2 . Basic idea of optimal subgradient scheme. This section summarizes the basic idea of OSGA [49] which is adapted for the considered unconstrained problem (9). Let us to consider an arbitrary convex function Ψ : X → R of the form (8). A vector gΨ (x) ∈ X ∗ is called a subgradient of Ψ at a given point x ∈ dom Ψ if for any z ∈ dom Ψ we have that Ψ (z) ≥ Ψ (x) + hgΨ (x), z − xi. The set of all subgradients of the function Ψ at x is called subdifferential of Ψ at x denoted by ∂Ψ (x). The subdifferential of a general convex function Ψ at any point x in its domain is nonempty, and it is unique if Ψ is differentiable. We now consider a minimization problem of the form (8) and suppose that the first-order black-box oracle O(x) = (Ψ (x), gΨ (x)), is available, where gΨ (x) ∈ ∂Ψ (x) at x ∈ X . The primary objective of the optimal subgradient algorithm (OSGA) is to construct a sequence of iterations such that the error bound Ψ (xb ) − Ψ ∗ is monotonically reducing to be smaller than an accuracy parameter ε, i.e. 0 ≤ Ψb − Ψ ∗ ≤ ε, where Ψb := Ψ (xb ) is the function value in the best known point xb . To derive a framework for OSGA, let us to define the following affine global underestimator , which is a linear relaxation of Ψ , Ψ (z) ≥ γ + hh, zi, (11) for all z ∈ dom Ψ with γ ∈ R and h ∈ X ∗ , which is clearly available due to the first-order condition for the convex function Ψ . We also consider a continuously differentiable prox-function Q : X → R, which is a strongly convex function with the convexity parameter set to σ satisfying σ Q(λx + (1 − λ)z) + λ(1 − λ) kx − zk2 ≤ λQ(x) + (1 − λ)Q(z), 2 (12) for all x, z ∈ X and λ ∈ (0, 1) where k.k is an arbitrary norm of the vector space X . This is equivalent to Q(z) ≥ Q(x) + hgQ (x), z − xi + σ kz − xk2 , 2 (13) for all x, z ∈ X where gQ (x) denotes the gradient of Q at x. Furthermore, we suppose that Q0 := inf Q(z) > 0. z∈X (14) The point z0 ∈ X minimizing the problem (14) is called the center of Q. At the end of this section, we will propose some prox-functions satisfying these conditions. Afterwards, the following minimization subproblem is defined γ + hh, zi E(γ, h) := − inf . (15) z∈ dom Ψ Q(z) By the definition of prox-functions, it is clear that the denominator Q(z) is positive and bounded away from zero on the feasible region dom Ψ . The solution of this subproblem (15) is denoted by u := U (γ, h), and it is supposed that e := E(γ, h) and u can be effectively computed. In particular, for the unconstrained problem (8), we will show that this subproblem can be solved explicitly and cheaply for some suitable prox-functions. Theoretical analysis. We here establish a bound on error Ψ (xb ) − Ψ ∗ for the proposed scheme. The definition of e = E(γ, h) in (15), for arbitrary γb ∈ R and h ∈ X ∗ , simply implies that γb + hh, zi ≥ −E(γb , h)Q(z), for all z ∈ X . From (11) and setting γb = γ − Ψ (xb ), we have Ψ (z) − Ψ (xb ) ≥ γb + hh, zi, (16) 6 Masoud Ahookhosh for all z ∈ X . By this fact, assuming E(γb , h) ≤ η, setting z = x∗ and (16), we conclude the following upper bound on the function value error 0 ≤ Ψ (xb ) − Ψ ∗ ≤ ηQ(x∗ ). (17) It is clear that if an upper bound for Q(x∗ ) is known or assumed, this bound translates into a computable error estimate for the minimal function value. But even in the absence of such an upper bound, the optimization problem (8) can be solved with a target accuracy 0 ≤ Ψ (xb ) − Ψ ∗ ≤ εQ(x∗ ) (18) if one manages to decrease the error factor η from its initial value until η ≤ ε, for some target accuracy tolerance ε. This bound clearly means that the rate of decrease in the error bound Ψ (xb ) − Ψ ∗ , is at least the same as convergence rate of the sequence η. Indeed, this fact is the main motivation of OSGA. Neumaier in [49] considers a problem of the form (6) and derive the following complexity bounds for OSGA matching to the optimal complexity bounds of first-order methods for smooth and nonsmooth problems, see [39, 42]. Theorem 1 Suppose that f is a convex function of the form (6) and let {xk } is generated by OSGA. Then complexity bounds for smooth and nonsmooth problems are as follows: (i) (Nonsmooth complexity bound) If the point generated by OSGA stay in bounded region of the
interior of C, or f is Lipschitz continuous in C, then
the total number of iterations needed is O 1/ε2 . Thus the asymptotic worst case complexity is O 1/ε2 . (ii) (Smooth complexity bound) √ If f has Lipschitz continuous gradient, the total number of iterations needed for the algorithm is O (1/ ε). Algorithmic framework. Before presenting the OSGA algorithm, we first look at the first-order oracle of the problem (8) which is important in the presence of linear mappings, and then we outline the adapted algorithm. On the one hand, practical problems are commonly involving high-dimensional data, i.e. function and subgradient evaluations are quite expensive. On the other hand, in many applications the most computational cost of function and subgradient evaluations relates to applying direct and adjoint linear operators, originated from an existence of affine terms in (8). This facts suggest that we should apply linear operators and their adjoints as less as possible during our implementations. For the sake of simplicity in the rest of the paper, we employ Ψx and gx instead of Ψ (x) and gΨ (x), respectively. Considering the structure of (8), we see that affine terms Ai x for i = 1, 2, · · · , n1 and Wj x for j = 1, 2, · · · , n2 are appearing in both function and subgradient evaluations. As a result, in the oracle, we set vxi = Ai x for i = 1, 2, · · · , n1 and wxj = Wj x for j = 1, 2, · · · , n2 . Therefore, the function value and subgradient of Ψ at x are computed in the following scheme: Algorithm 1: NFO-FG nonsmooth first-order oracle Input: global parameters: Ai for i = 1, · · · , n1 ; Wj for j = 1, · · · , n2 ; Output: Ψx ; gΨ (x); begin vxi ← Ax for i = 1, · · · , n1 ; wxj ← P Wx for j = 1, ·P · · , n2 ; n1 n2 Ψx ← i=1 fi (vxi ) + j=1 ϕj (wxj ); P n1 ∗ P n2 i gΨ (x) ← i=1 Ai ∂fi (vx ) + j=1 Wi∗ ∂ϕj (wxj ); end local parameters: x; The structure of NFO-FG immediately implies that any call of the oracle O(x) = (Ψx , gΨ (x)) requires n1 + n2 calls of each direct and adjoint linear operator. It is also clear that by this scheme, one can avoid of double applying of expensive linear operators in the computation of function values and subgradients. In cases that algorithms need only a function value or a subgradient, we consider two special cases of NFO-FG that just return a function value or a subgradient, which are respectively called NFO-F and NFO-G. We also emphasize that in the cases that total computational cost of the oracle dominated by applying linear mappings, the complexity of an algorithm can be measured by counting the number of direct and adjoint linear operators used to achieve a prescribed accuracy ε. Optimal subgradient algorithms with application to large-scale linear inverse problems 7 Considering what established in this section and those introduced in Section 2 of [49] about linear relaxation constructions and the case of strongly convex functions, we outline a version of OSGA for multi-term affine composite function (8) as follows: Algorithm 2: OSGA for multi-term affine composite functions Input: global tuning parameters: δ; αmax ∈ ]0, 1[; 0 < κ′ ≤ κ; local parameters: x0 ; µ ≥ 0; Ψtarget ; Output: xb ; Ψxb ; begin choose an initial best point xb ; compute Ψxb and gΨ (xb ) using NFO-FG; if Ψxb ≤ Ψtarget then stop; else h ← gxb − µgQ (xb ); γ ← Ψxb − µQ(xb ) − hh, xb i; γb ← γ − Ψxb ; u ← U (γb , h); η ← E(γb , h) − µ; end α ← αmax ; while stopping criteria do not hold do x ← xb + α(u − xb ); compute Ψx and gΨ (x) using NFO-FG; g ← gx − µgQ (x); h ← h + α(g − h); γ ← γ + α(Ψx − µQ(x) − hg, xi − γ); x′b ← argminz∈{xb ,x} Ψ (z, vz ); Ψx′b ← min{Ψxb , Ψx }; γb′ ← γ − fx′b ; u′ ← U (γb′ , h); x′ ← xb + α(u′ − xb ); compute Ψx′ using NFO-F; choose xb in such a way that Ψxb ≤ min{Ψx′b , Ψx′ }; γ b ← γ − Ψxb ; u ← U (γ b , h); η ← E(γ b , h) − µ; x b ← x b ; Ψ x b = Ψx b ; if Ψxb ≤ Ψtarget then stop; else update the parameters α, h, γ, η and u using PUS; end end end Pay attention to the structure of problem (8) and the OSGA framework, we only require to calculate the following simple operations to implement the algorithm. • • • • x + βy, where x, y ∈ Ui for i = 1, · · · , n1 and β ∈ R; r + βs, where r, s ∈ Vj for j = 1, · · · , n2 and β ∈ R; Ai x for i = 1, · · · , n1 and Wj x for j = 1, · · · , n2 where x ∈ X ; A∗i y with y ∈ Ui for i = 1, · · · , n1 and Wj∗ y with y ∈ Vj for j = 1, · · · , n2 . As a result of (18), a natural stopping criterion for the algorithm is the condition η ≤ ǫ, but our experiments shows that satisfying this condition is slow and takes many iterations. Then more sophisticated stopping criteria for the algorithm is needed. For example a bound on the number of iterations, function evaluations, subgradient evaluations or the running time or a combination of them can be used. As discussed in [49], OSGA uses the following scheme for updating the given parameters α, h, γ, η and u: 8 Masoud Ahookhosh Algorithm 3: PUS parameters updating scheme Input: global tuning parameters: δ; αmax ∈ ]0, 1[; 0 < κ′ ≤ κ; Output: α; h; γ; η; u; begin R ← (η − η)/(δαη); if R < 1 then h ← h; else ′ α ← min(αeκ (R−1) , αmax ); end α ← α; if η < η then h ← h; γ ← γ; η ← η; u ← u; end end local parameters: α; η; h̄; γ̄; η̄; ū; We here emphasize that the composite version of optimal subgradient algorithms inherits the main characteristics of the original OSGA framework. Thus, we expect the same theoretical analysis to be the same as that reported for OSGA in [49] which means that the current version remains optimal. Prox-functions and solving the auxiliary subproblem. In the rest of this section, we first propose some prox-functions and then derive a closed form for the corresponding subproblem (15). As mentioned, the prox-function Q(z) should be a strongly convex function that is analytically known. This means that it takes its unique global minimum which is positive by the definition of Q. Suppose z0 is the global minimizer of Q. By the definition of center of Q and the first-order optimality condition for (14), we have that gQ (z0 ) = 0. This fact along with (13) imply Q(z) ≥ Q0 + σ kz − z0 k2 , 2 (19) which simply means that Q(z) > 0 for all z ∈ X . In addition, it is interested that prox-functions be separable. Taking this fact into account and (19), appropriate choices of prox-functions for unconstrained problem (8) can be σ (20) Q1 (z) := Q0 + kz − z0 k22 2 and Q2 (z) := Q0 + n σX wi (zi − (z0 )i )2 , 2 i=1 (21) where Q0 > 0 and wi ≥ 1 for i = 1, 2, · · · , n. It is not hard to show that both (20) and (21) are strongly convex functions satisfying (19). To introduce a more general class of prox-functions, let us to define the quadratic norm on vector space X using p kzkD := hDz, zi by means of a preconditioner D, where D is symmetric and positive definite. The associated dual norm on X ∗ is then given by p khk∗D := kD−1 hkD = hh, D−1 hi, where D−1 denotes an inverse of D. Given a symmetric and positive definite preconditioner D, then it is natural to consider the quadratic function Q(z) := Q0 + σ kz − z0 k2D , 2 (22) Optimal subgradient algorithms with application to large-scale linear inverse problems 9 where Q0 is a positive number and z0 ∈ X is the center of Q. Now, we should prove that Q is a proxfunction and satisfies (19). The fact that gQ (x) = σD(x − z0 ) leads to Q(x) + hgQ (x), z − xi + σ σ kz − xk2D = Q(x) + hgQ (x), z − xi + hD(z − x), z − xi 2 2 σ σ = Q0 + hD(x − z0 ), x − z0 i + hgQ (x), z − xi + hD(z − x), z − xi 2 2 σ σ = Q0 + hD(x − z0 ), z − z0 i + hD(z − z0 ), z − xi 2 2 σ σ = Q0 + hD(x − z0 ), z − z0 i + hD(z − x), z − z0 i 2 2 σ = Q0 + hD(z − z0 ), z − z0 i 2 = Q(z). This clearly means that Q is strongly convex by the convexity parameter σ, i.e., (22) is a prox-function. Moreover, z0 is the center of Q and gQ (z0 ) = 0. This fact together with (13) and (14) imply that (19) holds. At this point, we emphasize that an efficient solving of the subproblem (15) is highly related to the selection of prox-functions. While using a suitable prox-function allows a very efficient way of computing u, employing other selections may significantly slow down the process of solving the auxiliary subproblem. Here, it is assumed that by appropriate selections of prox-functions the subproblem (15) can be solved explicitly. In the following result, we verify the solution of (15) using the prox-function (22). Proposition. 21 Suppose Q is determined by (22) and Q0 > 0. Then Q is a prox-function with the center z0 and satisfies (19). Moreover, the subproblem (15) corresponds to this Q is explicitly solved as follows u = z0 − E(γ, h)−1 σ −1 D−1 h (23) with p β12 + 4Q0 β2 2β p 2 , = E(γ, h) = 2Q0 β1 + β12 + 4Q0 β2
where β1 = γ + hh, z0 i and β2 = 21 σ −2 − σ −1 khk2∗ . −β1 + (24) Proof We already proved that Q is a prox-function and (19) holds. The rest of the proof is similar to that discussed in [49] for σ = 1. ⊔ ⊓ It is obvious that (20) and (21) are special cases of (22), so the solution (23) and (24) derived for (22) can be simply adapted for the situation (20) or (21) used as a prox-function. Furthermore, notice that the error bound (17) is proportional to Q(x∗ ) implying that an acceptable choosing of x0 make the term Q(x∗ ) = Q0 + 12 kx∗ − x0 k2D enough small. Hence, selecting a suitable starting point x0 close to the optimizer x∗ as much as possible has a positive effect on giving a better complexity bound. This also propose that a reasonable choice for Q0 can be Q0 ≈ 12 kx∗ − x0 k2 . This fact will be considered along numerical experiments. We conclude this section by considering cases that the variable domain is a set of all m × n matrices, X = Rm×n . In this case, we introduce a suitable prox-function such that the subproblem (15) is solved in a closed form. Details are described in the following proposition. Proposition. 22 Suppose Q0 > 0 and Q is determined by σ Q(Z) := Q0 + kZ − Z0 k2F , 2 (25) where k.kF is the Frobenius norm. Then Q is a prox-function with the center Z0 and satisfies (19). Moreover, the subproblem (15) corresponded to this Q is explicitly solved by U = Z0 − E(γ, H)−1 σ −1 H with p ξ12 + 4Q0 ξ2 2ξ p 2 = , 2Q0 ξ1 + ξ12 + 4Q0 ξ2
where ξ1 = γ + T r(H T Z0 ) and ξ2 = 12 σ −2 − σ −1 kHk2F . E(γ, H) = −ξ1 + (26) (27) 10 Masoud Ahookhosh Proof By the features of the Frobenius norm, it is clear that Z = Z0 is the minimizer of (14). Also, Q is continuously differentiable gQ (x) = σ(X − Z0 ). Thus, we have σ σ σ Q(X) + hgQ (X), Z − Xi + kZ − Xk2F = Q0 + hX − Z0 , X − Z0 i + σhX − Z0 , Z − Xi + kZ − Xk2F 2 2 2 σ σ = Q0 + h(X − Z0 ), Z − Z0 i + h(Z − Z0 ), Z − Xi 2 2 σ = Q0 + h(Z − Z0 ), Z − Z0 i 2 = Q(Z). This means that Q is strongly convex by the convexity parameter σ, i.e., (22) is a prox-function. Also by (13) and (14), the condition (19) is clearly hold. Now, we consider the subproblem (15) and drive its minimizer. Let us to define the function Eγ,H : X → R using γ + hH, Zi , Eγ,H (Z) := − Q(Z) where it is continuously differentiable and attains its supremum at Z = U . Therefore, in the point Z = U , we obtain E(γ, H)Q(U ) = −γ − hH, U i. (28) It follows from gQ (Z) = σ(Z − Z0 ), (28) and the first-order necessary optimality condition that E(γ, H)σ(U − Z0 ) + H = 0 or equivalently U = Z0 − E(γ, H)−1 σ −1 H. Setting e = E(γ, H) and substituting it into (28) straightforwardly lead to   1 −1 −1 2 eQ(U ) + γ + hH, U i = e Q0 + ke σ HkF + γ + hH, Z0 − e−1 σ −1 Hi 2   1 −2 −2 2 = e Q0 + e σ kHkF + γ + T r(H T Z0 ) + e−1 σ −1 kHk2F 2  
1 −2 2 T −1 = Q0 e + γ + T r(H Z0 ) e + σ −σ kHk2F 2 = Q0 e2 + ξ1 e + ξ2 = 0,
where ξ1 = γ + T r(H T Z0 ) and ξ2 = 12 σ −2 − σ −1 kHk2F . Solving this quadratic equation and selecting the larger solution lead to p −ξ1 + ξ12 + 4Q0 ξ2 2ξ p 2 = . E(γ, H) = 2Q0 ξ1 + ξ12 + 4Q0 ξ2 This completes the proof. ⊔ ⊓ 3 Numerical experiments This section reports extensive numerical experiments of the proposed algorithm dealing with some famous inverse problems in applications compared with some popular algorithms and software packages. We assess the convergence speed of OSGA and some widely used first-order methods for solving such practical problems. The experiment is divided into two parts. We first consider some imaging problems and give results of experiments with OSGA and some state-of-the-art solvers. In the second part, we adapt some smooth optimization algorithms for nonsmooth problems and report encouraging results and comparisons. All codes are written in MATLAB, and default values of parameters for the algorithms and packages are used. For OSGA, we used the prox-function (20) with Q0 = 21 kx0 k2 + ǫ, where ǫ is the machine precision in MATLAB. We also set the following parameters for OSGA: δ = 0.9; αmax = 0.7; κ = κ′ = 0.5; Ψtarget = −∞. All implementations are executed on a Toshiba Satellite Pro L750-176 laptop with Intel Core i7-2670QM Processor and 8 GB RAM. Optimal subgradient algorithms with application to large-scale linear inverse problems 11 3.1 Image restoration Image reconstruction, also called image restoration, is one of the classical linear inverse problems dating back to the 1960s, see for example [7, 14]. The goal is to reconstruct images from some observations. Let y ∈ U be a noisy indirect observation of an image x ∈ X with an unknown noise δ ∈ U . Suppose that A : X → U is a linear operator. To recover the unknown vector x, we use the linear inverse model (1). As discussed, according to the ill-conditioned or singular nature of the problem, some regularizations are needed. Recall (3), (4) or the more general model (8) to derive a suitable solution for the corresponding inverse problem and refer to [33, 34, 50] for more regularization options. The strategy employed for solving such penalized optimization problems depends on features of objective functions and regularization terms like differentiability or nondifferentiability and convexity or nonconvexity. The quality of the reconstructed image highly depends on these features. A typical model for image reconstruction consists of a smooth objective function penalized by smooth or nonsmooth regularization penalties. These penalties are selected based on expected features of the recovered image. Suppose that the known image x ∈ X is represented by x = Wτ in some basis domain, e.g., wavelet, curvelet, or time-frequency domains, where the operator W : U → V can be orthogonal or any dictionary. There are two main strategies for restoring the image x from an indirect observation y, namely synthesis and analysis. To recover a clean image, the synthesis approach reconstructs the image by solving the following minimization problem 1 kAWτ − yk22 + λϕ(τ ), (29) min τ ∈U 2 where ϕ is a convex regularizer. If τ ∗ is an optimizer of (29), then an approximation of the clean image x will be reconstructed using x∗ = Wτ ∗ . Alternatively, the analysis strategy recovers the image by solving min x∈X 1 kA(x) − yk22 + λϕ(x). 2 (30) Similarity of the problems (29) and (30) suggests that they are strongly related. These two approaches are equivalent in some assumptions, while variant in others, for more details see [29]. Although various regularizers like lp -norm are popular in image restoration, it is arguable that the best known and frequently employed regularizer in analysis approaches for image restoration is total variation (TV), which will be discussed later in this section. In the sequel, we report our numerical experiments based on the analysis strategy. The pioneering work on total variation as a regularizer for image denoising and restoration was proposed by Rudin, Osher, and Fatemi in [51]. Total variation regularizations are able to restore discontinuities of images and recover the edges. This makes TV appealing and widely used in applications. TV is originally defined in the infinite-dimensional Hilbert space, however for the digital image processing it is interested to consider it in a finite-dimensional space like Rn of pixel values on a two-dimensional lattice. In practice, some discrete versions of TV are favoured to be used. Two standard choices of discrete TV-based regularizers, namely isotropic total variation and anisotropic total variation, are popular in signal and image processing, defined by Pm−1 Pn−1 p (Xi+1,j − Xi,j )2 + (Xi,j+1 − Xi,j )2 kXkIT V = i Pm−1 j Pn−1 + i |Xi+1,n − Xi,n | + i |Xm,j+1 − Xm,j |, and Pm−1 Pn−1 kXkAT V = i {|Xi+1,j − Xi,j | + |Xi,j+1 − Xi,j |} Pm−1 j Pn−1 |Xi+1,n − Xi,n | + i |Xm,j+1 − Xm,j |, + i for X ∈ Rm×n , respectively. The image restoration problem with the discrete TV-based regularization can be formulated by 1 kA(X) − Y k2F + λ ϕ(X), (31) min X∈Rm×n 2 √ where A : Rm×n → Rr×m is a linear operator, ϕ = k.kIT V or ϕ = k.kAT V and kXkF = Tr XX T is the Frobenius norm in which Tr stands for the trace of a matrix. In what follows, we consider some prominent classes of imaging problems, more specifically denoising, inpainting and deblurring, and conduct our tests on numerous images. 12 Masoud Ahookhosh 3.1.1 Denoising with total variation Image denoising is one of the most fundamental image restoration problems aiming to remove noise from images, in which the noise comes from some resources like sensor imperfection, poor illumination, or communication errors. While this task itself is a prominent imaging problem, it is important as a component in other processes like deblurring with the proximal mapping. The main objective is to denoise images while the edges are preserved. To do so, various models and ideas have been proposed. Among all of the regularizers, total variation proposed a model preserving discontinuities (edges). The particular case of (31) when A is the identity operator, is called denoising and is very popular with total variation regularizers. Although the original total variation introduced in [51] is an infinitedimensional continuous constrained problem, we here consider the following unconstrained discrete finitedimensional version 1 (32) XY,λ = argmin kX − Y k2F + λ ϕ(X), X∈Rm×n 2 where ϕ = k.kIT V or ϕ = k.kAT V . It is clear that both k.kIT V and k.kAT V are nonsmooth semi-norms and the proximal operator (32) cannot be solved explicitly for a given Y . Rather than solving the problem (32) explicitly, it has been approximately solved by iterative processes similar to that proposed by Chambolle in [23]. We note that k.kIT V and k.kAT V are nonsmooth and convex, and their subdifferentials are available. Therefore, OSGA can be employed for solving (32). We consider the isotropic version of problem (32) where Y denotes a noisy image and XY,λ is denoised approximation of the original X corresponding to a regularization parameter λ. We run IST, TwIST [15] and FISTA [9, 10] in order to compare their performances with OSGA’s results. For IST, TwIST and FISTA, we use original codes and only adding more stopping criteria are needed in our comparisons. We set λ = 0.05 for all algorithms. As IST, TwIST and FISTA need to solve their subproblems iteratively, we consider three different cases by limiting the number of internal iterations to chit = 5, 10 and 20 iterations. Let us consider the recovery of the 1024 × 1024 Pirate image by IST, TwIST, FISTA and OSGA to verify the results visually and in more details. The noisy image is generated by using the MATLAB internal function awgn with SNR = 15, and the algorithms are stopped after 50 iterations. The results are summarized in Figures 1 and 2. To see details of this experiment, we compare the results of implementations in different measures in Figure 2. The performances are measured by some relative errors for points and function values and also so-called signal-to-noise improvement (ISNR). Here, ISNR and PSNR are defined by    √  mn kY − X0 kF and PSNR = 20 log10 , (33) ISNR = 20 log10 kX − X0 kF kX − X0 kF where X0 denotes the m × n clean image, Y is the observed image and pixel values are in [0, 1]. Generally, these ratios measure the quality of the restored image X relative to the blurred or noisy observation Y . From Figure 2, it is observed that OSGA outperforms IST, TwIST and FISTA with respect to all considered measures. Figure 2, shows that the denoised image looks good for all considered algorithms, where the best function value and the second best PSNR are attained by OSGA requiring 25.46 seconds to be stopped. However, it can be seen that by increasing the number of Chambolle’s iterations, the running time is dramatically increased for IST, TwIST and FISTA, however, a better function value or PSNR (compared with OSGA) is not attained. 3.1.2 Inpainting with isotropic total variation Image inpainting is a set of techniques to fill lost or damaged parts of an image such that these modifications are undetectable for viewers who are not aware of the original image. The basic idea is originated from reconstructing of degraded or damaged paintings conventionally carried out by expert artists. This is without no doubt very time consuming and risky. Thus, the goal is to produce a modified image in which the inpainted region is unified into the image and looks normal for ordinary observers. The pioneering work on digital image inpainting was proposed by Bertalmi et al. in [13] based on nonlinear PDE. Since then lots of applications of inpainting in image interpolation, photo restoration, zooming, super-resolution images, image compression and image transmission have emerged. In 2006, Chan in [24] proposed an efficient method to recover piecewise constant or smooth images by combining total variation regularizers and wavelet representations. Motivated by this work, we consider Optimal subgradient algorithms with application to large-scale linear inverse problems 13 0.16 0.16 0.16 IST TwIST FISTA OSGA 0.14 IST TwIST FISTA OSGA 0.14 0.12 0.12 0.12 0.1 0.1 0.08 0.08 0.06 0.06 0.06 0.04 0.04 0.04 0.02 0.02 0.02 0 0 0 5 10 15 20 25 30 35 40 45 50 (a) rel. 1 vs. iterations (chit = 5) IST TwIST FISTA OSGA 0.14 0.1 0.08 0 0 5 10 15 20 25 30 35 40 45 50 (b) rel. 1 vs. iterations (chit = 10) 0 5 10 15 20 25 30 35 40 45 50 (c) rel. 1 vs. iterations (chit = 20) 0 10 0 10 0 10 IST TwIST FISTA OSGA −1 10 IST TwIST FISTA OSGA −1 10 10 −2 −2 10 −2 10 −3 10 −4 10 10 −3 10 −3 10 −4 10 −4 10 −5 10 IST TwIST FISTA OSGA −1 −5 −5 0 5 10 15 20 25 30 35 40 45 50 (d) rel. 2 vs. time (chit = 5) 10 10 0 20 40 60 80 100 120 0 140 (e) rel. 2 vs. time (chit = 10) 50 100 150 200 250 (f) rel. 2 vs. time (chit = 20) 0 0 10 0 10 10 IST TwIST FISTA OSGA −1 10 −1 10 10 −2 10 −3 10 −4 10 10 −3 10 −3 10 −4 10 −4 10 −5 −5 −5 0 5 10 15 20 25 30 35 40 45 50 (g) rel. 2 vs. iterations (chit = 5) 10 10 0 5 10 15 20 25 30 35 40 45 0 50 (h) rel. 2 vs. iterations (chit = 10) 7 7 6 6 5 5 4 4 4 3 3 3 2 2 IST TwIST FISTA OSGA 0 5 10 15 20 25 30 35 40 45 (j) ISN R vs. iterations (chit = 5) 0 15 20 25 30 35 40 45 50 7 6 5 2 IST TwIST FISTA OSGA 1 50 10 8 8 1 5 (i) rel. 2 vs. iterations (chit = 20) 8 0 IST TwIST FISTA OSGA −1 −2 −2 10 10 IST TwIST FISTA OSGA 0 5 10 15 20 25 30 35 40 45 IST TwIST FISTA OSGA 1 50 (k) ISN R vs. iterations (chit = 10) 0 0 5 10 15 20 25 30 35 40 45 50 (l) ISN R vs. iterations (chit = 20) Fig. 1: A comparison among IST, TwIST, FISTA and OSGA for denoisinging the 1024 × 1024 Pirate image when they stopped after 50 iterations. Each column stands for comparisons with a different number of Chambolle’s iterations. While the first row denotes the relative error rel. 1 := kxk − x∗ k2 /kx∗ k2 of points versus iterations, the second row shows the relative error rel. 2 := (fk − f ∗ )/(f0 − f ∗ ) of function values versus time. The third row illustrates rel. 2 versus iterations and fourth row stands for ISNR (33) versus iterations. 14 Masoud Ahookhosh (a) Original image (b) Noisy image (c) OSGA: F = 3721.52, PSNR = 30.12, T = 25.46 (d) IST (chit = 5): F = 3791.57, PSNR = 30.05, T = 22.76 (e) IST (chit = 10): F = 3750.28, PSNR = 30.16, T = 43.05 (f) IST (chit = 20): F = 3741.17, PSNR = 30.08, T = 71.42 (g) TwIST (chit = 5): F = 3791.57, PSNR = 30.05, T = 47.02 (h) TwIST (chit = 10): F = 3729.11, PSNR = 30.06, T = 124.25 (i) TwIST (chit = 20): F = 3735.04, PSNR = 30.00, T = 224.44 (j) FISTA (chit = 5): F = 3791.57, PSNR = 30.05, T = 20.48 (k) FISTA (chit = 10): F = 3750.28, PSNR = 30.16, T = 42.04 (l) FISTA (chit = 20): F = 3741.17, PSNR = 30.08, T = 65.88; Fig. 2: Denoising of the 1024 × 1024 Pirate image by IST, TwIST, FISTA and OSGA. Here, F denotes the best function value, PSNR is defined by (33) and T stands for CPU time of the algorithms after 50 iterations. Optimal subgradient algorithms with application to large-scale linear inverse problems 15 inpainting with isotropic total variation to reconstruct images with missing data. We consider the 512×512 Head CT image with 40% missing sample. The same as previous section, IST, TwIST, FISTA and OSGA are employed to recover the image. The parameter λ sets to 0.09. We draw your attention to the different running time of these algorithms in each step, we stop them after 50 seconds. To have a fair comparison, three versions of each IST, TwIST and FISTA corresponding to 5, 10 and 20 iterations of Chambolle’s algorithm for solving the proximal subproblems are considered. The results of these implementations are illustrated in Figures 3 and 4. Similar to denoising, Figure 3 reports comparisons with some relative errors for step points and function values and also ISNR. These results display that IST is not comparable with others, and also OSGA is superior to IST, TwIST and FISTA regarding all considered measures. Since the data is not scaled in the implementation, PSNR is computed by   √ 255 mn , PSNR = 20 log10 kX − X0 kF where X0 is the clean image. The definition clearly proposes that a bigger PSNR value means a smaller kX −X0 kF , suggesting a better quality for the restored image. Conversely, a smaller PSNR value indicates a bigger kX −X0 kF implying the worse image quality. Figure 5 shows that IST cannot recover the missing data effectively compared with other algorithms. Moreover, the results suggest that OSGA achieves the best function value 187564.02 and the best PSNR 33.13 among the solvers considered. 3.1.3 Deblurring with isotropic total variation Image blur is one of the most common problem which happens in the photography and can often ruin the photograph. In digital photography, the motion blur is caused by camera shakes, which is difficult to avoid in many situations. Deblurring has been a fundamental research problem in the digital imaging and can be considered as an inverse problem of the form (1). The same as other inverse problems, this is an inherently ill-conditioned problem leading to the optimization problem of the form (5) or (8). Deblurring is generally much harder than denoising because most of algorithms that have been developed for solving this problem need to solve a subproblem of the form (32) in each iteration. Here, the deblurring problem of the form (31) with the isotropic total variation regularizer is considered. As discussed about denoising and inpainting, IST cannot recover images efficiently and is not comparable with the other algorithms considered. Hence, we instead take advantage of the more sophisticated package SpaRSA [55] in our implementations, so the comparison for deblurring involves TwIST, SpaRSA, FISTA and OSGA. Here, the regularization parameter set to λ = 0.05. We now consider deblurring of the 512 × 512 blurred/noisy Elaine image. Let Y be a blurred/noisy version of this image generated by 9 × 9 uniform blur and adding Gaussian noise with SNR = 40 dB. We recover it using (31) by the algorithms in 25 seconds of the running time. The results are demonstrated in Figures 5 and 6. In Figure 5, the algorithms are compared regarding relative errors for step points and function values and also ISNR. In Figures 5, subfigures (a)-(b) indicate that OSGA obtains a better level of accuracy regarding both of these measures, while Figure 5 (d) suggests that OSGA is superior to the others regarding ISNR. To see the constructed results visually, Figure 6 is considered, where it reports last function value and PSNR of the algorithms. It can be seen that OSGA obtains the best function value 79252.60 and the best PSNR 32.23 among the considered algorithms. We finally consider 72 widely used test images to do more comparisons, where blurred/noisy images are generated by 9 × 9 uniform blur and adding Gaussian noise with SNR = 40 dB. The images recovered by TwIST, SpaRSA, FISTA and OSGA, where the algorithm are stopped after 100 iterations. The results of reconstructions are summarized in Table 1, where we report the best function value (F), PSNR and CPU time (T) suggesting OSGA outperforms the others. The results are illustrated in Figure 7, where the comparisons are based on the performance profile of Dolan and Moré in [27]. Figure 7 shows that OSGA wins with 84% and 93% score for function values and PSNR among all others. It also indicates that all algorithms recover images successfully. 16 Masoud Ahookhosh 0.7 0.7 0.7 IST TwIST FISTA OSGA 0.6 IST TwIST FISTA OSGA 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 100 200 300 400 500 600 (a) rel. 1 vs. iterations (chit = 5) 0 0 50 100 150 200 250 300 350 400 (b) rel. 1 vs. iterations (chit = 10) −1 −2 −3 50 60 (d) rel. 2 vs. time (chit = 5) 0 10 −1 20 30 40 50 −3 10 −2 10 −3 0 50 100 150 200 250 300 350 400 450 (h) rel. 2 vs. iterations (chit = 10) 10 0 25 20 20 20 15 10 5 5 5 400 500 200 250 300 350 400 600 (j) ISN R vs. iterations (chit = 5) 0 450 500 IST TwIST FISTA OSGA 10 300 150 IST TwIST FISTA OSGA 10 200 100 15 IST TwIST FISTA OSGA 100 50 (i) rel. 2 vs. iterations (chit = 20) 25 15 60 (f) rel. 2 vs. time (chit = 20) 25 0 50 IST TwIST FISTA OSGA −3 600 (g) rel. 2 vs. iterations (chit = 5) 0 40 −1 −2 500 30 10 10 400 20 IST TwIST FISTA OSGA −2 300 10 0 10 10 500 10 IST TwIST FISTA OSGA 200 0 60 −1 10 450 −3 10 0 100 400 (c) rel. 1 vs. iterations (chit = 20) 10 0 10 0 350 −2 (e) rel. 2 vs. time (chit = 10) 10 10 300 10 −3 40 250 IST TwIST FISTA OSGA −2 30 200 10 10 20 150 IST TwIST FISTA OSGA 10 10 100 −1 −1 10 0 50 0 IST TwIST FISTA OSGA 10 0 10 10 10 0 450 0 10 IST TwIST FISTA OSGA 0.6 0 50 100 150 200 250 300 350 400 450 (k) ISN R vs. iterations (chit = 10) 0 0 50 100 150 200 250 300 350 400 450 500 (l) ISN R vs. iterations (chit = 20) Fig. 3: A comparison among IST, TwIST, FISTA and OSGA for inpainting the 512 × 512 Head CT image when 40% of its data is missing. The algorithms are stopped after 50 seconds of the running time. Each column stands for comparisons with a different number of Chambolle’s iterations. The first row denotes the relative error rel. 1 = kxk − x∗ k2 /kx∗ k2 of points versus iterations, the second row shows the relative error rel. 2 = (fk − f ∗ )/(f0 − f ∗ ) of function values versus time, the third row illustrates rel.2 versus iterations, and fourth row stands for ISNR (33) versus iterations. Optimal subgradient algorithms with application to large-scale linear inverse problems 17 (a) Original image (b) Blurred/noisy image (c) OSGA: F = 187564.02, PSNR = 33.13 (d) IST (chit = 5): F = 710603.30, PSNR = 16.15 (e) IST (chit = 10): F = 1046712.18, PSNR = 13.61 (f) IST (chit = 20): F = 1284434.37, PSNR = 12.08 (g) TwIST (chit = 5): F = 189347.03, PSNR = 29.84 (h) TwIST (chit = 10): F 189533.03, PSNR = 29.94 = (i) TwIST (chit = 20): F 192946.47, PSNR = 29.26 (j) FISTA (chit = 5): F = 188676.25, PSNR = 30.37 (k) FISTA (chit = 10): F 189145.34, PSNR = 30.36 = (l) FISTA (chit = 20): F = 198692.14 & PSNR = 30.24 = Fig. 4: Inpainting of the 512 × 512 Head CT image by IST, TwIST, FISTA, OSGA when they stopped after 50 seconds of the running time. 18 Masoud Ahookhosh 0 10 0.07 TwIST SpaRSA FISTA OSGA 0.06 10 relative error of function values 0.05 relative error of steps TwIST SpaRSA FISTA OSGA −1 0.04 0.03 −2 10 −3 10 −4 10 0.02 −5 10 0.01 0 −6 0 10 20 30 40 50 60 iterations 70 80 90 100 10 0 5 15 20 25 time (a) rel. 1 vs. iterations, chit = 5 (b) rel. 2 vs. time, chit = 5 0 5 10 TwIST SpaRSA FISTA OSGA −1 10 4.5 4 3.5 −2 10 SNR improvement relative error of function values 10 −3 10 3 2.5 2 −4 10 1.5 TwIST SpaRSA FISTA OSGA 1 −5 10 0.5 −6 10 0 10 20 30 40 50 60 iterations 70 80 (c) rel. 2 vs. iterations, chit = 5 90 100 0 0 10 20 30 40 50 60 iterations 70 80 90 100 (d) ISN R vs. iterations, chit = 5 Fig. 5: A comparison among TwIST, SpaRSA, FISTA and OSGA for deblurring the 512 × 512 Elaine image with 9×9 uniform blur and Gaussian noise with SNR = 40 dB. TwIST, SpaRSA and FISTA exploit chit = 5, and all algorithms are stopped after 100 iterations. In the figure: (a) stands for the relative error rel. 1 = kxk −x∗ k2 /kx∗ k2 of points versus iterations, (b) shows the relative error rel. 2 = (fk −f ∗ )/(f0 −f ∗ ) of function values versus time, (c) demonstrates rel.2 versus iterations, and (d) depicts ISNR (33) versus iterations. Optimal subgradient algorithms with application to large-scale linear inverse problems (a) Original image (b) Blurred/noisy image (c) TwIST: f = 80733.24, PSNR = 31.99, T = 20.07 (d) SpaRSA: f = 79698.36, PSNR = 32.12, T = 12.16 (e) FISTA: f = 80795.59, PSNR = 31.92, T = 11.84 (f) OSGA: f = 79252.60, PSNR = 32.23, T = 14.82 19 Fig. 6: Deblurring of the 512×512 Elaine image with 9×9 uniform blur and Gaussian noise with SNR = 40 dB by TwIST, SpaRSA, FISTA and OSGA. The algorithms stopped after 100 iterations, and the number of Chambolle iterations for TwIST, SpaRSA and FISTA is set to chit = 5. 20 Masoud Ahookhosh 1 0.9 0.8 P(rp,s ≤ τ : 1 ≤ s ≤ ns) 0.7 0.6 0.5 0.4 0.3 0.2 TwIST SpaRSA FISTA OSGA 0.1 0 1 1.05 τ 1.1 1.15 (a) performance profile for function values 1 0.9 0.8 P(rp,s ≤ τ : 1 ≤ s ≤ ns) 0.7 0.6 0.5 0.4 0.3 TwIST SpaRSA FISTA OSGA 0.2 0.1 0 1 1.05 τ 1.1 1.15 (b) performance profile for PSNR Fig. 7: Deblurring of 72 test images with isotropic total variation using TwIST, SpaRSA FISTA, OSGA in 25 seconds of the running time. TwIST SpaRSA FISTA OSGA Image name Dimension F PSNR T F PSNR T F PSNR T F PSNR T Aerial Airplane Cameraman Chemical plant Clock Fingerprint 1 Lena Moon surface Pattern 1 Pattern 2 Star Aerial Airfield Airplane 1 Airplane 2 APC Barbara Blobs Blonde woman Boat Cameraman Car & APCs 1 Car & APCs 2 Clown Crowd 1 Crowd 2 Darkhair woman Dollar Elaine Fingerprint Flinstones Girlface Goldhill Head CT Houses 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 256 × 256 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 53380.99 16527.67 37077.60 37890.81 31581.71 54289.58 30700.61 16648.60 231150.20 176.53 23482.48 155827.40 169428.35 105947.22 29318.44 45314.07 114725.54 191940.70 96330.41 118404.12 111579.37 75501.10 69193.23 137437.13 127463.77 213890.70 85530.72 249508.76 80406.78 256043.16 245122.90 82611.48 110067.92 178265.55 250441.66 24.99 33.75 27.39 26.80 31.07 20.98 28.33 30.22 18.07 48.90 24.31 25.94 26.06 32.00 35.86 33.18 25.45 20.63 29.63 29.87 32.98 32.67 32.52 31.71 30.55 28.16 37.31 21.49 32.01 26.09 25.64 33.00 30.69 30.91 25.79 5.80 6.50 4.94 6.64 6.04 7.06 7.98 8.76 4.48 7.27 6.06 23.55 23.87 23.43 19.22 24.09 27.10 24.31 30.11 27.41 29.51 22.48 30.35 23.04 25.33 18.55 24.19 18.39 20.28 22.58 17.50 24.73 23.78 16.10 17.79 52966.30 16563.06 37230.96 37902.64 31740.38 55070.52 30610.62 16552.11 232285.41 190.83 23596.24 166304.12 169301.74 105942.26 29609.95 45128.61 114235.16 191945.61 95947.22 118194.91 111655.52 74826.32 68645.53 137241.10 127563.32 214814.30 85441.92 273059.31 79667.61 256803.08 246272.55 82842.10 109700.41 179561.46 262200.38 25.14 33.72 27.11 26.87 30.95 21.12 28.31 30.29 17.85 48.92 24.26 25.93 25.97 32.06 35.14 33.20 25.47 20.59 29.71 29.82 32.87 32.90 32.65 31.77 30.72 28.35 37.49 21.42 32.13 26.19 25.45 33.06 30.78 30.05 25.69 3.56 3.39 3.28 3.48 4.12 3.68 4.12 4.39 4.32 3.86 4.04 17.01 16.91 14.29 14.83 14.60 14.39 15.92 15.40 16.12 15.86 13.88 15.55 13.31 13.58 12.29 12.47 12.34 12.45 12.34 12.31 12.18 12.49 12.63 12.47 53513.56 16549.09 37096.87 37934.72 31604.08 54374.17 30730.93 16718.75 231937.14 186.27 23514.17 156015.30 169528.18 106078.89 29329.10 45425.04 115004.02 192059.89 96567.05 118590.52 111656.63 75868.96 69636.88 137550.26 127708.39 214494.66 85854.12 249686.20 80765.32 256811.98 245410.87 82896.35 110273.58 178400.71 250589.16 25.10 33.69 27.44 26.80 31.18 20.98 28.39 30.27 18.09 48.90 24.37 26.01 26.14 32.06 35.90 33.22 25.43 20.63 29.65 29.90 33.00 32.59 32.43 31.75 30.61 28.16 36.96 21.50 31.94 25.85 25.66 32.89 30.71 30.93 25.83 3.33 3.28 3.28 3.32 3.33 3.46 3.37 3.68 3.61 3.41 3.81 12.86 13.53 13.57 13.44 13.48 13.84 14.03 13.66 14.02 13.30 13.65 13.84 12.79 12.21 12.56 13.09 12.52 14.08 12.66 11.99 12.50 12.27 11.94 11.62 52832.73 16630.05 37089.05 37562.79 31502.30 54256.89 30468.62 16505.01 232025.68 172.69 23370.80 154008.47 168303.91 105076.82 29285.75 44985.53 113697.78 192043.10 95219.10 117624.44 110900.70 74614.10 68233.64 136517.93 125715.43 210665.21 84848.24 248817.22 79206.17 253486.89 243375.64 81871.62 109353.67 177648.72 250375.51 25.66 34.36 27.82 27.25 31.64 21.22 28.73 30.54 18.17 48.89 24.49 26.51 26.63 32.70 35.79 33.48 25.60 20.85 30.02 30.33 33.53 32.98 32.85 32.12 31.33 29.00 37.77 21.68 32.24 26.49 26.47 33.44 31.06 31.52 25.97 3.72 3.62 3.80 3.72 4.28 4.23 4.19 4.67 4.13 4.26 4.37 16.53 17.68 18.18 18.10 17.97 18.17 17.21 17.74 16.84 17.56 17.24 17.23 16.65 15.61 15.37 15.31 15.34 14.84 15.71 15.50 15.63 15.66 15.58 15.74 Optimal subgradient algorithms with application to large-scale linear inverse problems Table 1. Deblurring with isotropic TV (F, T & PSNR) 21 22 Table 1. Deblurring results (F, T & PSNR) (continued) 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 512 × 512 600 × 600 600 × 600 600 × 600 600 × 600 600 × 600 92142.42 600 × 600 1024 × 1024 1024 × 1024 1024 × 1024 1024 × 1024 1024 × 1024 1024 × 1024 155925.23 154117.24 91492.93 258205.38 43626.15 169572.40 126804.63 149477.52 85213.65 101362.85 108105.61 110681.19 154177.76 63474.30 85867.65 93860.07 73437.12 106703.13 105112.40 255921.91 257153.15 261572.25 222243.70 71117.95 201150.09 197734.24 129529.71 208297.87 219117.86 600 × 600 83431.60 423652.87 324683.85 462035.52 289127.64 921586.14 816333.83 159823.01 25.82 28.88 31.80 20.82 37.40 27.50 29.83 25.90 35.59 32.89 29.83 33.14 50.49 31.45 29.39 29.68 31.78 28.69 28.74 20.11 20.06 20.01 20.81 35.27 23.93 23.85 33.80 29.22 28.00 40.12 35.32 28.25 28.48 30.22 37.12 21.89 23.60 29.36 16.35 24.72 24.05 21.23 24.98 19.19 24.84 23.63 21.32 19.18 20.98 29.36 17.53 30.85 23.93 30.64 23.54 27.71 28.27 23.93 27.90 24.69 22.17 26.02 29.71 27.26 28.46 35.19 32.76 31.82 21.89 75.10 97.22 88.40 86.41 56.31 94.63 26.54 156039.11 154535.79 92149.44 260738.06 43713.46 169527.68 126502.75 150370.73 85179.01 102312.47 107788.51 110432.93 154454.47 63215.81 85277.68 93348.41 73088.97 106087.95 104405.27 255818.43 256665.77 266971.93 220980.31 71246.57 201437.41 198225.38 129835.46 207720.39 220633.95 92885.81 85917.89 422352.54 323003.05 460696.14 289405.02 927968.42 812073.87 160619.72 25.76 28.86 31.78 20.57 37.24 27.44 29.84 26.01 35.61 32.86 29.91 33.17 48.95 31.48 29.49 29.77 31.93 28.79 28.84 20.03 20.02 20.03 20.85 35.60 23.91 23.83 33.82 29.30 28.08 40.82 31.69 28.25 28.54 30.24 37.14 21.68 23.67 29.30 12.83 13.30 12.50 15.70 14.06 14.50 13.65 12.48 13.12 13.45 14.87 13.95 13.19 14.99 14.39 15.40 12.35 14.49 14.09 12.68 15.31 12.64 12.40 12.69 21.04 19.98 17.79 19.67 17.94 17.67 18.52 53.99 51.01 45.67 45.00 44.94 49.28 15.68 156013.11 154292.14 91620.97 258430.66 43668.81 169664.95 127007.96 149889.46 85397.57 101591.74 108301.20 110858.04 154228.55 63643.84 86445.51 94147.81 73768.38 107119.20 105581.13 256223.57 257413.32 261772.59 222443.50 71547.26 201282.62 197847.53 129693.50 208547.46 219692.19 92779.78 83775.26 424082.24 325787.77 462855.60 289515.55 923132.76 817252.58 160106.63 25.89 28.95 31.80 20.88 37.47 27.55 29.89 25.88 35.51 32.89 29.90 33.16 50.21 31.44 29.38 29.66 31.77 28.69 28.71 20.10 20.08 20.04 20.82 35.08 23.95 23.86 33.83 29.18 27.90 39.56 34.38 28.32 28.51 30.23 36.99 21.88 23.61 29.34 11.88 13.10 12.02 12.07 11.44 11.82 11.82 12.82 12.20 13.44 12.26 12.12 12.22 11.73 11.89 12.91 12.25 12.18 12.44 12.06 12.52 11.78 11.82 12.10 18.03 17.43 17.41 17.09 17.45 17.34 17.22 59.41 47.50 47.90 48.80 48.89 51.87 14.86 155847.17 152700.95 90544.16 258452.85 43420.69 169081.03 126280.35 148035.40 84823.91 100470.74 107122.20 109880.54 154711.06 62876.00 84963.72 93238.90 72621.40 105751.32 104103.47 255585.39 256847.89 261064.29 222099.26 70038.80 201705.89 198123.72 129342.37 207095.69 218101.17 90891.80 85334.45 420651.39 320916.50 457747.96 288520.76 927122.49 811272.73 159106.67 26.01 29.42 32.27 20.92 37.69 27.72 30.18 26.30 36.10 33.44 30.31 33.65 47.79 31.70 29.70 29.96 32.17 29.00 29.06 20.39 20.36 20.33 21.06 35.77 24.09 24.00 33.92 29.40 28.30 41.15 32.38 28.63 28.81 30.68 37.32 21.84 23.95 29.65 15.81 15.49 14.87 15.51 15.68 15.55 15.78 15.32 15.75 15.50 15.69 15.66 15.72 15.71 15.55 15.79 15.44 15.74 15.46 15.78 15.50 15.74 15.33 15.65 21.04 21.49 21.52 21.63 21.35 21.54 21.27 59.71 59.67 59.85 59.31 61.62 63.54 18.47 Masoud Ahookhosh Kiel Lake Lena Lena numbers Liftingbody Lighthouse Livingroom Mandril MRI spine Peppers Pirate Smiling woman Squares Tank 1 Tank 2 Tank 3 Truck Truck & APCs 1 Truck & APCs 2 Washington DC 1 Washington DC 2 Washington DC 3 Washington DC 4 Zelda Dark blobs 1 Dark blobs 2 House Ordered matches Random matches Rice Shepp-logan phantom Airport Pentagon Pirate Rose Testpat Washington DC Average Optimal subgradient algorithms with application to large-scale linear inverse problems 23 3.2 Comparisons with first-order methods As in the first section discussed, there is an excessive gap between the optimal complexity of first-order √ methods for solving smooth and nonsmooth convex problems to get an ε-solution, where it is O(1/ ε) for smooth problems and O(1/ε2 ) for nonsmooth ones. Inspired by this fact, we consider problems of both cases in the form (8) and compare the best-known first-order algorithms with OSGA. In spite of the fact that some first-order methods use smoothing techniques to solve nonsmooth problems [16, 17, 43], we only consider those that directly solve the original problem. In the other words, we only consider algorithms exploiting the nonsmooth first-order black-box oracle or employing the proximity operators along with the smooth first-order black-box oracle. Hence, during our implementations, the following algorithms are considered • PGA : Proximal gradient algorithm [52]; • NSDSG: Nonsummable diminishing subgradient algorithm [19]; • FISTA: Beck and Tebolle’s fast proximal gradient algorithm [9]; • NESCO: Nesterov’s composite optimal algorithm [46]; • NESUN: Nesterov’s universal gradient algorithm [48]; • NES83: Nesterov’s 1983 optimal algorithm [40]; • NESCS: Nesterov’s constant step optimal algorithm [42]; • NES05: Nesterov’s 2005 optimal algorithm [43]. The algorithms NES83, NESCS and NES05 originally proposed to solve smooth convex problems, where√ they use the smooth first-order black-box oracle and attain the optimal complexity of order O(1/ ε). Although, obtaining the optimal complexity bound of smooth problems is computationally very interesting, the considered problem (8) is commonly nonsmooth. Recently, Lewis and Overton in [36] investigated the behaviour of the BFGS method for nonsmooth nonconvex problems, where their results are interesting. Furthermore, Nesterov in [46] shows that the subgradient of composite functions preserve all important properties of the gradient of smooth convex functions. These facts, motivate us to do some computational investigations on the behaviour of optimal smooth first-order methods in the presence of the nonsmooth first-order black-box oracle, where the oracle passes a subgradient of composite objective functions instead of the gradient of smooth parts. In particular, we here consider Nesterov’s 1983 algorithm, [40], and adapt it to solve (8) by simply preparing the subgradient of the functions as follows: Algorithm 4: NES83 Nesterov’s 1983 algorithm for multi-term affine composite functions Input: select z such that z 6= y0 and gy0 6= gz ; local parameters: y0 ; ε > 0; Output: xk ; Ψxk ; begin a0 ← 0; x−1 ← y0 ; compute gy0 = gΨ (y0 ) using NFO-G; compute gz = gΨ (z) using NFO-G; α−1 ← ky0 − zk/kgy0 − gz k; while stopping criteria do not hold do α̂k ← αk−1 ; x̂k ← yk − α̂k gyk ; compute Ψx̂k using NFO-F; while Ψx̂k < Ψyk − 12 α̂k kgyk k2 do α̂k ← ρα̂k ; x̂k ← yk − α̂k gyk ; compute Ψx̂k using NFO-F; end xk ← x̂k ; Ψxk ← Ψx̂k ;  αk ← α̂k ; p ak+1 ← 1 + 4a2k + 1 /2; yk+1 ← xk + (ak − 1)(xk − xk−1 )/ak+1 ; compute gyk+1 = gΨ (yk+1 ) using NFO-G; end end 24 Masoud Ahookhosh In a similar way considered for NES83, the algorithms NESCS and NES05 are respectively adapted from Nesterov’s constant step [42] and Nesterov’s 2005 algorithms [43]. Using this technique, NES83, NESCS and NES05 are able to solve nonsmooth problems as well, however there is no theory to support their convergence at the moment. Among the algorithms considered, PGA, NSDSG, FISTA, NESCO and NESUN are originally introduced to solve nonsmooth problems, where NSDSG only employs the nonsmooth first-order black-box oracle and the others use the smooth first-order black-box oracle together with the proximal mapping. In the sense of complexity theory, NSDSG is just optimal for nonsmooth problems with the order O(1/ε2 ), PGA is not optimal with the√complexity of the order O(1/ε), and FISTA, NESCO and NESUN are optimal with the order O(1/ ε). These complexity bounds clearly imply that NSDSG and PGA are not theoretically comparable with the others, however we still consider them in our comparisons to see how much the optimal first-order methods are computationally interesting than the classical first-order methods, especially for solving large-scale problems in applications. To get a clear assessment of OSGA, we compare it with the algorithms only using the nonsmooth firstorder oracle (NSDSG, NES83, NESCS, NES05) and with those employing the smooth first-order oracle together with the proximal mapping (PGA, NSDSG, FISTA, NESCO, NESUN) in separate comparisons. Most of these algorithms require the global Lipschitz constant L to determine a step size. While NESCS, NES05, PGA and FISTA use the constant L to determine a step size, NESCO and NESUN get a lower approximation of L and adapt it by backtracking line searches, for more details see [46] and [48]. In our implementations of NESCO and NESUN, similar to [11], we determine the initial estimate L0 by L0 = kg(x0 ) − g(z0 )k∗ , kx0 − z0 k where x0 and z0 are two arbitrary points such that x0 6= z0 and g(x0 ) 6= g(z0 ). In [40], Nesterov used an equivalent scheme to √ determine an initial step size that can be seen in NES83. For NSDSG, the step size is calculated by α0 / k, where k is the current iteration number and α0 is a positive constant that should be specified as big as possible such that the algorithm is not divergent, see [19]. All the other parameters of the algorithms are set to those reported in the corresponding papers. We here consider solving an underdetermined system Ax = y, where A is a m × n random matrix and y is a random m-vector. This problem is frequently appeared in applications, and the goal is to determine x by optimization techniques. Considering the ill-conditioned feature of this problem, the most popular optimization models are (3) and (4), where (3) is clearly smooth and (4) is nonsmooth. In each case, we consider both dense and sparse data. More precisely, the dense data is constructed by setting A = rand(m, n), y = rand(1, m) and the initial point x0 = rand(1, n), and the sparse data is created by A = sprand(m, n, 0.05), y = sprand(1, m, 0.05) and the initial point x0 = rand(1, n, 0.05). The problems use m = 5000 and n = 10000 and employ λ = 1 similar to [46]. For the dense problems, the algorithms are stopped after 60 seconds of the running time, and for the sparse problems, they are stopped after 45 seconds. Let to define L̂ = max kai k2 , 1≤i≤n where ai for i = 1, 2, · · · , n are columns of A. For the dense and sparse problems, NESCS, NES05, PGA and FISTA respectively use L = 104 L̂ and L = 102 L̂. Moreover, NSDSG employs α0 = 10−7 and α0 = 10−4 for the dense and sparse problems, respectively. To assess the considered algorithms, we now solve (3) and (4) for both dense and sparse data, where the results are respectively illustrated in Figures 8 and 9. In Figure 8, subfigures (a) and (b) illustrate the results of dense data, while subfigures (c) and (d) demonstrate the results of sparse data. From this figure, it can be seen that NSDSG and PGA are not competitive with the others confirming the theoretical results about their complexity. More specifically, Figure 8 (a) and (c) imply that NES83 and NES05 are competitive and produce better results than NESCS. In this case, OSGA achieves the best results, especially for sparse data. Figure 8 (b) and (d) suggest that NESCO, NESUN and FISTA are comparable, but still OSGA generates better results. Figure 9 depicts the results of reconstructing an approximation solution for the nonsmooth problem (4) using the considered algorithms. The same interpretation for Figure 9 indicates that most of the optimal methods produce better results than classical ones and are competitive with each other. However, OSGA commonly obtains the best results among them. Summarizing the results of both figures, we see that the adapted algorithms NES83, NESCS and NES05 are comparable or even better than the the other optimal algorithms, especially for NES83. Optimal subgradient algorithms with application to large-scale linear inverse problems 25 11 10 11 10 NSDSG NES83 NESCS NES05 OSGA 10 10 9 10 9 10 8 10 8 10 function values function values PGA FISTA NESCO NESUN OSGA 10 10 7 10 6 10 7 10 6 10 5 5 10 10 4 4 10 10 3 3 10 10 2 10 2 0 10 1 10 2 10 iterations 3 10 10 4 10 0 10 (a) Function values vs. iterations 1 10 2 10 iterations 3 10 4 10 (b) Function values vs. iterations 5 5 10 4 10 3 10 10 4 10 3 function values function values 10 2 10 1 10 1 10 NSDSG NES83 NESCS NES05 OSGA 0 10 PGA FISTA NESCO NESUN OSGA 0 10 −1 −1 10 2 10 0 10 1 10 2 10 iterations 3 10 (c) Function values vs. iterations 4 10 10 0 10 1 10 2 10 iterations 3 10 4 10 (d) Function values vs. iterations Fig. 8: A comparison among the considered algorithms for solving the l22 − l22 problem (3): subfigures (a) and (b) stand for the dense data, where the algorithms stopped after 60 seconds; subfigures (c) and (d) show the results of the sparse data, where the algorithms stopped after 45 seconds. 3.3 Sparse optimization In recent decades, taking advantage of sparse solutions by using the structure of problems has become a common task in various applications of applied mathematics and particularly optimization, where applications are mostly arising in the fields of signal and image processing, statistics, machine learning and data mining. In the most cases, the problem involves high-dimensional data, and a small number of measurements is accessible, where the core of these problems involves an optimization problem. Thanks to the sparsity of solutions and the structure of problems, these optimization problems can be solved in reasonable time even for the extremely high-dimensional data sets. Basis pursuit, lasso, wavelet-based deconvolution and compressed sensing are some examples, where the last one receives lots of attentions during the recent years. Among fields involving sparse optimization, compressed sensing , also called compressive sensing and compressive sampling , is a novel sensing/sampling framework for acquiring and recovering objects like a sparse image or signal in the most efficient way possible with employing an incoherent projecting basis. On the one hand, conventional processes in the image/signal acquisition from frequency data 26 Masoud Ahookhosh 11 11 10 10 NSDSG NES83 NESCS NES05 OSGA 10 10 9 10 PGA FISTA NESCO NESUN OSGA 10 10 9 10 8 10 function values function values 8 10 7 10 6 10 7 10 6 10 5 10 5 4 10 3 10 10 4 10 2 10 3 0 10 1 10 2 10 iterations 3 10 10 4 10 0 10 (a) Function values vs. iterations 1 10 2 10 iterations 3 10 4 10 (b) Function values vs. iterations 5 10 5 10 4 10 4 10 3 10 3 function values function values 10 2 10 1 10 1 10 NSDSG NES83 NESCS NES05 OSGA 0 10 PGA FISTA NESCO NESUN OSGA 0 10 −1 −1 10 2 10 10 0 10 1 10 2 10 iterations 3 10 (c) Function values vs. iterations 4 10 0 10 1 10 2 10 iterations 3 10 4 10 (d) Function values vs. iterations Fig. 9: A comparison among the considered algorithms for solving the l22 − l1 problem (4): subfigures (a) and (b) stand for the dense data, where the algorithms stopped after 60 seconds; subfigures (c) and (d) show the results of the sparse data, where the algorithms stopped after 45 seconds. follow the Nyquist-Shannon density sampling theorem declaring that the number of samples required for reconstructing an image matches the number of pixels in the image and for recovering a signal without error is devoted by its bandwidth. On the other hand, it is known that most of the data we acquire can be thrown away with almost no perceptual loss. Then, the question is: why all the data should be acquired while most of them will be thrown away? Also, can we only recover parts of the data that will be useful in the final reconstruction? In response to these questions, a novel sensing/sampling approach was introduced which goes against common wisdom in data acquisition, namely compressed sensing, for example see [21, 22, 28] and references therein. It is known that compressed sensing supposes that the considered object, image or signal, has a sparse representation in some bases or dictionaries, and the considered representation dictionary is incoherent with the sensing basis, see [20, 28]. The idea behind this framework is centred on how many measurements are necessary to reconstruct an image/signal and the related nonlinear reconstruction techniques needed to recover this image/signal. It means that the object is recovered reasonably well from highly undersampled data, in spite of violating the traditional Nyquist-Shannon sampling constraints. In the other words, once an object is known to be sparse in a specific dictionary, one should find a set of measurements or sensing bases supposed to be incoherent with the dictionary. Afterwards, an underdetermined system of linear Optimal subgradient algorithms with application to large-scale linear inverse problems 27 equations of the form (1) is emerged, which is usually ill-conditioned. Some theoretical results in [21, 22, 28] yield that the minimum number of measurements required to reconstruct the original object provides a specific pair of measurement matrices and optimization problems. Thus, the main challenge is how to reformulate this inverse problem as an appropriate minimization problem involving regularization terms and solve it by a suitable optimization technique. We now consider the inverse problem (1) in the matrix form with A ∈ Rm×n , x ∈ Rn and y ∈ Rm , where the observation is deficient, m < n, as a common interest in compressed sensing. Therefore, the object x can not be recovered from the observation y, even in noiseless system Ax = y, unless some additional assumptions like sparsity of x is presumed. We here consider a sparse signal x ∈ Rn that should be recovered by the observation y ∈ Rm , where A is obtained by first filling it with independent samples of the standard Gaussian distribution and then orthonormalizing the rows, see [30]. In our experiments, we consider m = 5000 and n = 10000, where the original signal x consist of 300 randomly placed ±1 spikes, and y is generated using (1) by σ 2 = 10−6 . We reconstruct x using (4), where the regularization parameter is set to λ = 0.1kAyk∞ or λ = 0.001kAyk∞ . The algorithms stopped after 20 seconds of the running time, and the results are illustrated in Figure 10. 4 1.4 10 NSDSG NES83 NESCS NES05 PGA FISTA NESCO NESUN OSGA 3 1 0.8 MSE function values 10 NSDSG NES83 NESCS NES05 PGA FISTA NESCO NESUN OSGA 1.2 0.6 2 10 0.4 0.2 1 10 0 10 1 2 10 0 0 10 3 10 10 1 3 10 10 iterations (a) Function values vs. iterations (b) MSE vs. iterations 4 10 1.4 NSDSG NES83 NESCS NES05 PGA FISTA NESCO NESUN OSGA 3 10 2 10 NSDSG NES83 NESCS NES05 PGA FISTA NESCO NESUN OSGA 1.2 1 0.8 MSE function values 2 10 iterations 0.6 1 10 0.4 0 10 0.2 −1 10 0 10 1 2 10 10 iterations (c) Function values vs. iterations 3 10 0 0 10 1 2 10 10 3 10 iterations (d) MSE vs. iterations Fig. 10: Recovering a noisy sparse signal using l22 − l1 problem (4), where the algorithms stopped after 20 seconds. Subfigures (a) and (b) show the results for the regularization parameter λ = 0.1kAyk∞ , while subfigures (c) and (d) indicate the results for the regularization parameter λ = 0.001kAyk∞ . 28 Masoud Ahookhosh Original signal(n = 10000, number of nonzeros = 300) 1 0 −1 0 1000 2000 3000 500 1000 1500 1000 2000 1000 4000 5000 Noisy signal 6000 7000 8000 9000 10000 2000 2500 3000 Direct solution using the pseudo inverse 3500 4000 4500 5000 3000 4000 5000 6000 NES83 (MSE = 0.000729386) 7000 8000 9000 10000 2000 3000 4000 5000 6000 NESCS (MSE = 0.000947418) 7000 8000 9000 10000 1000 2000 3000 4000 5000 6000 NES05 (MSE = 0.000870573) 7000 8000 9000 10000 1000 2000 3000 4000 5000 6000 FISTA (MSE = 0.000849244) 7000 8000 9000 10000 1000 2000 3000 4000 5000 6000 OSGA (MSE = 0.000796183) 7000 8000 9000 10000 1000 2000 3000 4000 7000 8000 9000 10000 1 0 −1 0 1 0 −1 0 1 0 −1 0 1 0 −1 0 1 0 −1 0 1 0 −1 0 1 0 −1 0 5000 6000 Fig. 11: Recovering a noisy sparse signal using l22 − l1 problem (4) by NES83, NESCS, NES05, FISTA and OSGA. The algorithms stopped after 20 seconds, and the quality measure MSE is defined by (34). Figure 10 illustrates the results of the sparse signal recovery by NSDSG, NSE83, NESCS, NES05, PGA, FISTA, NESCO, NESUN and OSGA, where the the results are compared in the sense of function values and MSE. Indeed, MSE is a usual measure of comparisons in signal reconstruction defined by MSE(x) = 1 kx − x0 k2 , n (34) where x0 is the original signal. The results of Figure 10 indicate that NSDSG, PGA, NESCO and NESUN are unsuccessful to recover the signal accurately, where NESCO and NESUN could not leave the backtracking line searches in a reasonable number of inner iterations which dramatically decrease the step sizes. In Figure 10 (a) and (b), we see that NSE83, NESCS, NES05, FISTA and OSGA can reconstruct the signal accurately, however NES83 and OSGA are much faster than the others regarding both function values and MSE. To see the sensitivity of the algorithms to the regularization parameter, Figure 10 (c) and (d) demonstrate the results of recovery for the smaller regularization parameter λ = 0.001kAyk∞ . Optimal subgradient algorithms with application to large-scale linear inverse problems 29 It is clear that NESCS, NES05 and FISTA are more sensitive than NES83 and OSGA to the regularization parameter, where they could not get an accurate approximation of the original signal when the regularization parameter decreased. To show the results of our test visually, we illustrate the recovered signal of the experiment reported in Figure 10 (a) and (b) in Figure 11. Since NSDSG, PGA, NESCO and NESUN could not recover the signal accurately, we do not consider them in Figure 11. We also add the direct solution x = AT (AAT )−1 y to our comparison indicating that this solution is very inaccurate. The results suggest that NSE83, NESCS, NES05, FISTA and OSGA recover the signal visually acceptable. Summarizing the results of Figures 10 and 11, it can be seen that NES83 and OSGA just need about 15 iterations to achieve acceptable results, while NESCS, NES05 and FISTA require about 100 iterations to reach the same accuracy. 4 OSGA software package The first version of the OSGA software package for solving unconstrained convex optimization problems is publicly available at http://homepage.univie.ac.at/masoud.ahookhosh/uploads/OSGA_v1.1.tar.gz where it is implemented in MATLAB programming language. Some examples are available to show how the user can implement OSGA. The interface to each subprogram in the package is fully documented in the corresponding file. Moreover, the OSGA user’s manual [1] describes the design of the package and how the user can solve his/her own problems. 5 Concluding remarks Simple gradient and subgradient procedures are historically first numerical schemes proposed for solving smooth and nonsmooth convex optimization, respectively. Theoretical and computational results indicate that they are very slow to be able to solve large-scale practical problems in applications. However, these algorithms have a simple structure and need a low amount of memory in implementation. These features cause that first-order algorithms have received much attention during last few decades. In particular, both theoretical and computational aspects of optimal first-order methods are very interesting. In this paper, we study and develop an optimal subgradient algorithm called OSGA for the class of multi-term affine composite functions. The algorithm is flexible enough to handle a broad range of structured convex optimization problems and has a simple structure. It also does not need global parameters like Lipschitz constants except for the strong convexity parameter, which can be set to zero if it is unknown. We consider examples of linear inverse problems and report extensive numerical results by comparing OSGA with some state-of-the-art algorithms and software packages. The numerical results suggest that OSGA is computationally competitive or even better than the considered state-of-the-art solvers designed specifically for solving these special problems. It can often reach the fast rate of convergence even beyond the theoretical rate suggesting that OSGA is a promising algorithm for solving a wide range of large-scale convex optimization problems. We also develop some optimal first-order methods originally proposed for solving smooth convex optimization to deal with nonsmooth problems as well. To this end, we simply pass a subgradient of the nonsmooth composite objective function to these algorithms in place of the gradient. This simply means that we change their original oracle by the nonsmooth black-box oracle. The numerical results show that the adapted algorithms are promising and can produce competitive or even better results than the other considered algorithms specifically proposed for nonsmooth problems. However, they are not supported by any theoretical results at the moment, so it might be an interesting topic for a future research. The current paper only considers unconstrained convex optimization problems, however, the constrained version for some simple constraints developed in [2, 3, 5]. Furthermore, for problems involving costly direct and adjoint operators, a version of OSGA using the multi-dimensional subspace search is established in [5, 6]. Finally, many more applications like sparse optimization with the l1 /l2 or l1 /l∞ norm could be considered, but this is beyond the scope of the current work and remains for a future research. 30 Masoud Ahookhosh Acknowledgement. We would like to thank Arnold Neumaier for a careful reading of the manuscript and his insightful comments and feedback. We are also grateful to the authors of TwIST, SpaRSA and FISTA software packages for generously making their codes available on the web. References 1. Ahookhosh, M.: User’s manual for OSGA (Optimal SubGradient Algorithm), (2014) http://homepage.univie.ac.at/ masoud.ahookhosh/uploads/User’s_manual_for_OSGA.pdf [29] 2. Ahookhosh, M., Neumaier, A.: An optimal subgradient method for large-scale bound-constrained convex optimization, Manuscript, University of Vienna, (2014) [29] 3. Ahookhosh, M., Neumaier, A.: Optimal subgradient methods for structured convex constrained optimization I: theoretical results, Manuscript, University of Vienna, (2014) [29] 4. Ahookhosh, M., Neumaier, A.: Optimal subgradient methods for structured convex constrained optimization II: numerical results, Manuscript, University of Vienna, (2014) [29] 5. Ahookhosh, M., Neumaier, A.: High-dimensional convex optimization via optimal affine subgradient algorithms, in ROKS workshop, 83-84 (2013) [29] 6. Ahookhosh, M., Neumaier, A.: Fast subgradient algorithms for costly convex optimization problems, Manuscript, University of Vienna, (2014) [29] 7. Andrews, H., Hunt, B.: Digital Image Restoration. Englewood Cliffs, NJ: Prentice-Hall, (1977) [11] 8. Beck, A., Teboulle, M.: Mirror descent and nonlinear projected subgradient methods for convex optimization, Operations Research Letters 31, 167–175 (2003) [2, 3] 9. Beck, A., Teboulle, M.: A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2, 183–202 (2009) [2, 12, 23] 10. Beck, A., Teboulle, M.: Fast gradient-based algorithms for constrained total variation image denoising and deblurring, IEEE Transactions on Image Processing 18(11), 2419–2434 (2009) [12] 11. Becker, S.R., Candès, E.J., Grant, M.C.: Templates for convex cone problems with applications to sparse signal recovery, Mathematical Programming Computation 3, 165–218 (2011) [2, 3, 24] 12. Ben-Tal, A., Nemirovski, A.: Non-Euclidean restricted memory level method for large-scale convex optimization, Mathematical Programming 102, 407–456 (2005) [2] 13. Bertalmio, M., Sapiro, G., Caselles, V., Ballester, C.: Image inpainting, in Proceedings of SIGGRAPH, 417–424 (2000) [12] 14. Bertero, M., Boccacci, P.: Introduction to Inverse Problems in Imaging, Bristol, U.K.: IOP, (1998) [11] 15. Bioucas-Dias, J., Figueiredo, M.: A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration, IEEE Transactions on Image Processing 16(12), 2992–3004 (2007) [12] 16. Boţ, R.I., Hendrich, C.: A double smoothing technique for solving unconstrained nondifferentiable convex optimization problems, Computational Optimization and Applications 54(2), 239–262 (2013) [23] 17. Boţ, R.I., Hendrich, C.: On the acceleration of the double smoothing technique for unconstrained convex optimization problems, Optimization, 1–24, iFirst (2012) [23] 18. Bruckstein, A.M., Donoho, D.L., Elad, M.: From sparse solutions of systems of equations to sparse modeling of signals and images, SIAM Review 51(1), 34–81 (2009) [2] 19. Boyd, S., Xiao, L., Mutapcic, A.: Subgradient methods, (2003). http://www.stanford.edu/class/ee392o/subgrad_ method.pdf [23, 24] 20. Candés, E.: Compressive sampling, in Proceedings of International Congress of Mathematics, Vol. 3, Madrid, Spain, 1433–1452 (2006) [26] 21. Candés, E., Tao, T.: Near-optimal signal recovery from random projections: universal encoding etrategies?, IEEE Transactions of Information Theory 52(12), 5406–5425 (2006) [26, 27] 22. Candés, E., Romberg, J., Tao, T.: Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Transactions of Information Theory 52(2), 489–509 (2006) [26, 27] 23. Chambolle, A.: An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision 20(1-2), 89–97 (2004) [12] 24. Chan, T.F., Shen, J., Zhou, H.M.: Total variation wavelet inpainting, Journal of Mathematical Imaging and Vision 25, 107–125 (2006) [12] 25. Chen, X., Zhou, W.: Smoothing nonlinear conjugate gradient method for image restoration using nonsmooth nonconvex minimization, SIAM Journal on Imaging Sciences 3(4), 765–790 (2010) [2] 26. Devolder, O., Glineur, F., Nesterov, Y.: First-order methods of smooth convex optimization with inexact oracle, Mathematical Programming, (2013) DOI 10.1007/s10107-013-0677-5. [3] 27. Dolan, E., Moré, J.J.: Benchmarking optimization software with performance profiles, Mathematical Programming 91, 201–213 (2002) [15] 28. Donoho, D.L.: Compressed sensing, IEEE Transactions of Information Theory 52(4), 1289–1306 (2006) [26, 27] 29. Elad, M., Milanfar, P., Rubinstein, R.: Analysis versus synthesis in signal priors, Inverse Problems 23, 947–968 (2007) [11] 30. Figueiredo, M.A.T., Nowak, R.D., Wright, S.J.: Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems, IEEE Journal of Selected Topics in Signal Processing 1(4), 586–597 (2007) [27] 31. Gonzaga, C.C., Karas, E.W.: Fine tuning Nesterov’s steepest descent algorithm for differentiable convex programming, Mathematical Programming 138, 141–166 (2013) [3] 32. Gonzaga, C.C., Karas, E.W., Rossetto, D.R.: An optimal algorithm for constrained differentiable convex optimization, Optimization-online preprint 3053, (2011) http://www.optimization-online.org/DB_FILE/2011/06/3053.pdf [3] 33. Kaufman, L., Neumaier, A.: PET regularization by envelope guided conjugate gradients, IEEE Transactions on Medical Imaging 15, 385–389 (1996) [11] Optimal subgradient algorithms with application to large-scale linear inverse problems 31 34. Kaufman, L., Neumaier, A.: Regularization of ill-posed problems by envelope guided conjugate gradients, Journal of Computational and Graphical Statistics 6(4), 451–463 (1997) [11] 35. Lan, G.: Bundle-level type methods uniformly optimal for smooth and nonsmooth convex optimization, Mathematical Programming, (2013) DOI 10.1007/s10107-013-0737-x. [2, 3] 36. Lewis, A.S., Overton, M.L.: Nonsmooth optimization via quasi-Newton methods, Mathematical Programming 141, 135-163 (2013) [23] 37. Meng, X., Chen, H.: Accelerating Nesterov’s method for strongly convex functions with Lipschitz gradient, (2011) http://arxiv.org/abs/1109.6058 [3] 38. Nemirovsky, A.S., Nesterov, Y.: Optimal methods for smooth convex minimization, Zh. Vichisl. Mat. Fiz. (In Russian) 25(3), 356–369 (1985) [3] 39. Nemirovsky, A.S., Yudin, D.B.: Problem complexity and method efficiency in optimization, Wiley, New York (1983) [2, 3, 6] 40. Nesterov, Y.: A method of solving a convex programming problem with convergence rate O(1/k2 ), Doklady AN SSSR (In Russian), 269 (1983), 543-547. English translation: Soviet Math. Dokl. 27, 372–376 (1983) [3, 23, 24] 41. Nesterov, Y.: On an approach to the construction of optimal methods of minimization of smooth convex function, Èkonom. i. Mat. Metody (In Russian) 24, 509–517 (1988) [3] 42. Nesterov, Y.: Introductory lectures on convex optimization: A basic course, Kluwer, Dordrecht (2004) [2, 3, 6, 23, 24] 43. Nesterov, Y.: Smooth minimization of non-smooth functions, Mathematical Programming 103, 127–152 (2005) [2, 3, 23, 24] 44. Nesterov, Y.: Smoothing technique and its applications in semidefinite optimization, Mathematical Programming 110, 245–259 (2007) [3] 45. Nesterov, Y.: Excessive gap technique in nonsmooth convex minimization, SIAM Journal on Optimization 16(1), 235–249 (2005) [2, 3] 46. Nesterov, Y.: Gradient methods for minimizing composite objective function, Mathematical Programming, 140 (2013), 125–161. [2, 3, 23, 24] 47. Nesterov, Y.: Primal-dual subgradient methods for convex problems, Mathematical Programming 120, 221–259 (2009) [2] 48. Nesterov, Y.: Universal gradient methods for convex optimization problems, Technical Report, Center for Operations Research and Econometrics (CORE), Catholic University of Louvain (UCL), (2013) [2, 3, 23, 24] 49. Neumaier, A.: OSGA: a fast subgradient algorithm with optimal complexity, Manuscript, University of Vienna, (2014) http://arxiv.org/abs/1402.1125 [1, 3, 4, 5, 6, 7, 8, 9] 50. Neumaier, A.: Solving ill-conditioned and singular linear systems: A tutorial on regularization, SIAM Review 40(3), 636–666 (1998) [2, 11] 51. Osher, S., Rudin, L., Fatemi, E.: Nonlinear total variation basednoise removal algorithms, Physica 60, 259–268 (1992) [11, 12] 52. Parikh, N., Boyd, S.: Proximal Algorithms, Foundations and Trends in Optimization 1(3), 123–231 (2013) [2, 23] 53. Tikhonov, A.N.: Solution of incorrectly formulated problems and the regularization method, Soviet Mathematics Doklady 4, 1035–1038 (1963) [2] 54. Tseng, P.: On accelerated proximal gradient methods for convex-concave optimization, Manuscript (2008) http:// pages.cs.wisc.edu/~brecht/cs726docs/Tseng.APG.pdf [3] 55. Wright, S.J., Nowak, R.D., Figueiredo, M.A.T.: Sparse reconstruction by separable approximation, IEEE Transactions on Signal Processing 57(7), 2479–2493 (2009) [15]