An Investigation of Newton-Sketch and Subsampled Newton Methods

TL;DR

This paper compares two sketching techniques, Hessian subsampling and randomized Hadamard transformations, for improving Newton's method in large-scale finite-sum optimization problems, analyzing their tradeoffs and practical performance.

Contribution

It provides a comparative analysis of Hessian subsampling and randomized Hadamard sketching in Newton's method, including complexity and practical implementation insights.

Findings

Conjugate gradient method outperforms stochastic gradient in this context.

Complexity analysis favors Hessian subsampling under certain conditions.

Numerical experiments demonstrate tradeoffs between the two sketching 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 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.

Tables

Table 1. Table 1: Synthetic Datasets [ 14 ] .

Dataset	$n$ (train; test)	$d$	$\kappa$
synthetic1	(9000; 1000)	100	${10}^2$
synthetic2	(9000; 1000)	100	${10}^4$
synthetic3	(90000; 10000)	100	${10}^4$
synthetic4	(90000; 10000)	100	${10}^5$

Equations

$$\underset{w\in {\mathbb{R}}^d}{\min }F(w)=\frac{1}{n}\sum _{i=1}^n{F}_i(w),$$

$$A_k{p}_k=-\nabla F(w_k),w_{k+1}=w_k+{\alpha }_k{p}_k,$$

$$A_{k}=\Big{(}(S^{k}\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}})^{T}S^{k}\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}}\Big{)}.$$

$$A_k={\nabla }^2{F}_{T_k}(w_k)=\frac{1}{\left|T_k\right|}\sum _{i\in T_k}{\nabla }^2{F}_i(w_k),$$

$$A_k{p}_k=-\nabla F(w_k).$$

$$\displaystyle v=\Big{(}(S\nabla^{2}F(w)^{\nicefrac{{1}}{{2}}})^{T}S\nabla^{2}F(w)^{\nicefrac{{1}}{{2}}}\Big{)}u$$

$$v_1=(S{\nabla }^{2}F(w{)}^{\text{nicefrac}12})u,\text{and }v=(S{\nabla }^{2}F(w{)}^{\text{nicefrac}12}{)}^T{v}_1.$$

$$v=A_{k}u={\nabla }^2{F}_{T_k}(w_k)u=\frac{1}{|T_k|}\sum _{i\in T_k}{\nabla }^2{F}_i(w_k)u.$$

$$Q_i(p)=F(w_k)+\nabla F(w_k{)}^{T}p+\frac{1}{2}{p}^T{\nabla }^2{F}_i(w_k)p.$$

$$p_k^{t+1}=p_k^t-\bar{\alpha }\nabla Q_i(p_k^t)=(I-\bar{\alpha }{\nabla }^2{F}_i(w_k))p_k^t-\bar{\alpha }\nabla F(w_k),\text{mbox}with{p}_k^0=-\nabla F(w_k).$$

$$\Vert {\nabla }^2{F}_{T_k}(w_k)p_k+\nabla F(w_k)\Vert <\zeta \Vert \nabla F(w_k)\Vert ,$$

$$F(w)=\frac{1}{n}\sum _{i=1}^n\log (1+e^{-y^i(w^T{x}^i)})+\frac{\lambda }{2}\Vert w{\Vert }^2.$$

$$||p^r-p^{\ast }|{|}_A^2\le {\left(\frac{{\lambda }_{d-r+1}-{\lambda }_1}{{\lambda }_{d-r+1}+{\lambda }_1}\right)}^2\Vert p^0-p^{\ast }|{|}_A^2,\text{mbox}whereA{p}^{\ast }=b.$$

$$\mu I⪯{\nabla }^2{F}_i(w)⪯LI,\forall w\in {\mathbb{R}}^d.$$

$$\Vert {\nabla }^{2}F(w)-{\nabla }^{2}F(z)\Vert \le M\Vert w-z\Vert ,\forall w,z\in {\mathbb{R}}^d.$$

$$\Vert \mathbb{E}[({\nabla }^2{F}_i(w)-{\nabla }^{2}F(w){)}^2]\Vert \le {\sigma }^2,\forall w\in {\mathbb{R}}^d.$$

$$\mathbb{E}[\Vert w_k-w^{\ast }{\Vert }^2]\le \gamma (\mathbb{E}[\Vert w_k-w^{\ast }\Vert ]{)}^2.$$

$$|T_k|=T\ge \frac{64{\sigma }^2}{{\mu }^2},$$

$$\left(\frac{{\lambda }_{d-r_k+1}^{T_k}-{\lambda }_1^{T_k}}{{\lambda }_{d-r_k+1}^{T_k}+{\lambda }_1^{T_k}}\right)\le \frac{1}{8{\kappa }^{3/2}}.$$

$$\mathbb{E}\left[||w_{k+1}-w^{\ast }||\right]\le \frac{1}{2}\mathbb{E}\left[||w_k-w^{\ast }||\right],k=0,1,\dots .$$

$${\mathbb{E}}_k[\Vert w_{k+1}-w^{\ast }\Vert ]$$

$$\le \underset{\text{Term 1}}{\underbrace{{\mathbb{E}}_k[\Vert w_k-w^{\ast }-{\nabla }^2{F}_{T_k}^{-1}(w_k)\nabla F(w_k)\Vert ]}}+\underset{\text{Term 2}}{\underbrace{{\mathbb{E}}_k[\Vert p_k^r+{\nabla }^2{F}_{T_k}^{-1}(w_k)\nabla F(w_k)\Vert ]}}.$$

$${\mathbb{E}}_k[\Vert w_k-w^{\ast }-{\nabla }^2{F}_{T_k}^{-1}(w_k)\nabla F(w_k)\Vert ]\le$$

$$+\frac{1}{\mu }{\mathbb{E}}_k\left[\Vert \left({\nabla }^2{F}_{T_k}(w_k)-{\nabla }^{2}F(w_k)\right)(w_k-w^{\ast })\Vert \right].$$

$${\mathbb{E}}_k[\Vert w_k-w^{\ast }-{\nabla }^2{F}_{T_k}^{-1}(w_k)\nabla F(w_k)\Vert ]$$

$$\le \frac{M}{2\mu }\Vert w_k-w^{\ast }{\Vert }^2+\frac{1}{8}\Vert w_k-w^{\ast }\Vert .$$

$$\Vert p_k^r+{\nabla }^2{F}_{T_k}^{-1}\nabla F(w_k){\Vert }_{A_k}$$

$${\lambda }_1\Vert a{\Vert }^2\le a^{T}Aa\le b^{T}Ab\le {\lambda }_d\Vert b{\Vert }^2⟹\Vert a\Vert \le \sqrt{\kappa (A)}\Vert b\Vert \le \sqrt{\kappa }\Vert b\Vert .$$

$$\Vert p_k^r+{\nabla }^2{F}_{T_k}^{-1}(w_k)\nabla F(w_k)\Vert$$

$$\le \sqrt{\kappa }\frac{1}{8{\kappa }^{3/2}}\Vert \nabla F(w_k)\Vert \Vert {\nabla }^2{F}_{T_k}^{-1}(w_k)\Vert$$

$$\le \frac{L}{\mu }\frac{1}{8\kappa }\Vert w_k-w^{\ast }\Vert$$

$$=\frac{1}{8}\Vert w_k-w^{\ast }\Vert ,$$

$${\mathbb{E}}_k[\Vert w_{k+1}-w^{\ast }\Vert ]$$

Full text

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,

[TABLE]

in which each component function $F_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ 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 $d$ and the number of data points $n$ 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 $d$ is not large but the number of data points $n$ 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 $\alpha_{k}$, 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,

[TABLE]

where $\alpha_{k}$ is the steplength parameter and $A_{k}\succ 0$ 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 $\nabla F(w_{k})$ is exact, and therefore, the stochastic nature of the iteration is entirely due to the the choice of $A_{k}$. 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 $Q$-linear convergence rate, which makes them directly comparable with first order methods.

We now discuss how the Hessian approximations $A_{k}$ are constructed.

Hessian-Sketch

As proposed by Pilanci and Wainwright [16], the matrix $A_{k}$ 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 $S^{k}\in\mathbb{R}^{m\times n}$ with the property that $\mathbb{E}\left[{(S^{k})^{T}S^{k}}/{m}\right]=I_{n}$, and $m<n$. Next, the algorithm computes a square-root decomposition of the Hessian $\nabla^{2}F(w_{k})$, which we denote by $\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}}\in\mathbb{R}^{n\times d}$, and defines the Hessian approximation as

[TABLE]

The search direction of the algorithm is obtained by solving the linear system in (2.1). When $m<<n$ forming $A_{k}$ is much less expensive than forming the Hessian $\nabla^{2}F(w_{k})$.

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 $F(w)=\sum_{i=1}^{n}g_{i}(w;x_{i})$. The square-root matrix can be expressed as $\nabla^{2}F(w)^{\nicefrac{{1}}{{2}}}=D^{\nicefrac{{1}}{{2}}}X$, where $X\in\mathbb{R}^{n\times d}$ is the data matrix, $x_{i}\in\mathbb{R}^{1\times d}$ denotes the $i$-th row of $X$, and $D$ is a diagonal matrix given by $D={\rm diag}\{g_{i}^{\prime\prime}(w;x_{i})^{\nicefrac{{1}}{{2}}}\}_{i=1}^{n}$. The multiplication of the sketch matrix $S^{k}$ with $\nabla^{2}F(w)^{\nicefrac{{1}}{{2}}}$ 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 $H\in\mathbb{R}^{n\times n}$ as having entries $|H_{i,j}|\in\left[\frac{-1}{\sqrt{n}},\frac{1}{\sqrt{n}}\right]$. The randomized Hadamard projections are performed using the sketch matrix $S^{k}\in\mathbb{R}^{m\times n}$ whose rows are chosen by randomly sampling $m$ rows of the matrix $\sqrt{n}HD$, where $D\in\mathbb{R}^{n\times n}$ is a diagonal matrix with random $\pm{1}$ entries. The properties of this approach have been analyzed in [16], where it is shown that each iteration has complexity $\mathcal{O}(nd\log(m)+dm^{2})$, which is lower than the complexity of Newton’s method ($\mathcal{O}(nd^{2})$), when $n>d$ and $m=\mathcal{O}(d)$.

Hessian Subsampling

Hessian subsampling is a special case of sketching where the Hessian approximation is constructed by randomly sampling component Hessians $\nabla^{2}F_{i}$ of (1.1), and averaging them. Specifically, at the $k$-th iteration, one defines

[TABLE]

where the index set $T_{k}\subset\{1,...,n\}$ 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 $A_{k}$ 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 $A_{k}$ 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 $F_{i}$ 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 $d$ 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 $d$ is not large. Since we do not wish to impose such a restriction on $d$, 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 $p_{k}$ by running $r\leq d$ iterations of CG on the system

[TABLE]

Under a favorable distribution of the eigenvalues of the matrix $A_{k}$, the CG method can generate a useful approximate solution in much fewer than $d$ 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 $n\times d$ matrix, denoted by $\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}}$, such that $[\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}}]^{T}\nabla^{2}F(w_{k})^{\nicefrac{{1}}{{2}}}=\nabla^{2}F(w_{k})$. 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 $d$ is large it is imperative to use an iterative method such as CG. The entire matrix $A_{k}$ 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 $u\in\mathbb{R}^{d}$ one can compute

[TABLE]

in two steps as follows

[TABLE]

Our implementation of the Newton-Sketch method includes two tuning parameters: the number of rows $m_{\text{NS}}$ in the sketch matrix $S$, and the maximum number of CG iterations, $\max_{\text{CG}}$, used to compute the step. We define the sketch matrix $S$ through the Hadamard basis that allows for a fast matrix-vector multiplication, to obtain $S\nabla^{2}F(w)^{\nicefrac{{1}}{{2}}}$. This does not require forming the sketch matrix explicitly and can be executed at a cost of $\mathcal{O}(m_{\text{NS}}d\log n)$ flops. The maximum cost of every iteration is given as the evaluation of the gradient $\nabla F(w_{k})$ plus $2m_{\text{NS}}\times max_{\text{CG}}$ 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 $A_{k}$ times vectors can be computed by simply summing the individual Hessian-vector products. That is, given a vector $u\in\mathbb{R}^{d}$,

[TABLE]

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, $T$, and the maximum number of CG iterations used to compute the step, $\max_{\text{CG}}$. Thus, the maximum cost of a SSN-CG iteration is given by the cost of evaluating the gradient $\nabla F(w_{k})$ plus $T\times\max_{\text{CG}}$ 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 $m_{\text{SGI}}$. At every iteration of the cycle, an index $i\in\{1,2,...,n\}$ is chosen at random, and a gradient step is computed for the following quadratic model at $w_{k}$:

[TABLE]

This yields the iteration

[TABLE]

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 $\nabla F(w_{k})$ plus $m_{\text{SGI}}$ Hessian-vector products of the form $\nabla^{2}F_{i}v$. This method has two tuning parameters: the inner iteration step length $\bar{\alpha}$ and the length $m_{\text{SGI}}$ 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 $\alpha_{k}$ 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,

[TABLE]

where $\zeta$ is a tuning parameter. For SSN-SGI, the SGI method was terminated after a fixed number, $m_{\text{SGI}}$, of iterations. The SVRG method operates in cycles. Each cycle begins with a full gradient step followed by $m_{\text{SVRG}}$ stochastic gradient-type iterations. SVRG has two tuning parameters: the (fixed) step length employed at every inner iteration and the length $m_{\text{SVRG}}$ of the cycle. The cost per cycle is the sum of a full gradient evaluation plus $2m_{\text{SVRG}}$ evaluations of the component gradients $\nabla F_{i}$. 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 $\ell_{2}$ regularization:

[TABLE]

Here $(x^{i},y^{i})$, $i=1,\ldots,n$, denote the training examples and $\lambda=\frac{1}{n}$. 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 $b$. For the second order methods, $b$ 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 $m_{\text{SVRG}}$. We tested 9111$b=\{n/100,n/50,n/10,n/5,n/2,n,2n,5n,10n\}$ values of $b$ for each method and independently tested all possible combinations of the tuning parameters. For example, for the SVRG method, we set $b=n$, which implies that $m_{\text{SVRG}}=\frac{b}{2}=\frac{n}{2}$, and tuned the step length parameter $\alpha$. For Newton-Sketch, we chose pairs of parameters $(m_{\text{NS}},max_{\text{CG}})$ such that $\frac{m_{\text{NS}}\times max_{\text{CG}}}{2}=b$; for SSN-SGI, we set $m_{\text{SGI}}=b=n$, and tuned for the step length parameter $\alpha$; and for SSN-CG, we chose pairs of parameters $(T,max_{\text{CG}})$ such that $T\times max_{\text{CG}}=b$.

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 ($F(w_{k})-F^{*}$) 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 $\nabla^{2}F_{i}$ 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 $m_{\text{NS}}$ in Newton-Sketch (using Hadamard matrices) is significantly lower than the number of subsamples $T$ 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 $Ap=b$, the progress made by the CG method is not uniform but depends on the gap between eigenvalues. Specifically, if we denote the eigenvalues of $A$ by $0<\lambda_{1}\leq\lambda_{2}\cdots\leq\lambda_{d}$, then the iterates $\{p^{r}\}$, $r=1,2,\ldots,d$, generated by the CG method satisfy [12]

[TABLE]

Here $\|p\|^{2}_{A}=p^{T}Ap.$ 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 $|T_{k}|=T$ used to define the Hessian approximation is constant throughout the run of the SSN-CG method (but the sample $T_{k}$ itself varies at every iteration).

Assumptions

1. A.1

(Bounded Eigenvalues of Hessians) Each function $F_{i}$ is twice continuously differentiable and each component Hessian is positive definite with eigenvalues lying in a positive interval. That is there exists positive constants $\mu$, $L$ such that

[TABLE]

We define $\kappa={L}/{\mu}$, which is an upper bound on the condition number of the true Hessian. 2. A.2

(Lipschitz Continuity of Hessian) The Hessian of the objective function $F$ is Lipschitz continuous, i.e., there is a constant $M>0$ such that

[TABLE] 3. A.3

(Bounded Variance of Hessian components) There is a constant $\sigma$ such that, for all component Hessians, we have

[TABLE] 4. A.4

(Bounded Moments of Iterates) There is a constant $\gamma>0$ such that for any iterate $w_{k}$ generated by the SSN-CG method satisfies

[TABLE]

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 $\nabla^{2}F$ is assumed positive definite, and not the Hessian of each individual component, $\nabla^{2}F_{i}$. 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, $0<\lambda^{T_{k}}_{1}\leq\lambda^{T_{k}}_{2}\leq\cdots\leq\lambda^{T_{k}}_{d}$ denote the eigenvalues of subsampled Hessian $\nabla^{2}F_{T_{k}}(w_{k})$. In what follows, $\mathbb{E}$ denotes total expectation and $\mathbb{E}_{k}$ conditional expectation at the $k$th iteration.

Theorem 5.1**.**

Suppose that Assumptions A.1–A.4 hold. Let $\{w_{k}\}$ be the iterates generated by the Subsampled Newton-CG method (SSN-CG) with $\alpha_{k}=1$, and

[TABLE]

and suppose that the number of CG steps performed at iteration $k$ is the smallest integer $r_{k}$ such that

[TABLE]

Then, if $||w_{0}-w^{*}||\leq\min\left\{\frac{\mu}{2M},\frac{\mu}{2\gamma M}\right\}$, we have

[TABLE]

Proof.

We write

[TABLE]

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

[TABLE]

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

[TABLE]

Now, we analyze Term 2, which is the residual after $r_{k}$ iterations of the CG method. Assuming for simplicity that the initial CG iterate is $p^{0}_{k}=0$, we obtain from (5.1)

[TABLE]

where $A_{k}=\nabla^{2}F_{T_{k}}(w_{k})$. To express this in terms of unweighted norms, note that if $\|a\|^{2}_{A}\leq\|b\|^{2}_{A}$, then

[TABLE]

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

[TABLE]

where the last equality follows from the definition $\kappa={L}/{\mu}$.

Combining (5.11) and (5.12), we get

[TABLE]

We now use induction to prove (5.8). For the base case we have, by the assumption on $\|w_{0}\|$ ,

[TABLE]

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

[TABLE]

∎

Thus, if the initial iterate is near the solution $w^{*}$, then the SSN-CG method converges to $w^{*}$ at a linear rate, with a convergence constant of $1/2$. As we discuss below, when the spectrum of the Hessian approximation is favorable, condition (5.7) can be satisfied after a small number ($r_{k}$) of CG iterations.

We now establish a complexity result for the SSN-CG method that bounds the number of component gradient evaluations ($\nabla F_{i}$) and component Hessian-vector products ($\nabla^{2}F_{i}v$) needed to obtain an $\epsilon$-accurate solution. In establishing this result we assume that the cost of the product $\nabla^{2}F_{i}v$ 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 $T$ and $r_{k}$ be defined as in Theorem 5.1. Then, the total work required to get an $\epsilon$-accurate solution is

[TABLE]

This result involves the quantity $\bar{r}$, 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, $\bar{r}$ can be much smaller than $d$. Such favorable eigenvalue distributions can be observed in many applications, including machine learning problems. Figure 2 depicts the spectrum of $\nabla^{2}F_{T}(w^{*})$ 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 $\epsilon$-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 $\epsilon$ should be $\mathcal{O}({1}/{\kappa})$. As a result, the sketch dimension $m$ is given by

[TABLE]

where the second equality follows from the fact that the square of the Gaussian width $W$ is at most $\min\{n,d\}$ [16]. The per-iteration cost of Newton-Sketch with a randomized Hadamard transform is

[TABLE]

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

[TABLE]

Note that in order to compare this with (5.14), we need to bound the parameters $\sigma$ and $\bar{r}$. Using the definitions given in the Assumptions, we have that $\sigma^{2}$ is bounded by a multiple of $L^{2}$, and we know that CG converges in at most $d$ iterations, thus $\bar{r}\leq d$. In the most pessimistic case CG requires $d$ 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 $n$; and in Newton-Sketch, the sketch dimension is assumed to be less than $n$. This implies that for these methods we made the implicit assumption that the number of data points is large enough so that $n>\kappa^{2}$.

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

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 ($n$) and different condition numbers ($\kappa$), 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 ($n$), different number of variables ($d$) and different condition numbers ($\kappa$). We used the testing data sets where provided. For data sets where a testing set was not provided, we randomly split the data sets ($90\%,10\%$) for training and testing, respectively.

We focus on logistic regression classification; the objective function is given by

[TABLE]

where $(x^{i},y^{i})_{i=1}^{n}$ denote the training examples and $\lambda=\frac{1}{n}$ 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 $10^{2}-10^{5}$, 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 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 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 ($n$) 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 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 $\alpha$ is not chosen appropriately in the SVRG and SSN-SGI methods, these methods can diverge ($\alpha$ too large) or make very slow progress ($\alpha$ 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 $m_{\text{NS}}$ and $T$, respectively, and $max_{\text{CG}}$ 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 $T$ denotes the number of subsamples used in the subsampled Hessian and $m_{\text{NS}}$ denotes the number of rows of the Sketch matrix. In order to make a fair comparison, we chose $T=m_{\text{NS}}$ for each figure.

To calculate the eigenvalues we used the following procedure. For each data set and subsample size we:

Computed the true Hessian

[TABLE]

of (4.2) at $w^{*}$, and computed the eigenvalues of the true Hessian matrix (red). 2. 2.

Computed 10 different subsampled Hessians of (4.2) at the optimal point $w^{*}$ using different samples $T_{j}\subset\{1,2,\ldots,n\}$ for $j=1,2,\ldots,10$ ($T_{j}=|T|$)

[TABLE]

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

Computed 10 different sketched Hessians of (4.2) at the optimal point $w^{*}$ using different sketch matrices $S^{j}\in\mathbb{R}^{m_{\text{NS}}\times n}$ for $j=1,2,\ldots,10$

[TABLE]

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.