Stochastic Optimization Using a Trust-Region Method and Random Models

Introduction

Derivative free optimization (DFO) has recently grown as a field of nonlinear optimization which addressed optimization of black-box functions, that is functions whose value can be (approximately) computed by some numerical procedure or an experiment, while their closed-form expressions and/or derivatives are not available and cannot be approximated accurately or efficiently. Although the role of derivative-free optimization is particularly important when objective functions are noisy, traditional DFO methods have been developed primarily for deterministic functions. The fields of stochastic optimization and stochastic approximation on the other hand focus on optimizing functions that are stochastic in nature. Much of the focus of these methods depend on the availability and use of stochastic derivatives, however, some work has addressed stochastic black box functions, typically by some sort of a finite differencing scheme .

In this paper, using methods developed for DFO, we aim to solve

where the noise $\varepsilon$ is a random variable.

In recent years, some DFO methods have been extended to and analyzed for stochastic functions . Additionally, stochastic approximation methodologies started to incorporate techniques from the DFO literature . The analysis in all that work assumes some particular structure of the noise, including the assumption that the noisy function values give an unbiased estimator of the true function value.

There are two main classes of methods in this setting of stochastic optimization: stochastic gradient (SG) methods (such as the well known Robbins-Monro method) and sample averaging (SA) methods. The former (SG) methods roughly work as follows: they obtain a realization of an unbiased estimator of the gradient at each iteration and take a step in the direction of the negative gradient. The step sizes progressively diminish and the iterates are averaged to form a sequence that converges to a solution. These methods typically have very inexpensive iterations, but exhibit slow convergence and are strongly dependent on the choice of algorithmic parameters, particularly the sequence of step sizes. Many variants exist that average the gradient information from past iterations and are able to accept sufficiently small, but nondecreasing step sizes . However, the convergence remains slow, and parameter tuning remains necessary in most cases. Moreover, the majority of these methods have been developed exclusively for convex functions and may not converge in non convex settings. An exception is the randomized stochastic (accelerated) gradient (RS(A)G) method presented in , . This method employs random stopping criteria to provide theoretical first-order convergence even in nonconvex cases, but does not appear to be very practical.

The second class of methods, (SA), is based on sample averaging of the function and gradient estimators,which is applied to reduce the variance of the noise. These methods repeatedly sample the function value at a set of points in hopes to ensure sufficient accuracy of the function and gradient estimates. For a thorough introduction and references therein, see . The optimization method and sampling process are usually tightly connected in these approaches, hence, again, algorithmic parameters need to be specially chosen and tuned. These methods tend to be more robust with respect to parameters and enjoy faster convergence at a cost of more expensive iterations. However, none of these methods are applicable in the case of biased noise and they suffer significantly in the presence of outliers.

The goal of this paper is to show that a standard efficient unconstrained optimization method, such as a trust region method, can be applied, with very small modifications, to stochastic nonlinear (not necessarily convex) functions and can be guaranteed to converge to first order stationary points as long as certain conditions are satisfied. We present a general framework, where we do not specify any particular sampling technique. The framework is based on the trust region DFO framework , and its extension to probabilistic models . In terms of this framework and the certain conditions that must be satisfied, we essentially assume that

the local models of the objective function constructed on each iteration satisfy some first order accuracy requirement with sufficiently high probability,

and that function estimates at the current iterate and at a potential next iterate are sufficiently accurate with sufficiently high probability.

The main novelty of this work is the analysis of the framework and the resulting weaker, more general, conditions compared to prior work. In particular,

we do not assume that the probabilities of obtaining sufficiently accurate models and estimates are increasing (they simply need to be above a certain constant) and

we do not assume any distribution of the random models and estimates. In other words, if a model or estimate is inaccurate, it can be arbitrarily inaccurate, i.e. the noise in the function values can have nonconstant bias.

It is also important to note that while our framework and model requirements are borrowed from prior work on DFO, this framework applies to derivative-based optimization as well. Later in the paper we will discuss different settings which will fit into the proposed framework.

This paper consists of two main parts. In the first part we propose and analyze a trust region framework, which utilizes random models of $f(x)$ at each iteration to compute the next potential iterate. It also relies on (random, noisy) estimates of the function values at the current iterate and the potential iterate to gauge the progress that is being made. The convergence analysis then relies on requirements that these models and these estimates are sufficiently accurate with sufficiently high probability. Beyond these conditions, no assumptions are made about how these models and estimates are generated. The resulting method is a stochastic process that is analyzed with the help of martingale theory. The method is shown to converge to first order stationary points with probability one.

There is a very large volume of work on SA and SG, most of which is quite different from our proposed analysis and method. However, we will mention a few works here that are most closely related to this paper and highlight the differences. The three methods existing in the literature we will compare with are by Deng and Ferris , SPSA (simultaneous perturbations stochastic approximation) , and SCSR (sampling controlled stochastic recursion) . These three settings and methods are most closely related to our work because they all rely on using models of the objective function that can both incorporate second-order information and whose accuracy with respect to a “true” model can be dynamically adjusted. In particular, Deng and Ferris apply the trust-region model-based derivative free optimization method UOBYQA in a setting of sample path optimization . In and , the author applies an approximate gradient descent and Newton method, respectively, with gradient and Hessian estimates computed from specially designed finite differencing techniques, with decaying finite differencing parameter. In a very general scheme is presented, where various deterministic optimization algorithms are generalized as stochastic counterparts, with the stochastic component arising from the stochasticity of the models and the resulting step of the optimization algorithm. We now compare some key components of these three methods with those of our framework, which we hereforth refer to as STORM (STochastic Optimization with Random Models).

Deng and Ferris: The assumptions of the sample path setting are roughly as follows: on each iteration $k$, given a collection of points $X^{k}=\{x_{1}^{k},\ldots,x_{p}^{k}\}$ one can compute noisy function values $f(x_{1}^{k},\varepsilon^{k}),\ldots,f(x_{p}^{k},\varepsilon^{k})$. The noisy function values are assumed to be realizations of an unbiased estimator of true values $f(x_{1}^{k}),\ldots,f(x_{p}^{k})$. Then, using multiple, say $N_{k}$, realizations of $\varepsilon^{k}$, average function values $f^{N_{k}}(x_{1}^{k}),\ldots,f^{N_{k}}(x_{p}^{k})$ can be computed. A quadratic model $m^{N_{k}}_{k}(x)$ is then fit into these function values, and so a sequence of models $\{m^{N_{k}}_{k}(x)\}$ is created using a nondecreasing sequence of sampling rates $\{N_{k}\}$. The assumption on this sequence of models is that each of them satisfies a sufficient decrease condition (with respect to the true model of the true function $f$) with probability $1-\alpha_{k}$, such that $\sum_{k=1}^{\infty}\alpha_{k}<\infty$. The trust region maintenance follows the usual scheme like that in UOBYQA, hence the steps taken by the algorithm can be increased or decreased depending on the observed improvement of the function estimates.

SPSA: The first order version of this method assumes that $f(x,\varepsilon)$ is an unbiased estimate of $f(x)$, and the second order version, 2SPSA, assumes that an unbiased estimate of $\nabla f(x)$, $g(x,\varepsilon)\in{\bf R}^{n}$, can be computed. Gradient (in the first order case) and Hessian (in the second order case) estimates are constructed using an interesting randomized finite differencing scheme. The finite difference step is assumed to be decaying to zero. An approximate steepest descent direction or approximate Newton direction are then constructed and a step of length $t_{k}$ is taken along this direction. The sequence $\{t_{k}\}$ is assumed to be decaying in the usual Robbins-Monro way, that is $t_{k}\to 0$, $\sum_{k}t_{k}=\infty$. Hence, while no increase in accuracy of the models is assumed (they only need to be accurate in expectation), the step size parameter and the finite differencing parameter need to be tuned. Decaying step sizes often lead to slow convergence, as has been observed often in stochastic optimization literature.

SCSR: This is a very general scheme which can include multiple optimization methods and sampling rates. The key ingredients of this scheme are a deterministic optimization method, and a stochastic variant that approximates it. The stochastic step (recursion) is assumed to be a sufficiently accurate approximation of the deterministic step with increasing probability (the probabilities of failure for each iteration are summable). This assumption is stronger than the one in this paper. In addition, another key assumption made for SCSR is that the iterates produced by the base deterministic algorithm converge to the unique optimal minimizer $x^{*}$. Not only we do not assume here that the minimizer/stationary point is unique, but we also do not assume a priori that the iterates form a convergent sequence, since they may not do so in a non convex setting, while every iterate subsequence converges to a stationary point.

STORM: Like the Deng and Ferris method, we utilize a trust-region, model-based framework, where the size of the trust region can be increased or decreased according to empirically observed function decrease and the size of the observed approximate gradients. The desired accuracy of the models is tied only to the trust region radius in our case, while for Deng and Ferris, it is tied to both the radius and the size of true model gradients (the second condition is harder to ensure). In either method, this desired accuracy is assumed to hold with some probability - in STORM, this probability remains constant throughout the progress of the algorithm, while for Deng and Ferris it has to converge to $1$ sufficiently rapidly.

There are three major advantages to our results. First of all, in the case of unbiased noise, the sampling rate is directly connected to the desired accuracy of the estimates and the probability with which this accuracy is achieved. Hence, for the STORM method, the sampling rate may increase or decrease according to the trust region radius, eventually increasing only when necessary, i.e. when the noise become dominating. For all the other methods listed here, the sampling rate is assumed to increase monotonically. Secondly, in the case of biased noise, we can still prove convergence of our method, as long as the desired accuracy is achieved with a fixed probability. In other words, we allow for the noise to be arbitrarily large with a small, but fixed probability, on each iteration. This allows us to consider new models of noise which cannot be handled by any of the other methods discussed here. In addition, STORM incorporates first and second order models without changing the algorithm - the step size parameter (i.e., the trust region radius) and other parameters of the method are chosen almost identically to the standard practices of the trust region methods, which have proved to be very effective in practice for unconstrained nonlinear optimization. In Section 6 we show that the STORM method is very effective in different noise settings and is very robust with respect to sampling strategies.

Finally, we want to point to an unpublished work , which proposes a very similar method to the one in this paper. Both methods were developed based on the trust region DFO method with random models for deterministic functions analyzed in and extended to the stochastic setting. Some of the assumptions in this paper were inspired by an early version of . However, the assumptions and the analysis in are quite different from the ones in this paper. In particular, they rely on the assumption that $f(x,\varepsilon)$ is an unbiased estimate of $f(x)$, hence their analysis does not extend to the biased case. Also they assume that the probability of having an accurate model at the $k$-th iteration is at least $1-\alpha_{k}$, such that $\alpha_{k}\to 0$, while for our method this probability can remain bounded away from zero. Similarly, they assume that the probability of having accurate function estimates at the $k$-th iteration also converges to $1$ sufficiently rapidly, while in our case it is again constant. Their analysis, as a result, is very different from ours, and does not generalize to various stochastic settings (they only focus on the derivative free setting with additive noise). The advantage of their method is that they do not need to put a restriction on acceptable step sizes, when the norm of the gradient of the model is small. We, on the other hand, impose such restriction in our method and use it in the proof of our main result. However, as we discuss later in the paper this restriction can be relaxed at the cost of more complex algorithm and analysis. In practice, we do not implement this restriction, hence our basic implementation is virtually identical to that in except that we implement a variety of model building strategies, while only one such strategy (regression models based on randomly rotated orthogonal samples sets) is implemented in . Thus we do not directly compare the empirical performance of our method with the method in since we view them as more or less the same method.

We conclude this section by introducing some frequently used notations and their meanings. The rest of the paper is organized as follows. In Section 2 we introduce the trust region framework, followed by Section 3, where we discuss the requirements on our random models and function estimates. The main convergence results are presented in Section 4. In Section 5 we discuss various noise scenarios and how sufficiently accurate models and estimates can be constructed in these cases. Finally, we present computational experiments based on these various noise scenarios in Section 6.

Let $\|\cdot\|$ denote the Euclidean norm and $B(x,\Delta)$ denote the ball of radius $\Delta$ around $x$, i.e., $B(x,\Delta):\{y:\|x-y\|\leq\Delta\}$. $\Omega$ denotes the probability sample space, according to the context, and a sample from that space is denoted by $\omega\in\Omega$. As a rule, when we describe a random process within the algorithmic framework, uppercase letters, e.g. the $k$-th iterate $X_{k}$, will denote random variables, while lowercase letters will denote realizations of the random variable, e.g. $x_{k}=X_{k}(\omega)$ is the $k$-th iterate for a particular realization of our algorithm.

Trust Region Method

We consider the trust-region class of methods for minimization of unconstrained functions. They operate as follows: at each iteration $k$, given the current iterate $x_{k}$ and a trust-region radius $\delta_{k}$, a (random) model $m_{k}(x)$ is built, which serves as an approximation of $f(x)$ in $B(x_{k},\delta_{k})$. The model is assumed to be of the form

It is possible to generalize our framework to other forms of models, as long as all conditions on the models, described below, hold. We consider quadratic models for simplicity of the presentation and because they are the most common. The model $m_{k}(x)$ is minimized (approximately) in $B(x_{k},\delta_{k})$ to produce a step $s_{k}$ and (random) estimates of $f(x_{k})$ and $f(x_{k}+s_{k})$ are obtained, denoted by $f_{k}^{0}$ and $f_{k}^{s}$ respectively. The achieved reduction is measured by comparing $f_{k}^{0}$ and $f_{k}^{s}$ and if reduction is deemed sufficient, then $x_{k}+s_{k}$ is chosen as the next iterate $x_{k+1}$. Otherwise the iterate remains $x_{k}$. The trust-region radius $\delta_{k+1}$ is then chosen by either increasing or decreasing $\delta_{k}$ according to the outcome of the iteration. The details of the algorithm are presented in Algorithm 1.

The trial step computed on each iteration has to provide sufficient decrease of the model; in other words it has to satisfy the following standard fraction of Cauchy decrease condition:

For every $k$, the step $s_{k}$ is computed so that

for some constant $\kappa_{fcd}\in(0,1].$

If progress is achieved and a new iterate is accepted in the $k$-th iteration then we call this a successful iteration. Otherwise, the iteration is unsuccessful (and no step is taken). Hence a successful iteration occurs when $\rho_{k}\geq\eta_{1}$ and $\|g_{k}\|\geq\eta_{2}\delta_{k}$. However, a successful iteration does not necessarily yield an actual reduction in the true function $f$. This is because the values of $f(x)$ are not accessible in our stochastic setting and the step acceptance decision is made merely based on the estimates of $f(x_{k})$ and $f(x_{k}+s_{k})$. If these estimates, $f_{k}^{0}$ and $f_{k}^{s}$, are not accurate enough, a successful iteration can result in an increase of the true function value. Hence we consider two types of successful iterations - those where $f(x)$ is in fact decreased proportionally to $f_{k}^{0}-f_{k}^{s}$, which we call true successful iterations, and all other successful iterations, where the decrease of $f(x)$ can be arbitrarily small or even negative, which we call false successful iterations. Our setting and algorithmic framework does not allow us to determine which successful iterations are true and which ones are false, however, we will be able to show that true successful iterations occur sufficiently often for convergence to hold, if the random estimates $f_{k}^{0}$ and $f_{k}^{s}$ are sufficiently accurate.

A trust region framework based on random models was introduced and analyzed in . In that paper, the authors introduced the concept of probabilistically fully-linear models to determine the conditions that random models should satisfy for convergence of the algorithm to hold. However, the randomness in the models in their setting arises from the the construction process, and not from the noisy objective function. It is assumed in that the function values at the current iterate and the trial point can be computed exactly and hence all successful iterations are true in that case. In our case, it is necessary to define a measure for the accuracy of the estimates $f_{k}^{0}$ and $f_{k}^{s}$ (which, as we will see, generally has to be tighter than the measure of accuracy of the model). We will use a modified version of the probabilistic estimates introduced in .

Probabilistic Models and Estimates

To formalize conditioning on the past, let $\mathcal{F}_{k-1}^{M\cdot F}$ denote the $\sigma$-algebra generated by $M_{0},\cdots,M_{k-1}$ and $F_{0},\cdots,F_{k-1}$ and let $\mathcal{F}_{k-{1}/{2}}^{M\cdot F}$ denote the $\sigma$-algebra generated by $M_{0},\cdots,M_{k}$ and $F_{0},\cdots,F_{k-1}$.

To formalize sufficient accuracy, let us recall a measure for the accuracy of deterministic models introduced in and (with the exact notation introduced in ).

Suppose $\nabla f$ is Lipschitz continuous. A function $m_{k}$ is a $\kappa$-fully linear model of $f$ on $B(x_{k},\delta_{k})$ provided, for $\kappa=(\kappa_{ef},\kappa_{eg})$ and $\forall y\in B$,

In this paper we extend the following concept of probabilistically fully-linear models which is proposed in .

A sequence of random models $\{M_{k}\}$ is said to be $\alpha$-probabilistically $\kappa$-fully linear with respect to the corresponding sequence $\{B(X_{k},\Delta_{k})\}$ if the events

where $\mathcal{F}^{M}_{k-1}$ is the $\sigma$-algebra generated by $M_{0},\cdots,M_{k-1}$.

These probabilistically fully-linear models have the very simple properties that they are fully-linear (i.e., accurate enough) with sufficiently high probability, conditioned on the past, and they can be arbitrarily inaccurate otherwise. This property is somewhat different from the properties of models typical to stochastic optimization (such as, for example, stochastic gradient based models), where assumptions on the expected value and the variance of the models is imposed. We will discuss this in more detail in Section 5.

In this paper, aside from sufficiently accurate models, we require estimates of the function values $f(x_{k})$, $f(x_{k}+s_{k})$ that are sufficiently accurate. This is needed in order to evaluate whether a step is successful, unlike the case in where the exact values $f(x_{k})$ and $f(x_{k}+s_{k})$ are assumed to be available. The following definition of accurate estimates is a modified version of that used in .

The estimates $f_{k}^{0}$ and $f_{k}^{s}$ are said to be $\epsilon_{F}$-accurate estimates of $f(x_{k})$ and $f(x_{k}+s_{k})$, respectively, for a given $\delta_{k}$ if

We now modify Definitions 3.2 and 3.3 and introduce definitions of probabilistically accurate models and estimates which we will use throughout the remainder of the paper.

A sequence of random models $\{M_{k}\}$ is said to be $\alpha$-probabilistically $\kappa$-fully linear with respect to the corresponding sequence $\{B(X_{k},\Delta_{k})\}$ if the events

where $\mathcal{F}^{M\cdot F}_{k-1}$ is the $\sigma$-algebra generated by $M_{0},\cdots,M_{k-1}$ and $F_{0},\cdots,F_{k-1}$.

A sequence of random estimates $\{F_{k}^{0},F_{k}^{s}\}$ is said to be $\beta$-probabilistically $\epsilon_{F}$-accurate with respect to the corresponding sequence $\{X_{k},\Delta_{k},S_{k}\}$ if the events

where $\epsilon_{F}$ is a fixed constant and $\mathcal{F}_{k-{1}/{2}}^{M\cdot F}$ is the $\sigma$-algebra generated by $M_{0},\cdots,M_{k}$ and $F_{0},\cdots,F_{k-1}$.

Using Definitions 3.4 and 3.5 we will require in our analysis that our method has access to $\alpha$-probabilistically $\kappa$-fully linear models, for some fixed $\kappa=(\kappa_{ef},\kappa_{eg})$ and to $\beta$-probabilistically $\epsilon_{F}$-accurate estimates, for some fixed, sufficiently small $\epsilon_{F}$. Thus, the model and the estimate accuracy will be assumed to be proportional to $\delta_{k}^{2}$ (with some probability), the condition on the estimates is somewhat tighter because of the upper bound on $\epsilon_{F}$. However, we will see that this upper bound is not very small.

Procedures for obtaining probabilistically fully-linear models and probabilistically accurate estimates under different models of noise are discussed in Section 5.

Convergence Analysis

We now present first-order convergence analysis for the general framework described in Algorithm 1. Towards that end, we assume that the function $f$ and its gradient are Lipschitz continuous in regions considered by the algorithm realizations.

The second assumption provides a uniform upper bound on the model Hessian.

There exists a positive constant $\kappa_{bhm}$ such that, for every $k$, the Hessian $H_{k}$ of all realizations $m_{k}$ of $M_{k}$ satisfy

Note that since we are concerned with convergence to a first order stationary point in this paper, the bound $\kappa_{bhm}$ can be chosen to be any nonnegative number, including zero. Allowing a larger bound will give more flexibility to the algorithm and may allow better Hessian approximations, but as we will see in the convergence analysis, this imposes restrictions on the trust region radius and some other algorithmic parameters.

We now state the following result from martingale literature (see Exercise 5.3.1) that will be useful later in our analysis.

Let $G_{k}$ be a submartingale, i.e., a sequence of random variables which, for every $k$,

where $\mathcal{F}_{k-1}^{G}=\sigma(G_{0},\ldots,G_{k-1})$ is the $\sigma$-algebra generated by $G_{0},\ldots,G_{k-1}$, and $E[G_{k}|\mathcal{F}_{k-1}^{G}]$ denotes the conditional expectation of $G_{k}$ given the past history of events $\mathcal{F}_{k-1}^{G}$.

Assume further that $G_{k}-G_{k-1}\leq M<\infty$, for every $k$. Then,

We now prove some auxiliary lemmas that provide conditions under which decrease of the true objective function $f(x)$ is guaranteed. The first lemma states that if the trust region radius is small enough relative to the size of the model gradient and if the model is fully linear, then the step $s_{k}$ provides a decrease in $f(x)$ proportional to the size of the model gradient. Note that the trial step may still be rejected if the estimates $f_{k}^{0}$ and $f_{k}^{s}$ are not accurate enough.

Suppose that a model $m_{k}$ of the form (2) is a $(\kappa_{ef},\kappa_{eg})$-fully linear model of $f$ on $B(x_{k},\delta_{k})$. If

then the trial step $s_{k}$ leads to an improvement in $f(x_{k}+s_{k})$ such that

Using the Cauchy decrease condition, the upper bound on model Hessian and the fact that $\|g_{k}\|\geq\kappa_{bhm}\delta_{k}$, we have

Since the model is $\kappa$-fully linear, one can express the improvement in $f$ achieved by $s_{k}$ as

where the last inequality is implied by $\delta_{k}\leq\frac{\kappa_{fcd}}{8\kappa_{ef}}\|g_{k}\|$. ∎

The next lemma shows that for $\delta_{k}$ small enough relative to the size of the true gradient $\nabla f(x_{k})$, the guaranteed decrease in the objective function, provided by $s_{k}$, is proportional to the size of the true gradient.

Under Assumption 4.3, suppose that a model is $(\kappa_{ef},\kappa_{eg})$-fully linear on $B(x_{k},\delta_{k})$. If

then the trial step $s_{k}$ leads to an improvement in $f(x_{k}+s_{k})$ such that

where $C_{1}=\frac{\kappa_{fcd}}{4}\cdot\max\left\{\frac{\kappa_{bhm}}{\kappa_{bhm}+\kappa_{eg}},\frac{8\kappa_{ef}}{8\kappa_{ef}+\kappa_{fcd}\kappa_{eg}}\right\}.$

The definition of a $\kappa$-fully-linear model yields that

Since condition (11) implies that $\|\nabla f(x_{k})\|\geq\max\left\{\kappa_{bhm}+\kappa_{eg},\frac{8\kappa_{ef}}{\kappa_{fcd}}+\kappa_{eg}\right\}\delta_{k}$, we have

Hence, the conditions of Lemma 4.5 hold and we have

Since $\|g_{k}\|\geq\|\nabla f(x)\|-\kappa_{eg}\delta_{k}$ in which $\delta_{k}$ satisfies (11), we also have

We now prove the lemma that states that, if a) the estimates are sufficiently accurate, b) the model is fully-linear and c) the trust-region radius is sufficiently small relatively to the size of the model gradient, then a successful step is guaranteed.

Under Assumption 4.3, suppose that $m_{k}$ is $(\kappa_{ef},\kappa_{eg})$-fully linear on $B(x_{k},\delta_{k})$ and the estimates $\{f_{k}^{0},f_{k}^{s}\}$ are $\epsilon_{F}$-accurate with $\epsilon_{F}\leq\kappa_{ef}$. If

Since $\delta_{k}\leq\frac{\|g_{k}\|}{\kappa_{bhm}}$, the Cauchy decrease condition and the uniform bound on $H_{k}$ immediately yield that

The model $m_{k}$ being $(\kappa_{ef},\kappa_{eg})$-fully linear implies that

Since the estimates are ${\epsilon}_{F}$-accurate with $\epsilon_{F}\leq\kappa_{ef}$, we obtain

where we have used the assumption $\delta_{k}\leq\frac{\kappa_{fcd}(1-\eta_{1})}{8\kappa_{ef}}\|g_{k}\|$ to deduce the last inequality. Hence, $\rho_{k}\geq\eta_{1}$. Moreover, since $\|g_{k}\|\geq\eta_{2}\delta_{k}$, the $k$-th iteration is successful. ∎

Finally, we state and prove the lemma which guarantees an amount of decrease of the objective function on a true successful iteration.

Under Assumption 4.3, suppose that the estimates $\{f_{k}^{0},f_{k}^{s}\}$ are ${\epsilon}_{F}$-accurate with ${\epsilon}_{F}<\frac{1}{4}\eta_{1}\eta_{2}\kappa_{fcd}\min\left\{\frac{\eta_{2}}{\kappa_{bhm}},1\right\}$. If a trial step $s_{k}$ is accepted (a successful iteration occurs), then the improvement in $f$ is bounded below as follows

where $C_{2}=\frac{1}{2}\eta_{1}\eta_{2}\kappa_{fcd}\min\left\{\frac{\eta_{2}}{\kappa_{bhm}},1\right\}-2\epsilon_{F}>0$.

An iteration being successful indicates that $\|g_{k}\|\geq\eta_{2}\delta_{k}$ and $\rho\geq\eta_{1}$. Thus,

Then, since the estimates are $\epsilon_{F}$-accurate, we have that the improvement in $f$ can be bounded as

where $C_{2}=\frac{1}{2}\eta_{1}\eta_{2}\kappa_{fcd}\min\left\{\frac{\eta_{2}}{\kappa_{bhm}},1\right\}-2\epsilon_{F}>0$. ∎

To prove convergence of Algorithm 1 we will need to assume that models $\{M_{k}\}$ and estimates $\{F_{k}^{0},F_{k}^{s}\}$ are sufficiently accurate with sufficiently high probability.

Given values of $\alpha,\beta\in(0,1)$ and $\epsilon_{F}>0$, there exist $\kappa_{eg}$ and $\kappa_{ef}$ such that the the sequence of models $\{M_{k}\}$ and estimates $\{F_{k}^{0},F_{k}^{s}\}$ generated by Algorithm 1 are, respectively, $\alpha$-probabilistically $(\kappa_{ef},\kappa_{eg})$- fully-linear and $\beta$-probabilistically $\epsilon_{F}$-accurate.

Note that this assumption is a statement about the existence of constants $\kappa=(\kappa_{ef},\kappa_{eg})$ given an $\alpha$, $\beta$ and $\epsilon_{F}$ - we will determine exact conditions on $\alpha$, $\beta$ and $\epsilon_{F}$ in Theorem 4.11 and Lemma 4.12 below.

The following theorem states that the trust-region radius converges to zero with probability $1$.

Let Assumptions 4.1 and 4.3 be satisfied and assume that in Algorithm 1 the following holds.

The step acceptance parameter $\eta_{2}$ is chosen so that

The accuracy parameter of the estimates satisfies

Then $\alpha$ and $\beta$ can be chosen so that, if Assumption 4.9 holds for these values, then the sequence of trust-region radii, $\{\Delta_{k}\}$, generated by Algorithm 1 satisfies

We base our proof on properties of the random function $\Phi_{k}=\nu f(X_{k})+(1-\nu)\Delta_{k}^{2}$, where $\nu\in(0,1)$ is a fixed constant, which is specified below. A similar function is used in the analysis in , but the analysis itself is different. The overall goal is to show that there exists a constant $\sigma>0$ such that for all $k$

Since $f$ is bounded from below and $\Delta_{k}>0$, we have that $\Phi_{k}$ is bounded from below for all $k$ and hence if (24) holds on every iteration, then by summing (24) over $k\in(1,\infty)$ and taking expectations on both sides we can conclude that (23) holds with probability $1$. Hence, to prove the theorem we need to show that (24) holds on each iteration.

Let us pick some constant $\zeta$ which satisfies

We now consider two possible cases: $\|\nabla f(x_{k})\|\geq\zeta\delta_{k}$ and $\|\nabla f(x_{k})\|<\zeta\delta_{k}$. We will show that (24) holds in both cases and hence it holds on every iteration. Given $\zeta$ we now select $\nu\in(0,1)$ such that

As usual, let $x_{k}$, $\delta_{k}$, $s_{k}$, $g_{k}$, and $\phi_{k}$ denote realizations of random quantities $X_{k}$, $\Delta_{k}$, $S_{k}$, $G_{k}$, and $\Phi_{k}$, respectively.

Let us consider some realization of Algorithm 1. Note that on all successful iterations, $x_{k+1}=x_{k}+s_{k}$ and $\delta_{k+1}=\min\{\gamma\delta_{k},\delta_{max}\}$ with $\gamma>1$, hence

On all unsuccessful iterations, $x_{k+1}=x_{k}$ and $\delta_{k+1}=\frac{1}{\gamma}\delta_{k}$, i.e.

For each iteration and each of the two cases we consider, we will analyze the four possible combined outcomes of the events $I_{k}$ and $J_{k}$ as defined in (7) and (8), respectively.

Before presenting the formal proof let us outline the key ideas. We will show that, unless both the model and the estimates are bad on iteration $k$, we select $\nu\in(0,1)$ sufficiently close to $1$, so that the decrease in $\phi_{k}$ on a successful iteration is greater than the decrease on an unsuccessful iteration (which is equal to $b_{1}$, according to (28)). When the model and the estimates are both bad, an increase in $\phi_{k}$ may occur. This increase is bounded by a value proportional to $\delta_{k}^{2}$ when $\|\nabla f(x_{k})\|<\zeta\delta_{k}$. When $\|\nabla f(x_{k})\|\geq\zeta\delta_{k}$, though, the increase in $\phi_{k}$ may be proportional to $\|\nabla f(x_{k})\|\delta_{k}$. However, since iterations with good models and good estimates will provide decrease in $\phi_{k}$ also proportional to $\|\nabla f(x_{k})\|\delta_{k}$, by choosing values of $\alpha$ and $\beta$ close enough to 1, we can ensure that in expectation $\phi_{k}$ decreases.

$I_{k}$ and $J_{k}$ are both true, i.e., both the model and the estimates are good on iteration $k$. From the definition of $\zeta$, we know

Then since the model $m_{k}$ is $\kappa$-fully linear and, from $\eta_{2}>\kappa_{bhm}$, ${\epsilon}_{F}\leq\kappa_{ef}$ and $0<\eta_{1}<1$, it is easy to show that the condition (11) in Lemma 4.6 holds. Therefore, the trial step $s_{k}$ leads to a decrease in $f$ as in (12).

and the estimates $\{f_{k}^{0},f_{k}^{s}\}$ are ${\epsilon}_{F}$-accurate, with ${\epsilon}_{F}\leq\kappa_{ef}$, the condition (15) in Lemma 4.7 holds. Hence, iteration $k$ is successful, i.e. $x_{k+1}=x_{k}+s_{k}$ and $\delta_{k+1}=\gamma\delta_{k}$.

with $C_{1}$ defined in Lemma 4.6. Since $\|\nabla f(x_{k})\|\geq\zeta\delta_{k}$ we have

$I_{k}$ is true and $J_{k}$ is false, i.e., we have a good model and bad estimates on iteration $k$.

In this case, Lemma 4.6 still holds, that is $s_{k}$ yields a sufficient decrease in $f$, hence, if the iteration is successful, we obtain (29) and (30). However, the step can be erroneously rejected, because of inaccurate probabilistic estimates, in which case we have an unsuccessful iteration and (28) holds. Since (26) holds the right hand side of the first relation in (30) is strictly smaller than the right hand side of the first relation in (28) and therefore, (28) holds whether the iteration is successful or not.

$I_{k}$ is false and $J_{k}$ is true, i.e., we have a bad model and good estimates on iteration $k$. In this case, iteration $k$ can be either successful or unsuccessful. In the unsuccessful case (28) holds. When the iteration is successful, since the estimates are ${\epsilon}_{F}$-accurate and (22) holds then by Lemma 4.8 (20) holds with $C_{2}\geq\frac{1}{4}\eta_{1}\eta_{2}\kappa_{fcd}$. Hence, in this case we have

Again, due to the choice of $\nu$ satisfying (26) we have that, as in case (b), (28) holds whether the iteration is successful or not.

$I_{k}$ and $J_{k}$ are both false, i.e., both the model and the estimates are bad on iteration $k$.

Inaccurate estimates can cause the algorithm to accept a bad step, which may lead to an increase both in $f$ and in $\delta_{k}$. Hence in this case $\phi_{k+1}-\phi_{k}$ may be positive. However, combining the Taylor expansion of $f(x_{k})$ at $x_{k}+s_{k}$ and the Lipschitz continuity of $\nabla f(x)$ we can bound the amount of increase in $f$, hence bounding $\phi_{k+1}-\phi_{k}$ from above. By adjusting the probability of outcome (d) to be sufficiently small, we can ensure that in expectation $\Phi_{k}$ is sufficiently reduced.

In particular, from Taylor’s Theorem and Lipschitz continuity of $\nabla f(x)$ we have, respectively,

From this we can derive that any increase of $f(x_{k})$ is bounded by

where $C_{3}=1+\frac{3L}{2\zeta}$. Hence, the change in function $\phi$ is bounded:

Now we are ready to take the expectation of $\Phi_{k+1}-\Phi_{k}$ for the case when $\|\nabla f(X_{k})\|\geq\zeta\Delta_{k}$. We know that case (a) occurs with a probability at least $\alpha\beta$ (conditioned on the past) and in that case $\phi_{k+1}-\phi_{k}=b_{2}<0$ with $b_{2}$ defined in (29). Case (d) occurs with probability at most $(1-\alpha)(1-\beta)$ and that case $\phi_{k+1}-\phi_{k}$ is bounded from above by $b_{3}>0$. Cases (b) and (c) occur otherwise and in those cases $\phi_{k+1}-\phi_{k}$ is bounded from above by $b_{1}<0$, with $b_{1}$ defined in (28). Finally we note that $b_{1}>b_{2}$ due to our choice of $\nu$.

Hence, we can combine (28), (29), (31), and (32), and use $B_{1}$, $B_{2}$, and $B_{3}$ as random counterparts of $b_{1}$, $b_{2}$, and $b_{3}$, to obtain the following bound

where the last inequality holds because $\alpha\beta-\frac{1}{\gamma^{2}}(\alpha(1-\beta)+(1-\alpha)\beta)+(1-\alpha)(1-\beta)\leq[\alpha+(1-\alpha)][(\beta+(1-\beta)]=1.$

Let us choose $0<\alpha\leq 1$ and $0<\beta\leq 1$ so that they satisfy

where the last inequality is the result of (26). It is important to note that the quantity $\frac{1}{2}$ in the numerator of (33) is chosen so that the first inequality of the above expression holds. This quantity can be made arbitrarily small, and with appropriate adjustment to (26) one can have

where $\theta_{1}$ is positive and arbitrarily close to zero and $\theta_{2}>1$ and arbitrarily close to one. We will return to this comment in a later remark, but for the sake of simplicity, we choose $\theta_{1}=\frac{1}{2}$ and $\theta_{2}=2$.

Recall that $\|\nabla f(X_{k})\|\geq\zeta\Delta_{k}$, hence

For the purposes of this lemma and the liminf-type convergence result, which will follow, bound (35) is sufficient. We will use bound (34) in the proof of the $\lim$-type convergence result.

First we note that if $\|g_{k}\|<\eta_{2}\delta_{k}$, then we have an unsuccessful step and (28) holds. Hence, we now assume that $\|g_{k}\|\geq\eta_{2}\delta_{k}$ and again consider four possible outcomes. We will show that in all situations, except when both the model and the estimates are bad, (28) holds. In the remaining case, because $\|\nabla f(x_{k})\|<\zeta\delta_{k}$ the increase in $\phi_{k}$ can be bounded from above by a multiple of $\delta_{k}^{2}$. Hence by selecting appropriate values for probabilities $\alpha$ and $\beta$ we will be able to establish the bound on expected decrease in $\Phi_{k}$ as in Case 1.

$I_{k}$ and $J_{k}$ are both true, i.e., both the model and the estimates are good on iteration $k$.

The iteration may or may not be successful, even though $I_{k}$ is true. On successful iteration good model ensures reduction in $f$. Applying the same argument as in the case 1(c) we establish (28).

$I_{k}$ is true and $J_{k}$ is false, i.e., we have a good model and bad estimates on iteration $k$.

On unsuccessful iterations, (28) holds. On successful iterations, $\|g_{k}\|\geq\eta_{2}\delta_{k}$ and $\eta_{2}\geq\kappa_{bhm}$ imply that

Since $I_{k}$ is true, the model is $\kappa$-fully-linear, and the function decrease can be bounded as

It follows that, if $k$-th iterate is successful, then

Again by choosing $\nu\in(0,1)$ so that (26) holds, we ensure that right hand side of (36) is strictly smaller than that of (28), hence (28) holds, whether the iteration is successful or not.

Remark: $\eta_{2}$ may need to be a relatively large constant to satisfy (21). This is due to the fact that the model has to be sufficiently accurate to ensure decrease in the function if a step is taken, since the step is accepted based on poor estimates. Note that $\eta_{2}$ restricts the size of $\Delta_{k}$, which is used both as a bound on the step size and the control of the accuracy. In general it is possible to have two separate quantities (related by a constant) - one to control the step size and another to control the accuracy. Hence, it is possible to modify our algorithm to accept steps larger than $\|g_{k}\|/\eta_{2}$. This will make the algorithm more practical, but the analysis much more complex. In this paper, we choose to stay with the simplest version, but keep in mind that the condition (26) is not terminally restrictive.

$I_{k}$ is false and $J_{k}$ is true, i.e., we have a bad model and good estimates on iteration $k$.

This case is analyzed identically to the case 1(c).

$I_{k}$ and $J_{k}$ are both false, i.e., both the model and the estimates are bad on iteration $k$.

Here we bound the maximum possible increase in $\phi_{k}$ using the Taylor expansion and the Lipschitz continuity of $\nabla f(x)$.

We are now ready to bound the expectation of $\phi_{k+1}-\phi_{k}$ as we did in Case 1, except that in Case 2 we simply combine (37), which holds with probability at most $(1-\alpha)(1-\beta)$ and (28) which holds otherwise.

if we choose probabilities $0<\alpha\leq 1$ and $0<\beta\leq 1$ so that the following holds,

In conclusion, combining (35) and (39), and noting that $1-\frac{1}{\gamma^{2}}<\gamma^{2}-1$ we have

which implies, that (24) holds with $\sigma=-\frac{1}{2}(1-\nu)(1-\frac{1}{\gamma^{2}}-1)<0$. This concludes the proof of the theorem.

To summarize the conditions on the probabilities involved in Theorem 4.11 to ensure that the theorem holds, we state the following additional lemma.

Let all assumptions of Theorem 4.11 hold. The statement of Theorem 4.11 holds if the $\alpha$ and $\beta$ are chosen to satisfy the following conditions:

with $C_{1}=\frac{\kappa_{fcd}}{4}\cdot\max\left\{\frac{\kappa_{bhm}}{\kappa_{bhm}+\kappa_{eg}},\frac{8\kappa_{ef}}{8\kappa_{ef}+\kappa_{fcd}\kappa_{eg}}\right\}$ and $\zeta=\kappa_{eg}+\eta_{2}.$

The proof follows simply from combining expression for $C_{3}$ and condition (26) with (33) and (38).

Clearly, choosing $\alpha$ and $\beta$ sufficiently close to $1$ will satisfy this condition.

We will briefly illustrate through a simple example how these algorithmic parameters scale with problem data.

Recall that $L$ is the Lipschitz constant of the gradient of $f$ and of $f$ over ${\cal L}_{enl}(x^{0})$. It is reasonable to expect that $\kappa_{ef}$ and $\kappa_{eg}$ are quantities that scale with $L$, since Taylor models satisfy this condition, as do polynomial interpolation and regression models based on well-poised data sets . Let us assume for the sake of an example that $\kappa_{ef}=\kappa_{eg}=10L$. The bound on model Hessians $\kappa_{bhm}$ can be chosen to be arbitrarily small, at the expense of limiting the class of models, however, it is clearly reasonable to choose it as something that scales with $L$, if this information is available. Let us assume that $\kappa_{bhm}=10L$, as well. In a standard trust region method, a common choice of algorithmic parameters would use $\kappa_{fcd}=\frac{1}{2}$, $\gamma=2$, and $\eta_{1}=\frac{1}{2}$.

The reader can easily verify that with these parameter choices and previous assumptions, Lemma 4.12 states that we must choose $\eta_{2}\geq 32L$. The intermediate constants satisfy $\zeta\geq 42L$ and $C_{1}=\frac{2}{17}$. Without loss of generality, we will simply accept $\zeta=42L$.

From observing that, given the above values of the constants,

We note that as $\eta_{1}$ and $\kappa_{fcd}$ constants are driven closer to $1$, the constant $440$ can be reduced by up to a factor of $4$.

Supposing that $\kappa_{ef},\kappa_{eg}$, and $\kappa_{bhm}$ scale linearly with $L$, then $\eta_{2}$, $\epsilon_{F}$, and the expressions relating to $\alpha$ and $\beta$ in Corollary 4.12 are all functions in $L$ satisfying

where we use the notation $\Theta(\cdot)$ to indicate the $O(\cdot)$ relationship with moderate constants.

Recall our remark made earlier in Theorem 4.11, case 2b, on how $\eta_{2}$ bounds our step sizes. Indeed, if $\eta_{2}\geq\Theta(L)$ has to be imposed this may force the algorithm to take small step sizes throughout. However, as mentioned earlier, the analysis of Theorem 4.11 can be modified by introducing a tradeoff between the size of $\eta_{2}$ and the accuracy parameters $\epsilon_{F}$ and $\kappa_{ef}$ (as both of these constant parameters can be made smaller). It also may be advantageous to choose $\eta_{2}$ dynamically. Exploring this is a subject for future work. In the practical implementations that we will discuss in Section 6, we do not make use of the algorithmic parameter $\eta_{2}$ at all, and so even though $\eta_{2}$ is effectively arbitrarily close to , the algorithm still works.

Note that if $\beta=0$, then $\Delta_{k}\to 0$ for any positive value of $\alpha$, which is the case shown in , since, as we pointed earlier, condition (33) can be restated with $\frac{1}{2}$ replaced by an arbitrary small positive value.

1 The liminf-type convergence

We are ready to prove a liminf-type first-order convergence result, i.e., that a subsequence of the iterates drive the gradient of the objective function to zero. The proof follows closely that in , the key difference being the assumption on the function estimates that are needed to ensure that a good step gets accepted by Algorithm 1.

Let the assumptions of Theorem 4.11 and Lemma 4.12 hold. Suppose additionally that $\alpha\beta\geq\frac{1}{2}$. Then the sequence of random iterates generated by Algorithm 1, $\{X_{k}\}$, almost surely satisfies

We prove this result by contradiction conditioned on the almost sure event $\Delta_{k}\rightarrow 0$. Let us thus assume that there exists $\epsilon^{\prime}$ such that, with positive probability, we have

Let $\{x_{k}\}$ and $\{\delta_{k}\}$ be realizations of $\{X_{k}\}$ and $\{\Delta_{k}\}$, respectively for which $\|\nabla f(x_{k})\|\geq\epsilon^{\prime},\quad\forall k.$ Since $\lim\limits_{k\to\infty}\delta_{k}=0$ (because we conditioned on $\Delta_{k}\rightarrow 0$), there exists $k_{0}$ such that for all $k\geq k_{0}$,

We define a random variable $R_{k}$ with realizations $r_{k}=\log_{\gamma}\left(\dfrac{\delta_{k}}{b}\right)$. Then for the realization $\{r_{k}\}$ of $\{R_{k}\}$, $r_{k}<0$ for $k\geq k_{0}$. The main idea of the proof is to show that such realizations occur only with probability zero, hence obtaining a contradiction with the initial assumption of $\|\nabla f(x_{k})\|\geq\epsilon^{\prime}$ $\forall k$.

We first show that $R_{k}$ is a submartingale. Recall the events $I_{k}$ and $J_{k}$ in Definitions 3.4 and 3.5. Consider some iterate $k\geq k_{0}$ for which $I_{k}$ and $J_{k}$ both occur, which happens with probability $P(I_{k}\cap J_{k})\geq\alpha\beta$. Since (48) holds we have exactly the same situation as in Case 1(a) in the proof of Theorem 4.11. In other words, we can apply Lemmas 4.6 and 4.7 to conclude that the $k$-th iteration is successful, hence, the trust-region radius is increased. In particular, since $\delta_{k}\leq\frac{\delta_{\max}}{\gamma}$, $\delta_{k+1}=\gamma\delta_{k}$. Consequently, $r_{k+1}=r_{k}+1$.

Let $\mathcal{F}_{k-1}^{I\cdot J}=\sigma({I_{0}},\cdots,{I_{k-1}})\cap\sigma({J_{0}},\cdots,J_{k-1})$. For all other outcomes of $I_{k}$ and $J_{k}$, which occur with total probability of at most $1-\alpha\beta$, we have $\delta_{k+1}\geq\gamma^{-1}\delta_{k}$. Hence

as long as $\alpha\beta\geq 1/2$, which implies that $R_{k}$ is a submartingale.

Now let us construct another submartingale $W_{k}$, on the same probability space as $R_{k}$ which will serve as a lower bound on $R_{k}$ and for which $\left\{\underset{k\to\infty}{\lim\sup}\ W_{k}=\infty\right\}$ holds almost surely. Define indicator random variables $\textbf{1}_{I_{k}}$ and $\textbf{1}_{J_{k}}$ such that $\textbf{1}_{I_{k}}=1$ if $I_{k}$ occurs, $\textbf{1}_{I_{k}}=0$ otherwise, and similarly, $\textbf{1}_{J_{k}}=1$ if $J_{k}$ occurs, $\textbf{1}_{J_{k}}=0$ otherwise. Then define

Notice that $W_{k}$ is a submartingale since

where the last inequality holds because $\alpha\beta\geq 1/2$. Since $W_{k}$ only has $\pm 1$ increments, it has no finite limit. Therefore, by Theorem 4.4, we have $\left\{\underset{k\to\infty}{\lim\sup}\ W_{k}=\infty\right\}$.

By the construction of $R_{k}$ and $W_{k}$, we know that $r_{k}-r_{k_{0}}\geq w_{k}-w_{k_{0}}$. Therefore, $R_{k}$ has to be positive infinitely often with probability one. This implies that the sequence of realizations $r_{k}$ such that $r_{k}<0$ for $k\geq k_{0}$ occurs with probability zero. Therefore our assumption that $\|\nabla f(X_{k})\|\geq\epsilon^{\prime}$ hold for all $k$ with positive probability is false and

2 The lim-type convergence

In this subsection we show that $\lim_{k\to\infty}\|\nabla f(X_{k})\|=0$ almost surely.

We now state an auxiliary lemma, which is similar to the one in , but requires a different proof because in our case the function values $f(X_{k})$ can increase with $k$, while in the case considered in , function values are monotonically nonincreasing.

Let the same assumptions that were made in Theorem 4.16 hold. Let $\{X_{k}\}$ and $\{\Delta_{k}\}$ be sequences of random iterates and random trust-region radii generated by Algorithm 1. Fix $\epsilon>0$ and define the sequence $\left\{K_{\epsilon}\right\}$ consisting of the natural numbers $k$ for which $\|\nabla f(X_{k})\|>\epsilon$ (note that $K_{\epsilon}$ is a sequence of random variables). Then,

From Theorem 4.11 we know that $\sum\Delta_{k}^{2}<\infty$ and hence $\Delta_{k}\to 0$ almost surely. For each realization of Algorithm 1 and a sequence $\{\delta_{k}\}$, there exists $k_{0}$ such that $\delta_{k}\leq\epsilon/\zeta$, $\forall k\geq k_{0}$, where $\zeta$ is defined as in Theorem 4.11. Let $K_{0}$ be the random variable with realization $k_{0}$ and let $K$ denote the sequence of indices $k$ such that $k\in K_{\epsilon}$ and $k\geq K_{0}$. Then for all $k\in K$, Case 1 of Theorem 4.11 holds, i.e., $\|\nabla f(X_{k})\|\geq\zeta\Delta_{k}$, since $\|\nabla f(X_{k})\|\geq\epsilon$ for all $k\in K$. From this and from (34) we have

Recall that $\Phi_{k}$ is bounded from below. Hence, summing up the above inequality for all $k\in K$ and taking the expectation, we have that

almost surely. Since $K_{\epsilon}\subseteq K\cup\{k\leq K_{0}\}$ and $K_{0}$ is finite almost surely then the statement of the lemma holds. ∎

We are now ready to state the $\lim$-type result.

Let the same assumptions as in Theorem 4.16 hold. Let $\{X_{k}\}$ be a sequence of random iterates generated by Algorithm 1. Then, almost surely,

The proof of this result, is almost identical to the proof of the same theorem in hence we will not present the proof here. The key idea of the proof is to show that if the theorem does not hold, then with positive probability

with $K_{\epsilon}$ defined as in Lemma 4.17. This result is shown using Lipschitz continuity of the gradient and does not depend on the stochastic nature of the algorithm. Since this result contradicts the almost sure result of Lemma 4.17, we can conclude that the statement of the theorem holds almost surely.

Constructing models and estimates in different stochastic settings.

We now discuss various settings of stochastic noise in the objective function and how $\alpha$-probabilistically $\kappa$-fully linear models and $\beta$-probabilistically $\epsilon_{F}$-accurate estimates can be obtained in these settings.

where $\omega$ is a random variable which induces the noise.

The example of noise which is typically considered in stochastic optimization is the “i.i.d.” noise, that is noise with distribution independent of the function values. Here we consider a somewhat more general setting, where the noise is unbiased for all $f$, i.e.,

This is the typical noise assumption in stochastic optimization literature. In the case of unbiased noise as above, constructing estimates and models that satisfy our assumptions is fairly straight-forward. First, let us consider the case when only the noisy function values are available (without any gradient information), where we want to construct a model that is $\kappa$-fully linear in a given trust region $B(x^{0},\delta)$ with some reasonably large probability, $\alpha$.

One can employ standard sample averaging approximation techniques to reduce the variance of the function evaluations. In particular, let $\bar{f}_{p}(x,\omega)=\frac{1}{p}\sum_{i=1}^{p}{f}(x,\omega_{i})$, where $\omega_{i}$ are the i.i.d. realizations of the noise $\omega$. Then, by Chebyshev inequality, for any $v>0$,

In particular, we want $v=\kappa^{\prime}_{ef}\delta^{2}$ for some $\kappa^{\prime}_{ef}>0$ and $\frac{V}{pv^{2}}\leq 1-\alpha^{\prime}$ for some $\alpha^{\prime}$, which can be ensured by choosing $p\geq\frac{V}{(\kappa^{\prime}_{ef})^{2}(1-\alpha^{\prime})\delta^{4}}$.

We now construct a fully linear model as follows: given a well-poised setSee for details on well-poised sets and how they can be obtained. $Y$ of $n+1$ points in $B(x^{0},\delta)$, at each point $y^{i}\in Y$, we compute $\bar{f}_{p}(y^{i},\omega)$ and build a linear interpolation model $m(x)$ such that $m(y^{i})=\bar{f}_{p}(y^{i},\omega)$, for all $i=1,\ldots,n+1$. Hence, for any $y^{i}\in Y$, we have

Moreover, the events $\{|m(y^{i})-f(y^{i})]|>\kappa^{\prime}_{ef}\delta^{2}\}$ are independent, hence

It is easy to show using, for example, techniques described in , that $m(x)$ is a $\kappa$-fully linear model of $E_{\omega}[{f}(x,\omega)]$ in $B(x^{0},\delta)$ for appropriately chosen $\kappa=(\kappa_{eg},\kappa_{ef})$, with probability at least $\alpha=(\alpha^{\prime})^{n+1}$.

Computing the $\beta$-probabilistically $\epsilon_{F}$-accurate estimates of $f(x,\omega)$ can be done analogously to the construction of the models described above.

The majority of stochastic optimization and sample average approximation methods focus on derivative based optimization where it is assumed that, in addition to $f(x,\omega)$, $\nabla_{x}f(x,\omega)$ is also available, and that the noise in the gradient computation is also independent of $x$, that is

(in general the variance of the gradient and the function value are not the same, but here for simplicity we bound both by $V$).

In the case when the noisy gradient values are available, the construction of fully linear models in $B(x^{0},\delta)$ is simpler. Let $\bar{\nabla}{f}_{p}(x,\omega)=\frac{1}{p}\sum_{i=1}^{p}{\nabla f}(x,\omega_{i})$. Again, by extension of Chebychev inequality, for $p$ such that

Hence the linear expansion $m(x)=\bar{f}_{p}(x^{0},\omega)+\bar{\nabla}f_{p}(x^{0},\omega)^{T}(x-x^{0})$ is a $\kappa$-fully linear model of $f(x)=E_{\omega}[{f}(x,\omega)]$ on $B(x^{0},\delta)$ for appropriately chosen $\kappa=(\kappa_{eg},\kappa_{ef})$, with probability at least $\alpha=(\alpha^{\prime})^{2}$.

In it is shown that least squares regression models based on sufficiently large strongly poised sample sets are $\alpha$-probabilistically $\kappa$-fully linear models.

There are many existing methods and convergence results using sample average approximations and stochastic gradients for stochastic optimization with i.i.d. or unbiased noise. Some of these methods have been shown to achieve optimal sampling rate , that is they converge to the optimal solution while sampling the gradient at the best possible rate. We do not provide convergence rates in this paper (it is a subject for future research), hence it remains to be seen if our algorithm can achieve the optimal rate. Our contribution here is the method which applies beyond the i.i.d. case, as we will discuss below. In Section 6, however, we demonstrate that our method can have superior numerical behavior compared to standard sample averaging even in the case of the i.i.d. noise, so it is at least competitive in practice.

Function computation failures.

Another example is solving a system of nonlinear black-box equations. Assume that we seek $x$ such that $\sum_{i}(f_{i}(x))^{2}=0$, for some functions $f_{i}(x)$, $i=1,\ldots,m$ that are computed by numerical simulation, with noise. As is often done in practice (and is supported by our theory) the noise in the function computation is reduced as the algorithm progresses, for example, by reducing the size of a discretization, step size, or convergence tolerance within the black-box computation. These adjustments for noise reduction usually increase the workload of the simulation. With the increase of the workload, there is an increased probability of failure of the code. Hence, the smaller the values of $f_{i}(x)$, the more likely the computation of $f_{i}(x)$ will fail and some inaccurate value is returned.

where $\sigma(x)$ is the probability with which the function $f(x)$ is computed inaccurately, and $\omega(x)$ is some random function of $x$, for which only an upper bound $V$ is known. This case is idealized, because we assume that with probability $1-\sigma(x)$, $f(x)$ is computed exactly. It is trivial to extend this example to the case when $f(x)$ is computed with an error, but this error can be made sufficiently small.

For this model of function computation failures we have

and it is clear, that for any $\sigma(x)>0$, unless $E[\omega(x)]\equiv\text{some constant}$, optimizing $E_{\omega}[{f}(x,\omega)]$ does not give the same result as optimizing $f(x)$. Hence applying Monte-Carlo sampling within an optimization algorithm solving this problem is not a correct approach.

We now observe that constructing $\alpha$-probabilistically $\kappa$-fully linear models and $\beta$-probabilistically $\epsilon_{F}$-accurate estimates is trivial in this case, assuming that $\sigma(x)\leq\sigma$ for all $x$, when $\sigma$ is small enough. In particular, given a trust region $B(x^{0},\delta)$, sampling a function $f(x)$ on a sample set $Y\subset B(x^{0},\delta)$ well-poised for linear interpolation will produce a $\kappa$-fully linear model in $B(x^{0},\delta)$ with probability at least $(1-\sigma)^{|Y|}$, since with this probability all of the function values are computed exactly. Similarly, for any $s\in B(x^{0},\delta)$, the function estimates $F^{0}$ and $F^{s}$ are both correct with probability at least $(1-\sigma)^{2}$. Assuming that $(1-\sigma)^{|Y|}\geq\alpha$ and $(1-\sigma)^{2}\geq\beta$, where $\alpha$ and $\beta$ satisfy the assumptions of Theorem 4.11 and Lemma 4.12 and $\alpha\beta\geq\frac{1}{2}$ as in Theorem 4.16, we observe that the resulting models satisfy our theory.