Catalyst Acceleration for Gradient-Based Non-Convex Optimization

We introduce a generic scheme to solve nonconvex optimization problems using gradient-based algorithms originally designed for minimizing convex functions. Even though these methods may originally require convexity to operate, the proposed approach allows one to use them on weakly convex objectives, which covers a large class of non-convex functions typically appearing in machine learning and signal processing. In general, the scheme is guaranteed to produce a stationary point with a worst-case efficiency typical of first-order methods, and when the objective turns out to be convex, it automatically accelerates in the sense of Nesterov and achieves near-optimal convergence rate in function values. These properties are achieved without assuming any knowledge about the convexity of the objective, by automatically adapting to the unknown weak convexity constant. We conclude the paper by showing promising experimental results obtained by applying our approach to incremental algorithms such as SVRG and SAGA for sparse matrix factorization and for learning neural networks.

Introduction
A large number of machine learning problems are formulated as the minimization of a regularized empirical risk where each function f i : R p → R is smooth and the regularization ψ : R p → R may be nonsmooth. Note that the composite setting also encompasses constrained minimization. Whereas a significant amount of work has been devoted to this composite setting for convex problems, leading in particular to fast incremental algorithms [see, e.g., 8,19,36,38,39], the question of minimizing efficiently (1) when the functions f i and ψ may be nonconvex is still largely open today. Yet, nonconvex problems in machine learning are of high interest. For instance, the variable x may represent the parameters of a neural network, where each term f i (x) measures the fit between x and a data point indexed by i, or (1) may correspond to a nonconvex matrix factorization problem (see Section 5).
Besides, even when the data-fitting functions f i are convex, it is also typical to consider nonconvex regularization functions ψ, for example for feature selection in signal processing [16].
Several pioneering works attempted to transfer ideas from the convex world to the nonconvex one. Our paper has a similar goal and studies the extension of Nesterov's acceleration for convex problems [25] to nonconvex composite ones. Unfortunately, the concept of acceleration for nonconvex problems is unclear from a worst-case complexity point of view: gradient descent requires O(ε −2 ) iterations to guarantee a residual gradient norm smaller than ε, whereas state-of-the-art methods [e.g., 6] only achieve marginal gain with complexity O(ε −7/4 log(1/ε)) by using second-order information.
For this reason, our work fits within a broader stream of recent research on methods that do not perform worse than gradient descent in the nonconvex case (in terms of worst-case complexity), while automatically accelerating for minimizing convex functions. The hope when applying such methods to nonconvex problems is to see an acceleration in practice, by heuristically exploiting convexity that is "hidden" in the objective (for instance, local convexity near the optimum, or convexity along trajectories).
The main contribution of this paper is a generic meta-algorithm, dubbed 4WD-Catalyst, which is able to use a gradient-based optimization method M, originally designed for convex problems, and turn it into an accelerated scheme that also applies to nonconvex objective functions. The proposed 4WD-Catalyst can be seen as a 4-Wheel-Drive extension of Catalyst [21] to all optimization "terrains" (convex and nonconvex), while Catalyst was originally proposed for convex optimization. As a result, we are able to accelerate first-order methods M designed to minimize convex objectives of form (1). 27], the first accelerated method performing universally well for nonconvex and convex problems was introduced by [12]. Specifically, it addresses composite problems such as (1) with n = 1, and, provided the iterates are bounded, it performs no worse than gradient descent on nonconvex instances with complexity O(ε −2 ) on the gradient norm residual. When the problem is convex, it accelerates with complexity O(ε −2/3 ). Extensions to accelerated Gauss-Newton type methods were also recently developed by [9].

Related work. Inspired by Nesterov's acceleration method for convex optimization [see
In a follow-up work [13], a different algorithm for the same problem class is described. This scheme monotonically interlaces proximal gradient descent steps and Nesterov's extrapolation. The boundedness assumption on the iterates can then be dropped. Extensions when the gradient of ψ is only Hölder continuous can be devised.
In [20], the authors follow a similar strategy, focusing instead on convergence guarantees under the so-called Kurdyka-Łojasiewicz inequality-a property corresponding to polynomial-like growth of the function, as shown by [5]. Our scheme is in the same spirit as these papers, since it monotonically interlaces proximal-point steps (instead of proximal-gradient) and acceleration steps. A fundamental difference is that our method is generic and accommodates inexact computations, since we allow the subproblems to be approximately solved by any method we wish to accelerate.
By considering C 2 -smooth nonconvex objective functions f with Lipschitz continuous gradient ∇f and Hessian ∇ 2 f , the authors of [6] propose an algorithm with complexity O(ε −7/4 log(1/ε)), based on iteratively solving convex subproblems closely related to the original problem. It is not clear if the complexity of their algorithm improves in the convex setting. Note also that the algorithm proposed in [6] is inherently for C 2 -smooth minimization. This implies in particular that incorporating nonsmooth regularizers could be challenging.
Finally, a stochastic method related to SVRG [18] for minimizing large sums while automatically adapting to the weak convexity constant of the objective function is proposed in [1]. When the weak convexity constant is small (i.e., the function is nearly convex), the proposed method enjoys an improved efficiency estimate. This algorithm, however, does not automatically accelerate for convex problems, in the sense that the overall rate is slower than O(ε −3/2 ) in terms of target accuracy ε on the gradient norm.
Organization of the paper. Section 2 presents mathematical tools for non-convex analysis; Sections 3 and 4 introduces the main algorithm and important extensions, resp. Finally, Section 5 presents experimental results on matrix factorization and training of neural nets.

Tools for non-convex optimization
First, we will use the concept of weak convexity, which we recall below. Definition 2.1 (Weak convexity). A function f : R p → R is ρ−weakly convex if for any points x, y ∈ R p and λ ∈ [0, 1], the approximate secant inequality holds: Notice that ρ-weak convexity with ρ = 0 is exactly the definition of a convex function. An elementary algebraic manipulation shows that f is ρ-weakly convex if and only if the function Useful characterizations of ρ-weakly convex functions rely on differential properties. Since the functions we consider in the paper are nonsmooth, we use a generalized derivative construction. We mostly follow the standard monograph on the subject [35].
Thus, a vector v lies in ∂f (x) whenever the linear function y → f (x) + v T (y − x) is a lower-model of f , up to first-order around x. In particular, the subdifferential ∂f (x) of a differentiable function f is the singleton {∇f (x)}, while for a convex function f it coincides with the subdifferential in the sense of convex analysis. It is useful to keep in mind that the sum rule, ∂(f + g)(x) = ∂f (x) + ∇g(x), holds for any differentiable function g.
We are interested in deriving complexity bounds on the number of iterations required by a method M to guarantee dist 0, ∂f (x) ≤ ε .
In our convergence analysis, we will also use the following differential characterization of ρ-weakly convex functions, which generalize classical properties of convex functions. A proof may be found in Theorem 12.17 of [35].

Theorem 2.3 (Differential characterization).
For any lower-semicontinuous function f : R p → R, the following properties are equivalent: holds for all x, y ∈ R p and v ∈ ∂f (x).

(hypo-monotonicity) The inequality
holds for all x, y ∈ R p and v ∈ ∂f (x), w ∈ ∂f (y).
Weakly convex functions have appeared in a wide variety of contexts, and under different names. Some notable examples are globally lower-C 2 [34], prox-regular [31], proximally smooth functions [7], and those functions whose epigraph has positive reach [10].

Convergence rates in composite minimization
We now briefly discuss convergence rates, which are typically given in different forms in the convex and non-convex cases. For that purpose, we consider a strongly-convex composite minimization problem where f 0 : R p → R is µ-strongly convex and smooth with L-Lipschitz continuous gradient ∇f 0 , and ψ : R p → R is a closed convex function with a computable proximal map prox βψ (y) := argmin Let x * be the minimizer of h and h * be the minimal value of h. In general, there are three types of measures of optimality that one can monitor: Since h is strongly convex, all three are equivalent in terms of convergence rates if one can take an extra prox-gradient step: [x] L := prox ψ/L (x − L −1 ∇f 0 (x)).
To see this, define the displacement vector g L (x) := L(x − [x] L ), and notice the inclusion g L (x) ∈ ∂h([x] L ). In particular g L (x) = 0 if and only if x is the minimizer. The following set of inequalities follow directly from Theorem 2.2.7 in [27]: Thus, an estimate of any one of the four quantities x − x * , h(x) − h * , g L (x) , or dist(0, ∂h(x)) directly implies an estimate of the other three evaluated either at x or at [x] L .

4WD-Catalyst acceleration: the basic scheme
We now present a generic acceleration scheme for minimizing where the function f is ρ-weakly convex. The proposed approach extends the Catalyst acceleration of [21], and comes with a simplified convergence analysis.
We shall make several technical assumptions on f . We assume that: } is bounded and we let B be its diameter.
(A2) We have available an upper bound C > ρ.
Note that, for the composite setting, the reader should just imagine setting C = 2L, where L is the Lipschitz constant of the smooth summand.
For the sake of clarity, in this section, we shall focus on a basic version of our approach. Indeed we assume that the regularization constant we choose, κ, is fixed throughout the iterations and κ > ρ. We shall show in Section 4 how to modify the basic scheme to proceed without the assumption κ > ρ. We outline the basic scheme in Algorithm 1.
Linear convergence and acceleration. As in the regular Catalyst algorithm of [21], our scheme wraps in an outer loop a minimization algorithm M used in an inner loop. Note that, as in the regular Catalyst, the minimization algorithm M is assumed to enjoy a linear convergence rate for minimizing strongly convex functions, yet is used as a sub-routine, even though the objective may be non-convex.
In each iteration, 4WD-Catalyst forms minimization subproblems of the form where f κ (· ; y) is strongly convex as soon as κ > ρ. We will precisely specify the function f κ (· ; y) shortly; for now, the reader should think of f κ (· ; y) as a model of our target function f centered around y, which is easier to minimize. Denote f * κ (y) = inf x∈R p f κ (x; y). We say that the minimization algorithm M, generating a sequence of iterates (z t ) t≥0 to solve the problem (P), has a linear convergence rate if there exist A κ ≥ 0 and τ κ ∈ (0, 1) such that where y is an arbitrary centering point, z 0 is the initial point, and ∂f κ (z t , y) denotes the subdifferential of f κ (·, y) evaluated at z t . In other words, we assume that the minimization algorithm M initiated at z 0 for (P) returns a point z satisfying dist(0, ∂f (z; y)) ≤ ε after at most T iterations where Remark 3.1. The linear convergence we assume here for M differs from the one considered by [21], which was given in terms of function values. However, if the problem is a composite one, both points of view are near-equivalent, as discussed in Section 2.1 and the precise relationship given in Appendix B.
Model (P) and approximate Moreau envelope. The model functions f κ (·; y) we use are simple quadratic perturbations of the target function f : The parameter κ ≥ 0 controls the amount of quadratic perturbation, and we call y the prox-center. The perturbed function f κ (x; y) has an interesting property. Indeed the function f κ (·; y) is (κ − ρ)-strongly convex for any κ > ρ. Therefore, if we have at hand an algorithm that enjoys linear convergence guarantees for strongly convex objectives, we can use it to minimize f κ (x; y) in place of f . We shall analyze this general strategy in Section 4. The simplest and most famous instance of this strategy is the proximal-point algorithm [23,33]. When applied to the function f , the algorithm proceeds as: for κ > ρ. The proximal point algorithm actually performs gradient descent on the function φ defined by φ κ (y) := inf x∈R p f κ (x; y) known as the Moreau-Yosida envelope of f .
The Moreau-Yosida envelope enjoys well-known properties, which we briefly recall here [see, e.g. 2,4,17]. For κ > ρ, the function φ κ is differentiable with Lipschitz continuous gradient, all minimizers of φ κ are minimizers of f , and the equation ∇φ κ (x t ) = κ(x t+1 − x t ) holds. Hence the proximal point algorithm is the gradient descent recursion where f * is any lower bound on f . Indeed, first-order optimality conditions yield, Note that the proximal point algorithm, as stated, is a rather conceptual algorithm since the proximal subproblems cannot be solved exactly in general. Instead, one must apply an iterative method to solve the subproblems. However, great care must be taken in analyzing the effect of inexactness associated with approximately solving these subproblems [see 15,21].
Accelerated first-order methods We shall also build upon Nesterov's accelerated gradient method [24,27] for minimizing a smooth convex function f with κ-Lipschitz continuous gradients. The method works as follows: where the weights are given by Note that x k is defined through the linearization of f at y k . If in the definition of x k , the linearization f (y k ) + ∇f (y k ) T (x − y k ) is replaced by f (x), then this is exactly the proximal point update x k = argmin x f κ (x; y k ). This accelerated proximal point method was investigated by [15], and has the convergence guarantee min i=1,...,k An important observation is that this conceptual algorithm works only for convex functions. Moreover, the proximal point steps cannot be exactly computed in general. See [11,21] for a discussion in the context of machine learning.
Main contribution We introduce in this paper a general approach, called 4WD-Catalyst, which allows to get the same worst-case guarantees as the proximal point algorithm for nonconvex problems and the same worst-case guarantees as the accelerated method for convex ones. The approach can be seen as an extension of the Catalyst acceleration scheme [11,21] to possibly nonconvex problems. The proposed approach is described in extenso in Algorithm 1. At each iteration, the scheme performs an approximate step of the proximal point method (7) yieldingx and then an approximate step of the accelerated method (8)-(11) yieldingx. The iterate x for the next round of the algorithm is declared in (12) to be the one achieving the minimal function value betweenx andx. This strategy is inspired by [13], who interlace a proximal-gradient method with an accelerated gradient method [26].
However, our scheme has two important differences with [13]. First, we proposed a generic scheme, that is able to accelerate a given optimization method M. Second, we rely a different type of acceleration scheme than the one considered in [13]. Indeed, the one considered by [13] turns out to be inappropriate for our purpose. Specifically, a naive proximal-point extension of their algorithm would require the strong convexity constants of the proximal subproblems to tend to zero, which does not allow us to control the complexity of solving these subproblems. Carefully controlling the complexity of solving the subproblems is in fact one of our key technical contributions.
Two steps. The proposed acceleration scheme builds two main sequences of iterates (x k ) k and (x k ) k , and two auxiliary sequences (y k ) k and (v k ) k . At each iteration, we proceed in two steps, corresponding to approximate minimization using M of the two different model functions f κ (·; x k−1 ) and f κ (·; y k ).
Proximal point step. Using M, we first perform an approximate minimization of a quadratic perturbation of f centered around x k−1 , namely an inexact proximal point step: We build where α k and v k−1 are defined by simple formulas we will detail in the next paragraph. The point y k is then used as the prox-center in f κ (x; y k ) for the following step.
Accelerated proximal point step. Using M, we next perform an approximate minimization of a quadratic perturbation of f now centered around y k Algorithm 1 4WD-Catalyst input Fix a point x 0 ∈ dom f , a real number κ > 0, two sequences of positive real numbers {ε k } k≥1 and {δ k } k≥1 , and an optimization method M. initialization: 2. Set 3. Choosex k using M such that 6. Choose x k to be any point satisfying until the stopping criterion dist 0, ∂f (x k ) < ε [27].
Picking the best. At the end of iteration k, we have at hand two iterates, resp.x k andx k . Following [12], we simply choose the best of the two in terms of their objective values, that is we choose The proposed scheme blends the two steps in a synergistic way, allowing us to recover the near-optimal rates of convergence in both worlds: convex and non-convex. Intuitively, whenx k is chosen, it means that Nesterov's extrapolation step "fails" to accelerate convergence.
Convergence analysis. We now present the theoretical properties of Algorithm 1 for minimization algorithms M with deterministic convergence rates of the form (3). When the rate is given as an expectation, see [21] for more details. For space limitation reasons, we shall mainly sketch the proof here, and defer the full proofs to the Appendix. Denote f * := liminf k→∞ f (x k ). We first establish the convergence rates for the sequence of best iterates generated by 4WD-Catalyst in terms of the gradient norm. We present the convergence rate for the outer loop, without taking into account the complexity required to solve the subproblems.
Define f * := liminf k→∞ f (x k ). Then for any N ≥ 1, the iterates generated by Algorithm 1 satisfy This result is based on the proof of inexact proximal algorithms [4,14,21], (see Appendix A), and is stated under the assumption that κ is set to a value that is larger than the true weak convexity constant ρ. The price we pay for the lack of knowledge of ρ and inexactness is an additional term This assumption can in fact be lifted, and we present an algorithm that automatically adapts κ in Section 4.
The choice of the sequences (ε k ) k and (δ k ) k is completely straightforward here and does not require any additional parameter tuning. The sequences (ε k ) k and (δ k ) k determine the stopping criterion of the minimization algorithm M respectively for the two major steps of 4WD-Catalyst. When the function turns out to be convex, we can characterize the convergence rate both in terms of decrease of the distance of subgradients to the origin and the decrease in function values: and where x * is any minimizer of f .
. Assume that we add to (7) the descent con- ; then, it is easy to see that the sequence {f (x k )} k≥1 is monotonically decreasing. This condition is reasonable if the iterative method M is initialized with x k−1 for minimizing f κ (.; x k−1 ). This assumption, as we later will see, simplifies (24) in the proof of Theorem 3.2 and Theorem 3.3.
We use a proof technique inspired by [3,37] combined with a technical argument from [21], which allows us to get an upper bound on function values This theorem establishes a rate of O(N −2 ) for suboptimality in function value and convergence in O(N −3/2 ) for the minimal norm of subgradients. The first rate is optimal in terms of information-based complexity for the minimization of a convex composite function [27,29]. The second can be improved to O(N −2 log(N )) through a regularization technique, if one has available an estimate on the distance of the initial point to an optimal solution [28].
A major deficiency of this analysis is that it does not take into account the number of iterations performed by M for approximately solving the subproblems. Therefore, a more careful analysis is needed to obtain the global computational complexity of the full algorithm. It turns out, for technical reasons which we discuss below, to obtain global complexity bounds, we must modify the algorithm. We now present an adaptive and more flexible version of 4WD-Catalyst, called 4WD-Catalyst-Automatic, which we first analyze without taking into account the complexity of the subproblems, and then, we establish its global complexity guarantees.

The 4WD-Catalyst-Automatic algorithm
In this section, we work towards understanding the global efficiency of Algorithm 1, which takes into account the cost of approximately solving the proximal subproblems (7) and (9), which leads to the intuition behind 4WD-Catalyst-Automatic; then, we present the complexity analysis in details.
We now state our assumptions on the algorithm M used to approximately solve the subproblems. We assume that for any κ > ρ, the method M applied to min x f κ (x; y) drives subgradients along the iterates to zero at a linear rate as in (3).
Linear convergence of M: For any κ > ρ, there exist A κ ≥ 0 and τ κ ∈ (0, 1) such that 1. For any prox-center y ∈ R p and initial point z 0 ∈ R p , the iterates {z t } t≥1 generated by M on the problem min z f κ (z; y) satisfy (3), namely

There exists a constant T 0 satisfying
3. The rates τ κ are increasing in κ.
In particular, the first two assumptions together imply Thus, we have in essence eliminated the dependence of A κ on κ. All interesting algorithms (e.g., gradient descent, SVRG) satisfy the three properties above with T 0 = 2.
Inexact subproblems and non-convexity. Consider a minimization method M satisfying the assumptions above. The number of iterations of M needed to evaluatex k in Line 7 of Algorithm 1, when initialized at x k−1 , is at most Similarly, the number of iterations of M needed to evaluatex k in Line 9 of Algorithm 1, when initialized at In order to get global complexity guarantees, we must upper-bound as a function of the iterate k. Herein lies a major challenge. While we can deduce reasonable upper bounds on f κ (x k−1 , y k ) − f * κ (y k ) when f is convex, it is much more challenging when f may be nonconvex. Furthermore, for the algorithm to be practical, we would like to try setting κ < ρ, so that the total number of the outer iterations is kept as small as possible. However, as we do so, we also lose all control on the number of iterations required by M to computex k andx k , since the linear convergence guarantees are no longer applicable for the possibly nonconvex subproblem.
Adaptive versions of 4WD-Catalyst. To solve the previous theoretical and practical challenges, we introduce an adaptive variant, called 4WD-Catalyst-Automatic, of 4WD-Catalyst. The adaptive algorithm is presented in Algorithm 2. There are two main modifications to Algorithm 1.
The first idea consists of using a pre-defined number of iterations for solving (9), corresponding to the maximal number of iterations of M that we are willing to tolerate for evaluatingx k . The main motivation is three-fold: (i) the complexity of such a strategy is known in advance; (ii) in the convex case, we will show that as soon as this number is sufficiently large, we will get acceleration because convexity allows us to control the terms f κ (x k−1 ; y k ) − f * κ (y k ); (iii) for nonconvex problems, we remark from the proof of Theorem 3.2 that the convergence guarantees are obtained regardless of the quality ofx k for solving (9). In other words, the conclusions of Theorem 3.2 are valid even when the iteratesx k are not δ k -approximate stationary solutions for non-convex functions. This allows us to choose the same strategy for computingx k in the convex and non-convex case-that is, we use the pre-defined number of iterations of M and we also choose a smoothing parameter κ cvx for the corresponding subproblem (18) that is set as if the objective was convex (meaning as in [21]).
The second idea lies in an adaptive procedure, Auto-adapt; see Algorithm 3. This procedure uses a pre-defined number of iterations of M to compute the iteratesx k while automatically choosing the smoothing parameter κ. We show that as soon as the number of iterations allocated to M for computingx k is large enough, acceleration is preserved in the convex case, while keeping the conclusions of Theorem 3.2 valid. Under our assumptions on M (see Section 4), we show that the procedure Autoadapt(x, κ, ε, T ) indeed terminates regardless of the number of iterations T chosen with the assumption that τ κ → 1 as κ → ∞ (see Lemma B.1). If the later is not the case, we can still choose T sufficiently large so that Auto-adapt(x, κ, ε, T ) terminates (See Proposition 4.3). When a formula for τ κ is explicitly given as a function of κ, the number of calls made by Auto-adapt(x, κ, ε, T ) to M is easy to compute.
Convergence analysis. Let us next postulate that T and S are chosen large enough to guarantee that x k andx k are ε k and δ k -stationary for the corresponding subproblems, and see how the outer algorithm complexity resembles the guarantees of Theorem 3.2 and Theorem 3.3. The main technical difference is Algorithm 2 4WD-Catalyst-Automatic input Fix a point x 0 ∈ dom f , real numbers κ min , κ cvx > 0 and T, S > 0, sequences of real numbers {ε k } k≥1 > 0 and {δ k } k≥1 > 0, and an optimization method M. initialization: repeat for k = 1, 2, . . .

Compute
3. Apply method M initialized at x k−1 for S log 1/δ 2 k iterations: 6. Choose x k to be any point satisfying until the stopping criterion dist 0, ∂f (x k ) < ε end while output z T . that we can no longer assume κ > ρ, which was used to show the key inequality (22). This is where the descent condition f κ (x k ; x k−1 ) ≤ f κ (x k−1 ; x k−1 ) comes in.
From now on, denote by κ k , the final constant κ used by Algorithm 3 in iteration k. Set κ max := max k≥1 κ k .
Then for any N ≥ 1, the iterates generated by Algorithm 2 satisfy, min j=1,...,N If in addition the function f is convex and S is chosen so thatx k is a δ k -stationary for f κcvx (·; y k ), then min j=1,...,2N where x * is any minimizer of the function f .

Inner-loop Complexity
In light of Theorem 4.1, we must now understand how to choose T and S as small as possible, while guaranteeing thatx k andx k are ε k and δ k -stationary solutions for f κ k (·; x k−1 ) and f κcvx (·; y k ), respectively, and the descent inequality f κ k (x k ; x k−1 ) ≤ f κ k (x k−1 ; x k−1 ) holds for each k. To simplify notation, we will often drop the subscript k in κ k ; the index will be clear from context. Henceforth, we assume that f satisfies assumptions (A1) and (A2). The quantities T and S must be set beforehand without knowing the true value of the weak convexity constant ρ. However, in general, we do know an upper bound C on ρ; for instance, in the composite setting ρ ≤ L, where L is the Lipschitz constant of the smooth part's gradient. To help determine the appropriate values of T and S, recall by (16) the number of iterations of the method M initialized at z 0 to produce an ε-stationary point for f κ (·; y) is at most There are two components controlling this value: (i) the term A κ (f κ (z 0 ; y) − f * κ (y)) and (ii) the term 1 τκ . Additionally, we must ensure that the iteratesx k satisfy a functional decrease, Note the values f (x k ) always decrease in Algorithm 2; hence we initialize the method M at x k−1 to produce bothx k andx k . Theorem 4.2 shows how to choose T and S and thus obtain the complexities of solving the subproblems. The idea is simply to choose T such that Algorithm 3 terminates when κ is not much greater than C + ρ for any index k. Let B be the diameter of the sublevel set, {x : f (x) ≤ f (x 0 )}.

and we choose T and S in Algorithm 2 to be the smallest numbers satisfying
for all k. In particular, Then κ max ≤ 6C and the following hold for any index k ≥ 1: The difficulty in proving Theorem 4.2 is twofold. First, Auto-adapt must terminate for all sufficiently large choices of T even when τ κ → 1 as κ → ∞ while simultaneously producing a point satisfying both ε k -stationarity of f κ (·; x k−1 ) and the descent condition. Secondly, the growth of the y k can not be controlled when the function f is non-convex. Consequently, the δ k -stationarity of f κmax (·, y k ) can only be guaranteed provided the function f is convex. We provide a work-around to control the complexity of computing the iteratesx k . Appendix C is devoted to proving Theorem 4.2, but we outline below the general procedure and state the two main propositions (see Proposition 4.3 and Proposition 4.6).
We summarize the proof of Theorem 4.2 as followed: 1. When κ > ρ + C, we compute the number of iterations of M to produce both an ε k -stationary solution and a point that satisfies the decrease condition. Such a point will becomex k .

2.
When the function f is convex, we compute the number of iterations of M to produce a point which satisfies the δ k -stationary condition. Such a point will become the pointx k .
3. We compute the number of times we must double κ min until it becomes larger than ρ + C. Thus eventually the condition κ > ρ + C will occur. 4. We always set the number of iterations of M to producex k andx k as in Step 1 and Step 2, respectively, regardless of whether f κ (·; x k ) is convex or f is convex.
The next proposition shows that Auto-adapt terminates with a suitable choice forx k after T log(1/ε k ) number of iterations.

Proposition 4.3 (Inner complexity forx k ). Suppose κ > ρ + C. If the method M is initialized at a point
x ∈ R p and is applied to the problem for at least the number of iterations Remark 4.4 ( Determining τ C ). Proposition 4.3 shows that by choosing T in such a way, Algorithm 3 terminates with an ε k -stationary point that also decreases the initial function value. The choice of T depends on τ C , which in turn relies on the method M. We see that τ C can be lower bounded by a value independent of the weak convexity of ρ. For instance, when ρ ≤ L, we can set C to equal 2L. If the method M is gradient-descent, then 1/τ Similarly, if the method M is SVRG then 1/τ C = n + C C−ρ ≤ n + 2. In most cases, the value of τ C can be estimated using known quantities of the problem and without the knowledge of ρ.
Remark 4.5. The constant A from our assumptions on M item (2) might depend on the weak convexity constant ρ. In the cases of SVRG and gradient descent, using T 0 = 2, the choice of A is simply ρ. We can apply our upper bound on ρ so that A is independent of both κ and ρ.
Under the additional assumption that the function f is convex, we produce a δ k -stationarity point for f κcvx (·; y k ) when the number of iterations, S, is chosen sufficently large.
where we define B = diam({x : f (x) ≤ f (x 0 )}), then, the outputz S k =x k is an δ k -stationary solution for all k ≥ 1.
We can now derive global complexity bounds by combining Theorem 4.1 and Theorem 4.2, and a good choice for the constant κ cvx .
Global complexity bounds for 4WD-Catalyst-Automatic : Set s = 1 and choose T and S according to Theorem 4.2. We letÕ hide universal constants and logarithmic dependencies in A, B, C, ε, κ min , κ cvx , τ C , τ κcvx and f (x 0 ) − f * . Then the following are true.
1. Algorithm 2 will generate a point x satisfying dist 0, ∂f (x) ≤ ε after at most iterations of the method M.
2. If f is convex, then Algorithm 2 will generate a point x satisfying dist 0, ∂f (x) ≤ ε after at mostÕ iterations of the method M.
3. If f is convex, then Algorithm 2 will generate a point x satisfying f (x) − f * ≤ ε after at most iterations of the method M.
Choosing C and κ cvx : We consider now choosing the parameters for the composite setting. The convergence rate of Algorithm 2 is driven by C (the upper bound on ρ) and κ cvx . For typical methods, the term τ κ takes the form τ κ = ϕ κ/(κ + L) for some function ϕ. Therefore, good choices for κ cvx and C are κ cvx , C ≈ L, since then the rate τ κ does not depend on L. For example, consider the composite setting. Then gradient descent satisfies τ κ = κ/(κ + L) while for SVRG (modulo taking expectations) we have τ κ = κ/ (L + (n + 1)κ). Therefore these choices for κ cvx and C, up to log terms, give the best known complexity bounds for methods in both convex and nonconvex setting. We present a summary of common choices for the parameters when minimizing the composite form, Variable Description Initialization  Figure 1: Dictionary learning experiment. We plot the function value of SVRG with different stepsizes η, where η is given by the user, before and after 4WD-Catalyst. For unaccelerated SVRG, one iteration corresponds to one pass over the data. The curves in red correspond to the largest step-sizes for which SVRG is stable. From left to right, we vary the size of the dataset. Dotted lines correspond to accelerated algorithms, related to the plain line of the same color representing the unaccelerated counterpart.

Experiments
As a proof of concept, we demonstrate the 4WD-Catalyst-Automatic acceleration for two non-convex machine learning problems. In all experiments, we consider the acceleration of the SVRG method originally designed for convex problems [39], and then studied for nonconvex ones by [32]. Unfortunately, theoretical step-sizes for SVRG in non-convex settings often depend on quantities that are unknown (e.g., Lipschitz constant L), and are pessimistic, leading to slow convergence. Therefore, we chose to run SVRG with various step-sizes that are potentially larger than the theoretical ones for non-convex problems, and report the corresponding performance when SVRG is (i) applied directly to the objective function, and (ii) used as a sub-routine M of the proposed approach 4WD-Catalyst-Automatic. In each case, this means guessing a Lipschitz constant L, using the step-size η = 1/L in SVRG, and using κ min = κ cvx = 2L/n, as advised by [21]. The number of iterations of M is always set to S = T = n, meaning performing one pass over the data for solving each sub-problem.
We will now present the two experiments, conducted on sparse matrix factorization and neural networks, and discuss the results, including cases where the acceleration is significant in practice and also cases where it fails.
Sparse matrix factorization a.k.a. dictionary learning. The sparse matrix factorization problem consists of representing a large dataset X = [x 1 , · · · , x n ] ∈ R m×n as a product X ≈ DA, where D in R m×p is called a dictionary, and A in R p×n is a sparse matrix. The classical formulation [see 22] is the non-convex problem where A = [α 1 · · · α n ] carries the decomposition coefficients of signals x 1 · · · x n , ψ is a sparsityinducing regularization and C is chosen as the set of matrices whose columns are in the ℓ 2 -ball. An equivalent point of view is the finite-sum problem min D∈C When we take elastic-net regularization ψ(α) = µ 2 α 2 + λ α 1 of [40], f i (D) is strongly convex, making the problem appropriate for the SVRG algorithm. We report the corresponding results in Figure 1,

Number of iterations Function value
Neural network, dataset covtype, n=100000 Figure 2: Simple experiment with a two-layer neural network on subsets of two datasets alpha (top) and covtype (bottom). The experimental setup is the same as in Figure 1. From left to right, we vary the size of the dataset's subset from n = 1 000 to n = 100 000.
when learning a dictionary in R m×p with p = 256 elements, on a set of whitened normalized image patches of size m = 8 × 8. Parameters are standard ones in this literature [22]-that is, a small value µ = 1e − 5, and λ = 0.25, leading to sparse matrices A (on average ≈ 4 non-zero coefficients per column of A).
Neural networks. We consider simple binary classification problems for learning neural networks. Assume that we are given a training set {a i , b i } n i=1 , where the variables b i in {−1, +1} represent class labels, and a i in R p are feature vectors. The estimator of a label class is now given by a two-layer neural networkb = sign(w ⊤ 2 σ(W ⊤ 1 a)), where W 1 in R p×d represents the weights of a hidden layer with d neurons, w 2 in R d carries the weight of the network's second layer, and σ(u) = log(1 + e u ) is a nonlinear function, applied pointwise to its arguments. We use the logistic loss to fit the estimators to the true labels and report our experimental results on two datasets alpha and covtype in Figure 2. The weights of the network are randomly initialized and we fix the number of hidden neurons to d = 100. Note that all parameter choices of the algorithm are also summarized at the end of Section 4.
Remark on the computational complexity. Since our scheme solves two sub-problems per iteration, the complexity per iteration is twice that of the unaccelerated SVRG. On the other hand, we have also noticed that in most experiments where acceleration was successful,x k was almost always preferred tō x k , meaning that half of the computation was in fact not used by the algorithm.
Experimental conclusions. When the dataset is small, the acceleration is significant in both experiments; when n increases, the effect is less important and our approach may even slightly hurt the convergence speed (see Figure 1 on the right). An intuitive explanation may be extrapolated from the properties of SVRG for minimizing convex functions. The complexity of SVRG for strongly-convex problems is indeed in O((n + Q) log(1/ε)), where Q is the condition number of f [39]. Assume that one may increase n without changing Q, then the complexity becomes O(n log(1/ε)) as soon as n ≫ Q, which is already optimal, preventing further acceleration. A similar phenomenon is also reported for Catalyst [21].
The second conclusion is particularly interesting since it suggests that Nesterov's extrapolation, or the Moreau-Yosida smoothing, may be helpful to to escape bad stationary points, e.g., saddle-points. The

Number of iterations Log of Subgradient Norm
Matrix factorization, n=100000 catalyst svrg η=0.1 svrg η=0.1 catalyst svrg η=0.01 svrg η=0.01 catalyst svrg η=0.001 svrg η=0.001 Figure 3: Dictionary learning experiment. We plot the norm of gradients of SVRG with different stepsizes η, where η is given by the user, before and after 4WD-Catalyst. For unaccelerated SVRG, one iteration corresponds to one pass over the data. The curves in red correspond to the largest step-sizes for which SVRG is stable. From left to right, we vary the size of the dataset. Dotted lines correspond to accelerated algorithms, related to the plain line of the same color representing the unaccelerated counterpart.   Figure 3. From left to right, we vary the size of the dataset's subset from n = 1 000 to n = 100 000.
success of our scheme for the neural network experiment is indeed intriguing, but most striking are the convergence curves in terms of residual gradient norm, which we see from Figure 4. With the random initialization, the norm of the gradient at the initial point is indeed small, and we see that our acceleration scheme seems very helpful to escape the (bad) initial point.

Conclusion.
In this work, we propose a new generic meta-algorithm, 4WD-Catalyst, and its variant, 4WD-Catalyst-Automatic, for solving (1) which performs no worse than gradient descent in the nonconvex case (at mostÕ(1/ε) number of iterations) while automatically accelerating when the function is convex (at mostÕ(1/ε 2/3 ) number of iterations). The proposed meta-algorithm wraps a gradient-based method M, designed for convex problems, with a proximal point method and Nesterov's extrapolation and turns it into an accelerated scheme that also works for nonconvex problems. Unlike previous work for solving (1) with nonconvex functions, 4WD-Catalyst-Automaticdoes not require Hessians, instead we approximately compute the weak convexity constant using a doubling procedure while simultaneously finding the smoothing parameter κ. This together with the algorithm working for both convex and nonconvex objective functions should heuristically speed up convergence. As proof of concept, we conducted two experiments using our approach applied to SVRG: neural networks and matrix factorization.
The results were promising showing acceleration occurs particularly when the dataset is small, but more interestingly the Nesterov's extrapolation may help escape bad stationary points. In this appendix, Section A is devoted to proving the outer convergence guarantees for 4WD-Catalyst as stated in Theorem 3.2 and Theorem 3.3. The beginning of Section B illustrates in the composite setting the relationship between convergence guarantees for M stated in terms of function values versus subdifferentials. Next in Lemma B.1, we state and prove that the adaptation to 4WD-Catalyst, Auto-adapt, terminates provided an additional assumption on M. Then Section B finishes with the proof for the outer convergence result for 4WD-Catalyst-Automatic (Theorem 4.1). The inner convergence results (Theorem 4.2) as well as the supporting Proposition 4.3 and Proposition 4.6 are proved in the final section of the appendix, Section C.

A Theoretical analysis of the basic algorithm
We present here proofs of the theoretical results of the paper. Althroughout the proofs, we shall work under the Assumptions on f stated in Section 3 and the Assumptions on M stated in Section 4.

A.1 Convergence guarantee of 4WD-Catalyst
In Theorem 3.2 and Theorem 3.3 under an appropriate tolerance policy on the proximal subproblems (7) and (9), 4WD-Catalyst performs no worse than an exact proximal point method in general, while automatically accelerating when f is convex. For this, we need the following observations. Lemma A.1 (Growth of (α k )). Suppose the sequence {α k } k≥1 is produced by Algorithm 1. Then, the following bounds hold for all k ≥ 1: Proof. This result is noted without proof in a remark of [37]. For completeness, we give below a simple proof using induction. Clearly, the statement holds for k = 1. Assume the inequality on the right-hand side holds for k. By using the induction hypothesis, we get as claimed and the expression for α k+1 is given by explicitly solving (11). To show the lower bound, we note that for all k ≥ 1, we have Using the established upper bound α k ≤ 2 k+1 yields The result follows.
Next we establish convergence guarantees of Theorem 3.2 and Theorem 3.3 for 4WD-Catalyst .
Proof of Theorem 3.3 and Theorem 3.3. The proof of (13) follows the analysis of the inexact proximal point method [21,14,4].
Completing the square on the right-hand-side, we obtain Using (12), we get an approximate descent guarantee on f along the iterates Summing up from j = 1, . . . , N , we get min j=1,...,N Using Lemma A.2, we finally get min j=1,...,N Using the identity N j=1 ε 2 j ≤ ∞ 0 8 (1+x) 3 dx = 4, the first result is shown. Next, suppose the function f is convex. Our analysis is similar to that of [37,3]. Fix ξ k ∈ ∂f κ (x k ; y k ) with ξ k < δ k . For any x ∈ R n , Equation (12) and the strong convexity of the function f κ (·; y k ) yields We substitute x = α k x * + (1 − α k )x k−1 where x * is any minimizer of f . Using the convexity of the function f and Equations (8) and (10), we deduce Subtracting f (x * ) from both sides, we have Iterating N times using the equality 1−α k and we obtain the overall estimate Using the lower bound on α j 's from Lemma A.1, we note N j=1 δ j The function value result follows by applying the upper bound in Lemma A.1. Finally, let us establish the bound (14). The idea is to combine the descent guarantee (22) with the rapid convergence of f (x k ) to f (x * ) in (15). Using (22) and summing up over N + 1, . . . , 2N , we have min j=N +1,...,2N The last inequality follows from 2N Combining this inequality with (15), we obtain the gradient bound.

B Analysis of 4WD-Catalyst-Automatic and Auto-adapt
Linear convergence interlude. Our assumption on the linear rate of convergence of M (see (3)) may look strange at first sight. Nevertheless, most linearly convergent first-order methods M for composite minimization either already satisfy this assumption or can be made to satisfy it by introducing an extra prox-gradient step. To see this, recall the convex composite minimization problem from Section 2.1 where 1. f 0 : R p → R is convex and C 1 -smooth with the gradient ∇f 0 that is L-Lipschitz, 2. ψ : R p → R is a closed convex function with a computable proximal map See [30] for a survey of proximal maps. Typical linear convergence guarantees of an optimization algorithm assert existence of constants A ∈ R and τ ∈ (0, 1) satisfying for each t = 0, 1, 2, . . . , ∞. To bring such convergence guarantees into the desired form (3), define the prox-gradient step [z] L := prox ψ/L (z − L −1 ∇f 0 (z)), and the displacement vector and notice the inclusion g L (z) ∈ ∂h([z] L ). The following inequality follows from [29]: Thus, the linear rate of convergence (29) implies which is exactly in the desired form (3).

B.1 Convergence analysis of the adaptive algorithm: 4WD-Catalyst-Automatic
First, under some reasonable assumptions on the method M (see Section 4), the sub-method Autoadapt terminates.
We prove the main result, Theorem 4.1, for 4WD-Catalyst-Automatic.
Proof of Theorem 4.1. The only difference in the proof, in comparison to the proof of Theorem 3.2, is that we cannot immediately conclude inequality (22). Instead, fix ζ k ∈ ∂f κ (x k ; x k−1 ) with ζ k < ε k and set κ = κ k . The decrease condition then implies With this inequality replacing the bound (22), the rest of the proof is identical to that of Theorem 3. Proof of Theorem 4.2. We consider two cases: (i) the function f is non-convex and (ii) the function f is convex. First, we consider the non-convex setting. To producex k , the method M is called number of times. This follows from Proposition 4.3 and Lemma C.1. Since the iteratex k is the output of Algorithm 3, then κ > ρ + C so Proposition 4.3 states thatx k is an ε k -stationary solution and satisfies f κ (x k ; x k−1 ) ≤ f κ (x k−1 ; x k−1 ), which are the two conditions required to apply Theorem 4.1.
To producex k , the method M is called S log(1/δ 2 k ) number of times. (Note: the proof of Theorem 4.1 does not needx k to be a δ k -stationary solution in the non-convex case).
Next, suppose the function f is convex. As before, to producex k the method M is called (32) times. (Note: the proof of Theorem 4.1 does not requirex k to be a ε k -stationary solution or satisfy the decrease condition in the convex case). To producex k , the method M is called S log(1/δ 2 k ) number of times. By Proposition 4.6, the iteratex k satisfies the δ k -stationarity condition; a key ingredient in the proof of Theorem 4.1. Plugging in the values of δ k and ε k , both the non-convex and convex cases follow.
C.1 Inner complexity forx k : proof of Proposition 4.3 Next, we supply the proof of Proposition 4.3 which shows that by choosing κ large enough, Algorithm 3 terminates.
Proof of Proposition 4.3. When the number of iterations t > T 0 , there exists a constant A > 0 independent of κ such that dist 2 0, ∂f κ (z t ; x) ≤ A(1 − τ κ ) t−T 0 (f κ (x; x) − f * κ (x)). By the linear convergence of M and the (κ-ρ)-strong convexity of f κ (·; x), we have as τ κ is an increasing function of κ. Since M ≥ 1 ). Moreover, we see where the last inequality uses both the fact that f (x) are decreasing and f * κ (x) ≥ f * . Again because Next, we compute the maximum number of times we must double κ until κ > ρ + C. Since κ is doubled (Algorithm 3) and T is chosen as in Proposition 4.3 , the maximum the value κ, κ max , takes is 2(ρ + C) ≤ 6C.
C.2 Inner complexity forx k : proof of Proposition 4.6 In this section, we prove Proposition 4.6, an inner complexity result for the iteratesx k . Recall that the inner-complexity analysis forx k is important only when f is convex (see Section 4). Therefore, we assume throughout this section that the function f is convex. If Proposition 4.3 holds, then the descent condition f (x k ) ≤ f (x k−1 ) is satisfied; hence we know x k ∈ {x : f (x) ≤ f (x 0 )} for all k ≥ 0. Let B be the diameter of the sublevel set {x : f (x) ≤ f (x 0 )}. The following lemmas are useful in proving Proposition 4.6.
Lemma C.2 (Controlling the iteratesx k ). Let x * be any minimizer of the function f . If the iteratesx k are generated by Algorithm 2 and ρ = 0 and the iteratesx j are δ j -stationary solutions for all 1 ≤ j ≤ k, then for any index k ≥ 1 and constant κ > 0, these iterates satisfy the inequalities: Next, suppose (33) holds. By (35), we deduce for k + 1: x k − y k+1 ≤ α k+1 x k − x k−1 + x k − y k + x k−1 − y k .
Using (33), we conclude the result hold for k + 1.


extfieldsofstudy:
Medicine
Computer Science

text:
Multi-categorical deep learning neural network to classify retinal images: A pilot study employing small database

Deep learning emerges as a powerful tool for analyzing medical images. Retinal disease detection by using computer-aided diagnosis from fundus image has emerged as a new method. We applied deep learning convolutional neural network by using MatConvNet for an automated detection of multiple retinal diseases with fundus photographs involved in STructured Analysis of the REtina (STARE) database. Dataset was built by expanding data on 10 categories, including normal retina and nine retinal diseases. The optimal outcomes were acquired by using a random forest transfer learning based on VGG-19 architecture. The classification results depended greatly on the number of categories. As the number of categories increased, the performance of deep learning models was diminished. When all 10 categories were included, we obtained results with an accuracy of 30.5%, relative classifier information (RCI) of 0.052, and Cohen’s kappa of 0.224. Considering three integrated normal, background diabetic retinopathy, and dry age-related macular degeneration, the multi-categorical classifier showed accuracy of 72.8%, 0.283 RCI, and 0.577 kappa. In addition, several ensemble classifiers enhanced the multi-categorical classification performance. The transfer learning incorporated with ensemble classifier of clustering and voting approach presented the best performance with accuracy of 36.7%, 0.053 RCI, and 0.225 kappa in the 10 retinal diseases classification problem. First, due to the small size of datasets, the deep learning techniques in this study were ineffective to be applied in clinics where numerous patients suffering from various types of retinal disorders visit for diagnosis and treatment. Second, we found that the transfer learning incorporated with ensemble classifiers can improve the classification performance in order to detect multi-categorical retinal diseases. Further studies should confirm the effectiveness of algorithms with large datasets obtained from hospitals.

Introduction
Retina is a photosensitive layer of optic nerve tissue lining in the inner surface of the eyeball. Retinal damages due to various diseases can eventually lead to irreversible vision loss. As population aging has emerged as a major demographic trend worldwide, patients suffering from chorioretinal diseases such as age-related macular degeneration (AMD) and diabetic retinopathy (DMR) are expected to increase in the future [1]. AMD is can cause blindness [2]. DMR, which a common lifestyle disease, is also the a major cause of blindness in patients with diabetes mellitus [3]. Other retinal diseases including retinal vessel occlusion, hypertensive retinopathy, and retinitis are significant causes of vision impairment. If early diagnosis and treatment are implemented prior to the initial stage of blindness progression, visual loss can be avoided in many cases. Hence, more precise screening program is required for early treatment in high-risk group in an effort to reduce socioeconomic burdens of visual loss caused by retinal diseases. DMR screening that uses fundus photograph is universally adopted for diabetes patients. Moreover, screening such as AMD is the most appropriate approach for early intervention in the asymptomatic stage [4]. Conducting AMD screening and DMR screening is cost-effective in a public health setting [5]. However, manual analysis for multiple fundus photographs for an accurate screening requires a great deal of efforts of ophthalmologists.
Many previous studies have focused on automated detection of retinal diseases by using machine learning algorithms in order to analyze a large number of fundus photographs taken from retinal screening programs [6,7]. Various machine learning algorithms-K-nearest neighbor algorithm, Naive Bayes classifier, artificial neural network (ANN), and support vector machine (SVM)-were applied to automated retinal disease detection [8]. However, only a few studies developed machine learning models for AMD detection; whereas, most studies devoted in identifying DMR [9].
Deep learning for analyzing medical images appeared in the field of machine learning technique [10]. There were several reports on introduced ANN models to mark the difference between glaucoma and non-glaucoma [11,12]. A glaucoma research group reported visual field analysis by using deep feed-forward neural network to discover preperimetric glaucoma [13]. An automated deep convolutional neural network (CNN) was applied in the grading severity of nuclear cataract [14]. Moreover, a similar technique model identifying retinopathy of prematurity by using babies' retinal images was developed [15]. According to the recent outcome from Abramoff's research team, this learning technique demonstrated better performance in terms of automated DMR detection than previous algorithms [16]. The Google research team has introduced the advanced deep learning model capable of diagnosing DMR as well as human ophthalmologists [17]. Using similar deep learning techniques, fundus photographs and optical coherence tomography were used for analyzing the AMD [18,19]. Yet all the studies for retinal image classification selected binary classification through which "one disease versus normal" problems were settled. Although studies performed in the past released the outcome that high performances of classification in controlled experimental settings, it is practically difficult to apply the binary classification model into the real clinical setting where visiting patients suffer from various retinal diseases. Nonetheless, studies about multi-categorical classification aiming at identifying ocular diseases have been very limited.
In this study, we applied deep learning using a state-of-the-art CNN for fundus photography analysis in multi-categorical disease settings. This paper articulates a pilot study designed for deep learning assessment on multi-categorical classification by using small open retinal image database.

Methods
We utilized publically available retinal image database at STructured Analysis of the REtina (STARE) project (available at http://www.ces.clemson.edu/~ahoover/stare) [20] in order to evaluate a multi-categorical deep learning model. Fig 1 shows a flow diagram of the proposed system. The experimental process complied with the Declaration of Helsinki. The Ethics committee approval was not required, because researchers instead used public database. The STARE project aimed to develop an image-understanding system to distinguish retinal diseases from fundus images. Database is comprised of retinal color images acquired by a TRV-50 fundus camera (Topcon Corp., Tokyo, Japan) at a 35 degrees field with a resolution of 605 x 700 pixels. The database contains 397 images in 14 disease categories including emboli, branch retinal artery occlusion (BRAO), cilio-retinal artery occlusion, branch retinal vein occlusion (BRVO), central retinal vein occlusion (CRVO), hemi-CRVO, background diabetic retinopathy (BDR), proliferative diabetic retinopathy (PDR), arteriosclerotic retinopathy, hypertensive retinopathy, Coat's disease, macroaneurism, choroidal neo-vascularization(CNV), and the other retinal status. However, there are several problems on original categorization. First, categorization was unequally distributed. For example, dataset includes only a single image of cilio-retinal artery occlusion, while more than 60 images were included into the BDR. Second, important disease groups, such as dry AMD (drusens in macula), wet AMD, central retinal arterial occlusion (CRAO), and retinitis were not classified. Therefore, two ophthalmologists (T.K.Y and J.G.S) reviewed all images and assigned new categories. Several categories were removed due to extremely small number of images such as emboli without vessel occlusion (6 images), cilioretinal artery occlusion (1 image), arteriosclerotic retinopathy without DMR (6 images) and macroaneurism (8 images). BRVO, CRVO, and hemi-CRVO categories were incorporated into a single category of retinal vein occlusion (RVO) and categories of BRAO. Furthermore, CRAO was also integrated into retinal arterial occlusion (RAO). We excluded 28 low-resolution images and 69 the ambiguity of retinal images that present overlapped multiple diseases. At the last stage, we involved remaining 279 images and classified them into 10 categories including normal (25 images), BDR (63 images), PDR (17 images), dry AMD (25 images), wet AMD (48 images), RVO (38 images), RAO (12 images), hypertensive retinopathy (19 images), Coat's disease (12 images) and retinitis (20 images).
Researches who noted that deep learning in medical image analysis used trained deep CNN models from scratch. The recent development of rapid parallel solvers with GPU promoted to train huge parameters in deeper CNN models. State-of-the-art deep learning algorithms presented in the ImageNet Large Scale Visual Recognition Competition (ILSVRC) originally centered on multi-categorical (or multiclass) classification problems since data provided by the ILSVRC contains more than a million training images from 1,000 object categories [21]. We used the deep learning CNN model with 19 layers (VGG-19) along with MatConvNet (available at http://www.vlfeat.org/matconvnet). Oxford University Visual Geometry Group developed MatConvNet, an open toolbox that implements new CNN including VGG-16, and VGG-19, placing second in the ILSVRC in 2014. VGG-19 has been widely adopted to solve image classification problems as the GoogLeNet, which won ILSVRC prize in 2014. There was no significant difference in terms of performance between VGG-19 and GoogLeNet. Moreover, VGG-19 is simple and efficient for individual user [22]. VGG-19 uses 224 x 224 pixels with RGB 3-channel images as input variables. MatConvNet also provides AlexNet, a standard CNN architecture as well as a winner of ILSVRC in 2012. There are three types of layers of all CNN models in this study: convolutional (computing the output of the connected local input neurons reading patterns), max pooling (sub-sampling the inputs) and fully connected layers (allocating final scores of each class). VGG-19 is composed of 16 convolutional layers and 3 fully connected layers.
In this study, we compared and analyzed four distinct deep learning models. VGG-19 and AlexNet were operated based on a stochastic gradient descent (SGD) method. Other models adopted a transfer learning technique based on a pre-trained VGG-19 model. If the models were trained ones by using SGD, we adopted a pre-trained model as a starting point for learning the network weights. The pre-trained model, previously trained on a subset of current Ima-geNet database (provided by MatConvNet), was further processed with fundus image dataset. VGG-19 and AlexNet were algorithms applied in training by using momentum 0.9, and a fixed learning rate of 10 −6 for 50 epochs.
A machine learning model can apply accumulated knowledge to a new task domain by applying transfer learning technique [23]. We retain original 279 retinal images from 10 categories, which are relatively small to deal with a CNN with millions of parameters. Thus, they should be adjusted optimally. A transfer learning by using pre-trained CNN can contribute to avoid the problem associated with a few dataset in medicine. Previous studies suggested that CNN intermediate layer outputs can function as input features to train other classifiers and that this skill demonstrated satisfactory accomplishment performance in dealing with various problems; thus, we used a pre-trained VGG-19 model [23,24]. Although this pre-trained model was optimally programmed to identify 1000 objects in ImageNet, our hypothesis is that several pre-trained texture features might be apt for analysis on photographs. When the CNN served as a feature extractor, multiclass random forest (RF) and SVM models were trained by utilizing the 4096 input features from the last covered layer of pre-trained VGG-19 model. SVM is a universally well-established technique based on mapping data in a higher dimensional space via a kernel function and selecting the maximum-margin hyper-plane dividing training data. Multiclass SVM adopted one-vs-one design that builds up binary SVM models for all pairs of classes. A decision function of one-vs-one design assigns an instance to a class that involves countless votes [25]. RF refers to a robust and powerful multiclass classification method that promotes many classification trees from random subsets of predictors and bootstrap samples. Previous studies revealed that these two multiclass classifier are the most robust techniques that outperforms other algorithms including decision trees, k-nearest neighbor, and back-propagation neural networks [26].
We applied ensemble classifiers to the transfer learning process to enhance the performance. Ensemble methods have been proved to be a potent tool to stabilize and improve the performance of machine learning classifiers [27]. By using the disease-labeled image data and 4096 input features from the last covered layer of pre-trained VGG-19 model, we trained ensemble SVM classifiers by complying with clustering and voting approaches (iRSpot-EL) [28], K-means clustering with dynamic selection strategy (D3C) [29], multiple kernel learning [30], and AdaBoost (deep SVM) [31]. We replaced the previous ensemble classifier with preserving the structure of iRSpot-EL (modified iRSpot-EL). Previous researchers defined the distance between two classifiers C(i) and C(j) as follows: where m represents the number of training samples, d ik refers to the misclassification probability of classifier C(i) on the kth sample, and d ik Δd jk can be calculated as follows: if CðiÞ and CðjÞ incorrectly predicts the kth sample ( Based on the distance, the affinity propagation clustering algorithm was estimated. 300 different multiclass SVM models were constructed by using following parameter combinations in order to set up different multi-categorical classifiers: where OVO represents one-versus-one classifier, OVA stands for one-versus-all classifier [32], DAG notes directed acyclic graph classifier [33], and RBF_SVM refers to radial basis function SVM with a penalty parameter C and scaling factor σ. After 300 different multiclass SVM models were acquired, they were classified into seven clusters in compliance with affinity propagation clustering [34]. The ensemble process was implemented via the following fractional votes: where P i represents the probability from the classifier C(i), and F i refers to its fraction used, which was optimized on the validation sets. This ensemble multi-categorical classifier is compared to single SVM, RF, and the above ensemble methods. We also applied feature selection methods with subsets of the top-ranked 1024, 2048, 3072, and 4096 input features in order to examine how feature dimensionality and change the performance. Features were selected by using Kruskal-Walis one-way ANOVA (KW), ratio of features between-categories to withincategory sum of squares (BW) [35], and Max-Relevance-Max-Distance (MRMD) [36]. All parameters of each ensemble method were highly tuned to promote the performance. We augmented data by oversampling images with translation, rotation, brightness change, and additive Gaussian noise due to reduced size of fundus image dataset for training CNN models [37]. Data augmentation is a widely used approach to boost the generalization of deep learning models. We randomly retrieved transformed 1000 fundus images per each disease class because of the imbalance of data problems. Specifically, we obtained samples with translation from the range [-10%, +10%] of the image width, with rotation from [-15˚, +15˚], and with brightness change from a range of [-10%, +10%]. Additive Gaussian noise has a uniformly sampled sigma from [0, 0.04]. All images were organized according to the input size of the pre-trained model (224 x 224 pixels) in the course of oversampling.
The measurement of multi-categorical problems was based on the accuracy, relative classifier information (RCI), and Cohen's kappa metric [38]. Accuracy is a standard metric for evaluation of a classifier. It is defined as follows: where the element q ij refers to the number of test times and test input actually labeled C i is C j noted by the classifier, and these elements organize the confusion matrix. Although it is easy to notice the accuracy, it cannot give full accounts on the actual performance in multi-categorical problems. The RCI is an entropy-based measure applicable to multi-categorical decision problems [39]. This quantifies how much uncertainty of classification had been reduced by a machine learning classifier [25]. It is defined as follows: where log refers to natural logarithm transformation. RCI represents the performance with unbalanced classes capable of distinguishing among different misclassification distributions. Cohen's kappa is an alternative to classification rate that compensates for random hits [40]. It is defined as follows: Kappa is a standard meter for a multi-categorical problem generally applied in several fields such as brain-computer interface. Matlab 2016a (Mathworks, Natick, MA, USA) was prepared to perform the algorithms. When we trained binary classifiers, we maximized the Youden's index to select cut-off points and granted equivalent portions to sensitivity and specificity [41]. We used MedCalc 12.3 (MedCalc, Mariakerke, Belgium) for Receiver Operating Characteristic (ROC) analysis. When we conducted ROC analysis, we divided all dataset (10,000 images) into training dataset (70%) and test dataset (30%). When training process was completed, test dataset was valid to create ROC curves. This process demonstrated that our training process did not derail an ordinary deep learning study process if we obtained the similar classification performance like the previous binary classification research.
We compared four deep learning models: transfer learning with random forest based on VGG-19 structure (VGG19-TL-RF), transfer learning with Gaussian kernel SVM based on VGG-19 structure (VGG19-TL-SVM), VGG-19 and AlexNet. The 5-fold cross validation scheme validated deep learning models. When there was a slight change in the number of categories, retinal diseases were classified in accordance with its importance (judged by T.K.Y) as follows: BDR, PDR, dry AMD, wet AMD, RVO, RAO, Hypertensive retinopathy, Coat's disease, and retinitis. We adopted a grid search where ranges of parameter values were tested in order to obtain the optimal result from RF and SVM. We conducted the grid search in all dataset of 10 categories. A radial basis kernel function with a penalty parameter C of 100 and scaling factor σ of 10 determined the optimal model of SVM. 1,000 trees and five predictors for each node for RF were optimal. Since transfer learning approaches did not train the huge deep learning structure, the only a short time was needed to perform the grid search process in order to decide parameters. We did not perform left-right alignment of fundus images to identify the generalized performance of deep learning algorithms. We employed the NVIDIA GEFORCE GTX1060 3GB GPU for transfer learning and GTX980 6GB for SGD with Intel core i7 processor to train deep learning models more rapid. As categories multiplied, performance of deep learning models underperformed. When only two categories (normal and BDR) were involved in the VGG19-TL-RF, the overall classification accuracy, RCI, and kappa were 87.4%, 0.453, and 0.747, respectively; whereas, the accuracy, RCI, and kappa were 30.5%, 0.052, and 0.224, respectively when all 10 categories were included.

Results
We analyzed a detailed binary classification in order to determine categories that reduced the multi-categorical classification performance. We trained each binary classification model by using VGG19-TL-RF in all pairs of categories to examine binary discriminative powers between retinal diseases. Fig 3 presents an accuracy of 5-fold cross validation of pair-wise binary classification for all diseases. All diseases except RAO were separated from normal retina in which accuracy was over 80.0%. Discrimination between BDR and hypertensive retinopathy found the worst among all pairs (accuracy 58.4%). The accuracy of binary classifier discriminating wet AMD and Coat's disease showed 66.1%, which was lower than the mean accuracy.
We developed more general screening classification models to detect retinal abnormality. Fig 4 shows the ROC curves of each deep learning model for binary classification between normal or any disease status (normal versus abnormal). The VGG-TL-RF predicted abnormal retinal disease status (including 9 retinal diseases) with an area under the curve (AUC) of 0.903, sensitivity of 80.3%, and specificity of 85.5%. This result showed outperformance among others.
Considering the most important two disease groups including DMR and AMD, we combined images of normal, BDR, PDR, dry AMD, and wet AMD in the experiment as introduced in Table 1. VGG19-TL-RF also succeeded in all situations. For screening for early cases without progressed stages (a scenario of early disease screening), accuracy of VGG19-TL-RF using normal, BDR, and dry AMD was 72.8%, RCI of 0.283, and kappa of 0.577. When all disease categories on DMR and AMD (five categories: normal, BDR, PDR, dry AMD, and wet AMD) were included for more clinical situation, accuracy of VGG19-TL-RF was 59.1%, RCI of 0.151, and kappa of 0.485. When only DMR groups (normal, BDR, and PDR) were involved (a scenario of DMR staging), accuracy of VGG19-TL-RF found 80.8%, RCI of 0.420, and kappa of 0.711. When AMD groups (normal, dry AMD, and wet AMD) were included (a scenario of AMD staging), accuracy of VGG19-TL-RF was 77.2%, RCI of 0.371, and kappa of 0.657.
We also analyzed the performance of ensemble learning in order to improve the deep learning classifier. Modified iRSpot-EL, multiple kernel learning, D3C, and deep SVM improved the classification performance in transfer learning setting ( Table 2). In particular, the modified iRSpot-EL model (accuracy of 36.7%; RCI of 0.053; kappa of 0.225) performed better than other methods. When we reduced numbers of input features, MRMD feature selection  We divided all data set (10,000 images) into training dataset (70%) and test dataset (30%). Retinal disease status includes diabetic retinopathy, age-related macular degeneration, retinal vein occlusion, retinal artery occlusion, hypertensive retinopathy, Coat's disease, and retinitis. methods raised the accuracy (Fig 5). The remaining algorithms (KW and BW) revealed a poor performance. Although feature selection was performed, transfer learning classifiers by using all input features achieved optimal performance.

Discussion
This study is an investigation on multi-categorical deep learning algorithms for automated detection of multiple retinal diseases. Findings from this study revealed that the current deep learning algorithms were ineffective to classify multi-class retinal images from small datasets. It did not show effective and practical outcomes for a computer-aided clinical application. However, this paper suggests that the automated classifier should aim to discriminate at least normal, DMR, and AMD with multi-categorical classification due to importance of each disease in order to construct an automated retinal disease classification model by using deep learning model for the general public screening (especially for elderly). The performance of deep learning models was dropped, as categories multiplied. In fact, this result is quite natural. When facing two categories (random performance: 50%), it is fairly common to have better performance than three-category situation (random performance: 33.3%). As categories multiplied, the expected accuracy at random distribution dropped. This result corresponded with previous studies [42]. A recent study that applied the GoogLeNet architecture for classifying skin cancer demonstrated that the growing number of classes proved underperformed (with an accuracy of 72.1% in a 3-class problem and 55.4% in a 9-class problem) [43]. An elaborate research should be conducted to augment this predictive power for multi-categorical classification. General clinicians observe that fundus image analysis on multi-categorical retinal disease diagnosis is difficult. For those unfamiliar with ophthalmology, they may find it more difficult than a binary classification problem on "one disease versus normal". Although we used the new deep learning technique, diagnostic performance was inadequate to apply in clinical practice (S1 File-example of misclassification). Since every single disease has its own pathophysiological characteristics and displays retinal image under different patterns of progression, designing only one machine technique classifying multiple retinal diseases appears challenging. For example, a recent research showed that a specific method adopting Radon transform and discrete wavelet transform can detect AMD with an accuracy of 100% [44]. Therefore, it is necessary to construct disease-specific algorithms marking a distinction between retinal diseases in order to upgrade the performance of multi-categorical classification.
Recent medical images including fundus photograph were obtained with increased resolution. Our final goal is to build up a machine learning model to classify multi-category retinal diseases. The CNN method can provide high-level feature extraction and multi-category classification solution by using a huge and intensive computation system from high-resolution image. Although previous machine learning techniques reached high and successful achievements, they were difficult to generalize due to the absence of high-level abstraction [45].
However, there are also still several challenges to apply deep learning to clinical practice. One previous paper addressed ethical and political issues in terms of establishing database [46]. Collecting large-sized multiple retinal diseases data has been difficult for this reason. Another obstacle is that actual clinical issues consist of multiclass classification problems. Previous research concentrated on binary classification for retinal disease prediction. Although Google developed the deep learning model that works better than ophthalmologists, their model 'Inception-v3' based on the GoogLeNet structure was optimized to binary classification for DMR identification [17]. This model was trained by adding large image database collected only for DMR screening from diabetes patients. In fact, there are a couple of binary classification methods that showed similar performance in Google's deep learning model. A majority of previous studies with higher accuracy for automated DMR detection have been based on SVM [47,48]. RF can be used to assess automated DMR evaluation [49]. Abramoff's research group applied deep learning CNN to detect DMR and reported significantly improved performance [16].
This study suggests that further studies on automated diagnosis by using retinal image should identify multi-categorical classification. Binary classification models for DMR patients and healthy subjects were restricted under the clinical circumstance due to the prevalence of the several diseases. This is an obvious issue, yet researchers ignored this. According to previous epidemiologic studies, the prevalence rate of DMR and AMD in U.S. were 3.4% and 6.5%, respectively [50,51]. RVO can in no way be negligible either since the prevalence of RVO is reported to be about 0.5% [52]. Facing unlearned disease is more likely to occur for binary classification models. Therefore, there are problems in binary classification models for DMR in terms of the application for the clinical practice. One study suggested that researchers should consider screening for intermediate AMD and DMR simultaneously [5]. Binary classifiers cannot be applied into this approach because of results. However, there were a few studies on the development of a multi-categorical classification model in the field of ophthalmology. Most of these studies focused on the grading severity of diseases. Previously published papers claim that a concept of multi-categorical classification was applied to predict AMD progression by using SVM and RF [53]. Multiclass SVM also worked well in classification of DMR severity [54]. The latest research practiced deep learning CNN technique into a 5-class grading of DMR [55]. This study used C4.5 and a random tree method for multi-categorical classification of DMR and glaucoma [56]. However, this study misunderstood glaucoma, because glaucoma is not a retinal disease and cannot be diagnosed by only utilizing fundus photos [57]. To the best of our knowledge, there has been no research that identified multiple retinal disease groups by utilizing machine learning techniques.
Deep transfer learning model with RF performed better than other well-trained models. We used pre-trained CNN (trained using the ImageNet) in order to add data into a new feature space by spreading the fundus image database to the CNN due to lack of fundus image database available for training deep learning models. We applied traditional multi-categorical machine learning technique including RF and SVM for transfer learning. Although fundus photograph data was actually increased, fully-trained models (VGG19 and AlexNet) failed to exceed a small number of training images. Previous research revealed that transfer learning showed excellent performance comparable to fully-trained deep learning (trained with scratch) [24,58]. Although they are untrained and used fundus images, our results also revealed that transfer learning might be appropriate for multi-categorical classification of fundus images since pre-trained models included a wide range of powerful pattern extractors such as color, texture, and shape. RF was identified as the most robust machine learning classifier in previous researches [26,59]. Although classifiers were used in transfer learning, findings were consistent with previous findings. RF do not trigger problems due to over-fitting; whereas, the SVM undergoes a fine tuning process in order to avoid over-fitting and further reveal a high performance [26].
The current study has several limitations. First, a small number of images from a single study database produced fundamental limits. Many deep learning researchers come to agreement that such a small number of each category is insufficient to test the effectiveness of the proposed method. Deep learning technique generally requires more than a million samples to train without overfitting [17]. We used data augmentation and transfer learning in order to overcome this challenge. Nevertheless, the finding from performance proved to be unfulfilled. A series of trials and errors identified in this study will develop adequate methods for further studies. In addition, there was no external validation dataset to confirm the performance of classification models. If additional studies are performed, researchers should consider more retinal images from multicenter by including multiple ethnicities for a concrete training process. Second, our dataset has limited categories of retinal disease. Although we confirmed high prevalence of major retinal diseases, we missed several important retinal diseases such as retinal detachment, chorioretinal melanoma, and myopic degeneration. A multicenter scale project is needed for gathering more detailed fundus image data on rare and crucial diseases. Third, the reference diagnosis of the STARE database was often ambivalent despite comments by an ophthalmologist on all fundus images. Accurate diagnosis of retinal diseases should be verified by optical coherence tomography or fluorescein angiography [60]. However, the STARE database did not provide more detailed diagnosis. Further research should include retinal images diagnosed by clinical standard methods.
In this paper, we investigated multi-categorical classification of deep learning for automated diagnosis by using fundus photograph. Prediction models in this pilot study failed to show the advantage of using the deep learning employing multi-class retinal image datasets due to small size of datasets regarding classification performance. However, this is the primary attempt to construct deep learning models for multi-categorical classification problem of multiple retinal diseases. Importantly, diagnosis with a new deep learning technique underperformed as numbers of disease categories increased. Ensemble classifiers such as clustering and voting approach, dynamic selection, and AdaBoost could boost the classification performance for detecting multi-categorical retinal diseases. Further studies should focus on the construction of an extended prediction model with a diverse range of retinal diseases by applying multi-categorical classification techniques and sufficient amounts of datasets collected from hospitals. This study will provide proper ways to ophthalmologists who continue researching on deep learning in terms for clinical use.
Supporting information S1 File. Example of misclassification. (PDF)