Combining Stochastic Adaptive Cubic Regularization with Negative Curvature for Nonconvex Optimization

TL;DR

This paper introduces SANC, a novel optimization algorithm that combines negative curvature with adaptive cubic regularization, improving convergence and practicality for large-scale nonconvex machine learning problems.

Contribution

It proposes the first method integrating negative curvature with adaptive cubic regularization, enhancing efficiency and applicability in large-scale nonconvex optimization.

Findings

SANC effectively escapes saddle points.

Experimental results show improved convergence on neural network problems.

The method is practical for large-scale machine learning applications.

Abstract

We focus on minimizing nonconvex finite-sum functions that typically arise in machine learning problems. In an attempt to solve this problem, the adaptive cubic regularized Newton method has shown its strong global convergence guarantees and ability to escape from strict saddle points. This method uses a trust region-like scheme to determine if an iteration is successful or not, and updates only when it is successful. In this paper, we suggest an algorithm combining negative curvature with the adaptive cubic regularized Newton method to update even at unsuccessful iterations. We call this new method Stochastic Adaptive cubic regularization with Negative Curvature (SANC). Unlike the previous method, in order to attain stochastic gradient and Hessian estimators, the SANC algorithm uses independent sets of data points of consistent size over all iterations. It makes the SANC algorithm more practical to apply for solving large-scale machine learning problems. To the best of our knowledge, this is the first approach that combines the negative curvature method with the adaptive cubic regularized Newton method. Finally, we provide experimental results including neural networks problems supporting the efficiency of our method.

Tables

Table 1. Table 1 : Datasets used for numerical experiments. n 𝑛 n stands for the number of datapoints, f 𝑓 f stands for the number of features, and c 𝑐 c stands for the number of classes. L 1 subscript 𝐿 1 L_{1} and L 2 subscript 𝐿 2 L_{2} represent the parameters used in the SANC and NCD method.

	$n$	$f$	$c$	learning model	$L_1$	$L_2$
w1a	2477	300	2	logistic regression	10.0	10.0
w8a	49749	300	2	logistic regression	10.0	10.0
a9a	32561	123	2	logistic regression	10.0	10.0
ijcnn1	49990	22	2	logistic regression	10.0	100.0
covtype	581012	54	2	logistic regression	10.0	100.0
higgs	11000000	28	2	logistic regression	10.0	100.0
segment	2310	19	3	multi-layer neural networks	100.0	100.0
seismic	78823	50	7	multi-layer neural networks	100.0	100.0
MNIST	60000	28$\times$28$\times$1	10	convolutional neural networks	100.0	100.0
CIFAR10	50000	32$\times$32$\times$3	10	convolutional neural networks	1000.0	1000.0

Equations

$$\underset{\mathbf{x}\in {\mathbb{R}}^d}{\min }f(\mathbf{x})=\mathbb{E}\left[F(\mathbf{x},\xi )\right]$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^d}{\min }f(\mathbf{x})\text{mathchar}58=\frac{1}{n}\sum _{i=1}^{n}f(\mathbf{x};{\mathbf{q}}_i,r_i)=\frac{1}{n}\sum _{i=1}^n{f}_i(\mathbf{x})$$

$$\mathbf{g}(\mathbf{x})\text{mathchar}58=\frac{1}{|{\mathcal{S}}_{\mathbf{g}}|}\sum _{i\in {\mathcal{S}}_{\mathbf{g}}}\nabla f_i(\mathbf{x})$$

$$\mathbf{B}(\mathbf{x})\text{mathchar}58=\frac{1}{|{\mathcal{S}}_{\mathbf{B}}|}\sum _{i\in {\mathcal{S}}_{\mathbf{B}}}{\nabla }^2{f}_i(\mathbf{x})$$

$$\Vert \nabla f({\mathbf{x}}^{\ast })\Vert \le ϵ\text{ and }{\lambda }_{\min }\left(H({\mathbf{x}}^{\ast })\right)\ge -ϵ$$

$${\rho }_t=\frac{f({\mathbf{x}}_t)-f({\mathbf{x}}_t+{\mathbf{s}}_t)}{f({\mathbf{x}}_t)-{\tilde{m}}_t({\mathbf{s}}_t)}$$

$${\sigma }_{t+1}=\{\begin{array}{ll}\max \{\min \{{\sigma }_t,\Vert {\mathbf{g}}_t\Vert \},{ϵ}_m\}& \text{if }{\rho }_t>{\eta }_2\text{ [very successful iteration]}\\ {\sigma }_t& \text{if }{\eta }_2\ge {\rho }_t\ge {\eta }_1\text{ [successful iteration]}\\ \gamma {\sigma }_t& \text{otherwise }\text{ [unsuccessful iteration]}\end{array}$$

$${\tilde{m}}_t(\mathbf{s})\text{mathchar}58=f({\mathbf{x}}_t)+{\mathbf{g}}_t^T\mathbf{s}+\frac{1}{2}{\mathbf{s}}^T{\mathbf{B}}_t\mathbf{s}+\frac{1}{3}{\sigma }_t\Vert \mathbf{s}{\Vert }^3$$

$${\tilde{m}}_t({\mathbf{s}}_t)\le {\tilde{m}}_t({\mathbf{s}}_t^c),{\mathbf{s}}_t^c=-{\alpha }_c{\mathbf{g}}_t$$

$${\mathcal{K}}_j({\mathbf{g}}_t,{\mathbf{B}}_t)=span\left\{{\mathbf{g}}_t,{\mathbf{B}}_t{\mathbf{g}}_t,{\mathbf{B}}_t^2{\mathbf{g}}_t,\dots ,{\mathbf{B}}_t^{j-1}{\mathbf{g}}_t\right\},$$

$${\mathbf{s}}_t={\mathbf{Q}}_t{\mathbf{u}}_t$$

$${\mathbf{u}}_t=\arg \underset{\mathbf{u}\in {\mathcal{K}}_j({\mathbf{g}}_t,{\mathbf{B}}_t)}{\min }{\gamma }_0{\mathbf{u}}^T{e}_1+\frac{1}{2}{\mathbf{u}}^T\mathbf{Tu}+\frac{1}{3}{\sigma }_t\Vert \mathbf{u}{\Vert }^3$$

$${\mathbf{g}}_t^T{\mathbf{s}}_t+{\mathbf{s}}_t^T{\mathbf{B}}_t{\mathbf{s}}_t+{\sigma }_t\Vert {\mathbf{s}}_t{\Vert }^3=0$$

$${\mathbf{s}}_t^T{\mathbf{B}}_t{\mathbf{s}}_t+{\sigma }_t\Vert {\mathbf{s}}_t{\Vert }^3\ge 0$$

$${\mathbf{Q}}_t^T{\mathbf{B}}_t{\mathbf{Q}}_t=\mathbf{T}.$$

$${\lambda }_{\min }\left({\mathbf{B}}_t\right)\ge {\mathbf{v}}^T{\mathbf{B}}_t\mathbf{v}-{ϵ}^{\prime }$$

$${\mathbf{d}}_t=\{\begin{array}{ll}-\frac{2|{\mathbf{v}}_t^T{\mathbf{B}}_t{\mathbf{v}}_t|}{L_2}z{\mathbf{v}}_t& \text{if }\frac{2(-{\mathbf{v}}_t^T{\mathbf{B}}_t{\mathbf{v}}_t{)}^3}{3L_2^2}-\frac{ϵ({\mathbf{v}}_t^T{\mathbf{B}}_t{\mathbf{v}}_t{)}^2}{6L_2^2}>\frac{\Vert {\mathbf{g}}_t{\Vert }^2}{4L_1}-\frac{{ϵ}_{\mathbf{g}}^2}{L_1}\\ -\frac{1}{L_1}{\mathbf{g}}_t& \text{otherwise}\end{array}$$

$$\Vert \nabla f({\mathbf{x}}_{{\mathbf{t}}^{\ast }})\Vert \le ϵ\text{ and }-{\lambda }_{\min }({\mathbf{Q}}_{t^{\ast }}^T{\mathbf{B}}_{t^{\ast }}{\mathbf{Q}}_{t^{\ast }})\le ϵ.$$

$$\Vert {\mathbf{g}}_t-\nabla f({\mathbf{x}}_t)\Vert \le {ϵ}_{\mathbf{g}}$$

$$\Vert {\mathbf{B}}_t-\mathbf{H}({\mathbf{x}}_t)\Vert \le {ϵ}_{\mathbf{B}}$$

$$\Vert {\nabla }_s{\tilde{m}}_t({\mathbf{s}}_{t,k})\Vert \le {ϵ}_{\mathbf{g}}$$

$${\mathbf{d}}_t=0\text{ and }\Vert {\mathbf{s}}_t\Vert \le {ϵ}_s$$

$$|{\mathcal{S}}_{\mathbf{g},t}|\ge \frac{4L_0^2{\left(1+2\sqrt{\log \frac{1}{\delta }}\right)}^2}{{ϵ}_{\mathbf{g}}^2},$$

$$|{\mathcal{S}}_{\mathbf{B},t}|\ge \frac{16L_1^2\log \frac{2d}{\delta }}{{ϵ}_{\mathbf{B}}^2},$$

$${\sigma }_{\min }\le {\sigma }_t\le {\sigma }_{\max }\text{ for all }t\ge 0$$

$$f({\mathbf{x}}_t)-\mathbb{E}[f({\mathbf{x}}_t+{\mathbf{d}}_t)]\ge \max \left\{\frac{{{ϵ}^{\prime }}^3}{24L_2^2},\frac{\Vert {\mathbf{g}}_t{\Vert }^2}{4L_1}-\frac{{ϵ}_{\mathbf{g}}^2}{L_1}\right\}$$

$$f({\mathbf{x}}_t)-{\tilde{m}}_t({\mathbf{s}}_t)\ge \frac{1}{6}{\sigma }_t\Vert {\mathbf{s}}_t{\Vert }^3$$

$$\Vert {\mathbf{s}}_t\Vert \ge \frac{1}{{\kappa }_l}\left(-{ϵ}_{\mathbf{B}}+\sqrt{{ϵ}_{\mathbf{B}}^2+2{\kappa }_l\left(\Vert \nabla f({\mathbf{x}}_t+{\mathbf{s}}_t)\Vert -2{ϵ}_{\mathbf{g}}\right)}\right)$$

$${\kappa }_l\text{mathchar}58=L_2+2{\sigma }_{\max }$$

$$1>ϵ>\max \{3{ϵ}_{\mathbf{g}},144{ϵ}_{\mathbf{B}}^2,\frac{2}{{\kappa }_l}{ϵ}_{\mathbf{B}}^2\},$$

$$l_1+l_2+1=\mathcal{O}({ϵ}^{-3})$$

$$\displaystyle l_{1}\mathrel{\mathop{\mathchar 58\relax}}=\bigg{\lceil}\frac{f(\mathbf{x}_{0})-f_{low}}{\min\{\kappa_{s1}\epsilon^{\frac{3}{2}},\kappa_{d}\epsilon^{3}\}}\bigg{\rceil},$$

Full text

In this paper, we suggest an algorithm combining negative curvature with the adaptive cubic regularized Newton method to update even at unsuccessful iterations. We call this new method Stochastic Adaptive cubic regularization with Negative Curvature (SANC). Unlike the previous method, in order to attain stochastic gradient and Hessian estimators, the SANC algorithm uses independent sets of data points of consistent size over all iterations. It makes the SANC algorithm more practical to apply for solving large-scale machine learning problems. To the best of our knowledge, this is the first approach that combines the negative curvature method with the adaptive cubic regularized Newton method. Finally, we provide experimental results including neural networks problems supporting the efficiency of our method.

Machine Learning, ICML, Stochastic Optimization, Adaptive Cubic Regularization, Negative Curvature

1 Introduction

We have focused on solving the nonconvex unconstrained problems. Especially, in machine learning problems with large dataset, we frequently exploit stochastic optimization schemes; when it comes to solving an empirical risk minimization(ERM) problem, instead of calculating exact gradient and Hessian at an iterate, their estimators are used to find an update step. In this context, in order to find second order critical point, two of the most popular optimization approaches using Hessian are considered:

(i) the negative curvature method and

(ii) the Newton method.

Given an iterate, the negative curvature methods (Curtis & Robinson, 2017; Liu et al., 2018; Cano et al., 2017) need to calculate the eigenvector corresponding to the left-most eigenvalue of Hessian. Using this eigenvector along with a novel approach to switch negative curvature direction and gradient descent direction, this method shows its convergence to a second order critical point. We remark that because calculating exact eigenvector and corresponding left-most eigenvalue are time-consuming, this method may use approximate approaches such as the Lanczos method or the Oja’s algorithm (Oja, 1982).

The second approach is the Newton method. (Martens, 2010; Martens & Sutskever, 2011; Agarwal et al., 2017b; Vinyals & Povey, 2012) Given an iterate, to find an update step, Newton method solves a linear system which denotes the first order necessary condition of optimality for a second-order Taylor approximation of an objective function. To prevent the storage of the Hessian matrix which can be prohibitive to store especially for large-scale problems, We usually use Krylov subspace methods to iteratively attain an approximated solution of the linear system over a Krylov subspace. The Krylov subspace methods also utilize Hessian-vector products (Pearlmutter, 1994) to construct a Krylov subspace whose dimension increases linearly over iterations.

In a similar vein, the cubic regularized Newton method (Nesterov & Polyak, 2006) has been used to escape from a strict saddle point and ultimately converge to a second-order stationary point for nonconvex problems. More recently, Adaptive Regularization algorithm using Cubics (ARC) and its stochastic variant (Cartis et al., 2011a, b; Kohler & Lucchi, 2017; Bergou et al., 2018) are proposed. The ARC algorithm adopts the trust-region concept to determine if the obtained Newton update step is successful or not. The next iterate is updated with the update step only when the update step is successful.

We focus our attention on integrating the negative curvature method into the ARC algorithm to update with a negative curvature direction when update Newton step is not successful . Also, it is noted that calculating a negative curvature direction is efficient because it is on the same Krylov subspace where the update step of ARC is obtained.

Thus, we make the following contributions in this paper.

• •

We provide an optimization algorithm which mainly uses the ARC algorithm and also exploit negative curvature when the obtained update step is not successful.

• •

We also provide the theoretical results which include worst-case iteration complexity for achieving approximate first- and second-order optimality and corresponding lower bounds on the cardinalities of the data subsets on which stochastic gradient and Hessian estimators are calculated.

• •

We verify the efficiency of our method numerically by comparing with other standard first- and second-order methods for nonconvex optimization to show its prowess especially for neural networks problems.

2 Related Works

Second-order methods In the stochastic optimization setting, second-order methods have been researched for decades. J. Martens (Martens, 2010; Martens & Sutskever, 2011) published papers about the Hessian-free optimization method to train neural networks including deep auto-encoder and recurrent neural networks. In these papers, it is argued that using a Gauss-Newton approximation instead of a Hessian approximation is more effective to train neural networks. Also, he employs a damping coefficient to maintain the positive definiteness of a Hessian approximation, which helps restrict a step size. Consequently, with a damping coefficient $\lambda$, a linear system is solved, $(\mathbf{H}+\lambda I)\mathbf{s}=-\mathbf{g}$, to attain an update step $\mathbf{s}$ by using the conjugate gradient algorithm, at every iterate. The idea to use a damping coefficient is highly similar to that of adopting an adaptive cubic regularization term. But it is noted that adaptive cubic regularization has a theoretical adaptive rule to guarantee its convergence whereas the Newton method with a damping coefficient does not provide.

To circumvent solving a linear system, many researchers have considered the quasi-Newton method (Wang et al., 2017b; Byrd et al., 2016) as an alternative. In stochastic optimization, they applied an L-BFGS update formula using stochastic gradient and Hessian estimators. However, since an assumption of positive definiteness of a Hessian is required to prove its convergence, it is not capable to provide convergence guarantees for nonconvex optimization problems. Experience has shown that some gains in performance in machine learning applications can be achieved, but the full potential of the stochastic quasi-Newton schemes is not yet known. (Bottou et al., 2018)

Cubic regularized Newton method The cubic regularized Newton method is proposed by Griewank (Griewank, 1981) firstly and is also studied by Nesterov et al. (Nesterov & Polyak, 2006). In the deterministic setting, the use of a cubic regularization as an over-estimator of an objective function to solve unconstrained optimization problems has been studied. A second-order Taylor approximation with a cubic regularization term is constructed and solved at each iteration. Nesterov and Polyak (Nesterov & Polyak, 2006) provided the theoretical analysis of its global convergence rate for nonconvex unconstrained problems.

Their theoretical results for the cubic regularization paved the way for the upcoming ARC algorithm (Cartis et al., 2011a, b). This algorithm uses an adaptive coefficient in the cubic regularization term. This coefficient varies according to the difference between the local cubic model value and actual function value at the next iterate. The way of adjusting a cubic coefficient in the ARC algorithm is similar to that of adjusting a radius of the trust region methods. They also showed the mathematical proof for the global worst-case complexity bounds.

More recently, Kohler and Lucchi (Kohler & Lucchi, 2017) provided stochastic ARC variant for solving finite-sum structure nonconvex objectives. They proved the lower bound on cardinalities of the subsampling subsets to calculate stochastic gradient and Hessian estimators. But in their method, the sizes of the subsets are increasing over iterations which hinders training large scale neural networks.

Wang et al. (Wang et al., 2017b, 2018) proposed approaches to incorporate a momentum acceleration, or a stochastic variance reduced gradient(SVRG)((Johnson & Zhang, 2013)) into a framework of the cubic regularized Newton method. However, they did not use an adaptive cubic coefficient in their paper.

Negative curvature approaches There have been relatively little works on the negative curvature methods. Curtis and Robinson (Curtis & Robinson, 2017) proposed several algorithms exploiting negative curvature for solving deterministic and stochastic optimization problems. In this work, a current iterate is updated with a direction of (stochastic) gradient descent and negative curvature. The ’dynamic method’ is also proposed which adaptively estimates a Lipschitz constant of gradient continuity to estimate a step length in the stochastic optimization framework. They showed that the ’dynamic method’ is efficient in performance in stochastic optimization problems, but they did not provide a proof of its convergence.

Most recently, in a similar vein, Liu et al. (Liu et al., 2018) proposed the adaptive negative curvature descent (NCD) method which gives some adaptability to the terminating criteria of the Lanczos procedure depending on the magnitude of the subsampled gradients. They also provided variants of the adaptive NCD method whose worst-case time complexity is $\mathcal{\widetilde{O}}(d/\epsilon^{3.5})$ for stochastic optimization, where $\mathcal{\widetilde{O}}$ hides logarithmic terms.

Carmon et al. (Carmon et al., 2018) used negative curvature directions at the first phase of iterates and then switched it to accelerated stochastic gradient descent method when an iterate reaches an almost convex region.

Also, a direction of negative curvature has been used for escaping from a strict saddle point. For the details on it, please refer to (Reddi et al., 2017; Xu et al., 2018).

3 Formulation of Stochastic Adaptive Cubic Regularization with Negative Curvature

3.1 Problem statement

We consider the following stochastic optimization problem:

[TABLE]

where $F$ is a twice continuously differentiable function and possibly nonconvex, $\mathbf{x}\in\mathbb{R}^{d}$ is a variable, and $\mathbb{E}$ denotes the expectation with respect to $\xi$, a random variable with the distribution $P$. If it is assumed that there is no prior knowledge about the distribution $P$, one may use the following finite-sum structure, as a proxy of the problem, of the form 1.

[TABLE]

where a function $f_{i}$ is an abbreviation of $f(\mathbf{x};\mathbf{q}_{i},r_{i})$, and $\{(\mathbf{q}_{1},r_{1}),\dots,(\mathbf{q}_{n},r_{n})\}$ are $n$ sampled dataset. This optimization problem 2 is usually referred to as an empirical risk minimization (ERM) problem. We frequently encounter this kind of optimization problems in supervised machine learning. For example, given a dataset, where $\mathbf{q}_{i}\in\mathbb{R}^{m}$ is a feature and $r_{i}\in\mathbb{R}$ is its corresponding label for all $i\in\{1,\dots,n\}$, $f_{i}(\mathbf{x})=1-\tanh(r_{i}\mathbf{x}^{T}\mathbf{q}_{i})+\lambda\lVert\mathbf{x}\rVert_{2}^{2}$ is a nonconvex support vector machine (SVM) (Wang et al., 2017a) with a convex regularization term.

Also, when it comes to online or stochastic optimization, it is assumed that it is not possible to access directly to an exact function value, gradient, or Hessian at $\mathbf{x}$, i.e., $f(\mathbf{x})$, $\nabla f(\mathbf{x})$, or $\mathbf{H}(\mathbf{x})$. Instead, we may resort to using unbiased estimators of the gradient and Hessian As a prototypical method to attain stochastic gradient and Hessian estimators, one frequently consider using the mini-batch approach. We construct sets of data points independently drawn from the dataset, $\mathcal{S}_{\mathbf{g}}$, $\mathcal{S}_{\mathbf{B}}$. The stochastic gradient and Hessian estimators, say, $\mathbf{g}(\mathbf{x})$ and $\mathbf{B}(\mathbf{x})$, can be defined as,

[TABLE]

In this context, we aim to find an $\epsilon$-approximate second-order stationary point $\mathbf{x}^{*}$ of the problem 2 which satisfies,

[TABLE]

where $\lambda_{\min}(\cdot)$ denotes the left-most eigenvalue.

For brevity, we use $\lVert\cdot\rVert$ as an abbreviation of the $l2$ norm, $\lVert\cdot\rVert_{2}$, in what follows. For the same reason, we also use $\mathbf{g}_{t}$, and $\mathbf{B}_{t}$ instead of $\mathbf{g}(\mathbf{x}_{t})$, and $\mathbf{B}(\mathbf{x}_{t})$, respectively.

3.2 Algorithm

In this section, we describe a general approach to our method which is presented as Algorithm 1. We call this method Stochastic Adaptive cubic regularization with Negative Curvature (SANC) in what follows. The main framework of Algorithm 1 is same as that of the ARC algorithm. Given an iterate $\mathbf{x}_{t}$, the algorithm determines if the Newton step $\mathbf{s}_{t}$ is successful or not. It depends on the prescribed parameter $\eta_{1}$ and $\eta_{2}$. It also depend on a measure $\rho_{t}$ in Algorithm 1 which calculates the fraction of the predicted model value decrease $f(\mathbf{x}_{t})-\widetilde{m}_{t}(\mathbf{s}_{t})$ and the actual decrease $f(\mathbf{x}_{t})-f(\mathbf{x}_{t}+\mathbf{s}_{t})$. If the difference between the predicted and the actual decrease is small enough, then the iteration may become successful, otherwise, it is unsuccessful.

If the iteration is successful, the next iterate is set to $\mathbf{x}_{t+1}=\mathbf{x}_{t}+\mathbf{s}_{t}$. If the iteration is unsuccessful, we try to find a negative curvature direction to update.

After updating the iterate, we need to update the cubic coefficient. When the iteration is very successful, an adaptive cubic coefficient $\sigma_{t}$ is decreasing according to Eq. 15 in Algorithm 1. Also, if the iteration is unsuccessful, $\sigma_{t}$ is increasing, which leads to the Newton step with a small magnitude. Thus, this adaptive cubic coefficient $\sigma_{t}>0$ plays a similar role as a reciprocal of radius in the trust region methods.

We describe the details on how we attain a Newton step and a negative curvature direction.

3.3 Finding Newton step from local cubic model

In this subsection, we elucidate a way to calculate a Newton step $\mathbf{s}_{t}$. Given an iterate $\mathbf{x}_{t}$, we can attain the stochastic gradient and Hessian estimator, $\mathbf{g}_{t}$, $\mathbf{B}_{t}$, over $\mathcal{S}_{\mathbf{g},t}$ and $\mathcal{S}_{\mathbf{B},t}$, respectively. With $\mathbf{g}_{t}$, $\mathbf{B}_{t}$, and a sufficiently large $\sigma_{t}>0$, we establish an approximate local cubic model $\widetilde{m}_{t}(\mathbf{s})$ of the form,

[TABLE]

We can attain a global minimizer $\mathbf{s}_{t}$ corresponding to Eq. 6 by $\mathbf{s}_{t}\mathrel{\mathop{\mathchar 58\relax}}=\arg\!\min_{\mathbf{s}\in\mathbb{R}^{d}}{\widetilde{m}_{t}(\mathbf{s})}$ regardless of the positive definiteness of $\mathbf{B}_{t}$. As Cartis et al. (Cartis et al., 2011a) pointed out, attaining the global minimizer of Eq. 6 is time-consuming. So they suggested that solving Eq. 6 is relaxed by satisfying the Cauchy condition,

[TABLE]

where $\alpha_{c}\in\arg\!\min_{\alpha>0}\widetilde{m}_{t}(-\alpha\mathbf{g}_{t})=f(\mathbf{x}_{t})-\alpha\lVert\mathbf{g}_{t}\rVert^{2}+\frac{1}{2}\alpha^{2}\mathbf{g}_{t}^{T}\mathbf{B}_{t}\mathbf{g}_{t}+\frac{1}{3}\sigma_{t}\alpha^{3}\lVert\mathbf{g}_{t}\rVert^{3}$. Also, they provided a convergence guarantee with the satisfaction of the Cauhy condition.

As a consequence, the use of the Lanczos method can be considered to seek an approximate minimizer of Eq. 6 over the Krylov subspace. The Lanczos method is a kind of iterative methods to solve a linear system. The Lanczos method constructs an orthogonal matrix, $\mathbf{Q}_{t}$, sequentially whose columns are a basis of a Krylov subspace $\mathcal{K}_{j}$ of $\mathbb{R}^{d}$. It is noted that if the Krylov subspace is formed by $\mathbf{g}_{t}$ and $\mathbf{B}_{t}$ successively as

[TABLE]

then $\mathcal{K}_{j}(\mathbf{g}_{t},\mathbf{B}_{t})$ always contain a subspace $\mathbf{g}_{t}$ even at Lanczos iteration $j=1$. Hence, an approximate minimizer $\mathbf{s}_{t}$ over the Krylov subspace $\mathcal{K}_{j\geq 1}(\mathbf{g}_{t},\mathbf{B}_{t})$ always holds the Cauchy condition 7 when preconditioning is not applied. As a consequence, the approximate minimizer $\mathbf{s}_{t}$ can be retrieved as

[TABLE]

where

[TABLE]

In Eq. 10, $\gamma_{0}=\lVert\mathbf{g}_{t}\rVert$, $e_{1}$ is $1st$ column vector of identity matrix and $\mathbf{T}$ is a $j$ by $j$ symmetric tridiagonal matrix. We would like to point out that solving Eq. 10 is more efficient than solving Eq. 6 because of its smaller dimension and a structured construction of $\mathbf{T}$. Carmon and Duchi (Carmon & Duchi, 2016) argued that it takes $\mathcal{O}(\log(1/\epsilon))$ iterations with a gradient descent method to find an $\epsilon$-approximate second-order stationary point for Eq. 6 with mild assumptions. Thus, it is obvious that solving Eq. 10 requires much less iterations than solving Eq. 6 because it uses a structured tridiagonal matrix and has smaller dimension.

As shown in Lemma 3.2 in (Cartis et al., 2011a), given an iterate $\mathbf{x}_{t}$, an global and approximate minimizer of $\widetilde{m}_{t}(\mathbf{s})$ over the Krylov subspace always satisfy the following requirements,

[TABLE]

These statements 11a,11b play a vital role to prove the proofs of the main theorem of our algorithm in Section 4.

It is not surprising that the Cauchy point $\mathbf{s}_{t}^{c}$ also satisfy Eq. 11a and Eq. 11b. If we use a larger Krylov subspace, we can attain more accurate $\mathbf{s}_{t}$ which improves the performance of the algorithm. However, at the same time, at each Lanczos iteration, we need $\mathcal{O}\left(|\mathcal{S}_{\mathbf{B}}|d\right)$ operations so we have to choose the termination condition of the Lanczos iteration carefully. In the next subsection, we describe details on how to attain a direction of negative curvature which can be utilized at unsuccessful iterations.

3.4 Finding negative curvature

Given an iterate $\mathbf{x}_{t}$, suppose that the stochastic Hessian estimator $\mathbf{B}_{t}$ is indefinite and $\mathbf{s}_{t}^{T}\mathbf{B}_{t}\mathbf{s}_{t}<0$. With a fixed $\sigma_{t}$, the magnitude of the calculated Newton step $\mathbf{s}_{t}$ should be relatively larger than that with an arbitrary positive definite Hessian because it is needed to satisfy Eq. 11b. So, we can conjecture that if $\sigma_{t}$ is not large enough, then $\mathbf{x}_{t}+\mathbf{s}_{t}$ is not in the vicinity of $\mathbf{x}_{t}$. As a consequence, it is likely that the Newton step $\mathbf{s}_{t}$ becomes unsuccessful.

At this point, if we know that the stochastic Hessian estimator $\mathbf{B}_{t}$ is indefinite, then we can utilize a negative curvature direction to update an iterate. Luckily, from the previous subsection, the Newton step $\mathbf{s}_{t}$ is achieved by the Lanczos method, and if the orthogonal matrix $\mathbf{Q}_{t}$ is saved or regenerated after terminating the Lanczos iterations, we can utilize the following relationship,

[TABLE]

It is noted that since it is structured, eigenpairs $(\lambda_{i}^{(j)},\mathbf{v}_{i}^{(j)})$ of $\mathbf{T}$ can be cheaply computed. Also the left-most eigenvalue $\lambda_{1}^{(j)}$ and corresponding eigenvector $\mathbf{v}_{1}^{(j)}$ of $\mathbf{T}$ are used to approximate the left-most eigenpair of $\mathbf{B}_{t}$. The approximate left-most eigenvector of $\mathbf{B}_{t}$ can be attained by, $\mathbf{v}_{t}=\mathbf{Q}_{t}\mathbf{v}_{1}^{(j)}$. We can approximate the left-most eigenpair of $\mathbf{B}_{t}$ by that of $\mathbf{T}$, i.e., $(\lambda_{t},\mathbf{v}_{t})=(\lambda_{1}^{(j)},\mathbf{Q}_{t}\mathbf{v}_{1}^{(j)})$. This approximation is often called the Ritz approximation.

Please note that with a small $\epsilon^{\prime}>0$, there exists $\mathbf{v}\in\mathbb{R}^{d}$ such that

[TABLE]

where $\mathbf{v}$ can be achieved within $\mathcal{\widetilde{O}}\left(\frac{d}{\sqrt{\epsilon^{\prime}}}\right)$ by the Lanczos method, as proved in (Kuczyński & Woźniakowski, 1992).

Usually, the Lanczos iteration is truncated at $j\ll d$ especially for large-scale problems in order to make the algorithm efficient. Also, computing eigenpairs of $\mathbf{T}$ takes only $\mathcal{O}(j^{2})$ operations because $\mathbf{T}$ is structured. So, it is computationally cheap to calculate the left-most eigenpair of $\mathbf{T}$ in large-scale problems. Also, one can use divide-and-conquer algorithm (Coakley & Rokhlin, 2013) to calculate the eigenpair of $\mathbf{T}$.

In order to prove the convergence guarantee of our algorithm, we maintain a framework of the ARC algorithm and use a negative curvature step $\mathbf{d}_{t}$ at an unsuccessful iteration (Liu et al., 2018) as

[TABLE]

where $z$ is the Rademacher random variable. In Eq. 14, $\frac{2(-\mathbf{v}_{t}^{T}\mathbf{B}_{t}\mathbf{v}_{t})^{3}}{3L_{2}^{2}}-\frac{\epsilon(\mathbf{v}_{t}^{T}\mathbf{B}_{t}\mathbf{v}_{t})^{2}}{6L_{2}^{2}}$ is the expected decrease when the direction is negative curvature, whereas $\frac{\lVert\mathbf{g}_{t}\rVert^{2}}{4L_{1}}-\frac{\epsilon_{\mathbf{g}}^{2}}{L_{1}}$ is the expected decrease when stochastic gradient descent is used. In general, $\mathbf{d}_{t}$ is a combination of negative curvature and stochastic gradient descent which leads to decreasing objective a lot as possible. Thus, the SANC updates at every iteration even if the Newton direction is not successful, which enables the algorithm to be more efficient than the ARC algorithm. Now, we present the theorem regarding the convergence and the sampling approximations of gradients and Hessian in the next section.

4 Theoretical Analysis

In this section, we provide the worst-case iteration complexity for achieving approximate first- and second-order optimality of the SANC algorithm. Firstly, we explain some assumptions regarding the gradient and Hessian approximations. Also, we describe some assumptions regarding the convergence of the algorithm and the terminating criteria of the Lanczos method, which are crucial to the proof of the convergence of the SANC algorithm.

4.1 Definition and assumptions

Definition 1** (Approximate Second-Order Stationary Point).**

With a small $\epsilon>0$, $x_{t^{*}}$ of iterates $\{x_{t}\}_{t\geq 0}$ is called an $\epsilon$-approximate second-order stationary point if it satisfies

[TABLE]

Notably, even though the $\epsilon$-approximate second-order stationary point does not necessarily denote that it is close to an exact second-order stationary point, but if strict saddle property (Lee et al., 2016) holds, then it is guaranteed that the $\epsilon$-approximate second-order stationary point is in the vicinity to the exact second-order stationary point for sufficiently small $\epsilon$.

To establish the convergence, we need to assume the following statements.

Assumption 1** (Function Boundness).**

Nonconvex function $f_{i}$ is twice differentiable and bounded below by $f_{inf}>-\infty$ for all $i$. Thus, nonconvex function $f$ is also twice differentiable and bounded below by $f_{inf}>-\infty$.

Assumption 1 plays an important role to impose an upper bound on the sum of the sufficient decreases of two consecutive function values.

Assumption 2** (Lipschitz Continuity).**

For all $i$, the function $f_{i}$, $\nabla f_{i}$, and $\mathbf{H}_{i}$ are Lipschitz continuous with Lipschitz constants $L_{0}$, $L_{1}$, and $L_{2}$, respectively.

The assumption that $f_{i}$ is Lipschitz continuous is relatively rare to use, but in our analysis, it is used to bound the cardinality of the set for stochastic gradient estimator. Instead, we can also assume the exponential tail behavior of gradients realization.

To attain stochastic gradient and Hessian estimators, we resort to the following assumptions to independently draw the sets from the dataset.

Assumption 3** (Gradient and Hessian Approximation Bounds).**

For all iteration $t\geq 0$,

[TABLE]

With sufficiently small $\epsilon_{\mathbf{g}}>0$ and $\epsilon_{\mathbf{B}}>0$

Assumption 3 is different with those in (Cartis et al., 2011a; Kohler & Lucchi, 2017). In their assumption, $\lVert\mathbf{s}_{t}\rVert$ was included on the RHS of the inequality. As an iteration point is approaching an optimum, the discrepancy between exact and approximate gradient/Hessian needs a tighter agreement because $\lVert\mathbf{s}_{t}\rVert$ tends to decrease. Smaller $\lVert\mathbf{s}_{t}\rVert$ makes the size of the sets increasing, which is harmful to apply for solving large-scale machine learning problems.

On the contrary, in our analysis, we can maintain the size of subsets consistent because of Assumption 3. It makes our algorithm more effective to solve large-scale problems.

The following assumption about the termination criteria of the Lanczos method is used to establish the proof of Lemma 6.

Assumption 4** (Lanczos Method Termination Criteria).**

For each iteration $t$, let us assume that the Lanczos procedure stops at the Lanczos iteration $k$ when the following criterion holds,

[TABLE]

The above assumption is used when we establish the proof of Lemma 6. The proof uses the norm of discrepancy between $\nabla f(\mathbf{x}_{t}+\mathbf{s}_{t})$ and $\nabla_{s}\widetilde{m}_{t}(\mathbf{s}_{t})$. Also, we need to assume the following statement to establish the upper bound on $\sigma_{t}$, which has a vital role in the convergence proof of our method.

Assumption 5** (Convergence).**

For Algorithm 1, the process $\{\mathbf{x}_{t}\}_{t\geq 0}$ stops when the following criteria hold,

[TABLE]

4.2 Gradient and Hessian sampling bound

Now, we establish the lower bound on the cardinalities of the sets for stochastic gradient and Hessian estimators. The following theorems are based on the Bernstein inequalities in (Gross, 2011).

Theorem 1** (Gradient Sampling Bound).**

If

[TABLE]

then $\mathbf{g}_{t}$ satisfies the gradient approximation bound Eq. 16a for all $t\geq 0$ with high probability $1-\delta$.

Theorem 2** (Hessian Sampling Bound).**

If

[TABLE]

then $\mathbf{B}_{t}$ satisfies the Hessian approximation bound Eq. 16b for all $t\geq 0$ with high probability $1-\delta$.

The proofs regarding Theorem 1, 2 are in the supplementary material. As shown in the Theorem 1, 2, the cardinalities of the sets $\mathcal{S}_{\mathbf{g}}$, $\mathcal{S}_{\mathbf{B}}$, do not depend on the norm of $\mathbf{s}$. Because these cardinalities bounds rely only on the constants, $L_{1}$, $d$, $\delta$, and $\epsilon_{\mathbf{B}}$, they can be consistent over all iterations. Now, with the sampling bounds, we can establish the worst-case iteration complexity for achieving approximate first- and second-order optimality.

4.3 Worst-case iteration complexity bound for achieving approximate first- and second-order optimality

Lemma 3** (Cubic Regularization Coefficient Bound).**

Let the A2, A3 and A5 hold. For finite positive values $\sigma_{\min}$ and $\sigma_{\max}$, it holds

[TABLE]

where $\sigma_{\max}\mathrel{\mathop{\mathchar 58\relax}}=\max\left\{\sigma_{0},\gamma\left(\frac{3}{2}L_{2}+\frac{3(\epsilon_{\mathbf{g}}+\frac{1}{2}\epsilon_{\mathbf{B}})}{\epsilon_{s}^{2}}\right)\right\}$.

The following lemma is about the sufficient decrease with the negative curvature step $\mathbf{d}_{t}$.

Lemma 4** (Sufficient Decrease with Negative Curvature Step, Lemma 3 in (Liu et al., 2018)).**

When $\mathbf{v}_{t}^{T}\mathbf{B}_{t}\mathbf{v}_{t}\leq-\epsilon^{\prime}/2$ and $\epsilon_{\mathbf{B}}\leq\epsilon^{\prime}/12$, the negative curvature step $\mathbf{d}_{t}$ satisfies that

[TABLE]

with high probability $1-\delta$.

In Lemma 4, we use the expectation of $f(\mathbf{x}_{t}+\mathbf{d}_{t})$ with respect to the Rademacher random variable $z$. This random variable takes the main role to simplify the lower bound of Eq. 22. Before describing the sufficient decrease with the Newton step $\mathbf{s}_{t}$, the following lemmas are needed.

Lemma 5** (Local Cubic Model Decrease, Lemma 3.3 in (Cartis et al., 2011a)).**

Suppose that $\mathbf{s}_{t}$ satisfies Eq. 11a and Eq. 11b. Then for all (very) successful iterations $t$,

[TABLE]

Lemma 6** (Lower Bound of $\mathbf{s}_{t}$).**

Suppose that A2, A3, A4, and A5 hold. For all (very) successful iterations $t$, the Newton step $\mathbf{s}_{t}$ satisfies that

[TABLE]

where

[TABLE]

Plugging Eq. 24 on the RHS of Eq. 23 yields the sufficient decrease with the Newton step $\mathbf{s}_{t}$ which is used in the proof of the following main theorem of our algorithm.

Theorem 7** (Worst-Case Iteration Complexity for Approximate First- and Second-Order Optimality).**

Let A1, A2, A3, A4, and A5 hold. For

[TABLE]

Algorithm 1 provides an iteration $t^{*}$ such that $\lVert\nabla f(\mathbf{x}_{t^{*}})\rVert\leq\epsilon$ within at most $\mathcal{O}(\epsilon^{-3/2})$ iterations.

Also it provides an iteration $t^{*}$ such that $\lVert\nabla f(\mathbf{x}_{t^{*}})\rVert\leq\epsilon$ and $-\lambda_{\min}(\mathbf{Q}_{t^{*}}^{T}\mathbf{B}_{t^{*}}\mathbf{Q}_{t^{*}})\leq\epsilon$ within at most

[TABLE]

iterations. Both are with high probability $1-\delta$. Here,

[TABLE]

Remark 1 It is noted that Eq 26 is intuitively reasonable because gradient and Hessian approximation bounds assumption 3 does not depend on the norm of $\mathbf{s}_{t}$. Using smaller $\epsilon_{\mathbf{g}}$ and $\epsilon_{\mathbf{B}}$, iterates of the SANC algorithm can converge to a closer point to an exact optimum.

Remark 2 Theorem 7 imply that the time complexity of the SANC algorithm is about $\mathcal{O}(\epsilon^{-3.5})$ when equipped with the Lanczos method for attaining steps, $\mathbf{s}_{t}$ or $\mathbf{d}_{t}$ to find a point satisfying second order optimality.

Remark 3 For the first order optimality, the iteration complexity $\mathcal{O}(\epsilon^{-1.5})$ of Algorithm 1 is same as that of ARC algorithm. If it is equipped with a method with per iteration complexity of $\mathcal{O}(\epsilon^{-0.25})$, the SANC find first order optimality in a time $\mathcal{O}(\epsilon^{-1.75})$ as (Agarwal et al., 2017a; Carmon et al., 2018).

All the proofs for the above theorems are provided in the appendix.

5 Empirical Studies

5.1 Practical implementation of SANC

As a default setting, we used $\gamma=2$, $\eta_{1}=0.2$, $\eta_{2}=0.8$, and $\sigma_{0}=1.$. These parameter settings are same as those in SCR (Kohler & Lucchi, 2017). For the neural networks problems, we used $\eta_{1}=0.1$, and $\eta_{2}=0.3$ tuned by some experiments.

$L_{1}$ and $L_{2}$ parameters in SANC have tuned. The searching range for $L_{1}$ and $L_{2}$ are $10^{-3\mathrel{\mathop{\mathchar 58\relax}}1\mathrel{\mathop{\mathchar 58\relax}}3}$. The used parameter values are presented in Table 1.

For the Lanczos method, we truncated it at Lanczos iteration $5$ for all experiments. Also, an approximate minimizer of an approximate local cubic model over a Krylov subspace was obtained by the ’conjugate gradient method’ in scipy.optimize.minimize. To calculate the left-most eigenpair of $\mathbf{T}$, we used eigh_tridiagonal in scipy.linalg.

For SANC and SCR, we have to calculate $f(\mathbf{x}_{t})$ and $f(\mathbf{x}_{t}+\mathbf{s}_{t})$ to measure $\rho_{t}$. For the neural networks problems, because of the shortage of the memory storage, we made another set of dataset whose size is same as $\mathcal{S}_{\mathbf{g}}$ and $\mathcal{S}_{\mathbf{B}}$ to estimate $f(\mathbf{x}_{t})$ and $f(\mathbf{x}_{t}+\mathbf{s}_{t})$.

5.2 Setup

In this section, we present some experimental results to show the effectiveness of the SANC algorithm for stochastic nonconvex optimization problems. In our numerical experiments, we considered three machine learning problems, namely,

(i) the logistic regression

(ii) the multi-layer neural networks and

(iii) the convolutional neural networks (CNN)

with real datasets.

For the logistic regression problems, with the binary classification dataset, i.e., $\{\mathbf{x}_{i},y_{i}\}$, $y_{i}\in\{0,1\}$, in order to find the optimum $\mathbf{w}^{*}$, we solved the following problem,

[TABLE]

where $\Omega(\mathbf{w})=\sum_{i=1}^{d}\frac{w_{i}^{2}}{1+w_{i}^{2}}$, and $\lambda$ is a fixed regularization coefficient. We initialized all variables by $\mathbf{w}_{0}=1$.

For multi-layer neural networks, we used two hidden layers. The number of neurons in the first hidden layers is $300$ and the number of neurons in the second hidden layers is $500$. At the outer layer, softmax functions were used. Cross entropy loss was used for obtaining an objective function. We also added $l2$ norm as a convex regularization term with a coefficient $\lambda=0.01$. Also, we used hyperbolic tangent functions as nonlinear activation functions.

For the CNN, we used two convolutional receptive filters. The first one is $5\times 5\times 32$ dimensional and the second one is $5\times 5\times 64$ dimensional. The fully connected layer was added at the output of the convolutional layers, which has $1000$ neurons. Also, we used hyperbolic tangent functions for nonlinear activation functions. We also added $l2$ norm as a convex regularization term with a coefficient $\lambda=0.01$.

All the variables including the weights and the biased in the multi-layer neural networks and the CNN were initialized by the Xavier initialization (Glorot & Bengio, 2010).

We used real datasets from libsvm (Chang & Lin, 2011) for the logistic regression problem and multi-layer neural networks. For CNN, we used MNIST (LeCun et al., 1998) and CIFAR10 datasets (Krizhevsky & Hinton, 2009).

The sizes of $\mathcal{S}_{\mathbf{g}}$ and $\mathcal{S}_{\mathbf{B}}$ are $\lceil\text{the number of datapoints}/20\rceil$ for the logistic regression problems. For the multi-layer neural networks and CNN, we computed stochastic gradient and Hessian estimates using independently drawn mini-batches of size 128. These batch size settings are used for SANC as well as all the baselines which are described below.

We compared our method, SANC, with various optimization methods described in Section 2. These include stochastic gradient descent(SGD), Sub-sampled Cubic Regularization (SCR) (Kohler & Lucchi, 2017), Cubic Regularization (CR) (Nesterov & Polyak, 2006), and Negative Curvature Direction (NCD) method (Liu et al., 2018). It is noted that the only difference between SANC and SCR is to use a direction of negative curvature. The details on the setting of the each algorithm is as follows,

• •

SGD: Fixed step-size is used. We used 0.01 for the logistic regression problems and 0.001 for the multi-layer neural networks and CNN problems, which is the best choice among $10^{-3\mathrel{\mathop{\mathchar 58\relax}}1\mathrel{\mathop{\mathchar 58\relax}}3}$, respectively.

• •

SCR (Kohler & Lucchi, 2017): We used the same parameter values what we used for SANC. We used same sets, $\mathcal{S}_{\mathbf{g}}$, $\mathcal{S}_{\mathbf{B}}$ in SANC, which is different with the experiments in (Kohler & Lucchi, 2017) where they used increasing sized sets.

• •

CR (Nesterov & Polyak, 2006): For the fixed cubic coefficient $\sigma$, we used $5.$ for all problems.

• •

NCD (Liu et al., 2018): We used the same parameter values what we used for SANC.

All baselines are our implementations using Tensorflow.

5.3 Results

Figure 1 shows that the SANC method has advantages to use a direction of negative curvature at unsuccessful iterations. Because of the existence of the nonconvex regularization term, initial iterates have indefinite Hessians. As pointed out in Section 3, with a smaller $\sigma_{0}$, we experienced that the SCR method suffers from unsuccessful iterations in initial iterations.

Figure 2 illustrates the training loss results for three problems. It shows that SANC and SCR are identical for the logistic regression problems. This is because SANC does not make the update with a direction of negative curvature in that problems. As shown in Figure 2(b), 2(c), SANC reveals competitive loss history against others. Because neural networks problems have so many local minima and saddle points, it is suitable to use a negative curvature step to find a better step at unsuccessful iteration. Please refer to the supplementary material for more numerical results.

6 Conclusion

In this work, we have developed a stochastic optimization framework which is called SANC for nonconvex optimization. The novelty of the proposed algorithms is that we combine the adaptive cubic regularized Newton method and the negative curvature method. This method not only maintain the trust-region like framework of the adaptive cubic regularized Newton method but also could reduce the number of oracle calls for the large-scale machine learning problems with real datasets. We would like to point out that a Krylov subspace where an approximate minimizer of an approximate local cubic model is established can also be efficiently used for approximating the left-most eigenvector of the Hessian. This eigenvector can be utilized as a direction of negative curvature. Also, as shown in empirical results, a subsampling technique, that we proposed, is suitable to be used for solving large-scale machine learning problems. As far as we know, this approach is the first work to integrate the negative curvature method into the Newton method. We think the accelerated method and the stochastic variance reduced gradient method can perform well with our algorithm to attain better theoretical and empirical results.

acknowledgements

This work was partially supported by the National Research Foundation (NRF) of Korea. (NRF-2018R1D1A1B07043406).

Appendix A Appendix

For the completeness, we rewrite all the lemmas and the theorems of the article.

A.1 Gradient and Hessian sampling bound

Before addressing Theorem 1 in the article, we introduce the vector Bernstein inequality and its following lemma. These proofs are based on the work from (Kohler & Lucchi, 2017).

Lemma 8** (Vector Bernstein Inequality, Theorem 12 in (Gross, 2011)).**

Let $\mathbf{X}_{1},\dots\mathbf{X}_{m}$ be independent zero-mean vector-valued random variables. Let

[TABLE]

then

[TABLE]

where $V=\sum_{i}\mathbb{E}\left[\lVert\mathbf{X}_{i}\rVert_{2}^{2}\right]$ and $t\leq V/max_{i}\lVert\mathbf{X}_{i}\rVert_{2}$.

Proof.

please refer to the proof of the theorem 12 in (Gross, 2011) ∎

Before describing the proof for Theorem 1 in the paper, we describe the gradient bound in terms of the sampling size to establish $\epsilon_{\mathbf{g}}$ error.

Lemma 9** (Approximate Gradient Deviation Bound).**

For $\epsilon_{\mathbf{g}}\leq 2L_{0}\left(1+\frac{1}{\sqrt{|\mathcal{S}_{\mathbf{g}}|}}\right)$, if

[TABLE]

then for small $\epsilon_{\mathbf{g}}>0$ we can establish the condition 16a with high probability $1-\delta$

Proof.

We use Lemma 8 for the vector Bernstein inequality. Before using it, first, we make independent zero-mean vector-valued random variables,

[TABLE]

Then,

[TABLE]

where the first inequality is from the triangular inequality and the second inequality is from the Lipschitz continuity of the function values. Hence,

[TABLE]

From Eq. 35, let $|\mathcal{S}_{\mathbf{g}}|\epsilon_{\mathbf{g}}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{V}+t$ then,

[TABLE]

Thus, the probability of $\epsilon_{\mathbf{g}}<\frac{\left(1+\sqrt{4\log(1/\delta)}\right)2L_{0}}{\sqrt{|\mathcal{S}_{\mathbf{g}}}|}$ is higher or equal to $1-\delta$. Rearranging this condition yields the conclusion. ∎

Theorem 1** (Gradient Sampling Bound).**

If

[TABLE]

then $\mathbf{g}_{t}$ satisfies the gradient approximation bound Eq. 16a for all $t\geq 0$ with high probability $1-\delta$.

Proof.

From Lemma 9,

[TABLE]

Squaring both sides finalizes the proof. ∎

For the Hessian sampling size bound, we use a similar approach. The proof is based on the Operator Bernstein Inequality. The big difference between the gradient sampling and Hessian sampling is that Hessian sampling size depends on the number of variables, $d$.

Lemma 10** (Operator-Bernstein Inequality, Theorem 6 in (Gross, 2011)).**

Let $\mathbf{X}_{i},\,i=1,\dots,m$ be i.i.d., $n$ by $n$, zero-mean, Hermitian matrix-valued random variables. Assume $V_{0},c\in\mathbb{R}$ are such that $\lVert\mathbb{E}[\mathbf{X}_{i}^{2}]\rVert\leq V_{0}^{2}$ and $\lVert\mathbf{X}_{i}\rVert\leq c$. Set $S=\sum_{i=1}^{m}\mathbf{X}_{i}$ and let $V=mV_{0}^{2}$. Then,

[TABLE]

for $t\leq 2V/c$.

Proof.

Please refer to the proof of the theorem 6 in (Gross, 2011). ∎

Let us denote $\mathbf{H}_{i}$ as an abbreviation of $\nabla^{2}f_{i}(\mathbf{x})$ for some $i$ and $\mathbf{H}$ denote the expectation of $\mathbf{H}_{i}$.

Lemma 11** (Approximate Hessian Deviation Bound).**

For $\epsilon_{\mathbf{B}}\leq 4L_{1}$, if

[TABLE]

then for small $\epsilon_{\mathbf{B}}>0$ we can establish the condition 16b with high probability $1-\delta$.

Proof.

To use the operator-Bernstein inequality 48, We introduce an independent zero mean Hermitian random variable,

[TABLE]

Then,

[TABLE]

where the inequality is from the triangular inequality and the fact that $f_{i}$ has Lipschitz continuous gradient. Let define $V_{0}^{2}\mathrel{\mathop{\mathchar 58\relax}}=4L_{1}^{2}$ and $V\mathrel{\mathop{\mathchar 58\relax}}=|\mathcal{S}_{\mathbf{B}}|V_{0}^{2}=|\mathcal{S}_{\mathbf{B}}|4L_{1}^{2}$ Then,

[TABLE]

Hence from the operator-Bernstein inequality 48,

[TABLE]

for $|\mathcal{S}_{\mathbf{B}}|\epsilon_{\mathbf{B}}\leq\frac{4|\mathcal{S}_{\mathbf{B}}|L_{1}^{2}}{L_{1}}$ that means $\epsilon_{\mathbf{B}}\leq 4L_{1}$. Thus,

[TABLE]

∎

From this result, we can prove the Hessian sampling bound.

Theorem 2** (Hessian Sampling Bound).**

If

[TABLE]

then $\mathbf{B}_{t}$ satisfies the Hessian approximation bound Eq. 16b for all $t\geq 0$ with high probability $1-\delta$.

Proof.

From the previous lemma 11,

[TABLE]

∎

A.2 Worst-case iteration complexity bound for achieving approximate first- and second-order optimality

Lemma 3** (Cubic Regularization Coefficient Bound).**

Let the A2, A3 and A5 hold. For finite positive values $\sigma_{\min}$ and $\sigma_{\max}$, it holds,

[TABLE]

where $\sigma_{\max}\mathrel{\mathop{\mathchar 58\relax}}=\max\left\{\sigma_{0},\gamma\left(\frac{3}{2}L_{2}+\frac{3(\epsilon_{\mathbf{g}}+\frac{1}{2}\epsilon_{\mathbf{B}})}{\epsilon_{s}^{2}}\right)\right\}$.

Proof.

For $\sigma_{\min}$, $\sigma_{t}$ is always bounded below by $\epsilon_{m}$, the machine precision. Further, we need to show that $\sigma_{t}\leq\sigma_{\max}$. For each iteration $t\geq 0$,

[TABLE]

where $\xi_{t}$ is a line segment of $\mathbf{x}_{t}$ and $\mathbf{x}_{t}+\mathbf{s}_{t}$.

To achieve $f(\mathbf{x}_{t}+\mathbf{s}_{t})\leq\widetilde{m}_{t}(\mathbf{s}_{t})$ from Eq. 65, we need to have the condition that

[TABLE]

Then, from the definition of $\rho_{t}$, the iteration $t$ is very successful and $\sigma_{t+1}\leq\sigma_{t}$. Let us think that $\sigma_{t}$ is slightly less than $\frac{3}{2}L_{2}+\frac{3(\epsilon_{\mathbf{g}}+\frac{1}{2}\epsilon_{\mathbf{B}})}{\epsilon_{s}^{2}}$, then it could be multiplied by $\gamma$ and the case iterations can be started with a large $\sigma_{0}$. This gives the desired result. ∎

For the two consecutive iterations, we have to establish a sufficient decrease. Hence we provide the following lemma for the negative curvature step.

Lemma 4** (Sufficient Decrease with Negative Curvature Step, Lemma 3 in (Liu et al., 2018)).**

When $\mathbf{v}_{t}^{T}\mathbf{B}_{t}\mathbf{v}_{t}\leq-\epsilon^{\prime}/2$ and $\epsilon_{\mathbf{B}}\leq\epsilon^{\prime}/12$, the negative curvature step $\mathbf{d}_{t}$ satisfies that

[TABLE]

with high probability $1-\delta$.

Proof.

Please refer to the proof of the Lemma 3 in (Liu et al., 2018). ∎

From the paper, we showed that the approximate minimizer $s_{t}$ over the Krylov subspace $\mathcal{K}_{j}$ satisfies Eq. 11a and Eq. 11b. The next lemma gives a lower bound on the approximate local cubic model decrease when Eq. 11a and Eq. 11b are satisfied for a successful iteration $t\in\mathcal{N}$.

Lemma 5** (Local Cubic Model Decrease, Lemma 3.3 in (Cartis et al., 2011a)).**

Suppose that $\mathbf{s}_{t}$ satisfies Eq. 11a and Eq. 11b. Then for all (very) successful iterations $t$,

[TABLE]

Proof.

Please refer to the proof of the Lemma 3.3 in (Cartis et al., 2011a). Please note that in (Cartis et al., 2011a) they established an approximate local cubic model with exact gradient $\nabla f$ and approximate Hessian $\mathbf{B}$, but in our approach we establish an approximate local cubic model with approximate gradient $\mathbf{g}$ and Hessian $\mathbf{B}$. However, the procedure of the proof is same except for using $\mathbf{g}$ instead of $\nabla f$. ∎

Also, we establish the relationship between the magnitude of $\mathbf{s}_{t}$ and $\mathbf{g}_{t}$ for a successful iteration $t$.

Lemma 6** (Lower Bound of $\mathbf{s}_{t}$).**

Suppose that A2, A3, A4, and A5 hold. For all (very) successful iterations $t$, the Newton step $\mathbf{s}_{t}$ satisfies that

[TABLE]

where

[TABLE]

Proof.

For a successful iteration $t$,

[TABLE]

where $\xi_{t}$ is a line segment between $\mathbf{x}_{t}$ and $\mathbf{x}_{t}+\mathbf{s}_{t}$. First inequality is from triangular inequality and second inequality is from Assumption 3 and Eq. 21.

Hence with the TC 4,

[TABLE]

∎

Lemma 12**.**

With real numbers $a>0$, $b>0$, $t\in(0,1)$, $-a+\sqrt{a^{2}+b^{2}}\geq tb$ holds if and only if $b\geq\frac{2t}{\sqrt{1-t^{2}}}a$.

Proof.

[TABLE]

∎

Theorem 7** (Worst-Case Iteration Complexity for Approximate First- and Second-Order Optimality).**

Let A1, A2, A3, and A4 hold. Algorithm 1 provides an iteration $t^{*}$ such that $\lVert\nabla f(\mathbf{x}_{t^{*}})\rVert\leq\epsilon$ and $-\lambda_{\min}(\mathbf{Q}_{t^{*}}^{T}\mathbf{B}_{t^{*}}\mathbf{Q}_{t^{*}})\leq\epsilon$ within at most

[TABLE]

iterations with high probability $1-\delta$ for

[TABLE]

where

[TABLE]

Proof.

1. (a)

First, we need to show that an first order stationary point exists within at most the certain number of iterations.

For a sufficiently large iteration $j$, suppose that $\lVert\nabla f(\mathbf{x}_{t+1})\rVert>\epsilon,\,\forall t\leq j$. From lemma 12, plugging $\epsilon_{\mathbf{B}}$, $\sqrt{\frac{2L_{2}+4\sigma_{\max}}{3}\epsilon}$ and $1/2$ into $a$, $b$ and $t$, respectively, yields

[TABLE]

if and only if $\epsilon\leq\frac{2}{\kappa_{l}}\epsilon_{\mathbf{B}}$ holds.

From the definition of $\epsilon$ Eq. 84, we assume that $\lVert\nabla f(\mathbf{x}_{t}+\mathbf{s}_{t})\rVert>\epsilon$, then

[TABLE]

Hence, from Eq. 77 and Eq. 89,

[TABLE]

For a successful iterate $\mathbf{x}_{t}$, plugging Eq. 92 into Eq. 23 and using Eq. 5 yield,

[TABLE]

From Lemma 4 with $\epsilon^{\prime}=\epsilon^{\frac{1}{2}}$ for an iteration $t\in\mathcal{U}_{j}$,

[TABLE]

That means

[TABLE]

with high probability 1-$\delta$.

Therefore, we can attain an upper bound on $j$,

[TABLE] 2. (b)

We also have to attain an upper bound for the second-order optimality as the same approach in a. For a sufficiently large iteration $j$, suppose that $-\lambda_{\min}(\mathbf{Q}_{t}^{T}\mathbf{B}_{t}\mathbf{Q}_{t})>\epsilon,\,\forall t\leq j$. Eq. 11b implies that the matrix $\mathbf{Q}_{t}^{T}\mathbf{B}_{t}\mathbf{Q}_{t}+\sigma_{t}\lVert\mathbf{s}_{t}\rVert\mathbf{I}$ is positive semidefinite which gives

[TABLE]

because $-\lambda_{\min}(\mathbf{Q}_{t}^{T}\mathbf{B}_{t}\mathbf{Q}_{t})>\epsilon,\,\forall t\leq j$.

Hence, given an iteration $t\leq j$, from Eq. 23 of Lemma 5 and Eq. 5 of $\rho_{t}$

[TABLE]

Summing up to iteration $j$ gives,

[TABLE]

Thus, the upper bound on $j$ is

[TABLE] 3. (c)

To get an upper bound on all the possible iterations that occur either with $\lVert\nabla f(\mathbf{x}_{t^{*}})\rVert>\epsilon$ or $-\lambda_{\min}(\mathbf{Q}_{t^{*}}^{T}\mathbf{B}_{t^{*}}\mathbf{Q}_{t^{*}})>\epsilon$ , the sum of the bounds Eq. 100 and Eq. 109 gives,

[TABLE]

That equivalently means that if $j$ exceeds $l_{1}+l_{2}+1$, there exists at least an iterations $t$ which necessarily satisfy the condition $\lVert\nabla f(\mathbf{x}_{t^{*}})\rVert\leq\epsilon$ and $-\lambda_{\min}(\mathbf{Q}_{t^{*}}^{T}\mathbf{B}_{t^{*}}\mathbf{Q}_{t^{*}})\leq\epsilon$.

∎

More numerical results