L Pronzato
email: pronzato@i3s.unice.fr

A Zhigljavsky
email: zhigljavskyaa@cf.ac.uk

E Bukina
email: bukina@i3s.unice.fr

Estimation of spectral bounds in gradient algorithms *

We consider the solution of linear systems of equations Ax = b, with A a symmetric positive-definite matrix in R n×n , through Richardson-type iterations or, equivalently, the minimization of convex quadratic functions (1/2)(Ax, x) -(b, x) with a gradient algorithm. The use of step-sizes asymptotically distributed with the arcsine distribution on the spectrum of A then yields an asymptotic rate of convergence after k < n iterations, k → ∞, that coincides with that of the conjugate-gradient algorithm in the worst case. However, the spectral bounds m and M are generally unknown and thus need to be estimated to allow the construction of simple and cost-effective gradient algorithms with fast convergence. It is the purpose of this paper to analyse the properties of estimators of m and M based on moments of probability measures ν k defined on the spectrum of A and generated by the algorithm on its way towards the optimal solution. A precise analysis of the behavior of the rate of convergence of the algorithm is also given. Two situations are considered: (i) the sequence of step-sizes corresponds to i.i.d. random variables, (ii) they are generated through a dynamical system (fractional parts of the golden ratio) producing a low-discrepancy sequence. In the first case, properties of random walk can be used to prove the convergence of simple spectral bound estimators based on the first moment of ν k . The second option requires a more careful choice of spectral bounds estimators but is shown to produce much less fluctuations for the rate of convergence of the algorithm.

keywords estimation of leading eigenvalues; arcsine distribution; gradient algorithms; conjugate gradient; Fibonacci numbers MSC 65F10; 65F15

Introduction and motivation

For {ν k } ∞ k=0 a sequence of probability measures supported on a real interval [m, M ], the sequence of first moments µ (1)

We consider the behavior of estimators [START_REF] Barzilai | Two-point step size gradient methods[END_REF] and their extensions defined below when ν k are probability measures associated with a gradient algorithm for the minimization of a convex quadratic function with matrix A ∈ R n×n (symmetric positive-definite) and µ (k) 1 = (Ag k , g k )/(g k , g k ), with g k , the gradient at step k, obeying the recurrence equations

g k+1 = g k -γ k Ag k , k = 0, 1, 2 . . . (2)
Here, γ k > 0 is the step-size at iteration k and is determined by some rule that characterizes the algorithm. For instance, γ k = 1/µ (k) 1

for the Steepest Descent (SD) algorithm, γ k = (Ag k , g k )/(A 2 g k , g k ) for the method of Minimum Residues (MR), see [START_REF] Krasnosel | An iteration process with minimal residues[END_REF], [18, p. 134], and γ k = 1/µ (k-1) 1

for the method of Barzilai and Borwein [START_REF] Barzilai | Two-point step size gradient methods[END_REF]. The algorithms considered in [START_REF] Fischer | A stable Richardson iteration method for complex linear systems[END_REF][START_REF] Podvigina | An optimized iterative method for numerical solution of large systems of equations based on the extremal property of zeros of Chebyshev polynomials[END_REF][START_REF] Tal-Ezer | Polynomial approximation of functions of matrices and applications[END_REF] rely on the generation of an infinite sequence of step-sizes γ k (possibly random), such that β k = 1/γ k ∈ [m, M ] for all k, with m and M respectively the minimum and maximum eigenvalues of A. As shown in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF], when the sequence {β k } is asymptotically distributed with the arcsine density in [m, M ], then the asymptotic rate of convergence is competitive compared to that of Conjugate Gradients (CG) [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], Conjugate Residuals (CR) [5, p.547], or other methods based on Krylov spaces, like e.g. MINRES [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF] (it coincides with that exhibited by CG and CR in the worst case, in terms of choice of the starting point and locations of the n -2 internal eigenvalues of A in (m, M ), when the algorithm is stopped before n iterations). Such algorithms, with step-sizes generated externally, are simpler than CR and CG and thus of particular interest in situations where n is so large that the algorithm is stopped well before n iterations (in particular, it is shown in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF] that for some sequences of step-sizes the number of scalar products to be computed out of k iterations only grows as O(log k), see also Sect. 3.3 and 5). The generation of suitable sequences of step-sizes γ k requires, however, the knowledge of the spectral bounds m and M . Since they are usually unknown, it is suggested in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF][START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] to estimate them through the evaluation of moments of probability measures generated by the algorithm itself. It is the purpose of this paper to analyse the asymptotic properties of the estimators (1) of the leading eigenvalues of A. In particular, two situations will be considered: (i) the γ k form a sequence of i.i.d. random variables, (ii) they are constructed from a low discrepancy sequence. In addition, a more precise analysis of the behavior of the rate of convergence of the algorithm than in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF][START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] will be provided (Sect. 4).

A sequence of probability measures associated with a gradient algorithm

Consider a linear system of equations Ax = b ,

where x ∈ R n is an unknown vector, A is a n × n symmetric positive-definite matrix such that 0

< m = inf (z,z)=1 (Az, z) < M = sup (z,z)=1 (Az, z) < ∞ and b is a given vector in R n .
Equivalently, one may consider the minimization of the convex quadratic function

f (x) = 1 2 (Ax, x) -(b, x) . (4)
The gradient g k = Ax k -b then corresponds to the (minus) residual of the system (3) at x k . Richardson-like methods for solving (3) correspond to gradient algorithms for minimizing (4) and obey to the following iterations

x k+1 = x k -γ k g k , k = 0, 1, 2 . . . , (5)
where x 0 ∈ R n is a starting vector and γ k > 0 is the step-size at iteration k. The methods to be considered are of particular interest in situations where the dimension n can be very large (it could even be infinite, with A a self-adjoint operator in a Hilbert space, but we shall restrict the presentation to the finite dimensional situation) and we shall assume that n is much larger than the number of iterations k * needed to achieve the precision required (formally,

k * = o(n) with k * → ∞ in asymptotic considerations).
The iterations (5) can be rewritten in terms of gradients g k , which gives (2), with g 0 = Ax 0 -b ∈ R n the initial gradient. We define the rate of convergence of the algorithm [START_REF] Golub | Matrix Computations[END_REF] towards the solution at iteration j as

r j = (g j+1 , g j+1 ) (g j , g j ) ;
the rate of convergence after k iterations is then

R k = k-1 j=0 r j = (g k , g k ) (g 0 , g 0 ) . (6)
Other convergence rates, which are asymptotically equivalent to r j (see [START_REF] Pronzato | Asymptotic behaviour of a family of gradient algorithms in R d and Hilbert spaces[END_REF]Th. 6]), can also considered. In particular,

r ′ j = f (x j+1 ) -f (x * ) f (x j ) -f (x * ) = (A -1 g j+1 , g j+1 ) (A -1 g j , g j )
is often used when minimizing a quadratic function (4), with x * = A -1 b its minimizer. The method of Steepest-Descent (SD) chooses γ k in (5) that minimizes r ′ k and the method of Minimum Residues (MR) chooses γ k that minimizes r k , both are myopic and only look one-step forward. The method of Conjugate Gradients (CG) minimizes

R ′ k = k-1 j=0 r ′ j = f (x k ) -f (x * ) f (x 0 ) -f (x * ) = (A -1 g k , g k ) (A -1 g 0 , g 0 )
with respect to the sequence γ 0 , γ 1 , . . . , γ k-1 and the method of Conjugate Residuals (CR) does the same with R k ; we shall denote R CR k = min γ 0 ,...,γ k-1 R k . Although CR minimizes R k for all k, for any k < n one may nevertheless have, in the worst case with respect to the starting point x 0 and eigenvalues λ

2 , . . . , λ n-1 , R CR k = R * k , where R * k = R k/2 ∞ + R -k/2 ∞ 2 -2 = C -2 k ρ + 1 ρ -1 , with ρ = M/m the condition number of A, C k (•) the k-th Chebyshev polynomial of the first kind, C k (t) = cos[k arccos(t)] = (1/2)[(t + √ t 2 -1) k + (t - √ t 2 -1) k ],
and

R ∞ = lim k→∞ (R * k ) 1/k = √ ρ -1 √ ρ + 1 2 , (7)
see [START_REF] Forsythe | On the asymptotic directions of the s-dimensional optimum gradient method[END_REF][START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF]. Hence, although one regularly observes values of R CR k that are significantly smaller than R * k for k < n (and although R CR n = 0, that is, the solution is found exactly in n iterations in the exact arithmetic), for any k < n one has max

x 0 ,A R CR k = R * k , with (R * k ) 1/k decreasing monotonically to R ∞ as k → ∞.
As shown in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF], the same asymptotic rate R ∞ can be obtained when the sequence {β k } = {1/γ k } is generated externally with some suitable distribution in [m, M ]; see also Sect. [START_REF] Forsythe | On the asymptotic directions of the s-dimensional optimum gradient method[END_REF].

From (2), the gradient g k after k iterations can be written as

g k = P k (A)g 0 , (8)
where P k denotes the polynomial

P k (A) = (I -γ k-1 A)(I -γ k-2 A) . . . (I -γ 0 A) .
Let m = λ 1 ≤ . . . ≤ λ n = M be the eigenvalues of A and {q 1 , . . . , q n } be the set of corresponding orthonormal eigenvectors (we assume that no information is available about the eigenvalues λ i and eigenvectors q i , i = 1, . . . , n, and that the condition number M/m may be large). When decomposing the initial vector g 0 in the basis {q 1 , . . . , q n } as g 0 = n i=1 α i q i , (8) implies

g k = n i=1 α i P k (λ i )q i (9) for all k ≥ 1. The squared L 2 -norm of g 0 is g 0 2 = (g 0 , g 0 ) = n i=1 α 2 i and the squared L 2 -norm of g k is thus g k 2 = (g k , g k ) = n i=1 α 2 i P 2 k (λ i ) .
The convergence rate (6) after k iterations is then given by

R k = n i=1 α 2 i P 2 k (λ i ) n i=1 α 2 i = n i=1 p (0) i P 2 k (λ i ) ,
where p

(0) i = α 2 i n j=1 α 2 j ≥ 0 and n i=1 p (0) i = 1.
Without loss of generality all α 2 i can be assumed to be strictly positive. Indeed, if α i = 0 for some i then the matrix A = n j=1 λ j q j q ⊤ j can be replaced with à = j =i λ j q j q ⊤ j ; the equality α i = 0 would mean that (A x 0 , q i ) = (b, q i ) (and therefore (A x k , q i ) = (b, q i ) for all k).

From ( 9), (g k , q i ) = α i P k (λ i ). We define

p (k) i = (g k , q i ) 2 (g k , g k ) = α 2 i P 2 k (λ i ) n j=1 α 2 j P 2 k (λ j )
and interpret this as a mass at λ i . Then, the measure ν k defined by its masses ν k (λ i ) = p (k) i at λ = λ i (i = 1, . . . , n) characterizes the normalized vector g k / g k (up to the signs of the (g k , q i ), which are irrelevant for analyzing the behavior of the algorithm). For any real α, define µ (k) α as the α-th moment of the probability measure ν k :

µ (k) α = µ α (ν k ) = n i=1 λ α i p (k) i = (A α g k , g k ) (g k , g k ) . ( 10
)
Using the basic iteration (2), we obtain the following updating formula which expresses the measure ν k+1 through the measure ν k :

p (k+1) i = ν k+1 (λ i ) = α 2 i P 2 k+1 (λ i ) (g k+1 , g k+1 ) = α 2 i (1 -γ k λ i ) 2 P 2 k (λ i ) (g k+1 , g k+1 ) = (1 -γ k λ i ) 2 p (k) i r k , i = 1, . . . , n , (11)
where k ≥ 0 and

r k = (g k+1 , g k+1 ) (g k , g k ) = (g k , g k ) -2γ k (Ag k , g k ) + γ 2 k (A 2 g k , g k ) (g k , g k ) = 1 -2γ k µ (k) 1 + γ 2 k µ (k) 2 . ( 12
)
3 Estimation of the leading eigenvalues of A

Defining the estimators

Take any probability measure ν on [m, M ] with 0 < m < M < ∞ and denote by µ α its moment of order α, µ α = µ α (ν) = M m t α ν(dt) (so that µ α (ν k ) is defined by [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). The Cauchy-Schwarz inequality implies µ α+2 µ α ≥ (µ α+1 ) 2 for any α . Moreover, t(M -t) ≥ 0 for all t ∈ [m, M ] so that

M m t α (M -t) ν(dt) = M µ α -µ α+1 ≥ 0; that is, µ α+1 /µ α ≤ M . Similarly, m ≤ µ α+1 /µ α .
We thus obtain the following chain of inequalities

m ≤ µ (k) 1 ≤ µ (k) 2 µ (k) 1 ≤ µ (k) 3 µ (k) 2 ≤ µ (k) 4 µ (k) 3 ≤ • • • ≤ M , (13)
which are valid for all k = 0, 1, . . . ; note that µ

(k) 0 = 1 for all k.
In what follows we shall restrict our attention to the estimators of m and M defined by

m k = min j=0,...,k µ (j) 1 , M (i) k = max j=0,...,k µ (j) i /µ (j) i-1 , i ≥ 1 . ( 14
)
According to [START_REF] Pronzato | Renormalised steepest descent in Hilbert space converges to a two-point attractor[END_REF], the larger i in

M (i)
k the more precise the estimation of M , which has a significant influence on the behavior of the algorithm, see [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF]. Calculating high order moments has some computational cost, however, and a compromise must be made. The algorithm presented in Sect. 5 uses i = 4.

Let {α k } denote a sequence in [-1, 1] with asymptotic distribution function F α (•) symmetric with respect to zero (different types of sequences will be considered below). We shall consider algorithms defined as follows: we initiate (5) by two SD iterations with

γ k = 1/µ (k) 1 , or two MR iterations with γ k = µ (k) 1 /µ (k) 2 ; for each subsequent iteration the inverse step-size β k = 1/γ k is obtained by rescaling the k-th element α k of the sequence {α k } into [ m k , M k ], that is, β k = α k ( M k -m k )/2 + ( M k + m k )/2 . ( 15
)
The assumption that the β k are generated by symmetric pairs in [m, M ], that is,

β 2j+1 = M + m -β 2j
, is used in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] to derive the expression for the asymptotic convergence rate R of the algorithm, R = lim k→∞ R 1/k k ; see also Sect. 4. This is why we shall also consider the case when ( 15) is replaced by

β 2j = α j ( M 2j -m 2j )/2 + ( M 2j + m 2j )/2 , β 2j+1 = -α j ( M 2j+1 -m 2j+1 )/2 + ( M 2j+1 + m 2j+1 )/2 . ( 16
)
As shown in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF], using the largest β first in a pair (β 2j , β 2j+1 ) permits to improve the monotonicity of the algorithm [START_REF] Golub | Matrix Computations[END_REF]. We shall thus also consider the case when

β 2j = |α j |( M 2j -m 2j )/2 + ( M 2j + m 2j )/2 , β 2j+1 = -|α j |( M 2j+1 -m 2j+1 )/2 + ( M 2j+1 + m 2j+1 )/2 , (17)
Lemma 1 below shows that ( m k + M k )/2 converges to (m+M )/2 so that the symmetry condition with respect to (m + M )/2 will be asymptotically satisfied.

Consistency of m k and M k

The estimators m k and M k satisfy the following asymptotic symmetry property.

Lemma 1 Assume that in algorithm (5) β k = 1/γ k is generated according to one of the rules (15), ( 16) or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], with m k and

M k = M (i)
k given by ( 14) for all k for some i ≥ 1 and {α k } having an asymptotic distribution function

F α (•) in [-1, 1] symmetric with respect to zero. Then we have M -M ∞ = m ∞ -m ≥ 0 , (18)
where

m ∞ = lim k→∞ m k and M ∞ = lim k→∞ M k .
The proof, based on establishing a contradiction if we assume that the symmetry condition ( 18) is violated, is given in Sect. 6.

To obtain a more precise characterization of the limiting behaviors of m k and M k we shall make use of the following property, shown in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF].

Theorem 1 Set β k = 1/γ k (k = 0, 1, . . .) and assume that β k > 0 and β k / ∈ {m, M } for all k and that the sequence {β k } has an asymptotic distribution function F β (•) which is supported on an interval [m ′ , M ′ ] with 0 < m ′ ≤ M ′ < ∞. Suppose, moreover, that this limiting distribution satisfies

log(t -λ) 2 dF β (t) < max log(M -t) 2 dF β (t) , log(t -m) 2 dF β (t) , (19)
for all λ ∈ (m, M ). Then the algorithm (5) associated with the sequence {β k } is such that lim k→∞ ν k (λ i ) = 0 for all i = 2, . . . , n -1. Furthermore, there exist constants C > 0, k 0 > 0 and

0 ≤ θ < 1 such that n-1 i=2 ν k (λ i ) ≤ Cθ k for k > k 0 .
The main condition in Th. 1 is [START_REF] Slater | Gaps and steps for the sequence n θ mod 1[END_REF]; it implies that the ratio P 2 k (λ)/(P 2 k (m) + P 2 k (M )) tends to 0 (as k → ∞) exponentially fast for any λ ∈ (m, M ). It means that once the attraction of the sequence {ν k } to the set of measures supported at m and M is obtained, i.e. ν k (m)+ν k (M ) → 1, it is roughly enough to consider the behavior obtained for two-point measures.

Remark 1 The attraction of the sequence {ν k } to the set of measures supported at m and M does not imply that m k → m and M k → M . Indeed, consider the case when the sequence of step-sizes is self-generated by the algorithm itself. For instance, for SD we have

β k = µ (k)
1 for all k and the limiting measure for {β k } is the two-point measure allocating weights 1/2 at z and M + m -z for some z ∈ (m, M ). The condition [START_REF] Slater | Gaps and steps for the sequence n θ mod 1[END_REF] is then equivalent to z belonging to the stability interval defined in [START_REF] Pronzato | Renormalised steepest descent in Hilbert space converges to a two-point attractor[END_REF][START_REF] Pronzato | Asymptotic behaviour of a family of gradient algorithms in R d and Hilbert spaces[END_REF], see [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF], and m ∞ and M ∞ satisfy (18) but do not coincide with m and M . Remark 2 When the β k are generated by [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], then, under the conditions of Lemma 1, they asymptotically satisfy β 2j+1 = M + m -β 2j . When ν 2j is a two-point measure supported at m and M , we then have ν 2j+2 = ν 2j . Additionally to Th. 1, p (2j) 1 tends to a constant p ∞ as j → ∞, with p ∞ depending on the starting measure ν 0 ( i.e., on the starting point x 0 ) and on the spectrum of A.

We first consider the case when the α k used to generate the step-sizes γ k = 1/β k via (15) are i.i.d. with a suitable distribution.

Theorem 2 Assume that in algorithm [START_REF] Golub | Matrix Computations[END_REF]

β k = 1/γ k satisfies (15) with m k and M k = M (i) k
given by ( 14) for all k for some i ≥ 1 and that the α k are i.i.d. in [-1, 1] with a distribution function F α (•) symmetric with respect to zero. Assume, moreover, that

log(t -u) 2 dF α (t) < log(1 -t) 2 dF α (t) < ∞ for all u ∈ (-1, 1) (20)
and that

F α (1 -x) < 1 for any x > 0. Then, m ∞ = lim k→∞ m k = m and M ∞ = lim k→∞ M k = M almost surely.
The proof is given in Section 6. The idea is roughly as follows. In view of Th. 1, the asymptotic behavior of the measures ν k is very similar to the behavior of measures νk which use the same updating formulas but are supported on the two-point set {m, M }. However, for this sequence of two-point measures νk the sequence of random variables log νk (m) -log(1 -νk (m)) is a random walk and therefore the values νk (m) approach 0 and 1 (with any fixed precision) infinitely often. This implies that the sequence of first moments of νk gets arbitrarily close to m and M infinitely often. The same occurs for the original sequence of measures ν k (proving this requires some technicalities).

As shown in the next theorem, when using

M k = M (i)
k with i ≥ 2 we do not have to use i.i.d. α k to obtain m ∞ = m and M ∞ = M . Moreover, the sequence {β k } can be generated by [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF].

Theorem 3 Assume that in algorithm (5) the β k = 1/γ k are generated according to one of the rules (15), ( 16) or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], with {α k } having an asymptotic distribution function F α (•) in [-1, 1] symmetric with respect to zero, and that m k and M k = M (i) k are given by ( 14) for all k for some i ≥ 2. Assume, moreover, that F α (•) satisfies [START_REF] Tal-Ezer | Polynomial approximation of functions of matrices and applications[END_REF] and that F α (1 -x) < 1 for any x > 0. Then,

m ∞ = lim k→∞ m k = m and M ∞ = lim k→∞ M k = M .
The proof is given in Section 6. The proof of Th. 2 must be modified since, when the α i are not randomly generated, we cannot be sure to have simultaneously a large value of β k and a small value of µ When the α k are i.i.d., the result of Th. 3 holds almost surely. The theorem also covers the case where {α k } is a deterministic sequence, for instance generated via a dynamical system, which permits to obtain sequences of rates r k much less erratic than when using random step-sizes, see Sect. 4. Notice that the proof of Th. 3 does not apply when

M k = M (1) k = max j=0,...,k µ (j) 1 (although simulations seem to indicate consistency of m k and M k = M (1)
k when [START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF] or ( 16) is used, see [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF]).

Controlling the number of updates for m k and M k

The calculations of µ

(k) 1 and µ (k) 2 /µ (k) 1
require the evaluation of several inner products in R n ; therefore, by minimizing the number of iterations where m k or M k are updated we can reduce the computational cost of the algorithm. Updates of m k and M k in the situation of Th. 2 or Th. 3 can be stimulated by the convergence of the measure ν k to the set of measures supported at m and M and by the fact that, on the route to this set of two-point measures, ν k can fluctuate between measures supported at m (when m k + M k > m + M ) and measures supported at M (when m k + M k < m + M ). On the other hand, iterations where m k (resp. M k ) has a good chance to get significantly updated are those for which the next measure ν k+1 will be close to the delta measure at m (resp. at M ). This may happen in particular when β k is the smallest (resp. largest) among all β j , j ≤ k. It can be related to record moments for the α j and we shall thus consider the situation where updates of m k or M k are allowed only at those iterations where α j is a new record.

For any sequence {z k } = z 0 , z 1 , z 2 . . . define the two sequences or record moments

{L min j } = {L min j }[{z k }] and {L max j } = {L max j }[{z k }] by L min 0 = L max 0
= 0 and, for all j ≥ 0,

L min j+1 = min{k > L min j : z k < z L min j } , L max j+1 = min{k > L max j : z k > z L max j } .
We also define the numbers of lower and upper record moments, respectively δ min

j = δ min j [{z k }] and δ max j = δ max j [{z k }], by δ min j = #{i ≥ 0 : L min i ≤ j} = 1 + max{i ≥ 0 : L min i ≤ j} , δ max j = #{i ≥ 0 : L max i ≤ j} = 1 + max{i ≥ 0 : L max i ≤ j} .
For any k, denote by βk the value obtained when m and M are substituted for m k and M k in (15), [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF]. The record moments for { βk } then coincide with those for {α j }, but those for β k may differ since in general m k + M k = m + M . However, due to Lemma 1, the dissimilarity is asymptotically negligible. When (15) is used, upper (resp. lower) record moments for βk coincide with upper (resp. lower) record moments for α k ; when [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] or ( 17) is used, records for βk arrive in pairs, records for β2j and β2j+1 being associated with a α j that becomes new record (either lower or upper).

{α k } is constructed from a low discrepancy sequence

Consider in particular the sequence given by α k = cos(πu k ) for all k ≥ 0, with u k = (k + 1)ϕ mod 1 (the fractional part of (k + 1)ϕ), where ϕ = ( √ 5 + 1)/2 ≃ 1.61803 . . . is the golden ratio (note that the sequence {u k } can equivalently be constructed through the dynamical system u 0 = ϕ-1, u k+1 = (u k +ϕ) mod 1, k ≥ 0). This construction is motivated by the associated rate of convergence of the algorithm, see Sect. [START_REF] Forsythe | On the asymptotic directions of the s-dimensional optimum gradient method[END_REF] } can then be expressed in terms of the Fibonacci numbers F j as follows: L min j = F 2j+1 -1, L max j = F 2j+2 -1 for j = 0, 1, . . . This directly follows from the following two classical results of the theory of Diophantine approximations: (i) for the sequence {kζ mod 1}, with any irrational ζ, the successive minimal and maximal values occur when k = q in the denominator of a convergent p/q for ζ in the standard continued fraction expansion of ζ, see [START_REF] Slater | Gaps and steps for the sequence n θ mod 1[END_REF]; (ii) the convergents of ζ = ϕ -1 are F j /F j+1 for j > 1. The number of upper record moments for {α j } is then δ max k [{α j }] = 1 + max{j : F 2j+2 -1 ≤ k} = 1 + max{j : F 2j+2 ≤ k + 1}. Similar to Proposition 7 in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF] we can show that for all k > 1 we have the inequalities

C 0 log k -1 < δ k < C 0 log k + 1 with C 0 = 1/[2 log(ϕ)] ≃ 1.039. This yields the asymptotic relation δ k = C 0 log k + O(1) as k → ∞.
When we use [START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF], the sequence of record moments for { βk } is the same as for {α k }. If we use [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], then the sequences of record moments for βk are {L given by [START_REF] Pronzato | Asymptotic behaviour of a family of gradient algorithms in R d and Hilbert spaces[END_REF], the second to

= 2(F j+1 -1) if j is even and L max j = 2F j+1 -1 if j is odd, L min j = 2(F j+2 -1) if j is even and L max j = 2F j+2 -1 if j is odd
m 2j+1 = min j∈Lα µ (2j+1) 1 and M 2j+2 = max j∈Lα µ (2j+2) 4 µ (2j+2) 3 , (21)
where

L α = {L min j }[{α k }] {L max j }[{α k }]
denotes the sequence of lower and upper record moments for {α k }. This construction is motivated by the fact that when α j becomes a new (lower or upper) record, then β 2j+1 is large and µ (2j+1) 1 has a good chance to be small, while β 2j+2 is small and µ

(2j+2) 4 /µ (2j+2) 3
has a good chance to be large. The precision of the estimation of m and M given by the two estimators is compared in Fig. 1a. Figure 1b indicates the number of record moments for {α k } together with the number of iterations where m k and M k are updated. One may notice that both m 2j+1 and M 2j+2 are updated each time α j is a new record.

Fluctuations of the sequence of convergence rates

When ν k is a two-point measure supported at m and M , applying two successive iterations [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF] with

β k+1 = M + m -β k yields ν k+2 = ν k ,
the particular measure ν k , r k r k+1 = R 2 2 (β k ), where R 2 (β) = (M -β)(β -m) β(M + m -β) .
This is the key-point used in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] to prove the following.

Theorem 4 Assume that the conditions of Th. 1 are satisfied and that, moreover, the β k are generated by symmetric pairs for large k; that is, β 2j+1 = M + m -β 2j for all j ≥ j 0 , with

β 2j ∈ [m + ε, M -ε] for some ε ∈ (0, (M -m)/2). Then, lim k→∞ 1 k log R k = log (M -t) (t -m) t (m + M -t) dF β (t) = log (t -m) 2 t 2 dF β (t) , (22)
where R k is defined by [START_REF] Hall | Martingale Limit Theory and Its Applications[END_REF].

Th. 4 applies in particular when F β (•) has the arcsine density

f ǫ (•) on [m + ǫ, M -ǫ], f ǫ (β) = 1 π (β -m -ǫ)(M -ǫ -β) , with ǫ < (M -m)/2. In that case, as k → ∞, R 1/k k → R arcsine,ǫ = exp M -ǫ m+ǫ log (β -m) 2 β 2 f ǫ (β) dβ = M -m + 2 ǫ(M -m -ǫ) M + m + 2 (M -ǫ)(m + ǫ) 2 (23) = R ∞ (1 + 4 ǫ(M -m)) + O(ǫ) , ǫ → 0 ,
with R ∞ given by [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], see [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF]. In the rest of the section we are interested in the extension of Th. 4 to the case where the β k are generated by ( 16) or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] with estimated m k and M k and to the fluctuations of R 1/k k along its way to its limiting value.

The fact that in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF][START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] m k and M k are estimated brings a slight difference with Th. 4 in terms of asymptotic rate of convergence. This difference is marginal, however, as shown in the next theorem.

Theorem 5 Assume that in algorithm (5) the β k are generated by pairs as in [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF] or [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], with m k and M k = M (i) k given by ( 14) for all k for some i ≥ 2, and that the α j have an asymptotic distribution function

F α (•) in [-1 + δ, 1 -δ], 0 < δ < 1/2,
symmetric with respect to zero and satisfying [START_REF] Tal-Ezer | Polynomial approximation of functions of matrices and applications[END_REF] and F α (1 -δ -x) < 1 for any x > 0. Then the result [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF] of Th. 4 remains valid, with

F β (•) a rescaled version of F α (•) in some interval [m + ǫ ′ , M -ǫ ′ ], that is, dF β (t) = 2 M -m -2ǫ ′ dF α 2t -m -M M -m -2ǫ ′ , (24)
where

ǫ ′ satisfies δ(M -m)/2 ≤ ǫ ′ ≤ √ M ( √ M + √ m)δ/2 + O(δ 2 ) . ( 25
)
The proof is given in Sect. 6. The fact that the α k now lie in [-1 + δ, 1 -δ] with δ > 0 makes it necessary to slightly modify the proof of Th. 3.

Consider now the fluctuations of the asymptotic convergence rate around its limiting value. Suppose that m and M are perfectly estimated, so that [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF][START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF], with the β 2j having a distribution F β (•) satisfying the condition in Th. 4. Since the β k are exactly symmetric in [m, M ], we can take δ 1 = 0 in the proof of Th. 5, so that

β 2j+1 = M + m -β 2j in
log[R 2 (β 2j )] -Bθ 2j < 1 2 log(r 2j r 2j+1 ) < log[R 2 (β 2j )] + Bθ 2j
for some B > 0 and j > j 1 large enough, which gives

1 2 log R 2i R 2j 1 - i-1 j=j 1 log[R 2 (β 2j )] < ∞ j=0 Bθ 2j = B 1 -θ ,
with θ as in Th. 1. Define

L α = log[R 2 (t)] dF β (t) = log 1 + t t + M +m M -m 2 dF α (t) , V α =   log 1 + t t + M +m M -m 2 -L α   2 dF α (t) .
(Note that both quantities are well defined when

F α (•) is concentrated on [-1 + δ, 1 -δ], δ > 0.)
When the α k are i.i.d. with the distribution F α (•), we have, for i → ∞,

1 i log R 2i a.s. → L α , √ i log √ R 2i i -L α d → ξ ∼ N (0, V α ) ,
and

R 1/2i 2i a.s. → exp(L α ) , √ i R 1/2i 2i -exp(L α ) d → ξ ∼ N (0, V α exp(2L α )) .
Moreover, from the law of the iterated logarithm, lim sup

i→∞ log √ R 2i -iL α √ 2iV α √ log log i = 1 a.s. and lim inf i→∞ log √ R 2i -iL α √ 2iV α √ log log i = -1 a.s. ,
implying that, for any ε > 0,

R 1/2i 2i > exp L α + (1 -ε) √ 2V α log log i √ i (26)
infinitely often (a.s.), and thus indicating that the fluctuations of the normalized convergence rate R

1/2i 2i are unavoidably large. Suppose now that {α k } is constructed from a low-discrepancy sequence, as in Sect. 3.3. Then

1 i -j 1 i-1 j=j 1 log[R 2 (β 2j )] -L α < C α log i i
for some large enough j 1 (see the proof of Th. 5) and some constant C α depending on the sequence considered, see, e.g., [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]. Therefore, in that case the normalized convergence rate satisfies R

1/2i 2i / exp(L α ) < i Cα/i and shows much less fluctuations on its route to its limiting value exp(L α ) than when the α k are i.i.d. random variables. Indeed, denoting D 2i the difference R

1/2i 2i -exp(L α ), we have (D 2i ) i.i.d. |(D 2i ) LDS | > exp (1-ε) √ 2Vα log log i √ i - 1
i Cα/i -1 i.o. a.s. for any ε > 0 where the right-hand side behaves like c √ i log log i/ log i as i → ∞ (c = (1 -ε) √ 2V α /C α ). This justifies the preference given to low discrepancy sequences over random sequences in the algorithms presented in [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF]. One of them is summarized in Sect. 5.

Example 2

We consider the same problem as in Sect. 3.3 (we take now x 0 = (10 5 , 1, 1, . . . , 1) ⊤ to slow down convergence and better illustrate the different behaviors for the two types of stepsize sequences). The β k are generated by [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] with m k and M k = M (4) k given by [START_REF] Pronzato | Asymptotic behaviour of a family of gradient algorithms in R d and Hilbert spaces[END_REF]. Figure 2 shows R 1/2i 2i as a function of i for the cases when {α k } is the low discrepancy sequence given by α k = cos(π[(k + 1)ϕ mod 1]) for all k ≥ 0 and when {α k } is a sequence of i.i.d. random variables having the arcsine distribution; both sequences are generated on [-1 + δ; 1 -δ] with δ = 0.005.

Prototype algorithm and simulation results

The estimation of the spectral bounds m and M via (21) permits to construct a gradient algorithm which is quite parsimonious in terms of number of computations of inner products. In order to avoid using multiplications by A when calculating µ Stage II (iterations)

µ (k) 1 = (Ag k , g k ) (g k , g k ) = β k 1 - (g k , g k+1 ) (g k , g k )
II.1 If M k > M k-1 then set β k = M k . Otherwise set β k = m k +( M k -m k )z j and j ← j +1. II.2 Set x k+1 = x k -(1/β k )g k and g k+1 = Ax k+1 -b. II.3 If j -2 ∈ L max then compute m k+1 = min{ m k , µ (k)
1 }, M k+1 = max{ M k , µ The stopping rule used by the algorithm is simply (g k , g k ) < ǫ for some given ǫ. The value of (g k , g k ) is available for k such that j -2 ∈ L max , at such iterations we can thus check the condition (g k , g k ) < ǫ directly. Since M k / m k provides an under-estimate for ρ, and hence and under-estimate for R ∞ given by ( 7), we can thus estimate the number of iterations still remaining to achieve the required precision. The proposed algorithm only requires one matrix-vector multiplication per iteration (used to calculate the gradient g k = Ax k -b), like other gradient methods and Krylov-space based algorithms, in particular CR and CG. When A is sparse, the computation of inner products also contributes significantly to the total computational cost; when using parallel computing with distributed memory machines, it may even yield the main contribution to the efficiency loss, see [START_REF] Van Der | Iterative Methods for Large Linear Systems[END_REF]Sect. 4.4]. The standard formulation of CG (and also CR) requires the computation of two inner products per iteration (some sophisticated versions of CG compute these two inner products in parallel, at the possible expense of a slight increase of storage and maybe reduced numerical stability, see for instance [START_REF] Meurant | The block preconditioned conjugate gradient method on vector computers[END_REF][START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF]). The prototype algorithm above requires the computation of four inner products in the initial two iterations and then four inner products (possibly computed in parallel) each time the estimates m k and M k are updated. This is done when j -2 ∈ L. Therefore, the total number of inner products computed within k + 1 steps of the proposed algorithm is equal to N k = 4 + 4δ j where j = j(k) is defined by the algorithm and δ j satisfies δ j = log j/ log(ϕ) + O(1) as j → ∞. This and the fact that k/j(k) → 1 as k → ∞ imply that the number of inner products computed within k + 1 steps is approximately 4 + 4 log k/ log ϕ ≃ 4 + 8.31 log k.

Simulation results

Example 3 In this artificial example, A is diagonal with m = 1, M = 1 000, n = 1 000 and b = Ac with c random (uniformly distributed on the unit n-dimensional sphere S n ). We consider two configurations for the eigenvalues of A and starting point x 0 . In the first case, the n eigenvalues are uniformly distributed in [m, M ] and x 0 is random (uniformly distributed on S n ). The second configuration corresponds to the worst-case situation for n -1 steps of the Conjugate Residual (CR) algorithm: the eigenvalues are

λ i = (M + m)/2 + (M -m)/2 cos[π(i -1)/(n - 1
)] for i = 1, . . . , n and x 0 is such that the α 2 i in the decomposition (9) are proportional to τ 2 1 = 1/2λ 1 , τ 2 j = 1/λ j for j = 2, . . . , n -1 and τ 2 n = 1/2λ n , see [START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF] for details. Figures 3(a) and 3(b) present the evolution of log 10 g k as a function of k for the algorithm above (Alg) and the CR algorithm in the two configurations respectively. Since CR is optimal for R k given by ( 6), our method is not competitive in this respect. On the other hand, the evolution of log 10 g k as a function of the number of inner products computed is plotted in Figures 3(c We set b = Ac, c and x 0 are uniformly distributed on S n . Figure 5 presents the evolution of log 10 g k against the number inner products in Examples 2 and 3, computed for the algorithm above (Alg) and the CR algorithm. The reduced complexity of the proposed algorithm compared to CR is manifest.

Proofs

Proof of Lemma 1.

Since { m k } forms a non-increasing sequence bounded from below by m, m k → m∞ as k → ∞ for some m∞ ≥ m.

Similarly, M k → M∞ for some M∞ ≤ M . Denote ǫ1 = m∞ -m, ǫ2 = M -M∞, ǫ1, ǫ2 ≥ 0. Suppose that 0 ≤ ǫ1 < ǫ2. Then, the asymptotic distribution of the sequence {β k } is biased towards m. From [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF], the sequence {ν k } tends to concentrate at M so that M k → M as k → ∞, implying ǫ2 = 0, which contradicts the assumption 0 ≤ ǫ1 < ǫ2. Similarly, the assumption ǫ1 > ǫ2 ≥ 0 leads to a contradiction; therefore, M -M∞ = m∞ -m.

Proof of Th. 2.

From Lemma 1, m∞ = m + ǫ and M∞ = M -ǫ for some ǫ ≥ 0. Assuming that ǫ > 0, we show that this leads to a contradiction. The proof is in two steps. In (i) we show that m k is repeatedly updated; in (ii) we show that this implies m∞ = m and thus ǫ = 0.

(i) We have from ( 11) p

  k (dt) gives obvious estimators of m and M through m k = min j=0,...,k µ

1 ,

 1 see part (ii) of the proof of Th. 2. As a consequence, i = 1 does not guarantee the convergence of m k and M k = M (i) k to respectively m and M and we now have to use a more precise estimator of M with i > 1.

3. 3 . 1

 31 {α k } forms an i.i.d. sequence of random variables When the α k are i.i.d., the numbers of lower and upper record moments δ min j [{α k }] and δ max j [{α k }] satisfy δ min j / log j → 1 and δ max j / log j → 1 almost surely, see [2, p. 258]. Therefore, when we use (15), δ min j [{ βk }]/ log j → 1 and δ max j [{ βk }]/ log j → 1, whereas δ min j [{ βk }]/ log j → 2 and δ max j [{ βk }]/ log j → 2 when we use (16) or[START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF].

  min j }[{ βk }] = {0, 1, 3, 5, 9, 15, . . .} and {L max j }[{ βk }] = {0, 2, 4, 8, 14, . . .}, with L min j+1 = L max j + 1 for j = 0, 1, . . . and, in terms of Fibonacci numbers, L max j = 2(F j+2 -1) for j = 0, 1, . . . The number δ k of upper record moments for { βj } thus satisfies δ k / log k → 2C 0 = 1/ log(ϕ) ≃ 2.078 as k → ∞; the same is obviously true for the number of lower record moments. If we use (16), the sequences of record moments for βk are {L min j }[{ βk }] = {0, 3, 4, 9, 14, 25 . . .} and {L max j }[{ βk }] = {0, 1, 2, 5, 8, 15, 24, . . .}; that is, in terms of Fibonacci numbers, L max j

Example 1 2 ,

 12 . The numbers of upper and lower record moments satisfy again δ k / log k → 1/ log(ϕ) ≃ 2.078 as k → ∞. We set n = 800, m = 1 and M = 1000, the eigenvalues of A are uniformly distributed in [m, M ] and b = Ac in[START_REF] Fischer | A stable Richardson iteration method for complex linear systems[END_REF][START_REF] Forsythe | On the asymptotic directions of the s-dimensional optimum gradient method[END_REF] with c uniformly distributed on the unit n-dimensional sphere S n . We apply the gradient iterations (5) with x 0 uniformly distributed on S n ; the first two iterations correspond to the method of minimum residues with γ k = µ k = 0, 1, and the subsequent iterations use {γ k } = {1/β k }, where the β k are generated via (17) using the low discrepancy sequence {α k } above: α k = cos(π[(k + 1)ϕ mod 1]) for all k ≥ 0. We compare the behaviors of two estimators of m and M . The first one corresponds to m k and M (4) k

Figure 1 :

 1 Figure 1: (a) Evolution of log 10 (M -M k ) and log 10 ( m k -m) with m k and M k = M (4) k given by (14), respectively in dotted and dash-dotted lines, and evolution of log 10 (M -M k ) and log 10 ( m k -m) with m k and M k given by (21), respectively in solid and dashes lines. (b) Number of record moments δ max k [{α j }] (hexagrams) and δ min k [{α j }] (circles) as functions of k, the number of iterations where m k and M k are updated are indicated respectively by stars and dots, the solid line corresponds to (log k)/[2 log(ϕ)] ≃ 1.039 log k.

Figure 2 :I. 5

 25 Figure 2: Rate R 1/2i 2i as a function of i for α k = cos(π[(k + 1)ϕ mod 1]) for all k ≥ 0 (solid line) and for {α k } a random i.i.d. sequence having the arcsine distribution (dash-dotted line). The dashed line indicates the bound R arcsine,ǫ ′ given by (23) with ǫ ′ = δ(M -m)/2 (δ = 5 • 10 -3 ), see (25); the dotted line denotes the right-hand side of (26).

  check the stopping rule (g k , g k ) ≤ ǫ. Otherwise set m k+1 = m k , M k+1 = M k . II.4 Set k ← k + 1 and return to Step II.1.

Figure 3 : 4 . 5 A

 345 Figure 3: log 10 g k as a function of k (top), log 10 g k against the number of inner products computed (bottom) in Example 3. The eigenvalues of A are uniformly distributed in [m, M ] = [1, 1 000] in (a) and (c); (b) and (d) correspond to the worst-case situation for the CR algorithm.

Figure 4 :

 4 Figure 4: Non-zero elements of A in Examples 4 (left) and 5 (right).

Figure 5 :

 5 Figure 5: log 10 g k against the number of inner products computed in Examples 4 (left) and 5 (right).

. 1 -

 1 When β k satisfies[START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF] with Fα(•) satisfying (20) and ǫ > 0, then[START_REF] Slater | Gaps and steps for the sequence n θ mod 1[END_REF] is satisfied so that p δ k with

  . The corresponding sequences of record moments are {L min {1, 1, 2, 3, 5, 8, 13, 21, 34 . . .} the sequence of Fibonacci numbers, with exact expression F N = (ϕ N -(-1/ϕ) N )/ √ 5. {L min j } and {L max j

j }[{α k }] = {0, 1, 4, 12, 33, . . .} and {L max j }[{α k }] = {0, 2, 7, 20, 54, . . .}. Denote {F N } ∞ N =1 =

* Part of this work was accomplished while the first two authors were invited at the Isaac Newton Institute for Mathematical Sciences, Cambridge, UK; the support of the INI and of CNRS is gratefully acknowledged. The work of E. Bukina was partially supported by the EU through a Marie-Curie Fellowship (EST-SIGNAL program: http://est-signal.i3s.unice.fr) under the contract Nb. MEST-CT-2005-021175.

which can easily be derived from [START_REF] Embrechts | Modelling Extremal Events[END_REF]. We generate the β k = 1/γ k according to [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] using the low discrepancy sequence α k = cos(π[(k + 1)ϕ mod 1]) for all k ≥ 0. The construction [START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] tends to favor the estimation of m against that of M , see [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF], which results in the concentration of ν k at M . From [START_REF] Podvigina | An optimized iterative method for numerical solution of large systems of equations based on the extremal property of zeros of Chebyshev polynomials[END_REF], the convergence is monotonic at step k (i.e., r k < 1) when

). When ν k gets close to the delta measure at M , this ratio becomes close to M/2 when M/m is large, and the monotonicity condition is violated frequently (approximately every second iteration). In the algorithm proposed below this is avoided by forcing ν k to become concentrated at m rather than M , the monotonicity property r k < 1 being always satisfied when µ

This can be achieved by using a step with large β k when ν k-1 becomes close to the delta measure at M . In practice, we simply use β k = M k when we observe M k > M k-1 . The algorithm is summarized below; its MATLAB implementation is available at http://www.i3s.unice.fr/ ~pronzato/Matlab/goldenArcsineQ.m. We define v j = ϕ(j + 1) mod 1 and, for j = 0, 1, . . . we set where

0 ≤ δ k ≤ Dθ k when k > k0 for some k0, D > 0 and 0 ≤ θ < 1. Denoting

Suppose that there is no update of m k and M k after some k2 and denote

Then, for j > k2, {ξj} forms a sequence of i.i.d. random variables and (28) indicates that, for k > k2, log Q k+1log Q k 2 behaves like a random walk. The random variables ξj have mean

From Lemma 1, ǫm = ǫM = ǫ. Since M (ǫ, ǫ) = 0, the random walk has no drift and we have lim sup k→∞ log Q k = lim inf k→∞ log Q k = ∞ a.s., which contradicts the assumption of no update of m k and M k after iteration k2. (One may notice that ξ2j + ξ2j+1 = 0 when the β k are generated according to [START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF], so that the argument cannot be used in that case.) Suppose now that there is no update of m k for k > k2 and denote ǫm = m k 2 -m. From the argument above, M k is repeatedly updated, which, from Lemma 1, is only possible if ǫ M,k = M -M k > ǫm. The ξj, for j > k2, are now neither independent nor identically distributed, but IE(ξj|Fj-1) = M (ǫm, ǫM,j) < 0, with {Fj} the sequence of σ-fields σ(α0, α1, . . . , αj). From (28), the sequence of log Q k+1log Q k 2 , k > k2, thus forms a supermartingale relative to {Fj}. Consider now S k = k j=k 2 ξj -IE(ξj|Fj-1), which forms a martingale sequence. Since we assume that ǫ > 0, the increments |ξj -IE(ξj|Fj-1)| are bounded and S k / √ k satisfies the central limit theorem (see, e.g., [START_REF] Hall | Martingale Limit Theory and Its Applications[END_REF]). This implies that lim inf k→∞ log Q k = -∞ and therefore M∞ = M , which contradicts the assumption of ǫ > 0. We have thus proved that m k is updated infinitely often.

(ii) Similarly to (i), m∞ = m + ǫ and M∞ = M -ǫ with ǫ > 0, β k satisfying [START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF] with Fα(•) satisfying [START_REF] Tal-Ezer | Polynomial approximation of functions of matrices and applications[END_REF] imply that ( 19) is satisfied. Denoting p k = p (k) 1 , direct calculations using [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF] and [START_REF] Podvigina | An optimized iterative method for numerical solution of large systems of equations based on the extremal property of zeros of Chebyshev polynomials[END_REF] give

when k > k0, for some constants D ≤ 0 and 0 ≤ θ < 1, where

] and β k ∈ (m + ǫ, M -ǫ) by construction. (Notice that the term within square brackets on the right-hand side of ( 29) is an increasing function of both p k and β k .) The fact that Fα(1 -x) < 1 for any x > 0 implies that lim sup k→∞ β k = M -ǫ. We have shown in (i) that m k is updated infinitely often. Therefore, for any δ1, δ2 > 0, there exists a subsequence {ji} such that βj i > M -ǫ -δ1 and µ

< m + ǫ + δ2. For a two-point measure supported at m and M , this second inequality implies pj i > (M -ǫ -m)/(M -m) -δ3, with δ3 → 0 as δ2 → 0. Due to Th. 1, we thus have

for some constant B ≤ 0 and all ji > k0. Together with (29), it gives

where δ4 can be made arbitrarily small by taking δ1, δ2 small enough and i large enough. This implies that µ

for δ1, δ2 small enough and i large enough, which contradicts m∞ = m + ǫ.

Proof of Th. 3.

The proof is similar to part (ii) of the proof of Th. 2. From Lemma 1, m∞ = m + ǫ and M∞ = M -ǫ for some ǫ ≥ 0. Suppose that ǫ > 0. When β k satisfies one of the rules [START_REF] Pronzato | A dynamical-system analysis of the optimum s-gradient algorithm[END_REF][START_REF] Pronzato | Gradient algorithms for quadratic optimization with fast convergence rates[END_REF][START_REF] Saad | Practical use of polynomial preconditionings for the conjugate gradient method[END_REF] with Fα(•) satisfying [START_REF] Tal-Ezer | Polynomial approximation of functions of matrices and applications[END_REF], then ( 19) is satisfied which implies (29) for k > k0 and some constants D ≤ 0 and 0 ≤ θ < 1, where

≤ M -ǫ. For a two-point measure supported at m and M , this implies

; in view of Th. 1, we thus have

for some constant B ≤ 0 and k > k0. Since Fα(1 -x) < 1 for any x > 0, lim sup k→∞ β k = M -ǫ and, for any δ1 > 0 and any k1, there exist some k > k1 such that β k > M -ǫ -δ1. In view of ( 29) and (30) this implies that

where δ2 can be made arbitrarily small by taking δ1 small enough and k1 large enough. This implies in turn that

for δ1 small enough and k1 large enough, which contradicts m∞ = m + ǫ with ǫ > 0.

Proof of Th. 5.

Following the same arguments as in the proof of Th. 3, (30) implies that µ

< m + ǫ for δ < δǫ and k large enough, with

From this we obtain that for small δ, 0

The conditions of Th. 1 are satisfied, so that n-1 i=2 ν k (λi) ≤ Cθ k for k > 2j0 for some constants C > 0, j0 > 0 and 0 ≤ θ < 1. Also, accounting for the fact that the distribution of β2j is not exactly symmetric in [m, M ], for any δ1 > 0 there exists some j1 such that for all j > j1, |r2jr2j+1 -R 2 2 (β2j)| < δ1. Altogether, for j > j1 large enough, R for some ǫ ′ satisfying (25), which can be written as [START_REF] Zhigljavsky | An asymptotically optimal gradient algorithm for quadratic optimization with low computational cost[END_REF] with F β (•) given by (24).