Sub-sampled Newton Methods with Non-uniform Sampling

Introduction

Consider the following optimization problem

There is a plethora of first-order optimization algorithms [Bub14, NW06] for solving (1). However, for ill-conditioned problems, it is often the case that first-order methods return a solution far from the minimizer, $\mathbf{w}^{\ast}$, albeit a low objective value. (See Figure 2 in Section 5 for example.) On the other hand, most second-order algorithms prove to be more robust to such ill conditioning. This is so since, using the curvature information, second-order methods properly rescale the gradient, such that it is a more appropriate direction to follow. For example, take the canonical second-order method, i.e., Newton’s method, which, in the unconstrained case, has updates of the form

where $\mathbf{g}(\mathbf{w}_{t})$ and $\mathbf{H}(\mathbf{w}_{t})$ denote the gradient and the Hessian of $F$ at $\mathbf{w}_{t}$, respectively. Classical results indicate that under certain assumptions, Newton’s method can achieve a locally super-linear convergence rate, which can be shown to be problem independent! Nevertheless, the cost of forming and inverting the Hessian is a major drawback in using Newton’s method in practice.

In this regard, there has been a long line of work that tries to provide sufficient second-order information with feasible computations. For example, among the class of quasi-Newton methods, the BFGS algorithm [NW06] and its limited memory version [LN89] are the most celebrated. However, the convergence guarantee of these methods can be much weaker than Newton’s methods. More recently, authors in [PW15, EM15, RKM16] considered using sketching and sampling techniques to construct an approximate Hessian matrix and using it in the update rule (2). They showed that such algorithms inherit a local linear-quadratic convergence rate with a substantial computational gain.

In this work, we propose novel, robust and highly efficient non-uniformly sub-sampled Newton methods (SSN) for a large sub-class of problem (1), where the Hessian of $F(\mathbf{w})$ in (1) can be written as

Second, when the dimension of the problem, i.e., $d$, is so large that solving the above linear system, i.e., (4), becomes infeasible, we consider solving (4) inexactly by using an iterative solver $\mathcal{A}$, e.g., Conjugate Gradient or Stochastic Gradient Descent, with a few iterations such that a high-quality approximate solution can be produced with a less complexity. Such inexact updates used in many second-order optimization algorithms have been well studied in [Byr+11, DES82].

Under certain conditions, it can be shown that this type of randomized Newton-type algorithm can achieve a local linear-quadratic convergence rate shown as below (formally stated in Lemma 7)

where $C_{q},C_{l}$ are some constants that can be controlled by the Hessian approximation quality, i.e., choice of sampling scheme $\mathcal{S}$,To be more precise, by a sampling scheme here, we mean the way we construct the sampling distribution $\{p_{i}\}_{i=1}^{n}$, e.g., uniform sampling distribution or leverage scores sampling distribution, and the value of sampling size $s$. and the solution quality of (4), i.e., choice of solver $\mathcal{A}$.By a solver here, we mean the choice of the specific algorithm with parameters specified, e.g., number of iterations, we use to obtain a high-quality approximation solution to the subproblem.

Different choices of sampling scheme $\mathcal{S}$ and solver $\mathcal{A}$ lead to different complexities in SSN. Below, we briefly discuss their effects. As can be seen, the total complexity of our algorithm can be characterized by the following three factors, each of which is affected by $\mathcal{S}$ or $\mathcal{A}$, or both.

Number of total iterations $T$ determined by the convergence rate, i.e., $C_{q}$ and $C_{l}$ in (5) which is affected by sampling scheme $\mathcal{S}$ and solver $\mathcal{A}$.

In each iteration, the time $t_{const}$ it needs to construct $\{p_{i}\}_{i=1}^{n}$ and sample $s$ terms, which is determined by sampling scheme $\mathcal{S}$.

With these, the total complexity can be expressed as

where $t_{grad}$ is the time it takes to compute the full gradient $\nabla F(\mathbf{w}_{t})$ which is not affected by the choice of $\mathcal{S}$ and $\mathcal{A}$ and will not be discussed in the rest of this paper.

As discussed above, the choice of sampling scheme $\mathcal{S}$ and solver $\mathcal{A}$ plays an important role in our algorithm. Below, we focus on $\mathcal{S}$ first and discuss a few concrete sampling schemes. A natural and simple approach is uniform sampling discussed in [EM15, RKM16, RKM16a]. The greatest advantage of uniform sampling is its simplicity of construction. However, in the presence of high non-uniformity among $\{\nabla^{2}f_{i}(\mathbf{w})\}_{i=1}^{n}$, the sampling size required to sufficiently capture the curvature information of the Hessian can be very large, which makes the resulting algorithm less beneficial.

In this paper, we consider two non-uniform sampling schemes, block norm squares and a new, and more general, notion of leverage scores named block partial leverage scores (Definition 6); see Section 3 for detailed construction. The motivation for using these schemes is that we can in fact view the sufficient conditions for achieving the local linear-quadratic convergence rate (5) as matrix approximation guarantees and these two sampling schemes can yield high-quality matrix approximations; see Section 4.1 for more details. Recall that, the choice of $\mathcal{S}$ affects the three terms, namely, $T$ (manifested in $C_{q}$ and $C_{l}$), $t_{const}$, $t_{solve}$, in (6). By leveraging and extending theories in randomized linear algebra, we can show how these terms may become when different sampling schemes are used as presented in Table 1. Detailed theory is elaborated in Section 4. As we can see, block norm squares sampling and leverage scores sampling require a smaller sampling size $s$ than uniform sampling does since $t_{solve}=sd^{2}$. Furthermore, the dependence of $C_{q}$ and $C_{l}$ on the condition number $\kappa$ reveals that leverage scores sampling is more robust to ill-conditioned problems; see Section 4.6.1 for more detailed discussions.

Next, we discuss the effect of the solver $\mathcal{A}$. Typically, a direct solver takes $\operatorname{\mathcal{O}}(sd^{2})$ time to solve the subproblem (4) where $s$ is the sampling size. This becomes prohibitive when $d$ is large. However, iterative solvers allow one to obtain a high quality approximation solution with a few iterations which may drastically drive down the complexity. For example, Conjugate Gradient (CG) takes $\operatorname{\mathcal{O}}(sd\sqrt{\kappa}\log(1/\epsilon_{0}))$ to return an approximate solution with relative error $\epsilon_{0}$ where $\kappa$ is the condition number of the problem. In Lemma 7 we show that such inexactness will not deteriorate the performance of SSN too much.

Indeed, based on (5), it is possible to choose the parameters, e.g., sampling size $s$ in the sampling scheme $\mathcal{S}$ and number of iterations to run in solver $\mathcal{A}$, so that SSN converges in a constant linear rate, e.g.,

By this way, the complexity per iteration of SSN can be explicitly given. In Table 2 we summarize these results with comparison to other stochastic second-order approaches such as LiSSA [ABH16]. It can be seen from Table 2 that compared to Newton’s methods, these stochastic second-order methods trade the coefficient of the leading term $\operatorname{\mathcal{O}}(nd)$ with some lower order terms that only depend on $d$ and condition numbers. Although SSN with non-uniform sampling has a quadratic dependence on $d$, its dependence on the condition number is better than the other methods. There are two main reasons. First, the total power of the condition number is lower, regardless of the versions of the condition number needed. Second, SSN (leverage scores) and SSN (block norm squares) only depend on $\kappa$ which can be significantly lower than the other two definitions of condition number, i.e., $\hat{\kappa}$ and $\bar{\kappa}$; see Section 4.6.2 for more details.

As we shall see (in Section 5), our algorithms converge much faster than other competing methods with ridge logistic regression. In particular, on several real datasets with a moderately high condition number and large $n$, our methods are at least twice as fast as Newton’s methods in finding a medium- or high-precision solution, while other methods including first-order methods converge slowly. Indeed, this phenomenon is well supported by our theoretical findings—the complexity of our algorithms has a lower dependence on the problem condition number and is immune to any non-uniformity among $\{\mathbf{A}_{i}(\mathbf{w})\}_{i=1}^{n}$ which may cause a factor of $n$ in the complexity (Table 2). In the following we present other prior work and details of our main contributions.

Recently, within the context of randomized second-order methods, many algorithms have been proposed that aim at reducing the computational costs involving pure Newton’s method. Among them, algorithms that employ uniform sub-sampling constitute a popular line of work [Byr+11, EM15, Mar10, VP11]. In particular [RKM16, RKM16a] consider a more general class of problems and, under a variety of conditions, thoroughly study the local and global convergence properties of sub-sampled Newton methods where the gradient and/or the Hessian are uniformly sub-sampled. Our work here, however, is more closely related to a recent work [PW15] (Newton Sketch), which considers a similar class of problems and proposes sketching the Hessian using random sub-Gaussian matrices or randomized orthonormal systems. Furthermore, [ABH16] proposes a stochastic algorithm (LiSSA) that, for solving the sub-problems, employs some unbiased estimators of the inverse of the Hessian.

The main technique used by [PW15, RKM16, RKM16a] and our work is sketching, which is a powerful technique in randomized linear algebra and many other applications [Woo14, Mah11, YMM16]. As mentioned above, Hessian approximation can be viewed as a matrix approximation problem. In terms of this, error analysis of matrix approximation based on leverage scores sampling has been well studied [DMM08, Dri+12, CMM15]. [HI15] show the lower bounds of sampling size for both block norm squares sampling and leverage scores sampling in approximating the Gram product matrix. Also [Coh+15] implies that uniform sampling cannot guarantee spectral approximation when the sampling size is only dependent on the lower dimension.

2 Main contributions

The contributions of this paper can be summarized as follows.

For the class of problems considered here, unlike the uniform sampling used in [Byr+11, EM15, RKM16, RKM16a], we employ two non-uniform sampling schemes based on block norm squares and a new, and more general, notion of leverage scores named block partial leverage scores (Definition 6). It can be shown that in the case of extreme non-uniformity among $\{\mathbf{A}_{i}(\mathbf{w})\}_{i=1}^{n}$, uniform sampling might require $\Omega(n)$ samples to capture the Hessian information appropriately. However, we show that our non-uniform sampling schemes result in sample sizes completely independent of $n$ and are immune to such non-uniformity.

Within the context of sub-sampled Newton-type algorithms, [Byr+11, RKM16a] incorporate inexact updates where the sub-problems are solved only approximately and study global convergence properties of their algorithms. We extend the study of inexactness to the finer level of local convergence analysis.

We provide a general structural result (Lemma 7) showing that, as in [EM15, PW15, RKM16], our main algorithm exhibits a linear-quadratic solution error recursion. However, we show that by using our non-uniform sampling strategies, the factors appearing in such error recursion enjoy a much better dependence on problem specific quantities, e.g., such as the condition number (Table 1′). For example, using block partial leverage score sampling, the factor for the linear term of the error recursion (14) is of order $\mathcal{O}(\sqrt{\kappa})$ as opposed to $\mathcal{O}(\kappa)$ for uniform sampling.

We numerically demonstrate the effectiveness and robustness of our algorithms in recovering the minimizer of ridge logistic regression on several real datasets with a moderately large condition number (Figures 1, 2 and 3). In particular, our algorithms are at least twice as fast as Newton’s methods in finding a medium- or high-precision solution, while other methods including first-order methods converge slowly.

The remainder of the paper is organized as follows. We begin in Section 2 with notation and assumptions that will be used throughout the paper. In Section 3, we describe our main algorithm and propose two non-uniform sampling schemes. Section 4 provides theoretical analysis. Finally, we present our numerical experiments in Section 5.

Background

Denote $\mathcal{K}$ be the tangent cone of constraints $\mathcal{C}$ at the optimum $\mathbf{w}^{\ast}$, i.e., $\mathcal{K}=\{\Delta|\mathbf{w}^{\ast}+t\Delta\in\mathcal{C}\text{ for some }t>0\}$.

Given a symmetric matrix $\mathbf{A}$ and a cone $\mathcal{K}$, we define the $\mathcal{K}$-restricted maximum and minimum eigenvalues as follows.

2 Assumptions

Throughout the paper, we use the following assumptions regarding the properties of the problem.

$F(\mathbf{w})$ is convex and twice differentiable. The Hessian is $L$-Lipschitz continuous, i.e.

$F(\mathbf{x})$ is locally strongly convex and smooth, i.e.,

Here we define the local condition number of the problem as $\kappa:=\nu/\mu$.

Main Algorithm: SSN with Non-uniform Sampling

Our proposed SSN method with non-uniform sampling is given in Algorithm 1. The core of our algorithm is based on choosing a sampling scheme $\mathcal{S}$ that, at every iteration, constructs a non-uniform sampling distribution $\{p_{i}\}_{i=1}^{n}$ over $\{\mathbf{A}_{i}(\mathbf{w}_{t})\}_{i=1}^{n}$ and then samples from $\{\mathbf{A}_{i}(\mathbf{w}_{t})\}_{i=1}^{n}$ to form the approximate Hessian, $\widetilde{\mathbf{H}}(\mathbf{w}_{t})$. The sampling sizes $s$ needed for different sampling distributions will be discussed in Sections 4.2 and 4.3. Since $\mathbf{H}(\mathbf{w})=\sum_{i=1}^{n}\mathbf{A}^{T}_{i}(\mathbf{w})\mathbf{A}_{i}(\mathbf{w})+\mathbf{Q}(\mathbf{w})$, the Hessian approximation essentially boils down to a matrix approximation problem. Here, we generalize the two popular non-uniform sampling strategies, i.e., leverage score sampling and block norm squares sampling, which are commonly used in the field of randomized linear algebra, particularly for matrix approximation problems [HI15, Mah11]. With an approximate Hessian constructed via non-uniform sampling, we may choose an appropriate solver $\mathcal{A}$ to the solve the sub-problem in Step 11 of Algorithm 1. Below we elaborate on the construction of the two non-uniform sampling schemes. Indeed, the sampling distribution is defined based on the matrix representation of $\{\mathbf{H}_{i}(\mathbf{w}_{t})\}_{i=1}^{n}$ — its augmented matrix defined as follows.

Define the augmented matrix of $\{\mathbf{H}_{i}(\mathbf{w}_{t})\}_{i=1}^{n}$ as

For the ease of presentation, throughout the rest of this section and next section, we use $\mathbf{A}$ and $\mathbf{Q}$ to denote $\mathbf{A}(\mathbf{w})$ and $\mathbf{Q}(\mathbf{w})$, respectively, as long as it is clear in the text.

The first option is to construct a sampling distribution based on the magnitude of $\mathbf{A}_{i}$. That is, define

This is an extension to the row norm squares sampling in which the intuition is to capture the importance of the blocks based on the “magnitudes” of the sub-Hessians.

The second option is to construct a sampling distribution based on leverage scores. Compared to the traditional matrix approximation problem, this problem has two major difficulties. First, here blocks are being sampled, not single rows. Second, the matrix being approximated involves not only $\mathbf{A}$ but also $\mathbf{Q}$.

To address the first difficulty, we follow the work by [CSHS11] in which a sparse sum of semi-definite matrices is found by sampling based on the trace of each semi-definite matrix after a proper transformation. By expressing $\mathbf{A}^{T}\mathbf{A}=\sum_{i=1}^{n}\mathbf{A}_{i}^{T}\mathbf{A}_{i}$, one can show that their approach is essentially sampling based on the sum of leverage scores that correspond to each block. For the second difficulty, inspired by the recently proposed ridge leverage scores [AM15, CMM15], we consider the leverage scores of a matrix that concatenates $\mathbf{A}$ and $\mathbf{Q}^{\frac{1}{2}}$. Combining these motivates us to define a new notion of leverage scores —- block partial leverage scores which is define as follows formally.

Then the sampling distribution is defined as

Remark. When each block of $\mathbf{A}$ has only one row and $\mathbf{Q}=\mathbf{0}$, the partial block leverage scores are equivalent to the ordinary leverage scores defined in Definition 4.

Theoretical Results

Below we study these contributing factors. Lemma 7 in Section 4.1 gives a general structural lemma that characterize the convergence rate which determines $T$. In Sections 4.2 and 4.3, Lemmas 8 and 10 discuss $t_{const}$ for the two sampling schemes respectively while Lemmas 9 and 11 give the required sampling size $s$ for the two sampling schemes respectively which directly affects the $t_{solve}$. Furthermore, $t_{solve}$ is also affected by the choice of solver which will be discussed in Section 4.4. Finally, the complexity results are summarized in Section 4.5 and a comparison with other methods is provided in Section 4.6.

Before diving into details of the complexity analysis, we state a structural lemma that characterizes the local convergence rate of our main algorithm, i.e., Algorithm 1. As discussed earlier, there are two layers of approximation in Algorithm 1, i.e., approximation of the Hessian by sub-sampling and inexactness of solving (8). For the first layer, we require the approximate Hessian to satisfy one of the following two conditions (in Sections 4.2 and 4.3 we shall see our construction of approximate Hessian via non-uniform sampling can achieve these conditions with a sampling size independent of $n$).

Note that (C1) and (C2) are two commonly seen guarantees for matrix approximation problems. In particular, (C2) is stronger in the sense that the spectrum of the approximated matrix $\mathbf{H}(\mathbf{w}_{t})$ is well preserved. Below in Lemma 7, we shall see such a stronger condition ensures a better dependence on the condition number in terms of the convergence rate. For the second layer of approximation, we require the solver to produce an $\epsilon_{0}$-approximate solution $\mathbf{w}_{t+1}$ satisfying

where $\mathbf{w}_{t+1}^{\ast}$ is the exact optimal solution to (8). Note that (12) implies an $\epsilon_{0}$-relative error approximation to the exact update direction, i.e., $\|\mathbf{v}-\mathbf{v}^{\ast}\|\leq\epsilon\|\mathbf{v}^{\ast}\|$ where $\mathbf{v}=\mathbf{w}_{t+1}-\mathbf{w}_{t},\>\mathbf{v}^{\ast}=\mathbf{w}_{t+1}^{\ast}-\mathbf{w}_{t}$.

Then requirement (12) is equivalent to finding an approximation solution $\mathbf{v}$ such that

Let $\{\mathbf{w}_{t}\}_{i=1}^{T}$ be the sequence generated based on update rule (8) with initial point $\mathbf{w}_{0}$ satisfying $\|\mathbf{w}_{0}-\mathbf{w}^{\ast}\|\leq\frac{\mu}{4L}$. Under Assumptions 1 and 2, if condition (C1) or (C2) is met, we have the following results.

If the subproblem is solved exactly, then the solution error satisfies the following recursion

where $C_{q}$ and $C_{l}$ are specified in (15) or (16) below.

If the subproblem is solved approximately and $\mathbf{w}_{t+1}$ satisfies (12), then the solution error satisfies the following recursion

where $C_{q}$ and $C_{l}$ are specified in (15) or (16) below.

Specifically, given any $\epsilon\in(0,1/2)$,

and due to Lemma 16 in Appendix A, (C2) is equivalent to

From this it is not hard to see that (C2) is strictly stronger than (C1). Also, in this case the Hessian approximation problem boils down to a matrix approximation problem. That is, given $\mathbf{A}$ and $\mathbf{Q}$, we want to construct a sampling matrix $\mathbf{S}$ efficiently such that the matrix $\mathbf{A}^{T}\mathbf{A}+\mathbf{Q}$ is well preserved. As we mentioned, leverage scores sampling and block norm squares sampling are two popular ways for this task. In the next two subsections we will focus on the theoretical properties of these two schemes.

2 Results for block partial leverage scores sampling

Since the block partial leverage scores are defined as the standard leverage scores of some matrix, we can make use of the fast approximation algorithm for standard leverage scores. Specifically, apply a variant of the algorithm in [Dri+12] by using the sparse subspace embedding [clarkson13sparse] as the underlying sketching method to further speed up the computation.

Given $\mathbf{w}$, under Assumption 3, with high probability, it takes $t_{const}=\operatorname{\mathcal{O}}(\operatorname{nnz}(\mathbf{A})\log n)$ time to construct a set of approximate leverage scores $\{\hat{\tau}^{\mathbf{Q}}_{i}(\mathbf{A})\}_{i=1}^{n}$ that satisfy $\tau^{\mathbf{Q}}_{i}(\mathbf{A})\leq\hat{\tau}^{\mathbf{Q}}_{i}(\mathbf{A})\leq\beta\cdot\tau^{\mathbf{Q}}_{i}(\mathbf{A})$ where $\{\tau_{i}\}_{i=1}^{n}$ are the block partial leverage scores of $\mathbf{H}(\mathbf{w})=\sum_{i=1}^{n}\mathbf{H}_{i}(\mathbf{w})+\mathbf{Q}(\mathbf{w})$ where $\mathbf{A}$ is the augmented matrix of $\{\mathbf{H}_{i}(\mathbf{w})\}_{i=1}^{n}$, and $\beta$ is a constant.

2.2 Sampling size

The following theorem indicates that if we sample the blocks of $\mathbf{A}$ based on block partial leverage scores with large enough sampling size, (18) holds with high probability.

Given $\mathbf{A}$ with $n$ blocks, $\mathbf{Q}\succeq\mathbf{0}$ and $\epsilon\in(0,1)$, let $\{\tau^{\mathbf{Q}}_{i}(\mathbf{A})\}_{i=1}^{n}$ be its block partial leverage scores and $\{\hat{\tau}^{\mathbf{Q}}_{i}(\mathbf{A})\}_{i=1}^{n}$ be their overestimates, i.e., $\hat{\tau}^{\mathbf{Q}}_{i}(\mathbf{A})\geq\tau^{\mathbf{Q}}_{i}(\mathbf{A}),i=1,...,n$. Let $p_{i}=\frac{\hat{\tau}^{\mathbf{Q}}_{i}(\mathbf{A})}{\sum_{j=1}^{n}\hat{\tau}^{\mathbf{Q}}_{j}(\mathbf{A})}$. Construct $\mathbf{S}\mathbf{A}$ by sampling the $i$-th block of $\mathbf{A}$ with probability $q_{i}=\min\{s\cdot p_{i},1\}$ and rescaling it by $1/\sqrt{q_{i}}$. Then if

with probability at least $1-\delta$, (18) holds, thus (C2) holds.

Remark. When $\{\tau^{\mathbf{Q}}_{i}(\mathbf{A})\}_{i=1}^{n}$ are the exact scores, since $\sum_{i=1}^{n}\tau_{i}^{\mathbf{Q}}(\mathbf{A})\leq\sum_{i=1}^{N+d}\tau_{i}(\bar{\mathbf{A}})=d$ where $\bar{\mathbf{A}}=\begin{pmatrix}\mathbf{A}\\ \mathbf{Q}^{\frac{1}{2}}\end{pmatrix}$, the above theorem indicates that less than $\operatorname{\mathcal{O}}(d\log d/\epsilon^{2})$ blocks are needed for (18) to hold.

3 Results for block norm squares sampling

To sample based on block norm squares, one has first compute the Frobenius norm of every block in the augmented matrix $\mathbf{A}$. This requires $\operatorname{\mathcal{O}}(\operatorname{nnz}(\mathbf{A}))$ time.

Given $\mathbf{w}$, under Assumption 3, it takes $t_{const}=\operatorname{\mathcal{O}}(\operatorname{nnz}(\mathbf{A}))$ time to construct a block norm squares sampling distribution for $\mathbf{H}(\mathbf{w})=\sum_{i=1}^{n}\mathbf{H}_{i}(\mathbf{w})+\mathbf{Q}(\mathbf{w})$ where $\mathbf{A}$ is the augmented matrix of $\{\mathbf{H}_{i}(\mathbf{w})\}_{i=1}^{n}$.

3.2 Sampling size

The following theorem [HI15] show the approximation error bound for Gram matrix. Here we extend it to our augmented matrix setting as follows,

Given $\mathbf{A}$ with $n$ blocks, $\mathbf{Q}\succeq\mathbf{0}$ and $\epsilon\in(0,1)$, for $i=1,\ldots,n$, let $r_{i}=\|\mathbf{A}_{i}\|_{F}^{2}$. Let $p_{i}=\frac{r_{i}}{\sum_{j=1}^{n}r_{j}}$. Construct $\mathbf{S}\mathbf{A}$ by sampling the $i$-th block of $\mathbf{A}$ with probability $q_{i}=\min\{s\cdot p_{i},1\}$ and rescaling it by $1/\sqrt{q_{i}}$. Then if

with probability at least $1-\delta$, (17) holds, thus (C1).

4 Discussion on the choice of solver

to denote the time it needs to solve the subproblem (8) using solver $\mathcal{A}$.

5 Complexities

Again, recall that in (11) the complexity of the sub-sampled Newton methods can be expressed as $T\cdot(t_{const}+t_{grad}+t_{solve})$. Combining the results from the previous few subsections, we have the following lemma characterizing the total complexity.

For Algorithm 1 with sampling scheme $\mathcal{S}$ and solver $\mathcal{A}$, the total complexity is

and the solution error is specified in Lemma 7. In the above, $t_{const}$ is specified in Theorem 8 and Theorem 10 and $s$ is specified in Theorem 9 and Theorem 11 depending on the choice of $\mathcal{S}$; $\mathcal{T}(\mathcal{A},\mathcal{C},s,d)$ is discussed in Section 4.4.

Indeed, Lemma 7 implies that the sub-sampled Newton method inherits a local constant linear convergence rate. This can be shown by choosing specific values for $\epsilon$ and $\epsilon_{0}$ in Lemma 7. The results are presented in the following corollary.

if block partial leverage scores sampling is used, the complexity per iteration in the local phase is

if block norm squares sampling is used, the complexity per iteration in the local phase is

6 Comparisons

As discussed above, the sampling scheme $\mathcal{S}$ plays a crucial role in sub-sampled Newton methods. Here, we compare the two proposed non-uniform sampling schemes, namely, block partial leverage scores sampling and block norm squares sampling, with uniform sampling. SSN with uniform sampling was discussed in [RKM16]. For completeness, we state the sampling size bound for uniform sampling. Note that, this upper bound for $s$ is tighter than the original analysis in [RKM16].

Given $\mathbf{A}$ with $n$ blocks, $\mathbf{Q}\succeq\mathbf{0}$ and $\epsilon\in(0,1)$, construct $\mathbf{S}\mathbf{A}$ by uniform sampling $s$ blocks from $\mathbf{A}$ and rescaling it by $\sqrt{n/s}$. Then if

with probability at least $1-\delta$, (17) holds, thus (C1) holds.

where constants $L,\mu,\nu$ are defined in Assumptions 1 and 2. Also, throughout this subsection, for randomized algorithms, we choose parameters such that the failure probability is a constant.

As can be seen in Table 1′, the greatest advantage of uniform sampling scheme comes from its simplicity of construction. On the other hand, as discussed in Sections 4.2.1 and 4.3.1, it takes nearly input sparsity time to construct the leverage scores sampling distribution or the block norm squares sampling distribution. When it comes to the sampling size $s$ for achieving (17) or (18), as suggested in (24), the one for uniform sampling can become $\Omega(n)$ when $\mathbf{A}$ is very non-uniform, i.e., $\max_{i}\|\mathbf{A}_{i}\|\approxeq\|\mathbf{A}\|$. It can be shown that for a given $\epsilon$, block norm squares sampling requires the smallest sampling size which leads to the smallest value of $t_{solve}$ in Table 1′.

It is worth pointing that, although either (17) or (18) is sufficient to yield a local linear-quadratic convergence rate, as (18) is essentially a stronger condition, it has better constants, i.e., $C_{q}$ and $C_{l}$. This fact is reflected in Table 1′. The constants $C_{q}$ and $C_{l}$ for leverage scores sampling have a better dependence on the local condition number $\kappa$ than the other two schemes since leverage scores sampling yields a sampling matrix that satisfies the spectral approximation guarantee (18). In fact, this difference can dramatically affect the performance of the algorithm when dealing with ill-conditioned problems. This is verified by numerical experiments; see Figure 1 in Section 5 for details.

Remark. Note that all the analysis including error recursion (Lemma 7) and the required sampling size $s$ for different sampling schemes are provided as upper bounds. There will be cases that the sampling size bound indicates a large value for $s$, in fact a much smaller sampling size $s$ suffices to yield good performance. For example, when the leverage scores are equal or close to a uniform distribution, the actual required sampling size for uniform sampling scheme is much less than (24).

6.2 Comparison between various methods

Next, we compare our main algorithm with other stochastic second-order methods including [RKM16, ABH16]. Since these essentially imply a constant linear convergence rate, i.e.,

An immediate relationship between the three condition numbers is $\kappa(\mathbf{w})\leq\hat{\kappa}(\mathbf{w})\leq\bar{\kappa}(\mathbf{w}).$ The connections between these condition numbers depend on the properties of $\mathbf{H}_{i}(\mathbf{w})$. Roughly speaking, when all $\mathbf{H}_{i}(\mathbf{w})$’s are “close” to each other, then $\lambda_{\max}^{\mathcal{K}}(\sum_{i=1}^{n}\mathbf{H}_{i}(\mathbf{w}))\approx\sum_{i=1}^{n}\lambda_{\max}^{\mathcal{K}}(\mathbf{H}_{i}(\mathbf{w}))\approx n\cdot\max_{i}\lambda_{\max}^{\mathcal{K}}(\mathbf{H}_{i}(\mathbf{w}))$, and thus $\kappa\approx\hat{\kappa}$. And similarly, $\kappa\approx\bar{\kappa}$. While in many cases, some $\mathbf{H}_{i}(\mathbf{w})$’s can be very different from the rest. For example, when solving linear regression, the Hessian $\mathbf{H}(\mathbf{w})=\mathbf{A}^{T}\mathbf{A}$, where $\mathbf{A}$ is the data matrix with each row as a data point. When the rows are not very uniform, it can be the case that $\kappa$ is smaller than $\hat{\kappa}$ and $\bar{\kappa}$ by a a factor of $n$.

Given the notation we defined, we summarize the complexities of different algorithms in Table 2′ (identical to Table 2 in Section 1) including Newton’s methods with CG solving the subproblem. One immediate conclusion we can draw is that compared to Newton’s methods, these stochastic second-order methods trade the coefficient of the leading term $\operatorname{\mathcal{O}}(nd)$ with some lower order terms that only depend on $d$ and condition numbers (assuming $\operatorname{nnz}(\mathbf{A})\approx nd$). Therefore, one should expect these algorithm to perform well when $n\gg d$ and the problem is fairly well-conditioned.

Although SSN with non-uniform sampling has a quadratic dependence on $d$, its dependence on the condition number is better than the other methods. There are two main reasons. First, the total power of the condition number is lower, regardless the versions of the condition number needed. Second, SSN (leverage scores) and SSN (block norm squares) only depend on $\kappa$ which can be significantly lower than the other two definitions of condition number according to the discussion above. Overall, SSN (leverage scores) is more robust to the ill-conditioned problems.

Numerical Experiments

For our numerical simulations, we consider a very popular instance of GLMs, namely, logistic regression, where $\psi(u,y)=\log(1+\exp(-uy))$ and $\mathcal{Y}=\{\pm 1\}$. Table 4 summarizes the datasets used in our experiments.

We compare the performance of the following five algorithms: (i) Newton: the standard Newton’s method, (ii) Uniform: SSN with uniform sampling, (iii) PLevSS: SSN with partial leverage scores sampling, (iv) RNormSS: SSN with block (row) norm squares sampling, and (v) LBFGS-k: standard L-BFGS method [LN89] with history size $k$, (vi) GD: Gradient Descent, (vii) AGD: Accelerated Gradient Descent(AGD) [Nes13]. Note that, despite all of our effort, we could not compare with methods introduced in [ABH16, EM15] as they seem to diverge in our experiments.

All algorithms are initialized with a zero vector.Theoretically, the suitable initial point for all the algorithms is the one with which the standard Newton’s method converges with a unit stepsize. Here, $\mathbf{w}_{0}=\mathbf{0}$ happens to be one such good starting point. We also use CG to solve the sub-problem approximately to within $10^{-6}$ relative residue error. In order to compute the relative error $\|\mathbf{w}_{t}-\mathbf{w}^{*}\|/\|\mathbf{w}^{*}\|$, an estimate of $\mathbf{w}^{\ast}$ is obtained by running the standard Newton’s method for sufficiently long time. Note here, in SSN with partial leverage score sampling, we recompute the leverage scores every $10$ iterations. Roughly speaking, these “stale” leverage scores can be viewed as approximate leverage scores for the current iteration with approximation quality that can be upper bounded by the change of the Hessian and such quantity is often small in practice. So reusing the leverage scores allows us to further drive down the running time.

We first investigate the effect of the condition number, controlled by varying $\lambda$, on the performance of different methods, and the results are depicted in Figure 1. It can be seen that in well-conditioned cases, all sampling schemes work equally well. However, as the condition number gets larger, the performance of uniform sampling deteriorates, while non-uniform sampling, in particular leverage score sampling, shows a great degree of robustness to such ill-conditioning effect. The experiments shown in Figure 1 are consistent with the theoretical results of Table 1′.

Next, we compare the performance of various methods as measured by relative-error of the solution vs. running time. First, we provide a set of empirical comparison between first-order and second-order methods in Figure 2.For each sub-sampled Newton method, the sampling size is determined by choosing the best value from $\{10d,20d,30d,...,100d,200d,300d,...,1000d\}$ in the sense that the objective value drops to $1/3$ of initial function value first. This is on dataset CT Slice with two different $\lambda$’s. As can be seen clearly in Figure 2, SSN with non-uniform sampling not only drives down the loss function $F(\mathbf{w})$ to an arbitrary precision much more quickly, but also recovers the minimizer $\mathbf{w}^{\ast}$ to a high precision while first-order methods such as Gradient Descent converge very slowly. More importantly, unlike SSN with uniform sampling and LBFGS, non-uniform SSN exhibits a better robustness to condition number as it performance doesn’t deteriorate much when the problem becomes more ill-conditioned (by setting the regularization $\lambda$ smaller in Figures 2(c) and Figure 2(d)). This robustness to condition number allows our approach to excel for a wider range of models.

A more comprehensive comparison among various second-order methods on the four datasets is presented in Figure 3. It can be seen that, in most cases, SSN with non-uniform sampling schemes, i.e., PLevSS and RNormSS, outperform the other algorithms, especially Newton’s method. In particular, they can be as twice faster as Newton’s method. This is because on datasets with large $n$, the computational gain of our sub-sampled Newton methods in forming the (approximate) Hessian is significant while their convergence rate is only slightly worse than Newton’s method (not shown here). Moreover, recall that in Section 4.6 we discussed that the convergence rate of SSN with uniform sampling relies on $\bar{\kappa}$. When the problem exhibits a high non-uniformity among data points, i.e., $\bar{\kappa}$ is much higher than $\kappa$ as shown in Table 4, uniform sampling scheme performs poorly, e.g., in Figure 3(b).

Conclusions

In this paper, we propose non-uniformly sub-sampled Newton methods with inexact update for a class of constrained problems. We show that our algorithms have a better dependence on the condition number and enjoy a lower per-iteration complexity, compared to other similar existing methods. Theoretical advantages are numerically demonstrated.

Acknowledgments. We would like to acknowledge the Army Research Office and the Defense Advanced Research Projects Agency for providing partial support for this work.

References

Appendix A Results of Block Partial Leverage Scores

In this work, we propose to use a new notion of leverage scores, namely, block partial leverage scores, for approximating matrix of the form $\mathbf{A}^{T}\mathbf{A}+\mathbf{Q}$. In this section, we give its theoretical guarantee, i.e., quality of approximation, which will be used in the proofs later.

Denote $\bar{\mathbf{A}}=\begin{pmatrix}\mathbf{A}\\ \mathbf{Q}^{\frac{1}{2}}\end{pmatrix}$. Let $\bar{\mathbf{A}}=\bar{\mathbf{U}}\mathbf{R}$ where $\bar{\mathbf{U}}$ has orthonormal columns. Then define $\mathbf{U}=\mathbf{A}\mathbf{R}^{-1}$ and $\mathbf{U}_{i}=\mathbf{A}_{i}\mathbf{R}^{-1}$ for $i=1,\ldots,n$. By definition, the true partial leverage scores $\tau_{i}^{\mathbf{Q}}(\mathbf{A})$’s are defined as $\tau_{i}=\textbf{tr}(\mathbf{U}_{i}\mathbf{U}_{i}^{T})$. For simplicity, we use $\tau_{i}$ to denote $\tau_{i}^{\mathbf{Q}}(\mathbf{A})$.

In the following we bound $\|\mathbf{U}^{T}\mathbf{S}^{T}\mathbf{S}\mathbf{U}-\mathbf{U}^{T}\mathbf{U}\|\leq\epsilon$ with high probability. For $i=1,\ldots,n$, define

Also define $\mathbf{Y}=\sum_{i=1}^{n}\mathbf{X}_{i}$. We have $\mbox{}{\mathbf{E}}\left[\mathbf{X}_{i}\right]=0$. In the following we bound $\|\mathbf{Y}\|$ using matrix Bernstein bound.

Next, we bound $\mbox{}{\mathbf{E}}\left[\mathbf{Y}^{2}\right]=\sum_{i=1}^{n}\mbox{}{\mathbf{E}}\left[\mathbf{X}_{i}^{2}\right]$. We have

Since $\mathbf{U}$ consists of a subset of rows of $\bar{\mathbf{U}}$, one can show that $\mathbf{U}^{T}\mathbf{U}\preceq\bar{\mathbf{U}}^{T}\bar{\mathbf{U}}=\mathbf{I}$.

Given these, by the matrix Bernstein bound[Tro15], we have when

Remark. In (28), since each element of $\mathbf{D}$ is no greater than $1$ and $\mathbf{U}^{T}\mathbf{U}\preceq\mathbf{I}$, we have

Appendix B Proofs

In the section, we will prove the results in Lemma 7. Specifically we have two parts of proof. Part 1 is to prove the case when the subproblem is solved exactly and Part 2 is for the case when the subproblem is solved approximately. And in Part 1, we also have two cases, one is recursion inequality under condition (C1) and the other is the recursion inequality under condition (C2).

Since $\|\Delta_{t}\|_{2}\leq\frac{\mu}{4L}$, then

In the part, we consider the subproblem is solved exactly every iteration in Algorithm 1. Specifically we show the following two results:

Under condition (C1), the error recursion (13) holds with factors in (15).

Under condition (C2), the error recursion (13) holds with factors in (16).

If the subproblem $\min_{\mathbf{w}\in\mathcal{C}}\Psi_{t}(\mathbf{w})$ is solved exactly in Algorithm 1, namely

Then $\Psi_{t}(\mathbf{w}_{t+1})\leq\Psi_{t}(\mathbf{w}^{*})$. By expanding both sides, we have

The third term on the right hand side $\nabla F(\mathbf{w}^{*})^{T}\Delta_{t+1}\geq 0$ because of the optimality condition.

Now we move on to prove the case when the subproblem is approximately solved. Specifically we want to show under condition (12), the error recursion (14) holds where $C_{l},C_{q}$ are the same constants in the case when the problem is solved exactly.

Consider at the iteration $t$ in Algorthm 1. First, $\mathbf{w}_{t+1}^{*}=\operatorname{argmin}_{\mathbf{w}\in\mathcal{C}}\Psi_{t}(\mathbf{w})$ (Note that here $\mathbf{w}_{t+1}$ is not the minimizer any more). Then based on the proof in Part 1, we have

B.2 Proof of Theorem 8

B.3 Proof of Theorem 9

Based on Lemma 16 (stated below), we convert condition (C2) to a standard matrix product approximation guarantee. A direct corollary of Theorem 15 completes the proof.

First, it is straightforward to see (b) $\Rightarrow$ (a) by setting $\mathbf{x}=\mathbf{y}$ in (82). So now we prove the other direction.

Denote $\bar{\mathbf{A}}=\begin{pmatrix}\mathbf{A}\\ \mathbf{Q}^{\frac{1}{2}}\end{pmatrix}$. Let $\bar{\mathbf{A}}=\bar{\mathbf{U}}\mathbf{R}$ where $\bar{\mathbf{U}}$ has orthonormal columns. Then define $\mathbf{U}=\mathbf{A}\mathbf{R}^{-1}$ and $\mathbf{U}_{i}=\mathbf{A}_{i}\mathbf{R}^{-1}$ for $i=1,\ldots,n$. Then (81) is equivalent to

Since $\mathbf{R}^{-1}$ is a full rank matrix, then

B.4 Proof of Corollary 13

According to Lemma 7, we have the follow error recursion

If leverage scores sampling is used, then

Now choose $\epsilon\leq\min\left\{\frac{1}{10\sqrt{\kappa}},0.1\right\}$ and $\epsilon_{0}\leq 0.01$, then we can get $\rho<0.9$.

Therefore, the complexity per iteration is

Similarly, if block norm squares sampling is used, then

Now choose $\epsilon\leq\min\left\{\frac{1}{10\kappa},0.1\right\}$ and $\epsilon_{0}\leq 0.01$, then we can get $\rho<0.9$. Similar to the case using leverage scores sampling, we get the total complexity per iteration which is

B.5 Proof of Theorem 14

And $\mathbf{Y}=\sum_{j=1}^{s}\mathbf{X}_{j}(=\mathbf{A}^{T}\mathbf{S}^{T}\mathbf{S}\mathbf{A}-\mathbf{A}^{T}\mathbf{A}$). Then $\mbox{}{\mathbf{E}}\left[\mathbf{X}_{j}\right]=0$. In the following we will bound $\|\mathbf{Y}\|$ through matrix Bernstein inequality. For convenience, let’s denote $K_{t}:=\max_{i}\|\mathbf{A}_{i}\|^{2}$.

By the matrix Bernstein bound [Tro15], we have when

Now choose scale the $\epsilon_{0}=\|A\|^{2}\epsilon$, where $\epsilon\in(0,1)$, then

with probability at least $1-\delta$, $\|\mathbf{A}^{T}\mathbf{S}^{T}\mathbf{S}\mathbf{A}-\mathbf{A}^{T}\mathbf{A}\|\leq\epsilon\cdot\|\mathbf{A}^{T}\mathbf{A}\|$ holds. Since $\mathbf{Q}\succeq\mathbf{0}$, then $\|\mathbf{A}^{T}\mathbf{A}\|\leq\|\mathbf{A}^{T}\mathbf{A}+\mathbf{Q}\|$. Therefore condition (C1) holds. And this completes the proof.