Serge Gratton
email: serge.gratton@enseeiht.fr

D Avid Titley-Peloquin

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible Stochastic Conditionin g of Matrix Functions

Keywords: sensitivity analysis, perturbation analysis, c onditioning, uncertainty propagation, matrix functions

We investigate the sensitivity of matrix fun ctions to random noise in their input. We propose the notion of a stochastic condition number, which determines, to first order, the sensitivity of a matrix function to random noise. We derive an upper bound on the stochastic condition number that can be estimated efficiently by using "small-sample" estimation techniques. The bound can be used to estimate the median, or any other quantile, of the error in a function's output when its input is subjected to random noise. We give numerical experiments illustrating the effectiveness of our stochastic error estimate.

1. Introduction. How sensitive are matrix functions to perturbations in their input? This is a fundamental question in numerical linear algebra. By "matrix function" we mean a general mapping F : Ω → R p×q , where Ω is an open subset of R m×n . As for the perturbations, these might represent uncertainties in the data or rounding errors arising from computations in finite precision arithmetic. The goal is to quantify the effect that such uncertainties might have on the computed function value.

Sensitivity analyses in numerical linear algebra are usually deterministic in nature. Worstcase error bounds are obtained that hold asymptotically in the small-perturbation limit. (For an excellent overview, we recommend [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF][START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF].) To the best of our knowledge, few attempts have been made in the literature to rigorously quantify the sensitivity of matrix functions to random noise. We believe that it is important to address this issue for the following reasons.

First, although rounding errors are known not to be random, some random models of rounding have been proposed, and, in fact, some automatic error analysis software does use random perturbations. (See, for example, the discussion in [10, sections 2.8, 26.5].) Therefore, it seems worthwhile to attempt to gain at least some insight into the effect of rounding errors on computed function values by using a stochastic analysis. Furthermore, in many applications, the uncertainty due to rounding errors is dominated by measurement noise in the data, which is explicitly modeled as random. In this case, it is a natural idea to take into account the random nature of the noise in the sensitivity analysis.

For example, there has been some recent interest in the data assimilation community in the sensitivity of Krylov method iterates to random noise [START_REF] Zhu | Observation sensitivity calculations using the adjoint of the gridpoint statistical interpolation (GSI) analysis system[END_REF][START_REF] Moore | Estimates of analysis and forecast errors derived from the adjoint of 4D-Var[END_REF]. In these applications, the Krylov iteration forms the inner loop of a nested outer-inner iterative scheme, and the uncertainty must be propagated through the inner loop. Specifically, let x (j) (b + h) denote the jth iterate of the method of conjugate gradients (CG) applied to solving Ax = b + h, where A ∈ R n×n is deterministic, symmetric, and positive definite, b ∈ R n is deterministic, and h is modeled as a normally distributed random vector, h ∼ N (0, Σ). Ideally, one would solve the linear system exactly. Of course, since the exact solution is a linear function of h, the distribution of x(b + h) = A -1 (b + h) is easily derived. However, in large-scale applications, one often has no choice but to terminate CG after very few iterations, long before convergence has occurred. In this setting, statistical information about x (j) (b + h), for a small fixed j, as opposed to x(b + h), is required. This is complicated by the fact that, in general, x (j) (b) is known to be very nonlinear in b. This example motivated us to study the problem in a more general context.

First, let us review the classical deterministic analysis. Let Ω represent an open subset of R m×n , F : Ω → R p×q be a matrix function defined everywhere on Ω, and A ∈ Ω be given data. We seek to estimate some measure of the difference F (A + H) -F (A), for instance, its norm or the magnitude of some of its elements, given some limited information about the perturbation H. For example, suppose we know that in some norm H is bounded by δ, and from this information we wish to bound F (A + H) -F (A) . In this case the quantity of interest would be (1.1) c δ = sup

H ≤δ F (A + H) -F (A) .
Unfortunately, it is usually unfeasible to compute c δ . Instead, the normwise absolute condition number of F at A, defined as

(1.2) κ = lim δ→0 c δ δ = lim δ→0 sup H ≤δ F (A + H) -F (A) δ ,
leads to the first-order estimate (FOE)

(1.3) F (A + H) -F (A) κ H .
Rice [START_REF] Rice | A theory of condition[END_REF] has shown that if F is Fréchet differentiable at A, then κ is the operator norm of the Fréchet derivative of F at A. There is a large body of literature dedicated to computing or bounding normwise and componentwise condition numbers of matrix functions. Despite its widespread use, sensitivity analysis using condition numbers has some drawbacks. First, taking the supremum over H ≤ δ as in (1.1) and (1.2) may lead to an unnecessarily pessimistic estimate of the typical sensitivity of F because of pathological values of H that are highly unlikely to occur in practice. Furthermore, the FOE (1.3) is only valid asymptotically for "sufficiently small" H . In practice it is usually hard to determine how small is sufficiently small, and noise present in real data is often not sufficiently small. This paper addresses the first above-mentioned drawback. Specifically, we model H as a random matrix and attempt to quantify some of the statistical properties of the random variable F (A + H) -F (A) F , where • F denotes the matrix Frobenius norm. One idea that has appeared in the literature [START_REF] Turing | Rounding-off errors in matrix processes[END_REF][START_REF] Fletcher | Expected conditioning[END_REF][START_REF] Stewart | Stochastic perturbation theory[END_REF] is estimating to first order the root-mean-squared (RMS) error

E F (A + H) -F (A) 2 F .
(Throughout, E{x} denotes the expected value of x and V{x} = E (x -E{x})(x -E{x}) T .) Unfortunately, the above RMS error is infinite for many common matrix functions and common distributions. We argue that much more insight into the sensitivity of F can be obtained by considering the quantiles of the random variable F (A + H) -F (A) F rather than its RMS value. In particular, a useful measure of central tendency is the 0.5-quantile, i.e., the median. As a first step, in analogy with (1.2), we only consider small perturbation asymptotics. We define a stochastic condition number κΣ in terms of the median, as follows. Let the elements of H be random variables following some distribution such that E{H} = 0 and V{vec(H)} = σ 2 Σ for a given positive semidefinite matrix Σ ∈ R mn×mn . Then

κΣ ≡ lim sup σ→0 Med F (A + H) -F (A) F σ .
We shall obtain an upper bound on the stochastic condition number that holds regardless of the distribution of the elements of H. With minor modifications the bound applies to other quantiles as well. The results are obtained by bounding the probability

Prob F (A + H) -F (A) F ≥ τ
for any given τ > 0, to first order in σ. Our hope is that this investigation will be the first step toward obtaining global bounds that hold for arbitrary σ > 0, at least for some classes of matrix functions. We shall make further comments about this point in the conclusion.

Our upper bound on the stochastic condition number can be estimated efficiently, for example, by using the small-sample estimates of Kenney et. al. [START_REF] Kenney | Small-sample statistical condition estimates for general matrix functions[END_REF][START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF]. The resulting error estimate is much more representative of the sensitivity of F to random perturbations than an FOE based on the deterministic condition number.

The rest of the paper is organized as follows. In section 2 we review some basic notions that are useful in later derivations. In section 3 we present some measures of sensitivity to random noise and show how these compare to the deterministic condition number (1.2) and FOE (1.3). Section 4 covers the efficient computation of our stochastic error estimate. Numerical experiments illustrating the theory are given in section 5, and we conclude with a short discussion in section 6.

2.

Preliminaries. First we establish our notation and review some basic notions that are used in later sections.

Notation. We use uppercase letters to denote matrices and lowercase Roman letters for vectors and indices. Subscripts are used to denote elements of a vector or matrix, while superscripts denote terms in a sequence. For example, A ij is the (i, j)th element of the matrix A, and h j is the jth element of the vector h, while u (i) ∈ R n is the ith term in a sequence of vectors. Otherwise, scalars are designated by lowercase Greek letters.

Fréchet derivative. Throughout we assume that F : Ω ⊆ R m×n → R p×q is Fréchet differentiable at A ∈ Ω. The Fréchet derivative of F at A is the unique bounded linear operator F ′ (A) defined by the relation (2.1)

F (A + H) = F (A) + F ′ (A)(H) + R(H), lim H→0 R(H) H = 0.
The matrix representation of F ′ (A) in the standard basis is the Jacobian matrix J A ∈ R pq×mn . In other words,

(2.2) vec(F ′ (A)(H)) = J A vec(H).
Multi-index notation. A multi-index is simply a vector of natural numbers,

v = [v 1 , v 2 , . . . , v n ] T ∈ N n 0 ,
where N 0 = {0, 1, 2, . . . }. The magnitude and factorial of an n-dimensional multi-index are defined as

|v| = v 1 + • • • + v n , v! = v 1 ! • • • v n !, while, for v ∈ N n 0 and h = [h 1 , . . . , h n ] T ∈ R n , h v = h v 1 1 h v 2 2 • • • h vn n .
Higher-order partial derivatives can be written compactly as

∂ v x = ∂ v 1 ∂x v 1 1 • • • ∂ vn ∂x vn n .
In this notation, the Taylor series of an analytic function g : R n → R can be expressed as

(2.3) g(x + h) = v∈N n 0 ∂ v x g(x) v! h v = g(x) + ∇g(x) T h + v∈N n 0 |v|≥2 ∂ v x g(x) v! h v .
Probability notation. We work with the standard probability space consisting of the sample space R, its σ-algebra of Borel sets, and Lebesgue measure, with straightforward extension to random vectors and matrices. We use Prob{X} to denote the probability of an event X. E{x} and V{x} stand for the expected value and covariance matrix of a random vector x. As shorthand for E{x} = u and V{x} = Σ we write x ∼ (u, Σ), and if additionally x follows a multivariate normal distribution, x ∼ N (u, Σ). A median of a random variable α is any scalar τ satisfying both

(2.4) Prob{α ≤ τ } ≥ 1 2 and Prob{α ≥ τ } ≥ 1 2 .
If α is a continuous random variable, then Prob{α = τ } = 0. In this case, (2.4)

is equivalent to Prob{α ≤ τ } = Prob{α ≥ τ } = 1 2 .
If additionally α has a strictly positive probability density function, then its median is unique.

In any case, when we write the median of α, or Med{α}, we are referring to Med{α} = sup τ : (2.4) holds .

More generally, for q ∈ (0, 1), the qth quantile of α is any scalar τ satisfying

(2.5) Prob{α ≤ τ } ≥ q, Prob{α ≥ τ } ≥ 1 -q.
When we write the qth quantile of α, we mean

Q q {α} = sup τ : (2.5) holds .
(The median is the 0.5-quantile.) In contrast to the expected value, the quantiles of a random variable are always finite. Basic probability inequalities. If α and β are random variables such that α ≤ β, then for any τ ∈ R,

(2.6) Prob{α ≥ τ } ≤ Prob{β ≥ τ } , Prob{β ≤ τ } ≤ Prob{α ≤ τ } .
For random variables α and β, for any τ, ǫ ∈ R,

(2.7) Prob{α + β ≥ τ } ≤ Prob{α ≥ τ (1 -ǫ)} + Prob{β ≥ τ ǫ} .
Markov's inequality. Given a nonnegative random variable α, then for any τ > 0,

(2.8) Prob{α ≥ τ } ≤ E{α} τ .
Quadratic forms. If M ∈ R m×n and x ∼ (0, Σ),

(2.9)

E{ M x 2 } ≤ E M x 2 2 = M Σ 1/2 F .
If H ∈ R m×n and vec(H) ∼ (0, Σ),

(2.10) E{ H F } ≤ E H 2 F = Σ 1/2 F .
3. Notions of stochastic conditioning.

3.1. "Expected conditioning," revisited. As mentioned in the introduction, there is not a large body of literature dealing with the sensitivity of matrix functions to random noise. We could trace the idea back to Turing [START_REF] Turing | Rounding-off errors in matrix processes[END_REF], who considers the function

F (A) = A -1 b for nonsingular A ∈ R n×n and b ∈ R n . It is straightforward to verify that (A + H) -1 b -A -1 b = -A -1 HA -1 b + o( H ).

Turing ignores the o( H ) term and gives an expression for the RMS value of

A -1 HA -1 b 2 .
Fletcher considers the same problem in [START_REF] Fletcher | Expected conditioning[END_REF]. (See also [10, page 136].) In Fletcher's approach, the elements H ij of H are independent random variables with mean 0 and variance σ 2 A 2 ij . Under these assumptions,

E A -1 HA -1 b 2 2 = σ [A -1 ][A][A -1 b] 1/2 1 ,
where [A] denotes the matrix whose entries are A 2 ij . Fletcher calls the above the expected condition number of F (A) = A -1 b. He also generalizes the above to other matrix functions. Once again, the idea is to take a first-order expansion,

F (A + H) -F (A) = F ′ (A)(H) + o( H ),
ignore the o( H ) term, and find an expression for

(3.1) E F ′ (A)(H) 2 F .
Stewart [START_REF] Stewart | Stochastic perturbation theory[END_REF] independently derives a similar result. In Stewart's case, the random perturbation has the form H = S c GS T r , where G is a random matrix whose elements are uncorrelated with mean 0 and variance 1, and S c and S r are deterministic scalings.

The quantity in (3.1) is easy to analyze because it involves a quadratic form in the elements of H. Suppose that vec(H) ∼ (0,

σ 2 Σ),
where σ ∈ R and Σ ∈ R mn×mn is symmetric nonnegative definite. Then from (2.2) and (2.9),

E F ′ (A)(H) 2 F = E vec(F ′ (A)(H)) 2 2 = E J A vec(H) 2 2 = σ 2 J A Σ 1/2 2 F .
Thus, assuming that one can indeed ignore the o( H ) term in the expansion (2.1), one obtains

E F (A + H) -F (A) 2 F ≈ E F ′ (A)(H) 2 F = σ J A Σ 1/2 F . (3.2)
This reduces to the previously mentioned results of Fletcher and Stewart when the covariance matrix Σ is chosen according to their respective approaches.

One might question whether there is a rigorous theoretical justification for dropping the o( H ) term in the above. For instance, because it was obtained from a first-order expansion in H, we might conclude that (3.2) is tight for perturbations with sufficiently small covariance matrix, i.e., for sufficiently small σ. In fact, as noted by Stewart [START_REF] Stewart | Stochastic perturbation theory[END_REF], this is not the case, and (3.2) is often meaningless. For example, with F (A) = A -1 b and vec(H) ∼ N (0, σ 2 I mn ), the left-hand side of (3.2) is infinite for any σ > 0, while the right-hand side tends to 0 as σ → 0. Stewart [START_REF] Stewart | Stochastic perturbation theory[END_REF]Theorem 2.8] gives some arguments as to why the lower-order terms can be ignored. Nevertheless, it is not immediately clear under which conditions the approximation (3.2) does indeed make sense.

In order to guarantee that (3.2) is tight, at least in the small σ limit, we found it necessary to make the following assumptions on F and the nature of the random noise H. (We abuse notation and interchangeably use F (A) and F (a), where a = vec(A).)

(3.3) • F : R m×n → R p×q is entire, i.e., each component F ij (a) of F has a Taylor series (2.
3) that is absolutely convergent for all a ∈ R mn . • vec(H) ∼ (0, σ 2 Σ), and the elements H ij of H have their kth moment bounded by cσ k for some constant c > 0.

In Theorem 3.2 we show that, under the above assumptions, (3.2) is indeed valid to first order in σ. Specifically,

lim σ→0 E F (A + H) -F (A) 2 F σ = J A Σ 1/2 F .
Thus, if the assumptions (3.3) hold, (3.2) is indeed a valid FOE of the RMS error:

(3.4) E F (A + H) -F (A) 2 F ≈ σ J A Σ 1/2 F .
This is a stochastic analogue of the deterministic FOE in (1.3).

In the derivation of this result we will need the following lemma, whose proof is given in the appendix.

Lemma 3.1. Let ξ(x) represent a power series in multi-index notation which is absolutely convergent for all x ∈ R n :

ξ(x) = v∈N n 0 |v|≥2 α v x v ∈ R, α v ∈ R, x ∈ R n .
If h is a random vector whose elements h i have their kth moment bounded by cσ k , then for any u ∈ R n ,

lim σ→0 E (u T h)ξ(h) σ 2 = 0 and lim σ→0 E ξ(h) 2 σ 2 = 0.
Theorem 3.2. Suppose that F : R m×n → R p×q is entire, and let J A denote the Jacobian (2.2) at A ∈ R m×n . Additionally, let the elements of H be random variables whose kth moments are bounded by cσ k , and

h = vec(H) ∼ (0, σ 2 Σ). Then (3.5) lim σ→0 E F (A + H) -F (A) 2 F σ = J A Σ 1/2 F . Proof. Denote a = vec(A) ∈ R mn . For each element F ij (a) of F , define ξ ij (h) = v∈N mn 0 |v|≥2 ∂ v a F ij (a) v! h v ∈ R.
Then, from (2.3),

F ij (a + h) -F ij (a) 2 = ∇F ij (a) T h + ξ ij (h) 2 = (∇F ij (a) T h) 2 + 2(∇F ij (a) T h)ξ ij (h) + ξ ij (h) 2 .
In the above,

E (∇F ij (a) T h) 2 = V ∇F ij (a) T h = σ 2 ∇F ij (a) T Σ∇F ij (a) = σ 2 Σ 1/2 ∇F ij (a) 2 2
, while, after dividing by σ 2 and taking the limit as σ → 0, from Lemma 3.1, the other two terms tend to 0. Therefore, the limit in (3.5) is

lim σ→0 p i=1 q j=1 E F ij (a + h) -F ij (a) 2 σ = lim σ→0 p i=1 q j=1 σ 2 Σ 1/2 ∇F ij (a) 2 2 σ =    ∇F 11 (a) T . . . ∇F pq (a) T    Σ 1/2 F = J A Σ 1/2 F .
Theorem 3.2 is applicable to many matrix functions, such as matrix polynomials and the matrix exponential. However, in order to generalize the result so that it is applicable to more general classes of functions and random perturbations, a different framework is required. This is the topic of the following section.

3.2.

A more general stochastic condition number. In this section we provide a generalization of Theorem 3.2. Our motivation is the following. First, recall that Theorem 3.2 was derived under the two assumptions (3.3), one of which being that F is entire. This assumption is rather restrictive: many interesting functions in numerical linear algebra are Fréchet differentiable but not entire. Furthermore, even if the assumptions (3.3) required by Theorem 3.2 do hold, more insight about the sensitivity of F can be gained from not only the RMS value of F (A + H) -F (A) F but also its quantiles. For example, we might require that the error remain below a certain threshold with probability 95%, in which case we would work with the 0.95-quantile.

Suppose that

F : Ω ⊆ R m×n → R p×q is Fréchet differentiable at A ∈ Ω and that vec(H) ∼ (0, σ 2 Σ),
without any other assumption on the distribution or the moments of the elements of H. In this more general setting, it is no longer possible to obtain a meaningful first-order RMS estimate as in Theorem 3.2. Consider the following quantity:

(3.6) κΣ ≡ lim sup σ→0 Med F (A + H) -F (A) F σ .
This is analogous to the limit in (3.5), but with the RMS replaced by the median. In analogy to (1.2), we call κΣ the Frobenius-norm stochastic condition number of F at A with respect to (0, Σ) perturbations, or simply the stochastic condition number of F . We show below that

κΣ ≤ 2 J A Σ 1/2 F .
With minor modifications the result is applicable not only to the median but to other quantiles as well; see Remark 3.2 below. The inequality is a consequence of the following theorem. Theorem 3.3. Suppose that F : Ω ⊆ R m×n → R p×q is Fréchet differentiable at A ∈ Ω, and let J A denote the Jacobian (2.2) at A. If vec(H) ∼ (0, σ 2 Σ), then for any τ > 0, lim sup

σ→0 Prob F (A + H) -F (A) F σ ≥ τ ≤ J A Σ 1/2 F τ .
Proof. Denote

P τ σ = Prob F (A + H) -F (A) F σ ≥ τ .
From (2.1) along with (2.6) and (2.7), for any ǫ ∈ (0, 1) we can decompose P τ σ as follows:

P τ σ = Prob F ′ (A)(H) + R(H) F σ ≥ τ ≤ Prob F ′ (A)(H) F σ ≥ τ (1 -ǫ) + Prob R(H) F σ ≥ τ ǫ . (3.7)
Under the assumption that vec(H) ∼ (0, σ 2 Σ), from (2.2), (2.8), and (2.9) we obtain

Prob F ′ (A)(H) F σ ≥ τ (1 -ǫ) ≤ E{ F ′ (A)(H) F } στ (1 -ǫ) = E{ J A vec(H) 2 } στ (1 -ǫ) ≤ J A Σ 1/2 F τ (1 -ǫ) . (3.8)
Next we need to show that for any τ > 0 and any ǫ ∈ (0, 1), the residual term in (3.7) is bounded above by ǫ if σ is sufficiently small. Actually, this was shown by Stewart in [20, Theorem 2.8]. For completeness, we include a short proof here. For any β > 0,

R τ σǫ = Prob R(H) F σ ≥ τ ǫ = Prob R(H) F σ ≥ τ ǫ ∩ H F < β + Prob R(H) F σ ≥ τ ǫ ∩ H F ≥ β ≤ Prob R(H) F σ ≥ τ ǫ ∩ H F < β + Prob{ H F ≥ β} .
From (2.1), for any α > 0, there exists β such that

(3.9) R(H) F H F ≤ α when H F ≤ β.
Therefore, for any α > 0, there is a corresponding β such that

R τ σǫ ≤ Prob α H F σ ≥ τ ǫ ∩ H F < β + Prob{ H F ≥ β} ≤ Prob α H F σ ≥ τ ǫ + Prob{ H F ≥ β} ≤ αE{ H F } στ ǫ + E{ H F } β ≤ α Σ 1/2 F τ ǫ + σ Σ 1/2 F β .
(The last two inequalities follow from (2.8) and (2.10).) Set α = τ ǫ 2 /(2 Σ 1/2 F ) and obtain the corresponding β such that (3.9) holds. Note that β depends on α but is independent of σ. Then, for all

σ ≤ σ * (ǫ) = ǫβ/(2 Σ 1/2 F ), R τ σǫ ≤ ǫ 2 + ǫ 2 = ǫ.
Substituting the above and (3.8) into (3.7) we see that for any τ > 0 and any ǫ ∈ (0, 1), (3.10)

P τ σ ≤ J A Σ 1/2 F τ (1 -ǫ) + ǫ = γ(ǫ) when σ ≤ σ * (ǫ).
It follows that lim sup

σ→0 P τ σ ≤ γ(0) = J A Σ 1/2 F τ .
We can use Theorem 3.3 to bound the condition number κΣ as follows. Setting τ to be the median in (3.10), so that P τ σ ≥ 1/2, we obtain

(3.11) 1 2 ≤ J A Σ 1/2 F Med F (A+H)-F (A) F σ (1 -ǫ) + ǫ.
Rearranging the above shows that for any ǫ ∈ (0, 1/2), when σ is sufficiently small, Med

F (A + H) -F (A) F σ ≤ 2 J A Σ 1/2 F (1 -2ǫ)(1 -ǫ) .
Therefore, as in the proof of Theorem 3.3, we have

(3.12) κΣ = lim sup σ→0 Med{ F (A + H) -F (A) F } σ ≤ 2 J A Σ 1/2 F .
This leads to the the FOE

(3.13) Med F (A + H) -F (A) F 2σ J A Σ 1/2 F .
We conclude this section with two additional generalizations of the stochastic condition number, given in the following remarks.

Remark 3.1. It might also be of interest to measure the sensitivity of a linear function of F , namely,

W T 1 F (A + H) -F (A) W 2 F , for given W 1 ∈ R p×r 1 and W 2 ∈ R q×r 2 .
For instance, we might be interested in the error in just one element of F , in which case W 1 and W 2 would be standard basis vectors. It can be verified that the results of Theorem 3.3 carry through in this case:

lim sup σ→0 Med W T 1 F (A + H) -F (A) W 2 F σ ≤ 2 (W T 2 ⊗ W T 1 )J A Σ 1/2 F ,
where

W T 2 ⊗ W T 1 ∈ R r 1 r 2 ×pq
denotes the Kronecker product. Remark 3.2. We have chosen to define κΣ in terms of the median for two reasons. The first is so that it may be easily comparable with the RMS approach described in section 3.1.

The second is that we think it is legitimate to be interested in a measure of "central tendency" or "typical error." Besides the median, we could bound in the same way any other quantile of the random variable F (A + H) -F (A) F . Replacing 1/2 in the left-hand side of (3.11) by an arbitrary 1q ∈ (0, 1), from (2.5) we obtain

lim sup σ→0 Q q { F (A + H) -F (A) F } σ ≤ 1 1 -q J A Σ 1/2 F .
For example, we might require that the error remain below a certain threshold with probability 95%, in which case we would set q = 0.95. The resulting error estimate is only a factor 10 larger than κΣ .

Sharpness of the bounds.

Compare the bound on the stochastic condition number (3.12) and the resulting FOE (3.13) with their RMS counterparts (3.5) and (3.4). Recall that (3.12) and (3.13) apply much more generally: the only required assumptions are that F be Fréchet differentiable at A and that vec(H) ∼ (0, σ 2 Σ), without any restriction on the distribution of the elements of H. Furthermore, (3.12) and (3.13) can easily be modified to apply to any other quantile-not just the median.

Being more generally applicable, (3.12) and (3.13) are less sharp than (3.5) and (3.4). The additional factor 2 in (3.12) and (3.13) is not significant. Typically one is not interested in the precise value of the error but only in its order of magnitude. Furthermore, our numerical tests indicate that the factor 2 can typically be ignored.

On the other hand, in (3.5) there is equality, whereas (3.12) is merely an upper bound. It would be useful to obtain a nontrivial lower bound for κΣ . In fact, this is impossible to achieve without additional assumptions on the distribution of the elements of H. The following simple example shows that κΣ can be arbitrarily smaller than

J A Σ 1/2 F .
Let H be a random scalar such that H = ±ǫ, each with probability 1/3, ± √ 3, each with probability 1/6, for ǫ > 0. As required in Theorem 3.3, H has mean 0 and finite variance Σ = 1 + 2ǫ 2 /3. Consider the scalar linear function F (A) = A. Then, for all A ∈ R, the stochastic condition number κΣ in (3.6) simplifies to κΣ = Med{|H|} = ǫ, while the upper bound in (3.12) is

2 J A Σ 1/2 F = 2 √ Σ ≥ 2.
Therefore, κΣ J A Σ 1/2 F → 0 as ǫ → 0, showing that no nontrivial lower bound in (3.12) is possible in general. The above example can easily be extended to the case in which the distribution of H is continuous.

The limiting behavior of P τ σ in Theorem 3.3 is given by

Prob F ′ (A)(H) F σ ≥ τ = Prob J A vec(H) 2 2 σ 2 ≥ τ 2 .
Even if one makes very strong assumptions on the nature of H, it seems difficult to derive a useful relationship between J A vec(H) 2 2 /σ 2 and its expected value J A Σ 1/2 2 F . For example, if vec(H) ∼ N (0, σ 2 Σ), then it is straightforward to verify that J A vec(H) 2 F is a linear combination of independent χ 2 random variables with one degree of freedom:

J A vec(H) 2 2 σ 2 = r i=1 σ 2 i χ 2 i ,
where r and σ i are, respectively, the rank and nonzero singular values of J A Σ 1/2 . Tail bounds for linear combination of independent χ 2 random variables do exist (e.g., [START_REF] Laurent | Adaptive estimation of a quadratic functional by model selection[END_REF]Lemma 1] or [4, Theorem 5.2]). For example, in our notation, [15, Lemma 1] shows that for any τ > 0,

Prob    J A vec(H) 2 2 J A Σ 1/2 2 F ≤ 1 -2 τ r i=1 σ 4 i r i=1 σ 2 i    ≤ exp(-τ ). Therefore, if r i=1 σ 4 i ≪ r i=1 σ 2 i ,
then J A vec(H) 2 is very unlikely to be much smaller than J A Σ 1/2 F , so (3.12) and (3.13) are likely to be sharp. Unfortunately, such bounds depend on the coefficients in the linear combination (here the unknown singular values σ i of J A Σ 1/2 ) and/or are generally quite pessimistic. The matter is even more complicated if vec(H) is not normally distributed. To summarize, it is possible to create distributions such that our bound (3.12) and FOE (3.13) are arbitrarily larger than the actual stochastic condition number. Furthermore, even in simple specific cases, it appears to be quite difficult to fully quantify the sharpness of (3.12) and (3.13). Nevertheless, our numerical experiments indicate that these give reasonable order-of-magnitude estimates of error quantiles, and we believe that they are suitable for practical use.

3.4.

Comparison with a deterministic error estimate. Using (2.1) and (2.2), we can write the Frobenius-norm deterministic condition number (1.2) as follows:

(3.14) κ = lim δ→0 sup H F ≤δ F (A + H) -F (A) F δ = sup H F ≤1 F ′ (A)(H) F = J A 2 .
Therefore, in the Frobenius norm, the deterministic FOE (1.3) is

(3.15) F (A + H) -F (A) F J A 2 H F .
In general, it is difficult to compare the above to the stochastic FOE from (3.13):

Med{ F (A + H) -F (A) F } 2σ J A Σ 1/2 F .
In fact, (3.14) and (3.15) are not particularly suitable for measuring the sensitivity of F to random noise since they are based on the condition H F ≤ δ. In general, from the distribution of the elements of H, not much is known about the distribution of H F . Nevertheless, some insight can be gained from some specific cases. For example, suppose that vec(H) ∼ N (0, σ 2 I mn ). In this case, it is highly unlikely that H F lies far from its RMS value, and (3.15) leads to

F (A + H) -F (A) F σ √ mn J A 2 ≡ δ det .
Compare the above to the stochastic FOE (3.13) of the median with Σ = I mn :

Med{ F (A + H) -F (A) F } 2σ J A F ≡ δ med . Because J A ∈ R pq×mn , J A F ≤ J A 2 min{ √ mn,
√ pq}, and it follows that

(3.16) δ med δ det = 2 J A F √ mn J A 2 ≤ 2 min{ √ mn, √ pq} √ mn .
In particular, if

• J A F ≪ J A 2 min{
√ mn, √ pq}, i.e., the Jacobian matrix has a few singular values that are very large relative to its remaining singular values; and/or • pq ≪ mn, i.e., the domain of F is in a much larger space than its range, then δ med /δ det ≪ 1, and there can be a very large difference between the worst-case sensitivity and the typical sensitivity of F to random noise, at least asymptotically in the small σ limit. Numerical examples comparing stochastic and deterministic FOEs are given in section 5.

Estimating the stochastic condition number.

Recall that the deterministic Frobeniusnorm condition number κ in (3.14) involves the operator norm of F ′ (A) or, equivalently, the spectral norm of the Jacobian matrix:

κ = sup H F ≤1 F ′ (A)(H) F = J A 2 .
The above can be computed using the power method; see, e.g., [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Algorithm 3.20]. This requires the evaluation of the Fréchet derivative F ′ (A)(H) and its adjoint at a few different values of H or, equivalently, the computation of a few matrix-vector products with the Jacobian matrix and with its transpose. These can be computed by standard automatic differentiation techniques. (See, for example, [18, section 7.2] or [START_REF] Griewank | Automatic Differentiation of Algorithms: Theory, Implementation, and Application[END_REF] for an introduction to automatic differentiation.) A number of specialized algorithms have also recently been developed for this purpose; see, e.g., [12, section 7], [START_REF] Al-Mohy | The complex step approximation to the Fréchet derivative of a matrix function[END_REF], and [START_REF] Al-Mohy | Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation[END_REF][START_REF] Al-Mohy | Computing the Fréchet derivative of the matrix logarithm and estimating the condition number[END_REF][START_REF] Higham | An improved Schur-Padé algorithm for fractional powers of a matrix and their Fréchet derivatives[END_REF] for methods specific to some important matrix functions.

On the other hand, the upper bound 2 J A Σ 1/2 F on the stochastic condition number in (3.12) involves the matrix Frobenius norm of the scaled Jacobian matrix J A Σ 1/2 ∈ R pq×mn . Several randomized estimators can be used for this purpose. For a recent survey, we recommend [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF]. An attractive feature of these methods is that they are adjoint-free: only matrixvector products with J A (and not its transpose) are required.

In our numerical experiments we use the "small-sample" estimator of Gudmundsson, Kenney, and Laub [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF][START_REF] Kenney | Small-sample statistical condition estimates for general matrix functions[END_REF], defined as

(4.1) η k (J A Σ 1/2 ) = mn k J A Σ 1/2 Q F = mn k k i=1 J A Σ 1/2 q (i) 2 2 ,
where Q = [q (1) , . . . , q (k) ] ∈ R mn×k is the "thin" Q factor in the QR decomposition of an mn × k matrix whose elements are mutually independent standard normal variables. Small values of k are required to obtain good order of magnitude estimates with a high probability.

In our tests we use k = 5. From (4.1) we see that computing η k (J A Σ 1/2 ) reduces to computing k matrix-vector products with J A Σ 1/2 . This entails the computation of u (i) = Σ 1/2 q (i) , which can be very expensive if an appropriate factorization of Σ is not known a priori. Nevertheless, this step can sometimes be performed easily, for example, if Σ is diagonal as in Fletcher's approach discussed in section 3.1. For large-scale problems involving noise following a multivariate normal distribution with general covariance matrices Σ, an iterative strategy for computing u (i) is proposed in [START_REF] Chow | Preconditioned Krylov subspace methods for sampling multivariate Gaussian distributions[END_REF]. Finally, several methods can be used to compute the matrix-vector products J A u (i) , as mentioned above in our discussion of the power method.

Numerical experiments.

To illustrate the theory we give numerical examples involving the solution of linear equations, least squares problems, and a Krylov subspace iterative method. Our tests are meant to compare the effectiveness of the FOEs of the error obtained from the deterministic condition number and the stochastic condition number when vec(H) ∼ (0, σ 2 Σ) for a wide range of values of σ.

Deterministic FOE

(see (3.15)) :

F (A + H) -F (A) F F (A) F J A 2 H F F (A) F . (5.1) Stochastic FOE (see (3.13)) : F (A + H) -F (A) F F (A) F σ J A Σ 1/2 F F (A) F . (5.2)
Note that, H being random, the deterministic FOE is random, while the stochastic FOE is a deterministic quantity.

5.1. Solution of linear systems of equations. First, we provide examples using the problem discussed in section 3.1 that originally motivated the investigations of Turing, Fletcher, and Stewart:

F (A) = A -1 b for nonsingular A ∈ R n×n and b ∈ R n .
In the following examples, n = 500 and A is created from its SVD, A = U SV T , where U ∈ R n×n and V ∈ R n×n are the Q factors in the QR factorization of random matrices and S ∈ R n×n is diagonal with the n singular values of A, logarithmically equally spaced between 10 -ℓ and 1, on its main diagonal. To highlight the effect of the smallest singular value of A, we give results with ℓ = 2 and ℓ = 8. The vector b is formed as follows: b = Ax, where x = cos(1), . . . , cos(n)

T .

We perturb the input A with random noise H satisfying h = vec(H) ∼ (0, σ 2 Σ) for various values of σ. In the first test, the elements H ij are independent normal variables with mean 0 and variance σ 2 ; in other words, h ∼ N (0, σ 2 Σ) with Σ = I n 2 . In the second test, the elements H ij are mutually independent and equal to ±σA ij each with probability 1/2; in other words, h follows a scaled Bernoulli distribution with Σ = diag(vec([A])). This last covariance matrix comes from Fletcher's approach discussed in section 3.1.

For each value of σ we compute 1000 samples of H and of the resulting normwise relative error F (A + H) -F (A) F / F (A) F . We plot the sample median, as well as error bars representing the 5th and 95th sample percentiles, versus σ. Similarly, for the same 1000 random samples of H, we also show the sample median, as well as the 5th and 95th sample percentiles of the deterministic FOE (5.1), versus σ. We also plot the stochastic FOE (5.2), which is linear in σ.

In this simple example, expressions for the Fréchet derivative F ′ (A)(H) = -A -1 Hx and the Jacobian matrix J A = -x T ⊗ A -1 are known. The spectral norm of the Jacobian is J A 2 = x 2 A -1 2 , which is known from the construction of A and b. As for the scaled norm J A Σ 1/2 F , for each given covariance matrix Σ, we use the estimate η k (J A Σ 1/2 ) with k = 5 samples.

Results are plotted in Figure 1. In all cases, for sufficiently small σ, the stochastic FOE is an excellent estimate of the median relative error. Here "sufficiently small" is roughly σ ≈ σ min (A), the smallest singular value of A. We have not found a rigorous explanation for this. The deterministic FOE is roughly two orders of magnitude larger than the stochastic FOE. As discussed in section 3.4, the difference between the two may increase with increasing n.

In these examples, the deterministic FOE, which involves the random variable H F , is very concentrated about its mean: the 5th and 95th sample percentiles essentially overlap on the loglog plot. (This is due to the central limit theorem, H 2 F being a sum of random variables.) For small σ the distribution of the relative error is also very concentrated; see the comments in section 3.3. For larger values of σ, the relative error seems to level off and concentrate at 1. Explaining the large σ asymptotics, as well as obtaining nonasymptotic error estimates, is the subject of ongoing research. 5.2. Linear least squares problems. Our next example involves a linear least squares problem,

F (b) = arg min x b -Ax 2 = A † b,
where A ∈ R m×n has full column rank and A † = (A T A) -1 A T denotes the Moore-Penrose generalized inverse.

We perturb the input b with random noise h ∼ N (0, σ 2 Σ). However, instead of fixing the problem dimensions and varying σ, in this example we fix σ = 10 -4 and vary m and n. Specifically, we test

• n = 10 1 , . . . , 10 5 and m = 100n as well as

• n = 10 and m = 10 1 , . . . , 10 7 .

Although it may seem artificial to create matrices with such a large aspect ratio, we purposefully do so to highlight the fact that when the domain of F is in a much larger space than its range, i.e., when m ≫ n in this case, then there can be a very big difference between the deterministic and stochastic FOEs. (See the discussion in section 3.4.) Once again, A is formed from its SVD: A = U SV T . In each case S is diagonal with entries logarithmically equally spaced between 10 -4 and 1, while U and V are Householder matrices:

u =    cos(1) . . . cos(m)    , v =    sin(1) . . . sin(n)    , U = I -2uu † , V = I -2vv † .
We do not form U , S, V , or A explicitly to avoid running out of computer memory. We let b = A 1, . . .

1 T + v, where v ∼ N (0, σ 2 Σ).
The Jacobian matrix is simply J b = A † , whose spectral norm is known from the construction of A and whose scaled Frobenius norm J A Σ 1/2 F we estimate using η k (J A Σ 1/2 ) with k = 5 samples. For each value of m, we compute 1000 samples of h ∼ N (0, σ 2 Σ). In Figure 2 we plot the sample median as well as the 5th and 95th sample percentiles of the normwise relative error, along with the deterministic and stochastic FOEs, versus m. Once again, in these examples, the stochastic FOE is an excellent example of the median relative error for all values of m tested. On the other hand, the deterministic FOE can be several orders of magnitude larger and sometimes increasingly so with increasing m/n. 5.3. The method of conjugate gradients. As mentioned in the introduction, there has been some recent interest in the data assimilation and optimization communities in the sensitivity of Krylov subspace iterates (e.g., [START_REF] Zhu | Observation sensitivity calculations using the adjoint of the gridpoint statistical interpolation (GSI) analysis system[END_REF][START_REF] Moore | Estimates of analysis and forecast errors derived from the adjoint of 4D-Var[END_REF][START_REF] Moré | Estimating computational noise[END_REF]). Let x (j) denote the jth iterate of the method of CG applied for solving Ax = b, where A ∈ R n×n is symmetric and positive definite. In previous work [START_REF] Gratton | Differentiating the method of conjugate gradients[END_REF], we derived expressions for the Jacobian matrix of x (j) with respect to b. This led to bounds on the (deterministic) condition number of x (j) with respect to perturbations in b. Here we perform a numerical experiment to investigate the sensitivity of CG iterates to random noise in the right-hand side vector b.

Let x (j) = F (b) denote the jth iterate of CG. In this experiment we set j = 10. We perturb the input b with random noise h, where h ∼ (0, σ 2 Σ) for various values of σ. As in the previous examples, we perform tests with h following a multivariate normal distribution and a scaled Bernoulli distribution, and we test both Σ = I n and Σ = diag(vec([b])).

In this example, n = 1000 and the matrix A ∈ R n×n is created via its spectral decomposition: A = U ΛU T . U ∈ R n×n is the Q factor in the QR factorization of a random matrix, and Λ contains the eigenvalues λ i of A on its main diagonal. The eigenvalues λ i are chosen logarithmically equally spaced between 0.1 and 1, with one extra eigenvalue 10 -6 small with respect to the rest of the spectrum. The vector b is taken as b = cos(1), . . . , cos(n)

T .

Because of the small problem dimensions, we can explicitly create the n × n Jacobian matrix J b using automatic differentiation. We then directly compute its spectral norm and scaled Frobenius norm. As in the previous example, for each value of σ we compute 1000 samples of h and of the resulting normwise relative error F (b + h) -F (b) 2 / F (b) 2 . We plot the sample median, as well as the 5th and 95th sample percentiles of the normwise relative error, versus σ, along with the deterministic and stochastic FOEs of the relative error.

Results are plotted in Figure 3. In all four cases, the behavior is very similar. The stochastic FOE is an excellent estimate of the median relative error. This seems to be the case even for large values of σ. In other words, the CG iterates seem to behave linearly with respect to random perturbations. This is very surprising as Krylov methods are known for being highly nonlinear. We have not found a theoretical explanation for this behavior. As before, the deterministic FOE is a few orders of magnitude larger than the stochastic FOE.

Conclusions.

We have defined a notion of stochastic condition number of a matrix function F to random noise H in its input. First we gave sufficient conditions for the firstorder RMS approach from [START_REF] Turing | Rounding-off errors in matrix processes[END_REF][START_REF] Fletcher | Expected conditioning[END_REF] to be applicable. Our sensitivity analysis based on the work of Stewart [START_REF] Stewart | Stochastic perturbation theory[END_REF] applies much more generally-for all Fréchet differentiable functions and any random noise H such that vec(H) ∼ (0, σ 2 Σ), regardless of the distribution of vec(H). Consequently, we proposed the stochastic condition number (3.12). Analogously to the deterministic condition number (1.2), this measures the asymptotic sensitivity of F in the small perturbation limit, but in this setting "small" refers to the norm of the covariance matrix of vec(H). The stochastic condition number can be computed very efficiently by random sampling: few samples are required to obtain excellent order-of-magnitude estimates with probability close to 1. Our numerical experiments indicate that the resulting error estimates are excellent measures of the median error when A is subjected to random noise.

One drawback of the stochastic condition number is that, like the deterministic condition number, it leads only to a first-order error estimate (valid for sufficiently small σ). The behavior of the error for large values of σ is problem dependent. In our examples involving linear systems, the error estimate was only valid for values of σ smaller than roughly the smallest singular value of A. Surprisingly, however, for perturbations to the right-hand side vector in CG, the small σ asymptotics seemed to be descriptive for arbitrary σ. Obtaining global bounds on Prob F (A + H) -F (A) F ≥ τ for these matrix functions is the subject of current research. Finally, in this work we have restricted ourselves to real matrices and real perturbations. An extension to the complex case seems possible, for example, by considering complex random perturbations H = H R + iH I (where i = √ -1 and H R , H I are real random matrices) of complex matrices A ∈ C m×n . This topic is also currently under investigation.

Appendix. Proof of Lemma 3.1.

Proof. For the first part, using linearity of expectation, we obtain

E (u T h)ξ(h) = E

Figure 1 .

 1 Figure 1. Relative error, and deterministic and stochastic FOEs plotted versus σ. The function is F (A) = A -1 b, and the input A is perturbed by H satisfying vec(H) ∼ (0, σ 2 Σ). The dotted vertical lines represent σmin(A), the smallest singular value of A. Top: σmin(A) = 10 -2 . Bottom: σmin(A) = 10 -8 . Left: Normal distribution, Σ = I n 2 . Right: Scaled Bernoulli distribution, Σ = diag(vec([A])).

Figure 2 .

 2 Figure 2. Relative error, and deterministic and stochastic FOEs plotted versus m. The function is F (b) = A † b, and the input b is perturbed by h ∼ N (0, σ 2 Σ). Top: Constant aspect ratio m/n. Bottom: Increasing aspect ratio m/n. Left: Σ = Im. Right: Σ = diag(vec([b])).

Figure 3 .

 3 Figure 3. Relative error, and deterministic and stochastic FOEs plotted versus σ. The function is x10 = F (b), and the input b is perturbed by h ∼ (0, σ 2 Σ) following two distributions. Top: Normal distribution. Bottom: Scaled Bernoulli distribution. Left: Σ = In. Right: Σ = diag(vec([b])).

αα v u j cσ |v|- 1 ≤ c u 1

 11 v u j E{h j h v } .Repeatedly using the Cauchy-Schwarz inequality| E{x 1 x 2 } | ≤ E x 2 1 E x 2 2 and the fact that E h k i ≤ cσ k , we obtain E{h j h v } = E{h j h v 1 1 , . . . , h vn n } ≤ cσ |v|+1 .Therefore, for σ < 1, E (u T h)ξ(h0 as σ → 0. (The last summation is bounded since ξ(x) is absolutely convergent at x = [1, . . . , 1] T .)Similarly, for the second part, E ξ(h) 2 = E converges to 0 as σ → 0.

This research was sup ported by the "Fondation Sciences et Technologies pour l'Ae ´ronautique et l'Espace (FCS STAE)," within the "Re ´sea u The ´matique de Recherche Avance ´e (RTRA),


S Gratton
email: serge.gratton@enseeiht.fr

C W Royer
email: croyer2@wisc.edu

L N Vicente

Z Zhang

Complexity and global rates of trust-region methods based on probabilistic models

Trust-region algorithms have been proved to globally converge with probability one when the accuracy of the trust-region models is imposed with a certain probability conditioning on the iteration history. In this paper, we study their complexity, providing global rates and worst case complexity bounds on the number of iterations (with overwhelmingly high probability), for both first and second order measures of optimality. Such results are essentially the same as the ones known for trust-region methods based on deterministic models. The derivation of the global rates and worst case complexity bounds follows closely from a study of direct-search methods based on the companion notion of probabilistic descent.

Introduction

Trust-region methods form a well established and understood class of methods for the minimization of a nonlinear (possibly nonsmooth) function subject or not to constraints on its variables (see the book [START_REF] Conn | Trust-Region Methods[END_REF] and the recent survey [START_REF] Yuan | Recent avances in trust region algorithms[END_REF]). They have also been comprehensively studied in the context of derivative-free optimization (DFO), where the derivatives of the objective or constraint functions cannot be computed or approximated (see the book [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] and the recent survey [START_REF] Custódio | Methodologies and software for derivativefree optimization[END_REF]). In this paper we focus on the unconstrained minimization of a smooth objective function f : R n → R without using its derivatives.

In the derivative-free setting, trust-region algorithms use sampled points to build a model of the objective function around the current iterate, typically by quadratic interpolation. The quality of these models is measured by the accuracy they provide relatively to a Taylor expansion in a ball B(x, δ) of center x and radius δ. Models that are as accurate as first-order Taylor ones are called fully linear [START_REF] Conn | Geometry of interpolation sets in derivativefree optimization[END_REF][START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. Definition 1.1 Given a function f ∈ C 1 and constants κ ef , κ eg > 0, a C 1 function m : R n → R is called a (κ eg , κ ef )-fully linear model of f on B(x, δ) if, for all s ∈ B(0, δ), |m(s) -f (x + s)| ≤ κ ef δ 2 , ∇m(s) -∇f (x + s) ≤ κ eg δ.

Fully linear models are not necessarily linear or affine functions. Models that are as accurate as second-order Taylor ones are called fully quadratic [START_REF] Conn | Geometry of interpolation sets in derivativefree optimization[END_REF][START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. Similarly, such models are not necessarily quadratic functions. Definition 1.2 Given a function f ∈ C 2 and constants κ ef , κ eg , κ eh > 0, a C 2 function m : R n → R is called a (κ eh , κ eg , κ ef )-fully quadratic model of f on B(x, δ) if, for all s ∈ B(0, δ),

|m(s) -f (x + s)| ≤ κ ef δ 3 , ∇m(s) -∇f (x + s) ≤ κ eg δ 2 , ∇ 2 m(s) -∇ 2 f (x + s) ≤ κ eh δ.
The construction of fully linear/quadratic models based on sampled sets raises a number of geometrical questions. Conn, Scheinberg, and Vicente [START_REF] Conn | Geometry of interpolation sets in derivativefree optimization[END_REF][START_REF] Conn | Geometry of sample sets in derivative-free optimization: Polynomial regression and underdetermined interpolation[END_REF][START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] provided a systematic approach to the subject of deterministic sampling geometry in DFO, establishing error bounds for polynomial interpolation and regression models in terms of a constant measuring the quality or well poisedness of the corresponding sample set (ensuring then the fully linear/quadratic properties). They also showed how to deterministically build or update such sets to ensure that such a constant remains moderate in size. Some numerical studies pioneered by [START_REF] Fasano | On the geometry phase in model-based algorithms for derivative-free optimization[END_REF] have however shown that trust-region methods can tolerate the use of models updated without strict geometry requirements, although it is also known [START_REF] Scheinberg | Self-correcting geometry in model-based algorithms for derivative-free unconstrained optimization[END_REF] that convergence cannot be ensured to first-order critical points without appealing to fully linear models when the size of the model gradient becomes small (a procedure known as the criticality step).

A DFO context of expensive function evaluations often makes it unaffordable to construct a deterministic model that is guaranteed to be fully quadratic, as such a process requires (n + 1)(n + 2)/2 function evaluations. Practical approaches rely on considerably less points (but at least n+ 1 to preserve fully linearity), and use the remaining degrees of freedom to minimize the norm of the model Hessian or its distance to the previous one. The most studied examples use minimum Frobenius type norms [START_REF] Conn | A derivative free optimization algorithm in practice[END_REF][START_REF] Powell | Least Frobenius norm updating of quadratic models that satisfy interpolation conditions[END_REF], yet in [START_REF] Bandeira | Computation of sparse low degree interpolating polynomials and their application to derivative-free optimization[END_REF] it was proposed to apply the theory of sparse ℓ 1 -recovery to build quadratic models based on random sampling. Such models were proved to be fully quadratic with high probability even when considerably less than (n+1)(n+2)/2 points were used, depending on the sparsity of the Hessian of the objective.

Such findings have then called for a probabilistic analysis of derivative-free trust-region algorithms [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF], where the accuracy of the models is only guaranteed with a certain probability. It was shown in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] that the resulting trust-region methods converge with probability one to first and second order critical points. The main purpose of this paper is to establish (with overwhelmingly high probability) the rate under which these methods drive to zero the corresponding criticality measures: the norm of the 2 gradient ∇f(x k ) (in the first-order case) and the minimum σ(x k ) = min{ ∇f(x k ) , -λ min (∇ f(x k ))} between the gradient and the symmetric of the smallest eigenvalue of the Hessian (in the second-order case). We will see that these results match the convergence rates known for deterministic trust-region algorithms. The proofs rely heavily on the technique developed in [START_REF] Gratton | Direct search based on probabilistic descent[END_REF] for establishing global rates and worst case complexity bounds for randomized algorithms in which the new iterate depends on some object (directions in [START_REF] Gratton | Direct search based on probabilistic descent[END_REF], models here) and the quality of the object is favorable with a certain probability. The technique is based on counting the number of iterations for which the quality is favorable and examining the probabilistic behavior of this number. Although the road map for our paper was described in [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Section 6], its actual concretization poses a few delicate issues and, in addition, we go beyond [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Section 6] in other aspects (bounds in expectation, coverage of the second-order case).

Let us now review what is known about the complexity of trust-region methods in the deterministic unconstrained case1 . Using first-order Taylor models, trust-region methods are known [START_REF] Gratton | Recursive trust-region methods for multiscale nonlinear optimization[END_REF] to take at most O(ǫ -2 ) iterations to reduce the ∇f (•) below ǫ ∈ (0, 1) (see also [START_REF] Grapiglia | On the worst-case complexity of nonlinear stepsize control algorithms for convex unconstrained optimization[END_REF] for the corresponding bounds under convexity O(ǫ -1 ) and strong convexity O(-log(ǫ))). Using second-order Taylor models, it has been proved [START_REF] Cartis | Complexity bounds for second-order optimality in unconstrained optimization[END_REF] that at most O(max{ǫ -2 g ǫ -1 h , ǫ -3 h }) iterations are needed to reduce simultaneously ∇f (•) below ǫ g ∈ (0, 1) and -λ min (∇ 2 f (•)) below ǫ h ∈ (0, 1). Later in [START_REF] Grapiglia | Nonlinear stepsize control algorithms: Complexity bounds for first-and second-order optimality[END_REF] it was also proved a complexity bound of the type O(max{ǫ -3 g , ǫ -3 h }). In a recent paper [START_REF] Curtis | A trust region algorithm with a worst-case iteration complexity of O ǫ -3/2 for nonconvex optimization[END_REF] it was suggested a modification to the classical trust-region approach to make it achieving the first-order O(ǫ -1.5 ) bound of cubic regularization methods.

In the derivative-free case, we are also interested in counting function evaluations and to understand the dependence of the complexity bounds in terms of the dimension n of the problem. Using fully linear models instead, the first-order bound [START_REF] Garmanjani | Trust-region methods without using derivatives: Worst case complexity and the non-smooth case[END_REF] is then of the form O(κ -2 eg ǫ -2 ), and since interpolation techniques can ensure κ eg = O( √ n) with at most O( √ n) evaluations per iteration, one recovers the bounds O(nǫ -2 ) for iterations and O(n 2 ǫ -2 ) for function evaluations, also known for direct search [START_REF] Vicente | Worst case complexity of direct search[END_REF]. Using fully quadratic models, the second-order complexity bounds [START_REF] Garmanjani | Trust-region methods without using derivatives: Worst case complexity and the non-smooth case[END_REF][START_REF] Júdice | Trust-Region Methods without using Derivatives: Worst Case Complexity and the Non-smooth Case[END_REF] do not entirely match the derivative-based case as a single tolerance ǫ must be used, being of the form O(n 3 ǫ -3 ) (resp. O(n 5 ǫ -3 )) for measuring the number of iterations (resp. function evaluations) needed to drive σ(•) below ǫ.

We are ready to start presenting our ideas on how to derive global rates and complexity bounds for trust-region methods based on probabilistic models. We will do so in Section 2 for first-order stationarity and in Section 3 for the second-order counterpart. In Section 4 we comment on the extension of our work to other settings.

Complexity of first-order trust-region methods based on probabilistic models

We consider now the scenario analyzed in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] where the models used in a trust-region method are randomly generated at each iteration. As a result, the iterates and trust-region radii produced by the algorithm will also be random. Upper case letters will be then used to designate random variables and lower case their realizations. Hence, m k , x k , δ k will denote respectively the realizations of the random model, iterate, and trust-region radius M k , X k , ∆ k at iteration k.

The random models are then asked to be fully linear with a certain favorable property regardless of the past iteration history. The following definition was proposed in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] to analyze global convergence of the corresponding trust-region methods to first-order critical points.

Definition 2.1 We say that a sequence of random models {M k } is (p)-probabilistically (κ eg , κ ef )fully linear for a corresponding sequence {B(X k , ∆ k )} if the events

S k = {M k is a (κ eg , κ ef )-fully linear model of f on B(X k , ∆ k )}
satisfy P(S 0 ) ≥ p and, for each k ≥ 1, the following submartingale-like condition

P (S k |M 0 , . . . , M k-1 ) ≥ p.
An example is given in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] where, using random matrix theory, it was shown that linear interpolation based on Gaussian sample sets of cardinality n + 1 (with a fixed point and the remaining n points being generated randomly from a standard Gaussian distribution) gives rise to fully linear models with a favorable probability (say p > 1/2). By using increasingly more sample points and building the models by linear regression it is possible to reach a probability as high as desired ([12, Chapter 4]; see also [START_REF] Larson | Stochastic derivative-free optimization using a trust region framework[END_REF]).

Algorithm and assumptions

To simplify the presentation, we describe the trust-region methods under consideration (later given in Algorithm 2.1) for each realization of the model randomness. A few components of these methods are classical, with or without derivatives. At each iteration k, one minimizes a quadratic model

m k (x k + s) = f (x k ) + g ⊤ k s + 1 2 s ⊤ H k s in a trust region of the form B(x k , δ k ).
For global convergence to first-order criticality, the Hessian models are assumed uniformly bounded and the step s k is asked to satisfy a fraction of the model decrease given by the negative model gradient within the trust region. These two assumptions are formalized next.

Assumption 2.1 There exists a positive constant κ bhm such that, for every k, the Hessians

H k of all realizations m k of M k satisfy H k ≤ κ bhm .
Assumption 2.2 For every k, and for all realizations m k of M k (and of X k and ∆ k ), we are able to compute a step s k so that it satisfies a fraction of Cauchy decrease, i.e.,

m(x k ) -m(x k + s k ) ≥ κ f cd 2 g k min g k H k , δ k , (2.1)
for some constant κ f cd ∈ (0, 1], and with the convention that

g k H k = ∞ if H k = 0.
Finally, the step acceptance and trust-region radius update are based on the ratio between the actual decrease in the objective function and the one predicted by the model, namely

ρ k = f (x k ) -f (x k + s k ) m k (x k ) -m k (x k + s k ) .
What is different now from the classical derivative-based case is that some form of criticality step has to be taken into account, where models are recomputed in regions small enough compared to the size of the model gradient. Following [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF], the presentation and the analysis can be Algorithm 2.1: A simple first-order derivative-free trust-region framework Fix parameters η 1 , η 2 , δ max > 0 and 0 < γ 1 < 1 < γ 2 . Select initial x 0 and δ 0 ∈ (0, δ max ).

for k = 0, 1, . . . do Build a quadratic model m k (x k + s) of f , and compute s k by approximately minimizing m k in B(x k , δ k ) so that it satisfies (2.1). If ρ k ≥ η 1 , set x k+1 = x k + s k and δ k+1 = min {γ 2 δ k , δ max } if g k ≥ η 2 δ k , γ 1 δ k otherwise.
Otherwise, set x k+1 = x k and δ k+1 = γ 1 δ k . end significantly simplified if this requirement is mitigated at each iteration. So, in Algorithm 2.1, the trust-region radius is reduced (at iterations where ρ k is large enough and the step is taken) provided δ k is too large compared to g k .

Note that we have slightly extended the framework in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] by using two different parameters (namely γ 1 and γ 2 ) to update the trust-region radius, instead of using a single one and its inverse. As we will see, these parameters are intimately connected to the minimum probability with which the models are required to be fully linear. Also, the safeguard δ max is not used in the analysis of the first-order methods and is only there for coherence with the second-order case (where it appears in the analysis) as well as with [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF].

The algorithm will be analyzed under the following two assumptions on f . As in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] it would be enough to assume continuously differentiability in an enlarged initial level set, but we skip this detail for keeping the presentation simple.

Assumption 2.3

The function f is continuously differentiable on R n , and its gradient is Lipschitz continuous.

Assumption 2.4

The objective function f is bounded from below on R n , and we denote by f low a lower bound.

At this point we can state a fundamental result for establishing the complexity of trust-region methods based on probabilistic models. As in the classical setting of trust-region methods, it ensures that the step is taken if the trust-region radius is small enough compared to the size of the true gradient. There are two differences compared to [2, Lemma 3.2]: first, the result is stated for the true gradient, in alignment to what is needed to establish complexity bounds; second, it is inferred additionally that the trust-region radius is increased under the same condition.

Lemma 2.1 If m k is (κ eg , κ ef )-fully linear on B(x k , δ k ) and δ k < κ eg + max η 2 , κ bhm , 4κ ef κ f cd (1 -η 1 ) -1 ∇f (x k ) , (2.2)
then at the k-th iteration the step is taken (x k+1 = x k + s k ) and δ k is increased.

Proof. From (2.2), one has

κ eg δ k + max η 2 , κ bhm , 4κ ef κ f cd (1 -η 1 ) δ k < ∇f (x k ) ,
and from this and Definition 1.1,

max η 2 , κ bhm , 4κ ef κ f cd (1 -η 1 ) δ k < ∇f (x k ) -g k -∇f (x k ) ≤ g k .
Hence,

δ k < min 1 η 2 , 1 κ bhm , κ f cd (1 -η 1 ) 4κ ef g k ,
and from [12, Proof of Lemma 10.6] we obtain ρ k ≥ η 1 . Since the first term in the minimum gives η 2 δ k < g k , the trust-region radius is increased.

Behavior of the trust-region radius

We will now prove that the sequence of the trust-region radii is square summable and establish an explicit upper bound for the sum. The proof makes use of the set of indices corresponding to iterations where the trust-region radius is increased, that is

K = {k ∈ N : ρ k ≥ η 1 and g k ≥ η 2 δ k } . (2.3)
Remark 2.1 In context of direct search based on probabilistic descent [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Lemma 4.1], it was proved a similar result for the sequence of the step size α k . There, an iteration attains a decrease of the order of α 2 k whenever the step is taken, in which case such direct-search algorithms always increase the step size. In the trust-region context, as we will see in the proof below, a decrease of the order of δ 2 k corresponds only to the iterations where the trust-region radius is increased. There are iterations where the step is taken and the trust-region is decreased, in which case the decrease obtained is not necessarily of this order. Lemma 2.2 For any realization of Algorithm 2.1,

∞ k=0 δ 2 k ≤ β := γ 2 2 1 -γ 2 1 δ 2 0 γ 2 2 + f 0 -f low η ,
where

f 0 = f (x 0 ) and η = η 1 η 2 κ f cd 2 min η 2 κ bhm , 1 .
Proof. We only need to address the case where there are infinitely many iterations at which the trust-region radius is increased (i.e., |K| = ∞). For any k ∈ K, from (2.1),

f (x k ) -f (x k + s k ) ≥ η 1 (m k (x k ) -m k (x k + s k )) ≥ η 1 κ f cd 2 η 2 min η 2 κ bhm , 1 δ 2 k = ηδ 2 k .
Consequently, if we sum a finite number of consecutive iterations in K, we obtain

η j∈K j≤k δ 2 j ≤ j∈K j≤k [f (x j ) -f (x j+1 )] ≤ j≤k [f (x j ) -f (x j+1 )] ≤ f 0 -f (x k+1 ) ≤ f 0 -f low , leading to k∈K δ 2 k ≤ f 0 -f low η .
From now on the proof is as in the proof of [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Lemma 4.1]. Let K = {k 1 , k 2 , . . .} and, for auxiliary reasons, k 0 = -1 and δ -1 = δ 0 /γ 2 . The sum ∞ k=0 δ2 k can thus be rewritten as

∞ k=0 δ 2 k = ∞ i=0 k i+1 k=k i +1 δ 2 k .
Besides, one has for each index i,

δ k ≤ γ 2 (γ 1 ) k-k i -1 δ k i for k = k i + 1, . . . , k i+1 . Hence, k i+1 k=k i +1 δ 2 k ≤ γ 2 2 1 -γ 2 1 δ 2 k i ,
and we finally obtain the desired result:

∞ k=0 δ 2 k ≤ γ 2 2 1 -γ 2 1 ∞ i=0 δ 2 k i ≤ γ 2 2 1 -γ 2 1 δ 2 0 γ 2 2 + f 0 -f low η .
This is a stronger version of the result lim k→∞ δ k = 0 proved in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] for the purpose of global convergence.

Number of iterations with fully linear models

Recall the set of indices (2.3) corresponding to iterations where the step is taken and the trustregion radius is increased, and let y k denote the corresponding indicator (y k = 1 if k ∈ K, y k = 0 otherwise). One sees that δ k and K play the same roles as the step size α k and the set of successful iterations for the analysis of the direct search based on probabilistic descent [START_REF] Gratton | Direct search based on probabilistic descent[END_REF].

Let now ∇f (x k ), with xk ∈ {x 0 , . . . , x k }, represent a minimum norm gradient among ∇f (x 0 ), . . . , ∇f (x k ). Let us also consider the minimum probability

p 0 = ln(γ 1 ) ln (γ 1 /γ 2 ) (2.4)
that will be assumed when applying Definition 2.1. When γ 1 = 1/2 and γ 2 = 2, one has p 0 = 1/2. The next step in the analysis is to show, similarly to [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Lemma 4.2], that if ∇f (x k ) is too large then necessarily not too many iterations benefited from a fully linear model 2 . The binary variable z k below indicates whether m k is (κ eg , κ ef )-fully linear or not.

k-1 l=0 z l ≤ β (min{δ 0 /γ 2 , κ ∇f (x k ) }) 2 + p 0 k, where κ = κ eg + max η 2 , κ bhm , 4κ ef κ f cd (1 -η 1 ) -1 . Proof. For each l ∈ {0, 1, . . . , k -1}, define v l = 1 if δ l < min{γ -1 2 δ 0 , κ ∇f (x k ) }, 0 otherwise.
(2.5)

The proof relies then on z l ≤ (1 -v l ) + v l y l , which is true because, when v l = 1, Lemma 2.1 implies that y l ≥ z l (since ∇f (x k ) ≤ ∇f (x k-1 ) ≤ ∇f (x l ) ). It suffices then to separately prove

k-1 l=0 (1 -v l ) ≤ β (min{δ 0 /γ 2 , κ ∇f (x k ) }) 2 (2.6) and k-1 l=0 v l y l ≤ p 0 k. (2.7)
Inequality (2.6) is justified by Lemma 2.2 and the fact that (2.5) implies

1 -v l ≤ δ 2 l (min{δ 0 /γ 2 , κ ∇f (x k ) }) 2 .
The proof of inequality (2.7) is verbatim as in the proof of [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Lemma 4.2]. It is in here that the particular choice (2.4) for p 0 comes into a play.

Worst case complexity and global rate

From Lemma 2.3, one then has the following inclusion of events

∇f ( Xk ) ≥ ǫ ⊂ k-1 l=0 Z l ≤ β κ 2 ǫ 2 + p 0 k , (2.8)
for any ǫ satisfying 0 < ǫ < δ 0 κγ 2 .

(2.9)

On the other hand the probabilistic behavior of the event on the right-hand side is known from the application of the Chernoff bound to the lower tail of k-1 l=0 Z l (see, e.g., [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Lemma 4.4]), and here is where the conditioning on the past in Assumption 2.1 comes to play a role. Lemma 2.4 Suppose that {M k } is (p)-probabilistically (κ eg , κ ef )-fully linear and λ ∈ (0, p). Then

π k (λ) := P k-1 l=0 Z l ≤ λk ≤ exp - (p -λ) 2 2p k .
Thus, when ǫ < δ 0 /(κγ 2 ) and

k ≥ 2β (p -p 0 )κ 2 ǫ 2 ,
(2.10) the inclusion (2.8) and the monotonicity of π k (•) will give us (setting λ = (p+p 0 )/2 in Lemma 2.4)

P ∇f ( Xk ) ≤ ǫ ≥ P ∇f ( Xk ) < ǫ ≥ 1 -π k β kκ 2 ǫ 2 + p 0 ≥ 1 -π k p -p 0 2 + p 0 ≥ 1 -exp - (p -p 0 ) 2 8p k , (2.11)
leading to the following global rate result.

Theorem 2.1 Suppose that {M k } is (p)-probabilistically (κ eg , κ ef )-fully linear with p > p 0 and

k ≥ 2βγ 2 2 (p -p 0 )δ 2 0 .
(2.12)

Then, the minimum gradient norm ∇f ( Xk ) satisfies

P ∇f ( Xk ) ≤ √ 2β 1 2 (p -p 0 ) 1 2 κ 1 2 1 √ k ≥ 1 -exp - (p -p 0 ) 2 8p k . Proof. Let ǫ = √ 2β 1 2 (p -p 0 ) 1 2 κ 1 2 1 √ k .
Then (2.10) holds with equality, and ǫ < δ 0 /(κγ 2 ) is guaranteed by (2.12). Hence (2.11) gives us the bound.

We have thus proved a global rate of 1/ √ k for the norm of the gradient with overwhelmingly high probability.

Similarly, one can prove a worst-case bound of the order of O(ǫ -2 ) for the first iteration index K ǫ for which ∇f ( XKǫ ) ≤ ǫ, also with overwhelmingly high probability (and we note that K ǫ is a random variable due to the randomness of the models). The proof relies on again applying (2.11), on the observation that P(K ǫ ≤ k) = P( ∇f ( Xk ) ≤ ǫ), and on taking k as

k = 2β (p -p 0 )κ 2 ǫ 2 .
Theorem 2.2 Suppose that {M k } is (p)-probabilistically (κ eg , κ ef )-fully linear with p > p 0 and that ǫ satisfies (2.9). Then, K ǫ satisfies

P K ǫ ≤ 2β (p -p 0 )κ 2 ǫ 2 ≥ 1 -exp - β(p -p 0 )δ 2 4pκ 2 ǫ 2 .
Results in expectation are a natural byproduct of our analysis. It can be shown that E( ∇f ( Xk ) ) is bounded by a function of the order of k -1 2 + exp(-k) up to multiplicative constants (see [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Proposition 5.2]). It is also possible to bound the expected value of K ǫ [START_REF] Royer | Derivative-free Optimization Methods based on Probabilistic and Deterministic Properties: Complexity Analysis and Numerical Relevance[END_REF].

Theorem 2.3 Suppose that {M k } is (p)-probabilistically (κ eg , κ ef )-fully linear with p > p 0 and that ǫ satisfies (2.9). Then,

E (K ǫ ) ≤ c 1 ǫ -2 + 1 1 -exp(-c 2 )
,

where

c 1 = 2β (p -p 0 )κ 2 , c 2 = (p -p 0 ) 2 8p .
Proof. Since K ǫ only takes non-negative integer values, the expected value of K ǫ satisfies:

E (K ǫ ) = ∞ j=0 jP (K ǫ = j) = j≥0 k≥0 k<j P (K ǫ = j) = k≥0 j≥0 k<j P (K ǫ = j) = k≥0 P (K ǫ > k) .
Hence,

E (K ǫ ) = 0≤k<c 1 ǫ -2 P (K ǫ > k) + k≥c 1 ǫ -2 P (K ǫ > k) ≤ c 1 ǫ -2 + 1 + k≥c 1 ǫ -2 P (K ǫ > k) = c 1 ǫ -2 + 1 + k≥c 1 ǫ -2 P ∇f ( Xk ) > ǫ .
For any index k ≥ c 1 ǫ -2 , similar to (2.11), we have

P ∇f ( Xk ) > ǫ ≤ exp(-c 2 k).
As a result,

E (K ǫ ) ≤ c 1 ǫ -2 + 1 + k≥c 1 ǫ -2 exp(-c 2 k) ≤ c 1 ǫ -2 + k≥0 exp(-c 2 k) ≤ c 1 ǫ -2 + 1 1 -exp(-c 2 )
, which proves the desired result.

The obtained bound is thus

O κ -2 ǫ -2 p -p 0 + O(1),
which matches the results of [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF] for line-search methods based on probabilistic models (where p 0 is taken as 1/2). We also emphasize that these expectation bounds exhibit a dependence on the inverse of p -p 0 .

A note on global convergence

Our complexity theory implies (see [START_REF] Gratton | Direct search based on probabilistic descent[END_REF]Proposition 5.1])

P inf k≥0 ∇f (X k ) = 0 = 1.
If we assume for all realizations of Algorithm 2.1 that the iterates never arrive at a stationary point in a finite number of iterations, then the events {lim inf k→∞ ∇f (X k ) = 0} and {inf k≥0 ∇f (X k ) = 0} are identical and we recover the liminf result in probability one of [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF]Theorem 4.2]. Note also that such a result could be derived even more directly by using the argument of [19, Lemma 3.2 and Theorem 3.1]. The paper [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] takes it one step further, establishing also a lim-type result.

3 Complexity of second-order trust-region methods based on probabilistic models

The same proof technology enables us to derive a similar complexity study for the class of trustregion algorithms under consideration but now with the goal of approaching or converging to second-order critical points. To do so, additional assumptions need to be enforced regarding both the quality of the models and the properties of the step resulting from the approximate solution of the trust-region subproblem. We start by the probabilistic counterpart to the concept of fully quadratic models of Definition 1.2.

Definition 3.1 We say that a sequence of random models

{M k } is (p)-probabilistically (κ eh , κ eg , κ ef )-fully quadratic for a corresponding sequence {B(X k , ∆ k )} if the events S k = {M k is a (κ eh , κ eg , κ ef )-fully quadratic model of f on B(X k , ∆ k )}
satisfy P(S 0 ) ≥ p and, for each k ≥ 1, the following submartingale-like condition

P (S k |M 0 , . . . , M k-1 ) ≥ p.
It was shown in [START_REF] Bandeira | Computation of sparse low degree interpolating polynomials and their application to derivative-free optimization[END_REF] how to build fully quadratic models with high probability from quadratic interpolation and uniformly generated sample sets. It is also proved there that such a procedure may recover such models with considerably less than (n + 1)(n + 2)/2 function evaluations when the Hessian of the function is sparse.

Algorithm and assumptions

As before, we consider quadratic models around the iterate x k , with the same definitions for m k , g k , and H k (and the imposition of Assumption 2.1). As curvature is now present in our analysis, we will make use of the notation τ k = λ min (H k ). The solution of the trust-region subproblem has now to be second-order accurate. Assumption 3.1 For every k, and for all realizations m k of M k (and of X k and ∆ k ), we are able to compute a step s k so that it satisfies both a fraction of Cauchy decrease and a fraction of eigendecrease, i.e.,

m(x k ) -m(x k + s k ) ≥ κ f od 2 max g k min g k H k , δ k , -τ k δ 2 k . (3.1)
for some constant κ f od ∈ (0, 1].

The first part of (3.1) can be satisfied by a Cauchy step, while considering a step of norm δ k along an eigenvector of the model Hessian associated with the eigenvalue τ k yields the secondorder decrease in δ 2 k . A step satisfying (3.1) can be obtained by taking the one corresponding to the largest decrease caused in the model value.

Such considerations lead us from Algorithm 2.1 to Algorithm 3.1, preserving the overall structure of the method (and, in particular, the updating rules for the trust-region radius). As mentioned in the Introduction we make use of the second-order criticality measure

σ(x) = max ∇f (x) , -λ min (∇ 2 f (x)) ,
for which a natural estimator at x k is

σ m (x k ) := max { g k , -τ k } .
In the case of a fully quadratic model, the two quantities are related as follows.

Lemma 3.1 [12, Lemma 10.15] If m k is (κ eh , κ eg , κ ef )-fully quadratic on B(x k , δ k ), then |σ(x k ) -σ m (x k )| ≤ κ σ δ k , (3.2)
where κ σ = max {κ eg δ max , κ eh }.

δ k+1 = min {γ 2 δ k , δ max } if σ m (x k ) ≥ η 2 δ k , γ 1 δ k otherwise.
Otherwise, set x k+1 = x k and δ k+1 = γ 1 δ k . end Lemma 3.2 If m k is (κ eh , κ eg , κ ef )-fully quadratic on B(x k , δ k ) and

δ k < κ σ + max η 2 , κ bhm , 4κ ef δ max κ f od (1 -η 1 ) , 4κ ef κ f od (1 -η 1 ) -1 σ(x k ),
then at the k-th iteration the step is taken (x k+1 = x k + s k ) and δ k is increased.

Proof. The proof is similar to the one of Lemma 2.1. Combining (3.2) with the error bound of Lemma 3.1 yields

δ k < min 1 η 2 , 1 κ bhm , κ f od (1 -η 1 ) 4κ ef , κ f od (1 -η 1 ) 4κ ef δ max σ m (x k ).
This allows to directly conclude m that ρ k ≥ η 1 (see [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]Lemma 10.17]). Also, since the first term in the minimum gives η 2 δ k ≤ σ (x k ), the trust-region radius is increased.

Behavior of the trust-region radius

Given the decrease properties now enforced by the algorithm, similar to Lemma 2.2, we can prove that the sequence of the trust-region radii is cube summable. Let K 2nd be the set of indexes corresponding to iterations where the step is taken and the trust-region radius is increased, i.e.,

K 2nd = {k ∈ N : ρ k ≥ η 1 and σ m (x k ) ≥ η 2 δ k } .
Lemma 3.3 For any realization of Algorithm 3.1,

∞ k=0 δ 3 k ≤ β 2nd := γ 3 2 1 -γ 3 1 δ 3 0 γ 3 2 + f 0 -f low η 2nd ,
where

η 2nd = η 1 η 2 κ f od 2 min min η 2 κ bhm , 1 1 δ max , 1 .
Proof. Similar to the proof of Lemma 2.2, we only need to consider the case where

|K 2nd | = ∞. For any k ∈ K 2nd , if σ m (x k ) = g k , we have by Assumption 3.1 that f (x k ) -f (x k+1 ) ≥ η 1 [m k (x k ) -m k (x k + s k )] = η 1 κ f od 2 g k min g k H k , δ k ≥ η 1 η 2 κ f od 2 min η 2 κ bhm , 1 δ 2 k ≥ η 1 η 2 κ f od 2 min η 2 κ bhm , 1 1 δ max δ 3 k .
Meanwhile, if σ m (x k ) = -τ k , a similar reasoning leads to:

f (x k ) -f (x k+1 ) ≥ η 1 η 2 κ f od 2 δ 3 k .
As a result, for any index k ∈ K 2nd , one has

f (x k ) -f (x k+1 ) ≥ η 1 η 2 κ f od 2 min min η 2 κ bhm , 1 1 δ max , 1 δ 3 k = η 2nd δ 3 k .
The rest of the proof follows the lines of Lemma 3.3 replacing the squares of the trust-region radii by cubes.

Number of iterations with fully quadratic models

An upper bound on the number of iterations where the models are fully quadratic is derived similarly as in the first-order case, and all that is required is to define x ˜k ∈ {x 0 , . . . , x k } such that σ(x ˜k) is a minimum value among {σ(x 0 ), . . . , σ(x k )} and z k as the binary variable indicating if m k is (κ eh , κ eg , κ ef )-fully quadratic or not.

Lemma 3.4 Given a realization of Algorithm 3.1 and a positive integer k,

k-1 l=0 z l ≤ β 2nd (min{δ 0 /γ 2 , κ 2nd σ(x k )}) 3 + p 0 k, where κ 2nd = κ σ + max η 2 , κ bhm , 4κ ef δ max κ f od (1 -η 1 ) , 4κ ef κ f od (1 -η 1 ) -1 .

Worst case complexity and global rate

The derivation of the second-order complexity theory is based on observing that Lemma 3.4 implies now

σ( Xk ) ≥ ǫ ⊂ k-1 l=0 Z l ≤ β 2nd κ 3 2nd ǫ 3 + p 0 k , for any ǫ satisfying 0 < ǫ < δ 0 κ 2nd γ 2 . ( 3.3)
Then a result identical to Lemma 2.4 can be ensured considering the definition of Z l based on fully quadratic models and replacing the probabilistic fully linear assumption by the probabilistic fully quadratic one. The rest of the analysis proceeds similarly with minor changes in constants. The global rate is now 1/ 3 √ k, as shown below. Then, the minimum second-order criticality measure σ( Xk ) satisfies

P   σ( Xk ) ≤ 3 √ 2β 1 3 2nd (p -p 0 ) 1 3 κ 1 3 2nd 1 3 √ k   ≥ 1 -exp - (p -p 0 ) 2 8p k .
Similar conclusions as those of Subsection 2.5 can be drawn here regarding global convergence results of the inf and liminf type that can be deduced from the above result, and regarding their interplay with the second-order convergence theory of [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF]. The only difference from the first-order case is the difficulty in obtaining a lim-type result [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF].

A worst-case complexity bound of the order of ǫ -3 is established similarly with overwhelmingly high probability. Theorem 3.2 Suppose that {M k } is (p)-probabilistically (κ eh , κ eg , κ ef )-fully quadratic with p > p 0 and that ǫ satisfies (3.3). Let K ǫ be the first iteration index for which σ( XKǫ ) ≤ ǫ. Then, K ǫ satisfies

P K ǫ ≤ 2β 2nd (p -p 0 )κ 3 2nd ǫ 3 ≥ 1 -exp - β 2nd (p -p 0 ) 4pκ 3 2nd ǫ 3 .
As in Theorem 2.3, we can also bound the expected number of iterations needed to reach the desired accuracy. Theorem 3.3 Suppose at {M k } is (p)-probabilistically (κ eh , κ eg , κ ef )-fully quadratic with p > p 0 and that ǫ satisfies (3.3). Then,

E (K ǫ ) ≤ c 3 ǫ -3 + 1 1 -exp(-c 2 )
, where c 3 = 2β 2nd (p -p 0 )κ 3 2nd , and c 2 is defined as in Theorem 2.3.

Remarks on open questions

Recently a number of papers have appeared proposing and analyzing derivative-free trust-region methods for the unconstrained optimization of a stochastic function. In this setting, what is observable is f (x, ε(ω)), where ε is a random variable. The objective function f (x) may be given by E( f (x, ε)). One approach [START_REF] Shashaani | ASTRO-DF: A class of adaptive sampling trust-region algorithms for derivative-free stochastic optimization[END_REF] extended the framework [START_REF] Conn | Global convergence of general derivativefree trust-region algorithms to first and second order critical points[END_REF] using Sample-Average Approximation (SAA). The number of observations in each Monte Carlo oracle may be up to O(δ -4 k ). First-order global convergence was proved with probability one but for algorithmic parameters that depend on unknown problem constants. Another approach [START_REF] Chen | Stochastic optimization using a trust-region method and random models[END_REF] extended trust-region methods based on probabilistic models [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] to cover also probabilistic estimates of the objective function. In the non-biased case with f (x) = E( f (x, ε)), the probabilistic assumptions can be ensured by SAA within O(δ -4 k ) observations. This approach can also handle biased cases like failures in function evaluations or even processor failures (thus accommodating gradient failures when using finite differences). First-order global convergence was also proved with probability one but again for algorithmic parameters that depend on unknown problem constants. A similar approach [START_REF] Larson | Stochastic derivative-free optimization using a trust region framework[END_REF] led to first-order global convergence in probability (weaker than with probability one), but under more practical assumptions. Very recently, a paper [START_REF] Blanchet | Convergence rate analysis of a stochastic trust region method for nonconvex optimization[END_REF] developed a complexity analysis for the approach in [START_REF] Chen | Stochastic optimization using a trust-region method and random models[END_REF] showing a complexity bound of O(ǫ -2 ) on the expected number of iterations needed to satisfy ∇f (X k ) ≤ ǫ (again under unverifiable assumptions on algorithmic parameters).

It is an open question whether our proof technology can improve upon these stochastic optimization approaches in the sense of establishing global convergence with probability one and global rates and complexity bounds with overwhelmingly high probability without unverifiable assumptions on algorithmic parameters. Another challenging prospect for future work is to develop better rates and bounds for the convex and strongly cases for either deterministic or stochastic functions.

Algorithm 3 . 1 :

 31 A simple second-order derivative-free trust-region framework Fix parameters η 1 , η 2 , δ max > 0 and 0 < γ 1 < 1 < γ 2 . Select initial x 0 and δ 0 ∈ (0, δ max ).for k = 0, 1, 2, . . . do Build a quadratic model m k (x k + s) of f ,and compute s k by approximately minimizing m k in B(x k , δ k ) so that it satisfies (3.1). If ρ k ≥ η 1 , set x k+1 = x k + s k and

Theorem 3 . 1 0 and k ≥ 2β 2nd γ 3 2 (p -p 0 )δ 3 0 .

 3120 Suppose that {M k } is (p)-probabilistically (κ eh , κ eg , κ ef )-fully quadratic with p > p

The notation O(A) will stand for a scalar times A, with this scalar depending solely on the problem considered or constants from the algorithm. The dependence on the problem dimension n will explicitly appear in A when considered appropriate.

As opposed to Lemma 2.2 this result follows more directly from the theory in[START_REF] Gratton | Direct search based on probabilistic descent[END_REF], given that the discrepancy mentioned in Remark 2.1 does not need special treatment.

Support for this author was partially provided by FCT under grants UID/MAT/00324/2013 and P2020 SAICT-PAC/0011/2015. Support for this author was provided by the Hong Kong Polytechnic University under the start-up grant 1-ZVHT.