Second-Order Stochastic Optimization for Machine Learning in Linear Time

Introduction

In recent literature stochastic first-order optimization has taken the stage as the primary workhorse for training learning models, due in large part to its affordable computational costs which are linear (in the data representation) per iteration. The main research effort devoted to improving the convergence rates of first-order methods have introduced elegant ideas and algorithms in recent years, including adaptive regularization [DHS11], variance reduction [JZ13, DBLJ14], dual coordinate ascent [SSZ13], and many more.

In contrast, second-order methods have typically been much less explored in large scale machine learning (ML) applications due to their prohibitive computational cost per iteration which requires computation of the Hessian in addition to a matrix inversion. These operations are infeasible for large scale problems in high dimensions.

In this paper we propose a family of novel second-order algorithms, LiSSA (Linear time Stochastic Second-Order Algorithm) for convex optimization that attain fast convergence rates while also allowing for an implementation with linear time per-iteration cost, matching the running time of the best known gradient-based methods. Moreover, in the setting where the number of training examples $m$ is much larger than the underlying dimension $d$, we show that our algorithm has provably faster running time than the best known gradient-based methods.

Formally, the main optimization problem we are concerned with is the empirical risk minimization (ERM) problem:

Our focus is second-order optimization methods (Newton’s method), where in each iteration, the underlying principle is to move to the minimizer of the second-order Taylor approximation at any point. Throughout the paper, we will let $\nabla^{-2}f(\mathbf{x})\triangleq\left[\nabla^{2}f(\mathbf{x})\right]^{-1}$. The update of Newton’s method at a point $\mathbf{x}_{t}$ is then given by

Certain desirable properties of Newton’s method include the fact that its updates are independent of the choice of coordinate system and that the Hessian provides the necessary regularization based on the curvature at the present point. Indeed, Newton’s method can be shown to eventually achieve quadratic convergence [Nes13]. Although Newton’s method comes with good theoretical guarantees, the complexity per step grows roughly as $\Omega(md^{2}+d^{\omega})$ (the former term for computing the Hessian and the latter for inversion, where $\omega\approx 2.37$ is the matrix multiplication constant), making it prohibitive in practice. Our main contribution is a suite of algorithms, each of which performs an approximate Newton update based on stochastic Hessian information and is implementable in linear $O(d)$ time. These algorithms match and improve over the performance of first-order methods in theory and give promising results as an optimization method on real world data sets. In the following we give a summary of our results. We propose two algorithms, LiSSA and LiSSA-Sample.

LiSSA: Algorithm 1 is a practical stochastic second-order algorithm based on a novel estimator of the Hessian inverse, leading to an efficient approximate Newton step (Equation 1). The estimator is based on the well known Taylor approximation of the inverse (Fact 2) and is described formally in Section 3.1. We prove the following informal theorem about LiSSA.

LiSSA returns a point $\mathbf{x}_{t}$ such that $f(\mathbf{x}_{t})\leq\min_{\mathbf{x}^{*}}f(\mathbf{x}^{*})+\varepsilon$ in total time

where $\kappa$ is the underlying condition number of the problem and $S_{1}$ is a bound on the variance of the estimator.

The precise version of the above theorem appears as Theorem 3.3. In theory, the best bound we can show for $S_{1}$ is $O(\kappa^{2})$; however, in our experiments we observe that setting $S_{1}$ to be a small constant (often 1) is sufficient. We conjecture that $S_{1}$ can be improved to $O(1)$ and leave this for future work. If indeed $S_{1}$ can be improved to $O(1)$ (as is indicated by our experiments), LiSSA enjoys a convergence rate comparable to first-order methods. We provide a detailed comparison of our results with existing first-order and second-order methods in Section 1.2. Moreover, in Section 7 we present experiments on real world data sets that demonstrate that LiSSA as an optimization method performs well as compared to popular first-order methods. We also show that LiSSA runs in time proportional to input sparsity, making it an attractive method for high-dimensional sparse data.

LiSSA-Sample: This variant brings together efficient first-order algorithms with matrix sampling techniques [LMP13, CLM+15] to achieve better runtime guarantees than the state-of-the-art in convex optimization for machine learning in the regime when $m>d$. Specifically, we prove the following theorem:

LiSSA-Sample returns a point $\mathbf{x}_{t}$ such that $f(\mathbf{x}_{t})\leq\min_{\mathbf{x}^{*}}f(\mathbf{x}^{*})+\varepsilon$ in total time

The above result improves upon the best known running time for first-order methods achieved by acceleration when we are in the setting where $\kappa>m>>d$. We discuss the implication of our bounds and further work in Section 1.3.

We also remark that all of our results focus on the very high accuracy regime. In general the benefits of linear convergence and second-order methods can be seen to be effective only when considerably small error is required. This is also the case for recent advances in fast first-order methods where their improvement over stochastic gradient descent becomes apparent only in the high accuracy regime. Our experiments also demonstrate that second-order methods can improve upon fast first-order methods in the regime of very high accuracy. While it is possible that this regime is less interesting for generalization, in this paper we focus on the optimization problem itself.

We further consider the special case when the function $f$ is self-concordant. Self-concordant functions are a sub-class of convex functions which have been extensively studied in convex optimization literature in the context of interior point methods [Nem04]. For self-concordant functions we propose an algorithm (Algorithm 5) which achieves linear convergence with running time guarantees independent of the condition number. We prove the formal running time guarantee as Theorem 6.2.

We believe our main contribution to be a demonstration of the fact that second-order methods are comparable to, or even better than, first-order methods in the large data regime, in both theory and practice.

LiSSA: The key idea underlying LiSSA is the use of the Taylor expansion to construct a natural estimator of the Hessian inverse. Indeed, as can be seen from the description of the estimator in Section 3.1, the estimator we construct becomes unbiased in the limit as we include additional terms in the series. We note that this is not the case with estimators that were considered in previous works such as that of [EM15], and so we therefore consider our estimator to be more natural. In the implementation of the algorithm we achieve the optimal bias/variance trade-off by truncating the series appropriately.

An important observation underlying our linear time $O(d)$ step is that for GLM functions, $\nabla^{2}f_{i}(\mathbf{x})$ has the form $\alpha\mathbf{v}_{i}\mathbf{v}_{i}^{\top}$ where $\alpha$ is a scalar dependent on $\mathbf{v}_{i}^{\top}\mathbf{x}$. A single step of LiSSA requires us to efficiently compute $\nabla^{2}f_{i}(\mathbf{x})\mathbf{b}$ for a given vector $\mathbf{b}$. In this case it can be seen that the matrix-vector product reduces to a vector-vector product, giving us an $O(d)$ time update.

LiSSA-Sample: LiSSA-Sample is based on Algorithm 2, which represents a general family of algorithms that couples the quadratic minimization view of Newton’s method with any efficient first-order method. In essence, Newton’s method allows us to reduce (up to $\log\log$ factors) the optimization of a general convex function to solving intermediate quadratic or ridge regression problems. Such a reduction is useful in two ways.

First, as we demonstrate through our algorithm LiSSA-Sample, the quadratic nature of ridge regression problems allows us to leverage powerful sampling techniques, leading to an improvement over the running time of the best known accelerated first-order method. On a high level this improvement comes from the fact that when solving a system of $m$ linear equations in $d$ dimensions, a constant number of passes through the data is enough to reduce the system to $O(d\log(d))$ equations. We carefully couple this principle and the computation required with accelerated first-order methods to achieve the running times for LiSSA-Sample. The result for the quadratic sub-problem (ridge regression) is stated in Theorem 5.1, and the result for convex optimization is stated in Theorem 5.2.

The second advantage of the reduction to quadratic sub-problems comes from the observation that the intermediate quadratic sub-problems can potentially be better conditioned than the function itself, allowing us a better choice of the step size in practice. We define these local notions of condition number formally in Section 2.1 and summarize the typical benefits for such algorithms in Theorem 4.1. In theory this is not a significant improvement; however, in practice we believe that this could be significant and lead to runtime improvements.While this is a difficult property to verify experimentally, we conjecture that this is a possible explanation for why LiSSA performs better than first-order methods on certain data sets and ranges of parameters.

To achieve the bound for LiSSA-Sample we extend the definition and procedure for sampling via leverage scores described by [CLM+15] to the case when the matrix is given as a sum of PSD matrices and not just rank one matrices. We reformulate and reprove the theorems proved by [CLM+15] in this context, which may be of independent interest.

2 Comparison with Related Work

In this section we aim to provide a short summary of the key ideas and results underlying optimization methods for large scale machine learning. We divide the summary into three high level principles: first-order gradient-based methods, second-order Hessian-based methods, and quasi-Newton methods. For the sake of brevity we will restrict our summary to results in the case when the objective is strongly convex, which as justified above is usually ensured by the addition of an appropriate regularizer. In such settings the main focus is often to obtain algorithms which have provably linear convergence and fast implementations.

First-Order Methods: First-order methods have dominated the space of optimization algorithms for machine learning owing largely to the fact that they can be implemented in time proportional to the underlying dimension (or sparsity). Gradient descent is known to converge linearly to the optimum with a rate of convergence that is dependent upon the condition number of the objective. In the large data regime, stochastic first-order methods, introduced and analyzed first by [RM51], have proven especially successful. Stochastic gradient descent (SGD), however, converges sub-linearly even in the strongly convex setting. A significant advancement in terms of the running time of first-order methods was achieved recently by a clever merging of stochastic gradient descent with its full version to provide variance reduction. The representative algorithms in this space are SAGA [RSB12, DBLJ14] and SVRG [JZ13, ZMJ13]. The key technical achievement of the above algorithms is to relax the running time dependence on $m$ (the number of training examples) and $\kappa$ (the condition number) from a product to a sum. Another algorithm which achieves similar running time guarantees is based on dual coordinate ascent, known as SDCA [SSZ13].

Further improvements over SAGA, SVRG and SDCA have been obtained by applying the classical idea of acceleration emerging from the seminal work of [Nes83]. The progression of work here includes an accelerated version of SDCA [SSZ16]; APCG [LLX14]; Catalyst [LMH15], which provides a generic framework to accelerate first order algorithms; and Katyusha [AZ16], which introduces the concept of negative momentum to extend acceleration for variance reduced algorithms beyond the strongly convex setting. The key technical achievement of accelerated methods in general is to reduce the dependence on condition number from linear to a square root. We summarize these results in Table 1.

LiSSA places itself naturally into the space of fast first-order methods by having a running time dependence that is comparable to SAGA/SVRG (ref. Table 1). In LiSSA-Sample we leverage the quadratic structure of the sub-problem for which efficient sampling techniques have been developed in the literature and use accelerated first-order methods to improve the running time in the case when the underlying dimension is much smaller than the number of training examples. Indeed, to the best of our knowledge LiSSA-Sample is the theoretically fastest known algorithm under the condition $m>>d$. Such an improvement seems out of hand for the present first-order methods as it seems to strongly leverage the quadratic nature of the sub-problem to reduce its size. We summarize these results in Table 1.

Second-Order Methods: Second-order methods such as Newton’s method have classically been used in optimization in many different settings including development of interior point methods [Nem04] for general convex programming. The key advantage of Newton’s method is that it achieves a linear-quadratic convergence rate. However, naive implementations of Newton’s method have two significant issues, namely that the standard analysis requires the full Hessian calculation which costs $O(md^{2})$, an expense not suitable for machine learning applications, and the matrix inversion typically requires $O(d^{3})$ time. These issues were addressed recently by the algorithm NewSamp [EM15] which tackles the first issue by subsampling and the second issue by low-rank projections. We improve upon the work of [EM15] by defining a more natural estimator for the Hessian inverse and by demonstrating that the estimator can be computed in time proportional to $O(d)$. We also point the reader to the works of [Mar10, BCNN11] which incorporate the idea of taking samples of the Hessian; however, these works do not provide precise running time guarantees on their proposed algorithm based on problem specific parameters. Second-order methods have also enjoyed success in the distributed setting [SSZ14].

Quasi-Newton Methods: The expensive computation of the Newton step has also been tackled via estimation of the curvature from the change in gradients. These algorithms are generally known as quasi-Newton methods stemming from the seminal BFGS algorithm [Bro70, Fle70, Gol70, Sha70]. The book of [NW06] is an excellent reference for the algorithm and its limited memory variant (L-BFGS). The more recent work in this area has focused on stochastic quasi-Newton methods which were proposed and analyzed in various settings by [SYG07, MR14, BHNS16]. These works typically achieve sub-linear convergence to the optimum. A significant advancement in this line of work was provided by [MNJ16] who propose an algorithm based on L-BFGS by incorporating ideas from variance reduction to achieve linear convergence to the optimum in the strongly convex setting. Although the algorithm achieves linear convergence, the running time of the algorithm depends poorly on the condition number (as acknowledged by the authors). Indeed, in applications that interest us, the condition number is not necessarily a constant as is typically assumed to be the case for the theoretical results in [MNJ16].

Our key observation of linear time Hessian-vector product computations for machine learning applications provides evidence that in such instances, obtaining true Hessian information is efficient enough to alleviate the need for quasi-Newton information via gradients.

3 Discussion and Subsequent Work

In this section we provide a brief survey of certain technical aspects of our bounds which have since been improved by subsequent work.

An immediate improvement in terms of $S_{2}\sim\kappa$ (in fact suggested in the original manuscript) was achieved by [BBN16] via conjugate gradient on a sub-sampled Hessian which reduces this to $\sqrt{\kappa}$. A similar improvement can also be achieved in theory through the extensions of LiSSA proposed in the paper. As we show in Section 7, the worse dependence on condition number has an effect on the running time when $\kappa$ is quite large.Equivalently, $\lambda$ is small. Accelerated first-order methods, such as APCG [LLX14], outperform LiSSA in this regime. To the best of our knowledge second-order stochastic methods have so far not exhibited an improvement in that regime experimentally. We believe a more practical version of LiSSA-Sample could lead to improvements in this regime, leaving this as future work.

To the best of our knowledge the factor of $S_{1}=\kappa^{2}$ that appears to reduce the variance of our estimator has yet not been improved despite it being $O(1)$ in our experiments. This is an interesting question to which partial answers have been provided in the analysis of [YLZ17].

Significant progress has been made in the space of inexact Newton methods based on matrix sketching techniques. We refer the reader to the works of [PW15, XYRK+16, Coh16, LACBL16, YLZ17] and the references therein.

We would also like to comment on the presence of a warm start parameter $\frac{1}{\kappa M}$ in our proofs of Theorems 3.3 and 5.1. In our experiments the warm start we required would be quite small (often a few steps of gradient descent would be sufficient) to make LiSSA converge. The warm start does not affect the asymptotic results proven in Theorems 3.3 and 5.1 because getting to such a warm start is independent of $\varepsilon$. However, improving this warm start, especially in the context of Theorem 5.1, is left as interesting future work.

On the complexity side, [AS16] proved lower bounds on the best running times achievable by second-order methods. In particular, they show that to get the faster rates achieved by LiSSA-Sample, it is necessary to use a non-uniform sampling based method as employed by LiSSA-Sample. We would like to remark that in theory, up to logarithmic factors, the running time of LiSSA-Sample is still the best achieved so far in the setting $m>>d$. Some of the techniques and motivations from this work were also generalized by the authors to provide faster rates for a large family of non-convex optimization problems [AAZB+17].

4 Organization of the Paper

The paper is organized as follows: we first present the necessary definitions, notations and conventions adopted throughout the paper in Section 2. We then describe our estimator for LiSSA, as well as state and prove the convergence guarantee for LiSSA in Section 3. After presenting a generic procedure to couple first-order methods with Newton’s method in 4, we present LiSSA-Sample and the associated fast quadratic solver in Section 5. We then present our results regarding self-concordant functions in Section 6. Finally, we present an experimental evaluation of LiSSA in Section 7.

Preliminaries

The following is a well known fact about the inverse of a matrix $A$ s.t. $\|A\|\leq 1$ and $A\succ 0$:

For an $\alpha$-strongly convex and $\beta$-smooth function $f$, the condition number of the function is defined as $\kappa(f)\triangleq\frac{\beta}{\alpha}$, or $\kappa$ when the function is clear from the context. Note that by definition this corresponds to the following notion:

We define a slightly relaxed notion of condition number where the $max$ moves out of the fraction above. We refer to this notion as a local condition number $\kappa_{l}$ as compared to the global condition number $\kappa$ defined above:

It follows that $\kappa_{l}\leq\kappa$. The above notions are defined for any general function $f$, but in the case of functions of the form $f(\mathbf{x})=\frac{1}{m}\sum\limits_{k=1}^{m}f_{k}(\mathbf{x})$, a further distinction is made with respect to the component functions. We refer to such definitions of the condition number by $\hat{\kappa}$. In such cases one typically assumes the each component is bounded by $\beta_{max}(\mathbf{x})\triangleq\max\limits_{k}\lambda_{\max}(\nabla^{2}f_{k}(\mathbf{x}))$. The running times of algorithms like SVRG depend on the following notion of condition number:

Similarly, we define a notion of local condition number for $\hat{\kappa}$, namely

and it again follows that $\hat{\kappa_{l}}\leq\hat{\kappa}$.

For our (admittedly pessimistic) bounds on the variance we also need a per-component strong convexity bound $\alpha_{\min}(\mathbf{x})\triangleq\min\limits_{k}\lambda_{\min}(\nabla^{2}f_{k}(\mathbf{x}))$. We can now define

Assumptions: In light of the previous definitions, we make the following assumptions about the given function $f(\mathbf{x})=\frac{1}{m}\sum\limits_{k=1}^{m}f_{k}(\mathbf{x})$ to make the analysis easier. We first assume that the regularization term has been divided equally and included in $f_{k}(\mathbf{x})$. We further assume that each $\nabla^{2}f_{k}(\mathbf{x})\preceq I$.The scaling is without loss of generality, even when looking at additive errors, as this gets picked up in the log-term due to the linear convergence. We also assume that $f$ is $\alpha$-strongly convex and $\beta$-smooth, $\hat{\kappa}_{l}$ is the associated local condition number and $\nabla^{2}f$ has a Lipschitz constant bounded by $M$.

We now collect key concepts and pre-existing results that we use for our analysis in the rest of the paper.

Matrix Concentration: The following lemma is a standard concentration of measure result for sums of independent matrices.The theorem in the reference states the inequality for the more nuanced bounded variance case. We only state the simpler bounded spectral norm case which suffices for our purposes. An excellent reference for this material is by [Tro12].

Consider a finite sequence $\{X_{k}\}$ of independent, random, Hermitian matrices with dimension $d$. Assume that

Define $Y=\sum_{k}X_{k}$. Then we have for all $t\geq 0$,

Accelerated SVRG: The following theorem was proved by [LMH15].

Given a function $f(\mathbf{x})=\frac{1}{m}\sum_{k=1}^{m}f_{k}(\mathbf{x})$ with condition number $\kappa$, the accelerated version of SVRG via Catalyst [LMH15] finds an $\varepsilon$-approximate minimum with probability $1-\delta$ in time

Sherman-Morrison Formula: The following is a well-known expression for writing the inverse of rank one perturbations of matrices:

LiSSA: Linear (time) Stochastic Second-Order Algorithm

In this section, we provide an overview of LiSSA (Algorithm 1) along with its main convergence results.

Based on a recursive reformulation of the Taylor expansion (Equation 2), we may describe an unbiased estimator of the Hessian. For a matrix $A$, define $A^{-1}_{j}$ as the first $j$ terms in the Taylor expansion, i.e.,

One can also define and analyze a simpler (non-recursive) estimator based on directly sampling terms from the series in Equation (2). Theoretically, one can get similar guarantees for the estimator; however, empirically our proposed estimator exhibited better performance.

2 Algorithm

Our algorithm runs in two phases: in the first phase it runs any efficient first-order method FO for $T_{1}$ steps to shrink the function value to the regime where we can then show linear convergence for our algorithm. It then takes approximate Newton steps based on the estimator from Definition 3.1 in place of the true Hessian inverse. We use two parameters, $S_{1}$ and $S_{2}$, to define the Newton step. $S_{1}$ represents the number of unbiased estimators of the Hessian inverse we average to get better concentration for our estimator, while $S_{2}$ represents the depth to which we capture the Taylor expansion. In the algorithm, we compute the average Newton step directly, which can be computed in linear time as observed earlier, instead of estimating the Hessian inverse.

3 Main Theorem

In this section we present our main theorem which analyzes the convergence properties of LiSSA. Define $FO(M,\hat{\kappa}_{l})$ to be the total time required by a first-order algorithm to achieve accuracy $\frac{1}{4\hat{\kappa}_{l}M}$.

Consider Algorithm 1, and set the parameters as follows: $T_{1}=FO(M,\hat{\kappa}_{l})$, $S_{1}=O(\left(\hat{\kappa}^{\max}_{l}\right)^{2}\ln(\frac{d}{\delta}))$, $S_{2}\geq 2\hat{\kappa}_{l}\ln(4\hat{\kappa}_{l}).$ The following guarantee holds for every $t\geq T_{1}$ with probability $1-\delta$,

As an immediate corollary, we obtain the following:

For a GLM function $f(\mathbf{x})$, Algorithm 1 returns a point $\mathbf{x}_{t}$ such that with probability at least $1-\delta$,

We now prove our main theorem about the convergence of LiSSA (Theorem 3.3).

Note that since we use a first-order algorithm to get a solution of accuracy at least $\frac{1}{4\hat{\kappa}_{l}M}$, we have that

where $\gamma=16\hat{\kappa}^{\max}_{l}\sqrt{\frac{\ln(d\delta^{-1})}{S_{1}}}+\frac{1}{16}$.

Substituting the values of $S_{1}$ and $S_{2}$, combining Equation (3) and Lemma 3.5, and noting that $\|\nabla^{-2}f(\mathbf{x}_{t})\|\leq\hat{\kappa}_{l}$, we have that at the start of the Newton phase the following inequality holds:

It can be shown via induction that the above property holds for all $t\geq T_{1}$, which concludes the proof. ∎

Define $\chi(\mathbf{x}_{t})=\int_{0}^{1}\nabla^{2}f(\mathbf{x}^{*}+\tau(\mathbf{x}_{t}-\mathbf{x}^{*}))d\tau$. Note that $\nabla f(\mathbf{x}_{t})=\chi(\mathbf{x}_{t})(\mathbf{x}_{t}-\mathbf{x}^{*})$. Following an analysis similar to that of [Nes13], we have that

Following from the previous equations, we have that

We now analyze the above two terms $a$ and $b$ separately:

The second inequality follows from the Lipschitz bound on the Hessian. The second term can be bounded as follows:

The previous claim follows from Lemma 3.6 which shows a concentration bound on the sampled estimator and by noting that due to our assumption on the function, we have that for all $\mathbf{x},\;\|\nabla^{2}f(\mathbf{x})\|\leq 1$, and hence $\|\chi(\mathbf{x})\|\leq 1$.

Putting the above two bounds together and using the triangle inequality, we have that

First note the following statement which is a straightforward implication of our construction of the estimator:

We also know from Equation (2) that for matrices $X$ such that $\|X\|\leq 1$ and $X\succ 0$,

Since we have scaled the function such that $\|\nabla^{2}f_{k}\|\leq 1$, it follows that

Also note that since $\nabla^{2}f(\mathbf{x}_{t})\succeq\frac{I}{\hat{\kappa}_{l}}$, it follows that $\|I-\nabla^{2}f(\mathbf{x}_{t})\|\leq 1-\frac{1}{\hat{\kappa}_{l}}$. Observing the second term in the above equation,

We can now apply Theorem 2.1, which gives the following:

Setting $\varepsilon=16\hat{\kappa}^{\max}_{l}\sqrt{\frac{ln(\frac{d}{\delta})}{S_{1}}}$ gives us that the probability above is bounded by $\delta$. Now putting together the bounds and Equation (4) we get the required result. ∎

4 Leveraging Sparsity

A key property of real-world data sets is that although the input is a high dimensional vector, the number of non-zero entries is typically very low. The following theorem shows that LiSSA can be implemented in a way to leverage the underlying sparsity of the data. Our key observation is that for GLM functions, the rank one Hessian-vector product can be performed in $O(s)$ time where $s$ is the sparsity of the input $\mathbf{x}_{k}$.

For GLM functions Algorithm 1 returns a point $\mathbf{x}_{t}$ such that with probability at least $1-\delta$

We will prove the following theorem, from which Theorem 3.7 will immediately follow.

Consider Algorithm 1, let $f$ be of the form described above, and let $s$ be such that the number of non zero entries in $\mathbf{x}_{i}$ is bounded by $s$. Then each step of the algorithm can be implemented in time $O(ms+(\kappa^{\max}_{l})^{2}\kappa_{l}s)$.

LiSSA: Extensions

In this section we first describe a family of algorithms which generically couple first-order methods as sub-routines with second-order methods. In particular, we formally describe the algorithm LiSSA-Quad (Algorithm 2) and provide its runtime guarantee (Theorem 4.1). The key idea underlying this algorithm is that Newton’s method essentially reduces a convex optimization problem to solving intermediate quadratic sub-problems given by the second-order Taylor approximation at a point, i.e., the sub-problem $Q_{t}$ given by

where $\mathbf{y}\triangleq\mathbf{x}-\mathbf{x}_{t-1}$. The above ideas provide an alternative implementation of our estimator for $\nabla^{-2}f(\mathbf{x})$ used in LiSSA. Consider running gradient descent on the above quadratic $Q_{t}$, and let $\mathbf{y}^{i}_{t}$ be the $i^{th}$ step in this process. By definition we have that

It can be seen that the above expression corresponds exactly to the steps taken in LiSSA (Algorithm 2, line 8), the difference being that we use a sample of the Hessian instead of the true Hessian. Therefore LiSSA can also be interpreted as doing a partial stochastic gradient descent on the quadratic $Q_{t}$. It is partial because we have a precise estimate of gradient of the function $f$ and a stochastic estimate for the Hessian. We note that this is essential for the linear convergence guarantees we show for LiSSA.

The above interpretation indicates the possibility of using any first-order linearly convergent scheme for approximating the minimizer of the quadratic $Q_{t}$. In particular, consider any algorithm $ALG$ that, given a convex quadratic function $Q_{t}$ and an error value $\varepsilon$, produces a point $\mathbf{y}$ such that

with probability at least $1-\delta_{ALG}$, where $\mathbf{y}^{*}_{t}=\operatorname*{argmin}Q_{t}$. Let the total time taken by the algorithm $ALG$ to produce the point be $T_{ALG}(\varepsilon,\delta_{ALG})$. For our applications we require $ALG$ to be linearly convergent, i.e. $T_{ALG}$ is proportional to $\log(\frac{1}{\varepsilon})$ with probability at least $1-\delta_{ALG}$.

Given such an algorithm $ALG$, LiSSA-Quad, as described in Algorithm 2, generically implements the above idea, modifying LiSSA by replacing the inner loop with the given algorithm $ALG$. The following is a meta-theorem about the convergence properties of LiSSA-Quad.

Given the function $f(\mathbf{x})=\sum f_{i}(\mathbf{x})$ which is $\alpha$-strongly convex, let $\mathbf{x}^{*}$ be the minimizer of the function and $\{\mathbf{x}_{t}\}$ be defined as in Algorithm 2. Suppose the algorithm $ALG$ satisfies condition (5) with probability $1-\delta_{ALG}$ under the appropriate setting of parameters $ALG_{params}$. Set the parameters in the algorithm as follows: $T_{1}=T_{ALG}(1/4\alpha M)$, $T=\log\log(1/\varepsilon)$, $\delta_{ALG}=\delta/T$, where $\varepsilon$ is the final error guarantee one wishes to achieve. Then we have that after $T$ steps, with probability at least $1-\delta$,

In particular, LiSSA-Quad(ALG) produces a point $\mathbf{x}$ such that

in total time $O(T_{ALG}(\varepsilon,\delta_{ALG})\log\log(1/\varepsilon))$ with probability at least $1-\delta$ for $\varepsilon\rightarrow 0$.

Note that for GLM functions, the $\nabla Q_{t}(\mathbf{y})$ at any point can be computed in time linear in $d$. In particular, a full gradient of $Q_{t}$ can be computed in time $O(md)$ and a stochastic gradient (corresponding to a stochastic estimate of the Hessian) in time $O(d)$. Therefore, a natural choice for the algorithm $ALG$ in the above are first-order algorithms which are linearly convergent, for example SVRG, SDCA, and Acc-SDCA. Choosing a first-order algorithm FO gives us a family of algorithms LiSSA-Quad(FO), each with running time comparable to the running time of the underlying FO, up to logarithmic factors. The following corollary summarizes the typical running time guarantees for LiSSA-Quad(FO) when FO is Acc-SVRG.

Given a GLM function $f(\mathbf{x})$, if $ALG$ is replaced by Acc-SVRG [LMH15], then under a suitable setting of parameters, LiSSA-Quad produces a point $\mathbf{x}$ such that

We run the algorithm $A$ to achieve accuracy $\varepsilon^{2}$ on each of the intermediate quadratic functions $Q_{t}$, and we set $\delta_{A}=\delta/T$ which implies via a union bound that for all $t\leq T$,

Assume that for all $t<T$, $\|\mathbf{x}_{t}-\mathbf{x}^{*}\|\geq\varepsilon$ (otherwise the theorem is trivially true). Using the analysis of Newton’s method as before, we get that for all $t\leq T$,

where the second inequality follows from the analysis in the proof of Theorem 3.3 and Equation (6). We know that $\|\mathbf{x}_{0}-\mathbf{x}_{t}\|\leq\sqrt{\frac{\alpha}{M}}$ from the initial run of the first-order algorithm $FO$. Applying the above inductively and using the value of $T$ prescribed by the theorem statement, we get that $\|\mathbf{x}_{T}-\mathbf{x}^{*}\|\leq\varepsilon$. ∎

Runtime Improvement through Fast Quadratic Solvers

The previous section establishes the reduction from general convex optimization to quadratic functions. In this section we show how we can leverage the fact that for quadratic functions the running time for accelerated first-order methods can be improved in the regime when $\kappa>m>>d$. In particular, we show the following theorem.

$\kappa_{sample}(A)$ is the condition number of an $O(d\log(d))$ sized sample of A and is formally defined in Equation (11). We can now use Algorithm 4 to compute an approximate Newton step by setting $A=\nabla^{2}f(x)$ and $\mathbf{b}=\nabla f(x)$. We therefore propose LiSSA-Sample to be a variant of LiSSA-Quad where Algorithm 4 is used as the subroutine ALG and any first-order algorithm can be used in the initial phase. The following theorem bounding the running time of LiSSA-Sample follows immediately from Theorem 4.1 and Theorem 5.1.

Given a GLM function $f(\mathbf{x})=\sum_{i}f_{i}(\mathbf{x})$, let $\mathbf{x}^{*}=\operatorname*{argmin}f(\mathbf{x})$. LiSSA-Sample produces a point $\mathbf{x}$ such that

with probability at least $1-\delta$ in total time

In this section we provide a short overview of Algorithm 4. To simplify the discussion, lets consider the case when we have to compute $A^{-1}\mathbf{b}$ for a $d\times d$ matrix $A$ given as $A=\sum_{i=1}^{m}\mathbf{v}_{i}\mathbf{v}_{i}^{T}=VV^{T}$ where the $i^{th}$ column of $V$ is $\mathbf{v}_{i}$. The computation can be recast as minimization of a convex quadratic function $Q(\mathbf{y})=\frac{\mathbf{y}^{T}A\mathbf{y}}{2}+\mathbf{b}^{T}y$ and can be solved up to accuracy $\varepsilon$ in total time $\left(m+\sqrt{\kappa(A)m}\right)d\log(1/\varepsilon)$ as can be seen from Theorem 2.2 Algorithm 4 improves upon the running time bound in the case when $m>d$. In the following we provide a high level outline of the procedure which is formally described as Algorithm 4.

Given $A$ we will compute a low complexity constant spectral approximation $B$ of $A$. Specifically $B=\sum_{i=1}^{O(d\log(d))}\mathbf{u}_{i}\mathbf{u}_{i}^{T}$ and $B\preceq A\preceq 2B$. This is achieved by techniques developed in matrix sampling/sketching literature, especially those of [CLM+15]. The procedure requires solving a constant number of $O(d\log(d))$ sized linear systems, which we do via Accelerated SVRG.

We use $B$ as a preconditioner and compute $BA^{-1}\mathbf{b}$ by minimizing the quadratic $\frac{\mathbf{y}^{T}AB^{-1}\mathbf{y}}{2}+\mathbf{b}^{T}\mathbf{y}$. Note that this quadratic is well conditioned and can be minimized using gradient descent. In order to compute the gradient of the quadratic which is given by $AB^{-1}\mathbf{y}$, we again use Accelerated SVRG to solve a linear system in B.

Finally, we compute $A^{-1}\mathbf{b}=B^{-1}BA^{-1}\mathbf{b}$ using Accelerated SVRG to solve a linear system in B.

In the rest of the section we formally describe the procedure outlined and provide the necessary definitions. One key nuance we must take into account is the fact that based on our assumption we have included the regularization term into the component functions. Due to this, the Hessian does not necessarily look like a sum of rank one matrices. Of course, one can decompose the identity matrix that appears due to the regularizer as a sum of rank one matrices. However, note that the procedure above requires that each of the sub-samples must have good condition number too in order to solve linear systems on them with Accelerated SVRG. Therefore, the sub-samples generated must look like sub-samples formed from the Hessians of component functions. For this purpose we extend the procedure and the definition for leverage scores described by [CLM+15] to the case when the matrix is given as a sum of PSD matrices and not just rank one matrices. We reformulate and reprove the basic theorems proved by [CLM+15] in this context. To maintain computational efficiency of the procedure, we then make use of the fact that each of the PSD matrices actually is a rank one matrix plus the Identity matrix. We now provide the necessary preliminaries for the description of the algorithm and its analysis.

2 Preliminaries for Fast Quadratic Solver

For all the definitions and preliminaries below assume we are given a $d\times d$ PSD matrix $A\triangleq\sum_{i=1}^{m}A_{i}$ where $A_{i}$ are also PSD matrices. Let $A\cdot B\triangleq Tr(B^{T}A)$ be the standard matrix dot product. Given two matrices $A$ and $B$ we say $B$ is a $\lambda$-spectral approximation of $A$ if $\frac{1}{\lambda}A\preceq B\preceq A$.

If $B$ is a $\lambda$-spectral approximation of $A$, then $\tau_{i}(A)\leq\tau_{i}^{B}(A)\leq\lambda\tau_{i}(A)$.

When the sample is unweighted, i.e., $\mathbf{w}=\mathbf{1}$, we will simply denote the above as $Sample(I)$. We can now define $\kappa_{sample}(A,r)$ to be

The following lemma is a generalization of the leverage score sampling lemma (Lemma 4, [CLM+15]). The proof is very similar to the original proof by [CLM+15] and is included in the Appendix for completeness.

Given an error parameter $0<\varepsilon<1$, let $\mathbf{u}$ be a vector of leverage score overestimates, i.e., $\tau_{i}(A)\leq\mathbf{u}_{i}$, for all $i\in[m]$. Let $\alpha$ be a sampling rate parameter and let c be a fixed positive constant. For each matrix $A_{i}$, we define a sampling probability $p_{i}(\alpha)=\min\{1,\alpha\mathbf{u}_{i}c\log d\}$. Let $I$ be a random sample of indices drawn from $[m]$ by sampling each index with probability $p_{i}(\alpha)$. Define the weight vector $\mathbf{w}(\alpha)$ to be the vector such that $\mathbf{w}(\alpha)_{i}=\frac{1}{p_{i}(\alpha)}$. By definition of weighted samples we have that

where $x_{i}$ is a Bernoulli random variable with probability $p_{i}(\alpha)$.

If we set $\alpha=\varepsilon^{-2}$, $S=Sample(\mathbf{w}(\alpha),I)$ is formed by at most $\sum_{i}\min\{1,\alpha\mathbf{u}_{i}c\log(d)\}\leq\alpha c\log(d)\|\mathbf{u}\|_{1}$ entries in the above sum, and $\frac{1}{1+\varepsilon}S$ is a $\frac{1+\varepsilon}{1-\varepsilon}$ spectral approximation for $A$ with probability at least $1-d^{-c/3}$.

The following theorem is an analogue of the key theorem regarding uniform sampling (Theorem 1, [CLM+15]). The proof is identical to the original proof and is also included in the appendix for completeness.

Given any $A=\sum_{i=1}^{m}A_{i}$ as defined above, let $S=\sum_{j=1}^{r}X_{j}$ be formed by uniformly sampling $r$ matrices $X_{1}\ldots X_{r}\sim\{A_{i}\}$ without repetition. Define

Suppose we are given any $A=\sum_{i=1}^{m}A_{i}$ where each $A_{i}$ is of the form $A_{i}=\mathbf{v}_{i}\mathbf{v}_{i}^{T}+\lambda I$. Let $S=\sum_{j=1}^{m}X_{j}$ be formed by uniformly sampling $r$ matrices $X_{1}\ldots X_{r}\sim\{A_{i}\}$ without repetition. Define

Then $\hat{\tau}^{S}_{i}(A)\geq\tau_{i}(A)$ for all $i$, and

In the other case a similar inequality can be shown by noting via the Sherman-Morrison formula that

This proves Equation (12). A direct application of Theorem 5.7 now finishes the proof. ∎

3 Algorithms

In the following we formally state the two sub-procedures: Algorithm 4, which solves the required system, and Algorithm 3, which is the sampling routine for reducing the size of the system.

We prove the following theorem regarding the above algorithm REPEATED HALVING (Algorithm 3).

We first provide the proof of Theorem 5.1 using Theorem 5.9 and then provide the proof for Theorem 5.9. For the purpose of clarity of discourse we hide the terms that appear due to the probabilistic part of the lemmas. We can take a union bound to bound the total probability of failure, and those terms show up in the log factor.

We will first prove correctness which follows from noting that $\operatorname*{argmin}_{\mathbf{y}}Q(\mathbf{y})=BA^{-1}b$. Therefore we have that

Running Time Analysis: We now prove that the algorithm can be implemented efficiently. We first note the following two simple facts: since $B$ is a $2$-approximation of $A$, $\kappa(B)\leq 2\kappa(A)$; furthermore, the quadratic $Q(\mathbf{y})$ is $1$-strongly convex and $2$-smooth. Next, we consider the running time of implementing step 5. For this purpose we will perform gradient descent over the quadratic $Q(\mathbf{y})$. Note that the gradient of $Q(\mathbf{y})$ is given by

Define $h_{t}=Q(\mathbf{y}_{t})-\min_{\mathbf{y}}Q(\mathbf{y})$. Using the standard analysis of gradient descent, we show that the following holds true for $t>0$:

This follows directly from the gradient descent analysis which we outline below. To make the analysis easier, we define a true gradient descent series:

where $\beta\geq 1$ and $\alpha\leq 2$ are the strong convexity and smoothness parameters of $Q(\mathbf{y})$. Therefore, we have that

The running time of the above sub-procedure is bounded by the time to multiply a vector with $A$, which takes $O(md)$ time, and the time required to compute $\mathbf{v}_{t}$, which involves solving a linear system in $B$ at each step. Finally, in step 6 we once again compute the solution of a linear system in $B$. Combining these we get that the total running time is

The correctness of the algorithm is immediate from Theorems 5.8 and 5.6. The key challenge in the proof lies in proving that there is an efficient way to implement the algorithm, which we describe next.

Runtime Analysis: To analyze the running time of the algorithm we consider the running time of the computation required in step 7 of the algorithm. We will show that there exists an efficient algorithm to compute numbers $\gamma_{i}$ such that

We also know that $A^{\prime}$ is an $O(d\log(d))$ sized weighted sub-sample of $A^{\prime}$, and so it is of the form

Consider the following procedure for the computation of $\gamma^{\prime}(i)$:

For each $v_{i}$ compute $\gamma_{i}^{\prime\prime}\triangleq\sum_{j=1}^{k}<G_{j}^{\prime\prime},\mathbf{v}_{i}>^{2}$. This takes a total time of $O(kmd)$.

Here, $LIN(S,\varepsilon)$ is the running time of solving a linear system in $S$ with error $\varepsilon$. Therefore, substituting $k=O(\log(md))$, the total running time of the above algorithm is

It is now easy from the definitions of $\gamma_{i}^{\prime\prime}$ and $\gamma_{i}^{\prime}$ to see that

Setting $\gamma_{i}=4(\gamma_{i}^{\prime\prime}+k\varepsilon\|\mathbf{v}_{i}\|^{2})$, it follows from the definitions that

Now setting $\varepsilon=\frac{1}{8k\|V\|_{F}^{2}}$ satisfies the required inequality for $\gamma_{i}$. This implies that when sampling from $\gamma_{i}$, we will have a $2$-approximation, and the number of matrices will be bounded by $O(d\log(d))$.

Putting the above arguments together we get that the total running time of the procedure is bounded by

Condition Number Independent Algorithms

In this section we state our main result regarding self-concordant functions and provide the proof. The following is the definition of a self-concordant function:

Our main theorem regarding self-concordant functions is as follows:

Let $0<r<1$, let $\gamma\geq 1$, and let $\nu$ be a constant depending on $\gamma,r$. Set $\eta=10(1-r)^{2}$, $S_{1}=c_{r}=\frac{50}{(1-r)^{4}}$, where $c_{r}$ is a constant depending on $r$, and $T=\frac{f(\mathbf{x}_{1})-f(\mathbf{x}^{*})}{\nu}$. Then, after $t>T$ steps, the following linear convergence guarantee holds between two full gradient steps for Algorithm 5:

Moreover, the time complexity of each full gradient step is $O(md+c_{r}d^{2})$, where $c_{r}$ is a constant depending on $r$ (independent of the condition number).

To describe the algorithm we first present the requisite definitions regarding self-concordant functions.

An excellent reference for this material is the lecture notes on this subject by [Nem04].

A key object in the analysis of self-concordant functions is the notion of a Dikin ellipsoid, which is the unit ball around a point $\mathbf{x}$ in the norm given by the Hessian $\|\cdot\|_{\nabla^{2}f(\mathbf{x})}$ at the point. We will refer to this norm as the local norm around a point and denote it as $\|\cdot\|_{\mathbf{x}}$. Formally, we have:

The Dikin ellipsoid of radius $r$ centered at a point $\mathbf{x}$ is defined as

One of the key properties of self-concordant functions is that inside the Dikin ellipsoid, the function is well conditioned with respect to the local norm at the center, and furthermore, the function is smooth. The following lemmas makes these notions formal, and the proofs of these lemmas can be found in the lecture notes of [Nem04].

For all $\mathbf{h}$ such that $\|\mathbf{h}\|_{\mathbf{x}}<1$ we have that

For all $\mathbf{h}$ such that $\|\mathbf{h}\|_{\mathbf{x}}<1$ we have that

Another key quantity which is used both as a potential function as well as a dampening for the step size in the analysis of Newton’s method in general is the Newton decrement which is defined as $\lambda(\mathbf{x})\triangleq\|\nabla f(\mathbf{x})\|_{\mathbf{x}}^{*}=\sqrt{\nabla f(\mathbf{x})^{\top}\nabla^{-2}f(\mathbf{x})\nabla f(\mathbf{x})}$. The following lemma quantifies how $\lambda(\mathbf{x})$ behaves as a potential by showing that once it drops below 1, it ensures that the minimum of the function lies in the current Dikin ellipsoid. This is the property which we use crucially in our analysis.

2 Condition Number Independent Algorithms

In this section we describe an efficient linearly convergent method (Algorithm 5) for optimization of self-concordant functions for which the running time is independent of the condition number. We have not tried to optimize the complexity of the algorithms in terms of $d$ as our main focus is to make it condition number independent.

The key idea here is the ellipsoidal view of Newton’s method, whereby we show that after making a constant number of full Newton steps, one can identify an ellipsoid and a norm such that the function is well conditioned with respect to the norm in the ellipsoid. This is depicted in Figure 1.

At this point one can run any desired first-order algorithm. In particular, we choose SVRG and prove its fast convergence. Algorithm 6 (described in the appendix) states the modified SVRG routine for general norms used in Algorithm 5.

We now state and prove the following theorem regarding the convergence of Algorithm 5.

It follows from Lemma 6.8 that at $t=\frac{f(\mathbf{x}_{1})-f(\mathbf{x}^{*})}{\nu}$, the minimizer is contained in the Dikin ellipsoid of radius $r$, where $\nu$ is a constant depending on $\gamma,r$. This fact, coupled with Lemma 6.7, shows that the function satisfies the following property with respect to $W_{\mathbf{x}_{t}}$:

Using the above fact along with Lemma 6.9, and substituting for the parameters, concludes the proof. ∎

We describe the Algorithm N-SVRG (Algorithm 6 in the appendix). Since the algorithm and the following lemmas are minor variations of their original versions, we include the proofs of the following lemmas in the appendix for completeness.

Let $f$ be a self-concordant function over $\mathcal{K}$, let $0<r<1$, and consider $\mathbf{x}\in\mathcal{K}$. Let $W_{r}(\mathbf{x})$ be the Dikin ellipsoid of radius $r$, and let $\alpha=(1-r)^{2}$ and $\beta=(1-r)^{-2}$. Then, for all $\mathbf{h}$ s.t. $\mathbf{x}+\mathbf{h}\in W_{r}(\mathbf{x})$,

Let $f$ be a self-concordant function over $\mathcal{K}$, let $0<r<1$, let $\gamma\geq 1$, and consider following the damped Newton step as described in Algorithm 5 with initial point $\mathbf{x}_{1}$. Then, the number of steps $t$ of the algorithm before the minimizer of $f$ is contained in the Dikin ellipsoid of radius $r$ of the current iterate, i.e. $\mathbf{x}^{*}\in W_{r}(\mathbf{x}_{t})$, is at most $t=\frac{f(\mathbf{x}_{1})-f(\mathbf{x}^{*})}{\nu}$, where $\nu$ is a constant depending on $\gamma,r$.

Let $f$ be a convex function. Suppose there exists a convex set $\mathcal{K}$ and a positive semidefinite matrix $A$ such that for all $\mathbf{x}\in\mathcal{K}$, $\alpha A\preceq\nabla^{2}f(\mathbf{x})\preceq\beta A$. Then the following holds between two full gradient steps of Algorithm 6:

Experiments

In this section we describe our experiments and choice of parameters in detail. Table 2 provides details regarding the data sets chosen for the experiments. To make sure our functions are scaled such that the norm of the Hessian is bounded, we scale the above data set points to unit norm.

2 Comparison with Standard Algorithms

In Figures 2 and 4 we present comparisons between the efficiency of our algorithm with different standard and popular algorithms. In both cases we plot $\log(CurrentValue-OptimumValue)$. We obtained the optimum value for each case by running our algorithm for a long enough time until it converged to the point of machine precision.

Epoch Comparison: In Figure 2, we compare LiSSA with SVRG and SAGA in terms of the accuracy achieved versus the number of passes over the data. To compute the number of passes in SVRG and SAGA, we make sure that the inner stochastic gradient iteration in both the algorithms counts as exactly one pass. This is done because although it accesses gradients at two different points, one of them can be stored from before in both cases. The outer full gradient in SVRG counts as one complete pass over the data. We set the number of inner iterations of SVRG to $2m$ for the case when $\lambda=1/m$, and we parameter tune the number of inner iterations when $\lambda=10/m$. The stepsizes for all of the algorithms are parameter tuned by an exhaustive search over the parameters.

Time Comparison: For the comparison with respect to time (Figure 4), we consider the following algorithms: AdaGrad [DHS11], BFGS [Bro70, Fle70, Gol70, Sha70], gradient descent, SGD, SVRG [JZ13] and SAGA [DBLJ14]. The $\log(Error)$ is plotted as a function of the time elapsed from the start of the run for each algorithm. We next describe our choice of parameters for the algorithms. For AdaGrad we used the faster diagonal scaling version proposed by [DHS11]. We implemented the basic version of BFGS with backtracking line search. In each experiment for gradient descent, we find a reasonable step size using parameter tuning. For stochastic gradient descent, we use the variable step size $\eta_{t}=\gamma/\sqrt{t}$ which is usually the prescribed step size, and we hand tune the parameter $\gamma$. The parameters for SVRG and SAGA were chosen in the same way as before.

Choice of Parameters for LiSSA: To pick the parameters for our algorithm, we observe that it exhibits smooth behavior even in the case of $S_{1}=1$, so this is used for the experiments. However, we observe that increasing $S_{2}$ has a positive effect on the convergence of the algorithm up to a certain point, as a higher $S_{2}$ leads to a larger per-iteration cost. This behavior is consistent with the theoretical analysis. We summarize the comparison between the per-iteration convergence and the value of $S_{2}$ in Figure 5. As the theory predicts $S_{2}$ to be of the order $\kappa\ln(\kappa)$, for our experiments we determine an estimate for $\kappa$ and set $S_{2}$ to around $\kappa\ln(\kappa)$. This value is typically equal to $m$ in our experiments. We observe that setting $S_{2}$ in this way resulted in the experimental results displayed in Figure 2.

Comparison with Second-Order Methods: Here we present details about the comparison between LiSSA, NewSamp [EM15], and standard Newton’s method, as displayed in Figure 3. We perform this experiment on the MNIST data set and show the convergence properties of all three algorithms over time as well as over iterations. We could not replicate the results of NewSamp on all of our data sets as it sometimes seems to diverge in our experiments. For logistic regression on the MNIST data set, we could get it to converge by setting the value of $S_{1}$ to be slightly higher. We observe as is predicted by the theory that when compared in terms of the number of iterations, NewSamp and LiSSA perform similarly, while Newton’s method performs the best as it attains a quadratic convergence rate. This can be seen in Figure 3. However, when we consider the performance in terms of time for these algorithms, we see that LiSSA has a significant advantage.

Comparison with Accelerated First-Order Methods: Here we present experimental results comparing LiSSA with a popular accelerated first-order method, APCG [LLX14], as seen in Figure 6. We ran the experiment on the RealSIM data set with three settings of $\lambda=10^{-5},\ 10^{-6},\ 10^{-7}$, to account for the high condition number setting. We observe a trend that can be expected from the runtime guarantees of the algorithms. When $\lambda$ is not too low, LiSSA performs better than APCG, but as $\lambda$ gets very low we see that APCG performs significantly better than LiSSA. This is not surprising when considering that the running time of APCG grows proportional to $\sqrt{\kappa m}$, whereas for LiSSA the running time can at best be proportional to $\kappa$. We note that for accelerated first-order methods to be useful, one needs the condition number to be quite large which is not often the case for applications. Nevertheless, we believe that an algorithm with running time guarantees similar to LiSSA-Sample can get significant gains in these settings, and we leave this exploration as future work.

Acknowledgements

The authors would like to thank Zeyuan Allen-Zhu, Dan Garber, Haipeng Luo and David McAllester for several illuminating discussions and suggestions. We would especially like to thank Aaron Sidford and Yin Tat Lee for discussions regarding the matrix sampling approach.

References

Appendix A Remaining Proofs

Here we provide proofs for the remaining lemmas.

The proof is almost identical to the proof of Lemma 4 by [CLM+15]. As in there we will use the following lemma on matrix concentration, which appears as Lemma 11 in the work of [CLM+15].

Let $Y_{1}\ldots Y_{k}$ be independent random positive semidefinite matrices of size $d\times d$. Let $Y=\sum Y_{i}$ and let $Z=E[Y]$. If $Y_{i}\preceq R\cdot Z$, then

For every matrix $A_{i}$ choose $Y_{i}=\frac{A_{i}}{p_{i}}$ with probability $p_{i}$ and otherwise. Therefore we need to bound $\sum Y_{i}$. Note that $E[Y_{i}]=A$. We will now show that

which will finish the proof by a direct application of Lemma A.1. First we will show that

We need to show that $\forall\mathbf{x}\;\mathbf{x}^{T}A_{i}\mathbf{x}\leq\tau_{i}(A)\mathbf{x}^{T}A\mathbf{x}$. By noting that the $A_{i}$ are PSD we have that if $\mathbf{x}^{T}A\mathbf{x}=0$ then $\forall i,\mathbf{x}^{T}A_{i}\mathbf{x}=0$. We can now consider $\mathbf{x}=A^{+/2}\mathbf{y}$. Therefore, we need to show that

This is true because the maximum eigenvalue of $A^{+/2}A_{i}A^{+/2}$ is bounded by its $Trace(A^{+/2}A_{i}A^{+/2})$ (due to positive semidefiniteness), which is equal to $\tau_{i}(A)$ by definition. Therefore when $p_{i}<1$, i.e., $\alpha c\mathbf{u}_{i}\log(d)<1$, the facts $\tau_{i}(A)\leq\mathbf{u}_{i}$ immediately provide

When $p_{i}=1$ the above does not hold, but we can essentially replace $Y_{i}=A_{i}$ with $c\log d\varepsilon^{-2}$ variables each equal to $\frac{A_{i}}{c\log d\varepsilon^{-2}}$, each being sampled with probability 1. This does not change $E[\sum Y_{i}]$ but proves concentration. Now a direct application of Lemma A.1 finishes the proof. Also note that a standard Chernoff bound proves the required bound on the sample size. ∎

A.2 Proof of Lemma 5.7

A.3 Proof of Ellipsoidal Cover Lemma

Let $\lambda(\mathbf{x})$ be the Newton decrement at $\mathbf{x}$. By Lemma 6.6, we know that if $\lambda(\mathbf{x})\leq\frac{r}{1+r}<1$, then

Consider a single iteration of the algorithm at time $t$. If $\lambda(\mathbf{x}_{t})\leq\frac{r}{1+r}$, then we may conclude that $\mathbf{x}^{*}\in W_{r}(\mathbf{x}_{t})$. Therefore, it is only when $\lambda(\mathbf{x}_{t})>\frac{r}{1+r}$ that $\mathbf{x}^{*}$ may not be contained within $W_{r}(\mathbf{x}_{t})$. Since our update is of the form

where the first inequality follows from Lemma 6.5. It now follows that

steps, we can guarantee that we have arrived at $\mathbf{x}_{t}$ such that $\mathbf{x}^{*}\in W_{r}(\mathbf{x}_{t})$. ∎

A.4 Proof of SVRG Lemma

The proof follows the original proof of SVRG [JZ13] with a few modifications to take the general norm into account. For any $i$, consider

We know that $g_{i}(\mathbf{x}^{*})=\min_{\mathbf{x}}g_{i}(\mathbf{x})$ since $\nabla g_{i}(\mathbf{x}^{*})=0$. Therefore,

The second inequality follows from smoothness of $g_{i}$, where we note that $g_{i}$ is as smooth as $f_{i}$. Summing the above inequality over $i$ and setting $\mathbf{x}^{*}$ to be the minimum of $f$ so that $\nabla f(\mathbf{x}^{*})=0$, we get

The above inequalities follow by noting the following three facts:

$\|a+b\|_{A^{-1}}^{2}\leq 2\|a\|_{A^{-1}}^{2}+2\|b\|_{A^{-1}}^{2}\enspace.$

$\mbox{\bf E}\left[\|X-\mbox{\bf E}X\|_{A^{-1}}^{2}\right]\leq\mbox{\bf E}\|X\|_{A^{-1}}^{2}\enspace.$

Now note that conditioned on $\mathbf{x}_{t-1}$, we have that $\mbox{\bf E}v_{t}=\nabla f(\mathbf{x}_{t-1})$, and so

The second inequality uses the strong convexity property. Therefore, we have that