{'fraction_non_alphanumeric': 0.04221484610742306, 'fraction_numerical': 0.019342184671092336, 'mean_word_length': 3.947156291983878, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimized convergence of stochastic gradient descent by weighted averaging

Under mild assumptions stochastic gradient methods asymptotically achieve an optimal rate of convergence if the arithmetic mean of all iterates is returned as an approximate optimal solution. However, in the absence of stochastic noise, the arithmetic mean of all iterates converges considerably slower to the optimal solution than the iterates themselves. And also in the presence of noise, when a finite termination of the stochastic gradient method is considered, the arithmetic mean is not necessarily the best possible approximation to the unknown optimal solution. This paper aims at identifying optimal strategies in a particularly simple case, the minimization of a strongly convex function with i. i. d. noise terms and finite termination. Explicit formulas for the stochastic error and the optimization error are derived in dependence of certain parameters of the SGD method. The aim was to choose parameters such that both stochastic error and optimization error are reduced compared to arithmetic averaging. This aim could not be achieved; however, by allowing a slight increase of the stochastic error it was possible to select the parameters such that a significant reduction of the optimization error could be achieved. This reduction of the optimization error has a strong effect on the approximate solution generated by the stochastic gradient method in case that only a moderate number of iterations is used or when the initial error is large. The numerical examples confirm the theoretical results and suggest that a generalization to non-quadratic objective functions may be possible.


Introduction
In Polyak and Juditsky [14] it is shown that the stochastic gradient descent algorithm asymptotically achieves optimal complexity if short but constant step lengths are chosen and if the average over all iterates is used as final output. The aim of the current paper is to explain this phenomenon with a slightly different approach and to optimize the results further by the use of weighted averages while considering finite termination and the nature of the function to be minimized. The main source of this paper, Polyak and Juditsky [14], also inspired many other scientists to investigate weighted averages in conjunction with the stochastic gradient method in more detail. In Neu and Rosasco [13], a variant of the weighted average SGD with geometrically decreasing weights is analyzed in context of linear least-squares regression. Likewise, Cohen and Nedić [2] deal with the least-square regression. They consider constrained problems and derive upper bounds for the convergence rate and for the asymptotic ratio between convergence rate and empirical risk minimizer depending on the dimension. More abstractly, papers like Izmailov et al. [6] and Guo et al. [5], which deal with stochastic weight averaging, also build on [14]. Recently Sebbouh et al. [15] have shown almost sure convergence rates for weighted average SGD with decreasing step size. This result has already been supplemented for stronglyconvex and non-convex objective functions by Liu and Yuan [10]. There certainly is also research in this area independent of Polyak and Juditsky [14]. For example, the sampling of Needell et al. [12] and the paper of Shamir and Zhang [16], in which the polynomial-decay averaging and the suffix averaging are examined. There are numerous further modifications of the stochastic gradient approach such using momentum or heavy-ball iterations, see, e.g., [11], variance reduction [8], stochastic gradient boosting [4], and modifications tailored to specific applications. The above is merely a short and incomplete selection of related work. It seems, however, that the focus of the present paper has not been considered in this form before. This paper returns to a simple general format as considered in [14], aims at optimizing two parameters in this approach, and considers the case of an infinitely large training set in the numerical examples. A brief outline is given next.

Outline
A consequence of the results by Polyak and Juditsky [14] is that asymptotically the optimal rate of convergence of a stochastic gradient descent method is obtained when averaging all iterates with the same weight. This may seem counter-intuitive as one would expect the later iterates to be closer to the optimal solution and would therefore allow higher weights for later iterates. A short and intuitive explanation of why averaging over all iterates is optimal is given by the observation that the square root of the function value is reduced at a linear rate when sufficiently short steps with constant step lengths are chosen for minimizing a smooth, strongly convex function, while the reduction of the variance has a much slower rate of convergence. Thus, asymptotically, the stochastic effects determine the overall rate of convergence, and not the condition number of the function to be minimized. And from a stochastic point of view, averaging over all iterates with the same weight is a simple but optimal strategy. However, if only a limited number of stochastic descent steps are taken, the optimization aspect and the aspect of stochastic convergence need to be balanced to each other. The present paper is an attempt to derive a simple strategy that does so in an optimized form.
To this end, an elementary derivation of the optimality result in [14] is attempted in Section 2. for a particularly simple situation, the minimization of a strongly convex quadratic function f : R n → R, by a stochastic gradient approach. Exact and computable formulas for the variance of the averaged iterates as a measure of the stochastic error and for the contraction constant of the descent steps as a measure of the optimization error of the weighted iterates are derived. Given these explicit formulas standard nonlinear minimization algorithms with different starting points were applied to reduce both errors simultaneously by adjusting certain parameters associated with the weights and the step lengths of a stochastic gradient method.
While it turned out that the goal could not be reached of reducing both, optimization error and stochastic error at the same time, a significant reduction of the optimization error was possible by allowing a slight increase (e.g. of 10%) of the stochastic error.
Generalizations to further classes of smooth convex functions are discussed in Section 3..

Notation
The condition number with respect to the 2-norm of a square matrix A is denoted by cond(A) and the smallest and largest eigenvalues are denoted by λ min (A) and λ max (A). The Hadamard product (componentwise product) of two matrices X, Y with the same dimensions is denoted by X • Y and powers of vectors are also defined componentwise, e.g.
The diagonal of a square matrix A is denoted by diag(A) and a diagonal matrix with diagonal x ∈ R n is denoted by Diag(x). The i-th canonical unit vector is denoted by e i and the all-one-vector is denoted by e, its dimension being given by the context. Finally, let 1 n×n denote the n × n matrix with all entries equal to one.

A Standard Stochastic Descent Method
For large values of m (or when m = ∞) a stochastic gradient descent method of the following form is considered: Assume that a batch S k is chosen i. i. d. from the uniform distribution of {1, . . . , m}. Then, as is well known, the expected value of the gradient of f S k (x k ) : This motivates the stochastic gradient descent that uses the approximation ∇f S k (x k ) instead of ∇f (x k ) to define a sequence of iterates. Let γ : R + → R + be a weakly monotonously decreasing function (used in the analysis below) and consider a-priori-defined 1 step lengths γ k = γ(k) in a stochastic gradient descent approach with random noise term whose expected value is zero.

Assumptions
For simplicity it is assumed that the batch size |S k | is constant for all k. As ξ k depends on the current iterate, and thus on previous noise terms ξ i for i < k it may seem unrealistic to assume that the terms ξ k are i.i.d. as well. However, setting the noise can be written as where both ∆A S k and ∆b S k are i.i.d. and the expected values satisfy Since |S k | is constant for all k, also σ 2 2 ) are independent of k. Moreover, when the iterates x k converge to some limit x * then, asymptotically, the error terms ξ k are i.i.d.
This is the situation that is also analyzed in a different context as part (c) of Theorem 1 in [14], and this assumption will be used for simplicity below. 1 In the presence of noise, a line search is difficult.
For the analysis it can be assumed further, without loss of generality, that x * = 0 and that the Hessian of f is a diagonal matrix D with diagonal elements 0 < D 1,1 ≤ . . . ≤ D n,n . Thus, ∇f (x) = Dx.
Finally, assume that some upper bound for the largest eigenvalue D n,n is known so that the step lengths γ k can be chosen in the half-open interval (0, 1/D n,n ]. (A positive lower bound for D 1,1 is not assumed to be known.) Up to a factor of 2, the bound on γ k is essentially the bound in [14]. (Note that there is a minor error in [14]: in Assumption 2.2 (page 839) they write 2/ min i Re(D i,i ) while the correct statement would be min i 2Re(D i,i )/|D i,i | 2 in order for the argument at the bottom of page 844 of [14] to be valid.)

A Simple Algorithm with Two Parameters
In Polayk and Juditsky [14], an algorithm is considered with short step lengths and with output given by the arithmetic mean of all iterates x k generated via (2). In the absence of noise, i.e. in case that all terms ξ k are zero, it is clear that the final iterate x k is a much better approximation to the optimal solution than the average. On the other hand, as shown in [14], averaging all iterates with equal weights reduces the variance, and is asymptotically optimal for large k. In the following a method is considered that uses weighted averages with higher weights on later iterates. When the noise is small, the weighted average also is a better approximation to the optimal solution than the average with equal weights. To reduce the variance of the iterates with higher weight, a possible step length reduction for the later iterates is considered.
In the following it is assumed that an initial iterate x 0 is given and that iterates x k are generated via (2). It is further assumed that the output of the algorithm is given by weighted averagesx k defined asx and that γ k ≡ γ(k) in (2), and w j ≡ w(j) where The asymptotic analysis in [14] (Assumption 2.2) covers the case α ∈ (0, 1), β = 0 and M = 1. Larger values of α lead to a faster reduction of the step length, and for M > 1 the step length reduction in the early iterations is slower. In the numerical experiments the choice M = 1 + δk max is considered where k max is the value of k at which the iteration (5) is stopped, and δ ∈ [0, 1] is a parameter to be determined. If δ > 0 is fixed, the step lengths γ do not converge to zero when k max → ∞. It is the aim of this paper to determine optimal values α, c, M of the step lengths and an optimal value β of the weights for finite values of k max .
For completeness, the algorithm outlined above is stated in detail:
Set the weighted averagex f inal :=x k max /σ and returnx f inal as approximate minimizer of f .
For α = β = 0, Algorithm 1 reduces to an algorithm with asymptotically optimal parameters in [14]. In order to identify optimal parameters depending on the number of iterations allocated in advance an elementary analysis of the iterates is presented next.

Analysis Without Averaging
With the above assumptions, For the next step, there is new noise "−γ 1 ξ 1 " while the noise of the previous iteration is reduced, Let the product of contractions be denoted by where the empty product C i,k for i > k is equal to I by convention, for example, C k+1,k = I.
Since the largest entry of all C i is at the first diagonal position, it follows that ||C i,k || = k j=i (1 − γ j D 1,1 ). Again, the empty product is 1, i.e. C k+1,k = 1. The process (2) then leads to the aggregated representation with the expected value of x k given by E( Observe that E(x k ) coincides with the "ideal" iterates, meaning the iterates without noise and denote the "ideal" iterates by

Analysis With Averaging
Now, consider weighted averaged iterates (5) for the weight function w : R + → R + given in (6) and w j = w(j). Denote the "ideal" averages (without noise) bȳ and the "accumulated" noise in x k byξ The noise satisfies the recurrence relationξ 0 := 0 and For indices i ∈ {0, . . . , k − 1} and j ∈ {1, . . . , k} with k ∈ N a reordering of the sum leads to a representation of the "total" noisex k −x k inx k denoted byξ

Parameter Selection Without Averaging
For a comparison, first an algorithm is considered that generates the "plain" iterate x k max as output rather than the weighted averagex f inal . In this case, based on (8), the quantity to be minimized in the design of a method with optimal convergence would be a weighted sum of C 0,k−1 and of a bound of the variance of the noise term. Due to stochastic independence of the noise terms (4), the covariance matrix of x k satisfies Furthermore, by (4) the covariance matrices Cov(ξ i ) are all equal, Cov(ξ i ) ≡ Σ for some positive semidefinite Σ. In the case that Σ is a multiple of the identity matrix, one obtains an exact reformulation where the inequality follows from a repeated application of the triangle inequality and from sub-multiplicativity of the norm. Let Optimizing the expected norm of x k thus leads to the aim of choosing the step length function γ( . ) such that is minimized where µ ≥ 0 is a fixed weight. For µ := Σ 1/2 / x 0 the minimizers of (16) and (17) coincide; but unfortunately, the ratio of "noise" Σ 1/2 to "starting error" x 0 generally is not known. Nevertheless, the separation of "optimization error" (here C 0,k−1 x 0 ) and "stochastic error" (here ϑ(γ) Σ 1/2 ) can be extended to weighted averages in the next subsection and will then be exploited with the aim of identifying suitable parameters α, β, δ, M.

Parameter Selection With Averaging
Since the step lengths and weights are pre-defined the same is true for G i,k , and thus, since the noise terms ξ i are assumed to be i.i.d., also G i,k ξ i are independently distributed. The To reduce the expected value of ||ξ k || 2 , it is therefore the aim to define the weight function w and the step length function γ such that is small. As in the case of ϑ in (15), the upper bound κ 2 is the exact norm of the covariance matrix ofξ k when Σ is a multiple of the identity matrix. Simultaneously, the second goal is that also the norm ofx k should be small, i.e.
should be small. The above equation again uses the fact that all C i are diagonal matrices with the largest entry defining the norm at the (1, 1) position. Thus, in place of (17) it is the aim of choosing the step length function γ( . ) and the weight function w( . ) such that is minimized for a given fixed µ ≥ 0. Again, an "appropriate" choice of µ is not evident. For the special case (6) the evaluation of κ and τ with order k arithmetic operations is considered next.

Evaluation of All Norms
The choice of γ ∈ (0, 1/D n,n ] implies that . where the second equality again follows from the diagonal structure of C ℓ = I − γ ℓ D with the largest entry in absolute value always at the (1,1)-position.
Thus, starting from C 0,0 = 1 −c, all C 0,j can be computed for j = 1, 2, 3 . . . , k with a total of order k arithmetic operations. Likewise, starting with C k,k−1 , which is 1 by convention, the predecessors C j,k−1 can be iteratively determined by (15) can be evaluated with order k arithmetic operations.

Evaluation of All Norms
The aim of this section is to derive a scheme for evaluating all norms G i,k for 0 ≤ i ≤ k − 1 also with order k arithmetic operations, so that the parameters α, β, c, and M can easily be optimized for maximum iteration numbers k max up to the order of about 10 8 .
Note that for 0 ≤ i ≤ k − 1: where the last equality follows again since all C i are diagonal matrices with the largest entry defining the norm at the (1,1) position. Since C i+1,i = 1 it follows that In the last equation relation (22) has been used for G i+1,k in place of G i,k . Hence, For the case α = β = 0 as considered in [14], the above simplifies to G k−j,k = c(1 − (1 −c) j )/c, and the quantities κ and τ allow a closed form representation, The straightforward derivation of these formulas is omitted for brevity. For fixed values ofc > 0 it follows that κ is of the order 1/ √ k max and τ is of the order 1/k max so that convergence of Algorithm 1 to the optimal solution when k max → ∞ follows also for the case α = β = 0.
With these preparations it is now possible to evaluate τ (w, γ) and κ(w, γ) in (19) and (18) with order k arithmetic operations and thus to minimize the function r with w, γ as in (6). Here, larger values of µ are meaningful, when the noise is large compared to the distance of the initial iterate from optimality.
Note that when f is multiplied by some constant η > 0 thenc and all C i,j remain invariant, but c and all G i,k are multiplied by 1/η. Hence, when multiplying µ with η, the minimizer of (20) remains invariant.

Optimizing the Parameters of Weighted Averaging
In this sub-section the selection of the parameters α, β, c, and M in algorithm 1 is considered minimizing the function r in (20). Here, r is a weighted sum of τ and κ where τ ≥ 0 always is less than 1 and decreases with increasing values of β. On the other hand, κ is minimized for β = 0. It is the aim of the considerations below to balance the conflicting goals of minimizing both τ and κ.
To standardize the results in this sub-section the largest eigenvalue of f is fixed to D n,n = 1 throughout. Table 1 and Table 2 give some intuition about the values of τ and κ when α = β = 0 and c is also fixed to c = 1. The values of τ (0, 0) in Table 1 are not surprising, and are listed only as a comparison to Table 2 below. When k max is larger than the condition number then τ (0, 0) is roughly of the order condition number k max .
In Table 2 below it is interesting to observe that for large condition numbers such as condition number 10 4 in column 7, the bound κ(0, 0) for the variance increases first when k max increases, and starting from k max = 100 it reaches a maximum of κ(0, 0) ≈ 41.7920 for k max = 10 4.5 in row 6 before decreasing again for larger values of k max . (For the dimension, and below also for k max , square roots are rounded to integer values in Table 1 and Table 2. While condition numbers as large as 10 4 might not be typical for stochastic applications, and while moderate iteration numbers such as k max = 10 4.5 cannot render significant progress for such large condition numbers, it is interesting to consider possible improvements of κ for k max = 10 4.5 and D 1,1 = 10 −4 by optimizing τ + µ κ with respect to α, β, c, and M for different values of µ. For minimizing the function r = (τ +µ κ)/(1+µ) in (20) the descent algorithm "min f.m" from [9] was used that aims for a local minimizer near the starting point: The variables α and β were constrained to the intervals [0, 2] and [0, 5]. (The upper bounds α ≤ 2 and β ≤ 5 were chosen at will to limit the search space to a compact domain.) For c ≤ 1 a lower bound of 0.1 was chosen. (The reduction of τ is considered as too slow when c < 0.1.) Finally M was set as 1 + δk max with δ ∈ [0, 1]. In Table 3, the (approximately) optimal parameters identified with min f are listed for different weights µ > 0 and k max = 10 4.5 , D 1,1 = 10 −4 .
The numbers in the first column refer to the situation of Polyak and Juditsky [14] and are not the result of an optimization process. The two columns for µ = 0.017 show the situation that two approximate local minimizers were found with similar objective value "r" but rather different input arguments. The optimal parameters appear to be discontinuous near µ = 0.017. Since the "appropriate" weight µ > 0 depends on the unknown distance of x 0 to the optimal solution and on the unknown magnitude of the noise, the selection of α, β, c, and δ cannot be extracted from data such as Table 3, even if D 1,1 and D n,n are known.
Given an example that is not as poorly conditioned as in Table 3, namely D 1,1 = 0.03, different values of k max and weighting terms µ were considered for Table 4 such that at least two approximate local optimal solutions α, β, c, δ could be identified.  Table 4: Nearly optimal parameters α, β, c, δ for different values of µ and k max .
Because of the discontinuous dependence of α and β (and also of τ and κ) on µ as observed in Table 3 and Table 4, another selection process for optimizing the parameters α, β, c, δ was considered, namely Problem (23) uses a compact notation highlighting the changes compared to (20) in Table 3 and Table 4. As in Table 3 and Table 4, w depends on β ∈ [0, 5] and γ depends on α ∈ [0, 2], c ∈ [0.1, 1], and δ ∈ [0, 1], and w 0 , γ 0 refer to α = β = δ = 0, and c = 1. The restriction v 1 , v 2 ≤ 0.1 implies that neither τ nor κ are allowed to be more than 10% worse than the choice w 0 , γ 0 . Again, µ > 0 is a weight balancing optimization error and stochastic error; however, due to the upper bound on v 1 , v 2 the dependence on µ of the optimal solution turns out to be less pronounced.
Using "min fc.m" from [7] with different starting points, approximate local optimal solutions were computed for Problem (23). For upper bounds v 1 , v 2 ≤ 0 aiming at a simultaneous reduction of both τ and κ, no point apart from α = β = 0 could be identified via "min fc.m". However, allowing a small slack of v 1 , v 2 ≤ 0.1 as in (23)  For large initial errors x 0 this reduction of τ will outweigh the 10%-increase of κ. In this situation a significant gain in the optimization error could be achieved when allowing a small increase of the stochastic error.
To test the robustness of this solution, the other parameters k max and D 1,1 of Table 1 and Table 2 were evaluated as well for α = δ = 0, β = 0.7116: For higher condition numbers up to 10 4 the increase of κ is less than 21% and it is less than 10% for smaller condition numbers or larger values of k max . Likewise, the reduction of τ deteriorates for increasing condition numbers but improves for smaller condition numbers or for larger values of k max .

Practical Parameters
Summarizing, while a step length reduction generally does not lead to a significant improvement of convergence, moderately growing weights such as w j = j 0.7 do lead to faster reduction of the optimization error without deteriorating the stochastic error. Thus, the parameter setting (24) is used in the numerical examples in Subsection 3.2 with the slight modification that β is set to β = 0.7 (for simplicity).

Test Examples
Algorithm 1 is rather simple and the results in Section 2. are not very strong but "robust" in the sense that they are independent of the dimension n and of the number m in the definition (1) of f .

Simple Examples
For the first set of examples a situation is considered that does satisfy the strong assumptions of Section 2. and where m = ∞ so that variance reduction by periodic full gradient evaluations is not possible.
Here, the convergence is rather slow with a growing number of iterations -it is of the order 1/ √ k max , i.e., as in the analysis of Polyak and Juditsky [14], the stochastic effects do dominate the convergence. The reductions of the initial error observed in Table 6 are better by a factor of 2.5 or 2.6 than the theoretical bounds κ(0, 0) listed in the second column of Table 2.
For all test runs listed in Table 5 and Table 6, the final iterate had a norm x k max ∈ [1.21, 1.39] -since α = 0, the step length was constant and the final iterate was largely determined by the noise.
For a large initial error λ there is a significant gain when choosing β = 0.7 compared to β = 0 while there is not much loss of the choice β = 0.7 for small initial errors. (For λ = 10 8 the final iterate x 100000 has norm 1.26, i.e. it is much closer to the optimal solution than the average x f inal 2 with β = 0.) For large initial errors, the value of τ will dominate the convergence behavior. Indeed, for the entry with β = 0.7 in the last column of Table 7 a reduction of the initial error by 10 8 /6.1 ≈ 1.6 · 10 7 can be observed which is larger by a (moderate) factor of about 3.5 than the theoretical bound 1/τ (0, 0.7) ≈ 4.5 · 10 6 and close to bound given by the last entry of Column 2 in Table 1 with 1/τ (0, 0) = 10 7.0458 ≈ 1.1 · 10 7 for β = 0.
Summarizing, the above simple examples confirm the independence of the theoretical results with respect to the dimension n and the number m of functions in the definition of f , as long as the variance ∇f (x) − ∇f k (x) 2 is bounded independent of x and n. The observed results also indicate that the theoretical results in Table 1 and Table 2 are not overly pessimistic.

Relaxing the Assumptions
To test the limits of Algorithm 1 the next example concerns a somewhat more difficult case where again m = ∞, and where each single function f k carries rather little information about the function f to be minimized.
For this set of examples a non-singular randomly generated fixed matrix A ∈ R n×n is chosen and a scaling factor ρ := 1/ √ n. For k = 1, 2, . . . random vectors r k and b k are drawn independently from n-variate normal distributions with expected value 0 n and covariance matrices I n respectively ρ 2 I n , i.e. r k ∼ N n (0 n , I n ) and b k ∼ N n (0 n , ρ 2 I n ). The choice of ρ implies that E( b k 2 2 ) = 1. Define a k := Ar k (so that a k ∼ N n (0 n , AA T )) and for given where the last equation uses the linearity of the expected value for fixed x. Since And since the fourth momenta of r k and b k exist and f k is a quadratic function of r k and b k it follows that f k has bounded variance and exists almost surely. Moreover, ∇f (x) = AA T x, which coincides with the expected value of ∇f k (x): The noise defined by and the covariance matrix is given by To keep the presentation self-contained a short proof of (25) is given in the appendix. If x → 0, the noise terms ξ k indeed are i.i.d. in the limit, as assumed in (4), but for larger starting errors x 0 this assumption is violated, more so when the dimension n grows large.
Note that (since the matrices on the right hand side of (25) all are positive semidefinite) where AA T xx T AA T 1/2 = ∇f (x) . Hence, (25) implies that standard deviation of the noise ξ k (that is added to the gradient at each iteration of Algorithm 1) always dominates the norm of the gradient itself.
Similar to the analysis in Section 2., again, for the numerical experiments it is sufficient to use a diagonal matrix D in place of A. Indeed, consider the singular value decomposition A = UDV T with a diagonal matrix D containing the singular values of A and orthogonal matrices U, V . Since a k = Ar k = UDV T r k ∼ N n (0 n , UDDU T ) one obtains with the transformation z := U T x. Likewise, (25) translates to an equivalent formula for z andξ k := U T ξ k . The examples in Table 8 and Table 9 refer to n = 100, randomly generated D (uniform distribution scaled to D n,n = n −1/2 and D 1,1 = (10n) −1/2 ) so that the condition number of D 2 is 10. The scaling by n −1/2 was chosen to compensate for the norm r k which is of the order n 1/2 . For the choice k max = 1000 and β = 0.7 (for comparison also β = 0) the following average final errors were obtained:  Table 9: Final error in dependence of k max when x 0 2 = 10 3 .
Above, it takes about 10 times longer for the algorithm with β = 0 to reach an error of 0.0048 or 0.0049 than it takes with the choice β = 0.7. When multiplying D with a positive constant greater than one, convergence actually improves (as the noise ratio b k / a k decreases), but overall, this is an example where the stochastic effects dominate. This also implies that it may be difficult to improve over the simple scheme of Algorithm 1 for this type of example.

Non-Quadratic Test Examples
The intent of this paper of course is to motivate a rule that is more generally applicable for stochastic gradient descent approaches. This motivation is supported by the observation that if the functions f i are convex and twice continuously differentiable, then locally, the functions can be closely approximated by quadratic functions for which the analysis holds.
The MNIST database [3] provides 70000 handwritten digit images with corresponding labels 0,1,. . . ,9. Each image consists of n = 28 · 28 = 784 pixels with entries between 0 and 1. To test the performance of Algorithm 1, m = 60000 of the labeled images are used as training data to find a rule that predicts whether the digits of the remaining 10000 test images are 0 or not. This prediction can be compared with the labels of the test images to determine the false classification rate "FCR" i.e. the percentage of incorrectly classified images.
For the numerical experiments the data for the i-th image is put into the following form: the vector a i ∈ [0, 1] 784 contains the information about the image itself and b i ∈ {−1, +1} is the corresponding label, where "+1" means that the image shows the digit 0 and "-1" means this is not the case.
Let σ : R → ]0, 1[ , σ(t) := 1/(1 + e −t ) be the logistic function with derivative σ ′ (t) = σ(t)(1 − σ(t)) and consider f : R n → R defined by Note that t → log(σ(t)) is a smooth (and asymptotically exact) approximation of the function t → max{t, 0}. Thus minimizing f , approximates the maximization of the terms This motivates the following classification rule: Let an approximate minimizerx f inal of f be given and a test image a new . Then a new is classified to represent the digit "0" if a T newx f inal > 0, and a new is classified not to represent the digit "0" otherwise.
The derivatives of f are given by , Here, f is convex since the Hessian matrix is positive semidefinite, i. e. ∇ 2 f (x) 0.
To estimate a value for the maximum step length c in Algorithm 1, consider the case that a i is a random vector whose components are continuously uniformly distributed on ]0, 1[ . Thus, every component of a i has expected value 1/2 and second momentum 1/3 and the components of a i are independent of each other. Since the Hessian matrix can be estimated by  Table 10.
For this example it was also possible to apply Newton's method with line search which generated an approximation x opt with ∇f (x opt ) 2 ≈ 2.3 · 10 −10 and an associated false classification rate of FCR≈ 0.80%. The Hessian of f at x opt was numerically singular and the spectrum of ∇ 2 f (x opt ) was quite dense near zero and did not allow a clear identification of the null space. It was not possible to give a reliable estimate for the condition number of the Hessian of f , even when restricted to the range space of ∇ 2 f (x opt ). The Hessian certainly was very far from being well-conditioned. For all approximate solutionsx f inal generated in Table 10, the distance x f inal − x opt was about the same, namely close to 1100. In particular, a convergence of the iterates of Algorithm 1 to a minimizer of f could not be observed. Nevertheless a small gain in the rate of convergence with β = 0.7 compared to β = 0 could be observed for smaller values of k max . For larger values of k max the asymptotic optimality of β = 0 in [14] can be confirmed in the sense that there is not much difference of β = 0.7 and β = 0.
It is stressed that the intention of this example was not to propose a new classification scheme for MNIST -other classification schemes are certainly better -but to test Algorithm expected value returns a sum of 9 terms that are considered one by one: The expectation of the constant term AA T xx T AA T , of course, is the term itself. The expectations of −AA T xb T and −bx T AA T are both zero (since AA T x is a fixed vector). The expectation of bb T is ρ 2 I. The expectation of the two terms −AA T xx T aa T and −aa T xx T AA T is −AA T xx T AA T each. The expectation of aa T xx T aa T is given by Summing up all 9 expectations above leads to (25).

Conclusion
The optimization error and the stochastic error are analyzed for the simple case of a stochastic gradient method applied to a strongly convex quadratic function with stochastic gradients that satisfy the assumption of being i.i.d. By optimizing both errors for an algorithm with finite termination it was possible to modify the known asymptotically optimal parameter selection of a stochastic gradient method. As predicted by the analysis in [14], for large numbers of iterations k max the results cannot be improved, but for moderate numbers of iterations the numerical experiments confirm a gain in accuracy. This gain can be achieved by a simple parameter selection and without additional computational cost.
The extension to smooth non-convex functions or to limited memory Quasi-Newton approaches such as presented in [1] are the subject of future research.