Harder, Better, Faster, Stronger Convergence Rates for Least-Squares Regression

Introduction

Many supervised machine learning problems are naturally cast as the minimization of a smooth function defined on a Euclidean space. This includes least-squares regression, logistic regression (see, e.g., Hastie et al., 2009) or generalized linear models (McCullagh and Nelder, 1989). While small problems with few or low-dimensional input features may be solved precisely by many potential optimization algorithms (e.g., Newton™ method), large-scale problems with many high-dimensional features are typically solved with simple gradient-based iterative techniques whose per-iteration cost is small.

In this paper, we consider a quadratic objective function $f$ whose gradients are only accessible through a stochastic oracle that returns the gradient at any given point plus a zero-mean finite variance random error. In this stochastic approximation framework (Robbins and Monro, 1951), it is known that two quantities dictate the behavior of various algorithms, namely the covariance matrix $V$ of the noise in the gradients, and the deviation $\theta_{0}-\theta_{\ast}$ between the initial point of the algorithm $\theta_{0}$ and any of the global minimizer $\theta_{\ast}$ of $f$. This leads to a “bias/variance”€™ decomposition (Bach and Moulines, 2013; Hsu et al., 2014) of the performance of most algorithms as the sum of two terms: (a) the bias term characterizes how fast initial conditions are forgotten and thus is increasing in a well-chosen norm of $\theta_{0}-\theta_{\ast}$; while (b) the variance term characterizes the effect of the noise in the gradients, independently of the starting point, and with a term that is increasing in the covariance of the noise.

For quadratic functions with (a) a noise covariance matrix $V$ which is proportional (with constant $\sigma^{2}$) to the Hessian of $f$ (a situation which corresponds to least-squares regression) and (b) an initial point characterized by the norm $\|\theta_{0}-\theta_{\ast}\|^{2}$, the optimal bias and variance terms are known separately. On the one hand, the optimal bias term after $n$ iterations is proportional to $\frac{L\|\theta_{0}-\theta_{\ast}\|^{2}}{n^{2}}$, where $L$ is the largest eigenvalue of the Hessian of $f$. This rate is achieved by accelerated gradient descent (Nesterov, 1983, 2004), and is known to be optimal if the number of iterations $n$ is less than the dimension $d$ of the underlying predictors, but the algorithm is not robust to random or deterministic noise in the gradients (d’Aspremont, 2008; Devolder et al., 2014). On the other hand, the optimal variance term is proportional to $\frac{\sigma^{2}d}{n}$ (Tsybakov, 2003); it is known to be achieved by averaged gradient descent (Bach and Moulines, 2013), which for the bias term only achieves $\frac{L\|\theta_{0}-\theta_{\ast}\|^{2}}{n}$ instead of $\frac{L\|\theta_{0}-\theta_{\ast}\|^{2}}{n^{2}}$.

Our first contribution in this paper is to present a novel algorithm which attains optimal rates for both the variance and the bias terms. This algorithm analyzed in Section 4 is averaged accelerated gradient descent; beyond obtaining jointly optimal rates, our result shows that averaging is beneficial for accelerated techniques and provides a provable robustness to noise.

While optimal when measuring performance in terms of the dimension $d$ and the initial distance to optimum $\|\theta_{0}-\theta_{\ast}\|^{2}$, these rates are not adapted in many situations where either $d$ is larger than the number of iterations $n$ (i.e., the number of observations for regular stochastic gradient descent) or $L\|\theta_{0}-\theta_{\ast}\|^{2}$ is much larger than $n^{2}$. Our second contribution is to provide in Section 5 an analysis of a new algorithm (based on some additional regularization) that can adapt our bounds to finer assumptions on $\theta_{0}-\theta_{\ast}$ and the Hessian of the problem, leading in particular to dimension-free quantities that can thus be extended to the Hilbert space setting (in particular for non-parametric estimation).

In order to characterize the optimality of these new bounds, our third contribution is to consider an application to non-parametric regression in Section 6 and use the known lower bounds on the statistical performance (without computational limits), which happen to match our bounds obtained from a single pass on the data and thus show optimality of our algorithm in a wide variety of particular trade-offs between bias and variance.

Our paper is organized as follows: in Section 2, we present the main problem we tackle, namely least-squares regression; then, in Section 3, we present new results for averaged stochastic gradient descent that set the stage for Section 4, where we present our main novel result leading to an accelerated algorithm which is robust to noise. Our tighter analysis of convergence rates based on finer dimension-free quantities is presented in Section 5, and their optimality for kernel-based non-parametric regression is studied in Section 6.

Least-Squares Regression

In this section, we present our least-squares regression framework, which is risk minimization with the square loss, together with the main assumptions regarding our model and our algorithms. These algorithms will rely on stochastic gradient oracles, which will come in two kinds, an additive noise which does not depend on the current iterate, which will correspond in practice to the full knowledge of the covariance matrix, and a “multiplicative/additive” noise, which corresponds to the regular stochastic gradient obtained from a single pair of observations. This second oracle is much harder to analyze.

We make the following general assumptions:

$\mathcal{H}$ is a $d$-dimensional Euclidean space with $d\geq 1$. The (temporary) restriction to finite dimension will be relaxed in Section 6.

2 Averaged Gradient Methods and Acceleration

We focus in this paper on stochastic gradient methods with and without acceleration for a quadratic function regularized by $\frac{\lambda}{2}\|\theta-\theta_{0}\|^{2}$. Stochastic gradient descent (referred to from now on as “SGD”) can be described for $n\geq 1$ as

Accelerated stochastic gradient descent is defined by an iterative system with two parameters $(\theta_{n},\nu_{n})$ satisfying for $n\geq 1$

and we note that it can be computed online as $\bar{\theta}_{n}=\frac{n}{n+1}\bar{\theta}_{n-1}+\frac{1}{n+1}\theta_{n}$.

The key ingredient in the algorithms presented above is the unbiased estimate on the gradient $f_{n}^{\prime}(\theta)$, which we now describe.

3 Stochastic Oracles on the Gradient

The first oracle is the sum of the true gradient $f^{\prime}(\theta)$ and an uncorrelated zero-mean noise that does not depend on $\theta$. Consequently it is of the form

Averaged Stochastic Gradient Descent

In this section, we provide convergence bounds for regularized averaged stochastic gradient descent. The main novelty compared to the work of Bach and Moulines (2013) is (a) the presence of regularization, which will be useful when deriving tighter convergence rates in Section 5 and (b) a simpler more direct proof. We first consider the additive noise in Section 3.1 before considering the multiplicative/additive noise in Section 3.2.

We study here the convergence of the averaged SGD recursion defined by Eq. (1) under the simple oracle from Assumption ($\mathcal{A}_{4}$). For least-squares regression, it takes the form:

This is an easy adaptation of the work of Bach and Moulines (2013, Lemma 2) for the regularized case.

The proof (see Appendix A) relies on the fact that $\theta_{n}-\theta_{*}$ is obtainable in closed form since the cost function is quadratic and thus the recursions are linear, and follows from Polyak and Juditsky (1992).

The constraint on the step-size $\gamma$ is equivalent to $\gamma(L+\lambda)\leqslant 1$ where $L$ is the largest eigenvalue of $\Sigma$ and we thus recover the usual step-size from deterministic gradient descent (Nesterov, 2004).

When $n$ tends to infinity, the algorithm converges to the minimum of $f(\theta)+\frac{\lambda}{2}\|\theta-\theta_{0}\|^{2}$ and our performance guarantee becomes $\lambda^{2}\|\Sigma^{1/2}(\Sigma+\lambda I)^{-1}(\theta_{0}-\theta_{\ast})\|^{2}$. This is the standard “bias term” from regularized ridge regression (Hsu et al., 2014) which we naturally recover here. The term \frac{\tau^{2}}{n}\mathop{\rm tr}\big{[}\Sigma^{2}(\Sigma+\lambda I)^{-2}\big{]} is usually referred to as the “variance term” (Hsu et al., 2014), and is equal to $\frac{\tau^{2}}{n}$ times the quantity \mathop{\rm tr}\big{[}\Sigma^{2}(\Sigma+\lambda I)^{-2}\big{]}, which is often called the degrees of freedom of the ridge regression problem (Gu, 2013).

For finite $n$, the first term is the usual bias term which depends on the distance from the initial point $\theta_{0}$ to the objective point $\theta_{*}$ with an appropriate norm. It includes a regularization-based component which is function of $\lambda^{2}$ and optimization-based component which depends on $(\gamma n)^{-2}$. The regularization-based bias appears because the algorithm tends to minimize the regularized function instead of the true function $f$.

Given Eq. (7), it is natural to set $\lambda\gamma=\frac{1}{n}$, and the two components of the bias term are exactly of the same order $\frac{4}{\gamma^{2}n^{2}}{\|\Sigma^{1/2}(\Sigma+\lambda I)^{-1}(\theta_{0}-\theta_{*})\|^{2}}$. It corresponds up to a constant factor to the bias term of regularized least-squares (Hsu et al., 2014), but it is achieved by an algorithm accessing only $n$ stochastic gradients. Note that here as in the rest of the paper, we only prove results in the finite horizon setting, meaning that the number of samples is known in advance and the parameters $\gamma,\lambda$ may be chosen as functions of $n$, but remain constant along the iterations (when $\lambda$ or $\gamma$ depend on $n$, our bounds only hold for the last iterate).

Note that the bias term can also be bounded by $\frac{1}{\gamma n}\|\Sigma^{1/2}(\Sigma+\lambda I)^{-1/2}(\theta_{0}-\theta_{*})\|^{2}$ when only $\|\theta_{0}-\theta_{*}\|$ is finite. See the proof in Appendix A.2 for details.

Overall we get the same performance as the empirical risk minimizer with fixed design, but with an algorithm that performs a single pass over the data.

When $\lambda=0$ we recover Lemma 2 of Bach and Moulines (2013). In this case the variance term $\frac{\tau^{2}d}{n}$ is optimal over all estimators in $\mathcal{H}$ (Tsybakov, 2008) even without computational limits, in the sense that no estimator that uses the same information can improve upon this rate.

2 Multiplicative/Additive Noise

When the general stochastic oracle in Eq. (5) is considered, the regularized LMS algorithm defined by Eq. (1) takes the form:

We have a very similar result with an additional corrective term (second line below) compared to Lemma 1.

Assume $(\mathcal{A}_{1,3}$). Consider the recursion in Eq. (8). For any regularization parameter $\lambda\leq R^{2}/2$ and for any constant step-size $\gamma\leq\frac{1}{2R^{2}}$ we have

The proof (see Appendix B) relies on a bias-variance decomposition, each term being treated separately. We adapt a proof technique from Bach and Moulines (2013) which considers the difference between the recursions in Eq. (8) and in Eq. (6).

As in Lemma 1, the bias term can also be bounded by $\frac{1}{\gamma n}\|\Sigma^{1/2}(\Sigma+\lambda I)^{-1/2}(\theta_{0}-\theta_{*})\|^{2}$ and the variance term by $\gamma\mathop{\rm tr}[\Sigma(\Sigma+\lambda I)^{-1}\xi_{n}\otimes\xi_{n}]$ (see proof in Appendices B.4 and B.5). This is useful in particular when considering unstructured noise.

The variance term is the same than in the previous case. However there is a residual term that now appears when we go to the fully stochastic oracle (second line). This term will go to zero when $\gamma$ tends to zero and can be compared to the corrective term which also appears when Hsu et al. (2014) go from fixed to random design. Nevertheless our bounds are more concise than theirs, make significantly fewer assumptions and rely on an efficient single-pass algorithm.

In this setting, the step-size may not exceed $1/(2R^{2})$, whereas with an additive noise in Lemma 1 the condition is $\gamma\leq 1/(L+\lambda)$, a quantity which can be much bigger than $1/(2R^{2})$, as $L$ is the spectral radius of ${\Sigma}$ whereas $R^{2}$ is of the order of $\mathop{\rm tr}(\Sigma)$. Note that in practice, computing $L$ is as hard as computing $\theta_{\ast}$ so that the step-size $\gamma\propto 1/R^{2}$ is a good practical choice.

For $\lambda=0$ we recover results from Défossez and Bach (2015) with a non-asymptotic bound but we lose the advantage of having an asymptotic equivalent.

Accelerated Stochastic Averaged Gradient Descent

We study the convergence under the stochastic oracle from Assumption ($\mathcal{A}_{4}$) of averaged accelerated stochastic gradient descent defined by Eq. (2) which can be rewritten for the quadratic function $f$ as a second-order iterative system with constant coefficients:

Assume ($\mathcal{A}_{4,5}$). For any constant step-size $\gamma\Sigma\preccurlyeq I$, we have for $\delta=1$,

The proof technique relies on direct moment computations in each eigensubspace obtained by O’Donoghue and Candès (2013) in the deterministic case. Indeed as $\Sigma$ is a symmetric matrix, the space can be decomposed on an orthonormal eigenbasis of $\Sigma$, and the iterations are decoupled in such an eigenbasis. Although we only provide an upper-bound, this is in fact an equality plus other exponentially small terms as shown in the proof which relies on linear algebra, with difficulties arising from the fact that this second-order system can be expressed as a linear stochastic dynamical system with non-symmetric matrices. We only provide a result for additive noise.

The first bound $\frac{1}{\gamma n^{2}}\|\theta_{0}-\theta_{*}\|^{2}$ corresponds to the usual accelerated rate. It has been shown by Nesterov (2004) to be the optimal rate of convergence for optimizing a quadratic function with a first-order method that can access only to sequences of gradients when $n\leq d$. We recover by averaging an algorithm dedicated to strongly-convex function the traditional convergence rate for non-strongly convex functions. Even if it seems surprising, the algorithm works also for $\lambda=0$ and $\delta=1$ (see also simulations in Section 7).

The second bound also matches the optimal statistical performance $\frac{\tau^{2}d}{n}$ described in the observations following Lemma 1. Accordingly this algorithm achieves joint bias/variance optimality (when measured in terms of $\tau^{2}$ and $\|\theta_{0}-\theta_{\ast}\|^{2}$).

Overall, the bias term is improved whereas the variance term is not degraded and acceleration is thus robust to noise in the gradients. Thereby, while second-order methods for optimizing quadratic functions in the singular case, such as conjugate gradient (Polyak, 1987, Section 6.1) are notoriously highly sensitive to noise, we are able to propose a version which is robust to stochastic noise.

Note that when there is no assumption on the covariance of the noise we still have the variance bounded by \frac{\gamma n}{2}\mathop{\rm tr}\big{[}\Sigma(\Sigma+\lambda I)^{-1}V\big{]}; setting $\gamma=1/n^{3/2}$ and $\lambda=0$ leads to the bound $\frac{\|\theta_{0}-\theta_{*}\|^{2}}{\sqrt{n}}+\frac{\mathop{\rm tr}V}{\sqrt{n}}$. We recover the usual rate for accelerated stochastic gradient in the non-strongly-convex case (Xiao, 2010). When the value of the bias and the variance are known, we can achieve the optimal trade-off of Lan (2012) $\frac{R^{2}\|\theta_{0}-\theta_{*}\|^{2}}{n^{2}}+\frac{\|\theta_{0}-\theta_{*}\|\sqrt{\mathop{\rm tr}V}}{\sqrt{n}}$ for \gamma=\min\Big{\{}1/R^{2},\frac{\|\theta_{0}-\theta_{*}\|}{\sqrt{\mathop{\rm tr}V}n^{3/2}}\Big{\}}.

Tighter Convergence Rates

We have seen in Corollary 1 above that the averaged accelerated gradient algorithm matches the lower bounds $\tau^{2}d/n$ and $\frac{L}{n^{2}}\|\theta_{0}-\theta_{\ast}\|^{2}$ for the prediction error. However the algorithm performs better in almost all cases except the worst-case scenarios corresponding to the lower bounds. For example the algorithm may still predict well when the dimension $d$ is much bigger than $n$. Similarly the norm of the optimal predictor $\|\theta_{*}\|^{2}$ may be huge and the prediction still good, as gradients algorithms happen to be adaptive to the difficulty of the problem. In this section, we provide such a theoretical guarantee.

The following bound stands for the averaged accelerated algorithm. It extends previously known bounds in the kernel least-mean-squares setting (Dieuleveut and Bach, 2015).

The proof is straightforward by upper bounding the terms coming from regularization, depending on $\Sigma(\Sigma+\lambda I)^{-1}$, by a power of $\lambda$ times the considered quantities. More precisely, the quantity $\mathop{\rm tr}(\Sigma(\Sigma+\lambda I)^{-1})$ can be seen as an effective dimension of the problem (Gu, 2013), and is upper bounded by $\lambda^{-b}\mathop{\rm tr}(\Sigma^{b})$ for any $b\in[0;1]$. Similarly, $\|\Sigma^{1/2}(\Sigma+\lambda I)^{-1/2}\theta_{*}\|^{2}$ can be upper bounded by $\lambda^{-r}\|\Sigma^{r/2}(\theta_{0}-\theta_{*})\|^{2}$. A detailed proof of these results is given in Appendix D.

In order to benefit from the acceleration, we choose $\lambda=(\gamma n^{2})^{-1}$. With such a choice we have the following corollary:

Assume ($\mathcal{A}_{4,5}$), for any constant step-size $\gamma(\Sigma+\lambda I)\preccurlyeq I$, we have for $\lambda=\frac{1}{\gamma(n+1)^{2}}$ and \delta\in\big{[}1-\frac{2}{n+2},1\big{]}, for the recursion in Eq. (9):

The algorithm is independent of $r$ and $b$, thus all the bounds for different values of $(r,b)$ are valid. This is a strong property of the algorithm, which is indeed adaptative to the regularity and the effective dimension of the problem (once $\gamma$ is chosen). In situations in which either $d$ is larger than $n$ or $L\|\theta_{0}-\theta_{\ast}\|^{2}$ is larger than $n^{2}$, the algorithm can still enjoy good convergence properties, by adapting to the best values of $b$ and $r$.

For $b=0$ we recover the variance term of Corollary 1, but for $b>0$ and fast decays of eigenvalues of $\Sigma$, the bound may be much smaller; note that we lose in the dependency in $n$, but typically, for large $d$, this can be advantageous.

For $r=0$ we recover the bias term of Corollary 1 and for $r=1$ (no assumption at all) the bias is bounded by $\|\Sigma^{1/2}\theta_{*}\|^{2}\leq 4R^{2}$, which is not going to zero. The smaller $r$ is, the stronger the decrease of the bias with respect to $n$ is (which is coherent with the fact that we have a stronger assumption). Moreover, $r$ is only considered between 0 and 1: indeed, if $r<0$, the constant$\|(\gamma\Sigma)^{r/2}(\theta_{0}-\theta_{*})\|$ is bigger than $\|\theta_{0}-\theta_{*}\|$, but the dependence on $n$ cannot improve beyond $(\gamma n^{2})^{-1}$. This is a classical phenomenon called “saturation” (Engl et al., 1996). It is linked with the uniform averaging scheme: here, the bias term cannot forget the initial condition faster than $n^{-2}$.

A similar result happens to hold, for averaged gradient descent, with $\lambda=(\gamma n)^{-1}$ :

where $\text{Res}(b,r,n,\gamma))$ corresponds to a residual term, which is smaller than $\mathop{\rm tr}(\Sigma^{b})n^{b}\gamma^{1+b}$ if $r\geq 0$ and does not exist otherwise. The bias term’s dependence on $n$ is degraded, thus the “saturation” limit is logically pushed down to $r=-1$, which explains the $[-1;1]$ interval for $r$. The choice $\lambda=(\gamma n)^{-1}$ arises from Th. 1, in order to balance both components of the bias term $\lambda+(\gamma n)^{-1}$. This result is proved in Appendix D.

Considering a non-uniform averaging, as proposed as after Theorem 1 the $\min_{0\leq r\leq 1}$ in Th. 3 and Corollary 2 can be extended to $\min_{-1\leq r\leq 1}$. Indeed, considering a non-uniform averaging allows to have a faster decreasing bias, pushing the saturation limit observed below.

In finite dimension these bounds for the bias and the variance cannot be said to be optimal independently in any sense we are aware of. Indeed, in finite dimension, the asymptotic rate of convergence for the bias (respectively the variance), when $n$ goes to $\infty$ is governed by $L\|\theta_{0}-\theta_{\ast}\|^{2}/n^{2}$ (resp. $\tau^{2}d/n$). However, we show in the next section that in the setting of non parametric learning in kernel spaces, these bounds lead to the optimal statistical rate of convergence among all estimators (independently of their computational cost). Moving to the infinite-dimensional setting allows to characterize the optimality of the bounds by showing that they achieve the statistical rate when optimizing the bias/variance tradeoff in Corollary 2.

Rates of Convergence for Kernel Regression

Computational convergence rates give the speed at which an objective function can decrease depending on the amount of computation which is allowed. Typically, they show how the error decreases with respect to the number of iterations, as in Theorem 1. Statistical rates, however, show how close one can get to some objective given some amount of information which is provided. Statistical rates do not depend on some chosen algorithm: these bounds do not involve computation, on the contrary, they state the best performance that no algorithm can beat, given the information, and without computational limits. In particular, any lower bound on the statistical rate implies a lower bound on the computational rates, if each iteration corresponds to access to some new information, here pairs of observations. Interestingly, many algorithms these past few years have proved to match, with minimal computations (in general one pass through the data), the statistical rate, emphasizing the importance of carrying together optimization and approximation in large scale learning, as described by Bottou and Bousquet (2008). In a similar flavor, it also appears that regularization can be accomplished through early stopping (Yao et al., 2007; Rudi et al., 2015), highlighting this interplay between computation and statistics.

To characterize the optimality of our bounds, we will show that accelerated-SGD matches the statistical lower bound in the context of non-parametric estimation. Even if it may be computationally hard or impossible to implement accelerated-SGD with additive noise in the kernel-based framework below (see remarks following Theorem 5), it leads to the optimal statistical rate for a broader class of problems than averaged-SGD, showing that for a wider set of trade-offs, acceleration is optimal.

However, in such a setting, both quantities $\|\Sigma^{r/2}\theta_{*}\|_{\mathcal{H}}$ and $\mathop{\rm tr}(\Sigma^{b})$ may exist or not. It thus arises as a natural assumption to consider the smaller $r\in[-1;1]$ and the smaller $b\in[0;1]$ such that

$\|\Sigma^{r/2}\theta_{*}\|_{\mathcal{H}}<\infty$ (meaning that $\Sigma^{r/2}\theta_{*}\in\mathcal{H}$), ($\mathcal{A}_{6}$)

$\mathop{\rm tr}(\Sigma^{b})<\infty$. ($\mathcal{A}_{7}$)

The quantities considered in Sections 2 and 5 are the natural finite-dimensional twins of these assumptions. However in infinite dimension a quantity may exist or not and it is thus an assumption to consider its existence, whereas it can only be characterized by its value, big or small, in finite dimension.

In the last decade, De Vito et al. (2005); Cucker and Smale (2002) studied non-parametric least-squares regression in the RKHS framework. These works were extended to derive rates of convergence depending on assumption $(\mathcal{A}_{6})$: Ying and Pontil (2008) studied un-regularized stochastic gradient descent and derived asymptotic rate of convergence $O(n^{-\frac{1-r}{2-r}})$, for $-1\leq r\leq 1$; Zhang (2004) studies stochastic gradient descent with averaging, deriving similar rates of convergence for $0\leq r\leq 1$; whereas Tarrès and Yao (2011) give similar performance for $-1\leq r\leq 0$. This rate is optimal without assumption on the spectrum of the covariance matrix, but comes from a worst-case analysis: we show in the next paragraphs that we can derive a tighter and optimal rate for both averaged-SGD (recovering results from Dieuleveut and Bach (2015)) and accelerated-SGD, for a larger class of kernels for the latter.

We will first describe results for averaged-SGD, then increase the validity region of these rates (which depends on $r,b$) using averaged accelerated SGD. We show that the derived rates match statistical rates for our setting and thus our algorithms reach the optimal prediction performance for certain $b$ and $r$.

We have the following result, proved in Appendix D and following from Theorem 1: for some fixed $b,r$, we choose the best step-size $\gamma$, that optimizes the bias-variance trade-off, while still satisfying the constraint $\gamma\leq 1/(2R^{2})$. We get a result for the stochastic oracle (multiplicative/additive noise).

With $\lambda=\frac{1}{\gamma n}$, we have, if $r\leq b$, under Assumptions $(\mathcal{A}_{1,3,6,7}$) and the stochastic oracle Eq. (5), for any constant step-size $\gamma\leq\frac{1}{2R^{2}}$, with $\gamma\varpropto n^{\frac{-b+r}{b+1-r}}$, for the recursion in Eq. (8):

The term $o(1)$ stands for a quantity which is decreasing to 0 when $n\rightarrow\infty$. More specifically, this constant is smaller than $3\mathop{\rm tr}(\Sigma^{b})$ divided by $n^{\chi}$, where $\chi$ is bigger than 0 (see Appendix D). The result comes from Eq. (12), with the choice of the optimal step-size.

We recover results from Dieuleveut and Bach (2015), but with a simpler analysis resulting from the consideration of the regularized version of the problem associated with a choice of $\lambda$. However, we only recover rates in the finite horizon setting.

This result shows that we get the optimal rate of convergence under Assumptions $(\mathcal{A}_{6,7})$, for $r\leq b$. This point will be discussed in more details after Theorem 5.

We now turn to the averaged accelerated SGD algorithm. We prove that it enjoys the optimal rate of convergence for a larger class of problems, but only for the additive noise which corresponds to knowing the distribution of $x_{n}$.

2 Accelerated SGD

Similarly, choosing the best step-size $\gamma$, it comes from Theorem 3, that in the RKHS setting, under additional Assumptions $(\mathcal{A}_{6,7})$, we have for the the averaged accelerated algorithm the following result:

With $\lambda=\frac{1}{\gamma n^{2}}$, we have, if $r\leq b+1/2$, under Assumptions $(\mathcal{A}_{4,5,6,7})$, for any constant step-size $\gamma\leq\frac{1}{L+\lambda}$, with $\gamma\varpropto n^{\frac{-2b+2r-1}{b+1-r}}$, for the recursion in Eq. (9):

The rate $\frac{1-r}{b+1-r}$ is always between 0 and 1, and improves when our assumptions gets stronger ($r$ getting smaller, $b$ getting smaller). Ultimately, with $b\rightarrow 0$, and $r\rightarrow-1$, we recover the finite-dimensional $n^{-1}$ rate.

We can achieve this optimal rate when $r\leq b+1/2$. Beyond, if $r>b+1/2$, the rate is only $n^{-2(1-r)}.$ Indeed, the bias term cannot decrease faster than $n^{-2(1-r)}$, as $\gamma$ is compelled to be upper bounded.

The same phenomenon appears in the un-accelerated averaged situation, as shown by Theorem 4, but the critical value was then $r\leq b$. There is thus a region (precisely $b<r\leq b+1/2$) in which only the accelerated algorithm gets the optimal rate of convergence. Note that we increase the optimality region towards optimization problems which are more ill-conditioned, naturally benefiting from acceleration.

This algorithm cannot be computed in practice (at least with computational limits). Indeed, without any further assumption on the kernel $K$, it is not possible to compute images of vectors by the covariance operator $\Sigma$ in the RKHS. However, as explained in the following remark, this is enough to show optimality of our algorithm.

These rates happen to be optimal from a statistical perspective, meaning that no algorithm which is given access to the sample points and the distribution of $x_{n}$ can perform better for all functions that satisfy assumption $(\mathcal{A}_{7})$, for a kernel satisfying $(\mathcal{A}_{6}$). Indeed it is equivalent to assuming that the function lives in some ellipsoid in the space of squared integrable functions. Note that the statistical minimization problem (and thus the lower bound) does not depend on the kernel, and is valid without computational limits. The case of learning with kernels is studied by Caponnetto and De Vito (2007) which shows these minimax convergence rates under $(\mathcal{A}_{6,7})$, under assumption that $-1\leq r\leq 0$ (but state that it can be easily extended to $0\leq r\leq 1$). They do not assume knowledge of the distribution of the inputs; however, Massart (2007) and Tsybakov (2008) discuss optimal rates on ellipsoids, and Györfi et al. (2006) proves similar results for certain class of functions under a known distribution for the input data, showing that the knowledge of the distribution does not make any difference. This minimax statistical rate stands without computational limits and is thus valid for both algorithms (additive noise that corresponds to knowing $\Sigma$, and multiplicative/additive noise). The optimal tradeoff is derived for an extended region of $b,r$ (namely $r\leq b+1/2$ instead of $r\leq b$) in the accelerated case which shows the improvement upon non-accelerated averaged SGD.

The choice of the optimal $\gamma$ is difficult in practice, as the parameters $b,r$ are unknown, and this remains an open problem (see, e.g., Birgé, 2001, for some methods for non-parametric regression).

Experiments

We illustrate now our theoretical results on synthetic examples. For $d=25$ we consider normally distributed inputs $x_{n}$ with random covariance matrix $\Sigma$ which has eigenvalues $1/i^{3}$ , for $i=1,\dots,d$, and random optimum $\theta_{*}$ and starting point $\theta_{0}$ such that $\|\theta_{0}-\theta_{*}\|=1$. The outputs $y_{n}$ are generated from a linear function with homoscedastic noise with unit signal to noise-ratio ($\sigma^{2}=1$), we take $R^{2}=\mathop{\rm tr}\Sigma$ the average radius of the data and a step-size $\gamma=1/R^{2}$ and $\lambda=0$. The additive noise oracle is used. We show results averaged over $10$ replications.

We compare the performance of averaged SGD (AvSGD), AccSGD (usual Nesterov acceleration for convex functions) and our novel averaged accelerated SGD from Section 4 (AvAccSGD, which is not the averaging of AccSGD) on two different problems: one deterministic ($\|\theta_{0}-\theta_{*}\|=1$, $\sigma^{2}=0$) which will illustrate how the bias term behaves, and one purely stochastic ($\|\theta_{0}-\theta_{*}\|=0$, $\sigma^{2}=1$) which will illustrate how the variance term behaves.

For the bias (left plot of Figure 1), AvSGD converges at speed $O(1/n)$, while AvAccSGD and AccSGD converge both at speed $O(1/n^{2})$. However, as mentioned in the observations following Corollary 1, AccSGD takes advantage of the hidden strong convexity of the quadratic function and starts converging linearly at the end. For the variance (right plot of Figure 1), AccSGD is not converging to the optimum and keeps oscillating whereas AvSGD and AvAccSGD both converge to the optimum at a speed $O(1/n)$. However AvSGD remains slightly faster in the beginning.

Note that for small $n$, or when the bias $L\|\theta_{0}-\theta_{*}\|^{2}/n^{2}$ is much bigger than the variance $\sigma^{2}d/n$, the bias may have a stronger effect, although asymptotically, the variance always dominates. It is thus essential to have an algorithm which is optimal in both regimes, what is achieved by AvAccSGD.

Conclusion

In this paper, we showed that stochastic averaged accelerated gradient descent was robust to structured noise in the gradients present in least-squares regression. Beyond being the first algorithm which is jointly optimal in terms of both bias and finite-dimensional variance, it is also adapted to finer assumptions such as fast decays of the covariance matrices or optimal predictors with large norms.

Our current analysis is performed for least-squares regression. While it could be directly extended to smooth losses through efficient online Newton methods (Bach and Moulines, 2013), an extension to all smooth or self-concordant-like functions (Bach, 2014) would widen its applicability. Moreover, our accelerated gradient analysis is performed for additive noise (i.e., for least-squares regression, with knowledge of the population covariance matrix) and it would be interesting to study the robustness of our results in the context of least-mean squares. Finally, our analysis relies on single observations per iteration and could be made finer by using mini-batches (Cotter et al., 2011; Dekel et al., 2012), which should not change the variance term but could impact the bias term.

The authors would like to thank Damien Garreau for helpful discussions.

References

Appendix A Proof of Section 3

We proof here Lemma 1 which is the extension of Lemma 2 of Bach and Moulines for the regularized case. The proof technique relies on the fact that recursions in Eq. (6) are linear since the cost function is quadratic which allows us to obtain $\theta_{n}-\theta_{*}$ in closed form.

We then have using the definition of the average

For which we will compute the two sums separately

Gathering the three terms together, we thus have

Since all the matrices in this equality are symmetric positive-definite we are allowed to bound

Unfortunately $\|\Sigma^{-1}(\theta_{0}-\theta_{*})\|$ may not be finite. However we can use that for all $u\in$ we have $\frac{1-(1-u)^{n}}{nu}\leq 1$since $\frac{1-(1-u)^{n}}{u}=\sum_{k=0}^{n}(1-u)^{k}\leq n$ and have therefore the bound

which is interesting when only $\|\theta_{0}-\theta_{*}\|$ is finite.

A.3 Proof when the noise is not structured

The bound in Eq. (17) becomes less interesting when the noise is not structured. However using the same technique we have that \big{[}I-(I-\gamma\Sigma-\gamma\lambda I)^{n-k}\big{]}^{2}(\Sigma+\lambda I)^{-1}\preccurlyeq(n-k)\gamma I and we get the following upper-bound on the variance

which is meaningful when the noise is not structured.

Appendix B Proof of Theorem 1

In this section, we will prove Theorem 1. The proof relies on a decomposition of the error as the sum of three main terms which will be studied separately. We state decomposition in Section B.1 then prove upper bounds for the different terms in Sections B.2 and B.3.

We may rewrite the regularized stochastic gradient recursion as:

be an operator from $\mathcal{H}$ to $\mathcal{H}$. We have the expansion

Our goal is to study these three terms separately and bound $\|\Sigma^{1/2}(\bar{\theta}_{n}-\theta_{\ast})\|$ for each of them.

B.2 Regularization-based bias term

This is the term: $\theta_{n}-\theta_{\ast}=\gamma\sum_{k=1}^{n}M(n,k+1)\lambda(\theta_{0}-\theta_{\ast})$, which corresponds to the recursion

initialized with $\theta_{0}=\theta_{\ast}$, and no noise.

Following the proof technique of Bach and Moulines , we are going to consider a related recursion by replacing in Equation (21) the operator $x_{n}\otimes x_{n}$ by its expectation $\Sigma$. Thus, we consider $\eta_{n}$ defined as

which satisfies the recursion (with initialization $\eta_{0}=\theta_{\ast}$) and

In order to bound $\|\Sigma^{1/2}(\theta_{n}-\theta_{*})\|$, we will independently bound $\|\Sigma^{1/2}(\eta_{n}-\theta_{*})\|$ and $\|\Sigma^{1/2}(\theta_{n}-\eta_{n})\|$ using Minkowski’s inequality.

We can now bound the recursion for $\theta_{n}-\eta_{n}$ as follows, using standard online learning proofs [Nemirovski et al., 2009]:

This leads by taking full expectations and moving terms to

Thus, if $\gamma(R^{2}+2\lambda)\leqslant\frac{1}{2}$

This leads to, summing and using initial conditions $\theta_{0}-\eta_{0}=0$, then using convexity to upper bound \big{\langle}\bar{\theta}_{n}-\bar{\eta}_{n},\Sigma(\bar{\theta}_{n}-\bar{\eta}_{n})\big{\rangle}\leq\frac{1}{n+1}\sum_{k=0}^{n}\big{\langle}\theta_{n}-\eta_{n},\Sigma(\theta_{n}-\eta_{n})\big{\rangle},

that gives the first bound on the regularization-based bias

B.3 Expansion without the regularization term

We will follow here the outline of the proof of Györfi and Walk which considers a full expansion of the function value $\|\Sigma^{1/2}(\bar{\theta}_{n}-\theta_{*})\|^{2}$. This corresponds to

Using the operator $T$ on matrices defined below, this corresponds to showing

For $\theta_{0}-\theta_{\ast}=0$, we have:

Note that when $\gamma$ tends to zero, we recover the optimal variance term.

with $\theta_{i}-\theta_{\ast}=M(i,1)(\theta_{0}-\theta_{\ast})$, that is

We follow here the proof of Défossez and Bach and consider the operator $T$ from symmetric matrices to symmetric matrices defined as

of the form $TA=(\Sigma+\lambda I)A+(\Sigma+\lambda I)A-\gamma SA$.

The operator $S$ is self-adjoint and positive. Moreover:

with $E_{0}=(\theta_{0}-\theta_{\ast})(\theta_{0}-\theta_{\ast})^{\ast}$ and $A=\Sigma(\Sigma+\lambda I)^{-1}$. This leads to

where $\langle\langle\cdot,\cdot\rangle\rangle$ denote the dot-product between self-adjoint operators.

The sum is less than its limit for $n\to\infty$, and thus, we can get rid of the term $(I-\gamma T)^{n+1}$, and we need to bound

with M:=T^{-1}\big{[}\Sigma(\Sigma+\lambda I)^{-1}\big{]}, i.e., such that

The operator $(\Sigma+\lambda I)\otimes I+I\otimes(\Sigma+\lambda I)$ is self adjoint, and so is its inverse, thus:

Moreover we can upper bound $\mathop{\rm tr}(SM):$ using Equation (24) we have

then, using Assumption ($\mathcal{A}_{1}$) :

When $\lambda=0$, without noise, we then need to bound:

with $\theta_{i}-\theta_{\ast}=M(i,1)(\theta_{0}-\theta_{\ast})$, that is

For the regularization-based bias we also have

B.5 Proof when the noise is not structured

For $\|\theta_{0}-\theta_{*}\|=0$ we have $\theta_{n}-\theta_{*}=\gamma\sum_{k=1}^{n}M(n,k+1)\varepsilon_{k}x_{k}$ which leads to

Appendix C Convergence of Accelerated Averaged Stochastic Gradient Descent

We now prove Theorem 2. We thus consider iterates satisfying Eq. (9), under Assumptions ($\mathcal{A}_{4}$), ($\mathcal{A}_{5}$). We consider a fixed step size $\gamma$ such that $\gamma(\Sigma+\lambda I)\preccurlyeq I$. Seing Eq. (9) as a linear second order for $\theta_{n}$, we will derive from exact calculations a decomposition of the errors a sum of three terms that will be studied independently. The proof is organized as follows: in Section C.1, we state the formulation as a second order linear system and derive the three main terms that have to be studied (see Lemma 2). Section C.2 studies asymptotic behaviors of the three terms, ignoring some exponentially decreasing terms, in order to give insight of how they behave. This section is not necessary for the proof, indeed a direct and exact calculation in the eigenbasis of $\Sigma$, following O’Donoghue and Candès , is provided in Section C.3. Results are summed up in Section C.4.

We study the regularized stochastic accelerated gradient descent recursion defined for $n\geq 1$ by

starting from $\theta_{0}=\nu_{0}\in\mathcal{H}$. We may rewrite it for a quadratic function $f:\theta\mapsto\frac{1}{2}\langle\theta-\theta_{*},\Sigma(\theta-\theta_{*})\rangle$ for $n\geq 2$ as

with $\theta_{0}\in\mathcal{H}$ and \theta_{1}=\big{[}I-\gamma\Sigma-\gamma\lambda I\big{]}\theta_{0}+\gamma\xi_{1}+\gamma\lambda\theta_{0}+\gamma\Sigma\theta_{\ast}.

And by centering around the optimum, we get:

Thus this is a second order iterative system which is standard to cast in a linear form

with $T=I-\gamma\Sigma-\gamma\lambda I$, $F=\begin{pmatrix}(1+\delta)T&-\delta T\\ I&0\end{pmatrix}$, $\Theta_{n}=\begin{pmatrix}\theta_{n}-\theta_{\ast}\\ \theta_{n-1}-\theta_{\ast}\end{pmatrix}$, $\Theta_{0}=\begin{pmatrix}\theta_{0}-\theta_{*}\\ \theta_{0}-\theta_{*}\end{pmatrix}$, $\Xi_{n}=\begin{pmatrix}\xi_{n}\\ 0\end{pmatrix}$ and $\Theta_{\lambda}=\begin{pmatrix}\theta_{0}-\theta_{*}\\ 0\end{pmatrix}$.

We are interested in the behavior of the average $\bar{\Theta}_{n}=\frac{1}{n+1}\sum_{k=0}^{n}\Theta_{k}$ for which we have the following general convergence result:

Error thus decomposes as the sum of three main terms:

the two first ones are bias terms, one arising from the regularization (the first one), and one arising computation (the second one),

We remark that as we have assumed that $\Sigma$ is invertible, the matrix $I-F$ can be shown to be invertible for all the considered $\delta$.

The regularization-based term will be studied directly whereas the two others will be studied in two stages. First an heuristic will lead to an asymptotic bound then an exact computation will give a non-asymptotic bound. Then using $C=H=\begin{pmatrix}\Sigma&0\\ 0&0\end{pmatrix}$ would give a convergence result on the function value and $C=\begin{pmatrix}I&0\\ 0&0\end{pmatrix}$ a result on the iterate. The end of the section is devoted to the proof of this lemma.

The sequence $\Theta_{n}$ satisfies a linear recursion, from which we get, for all $n\geq 1$:

We study the averaged sequence: $\bar{\Theta}_{n}=\frac{1}{n+1}\sum_{k=0}^{n}\Theta_{k}$ . Using the identity $\sum_{k=0}^{n-1}F^{k}=(I-F^{n})(I-F)^{-1}$, we get

and $\sum_{k=1}^{n}(I-F^{k})=\sum_{k=0}^{n}(I-F^{k})=[n+1-(I-F^{n+1})(I-F)^{-1}]$.

Using summation formulas for geometric series, we derive:

C.2 Asymptotic expansion

To give the main terms that we expect, we first provide an asymptotic analysis, which shall only be understood as an insight and is not necessary for the proof. Operator $F$ will have only eigenvalues smaller than 1, thus ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|F^{j}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ will decrease exponentially to as $j\rightarrow\infty$ (even if ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|F\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|F\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ denotes the operator norm of $F$, i.e., $\sup_{\|x\|\leq 1}\|Fx\|.$ might be bigger than 1). The asymptotic analysis relies on ignoring all terms in which $F^{j}$ appears. We thus approximately have:

where, as it has been explained $\approx$ stands for an equality up to terms that will decay exponentially. However, these terms have to be studied very carefully, what will be done in the Section C.3.

Using the matrix inversion lemma we have for $C=\begin{pmatrix}c&0\\ 0&0\end{pmatrix}$,

This gives for the regularization based term

The computation of this term is exact (not asymptotic).

And if $c$ commutes with $\Sigma$ we have the bound for $\delta\in[\frac{1-\sqrt{\gamma\lambda}}{1+\sqrt{\gamma\lambda}},1]$

And for the variance term with $V=\begin{pmatrix}v&0\\ 0&0\end{pmatrix}$, we have $C^{1/2}(I-F)^{-1}V^{1/2}=\begin{pmatrix}c^{1/2}(\gamma\Sigma+\gamma\lambda I)^{-1}v^{1/2}&0\\ 0&0\end{pmatrix}$, and

This gives the three dominant terms. However in order to control the remainders we have to compute the eigenvalues more carefully, as done in the next section.

C.3 Direct computation without the regularization based term

The computation of the regularization based term being exact we derive now direct computation for the remainders. Following O’Donoghue and Candès we consider an eigen-decomposition of the matrix $F$, in order to study independently the recursion on eigenspaces. We assume $\Sigma$ has eigenvalues $(s_{i})$ and we decompose vectors in an eigenvector basis of $\Sigma$ with $\theta_{n}^{i}=p_{i}^{\top}\theta_{n}$ and $\xi_{n}^{i}=p_{i}^{\top}\xi_{n}$ and we have the reduced equation:

Depending on $\delta$, $F_{i}$ may have two distinct complex eigenvalues of same modulus, only one (double) eigenvalue, or two real eigenvalues. We only consider the two former cases, which we detail bellow.

has discriminant $\Delta_{i}=(1-\gamma(s_{i}+\lambda))((1+\delta)^{2}(1-\gamma(s_{i}+\lambda))-4\delta)$ which is non positive as far as $\delta\in[\delta_{-};\delta_{+}]$, with $\delta_{-}=\frac{1-\sqrt{\gamma(s_{i}+\lambda)}}{1+\sqrt{\gamma(s_{i}+\lambda)}}$, $\delta_{+}=\frac{1+\sqrt{\gamma(s_{i}+\lambda)}}{1-\sqrt{\gamma(s_{i}+\lambda)}}$.

C.3.1 Two distinct eigenvalues

We first assume that $F_{i}$ has two distinct complex eigenvalues $r_{\pm}=\frac{(1+\delta)(1-\gamma(s_{i}+\lambda))\pm\sqrt{-1}\sqrt{-\Delta_{i}}}{2}$ which are conjugate. Thus the roots are of the form $\rho_{i}e^{\pm i\omega_{i}}$ with $\rho_{i}=\sqrt{\delta(1-\gamma(s_{i}+\lambda))}$, $\cos(\omega_{i})=\frac{(1+\delta)(1-\gamma(s_{i}+\lambda))}{2\rho_{i}}$, $\omega_{i}\in[-\pi/2;\pi/2]$ and $\sin(\omega_{i})=\frac{\sqrt{-\Delta_{i}}}{2\rho_{i}}$.

Let $Q_{i}=\begin{pmatrix}r_{i}^{-}&r_{i}^{+}\\ 1&1\end{pmatrix}$ be the transfer matrix into an eigenbasis of $F_{i}$, i.e., $F_{i}=Q_{i}D_{i}Q_{i}^{-1}$ with $D_{i}=\begin{pmatrix}r_{i}^{-}&0\\ 0&r_{i}^{+}\end{pmatrix}$ and $Q_{i}^{-1}=\frac{1}{r_{i}^{-}-r_{i}^{+}}\begin{pmatrix}1&-r_{i}^{+}\\ -1&r_{i}^{-}\end{pmatrix}$.

We first compute the matrix $P_{i,k}$: With

and, when developing and regrouping terms which depend on $k$, we get :

We also have $P_{i,k}=C_{i}^{1/2}Q_{i}(I-D_{i}^{k})(I-D_{i})^{-1}Q_{i}^{-1}=\sum_{j=0}^{k-1}R_{i,j}$ with

but computing error terms based in $R_{i,j}$ before summing these errors gives a looser error bound than a tight calculation using $P_{i,k}$. More precisely, if we use $P_{{i,k}}\Theta^{i}_{0}=\sum_{j=0}^{k-1}R_{i,j}\Theta^{i}_{0}$ to upper bound $\|P_{{i,k}}\Theta^{i}_{0}\|\leq\sum_{j=0}^{k-1}\|R_{i,j}\Theta^{i}_{0}\|$, we end up with a worse bound.

This can be bound with the following lemma

For all $\rho\in(0,1)$ and $\omega\in[-\pi/2;\pi/2]$ and $r^{\pm}=\rho(\cos(\omega)\pm\sqrt{-1}\sin(\omega))$ we have:

We note that the exact constant seems empirically to be $2$. This lemma is proved as Lemma 8 in Appendix E. This gives for the bias term

We also have a looser bound using $P_{{i,k}}\Theta^{i}_{0}=\sum_{j=0}^{k-1}R_{i,j}\Theta^{i}_{0}$.

As for the variance term, with $V_{i}=\begin{pmatrix}v_{i}&0\\ 0&0\end{pmatrix}$, we have \mathop{\rm tr}P_{i,k}V_{i}P_{i,k}=\Big{\|}P_{i,k}\begin{pmatrix}\sqrt{v_{i}}\\ 0\end{pmatrix}\Big{\|}^{2}.

which we can bound using the following Lemma:

For all $\rho\in(0,1)$ and $\omega\in[-\pi/2;\pi/2]$ and $r^{\pm}=\rho(\cos(\omega)\pm\sqrt{-1}\sin(\omega))$ we have:

Where we note that the exact majoration seems to be 1.3. This Lemma is proved as Lemma 10 in Appendix E.

We can also have a looser bound using $P_{{i,k}}\begin{pmatrix}v_{i}^{1/2}\\ 0\end{pmatrix}=\sum_{j=0}^{k-1}R_{i,j}\begin{pmatrix}v_{i}^{1/2}\\ 0\end{pmatrix}$ and

and \big{\|}P_{{i,k}}\begin{pmatrix}v_{i}^{1/2}\\ 0\end{pmatrix}\big{\|}\leq\frac{\sqrt{c_{i}v_{i}}(k+1)k}{2}. And this gives for the Variance term

C.3.2 One coalescent eigenvalue

We now turn to the case where $F$ has two coalescent eigenvalues, which happens when the discriminant $\Delta=0$. We assume that $F_{i}$ has one coalescent eigenvalue $r_{i}=\frac{(1+\delta)(1-\gamma(s_{i}+\lambda))}{2}$. Then, with $\delta=\frac{1-\sqrt{\gamma(s_{i}+\lambda)}}{1+\sqrt{\gamma(s_{i}+\lambda)}}$, $r_{i}=\frac{(1+\delta)(1-\gamma(s_{i}+\lambda))}{2}=1-\sqrt{\gamma(s_{i}+\lambda)}$. Then $F_{i}$ can be trigonalized as $F_{i}=Q_{i}D_{i}Q_{i}^{-1}$ with $Q_{i}=\begin{pmatrix}r_{i}&1\\ 1&0\end{pmatrix}$, $D_{i}=\begin{pmatrix}r_{i}&1\\ 0&r_{i}\end{pmatrix}$ and $Q_{i}^{-1}=\begin{pmatrix}0&1\\ 1&-r_{i}\end{pmatrix}$. We note that for all $k\geq 0$, then $D_{i}^{k}=r_{i}^{k-1}\begin{pmatrix}r_{i}&k\\ 0&r_{i}\end{pmatrix}$.

Thus with $C_{i}^{1/2}Q_{i}=\begin{pmatrix}\sqrt{c_{i}}r_{i}&\sqrt{c_{i}}\\ 0&0\end{pmatrix}$ we have

And, computing as previously the matrices products, we derive:

With $V=\begin{pmatrix}v_{i}&0\\ 0&0\end{pmatrix}$,

Alternative bounds for the bias and the variance term, as in Equations(27), (30) may be derived as well. Combining all these results, we are now able to state Theorem 2.

C.4 Conclusion

Combining results from Lemma 2, and Equations (27), (30), (31), with $c=\Sigma$, and using the following simple facts:

Under assumption $\mathcal{A}_{4}$, $\mathcal{A}_{5}$, we have $V\preccurlyeq\tau^{2}\Sigma$.

The squared norm of a vector is the sum of its squared components on the orthonormal eigenbasis. For example $\|P_{n+1}\Theta_{0}\|^{2}=\sum_{i=1}^{d}\|P_{i,n+1}\Theta_{0}^{i}\|^{2}$.

This implies, using the Equation (28) for the initial point, using $c_{i}=\sigma_{i}$ and regrouping sums as traces or norms:

which gives exactly Theorem 2 using $V\preccurlyeq\tau^{2}\Sigma$ in the Variance term, and $\lambda^{1/2}(\Sigma+\lambda I)^{-1/2}\preccurlyeq I$ in the first term.

Appendix D Tighter bounds

In this section, we chow how tighter bounds naturally appear from the regularized quantities appearing in Theorems. It only relies on simple algebraic majorations, even if one has to be careful with the allowed intervals for $r,b$.

For any $\lambda\geq 0$, for any $b\in[0;1]$, if $\mathop{\rm tr}(\Sigma^{b})$ exists, we have :

As all operators can be diagonalized in a same eigenbasis with positive eigenvalues, we have,

And the calculations are exactly the same for $\mathop{\rm tr}(\Sigma^{-2}(\Sigma+\lambda I)^{-2})\leq\frac{\mathop{\rm tr}(\Sigma^{b})}{\lambda^{b}}$. ∎

As for the bias term, we need to bound the following quantities :

For any $\lambda\geq 0$, for any $r\in[-1;1]$, we have :

For any $\lambda\geq 0$, for any $r\in[-1;0]$, we have :

For any $\lambda\geq 0$, for any $r\in[0;1]$, we have :

(No result when $r\leq 0$ because of saturation effect)

Proof relies of simple following calculations:

D.2 Theorem 3 and Equation (12)

Theorem 3 and Equation (12) are directly derived from Theorem 1 and Theorem 2, using Lemmas 5 and 6.

To derive corollaries for the optimal $\gamma$, one has to find the $\gamma$ that balances the bias and variance term and to compute the products for such a step size.

We derive from Theorem 1, when choosing $\gamma=(\lambda n)^{-1}$, and using Lemmas 5 and 6, the following bound, under assumptions of Theorem 1 :

Where $\text{Res}(n,b,r,\gamma):=3\gamma^{1+b}n^{b}\mathop{\rm tr}(\Sigma^{b})$ if $-1\leq r\leq 0$ and $\text{Res}(n,b,r,\gamma):=0$ if $0\leq r\leq 1$. When choosing the optimal $\gamma\varpropto n^{\frac{-b+r}{b+1-r}}$, we have that $\gamma^{1+b}n^{b}=n^{-1+\frac{1+b}{1+b-r}}=n^{\chi}$, with $\chi=\frac{-r}{1+b-r}\geq 0$ if $r\leq 0$. Thus the residual term is always vanishing for $r\leq 0$ and does not exist for $r\geq 0$.

D.2.2 Theorem 3

Theorem 3 directly follows from Lemmas 5 and 6 and the choice of $\gamma\varpropto n^{\frac{-2b+2r-1}{b+1-r}}$.

Appendix E Technical Lemmas

The following sequence of Lemmas appear in the proof. They are mostly independent and rely on simple calculations.

The operator \big{[}(\Sigma+\lambda I)\otimes I+I\otimes(\Sigma+\lambda I)\big{]}^{-1} is a non-decreasing operator on $(S_{n},\preccurlyeq)$

It is equivalent to show that for any symmetric positive matrix $A\in S_{n}^{+}$,

Let $\Sigma=\sum_{i\geqslant 0}\mu_{i}e_{i}\otimes e_{i}$ is the eigenvalue decomposition of $\Sigma$, then

Thus, in the orthonormal basis of eigenvectors, this is thus Hadamard product between

and the matrix C=\left(\big{(}\frac{1}{\mu_{i}+\mu_{j}+2\lambda}\big{)}_{i,j\geqslant 0}\right). Matrix $C$ is a Cauchy matrix and is thus positive. Moreover the Hadamard product of two positive matrices is positive, which concludes the proof. ∎

Remark: surprisingly, the inverse operator $(\Sigma+\lambda I)\otimes I+I\otimes(\Sigma+\lambda I)$ is not non-decreasing. Indeed, $\preccurlyeq$ is not a total order on $S_{n}$ so we may have that an operator is non-decreasing and its inverse is not.

For all $\rho\in(0,1)$ and $\omega\in[-\pi/2;\pi/2]$ and $r^{\pm}=\rho(\cos(\omega)\pm\sqrt{-1}\sin(\omega))$ we have:

We note that $\rho_{i}^{k}A_{1}$ is a real number as is is a quotient of pure complex numbers, which come from the difference between a complex and its conjugate. We first write $A_{1}$ as a combination of sine and cosine functions:

This quantity can be simplified when $\rho\rightarrow 1$ or $\omega\rightarrow 0$. We thus modify the expression of $A_{1}$ to make these dependencies clearer:

So that in that final expression all the terms behave relatively simply when $\rho\rightarrow 1$ or $\omega\rightarrow 0$. We want to upper bound:

We thus consider separately the first and second term.

And considering separately the three terms in the numerator, using numerous times that for any $a,b\in[0;1]$, $|a-b|\leq 1-ab$:

We can also change $3\rho^{k}$ into $\sqrt{5}\rho^{k}$ We have used that $|(\rho-\cos(\omega))|\leq(1-\rho\cos(\omega))$. ∎

For any $\rho_{i}\in(0;1)$, for any $\omega_{i}\in[-\pi/2;\pi/2]$

For all $\rho\in(0,1)$ and $\omega\in[-\pi/2;\pi/2]$ and $r^{\pm}=\rho(\cos(\omega)\pm\sqrt{-1}\sin(\omega))$ we have:

Once again, as the considered quantity is real, we first express it as a combination of sine and cosine functions. We then use some simple trigonometric trics to upper bound the quantity.

Let’s turn our interest to the second part of the quantity: