An Investigation of Newton-Sketch and Subsampled Newton Methods

# An Investigation of Newton-Sketch and Subsampled Newton Methods

Albert S. Berahas Department of Industrial and Systems Engineering, Lehigh University.
albertberahas@gmail.com This author was supported by the Defense Advanced Research Projects Agency (DARPA).
Raghu Bollapragada Department of Industrial Engineering and Management Sciences, Northwestern University.
raghu.bollapragada@u.northwestern.edu This author was supported by the Defense Advanced Research Projects Agency (DARPA).
Jorge Nocedal Department of Industrial Engineering and Management Sciences, Northwestern University.
j-nocedal@northwestern.edu This author was supported by the Office of Naval Research grant N00014-14-1-0313 P00003, and by National Science Foundation grant DMS-1620022
###### 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 finite-sum 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,

 minw∈RdF(w)=1nn∑i=1Fi(w), (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 trade-offs 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]. Roosta-Khorasani 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 non-convex problems. Agarwal et al. [1] proposed to use a stochastic gradient-like 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,

 Akpk=−∇F(wk),wk+1=wk+αkpk, (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.

### Hessian-Sketch

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 square-root decomposition of the Hessian , which we denote by , and defines the Hessian approximation as

 Ak=((Sk∇2F(wk)\sfrac12)TSk∇2F(wk)\sfrac12). (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 square-root 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 Johnson-Lindenstrauss (JL) transform, e.g., Hadamard matrices, which allow for a fast matrix-matrix 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

 Ak=∇2FTk(wk)=1|Tk|∑i∈Tk∇2Fi(wk), (2.3)

where the index set is chosen either uniformly at random (with or without replacement) or in a non-uniform 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

In the rest of the paper, we use the acronym SSN (Sub-Sampled-Newton) to denote algorithm (2.1) where is chosen by Hessian subsampling (2.3), and use the term Newton-Sketch when other choices of the sketch matrix are employed, as in (2.2).

### Discussion

We thus have two distinct classes of methods at our disposal. On the one hand, Newton-Sketch 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 Newton-Sketch 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 finite-sum problems involving many variables and large data sets. In order to undertake this investigation we need to give careful consideration to the large-scale implementation of these methods.

## 3 Practical Implementations of the Algorithms

We can design efficient implementations of the Newton-Sketch 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 Newton-Sketch 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 Hessian-free approach, and compute an inexact solution of linear systems with coefficient matrices (2.2) or (2.3) using an iterative method that only requires matrix-vector 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

 Akpk=−∇F(wk). (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.

#### Newton-Sketch

The linear system (2.2) arising in the Newton-Sketch method requires access to a square-root 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 matrix-vector products with respect to the sketched square-root Hessian. More specifically, for any vector one can compute

 v=((S∇2F(w)\sfrac12)TS∇2F(w)\sfrac12)u

in two steps as follows

 v1=(S∇2F(w)\sfrac12)u,and v=(S∇2F(w)\sfrac12)Tv1. (3.2)

Our implementation of the Newton-Sketch 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 matrix-vector 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 matrix-vector products, where the factor 2 arises due to (3.2). We ignore the cost of forming the square-root 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 Newton-Sketch.

#### 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 Hessian-free approach. The product of times vectors can be computed by simply summing the individual Hessian-vector products. That is, given a vector ,

 v=Aku=∇2FTk(wk)u=1|Tk|∑i∈Tk∇2Fi(wk)u. (3.3)

We refer to this method as SSN-CG method. Unlike Newton-Sketch that requires the projection of the square-root matrix onto the lower dimensional manifold, the SSN-CG method can be implemented without the need of projections, and thus does not require additional computation beyond the matrix-vector 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 SSN-CG iteration is given by the cost of evaluating the gradient plus matrix-vector 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 Newton-Sketch there is no efficient way to compute the stochastic gradient as the sketched Hessian cannot be decomposed into individual components. The resulting SSN-SGI 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 :

 Qi(p)=F(wk)+∇F(wk)Tp+12pT∇2Fi(wk)p. (3.4)

This yields the iteration

 pt+1k=ptk−¯α∇Qi(ptk)=(I−¯α∇2Fi(wk))ptk−¯α∇F(wk), with p0k=−∇F(wk). (3.5)

A method similar to SSN-SGI has been analyzed in [1], where it is referred to as LiSSA, and [2] compares the complexity of LiSSA and SSN-CG. 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 SSN-SGI method is given by the evaluation of the gradient plus Hessian-vector 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 Newton-Sketch 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 Newton-Sketch 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 SSN-CG and SSN-SGI. The steplength in (2.1) for the Newton-Sketch and the SSN methods was determined by an Armijo back-tracking 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,

 ∥∇2FTk(wk)pk+∇F(wk)∥<ζ∥∇F(wk)∥, (4.1)

where is a tuning parameter. For SSN-SGI, 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 gradient-type 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 Newton-Sketch and SSN methods.

We test the methods on binary classification problems where the training function is given by a logistic loss with regularization:

 F(w)=1nn∑i=1log(1+e−yi(wTxi))+λ2∥w∥2. (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 Hessian-vector 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 9111 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 Newton-Sketch, we chose pairs of parameters such that ; for SSN-SGI, we set , and tuned for the step length parameter ; and for SSN-CG, 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 Hessian-vector products). For SVRG one iteration denotes one complete cycle.

Our first observation, from the second column of Figure 1, is that the SSN-SGI method is not competitive with the SSN-CG 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 Newton-Sketch and SSN-CG are much faster; on problem australian-scale, which is well-conditioned, all methods perform on par; on the simple problem rcv1 the benefits of using curvature information do not outweigh the increase in cost.

The Newton-Sketch and SSN-CG methods perform on par on most of the problems in terms of effective gradient evaluations. However, if we had reported CPU time, the SSN-CG method would show a dramatic advantage over Newton-Sketch. 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 Newton-Sketch (using Hadamard matrices) is significantly lower than the number of subsamples used in the SSN-CG. 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 SSN-CG method is very efficient in practice. In this section, we discuss some of its theoretical properties.

A complexity analysis of the SSN-CG method based on the worst-case performance of CG is given in [2]. That analysis, however, underestimates the power of the SSN-CG 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 SSN-CG 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]

 ||pr−p∗||2A≤(λd−r+1−λ1λd−r+1+λ1)2∥p0−p∗||2A,% where  Ap∗=b. (5.1)

Here We make use of this property in the following local linear convergence result for the SSN-CG 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 SSN-CG method (but the sample itself varies at every iteration).

### Assumptions

1. (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

 μI⪯∇2Fi(w)⪯LI,∀w∈Rd. (5.2)

We define , which is an upper bound on the condition number of the true Hessian.

2. (Lipschitz Continuity of Hessian) The Hessian of the objective function is Lipschitz continuous, i.e., there is a constant such that

 ∥∇2F(w)−∇2F(z)∥≤M∥w−z∥,∀w,z∈Rd. (5.3)
3. (Bounded Variance of Hessian components) There is a constant such that, for all component Hessians, we have

 ∥E[(∇2Fi(w)−∇2F(w))2]∥≤σ2,∀w∈Rd. (5.4)
4. (Bounded Moments of Iterates) There is a constant such that for any iterate generated by the SSN-CG method satisfies

 E[∥wk−w∗∥2]≤γ(E[∥wk−w∗∥])2. (5.5)

In practice, most finite-sum problems are regularized and therefore Assumption A.1 is satisfied. Assumption A.2 is a standard (and convenient) in the analysis of Newton-type methods. Concerning Assumption A.3, variances are always bounded for the finite-sum problem. Assumption A.4 is imposed on non-negative 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 per-iteration 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 worst-case 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 Newton-CG method (SSN-CG) with , and

 |Tk|=T≥64σ2μ2, (5.6)

and suppose that the number of CG steps performed at iteration is the smallest integer such that

 ⎛⎜⎝λTkd−rk+1−λTk1λTkd−rk+1+λTk1⎞⎟⎠≤18κ3/2. (5.7)

Then, if , we have

 E[||wk+1−w∗||]≤12E[||wk−w∗||],k=0,1,…. (5.8)
###### Proof.

We write

 Ek[∥wk+1−w∗∥] =Ek[∥wk−w∗+prk∥] ≤Ek[∥wk−w∗−∇2F−1Tk(wk)∇F(wk)∥]Term 1+Ek[∥prk+∇2F−1Tk(wk)∇F(wk)∥]Term 2. (5.9)

Term 1 was analyzed in [2], which establishes the bound

 Ek[∥wk−w∗−∇2F−1Tk(wk)∇F(wk)∥]≤ M2μ∥wk−w∗∥2 +1μEk[∥∥(∇2FTk(wk)−∇2F(wk))(wk−w∗)∥∥]. (5.10)

Using the result in Lemma 2.3 from [2] and (5.6), we have

 Ek[∥wk−w∗−∇2F−1Tk(wk)∇F(wk)∥] ≤M2μ∥wk−w∗∥2+σμ√T∥wk−w∗∥ ≤M2μ∥wk−w∗∥2+18∥wk−w∗∥. (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)

 ∥prk+∇2F−1Tk∇F(wk)∥Ak ≤⎛⎜⎝λTkd−rk+1−λTk1λTkd−rk+1+λTk1(wk)⎞⎟⎠∥∇2F−1Tk∇F(wk)∥Ak,

where . To express this in terms of unweighted norms, note that if , then

 λ1∥a∥2≤aTAa≤bTAb≤λd∥b∥2 ⟹ ∥a∥≤√κ(A)∥b∥≤√κ∥b∥.

Therefore, from Assumption A1 and due to the condition (5.7) on the number of CG iterations, we have

 ∥prk+∇2F−1Tk(wk)∇F(wk)∥ ≤√κ⎛⎜⎝λTkd−rk+1(wk)−λTk1(wk)λTkd−rk+1(wk)+λTk1(wk)⎞⎟⎠∥∇2F−1Tk(wk)∇F(wk)∥ ≤√κ18κ3/2∥∇F(wk)∥∥∇2F−1Tk(wk)∥ ≤Lμ18κ∥wk−w∗∥ =18∥wk−w∗∥, (5.12)

where the last equality follows from the definition .

Combining (5.11) and (5.12), we get

 Ek[∥wk+1−w∗∥] ≤M2μ∥wk−w∗∥2+14∥wk−w∗∥. (5.13)

We now use induction to prove (5.8). For the base case we have, by the assumption on ,

 E[∥w1−w∗∥] ≤M2μ∥w0−w∗∥2+14∥w0−w∗∥ ≤14∥w0−w∗∥+14∥w0−w∗∥ =12∥w0−w∗∥.

Now suppose that (5.8) is true for -th iteration. Then, by Assumption A.4

 E[Ek[∥wk+1−w∗∥]] ≤M2μE[∥wk−w∗∥2]+14E[∥wk−w∗∥] ≤γM2μE[∥wk−w∗∥]E[∥wk−w∗∥]+14E[∥wk−w∗∥] ≤14E[∥wk−w∗∥]+14E[∥wk−w∗∥] ≤12E[∥wk−w∗∥].

Thus, if the initial iterate is near the solution , then the SSN-CG 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 SSN-CG method that bounds the number of component gradient evaluations () and component Hessian-vector 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 SSN-CG 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 20-24. 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 worst-case analysis. Additional plots are given Appendix B.

Let us now consider the Newton-Sketch algorithm, evaluate the total work complexity to get an -accurate solution and compare with the work complexity of the SSN-CG method. Pilanci and Wainwright [16, Theorem 1] provide a linear-quadratic convergence rate (stated in probability) for the exact Newton-Sketch 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 user-supplied parameter should be . As a result, the sketch dimension is given by

 m=O(\sfracW2ϵ2)=O(κ2min{n,d}),

where the second equality follows from the fact that the square of the Gaussian width is at most [16]. The per-iteration cost of Newton-Sketch with a randomized Hadamard transform is

 O(ndlog(m)+dm2)=O(ndlog(κd)+d3κ4).

Hence, the total work complexity required to get an -accurate solution is

 O((n+κ4d2)dlog(\sfrac1ϵ)). (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 Newton-Sketch having a slightly higher cost compared to SSN-CG.

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 Newton-Sketch, 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 Newton-Sketch and subsampled Hessian Newton methods have recently received much attention, but their relative performance on large-scale applications has not been studied in the literature. In this paper, we focus on the finite-sum 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 Newton-Sketch method, which has a high per-iteration 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 SSN-CG 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. Second-order 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 large-scale machine learning. stat, 1050:2, 2017.
• [4] Chih-Chung Chang and Chih-Jen 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 sub-sampled 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, Jean-Philippe 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. Addison-Wesley, New Jersey, second edition, 1984.
• [13] Haipeng Luo, Alekh Agarwal, Nicolò Cesa-Bianchi, 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 linear-time optimization algorithm with linear-quadratic 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 Roosta-Khorasani and Michael W Mahoney. Sub-sampled 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 high-dimensional 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 Roosta-Khorasan, and Michael W Mahoney. Second-order optimization for non-convex machine learning: An empirical study. arXiv preprint arXiv:1708.07827, 2017.
• [25] Peng Xu, Farbod Roosta-Khorasani, and Michael W Mahoney. Newton-type methods for non-convex optimization under inexact hessian information. arXiv preprint arXiv:1708.07164, 2017.
• [26] Peng Xu, Jiyan Yang, Farbod Roosta-Khorasani, Christopher Ré, and Michael W Mahoney. Sub-sampled newton methods with non-uniform sampling. In Advances in Neural Information Processing Systems 29, pages 3000–3008, 2016.
• [27] Zhewei Yao, Peng Xu, Farbod Roosta-Khorasani, and Michael W Mahoney. Inexact non-convex newton-type methods. arXiv preprint arXiv:1802.06925, 2018.

## Appendix A Complete Numerical Results

In this Section, we present further numerical results, on the data sets listed in Tables 2 and 2.

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 11-14). 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

 minw∈RdF(w)=1nn∑i=1log(1+e−yi(wTxi))+λ2∥w∥2, (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 3-10), we ran four methods: SVRG, Newton-Sketch, SSN-SGI and SSN-CG. 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 hyper-parameters. We report training error vs. iterations, and training error vs. effective gradient evaluations (defined as the sum of all function evaluations, gradient evaluations and Hessian-vector 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 Newton-Sketch and SSN-CG methods outperform the SVRG and SSN-SGI methods. This is expected due to the ill-conditioning of the problems. The Newton-Sketch method performs slightly better than the SSN-CG 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 Newton-Sketch 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 46 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 sweet-spot, 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 ill-conditioning. 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 Newton-Sketch 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 710 we illustrate the sensitivity of the methods to the choice of hyper-parameters for 4 different data sets: synthetic1, australian-scale, mushroom and ijcnn. One needs to be careful when interpreting the sensitivity results. It appears that all methods have similar sensitivity to the chosen hyper-parameters, however, we should note that this is not the case. More specifically, if the step length is not chosen appropriately in the SVRG and SSN-SGI 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 Newton-Sketch and SSN-CG methods are more robust to the choice of the hyper-parameters, and easier to tune. For almost all choices of and , respectively, and the methods converge, with the choice of hyper-parameters 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:

1. Computed the true Hessian

 ∇2F(w∗)=1nn∑i∇2Fi(w∗),

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

2. Computed 10 different subsampled Hessians of (4.2) at the optimal point using different samples for ()

 ∇2FTj(w∗)=1∣∣Tj∣∣∑i∈Tj∇2Fi(w∗),

computed the eigenvalues of each subsampled Hessian and took the average of the eigenvalues across the 10 replications (blue).

3. Computed 10 different sketched Hessians of (4.2) at the optimal point using different sketch matrices for

 ∇2FSj(w∗)=(((Sj)∇2F(w∗)\sfrac12)TSj∇2F(w∗)\sfrac12),

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 1122 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.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters