An Investigation of NewtonSketch and Subsampled Newton Methods
Abstract
Sketching, a dimensionality reduction technique, has received much attention in the statistics community. In this paper, we study sketching in the context of Newton’s method for solving finitesum optimization problems in which the number of variables and data points are both large. We study two forms of sketching that perform dimensionality reduction in data space: Hessian subsampling and randomized Hadamard transformations. Each has its own advantages, and their relative tradeoffs have not been investigated in the optimization literature. Our study focuses on practical versions of the two methods in which the resulting linear systems of equations are solved approximately, at every iteration, using an iterative solver. The advantages of using the conjugate gradient method vs. a stochastic gradient iteration are revealed through a set of numerical experiments, and a complexity analysis of the Hessian subsampling method is presented.
Keywords: sketching, subsampling, Newton’s method, machine learning, stochastic optimization
1 Introduction
In this paper, we consider second order methods for solving the finite sum problem,
(1.1) 
in which each component function is smooth and convex. Second order information is incorporated into the algorithms using sketching techniques, which construct an approximate Hessian by performing dimensionality reduction in data space. Two forms of sketching are of interest in the context of problem (1.1). The simplest form, referred to as Hessian subsampling, selects individual Hessians at random and averages them. The second variant consists of sketching via randomized Hadamard transforms, which may yield a better approximation of the Hessian but can be computationally expensive. A study of the tradeoffs between these two approaches is the subject of this paper.
Our focus is on practical implementations designed for problems in which the number of variables and the number of data points are both large. A key ingredient is to avoid the inversion of large matrices by solving linear systems (inaccurately) at each iteration using an iterative linear solver. This leads us to the question of which is the most effective linear solver for the type of algorithms and problems studied in this paper. In the light of the recent popularity of the stochastic gradient method, it is natural to investigate whether it could be preferable (as an iterative linear solver) to the conjugate gradient method, which has been a workhorse in numerical linear algebra and optimization. On the one hand, a stochastic gradient iteration can exploit the stochastic nature of the Hessian approximations; on the other hand the conjugate gradient (CG) method enjoys attractive convergence properties. We present theoretical and computational results contrasting these two approaches and highlighting the superior properties of the conjugate gradient method. Given the appeal of the Hessian subsampling method that utilizes a CG iterative solver, we derive computational complexity bounds for it.
A Newton method employing sketching techniques was proposed by Pilanci and Wainwright [16] for problems in which the number of variables is not large but the number of data points can be extremely large. They assume that the objective function is such that a square root of the Hessian can be computed at reasonable cost (as is the case with generalized linear models). They analyze the relative advantages of various sketching techniques and provide computational complexity results under the assumption that linear systems are solved exactly. However, they do not present extensive numerical results, nor do they compare their approaches with Hessian subsampling.
Second order optimization methods that employ Hessian subsampling techniques have been studied in [2, 5, 6, 13, 16, 18, 20, 21, 22]. RoostaKhorasani and Mahoney [18] derived convergence rates of subsampled Newton methods in probability for different sampling rates. Bollapragada et al. [2] provided convergence and complexity results in expectation for strongly convex functions, and employed CG to solve the linear system inexactly. Recently, several works [25, 24, 27] have studied the theoretical and empirical performance of subsampled Newton methods for solving nonconvex problems. Agarwal et al. [1] proposed to use a stochastic gradientlike method for solving the linear systems arising in subsampled Newton methods.
Given the popularity of first order methods (see e.g., [19, 3, 9, 15, 23] and the references therein), one may question whether algorithms that incorporate stochastic second order information have much to offer in practice for solving large scale instances of problem (1.1). In this paper, we argue that this is indeed the case in many important settings. Second order information can provide stability to the algorithm by facilitating the selection of the steplength parameter , and by countering the effects of ill conditioning, without greatly increasing the cost of the iteration.
The paper is organized as follows. In Section 2, we describe sketching techniques for second order optimization methods. Practical implementations of the algorithms are presented in Section 3, and the results of our numerical investigation are given in Section 4. We present a convergence analysis for the Hessian subsampling method in Section 5, where we also compare its work complexity with that of Newton methods using other sketching techniques such as the Hadamard transforms. Section 6 makes some final remarks about the contributions of the paper.
2 Stochastic Second Order Methods
The methods studied in this paper are of the form,
(2.1) 
where is the steplength parameter and is a stochastic approximation of the Hessian obtained via the sketching techniques described below. We show that a useful estimate of the Hessian can be obtained by reducing the dimensionality in data space, making it possible to tackle problems with extremely large datasets. We assume throughout the paper that the gradient is exact, and therefore, the stochastic nature of the iteration is entirely due to the the choice of . We could have employed a stochastic approximation to the gradient but opted not to do so because the use of an exact gradient allows us to study Hessian approximations in isolation and endows the algorithms with a linear convergence rate, which makes them directly comparable with first order methods.
We now discuss how the Hessian approximations are constructed.
HessianSketch
As proposed by Pilanci and Wainwright [16], the matrix can be defined as a randomized approximation to the true Hessian that utilizes a decomposition of the matrix and approximates the Hessian via random projections in lower dimensions. Specifically, at every iteration, the algorithm defines a random sketch matrix with the property that , and . Next, the algorithm computes a squareroot decomposition of the Hessian , which we denote by , and defines the Hessian approximation as
(2.2) 
The search direction of the algorithm is obtained by solving the linear system in (2.1). When forming is much less expensive than forming the Hessian .
This approach is not always viable because a square root Hessian matrix may not be easily computable, but there are important classes of problems when it is. Consider for example, generalized linear models where . The squareroot matrix can be expressed as , where is the data matrix, denotes the th row of , and is a diagonal matrix given by . The multiplication of the sketch matrix with leads to the projection of the latter matrix on to a smaller dimensional subspace.
Randomized projections can be performed based on the fast JohnsonLindenstrauss (JL) transform, e.g., Hadamard matrices, which allow for a fast matrixmatrix multiplication and have been shown to be efficient in least squares settings. We follow [10, 16] and define the Hadamard matrices as having entries . The randomized Hadamard projections are performed using the sketch matrix whose rows are chosen by randomly sampling rows of the matrix , where is a diagonal matrix with random entries. The properties of this approach have been analyzed in [16], where it is shown that each iteration has complexity , which is lower than the complexity of Newton’s method (), when and .
Hessian Subsampling
Hessian subsampling is a special case of sketching where the Hessian approximation is constructed by randomly sampling component Hessians of (1.1), and averaging them. Specifically, at the th iteration, one defines
(2.3) 
where the index set is chosen either uniformly at random (with or without replacement) or in a nonuniform manner as described in [26]. The search direction is defined as the (inexact) solution of the system of equations in (2.1) with given by (2.3). By computing an approximate solution of the linear system that is only as accurate as needed to ensure progress toward the solution, one achieves computational savings, as discussed in the next section. Newton methods that use Hessian subsampling have recently received considerable attention, as mentioned above, and some of their theoretical properties are well understood [2, 6, 18, 26]. Although Hessian subsampling can be thought of as a special case of sketching, we view it here as a distinct method because it can be motivated directly from the finite sum structure (1.1) without referring to randomized projections.
Notation and Nomenclature
Discussion
We thus have two distinct classes of methods at our disposal. On the one hand, NewtonSketch enjoys superior statistical properties (as discussed below) but has a high computational cost. On the other hand, the SSN method is simple to implement and is inexpensive, but could have inferior statistical properties. For example, when individual component functions are highly dissimilar, taking a small sample of individual Hessians at each iteration of the SSN method may not yield a useful approximation to the true Hessian. For such problems, the NewtonSketch method that forms the matrix based on a combination of all component functions may be more effective.
In the rest of the paper, we compare the computational and theoretical properties of the two approaches when solving finitesum problems involving many variables and large data sets. In order to undertake this investigation we need to give careful consideration to the largescale implementation of these methods.
3 Practical Implementations of the Algorithms
We can design efficient implementations of the NewtonSketch and subsampled Hessian method for the case when the number of variables is very large by giving careful consideration to how the linear system in (2.1) is formed and solved. Pilanci and Wainwright [16] proposed to implement the NewtonSketch method using direct solvers, but this is only practical when is not large. Since we do not wish to impose such a restriction on , we utilize iterative linear solvers tailored to each specific method. We consider two iterative methods: the conjugate gradient method, and a stochastic gradient iteration [1] that is inspired by (but it not equivalent to) the stochastic gradient method of Robbins and Monro [17].
3.1 The Conjugate Gradient (CG) Method
For large scale problems, it can be very effective to employ a Hessianfree approach, and compute an inexact solution of linear systems with coefficient matrices (2.2) or (2.3) using an iterative method that only requires matrixvector products instead of the matrix itself. The conjugate gradient method is the best known method in this category [7]. One can compute a search direction by running iterations of CG on the system
(3.1) 
Under a favorable distribution of the eigenvalues of the matrix , the CG method can generate a useful approximate solution in much fewer than steps; we exploit this result in the analysis presented in Section 5. We now discuss how to apply the CG method to the two classes of methods studied in this paper.
NewtonSketch
The linear system (2.2) arising in the NewtonSketch method requires access to a squareroot Hessian matrix. This is a rectangular matrix, denoted by , such that . Note that this is not the kind of decomposition (such as LU or QR) that facilitates the formation and solution of of linear systems involving the matrix (2.2). Therefore, when is large it is imperative to use an iterative method such as CG. The entire matrix need not be formed, and we only need to perform two matrixvector products with respect to the sketched squareroot Hessian. More specifically, for any vector one can compute
in two steps as follows
(3.2) 
Our implementation of the NewtonSketch method includes two tuning parameters: the number of rows in the sketch matrix , and the maximum number of CG iterations, , used to compute the step. We define the sketch matrix through the Hadamard basis that allows for a fast matrixvector multiplication, to obtain . This does not require forming the sketch matrix explicitly and can be executed at a cost of flops. The maximum cost of every iteration is given as the evaluation of the gradient plus component matrixvector products, where the factor 2 arises due to (3.2). We ignore the cost of forming the squareroot matrix, which is not large in the case of generalized linear models; our test problems are of this form. We also ignore the cost of multiplying by the sketch matrix, which can be large. Therefore, the results we report in the next section are the most optimistic for NewtonSketch.
Subsampled Newton
The subsampled Hessian used in the SSN method is significantly simpler to define compared to sketched Hessian, and the CG solver is naturally suited to a Hessianfree approach. The product of times vectors can be computed by simply summing the individual Hessianvector products. That is, given a vector ,
(3.3) 
We refer to this method as SSNCG method. Unlike NewtonSketch that requires the projection of the squareroot matrix onto the lower dimensional manifold, the SSNCG method can be implemented without the need of projections, and thus does not require additional computation beyond the matrixvector products (3.3).
The tuning parameters involved in the implementation of this method are the maximum size of the subsample, , and the maximum number of CG iterations used to compute the step, . Thus, the maximum cost of a SSNCG iteration is given by the cost of evaluating the gradient plus matrixvector products (3.3).
3.2 A Special Stochastic Gradient Iteration
We can also compute an approximate solution of the linear system in (2.1) using a stochastic gradient iteration similar to that proposed in [1] and studied in [2]. This stochastic gradient iteration, that we refer to as SGI, is appropriate only in the context of the SSN method because for NewtonSketch there is no efficient way to compute the stochastic gradient as the sketched Hessian cannot be decomposed into individual components. The resulting SSNSGI method operates in cycles of length . At every iteration of the cycle, an index is chosen at random, and a gradient step is computed for the following quadratic model at :
(3.4) 
This yields the iteration
(3.5) 
A method similar to SSNSGI has been analyzed in [1], where it is referred to as LiSSA, and [2] compares the complexity of LiSSA and SSNCG. None of those papers report detailed numerical results, and the effectiveness of the SGI approach was not known to us at the outset of this investigation.
The total cost of an (outer) iteration of the SSNSGI method is given by the evaluation of the gradient plus Hessianvector products of the form . This method has two tuning parameters: the inner iteration step length and the length of the inner SGI cycle.
4 Numerical Study
In this section, we study the practical performance of the NewtonSketch and Subsampled Newton (SSN) methods described in the previous sections. As a benchmark, we compare them against a first order method, and for this purpose we chose SVRG [9], a method that has gained much popularity in recent years. The main goals of our numerical investigation are to determine whether the second order methods have advantages over first order methods, to identify the classes of problems for which this is the case, and to evaluate the relative strengths of the NewtonSketch and SSN approaches.
We test two versions of the SSN method that differ in the choice of iterative linear solver, as mentioned above; we refer to them as SSNCG and SSNSGI. The steplength in (2.1) for the NewtonSketch and the SSN methods was determined by an Armijo backtracking line search, starting from a unit step length. For methods that use CG to solve the linear system, the CG iteration was terminated as soon as the following condition is satisfied,
(4.1) 
where is a tuning parameter. For SSNSGI, the SGI method was terminated after a fixed number, , of iterations. The SVRG method operates in cycles. Each cycle begins with a full gradient step followed by stochastic gradienttype iterations. SVRG has two tuning parameters: the (fixed) step length employed at every inner iteration and the length of the cycle. The cost per cycle is the sum of a full gradient evaluation plus evaluations of the component gradients . Therefore, each cycle of SVRG is comparable in cost to the NewtonSketch and SSN methods.
We test the methods on binary classification problems where the training function is given by a logistic loss with regularization:
(4.2) 
Here , , denote the training examples and . For this objective function, the cost of computing the gradient is the same as the cost of computing a Hessianvector product. When reporting the numerical performance of the methods, we make use of this fact. The data sets employed in our experiments are listed in Appendix A. They were selected to have different dimensions and condition numbers, and consist of both synthetic and public data sets.
Experimental Setup
To determine the best implementation of each method, we experimented with settings of parameters to achieve best practical performance. To do so, we first determined, for each method, the total amount of work to be performed between gradient evaluations. We call this the iteration budget and denote it by . For the second order methods, controls the cost of forming and solving the linear system in (2.1), while for the SVRG method it controls the amount of work of a complete cycle of length . We tested 9^{1}^{1}1 values of for each method and independently tested all possible combinations of the tuning parameters. For example, for the SVRG method, we set , which implies that , and tuned the step length parameter . For NewtonSketch, we chose pairs of parameters such that ; for SSNSGI, we set , and tuned for the step length parameter ; and for SSNCG, we chose pairs of parameters such that .
4.1 Numerical Results
In Figure 1 we present a sample of results using four data sets (the rest of the results are given in Appendix A). We report training error () vs. iterations, as well as training error vs. effective gradient evaluations (defined as the sum of all function evaluations, gradient evaluations and Hessianvector products). For SVRG one iteration denotes one complete cycle.
Our first observation, from the second column of Figure 1, is that the SSNSGI method is not competitive with the SSNCG method. This was not easily predictable at the outset because the ability of the SGI method to use a new Hessian sample at each iteration seems advantageous. But the results show that this is no match to the superior convergence properties of the CG method. Comparing with SVRG, we note that on the more challenging problems gisette and australian NewtonSketch and SSNCG are much faster; on problem australianscale, which is wellconditioned, all methods perform on par; on the simple problem rcv1 the benefits of using curvature information do not outweigh the increase in cost.
The NewtonSketch and SSNCG methods perform on par on most of the problems in terms of effective gradient evaluations. However, if we had reported CPU time, the SSNCG method would show a dramatic advantage over NewtonSketch. Interestingly, we observed that the optimal number of CG iterations for the best performance of the two methods is very similar. On the other hand, the best sample sizes differ significantly: the optimal choice of the number of rows of the sketch matrix in NewtonSketch (using Hadamard matrices) is significantly lower than the number of subsamples used in the SSNCG. This is explained by the eigenvalue figures in Appendix B, which show that a sketched Hessian with small number of rows is able to adequately capture the eigenvalues of the true Hessian on most test problems.
More extensive results and commentary are provided in the Appendix.
5 Analysis
We have seen in the previous section that the SSNCG method is very efficient in practice. In this section, we discuss some of its theoretical properties.
A complexity analysis of the SSNCG method based on the worstcase performance of CG is given in [2]. That analysis, however, underestimates the power of the SSNCG method on problems with favorable eigenvalue distributions and yields complexity estimates that are not indicative of its practical performance. We now present a more realistic analysis of SSNCG that exploits the optimal interpolatory properties of CG on convex quadratic functions.
When applied to a symmetric positive definite linear system , the progress made by the CG method is not uniform but depends on the gap between eigenvalues. Specifically, if we denote the eigenvalues of by , then the iterates , , generated by the CG method satisfy [12]
(5.1) 
Here We make use of this property in the following local linear convergence result for the SSNCG method, which is established under standard assumptions. For simplicity, we assume that the size of the sample used to define the Hessian approximation is constant throughout the run of the SSNCG method (but the sample itself varies at every iteration).
Assumptions

(Bounded Eigenvalues of Hessians) Each function is twice continuously differentiable and each component Hessian is positive definite with eigenvalues lying in a positive interval. That is there exists positive constants , such that
(5.2) We define , which is an upper bound on the condition number of the true Hessian.

(Lipschitz Continuity of Hessian) The Hessian of the objective function is Lipschitz continuous, i.e., there is a constant such that
(5.3) 
(Bounded Variance of Hessian components) There is a constant such that, for all component Hessians, we have
(5.4) 
(Bounded Moments of Iterates) There is a constant such that for any iterate generated by the SSNCG method satisfies
(5.5)
In practice, most finitesum problems are regularized and therefore Assumption A.1 is satisfied. Assumption A.2 is a standard (and convenient) in the analysis of Newtontype methods. Concerning Assumption A.3, variances are always bounded for the finitesum problem. Assumption A.4 is imposed on nonnegative numbers and is less restrictive than assuming that the iterates are bounded.
Assumption A1 is relaxed in [18] so that only the full Hessian is assumed positive definite, and not the Hessian of each individual component, . However, [18] only provides periteration improvement in probability, and hence, convergence of the overall method is not established. Rather than dealing with these intricacies, we rely on Assumption A1 to focus instead on the interaction between the CG linear solver and the efficiency of method.
We now state the main theorem of the paper. We follow the proof technique in [2, Theorem 3.2] but rather than assuming the worstcase behavior of the CG method, we give a more realistic bound that exploits the properties of the CG iteration. Let, denote the eigenvalues of subsampled Hessian . In what follows, denotes total expectation and conditional expectation at the th iteration.
Theorem 5.1.
Suppose that Assumptions A.1–A.4 hold. Let be the iterates generated by the Subsampled NewtonCG method (SSNCG) with , and
(5.6) 
and suppose that the number of CG steps performed at iteration is the smallest integer such that
(5.7) 
Then, if , we have
(5.8) 
Proof.
We write
(5.9) 
Term 1 was analyzed in [2], which establishes the bound
(5.10) 
Using the result in Lemma 2.3 from [2] and (5.6), we have
(5.11) 
Now, we analyze Term 2, which is the residual after iterations of the CG method. Assuming for simplicity that the initial CG iterate is , we obtain from (5.1)
where . To express this in terms of unweighted norms, note that if , then
Therefore, from Assumption A1 and due to the condition (5.7) on the number of CG iterations, we have
(5.12) 
where the last equality follows from the definition .
Thus, if the initial iterate is near the solution , then the SSNCG method converges to at a linear rate, with a convergence constant of . As we discuss below, when the spectrum of the Hessian approximation is favorable, condition (5.7) can be satisfied after a small number () of CG iterations.
We now establish a complexity result for the SSNCG method that bounds the number of component gradient evaluations () and component Hessianvector products () needed to obtain an accurate solution. In establishing this result we assume that the cost of the product is the same as the cost of evaluating a component gradient; this is the case in many practical applications such as the problem tested in the previous section.
Corollary 5.2.
Suppose that Assumptions A.1–A.4 hold, and let and be defined as in Theorem 5.1. Then, the total work required to get an accurate solution is
(5.14) 
This result involves the quantity , which is the maximum number of CG iterations performed at any iteration of the SSNCG method. This quantity depends on the eigenvalue distribution of the subsampled Hessians. For various eigenvalue distributions, such as clustered eigenvalues, can be much smaller than . Such favorable eigenvalue distributions can be observed in many applications, including machine learning problems. Figure 2 depicts the spectrum of for 4 data sets used in the logistic regression experiments discussed in the previous section.
For gisette, MNIST and cina one can observe that a small number of CG iterations can give rise to a rapid improvement in accuracy in the solution of the linear system. For example, for cina, the CG method improves the accuracy in the solution by 2 orders of magnitude during CG iterations 2024. For gisette, CG exhibits a fast initial increase in accuracy, after which it is not cost effective to continue performing CG steps due to the spectral distribution of this problem. The least favorable case is given by the synthetic3 data set, which was constructed to have a high condition number and a wide spectrum of eigenvalues. The behavior of CG on this problem is similar to that predicted by its worstcase analysis. Additional plots are given Appendix B.
Let us now consider the NewtonSketch algorithm, evaluate the total work complexity to get an accurate solution and compare with the work complexity of the SSNCG method. Pilanci and Wainwright [16, Theorem 1] provide a linearquadratic convergence rate (stated in probability) for the exact NewtonSketch method, similar in form to (5.9)(5), with Term 2 in (5.9) equal to zero since linear systems are solved exactly. In order to achieve local linear convergence, the usersupplied parameter should be . As a result, the sketch dimension is given by
where the second equality follows from the fact that the square of the Gaussian width is at most [16]. The periteration cost of NewtonSketch with a randomized Hadamard transform is
Hence, the total work complexity required to get an accurate solution is
(5.15) 
Note that in order to compare this with (5.14), we need to bound the parameters and . Using the definitions given in the Assumptions, we have that is bounded by a multiple of , and we know that CG converges in at most iterations, thus . In the most pessimistic case CG requires iterations, and then (5.15) and (5.14) are similar, with NewtonSketch having a slightly higher cost compared to SSNCG.
In the analysis stated above, we made certain implicit assumptions the algorithms and the problems. In the SSN method, it was assumed that the number of subsamples is less than the number of examples ; and in NewtonSketch, the sketch dimension is assumed to be less than . This implies that for these methods we made the implicit assumption that the number of data points is large enough so that .
6 Final Remarks
The theoretical properties of NewtonSketch and subsampled Hessian Newton methods have recently received much attention, but their relative performance on largescale applications has not been studied in the literature. In this paper, we focus on the finitesum problem and observed that, although a sketched matrix based on Hadamard matrices has more information than a subsampled Hessian, and approximates the spectrum of the true Hessian more accurately, the two methods perform on par in our tests, as measured by effective gradient evaluations. However, in terms of CPU time, the subsampled Hessian Newton method is much faster than the NewtonSketch method, which has a high periteration cost associated with the construction of the linear system (2.2). A Subsampled Newton method that employs the CG algorithm to solve the linear system (the SSNCG method) is particularly appealing as one can coordinate the amount of second order information employed at every iteration with the accuracy of the linear solver, and find an efficient implementation for a given problem. The power of the CG algorithm as an inner iterative solver is evident in our tests and is quantified in the complexity analysis presented in this paper. In particular, we find that CG is a more efficient solver than a first order stochastic gradient iteration.
Our numerical tests show that the two second order methods studied in this paper are generally more efficient than SVRG (a popular first order method), on some logistic regression problems
References
 [1] Naman Agarwal, Brian Bullins, and Elad Hazan. Secondorder stochastic optimization for machine learning in linear time. Journal of Machine Learning Research, 18(116):1–40, 2017.
 [2] Raghu Bollapragada, Richard H Byrd, and Jorge Nocedal. Exact and inexact subsampled newton methods for optimization. IMA Journal of Numerical Analysis, 39(2):545–578, 2018.
 [3] Léon Bottou, Frank E Curtis, and Jorge Nocedal. Optimization methods for largescale machine learning. stat, 1050:2, 2017.
 [4] ChihChung Chang and ChihJen Lin. LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2:27:1–27:27, 2011. Software available at http://www.csie.ntu.edu.tw/ cjlin/libsvm.
 [5] Petros Drineas, Michael W Mahoney, S Muthukrishnan, and Tamás Sarlós. Faster least squares approximation. Numerische Mathematik, 117(2):219–249, 2011.
 [6] Murat A Erdogdu and Andrea Montanari. Convergence rates of subsampled newton methods. In Advances in Neural Information Processing Systems 28, pages 3034–3042, 2015.
 [7] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore, second edition, 1989.
 [8] Isabelle Guyon, Constantin F Aliferis, Gregory F Cooper, André Elisseeff, JeanPhilippe Pellet, Peter Spirtes, and Alexander R Statnikov. Design and analysis of the causation and prediction challenge. In WCCI Causation and Prediction Challenge, pages 1–33, 2008.
 [9] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems 26, pages 315–323, 2013.
 [10] Felix Krahmer and Rachel Ward. New and improved johnson–lindenstrauss embeddings via the restricted isometry property. SIAM Journal on Mathematical Analysis, 43(3):1269–1281, 2011.
 [11] Yann LeCun, LD Jackel, Léon Bottou, Corinna Cortes, John S Denker, Harris Drucker, Isabelle Guyon, UA Muller, E Sackinger, Patrice Simard, et al. Learning algorithms for classification: A comparison on handwritten digit recognition. Neural networks: the statistical mechanics perspective, 261:276, 1995.
 [12] David G. Luenberger. Linear and Nonlinear Programming. AddisonWesley, New Jersey, second edition, 1984.
 [13] Haipeng Luo, Alekh Agarwal, Nicolò CesaBianchi, and John Langford. Efficient second order online learning by sketching. In Advances in Neural Information Processing Systems 29, pages 902–910, 2016.
 [14] Indraneel Mukherjee, Kevin Canini, Rafael Frongillo, and Yoram Singer. Parallel boosting with momentum. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 17–32. Springer Berlin Heidelberg, 2013.
 [15] Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013.
 [16] Mert Pilanci and Martin J Wainwright. Newton sketch: A near lineartime optimization algorithm with linearquadratic convergence. SIAM Journal on Optimization, 27(1):205–245, 2017.
 [17] H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical Statistics, pages 400–407, 1951.
 [18] Farbod RoostaKhorasani and Michael W Mahoney. Subsampled newton methods. Mathematical Programming, pages 1–34.
 [19] Suvrit Sra, Sebastian Nowozin, and Stephen J Wright. Optimization for machine learning. 2012.
 [20] Jialei Wang, Jason D Lee, Mehrdad Mahdavi, Mladen Kolar, Nathan Srebro, et al. Sketching meets random projection in the dual: A provable recovery algorithm for big and highdimensional data. Electronic Journal of Statistics, 11(2):4896–4944, 2017.
 [21] Shusen Wang, Alex Gittens, and Michael W Mahoney. Sketched ridge regression: Optimization perspective, statistical perspective, and model averaging. Journal of Machine Learning Research, 18(218):1–50, 2018.
 [22] David P Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends® in Theoretical Computer Science, 10(1–2):1–157, 2014.
 [23] Stephen J Wright. Coordinate descent algorithms. Mathematical Programming, 151(1):3–34, 2015.
 [24] Peng Xu, Farbod RoostaKhorasan, and Michael W Mahoney. Secondorder optimization for nonconvex machine learning: An empirical study. arXiv preprint arXiv:1708.07827, 2017.
 [25] Peng Xu, Farbod RoostaKhorasani, and Michael W Mahoney. Newtontype methods for nonconvex optimization under inexact hessian information. arXiv preprint arXiv:1708.07164, 2017.
 [26] Peng Xu, Jiyan Yang, Farbod RoostaKhorasani, Christopher Ré, and Michael W Mahoney. Subsampled newton methods with nonuniform sampling. In Advances in Neural Information Processing Systems 29, pages 3000–3008, 2016.
 [27] Zhewei Yao, Peng Xu, Farbod RoostaKhorasani, and Michael W Mahoney. Inexact nonconvex newtontype methods. arXiv preprint arXiv:1802.06925, 2018.
Appendix A Complete Numerical Results
Dataset  (train; test)  

synthetic1  (9000; 1000)  100  
synthetic2  (9000; 1000)  100  
synthetic3  (90000; 10000)  100  
synthetic4  (90000; 10000)  100 
Dataset  (train; test)  Source  

australian  (621; 69)  14  ,  [4] 
reged  (450; 50)  999  [8]  
mushroom  (5,500; 2,624)  112  [4]  
ijcnn1  (35,000; 91,701)  22  [4]  
cina  (14,430; 1603)  132  [8]  
ISOLET  (7,017; 780)  617  [4]  
gisette  (6,000; 1,000)  5,000  [4]  
cov  (522,911; 58101)  54  [4]  
MNIST  (60,000; 10,000)  748  [11]  
rcv1  (20,242; 677,399)  47,236  [4]  
realsim  (65,078; 7,231)  20,958  [4] 
The synthetic data sets were generated randomly as described in [14]. These data sets were created to have different number of samples () and different condition numbers (), as well as a wide spectrum of eigenvalues (see Figures 1114). We also explored the performance of the methods on popular machine learning data sets chosen to have different number of samples (), different number of variables () and different condition numbers (). We used the testing data sets where provided. For data sets where a testing set was not provided, we randomly split the data sets () for training and testing, respectively.
We focus on logistic regression classification; the objective function is given by
(A.1) 
where denote the training examples and is the regularization parameter.
For the experiments in Sections A.1, A.2 and A.3 (Figures 310), we ran four methods: SVRG, NewtonSketch, SSNSGI and SSNCG. The implementation details of the methods are given in Sections 3 and 4. In Section A.1 we show the performance of the methods on the synthetic data sets, in Section A.2 we show the performance of the methods on popular machine learning data sets and in Section A.3 we examine the sensitivity of the methods on the choice of the hyperparameters. We report training error vs. iterations, and training error vs. effective gradient evaluations (defined as the sum of all function evaluations, gradient evaluations and Hessianvector products). In Sections A.1 and A.2 we also report testing loss vs. effective gradient evaluations.
a.1 Performance of methods  Synthetic Data sets
Figure 3 illustrates the performance of the methods on four synthetic data sets. These problems were constructed to have high condition numbers , which is the setting where Newton methods show their strength compared to a first order methods such as SVRG, and to have wide spectrum of eigenvalues, which is the setting that is unfavorable for the CG method.
As is clear from Figure 3, the NewtonSketch and SSNCG methods outperform the SVRG and SSNSGI methods. This is expected due to the illconditioning of the problems. The NewtonSketch method performs slightly better than the SSNCG method; this can be attributed, primarily, to the fact that the component function in these problems are highly dissimilar. The Hessian approximations constructed by the NewtonSketch method, use information from all component functions, and thus, better approximate the true Hessian. It is interesting to observe that in the initial stage of these experiments, the SVRG method outperforms the second order methods, however, the progress made by the first order method quickly stagnates.
a.2 Performance of methods  Machine Learning Data sets
Figures 4–6 illustrate the performance of the methods on 12 popular machine learning data sets. As mentioned in Section 4.1, the performance of the methods is highly dependent on the problem characteristics. On the one hand, these figures show that there exists an SVRG sweetspot, a regime in which the SVRG method outperforms all stochastic second order methods investigated in this paper; however, there are other problems for which it is efficient to incorporate some form of curvature information in order to avoid stagnation due to illconditioning. For such problems, the SVRG method makes slow progress due to the need for a very small step length, which is required to ensure that the method does not diverge.
Note: We did not run NewtonSketch on the cov data set (Figure 6). This was due to the fact that the number of data points () in this data set is large, and so it was prohibitive (in terms of memory) to create the padded sketch matrices.
a.3 Sensitivity of methods
In Figures 7–10 we illustrate the sensitivity of the methods to the choice of hyperparameters for 4 different data sets: synthetic1, australianscale, mushroom and ijcnn. One needs to be careful when interpreting the sensitivity results. It appears that all methods have similar sensitivity to the chosen hyperparameters, however, we should note that this is not the case. More specifically, if the step length is not chosen appropriately in the SVRG and SSNSGI methods, these methods can diverge ( too large) or make very slow progress ( too small). We have excluded these runs from the figures presented in this section. Empirically, we observe that the NewtonSketch and SSNCG methods are more robust to the choice of the hyperparameters, and easier to tune. For almost all choices of and , respectively, and the methods converge, with the choice of hyperparameters only affecting the speed of convergence.
Appendix B Eigenvalues
In this Section, we present the eigenvalue spectrum for different data sets and different subsample sizes. Here denotes the number of subsamples used in the subsampled Hessian and denotes the number of rows of the Sketch matrix. In order to make a fair comparison, we chose for each figure.
To calculate the eigenvalues we used the following procedure. For each data set and subsample size we:

Computed the true Hessian
of (4.2) at , and computed the eigenvalues of the true Hessian matrix (red).

Computed 10 different subsampled Hessians of (4.2) at the optimal point using different samples for ()
computed the eigenvalues of each subsampled Hessian and took the average of the eigenvalues across the 10 replications (blue).

Computed 10 different sketched Hessians of (4.2) at the optimal point using different sketch matrices for
computed the eigenvalues of each sketched Hessian and took the average of the eigenvalues across the 10 replications (green).
All eigenvalues were computed in MATLAB using the function eig. Figures 11–22 show the distribution of the eigenvalues for different data sets. Each figure represents one data set and depicts the eigenvalue distribution for three different subsample and sketch sizes. The red marks, blue squares and green circles represent the eigenvalues of the true Hessian, the average eigenvalues of the subsampled Hessians and the average eigenvalues of the sketched Hessians, respectively. Since we computed 10 different subsampled and sketched Hessians, we also show the upper and lower bounds of each eigenvalue with error bars.
As is clear from the figures below, contingent on a reasonable subsample size and sketch size, the eigenvalue spectra of the subsampled and sketched Hessians are able to capture the eigenvalue spectrum of the true Hessian. It appears however, that the eigenvalue spectra of the sketched Hessians are closer to the spectra of the true Hessian. By this we mean both that with fewer rows of the sketch matrix, the approximations of the average eigenvalues of the sketched Hessians are closer to the true eigenvalues of the Hessian, and that the error bars of the sketched Hessians are smaller than those of the subsampled Hessians. This is not surprising as the sketched Hessians use information from all components functions whereas the subsampled Hessians are constructed from a subset of the component functions. This is interesting in itself, however, the main reasons we present these results is to demonstrate that the eigenvalue distributions of machine learning problems appear to have favorable distributions for CG.