{'fraction_non_alphanumeric': 0.050032182997474875, 'fraction_numerical': 0.016066742585532506, 'mean_word_length': 4.077300150829562, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust Block Coordinate Descent

In this paper we present a novel randomized block coordinate descent method for the minimization of a convex composite objective function. The method uses (approximate) partial second-order (curvature) information, so that the algorithm performance is more robust when applied to highly nonseparable or ill conditioned problems. We call the method Robust Coordinate Descent (RCD). At each iteration of RCD, a block of coordinates is sampled randomly, a quadratic model is formed about that block and the model is minimized approximately/inexactly to determine the search direction. An inexpensive line search is then employed to ensure a monotonic decrease in the objective function and acceptance of large step sizes. We prove global convergence of the RCD algorithm, and we also present several results on the local convergence of RCD for strongly convex functions. Finally, we present numerical results on large-scale problems to demonstrate the practical performance of the method.

1. Introduction. In this work we are interested in solving the following convex composite optimization problem where f (x) is a smooth convex function and Ψ(x) is a (possibly) nonsmooth, block separable, extended real valued convex function (this will be defined precisely in Section 2.5). Problems of the form of (1.1) arise in many important scientific fields, and applications include machine learning [29], regression [27] and compressed sensing [4,3,5]. Often the term f (x) is a data fidelity term, and the term Ψ(x) represents some kind of regularization. Frequently, problems of the form of (1.1) are large-scale problems, i.e., the size of N is of the order of a million or a billion. Large-scale problems impose restrictions on the types of methods that can be employed for the solution of (1.1). In particular, the methods should have low per iteration computational cost, otherwise completing even a single iteration of the method might require unreasonable time. The methods must also rely only on simple operations such as inner products or matrix vector products, and ideally, they should offer fast progress towards optimality.
First order methods, and in particular randomized coordinate descent methods, have found great success in this area because they can take advantage of the underlying problem structure (separability and block structure), and satisfy the requirements of low computational cost and low storage requirements. For example, in [20] the authors show that their randomized coordinate descent method was able to solve sparse problems with millions of variables in a reasonable amount of time.
Unfortunately, randomized coordinate descent methods have two significant drawbacks. First, due to its coordinate nature, it is efficient only on problems with high degree of separability, and performance suffers when there is a high dependency between variables. Second, as a first-order method, coordinate descent methods do not usually capture essential curvature information of the problem and have been shown to struggle on complicated sparse problems [8].
The purpose of this work is to overcome these drawbacks by equipping a randomized block coordinate descent method with approximate partial second-order information. In particular, at every iteration of RCD the direction is obtained by solving approximately a block piece-wise quadratic model, where the model includes a matrix representing approximate second order information. Then, a line search is employed in order to guarantee a monotonic decrease of the objective function.
RCD randomly selects a block of coordinates at every iteration, which is inexpensive. Although the per iteration computational cost of the method may be higher than other randomized coordinate descent methods, we show that in practise the method is more robust and the total number of iterations decreases. In particular we show that RCD is able to solve difficult problems, on which other coordinate descent method may struggle. RCD uses an inexact search direction, (the termination condition for the block piecewise quadratic subproblem is inspired by [2]), coupled with a line search to ensure a monotonic decrease in the function values, and we prove global convergence of RCD and study its local convergence properties.
1.1. Literature review. Coordinate descent methods are some of the oldest iterative methods, and they are often better known in the literature under various names such as Jacobi methods, Gauss-Seidel methods, among others. It has been observed that these methods suffer from poor practical performance, particularly on ill-conditioned problems. However, as we enter the era of big data, coordinate descent methods are coming back into favour, because of their ability to provide approximate solutions of some realistic very large/huge scale problems in a reasonable amount of time.
Currently, randomized coordinate descent methods include that of Richtàrik and Takàč [20], where the method can be applied to unconstrained convex composite optimization problems of the form (1.1). The algorithm is supported by theoretical convergence guarantees in the form of high probability iteration complexity results, and [20] also reports very impressive practical performance on highly separable large scale problems. The work has also been extended to the parallel case [21], to include acceleration techniques [7], and to include the use of inexact updates [26].
Other important works on randomized coordinate descent methods include methods for huge-scale problems [16], work in [13] that improves the complexity analysis of [21], coordinate descent methods for group lasso [19,25] and general regularizers [24,28] and coordinate descent for constrained optimization problems [14].
Unfortunately, on ill-conditioned problems, or problems that are highly nonseparable, first order methods can display very poor practical performance, and this has prompted the study of methods that employ second order information. To this end, recently there has been a flurry of research on Newton-type methods for problems of the general form (1.1), or a special case where Ψ(x) = x 1 . For example, Karimi and Vavasis [10] have developed a proximal quasi-Newton method for l 1 -regularized least squares problems, Lee, Sun and Saunders [11,12] have proposed a family of Newtontype methods for solving problems of the form (1.1) and Scheinberg and Tang [23] present iteration complexity results for a proximal Newton-type method. Moreover, the authors in [2] extended standard inexact Newton-type methods to the case of minimization of a composite objective involving a smooth convex term plus an l 1regularizer term. Finally, there exists parallel deterministic [6] and sequential active set [22] block coordinate descent methods, where the authors incorporate some block second-order information in the algorithmic process.
We believe that the works in [2,6] can lead the way to allow general randomized coordinate descent to practically incorporate second-order information and in this paper we propose a method which combines the ideas in [2,6,20].
1.2. Core ideas and Major Contributions. In this section we list several of the core ideas and major contributions of this work on randomized block coordinate descent methods. The first two points briefly describe the idea of incorporating (approximate) partial second-order (curvature) information (which have been also presented in a similar way in [6]), whilst the last three are contributions of this paper.
1. Incorporation of some second order information. RCD uses a quadratic model to determine the search direction, which incorporates a user defined positive definite matrix H (i) (x k ). If H (i) (x k ) approximates the Hessian, then second order information is incorporated into the quadratic model, and the search direction obtained by minimizing the model is an approximate Newtontype direction. We stress that H (i) (x) can change at every iteration and this is an advantage over the method in [20] where the matrix is fixed at the start of the algorithm for each block of coordinates. 2. Inexact updates. To ensure that this method is computationally practical, it is imperative that the iterates are inexpensive, and RCD achieves this through the use of inexact updates. Any algorithm can be used to approximately minimize the quadratic model. Moreover, the stopping conditions for inner solve are easy to verify; they depend upon quantities that are easy/inexpensive to obtain, or may be available as a byproduct of the inner search direction solver. 3. Blocks can vary throughout iterations. If Ψ(x) is completely separable into coordinates then we do not restrict ourselves to a fixed block structure; rather we allow the blocks of coordinates to change at any iteration. This is important because every element of the Hessian can be accessed (this is discussed further in Section 3.2.2). 4. Line search. The algorithm includes a line search step to ensure a monotonic decrease of the objective function as iterates progress. The line search is inexpensive to perform because, at each iteration, it depends on a single block of coordinates only. One of the major advantages of incorporating second-order information combined with line search is to allow in practice the selection of large step sizes (close to one). This is because unit step sizes can substantially improve the practical efficiency of a method. We prove that if f is strongly convex, then close to the optimal solution unit step sizes are selected. In fact, for all experiments that we performed, unit step sizes were accepted by line search for the majority of the iterations. 5. Convergence theory. We provide global convergence results to show that the RCD algorithm is guaranteed to converge in the limit. We also provide local convergence theory for strongly convex functions f . In particular, depending on the choice of stopping condition for the inner search direction solve and the matrix H (i) (x k ), we show that close to the optimal solution RCD has on expectation block quadratic or superlinear rate of convergence.
1.3. Format of the paper. The paper is organised as follows. In Section 2 we introduce the notation and definitions that are used throughout this paper, as well as giving several technical results. We also define the quadratic model that is used in the algorithm, prove the equivalence of some stationarity conditions for problem (1.1), and define a continuous measure of the distance of the current point from the set of solutions of (1.1). A thorough description of the RCD algorithm is presented in Section 3, including how the blocks are selected/sampled at each iteration, a description of the search direction and line search, several suggestions for the matrices H (i) (x k ), and we also present several concrete examples.
The second half of the paper is devoted to providing convergence results and numerical experiments. In Sections 4, global convergence results are presented, which do not require f to be convex. Local convergence theory for RCD is presented in Section 5. There we show that, close to optimality line search accepts unit step sizes. Moreover, if both the stopping conditions for the inner search direction solve and the matrix H (i) (x k ) are chosen appropriately, then RCD has on expectation block quadratic or superlinear rate of convergence. Finally, several numerical experiments are presented in Section 6, which show that the algorithm performs very well in practice.

Preliminaries.
In this section we introduce the notation and definitions that are used in this paper, and we also present some important technical results. Throughout the paper · ≡ ·, · and · A ≡ ·, A· , where A is a positive definite matrix. Moreover, λ min (·) and λ max (·) denote the smallest and largest eigenvalue of ·, respectively.

Subgradient and subdifferential. For a function Φ :
are called the subgradients of Φ at point x. In words, all elements defining a linear function that supports the function Φ at point x are subgradients. The set of all s at a point x is called the subdifferential of Φ and it is denoted by ∂Φ(x).

Convex conjugate and proximal mapping. For a convex function Φ :
The proximal mapping of a convex function Ψ at x is and the proximal mapping of its convex conjugate Ψ * is The following relation holds between the two proximal mappings.
From Chapter 1 of [1], we also see that prox Ψ (·) and prox Ψ * (·) are nonexpansive Finally, from Chapter 1 of [1] we have 2.4. Block decomposition of R N . Let U ∈ R N ×N be a column permutation of the N × N identity matrix and further let U = [U 1 , U 2 , . . . , U n ] be a decomposition of U into n submatrices, where U i is N × N i and n i=1 N i = N . It is clear that any vector x ∈ R N can be written uniquely as , where x (i) ∈ R Ni and block i denotes a subset of {1, 2, . . . , N }. Moreover, these vectors are given by (2.5) x (i) := U T i x. 2.5. Block decomposition of Ψ. The function Ψ : R N → R∪{+∞} is assumed to be block separable. That is, we assume that Ψ(x) can be decomposed as: where the functions Ψ i : Notice that if n = N , Ψ(x) is said to be separable (into coordinates), whereas if n < N , then Ψ(x) is said to be block separable (separable into blocks of coordinates).
The following relationship will be used repeatedly in this work: 2.6. Block Lipschitz continuity of f . Throughout the paper we assume that the gradient of f is block Lipschitz, uniformly in x. This means that, for all x ∈ R N , i ⊆ {1, 2, . . . , n} and t (i) ∈ R Ni we have where ∇ i f (x) (2.5) = U T i ∇f (x). An important consequence of (2.8) is the following standard inequality [15, p.57]: For fixed x ∈ R N , we define a piecewise quadratic approximation of F around the point (x + t) ∈ R N as follows: and H (i) (x) ∈ R Ni×Ni is any positive definite matrix, which possibly depends on x.
Notice that Q(x; 0) = F (x) and that Q i (x, t (i) ) is the quadratic model for block i.

Stationarity conditions.
The following theorem gives the equivalence of some stationarity conditions of problem (1.1).
Theorem 2.2. The following are equivalent first order optimality conditions of problem (1.1).
Proof. It is easy to see that (i) are first-order optimality conditions of problem (1.1), which can be obtained by using the definition of subgradient. It is trivial to show that (i) ⇐⇒ (ii). By Lemma 2.1, we have that (iii) ⇐⇒ (iv). We now show that (iii) ⇐⇒ (ii). We rewrite (ii) as which is satisfied if and only if (iv) holds, hence, if and only if (iii) holds.
Let us define the continuous function where β is a positive constant, which is used in the local convergence analysis (Section 5). By Theorem 2.2, the points that satisfy g(x; 0) = ∇f (x)+ 1 β prox (βΨ) * (x − β∇f (x)) = 0 are stationary points for problem (1.1). Hence, g(x; 0) is a continuous measure of the distance from the set of stationary points of problem (1.1).
Furthermore, let us define which will be used as a continuous measure for the distance from stationarity of the block piecewise quadratic function Q i (x k ; t (i) ).
3. The Algorithm. In this section we present the Robust Coordinate Descent (RCD) algorithm for solving problems of the form (1.1). There are three key steps in the algorithm: (step 4) the coordinates are sampled randomly; (step 5) the quadratic model (2.11) is solved approximately until the stopping conditions (3.2) are satisfied to give a search direction; (step 6) a line search is performed to find a step size that ensures a sufficient reduction in the objective value. Once these key steps have been performed, the current point x k is updated to give a new point x k+1 , and the process is repeated.
The following assumption is used in RCD. The reason this assumption is used will be made clear in Section 3.1.
Assumption 3.1. The block decomposition of R N used within RCD, and the associated probability distribution, adhere to the block structure of Ψ(x).
We now present pseudocode for the algorithm, while a thorough description of each of the key steps in the algorithm will follow in the rest of this section.

Algorithm 1 Robust Coordinate Descent (RCD)
1: Input Choose x 0 ∈ R N , θ ∈ (0, 1/2) and β > 0. 2: Initialize a decomposition of R N and a probability distribution following Assumption 3.1 Sample a block of coordinates i with probability p i > 0.

6:
Perform a backtracking line search along the direction t . One of the crucial ideas of this algorithm is that the block of coordinates to be updated at each iteration is chosen randomly. This allows the coordinates to be selected very quickly. In this section we explain in detail, how the blocks are selected/sampled at each iteration. We also give examples of how coordinates can be randomly sampled such that Assumption 3.1 is satisfied.
3.1.1. Ψ is block separable with n < N . When Ψ has a fixed block structure (i.e., n < N ), the block decomposition of R N (via the matrix U = [U 1 , . . . , U n ]) described in Section 2.4 is fixed at the start of the algorithm to coincide with the block structure of Ψ, and does not change as iterations progress. There are several ways to initialize a sampling scheme to use in RCD that follow Assumption 3.1.
1. Fix the n blocks of coordinates according to the decomposition of R N defined by U , and in the algorithm, select each block of coordinates with some probability p i . (e.g., uniform probabilities p i = 1/n > 0 for all i = 1, . . . , n). 2. Perform (single) coordinate descent, where at each iteration of RCD, the coordinate i is selected with some probability p i (e.g., uniform probabilities p i = 1/N for all i). 3. Perform block coordinate descent, where each block of coordinates has cardinality N min := min{N 1 , . . . , N n }. The restriction is that, at any iteration k, the sampled coordinates forming block i, must all belong to the same block of N j coordinates defined by submatrix U j . (i.e., Assumption 3.1 is satisfied because the decomposition of U is obeyed.) Recall that Ψ is separable into n blocks. Let the total number of subblocks be l(> n), where we assume that each coordinate 1, . . . , N appears in at least one of the l blocks. Then each subblock is selected with probability p i .
When Ψ is separable into coordinates, we have complete control over the indices that are updated at each iteration.
Let τ denote the block size (number of coordinates that are updated at any iteration k), where 1 ≤ τ ≤ N . Note that there are N C τ subsets 1 of τ coordinates that can be made from the set {1, . . . , N }. At any iteration k of RCD, a subset of coordinates i k with |i k | = τ is sampled with some probability p i (e.g., uniform probabilities p i = 1/ N C τ > 0 for all i). Note that in practice, one never explicitly forms the N C τ different blocks in order to randomly pick one with some probability p i . Instead, τ coordinates are sampled randomly without replacement.

The search direction and Hessian approximation (Step 5).
In this section we describe how RCD determines the search direction. In particular, RCD forms a quadratic model for block i, and minimizes the model approximately until the stopping conditions (3.2) are satisfied, giving an 'inexact' search direction.
We also describe the importance of the choice of matrix H, which is an approximate second order information term. From now on, we will often use the shorthand H 3.2.1. The search direction. At each iteration the update/search direction is found as follows. The subproblem (3.1), (where Q i (x k ; t (i) ) is defined in (2.11)) is approximately solved, and the search direction t (i) k is accepted when the stopping conditions (3.2) are satisfied, for some η i k ∈ [0, 1). Notice that Hence, from (3.5), the stopping conditions (3.2) depend on block i only, and are therefore inexpensive to verify, meaning that they are implementable.

Remark 3.2.
(i) At some iteration k, it is possible that g i (x k ; 0) = 0. In this case, it is easy to verify that the optimal solution of subproblem (ii) Notice that, unless at optimality (i.e., g(x k ; 0) = 0), there will always be blocks i such that g i (x k ; 0) = 0, which implies that t We now provide a few suggestions for the choice of H (i) k . (This list is not intended to be exhaustive.) Notice that in each case there is a trade off between a matrix that is inexpensive to work with, and one that is a more accurate representation of the true block Hessian.
1. Clearly, the simplest option is to set H (i) k = I for all i and k. In this case no second order information is employed by the method.

A second option is to let H
can be computed and stored at the start of the algorithm and elements can be accessed throughout the algorithm as necessary. This is very effective if diag(∇ 2 f (x)) is a good approximation to ∇ 2 f (x).

A third option is to let H
k is a principal minor of the Hessian). In this case, H (i) k provides the most accurate second order information, but it is (potentially) more computationally expensive to work with. In practice the matrix ∇ 2 i f (x k ) is used in a matrix-free way and is not explicitly stored. For example, there may be an analytic formula for performing matrixvector products with ∇ 2 i f (x k ), or techniques from automatic differentiation could be employed, see [18,Section 7]. 4. Another option is to use a quasi-Newton type approach where H (i) k is an approximation to ∇ 2 i f (x k ) based on the limited-memory BFGS update scheme, see [18,Section 8]. This approach might be more suitable in cases that the problem is not very ill-conditioned and additionally performing matrix-vector products with If any of the matrices above are not positive definite, then they can be altered to make them so. For example, if H (i) k is diagonal, any zero that appears on the diagonal can be replaced with a positive number. Moreover, if ∇ 2 i f (x k ) is not positive definite, a multiple of the identity can be added to it.

An advantage of the RCD algorithm (if Option 3 is used for H
is that all elements of the Hessian can be accessed. This is because the blocks of coordinates can change at every iteration, and so too can matrix H   1) is not guaranteed. To this end, we include a line search step in our algorithm in order to guarantee monotonic decrease of function F . Essentially, the line search guarantees the sufficient decrease of F at every iteration, where sufficient decrease is measured by the loss function (3.4).
In particular, for fixed θ ∈ (0, 1/2), we require that for some α ∈ (0, 1], (3.3) is satisfied. (In Lemma 4.3 we prove that there exists a subinterval (0,α] of (0, 1] in (3.6) which shows that the calculation of the right hand side of (3.3) only depends upon block i, so it is inexpensive. Moreover, the line search condition (Step 5) involves the difference between function values F (x k ) − F (x k + αU i t (i) k ). Fortunately, while function values can be expensive to compute, the difference in the objective value between iterates need not be (this is discussed in more detail in Section 3.4).
k ) as a function of α only depends on block i, hence, it is inexpensive. Moreover, in some cases some of the quantities in (3.7) are already needed in the computation of the search direction t, so regarding the line search step, they essentially come "for free". where a T j is the jth row of a matrix A ∈ R m×n and b j is the jth component of vector b ∈ R m . As before, we need to evaluate (3.7). Let us split calculation of is more expensive because is depends upon the logarithm. In this case, one can calculate f (x 0 ) once at the beginning of the algorithm and then update f (x k + αU i t (i) k ) ∀k ≥ 1 less expensively. In particular, let us assume that the following terms: At the end of the process, f (x 1 ) and e −bj a T j x1 ∀j will be given for free, and the same process can be followed for the calculation of f ( 4. Global convergence theory without convexity of f . In this section we provide global convergence theory for the RCD algorithm. Note that we do not assume that f is convex. Throughout this section we denote H   Before proving global convergence of RCD, we present several technical results. The following lemma shows that if t (i) k is nonzero, then F is decreased. Lemma 4.3. Let Assumptions 4.1 and 4.2 hold. Let θ ∈ (0, 1/2) and let x k and i be generated by RCD. Then Step 6 of RCD will accept a step-size α that satisfies Furthermore, Proof. The proof closely follows that of [2, Theorem 3.1]. From (3.2), Rearranging gives By Assumption 4.2, for α ∈ (0, 1), we have Adding Ψ(x k ) to both sides of the above and rearranging gives By convexity of Ψ(x) we have that Then From the previous observe that if α satisfies 0 ≤ α ≤ (1 − θ) λi Li , then α also satisfies the backtracking line search step of RCD. Suppose that any α that is rejected by the line search is halved for the next line search trial. Then, it is guaranteed that the α that is accepted satisfies (4.2).
Since the line search condition (3.3) is guaranteed to be satisfied for some step size α, from (3.3) and convexity of we obtain F ( ). Using (4.5) and (4.2) in the last inequality we obtain (4.3). The following lemma bounds the norm of the direction t Proof. This proof closely follows that of [2, Theorem 3.1]. Using the reverse triangular inequality and the fact that prox Ψ * i (·) is nonexpansive, we have that Rearranging gives (4.8), and combining (4.3) and (4.8) gives (4.9).
We now have all the tools to prove global convergence of RCD. Proof. Taking the expectation of (4.3) gives From Remark 3.2(ii) it is guaranteed that t Note that λ 2 i /L i > 0 is constant for every i, so (4.10) implies that t (i) k → 0 for all i as k → ∞. Combining this with Lemma 4.4 gives that g i (x k ; 0) → 0 for all i as k → ∞. Since this holds for an arbitrary block i, the result follows.

Local convergence theory.
In this section we present local convergence theory for RCD. First we discuss some common assumptions that are needed. Throughout the section we set H denotes the principal minor of the Hessian ∇ 2 f (x) with row (equivalently column) indices in the subset i.
Assumption 5.1. Function f is strongly convex with strong convexity parameter µ f > 0.
By continuity of f we have that ∇ 2 f (x) is symmetric, and by strong convexity of f we have that µ f I ∇ 2 f (x).
The next theorem explains that, if the Hessian H ≡ ∇ 2 f (x) is positive definite, then every principal submatrix of it is also positive definite.
Theorem 5.2 (Theorem 4.3.15 in [9]). Let A ∈ R N ×N be a Hermitian matrix, let r be an integer with 1 ≤ r ≤ N , and let A r denote any r × r principal submatrix of A (obtained by deleting N − r rows and the corresponding columns from A). For each integer k such that 1 ≤ k ≤ r we have λ k (A) ≤ λ k (A r ) ≤ λ k+N −r (A), where λ k (·) denotes the kth eigenvalue of matrix ·, and the eigenvalues are ordered λ 1 ≤ . . . , λ N . Corollary 5.
3. If f is strongly convex with strong convexity parameter µ f > 0, then Assumption 5.4. We assume that the blocks of the Hessian of f are Lipschitz continuous. This means that for all x ∈ R N , i ⊆ {1, 2, . . . , n} and t (i) ∈ R Ni we have The following theorem is used to show that the backtracking line search accepts units step sizes close to optimality for any block i. Similarly to [2], in order to prove the previous statement we have to impose sufficient decrease of the quadratic model (3.1) at every iteration. This means that the inexactness conditions in (3.2) are replaced by where ξ ∈ (θ, 1/2) and η i k ∈ [0, 1). Notice that, substituting the equality in (4.4) into (5.3), gives Proof. The proof closely follows that of [12,Lemma 3.3]. Using Lipschitz continuity of H Using this with (5.4) in (5.5) we get ), which implies that RCD accepts a step α = 1. Corollary 5.6. By Theorem 4.5, t (i) k → 0 ∀i as k → ∞. Thus, there will be a region close to the optimal solution x * such that t (i) k ≤ 3λ i (ξ − θ)/M i for all i and for all k. Hence, by Theorem 5.5, in this region, the backtracking line search step in RCD accepts unit step sizes for any i.
The following assumption is mild since it is guaranteed to be satisfied by Corollary 5.6.
Assumption 5.7. Iteration x k is close to the optimal solution x * of (1.1) such that unit step sizes are accepted by the backtracking line search algorithm of RCD.
The next lemma is a technical result that will be used in Theorem 5.9. Lemma 5.8. Let Assumptions 4.2 and 5.1 hold. Let L max = max i {L 1 , . . . , L n }. Let β < 1/L max in (2.13). Then g i (x; t (i) ) inherits strong monotonicity of ∇ i f (x) ∀i: Proof. The proof is the same as [12, Lemma 3.9], but restricted to the ith block, so is omitted.
In the following theorem we demonstrate that RCD has on expectation block quadratic or superlinear local rate of convergence. k , η i k = min{1/2, g i (x k ; 0) }, and β < 1/L max in (2.13). Then, g i (x k ; 0) has a quadratic rate of convergence in expectation: where c is a positive constant. If η i k → 0 for k → ∞, then g i (x k ; 0) has a superlinear rate of convergence in expectation: Proof. For a given x k we define k . Using the Fundamental Theorem of Calculus (F.T.o.C.), we have Also, from the definition of a derivative we have Now, adding and subtracting 1 k )) from (5.7), followed by taking norms, applying the triangle inequality, and using (5.8), gives By Assumption 5.7, RCD accepts unit step sizes. Hence, setting δ = 1 in (5.9) gives Using the same trick as before with the F.T.o.C. we have the bound Then, using Cauchy-Schwarz we have By the triangular inequality and stopping conditions (3.2) we have Replacing (5.11) in (5.10) we have Moreover, by setting η i k = min{1/2, g i (x k ; 0) } we obtain Rearranging and taking expectation gives The right hand side of (5.14) is constant for all i, which implies quadratic convergence of g i (x k ; 0) in expectation. Moreover, if η i k → 0 as k → ∞, then from (5.12), Theorem 4.5 and Lemma 4.4, which implies superlinear convergence of g i (x k ; 0) in expectation.
6. Numerical Experiments. In this section we examine the performance of two versions of RCD and two versions of a Uniform Coordinate Descent method (UCDC) [20] on two common optimization problems. The first problem is an 1regularized least squares problem of the form (1.1) with where c > 0, x ∈ R N , A ∈ R m×N and b ∈ R m . The second problem is an 1 -regularized logistic regression problems of the form (1.1) with where c > 0, a j ∈ R N ∀j = 1, 2, . . . , m are the training samples and b j ∈ {−1, +1} are the corresponding labels. For (6.1) a synthetic sparse large scale experiment is performed and for (6.2) we compare the methods on two real world large scale problems from machine learning. Notice that for both (6.1) and (6.2), Ψ(x) = c x 1 , which is fully separable into coordinates. This means that, for RCD, we have complete control over the block decomposition, and the indices making up each block can change at every iteration.
All algorithms are coded in MATLAB, and for fairness, MATLAB is limited to a single computational thread for each test run. All experiments are performed on a Dell PowerEdge R920 running Redhat Enterprise Linux with four Intel Xeon E7-4830 v2 2.2GHz, 20M Cache, 7.2 GT/s QPI, Turbo (4x10Cores).
6.1. Implementations of RCD and UCDC. In this section we discuss some details of the implementations of methods RCD and UCDC.
6.1.1. RCD. For the RCD method, we fix the size of blocks τ > 1 (to be given in the numerical experiments subsections), and at every iteration of RCD, τ coordinates are sampled uniformly at random without replacement.
We implement two versions of RCD, which we denote by RCD v.1 and RCD v.2. The two versions only differ in how matrix H (i) k is chosen. In particular, for RCD v.1 we set H (i) ) for all i and k. In this case subproblem (3.1) is separable and it has a closed form solution is the well-known soft-thresholding operator which is applied component wise when u and v are vectors. Notice that since the subproblem is solved exactly there is no need to verify the stopping conditions (3.2). For RCD v.2, we set where ρ > 0 guarantees that H (i) k is positive definite for all i, k. Hence, the subproblem (3.1) is well defined. The larger ρ is the smaller the condition number of matrix H k becomes, hence, the faster subproblem (3.1) will be solved by an iterative solver. However, we do not want ρ to dominate matrix H k because the essential second order information from ∇ 2 f (x k ) will be lost.
In this setting of matrix H (i) k we solve subproblems (3.1) iteratively using an Orthant Wise Limited-memory Quasi-Newton (OWL) method, which can be downloaded from http://www.di.ens.fr/~mschmidt/Software/L1General.html. We chose OWL because it has been shown in [2] to result in a robust and efficient deterministic version of RCD, i.e. τ = N (one block of size N ). Note that we never explicitly form matrix H (i) k , we only perform matrix-vector products with it in a matrix-free manner.
6.1.2. UCDC. We also implement two versions of a uniform coordinate descent method as it is described in Algorithm 2 in [20]. For both versions the size τ of the blocks and the decomposition of R N into N/τ blocks are fixed a-priori and all blocks are selected by UCDC with uniform probability. We compare two versions of UCDC, denoted by UCDC v.1 and UCDC v.2 respectively. For UCDC v.1 we set τ = 1 and for UCDC v.2 we set τ > 1 (the exact τ is given later).
One of the key ingredients of UCDC are the block Lipschitz constants, which are explicitly required in the algorithm. For single coordinate blocks, the Lipschitz constants can be computed with relative ease. However, for blocks of size greater than 1, the block Lipschitz constants can be far more expensive to compute. (For example, for problem (6.1), the block Lipschitz constants correspond to the maximum eigenvalue of A T i A i , where A i := AU i .) For this reason, we do not compute the actual block Lipschitz constants, rather, we use an overapproximation.
To this end, let L j > 0 ∀j = 1, 2, · · · , N denote the coordinate Lipschitz constants of function f . Then the direction t (i) at every iteration is obtained by solving exactly subproblem (3.1) with using operator (6.3). Notice that for problem (6.1), (6.5) is equivalent to H is an overapproximation of the maximum eigenvalue of A T i A i . Moreover, notice that Algorithm 2 in [20] is a special case of RCD where the subproblem (3.1) is solved exactly and line search is unnecessary. This is because, by setting H (i) k as in (6.5), then subproblem (3.1) is an over estimator of function F along block coordinate direction t (i) (for details we refer the reader to [20]).

Termination Criteria and Parameter
Tuning. The only termination criteria that RCD and UCDC should have are maximum number of iterations or maximum running time. This is because using subgradients as a measure of distance from optimality or any other operation of similar cost are considered as expensive tasks for large scale problems. In our experiments RCD and UCDC are terminated when their running time exceeds the maximum allowed running time. Furthermore, for RCD we set parameter η (i) k in (3.2) equal to 0.9 ∀i, k and ρ = 10 −6 in (6.4). The maximum number of backtracking line search iterations is set to 10 and θ = 10 −3 . For UCDC the coordinate Lipschitz constants L j ∀j are calculated once at the beginning of the algorithm and this task is included in the overall running time. Finally, all methods are initialized with the zero solution.
6.3. 1 -Regularized Least Squares. In this subsection we present the performance of RCD and UCDC on the 1 -regularized least squares problem (6.1). For this problem the data A and b were synthetically constructed using a generator proposed in [17,Section 6], and we set c = 1. The advantage of this generator is that it produces data A and b with a known minimizer x * . We slightly modified the generator so that we could control the density of A.
The dimensions of the problem are N = 2 21 and m = N/4 and the generated matrix A is full rank (with at least one non zero component per column) and a density of ≈ 10 −4 mN . The optimal solution is set to have 0.01N non zero components with values uniformly at random in the interval [−1, 1]. For UCDC, the coordinate Lipschitz constants are L j := A j 2 2 j = 1, 2, · · · , N, and for RCD v.1, RCD v.2 and UCDC v.2, we set τ = 0.01N .
The result of this experiment is shown in Figure 6.1. In this figure notice that all methods were terminated after 10 4 seconds. For practical purposes, for UCDC v.1 results are shown every 10 4 iterations. For all other methods results are shown after the first iteration takes place and then at every iteration. Observe in sub Figure  6.1a that block methods RCD v.1, RCD v.2 and UCDC v.2 performed fewer iterations compared to the single coordinate UCDC v.1. This is due to much larger per iteration computational complexity of the former methods compared to the latter. RCD v.2 despite its larger per iteration computational complexity among all methods it was the only one that solved the problem to higher accuracy within the required maximum time. Moreover, observe in sub Figure 6.1b that for purely practical purposes it might be better to have a combination of methods RCD v.1 and v.2. The former could be used at the beginning of the process while the latter could be used at later stages in order to guarantee robustness and speed closer to the optimal solution. Finally, it is important to mention that on this problem for both RCD versions unit step sizes α were accepted by the backtracking line search for a major part of the process. Hence, backtracking line search was inexpensive.
6.4. 1 -Regularized Logistic Regression. In this section we present the performance of RCD and UCDC on the 1 -regularized logistic regression problems (6.2).  is not included in running time of the methods. Fig.1a shows how the objective function F (x) decreases as a function of the number of iterations. Fig.1b shows how the objective function F (x) decreases as a function of wall-clock time measured in seconds.
Such problems are important in machine learning and are used for training a linear classifier x ∈ R N that separates input data into two distinct clusters, for example, see [29] for further details.
We present the performance of the methods on two sparse large scale data sets. Problem details are given in Table 6.1, where A ∈ R m×N is a matrix whose rows are training samples. The data sets can be downloaded from the collection of LSVM problems in http: //www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/. For both experiments we set c = 10, which resulted in more than 99% classification accuracy of the used data sets.
By [20,   methods were terminated after one hour of running time. Notice that RCD versions were more efficient than both UCDC versions, with RCD v.1 being the fastest among all. An interesting observation in Figures 6.2a and 6.2c is that RCD versions had similar per iteration computational complexity since they performed similar number of iterations within the maximum allowed running time. However, for RCD v.1, it seems that diagonal information from the second order derivatives of f was enough to decrease faster the objective function for all iterations compared to RCD v.2. Finally, in this experiment we observed that both RCD versions accepted unit step sizes for a major part of the process.

Conclusion.
We presented a robust randomized block coordinate descent method for composite function problems (1.1), which we name Robust Coordinate Descent (RCD), that can properly handle second-order (curvature) information. The proposed method can vary from first-to second-order; depending on how large the block updates are set, how accurate second-order information are used and how inexactly the arising subproblems are solved. Although the per iteration computational complexity might be higher for RCD, we present synthetic and real world large scale examples where the number of iterations substantially decreases, as well as the overall time.
From the theoretical point of view, we prove global convergence of RCD and under standard assumptions we show that RCD exhibits on expectation block quadratic or superlinear rate of convergence.