Newton-Type Methods for Non-Convex Optimization Under Inexact Hessian Information

Introduction

Consider the generic unconstrained optimization problem

In doing so, we first consider (P0), and study the theoretical convergence properties of variants of these two algorithms in which, under favorable conditions, Hessian is suitably approximated. We show that our Hessian approximation conditions, in many cases, are weaker than the existing ones in the literature. In addition, and in contrast to some prior works, our conditions allow for efficient constructions of the inexact Hessian with a priori guarantees via various approximation methods, of which Randomized Numerical Linear Algebra (RandNLA), , techniques are shown to be highly effective.

Subsequently, to showcase the application of randomized techniques for construction of the approximate Hessian, we consider an important instance of (P0), i.e., large-scale finite-sum minimization, of the form

The rest of this paper is organized as follows. In Section 1.1, we first introduce the notation and definitions used throughout the paper. For completeness, in Section 1.2, we give a brief review of trust region (Section 1.2.1) and cubic regularization (Section 1.2.2) along with related prior works. Our main contributions are summarized in Section 1.3. Theoretical analysis of the proposed algorithms for solving generic non-convex problem (P0) are presented in Section 2. Various randomized sub-sampling strategies as well as theoretical properties of the proposed algorithms for finite-sum minimization problems (P1) and (P2) are given in Section 3. Conclusions and further thoughts are gathered in Section 4.

Unlike convex functions for which “local optimality” and “global optimality” are in fact the same, in non-convex settings, we are often left with designing algorithms that can guarantee convergence to approximate local optimality. In this light, throughout this paper, we make use of the following definition of $(\epsilon_{g},\epsilon_{H})$-Optimality:

We note that $(\epsilon_{g},\epsilon_{H})$-Optimality (even with $\epsilon_{g}=\epsilon_{H}=0$) does not necessarily imply closeness to any local minimum, neither in iterate nor in the objective value. However, if the saddle points satisfy the strict-saddle property , then an $(\epsilon_{g},\epsilon_{H})$-optimality guarantees vicinity to a local minimum for sufficiently small $\epsilon_{g}$ and $\epsilon_{H}$.

2 Background and Related Work

Arguably, the most straightforward approach for globalization of many Newton-type algorithms is the application of line-search. However, near saddle points where the gradient magnitude can be small, traditional line search methods can be very ineffective and in fact produce iterates that can get stuck at a saddle point . Trust region and cubic regularization methods are two elegant globalization alternatives that, specially recently, have attracted much attention. The main advantage of these methods is that they are reliably able to take advantage of the direction of negative curvature and escape saddle points. In this section we briefly review these algorithms as they pertain to the present paper and mention the relevant prior works.

TR methods encompass a general class of iterative methods which specifically define a region around the current iterate within which they trust the model to be a reasonable approximation of the true objective function. The most widely used approximating model, which we consider here, is done via a quadratic function. More specifically, using the current iterate $\mathbf{x}_{t}$, the quadratic variant of TR algorithm finds the next iterate as $\mathbf{x}_{t+1}=\mathbf{x}_{t}+\mathbf{s}_{t}$ where $\mathbf{s}_{t}$ is a solution of the constrained sub-problem

For a smooth non-convex objective and in order to obtain approximate first-order criticality, i.e., $\|\nabla F(\mathbf{x}_{t})\|\leq\epsilon_{g}$ for some $\epsilon_{g}\in(0,1)$, the complexity of an (inexact) trust-region method, which ensures at least a Cauchy (steepest-descent-like) decrease at each iteration, is shown to be of the same order as that of steepest descent, i.e., $\mathcal{O}(\epsilon_{g}^{-2})$; e.g., . Recent non-trivial modifications of the classical TR methods have also been proposed which improve upon the complexity to $\mathcal{O}(\epsilon_{g}^{-3/2})$; see and further extensions to a more general framework in . These bounds can be shown to be tight in the worst case. Under a more general algorithmic framework and in terms of objective function sub-optimality, i.e., $F(\mathbf{x})-F^{*}\leq\epsilon$, better complexity bounds, in the convex and strongly-convex settings, have been obtained which are of the orders of $\mathcal{O}(\epsilon_{g}^{-1})$ and $\mathcal{O}(\log({1}/{\epsilon_{g}}))$, respectively .

For non-convex problems, however, it is more desired to obtain complexity bounds for achieving approximate second-order criticality, i.e., Definition 1. For this, bounds in the orders of $\mathcal{O}(\max\{\epsilon_{H}^{-1}\epsilon_{g}^{-2},\epsilon_{H}^{-3}\})$ and $\mathcal{O}(\max\{\epsilon_{g}^{-3},\epsilon_{H}^{-3}\})$ have been obtained in and , respectively. Similar bounds were also given in under probabilistic model. Bounds of this order have shown to be optimal in certain cases .

More closely related to the present paper, there have been several results which study the role of derivative-free and probabilistic models in general, and Hessian approximation in particular, e.g., see and references therein.

2.2 Cubic Regularization

An alternative to the traditional line-search and TR for globalization of Newton-type methods is the application of cubic regularization. Such class of methods is characterized by generating iterates as $\mathbf{x}_{t+1}=\mathbf{x}_{t}+\mathbf{s}_{t}$ where $\mathbf{s}_{t}$ is a solution of the following unconstrained sub-problem

where $\sigma_{t}$ is the cubic regularization parameter chosen for the current iteration. As in the case of TR, the major bottleneck of CR involved solving the sub-problem (2b), for which various techniques have been proposed, e.g., .

To the best of our knowledge, the use of such regularization, was first introduced in the pioneering work of , and subsequently further studied in the seminal works of .From the worst-case complexity point of view, CR has a better dependence on $\epsilon_{g}$ compared to TR. More specifically, showed that, under global Lipschitz continuity assumption on the Hessian, if the sub-problem (2b) is solved exactly, then the resulting CR algorithm achieves the approximate first-order criticality with complexity of $\mathcal{O}(\epsilon_{g}^{-3/2})$. These results were extended by the pioneering and seminal works of to an adaptive variant, which is often referred to as ARC (Adaptive Regularization with Cubics). In particular, the authors showed that the worst case complexity of $\mathcal{O}(\epsilon_{g}^{-3/2})$ can be achieved without requiring the knowledge of the Hessian’s Lipschitz constant, access to the exact Hessian, or multi-dimensional global optimization of the sub-problem (2b). These results were further refined in where it was shown that, not only, multi-dimensional global minimization of (2b) is unnecessary, but also the same complexity can be achieved with mere one or two dimensional search. This $\mathcal{O}(\epsilon^{-3/2})$ bound has been shown to be tight . As for the approximate second-order criticality, showed that at least $\mathcal{O}(\max\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\})$ is required. With further assumptions on the inexactness of sub-problem solution, also show that one can achieve $\mathcal{O}(\max\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\})$, which is shown to be tight . Better dependence on $\epsilon_{g}$ can be obtained if one assumes additional structure, such as convexity, e.g., see as well as the acceleration scheme of .

Recently, for (strongly) convex problems, obtained sub-optimal complexity for ARC and its accelerated variants using Hessian approximations. In the context of stochastic optimization problems, considers cubic regularization with a priori chosen fixed regularization parameter using both approximations of the gradients and Hessian. Specific to the finite-sum problem (P1), and by a direct application of the theoretical results of , presents a sub-sampled variant of ARC, in which the exact Hessian and the gradient are replaced by sub-samples. However, unfortunately, their analysis suffers from a rather vicious circle: the approximate Hessian and gradient are formed based on an a priori unknown step which can only be determined after such approximations are formed.

3 Contributions

In this section, we summarize the key aspects of our contributions. In Section 2, we consider (P0) and establish the worst-case iteration complexities for variants of trust-region and adaptive cubic regularization methods in which the Hessian is suitably approximated. More specifically, our entire analysis is based on the following key condition on the approximate Hessian $\mathbf{H}(\mathbf{x})$:

For some $0<K_{H}<\infty$, $\epsilon>0$, the approximating Hessian, $\mathbf{H}(\mathbf{x}_{t})$, satisfies

where $\mathbf{x}_{t}$ and $\mathbf{s}_{t}$ are, respectively, the iterate and the update at iteration $t$.

Under Condition 1, we show that our proposed algorithms (Algorithms 1 and 2) achieve the same worst-case iteration complexity to obtain approximate second order critical solution as that of the exact variants (Theorems 1, 2, and 3).

In Section 3, we describe schemes for constructing $\mathbf{H}(\mathbf{x}_{t})$ to satisfy Condition 1. Specifically, in the context of finite-sum optimization framework, i.e., problems (P1) and (P2), we present various sub-sampling schemes to probabilistically ensure Condition 1 (Lemmas 16 and 17). Our proposed randomized sub-sampling strategies guarantee, with high probability, a stronger condition than (3a), namely

It is clear that (4) implies (3a). We then give optimal iteration complexities for Algorithms 1 and 2 for optimization of non-convex finite-sum problems where the Hessian is approximated by means of appropriate sub-sampling (Theorems 4, 5 and 6).

To establish optimal second-order iteration complexity, many previous works considered Hessian approximation conditions that, while enjoying many advantages, come with certain disadvantages. Our proposed Condition 1 aims to remedy some of these disadvantages. We first briefly review the conditions used in the prior works, and subsequently highlight the merits of Condition 1 in comparison.

For the analysis of trust-region, many authors have considered the following condition

For cubic regularization, the condition imposed on the inexact Hessian is often considered as

for some $0<C_{3}<\infty$, e.g., and other follow-up works. In fact, has also established optimal iteration complexity for trust-region algorithm under (5c). Both of (5a) and (5c), are stronger than the celebrated Dennis-Moré condition, i.e.,

Indeed, under certain assumptions, Dennis-Moré condition is satisfied by a number of quasi-Newton methods, although the same cannot be said about (5a) and (5c) .

3.2 Merits of Condition 1

For our trust-region analysis, we require Condition 1 with $\epsilon\in\operatorname{\mathcal{O}}(\max\left\{\epsilon_{H},\Delta_{t}\right\})$; see (11) in Theorem 1. Hence, when $\Delta_{t}$ is large, e.g., at the beginning of iterations, all the conditions (3a), (5a), and (5b) are equivalent, up to some constants. However, the constants in (5a) and (5b) can be larger than what is implied by (3a), amounting to cruder approximations in practice for when $\Delta_{t}$ is large. As iterations progress, the trust-region radius will get smaller, and in fact it is expected that $\Delta_{t}$ will eventually shrink to be $\Delta_{t}\in\Theta\left(\min\{\epsilon_{g},\epsilon_{H}\}\right)$. In prior works, e.g., , the convergence analysis is derived using $\epsilon_{H}=\epsilon_{g}$, whereas here we allow $\epsilon_{H}=\sqrt{\epsilon_{g}}$. As a result, the requirements (5a) and (5b) can eventually amount to stricter conditions than (3a).

As for (5c), the main drawback lies in the difficulty of enforcing it. Despite the fact that for certain values of $\|\mathbf{s}_{t}\|$ and $\epsilon$, e.g., $\epsilon\ll\|\mathbf{s}_{t}\|$, (5c) can be less restrictive than (3a), a priori enforcing (5c) requires one to have already computed the search direction $\mathbf{s}_{t}$, which itself can be done only after $\mathbf{H}(\mathbf{x}_{t})$ is constructed, hence creating a vicious circle. A posteriori guarantees can be given if one obtains a lower-bound estimate on the yet-to-be-computed step-size, i.e., to have $s_{0}>0$ such that $s_{0}\leq\|\mathbf{s}_{t}\|$. This allows one to consider a stronger condition as ${\left\|\left(\mathbf{H}(\mathbf{x}_{t})-\nabla^{2}F(\mathbf{x}_{t})\right)\right\|}\leq C_{3}s_{0}$, which can be enforced using a variety of methods such as those described in Section 3. However, to obtain such a lower-bound estimate on the next step-size, one has to resort to a recursive procedure, which necessitates repeated constructions of the approximate Hessian and subsequent solutions of the corresponding subproblems. Consequently, this procedure may result in a significant computational overhead and will lead to undesirable theoretical complexities.

In sharp contrast to (5c), the condition (3a) allows for theoretically principled use of many practical techniques to construct $\mathbf{H}_{t}$. For example, under (3a), the use of quasi-Newton methods to approximate the Hessian is theoretically justified. Further, by considering the stronger condition (4), many randomized matrix approximation techniques can be readily applied, e.g., ; see Section 3. To the best of our knowledge, the only successful attempt at guaranteed a priori construction of $\mathbf{H}_{t}$ using (5c) is done in . Specifically, by considering probabilistic models, which are “sufficiently accurate” in that they are partly based on (5c), studies first-order complexity of a large class of methods, including ARC, and discusses ways to construct such probabilistic models as long as the gradient is large enough, i.e., before first-order approximate-optimality is achieved. Here, by considering (3a), we are able to provide an alternative analysis, which allows us to obtain second-order complexity results.

Requiring (4), as a way of enforcing (3a), offers a variety of other practical advantages, which are not readily available with other conditions. For example, consider distributed/parallel environments where the data is distributed across a network and the main bottleneck of computations is the communications across the nodes. In such settings, since (4) allows for the Hessian accuracy to be set a priori and to remain fixed across all iterations, the number of samples in each node can stay the same throughout iterations. This prevents unnecessary communications to re-distribute the data at every iteration.

Algorithms and Convergence Analysis

We are now ready to present our main algorithms for solving the generic non-convex optimization (P0) along with their corresponding iteration complexity results to obtain a $(\epsilon_{g},\epsilon_{H})$-optimal solution as in (1). More precisely, in Section 2.1 and 2.2, respectively, we present modifications of the TR and ARC methods which incorporate inexact Hessian information, according to Condition 1.

We remind that, though not specifically mentioned in the statement of the theorems or the algorithms, when the computed steps are rejected and an iteration needs to be repeated with different $\Delta_{t}$ or $\sigma_{t}$, the previous $\mathbf{H}_{t}$ may seamlessly be used in the next iteration. This can be a desirable feature in many practical situations and is directly the result of enforcing (4); see also the discussion in Section 1.3.2.

For our analysis throughout the paper, we make the following standard assumption regarding the regularity of the exact Hessian of the objective function $F$.

$F(\mathbf{x})$ is twice differentiable and has bounded and Lipschitz continuous Hessian on the piece-wise linear path generated by the iterates, i.e. for some $0<K,L<\infty$ and all iterations

where $\mathbf{x}_{t}$ and $\mathbf{s}_{t}$ are, respectively, the iterate and the update step at iteration $t$.

Although, we do not know of a particular way to, a priori, verify (6a), it is clear that Assumption (6a) is weaker than Lipschitz continuity of the Hessian for all $\mathbf{x}$, i.e.,

Despite the fact that theoretically (6a) is weaker than (7), to the best of our knowledge as of yet, (7) is the only practical sufficient condition for verifying (6a).

Algorithm 1 depicts a trust-region algorithm where at each iteration $t$, instead of the true Hessian $\nabla^{2}F(\mathbf{x}_{t})$, only an inexact approximation, $\mathbf{H}_{t}$, is used. For Algorithm 1, the accuracy tolerance in (3a) is adaptively chosen as $\epsilon_{t}\leq\max\left\{\epsilon_{0},\Delta_{t}\right\}$, where $\Delta_{t}$ is the trust region in the t-th iteration and $\epsilon_{0}\in\mathcal{O}(\epsilon_{H})$ is some fixed threshold. This allows for a very crude approximation at the beginning of iterations, when $\Delta_{t}$ is large. As iterations progress towards optimality and $\Delta_{t}$ gets small, the threshold $\epsilon_{0}$ can prevent $\epsilon$ from getting unnecessarily too small.

In Algorithm 1, we require that the sub-problem (8) is solved only approximately. Indeed, in large-scale problems, where the exact solution of the sub-problem is the main bottleneck of the computations, this is a very crucial relaxation. Such approximate solution of the sub-problem (8) has been adopted in many previous work. Here, we follow the inexactness conditions discussed in , which are widely known as Cauchy and Eigenpoint conditions. Recall that the Cauchy and Eigen directions correspond, respectively, to one dimensional minimization of the sub-problem (8) along the directions given by the gradient and negative curvature.

Assume that we solve the sub-problem (8) approximately to find $\mathbf{s}_{t}$ such that

Here, $m_{t}(\cdot)$ is defined in (8), $\mathbf{s}_{t}^{C}$ (Cauchy point) is along negative gradient direction and $\mathbf{s}_{t}^{E}$ is along approximate negative curvature direction such that $\langle\mathbf{s}_{t}^{E},\mathbf{H}_{t}\mathbf{s}_{t}^{E}\rangle\leq\nu\lambda_{\min}(\mathbf{H}_{t})\|\mathbf{s}_{t}^{E}\|^{2}<0$, for some $\nu\in(0,1]$ (see Appendix B for a way to efficiently compute $\mathbf{s}_{t}^{E}$).

One way to ensure that an approximate solution to the sub-problem (8) satisfies (9), is by replacing (8) with the following reduced-dimension problem, in which the search space is a two-dimensional sub-space containing vectors $\mathbf{s}_{t}^{C}$, and $\mathbf{s}_{t}^{E}$, i.e.,

Of course, any larger dimensional sub-space $\mathcal{P}$ for which we have $\text{Span}\{\mathbf{s}_{t}^{C},\mathbf{s}_{t}^{E}\}\subseteq\mathcal{P}$ would also guarantee (9). In fact, a larger dimensional sub-space implies a more accurate solution to our original sub-problem (8).

We now set out to provide iteration complexity for Algorithm 1. Our analysis follows similar line of reasoning as that in . First, we show the discrepancy between the quadratic model and objective function in Lemma 1.

Given Assumption 1 and Condition (3a) with any $\epsilon_{t}>0$, we have

Applying Mean Value Theorem on $F$ at $\mathbf{x}_{t}$ gives $F(\mathbf{x}_{t}+\mathbf{s}_{t})=F(\mathbf{x}_{t})+\nabla F(\mathbf{x}_{t})^{T}\mathbf{s}_{t}+\frac{1}{2}\mathbf{s}_{t}^{T}\nabla^{2}F(\xi_{t})\mathbf{s}_{t}$, for some $\xi_{t}$ in the segment of $[\mathbf{x}_{t},\mathbf{x}_{t}+\mathbf{s}_{t}]$. We have

Combining with Conditions 1 and 2, we get the following two lemmas that characterize sufficient conditions for successful iterations.

Consider any $\epsilon_{H}>0$, let $\epsilon_{0}\triangleq\alpha(1-\eta)\nu\epsilon_{H}$ for some $\alpha\in(0,1)$, and suppose Condition 1 is satisfied with $\epsilon_{t}\leq\max\{\epsilon_{0},\Delta_{t}\}$, where $\Delta_{t}$ is the trust region at the t-th iteration. Given Assumption 1 and Condition 2, if $\lambda_{\min}(\mathbf{H}_{t})<-\epsilon_{H}$ and $\Delta_{t}\leq(1-\alpha)(1-\eta)\nu\left|\lambda_{\min}(\mathbf{H}_{t})\right|/(L+1)$, then the t-th iteration is successful, i.e. $\Delta_{t+1}=\gamma\Delta_{t}$.

Suppose $\Delta_{t}\leq\epsilon_{0}$. From (9b) and (10), we have

By the assumption on $\Delta_{t}$, we get $\rho_{t}\geq\eta$ and the iteration is successful. Now consider $\Delta_{t}\geq\epsilon_{0}$. Similar to the above, we have

which again by assumption on $\Delta_{t}$ and noting $\alpha<1$, we get $\rho_{t}\geq\eta$ and the iteration is successful. ∎∎

Suppose Condition 1 is satisfied with any $\epsilon_{t}>0$. Given Assumption 1 and Condition 2, if $\|\nabla F(\mathbf{x}_{t})\|>\epsilon_{g}$ and

then, the t-th iteration is successful, i.e. $\Delta_{t+1}=\gamma\Delta_{t}$.

By assumption on $\Delta_{t}$, (9a), and since $\|\nabla F(\mathbf{x}_{t})\|>\epsilon_{g}$, we have

where the last inequality follows by assumption on $\Delta_{t}$. So $\rho_{t}\geq\eta$, which means the iteration is successful. ∎∎

Lemma 4 gives a lower bound for the trust region radius before the algorithm terminates, i.e., this ensures that the trust region never shrinks to become too small.

Consider any $\epsilon_{g},\epsilon_{H}>0$ such that $\epsilon_{H}\leq\sqrt{\epsilon_{g}}$ and let $\epsilon_{0}\triangleq\alpha(1-\eta)\nu\epsilon_{H}$ for some $\alpha\in(0,1)$. Further, suppose Condition 1 is satisfied with $\epsilon_{t}\leq\max\{\epsilon_{0},\Delta_{t}\}$, where $\Delta_{t}$ is the trust region at the t-th iteration. For Algorithm 1, under Assumption 1 and Condition 2, we have $\Delta_{t}\geq\kappa_{\Delta}\min\{\epsilon_{g},\epsilon_{H}\},\forall t\geq 0$, where

We prove by contradiction. Assume that the t-th iteration is the first unsuccessful iteration such that $\Delta_{t+1}=\Delta_{t}/\gamma\leq\kappa_{\Delta}\min\{\epsilon_{g},\epsilon_{H}\}$, i.e., we have

Suppose $\lambda_{\min}(\mathbf{H}_{t})<-\epsilon_{H}$. By Lemma 2, since $\Delta_{t}\leq(1-\alpha)(1-\eta)\nu\left|\lambda_{\min}(\mathbf{H}_{t})\right|/{(L+1)}$, iteration $t$ must have been accepted and we must have $\Delta_{t+1}=\gamma\Delta_{t}>\Delta_{t}$, which is a contradiction. Now suppose $\|\nabla F(\mathbf{x}_{t})\|\geq\epsilon_{g}$. By assumption on $\Delta_{t}$, we have that $\Delta_{t}\leq\kappa_{2}\epsilon_{H}=\epsilon_{0}$, which implies that $\epsilon_{t}\leq\epsilon_{0}$. Since the function $h(a,b)\triangleq{-a+\sqrt{a^{2}+b}}$, for any fixed $b>0$, is decreasing in $a$, and for any fixed $a$, is increasing in $b\geq 0$, we have

As a result, since $\Delta_{t}\leq\min\{\kappa_{3}\epsilon_{g},\kappa_{4}\epsilon_{H}\}$, it must satisfy the condition of Lemma 3. This implies that iteration $t$ must have been accepted, which is a contradiction. ∎∎

The following lemma follows closely the line of reasoning in [11, Lemma 4.5].

Consider any $\epsilon_{g},\epsilon_{H}>0$ such that $\epsilon_{H}\leq\sqrt{\epsilon_{g}}$ and let $\epsilon_{0}\triangleq\alpha(1-\eta)\nu\epsilon_{H}$ for some $\alpha\in(0,1)$. Further, suppose Condition 1 is satisfied with $\epsilon_{t}\leq\max\{\epsilon_{0},\Delta_{t}\}$, where $\Delta_{t}$ is the trust region at the t-th iteration. Let $\mathcal{T}_{succ}$ denote the set of all the successful iterations before Algorithm 1 stops. Then, under Assumption 1, Condition 2, the number of successful iterations is upper bounded by,

where $\widehat{\kappa}_{\Delta}\triangleq\kappa_{\Delta}/2,\widetilde{\kappa}_{\Delta}\triangleq\nu\kappa_{\Delta}^{2}/2$, and $\kappa_{\Delta}$ is as defined in Lemma 4.

Suppose Algorithm 1 doesn’t terminate at the t-th iteration. Then we have either ${\left\|\nabla F(\mathbf{x}_{t})\right\|}\geq\epsilon_{g}$ or $\lambda_{\min}(\Delta^{2}F(\mathbf{x}_{t}))\leq-\epsilon_{H}$. In the first case, from (9a), we have

where $\kappa_{\Delta}$ is as defined in Lemma 4. Similarly, in the second case, from (9b), we obtain

Since $F(\mathbf{x}_{t})$ is monotonically decreasing, we have

Hence, we have $\left|\mathcal{T}_{succ}\right|\leq(F(\mathbf{x}_{0})-F_{\min})\max\{\epsilon_{g}^{-2}\epsilon_{H}^{-1},\epsilon_{H}^{-3}\}/(\eta\min\left\{\widehat{\kappa}_{\Delta},\widetilde{\kappa}_{\Delta}\right\})$. ∎∎

Now we are ready to present the final complexity in Theorem 1.

Consider any $\epsilon_{g},\epsilon_{H}>0$ such that $\epsilon_{H}\leq\sqrt{\epsilon_{g}}$ and let $\epsilon_{0}\triangleq\alpha(1-\eta)\nu\epsilon_{H}$ for some $\alpha\in(0,1)$ where $\eta$ is a hyper-parameter in Algorithm 1, and $\nu$ is as in (9b). Suppose the inexact Hessian, $\mathbf{H}(\mathbf{x})$, satisfies Condition 1 with the approximation tolerance, $\epsilon_{t}$, in (3a) as

where $\Delta_{t}$ is the trust region at the t-th iteration. For Problem (P0), under Assumption 1 and Condition 2, Algorithm 1 terminates after at most $T\in\operatorname{\mathcal{O}}\left(\max\{\epsilon_{g}^{-2}\epsilon_{H}^{-1},\epsilon_{H}^{-3}\}\right)$ iterations.

Suppose Algorithm 1 terminates at the t-th iteration. Let $\mathcal{T}_{succ}$ and $\mathcal{T}_{fail}$ denote the sets of all the successful and unsuccessful iterations, respectively. Then $T=\left|\mathcal{T}_{succ}\right|+\left|\mathcal{T}_{fail}\right|$ and $\Delta_{T}=\Delta_{0}\gamma^{\left|\mathcal{T}_{succ}\right|-\left|\mathcal{T}_{fail}\right|}$, where $\gamma$ is a hyper-parameter of Algorithm 1. From Lemma 4, we have $\Delta_{T}\geq\kappa_{\Delta}\min\{\epsilon_{g},\epsilon_{H}\}$. Hence, $\left(\left|\mathcal{T}_{succ}\right|-\left|\mathcal{T}_{fail}\right|\right)\log\gamma\geq\log\left(\kappa_{\Delta}\cdot\min\{\epsilon_{g},\epsilon_{H}\}/\Delta_{0}\right)$, which implies $\left|\mathcal{T}_{fail}\right|\leq\log\left(\Delta_{0}/\left(\kappa_{\Delta}\cdot\min\{\epsilon_{g},\epsilon_{H}\}\right)\right)/\log\gamma+\left|\mathcal{T}_{succ}\right|$. Combine the result from Lemma 5, we have the total iteration complexity as

where $\kappa_{\Delta},\widehat{\kappa}_{\Delta},\widetilde{\kappa}_{\Delta}$ are defined in the proofs of Lemmas 4 and 5, respectively. ∎∎

As it can be seen, the worst-case total number of iterations required by Algorithm 1 before termination, matches the optimal iteration complexity obtained in . Furthermore, from (3a), it follows that upon termination of Algorithm 1 after $T$ iterations, in addition to ${\left\|\nabla F(\mathbf{x}_{T})\right\|}\leq\epsilon_{g}$, we have $\lambda_{\min}\left(\nabla^{2}F(\mathbf{x}_{T})\right)\geq-(\epsilon_{H}+\epsilon_{T})$, i.e., the obtained solution satisfies $(\epsilon_{g},\epsilon_{T}+\epsilon_{H})$-Optimality as in (1).

For Algorithm 1, the Hessian approximation tolerance $\epsilon_{t}$ is allowed to be chosen per-iteration as $\epsilon_{t}\leq\mathcal{O}\left(\max\{\epsilon_{H},\Delta_{t}\}\right)$. This way, when $\Delta_{t}$ is large (e.g., at the beginning of iterations), one can employ crude Hessian approximations. As iterations progress towards optimality, $\Delta_{t}$ can get very small, in which case Hessian accuracy is set in the order of $\epsilon_{H}$. Note that by Lemma 4, we are always guaranteed to have $\Delta_{t}\in\Omega\left(\min\left\{\epsilon_{g},\epsilon_{H}\right\}\right)$. As a result, when $\epsilon_{g}\ll\epsilon_{H}$, e.g., $\epsilon_{H}^{2}=\epsilon_{g}=\epsilon$, we can have that $\Delta_{t}\ll\epsilon_{H}$. In such cases, the choice $\epsilon_{t}\leq\mathcal{O}\left(\max\{\epsilon_{H},\Delta_{t}\}\right)$ ensures that the Hessian approximation tolerance never gets unnecessarily too small.

2 Adaptive Cubic Regularization with Inexact Hessian

Similar to Section 2.1, in this section, we present the algorithm and its corresponding convergence results for the case of adaptive cubic regularization with inexact Hessian. In particular, Algorithm 2 depicts a variant of ARC algorithm where at each iteration $t$, the inexact approximation, $\mathbf{H}_{t}$, is constructed according to Condition 1. Here, unlike Section 2.1, we were unable to provide convergence guarantees with adaptive tolerance in (3a) and as result, $\epsilon$ is set fixed a priori to a sufficiently small value, i.e., $\epsilon\in\mathcal{O}(\sqrt{\epsilon_{g}},\epsilon_{H})$ to guarantee $(\epsilon_{g},\epsilon_{H})$-optimality.

Similar to Algorithm 1, here we also require that the sub-problem (12) in Algorithm 2 is solved only approximately. Although similar inexact solutions to the sub-problem (12) by using Cauchy and Eigenpoint has been considered in several previous work, e.g., , here we provide refined conditions which prove to be instrumental in obtaining iteration complexities with the relaxed Hessian approximation (3a), as opposed to the stronger Condition (5c).

Assume that we solve the sub-problem (12) approximately to find $\mathbf{s}_{t}$ such that

Here $m_{t}(\cdot)$ is defined in (12), $\mathbf{s}_{t}^{C}$ (Cauchy point) is along negative gradient direction and $\mathbf{s}_{t}^{E}$ is along approximate negative curvature direction such that $\langle\mathbf{s}_{t}^{E},\mathbf{H}_{t}\mathbf{s}_{t}^{E}\rangle\leq\nu\lambda_{\min}(\mathbf{H}_{t})\|\mathbf{s}_{t}^{E}\|^{2}<0$ for some $\nu\in(0,1]$ (see Appendix B for a way to efficiently compute $\mathbf{s}_{t}^{E}$).

Note that Condition (13) describes the quality of the descent obtained by Cauchy and Eigen directions more accurately than is usually found in similar literature. A natural way to ensure that the approximate solution to the sub-problem (12) satisfies (13), is by replacing the unconstrained high-dimensional sub-problem (12) with the following constrained but lower-dimensional problem, in which the search space is reduced to a two-dimensional sub-space containing vectors $\mathbf{s}_{t}^{C}$, and $\mathbf{s}_{t}^{E}$, i.e.,

and set $\mathbf{s}_{t}=\mathbf{U}\mathbf{v}_{t}$. As before, any larger dimensional sub-space $\mathcal{P}$ for which we have $\text{Span}\{\mathbf{s}_{t}^{C},\mathbf{s}_{t}^{E}\}\subseteq\mathcal{P}$ would also ensure (13), and, indeed, implies a more accurate solution to our original sub-problem (12).

Lemmas 6 and 7 describe the model reduction obtained by Cauchy and eigen points as required by Condition (3).

Consider the Cauchy direction as $\mathbf{s}_{t}^{C}=-\alpha\nabla F(\mathbf{x}_{t})$ where $\alpha=\arg\min_{\widehat{\alpha}\geq 0}m_{t}(-\widehat{\alpha}\nabla F(\mathbf{x}_{t}))$. We have

Consider the function $h(x;\beta)=\sqrt{x^{2}+\beta}-x$. It is easy to verify that, for $\beta\geq 0$, $h(x)$ is decreasing function of $x$. Now since $\langle\nabla F(\mathbf{x}_{t}),\mathbf{H}_{t}\nabla F(\mathbf{x}_{t})\rangle\leq K_{H}\|\nabla F(\mathbf{x}_{t})\|^{2}$, we get

Alternatively, following the proof of [9, Lemma 2.1], for any $\alpha\geq 0$, we get

Consider the quadratic polynomial $r(\alpha)=2\alpha^{2}\sigma_{t}\|\nabla F(\mathbf{x}_{t})\|+3\alpha K_{H}-6$. We have $r(\alpha)\leq 0$ for $\alpha\in[0,\bar{\alpha}]$, where

Note that $\sqrt{9K_{H}^{2}+48\sigma_{t}\|\nabla F(\mathbf{x}_{t})\|}\leq 8\sqrt{3}\max\left\{K_{H},\sqrt{\sigma_{t}\|\nabla F(\mathbf{x}_{t})\|}\right\}$ and trivially $3K_{H}\leq 4\sqrt{3}\max\left\{K_{H},\sqrt{\sigma_{t}\|\nabla F(\mathbf{x}_{t})\|}\right\}$. Hence, defining $\alpha_{0}\triangleq 1/(\sqrt{3}\max\{K_{H},\sqrt{\sigma_{t}\|\nabla F(\mathbf{x}_{t})\|}\})$, it is easy to see that $0<\alpha_{0}\leq\bar{\alpha}$. With this $\alpha_{0}$, we get $r(\alpha_{0})\leq 2/9+3/\sqrt{3}-6\leq-3.$ Therefore

By the first-order necessary optimality condition of $\alpha$, we get $\langle\nabla F(\mathbf{x}_{t}),\mathbf{u}_{t}\rangle+\alpha\langle\mathbf{u}_{t},\mathbf{H}_{t}\mathbf{u}_{t}\rangle+\sigma_{t}\alpha^{2}\|\mathbf{u}_{t}\|^{3}=0$, which implies $\langle\nabla F(\mathbf{x}_{t}),\mathbf{s}_{t}^{E}\rangle+\langle\mathbf{s}_{t}^{E},\mathbf{H}_{t}\mathbf{s}_{t}^{E}\rangle+\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}=0$. Next, since $\alpha$ is a minimizer of $m_{t}(\widehat{\alpha}\mathbf{u}_{t})$, we have $m_{t}(\alpha\mathbf{u}_{t})\leq m_{t}(-\alpha\mathbf{u}_{t})$, which implies $\langle\nabla F(\mathbf{x}_{t}),\mathbf{s}_{t}^{E}\rangle\leq 0$. Hence, we also obtain $\langle\mathbf{s}_{t}^{E},\mathbf{H}_{t}\mathbf{s}_{t}^{E}\rangle+\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}\geq 0$. From [11, Lemma 2.1], we get $-m_{t}(s_{t}^{E})\geq\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}/6=\left(-\langle\nabla F(\mathbf{x}_{t}),\mathbf{s}_{t}^{E}\rangle-\langle\mathbf{s}_{t}^{E},\mathbf{H}_{t}\mathbf{s}_{t}^{E}\rangle\right)/6\geq\nu|\lambda_{\min}(\mathbf{H}_{t})|\|\mathbf{s}_{t}^{E}\|^{2}/6$. Now, we have

which gives $\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}\geq\nu|\lambda_{\min}(\mathbf{H}_{t})|\|\mathbf{s}_{t}^{E}\|^{2}$ and $\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}\geq\nu^{3}\sigma_{t}^{-2}|\lambda_{\min}(\mathbf{H}_{t})|^{3}$. Hence, we have $-m_{t}(\mathbf{s}_{t}^{E})\geq\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}/6\geq\nu|\lambda_{\min}(\mathbf{H}_{t})|\|\mathbf{s}_{t}^{E}\|^{2}/6$ and $-m_{t}(\mathbf{s}_{t}^{E})\geq\sigma_{t}\|\mathbf{s}_{t}^{E}\|^{3}/6\geq\nu^{3}\sigma_{t}^{-2}|\lambda_{\min}(\mathbf{H}_{t})|^{3}/6$. ∎∎

The next lemma is used to show sufficient decrease in the objective function using the approximate solution of the sub-problem (12).

Given Assumption 1 and Condition 1, we have

Apply Mean Value Theorem on $F$ at $\mathbf{x}_{t}$ gives $F(\mathbf{x}_{t}+\mathbf{s}_{t})=F(\mathbf{x}_{t})+\nabla F(\mathbf{x}_{t})^{T}\mathbf{s}_{t}+\frac{1}{2}\mathbf{s}_{t}^{T}\nabla^{2}F(\xi_{t})\mathbf{s}_{t}$, for some $\xi_{t}$ in the segment of $[\mathbf{x}_{t},\mathbf{x}_{t}+\mathbf{s}_{t}]$. Now, it follows that

Given Assumption 1, Conditions 1 and 3, suppose

First consider $\|\mathbf{s}_{t}^{C}\|$ for which we have two cases.

If $\|\mathbf{s}_{t}\|\leq\|\mathbf{s}_{t}^{C}\|$, then from assumption on $\sigma_{t}$, it immediately follows that

If $\|\mathbf{s}_{t}\|\geq\|\mathbf{s}_{t}^{C}\|$, since $L\leq\sigma_{t}/2$, then

The second last inequality follows from (14).

Similarly, for $\|\mathbf{s}_{t}^{E}\|$, we have two cases.

If $\|\mathbf{s}_{t}\|\leq\|\mathbf{s}_{t}^{E}\|$, then from assumption on $\sigma_{t}$, it immediately follows that

If $\|\mathbf{s}_{t}\|\geq\|\mathbf{s}_{t}^{E}\|$, since $L\leq\sigma_{t}/2$, then

The second last inequality follows from (15) and the last line follows from $\epsilon\leq\frac{\nu\epsilon_{H}}{6}$. ∎

Given Assumption 1, Conditions 1 and 3, suppose at the t-th iteration, $\lambda_{\min}(\mathbf{H}_{t})<-\epsilon_{H}$, $\sigma_{t}\geq 2L$, and $\epsilon\leq\min\{1/6,(1-\eta)/3\}\nu\epsilon_{H}$. Then, the t-th iteration is successful, i.e. $\sigma_{t+1}=\sigma_{t}/\gamma$.

From (13b), Lemma 8, Lemma 9, as well as assumptions on $\sigma_{t}$ and $\epsilon$, we have

Hence, $\rho_{t}\geq\eta$, and the iteration is successful. ∎∎

Given Assumption 1, Conditions 1 and 3, suppose at the t-th iteration, $\|\nabla F(\mathbf{x}_{t})\|\geq\epsilon_{g}$, $\sigma_{t}\geq 2L$, and

Then, the t-th iteration is successful, i.e. $\sigma_{t+1}=\sigma_{t}/\gamma$.

Hence, again, by (13a), Lemma 8 and 9, it follows that

Hence, $\rho_{t}\geq\eta$, and the iteration is successful. ∎∎

Now we can upper bound the cubic regularization parameter before the algorithm terminates, as in Lemma 12.

Consider Assumption 1, Conditions 1 and 3, and

where $\nu,L,K_{H}$ are, respectively, defined as in (13b), (6a), (3b), and $\eta$ is a hyper-parameter of Algorithm 2. For Algorithm 2 we have for all $t$, $\sigma_{t}\leq\max\{\sigma_{0},2\gamma L\}$.

We prove by contradiction. Assume the t-th iteration is the first unsuccessful iteration such that $\sigma_{t+1}=\gamma\sigma_{t}\geq 2\gamma L$, which implies that $\sigma_{t}\geq 2L$. However, according to Lemmas 10 and 11, respectively, if $\lambda_{\min}(H_{t})<-\epsilon_{H}$ or $\|\nabla F(\mathbf{x}_{t})\|\geq\epsilon_{g}$, then the iteration is successful and hence we must have $\sigma_{t+1}=\sigma_{t}/\gamma\leq\sigma_{t}$, which is a contradiction. ∎∎

Now, similar to [11, Lemma 2.8], we can get the following result about the estimate of the total number of successful iterations before algorithm terminates.

Given Assumption 1, Conditions 1 and 3, let $\mathcal{T}_{succ}$ denote the set of all the successful iterations before Algorithm 2 stops. The number of successful iterations is upper bounded by,

where $\kappa_{\sigma}\triangleq\min\left\{\nu^{3}/(24\gamma^{2}L^{2}),\min\left\{{1}/{K_{H}},\sqrt{{1}/{(2\gamma L)}}\right\}/(2\sqrt{3})\right\}$.

Suppose Algorithm 2 doesn’t terminate at the t-th iteration. Then either we have ${\left\|\nabla F(\mathbf{x}_{t})\right\|}\geq\epsilon_{g}$ or $\lambda_{\min}(\nabla^{2}\mathbf{H}_{t})\leq-\epsilon_{H}$. In the first case, (13a) and Lemma 12 gives

Similarly, in the case where $\lambda_{\min}(\nabla^{2}\mathbf{H}_{t})\leq-\epsilon_{H}$, from (13b) and Lemma 12, we obtain $-m_{t}(\mathbf{s}_{t})\geq\nu^{3}|\lambda_{\min}(\mathbf{H}_{t})|^{3}/(6\sigma_{t}^{2})\geq\nu^{3}\epsilon_{H}^{3}/(24\gamma^{2}L^{2})$.

Since $F(\mathbf{x}_{t})$ is monotonically decreasing, we have

Now we show the final complexity bounds of Algorithm 2 in Theorem 2.

Consider any $0<\epsilon_{g},\epsilon_{H}<1$. Suppose the inexact Hessian, $\mathbf{H}(\mathbf{x})$, satisfies Condition 1 with the approximation tolerance, $\epsilon$, in (3a) as (16). For Problem (P0), under Assumption 1 and Condition 3, Algorithm 2 terminates after at most $T\in\mathcal{O}\left(\max\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\}\right)$ iterations.

Suppose Algorithm 2 terminates at the t-th iteration. Let $\mathcal{T}_{succ}$ and $\mathcal{T}_{fail}$ denote the sets of all the successful and unsuccessful iterations, respectively. Then $T=\left|\mathcal{T}_{succ}\right|+\left|\mathcal{T}_{fail}\right|$ and $\sigma_{T}=\sigma_{0}\gamma^{\left|\mathcal{T}_{fail}\right|-\left|\mathcal{T}_{succ}\right|}$. From Lemma 12, we have $\sigma_{T}\leq 2\gamma L$. Hence, $\left|\mathcal{T}_{fail}\right|\leq\log\left({2\gamma L}/{\sigma_{0}}\right)/\log\gamma+\left|\mathcal{T}_{succ}\right|$, which, using Lemma 13 gives the total iteration complexity as

where $\kappa_{\sigma}$ is defined in Lemma 13. ∎∎

In Theorem 2 (as well as Theorem 3 below), we require $\epsilon\in\mathcal{O}(\sqrt{\epsilon_{g}},\epsilon_{H})$. This can be rather strict and computationally unattractive, unless either crude solutions are required (e.g., in most machine learning applications very rough solutions are encouraged to avoid over-fitting), or the inexact Hessian is formed from a sub-set of data that is significantly smaller than the original dataset (e.g., see Section 3 in the context of big-data regimes where $n\gg 1$ and $|\mathcal{S}|\ll n$). Nonetheless, the theoretical existence of such tolerance, though small, implies a certain level of robustness of the algorithm, i.e., the complexity of the algorithm is not adversely affected by small errors in Hessian computations.

We note that, for iterations where $\epsilon\ll\|\mathbf{s}_{t}\|$, (3a) is indeed a more stringent condition than (5c). As iterations progress towards optimality, step-size can become small, in which case (3a) might be theoretically more preferable. Nonetheless, beyond a direct theoretical comparison among various Hessian approximation bounds in terms of their tightness, the main advantage of (3a) should be regarded in light of its simplicity, which allows for direct constructions of $\mathbf{H}_{t}$ with a priori guarantees.

Condition 3 seems to be the bare minimum required to guarantee convergence to an approximate second-order criticality. Intuitively, however, if an approximate solution to the sub-problem (12) satisfies more than (13), i.e., if we solve (12) more exactly than just requiring (13), one could expect to be able to improve upon the iteration complexity of Theorem 2. Indeed, suppose we solve the reduced sub-problem on progressively embedded sub-spaces with increasingly higher dimensions, all of which including “$\text{Span}\{\mathbf{s}_{t}^{C},\mathbf{s}_{t}^{E}\}$”, and stop when the corresponding solution $\mathbf{s}_{t}$ satisfies the following conditions.

Assume that we solve the sub-problem (12) approximately to find $\mathbf{s}_{t}$ such that, in addition to (13), we have

for some prescribed $\zeta\in(0,1)$. Here, $m_{t}(\cdot)$ is defined in (12).

Conditions on the inexactness of the sub-problems were initially pioneered in . However, the main drawback for these conditions is that the inexactness tolerance is closely tied with the magnitude of the gradient. More specifically, when gradient is small, e.g., near saddle points, the sub-problems are required to be solved exceedingly more accurately. In fact, at a saddle point where $\|\nabla F(\mathbf{x}_{t})\|=0$, these conditions imply an exact solution to the sub-problem. To the best of our knowledge, Condition 4 represents a novel criterion, which offers the best of both worlds: when gradient is large, we allow for crude solutions to the sub-problem, but near saddle-points where the gradient is small, inexactness will be determined by the step length, which can be significantly larger than the gradient. Using Condition 4, we can obtain the optimal iteration complexity for Algorithm 2, as shown in Theorem 3. First, we prove the following two lemmas which will be used later for the proof of Theorem 3.

Suppose $\|\nabla F(\mathbf{x}_{t})\|\geq\epsilon_{g}$. Given Assumption 1 and Condition 3, let (3a) hold with $\epsilon_{t}=\min\{\epsilon,\zeta\|\nabla F(\mathbf{x}_{t})\|\}$ where $\epsilon$ is as in (16) and $\zeta\in(0,1/2)$. Furthermore, suppose (12) is solved such that Condition 4 eventually holds. Then, we have $\|\mathbf{s}_{t}\|\geq\kappa_{g}\sqrt{\|\nabla F(\mathbf{x}_{t+1})\|}$, where

First, suppose $\|\mathbf{s}_{t}\|^{2}\leq\theta_{t}\|\nabla F(\mathbf{x}_{t})\|$. Using Condition 4, we get $\|\nabla F(\mathbf{x}_{t+1})\|\leq\|\nabla F(\mathbf{x}_{t+1})-\nabla m_{t}(\mathbf{s}_{t})\|+\|\nabla m_{t}(\mathbf{s}_{t})\|\leq\|\nabla F(\mathbf{x}_{t+1})-\nabla m_{t}(\mathbf{s}_{t})\|+\theta_{t}\|\nabla F(\mathbf{x}_{t})\|$. Noting that $\nabla m_{t}(\mathbf{s}_{t})=\nabla F(\mathbf{x}_{t})+\mathbf{H}_{t}\mathbf{s}_{t}+\sigma_{t}\|\mathbf{s}_{t}\|\mathbf{s}_{t}$, and using Mean Value Theorem for vector-valued functions, (6a) and (3a), we get

where the last equality follows from Lemma 12. From (6b), it follows that

As such, using $\theta_{t}\leq\zeta$ from Condition 4 as well as the assumption on $\epsilon_{t}$, we get

which implies that $(1-\zeta)\|\nabla F(\mathbf{x}_{t+1})\|\leq\left(L/2+2\gamma L\right)\|\mathbf{s}_{t}\|^{2}+\left(\epsilon_{t}+\theta_{t}K\right)\|\mathbf{s}_{t}\|$. Now using Condition 4, we consider two cases:

If $\|\mathbf{s}_{t}\|\geq 1$, then we get $\left(\epsilon_{t}+\theta_{t}K\right)\|\mathbf{s}_{t}\|\leq\left(\epsilon_{t}+\theta_{t}K\right)\|\mathbf{s}_{t}\|^{2}\leq(\epsilon+\zeta K)\|\mathbf{s}_{t}\|^{2}$. Hence, it follows that $(1-\zeta)\|\nabla F(\mathbf{x}_{t+1})\|\leq\left(L/2+2\gamma L+(\epsilon+\zeta K)\right)\|\mathbf{s}_{t}\|^{2}$.

If $\|\mathbf{s}_{t}\|\leq 1$, then from assumption on $\epsilon_{t}$ and (18) , we have $\epsilon_{t}\|\mathbf{s}_{t}\|\leq\zeta\|\nabla F(\mathbf{x}_{t})\|\|\mathbf{s}_{t}\|\leq\zeta(K\|\mathbf{s}_{t}\|^{2}+\|\nabla F(\mathbf{x}_{t+1})\|\|\mathbf{s}_{t}\|)\leq\zeta(K\|\mathbf{s}_{t}\|^{2}+\|\nabla F(\mathbf{x}_{t+1})\|)$. Now by assumption on $\theta_{t}$, we get $\left(\epsilon_{t}+\theta_{t}K\right)\|\mathbf{s}_{t}\|=\epsilon_{t}\|\mathbf{s}_{t}\|+\theta_{t}K\|\mathbf{s}_{t}\|\leq 2\zeta K\|\mathbf{s}_{t}\|^{2}+\zeta\|\nabla F(\mathbf{x}_{t+1})\|$, which, in turn, implies that $(1-2\zeta)\|\nabla F(\mathbf{x}_{t+1})\|\leq\left(L/2+2\gamma L+2\zeta K\right)\|\mathbf{s}_{t}\|^{2}$.

Now suppose, $\|\mathbf{s}_{t}\|^{2}\geq\theta_{t}\|\nabla F(\mathbf{x}_{t})\|$. As above, we have $\|\nabla F(\mathbf{x}_{t+1})\|\leq\|\nabla F(\mathbf{x}_{t+1})-\nabla m_{t}(\mathbf{s}_{t})\|+\|\nabla m_{t}(\mathbf{s}_{t})\|\leq\left({L}/{2}+2\gamma L+\zeta\right)\|\mathbf{s}_{t}\|^{2}+\epsilon_{t}\|\mathbf{s}_{t}\|$. If $\|\mathbf{s}_{t}\|\geq 1$, we have $\epsilon_{t}\|\mathbf{s}_{t}\|\leq\epsilon\|\mathbf{s}_{t}\|^{2}$, which gives $\|\nabla F(\mathbf{x}_{t+1})\|\leq\left({L}/{2}+2\gamma L+\zeta+\epsilon\right)\|\mathbf{s}_{t}\|^{2}$. Otherwise, if $\|\mathbf{s}_{t}\|\leq 1$, then $\|\mathbf{s}_{t}\|^{2}\geq\theta_{t}\|\nabla F(\mathbf{x}_{t})\|$ implies that $\|\mathbf{s}_{t}\|\geq\|\nabla F(\mathbf{x}_{t})\|$. From assumption on $\epsilon_{t}$, it follows that $\epsilon_{t}\|\mathbf{s}_{t}\|\leq\zeta\|\nabla F(\mathbf{x}_{t})\|\|\mathbf{s}_{t}\|\leq\zeta\|\mathbf{s}_{t}\|^{2}$, which in turn gives $\|\nabla F(\mathbf{x}_{t+1})\|\leq\left({L}/{2}+2\gamma L+2\zeta\right)\|\mathbf{s}_{t}\|^{2}$. ∎∎

be the set of all successful iterations, before Algorithm 2 terminates. Under the conditions of Lemma 14, we must have $\left|\mathcal{T}_{\text{succ}}\right|\in\mathcal{O}(\max\{\epsilon_{H}^{-3},\epsilon_{g}^{-3/2}\})$.

From (13b) and Lemma 12, if $\lambda_{\min}(\nabla^{2}\mathbf{H}_{t})\leq-\epsilon_{H}$, it follows that $-m_{t}(\mathbf{s}_{t})\geq\nu^{3}|\lambda_{\min}(\mathbf{H}_{t})|^{3}/(6\sigma_{t}^{2})\geq\nu^{3}\epsilon_{H}^{3}/(24\gamma^{2}L^{2})$. Note that $\mathcal{T}_{\text{succ}}=\mathcal{T}^{1}_{\text{succ}}\bigcup\mathcal{T}^{2}_{\text{succ}}\bigcup\mathcal{T}^{3}_{\text{succ}},$where

We bound each of these sets individually. Since $F(\mathbf{x}_{t})$ is monotonically decreasing, from [9, Lemma 3.3], $\sigma_{t}\geq\sigma_{\min}$, and Lemmas 12 and 14, we have

Hence, $\left|\mathcal{T}^{1}_{\text{succ}}\right|\leq\kappa^{1}_{\mathcal{T}_{\text{succ}}}\max\{\epsilon_{H}^{-3},\epsilon_{g}^{-3/2}\}$, where

As for $\mathcal{T}^{2}_{\text{succ}}$, we have

Hence, $\left|\mathcal{T}^{2}_{\text{succ}}\right|\leq\kappa^{2}_{\mathcal{T}_{\text{succ}}}\epsilon_{H}^{-3}$, where $\kappa^{2}_{\mathcal{T}_{\text{succ}}}\triangleq(F(\mathbf{x}_{1})-F_{\min}){24\gamma^{2}L^{2}}/(\eta\nu^{3})$. Finally, we have $\left|\mathcal{T}^{3}_{\text{succ}}\right|=1$, because in such a case, the algorithm stops in one iteration. Putting these bounds all together, we get $\left|\mathcal{T}_{\text{succ}}\right|\leq\max\{1,\kappa^{1}_{\mathcal{T}_{\text{succ}}},\kappa^{2}_{\mathcal{T}_{\text{succ}}}\}\max\{\epsilon_{H}^{-3},\epsilon_{g}^{-3/2}\}$. ∎∎

Now we can obtain the optimal complexity bound of Algorithm 2 in Theorem 3. The proof follows similarly as that of Theorem 2, and hence is omitted here.

Consider any $0<\epsilon_{g},\epsilon_{H}<1$. Suppose the inexact Hessian, $\mathbf{H}(\mathbf{x})$, satisfies Conditions (3) with the approximation tolerance, $\epsilon$, in (3a) as $\epsilon=\min\{\epsilon_{0},\zeta\epsilon_{g}\}$ where $\epsilon_{0}$ is as in (16), and $\zeta\in(0,1/2)$. For Problem (P0) and under Assumption 1, if the approximate solution to the sub-problem (12) satisfies Conditions 3 and 4, then Algorithm 2 terminates after at most $T\in\operatorname{\mathcal{O}}\left(\max\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\}\right)$ iterations.

From (3a), upon termination of Algorithm 2, the obtained solution satisfies $(\epsilon_{g},\epsilon+\epsilon_{H})$-Optimality as in (1), i.e., ${\left\|\nabla F(\mathbf{x}_{T})\right\|}\leq\epsilon_{g}$ and $\lambda_{\min}\left(\nabla^{2}F(\mathbf{x}_{T})\right)\geq-(\epsilon_{H}+\epsilon)$.

Finite-Sum Minimization

In this section, we give concrete and practical examples to demonstrate ways to construct the approximate Hessian, which satisfies Condition 1. By considering finite-sum minimization, a ubiquitous problem arising frequently in machine learning, we showcase the practical benefits of the proposed relaxed requirement (3a) for approximating Hessian, compared to the stronger alternative (5c). In Section 3.1, we describe randomized techniques to appropriately construct the approximate Hessian, followed by the convergence analysis of Algorithms 1 and 2 with such Hessian approximations in Section 3.2.

Indeed, a major advantage of (3a) over (5c) is that there are many approximation techniques that can produce an inexact Hessian satisfying (3a). Of particular interest in our present paper is the application of randomized matrix approximation techniques, which have recently shown great success in the area of RandNLA at solving various numerical linear algebra tasks . For this, we consider the highly prevalent finite-sum minimization problem (P1) and employ random sampling as a way to construct approximations to the exact Hessian, which are, probabilistically, ensured to satisfy (3a). Many machine learning and scientific computing applications involve finite-sum optimization problems of the form (P1) where each $f_{i}$ is a loss (or misfit) function corresponding to $i^{th}$ observation (or measurement), e.g., see and references therein.

Here, we consider (P1) in large-scale regime where $n,d\gg 1$. In such settings, the mere evaluations of the Hessian and the gradient increase linearly in $n$. Indeed, for big-data problems, the operations with the Hessian, e.g., matrix-vector products involved in the (approximate) solution of the sub-problems (8) and (12), typically constitute the main bottleneck of computations, and in particular when $n\gg 1$, are computationally prohibitive. For the special case of (P1) in which each $f_{i}$ is convex, randomized sub-sampling has shown to be effective in reducing such costs, e.g., . We now show that such randomized approximation techniques can indeed be effectively employed for the non-convex settings considered in this paper.

In this light, suppose we have a probability distribution, $\mathbf{p}=\{p_{i}\}_{i=1}^{n}$, over the set $\{1,2,\ldots,n\}$, such that for each index $i=1,2\ldots,n$, we have $\Pr(i)=p_{i}>0$ and $\sum_{i=1}^{n}p_{i}=1$. Consider picking a sample of indices from $\{1,2,\ldots,n\}$, at each iteration, randomly according to the distribution $\mathbf{p}$. Let $\mathcal{S}$ and $|\mathcal{S}|$ denote the sample collection and its cardinality, respectively and define

to be the sub-sampled Hessian. In big-data regime when $n\gg 1$, if $|\mathcal{S}|\ll n$, such sub-sampling can offer significant computational savings.

In this case, we can naturally consider uniform distribution over $\{1,2,\ldots,n\}$, i.e., $p_{i}=1/n,;\forall i$. Lemma 16 gives the sample size required for the inexact Hessian, $\mathbf{H}(\mathbf{x})$, to probabilistically satisfy (3), for when the indices are picked uniformly at random with or without replacement.

Given (20a), (20b) , and $0<\epsilon,\delta<1$, let

Hence, we can apply Operator-Bernstein inequality [35, Theorem 1] to get

Now (21) ensure that $2d\exp\{-\epsilon^{2}|\mathcal{S}|/(16K_{\max}^{2})\}\leq\delta$, which gives (22). ∎∎

Indeed, if (22) holds, then (3a) follows with the same probability. In addition, if $H$ is constructed according to Lemma 16, it is easy to see that (3b) is satisfied with $K_{H}=K_{\max}$ (in fact this is a deterministic statement). These two, together, imply that $H$ satisfies Condition 1, with probability $1-\delta$.

In certain settings, one might be able to construct a more “informative” distribution, $\mathbf{p}$, over the indices in the set $\{1,2,\ldots,n\}$, as opposed to oblivious uniform sampling. In particular, it might be advantageous to bias the probability distribution towards picking indices corresponding to those $f_{i}$’s which are more relevant, in certain sense, in forming the Hessian. If this is possible, then we can only expect to require smaller sample size as compared with oblivious uniform sampling. One such setting where this is possible is the finite-sum optimization of the form (P2), which is indeed a special case of (P1) and arise often in many machine learning problems .

It is easy to see that, the Hessian of $F$ in this case can be written as $\nabla^{2}F(\mathbf{x})=\mathbf{A}^{T}\mathbf{B}\mathbf{A}=\sum_{i=1}^{n}f_{i}^{\prime\prime}(\mathbf{a}_{i}^{T}\mathbf{x})\mathbf{a}_{i}\mathbf{a}_{i}^{T}/n$, where

Note that the absolute values are needed since for non-convex $f_{i}$, we might have $f_{j}^{\prime\prime}(\mathbf{a}_{j}^{T}\mathbf{x})<0$ (for the convex case where all $f_{j}^{\prime\prime}(\mathbf{a}_{j}^{T}\mathbf{x})\geq 0$, one can obtain stronger guarantees than Lemmas 16 and 17; see ). Using non-uniform sampling distribution (23), Lemma 17 gives sampling complexity for the approximate Hessian of (P2) to, probabilistically, satisfy (3).

Given (20a), (20c) and $0<\epsilon,\delta<1$, let

Since $c=\sum_{i=1}^{n}\left|b_{i}\right|{\left\|\mathbf{a}_{i}\right\|}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left|f_{i}^{\prime\prime}\right|{\left\|\mathbf{a}_{i}\right\|}^{2}\leq\frac{1}{n}\sum_{i=1}^{n}K_{i}=\widehat{K}$ and

The bound in (24) can be improved by replacing the dimension $d$ with a smaller quantity, known as intrinsic dimension; see Appendix A. As it can be seen from (20b) and (20c), since $\widehat{K}\leq K_{\max}$, the sampling complexity given by Lemma 17 always provides a smaller sample-size compared with that prescribed by Lemma 16. Indeed, the advantage of non-uniform sampling is more pronounced in cases where the distribution of $K_{i}$’s are highly skewed, i.e., a few large ones and many small ones, in which case we can have $\widehat{K}\ll K_{\max}$; see numerical experiments in . Also, from (25), it follows that the approximate matrix $\mathbf{H}$, constructed according to Lemma 17 satisfies (3b) with $K_{H}=\widehat{K}+\epsilon$, with probability $1-\delta$, which in turn, implies that Condition 1 is ensured, with probability $1-\delta$.

As concrete examples of the problems in the form (P2) where Lemma 17 can be readily used, Table 1 gives estimates for $K_{i}$ in (20a) for robust linear regression with smooth non-convex bi-weight loss, , as well as non-convex binary-classification using logistic regression with least squares loss, .

2 Probabilistic Convergence Analysis

Now, we are in the position to give iteration complexity for Algorithms 1 and 2 where the inexact Hessian matrix $\mathbf{H}_{t}$ is constructed according to Lemmas 16 or 17. Since the approximation is a probabilistic construction, in order to guarantee success, we need to ensure that we require a small failure probability across all iterations. In particular, in order to get an overall and accumulative success probability of $1-\delta$ for the entire $T$ iterations, the per-iteration failure probability is set as $(1-\sqrt[T]{(1-\delta)})\in\mathcal{O}(\delta/T)$. This failure probability appears only in the “log factor” for sample size in all of our results, and so it is not the dominating cost. Hence, requiring that all $T$ iterations are successful for a large $T$, only necessitates a small (logarithmic) increase in the sample size. For example, for $T\in\operatorname{\mathcal{O}}(\max\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\})$, as in Theorem 2, we can set the per-iteration failure probability to $\delta\min\{\epsilon_{g}^{2},\epsilon_{H}^{3}\}$, and ensure that when Algorithm 2 terminates, all Hessian approximations have been, accumulatively, successful with probability of $1-\delta$.

Using these results, we can have the following probabilistic, but optimal, guarantee on the worst-case iteration complexity of Algorithm 1 for solving finite-sum problem (P1) (or (P2)) and in the case where the inexact Hessian is formed by sub-sampling. Their proofs follow very similar line of reasoning as that used for obtaining the results of Section 2, and hence are omitted.

Consider any $0<\epsilon_{g},\epsilon_{H},\delta<1$. Let $\epsilon$ be as in (11) and set $\delta_{0}=\delta\min\{\epsilon_{g}^{2}\epsilon_{H},\epsilon_{H}^{3}\}$. Furthermore, for such $(\epsilon,\delta_{0})$, let the sample-size $|\mathcal{S}|$ be as in (21) (or (24)) and form the sub-sampled matrix $\mathbf{H}$ as in (19). For Problem (P1) (or (P2)), under Assumption 1 and Condition 2, Algorithm 1 terminates in at most $T\in\operatorname{\mathcal{O}}(\max\{\epsilon_{g}^{-2}\epsilon_{H}^{-1},\epsilon_{H}^{-3}\})$ iterations, upon which, with probability $1-\delta$, we have that $\|\nabla F(\mathbf{x})\|\leq\epsilon_{g}$, and $\lambda_{\min}(\nabla^{2}F(\mathbf{x}))\geq-\left(\epsilon+\epsilon_{H}\right)$.

Similarly, in the setting of optimization problems (P1) and (P2), with appropriate sub-sampling of the Hessian as in Lemmas 16 and 17, we can also obtain probabilistic worst-case iteration complexities for Algorithm 2 as in the deterministic case. Again, the proofs are similar to those in Section 2, and hence are omitted.

Consider any $0<\epsilon_{g},\epsilon_{H},\delta<1$. Let $\epsilon$ be as in (16) and set $\delta_{0}=\delta\min\{\epsilon_{g}^{2},\epsilon_{H}^{3}\}$. Furthermore, for such $(\epsilon,\delta_{0})$, let the sample-size $|\mathcal{S}|$ be as in (21) (or (24)) and form the sub-sampled matrix $\mathbf{H}$ as in (19). For Problem (P1) (or (P2)), under Assumption 1 and Condition 3, Algorithm 2 terminates in at most $T\in\operatorname{\mathcal{O}}(\max\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\})$ iterations, upon which, with probability $1-\delta$, we have that $\|\nabla F(\mathbf{x})\|\leq\epsilon_{g}$, and $\lambda_{\min}(\nabla^{2}F(\mathbf{x}))\geq-\left(\epsilon+\epsilon_{H}\right)$.

Consider any $0<\epsilon_{g},\epsilon_{H},\delta<1$. Let $\epsilon$ be as in Theorem 3 and set $\delta_{0}=\delta\min\{\epsilon_{g}^{3/2},\epsilon_{H}^{3}\}$. Furthermore, for such $(\epsilon,\delta_{0})$, let the sample-size $|\mathcal{S}|$ be as in (21) (or (24)) and form the sub-sampled matrix $\mathbf{H}$ as in (19). For Problem (P1) (or (P2)), under Assumption 1, Conditions 3 and 4, Algorithm 2 terminates in at most $T\in\operatorname{\mathcal{O}}(\max\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\})$ iterations, upon which, with probability $1-\delta$, we have that $\|\nabla F(\mathbf{x})\|\leq\epsilon_{g}$, and $\lambda_{\min}(\nabla^{2}F(\mathbf{x}))\geq-\left(\epsilon+\epsilon_{H}\right)$.