Title: A STOCHASTIC TRUST REGION METHOD FOR NON-CONVEX MINIMIZATION

Abstract: We target the problem of finding a local minimum in non-convex finite-sum minimization. Towards this goal, we first prove that the trust region method with inexact gradient and Hessian estimation can achieve a convergence rate of order O(1/k 2/3 ) as long as those differential estimations are sufficiently accurate. Combining such result with a novel Hessian estimator, we propose a sample-efficient stochastic trust region (STR) algorithm which finds an ( , √ )-approximate local minimum within Õ( √ n/ 1.5 ) stochastic Hessian oracle queries. This improves the state-of-the-art result by a factor of O(n 1/6 ). Finally, we also develop Hessian-free STR algorithms which achieve the lowest runtime complexity. Experiments verify theoretical conclusions and the efficiency of the proposed algorithms.

1. INTRODUCTION:

We consider the following finite-sum non-convex minimization problem

EQUATION

where each (non-convex) component function f i : R d → R is assumed to have L 1 -Lipschitz continuous gradient and L 2 -Lipschitz continuous Hessian. Since first-order stationary points could be saddle points with inferior generalization performance (Dauphin et al., 2014) , in this work we are particularly interested in computing ( , √ )-approximate second-order stationary points, -SOSP:

EQUATION

To find a local minimum of problem (1) , the cubic regularization approach (Nesterov & Polyak, 2006) and the trust region algorithm (Conn et al., 2000; Curtis et al., 2017) are two classical methods. Specifically, cubic regularization forms a cubic surrogate function for the objective F (x) by adding a third-order regularization term to the second-order Taylor expansion, and minimizes it iteratively. Such a method is proved to achieve an O(1/k 2/3 ) global convergence rate and thus needs O(n/ 1.5 ) stochastic first-and second-order oracle queries, namely the evaluation number of stochastic gradient and Hessian, to achieve a point that satisfies (2). On the other hand, trust region algorithms estimate the objective with its second-order Taylor expansion but minimize it only within a local region.

Recently, Curtis et al. (2017) proposes a trust region variant to achieve the same convergence rate as the cubic regularization approach. But both methods require computing full gradients and Hessians of F (x) and thus suffer from high computational cost in large-scale problems.

To avoid costly exact differential evaluations, many works explore the finite-sum structure of problem (1) and develop stochastic cubic regularization approaches. Both Kohler & Lucchi (2017b) and Xu et al. (2017) propose to directly subsample the gradient and Hessian in the cubic surrogate function, and achieve Õ(1/ 3.5 ) and Õ(1/ 2.5 ) stochastic first-and second-order oracle complexities respectively. By plugging a stochastic variance reduced estimator (Johnson & Zhang, 2013) and the Hessian tracking technique (Gower et al., 2018) into the gradient and Hessian estimation, the approach in (Zhou et al., 2018a) improves both the stochastic first-and second-order oracle complexities to Õ(n 0.8 / 1.5 ). Recently, Zhang et al. (2018) and Zhou et al. (2018b) develop more efficient stochastic cubic regularization variants, which further reduce the stochastic second-order oracle complexity to Õ(n 2/3 / 1.5 ) at the cost of increasing the stochastic first-order oracle complexity to Õ(n 2/3 / 2.5 ).

Table 1 : Stochastic first-and second-order oracle complexities, SFO and SSO for short respectively, of the proposed STR approaches and other state-of-the-art methods. When SSO is prioritized, our STR 1 has strictly better complexity than both SCR and Lite-SVRC. When SFO and SSO are treated equally, STR 2 improves the existing result in SVRC.

Algorithm SFO SSO TR (Conn et al., 2000; Curtis et al., 2017) O( n 1.5 ) O( n 1.5 ) CR (Nesterov & Polyak, 2006) O( n 1.5 ) O( n 1.5 ) SCR (Kohler & Lucchi, 2017a) Õ( 1 3.5 ) Õ( 1 2.5 ) SVRC (Zhou et al., 2018c) Õ( n 4/5 1.5 ) Õ( n 4/5 1.5 ) Lite-SVRC (Zhou et al., 2018b) Õ( n 2/3 2.5 ) Õ( n 2/3 1.5 ) STR 1 (this paper) Õ(min{ n 1.5 ,

EQUATION

1.5 }) STR 2 (this paper) Õ( n 3/4 1.5 ) Õ( n 3/4 1.5 )

Contributions::

In this paper we propose and exploit a formulation in which we make explicit control of the step size in the trust region method. This idea is leveraged to develop two efficient stochastic trust region (STR) approaches. We tailor our methods to achieve state-of-the-art oracle complexities under the following two measurements: (i) the stochastic second-order oracle complexity is prioritized;

(ii) the stochastic first-and second-order oracle complexities are treated equally. Specifically, in Setting (i), our method STR 1 employs a newly proposed estimator to approximate the Hessian and adopts the estimator in (Fang et al., 2018) for gradient approximation. Our novel Hessian estimator maintains an accurate second-order differential approximation with lower amortized oracle complexity. In this way, STR 1 achieves Õ(min{1/ 2 , √ n/ 1.5 }) stochastic second-order oracle complexity. This is lower than existing results for solving problem (1) . In Setting (ii), our method STR 2 substitutes the gradient estimator in STR 1 with one that integrates stochastic gradient and Hessian together to maintain an accurate gradient approximation. As a result, STR 2 achieves convergence in Õ(n 3/4 / 1.5 ) overall stochastic first-and second-order oracle queries. Finally, based on STR, we further develop Hessian-free STR algorithms, namely STR free and STR free +, which outperform existing Hessian-free algorithms theoretically.

1.1. RELATED WORK:

Computing local minimum to a non-convex optimization problem is gaining considerable amount of attentions in recent years. Both cubic regularization (CR) approaches (Nesterov & Polyak, 2006) and trust region (TR) algorithms (Conn et al., 2000; Curtis et al., 2017) can escape saddle points and find a local minimum by iterating the variable along the direction related to the eigenvector of the Hessian with the most negative eigenvalue. As the CR heavily depends on the regularization parameter for the cubic term, Cartis et al. (2011) propose an adaptive cubic regularization (ARC) approach to boost the efficiency by adaptively tuning the regularization parameter according to the current objective decrease. Noting the high cost of full gradient and Hessian computation in ARC, sub-sampled cubic regularization (SCR) (Kohler & Lucchi, 2017a) is developed for sampling partial data points to estimate the full gradient and Hessian. Recently, by exploring the finite-sum structure of the target problem, many works incorporate variance reduction technique (Johnson & Zhang, 2013) into CR and propose stochastic variance-reduced methods. For example, Zhou et al. (2018c) propose stochastic variance-reduced cubic (SVRC) in which they integrate the stochastic variance-reduced gradient estimator (Johnson & Zhang, 2013) and the Hessian tracking technique (Gower et al., 2018) with CR. Such a method is proved to be at least O(n 1/5 ) faster than CR and TR. Then Zhou et al. (2018b) use adaptive gradient batch size and constant Hessian batch size, and develop Lite-SVRC to further reduce the stochastic second-order oracle Õ(n 4/5 / 1.5 ) of SVRC to Õ(n 2/3 / 1.5 ) at the cost of higher gradient computation cost. Similarly, except turning the gradient batch size, Zhang et al. (2018) further adaptively sample a certain number of data points to estimate the Hessian and prove the proposed method to have the same stochastic second-order oracle complexity as Lite-SVRC.

2. PRELIMINARY:

Notation. We use v to denote the Euclidean norm of vector v and use A to denote the spectral norm of matrix A. Let S be the set of component indices. We define the minibatch average of component functions by f (x; S)

def = 1 |S| i∈S f i (x).

Then we specify the assumptions that are necessary to the analysis of our methods.

MetaAlgorithm 1 Inexact Trust Region Method:

Input: initial point x 0 , step size r, number of iterations K, construction of differential estimators g k and H k 1: for k = 0 to K -1 do 2:

Compute h k and λ k by solving (8); 3:

x k+1 := x k + h k ; 4: if λ k ≤ 3 /L 2 then 5: Output x = x k+1 ; 6:

end if 7: end for Assumption 2.1. F is bounded from below and its global optimal is achieved at x * . We denote

∆ = F (x 0 ) -F (x * ). Assumption 2.2. Each f i : R d → R has L 1 -Lipschitz continuous gradient: for any x, y ∈ R d ∇f i (x) -∇f i (y) ≤ L 1 x -y . (3) Assumption 2.3. Each f i : R d → R has L 2 -Lipschitz continuous Hessian: for any x, y ∈ R d ∇ 2 f i (x) -∇ 2 f i (y) ≤ L 2 x -y .

(4)

2.1. TRUST REGION METHOD:

Here we briefly introduce the trust region method (Conn et al., 2000) . In each step, it first solves the Quadratic Constraint Quadratic Program (QCQP) defined as

EQUATION

where r is the trust-region radius. Then it updates the new variable as

EQUATION

)

Since ∇ 2 F (x k ) is indefinite, the trust-region subproblem ( 5) is non-convex. But its global optimizer can be characterized by the following lemma (Corollary 7.2.2 in (Conn et al., 2000) ).

Lemma 2.1. Any global minimizer of problem ( 5) satisfies the equation

EQUATION

where the dual variable λ ≥ 0 should satisfy ∇ 2 F (x k ) + λI 0 and λ( h k -r) = 0.

In particular, the standard QCQP solver returns both the minimizer h k as well as the corresponding dual variable λ of subproblem (5) . In the following section, we first prove that the deterministic trust-region update ( 5) and ( 6) converges at the rate of O(1/k 2/3 ), much sharper than existing provable convergence rate O(1/ √ k) (Conn et al., 2000) , and then develop a more efficient stochastic trust-region approach.

3. METHODOLOGY:

Here we first introduce a general inexact trust region method which is summarized in MetaAlgorithm 1. It accepts inexact gradient estimation g k and Hessian estimation H k as input to the QCQP subproblem

EQUATION

Similar to (5), Lemma 2.1 characterizes the global optimum to problem (8) which can be efficiently solved by Lanczos method (Gould et al., 1999) . Assume the dual variable of the minimizer h k is λ k .

We prove that such inexact trust region method achieves the optimal O(1/k 2/3 ) convergence rate when the estimation g k and H k at each iteration are sufficiently close to their full (exact) counterparts ∇F (x k ) and ∇ 2 F (x k ) respectively:

EQUATION

)

Algorithm 2 STR 1

Input: initial point x 0 , step size r, number of iterations K

1: for k = 1 to K do 2:

Construct gradient estimator g k by Estimator 4;

3:

Construct Hessian estimator H k by Estimator 3; 4:

Compute h k and λ k by solving (8);

5:

x k+1 := x k + h k ; 6: if λ k ≤ 3 /L 2 then 7:

Output x = x k+1 ; 8:

end if 9: end for Such result allows us to derive stochastic trust-region variants with novel differential estimators that are tailored to ensure the optimal convergence rate. We state our formal results in Theorem 3.1, whose proof is deferred to Appendix B.1 due to the space limit. Theorem 3.1 (Main Result). Consider problem (1) under Assumption 2.1-2.3. If the differential estimators g k and H k satisfy Eqn. (9) for all k, MetaAlgorithm 1 finds an O( )-SOSP in less than 5 ) iterations by setting the trust-region radius as r = /L 2 .

K = O( √ L 2 ∆/ 1.

Remark 3.1. We emphasize that MetaAlgorithm 1 degenerates to the exact trust region method by taking g k = ∇F (x k ) and

H k = ∇ 2 F (x k ).

Such result is of its own interest because this is the first proof to show that the vanilla trust region method has the optimal O(1/k 2/3 ) convergence rate. Similar rate is achieved by Curtis et al. ( 2017) but with a much more complicated trust region variant.

Theorem 3.1 shows the explicit step size control of the trust region method: Since the dual variable satisfies λ k > 3 0.5 / √ L 2 > 0 for all but the last iteration, we always find the solution to the trustregion subproblem (8) in the boundary, i.e. h k = r, according to the complementary condition (15) in Appendix B.1. Such exact step size control property is missing in the cubic-regularization method where the step size is implicitly decided by the cubic regularization parameter.

More importantly, we emphasize that such explicit step size control is crucial to the sample efficiency of our variance reduced differential estimators. The essence of variance reduction is to exploit the correlations between the differentials in consecutive iterations. Intuitively, when two neighboring iterates are close, so are their differentials due to the Lipschitz continuity, and hence a smaller number of samples suffice to maintain the accuracy of the estimators. On the other hand, smaller step size reduces the per-iteration objective decrease which harms the convergence rate of the algorithm (see proof of Theorem 3.1). Therefore, the explicit step size control in trust region method allows us to well trade-off the per-iteration sample complexity and convergence rate, from which we can derive stochastic trust region approaches with the state-of-the-art sample efficiency. In contrast, existing trust region methods change the step size at every iteration according to progress made, which requires loss evaluations that can be as expensive as gradient computations (e.g. the non-convex linear model in Section 7) and is thus prohibitive for large-scale problems.

4. STOCHASTIC TRUST REGION METHOD: TYPE I:

Having the inexact trust region method as prototype, we now present our first sample-efficient stochastic trust region method, namely STR 1 , in Algorithm 2 which emphasizes cheaper stochastic second-order oracle complexity. As Theorem 3.1 already guarantees the optimal convergence rate of MetaAlgorithm 1 when the gradient estimator g k and the Hessian estimator H k meet the requirement ( 9), here we focus on constructing such novel differential estimators. Specifically, we first present our Hessian estimator in Estimator 3 and our first gradient estimator in Estimator 4, both of which exploit the trust region radius r = /L 2 to reduce their variances.

4. 

k := ∇ 2 f (x k ; H) -∇ 2 f (x k-1 ; H) + H k-1 ; 8: end if Estimator 4 Gradient Estimator: Case (1) 1: if mod(k, p 1 )= 0 then 2: g k := ∇F (x k ); 3: else 4: Draw s 1 samples indexed by G and g k = ∇f (x k ; G) -∇f (x k-1 ; G) + g k-1 ; 5: end if case ( < O(1/n))

where we compute the full Hessian to avoid approximation error, and Option II is designed for the moderate accuracy case ( > O(1/n)) where we only need an approximate Hessian estimator. Then, p 2 iterations follow with H k defined in a recurrent manner. These recurrent estimators exist for the first-order case (Nguyen et al., 2017; Fang et al., 2018) , but their bound only holds under the vector 2 norm. Here we generalize them into Hessian estimation with matrix spectral norm bound.

The following lemma analyzes the amortized stochastic second-order oracle (Hessian) complexity for Algorithm 3 to meet the requirement in Theorem 3.1. As we need an error bound under the spectral norm, we will appeal to the matrix Azuma's inequality (Tropp, 2012) . The proof is deferred to Appendix B.2. Lemma 4.1. Assume Algorithm 2 takes the trust region radius r = /L 2 as in Theorem 3.1. For any k ≥ 0, Estimator 3 produces estimator H k for the second order differential

EQUATION

s 2 = 32L 1 /( √ L 2 ) log(dK 0 /δ) in option II.

Here K 0 is a constant to be determined later. Consequently the amortized per-iteration stochastic second-order oracle complexity to construct

H k is no more than 2s 2 = min{64 √ n log dK0 δ , 64L1 √ L2 log dK0 δ }.

4.2. GRADIENT ESTIMATOR: CASE (1):

When the stochastic second-order oracle complexity is prioritized, we directly employ the SPIDER gradient estimator to construct g k (Fang et al., 2018) . Similar to the construction for H k , the estimator g k is also constructed in an epoch-wise manner as presented in Estimator 4, where p 1 controls the epoch length and s 1 controls the minibatch size.

We now analyze the stochastic first-order oracle complexity to meet the requirement in Theorem 3.1.

Lemma 4.2. Assume Algorithm 2 takes the trust region radius r = /L 2 . Estimator 4 produces estimator g k of the first order differential ∇F (

x k ) such that g k -∇F (x k ) ≤ /6 with probability at least 1 -δ/K 0 for any k ≥ 0, if we set p 1 = max{1, n L 2 /(cL 2 1 log K0 δ )} and s 1 = min{n, cnL 2 1 log(K 0 /δ)/( L 2 )}

, where the constant c = 1152 and K 0 is a constant to be determined later. Consequently, the amortized per-iteration stochastic first-order oracle complexity to construct

g k is min{n, 4cnL 2 1 log (K 0 /δ)/( L 2 )}.

The proof of Lemma 4.2 is similar to the one of Lemma 4.1 and is deferred to Appendix B.3. These two lemmas only guarantee that the differential estimators satisfy the requirement ( 9) in a single iteration and can be extended to hold for all k by using the union bound with K 0 = 2K, where K denotes the number of iterations. Combining such lifted result with Theorem 3.1, we can establish the computational complexity bound as follows.

Under review as a conference paper at ICLR 2020

Algorithm 5 STR 2 Input: initial point x 0 , step size r, number of iterations K 1: for k = 1 to K do 2:

Construct gradient estimator g k by Estimator 6;

3:

Construct Hessian estimator H k by Estimator 3; 4:

Compute h k and λ k by solving (8);

5:

x k+1 := x k + h k ; 6: if λ k ≤ 3 /L 2 then 7: Output x = x k+1 ; 8:

end if 9: end for Estimator 6 Gradient Estimator: Case (2)

1: if mod(k, p 1 )= 0 then 2: Let x := x k , g k := ∇F (x) 3: else 4:

Draw s 1 samples indexed by G; 5:

g k = ∇f (x k ; G) -∇f (x k-1 ; G) + g k-1 + [∇ 2 F (x) -∇ 2 f (x; G)](x k -x k-1 ); 6: end if Corollary 4.1 (Computational Complexity of STR 1 )

. Assume Algorithm 2 will use Estimator 4 to construct the first-order differential estimator g k and use Estimator 3 to construct the second-order differential estimator H k . To find a 12 -SOSP with probability at least 1 -δ, the overall stochastic first-order oracle complexity is

O(min{ n √ L2∆ 1.5 , √ nL1 2 log( L2∆ δ 1.5 )}) and the overall stochastic second- order oracle complexity is O(min{ √ nL2∆ 1.5 , L1∆ 2 } log( d √ L2∆ δ 1.5 )).

From Corollary 4.1 we see that Õ(min{ √ n/ 1.5 , 1/ 2 }) stochastic second-order oracle queries are sufficient for STR 1 to find an -SOSP which is significantly better than both the subsampled cubic regularization method Õ(1/ 2.5 ) (Kohler & Lucchi, 2017a) and the variance reduction based ones Õ(n 2/3 / 1.5 ) (Zhou et al., 2018b; Zhang et al., 2018) . Recently, Zhou & Gu (2019) developed a stochastic recursive variance-reduced cubic regularization (SRVRC) method which finds an ( , √ )approximate local minimum with Õ(n/ 1.5 , 1/ 3 ) SFO and Õ( √ n/ 1.5 , 1/ 2 ) SSO. But the result of SRVRC needs to assume stochastic gradient to be bounded, i.e., ∇f i (x) -∇F (x) ≤ σ. With this extra assumption, STR1 enjoys Õ(n/ 1.5 , n/ 2 , 1/ 3 ) SFO and Õ( √ n/ 1.5 , 1/ 2 ) SSO. Thus, if 1/ ≤ n ≤ 1/ 2 , STR1 outperforms SRVRC; otherwise they have the same complexities.

5. STOCHASTIC TRUST REGION METHOD: TYPE II:

In the above section, we focus on the setting where the stochastic second-order oracle complexity is prioritized over the stochastic first-order oracle complexity. In this setting, STR 1 achieves the stateof-the-art efficiency. In this section, we consider a different complexity measure where first-order and second-order oracle complexities are treated equally and our goal is to minimize the maximum of them. We note that, currently the best result is Õ(n 4/5 / 1.5 ) of the SVRC method (Zhou et al., 2018c) .

Since the Hessian estimator H k of STR 1 already delivers the superior Õ( √ n/ 1.5 ) stochastic Hessian complexity, in STR 2 (see Algorithm 5), we retain Estimator 3 for second-order differential estimation and use Estimator 6 to further reduce the stochastic gradient complexity.

5.1. GRADIENT ESTIMATOR: CASE (2):

When stochastic gradient and Hessian complexities are equally important, we use Hessian to improve the gradient estimation. Denote x(a) = ax t + (1 -a)x. From Assumption 2.3, we have

∇f i (x t ) -∇f i (x)-∇ 2 f i (x)(x t -x) = 1 0 [∇ 2 f i (x(a)) -∇ 2 f i (x)](x t -x)da ≤ L 2 2 x t -x 2 .

Such property can be used to improve Lemma 4.2 of Estimator 4. Specifically, define the correction

c k = [∇ 2 F (x) -∇ 2 f (x; G)](x k -x k-1 ),

where x is some reference point updated in an epoch-wise manner. Estimator 6 adds c k to the estimator in Estimator 4. Note that in Estimator 6, the first-and second-order oracle complexities are the same. We now analyze the first-order (and second-order) oracle complexity to meet requirement (9) .

Lemma 5.1. Assume Algorithm 5 takes the trust region radius r = /L 2 as in Theorem 3.1. For any k ≥ 0, Estimator 6 produces estimator g k for the first order differential ∇F (x k ) such that g k -∇F (x k ) ≤ /6 with probability at least 1 -δ/K 0 , if we set p 1 = n 0.25 and s 1 = n 0.75 c log(K 0 /δ), where c = 1152 and K 0 is a constant to be determined. Consequently, the amortized per-iteration stochastic first-order oracle complexity to construct g k is 2s 1 = 2n 0.75 c log(K 0 /δ).

The proof of Lemma 5.1 is similar to the one of Lemma 4.1 and is deferred to Appendix B.4. Similar to the previous section, Lemma 5.1 only guarantees that the gradient estimator satisfies the requirement ( 9) in a single iteration. Such result can be extended to hold for all k by using the union bound with K 0 = 2K, which together with Theorem 3.1 gives the following corollary. )) overall stochastic second-order oracle queries.

Corollary 5.1 shows that to find an -SOSP, both SFO and SSO of STR 2 are Õ(n 3/4 / 1.5 ) which surpasses the best existing one Õ(n 4/5 / 1.5 ) in (Zhou et al., 2018c) .

6 PRACTICAL STOCHASTIC TRUST REGION VARIANTS for any desired δ. This repeating technique has been studied by, e.g., (Allen-Zhu & Li, 2018; Allen-Zhu, 2018b) . To test whether the t-th run outputs an O( )-SOSP, we need to compute ∇F (x t ) and the smallest eigenvalue of ∇ 2 F (x t ). The latter one can be approximated by solving the QCQP

EQUATION

where H t is the full Hessian ∇ 2 F (x t ) or its estimation. One can show that MetaAlgorithm 7 finds an O( )-SOSP w.p. at least (1 -δ) in Õ(1/ 1.5 ) iterations. We defer the detailed analysis to Appendix C.

6.2. HESSIAN-FREE IMPLEMENTATION:

Based on MetaAlgorithm 7, we propose a Hessian-free method named STR free , which is summarized in Algorithm 8. STR free leverages the full/stochastic Hessian and Estimator 4 to construct H k and g k , respectively. Besides, it uses Lanczos method (Gould et al., 1999; Carmon & Duchi, 2018) as the QCQP solver, which can be implemented in a Hessian-free manner (i.e., using only Hessian-vector products without explicit Hessian matrix evaluations). Thus, H k is only accessed through Hessianvector products and is never explicitly constructed. Since Hessian-vector products can be computed in linear time (in terms of the dimension d) for many machine learning problems (Allen-Zhu, 2018b; Agarwal et al., 2017) , Hessian-free methods are usually more practical than Hessian based ones for high dimensional problems. The following theorem, whose proof can be found in Appendix D, establishes the runtime complexity (i.e., the total complexity of stochastic gradient and Hessian-vector product evaluations (Zhou & Gu, 2019 )) of STR free .

MetaAlgorithm 7 Inexact Trust Region Method II:

Input: initial point x 0 , step size r, number of inner iterations K, constants δ, ζ ∈ (0, 1), c 1 , c 2 , number of outer iterations T = Θ(log(1/δ)), sample size s (optional) 1: for t = 1 to T do 2:

x t ← INEXACTTR WEAK (x 0 , r, K, ζ);

3::

Option I: high accuracy case (small ) 4:

H t := ∇ 2 F (x t ); Randomly select k from {0, . . . , K -1}; L2 log(4d/δ). The hyper-parameters in Estimator 4 are set to the same values as those in Lemma 4.2. The number of iterations of Lanczos method is set to Õ(1/(L 2 ) 0.25 ). To find an O( )-SOSP w.p. at least 1 -δ, the runtime complexity is Õ(d min{n/ 1.75 , 1/ 2.75 + √ n/ 2 }log(1/δ)).

To solve the QCQP more efficiently, here we develop a faster solver which is based on the AppxPCA method (Allen-Zhu & Li, 2016) and KatyushaX W (Allen-Zhu, 2018a) . See details in Appendix E. By replacing Lanczos method with this solver in STR free , we further improve the runtime complexity to Õ(d min{n/ 1.5 + n 0.75 / 1.75 , 1/ 2.5 + √ n/ 2 }). We call this new algorithm STR free + whose details can be found in Appendix E. Table 2 shows that for the runtime complexities, both STR free and STR free + outperform existing methods. See more comparison and discussion in Appendix E.4.

7. EXPERIMENTS:

Here we compare the proposed STR with several state-of-the-art (stochastic) cubic regularized algorithms and trust region approaches, including trust region (TR) algorithm (Conn et al., 2000) , adaptive cubic regularization (ARC) (Cartis et al., 2011) , sub-sampled cubic regularization (SCR) (Kohler & Lucchi, 2017a), stochastic variance-reduced cubic (SVRC) (Zhou et al., 2018c) , Lite-SVRC (Zhou et al., 2018b) , and SRVRC (Zhou & Gu, 2019) . For STR, we estimate the gradient as the way in case (1) . This is because such a method enjoys lower Hessian computational complexity over the way in case (2) and for most problems, computing their Hessian matrices is much more time-consuming than computing their gradients. For the subproblems in these compared methods, we use Lanczos Table 2 : Runtime complexities of STR free , STR free +, and other state-of-the-art methods.

Algorithm Runtime:

Hessian-free Cubic (Carmon & Duchi, 2016 ) Õ( dn 2 ) Fast-Cubic (Agarwal et al., 2017) Õ( dn 1.5 + dn 0.75 1.75 ) Stochastic Cubic (Tripuraneni et al., 2018) Õ( d 3.5 ) * SRVRC free (Zhou & Gu, 2019) Õ(min{ dn 2 , d 3 }) * STR free (this paper) Õ(min{ dn 1.75 , d 2.75 + dn 0.5 2 }) STR free + (this paper) Õ(min{ dn 1.5 + dn 0.75 1.75 , d 2.5 + dn 0.5 2 }) * These entries rely on an additional assumption: ∇f i (x) -∇F (x) ≤ σ a.s. method (Gould et al., 1999; Kohler & Lucchi, 2017a) to solve the subproblem approximately in a Hessian-related Krylov subspace. We run simulations on seven datasets from LibSVM (a9a, ijcnn, codrna, phishing, w8a, epsilon and mnist). We run our algorithm for 40 epochs and use the output as the optimal value f * for sub-optimality estimation. Note the output has very small gradient already verified by Figure 2 and 4 in appendix. For all the considered algorithms, we set their initializations as zeros and tune their hyper-parameters optimally. For more experimental settings, e.g. details of testing datasets and algorithm parameter settings, please refer to Appendix F.

Two evaluation non-convex problems. Following (Kohler & Lucchi, 2017a; Zhou et al., 2018c) , we evaluate all considered algorithms on two learning tasks: the logistic regression with nonconvex regularizer and the nonlinear least square. Given n data points (x i , y i ) where x i ∈ R d is the sample vector and y i ∈ {-1, 1} is the label, logistic regression with non-convex regularizer aims at distinguishing these two kinds of samples by solving the following problem min w

1 n n i=1 log(1 + exp(-y i w T x i )) + λR(w; α), where the non-convex regularizer R(w; α) is defined as R(w; α) = d i=1 αw 2 i /(1 + αw 2 i ).

The nonlinear least square problem fits the nonlinear data by minimizing min w

1 2n n i=1 y i -φ(w T x i )

2 + λR(w, α). For these two kinds of problems, we set the parameters λ = 10 -3 and α = 10 for all testing datasets.

Comparison of Hessian based algorithms. Figure 1 summarizes testing results on the non-convex logistic regression problem. For each dataset, we report the function value gap v.s. the overall algorithm running time which can reflect the overall computational complexity of an algorithm, and also show the function value gap v.s. Hessian sample complexity which reveals the complexity of Hessian computation. From Figure 1 , one can observe that our proposed STR algorithm runs faster than the compared algorithms in terms of the algorithm running time, showing the overall superiority of STR. Furthermore, STR also reveals much sharper convergence curves in terms of the Hessian sample complexity which is consistent with our theory. This is because to achieve an -accuracy local minimum, the Hessian sample complexity of the proposed STR is Õ(n 0.5 / 1.5 ) and is superior over the complexity of the compared methods (see the comparison in Sec. 4.2). Indeed, this also explains why our algorithm is also faster in terms of algorithm running time, since for most optimization problems, Hessian matrix is much more computationally expensive than the gradient and thus more efficient Hessian sample complexity means faster overall convergence speed. Note, as all compared methods need to compute the Hessian and gradient, their memory complexity are all O(d 2 + d). Comparison of Hessian-free algorithms.

Here we compare our proposed Hessian-free STR, namely STR free , with other state-ofthe-art Hessian-free algorithms on the two high-dimensional datasets, including epsilon and mnist (see details in Appendix F). Here we do not compare STR free +, as it is based on AppxPCA method (Allen-Zhu & Li, 2016) and KatyushaX W (Allen-Zhu, 2018a) which require tuning a lot of hyper-parameters. From the results in Figure 3 , one can observe that compared with other algorithms, our STR free achieves the best convergence speed which demonstrates its high efficiency in realistic applications. Besides, one also can find that STR free is much faster than Hessian based STR since computing full Hessian is actually much computationally expensive than the computation of the Hessian vector.

8. CONCLUSION:

We proposed two stochastic trust region variants. Under two settings (whether stochastic first-and second-order oracle complexities are treated equally), the proposed methods achieve state-of-theart oracle complexities. We also propose Hessian-free variants with lowest runtime complexity. Experimental results testify our theoretical implications and the efficiency of the proposed algorithms.

A APPENDIX:

In this appendix, Sec. B first provides the proofs for the results in the manuscript. Then, we analyze MetaAlgorithm 7 and STR free in Sec. C and Sec. D, respectively. Next, in Sec. E, we develop a fast QCQP solver to further improve the computational complexity of STR free . Finally, more experimental details and results are presented in Sec. F.

B DEFERRED PROOFS:

B.1 PROOF OF THEOREM 3.1

Proof. For simplicity of notation, we denote

∇ k def = ∇F (x k ) -g k and ∇ 2 k def = ∇ 2 F (x k ) -H k .

From Assumption 2.3 we have

F (x k+1 ) ≤F (x k )+ ∇F (x k ), h k + 1 2 ∇ 2 F (x k )h k , h k + L 2 6 h k 3 =F (x k )+ ∇ k +g k , h k + 1 2 [∇ 2 k + H k ]h k , h k + L 2 6 h k 3 .

Use the CauchySchwarz inequality to obtain

EQUATION

)

The requirement ( 9) together with the trust region radius h ≤ r = /L 2 allows us to bound

EQUATION

)

The optimality of (5) indicates that there exists a dual variable λ k ≥ 0 so that (Corollary 7.2.2 in (Conn et al., 2000) )

First Order :

EQUATION

Second Order :

EQUATION

Complementary :

EQUATION

Multiplying ( 13) by h k , we have

EQUATION

)

Additionally, using ( 14) we have

(H k + λ k L 2 2 I)h k , h k ≥ 0,

which together with (16) gives

EQUATION

)

Moreover, the complementary property (15) indicates h k = /L 2 as we have λ k > 3 /L 2 > 0 before MetaAlgorithm 1 terminates. Plug ( 12), ( 16), and ( 17) into (11) and use h k = /L 2 :

EQUATION

)

Therefore, if we have

EQUATION

)

Using Assumption 2.1, we find

λ k ≤ 3 0.5 / √ L 2 in no more than 4 √ L 2 • (F (x 0 ) -F (x * ))/ 1.5 iterations.

We now show that once λ k ≤ 3 0.5 / √ L 2 , then x k+1 is already an O( )-SOSP: From (13), we have

EQUATION

)

The assumptions ∇ k ≤ /6 and ∇ 2 k ≤ √ L 2 /3 together with the trust region radius h ≤ /L 2 imply

EQUATION

)

On the other hand use Assumption 2.3 to bound

∇F (x k+1 ) -∇F (x k ) -∇ 2 F (x k )h k ≤ L 2 2 h k 2 ≤ 2 .

Combining these two results gives ∇F (x k+1 ) ≤ 3 . Besides, using Assumption 2.3, ∇ 2 k ≤ √ L 2 /3, and ( 14), we derive the Hessian lower bound

∇ 2 F (x k+1 ) ∇ 2 F (x k ) -L 2 • h k I H k -L 2 /3I-L 2 h k I -12 L 2 I.

Hence x k+1 is a 12 -stationary point. Additionally, we have h k = r according to the complementary condition ( 15) for all but the last iteration.

B.2 PROOF OF LEMMA 4.1

Proof. Without loss of generality, we analyze the case 0 ≤ k < q 2 for ease of notation. We first focus on Option II. The proof for Option I follows the similar argument.

Option II: Define for k = 0 and i ∈ [s 2 ]

B 0 i def = ∇ 2 f i (x 0 ) -∇ 2 F (x 0 ),

and define for k ≥ 1 and i ∈ [s 2 ]

B k i def = ∇ 2 f i (x k ) -∇ 2 f i (x k-1 ) -(∇ 2 F (x k ) -∇ 2 F (x k-1 )).

{B k i } is a martingale difference sequence. We have for all k and i,

EQUATION

and use Assumption 2.3 for k ≥ 1 to bound

B k i ≤ ∇ 2 f i (x k ) -∇ 2 f i (x k-1 ) + ∇ 2 F (x k ) -∇ 2 F (x k-1 ) ≤ 2 L 2 . From the construction of H k , we have H k -∇ 2 F (x k ) = s 2 i=1 B 0 i s 2 + k j=1 s2 i=1 B j i s 2 .

Thus using the matrix Azuma's Inequality in Theorem 7.1 of (Tropp, 2012) and k ≤ p 2 , we have

P r{ H k -∇ 2 F (x k ) ≥ t} ≤d • exp{- t 2 /8 s 2 i=1 4L 2 1 /s 2 2 + k j=1 s2 i=1 4 L 2 /s 2 2 } ≤d • exp{- t 2 /8 4L 2 1 /s 2 + 4p 2 L 2 /s 2 }.

Consequently, we have

P r{ H k -∇ 2 F (x k ) ≤ L 2 } ≥ 1 -δ/K 0 . by taking t = √ L 2 , s 2 = 16L 2 1 /( L 2 ) log(dK 0 /δ), s 2 = 32L 1 /( √ L 2 ) log(dK 0 /δ), and p 2 = L 1 /(2 √ L 2 ).

Option I: The proof is similar to the one of Option II except that we replace B 0 i with zero matrix. In such case, the matrix Azuma's Inequality implies

P r{ H k -∇ 2 F (x k ) ≥ t} ≤ d • exp{- t 2 /8 k j=1 s2 i=1 4 L 2 /s 2 2 } ≤ d • exp{- t 2 /8 4p 2 L 2 /s 2 }.

Thus by taking t = √ L 2 , s 2 = 32 √ n log(dK 0 /δ), and p 2 = √ n, we have the result.

Amortized Complexity: In option I, the choice of parameters ensures that: s 2 ≤ p 2 × s 2 and in option II: n ≤ p 2 × s 2 . Consequently, the amortized stochastic second-order oracle complexity is bounded from above by 2s 2 .

B.3 PROOF OF LEMMA 4.2

Without loss of generality, we analyze the case 0 ≤ k < q 1 for ease of notation. Define for k ≥ 1

and i ∈ [s 1 ] a k i def = ∇f i (x k ) -∇f i (x k-1 ) -(∇F (x k ) -∇F (x k-1 )).

{a k i } is a martingale difference sequence: for all k and i

EQUATION

From the construction of g k , we have

g k -∇F (x k ) = k j=1 s1 i=1 a j i s 1 .

Recall the Azuma's Inequality. Using k ≤ p 1 , we have

P r{ g k -∇F (x k ) ≥ t} ≤ exp{- t 2 /8 k j=1 s1 i=1 4 L 2 1 L2s 2 1 } ≤ exp{- t 2 /8 4 L 2 1 p 1 /(s 1 L 2 )

}.

Take t = /6 and denote c = 1152. To ensure that

P r{ g k -∇F (x k ) ≥ /6} ≤ δ/K 0 , we need cL 2 1 L2 log K0 δ ≤ s1 p1 .

The best amortized stochastic first-order oracle complexity can be obtained by solving the following two-dimensional programming:

min p1≥1,s1≥1 (n + s 1 (p 1 -1))/p 1 s.t. cL 2 1 L 2 log K 0 δ ≤ s 1 p 1 ,

which has the solution

s 1 = min{n, n • cL 2 1 log K 0 δ L2 }, and p 1 = max{1, n • L2 cL 2 1 log K 0 δ }.

Note that when we take s 1 = n, we directly compute g k = ∇F (x k ) without sampling.

The amortized stochastic first-order oracle complexity is obtained by plugging in the choice of s 1 and p 1 , which completes the proof.

B.4 PROOF OF LEMMA 5.1

Without loss of generality, we analyze the case 0 ≤ k < q 1 for ease of notation. Define for k ≥ 1

and i ∈ [s 1 ] b k i def = ∇f i (x k ) -∇f i (x k-1 ) -∇ 2 f i (x)(x k -x k-1 ) -[∇F (x k ) -∇F (x k-1 ) -∇ 2 F (x)(x k -x k-1 )]. {b k i } is a martingale difference sequence: for all k and i E[b k i |x k ] = 0. Besides, b k i has bounded norm: b k i ≤ ∇f i (x k )-∇f i (x k-1 ) -∇ 2 f i (x)(x k -x k-1 ) + ∇F (x k )-∇F (x k-1 ) -∇ 2 F (x)(x k -x k-1 ) . We can bound ∇f i (x k )-∇f i (x k-1 ) -∇ 2 f i (x)(x k -x k-1 ) as follows. ∇f i (x k ) -∇f i (x k-1 ) -∇ 2 f i (x)(x k -x k-1 ) = 1 0 ∇ 2 f i (x k-1 + t(x k -x k-1 )) -∇ 2 f i (x) (x k -x k-1 )dt ≤ 1 0 L 2 tx k + (1 -t)x k-1 -x dt • x k -x k-1 ≤ 1 0 t x k -x + (1 -t) x k-1 -x dt • L 2 r ≤ L 2 kr 2 ,

where the first inequality follows from Assumption 2.3 and the last inequality holds because x kx ≤ kr and x k-1 -x ≤ kr, where r is the trust region radius. Similarly, we have

∇F (x k ) - ∇F (x k-1 ) -∇ 2 F (x)(x k -x k-1 ) ≤ L 2 kr 2 . Thus, we bound b k i ≤ 2L 2 kr 2 ≤ 2p 1 From the construction of g k , we have g k -∇F (x k ) = k j=1 s1 i=1 b j i s 1 .

We use k ≤ p 1 and the Azuma's inequality to bound

P r{ g k -∇F (x k ) ≥ t} ≤ exp{- t 2 /8 k j=1 s1 i=1 4p 2 1 2 s 2 1 } ≤ exp{- t 2 /8 4 2 p 3 1 /s 1 }.

Thus, by taking t = /6 and c = 1152, we need s1

p 3 1 ≥ c log K0 δ . Further we want s 1 p 1 O(n)

and hence we take p 1 = n 0.25 and s 1 = n 0.75 c log K0 δ . The amortized stochastic first-order oracle complexity is bounded by 2s 1 .

C ANALYSIS OF METAALGORITHM 7:

We first show that INEXACTTR WEAK finds an O( )-SOSP in O(1/ 1.5 ) iterations with probability at least 2/3 as stated in the following lemma. Lemma C.1. Consider problem (1) under Assumptions 2.1-2.3. Suppose that the differential estimators g k and H k satisfy Eqn. ( 9) with probability at least (1 -ζ 4K ). Besides, suppose that hk is an approximate solution to (8) such that w.p. (1

EQUATION

)

where h k is a global solution to (8). By setting ζ = 1/3, r = /L 2 , and K = 4 √ L 2 ∆/ 1.5 , INEXACTTR WEAK outputs a 500 -SOSP w.p. at least 2/3. Proof. Combining ( 11) and ( 24), we have w.p. (1

EQUATION

where h k is a global solution to the QCQP (8) and hk is an approximate solution satisfying (24). We let λ k denote the dual variable corresponding to the global solution h k as defined in Lemma 2.1. We note that h k and λ k are used only in our analysis. The INEXACT WEAK algorithm only requires the approximate solution xk without knowledge of h k or λ k .

By the assumption that ( 9) holds with probability (1

EQUATION

)

Plugging ( 16), ( 17), and ( 26) into ( 25) and applying the union bound, we have w.p. at least (1

EQUATION

)

where the second inequality follows from (15):

EQUATION

Summing inequality ( 27) from k = 0 to K -1 and applying the union bound, we have w.p. at least

EQUATION

)

where the second inequality follows from Assumption 2.1 and our choice of the trust region radius.

By sampling k uniformly from {0, . . . , K -1}, we obtain

EQUATION

where the expectation is taken over the randomness of k. Combining ( 29) and (30) and taking K = 4∆ √ L 2 / 1.5 , we have w.p. at least (1 -ζ/2)

EQUATION

)

Since λ k is always no-negative, by Markov's inequality and the union bound, with probability at least 1 -ζ, we have

EQUATION

)

By taking ζ = 1/3, we have w.p. at least 2/3, λ k ≤ 42 /L 2 . The rest of the proof is similar to Theorem 3.1 and we have the result.

The following theorem shows that MetaAlgorithm 7 finds an O( )-SOSP w.p. L2 log(4d/δ). Then MetaAlgorithm 7 finds a 600 -SOSP with probability at least (1 -δ).

Proof. By Lemma C.1 and our choice of T , with probability (1 -δ/2), at least one of x t is a 500 -SOSP. On the other hand, since ψ t (ṽ

EQUATION

where the last inequality follows from our choice of c 2 .

Option I: Since H t = ∇ 2 F (x t ) is the full Hessian, ψ t (v t ) is the smallest eigenvalue of ∇ 2 F (x t ).

Applying the union bound, we conclude that MetaAlgorithm 7 outputs a 600 -SOSP w.p. (1 -δ).

Option II:

EQUATION

By Assumption 2.2, we have

EQUATION

Applying the matrix Azuma's Inequality in Theorem 7.1 of Tropp (2012) leads to

EQUATION

).

(36)

By taking s = 32L 2 1 L2 log(4d/δ) and applying the union bound, we have with probability 1 -δ,

EQUATION

where the last inequality follows from (33). This completes the proof.

D PROOF OF THEOREM 6.1

Proof. We first analyze the computational cost of Lanczos method. By Corollary 2 in (Carmon & Duchi, 2018) , for any desired accuracy ˜ , Lanczos method achieves this accuracy in O( r √ ˜ log r √ d ˜ p ) Lanczos iterations w.p. at least (1 -p) . Without loss of generality, we assume that the number of Lanczos iterations is strictly smaller than the dimension d, otherwise the QCQP subproblem can be solved exactly. We note that each Lanczos iteration involves computation of one matrix-vector product. Therefore, to satisfy the condition (24) in Lemma C.1, one needs to evaluate Õ(1/(L 2 ) 0.25 ) Hessian-vector products of the form H k v. Similarly, to solve (10) up to accuracy √ L 2 w.h.p., one needs to evaluate Õ(1/(L 2 ) 0.25 ) Hessian-vector products of the form H t v.

In MetaAlgorithm 7, to verify whether the candidate solution v t is indeed an O( )-SOSP, one needs at most O(n) stochastic gradient evaluations and O(min{n, log(4d/δ)L 2 1 /(L 2 )}/(L 2 ) 0.25 )) stochastic Hessian-vector product evaluations, where the latter one follows from the proof of Theorem 7. We proceed to analyze the computational complexity of the INEXACTTR WEAK procedure. Recall that the iteration complexity of MetaAlgorithm 7 is O(log(1/δ)/ 1.5 ). Following Lemma 4.2 and Corollary 4.1, the stochastic first-order oracle complexity is Õ(min{n/ 1.5 , √ n/ 2 }log(1/δ)). Following the proof of Lemma 4.1 and Corollary 4.1, when p 2 = 1, the overall stochastic Hessian sample complexity is Õ(min{n/ 1.5 , 1/ 2.5 } log(1/δ)). Since it takes Õ(1/ 0.25 ) Lanczos iterations to meet the condition (24) as stated above, the overall stochastic Hessian-vector product oracle complexity is Õ(min{n/ 1.75 , 1/ 2.75 }log(1/δ)). Combining the stochastic first-order and Hessian-vector product complexities, the overall runtime is Õ(d min{n/ 1.75 , 1/ 2.75 + √ n/ 2 }log(1/δ)).

E A FASTER HESSIAN-VECTOR BASED QCQP SOLVER:

We recall from the previous section that, to approximately solve a quadratic subproblem in STR free , Lanczos method requires Õ(min{n/ 0.25 , 1/ 1.25 }) stochastic Hessian-vector product evaluations.

In this section, we propose a faster QCQP solver with an Õ(min{n + n 0.75 / 0.25 , 1/ }) complexity.

Replacing Lanczos method with this QCQP solver in STR free results in a faster Hessian-free method, which we refer to as STR free +.

E.1 CONVEX REFORMULATION OF QCQP:

To begin with, we present a known result that is key to achieve faster algorithm than Lanczos method. We summarize this result in Lemma E.1 which shows that the trust region subproblem is equivalent to a convex QCQP. Lemma E. 1. (Convex Reformulation of QCQP (Flippo & Jansen, 1996; Wang & Xia, 2017 )) Denote λ min as the smallest eigenvalue of H k . Let u min be a corresponding eigenvector. W.l.o.g., we assume that g k , u min ≤ 0. Let µ = min{λ min , 0}. Then the QCQP (8) is equivalent to the convex problem

min h∈R d , h ≤r q k (h) = g k , h + 1 2 (H k -µI)h, h + 1 2 µr 2 (38)

in the sense that (8) and (38) have the same minimum function value. Moreover, when λ min < 0, for any optimal solution of (38), denoted by h k c ,

h k c + h k c , u min 2 -u min 2 ( h k c 2 -r 2 ) -h k c , u min u min 2 u min (39)

is a global minimizer of the original QCQP (8).

To perform the above reformulation, one needs to compute the exact eigenpair (λ min , u min ). Nevertheless, as we shall see, it is sufficient to compute an approximate eigenpair ( λ, ũ) such that

EQUATION

) where ˜ is a target accuracy to be determined later. We note that ˜ ≤ 2L 2 w.l.o.g. since H k ≤ L 2 . With this approximate eigenpair, it remains to solve the following convex problem

EQUATION

where μ = min{0, λ -˜ }. One can check that the problem (41) well approximates (38). Corollary E.1. Let q k * and qk * be the minimum function value of (38) and (41), respectively. Assume

EQUATION

)

We note that the above convex reformulation approach divides an indefinite QCQP into two subproblems: (i) computation of an approximate eigenpair ( λ, ũ); (ii) solving the convex problem (41). As we shall see, by exploiting the finite-sum structure of the Hessian H k , these two subproblems can be efficiently solved. We treat these two subproblems in the following two subsections, respectively.

E.2 FINDING THE SMALLEST EIGENVECTOR:

To find a unit vector that satisfies requirement (40), we resort to the AppxPCA method (Allen-Zhu & Li, 2016) , which first finds an approximate eigenvalue λ = λ min -˜ via binary search and then applies Power method to the positive definite matrix (H k -λI) -1 for a logarithmic number of iterations. Computing (H k -λI) -1 v for any vector v is equivalent to solving the ˜ -strongly convex problem (Allen-Zhu & Li, 2018 )

EQUATION

We note that

H k = 1 |S| i∈S ∇ 2 f i (x k ).

Specifically, in STR free , either |S| = n (i.e., H k is the full Hessian) or |S| = Õ(L 2 1 /(L 2 )) by Lemma 4.1. Therefore, φ k (•) can be expressed as sum of non-convex functions

EQUATION

By observing that each φ k i is non-convex and has (4L 2 )-Lipschitz gradient, we can use KatyushaX S (Allen-Zhu, 2018a) to solve problem (43) in Õ(|S| + |S| 3/4 L 2 /˜ ) stochastic Hessianvector product (i.e., ∇ 2 f i (x k )u) evaluations. The following result is taken from (Agarwal et al., 2017, Section G.3) , which gives the overall computation complexity of AppxPCA.

Algorithm 9 Fast QCQP Solver:

Input: H k , g k , r, ˜ , ˜ 1 1: Use AppxPCA to find ( λ, ũ) satisfying ( 40), in which the matrix inverse is solved by KatyushaX S ; 2: Use KatyushaX W to solve (41) up to accuracy ˜ 1 , i.e., find a vector h such that qk (h) -qk * ≤ ˜ 1 with high probability;

3: Return h + ( h, ũ 2 -ũ 2 ( h 2 -r 2 ) -h, ũ )ũ/ ũ 2 .

Algorithm 10 STR free + Option I:

H k := ∇ 2 F (x k ); 5:

Option II: Draw s samples indexed by H and let

H k := ∇ 2 f (x k ; H); 6: use Algorithm 9 to solve QCQP subproblems. Lemma E.2. Let H k = 1 |S| i∈S ∇ 2 f i (x k ) ∈ R d×d , where ∇ 2 f i (x k ) ≤ L 2 .

With probability at least 1 -p, AppxPCA produces a unit vector u satisfying u T H k u ≤ λ min + ˜ . The total stochastic Hessian-vector product oracle complexity is Õ(|S| + |S| 3/4 L 2 /˜ ).

E.3 SOLVING THE CONVEX QCQP:

In what follows, we show that the convex problem (41) can be solved efficiently. We first observe that problem (41) has a finite-sum structure and can be rewritten as an unconstrained problem of the form

EQUATION

where Ψ(h) = 0 if h ≤ r, otherwise Ψ(h) = +∞. We note that each qk i (h) in ( 45) has (4L 2 )-Lipschitz continuous gradient since ∇ 2 f i (x k ) ≤ L 2 and ˜ ≤ 2L 2 . Therefore, we can use KatyushaX W (Allen-Zhu, 2018a) to solve (45). By (Allen-Zhu, 2018a, Theorem 4.6) 

, KatyushaX W finds a point h such that E[q k (h) -qk * ] ≤ ˜ 1 using Õ(|S| + |S| 3/4 √ L 2 • r/ √ ˜ 1

) stochastic Hessianvector products, where ˜ 1 is the target accuracy to be determined later.

E.4 PUTTING IT ALL TOGETHER:

The complete procedure of our fast QCQP solver is summarized in Algorithm 9. Combining all the above results and setting r = /L 2 , ˜ = √ L 2 /2, and ˜ 1 = ˜ r 2 , one can find an approximate solution to QCQP (8) satisfying requirement (24) in Õ(|S| + |S| 3/4 L 0.25 2 / 0.25 ) stochastic Hessianvector product evaluations. By replacing Lanczos method with this solver in STR free , we derive a new Hessian-free method called STR free +, which is summarized in Algorithm 10. The following theorem establishes the overall runtime complexity of STR free + for finding an -SOSP. L2 log(4d/δ). The hyper-parameters in Estimator 4 are set to the same values as those in Lemma 4.2. Besides, let ˜ = √ L 2 /2 and ˜ 1 = ˜ r 2 in Algorithm 9. To find an O( )-SOSP w.p. at least 1 -δ, the runtime complexity is Õ(d min{n/ 1.5 + n 0.75 / 1.75 , 1/ 2.5 + √ n/ 2 }log(1/δ)).

Proof. The proof directly follows from that in Sec. D.

We compare the runtime complexity of STR free and STR free + with existing Hessian free methods in Table 2 . One can see that STR free strictly outperforms Hessian-free Cubic. Besides, STR free outperforms Fast-Cubic if n ≥ Ω(1/ 4/3 ), which is a mild condition for large-scale problems in the moderate accuracy case. STR free + strictly outperforms both Hessian-free Cubic and Fast-Cubic. We note that the runtime analyses in (Tripuraneni et al., 2018; Zhou & Gu, 2019) rely on an additional assumption which states that for all x, with probability 1,

EQUATION

)

Under this additional assumption, one can use the same argument as in Section B.2 and B.3 to prove that STR free achieves a runtime complexity of Õ(d min{n/1.75 , n 0.5 /2 + 1/ 2.75 , 1/ 3 }) 1 .

Similarly, the runtime complexity of STR free + would be Õ(d min{n/ 1.5 + n 0.75 / 1.75 , 1/ 2.5 + √ n/ 2 , 1/ 3 }). In this sense, both STR free and STR free + outperform Stochastic Cubic and SRVRC free . F ADDITIONAL EXPERIMENTAL RESULTS

F.1 MORE EXPERIMENTAL DETAILS:

Descriptions of Testing Datasets. We briefly introduce the seven testing datasets in the manuscript. Among them, six datasets are provided by the LibSVM website 2 , including (a9a, ijcnn, codrna, phishing, w8a and epsilon). The detailed information is summarized in Table 3 . We can observe that these datasets are different from each other in feature dimension, training samples, etc.

Experimental Settings. In the manuscript, following SVRC (Zhou et al., 2018c) and Lite-SVRC (Zhou et al., 2018b) , we select hyper parameters from a set, namely s 1 from {0.2n, 0.6n, n}, s 2 from {0.01n, 0.1n, 0.2n}, p 1 and p 2 from {0.01n 0.5 , 0.05n 0.5 , 0.1n 0.5 }. For the Hessian estimation at the beginning of each p 2 iterations, we use full Hessian. Similarly, for the gradient estimation at the beginning of each p 1 iterations, we adopt the full gradient as the gradient estimation.

Memory Analysis. SVRC, Lite-SVRC, and our method need to store the previous and current gradient and Hessian and thus their memory complexity is 2(d 2 + d). TR, CR (ARC) and SCR need to compute current and Hessian and thus has complexity d 2 + d. So these memory is of the same order but our method is much faster than TR, CR and SCR both validated by theory and experiments.

F.2 MORE EXPERIMENTS:

Here we give more experimental results on the gradient norm v.s. the algorithm running time and the Hessian sample complexity. Due to the space limit, in the manuscript we only provide the gradient-norm related results on the codrna dataset. Here we provide results of a9a and ijcnn datasets in Figure 4 . One can observe that on both the logistic regression with non-convex regularizer and the nonlinear least square problems, the proposed algorithm always shows sharper convergence behavior in terms of both the running time and the Hessian sample complexity. These observations are consistent with the results in Figure 2 in the manuscript. All these results demonstrate the high efficiency of our proposed algorithm and also confirm our theoretical implication.

:

To obtain this complexity, one needs to replace the full gradient in line

of Estimator 4 (and line 6 of MetaAlgorithm 7) with a mini-batch stochastic gradient when n ≥ Ω(1/ 2 ).2 https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/

Figures and Tables:

Corollary 5.1 (Computational Complexity of STR 2 ). Algorithm 5 finds a 12 -SOSP with probability at least 1 -δ, within O( n 0.75 √ L2∆ 5 )) overall stochastic first-order oracle queries and O( n 0.75 √ L2∆

H k := ∇ 2 F (x k ); 5: Option II: Draw s samples indexed by H and let H k := ∇ 2 f (x k ; H); 6: use Lanczos method to solve QCQP subproblems. Theorem 6.1. Consider Algorithm 8 for solving problem (1). Let ζ = 1/3, r = /L 2 , K = 4 √ L 2 ∆/ 1.5 , T = 3 2 log(2/δ), c 1 = 600, c 2 = 500, and s = 32L 2 1

Figure 1: Comparison on the logistic regression with non-convex regularizer.

Figure 2: Comparison on the nonlinear least square problem.

Figure2displays results of the compared algorithms on the nonlinear least square problem. STR shows very similar behaviors as those in Figure1. Specifically, STR achieves fastest convergence rate in terms of both algorithm running time and Hessian sample complexity. On the codrna dataset we further plot the gradient norm versus running time and Hessian sample complexity. One can obverse that the gradient in STR vanishes significantly faster than other algorithms which means that STR can find the stationary point with high efficiency. See Figure4in Appendix F.2 for more experimental results on gradient norm comparison. All these results confirm the superiority of the proposed STR.

Figure 3: Comparison of Hessian-free algorithms on non-convex logistic and soft-max regressions.

(1 -δ) after running INEXACTTR WEAK for Θ(log(1/δ)) times. Theorem C.1 (Iteration Complexity of MetaAlgorithm 7). In the same setting as Lemma C.1, let T = 3 2 log(2/δ), c 1 = 600, c 2 = 500, and s = 32L 2 1

Consider Algorithm 10 for solving problem(1).Let ζ = 1/3, r = /L 2 , K = 4 √ L 2 ∆/ 1.5 , T = 3 2 log(2/δ), c 1 =600, c 2 = 500, and s = 32L 2 1

Figure 4: Comparison of gradient norm on both the non-convex logistic regression and nonlinear least square problems.

1 HESSIAN ESTIMATOR Our epoch-wise Hessian estimator H k is given in Estimator 3, where p 2 controls the epoch length and s 2 (and optionally s 2 ) controls the minibatch size. At the beginning of each epoch, Estimator 3 has two options, designed for different target accuracy: Option I is preferable for the high accuracy

Table content: <table><tr><td colspan="3">Estimator 3 Hessian Estimator</td></tr><tr><td colspan="3">Input: Epoch length p 2 , sample size s 2 , s 2 (optional)</td></tr><tr><td colspan="3">1: if mod(k, p 2 )= 0 then</td></tr><tr><td>2:</td><td>Option I:</td><td>high accuracy case (small )</td></tr><tr><td>3:</td><td>H k := ∇ 2 F (x k );</td><td/></tr><tr><td>4:</td><td>Option II:</td><td>low accuracy case (moderate )</td></tr><tr><td>5:</td><td colspan="2">Draw s 2 samples indexed by H and let H k := ∇ 2 f (x k ; H );</td></tr><tr><td colspan="2">6: else</td><td/></tr><tr><td>7:</td><td colspan="2">Draw s 2 samples indexed by H and let H</td></tr></table>

Descriptions of the five testing datasets.

Table content: <table><tr><td/><td colspan="2">#sample #feature</td><td/><td colspan="2">#sample #feature</td></tr><tr><td>a9a</td><td>32,561</td><td>123</td><td>w8a</td><td>49,749</td><td>300</td></tr><tr><td>ijcnn</td><td>49,990</td><td>22</td><td>phishing</td><td>7,604</td><td>68</td></tr><tr><td>codrna</td><td>28,305</td><td>8</td><td>mnist</td><td>60,000</td><td>784</td></tr><tr><td>epsilon</td><td>40,000</td><td>2,000</td><td/><td/><td/></tr></table>

Peer review (iclr2020):
In this paper, the authors improved the state-of-the-art result by using inexact gradient and Hessian estimation in the training and proving under this case, it still have good performance. In order to control the difference between gradient and approximation gradient, the author use variance reduced estimator to exploit the correlation between consecutive iterations. In addition, the author consider for two cases with the importance of second-order oracle and use Hessian to improve the approximation of gradient when first-order and second-order oracle are equally important. Furthermore, the authors proposed a refined algorithm in practice which only uses stochastic gradient and Hessian-vector product information and showed better experiment results. In general, the paper is well written and easy to follow, but I still have some questions about this paper. First, there lacks explanation about why using biased estimators can still guarantee the convergence. For estimator 3 and estimator 4, it seems like that the stochastic gradient and Hessian are not unbiased approximation of the true value. That is a bit strange, since most popular stochastic estimators including stochastic gradient and stochastic variance-reduced gradient are unbiased approximation, which could guarantee that when the number of samples is large enough, the estimators can approximate the true quantities very well. More discussion about this issue is needed. . Second, to the best of my knowledge, the trust region radius is changing during the iteration for basic trust region algorithm. However, in meta-algorithm 1, the radius is fixed to a very small quantity which is related to the accuracy epsilon. I am very surprised that with a fixed small radius, STR can still achieve the best result among all baseline algorithms, especially compared with trust region algorithm with adaptive radius. Can the authors explain this phenomenon? Third, there is no discussion on space complexity. In practice, it is important to consider the space complexity. However, this work did not provide any space complexity analysis, especially to compare with first-order algorithms. It is interesting to see the trade-off between space complexity and time complexity among both first-order algorithms and second-order algorithms. Fourth, the results of Lemma 4.1 and 4.2 are all in high probability. When the probability is small, then the number of sample is large. Thus, it is not fair that the authors simply ignored them when they compared their algorithms with deterministic algorithm. Furthermore, in practice, to choose small may cause the total number of samples very big. The authors need to make more comments on this issue. Finally, I think the experiment results contradict the theoretical results. From figure 1, it can be seen that the STR method including many other baseline algorithms converge to the global minimum with linear convergence speed (a9a, epoch). However, the convergence analysis provided by the authors claim that the convergence speed is sublinear. I believe such a difference is due to a fact that the initialized parameter is near to the global minimum, thus the optimization landscape here is actually convex but not non-convex. That makes the experiment results vacuous.