Sub-sampled Cubic Regularization for Non-convex Optimization

Introduction

In this paper we address the problem of minimizing an objective function of the form

In this work, we focus our attention on trust region methods to optimize Eq. 1. These methods construct and optimize a local model of the objective function within a region whose radius depends on how well the model approximates the real objective. One of the keys for efficiency of these methods is to pick a model that is comparably easy to optimize, such as a quadratic function (Conn et al., 2000). Following the trust region paradigm, cubic regularization methods (Nesterov & Polyak, 2006; Cartis et al., 2011a) suggest finding the step ${\bf s}_{k}$ that minimizes a cubic model of the form

(Nesterov & Polyak, 2006) were able to show that, if the step is computed by globally minimizing the cubic model and if the Hessian ${\bf H}_{k}$ is globally Lipschitz continuous, Cubic regularization methods possess the best known worst case complexity to solve Eq. 1: an overall worst-case iteration count of order $\epsilon^{-3/2}$ for generating $\|\nabla f({\bf x}_{k})\|\leq\epsilon$, and of order $\epsilon^{-3}$ for achieving approximate nonnegative curvature. However, minimizing Eq. 2 in an exact manner impedes the performance of this method for machine learning applications as it requires access to the full Hessian matrix. More recently, (Cartis et al., 2011a) presented a method (hereafter referred to as ARC) which relaxed this requirement by assuming that one can construct an approximate Hessian ${\bf B}_{k}$ that is sufficiently close to ${\bf H}_{k}$ in the following way:

Furthermore, they showed that it is sufficient to find an approximate minimizer by applying a Lanczos method to build up evolving Krylov spaces, which can be constructed in a Hessian-free manner, i.e. by accessing the Hessians only indirectly via matrix-vector products. However there are still two obstacles for the application of ARC in the field of machine learning: (1) The cost of the Lanczos process increases linearly in $n$ and is thus not suitable for large datasets and (2) there is no theoretical guarantee that quasi-Newton approaches satisfy Eq. 3 and (Cartis et al., 2011a) do not provide any alternative approximation technique.

In this work, we make explicit use of the finite-sum structure of Eq. 1 by applying a sub-sampling technique in order to provide guarantees for machine learning applications. Towards this goal, we make the following contributions:

We provide a theoretical Hessian sampling scheme that is guaranteed to satisfy Eq. 3 with high probability.

We extend the analysis to inexact gradients and prove that the convergence guarantees of (Nesterov & Polyak, 2006; Cartis et al., 2011a) can be retained.

Since the dominant iteration cost lie in the construction of the Lanczos process and increase linearly in $n$, we lower the computational cost significantly by reducing the number of samples used in each iteration.

Finally, we provide experimental results demonstrating significant speed-ups compared to standard first and second-order optimization methods for various convex and non-convex objectives.

Related work

In large-scale learning, when $n\gg d$ most of the computational cost of traditional deterministic optimization methods is spent in computing the exact gradient information. A common technique to address this issue is to use sub-sampling in order to compute an unbiased estimate of the gradient. The simplest instance is Stochastic Gradient Descent (SGD) whose convergence does not depend on the number of datapoints $n$. However, the variance in the stochastic gradient estimates slows its convergence down. The work of (Friedlander & Schmidt, 2012) explored a sub-sampling technique for gradient descent in the case of convex functions, showing that it is possible to maintain the same convergence rate as full-gradient descent by carefully increasing the sample size over time. Another way to recover a linear rate of convergence for strongly-convex functions is to use variance-reduced methods (Johnson & Zhang, 2013; Defazio et al., 2014; Roux et al., 2012; Hofmann et al., 2015; Daneshmand et al., 2016). Recently, the convergence of SGD and its variance-reduced counterparts has also been extended to non-convex functions (Ghadimi & Lan, 2013; Reddi et al., 2016a) but the techniques used in these papers require using a randomized sampling scheme which is different from what is typically used in practice. Furthermore, the guarantees these methods provide are only in terms of convergence to critical points. However, the work of (Ge et al., 2015; Sun et al., 2015) recently showed that SGD can achieve stronger guarantees in the case of strict saddle functions. Yet, the convergence rate has a polynomial dependency to the dimension $d$ and the smallest eigenvalue of the Hessian which can make this method fairly impractical.

For second-order methods, the problem of avoiding saddle points is even worse as they might be attracted by saddle points or even points of local maximizers (Dauphin et al., 2014). Another predominant issue is the computation (and perhaps storage) of the Hessian matrix, which requires $O(nd^{2})$ operations as well as computing the inverse of the Hessian, which requires $O(d^{3})$ computations. Quasi-Newton methods such as the well-known (L-)BFGS algorithm partially address this issue by requiring $O(nd+d^{2})$ per-iteration cost (Nesterov, 2004) instead of $O(nd^{2}+d^{3})$. An increasingly popular alternative is to use sub-sampling techniques to approximate the Hessian matrix, such as done for example in (Byrd et al., 2011) and (Erdogdu & Montanari, 2015). The latter method, named NewSamp, approximates the Hessian with a low-rank approximation which reduces the complexity per iteration to $O(nd+|S|d^{2})$ with $|S|$ being the sample size Note that this method still requires $O(nd)$ computation for the gradient as it only subsamples the Hessian.. Although this is a significant reduction in terms of complexity, NewSamp yields a composite convergence rate: quadratic at first but only linear near the minimizer. Unlike NewSamp, our sampling scheme yields a locally quadratic rate of convergence (as well as faster global convergence). Our analysis also does not require using exact gradients and can thus further reduce the complexity per iteration.

Trust region methods are among the most effective algorithmic frameworks to avoid pitfalls such as local saddle points in non-convex optimization. Classical versions iteratively construct a local quadratic model and minimize it within a certain radius wherein the model is trusted to be sufficiently similar to the actual objective function. This is equivalent to minimizing the model function with a suitable quadratic penalty term on the stepsize. Thus, a natural extension is the cubic regularization method introduced by (Nesterov & Polyak, 2006) that uses a cubic over-estimator of the objective function as a regularization technique for the computation of a step to minimize the objective function. The drawback of their method is that it requires computing the exact minimizer of Eq. 2, thus requiring the exact gradient and Hessian matrix. However finding a global minimizer of the cubic model $m_{k}({\bf s})$ may not be essential in practice and doing so might be prohibitively expensive from a computational point of view. (Cartis et al., 2011a) introduced a method named ARC which relaxed this requirement by letting ${\bf s}_{k}$ be an approximation to the minimizer. The model defined by the adaptive cubic regularization method introduced two further changes. First, instead of computing the exact Hessian ${\bf H}_{k}$ it allows for a symmetric approximation ${\bf B}_{k}$. Second, it introduces a dynamic positive parameter $\sigma_{k}$ instead of using the global Lipschitz constant $L$.

There have been efforts to further reduce the computational complexity of this problem. For example, (Agarwal et al., 2016) refined the approach of (Nesterov & Polyak, 2006) to return an approximate local minimum in time which is linear in the input representation. Similar improvements have been made by (Carmon & Duchi, 2016) and (Hazan & Koren, 2016). These methods provide alternatives to minimize the cubic model and can thus be seen as complementary to our approach. Finally, the work of (Blanchet et al., 2016) proposed a stochastic variant of a trust region method but their analysis does not specify any accuracy level required for the estimation of the stochastic Hessian. (Cartis & Scheinberg, 2017) also analyzed a probabilistic cubic regularization variant that allows approximate second-order models but they did not provide an explicit derivation of sampling conditions.

Formulation

We are interested in optimizing Eq. 1 in a large-scale setting when the number of datapoints $n$ is very large such that the cost of solving Eq. 2 exactly becomes prohibitive. In this regard we identify a sampling scheme that allows us to retain the convergence results of deterministic trust region and cubic regularization methods, including quadratic local convergence rates and global second-order convergence guarantees as well as worst-case complexity bounds. A detailed theoretical analysis is given in Section 4. Here we shall first state the algorithm itself and elaborate further on the type of local nonlinear models we employ as well as how these can be solved efficiently.

Instead of using deterministic gradient and Hessian information as in Eq. 2, we use unbiased estimates of the gradient and Hessian constructed from two independent sets of points denoted by $S_{g}$ and $S_{B}$. We then construct a local cubic model that is (approximately) minimized in each iteration:

where ${\bf g}_{k}:=\frac{1}{|S_{g}|}\sum_{i\in S_{g}}\nabla f_{i}({\bf x}_{k})$ and ${\bf B}_{k}:=\frac{1}{|S_{B}|}\sum_{i\in S_{B}}\nabla^{2}f_{i}({\bf x}_{k})$.

The model derivative with respect to ${\bf s}_{k}$ is defined as:

2 Algorithm

Our Sub-sampled Cubic Regularization approach (SCR) is presented in Algorithm 1. At iteration step $k$, we sub-sample two sets of datapoints from which we compute a stochastic estimate of the gradient and the Hessian. We then solve the problem in Eq. 4 approximately using the method described in Section 3.4 and update the regularization parameter $\sigma_{k}$ depending on how well the model approximates the real objective. In particular, very successful steps indicate that the model is (at least locally) an adequate approximation of the objective such that the penalty parameter is decreased in order to allow for longer steps. For unsuccessful iterations we proceed exactly the opposite way. Readers familiar with trust region methods might see that one can interpret the penalty parameter $\sigma_{k}$ as inversely proportional to the trust region radius $\delta_{k}$.

3 Exact model minimization

Solving Eq. 4 requires minimizing an unconstrained non-convex problem that may have isolated local minima. As shown in (Cartis et al., 2011a) the global model minimizer ${\bf s}_{k}^{*}$ is characterized by following systems of equations,

In order to find a solution we can express ${\bf s}_{k}^{*}:={\bf s}_{k}(\lambda_{k}^{*})=-({\bf B}_{k}+\lambda^{*}_{k}{\bf I})^{-1}{\bf g}_{k}$, apply this in the second equation of (9) and obtain a univariate, nonlinear equation in $\lambda_{k}$

Furthermore, we need $\lambda^{*}_{k}\geq\max\{-\lambda_{1}({\bf B}_{k}),0\}$, where $\lambda_{1}({\bf B}_{k})$ is the leftmost eigenvalue of ${\bf B}_{k}$, in order to guarantee the semi-positive definiteness of $({\bf B}_{k}+\lambda^{*}_{k}{\bf I})$.

Thus, computing the global solution of $m_{k}$ boils down to finding the root of Eq. 10 in the above specified range of $\lambda_{k}$. The problem can be solved by Newton’s method, which involves factorizing ${\bf B}_{k}+\lambda_{k}{\bf I}$ for various $\lambda_{k}$ and is thus prohibitively expensive for large problem dimensions $d$. See Section 6.2 in (Cartis et al., 2011a) for more details. In the following Section we instead explore an approach to approximately minimize the model while retaining the convergence guarantees of the exact minimization.

4 Approximate model minimization

(Cartis et al., 2011a) showed that it is possible to retain the remarkable properties of the cubic regularization algorithm with an inexact model minimizer. A necessary condition is that ${\bf s}_{k}$ satisfies the two requirements stated in A1.

Note that the first equation is equal to $\nabla_{s}m_{k}({\bf s}_{k})^{\intercal}{\bf s}_{k}=0$ and the second to ${\bf s}_{k}^{\intercal}\nabla^{2}_{s}m_{k}({\bf s}_{k}){\bf s}_{k}\geq 0$.

As shown in (Cartis et al., 2011a) Lemma 3.2, the global minimizer of $m_{k}({\bf s}_{k})$ in a Krylov subspace $\mathcal{K}_{k}:=\text{span}\{{\bf g}_{k},{\bf H}_{k}{\bf g}_{k},{\bf H}_{k}^{2}{\bf g}_{k},...\}$ satisfies this assumption independent of the subspace dimension. This comes in handy, as minimizing $m_{k}$ in the Krylov subspace only involves factorizing a tri-diagonal matrix, which can be done at the cost of $O(d)$. However, a Lanczos-type method must be used in order to build up an orthogonal basis of this subspace which typically involves one matrix-vector product ($O((2d-1)n)$) per additional subspace dimension (see Chapter 5 in (Conn et al., 2000) for more details).

Thus, in order to keep the per iteration cost of SCR low and in accordance to ARC, we apply the following termination criterion to the Lanczos process in the hope to find a suitable trial step before $\mathcal{K}_{k}$ is of dimensionality $d$.

For each outer iteration $k$, assume that the Lanczos process stops as soon as some Lanczos iteration $i$ satisfies the criterion

where $\theta_{k}=\kappa_{\theta}\min(1,{}\left\|{\bf s}_{i,k}\right\|),\>\kappa_{\theta}\in(0,1)$.

However, we argue that especially for high dimensional problems, the cost of the Lanczos process may significantly slow down cubically regularized methods and since this cost increases linearly in $n$, carefully sub-sampled versions are an attractive alternative.

Theoretical analysis

In this section, we provide the convergence analysis of SCR. For the sake of brevity, we assume Lipschitz continuous Hessians right away but note that a superlinear local convergence result as well as the global first-order convergence theorem can both be obtained without the former assumption.

First, we lay out some critical assumptions regarding the gradient and Hessian approximations. Second, we show that one can theoretically satisfy these assumptions with high probability by sub-sampling first- and second-order information. Third, we give a condensed convergence analysis of SCR which is widely based on (Cartis et al., 2011a), but adapted for the case of stochastic gradients. There, we show that the local and global convergence properties of ARC can be retained by sub-sampled versions at the price of slightly worse constants.

By use of the triangle inequality, it follows that these assumptions hold for all ${\bf g}$ and ${\bf H}$, independent of the sample size. Furthermore, note that the Hessian and gradient norms are uniformly bounded as a consequence of A3.

In each iteration, the Hessian approximation ${\bf B}_{k}$ shall satisfy condition AM.4 from (Cartis et al., 2011a), which we restate here for the sake of completeness.

We explicitly stress the fact that this condition is stronger than the well-known Dennis Moré Condition. While quasi-Newton approximations satisfy the latter, there is no theoretical guarantee that they also satisfy the former (Cartis et al., 2011a). Furthermore, any sub-sampled gradient shall satisfy the following condition.

2 Sampling Conditions

As detailed in the Appendix, the following Lemma arises from the Vector Bernstein Inequality.

Let the sub-sampled gradient ${\bf g}_{k}$ be defined as in Eq. 4. For $\epsilon\leq 2\kappa_{f}$ we have with probability $(1-\delta)$ that

It constitutes a non-asymptotic bound on the deviation of the gradient norms that holds with high probability. Note how the accuracy of the gradients increases in the sample size. This bound yields the following condition.

then ${\bf g}_{k}$ satisfies the sufficient agreement condition A5 with probability $(1-\delta)$.

2.2 Hessian Sampling

In analogy to the gradient case, we use the matrix version of Bernstein’s Inequality to derive the following Lemma.

Let the sub-sampled Hessian ${\bf B}$ be defined as in Eq. 4. For $\epsilon\leq 4\kappa_{g}$ we have with probability $(1-\delta)$ that

This, in turn, can be used to derive a Hessian sampling condition that is guaranteed to satisfy the sufficient agreement condition (A4) with high probability.

then ${\bf B}_{k}$ satisfies the strong agreement condition A4 with probability $(1-\delta)$.

As expected, the required sample size grows in the problem dimensionality $d$ and in the Lipschitz constants $\kappa_{f}$ and $\kappa_{g}$. Finally, as outlined in the Appendix (Lemma 24), the samples size is eventually equal to the full sample size $n$ as SCR converges and thus we have

3 Convergence Analysis

The entire analysis of cubically regularized methods is prohibitively lengthy and we shall thus establish only the crucial properties that ensure global, as well as fast local convergence and improve the worst-case complexity of these methods over standard trust region approaches. Next to the cubic regularization term itself, these properties arise mainly from the penalty parameter updates and step acceptance criteria of the ARC framework, which give rise to a good relation between regularization and stepsize. Further details can be found in (Cartis et al., 2011a).

First, we note that the penalty parameter sequence $\{\sigma_{k}\}$ is guaranteed to stay within some bounded positive range, which is essentially due to the fact that SCR is guaranteed to find a successful step as soon as the penalty parameter exceeds some critical value $\sigma_{sup}$.

where $\sigma_{\inf}$ is defined in Step 7 of Algorithm 1 and

Furthermore, for any successful iteration the objective decrease can be directly linked to the model decrease via the step acceptance criterion in Eq. 8. The latter, in turn, can be shown to be lower bounded by the stepsize which combined gives the following result.

Suppose that ${\bf s}_{k}$ satisfies A1. Then, for all successful iterations $k\geq 0$

Finally, the termination criterion (13) also guarantees step sizes that do not become too small compared to the respective gradient norm which leads to the following Lemma.

Let A3, A4 and A5 hold. Furthermore, assume the termination criterion TC (A2) and suppose that ${\bf x}_{k}\rightarrow{\bf x}^{*},\text{as }k\rightarrow\infty$. Then, for all sufficiently large successful iterations, ${\bf s}_{k}$ satisfies

where $\kappa_{s}$ is the positive constant

3.2 Local convergence

We here provide a proof of local convergence for any sampling scheme that satisfies the conditions presented in Theorem 7 and Theorem 9 as well as the additional condition that the sample size does not decrease in unsuccessful iterations. We show that such sampling schemes eventually yield exact gradient and Hessian information. Based upon this observation, we obtain the following local convergence result (as derived in the Appendix).

Let A3 hold and assume that ${\bf g}_{k}$ and ${\bf B}_{k}$ are sampled such that 17 and 19 hold and $|S_{g,k}|$ and $|S_{B,k}|$ are not decreased in unsuccessful iterations. Furthermore, let $s_{k}$ satisfy A1 and

where ${\bf H}({\bf x}^{*})$ is positive definite. Moreover, assume the stopping criterion TC (A2). Then,

That is, ${\bf x}_{k}$ converges in q-quadratically to ${\bf x}^{*}$ as $k\rightarrow\infty$ with high probability.

3.3 Global convergence to first-order critical point

Lemma 10 and 11 allow us to lower bound the function decrease of a successful step in terms of the full gradient $\nabla f_{k}$ (as we will shorty detail in Eq. 31). Combined with Lemma 10, this allows us to give deterministic global convergence guarantees using only stochastic first order information.

Let A1, A3, A4 and A5 hold. Furthermore, let $\{f({\bf x}_{k})\}$ be bounded below by some $f_{\inf}>-\infty$. Then

3.4 Global convergence to second-order critical point

Unsurprisingly, the second-order convergence guarantee relies mainly on the use of second-order information so that the stochastic gradients do neither alter the result nor the proof as it can be found in Section 5 of (Cartis et al., 2011a). We shall restate it here for the sake of completeness.

Let A3, A4 and A5 hold. Furthermore, let $\{f({\bf x}_{k})\}$ be bounded below by $f_{\inf}$ and ${\bf s}_{k}$ be a global minimizer of $m_{k}$ over a subspace $\mathcal{L}_{k}$ that is spanned by the columns of the $d\times l$ orthogonal matrix ${\bf Q}_{k}$. As ${\bf B}\rightarrow{\bf H}$ asymptotically (Eq. 20), any subsequence of negative leftmost eigenvalues $\{\lambda_{\min}({\bf Q}_{k}^{\intercal}{\bf H}({\bf x}_{k}){\bf Q}_{k})\}$ converges to zero for sufficiently large, successful iterations. Hence

3.5 Worst-case iteration complexity

For the worst-case analysis we shall establish the two disjoint index sets $\mathcal{U}_{j}$ and $\mathcal{S}_{j}$, which represent the un- and successful SCR iterations that have occurred up to some iteration $j>0$, respectively. As stated in Lemma 10 the penalty parameter $\sigma_{k}$ is bounded above and hence SCR may only take a limited number of consecutive unsuccessful steps. As a consequence, the total number of unsuccessful iterations is at most a problem dependent constant times the number of successful iterations.

For any fixed $j\geq 0$, let Lemma 10 hold. Then we have that

Regarding the number of successful iterations we have already established the two key ingredients: (i) a sufficient function decrease in each successful iteration (Lemma 11) and (ii) a step size that does not become too small compared to the respective gradient norm (Lemma 12), which is essential to driving the latter below $\epsilon$ at a fast rate. Combined they give rise to the guaranteed function decrease for successful iterations

which already contains the power of 3/2 that appears in the complexity bound. Finally, by summing over all successful iterations one obtains the following, so far best know, worst-case iteration bound to reach $\epsilon$ first-order criticality.

Let A1, A3, A4 and A5 hold. Furthermore, be $\{f({\bf x}_{k})\}$ bounded below by $f_{\inf}$ and TC applied (A2). Then, for $\epsilon>0$ the total number of iterations SCR takes to generate the first iterate $j$ with ${}\left\|\nabla f({\bf x}_{j+1})\right\|\leq\epsilon$, and assuming $\epsilon\leq 1$, is

Note that the constants $\kappa_{i}$ and $\kappa_{j}$ involved in this upper bound both increase in the gradient inaccuracy $M$ and the Hessian inaccuracy $C$ (via $\kappa_{s}$ and $\sigma_{\sup}$), such that more inaccuracy in the sub-sampled quantities may well lead to an increased overall number of iterations.

Finally, we want to point out that similar results can be established regarding a second-order worst-case complexity bound similar to Corollary 5.5 in (Cartis et al., 2011b), which we do not prove here for the sake of brevity.

Experimental results

In this section we present experimental results on real-world datasets where $n\gg d\gg 1$. They largely confirm the analysis derived in the previous section. Please refer to the Appendix for more detailed results and experiments on higher dimensional problems.

We implement SCR as stated in Algorithm 1 and note the following details. Following (Erdogdu & Montanari, 2015), we require the sampling conditions derived in Section 4 to hold with probability $O(1-1/d)$, which yields the following practically applicable sampling schemes

The positive constants $C$ and $M$ can be used to scale the sample size to a reasonable portion of the entire dataset and can furthermore be used to offset $\kappa_{g}$ and $\kappa_{f}$, which are generally expensive to obtain.

However, when choosing $|S|$ for the current iteration $k$, the stepsize ${\bf s}_{k}$ is yet to be determined. Based on the Lipschitz continuity of the involved functions, we argue that the previous stepsize is a fair estimator of the current one and this is confirmed by experimental results. Finally, we would like to point out that the sampling schemes derived in Eq. 34 gives our method a clear edge over sampling schemes that do not take any iteration information into account, e.g. linearly or geometrically increased samples.

2 Baselines and datasets

3 Results

The results in Figure 1 confirm our intuition that SCR can reduce ARCs computation time without losing its global convergence property. Newton’s method is the closest in terms of performance. However, it suffer heavily from an increase in $d$ as can be seen by additional results provided in the appendix. Furthermore, it cannot optimize the non-convex version of covtype due to a singular Hessian. Notably, BFGS terminates early on the non-convex higgs dataset due to a local saddle point. Finally, the high condition number of covtype has a significant effect on the performance of SGD, SAGA and L-BFGS.

Conclusion

In this paper we proposed a sub-sampling technique to estimate the gradient and Hessian in order to construct a cubic model analogue to trust region methods. We show that this method exhibits the same convergence properties as its deterministic counterpart, which are the best known worst-case convergence properties on non-convex functions. Our proposed method is especially interesting in the large scale regime when $n\gg d$. Numerical experiments on both real and synthetic datasets demonstrate the performance of the proposed algorithm which we compared to its deterministic variant as well as more classical optimization methods. As future work we would like to explore the adequacy of our method to train neural networks which are known to be hard to optimize due to the presence of saddle points.

References

Appendix A Appendix

First, we extend the Vector Bernstein inequality as it can be found in (Gross, 2011) to the average of independent, zero-mean vector-valued random variables.

Let ${\bf x}_{1},\ldots,{\bf x}_{n}$ be independent vector-valued random variables with common dimension $d$ and assume that each one is centered, uniformly bounded and also the variance is bounded above:

Then we have for $0<\epsilon<\sigma^{2}/\mu$

Proof: Theorem 6 in (Gross, 2011) gives the following Vector Bernstein inequality for independent, zero-mean vector-valued random variables

First, we shall define $\epsilon=t+\sqrt{V}$, which allows us to rewrite the above equation as

always holds, we can formulate a slightly weaker Vector Bernstein version as follows

Since the individual variance is assumed to be bounded above, we can write

Now, since $n>1$ and $\epsilon>0$, as well as $P({\bf z}>a)$ is falling in $a$ and $\exp(-x)$ falling in $x$, we can use this upper bound on the variance of ${\bf z}$ in (39), which gives the desired inequality

This result was applied in order to find the probabilistic bound on the deviation of the sub-sampled gradient from the full gradient as stated in Lemma 6, for which we will give the proof next.

To apply vector Bernstein’s inequality (35) we need to center the gradients. Thus we define

and note that from the Lipschitz continuity of $f$ (A3), we have

in equation (35), we can require the probability of a deviation larger or equal to $\epsilon$ to be lower than some $\delta\in(0,1]$

Conversely, the probability of a deviation of

Of course, any sampling scheme that guarantees the right hand side of (16) to be smaller or equal to $M$ times the squared step size, directly satisfies the sufficient gradient agreement condition (A5). Consequently, plugging the former into the latter and rearranging for the sample size gives Theorem 7 as we shall prove now.

A.1.2 Hessian Sampling

Let ${\bf A}_{1},..,{\bf A}_{n}$ be independent random Hermitian matrices with common dimension $d\times d$ and assume that each one is centered, uniformly bounded and also the variance is bounded above:

Proof: Theorem 12 in (Gross, 2011) gives the following Operator-Bernstein inequality

However, ${\bf Z}=\frac{1}{n}\sum_{i=1}^{n}{\bf A}_{i}$ and thus

Together with the Operator-Bernstein inequality, (52) and (53) give the desired inequality (49).

This result exhibits that sums of independent random matrices provide normal concentration near its mean in a range determined by the variance of the sum. We apply it in order to derive the bound on the deviation of the sub-sampled Hessian from the full Hessian as stated in Lemma 8, which we shall prove next.

Proof of Lemma 8: Bernstein’s Inequality holds as $f\in C^{2}$ and thus the Hessian is symmetric by Schwarz’s Theorem. Since the expectation of the random matrix needs to be zero, we center the individual Hessians,

and note that now from the Lipschitz continuity of ${\bf g}$ (A3):

Hence, for $\epsilon\leq 4\kappa_{g}$, we are in the small deviation regime of Bernstein’s bound with a sub-gaussian tail. Then, we may plug

Finally, we shall require the probability of a deviation of $\epsilon$ or higher to be lower than some $\delta\in(0,1]$

which is equivalent to ${}\left\|{\bf B}({\bf x})-{\bf H}({\bf x})\right\|$ staying within this particular choice of $\epsilon$ with probability $(1-\delta)$, generally perceived as high probability.

Proof of Theorem 9: Since $\|{\bf A}{\bf v}\|\leq\|{\bf A}\|_{op}\|{\bf v}\|\>\mbox{ for every }{\bf v}\in V$ we have for the choice of the spectral matrix norm and euclidean vector norm that any ${\bf B}$ that fulfils ${}\left\|({\bf B}({\bf x})-{\bf H}({\bf x}))\right\|\leq C{}\left\|{\bf s}\right\|$ also satisfies condition A4. Furthermore

Note that there may be a less restrictive sampling conditions that satisfy A4 since condition (56) is based on the worst case bound $\|{\bf A}{\bf v}\|\leq\|{\bf A}\|_{op}\|{\bf v}\|$ which indeed only holds with equality if ${\bf v}$ happens to be (exactly in the direction of) the largest eigenvector of $A$.

Finally, we shall state a Lemma which illustrates that the stepsize goes to zero as the algorithm converges. The proof can be found in Section 5 of (Cartis et al., 2011a).

Let $\{f({\bf x}_{k})\}$ be bounded below by some $f_{\inf}>-\infty$. Also, let ${\bf s}_{k}$ satisfy A1 and $\sigma_{k}$ be bounded below by some $\sigma_{\inf}>0$. Then we have for all successful iterations that

A.1.3 Illustration

In the top row of Figure 2 we illustrate the Hessian sample sizes that result when applying SCR with a practical version of Theorem 9 to the datasets used in our experiments see Section A.3 for details. In the bottom row of Figure 2, we benchmark our algorithm to the deterministic as well as two naive stochastic versions of ARC with linearly and exponentially increasing sample sizes.

Note that both the linear and the exponential sampling schemes do not quite reach the same performance as SCR even though they were carefully fine tuned to achieve the best possible performance. Furthermore, the sampling size was manually set to reach the full sample size at the very last iteration. This highlights another advantage of the automatic sampling scheme that does not require knowledge of the total number of iterations.

A.2 Convergence Analysis

The lower bound $\sigma_{\inf}$ follows directly from Step 7 in the algorithm design (see Algorithm 1). Within the upper bound, the constant $\sigma_{0}$ accounts for the start value of the penalty parameter. Now, we show that as soon as some $\sigma_{k}>3(\frac{2M+C+\kappa_{g}}{2})$, the iteration is very successful and $\sigma_{k+1}<\sigma_{k}$. Finally, $\gamma_{2}$ allows for $\sigma_{k}$ being ’close to’ the successful threshold, but increased ’one last time’.

Any iteration with $f({\bf x}_{k}+{\bf s}_{k})\leq m({\bf s}_{k})$ yields a $\rho_{k}\geq 1\geq\eta_{2}$ and is thus very successful. From a 2nd-order Taylor approximation of $f({\bf x}_{k}+{\bf s}_{k})$ around ${\bf x}_{k}$ we have:

Requiring the right hand side to be non-positive and solving for $\sigma_{k}$ gives the desired result.

Proof of Lemma 11 : By definition of the stochastic model $m_{k}({\bf s}_{k})$ we have

where we applied equation (11) first and equation (12) secondly.

Before proving the lower bound on the stepsize $\|{\bf s}_{k}\|$ we first transfer the rather technical result from Lemma 4.6 in (Cartis et al., 2011a) to our framework of stochastic gradients. For this purpose, let ${\bf e}_{k}$ be the gradient approximation error, i.e. ${\bf e}_{k}:={\bf g}_{k}-\nabla f({\bf x}_{k})$.

Let $f\in C^{2}$, Lipschitz continuous gradients (A3) and TC (A2) hold. Then, for each (very-)successful $k$, we have

with $\kappa_{\theta}\in(0,1)$ as in TC (13).

where the last inequality results from TC (Eq. (13)). Now, we can find the following bounds on the individual terms:

We can rewrite the right-hand side by a Taylor expansion of $\nabla f_{k+1}({\bf x}_{k}+{\bf s}_{k})$ around ${\bf x}_{k}$ to get

Contrary to the case of deterministic gradients, the first and third summand no longer cancel out. Applying the triangle inequality repeatedly, we thus get an error term in the final bound on (a):

(b) To bound the second summand, we can write

Finally, using the definition of $\theta_{k}$ as in (13) (which also gives $\theta_{k}\leq\kappa_{\theta}$ and $\theta_{k}\leq\kappa_{\theta}h_{k}$) and combining (a) and (b) we get the above result.

Proof of Lemma 12: The conditions of Lemma 21 are satisfied. By multiplying $d_{k}{}\left\|{\bf s}_{k}\right\|$ out in equation (60), we get

Now, applying the strong agreement conditions (A4) and (A5), as well as the Lipschitz continuity of H, we can rewrite this as

for all sufficiently large, successful $k$. Solving for the stepsize ${}\left\|{\bf s}_{k}\right\|$ give the desired result.

A.2.2 Local convergence

Before we can study the convergence rate of SCR in a locally convex neighbourhood of a local minimizer $w_{*}$ we first need to establish three crucial properties:

a lower bound on ${}\left\|{\bf s}_{k}\right\|$ that depends on ${}\left\|{\bf g}_{k}\right\|$.

an upper bound on ${}\left\|{\bf s}_{k}\right\|$ that depends on ${}\left\|{\bf g}_{k+1}\right\|$.

conditions under which all steps are eventually very successful.

With this at hand we will be able to relate ${}\left\|{\bf g}_{k+1}\right\|$ to ${}\left\|{\bf g}_{k}\right\|$, show that this ratio eventually goes to zero at a quadratic rate and conclude from a Taylor expansion around ${\bf g}_{k}$ that the iterates themselves converge as well.

Let ${\bf g}_{k}$ and ${\bf B}_{k}$ be sampled such that 17 and 19 hold in each iteration $k$. Furthermore, for unsuccessful iterations, assume that the sample size is not decreasing.

We have already established a lower stepsize bound in Lemma 12 so let us turn our attention directly to 2.:

Suppose that ${\bf s}_{k}$ satisfies (11) and that the Rayleigh coefficient

Proof: Given the above assumptions we can rewrite (11) as follows

where we used Cauchy-Schwarz inequality as well as the fact that $\sigma_{k}>0,\>\forall k$. Solving (70) for ${}\left\|{\bf s}_{k}\right\|$ gives (69).

Let $\{f({\bf x}_{k})\}$ be bounded below by some $f_{\inf}>-\infty$. Also, let A1, A3 hold and let ${\bf g}_{k}$ and ${\bf B}_{k}$ be sampled according to A22. Then we have w.h.p. that

The sampling schemes from Theorem 7 and Theorem 9 imply that the sufficient agreement assumptions A5 and A4 hold with high probability. Thus, we can deduce from Lemma 10 that after a certain number of consecutive unsuccessful iterates the penalty parameter is so high ($\sigma_{k}\geq\sigma_{sup}$) that we are guaranteed to find a successful step. Consequently, the number of successful iterations must be infinite ($|\mathcal{S}|=\infty$) when we consider the asymptotic convergence properties of SCR. We are left with two possible scenarios:

(i) If the number of unsuccessful iterations is finite ($|\mathcal{U}|\leq\infty$ & $|\mathcal{S}|=\infty$) we have that $\exists\;\hat{k}$ after which all iterates are successful, i.e. $k\in\mathcal{S},\forall\;k>\hat{k}$. From Lemma 20 we know that for all successful iterations ${}\left\|{\bf s}_{k}\right\|\rightarrow 0\text{ as }k\rightarrow\infty$. Consequently, due to the sampling scheme as specified in Theorem 7 and Theorem 9, $\exists\;\bar{k}\geq\hat{k}$ with $|S_{g,k}|=|S_{B,k}|=n,\;\forall\;k\geq\bar{k}$.

As a result the sample sizes eventually equal $n$ with high probability in all conceivable scenarios which proves the assertionWe shall see that, as a result of Lemma 25, the case of an infinite number of unsuccessful steps can actually not happen.

Now that we have (asymptotic) stepsize bounds and gradient (Hessian) agreement we are going to establish that, when converging, all SCR iterations are indeed very successful asymptotically.

Let $f\in C^{2}$, $\nabla f$ uniformly continuous and ${\bf B}_{k}$ bounded above. Let ${\bf B}_{k}$ and ${\bf g}_{k}$ be sampled according to A22, as well as ${\bf s}_{k}$ satisfy (11). Furthermore, let

with $\nabla f({\bf w}_{*})=0$ and ${\bf H}({\bf w}_{*})$ positive definite. Then there exists a constant $R_{min}>0$ such that for all $k$ sufficiently large

Furthermore, all iterations are eventually very successful w.h.p.

Proof: Since $f$ is continuous, the limit (72) implies that $\left\{f({\bf w}_{k})\right\}$ is bounded below. Since ${\bf H}({\bf w}_{*})$ is positive definite per assumption, so is ${\bf H}({\bf w}_{k})$ for all $k$ sufficiently large. Therefore, there exists a constant $R_{\min}$ such that

As a result of Lemma 24 we have that ${}\left\|{\bf e}_{k}\right\|\rightarrow 0$ as $k\rightarrow\infty$. Hence, Lemma 23 yields ${}\left\|{\bf s}_{k}\right\|\leq 1/R_{\min}{}\left\|\nabla f_{k}\right\|$ which implies that the step size converges to zero as we approximate $w^{*}$. Consequently, we are able to show that eventually all iterations are indeed very successful. Towards this end we need to ensure that the following quantity $r_{k}$ becomes negative for sufficiently large $k$:

where $\eta_{2}\in(0,1)$ is the ”very successful” threshold.

(i) By a (second-order) Taylor approximation around $f({\bf w}_{k})$ and applying the Cauchy-Schwarz inequality, we have:

where the term ${}\left\|{\bf e}_{k}\right\|{}\left\|{\bf s}_{k}\right\|$ is extra compared to the case of deterministic gradients.

(ii) Regarding the second part we note that if ${\bf s}_{k}$ satisfies (11), we have by the definition of $R_{k}$ and equation (73) that

which negated gives the desired bound on (ii). All together, the upper bound on $r_{k}$ is written as

Let us add and subtract ${\bf H}({\bf w}_{k})$ to the second summand and apply the triangle inequality

Now applying ${}\left\|{\bf A}{\bf v}\right\|\leq{}\left\|{\bf A}\right\|{}\left\|{\bf v}\right\|$ we get

We have already established in Lemma 24 that ${}\left\|{\bf e}_{k}\right\|\rightarrow 0$ and ${}\left\|({\bf H}_{k}-{\bf B}_{k})\right\|\rightarrow 0$. Together with Lemma 23 and the assumption ${}\left\|\nabla f_{k}\right\|\rightarrow 0$ this implies $\|{\bf s}_{k}\|\rightarrow 0$. Furthermore, since $\tau\in$ we have that ${}\left\|{\bf w}_{k}+\tau{\bf s}_{k}\right\|\leq{}\left\|{\bf w}_{k}+{\bf s}_{k}\right\|\leq{}\left\|{\bf s}_{k}\right\|$. Hence, ${\bf H}({\bf w}_{k}+\tau{\bf s}_{k})$ and ${\bf H}({\bf w}_{k})$ eventually agree. Finally, $\eta_{2}<1$ and $R_{\min}>0$ such that $r_{k}$ is negative for all $k$ sufficiently large, which implies that every such iteration is very successful.

From Lemma 10 we have $\sigma_{k}\leq\sigma_{sup}$. Furthermore, all assumptions needed for the step size bounds of Lemma 12 and 23 hold. Finally, Lemma 25 gives that all iterations are eventually successful. Thus, we can combine the upper (69) and lower (24) bound on the stepsize for all $k$ sufficiently large to obtain

which we can solve for the gradient norm ratio

Consequently, as long as the right hand side of (82) stays below infinity, i.e. ${}\left\|{\bf e}_{k}\right\|/{}\left\|\nabla f({\bf w}_{k})\right\|\not\rightarrow\infty$, we have quadratic convergence of the gradient norms. From Lemma 24 we have that ${}\left\|{\bf e}_{k}\right\|\rightarrow 0$ as $k\rightarrow\infty$ w.h.p. and furthermore $\kappa_{s}$ is bounded above by a constant and $R_{\min}$ is a positive constant itself which gives quadratic convergence of the gradient norm ratio with high probability. Finally, the convergence rate of the iterates follows from a Taylor expansion around ${\bf g}_{k}$.

A.2.3 First order global convergence

Note that the preliminary results Lemma 11 and 12 allow us to lower bound the function decrease of a successful step in terms of the full gradient $\nabla f_{k+1}$. Combined with Lemma 10, this enables us to give a deterministic global convergence guarantee while using only stochastic first order informationNote that this result can also be proven without Lipschitz continuity of $H$ and less strong agreement conditions as done in Corollary 2.6 in (Cartis et al., 2011a)..

We will consider two cases regarding the number of successful steps for this proof.

Case (i): SCR takes only finitely many successful steps. Hence, we have some index $k_{0}$ which yields the very last successful iteration and all further iterates stay at the same point ${\bf x}_{k_{0}+1}$. That is ${\bf x}_{k_{0}+1}={\bf x}_{k_{0}+i},\>\forall\>i\geq 1$. Let us assume that ${}\left\|\nabla f({\bf x}_{k_{0}+1})\right\|=\epsilon>0$, then

Since, furthermore, all iterations $k\geq k_{0}+1$ are unsuccessful $\sigma_{k}$ increases by $\gamma$, such that

However, this is in contradiction with Lemma 10, which states that $\sigma_{k}$ is bounded above. Hence, the above assumption cannot hold and we have ${}\left\|\nabla f({\bf x}_{k_{0}+1})\right\|={}\left\|\nabla f({\bf x}^{*})\right\|=0$.

Case (ii): sARC takes infinitely many successful steps. While unsuccessful steps keep $f({\bf x}_{k})$ constant, (very) successful steps strictly decrease $f({\bf x}_{k})$ and thus the sequence $\{f({\bf x}_{k})\}$ is monotonically decreasing. Furthermore, it is bounded below per assumption and thus the objective values converge

All requirements of Lemma 11 and Lemma 12 hold and we thus can use the sufficient function decrease equation (31) to write

Since $(f({\bf x}_{k})-f_{\inf})\rightarrow 0\text{ as }k\rightarrow\infty$, $\sigma_{\inf}>0,\eta_{1}>0$ and $\kappa_{s}^{3}\geq 0$ (as $\sigma_{\sup}<\infty$), we must have ${}\left\|\nabla f({\bf x}_{k})\right\|\rightarrow 0$, giving the result.

A.2.4 Second order global convergence and worst case iteration complexity

For the proofs of Theorem 15 and Theorem 17 we refer the reader to Theorem 5.4 in (Cartis et al., 2011a) and Corollary 5.3 in (Cartis et al., 2011b). Note that, as already laid out above in the proofs of Lemma 10 and Lemma 11, the constants involved in the convergence Theorems change due to the stochastic gradients used in our framework.

A.3 Details concerning experimental section

We here provide additional results and briefly describe the baseline algorithms used in the experiments as well as the choice of hyper-parameters. All experiments were run on a CPU with a 2.4 GHz nominal clock rate.

The real-world datasets we use represent very common instances of Machine Learning problems and are part of the libsvm library (Chang & Lin, 2011), except for cifar which is from Krizhevsky & Hinton (2009). A summary of their main characteristic can be found in table 1. The multiclass datasets are both instances of so-called image classification problems. The mnist images are greyscale and of size $28\times 28$. The original cifar images are $32\times 32\times 3$ but we converted them to greyscale so that the problem dimensionality is comparable to mnist. Both datasets have $10$ different classes, which multiplies the problem dimensionality of the multinomial regression by $10$.

Stochastic Gradient Descent (SGD): To bring in some variation, we select a mini-batch of the size $\lceil n/10\rceil$ on the real world classification- and $\lceil n/100\rceil$ on the multiclass problems. On the artificial datasets we only sample $1$ datapoint per iteration and update the parameters with respect to this point. We use a problem-dependent, constant step-size as this yields faster initial convergence (Hofmann et al., 2015),(Roux et al., 2012).

SAGA: is a variance-reduced variant of SGD that only samples 1 datapoint per iteration and uses a constant step-size.

Broyden-Fletcher-Goldfarb-Shanno (BFGS) is the most popular and stable Quasi-Newton method.

Limited-memory BFGS is a variant of BFGS which uses only the recent $K$ iterates and gradients to construct an approximate Hessian. We used $K=20$ in our experiments. Both methods employs a line-search technique that satisfies the strong Wolfe condition to select the step size.

NEWTON is the classic version of Newton’s method which we apply with a backtracking line search.

For L-BFGS and BFGS we used the implementation available in the optimization library of scipy. All other methods are our own implementation. The code for our implementation of SCR is publicly available on the authors’ webpage.

All of our experiments were started from the initial weight vector ${\bf w}_{0}:=(0,\ldots,0)$.

The regularization parameter updating is analog to the rule used in the reported experiments of (Cartis et al., 2011a), where $\gamma=2$. Its goal is to reduce the penalty rapidly as soon as convergence sets in, while keeping some regularization in the non asymptotic regime. A more sophisticated approach can be found in (Gould et al., 2012). In our experiments we start with $\sigma_{0}=1,\eta_{1}=0.2,\text{ and }\eta_{2}=0.8$ as well as an initial sample size of $5\%$.

To test the influence of the dimensionality on the progress of the above applied methods we created artificial datasets of three different sizes, labeled as gaussian s, gaussian m and gaussian l.

with a mean of zero $\mu=(0,\ldots,0)$ and a covariance matrix that has reasonably uniformly distributed off-diagonal elements in the interval $(-1,1)$.

As expected, the classic Newton methods suffers heavily from an increase in the dimension. The regularized Newton methods on the other hand scale comparably very well since they only need indirect access to the Hessian via matrix-vector products. Evidently, these methods outperform the quasi-newton approaches even in high dimensions. Among these, the limited memory version of BFGS is significantly faster than its original variant.

In this section we leave the trust region method out because our implementation is not optimized towards solving multi-class problems. We do not run Newton’s method or BFGS either as the above results suggests that they are unlikely to be competitive. Furthermore, Figure 5 does not show logarithmic but linear suboptimality because optimizing these problems to high precision takes very long and yields few additional benefits. For example, the 25th SCR iteration drove the gradient norm from $3.8\cdot 10^{-5}$ to $5.6\cdot 10^{-8}$ after building up a Krylov space of dimensionality $7800$. It took 9.47 hours and did not change any of the first $13$ digits of the loss. As can be seen, SCR provides early progress at a comparable rate to other methods but gives the opportunity to solve the problem to high precision if needed.