Performance of first- and second-order methods for 1-regularized least squares problems

We study the performance of ﬁrst- and second-order optimization methods for (cid:96) 1 -regularized sparse least-squares problems as the conditioning of the problem changes and the dimensions of the problem increase up to one trillion. A rigorously deﬁned generator is presented which allows control of the dimensions, the conditioning and the sparsity of the problem. The generator has very low memory requirements and scales well with the dimensions of the problem.


Introduction
We consider the problem minimize f τ (x) := τ x 1 + 1 2 where x ∈ R n , · 1 denotes the 1 -norm, · 2 denotes the Euclidean norm, τ > 0, A ∈ R m×n and b ∈ R m . An application that is formulated as in (1) is sparse data fitting, where the aim is to approximate n-dimensional sampled points (rows of matrix A) using a linear function, which depends on less than n variables, i.e., its slope is a sparse vector. Let us assume that we sample m data points (a i , b i ), where a i ∈ R n and b i ∈ R ∀i = 1, 2, · · · , m. We assume linear dependence of b i on a i : where e i is an error term due to the sampling process being innacurate. Depending on the application some statistical information is assumed about vector e. In matrix form the previous relationship is: where A ∈ R m×n is a matrix with a i 's as its rows and b ∈ R m is a vector with b i 's as its components. The goal is to find a sparse vector x (with many zero components) such that the error Ax − b 2 is minimized. To find x one can solve problem (1). The purpose of the 1 norm in (1) is to promote sparsity in the optimal solution [1]. An example that demonstrates the purpose of the 1 norm is presented in Figure 1. Figure 1 shows a two dimensional instance where n = 2, m = 1000 and matrix A is full-rank. Notice that the data points a i ∀i = 1, 2, · · · , m have large variations with respect to feature [a i ] 1 ∀i, where [·] j is the jth component of the input vector, while there is only a small variation with respect to feature [a i ] 2 ∀i. This property is captured when problem (1) is solved with τ = 30. The fitted plane in Figure 1a depends only on the first feature [a] 1 , while the second feature [a] 2 is ignored because [x * ] 2 = 0, where x * is the optimal solution of (1). This can be observed through the level sets of the plane shown with the colored map; for each value of [a] 1 the level sets remain constant for all values of [a] 2 . On the contrary, this is not the case when one solves a simple least squares problem (τ = 0 in (1)). Observe in Figure 1a that the fitted plane depends on both features [a] 1 and [a] 2 . A variety of sparse data fitting applications originate from the fields of signal processing and statistics. Five representative examples are briefly described below.
-Magnetic Resonance Imaging (MRI): A medical imaging tool used to scan the anatomy and the physiology of a body [27]. -Image inpainting: A technique for reconstructing degraded parts of an image [7]. -Image deblurring: Image processing tool for removing the blurriness of a photo caused by natural phenomena, such as motion [21]. -Genome-Wide Association study (GWA): DNA comparison between two groups of people (with/without a disease) in order to investigate factors that a disease depends on [42]. -Estimation of global temperature based on historic data [22].
Data fitting problems frequently require the analysis of large scale data sets, i.e., gigabytes or terabytes of data. In order to address large scale problems  there has been a resurgence in methods with computationally inexpensive iterations. For example many first-order methods were recovered and refined, such as coordinate descent [17, 24, 34, 39-41, 44, 45], alternating direction method of multipliers [9,15,18,23,43], proximal first-order methods [2,12,33] and firstorder smoothing methods [5,6,30]. The previous are just few representative examples, the list is too long for a complete demonstration, many other examples can be found in [1,11]. Often the goal of modern first-order methods is to reduce the computational complexity per iteration, while preserving the theoretical worst case iteration complexity of classic first-order methods [29]. Many modern first order methods meet the previous goal. For instance, coordinate descent methods can have up to n times less computational complexity per iteration [35,34].
First-order methods have been very successful in various scientific fields, such as support vector machine [46], compressed sensing [14], image processing [12] and data fitting [22]. Several new first-order type approaches have recently been proposed for various imaging problems in the special issue edited by M. Bertero, V. Ruggiero and L. Zanni [8]. However, even for the simple unconstrained problems that arise in the previous fields there exist more challenging instances. Since first-order methods do not capture sufficient second-order information, their performance might degrade unless the problems are well conditioned [16]. On the other hand, the second-order methods capture the curvature of the objective function sufficiently well, but by consensus they are usually applied only on medium scale problems or when high precision accuracy is required. In particular, it is frequently claimed [2,5,6,20,38] that the second-order methods do not scale favourably as the dimensions of the problem increase because of their high computational complexity per iteration. Such claims are based on an assumption that a full second-order information has to be used. However, there is evidence [16,19] that for non-trivial problems, inexact second-order methods can be very efficient.
In this paper we will exhaustively study the performance of first-and second-order methods. We will perform numerical experiments for large-scale problems with sizes up to one trillion of variables. We will examine conditions under which certain methods are favoured or not. We hope that by the end of this paper the reader will have a clear view about the performance of firstand second-order methods.
Another contribution of the paper is the development of a rigorously defined instance generator for problems of the form of (1). The most important feature of the generator is that it scales well with the size of the problem and can inexpensively create instances where the user controls the sparsity and the conditioning of the problem. For example see Subsection 8.9, where an instance of one trillion variables is created using the proposed generator. We believe that the flexibility of the proposed generator will cover the need for generation of various good test problems.
This paper is organised as follows. In Section 2 we briefly discuss the structure of first-and second-order methods. In Section 3 we give the details of the instance generator. In Section 4 we provide examples for constructing matrix A. In Section 5, we present some measures of the conditioning of problem (1). These measures will be used to examine the performance of the methods in the numerical experiments. In Section 6 we discuss how the optimal solution of the problem is selected. In Section 7 we briefly describe known problem generators and explain how our propositions add value to the existing approaches. In Section 8 we present the practical performance of first-and second-order methods as the conditioning and the size of the problems vary. Finally, in Section 9 we give our conclusions.
2 Brief discussion on first-and second-order methods We are concerned with the performance of unconstrained optimization methods which have the following intuitive setting. At every iteration a convex function Q τ (y; x) : R n → R is created that locally approximates f τ at a given point x. Then, function Q τ is minimized to obtain the next point. An example that covers the previous setting is the Generic Algorithmic Framework (GFrame) which is given below. Details of GFrame for each method used in this paper are presented in Section 8.
Loosely speaking, close to the optimal solution of problem (1), the better the approximation Q τ of f τ at any point x the fewer iterations are required to solve (1). On the other hand, the practical performance of such methods is a trade-off between careful incorporation of the curvature of f τ , i.e. secondorder derivative information in Q τ and the cost of solving subproblem (3) in GFrame.
Discussion on two examples of Q τ which consider different trade-off follows. First, let us fix the structure of Q τ for problem (1) to be 1: Initialize x 0 ∈ R n and y 0 ∈ R n For k = 0, 1, 2, . . . until some termination criteria are satisfied 2: Create a convex function Qτ (y; y k ) that approximates fτ in a neighbourhood of y k 3: Approximately (or exactly) solve the subproblem 4: Find a step-size α > 0 based on some criteria and set end-for 5: Return approximate solution x k+1 Algorithm : Generic Framework (GFrame) where H ∈ R n×n is a positive definite matrix. Notice that the decision of creating Q τ has been reduced to a decision of selecting H. Ideally, matrix H should be chosen such that it represents curvature information of 1/2 Ax−b 2 at point x, i.e. matrix H should have similar spectral decomposition to A A.
2 ≤ 1} be a unit ball centered at x. Then, H should be selected in an optimal way: The previous problem simply states that H should minimize the sum of the absolute values of the residual f τ − Q τ over B. Using twice the fundamental theorem of calculus on f τ from x to y we have that (5) is equivalent to It is trivial to see that the best possible H is simply H = A A. However, this makes every subproblem (3) as difficult to be minimized as the original problem (1). One has to reevaluate the trade-off between a matrix H that sufficiently well represents curvature information of 1/2 Ax − b 2 at a point x compared to a simple matrix H that is not as good approximation but offers an inexpensive solution of subproblem (3). An example can be obtained by setting H to be a positively scaled identity, which gives a solution to problem (5) H = λ max (A A)I n , where λ max (·) denotes the largest eigenvalue of the input matrix and I n is the identity matrix of size n × n. The contours of such a function Q τ compared to those of function f τ are presented in Subfigure 2a.
Notice that the curvature information of function f τ is lost, this is because nearly all spectral properties of A A are lost. However, for such a function Q τ the subproblem (3) has an inexpensive closed form solution known as iterative shrinkage-thresholding [12,27]. The computational complexity per iteration is so low that one hopes that this will compensate for the losses of curvature information. Such methods, are called first-order methods and have been shown to be efficient for some large scale problems of the form of (1) [2].   (4), that is frequently used in first-order methods. In the right subfigure, function Q τ is a non separable quadratic (9) which is used in some of the second-order methods.
Another approach of constructing Q τ involves the approximation of the 1 -norm with the pseudo-Huber function where µ > 0 is an approximation parameter. This approach is frequently used by methods that aim in using at every iteration full information from the Hessian matrix A A, see for example [16,19]. Using (7), problem (1) is replaced with The smaller µ is the better the approximation of problem (8) to (1). The advantage is that f µ τ in (8) is a smooth function which has derivatives of all degrees. Hence, smoothing will allow access to second-order information and essential curvature information will be exploited. However, for very small µ certain problems arise for optimization methods of the form of GFrame, see [16]. For example, the optimal solution of (1) is expected to have many zero components, on the other hand, the optimal solution of (8) is expected to have many nearly zero components. However, for small µ one can expect to obtain a good approximation of the optimal solution of (1) by solving (8). For the smooth problem (8), the convex approximation Q τ at x is: The contours of such a function Q τ compared to function f τ are presented in Subfigure 2b. Notice that Q τ captures the curvature information of function f τ . However, minimizing the subproblem (3) might be a more expensive operation. Therefore, we rely on an approximate solution of (3) using some iterative method which requires only simple matrix-vector product operations with matrices A and A . In other words we use only an approximate second-order information. It is frequently claimed [2,5,6,20,38] that second-order methods do not scale favourably with the dimensions of the problem because of the more costly task of solving approximately the subproblems in (3), instead of having an inexpensive closed form solution. Such claims are based on an assumption that full second-order information has to be used when solving subproblem (3). Clearly, this is not necessary: an approximate second-order information suffices. Studies in [16,19] provided theoretical background as well as the preliminary computational results to illustrate the issue. In this paper, we provide rich computational evidence which demonstrates that second-order methods can be very efficient.

Instance Generator
In this section we discuss an instance generator for (1) for the cases m ≥ n and m < n. The generator is inspired by the one presented in Section 6 of [32]. The advantage of our modified version is that it allows to control the properties of matrix A and the optimal solution x * of (1). For example, the sparsity of matrix A, its spectral decomposition, the sparsity and the norm of x * , since A and x * are defined by the user. Throughout the paper we will denote the i th component of a vector, by the name of the vector with subscript i. Whilst, the i th column of a matrix is denoted by the name of the matrix with subscript i.

Instance Generator for m ≥ n
Given τ > 0, A ∈ R m×n and x * ∈ R n the generator returns a vector b ∈ R m such that x * := arg min x f τ (x). For simplicity we assume that the given matrix A has rank n. The generator is described in Procedure IGen below.
In procedure IGen, given τ , A and x * we are aiming in finding a vector b such that x * satisfies the optimality conditions of problem (1) where ∂ x 1 = [−1, 1] n is the subdifferential of the 1 -norm at point x. By fixing a subradient g ∈ ∂ x * 1 as defined in (10) and setting e = b − Ax * , the previous optimality conditions can be written as 1: Initialize τ > 0, A ∈ R m×n with m ≥ n and rank n, x * ∈ R n 2: Construct g ∈ R n such that g ∈ ∂ x * 1 : Procedure : Instance Generator (IGen) The solution to the underdetermined system (11) is set to e = τ A(A A) −1 g and then we simply obtain b = Ax * + e; Steps 3 and 4 in IGen, respectively. Notice that for a general matrix A, Step 3 of IGen can be very expensive. Fortunately, using elementary linear transformations, such as Givens rotations, we can iteratively construct a sparse matrix A with a known singular value decomposition and guarantee that the inversion of matrix A A in Step 3 of IGen is trivial. We provide a more detailed argument in Section 4.

Instance Generator for m < n
In this subsection we extend the instance generator that was proposed in Subsection 3.1 to the case of matrix A ∈ R m×n with more columns than rows, i.e. m < n. Given τ > 0, B ∈ R m×m , N ∈ R m×n−m and x * ∈ R n the generator returns a vector b ∈ R m and a matrix A ∈ R m×n such that For this generator we need to discuss first some restrictions on matrix A and the optimal solution x * . Let with |S| = s and A S ∈ R m×s be a collection of columns from matrix A which correspond to indices in S. Matrix A S must have rank s otherwise problem (1) is not well-defined. To see this, let sign(x * S ) ∈ R s be the sign function applied component-wise to x * S , where x * S is a vector with components of x * that correspond to indices in S. Then problem (1) reduces to the following: where x s ∈ R s . The first-order stationary points of problem (13) satisfy If rank(A S ) < s, the previous linear system does not have a unique solution and problem (1) does not have a unique minimizer. Having this restriction in mind, let us now present the instance generator for m < n in Procedure IGen2 below.
1: Initialize τ > 0, B ∈ R m×m with rank m, N ∈ R m×n−m , x * ∈ R n with S := {1, 2, · · · , s} and s ≤ m 2: Construct g ∈ R m such that g ∈ ∂ x * (1, 2, · · · , m) 1 : 3: Set e = τ B − g 4: Construct matrixÑ ∈ R m×n−m with the following loop: Without loss of generality it is assumed that all nonzero components of x * correspond to indices in S = {1, 2, · · · , s}. By fixing a partial subradient g ∈ ∂ x * (1, 2, · · · , m) 1 as in (14), where x * (1, 2, · · · , m) ∈ R m is a vector which consists of the first m components of x * , and defining a vector e = b − Ax * , the previous optimality conditions can be written as: It is easy to check that by definingÑ as in Step 4 of IGen2 conditions (15) are satisfied. Finally, we obtain b = Ax * + e. Similarly to IGen in Subsection (3.1), for Step 3 in IGen2 we have to perform a matrix inversion, which generally can be an expensive operation. However, in the next section we discuss techniques how this matrix inversion can be executed using a sequence of elementary orthogonal transformations.

Construction of matrix A
In this subsection we provide a paradigm on how matrix A can be inexpensively constructed such that its singular value decomposition is known and its sparsity is controlled. We examine the case of instance generator IGen where m ≥ n. The paradigm can be easily extended to the case of IGen2, where m < n.
Let Σ ∈ R m×n be a rectangular matrix with the singular values σ 1 , σ 2 , · · · , σ n on its diagonal and zeros elsewhere: where O m−n×n ∈ R m−n×n is a matrix of zeros, and let G(i, j, θ) ∈ R n×n be a Givens rotation matrix, which rotates plane i-j by an angle θ: where i, j ∈ {1, 2, · · · , n}, c = cos θ and s = sin θ. Given a sequence of Givens rotations {G(i k , j k , θ k )} K k=1 we define the following composition of them: Similarly, letG(l, p, ϑ) ∈ R m×m be a Givens rotation matrix where l, p ∈ {1, 2, · · · , m} and be a composition ofK Givens rotations. Using G andG we define matrix A as where P 1 , P 2 ∈ R m×m are permutation matrices. Since the matrices P 1G P 2 and G are orthonormal it is clear that the left singular vectors of matrix A are the columns of P 1G P 2 , Σ is the matrix of singular values and the right singular vectors are the columns of G. Hence, in Step 3 of IGen we simply set (A A) −1 = G(Σ Σ) −1 G , which means that Step 3 in IGen costs two matrixvector products with G and a diagonal scaling with (Σ Σ) −1 . Moreover, the sparsity of matrix A is controlled by the numbers K andK of Givens rotations, the type, i.e. (i, j, θ) and (l, p, ϑ), and the order of Givens rotations. Also, notice that the sparsity of matrix A A is controlled only by matrix G. Examples are given in Subsection 4.1.
It is important to mention that other settings of matrix A in (16) could be used, for example different combinations of permutation matrices and Givens rotations. The setting chosen in (16) is flexible, it allows for an inexpensive construction of matrix A and makes the control of the singular value decomposition and the sparsity of matrices A and A A easy.
Notice that matrix A does not have to be calculated and stored. In particular, in case that the method which is applied to solve problem (1) requires only matrix-vector product operations using matrices A and A , one can simply consider matrix A as an operator. It is only required to predefine the triplets (i k , j k , θ k ) ∀k = 1, 2, · · · , K for matrix G, the triplets (l k , p k , θ k ) ∀k = 1, 2, · · · ,K for matrixG and the permutation matrices P 1 and P 2 . The previous implies that the generator is inexpensive in terms of memory requirements. Examples of matrix-vector product operations with matrices A and A in case of (16) are given below in Algorithms MvPA and MvPAt, respectively. 1: Given a matrix A defined as in (16) and an input vector x ∈ R n , do 2: Set y 0 = x For k = 1, 2, . . . , K 3: Algorithm : Matrix-vector product with A (MvPA) 1: Given a matrix A defined as in (16) and input vector y ∈ R m , do 2: For k = 1, 2, . . . , K 5: x k = G K−k+1 x k−1 end-for 6: Return x K Algorithm : Matrix-vector product with A (MvPAt)

An example using Givens rotation
Let us assume that m, n are divisible by two and m ≥ n. Given the singular values matrix Σ and rotation angles θ and ϑ, we construct matrix A as where P is a random permutation of the identity matrix, G is a composition of n/2 Givens rotations: andG is a composition of m/2 Givens rotations: Notice that the angle θ is the same for all Givens rotations G k , this means that the total memory requirement for matrix G is low. In particular, it consists only of the storage of a 2 × 2 rotation matrix. Similarly, the memory requirement for matrixG is also low.

Control of sparsity of matrix A and A A
We now present examples in which we demonstrate how sparsity of matrix A can be controlled through Givens rotations.
In the example of Subsection 4.1, two compositions of n/2 and m/2 Givens rotations, denoted by G andG, are applied on an initial diagonal rectangular matrix Σ. If n = 2 3 and m = 2n the sparsity pattern of the resulting matrix A = (PGP )ΣG is given in Subfigure 3a and has 28 nonzero elements, while the sparsity pattern of matrix A A is given in Subfigure 4a and has 16 nonzero elements. Notice in this subfigure that the coordinates can be clustered in pairs of coordinates (1, 2), (3,4), (5, 6) and (7,8). One could apply another stage of Givens rotations. For example, one could construct matrix

Conditioning of the problem
Let us now precisely define how we measure the conditioning of problem (1). For simplicity, throughout this section we assume that matrix A has more rows than columns, m ≥ n, and it is full-rank. Extension to the case of matrix A with more columns than rows is easy and we briefly discuss this at the end of this section. We denote with span(·) the span of the columns of the input matrix. Moreover, S is defined in (12), S c is its complement.
Two factors are considered that affect the conditioning of the problem. First, the usual condition number of the second-order derivative of 1/2 Ax − b 2 2 in (1), which is simply κ(A A) = λ 1 (A A)/λ n (A A), where 0 < λ n ≤ λ n−1 ≤ · · · ≤ λ 1 are the eigenvalues of matrix A A. It is well-known that the larger κ(A A) is, the more difficult problem (1) becomes.
Second, the conditioning of the optimal solution x * of problem (1). Let us explain what we mean by the conditioning of x * . We define a constant ρ > 0 and the index set I ρ := {i ∈ {1, 2, · · · , n} | λ i (A A) ≥ ρ}. Furthermore, we define the projection P ρ = G ρ G ρ , where G ρ ∈ R n×r , r = |I ρ | and matrix G ρ has as columns the eigenvectors of matrix A A which correspond to eigenvalues with indices in I ρ . Then, the conditioning of x * is defined as For the case P ρ x * = 0, the denominator of (17) is the mass of x * which exists in the space spanned by eigenvectors of A A which correspond to eigenvalues that are larger than or equal to ρ. Let us assume that there exists some ρ which satisfies λ n (A A) ≤ ρ λ 1 (A A). If κ ρ (x * ) is large, i.e., P ρ x * 2 is close to zero, then the majority of the mass of x * is "hidden" in the space spanned by eigenvectors which correspond to eigenvalues that are smaller than ρ, i.e., the orthogonal space of span(G ρ ). In Section 1 we referred to methods that do not incorporate information which correspond to small eigenvalues of A A. Therefore, if the previous scenario holds, then we expect the performance of such methods to degrade. In Section 8 we empirically verify the previous arguments.
If matrix A has more columns than rows then the previous definitions of conditioning of problem (1) are incorrect and need to be adjusted. Indeed, if m < n and rank(A) = min(m, n) = m, then A A is a rank deficient matrix which has m nonzero eigenvalues and n − m zero eigenvalues. However, we can restrict the conditioning of the problem to a neighbourhood of the optimal solution of x * . In particular, let us define a neighbourhood of x * so that all points in this neighbourhood have nonzeros at the same indices as x * and zeros elsewhere, i.e. N := {x ∈ R n | x i = 0 ∀i ∈ S, x i = 0 ∀i ∈ S c }. In this case an important feature to determine the conditioning of the problem is the ratio of the largest and the smallest nonzero eigenvalues of A S A S , where A S is a submatrix of A built of columns of A which belong to set S.

Construction of the optimal solution
Two different techniques are employed to generate the optimal solution x * for the experiments presented in Section 8. The first procedure suggests a simple random generation of x * , see Procedure OsGen below.
Procedure : Optimal solution Generator (OsGen) The second and more complicated approach is given in Procedure OsGen2. This procedure is applicable only in the case that m ≥ n, however, it can be easily extended to the case of m < n. We focus in the former scenario since all experiments in Section 8 are generated by setting m ≥ n.
where 1n ∈ R n is a vector of ones and · 0 is the zero norm which returns the number of nonzero components of the input vector. Problem (18) can be solved approximately using an Orthogonal Matching Pursuit (OMP) [28] solver implemented in [3].

Procedure : Optimal solution Generator 2 (OsGen2)
The aim of Procedure OsGen2 is to find a sparse x * with κ ρ (x * ) arbitrarily large for some ρ in the interval λ n (A A) ≤ ρ λ 1 (A A). In particular, OsGen2 will return a sparse x * which can be expressed as x * = Gv. The coefficients v are close to the inverse of the eigenvalues of matrix A A. Intuitively, this technique will create an x * which has strong dependence on subspaces which correspond to small eigenvalues of A A. The constant γ is used in order to control the norm of x * .
The sparsity constraint in problem (18), i.e., x 0 ≤ s, makes the approximate solution of this problem difficult when we use OMP, especially in the case that s and n are large. To avoid this expensive task we can ignore the sparsity constraint in (18). Then we can solve exactly and inexpensively the unconstrained problem and finally we can project the obtained solution in the feasible set defined by the sparsity constraint. Obviously, there is no guarantee that the projected solution is a good approximation to the one obtained in Step 2 of Procedure OsGen2. However, for all experiments in Section 8 that we applied this modification we obtained sufficiently large κ ρ (x * ). This means that our objective to produce ill-conditioned optimal solutions was met, while we kept the computational costs low. The modified version of Procedure OsGen2 is given in Procedure OsGen3.
1: Given the required number s ≤ min(m, n) of nonzeros in x * , two non negative integers s 1 and s 2 such that s 1 + s 2 = s, a positive constant γ > 0, the right singular vectors G and singular values Σ of matrix A do: 2: Solve exactly x * := arg min where 1n ∈ R n is a vector of ones. Problem (19) can be solved exactly and inexpensively because G is an orthonormal matrix.

Existing Problem Generators
So far in Section 3.1 we have described in details our proposed problem generator. Moreover, in Section 4 we have described how to construct matrices A such that the proposed generator is scalable with respect to the number of unknown variables. We now briefly describe existing problem generators and explain how our propositions add value to the existing approaches.
Given a regularization parameter τ existing problem generators are looking for A, b and x * such that the optimality conditions of problem (1): are satisfied. For example, in [32] the author fixes a vector of noise e and an optimal solution x * and then finds A and b such that (20) is satisfied. In particular, in [32] matrix A = BD is used, where B is a fixed matrix and D is a scaling matrix such that the following holds.
(BD) e ∈ τ ∂ x * 1 , Matrix D is trivial to calculate, see Section 6 in [32] for details. Then by setting b = Ax * + e (20) is satisfied. The advantage of this generator is that it allows control of the noise vector e, in comparison to our approach where the vector noise e has to be determined by solving a linear system. On the other hand, one does not have direct control over the singular value decomposition of matrix A, since this depends on matrix D, which is determined based on the fixed vectors e and x * . Another representative example is proposed in [26]. This generator, which we discovered during the revision of our paper, proposes the same setting as in our paper. In particular, given A, x * and τ one can construct a vector b (or a noise vector e) such that (20) is satisfied. However, in [26] the author suggests that b can be found using a simple iterative procedure. Depending on matrix A and how ill-conditioned it is, this procedure might be slow. In this paper, we suggest that one can rely on numerical linear algebra tools, such as Givens rotation, in order to inexpensively construct b (or a noise vector e) using straightforwardly scalable operations. Additionally, we show in Section (8) that a simple construction of matrix A is sufficient to extensively test the performance of methods.

Numerical Experiments
In this section we study the performance of state-of-the-art first-and secondorder methods as the conditioning and the dimensions of the problem increase. The scripts that reproduce the experiments in this section as well as the problem generators that are described in Section 3 can be downloaded from: http://www.maths.ed.ac.uk/ERGO/trillion/.
-FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) [2] is an optimal first-order method for problem (1), which adheres to the structure of GFrame. At a point x, FISTA builds a convex function: where L is an upper bound of λ max (A A), and solves subproblem (3) exactly using shringkage-thresholding [12,27]. An efficient implementation of this algorithm can be found as part of TFOCS (Templates for First-Order Conic Solvers) package [6] under the name N83. In this implementation the parameter L is calculated dynamically. -PCDM (Parallel Coordinate Descent Method) [35] is a randomized parallel coordinate descent method. The parallel updates are performed asynchronously and the coordinates to be updated are chosen uniformly at random. Let be the number of processors that are employed by PCDM. Then, at a point x, PCDM builds convex approximations: is the ith diagonal element of matrix A A and β is a positive constant which is defined in Subsection 8.3. The Q i τ functions are minimized exactly using shrinkage-thresholding.
-PSSgb (Projected Scaled Subgradient, Gafni-Bertsekas variant) [37] is a second-order method. At each iteration of PSSgb the coordinates are separated into two sets, the working set W and the active set A. The working set consists of all coordinates for which, the current point x is nonzero. The active set is the complement of the working set W. The following local quadratic model is build at each iteration  [16] is also a secondorder method. At every point x pdNCG constructs a convex function Q τ exactly as described for (9). The subproblem (3) is solved inexactly by reducing it to the linear system: , which is solved approximately using preconditioned Conjugate Gradients (PCG). A simple diagonal preconditioner is used for all experiments. The preconditioner is the inverse of the diagonal of matrix ∇ 2 f µ τ (x).

Implementation details
Solvers pdNCG, FISTA and PSSgb are implemented in MATLAB, while solver PCDM is a C++ implementation. We expect that the programming language will not be an obstacle for pdNCG, FISTA and PSSgb. This is because these methods rely only on basic linear algebra operations, such as the dot product, which are implemented in C++ in MATLAB by default. According to the most recent list of commercial supercomputers, which is published in TOP500 list (http://www.top500.org), ARCHER is currently the 25th fastest supercomputer worldwide out of 500 supercomputers. ARCHER has a total of 118, 080 cores with performance 1, 642.54 TFlops/s on LIN-PACK benchmark and 2, 550.53 TFlops/s theoretical peak perfomance. The most computationally demanding experiments which are presented in Subsection 8.9 required more than half of the cores of ARCHER, i.e., 65, 536 cores out of 118, 080.

Parameter tuning
We describe the most important parameters for each solver, any other parameters are set to their default values. For pdNCG we set the smoothing parameter to µ = 10 −5 , this setting allows accurate solution of the original problem with an error of order O(µ) [31]. For pdNCG, PCG is terminated when the relative residual is less that 10 −1 and the backtracking line-search is terminated if it exceeds 50 iterations. Regarding FISTA the most important parameter is the calculation of the Lipschitz constant L, which is handled dynamically by TFOCS. For PCDM the coordinate Lipschitz constants L i ∀i = 1, 2, · · · , n are calculated exactly and parameter β = 1 + (ω − 1)( − 1)/(n − 1), where ω changes for every problem since it is the degree of partial separability of the fidelity function in (1), which is easily calculated (see [35]), and = 40 is the number of cores that are used. For PSSgb we set the number of L-BFGS corrections to 10.
We set the regularization parameter τ = 1, unless stated otherwise. We run pdNCG for sufficient time such that the problems are adequately solved. Then, the rest of the methods are terminated when the objective function f τ in (1) is below the one obtained by pdNCG or when a predefined maximum number of iterations limit is reached. All comparisons are presented in figures which show the progress of the objective function against the wall clock time. This way the reader can compare the performance of the solvers for various levels of accuracy. We use logarithmic scales for the wall clock time and terminate runs which do not converge in about 10 5 seconds, i.e., approximately 27 hours.

Increasing condition number of A A
In this experiment we present the performance of FISTA, PCDM, PSSgb and pdNCG for increasing condition number of matrix A T A when Procedure Os-Gen is used to construct the optimal solution x * . We generate six matrices A and two instances of x * for every matrix A; twelve instances in total.
The singular value decomposition of matrix A is A = ΣG , where Σ is the matrix of singular values, the columns of matrices I m and G are the left and right singular vectors, respectively, see Subsection 4.1 for details about the construction of matrix G. The singular values of matrices A are chosen uniformly at random in the intervals [0, 10 q ], where q = 0, 1, · · · , 5, for each of the six matrices A. Then, all singular values are shifted by 10 −1 . The previous resulted in a condition number of matrix A A which varies from 10 2 to 10 12 with a step of times 10 2 . The rotation angle θ of matrix G is set to 2π/3 radians. Matrices A have n = 2 22 columns, m = 2n rows and rank n. The optimal solutions x * have s = n/2 7 nonzero components for all twelve instances.
For the first set of six instances we set γ = 10 in OsGen, which resulted in κ 0.1 (x * ) ≈ 1 for all experiments. The results are presented in Figure 5. For these instances PCDM is clearly the fastest for κ(A A) ≤ 10 4 , while for κ(A A) ≥ 10 6 pdNCG is the most efficient.
For the second set of six instances we set γ = 10 3 in Procedure OsGen, which resulted in the same κ 0.1 (x * ) as before for every matrix A. The results are presented in Figure 6. For these instances PCDM is the fastest for very well conditioned problems with κ(A A) ≤ 10 2 , while pdNCG is the fastest for κ(A A) ≥ 10 4 .
We observed that pdNCG required at most 30 iterations to converge for all experiments. For FISTA, PCDM and PSSgb the number of iterations was varying between thousands and tens of thousands iterations depending on the condition number of matrix A A; the larger the condition number the more the iterations. However, the number of iterations is not a fair metric to compare solvers because every solver has different computational cost per iteration. In particular, FISTA, PCDM and PSSgb perform few inner products per iteration, which makes every iteration inexpensive, but the number of iterations is sensitive to the condition number of matrix A A. On the other hand, for pdNCG the empirical iteration complexity is fairly stable, however, the number of inner products per iteration (mainly matrix-vector products with matrix A) may increase as the condition number of matrix A A increases. Inner products are the major computational burden at every iteration for all solvers, therefore, the faster an algorithm converged in terms of wall-clock time the less inner products that are calculated. In Figures 5 and 6 we display the objective evaluation against wall-clock time (log-scale) to facilitate the comparison of different algorithms.

Increasing condition number of A A: non-trivial construction of x *
In this experiment we examine the performance of the methods as the condition number of matrix A A increases, while the optimal solution x * is generated using Procedure OsGen3 (instead of OsGen) with γ = 100 and   The two classes of experiments are distinguished based on the rotation angle θ that is used for the composition of Givens rotations G. In particular, for the first class of experiments the angle is θ = 2π/10 radians, while for the second class of experiments the rotation angle is θ = 2π/10 3 radians. The difference between the two classes is that the second class consists of matrices A A for which, a major part of their mass is concentrated in the diagonal. This setting is beneficial for PCDM since it uses information only from the diagonal of matrices A A. This setting is also beneficial for pdNCG since it uses a diagonal preconditioner for the inexact solution of linear systems at every iteration.
The results for the first class of experiments are presented in Figure 7. For instances with κ(A A) ≥ 10 6 PCDM was terminated after 1, 000, 000 iterations, which corresponded to more than 27 hours of wall-clock time.
The results for the second class of experiments are presented in Figure 8. Notice in this figure that the objective function is only slightly reduced. This does not mean that the initial solution, which was the zero vector, was nearly optimal. This is because noise with large norm, i.e., Ax * − b is large, was used in these experiments, therefore, changes in the optimal solution did not have large affect on the objective function.

Increasing dimensions
In this experiment we present the performance of pdNCG, FISTA, PCDM and PSSgb as the number of variables n increases. We generate four instances where the number of variables n takes values 2 20 , 2 22 , 2 24 and 2 26 , respectively. The singular value decomposition of matrix A is A = ΣG . The singular values in matrix Σ are chosen uniformly at random in the interval [0, 10] and then are shifted by 10 −1 , which resulted in κ(A A) ≈ 10 4 . The rotation angle θ of matrix G is set to 2π/10 radians. Moreover, matrices A have m = 2n rows and rank n. The optimal solutions x * have s = n/2 7 nonzero components for each generated instance. For the construction of the optimal solutions x * we use Procedure OsGen3 with γ = 100 and s 1 = s 2 = s/2, which resulted in κ 0.1 (x * ) ≈ 3 on average.
The results of this experiment are presented in Figure 9. Notice that all methods have a linear-like scaling with respect to the size of the problem.

Increasing density of matrix A A
In this experiment we demonstrate the performance of pdNCG, FISTA, PCDM and PSSgb as the density of matrix A A increases. We generate four instances (A, x * ). For the first experiment we generate matrix A = ΣG , where Σ is the matrix of singular values, the columns of matrices I m and G are the   left and right singular vectors, respectively. For the second experiment we generate matrix A = Σ(G 2 G) , where the columns of matrices I m and G 2 G are the left and right singular vectors of matrix A, respectively; G 2 has been defined in Subsection 4.2. Finally, for the third and fourth experiments we have A = Σ(GG 2 G) and A = Σ(G 2 GG 2 G) , respectively. For each experiment the singular values of matrix A are chosen uniformly at random in the interval [0, 10] and then are shifted by 10 −1 , which resulted in κ(A A) ≈ 10 4 . The rotation angle θ of matrices G and G 2 is set to 2π/10 radians. Matrices A have m = 2n rows, rank n and n = 2 22 . The optimal solutions x * have s = n/2 7 nonzero components for each experiment. Moreover, Procedure OsGen3 is used with γ = 100 and s 1 = s 2 = s/2 for the construction of x * for each experiment, which resulted in κ 0.1 (x * ) ≈ 2 on average.   The results of this experiment are presented in Figure 10. Observe, that all methods had a robust performance with respect to the density of matrix A.

Varying parameter τ
In this experiment we present the performance of pdNCG, FISTA, PCDM and PSSgb as parameter τ varies from 10 −4 to 10 4 with a step of times 10 2 . We generate four instances (A, x * ), where matrix A = ΣG has m = 2n rows, rank n and n = 2 22 . The singular values of matrices A are chosen uniformly at random in the interval [0, 10] and then are shifted by 10 −1 , which resulted in κ(A A) ≈ 10 4 for each experiment. The rotation angles θ for matrix G in A is set to 2π/10 radians. The optimal solution x * has s = n/2 7 nonzero components for all instances. Moreover, the optimal solutions are generated  using Procedure OsGen3 with γ = 100, which resulted in κ 0.1 (x * ) ≈ 3 for all four instances.
The performance of the methods is presented in Figure 11. Notice in Subfigure 11d that for pdNCG the objective function f τ is not always decreasing monotonically. A possible explanation is that the backtracking line-search of pdNCG, which guarantees monotonic decrease of the objective function [16], terminates in case that 50 backtracking iterations are exceeded, regardless if certain termination criteria are satisfied.

Performance of a second-order method on huge scale problems
We now present the performance of pdNCG on synthetic huge scale (up to one trillion variables) S-LS problems as the number of variables and the number of processors increase.    We generate six instances (A, x * ), where the number of variables n takes values 2 30 , 2 32 , 2 34 , 2 36 , 2 38 and 2 40 . Matrices A = ΣG have m = 2n rows and rank n. The singular values σ i for i = 1, 2, · · · , n of matrices A are set to 10 −1 for odd i's and 10 2 for even i's. The rotation angle θ of matrix G is set to 2π/3 radians. The optimal solutions x * have s = n/2 10 nonzero components for each experiment. In order to simplify the practical generation of this problem the optimal solutions x * are set to have s/2 components equal to −10 4 and the rest of nonzero components are set equal to 10 −1 .
Details of the performance of pdNCG are given in Table 1. Observe the nearly linear scaling of pdNCG with respect to the number of variables n and the number of processors. For all experiments in Table 1 pdNCG required 8 Newton steps to converge, 100 PCG iterations per Newton step on average, where every PCG iteration requires two matrix-vector products with matrix A. Observe in Subfigure 11d that for pdNCG the objective function f τ is not always decreasing monotonically. This is due to the backtracking line-search of pdNCG, which terminates in case that the maximum number of backtracking iterations is exceeded regardless if certain termination criteria are satisfied.  Table 1: Performance of pdNCG for synthetic huge scale S-LS problems. All problems have been solved to a relative error of order 10 −4 of the obtained solution

Conclusion
In this paper we developed an instance generator for 1 -regularized sparse least-squares problems. The generator is aimed for the construction of very large-scale instances. Therefore it scales well as the number of variables increases, both in terms of memory requirements and time. Additionally, the generator allows control of the conditioning and the sparsity of the problem. Examples are provided on how to exploit the previous advantages of the proposed generator. We believe that the optimization community needs such a generator to be able to perform fair assessment of new algorithms.
Using the proposed generator we constructed very large-scale sparse instances (up to one trillion variables), which vary from very well-conditioned to moderately ill-conditioned. We examined the performance of several representative first-and second-order optimization methods. The experiments revealed that regardless of the size of the problem, the performance of the methods crucially depends on the conditioning of the problem. In particular, the first-order methods PCDM and FISTA are faster for problems with small or moderate condition number, whilst, the second-order method pdNCG is much more efficient for ill-conditioned problems.