Stochastic Trust Region Methods with Trust Region Radius Depending on Probabilistic Models

TL;DR

This paper introduces a stochastic trust-region method called STRME, where the trust-region radius depends on probabilistic models and the latest gradient, with analyzed complexity and demonstrated benefits through numerical experiments.

Contribution

The paper proposes a novel stochastic trust-region algorithm with a radius linearly dependent on the latest gradient, and provides complexity analysis and empirical validation.

Findings

STRME matches existing algorithms in complexity analysis.

Numerical experiments show STRME outperforms existing stochastic trust-region methods.

The trust-region radius dependence improves convergence behavior.

Abstract

We present a stochastic trust-region model-based framework in which its radius is related to the probabilistic models. Especially, we propose a specific algorithm, termed STRME, in which the trust-region radius depends linearly on the latest model gradient. The complexity of STRME method in non-convex, convex and strongly convex settings has all been analyzed, which matches the existing algorithms based on probabilistic properties. In addition, several numerical experiments are carried out to reveal the benefits of the proposed methods compared to the existing stochastic trust-region methods and other relevant stochastic gradient methods.

Tables

Table 1. Table 1: Datasets for regularized logisitic regression

dataset	$n$	$d$	$\lambda$
a9a	32561	123	${10}^{-4}$
ijcnn1	49990	22	${10}^{-4}$

Equations

$$\underset{x\in {\mathbb{R}}^d}{\min }f(x),$$

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

$$\begin{array}{rl}\Vert \nabla f(y)-\nabla m_k(y)\Vert & \le {\kappa }_{eg}{\delta }_k,\text{and}\\ \left|f(y)-m_k(y)\right|& \le {\kappa }_{ef}{\delta }_k^2.\end{array}$$

$$I_k=\mathbb{I}\left\{M_k\text{ is a }\kappa \text{-fully linear model of}f\text{on}B(x_k,\delta ({\mu }_k))\right\}$$

$$P(I_k=1|{\mathcal{F}}_{k-1}^{M\cdot F})\ge \alpha .$$

$$\begin{array}{rl}\left|f_k^0-f(x_k)\right|& \le {ϵ}_F{\delta }_k^2,\text{and}\\ \left|f_k^d-f(x_k+d_k)\right|& \le {ϵ}_F{\delta }_k^2.\end{array}$$

$$J_k=\mathbb{I}\left\{F_k^0\text{and}{F}_k^d\text{are}{ϵ}_F\text{-accurate estimates of}f(x_k)\text{and}f(x_k+d_k),\text{respectively}\right\}$$

$$P(J_k=1|{\mathcal{F}}_{k-1/2}^{M\cdot F})\ge \beta ,$$

$$\begin{array}{r}\mathbb{E}[\left|F_k^0-f(X_k)\right||{\mathcal{F}}_{k-1/2}^{M,F}]\le {\kappa }_f{\Delta }_k^2,\text{and}\\ \mathbb{E}[\left|F_k^d-f(X_k+D_k)\right||{\mathcal{F}}_{k-1/2}^{M,F}]\le {\kappa }_f{\Delta }_k^2,\end{array}$$

$$P(W_k=1|{\mathcal{F}}_{k-1}^{M,F})=p,P(W_k=-1|{\mathcal{F}}_{k-1}^{M,F})=1-p,$$

$${\Lambda }_{k+1}\ge \min ({\Lambda }_k{e}^{{\lambda }_1{W}_k},\hat{\Lambda }),$$

$$\mathbb{E}[{\Phi }_{k+1}|{\mathcal{F}}_{k-1}^{M,F}]\le {\Phi }_k-Ch({\Lambda }_k).$$

$$\mathbb{E}[T_{ϵ}]\le \frac{p}{2p-1}(\frac{{\Phi }_0}{h(\hat{\Lambda })}+\frac{h({\Lambda }_0)}{h(\hat{\Lambda })}+1).$$

$$m_k(x_k+d)=f_k+g_k^{T}d+\frac{1}{2}{d}^T{B}_{k}d,$$

$$m_k(x_k)-m_k(x_k+d_k)\ge {\kappa }_{fcd}\Vert g_k\Vert \min \left\{\frac{\Vert g_k\Vert }{\Vert B_k\Vert },{\delta }_k\right\}.$$

$$\Vert B_k\Vert \le {\kappa }_{bhm}.$$

$${\mu }_k\le \min \left\{\frac{1}{{\kappa }_{bhm}},\frac{{\kappa }_{fcd}}{8{\kappa }_{ef}}\right\},$$

$$f(x_k)-f(x_k+d_k)\ge \frac{{\kappa }_{fcd}}{4}{\Vert g_k\Vert }^2{\mu }_k.$$

$$f(x_k)-f(x_k+d_k)\ge \frac{{\eta }_1{\kappa }_{fcd}}{4}{\Vert g_k\Vert }^2{\mu }_k.$$

$${\mu }_k\le \min \left\{\frac{1}{{\kappa }_{bhm}},\frac{{\kappa }_{fcd}(1-{\eta }_1)}{8{\kappa }_{ef}}\right\},$$

$${\Phi }_k=\nu (f(X_k)-f^{\ast })+(1-\nu )\frac{1}{L^2}{\Lambda }_k{\Vert \nabla f(X_k)\Vert }^2,$$

$$\hat{\Lambda }=\zeta ,\text{where ζ is a constant such that}\zeta \le \min \left\{{\mu }_{\max },\frac{{\kappa }_{fcd}(1-{\eta }_1)}{8{\kappa }_{ef}}\right\}.$$

$$\frac{1-\nu }{\nu }\le \min \left\{\frac{{\kappa }_{fcd}}{16\gamma {\mu }_{\max }^2},\frac{{\kappa }_{fcd}{L}^2}{32\gamma (1+{\kappa }_{eg}{\mu }_{\max }{)}^2},\frac{{\eta }_1{\kappa }_{fcd}}{8\gamma {\mu }_{\max }^2}\right\},$$

$$\alpha \beta \ge \frac{4{\gamma }^2}{4{\gamma }^2+(\gamma -1)},$$

$$\mathbb{E}[{\Phi }_{k+1}-{\Phi }_k|{\mathcal{F}}_{k-1}^{M,F}]\le -\frac{1}{2}\alpha \beta (1-\nu )(1-\frac{1}{\gamma })\frac{1}{L^2}{\Lambda }_k{\Vert \nabla f(X_k)\Vert }^2.$$

$$T_{ϵ}=\inf \left\{k\ge 0:\Vert \nabla f(X_k)\Vert \le ϵ\right\}.$$

$${\Phi }_k=\nu (f(X_k)-f^{\ast })+(1-\nu )\frac{1}{L^2}{\Lambda }_k{\Vert \nabla f(X_k)\Vert }^2.$$

$$h({\Lambda }_k)=C{\Lambda }_k{ϵ}^2$$

$$\mathbb{E}[T_{ϵ}]\le \frac{\alpha \beta }{(2\alpha \beta -1)}(\frac{M}{{ϵ}^2}+\mathcal{O}(1)),$$

$$\Vert x-x^{\ast }\Vert \le D,\forall x\in \mathcal{L}.$$

Full text

Mathematics Subject Classification: 65K05, 65K10, 90C60

1 Introduction

In this paper, we are concerned with the following unconstrained optimization problem

[TABLE]

where the objective function $f$ is assumed to be smooth and bounded from below. But we only have access to the value of $f$ and its derivative information with some noise. In recent years, the expected risk minimization (ERM) problem, which is fundamental in the field of machine learning and statistic, has become the focus of many researchers. The ERM problems can be formulated as follows:

[TABLE]

where $\mathbb{E}[\cdot]$ denotes the expectation taken with respect to the random variable $\xi\in\mathbb{R}^{d}$. However, because the probability distribution of $\xi$ is unknown in advance, solving (1.2) is intractable directly. Usually only noisy information about the gradient of $f$ is available. The empirical risk problem with a fixed amount of data (possibly very large) or the on-line setting problem where the data is flowing in sequentially, which involves an estimate of problem (1.2), is more often considered in practice. Through the whole paper, we mainly consider stochastic optimization methods to solve such kind of problems.

The classic stochastic optimization method is stochastic gradient descent (SGD) method, which dates back to the work by Robbins and Monro SGD-1951 . The method is prominent and adorable in large-scale machine learning due to simpleness and low-cost computing. However, because of the variance introduced by random sampling, the sequence of learning rate (step-size) progressively diminish both in theoretical analysis and practical implementation, which leads to slow convergence. Thus finding an appropriate learning rate is critical for the performance of SGD method, but it is not easy in practice. To deal with aforementioned issues, various adaptive gradient algorithms have emerged, for instance AdaGradAdaGrad , RMSPropRMSProp , AdamAdam , which are very popular in deep learning. Besides, variance reduction (VR) methods, to improve the performance of SGD method, are proposed, such as SVRGSVRG , SAGASAGA and SARAH SARAH . Especially, they have achieved linear convergence rate when solving strongly convex problems, which is a stronger result than that of SGD method. Furthermore, these methods have also been extended to solve non-convex problems such as the deep neural networks, and achieve good performance SAGA-nonconvex ; SVRG-nonconvex ; SARAH_nonconvex . The VR technique is applicable to the problem with a large but fixed sample set, for which the full gradient has to be calculated as a compromise to achieve the significant variance reduction. Hence they are not easy to fit on the on-line setting like SGD method and adaptive gradient methods.

Besides, many second order methods are proposed, which are known to perform better than the first-order methods on various highly nonlinear and ill-conditioned problemsoLBFGS ; SQN ; SdLBFGS ; NIM ; IQN . More recently, cubic regularization methods as a class of Newton-type variants have attracted a lot of interestXu_Peng ; CR_Kohler ; CR_Ghadimi ; CR_Jordan ; CR_VR_Lan ; CR_VR_Zhou . Especially, the VR technique is applied to improve the performance of cubic regularization methods CR_VR_Lan ; CR_VR_Zhou .

Sample averaging is a natural and well-known technique to reduce the variance of gradient or the noiseSample_average . Generally speaking, to guarantee the accuracy of the function and gradient estimators, the sample size has to increase when the algorithm goes to optimality. And the training sample is not only regarded as a fixed and finite set, that is to say that the sample averaging technique can be employed to the on-line setting.

Recently, with the success of deep neural networks, the development and analysis of methods for non-convex problems have attracted tremendous attention. As we know, traditional trust-region methods is a class of well-established and effective methods in nonlinear optimizationPowell_TR_a ; Yuan_TR_review . For details interested readers can refer to the review by YuanYuan_TR_review . With such a framework, we can utilize second order information when building the trust-region subproblem. Besides, due to the boundedness of the trust-region, the Hessian approximation matrix is not required to be positive definite. An advantage of trust-region methods is that they can be applied to non-convex and ill-conditioned problems. Although there are various effective methods as mentioned before, trust-region algorithms deserve more attention in stochastic optimization.

Actually, the traditional trust-region framework has already been considered to solve machine learning problems TR_logistic ; TR_Newton_linear_classification ; SFN_saddle ; Two_stage_TR ; Xu_Peng . Dauphin et al. SFN_saddle proposed a saddle free Newton (SFN) which exploits the exact Hessian information to escape saddle points. However its computation is high cost for large-scale and high-dimension problems. A two-stage subspace trust-region approach Two_stage_TR was proposed to train deep neural networks, in which the local second-order model is conducted based on the partial information computed from a subset of the data. But the approach lacks theoretical guarantees. TR_Newton_linear_classification ; TR_logistic are designed to solve a specific class of machine learning problems, which need to utilize the accurate derivative information to construct good models, and the accurate function values to obtain good estimators, of which the computational costs is too expensive to afford for general large-scale machine learning problems. Xu_Peng incorporates inexact Hessian information into the trust-region framework but the exact gradient and function values are required to be computed per iteration.

In Bandeira_TR_2014 ; Gratton_TR , the authors construct the inexact models to satisfy some first-order accurate conditions with sufficiently high probability when building the trust-region subproblem. Cartis and Scheinberg Cartis_Linesearch analyzes the complexity of line search and cubic regularization algorithms, which is based on the random models with certain probability but their function estimators are accurate, in non-convex, convex and strongly convex settings. These mainly focus on derivative free optimization (DFO) problems.

A stochastic trust-region algorithm named STORM for stochastic optimization setting has been introduced in STORM . Not only the model is conducted to satisfy some first-order accuracy requirements, but also the function values both at current iterate and next potential iterate are estimated with some probability, instead of the exact function values. And the liminf-type and lim-type first-order convergence has been analyzed. In addition, Blanchet et al. STORM_nonconvex has bounded the expected convergence complexity of STORM for non-convex problems. More recently, Paquette and Scheinberg S_line_search analyzes the complexity of a stochastic line search algorithm, of which the gradient and function estimators are both randomly sampled with some probability.

Previous works have established the convergence and complexity properties for such a trust-region framework. In this paper, we are particularly interested in introducing a new trust-region radius formula which depends on the latest probabilistic model, to improve the practical performance of such trust-region methods. Besides, we present an algorithm termed STRME in which the trust-region radius depends linearly on the model gradient just updated, following a piece of work initially proposed for deterministic optimizationFan_Yuan . The idea is meaningful and attractive. Because the trust region can be tailored according to the newly generated model, not just the success of the trial steps as that of STORM. Note that the trust-region radius of STORM at iteration $k$ is completely determined by the past iterations [math] to $k-1$. However, the framework we proposed is no longer that case. That is to say the trust-region radius in STRME, which related to the current model, is not measurable with respect to all the information generated from the past iterations [math] to $k-1$.

It will bring new challenges both in theoretical analysis and numerical experiments. The trust-region radius in STRME depends on the newly updated model makes it more complicated to analyze the complexity of STRME, compared to that of STORM. In S_line_search , Paquette and Scheinberg tackles this issue, of which the quantity $\Delta_{k}$ is also related to the currently updated gradient. Note that the quantity $\Delta_{k}$ in S_line_search is not really the trust-region radius of trust-region algorithms. To obtain the complexity results of STRME, we analyze the parameter $\Lambda_{k}$ (called relative trust-region radius), which depends on the past iterative information, instead of the trust-region radius itself. Our approach avoids the difficulty of analyzing the trust-region radius directly. The convergence analysis for such trust-region algorithms relies on the requirements that these quantities such as models and function estimators are sufficiently accurate with sufficiently high probability. And the accuracy of these quantities are controlled by the trust-region radius. However the trust-region radius is unknown before the model is updated. We will elaborate this at the beginning of Section 4 where we present numerical experiments.

These changes indeed bring forth some advantages. The choice of trust-region radius makes STRME algorithm scale invariant on problems. Besides, the trust-region radius which depends on the model gradient can capture more new information. We have to say that it is advisable to adjust the trust-region radius according to the latest probabilistic model. Our numerical experiments illustrate this viewpoint. We have tested on regularized logistic regression and a simple deep neural network on real datasets. The numerical experiments show that the trust-region radius of the proposed STRME method reduces asymptotically. And the oscillation in the trust-region radius is less severe in contrast to that of STORM method. In addition, we observe that the proposed algorithm can get more successful iterates after a long time training.

The major contributions of this paper are summarized as follows:

• (1)

The trust-region radius can be defined using the probabilistic model, in particular its gradient.

• (2)

We propose a specific algorithm termed STRME, in which the trust-region radius $\delta_{k}$ is linearly dependent on the norm of the model gradient. The complexity of STRME in non-convex, convex and strongly convex cases are analyzed, respectively. The expected number of iterations of STRME algorithm for non-convex problem is $O(\epsilon^{-2})$ by reaching $\left\|\nabla f(x)\right\|\leq\epsilon$, which is similar to the result in STORM_nonconvex . In addition, the expected convergence rates for general convex and strongly convex problems, which are $O(1/\epsilon)$ and $O(\log(1/\epsilon))$ for reaching $f(x)-f^{\ast}\leq\epsilon$, respectively.

• (3)

In numerical experiments, sample averaging technique is utilized to construct probabilistic models. Besides, we adopt the dogleg method to solve the trust-region subproblem for the regularized logistic regression problem. In the same way, the limited memory symmetric rank one (L-SR1) is employed to approximate the Hessian matrix, and then incorporate them into STRME algorithm to train a deep neural network problem. The results indicate that the proposed algorithm compares favorably to other stochastic optimization algorithms.

The outline of this paper is as follows. In Section 2 we give some definitions about the probabilistic models and estimates, and present a generic analysis framework based on the random models and estimates; In Section 3 we propose a specific algorithm named STRME and analyze the complexity in non-convex, convex and strongly convex cases; In Section 4 we report some numerical results on regularized logistic regression problem and a simple deep neural network to show the efficiency of STRME in different settings; In the end, we draw some conclusions in Section 5.

Notations.

Throughout this paper, we use $x^{\ast}$ to denote the global minimizer, $f^{\ast}=f(x^{\ast})$. Let $\left\|\cdot\right\|$ denote the Euclidean norm, i.e. $\left\|\cdot\right\|_{2}$, unless otherwise specified. Let $B(x,\Delta)$ denote the ball of the radius $\Delta$ around $x$. Let $\mathbb{I}\left\{A\right\}$ denote the indicator function of the event $A$, that is: if $A$ occurs, $\mathbb{I}\left\{A\right\}=1$; else, $\mathbb{I}\left\{A\right\}=0.$ A function $f\in\mathcal{C}^{1}(\mathbb{R}^{d})$, if the first derivation of $f$ exists and continuous. A function $f$ is $L$-smooth, if there is a constant $L>0$ such that $\left\|\nabla f(x)-\nabla f(y)\right\|\leq L\left\|x-y\right\|,\quad\forall\,x,y\in\mathbb{R}^{d}.$

2 A generic analysis framework based on random models and estimates

Let us first introduce a generic stochastic trust-region framework. The analysis for the framework can particularize to the specific algorithm for example the algorithm STRME proposed in Section 3, and the objective function, whether it is non-convex or convex, provided that the assumptions are satisfied.

The proposed Algorithm 2.1 covers the framework of STORM STORM . We can see that if $\delta_{k}=\mu_{k}$, the above algorithm will reduce to STORM algorithm. The main difference lies in the trust-region radius in Algorithm 2.1 which depends on the current model $m_{k}$. We directly update the parameter $\mu_{k}$ (called relative trust-region radius), not the trust-region radius $\delta_{k}$ itself, which avoids the troubles by the randomness rise to the currently random model. Of course, except the current random model, there may be other factors, for instance the iterative models of previous steps, i.e. $\left\{m_{t}\right\}_{t\leq k}$, and the difference of previous iterates $\left\{x_{t}-x_{t-1}\right\}_{t\leq k}$ and so on. The framework we present here does not involve the specific forms of the random models and the sufficient reduction condition. We will discuss them in the next part.

Note that Algorithm 2.1 generates a random process. Obviously, the randomness of the algorithm comes from the randomness of the models and estimates we have constructed per iteration. At a deep level, it is determined by the inexact or random information obtained from the problem we are trying to solve. To formalize the random process, we introduce some notations to describe the quantities of them and their realizations. Let $M_{k}$ denote the random model in $k$-th iteration, while $m_{k}=M_{k}(\omega)$ for its realization, where $\omega$ is a random variable. We know that the randomness of the models gives rise to the randomness of the iterates, relative trust-region radius and trial step produced by algorithm 2.1. These random variables are denoted by $X_{k}$, $\Lambda_{k}$ and $D_{k}$, respectively, while let $x_{k}=X_{k}(\omega)$, $\mu_{k}=\Lambda_{k}(\omega)$, and $d_{k}=D_{k}(\omega)$ to denote their realizations. Similarly, we use $\left\{F_{k}^{0},F_{k}^{d}\right\}$ to denote the random estimates of $f(X_{k})$ and $f(X_{k}+D_{k})$, while their realizations are denoted by $f_{k}^{0}=F_{k}^{0}(\omega)$ and $f_{k}^{d}=F_{k}^{d}(\omega)$. We will utilize those notations to analyze the random process later in this section, which is under some assumptions that model $M_{k}$ and estimates $F_{k}^{0},F_{k}^{d}$ are sufficiently accurate with some probability conditioned on the past. In order to formalize all the randomized information before $k$-th iteration, let $\mathcal{F}_{k-1}^{M\cdot F}$ denote the $\sigma$-algebra generated by $\left\{M_{0},\cdots,M_{k-1}\right\}$ and $\left\{F_{0}^{0},F_{0}^{d},\cdots,F_{k-1}^{0},F_{k-1}^{d}\right\}$. After the current model $m_{k}$ is constructed, let $\mathcal{F}_{k-1/2}^{M\cdot F}$ denote the $\sigma$-algebra generated by $\left\{M_{0},\cdots,M_{k}\right\}$, and $\left\{F_{0}^{0},F_{0}^{d},\cdots,F_{k-1}^{0},F_{k-1}^{d}\right\}$.

Next, we will introduce some definitions to precise our requirements on the probabilistic models and estimates.

2.1 Probabilistic models and estimates

First we recall the measure for accuracy of deterministic models, which is introduced in Random_model_1 ; Random_model_2 .

Definition 2.1**.**

We say a model $m_{k}$ is $\kappa$-fully linear model of f on $B(x_{k},\delta_{k})$, for $\kappa=(\kappa_{ef},\kappa_{eg})$, if   $\forall y\in B(x_{k},\delta_{k})$,

[TABLE]

The extending concept of the above definition is probabilistically fully-linear model which is described in STORM .

Definition 2.2**.**

A sequence of random model $M_{k}$ is said to be $\alpha$-probabilistically $\kappa$-fully linear with respect to the corresponding sequence $\left\{X_{k},\Lambda_{k}\right\}$, if the events

[TABLE]

satisfy the condition:

[TABLE]

The above definition states that the model $M_{k}$ is a locally good approximation of the first-order Taylor expansion of the objective function with probability at least $\alpha$, conditioned on $\mathcal{F}_{k-1}^{M\cdot F}$. However, there is still some possibility such that the model is inaccurate, even very bad. To guarantee the quality of the trial step, we hope the random model $M_{k}$ closer to the first-order Taylor expansion. However, the corresponding computation cost will increase. Thus there is a trade-off between the accuracy of the model and the computation cost.

Taking aside of the accurate model, the estimates of $f(x_{k})$ and $f(x_{k}+d_{k})$ are also required to be sufficiently accurate. The deterministic version of accurate estimates is formally stated as follows.

Definition 2.3**.**

The estimates $f_{k}^{0}$ and $f_{k}^{d}$ are $\epsilon_{F}$-accurate estimates of $f(x_{k})$ and $f(x_{k}+d_{k})$, if

[TABLE]

The definition of probabilistically accurate estimates is shown as follows which is a modified version of that in Random_estimate .

Definition 2.4**.**

A sequence of random estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ is said to be $\beta$-probabilistically $\epsilon_{F}$-accurate with respect to the corresponding sequence $\left\{X_{k},\Lambda_{k},D_{k}\right\}$, if the events

[TABLE]

satisfy the condition:

[TABLE]

where $\epsilon_{F}$ is a fixed constant.

Using Definitions 2.2 and 2.4, we assume that model $M_{k}$ and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ satisfy the following assumption in our analysis.

Assumption 1**.**

The followings hold for the quantities used in Algorithm 2.1

1. (i)

There exist $\kappa_{ef}$, $\kappa_{eg}>0$ such that the sequence of random models $M_{k}$ is $\alpha$-probabilistically $(\kappa_{ef},\kappa_{eg})$-fully linear, for a sufficiently large $\alpha\in(0,1)$. 2. (ii)

There exist $\epsilon_{F}>0$ such that the sequence of estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ is $\beta$-probabilistically $\epsilon_{F}$-accurate, for a sufficiently large $\beta\in(0,1)$. 3. (iii)

The sequence of estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ generated by Algorithm 2.1 satisfies the following condition that

[TABLE]

with $\kappa_{f}>0$.

Remark 2.1**.**

The above assumption is significant in the following convergence rates analysis. Compared to the Assumption 3.1 in STORM_nonconvex , we know that Assumption 1(iii) is additional. Nevertheless it is essential, which states that estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ can not be too worse in expectation in contrast to the true function values $\left\{f(X_{k}),f(X_{k}+D_{k})\right\}$ when estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are not sufficiently accurate. One may doubt that this additional condition is a little stronger than that of stochastic line search S_line_search . One would agree that (2.5) given above and (2.3) of S_line_search are close in spirit. However, it does not possible to obtain one from the another. Actually, the term $\Lambda_{k}\left\|\nabla f(X_{k})\right\|^{2}$ in (2.3) of S_line_search is hard to be calculated directly in the construction of estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$. From later analysis of Theorem 3.1, we can see that if Assumption 1(iii) is relaxed as Assumption 2.4 in S_line_search , the bound for $\mathbb{E}[\Phi_{k+1}-\Phi_{k}]$ will be worse than the current results. Moreover the estimates for $\nu$ and other parameters will be more complicated. Note that Assumption 1(iii) can easily be satisfied in the practical implementation as long as Assumption1(ii) holds. We will explain this later in our numerical experiments.

To make the analysis simple and easy to understand, we use the following statements.

• •

If $I_{k}=1$, we say that the model is true; otherwise, we say that the model is false.

• •

If $J_{k}=1$, we say that the estimates are tight; otherwise, we say that the estimates are loose.

• •

If an iteration $k$ is accepted, we say that the iteration is successful; otherwise, we say that the iteration is failed.

In the end of this subsection, we would like to introduce a definition of convergence criterion for analysis, which is named $\epsilon$-solution. When $f$ is unknown to be convex, we say that $X_{k}$ is an $\epsilon$-solution if $\left\|\nabla f(X_{k})\right\|\leq\epsilon$. However, when $f$ is convex or strongly convex, we say that $X_{k}$ is an $\epsilon$-solution if $f(X_{k})-f^{\ast}\leq\epsilon$.

Remark 2.2**.**

There are three cases of non-convex, convex and strongly convex to be discussed in this article. Due to the intractability of the general non-convex problem, it is unreasonable to use the same criterion as the convex problem. Thus the definition of $\epsilon$-solution for non-convex case is different from that for convex case.

2.2 Analysis of the stochastic process

In this part, we aim to estimate the upper bound of an expected stopping time by observing the behavior of stochastic process generated by Algorithm 2.1. The results can be applied to analyze the convergence rates of the proposed algorithm in different settings.

We first give some basic definitions before theoretical analysis.

Definition 2.5**.**

Let $X=\left\{X_{k},k\geq 0\right\}$ be a stochastic process. We say $T$ is a stopping time with respect to $X$ if for each $k>0$, the event $\left\{T=k\right\}$ is completely determined by the total information up to time $k$, that is $\left\{X_{0},X_{1},\cdots,X_{k}\right\}$.

Here we give a random variable $T_{\epsilon}$, which is the total number of iterations until an $\epsilon$-solution is achieved.

Remark 2.3**.**

$T_{\epsilon}$* is a special stopping time and dependent on randomness of the proposed algorithm and the $\epsilon$-solution we have defined.*

Next, we consider a stochastic process $\left\{\Lambda_{k},\Phi_{k}\right\}$ such that $\Lambda_{k}\in\left[0,\infty\right)$ and $\Phi_{k}\in\left[0,\infty\right)$ for all $k>0$. Let us give the definition of a special random event $W_{k}$ as follows

[TABLE]

where $p\in[0,1]$. We now assume that $\Lambda_{k}$ and $\Phi_{k}$ satisfy the assumption mentioned below, for all $k<T_{\epsilon}$.

Assumption 2**.**

2. (i)

There exist constants $\Phi_{\max}>0$ and $\mu_{\max}>0$ such that $\Phi_{k}\leq\Phi_{\max}$ and $\Lambda_{k}\leq\mu_{\max}$, respectively. 3. (ii)

There exists a constant $\hat{\Lambda}>0$ such that for all $k\leq T_{\epsilon}$, the following properties hold

[TABLE]

where $\lambda_{1}\in\mathbb{R}$ and $W_{k}$ satisfies (2.6) with $p>\frac{1}{2}$. 4. (iii)

There exists a constant $C>0$, and a non-decreasing function $h(\cdot)$ which is positive on any positive domain, such that for all $k<T_{\epsilon}$,

[TABLE]

Based on Assumptions 2, the following theorem (see STORM_nonconvex ) illustrates the upper bound on the expected number of iterations $T_{\epsilon}$ for obtaining an $\epsilon$-solution.

Theorem 2.1**.**

Let Assumption 2 holds. Then

[TABLE]

The analysis of the renewal-reward process in STORM_nonconvex is appropriate for the stochastic process generated by Algorithm 2.1. So here we omit the proof of Theorem 2.1. For more details, we refer the readers to Theorem 2.2 in STORM_nonconvex . Theorem 2.1 is very important to the following analysis of the complexity of Algorithm 2.1. The difficulties lie in finding the non-decreasing function $h(\cdot)$ and the constant $\hat{\Lambda}$. The choice of trust-region radius indeed introduces some differences and difficulties, compared to the analysis of STORM_nonconvex and S_line_search . In the next part, we will give more analysis and discussions.

3 Stochastic trust-region with probabilistic models and estimates

In this section, we propose a specific trust-region framework named STRME based on probabilistic models and estimates. The main steps are described as follows.

At each iteration $k$, given a current point $x_{k}$ and trust-region radius $\delta_{k}$, the model is built as

[TABLE]

to approximate $f(x)$ in $B(x_{k},\delta_{k})$. The quadratic model is simple and widely used in many trust-region algorithms. Of course, other models, for example the conic model (see Yuan_TR_review ), can also be applied to the framework as long as some requirements we stated are met.

The trust-region radius is defined as $\delta_{k}=\mu_{k}\left\|g_{k}\right\|$. Actually, one can try other more general choices that $\delta_{k}=\mu_{k}^{r_{1}}\left\|g_{k}\right\|^{r_{2}}$ with $r_{1},r_{2}\geq 0$ TR_Yuan_Trustregion . For simplicity, we only consider the case that $r_{1},r_{2}=1$. In the following steps, we choose to update the parameter $\mu_{k}$ (relative trust-region radius). Note that due to the randomness of the model, $\left\|g_{k}\right\|$ can be very small even zero even though the algorithm does not converge yet. In this case, it does not make any sense to continue the following process. Thus we add steps in Algorithm 3.1 to check if $\left\|g_{k}\right\|>\epsilon$.

The trial step $d_{k}$ is produced by minimizing the model $m_{k}(x_{k}+d)$ in a neighborhood of $x_{k}$ exactly or inexactly. Then we compute the random estimates $f_{k}^{0}$ and $f_{k}^{d}$ of $f(x_{k})$ and $f(x_{k}+d_{k})$ respectively to measure the actual function reduction. Once the trial step is obtained, we can use the ratio $\rho_{k}$, which is defined below, to judge how good the trial step $d_{k}$ is. Based on this criterion, if the trial point $x_{k}+d_{k}$ yields sufficient reduction, we accept the trial step $d_{k}$; otherwise, we reject it. At the end of each iteration, the relative trust-region $\mu_{k}$ is chosen according to the outcome of the iterates. The details of the algorithm are described as follows.

At each iteration, the trial step $d_{k}$ is computed to satisfy the well-known $Cauchy\,decrease$ condition, which is given as follows.

Assumption 3**.**

[TABLE]

Besides, in the case that $\delta_{k}=\mu_{k}\left\|g_{k}\right\|$, the condition $\left\|g_{k}\right\|\geq\eta_{2}\delta_{k}$ in Algorithm 2.1 is equivalent to the condition $\mu_{k}\leq\frac{1}{\eta_{2}}$ . We notice that the parameter $\eta_{2}$ in algorithm 2.1 is usually very small. When $\eta_{2}$ is small, we can see that the condition $\mu_{k}\leq\frac{1}{\eta_{2}}$ is easily satisfied. We might as well let $\mu_{\max}\leq\frac{1}{\eta_{2}}$, thus the condition $\left\|g_{k}\right\|\geq\eta_{2}\delta_{k}$ in Algorithm 2.1 can be satisfied automatically.

3.1 Theoretical properties of STRME

We are ready to present the theoretical properties of the framework described in Algorithm 3.1. First, we give an assumption which states that the Hessian approximation matrix $B_{k}$ in model $m_{k}$ is uniformly upper bounded.

Assumption 4**.**

There exists a constant $\kappa_{bhm}>0$ such that, for all $k\geq 0$,

[TABLE]

We now provide some auxiliary lemmas to show that the decrease of the objective function $f(x)$ is guaranteed under some conditions. The following lemma states that if the model $m_{k}$ is true (fully linear) and the relative trust-region radius $\mu_{k}$ is upper bounded by a given number, then the actual reduction of the objective function is achieved. Although the theoretical content of Lemmas 3.1 to 3.4 closely relies on the existing arguments related to the STORM algorithm, for the integrity of analysis, the proofs of these lemmas will still be attached in the Appendix.

Lemma 3.1**.**

Suppose that model $m_{k}$ is true. If

[TABLE]

then

[TABLE]

The next lemma states that the decrease of the objective function is achieved if the estimates $\left\{f_{k}^{0},f_{k}^{d}\right\}$ are tight and iteration $k$ is successful.

Lemma 3.2**.**

Suppose that estimates $f_{k}^{0}$ and $f_{k}^{d}$ are tight with $\epsilon_{F}\leq\frac{\eta_{1}\kappa_{fcd}}{8\mu_{\max}}$ and $\mu_{k}\leq\frac{1}{\kappa_{bhm}}$. If $d_{k}$ is accepted, then

[TABLE]

The following lemma shows when model $m_{k}$ and estimates $\left\{f_{k}^{0},f_{k}^{d}\right\}$ are both sufficiently accurate, if $\mu_{k}$ is not too large, then the iteration will be successful.

Lemma 3.3**.**

Suppose that model $m_{k}$ is true, and estimates $f_{k}^{0}$ and $f_{k}^{d}$ are tight with $\epsilon_{F}\leq\kappa_{ef}$. If

[TABLE]

then the k-th iteration is successful.

We now turn to consider the random process $\left\{\Phi_{k},\Lambda_{k}\right\}$ derived from the process generated from Algorithm 3.1. The following analysis is based on the function

[TABLE]

for some $\nu\in(0,1)$. It is obvious that $\Phi_{k}\geq 0$. Actually, the random variable $\Phi_{k}$ can be regarded as a kind of measure of progress to optimality. It plays an important role in the analysis of such trust-region algorithms. In Algorithm 3.1, we update the trust-region radius as $\Delta_{k}=\Lambda_{k}\left\|G_{k}\right\|$. One may find that the trust-region radius $\Delta_{k}$ depends on the randomness introduced by the current model $M_{k}$, that is to say, $\Delta_{k}$ is not measurable with respect to $\mathcal{F}_{k-1}^{M,F}$. However, in STRME, we update the parameter $\Lambda_{k}$, which is completely determined by $\mathcal{F}_{k-1}^{M,F}$, instead of $\Delta_{k}$. Our approach avoids the difficulty of using the trust-region radius directly to analyze the complexity of STRME.

We have to admit that the theoretical analysis of STRME algorithm is similar to that of S_line_search in spirit, but there are some distinctions between them. The most distinctive lies in the measure function $\Phi_{k}$. In S_line_search , $\Phi_{k}$ consists of the three terms $f(X_{k})-f^{\ast}$, $\frac{1}{L^{2}}\Lambda_{k}\left\|\nabla f(X_{k})\right\|^{2}$ and $\Delta_{k}^{2}$. However, in our analysis, only two of them, that is $f(X_{k})-f^{\ast}$ and $\frac{1}{L^{2}}\Lambda_{k}\left\|\nabla f(X_{k})\right\|^{2}$, are used to evaluate the reduction of the algorithm. Actually, for STRME algorithm, trust-region $\Delta_{k}$ does not necessarily increase even if the iteration is successful.

We aim to bound the expected number of iterations $\mathbb{E}[T_{\epsilon}]$. Before that, we have to prove Assumption 2 holds for the process $\left\{\Phi_{k},\Lambda_{k}\right\}$. It is apparent that Assumption 2(i) holds with the definition of $\Phi_{k}$ and $\mu_{k}\leq\mu_{\max}$ in Algorithm 3.1. This assumption is not related with the convexity of the objective function, so it holds in all three cases we will consider later. Let us define the constant $\hat{\Lambda}$ in Assumption 2(ii) as follows:

[TABLE]

In our analysis, we might as well claim that $\mu_{\max}\leq\min\left\{\frac{1}{\kappa_{bhm}},\frac{\kappa_{fcd}}{8\kappa_{ef}}\right\}.$ For simplicity, we assume that $\Lambda_{0}=\gamma^{i}\hat{\Lambda}$ and $\mu_{\max}=\gamma^{j}\hat{\Lambda}$ for some integers $i,j>0$. As a result, for any $k>0$, we have $\Lambda_{k}=\gamma^{i}\hat{\Lambda}$ for some integer $i$. Next, we will show that Assumption 2(ii) holds provided the constant $\hat{\Lambda}$ is defined as above.

Lemma 3.4**.**

Let $\alpha$ and $\beta$ safisfy that $\alpha\beta>\frac{1}{2}$, then we have Assumption 2(ii) holds with $W_{k}=2(I_{k}J_{k}-\frac{1}{2})$, $\lambda_{1}=\log(\gamma)$, and p = $\alpha\beta$.

Using the above Lemmas 3.1 to 3.3, we can derive the following result.

Theorem 3.1**.**

Let Assumptions 1 hold with $\epsilon_{F}\leq\min\left\{\frac{\eta_{1}\kappa_{fcd}}{8\mu_{\max}},\kappa_{ef}\right\}$ and $\kappa_{f}\leq 2\eta_{1}\kappa_{ef}$. Besides we assume that $f$ is $L$-smooth and Assumption 4 is satisfied. Then there exist a constant $\nu$ and sufficiently large $\alpha,\beta$ satisfying the following conditions

[TABLE]

[TABLE]

such that

[TABLE]

Remark 3.1**.**

Due to the possibility of inaccurate model $M_{k}$ and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$, the function value of $f$ may increase. So $\Phi_{k}$ is designed to balance the decrease and increase of $f(X_{k})$. The above theorem shows that the decrease of expected $\Phi_{k}$ can be achieved by carefully choosing $\nu$ and the probability $\alpha$ and $\beta$.

3.2 Convergence rates for non-convex problems

We now show the global convergence rate of Algorithm 3.1 when $f$ is unknown to be convex, that is the following assumption holds.

Assumption 5**.**

$f\in\mathcal{C}^{1}(\mathbb{R}^{d})$* is bounded below by $f^{\ast}$ and $L$-smooth.*

Our goal is to bound the expected number of iterations until an $\epsilon$-solution occurs, i.e. $\mathbb{E}[T_{\epsilon}]$. The definition of $T_{\epsilon}$ is described as follows

[TABLE]

For all $k\leq T_{\epsilon}$, we know $\left\|\nabla f(X_{k})\right\|\geq\epsilon$. Let us recall the definition of $\Phi_{k}$, i.e.

[TABLE]

In this case, $h(\cdot)$ can be defined as

[TABLE]

where $C=\frac{1}{2}\alpha\beta(1-\nu)(1-\frac{1}{\gamma})$, which is non-decreasing on any positive domain. From Theorem 3.1, we know that Assumption 2(iii) will hold if the conditions in Theorem 3.1 are satisfied. Applying the result in Lemma 3.4, Assumption 2(ii) holds under the conditions that $\alpha\beta>\frac{1}{2}$ and $\hat{\Lambda}$ is defined as (3.6). Clearly Assumption 2(i) holds. Thus, we can conclude that Assumption 2 holds under certain conditions. Thus the conclusion of Theorem 2.1 is true in this case. Then the following complexity result for Algorithm 3.1 can be achieved by simple analysis.

Theorem 3.2**.**

Suppose that $f$ satisfies Assumption 5. Under the conditions in Theorem 3.1, if $\alpha\beta>\frac{1}{2}$, then for Algorithm 3.1, to achieve $\epsilon$-solution, the expected number of iterations is bounded as follows

[TABLE]

where $M=\frac{2\nu(f(x_{0})-f^{\ast})+2(1-\nu)\frac{1}{L^{2}}\mu_{0}\left\|\nabla f(x_{0})\right\|^{2}}{\alpha\beta(1-\nu)(1-\frac{1}{\gamma})\hat{\Lambda}}$.

Remark 3.2**.**

Note that the dependency of the expected number of iterations for obtaining an $\epsilon$-solution is $\mathcal{O}(1/\epsilon^{2})$, which is similar to that in STORM_nonconvex for nonconvex problems.

3.3 Convergence rates for convex problems

In this part, we will analyze the expected complexity for STRME when $f$ is convex. First we give the following assumption.

Assumption 6**.**

$f\in\mathcal{C}^{1}(\mathbb{R}^{d})$* is convex. The level set $\mathcal{L}=\left\{x\in\mathbb{R}^{d}:f(x)\leq f(x_{0})\right\}$ is bounded, and there exists a constant $D>0$ such that*

[TABLE]

In convex setting, we aim to bound the expected $T_{\epsilon}$, which is defined as below

[TABLE]

for an $\epsilon$-solution. In this case, we define a function

[TABLE]

to replace $\Phi_{k}$ to measure the progress of the iterations. For $k\leq T_{\epsilon}$, $f(X_{k})-f^{\ast}\geq\epsilon$, then we have $\Phi_{k}\geq\nu\epsilon$. Thus $\Psi_{k}\geq 0$, which is well-defined. In the later analysis, we can demonstrate that the random process $\left\{\Lambda_{k},\Psi_{k}\right\}$ satisfies Assumption 2. Therefore, Theorem 2.1 can be applied to derive the upper bound of $\mathbb{E}[T_{\epsilon}]$. For more details, please refer to the proof of Theorem 3.3 in the Appendix.

Theorem 3.3**.**

We assume that $f$ is $L$-smooth and satisfies Assumption 6. Under the conditions in Theorem 3.1 and $\alpha\beta>\frac{1}{2}$, for Algorithm 3.1, in order to achieve an $\epsilon$-solution, we have

[TABLE]

where $M=\frac{2(\nu L+(1-\nu)\mu_{\max})^{2}D^{2}}{\alpha\beta\nu(1-\nu)(1-\frac{1}{\gamma})\hat{\Lambda}}.$

Remark 3.3**.**

The above theorem states that if $f$ is convex, the dependency on $\epsilon$ for the complexity bound is $\mathcal{O}(1/\epsilon)$, which is a stronger result than that in Theorem 3.2. Note that the above result is the same as that of traditional trust-region for general convex problems.

3.4 Convergence rates for strongly convex problems

In this part, we will derive the complexity bound for Algorithm 3.1 in the strongly convex setting. In this case, we assume the following assumption holds.

Assumption 7**.**

$f\in\mathcal{C}^{1}(\mathbb{R}^{d})$* is strongly convex, i.e. there exists a constant $\sigma>0$ such that*

[TABLE]

First we use the definition of $T_{\epsilon}$ as given in (3.14):

[TABLE]

Our aim is to bound the expected number of iteration i.e. $\mathbb{E}[T_{\epsilon}]$, to obtain an $\epsilon$-solution for strongly convex problems. Instead of using $\Phi_{k}$ to measure the progress towards optimality, we define a function

[TABLE]

For $k\leq T_{\epsilon}$, $f(X_{k})-f^{\ast}\geq\epsilon$, then we have $\Phi_{k}\geq\nu\epsilon$, which implies that $\Psi_{k}\geq 0$. So the definition of $\Psi_{k}$ is reasonable. In the later part, we will obtain the upper bound of $\mathbb{E}[T_{\epsilon}]$ with the help of random process $\left\{\Lambda_{k},\Psi_{k}\right\}$. The details will be shown in the proof of the following theorem.

Theorem 3.4**.**

Assume that $f$ is $L$-smooth and satisfies Assumption 7. If the conditions in Theorem 3.1 and $\alpha\beta>\frac{1}{2}$ hold, then for Algorithm 3.1, in order to achieve an $\epsilon$-solution, we have

[TABLE]

where $M=\frac{(\frac{\nu}{2\sigma}+(1-\nu)\frac{1}{L^{2}}\mu_{\max})}{\alpha\beta(1-\nu)(1-\frac{1}{\gamma})\hat{\Lambda}}.$

Similarly, the proofs of the above theorem will be given in the Appendix.

Remark 3.4**.**

The result in Theorem 3.4 show Algorithm 3.1 takes at most $\mathcal{O}(\log(1/\epsilon))$ iterations in expectation, to achieve an $\epsilon$-solution. The result coincides with that for trust-region methods in deterministic setting.

4 Numerical experiments

In this section, we empirically test our STRME algorithm and compare its performance with STORMSTORM and some related algorithms.

We test on two type of problems: (i) regularized logistic regression problem, which is strongly convex; (ii) deep neural networks, which is highly non-linear and non-convex. The function value of training data (called training loss) and accuracy (percentage of correctly classified testing data) are adopted as criteria to measure the performance of all the algorithms that are tested. For all those algorithms, we compare these criteria against the number of effective pass through the data, that is total gradient calls divided by $N$(training data size). All algorithms were terminated when the maximum budget of the gradient evaluations is larger than the maximum value $SFO_{\max}$ we have set.

All algorithms are implemented in Anaconda3 (python 3.6.2) under Windows 7 operating system on Dell desktop with Intel(R) Core(TM) i7-4790U CPU @3.6GHz, 8GB Memory.

4.1 How to implement the probabilistic models and estimates to satisfy Assumption 1

In these numerical experiments, we focus on the derivative-based problems where $f(x,\omega)$ and $\nabla f(x,\omega)$ are available. At each iteration point $x_{k}$, we assume that the noise $\omega$ in function value and gradient computation is conditional unbiased and the corresponding variance is conditional bounded for all $f$, i.e.

[TABLE]

We now discuss how to obtain $\alpha$-probabilistically $\kappa$-fully linear models and $\beta$-probabilistically $\epsilon_{F}$-accurate estimates. In the later analysis, we employ the standard sample averaging approximation technique Sample_average to construct the sufficient accurate model and estimates. Let

[TABLE]

where $\omega_{i}$ are the i.i.d. and finite realizations of the noise $\omega$ and the sample set $S_{k}\subseteq\left\{1,2,\cdots,N\right\}$ of size $\left|S_{k}\right|=p_{k}$. The local approximation model can be constructed as $m_{k}(x)=f_{k}+g_{k}^{T}(x-x_{k})+\frac{1}{2}(x-x_{k})^{T}B_{k}(x-x_{k})$, where $x\in B(x_{k},\delta_{k})$. We recall the definitions in Section 2.1. In order to satisfy Assumption 1(i),we need the following conditions hold

[TABLE]

for all $x\in B(x_{k},\delta_{k})$. By Chebyshev’s inequality in Lemma 6.2, at current point $x_{k}$, for any $v>0$, we have

[TABLE]

By (4.4), the conditions (4.3) hold at $x=x_{k}$, as long as the sample rate $p_{k}$ satisfies the following condition that

[TABLE]

For $\forall x\in B(x_{k},\delta_{k})/\left\{x_{k}\right\}$, if $\left\|g_{k}-\nabla f(x_{k})\right\|\leq\kappa_{eg}\delta_{k}$, we have

[TABLE]

The second inequality relies on the fact that $\left\|x-x_{k}\right\|\leq\delta_{k}$, $f$ is $L-$smooth, and the second-order matrix $B_{k}$ satisfies Assumption 4. Thus we can conclude that if the sample rate $p_{k}$ satisfies (4.5), condition (4.3) will hold for any $x\in B(x_{k},\delta_{k})$. Consequently, $m_{k}(x)$ is $\alpha$-probabilistically $\kappa$-fully linear models. In practice, the function value is not explicitly computed at model $m_{k}$, so we only require that $p_{k}\geq\frac{V_{g}}{(1-\alpha^{\prime})\kappa_{eg}^{2}\delta_{k}^{2}}.$

Similarly, to obtain $\beta$-probabilistically $\epsilon_{F}$-accurate estimates, let estimates $\left\{f_{k}^{0},f_{k}^{d}\right\}$ be

[TABLE]

where $S_{k}^{0}\subseteq\left\{1,2,\cdots,N\right\}$ and $\left|S_{k}^{0}\right|=q_{k}$. In order to satisfy Assumption 1(ii) such that

[TABLE]

we require

[TABLE]

By the Hölder’s inequality for expectation that $\mathbb{E}[ab]\leq(\mathbb{E}[a^{2}])^{1/2}(\mathbb{E}[b^{2}])^{1/2}$ for $a,b>0$, we have

[TABLE]

Since $q_{k}\geq\frac{V_{f}}{(1-\beta)\epsilon_{F}^{2}\delta_{k}^{4}}$, we have $\mathbb{E}[\frac{\left|f_{k}^{0}-f(x_{k})\right|^{2}}{\epsilon_{F}^{2}\delta_{k}^{4}}|\mathcal{F}_{k-1/2}^{M,F}]\leq\frac{V_{f}}{q_{k}\epsilon_{F}^{2}\delta_{k}^{4}}\leq\mathcal{O}(1)$. Thus there exists a constant $\kappa_{f}>0$ such that $\mathbb{E}[\left|f_{k}^{0}-f(x_{k})\right||\mathcal{F}_{k-1/2}^{M,F}]$ is bounded by $\kappa_{f}\delta_{k}^{2}$. In the same way, similar result can be applied to $f_{k}^{d}-f(x_{k}+d_{k})$. We can claim that Assumption 1(iii) holds. Thus, we have shown how to construct the random model $m_{k}$ and estimates $\left\{f_{k}^{0},f_{k}^{d}\right\}$ to satisfy Assumption 1.

However, we can see that the trust-region radius $\delta_{k}$ is unknown before the sample rate is chosen as in (4.5) and (4.9). It brings a challenge to construct a fully linear model in practice. First of all, obviously, in the case that the total sample set is large but limited, such sample rate $p_{k}$ must exist, such as $p_{k}=N$, to make the model fully linear. But it is really difficult to give a computable condition of $p_{k}$ in theory, and prove that under such condition the model is fully linear with probability. Thus we settle for the minimal sample rate of increase in $p_{k}$ as byrd2012sample . As previous discussed in (4.5), the sample size $p_{k}$ satisfies

[TABLE]

In order to guarantee the fully linear property, from the analysis in Section 3, the critical is that there exists a constant $G_{1}>0$ such that the following condition holds

[TABLE]

The above condition yields

[TABLE]

In byrd2012sample , under the condition (4.12), they have analyzed gradient-based mini-batch optimization algorithm for the strongly convex problem that the sample size should grow geometrically with iteration $k$. Intuitively, it make sense because if condition (4.12) holds, $\left\|\nabla F(x_{k})\right\|^{2}$ is geometrically decreasing in this case, so does $\left\|g_{k}\right\|^{2}$. Theorem 3.4 shows that our algorithm achieves the same complexity rate as Theorem 4.1 in byrd2012sample . By (4.11), in the strongly convex case, we impose the sample rate $p_{k}$ to be exponentially increased, i.e. $p=a^{k}$ for some $a>1$. However, for the non-convex problem, the situation is different. From the analysis in Section 3.2, to obtain an $\epsilon$-solution ($\left\|\nabla F(x_{k})\right\|\leq\epsilon$), we have $T_{\epsilon}\approx\mathcal{O}(\frac{1}{\epsilon^{2}})$. From (4.11) and (4.13), the sample rate $p_{k}$ should be linearly increasing with iteration $k$.

As we know, the situation in practice is very complicated. Thus we require that the sample rate grows exponentially with iteration $k$ for regularized logistic regression problem and linearly increases with iteration $k$ for the deep neural networks, respectively.

4.2 Experimental results on regularized logistic regression problem

In this subsection, we consider the following smooth (strongly convex) regularized logistic loss problem considered in STORM :

[TABLE]

where $\left\{(a_{i},b_{i})\right\}_{i=1}^{N}$ is a training sample set with $a_{i}\in\mathbb{R}^{d}$ being the feature vector and $b_{i}\in\left\{-1,+1\right\}$ being the corresponding label. And $\lambda\geq 0$ is the regularization parameter. As in the typical machine learning setting, we assume that $n$ is very large and $N\gg d$. So computing $F(x)$, as well as $\nabla F(x)$ and $\nabla^{2}F(x)$ are very expensive. In our work, we randomly (without replacement) choose a subset $I_{k}\subseteq\left\{1,2,\cdots,N\right\}$ to estimate the quantities in our algorithms. For the algorithms that need to compute Hessian matrix, the same sample is drawn for gradient and Hessian evaluations. In this setting, we re-sample the sample set for $f_{k}^{0}$ and $f_{k}^{d}$.

We compare our algorithm STRME with STORM which is implemented as algorithm 5 in STORM . In our numerical experiments, we construct two versions of STRME: one which only computes the stochastic gradients and sets $B_{k}=0$ is the first-order version (called STRME-1st), the other one which in addition to the stochastic gradients, computes stochastic Hessian estimators is the second-order version (called STRME-2st). We use the classic dogleg method in Numerical-Op to solve the second order subproblem , and the corresponding algorithms, we call STRME-dogleg and STORM-dogleg. For the details of the implementation, one can refer to Algorithm 6.1 in Appendix. Besides, We compare against a special adaptive solver AdaGrad AdaGrad , which takes the adaptive step size but does not have to compute function value, and only computes the average stochastic gradients.

For problem (4.14), all algorithms were tested with different input parameters. We set $x_{0}=0$ as the starting point for all algorithms. For the three quantities gradient, Hessian and ($f_{k}^{0}$, $f_{k}^{d}$) evaluations, we adopt the linearly increased sample rule that $b_{k}=\min\left\{b_{\max},\max\left\{t_{0}k+b_{0},\frac{1}{\delta^{2}}\right\}\right\}$, where $b_{0}=d+1$, $b_{\max}=N$, $t_{0}=100$ for STORM and our STRME algorithm as in STORM . For AdaGrad, the mini-batch size $b=d+1.$ We set the same random seed to generate random sample sequence for all the algorithms.

The regularized logistic regression problem we consider in this subsection is strongly convex. As we discussed in section 4.1, to make sure that the model is fully linear, the sample size should be exponentially increased. Thus we also tested the mini-batch size $b_{k}=\min\left\{b_{\max},\max\left\{b_{0}a^{k},\frac{1}{\delta^{2}}\right\}\right\}$, where $b_{0}=d+1$, $b_{\max}=N$, $t_{0}=100$, and $a\in\left\{1.1,1.3,1.5,2\right\}$ for our algorithm STRME and STORM. And we compare the two sample rules as well at the end of this part.

In this subsection, we test on two datasets a9a and ijcnn1 from the LIBSVM website 222https://www.csie.ntu.edu.tw/cjlin/libsvmtools/datasets/. We list the datasets in Table 1, in which $n$ denotes the total sample size, and $d$ is the dimension of the dataset, and $\lambda$ is the regularization parameter. We use 0.95 partition of the data as the training set, and the remaining as the testing set, just like in STORM .

In Figure 1, we present the results on the a9a dataset. For STORM-1st and STORM-2st, we use the same parameters in STORM : $\delta_{\max}=10,\delta_{0}=1,\gamma=2,\eta_{1}=0.1,\eta_{2}=0.001$. For STRME-1st and STRME-2st, the following parameters are used: $\mu_{\max}=10^{3},\mu_{0}=1,\eta_{1}=0.1,\gamma=2$. We set step size $\eta=1$ for AdaGrad.

In Figure 2, we report the results on the ijcnn1 dataset. The parameters $\mu_{0}=10,\mu_{\max}=10^{3},\eta_{1}=0.1,\gamma=2$ are set for STRME-1st and $\mu_{0}=1$, others are the same for STRME-dogleg. For STORM-1st and STORM-dogleg, we set $\delta_{0}=1,\delta_{\max}=10,\eta_{1}=0.1,\gamma=2$. For AdaGrad, we set step size $\eta=1$.

From Figure 1 and 2, we can find that the proposed STRME is comparable to STORM in this setting, both in terms of the training function value and accuracy.

At the end of this part, we compare the exponentially increased sample rule with the linearly increased case in Figure 3. We choose the exponential ratio $a=1.1$ from $\left\{1.1,1.3,1.5,2\right\}$ for all the experiments. The other parameters are the same as the linearly increased sample. For the sake of distinction, we use STRME-L and STRME-E to denote the linearly increased and exponentially increased sample rules respectively, so does STORM. As a whole, from Figure 3, we can see that the performance on the two sample rules is very similar.

4.3 Experimental results on a simple deep neural network(DNN)

In this subsection, we consider to train a fully-connected 2-layer net with 50 hidden units (784-50-10) neural networks with MNIST333http://deeplearning.net/data/mnist/, a benchmark dataset of handwritten digits. We used softmax output, sigmoid hidden functions, and the cross-entropy error function. The $l_{2}$ regularization parameter $\lambda=10^{-3}$, suggested in SVRG-nonconvex .

We compare our algorithm STRME with STORMSTORM . As in the previous subsection, we construct two versions of STRME: one which only computes the stochastic gradients and set $B_{k}=0$ is the first-order version, the other one which in addition to the stochastic gradient, computes the quasi-Newton matrix $B_{k}$ to approximate the true Hessian matrix is the second-order version. Besides, we make an implementation of the SdLBFGS SdLBFGS , which is an efficient second-order algorithm for the non-convex problem, to compare with STRME and STORM.

In our work, we run one epoch mini-batch SGD algorithm to obtain an initial point for all algorithms. In our implementation, we randomly (without replacement) choose the subset $I_{k}\subseteq\left\{1,2,\cdots,n\right\}$ to estimate the gradient and Hessian pair $(g_{k},B_{k})$ and objective function value pair ($f_{k}^{0}$, $f_{k}^{d}$). We have attempted the three cases: (i) sample gradient and Hessian pair, and re-sample $f_{k}^{0}$ and $f_{k}^{d}$ independently; (ii) sample the two pair independently; (iii) sample the two pairs with the same subset. However, for the first two cases, the results are not satisfactory. Therefore, in our numerical experiments, we only consider the last one. For STORM and STRME, we set the mini-batch size $b_{k}=\min\left\{b_{\max},\max\left\{t_{0}k+b_{0},\frac{1}{\delta^{2}}\right\}\right\}$, where $b_{0}=d+1$, $b_{\max}=N$, $t_{0}=10$.

4.3.1 Experimental results on the first-order probabilistic model

In this part, we first construct a first-order probabilistic model $m_{k}$ at $x_{k}$, i.e.

[TABLE]

where $f_{k}$ and $g_{k}$ are computed as (4.2), to test the STRME (STRME-1st) framework. Actually, in practice, we do not need to compute $f_{k}$. In our numerical experiments, we compare the numerical performace of STRME-1st with the related STORM-1st.

We now give details of parameters in the proposed STRME-1st and STORM-1st. For STRME-1st, $\mu_{0}$ and $\mu_{\max}$ are two important parameters. The parameters $\mu_{0}$ is chosen from $\left\{0.01,0.1,1\right\}$, and the best $\mu_{0}$ is achieved at $\mu_{0}=0.1$. Compared to $\mu_{0}$, the parameter $\mu_{\max}$ is more important in the later iteration process. Therefore the range of $\mu_{\max}$ is more elaborate. In our implementation, let $\mu_{\max}\in\left\{1/2,1,2,2^{2},2^{3}\right\}$, and the best tuned $\mu_{\max}$ is achieved at $\mu_{\max}=2$. For STORM-1st, we test on $\delta_{0}\in\left\{0.01,0.1,1\right\}$ and $\delta_{\max}\in\left\{1/2,1,2,2^{2},2^{3}\right\}$, the best performance is achieved with the inputs $(\delta_{0},\delta_{\max})=(0.1,1)$. The results are showing in the Figure 4.

From Figure 4, we can see that STRME-1st is better than STORM-1st with the model constructed as the beginning of this section. Not surprisingly, the trust-region radius constructed by STRME can make better use of gradient information. However, we can construct a more efficient second-order model to make the STRME algorithm perform better.

4.3.2 Experiments on second-order probabilistic models

In this part, we construct a specific second-order probabilistic model and implement STRME algorithm framework in deep neural network. In our work, we consider the limited memory symmetric rank one method (L-SR1) to generate the second order quasi-Newton matrix $B_{k}$, and build up the second-order model $m_{k}(x_{k}+d)$ as follows:

[TABLE]

where $f_{k}$ and $g_{k}$ are defined as the beginning of this section. Next, we will show how to update the quasi-Newton matrix $B_{k}$.

Let $s_{k}\stackrel{{\scriptstyle\Delta}}{{=}}x_{k+1}-x_{k}$ and $y_{k}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{b_{k}}(\sum_{j=1}^{b_{k}}\nabla f_{i_{j}}(x_{k+1})-\nabla f_{i_{j}}(x_{k})).$ Given a initial matrix $B_{0}$, provided that $(y_{k}-B_{k}s_{k})^{T}s_{k}\neq 0$, then $B_{k+1}$ can be defined as

[TABLE]

Limited memory symmetric rank one (SR1) method stores and uses the $m$ most recently computed pair $\left\{(s_{k},y_{k})\right\}$ ($m\ll d$). To describe the compact representation of a L-SR1 matrix, we need to define, for $k\geq m$:

[TABLE]

Moreover, we need the following decomposition of $S_{k}^{T}Y_{k}:$

[TABLE]

where $L_{k}$ is strictly lower triangular, $D_{k}$ is diagonal, and $R_{k}$ is strictly upper triangular. We assume that all the updates are well-defined, that is $s_{k}^{T}(y_{k}-B_{k}s_{k})\neq 0$, otherwise we skip the update. The compact form of L-SR1 Compact_LSR1 can be written as

[TABLE]

where $U_{k}\in\mathbb{R}^{d\times(k+1)}$, $V_{k}\in\mathbb{R}^{(k+1)\times(k+1)}$, and $B_{0}$ is a diagonal matrix. $U_{k}$ and $V_{k}$ are given by

[TABLE]

Now the quadratic probabilistic model defined by the L-SR1 method is constructed. The trust-region subproblem will be

[TABLE]

Applying this model into STRME, we can obtain STRME-Lsr1 algorithm. In the same way, we can obtain STORM-Lsr1 algorithm. With respect to the specific implementation of the subproblem and how to solve the trust-region subproblem efficiently, one can refer to the OBS method in OBS .

In this part, we compare our algorithm STRME-Lsr1 with STORM-Lsr1 and its first-order form. The result is shown in Figure 5. We test on different choice of matrix $B_{0}$: (i) $B_{0}=\frac{s_{0}^{T}y_{0}}{s_{0}^{T}s_{0}}$; (ii) $B_{0}=\frac{y_{0}^{T}y_{0}}{s_{0}^{T}y_{0}}$; (iii) $B_{0}=\tau_{0}\mathbb{I}_{d}$, where $\tau_{0}$ is adjusted with the same scale as the initial value of $\frac{y_{0}^{T}y_{0}}{s_{0}^{T}y_{0}}$. Unfortunately, the first two choices are not satisfactory. Thus, we use (iii) to define the initial matrix $B_{0}$.

In our numerical experiement, $B_{0}=\tau_{0}\mathbb{I}_{d}$, where $\tau_{0}\in\left\{0.5,1,1.5\right\}$, and the best $\tau_{0}$ is achieved at $\tau_{0}=1$. We have tested the limited memory size $m\in\left\{10,20,30,40\right\}$. We find that $m=30$ is performed relatively better. The parameters $\mu_{0}$ is chosen from $\left\{0.01,0.1,1\right\}$, and the best $\mu_{0}$ is achieved at $\mu_{0}=0.1$. For $\mu_{\max}$, we test on the range $\left\{0.1,1,10,100\right\}$, the best one is $\mu_{\max}=10$, which implies that unlike the first-order STRME, the second-order STRME is not so sensitive to $\mu_{\max}$ and allow some larger values. For STORM-Lsr1, we test the $\delta_{0}$ and $\delta_{\max}$ with the same range as $\mu_{0}$ and $\mu_{\max}$, the best choices are achieved at $\delta_{0}=0.1,\delta_{\max}=1$, respectively. From Figure 5, we can see that our algorithm STRME-Lsr1 performs better than STORM-Lsr1. And it is not hard to find that STRME-1st is not worse than STORM-Lsr1.

Moreover, in Figure 6, we show the behavior of successful and failed ratio $\varsigma$ which denotes the total number of successful iterations divided by the total number of failed iterations. In this case, we set the maximum number of gradient $SFO_{\max}=100N.$ For the first-order methods, we can see that the ratio $\varsigma$ of STORM-1st is basically stable around 1, and the ratio of STRME-1st is higher and increasing in the later period. Besides, we can see that in the middle and later period of the algorithm STRME-Lsr1 and STORM-Lsr1, the value of $\varsigma$ is basically stable. However, for STRME-Lsr1, the value of $\varsigma$ in stable condition is still larger than that for STORM-Lsr1.

Beyond that we also compare the trust-region radius $\delta_{k}$ of the two algorithms with a long time training in Figure 7, where $SFO_{\max}=100N$. We can see that for STRME, whether STRME-1st or STRME-Lsr1, the bandwidth of trust-region radius is narrower than that of STORM. This means the oscillation of the trust-region radius for STRME is less severe in contrast to that of STORM methods. Besides, for STRME-Lsr1, $\delta_{k}$ is overall declining and smaller than that in STORM-Lsr1. At the same time, we observe that the second-order methods permit larger trust-region radius from the Figure 7.

Moreover, we implement the well-known second-order algorithm SdLBFGS SdLBFGS to test the performance of our algorithm with a long time running ($SFO_{\max}=100N$). In this case, we set batch size $b=d+1$ and step size $\eta_{k}=\frac{\beta_{0}}{(k/10+1)}$ where $\beta_{0}\in\left\{0.1,1,10\right\}$ for SdLBFGS. The best tuned step size is obtained at $\beta_{0}=10$. In SdLBFGS , they set $\eta_{k}=\frac{10}{k+1}$. By numerical comparison, we find that the result of $\eta_{k}=\frac{10}{(k/10+1)}$ is better than that for the choice $\eta_{k}=\frac{10}{k+1}$. The results in Figure 8 illustrate that STRME-Lsr1 performs better than the best tuned SdLBFGS.

5 Conclusion

We have presented a stochastic trust-region framework in which the trust-region radius depends on the currently probabilistic model. To verify the effectiveness of the framework, we have proposed a specific algorithm named STRME in which the trust-region radius is linearly associated with the model gradient. We have analyzed the expected number of iterations of STRME for three different cases: non-convex, convex, and strongly convex. We can see that our algorithm enjoys the same complexity properties as the existing schemes. Moreover, our algorithm compares favorably to STORM algorithm and other stochastic algorithms on several testing problems involving the real datasets. Actually, in addition to STRME, there are many other approaches to explore the trust-region radius related with random models. We point out that the work in this paper is limited to the case that the objective function is smooth. There are some important and latest works for non-smooth problems, for example RSPG ; prox-SVRG-nonconvex ; QuickeNing ; StSR1 . It is worthwhile to extend the stochastic trust-region framework to the non-smooth setting. Moreover, the effectiveness of the stochastic trust-region method is relevant to the model. How to construct a more efficient model is an interesting subject for future research.

6 Appendix

A: Proofs of lemmas and theorems in Section 3

Proof of Lemma 3.1

Proof.

From Assumption 3, the trial step $d_{k}$ will lead a sufficient reduction on $m_{k}$ such that

[TABLE]

Since $\mu_{k}\leq\frac{1}{\kappa_{bhm}}$, we have $\delta_{k}=\mu_{k}\left\|g_{k}\right\|\leq\frac{\left\|g_{k}\right\|}{\kappa_{bhm}}\leq\frac{\left\|g_{k}\right\|}{\left\|B_{k}\right\|}$ which together with (6.1), implies that

[TABLE]

Suppose that model $m_{k}$ is true, we can obtain that

[TABLE]

Because of the condition that $\mu_{k}\leq\frac{\kappa_{fcd}}{8\kappa_{ef}}$, we have

[TABLE]

Thus, the desired result is proved. ∎

Proof of Lemma 3.2

Proof.

If $d_{k}$ is accepted, which implies that $\rho_{k}\geq\eta_{1}$, then

[TABLE]

where the last inequality follows from $\mu_{k}\leq\frac{1}{\kappa_{bhm}}$.

If the estimates $\left\{f_{k}^{0},f_{k}^{d}\right\}$ are tight, the improvement in $f$ can be bounded by

[TABLE]

Since $\epsilon_{F}\leq\frac{\eta_{1}\kappa_{fcd}}{8\mu_{\max}}$, we know that $\epsilon_{F}\leq\frac{\eta_{1}\kappa_{fcd}}{8\mu_{k}}$, which deduces the last inequality.

∎

Proof of Lemma 3.3

Proof.

Because $\mu_{k}\leq\frac{1}{\kappa_{bhm}}$, the $Cauchy\,decrease$ condition yields

[TABLE]

Assume that model $m_{k}$ are true, which means, for all $y\in B(x_{k},\delta_{k})$, we have

[TABLE]

And the estimates $f_{k}^{0}$ and $f_{k}^{d}$ are tight with $\epsilon_{F}\leq\kappa_{ef}$, we have

[TABLE]

The ratio $\rho_{k}$ can be rewritten as

[TABLE]

Applying the inequlaities (6.7), (6.8) and (6.9) to the above equality, and then using $\delta_{k}=\mu_{k}\left\|g_{k}\right\|$, we can obtain

[TABLE]

Since $\mu_{k}\leq\frac{\kappa_{fcd}(1-\eta_{1})}{8\kappa_{ef}}$, we have $\frac{8\kappa_{ef}}{\kappa_{fcd}}\mu_{k}\leq 1-\eta_{1}$. Thus, we conclude that $\rho_{k}\geq\eta_{1}$, which means that the $k$-th iteration is successful.

∎

Proof of Lemma 3.4

Proof.

First, we show that for all $k<T_{\epsilon}$, the following inequality holds

[TABLE]

If $\Lambda_{k}>\hat{\Lambda}$, we have $\Lambda_{k+1}\geq\gamma\hat{\Lambda}$ by the update process of the sequence $\Lambda_{k}$. Hence, $\Lambda_{k+1}\geq\hat{\Lambda}$. Now we assume that $\Lambda_{k}\leq\hat{\Lambda}$, by the definition of $\hat{\Lambda}$, we have

[TABLE]

If $I_{k}=1$ and $J_{k}=1$, i.e. model $M_{k}$ and estimates are all sufficiently accurate, from Lemma 3.3, we know that the iteration $k$ is successful. Thus $X_{k+1}=X_{k}+D_{k}$ and $\Lambda_{k+1}=\gamma\Lambda_{k}$. If $I_{k}J_{k}=0$, whether the iteration $k$ is successful or failed, we all have $\Lambda_{k+1}\geq\frac{1}{\gamma}\Lambda_{k}$.

From the above analysis, we conclude that (6.11) holds. Moreover we have observed that $p=P(W_{k}=1|\mathcal{F}_{k-1}^{M,F})=P(I_{k}J_{k}=1|\mathcal{F}_{k-1}^{M,F})\geq\alpha\beta>\frac{1}{2}$, which implies that Assumption 2(ii) holds.

∎

Proof of Theorem 3.1

Proof.

First, we recall the definition of $\Phi_{k}$, that is

[TABLE]

In the following proof, we consider three separate cases: (i) model $M_{k}$ is true and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight ($I_{k}=J_{k}=1$); (ii) model $M_{k}$ is false and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight ($J_{k}=1,I_{k}=0$); (iii) estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are loose ($J_{k}=0$). For each of these cases, we will analyze two possible outcomes of the iteration process, that is the iteration $k$ is successful or failed. Based on the above classifications, we rewrite the expected decrease of $\Phi_{k}$ as

[TABLE]

Before presenting the formal proof, we brief describe the key ideas. By choosing a suitable constant $\nu$, we can derive an upper bound on the expected decrease of $\Phi_{k}$ for each of the cases. When model $M_{k}$ is true and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight, no matter whether the iteration $k$ is successful or not, it will give rise to the decrease of $\Phi_{k}$ which is in proportion to $\Lambda_{k}\left\|\nabla f(X_{k})\right\|^{2}$. For the other two cases, due to the model error or inaccurate estimates, $\Phi_{k}$ may increase. However, the increment of $\Phi_{k}$ is still in proportion to $\Lambda_{k}\left\|\nabla f(X_{k})\right\|^{2}$. Therefore, by choosing sufficiently large $\alpha$ and $\beta$, the expectation of $\Phi_{k}$ can be guaranteed to decrease.

• (i).

Model $M_{k}$ is true and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight ($I_{k}=J_{k}=1$).

• a.

Iteration $k$ is successful. In this case, we have $X_{k+1}=X_{k}+D_{k}$, and $\Lambda_{k+1}=\gamma\Lambda_{k}.$

Because model $M_{k}$ is true. From Lemma 3.1, we know that if

[TABLE]

we have

[TABLE]

Besides, we can easily find the relation between $\left\|G_{k}\right\|$ and $\left\|\nabla f(X_{k})\right\|$, i.e.

[TABLE]

Using the above inequality, the triangle inequality and the fact that $\Lambda_{k}\leq\mu_{\max}$, we obtain

[TABLE]

Since $f$ is $L$-smooth and iteration $k$ is successful, we have that

[TABLE]

Using the fact that $(a+b)^{2}\leq 2(a^{2}+b^{2})$ and by a simple calculation, we get

[TABLE]

Particularly, the following holds with $\Lambda_{k+1}=\gamma\Lambda_{k}$ and $\,\Lambda_{k}\leq\mu_{\max}$

[TABLE]

Applying (6.15) and (6.19), we get

[TABLE]

We choose $\nu\in(0,1)$ to satisfy

[TABLE]

which implies that

[TABLE]

Consequently, (6.20) can be written as

[TABLE]

Applying (6.16) to (6.22), we have

[TABLE]

Furthermore, we assume that $\nu$ satisfies

[TABLE]

which yields

[TABLE]

(6.23) and (6.24) show that

[TABLE]

• b.

Iteration k is failed. In this case, we know $X_{k+1}=X_{k}$ and $\Lambda_{k+1}=\frac{1}{\gamma}\Lambda_{k}$, which deduce that

[TABLE]

We can choose a suitable $\nu\in(0,1)$ such that

[TABLE]

Then,

[TABLE]

Taking conditional expectation on $\mathcal{F}_{k-1}^{M,F}$, we have

[TABLE]

• (ii).

Model $M_{k}$ is false and estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight ($J_{k}=1,I_{k}=0$).

• a.

Iteration $k$ is successful. Because the iteration is successful, we have $X_{k+1}=X_{k}+D_{k}$, and $\Lambda_{k+1}=\gamma\Lambda_{k}.$ In this case, the estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are tight. Applying lemma 3.2, we have

[TABLE]

with $\epsilon_{F}\leq\frac{\eta_{1}\kappa_{fcd}}{8\mu_{\max}}$ and $\mu_{\max}\leq\frac{1}{\kappa_{bhm}}$. Due to the fact that $k$ is successful and $f$ is $L$-smooth, so (6.19) holds. Combining (6.19) and (6.30), we have

[TABLE]

Then we choose a suitable $\nu\in(0,1)$ such that $\frac{1-\nu}{\nu}\leq\frac{\eta_{1}\kappa_{fcd}}{8\gamma\mu_{\max}^{2}}$, which implies that

[TABLE]

Thus, it follows that

[TABLE]

• b.

Iteration $k$ is failed. Here, we have $X_{k+1}=X_{k}$ and $\Lambda_{k+1}=\frac{1}{\gamma}\Lambda_{k}$. In this case, (6.26) holds.

No matter whether the iteration $k$ is successful or failed, we always have

[TABLE]

Taking conditional expectation on the above inequality, we obtain

[TABLE]

• (iii).

Estimates $\left\{F_{k}^{0},F_{k}^{d}\right\}$ are loose ($J_{k}=0$).

• a.

Iteration $k$ is successful. In this case, we have $X_{k+1}=X_{k}+D_{k}$, and $\Lambda_{k+1}=\gamma\Lambda_{k}$. Because iteration $k$ is accepted and the trial step $D_{k}$ satisfies Assumption 3, we have

[TABLE]

where the last inequality reuses the fact that $\Delta_{k}=\Lambda_{k}\left\|G_{k}\right\|$ and $\mu_{\max}\leq\frac{1}{\kappa_{bhm}}$. Applying (6.36), a successful iteration deduces the following bound for the increment of $f$

[TABLE]

Then, the upper bound for $\Phi_{k+1}-\Phi_{k}$ can be obtained as follows

[TABLE]

where the first inequality uses (6.19), which is still true in this setting. We can choose $\nu\in(0,1)$ to satisfy $\frac{1-\nu}{\nu}\leq\frac{\eta_{1}\kappa_{fcd}}{4\gamma\mu_{\max}^{2}}$ so that

[TABLE]

Then (6.38) can be rewritten as

[TABLE]

• b.

Iteration $k$ is failed. In this case, we have $X_{k+1}=X_{k}$ and $\Lambda_{k+1}=\frac{1}{\gamma}\Lambda_{k}$. Then (6.26) holds.

Compared to (6.26), we see that (6.40) dominates the upper bound of $\Phi_{k+1}-\Phi_{k}$. Taking conditional expectations on inequality (6.40) and applying Assumption 1(iii), we have

[TABLE]

From relations $\kappa_{f}\leq 2\eta_{1}\kappa_{ef}$ and $\Lambda_{k}\leq\mu_{\max}\leq\min\left\{\frac{1}{\kappa_{bhm}},\frac{\kappa_{fcd}}{8\kappa_{ef}}\right\}$, it follows that

[TABLE]

Thus

[TABLE]

Now combining (6.29), (6.35) and (6.43), we can show that

[TABLE]

By a simple caculation, we have $P(\mathbb{I}_{\left\{I_{k}J_{k}=1\right\}}|\mathcal{F}_{k}^{M,F})\geq\alpha\beta$. The sixth line of above inequality uses the fact that event $\left\{(1-I_{k})J_{k}=1\right\}$ and $\left\{(1-J_{k})=1\right\}$ are disjoint, which implies that

[TABLE]

Thus, we have $P(\left\{(1-I_{k})J_{k}=1\right\}\bigcup\left\{(1-J_{k})=1\right\}|\mathcal{F}_{k-1}^{M,F})=P(1-I_{k}J_{k}=1|\mathcal{F}_{k-1}^{M,F})\leq 1-\alpha\beta$. Choosing suitable $\alpha$ and $\beta$ such that

[TABLE]

which implies that

[TABLE]

then we get the last inequality of (6.44).

Finally, the proof is completed.

∎

Proof of Theorem 3.3

Proof.

In this case, $f$ satisfies Assumption 6. Let us consider the stochastic process $\left\{\Lambda_{k},\Psi_{k}\right\}$ with

[TABLE]

The convexity of $f$ implies that

[TABLE]

Let $x=x^{\ast}$, $y=X_{k}$, it follows from the above inequality that

[TABLE]

Because $f$ is $L$-smooth, we have $\left\|\nabla f(X_{k})-\nabla f(x^{\ast})\right\|\leq L\left\|X_{k}-x^{\ast}\right\|$. Due to Assumption 6, we know the level set $\mathcal{L}$ is bounded, then

[TABLE]

Combining (6.49) and (6.50), we have

[TABLE]

From the above inequality and the result of Theorem 3.1, we have

[TABLE]

(6.52) implies that $\mathbb{E}[\Phi_{k+1}|\mathcal{F}_{k}^{M,F}]\leq\Phi_{k}$. Recalling the definition of $\Psi_{k}$, for all $k<T_{\epsilon}$, we have

[TABLE]

The first inequality of (6.53) follows from Jensen’s inequality which will be given in Lemma 6.3. The second inequality uses (6.52). The last inequality is due to the fact that $\mathbb{E}[\Phi_{k+1}|\mathcal{F}_{k}^{M,F}]\leq\Phi_{k}$. Here, we define an non-decreasing function $h(\cdot)$ as follows

[TABLE]

where $C_{1}=\frac{\alpha\beta(1-\nu)(1-\frac{1}{\gamma})}{2(\nu L+(1-\nu)\mu_{\max})^{2}D^{2}}$. Then we know that Assumption 2(iii) holds.

From Lemma 3.4, we can easily obtain that if $\alpha\beta>\frac{1}{2}$ and $\hat{\Lambda}$ is defined as (3.6), Assumption 2(ii) satisfies. Then we have Assumption 2 holds. Thus, the conclusion of Theorem 2.1 is true in this case.

Finally, substituting the expression of $\Psi_{0}$, $\hat{\Lambda}$ and $h$ into Theorem 2.1, we have

[TABLE]

where $M=\frac{2(\nu L+(1-\nu)\mu_{\max})^{2}D^{2}}{\alpha\beta\nu(1-\nu)(1-\frac{1}{\gamma})\hat{\Lambda}}.$ Here, we simplify the constant term as $\mathcal{O}(1)$.

Now, this completes the proof of Theorem 3.3. ∎

Proof of Theorem 3.4

Proof.

In this setting, $f$ is strongly convex with $\sigma>0$. We will consider the measure $\Psi_{k}$ as follows

[TABLE]

to analyze the theoretical complexity. Due to the strongly convexity, we have

[TABLE]

Then,

[TABLE]

It follows from (6.57) and Theorem 3.1 that

[TABLE]

The above inquality implies

[TABLE]

Recalling the definition $\Psi_{k}=\log(\Phi_{k})+\log(\frac{1}{\nu\epsilon})$, we have

[TABLE]

We can define

[TABLE]

where $C_{2}=\frac{\alpha\beta(1-\nu)(1-\frac{1}{\gamma})}{(\frac{\nu}{2\sigma}+(1-\nu)\frac{1}{L^{2}}\mu_{\max})}$.

From Lemma 3.4, we can easily see that Assumption 2(ii) satisfies if $\alpha\beta>\frac{1}{2}$ and $\hat{\Lambda}$ is defined as (3.6). Thus Assumption 2 holds. So the conclusion of Theorem 2.1 is true in strongly convex setting.

By substituting the expression of $\Psi_{0}$, $\hat{\Lambda}$ and $h$ into Theorem 2.1, we have

[TABLE]

where $M=\frac{(\frac{\nu}{2\sigma}+(1-\nu)\frac{1}{L^{2}}\mu_{\max})}{\alpha\beta(1-\nu)(1-\frac{1}{\gamma})\hat{\Lambda}}.$

Now the proof is finished.

∎

B: Related lemmas and algorithms for second-order STRME

Lemma 6.1** (Chebyshev Inequalitydurrett2019probability ).**

If $X$ is a random variable with mean $\mathbb{E}[X]$ and variance $Var(X)$, then

[TABLE]

Based on the exercise 4.1.2 in durrett2019probability , we prove the Chebyshev Inequality with conditional expectation.

Lemma 6.2**.**

If $X$ is a random variable given the $\sigma$-field $\mathcal{F}$, then

[TABLE]

Proof.

Let $A=\left\{X\,|\,|X-\mathbb{E}[X]|\geq v\right\}$, for $\forall v>0$, then

[TABLE]

Thus, the proof is finished. ∎

Lemma 6.3** (Jensen InequalityProbability_Bertsekas ).**

Assume that $f$ is continuous and convex. If X is a random variable, then

[TABLE]

Remark 6.1**.**

If $f$ is concave in Lemma 6.3, then we get the opposite result, i.e. $\mathbb{E}[f(X)]\leq f(\mathbb{E}[X])$.