Adaptive Stochastic Variance Reduction for Subsampled Newton Method with Cubic Regularization

Keywords:

cubic-regularized Newton method, subsampling, stochastic variance reduction, randomized algorithm, iteration complexity, sample complexity.

Introduction

We consider the problem of minimizing the average of a large number of loss functions:

Our goal is to find an approximate local solution $x$ that satisfies

For minimizing a general smooth and nonconvex function $F$, Nesterov and Polyak introduced a modified Newton method with cubic regularization (CR). Specifically, each iteration of the CR method consists of the following updates:

Assuming the Hessian $\nabla^{2}F$ to be Lipschitz continuous, it is shown in that the CR method finds an approximate solution satisfying (2) within $\mathcal{O}(\epsilon^{-3/2})$ iterations. This is better than purely gradient-based methods, which need $\mathcal{O}(\epsilon^{-2})$ iterations to reach a point $x$ satisfying $\|\nabla F(x)\|\leq\epsilon$ [12, Section 1.2.3]. However, the computational cost per iteration of CR can be much higher than gradient-based methods.

Much recent efforts have been devoted to improving the efficiency of CR by exploiting the finite-sum structure in (1); see, e.g., . An natural approach is to replace $\nabla F(x^{k})$ and $\nabla^{2}F(x^{k})$ by subsampled approximations:

where $\mathcal{S}_{k},\mathcal{B}_{k}\subseteq\{1,\ldots,N\}$ are two sets (or multisets for sampling with replacement) of random indices at the $k$th iteration. The cost of computing the Hessians $\nabla^{2}f_{i}$ usually dominates that of the gradients. Moreover, the cost of solving the CR subproblem (3) may grow fast when the batch size $|\mathcal{B}_{k}|$ increases, especially when using iterative methods such as gradient descent or the Lanczos method . Therefore, an important measure of efficiency is the number of second-order oracle calls for $\nabla^{2}f_{i}$, i.e., the Hessian sample complexity.

In this paper, we develop an adaptive subsampling CR method that requires $\mathcal{O}(N+N^{2/3}\epsilon^{-3/2})$ second-order oracle calls in expectation. Assuming that $\epsilon$ is small enough, we often simply refer to it as $\mathcal{O}(N^{2/3}\epsilon^{-3/2})$. Notice that using the choices in (4) would require $\mathcal{O}(N\epsilon^{-3/2})$ Hessian samples. Thus this is a significant improvement especially when $N$ is very large. In the rest of this section, we discuss several related work, and then outline our contributions.

It is shown in that the order of convergence rate of the CR method remains the same as long as $g_{k}$ and $H_{k}$ in (3) satisfy

with $\xi_{k}$ being defined in (3) and $C_{1},C_{2}$ being some positive constants. Here $\|\cdot\|$ for vectors denotes their Euclidean norm, and for matrices denotes their spectral norm. In order to exploit the averaging structure in (1), a subsampled cubic regularization (SCR) method was proposed in where $g^{k}$ and $H^{k}$ are calculated as in (5) and (6) (here we omit additional accept/reject steps based on trust-region methods). Matrix concentration inequalities are used in to derive appropriate sample sizes $|\mathcal{S}_{k}|$ and $|\mathcal{B}_{k}|$ such that the conditions (7) and (8) hold with high probability. In particular, matrix Bernstein inequality (e.g., ) implies that with probability at least $1-\delta$,

where $L$ is a uniform Lipschitz constant of $\nabla f_{i}$ for all $i$. Therefore, if we upper bound the right-hand side above by $C_{2}\|\xi_{k}\|$, then (8) holds with probability at least $1-\delta$ provided that

A similar condition on $|\mathcal{S}_{k}|$ is also derived in . The overall gradient and Hessian sample complexities for SCR, with or without replacement, are summarized in Table 1.1 (based on the analysis in ). When $\epsilon\leq 1/N$, SCR can be much worse than the deterministic CR method.

In order to further reduce the Hessian sample complexity, several recent works combine CR with stochastic variance reduction techniques. Stochastic variance-reduced gradient (SVRG) method was first proposed to reduce gradient sample complexity of randomized first-order algorithms (see for convex optimization and for nonconvex optimization). Two different SVRC (stochastic variance-reduced cubic regularization) methods were proposed in and respectively. Both of them employed the same variance-reduction technique to reduce the Hessian sample complexity. In particular, incorporated additional second-order corrections in gradient variance-reduction, therefore it obtained better gradient sample complexity, but with slightly worse Hessian sample complexity than . See Table 1.1 for a summary of their results. All these methods require $\mathcal{O}(\epsilon^{-3/2})$ calls to solve the cubic regularized sub-problem (3).

Subsampled Newton methods without CR have been studied in, e.g., , but their convergence rates are worse than the ones that are based on CR.

Among the works on CR with stochastic variance reduction, most of them (e.g., ) rely on matrix concentration bounds such as (10) to set the sample sizes $|\mathcal{B}_{k}|$ and $|\mathcal{S}_{k}|$ according to $\|\xi^{k}\|$ at each iteration $k$, The problem is that $\xi^{k}$ is the solution to the CR sub-problem in (3), where $g^{k}$ and $H^{k}$ need to be obtained from $|\mathcal{S}_{k}|$ and $|\mathcal{B}_{k}|$ samples in the first place. While one can assume bounds like (10) hold in the complexity analysis, they do not provide practically implementable variance-reduction schemes.

In contrast, the sample sizes in are set as constants across all iterations, which only depend on $N$ and $\log(d)$. The disadvantage of this approach is that it can be very conservative. According to concentration bounds such as (10), the sample sizes at the initial stage of the algorithm (when $\|\xi^{k}\|$ is large) can be set much smaller than the ones required in the later stage (when $\|\xi^{k}\|$ is very small). Therefore when using a constant sample size, much of the samples in the initial stage can be wasteful.

There is another technicality of using matrix concentration bounds: we need inequalities such as (9) to hold for all iterations with high probability, say with probability at least $1-\delta_{0}$. Then the probability margin $\delta$ per iteration and $\delta_{0}$ should satisfy $(1-\delta)^{T}\geq 1-\delta_{0}$, where $T=\mathcal{O}(\epsilon^{-3/2})$ is the number of iterations for CR methods. Therefore, $\delta=\mathcal{O}(\delta_{0}\epsilon^{3/2})$ and there is an additional $\mathcal{O}\left(\log(d/\epsilon\delta_{0})\right)$ factor in the sampling complexities (see Table 1.1). The results in are for convergence in expectation, nevertheless they still have a $\log(d)$ factor and unusually large constants in their bounds.

2 Contributions and outline

In this paper, we develop an adaptive sampling scheme for stochastic variance reduction in the subsampled Newton method with cubic regularization. In particular, the gradient sample size $|\mathcal{S}_{k}|$ and Hessian sample size $|\mathcal{B}_{k}|$ are chosen adaptively in each iteration to ensure the following conditions hold in expectation (conditioned on $x^{k}$):

where $C^{\prime}_{1}$ and $C^{\prime}_{2}$ are some positive constants. The major difference from (7) and (8) is that here $\|\xi^{k-1}\|$ is a known quantity conditioned on $x^{k}$, before choosing the sample sizes to form the approximations $g^{k}$ and $H^{k}$. Indeed we choose $|\mathcal{S}_{k}|$ and $|\mathcal{B}_{k}|$ based on $\|\xi^{k-1}\|$. Such an adaptive scheme is readily implementable in practiceRight before submitting this paper, we discovered a recent note which independently showed that the conditions (7) and (8) are sufficient to retain the same convergence rate of the exact CR method. However it does not provide improved sample complexity over previous work listed in Table 1.1.. Moreover, it does not waste samples in the early stage of the algorithm as constant sample sizes do.

We show that our adaptive subsampled CR method has an expected iteration complexity $\mathcal{O}(\epsilon^{-3/2})$, which is the number of times the CR subproblem in (3) needs to be solved in order to find a point $x$ satisfying (2), which is the same as that of the deterministic CR method. However, the total Hessian sample complexity of our method is $\mathcal{O}(N^{2/3}\epsilon^{-3/2})$, which is much better than $\mathcal{O}(N\epsilon^{-3/2})$ of the full CR method, and is indeed better than all previous works listed in Table 1.1.

In addition to the improved Hessian sample complexity, the techniques in our analysis are quite different from those adopted in previous works. In particular, we avoid using any high probability bounds based on matrix or vector concentration inequalities. Instead, our analysis is based on novel bounds on the 3rd and 4th order moments of the average of independent random matrices, which are of independent interest on their own. The type of convergence studied in is also in expectation. However, their analysis still relies on some matrix concentration inequalities, thus their results contain the $\log(d)$ factor and excessively large constants.

The rest of this paper is organized as follows. In Section 2, we present the adaptive SVRC method and its convergence analysis. We show that it retains the ${\mathcal{O}}(\epsilon^{-3/2})$ iteration complexity of the exact cubic regularization method, but with only ${\mathcal{O}}(N^{2/3})$ Hessian samples per iteration. In addition, we show that sampling with or without replacement have the same order of sample complexity. In Section 3, we study a non-adaptive SVRC method with fixed sample size at each iteration. We show that if exact gradients are available, then it attains the same total Hessian sample complexity ${\mathcal{O}}(N^{2/3}\epsilon^{-3/2})$. If both gradient and Hessian need to be subsampled, we examine the SVRC method of using the higher moments bounds developed in this paper and obtain refined analysis.

Adaptive variance reduction for cubic regularization

In this section we first present the adaptive SVRC method, then analyze its convergence rate and Hessian sample complexity.

where $\epsilon_{g}$ and $\epsilon_{H}$ are determined by $\xi^{k}_{t-1}$ and the desired tolerance $\epsilon$:

Then we compute the subsampled gradient and Hessian as

This construction follows the stochastic variance reduction scheme proposed in , which has been adopted by recent works on subsampled Newton method with cubic regularization .

We make several remarks regarding the subsampling scheme in Algorithm 1. First, compared with the algorithms in , which use large enough sample sizes to ensure (7) and (8) with high probability, Algorithm 1 aims to ensure (11) and (12) in expectation. Moreover, instead of depending on $\xi^{k}_{t}$ which is not available until after the current iteration, our sample sizes are determined by $\xi^{k}_{t-1}$, which is computed in the previous iteration. Second, compared with the constant sample size used in , our adaptive sampling scheme may use much less samples in the early stages when $\|\xi^{k}_{t-1}\|$ is relatively large. Our overall Hessian sample complexity is better than either of these two previous approaches.

An alternative approach to avoid the dependence of sample sizes on $\xi^{k}_{t}$ is to use the full gradient, i.e., let $g_{t}^{k}=\nabla F(x_{t}^{k})$, and use (15) or (6) to compute the Hessian approximation . Then condition (7) holds automatically, and one can replace (8) with

because it can be be shown that $\|\xi^{k}_{t}\|\leq c\|\nabla F(x^{k}_{t})\|$ for some large enough constant $c$. Since the full gradient $\nabla F(x^{k}_{t})$ can be computed before sampling $\mathcal{B}^{k}_{t}$, this condition can be used to determine sample size $|\mathcal{B}^{k}_{t}|$. However, it can be shown that one have roughly $\|\nabla F(x^{k}_{t})\|\leq{\mathcal{O}}(\|\xi^{k}_{t-1}\|^{2})$ plus some random noise, thus the mini-batch size required for (16) to hold can be much larger than that for (12).

We make the following assumption regarding the objective function in (1):

Consequently, $\nabla F$ and $\nabla^{2}F$ are $L$ and $\rho$-Lipschitz continuous respectively.

This assumption is very similar to those adopted in subsampled Newton methods with cubic regularization . The only difference is that here we use the Frobenius norm, instead of the spectral norm, for the Hessian smoothness assumption. The advantage of using the Frobenius norm is that it works well with matrix inner product, which allows us to derive simple bounds on the 3rd and 4th order moments for the average of random matrices. On the other hand, the Lipschitz constant $\rho$ is always larger than the one corresponding to the spectral norm, up to a factor of $\sqrt{d}$ in the worst case. However, when the Hessians $\nabla^{2}f_{i}$ are of low-rank or ill-conditioned, which is often the case in practice, the Lipschitz constants for different norms are very close.

Before presenting the main results, we first provide a few supporting lemmas. The first one gives a simple bound on the 4th order moment of the average of i.i.d. random matrices with zero mean. Its proof is given in Appendix A.

Based on the above lemma, we can bound the 2nd and 4th order variances of the gradient and Hessian approximations. The following lemma is proved in Appendix B.

Let the variance reduced gradient $g^{k}_{t}$ and Hessian $H^{k}_{t}$ be constructed according to (14) and (15). Then they satisfy the following equalities and inequalities

As a result, we have the following corollary, whose proof is given in Appendix C.

Let $H^{k}_{t}$, $g^{k}_{t}$, $\xi^{k}_{t}$ and the mini-batch index sets $\mathcal{B}^{k}_{t}$ and ${\mathcal{S}}^{k}_{t}$ be generated according to Algorithm 1. Then we have

Now we are ready to present the descent property of Algorithm 1.

Suppose the sequence $\{x^{k}_{t}\}^{k=1,...,K}_{i=1,...,m}$ is generated by Algorithm 1. Then the following descent property holds

As a remark, as long as $\frac{\sigma}{4}-\frac{\rho}{2}-\frac{L}{3}>\frac{5\rho}{2}+\frac{2L}{3}$, the variance term associated with $\|\xi^{k}_{t-1}\|^{3}$ in current step can be dominated by the descent associated with $\|\xi^{k}_{t-1}\|^{3}$ in the previous step. At the first step of each stage, the variance term is 0, hence we define $\xi^{k}_{-1}=0$ for a unified expression.

Consider the cubic regularization subproblem in Algorithm 1:

Then by the Lipschitz continuous condition of the objective function, we have

where the second line is due to (21) and the fifth line is due to (22). By the following variant of Young’s inequality

Since $g^{k}_{0}=\nabla F(x^{k}_{0})$ and $H^{k}_{0}=\nabla^{2}F(x^{k}_{0})$ at the beginning of each stage, we have

By substituting the variance bounds in Corollary 2.4 into (24), we get

Let $x^{k}_{t+1}$, $\xi^{k}_{t-1}$ and $\xi^{k}_{t}$ be generated by Algorithm 1. Then the following relations hold

Finally we are ready to present the main result on iteration complexity.

Choose the parameter $\sigma>13\rho+4L$ and let $k^{*},i^{*}$ be given by either output option in Algorithm 1 after running for $K$ stages, then

where $F^{\star}=\min_{x}F(x)$. As a result,

If we choose $K$ and $m$ such that Km={\mathcal{O}}\bigl{(}\epsilon^{-3/2}\bigr{)}, then within {\mathcal{O}}\bigl{(}\epsilon^{-3/2}\bigr{)} iterations, Algorithm 1 would reach a point $x^{k^{*}}_{t^{*}+1}$ such that

In other words, the approximate optimality conditions in (2) are satisfied in expectation.

First, taking expectation over the whole history of random samples for (20) and (25) and summing them up for $t=0,...,m-1$ yield

Under the assumption $\sigma>13\rho+4L$, we have $\sigma/4-3\rho-L>0$. Further summing over the stages $k=1,\ldots,K$, we obtain

For option 1 in the output rule, due to the concavity of $\min$ function, Jensen’s inequality gives

which is precisely (26). For option 2, since $k^{*}$ and $t^{*}$ are randomly chosen,

Therefore, for both options, inequality (26) holds.

Next, we derive guarantees for approximating the first and second-order stationary conditions. According to Lemma 2.6,

where the first inequality is due to Hölder’s inequality, the second inequality is due to Jensen’s inequality, the third inequality is due to (26) and the last inequality is due to the fact that $(a+b)^{\theta}\leq a^{\theta}+b^{\theta}$ for all $a,b>0$ and $0\leq\theta\leq 1$. Finally, substituting the above inequality into (30) yields the desired result in (27).

In order to bound the minimum eigenvalue of the Hessian, we have from Lemma 2.6,

Similar to the arguments used for proving the first-order bound,

Combining the inequality above with (31) gives the desired bound in (28). ∎

2 Bounding the Hessian sample complexity

Due to the adaptive mini-batch size rule, the total Hessian sample complexity is not given explicitly. In this subsection, we provide a bound on the complexity of Hessian sampling.

Let the total number of Hessian samples in Algorithm 1 be $B_{H}$. If we set the length of each stage and number of stages to be

then the expectation of $B_{H}$ taken to reach a second-order $\epsilon$-solution will be

According to (13), it suffices to use the following sample size for approximating the Hessian,

where the first inequality is due to the triangle inequality and the second one is due to the Cauchy-Schwarz inequality. Summing up for all $k=1,\ldots,K$ and $t=0,\ldots,m$, we get

where we used $Km=\epsilon^{-3/2}$ as specified in (32). Taking expectation on both sides gives

where the second inequality is due to Jensen’s inequality and the third inequality is due to (29).

where we used $Km=\epsilon^{-3/2}$ and the constant $C$ is defined as

Note that the above bound holds only when $\epsilon$ is small enough so that (32) makes sense. In order to cover the case when $\epsilon$ is large, we can write

Through similar arguments, one can bound the total gradient sample complexity $B_{G}$ by

3 Analysis of sampling without replacement

In this section, we show that sampling without replacement will not change the sample complexity of Algorithm 1. Different Hessian sample complexity bounds for sampling with and without replacement are derived in (see Table 1.1), and both of them are worse than the $\mathcal{O}(N^{2/3}\epsilon^{-2/3})$ complexity obtained in this paper. Again, we start from the variance estimation of the subsampled gradient and Hessian.

The proof of this lemma is given in Appendix E. As a result of this lemma, we have the following corollary. The proof is similar to that of Lemma 2.3, which we omit for simplicity.

Let $g^{k}_{t}$ and $H^{k}_{t}$ be constructed by (14) and (15) respectively, with the mini-batches ${\mathcal{S}}^{k}_{t}$ and $\mathcal{B}^{k}_{t}$ be sampled without replacement. Then

Thus by Corollary 2.11, the following inequality holds approximately

We can derive the desired mini-batch sizes for Hessian and gradient smapling by requiring

For the Hessian sample complexity, if we want to have a sublinear sample size $N^{\alpha}$ with $\alpha<1$, then we should expect $|\mathcal{B}^{k}_{t}|\leq\mathcal{O}(N^{\alpha})$ holds during the iterations. This requires

Consider Algorithm 1, where we sample ${\mathcal{S}}^{k}_{t}$ and $\mathcal{B}^{k}_{t}$ without replacement and set the length of each epoch and number of epochs to be $m=\mathcal{O}(N^{1/3})$ and $K=\epsilon^{-3/2}/m$ respectively. Then the totoal number of Hessian samples $B_{H}$ and total number of gradient samples $B_{G}$ required to reach a second-order $\epsilon$-solution is

Analysis of non-adaptive SVRC Schemes

In this section, we consider variants of the SVRC methods that use fixed gradient and Hessian sample sizes across all iterations. In the first variant, we use the full gradients but subsampled Hessians. In this case, we show that the total Hessian sample complexity is still $\mathcal{O}(N^{2/3}\epsilon^{-3/2})$. The second variant uses both approximate gradients and approximate Hessians, and adds a correction term to the gradient approximation based on the second-order information. This variant is proposed in , where an $\widetilde{\mathcal{O}}(N^{4/5}\epsilon^{-3/2})$ sample complexity is proved for both the gradient and Hessian approximations. We obtain the same order of sample complexity using the higher moment bounds developed in this paper (instead of using concentration inequalities), which avoid the $poly(\log d)$ factor and excessively large constant in the results of .

Algorithm 2 describes the SVRC method using full gradient in each iteration, i.e., $g^{k}_{t}=\nabla F(x^{k}_{t})$. One remark regarding output option 1 is that the best choice of $\gamma$, depending on the parameters $B,m,\sigma$ and $\rho$, is given in Lemma 3.3. However, one does not need to know the value exactly, and any choice of $\gamma=\Theta(\rho)$ will not affect our complexity result.

Similar to the derivation of (23), we have the following result,

Note that the expectation is not taken over $|\mathcal{B}^{k}_{t}|$ since it is a predetermined constant, which we denote as $B$ from now on. Let us define a Lyapunov function

where the coefficients $c_{t}$ are constructed recursively by setting $c_{m}=0$ and

Here $\theta_{1},\theta_{2}$ are some constant to be determined later. Next, we prove a monotone decreasing property of this Lyapunov function over one epoch of Algorithm 2.

First, we note the following simple fact, which is proved in Appendix F.

For any $a,b,\theta_{1},\theta_{2}>0$, the following inequality holds:

Now, using the definition of $R^{k}_{i}$ in (37) and expression of $c_{i}$ (38), we see that the above inequality leads to the following descent property of the Lyapunov function.

Let the sequence $\{x^{k}_{t}\}$ be generated by Algorithm 2, then for all $k$ and $0\leq t\leq m-1$,

Next, we prove that when the parameters are properly choosen, then $\gamma={\mathcal{O}}(\rho)$.

Suppose we set the batch size $B=\alpha N^{2/3}$, where $\alpha\geq 8$ is a constant. If we set $\sigma\geq 3\rho$, $m=(1/3)N^{1/3}$, $\theta_{1}=N^{1/9}$ and $\theta_{2}=N^{1/18}$, then

By (38) and the values of $\theta_{1}$, $\theta_{2}$, $B$ and $m$, we have

Adding $\frac{\rho}{\alpha^{3/2}N^{2/3}}$ to both sides of the above equality, we obtain

where we used $c_{m}=0$ and by choosing $m={(1/3)N^{1/3}}$ we have $3N^{-1/3}=1/m$ and $(1+1/m)^{m}\leq e$. In addition,

Therefore if $\alpha\geq 8>(8e)^{2/3}$, then we have $\gamma={\mathcal{O}}(\rho)$. ∎

Suppose in Algorithm 2 we set $m=(1/3)N^{1/3}$, $|\mathcal{B}^{k}_{t}|=B=8N^{2/3}$, $\sigma\geq 3\rho$, and $\gamma$ is defined in (42). Let $k^{*},t^{*}$ be given by either of the two output options in the algorithm, then

and the expected total Hessian samples required to reach a second-order $\epsilon$-solution is ${\mathcal{O}}(N^{2/3}\epsilon^{-3/2})$.

Then Lemma 3.2 and Lemma 3.3 immediately give

Using similar arguments for the output option 1 and 2 in the proof of Theorem 2.7, the first bound (43) follows.

From (43), by Jensen’s inequality, one simply gets

Then following the proof of Lemma 2.6 in Appendix D, more specifically (53), we have

Similarly, using (53), one can bound the minimum eigenvalue of the Hessian as

Finally, by setting $B=8N^{2/3}$, $m=(1/3)N^{1/3}$ and $K=\epsilon^{-3/2}/m$, the total number of Hessian samples is

2 The case of using subsampled gradient

In this subsection, we analyze a non-adaptive SVRC scheme proposed in , which is shown in Algorithm 3. Similar to Algorithm 2, one need not knwon the best value of $\gamma$ given in Lemma 3.7, any choice of $\gamma=\Theta(\rho)$ will be acceptable. In this scheme, the approximate gradients and Hessians are constructed as follows:

The following lemma is proved in Appendix G.

Let $g^{k}_{t}$ be constructed by (44) and $|{\mathcal{S}}^{k}_{t}|=S$, then the following inequalities hold:

By the above lemma and a discussion similar to that of (23) and (36), we have

where the last inequality is due to Lemmas 2.3 and 3.5 and Jensen’s inequality. In the rest of the analysis, we use the Lyapunov function $R_{t}^{k}$ defined in (37) but with a new set of parameters by setting $c_{m}=0$ and

Again the constants $\theta_{1}$ and $\theta_{2}$ will be determined later. We present the following results without proof due to the similarity to their counterparts in previous subsection.

Suppose the sequence $\{x^{k}_{t}\}_{t=1,...,m}^{k=1,...,K}$ is generated by Algorithm 3, then we have

Define \gamma=\displaystyle\min_{1\leq t\leq m}\textstyle\left\{\frac{\sigma}{4}-\frac{5\rho}{6}-c_{t+1}\bigl{(}1+\theta_{1}^{6}+2\theta_{2}^{3}\bigr{)}\right\}. Then for $k=1,...,K$,

If we set $\sigma\geq 4\rho$, $m=(1/3)N^{1/5}$, $B=\alpha N^{2/5}$, $S=\alpha^{2}N^{4/5}$, $\theta_{1}=N^{1/15}$ and $\theta_{2}=N^{1/30}$, then we have $\rho>0$ as long as $\alpha\geq 12$ and

Let $k^{*},t^{*}$ be chosen by either of the two output options in the algorithm, then

and the total Hessian and gradient sample complexity for finding a second-order $\epsilon$-solution is {\mathcal{O}}\bigl{(}N^{4/5}\epsilon^{-3/2}\bigr{)}.

Discussion

We considered the problem of minimizing the average of a large number of smooth and possibly nonconvex functions, $F(x)=(1/N)\sum_{i=1}^{N}f_{i}(x)$, using subsampled Newton method with cubic regularization. We presented an adaptive variance reduction method that requires ${\mathcal{O}}(N+N^{2/3}\epsilon^{-3/2})$ Hessian samples for finding an approximate solution satisfying $\|\nabla F(x)\|\leq\epsilon$ and $\nabla^{2}F(x)\succeq-\sqrt{\epsilon}I$. This result holds for both sampling with and without replacement. Our analysis do not rely on high probability bounds from matrix concentration inequalities, instead, we use bounds on 3rd and 4th moments of the average of random matrices.

We have focused on the Hessian sample complexity by assuming that the solution to the cubic regularization (CR) subproblem (3) is available at each iteration. Nesterov and Polyak showed that the CR subproblem is equivalent to a convex one-dimensional optimization problem, but this approach requires eigenvalue decomposion of the approximate Hessian $H^{k}$, which can be very costly for large-dimensional problems. Several recent works propose to solve the CR subproblem using iterative algorithms such as gradient descent or Lanczos method , and approximate trust-region solver can also be used. However, the overall complexity of the combined methods may still be high.

An efficient approximate solver for the CR subproblem has been developed in , which leads to a total computational complexity of ${\mathcal{O}}(N\epsilon^{-3/2}+N^{3/4}\epsilon^{-7/4})$ for minimizing the finite average problem (1). This complexity is measured in terms of the total number of Hessian-vector products, i.e., multiplications by the component Hessians $\nabla^{2}f_{i}(x^{k})$. Similar results have also been obtained by . For many mahcine learning problems, including generalized linear models and training neural networks, such Hessian-vector products can be computed in ${\mathcal{O}}(d)$ time . We notice that the CR subproblem solved in includes all $N$ component Hessians. Thus it is possible to further reduce the overall computational complexity by combining the efficient CR solver developed in and the lower Hessian sample complexity obtained in this paper.

Acknowledgments

We thank Zeyuan Allen-Zhu for helpful discussions on the complexities of several first and second-order methods for nonconvex optimization and the approximate cubic regularization solver in .

Appendix A Proof of Lemma 2.2

Appendix B Proof of Lemma 2.3

where $j$ is a uniform sample from $\{1,...,N\}$. Note that

Similarly, one can get the variance bounds for the gradient estimate. This part of the proof is omitted for simplicity. ∎

Appendix C Proof of Corollary 2.4

By Lemma 2.3 and the mini-batch size rule in Algorithm 1,

Due to the concavity of the square root function $\sqrt{\cdot}$, we can apply Jensen’s inequality to obtain

Again, with the concavity of the function $(\cdot)^{3/4}$, applying Jensen’s inequality yields

Following a similar line of arguments, one can get the bounds for the gradient variances. We omit the details for simplicity. ∎

Appendix D Proof of Lemma 2.6

By the Lipschitz continuity of $\nabla^{2}F$ and the fact $x^{k}_{t+1}=x^{k}_{t}+\xi^{k}_{t}$,

In addition, the optimality condition (21) implies

Taking expectation on both sides of the above inequality and applying Corollary 2.4 yield

For the secnod inequality, the optimality condition (22) implies

Therefore, according to the Lipschitz continuity of $\nabla^{2}F$,

Taking expectation on both sides and applying the Corollary 2.4 results in

which is the second desired inequality. ∎

Appendix E Proof of Lemma 2.10

Equation (34) is a standard result for variance analysis of sampling without replacement scheme. Hence we omit the proof of this equation. To prove (35), we start with the expansions in (51) and calculate the terms $T_{1},...,T_{7}$ for sampling without replacement. First, the following three terms do not change:

where in the last equality we used result for $T_{2}$. In similar ways, one can find the expressions for $T_{6}$ and $T_{7}$:

Summing these terms up gives the desired equation (35). ∎

Appendix F Proof of Lemma 3.1

We expand $(a+b)^{3}$ and then use Young’s inequality,

Appendix G Proof of Lemma 3.5

where the last inequality is due to the Lipschitz continuity of $\nabla^{2}f$. ∎

References