Quadratic and Cubic Regularisation Methods with Inexact function and Random Derivatives for Finite-Sum Minimisation

Given $1\leq q\leq p\leq 2$, we make the following assumptions on (1).

AS.1

There exists a constant $f_{\rm low}$ such that $f(x)\geq f_{\rm low}$ for all $x\in\mathbb{R}^{n}$.

AS.2

$f\in\mathcal{C}^{p}(\mathcal{B})$ with $\mathcal{B}$ a convex neighbourhood of $\mathbb{R}^{n}$. Moreover, there exists a nonnegative constant $L_{p}$, such that, for all $x,y\in\mathcal{B}$:
$\|\nabla_{x}^{p}f(x)-\nabla_{x}^{p}f(y)\|\leq L_{p}\|x-y\|$.

To obtain complexity results for approximating first- and second-order minimisers, we use a compact formulation to characterise the optimality conditions. As in [16], given the order of accuracy $q\in\{1,2\}$, the tolerance vector $\epsilon=\epsilon_{1}$, if $q=1$, or $\epsilon=(\epsilon_{1},\epsilon_{2})$, if $q=2$, we say that $x\in\mathbb{R}^{n}$ is a $q$-th order $\epsilon$-approximate minimiser for (1) if

where

The optimaliy measure $\phi_{f,j}$ is a nonnegative (continuous) function that can be used as a measure of closeness to $q$-th order stationary points [16]. For $\epsilon_{1}=\epsilon_{2}=0$, it reduces to the known first- and second-order optimality conditions, respectively. Indeed, assuming that $q=1$ we have from (3)–(4) with $j=1$ that $\phi_{f,1}(x)=\|\nabla_{x}f(x)\|=0$. If $q=2$, (3)–(4) further imply that $\phi_{f,2}(x)=\max_{\|d\|\leq 1}\left(-\frac{1}{2}d^{\top}\nabla_{x}^{2}f(x)d\right)=0$, which is the same as requiring the semi-positive definiteness of $\nabla_{x}^{2}f(x)$.

We now define our Inexact Adaptive Regularisation ($IAR_{qp}$) scheme, whose purpose is to find a $q$-th $\epsilon$-approximate minimiser of (1) (see (3)) using a $p+1$ order model, for $1\leq q\leq p\leq 2$.

The scheme, sketched in Algorithm 1 on page 1, is defined in analogy with the Adaptive Regularization Approach (see, [17]), but now uses the inexact values $\overline{f}(x_{k})$, $\overline{f}(x_{k}+s)$, $\overline{\nabla_{x}^{j}f}(x_{k})$ instead of $f(x_{k})$, $f(x_{k}+s)$, $\nabla_{x}^{j}f(x_{k})$, $j=1,...,q$, respectively. At iteration $k$, the model

is built at Step 1 and approximately minimised, at Step 2, finding a step $s_{k}$ such that

and

for $j=1,...,q$, and some $\theta\in(0,\frac{1}{2})$. Taking into account that,

we have

for $p=1,2$ and

for $q=p=2$.

The existence of such a step can be proved as in Lemma $4.4$ of [16] (see also [3]). We note that the model definition in (7) does not depend on the approximate value of $f$ at $x_{k}$. The ratio $\rho_{k}$, depending on inexact function values and model, is then computed at Step $4$ and affects the acceptance of the trial point $x_{k}+s_{k}$. Its magnitude also influences the regularisation parameter $\sigma_{k}$ update at Step $5$. While the gradient $\overline{\nabla_{x}f}(x_{k})$ and the Hessian $\overline{\nabla_{x}^{2}f}(x_{k})$ approximations can be seen as random estimates, the values $\overline{f}(x_{k})$, $\overline{f}(x_{k}+s_{k})$ are required to be deterministically computed to satisfy

in which $\omega_{k}$ is iteratively defined at Step $6$. As for the implementation of the algorithm, we note that $\overline{\phi}_{m_{k},1}(s_{k})=\|\nabla_{s}m_{k}(s_{k})\|$, while $\overline{\phi}_{m_{k},2}(s_{k})$ can be computed via a standard trust-region method at a cost which is comparable to that of computing the Hessian left-most eigenvalue. The approximate minimisation of the cubic model (7) at each iteration can be seen, for $p=2$, as an issue in the AR framework. However, an approximate minimiser can be computed via matrix-free approaches accessing the Hessian only through matrix-vector products. A number of procedures have been proposed, ranging from Lanczos-type iterations where the minimisation is done via nested, lower dimensional, Krylov subspaces [21], up to minimisation via gradient descent (see, e.g., [1, 14, 15]) or the Barzilai-Borwein gradient method [10]. Hessian-vector products can be approximated by the finite difference approximation, with only two gradient evaluations [10, 15]. All these matrix-free implementations remain relevant if $\nabla_{x}^{2}f(x_{k})$ is defined via subsampling, proceeding as in Section IV. Interestingly, back-propagation-like methods in deep learning also allow computations of Hessian-vector products at a similar cost [33, 38].

In what follows, all random quantities are denoted by capital letters, while the use of small letters denotes their realisations. We refer to the random model $M_{k}$ at iteration $k$, while $m_{k}$ = $M_{k}(\zeta_{k})$ is its realisation, with $\zeta_{k}$ being a random sample taken from a context-dependent probability space. As a consequence, the iterates $X_{k}$, as well as the regularisers $\Sigma_{k}$, the steps $S_{k}$ and $\Omega_{k}$, are the random variables such that $x_{k}=X_{k}(\zeta_{k})$, $\sigma_{k}=\Sigma_{k}(\zeta_{k})$, $s_{k}=S_{k}(\zeta_{k})$ and $\omega_{k}=\Omega_{k}(\zeta_{k})$. For the sake of brevity, we will omit $\zeta_{k}$ in what follows. Due to the randomness of the model construction at Step $1$, the $IAR_{qp}$ algorithm induces a random process formalised by $\{X_{k},S_{k},M_{k},\Sigma_{k},\Omega_{k}\}$ ($x_{0}$ and $\sigma_{0}$ are deterministic quantities). For $k\geq 0$, we formalise the conditioning on the past by using $\mathcal{A}_{k-1}^{M}$, the $\hat{\sigma}$-algebra induced by the random variables $M_{0}$, $M_{1}$,…, $M_{k-1}$, with $\mathcal{A}_{-1}^{M}=\hat{\sigma}(x_{0})$. We also denote by $d_{k,j}$ and $\overline{d}_{k,j}$ the arguments in the maximum in the definitions of $\phi_{m_{k},j}(s_{k})$ and $\overline{\phi}_{m_{k},j}(s_{k})$, respectively. We say that iteration $k$ is accurate if

where

We emphasize that the above accuracy requirements are adaptive. At variance with the trust-region methods of [12, 18, 24], the above conditions do not require the model to be fully linear or quadratic in a ball centered at $x_{k}$ of radius at least $\|s_{k}\|$. As standard in related papers [12, 18, 24, 32], we assume a lower bound on the probability of the model to be accurate at the $k$-th iteration.

AS.3

For all $k\geq 0$, conditioned to $\mathcal{A}_{k-1}^{M}$, we assume that (14) is satisfied with probability at least $p_{*}\in(\frac{1}{2},1]$ independent of $k$.

We stress that the inequalities in (14) can be satisfied via uniform subsampling with probability of success at least $(1-t_{\ell})$, using the operator Bernstein inequality (see, Subsection IV-A and Section $7.2$ in [6]), and this provides AS.3 with $p_{*}=\prod_{\ell=1}^{p}(1-t_{\ell})$.

Before starting with the worst-case evaluation complexity of the $IAR_{qp}$ algorithm, the following clarifications are important. For each $k\geq 1$ we assume that the computation of $s_{k-1}$ and thus of the trial point $x_{k-1}+s_{k-1}$ ($k\geq 1$) are deterministic, once the inexact model $m_{k-1}(s)$ is known, and that (12)-(13) at Step $3$ of the algorithm are enforced deterministically; therefore, $\rho_{k-1}$ and the fact that iteration $k-1$ is successful are deterministic outcomes of the realisation of the (random) inexact model. Hence,

can be seen as a family of hitting times depending of $\epsilon$ and corresponding to the number of iterations required until (3) is met for the first time. The assumptions made on top of this subsection imply that the variables $X_{k-1}+S_{k-1}$ and the event $\{X_{k}=X_{k-1}+S_{k-1}\}$, occurring when iteration $k-1$ is successful, are measurable with respect to $\mathcal{A}_{k-1}^{M}$. Our aim is then to derive an upper bound on the expected number of steps $\mathbb{E}[N_{\epsilon}]$ needed by the $IAR_{qp}$ algorithm, in the worst-case, to reach an $\epsilon$-approximate $q$-th-order-necessary minimiser, as in (3).

A crucial property for our complexity analysis is to show that, when the model (7) is accurate, iteration $k$ is successful but $\|\nabla_{x}f(x_{k+1})\|>\epsilon_{1}$, then $\|s_{k}\|^{p+1}$ is bounded below by $\psi(\sigma_{k})\epsilon_{1}^{\frac{p+1}{p-q+1}}$ in which $\psi(\sigma)$ is a decreasing function of $\sigma$. This is central in [18] for proving the worst-case complexity bound for first-order optimality. By virtue of the compact formulation (4) of the optimality measure, for $p=2$, the lower bound on $\|s_{k}\|^{3}$ also holds for second-order optimality and a suitable power of $\epsilon_{2}$.

Theorem 4 shows that for $q=p=1$ the $IAR_{1}$ algorithm reaches the hitting time in at most $\mathcal{O}(\epsilon_{1}^{-2})$ iterations. Thus, it matches the complexity bound of gradient-type methods for nonconvex problems employing exact gradients. Moreover, $IAR_{qp}$ needs $\mathcal{O}(\epsilon_{1}^{-3/2})$ iterations when $q=1$ and $p=2$ are considered, while to approximate a second-order optimality point (case $p=q=2$) $\mathcal{O}\bigl{(}\max[\epsilon_{1}^{-{3/2}},\epsilon_{2}^{-3}]\bigr{)}$ are needed in the worst-case. Therefore, for $1\leq q\leq p\leq 2$, Theorem 4 generalises the complexity bounds stated in [16, Theorem 5.5] to the case where evaluations of $f$ and its derivatives are inexact, under probabilistic assumptions on the accuracies of the latter. Since in [16, Theorems 6.1 and 6.4] it is shown that the complexity bounds are sharp in order of the tolerances $\epsilon_{1}$ and $\epsilon_{2}$, for exact evaluations and Lipschitz continuous derivatives of f, this is also the case for the bounds of Theorem 4.

We stress that (12)-(14) ensure that $\nu_{0,k}=\omega_{k}\overline{\Delta T}_{f,p}(x_{k},s_{k})$, $\nu_{\ell,k}=\Omega_{k}\frac{\overline{\Delta T}_{k,\min}}{6\tau_{k}^{\ell}},\textrm{ for }{{\ell\in\{1,...,p\}}}.$ The computation of both $\overline{\Delta T}_{f,p}(x_{k},s_{k})$ and $\overline{\Delta T}_{k,\min}$ requires the knowledge of the step $s_{k}$, that is available at Step 3 when function estimators are built, but it is not yet available at Step 1 when approximations of the gradient (and the Hessian) have to be computed. Thus, practical implementation of Step 1 calls for an inner loop where the accuracy requirements is progressively reduced in order to meet (14). This process terminates as for sufficiently small accuracy requirements the right-hand side in (20)-(21) reaches the full sample $N$. To make clear this point a detailed description of Step 1 and Step 2 of the $IAR_{qp}$ algorithm with $q=p=1$, renamed for brevity as $IAR_{1}$, is given below.

Regarding algorithm $IAR_{qp}$ with $q=1$ and $p=2$, hereafter called $IAR_{2}$, a similar loop as in Step 1 of Algorithm2 is needed that involves the step computation. The approximated gradient and Hessian are computed via (20) and (21), using predicted accuracy requirements $\varepsilon_{1,0}$ and $\varepsilon_{2,0}$, respectively. Then, the step $s_{k}$ is computed and if (14) is not satisfied then the accuracy requirements are progressively decreased and the step is recomputed until (14) is satisfied for the first time. The values $\kappa_{\epsilon}=\gamma_{\epsilon}=0.5$ are considered to initialize and decrease the the accuracy requirements in both $IAR_{1}$ and $IAR_{2}$ Algorithms.

As one may easily see, the step $s_{k}$ taken in $IAR_{1}$ Algorithm (see the form of the model (7) and Step 2 in Algorithm 2) is simply the global minimiser the model. Therefore, for any successful iteration $k$:

that can be viewed as the general iteration of a gradient descent method under inexact gradient evaluations with an adaptive choice of the learning rate (it corresponds to the opposite of the coefficient of $\overline{\nabla f}(x_{k})$ in (28)), taken as the inverse of the quadratic regulariser $\sigma_{k}>0$.

On the other hand, in the $IAR_{2}$ procedure, an approximated minimizer $s_{k}$ of the third model $m_{k}$ (7) has to be computed. Such approximate minimizer has to satisfy $m_{k}(s_{k})\leq m_{k}(0)$ and (9), i.e. $\|\nabla_{s}m_{k}(s_{k})\|\leq\theta{\epsilon_{1}}$. In our implementation, we used $\theta=0.5$, $\epsilon_{1}=10^{-3}$ and the computation of $s_{k}$ has been carried out using the Barzilai-Borwein method [34] with a non-monotone line-search following [25, Algorithm 1] and using parameter values as suggested therein. The Hessian-vector product required at each iteration of the Barzilai-Borwein procedure is approximated via finite-differences in the model gradient.

To assess the computational cost of the $IAR_{1}$ and $IAR_{2}$ Algorithms, we follow the approach in [35] for what concerns the cost measure definition. Specifically, at the generic iteration $k$, we count the $N$ forward propagations needed to evaluate the objective in (18) at $x_{k}$ has a unit Cost Measure (CM), while the evaluation of the approximated gradient at the same point requires an extra cost of $\left(2|\mathcal{G}_{k}\smallsetminus(\mathcal{G}_{k}\cap\mathcal{D}_{k,1})|+|\mathcal{G}_{k}\cap\mathcal{D}_{k,1}|\right)/N$ CM. Moreover, at each iteration of the Barzilai-Borwein procedure of the $IAR_{2}$ Algorithm, the approximation of the Hessian-vector product by finite differences calls for the computation of $\overline{\nabla f}(x_{k}+hv)$, ($h\in\mathbb{R}^{+}$), leading to additional $|\mathcal{H}_{k}|$ forward and backward propagations, at the price of the weighted cost $2|\mathcal{H}_{k}|/N$ CM, and a potential extra-cost of $|\mathcal{H}_{k}\smallsetminus(\mathcal{G}_{k}\cap\mathcal{H}_{k})|/N$ CM to compute $\nabla f_{i}(x_{k})$, for $i\in\mathcal{H}_{k}\smallsetminus(\mathcal{G}_{k}\cap\mathcal{H}_{k})$. A budget of $80$ and $100$ CM have been considered as the stopping rule for our runs for the MNIST-B and the GISETTE dataset, respectively.
We now solve the binary classification task described at the beginning of Section IV on the MNIST-B and GISETTE datasets, described in Table I. We test the $IAR_{1}$ algorithm on both datasets without net, while different ANN structures are considered for the MNIST-B via $IAR_{1}$. A comparison between $IAR_{2}$ and $IAR_{1}$ is performed on the GISETTE dataset. The accuracy in terms of the binary classification rate on the testing set, for each of the performed methods and ANN structures, is reported in Table II, averaging over 20 runs.