Momentum and Stochastic Momentum for Stochastic Gradient, Newton, Proximal Point and Subspace Descent Methods

Introduction

Two of the most popular algorithmic ideas for solving optimization problems involving big volumes of data are stochastic approximation and momentum. By stochastic approximation we refer to the practice pioneered by Robins and Monro of replacement of costly-to-compute quantities (e.g., gradient of the objective function) by cheaply-to-compute stochastic approximations thereof (e.g., unbiased estimate of the gradient). By momentum we refer to the heavy ball technique originally developed by Polyak to accelerate the convergence rate of gradient-type methods.

While much is known about the effects of stochastic approximation and momentum in isolation, surprisingly little is known about the combined effect of these two popular algorithmic techniques. For instance, to the best of our knowledge, there is no context in which a method combining stochastic approximation with momentum is known to have a linear convergence rate. One of the contributions of this work is to show that there are important problem classes for which a linear rate can indeed be established for a range of stepsize and momentum parameters.

In this paper we study three closely related problems:

(bounded) concave quadratic maximization.

These problems and the relationships between them are described in detail in Section 3. Here we only briefly outline some of the key relationships. By stochastic optimization we refer to the problem of the form

2 Three stochastic methods

Method 1: stochastic gradient descent (SGD),

Method 2: stochastic Newton method (SN), and

Method 3: stochastic proximal point method (SPP);

all with a fixed stepsize $\omega>0$. The methods will be described in detail in Section 3; see also Table 2 for a quick summary.

The equivalence of these methods is useful for the purposes of this paper as it allows us to study their variants with momentum by studying a single algorithm only. We are not aware of any successful attempts to analyze momentum variants of SN and SPP and as we said before, there are no linearly convergent variants of SGD with momentum in any setting.

3 Best approximation, duality and stochastic dual subspace ascent

It was shown in in the $\omega=1$ case and in in the general $\omega>0$ case that SGD, SN and SPP converge to a very particular minimizer of $f$: the projection of the starting point $x_{0}$ onto the solution set of the linear system (2). This naturally leads to the best approximation problem, which is the problem of projectingIn the rest of the paper we consider projection with respect to an arbitrary Euclidean norm. a given vector onto the solution space of the linear system (2):

The dual of the best approximation problem is an unconstrained concave quadratic maximization problem . Consistency of ${\bf A}x=b$ implies that the dual is bounded. It follows from the results of that for $\omega=1$, the random iterates of SGD, SN and SPP arise as affine images of the random iterates produced by an algorithm for solving the dual of the best approximation problem (3), known as

Method 4: stochastic dual subspace ascent (SDSA).

In this paper we show that this equivalence extends beyond the $\omega=1$ case, specifically for $0<\omega<2$, and further study SDSA with momentum. We then show that SGD, SN and SPP with momentum arise as affine images of SDSA with momentum. SDSA proceeds by taking steps in a random subspace spanned by the columns of ${\bf S}$ randomly drawn in each iteration from ${\cal D}$. In this subspace, the method moves to the point which maximizes the dual objective, $D(y)$. Since ${\cal D}$ is an arbitrary distribution of random matrices, SDSA moves in arbitrary random subspaces, and as such, can be seen as a vast generalization of randomized coordinate descent methods and their minibatch variants .

4 Structure of the paper

The remainder of this work is organized as follows. In Section 2 we summarize our contributions in the context of existing literature. In Section 3 we provide a detailed account of the stochastic optimization problem, the best approximation and its dual. Here we also describe the SGD, SN and SPP methods. In Section 4 we describe and analyze primal methods with momentum (mSGD, mSN and mSPP), and in Section 5 we describe and analyze the dual method with momentum (mSDSA). In Section 6 we describe and analyze primal methods with stochastic momentum (smSGD, smSN and smSPP). Numerical experiments are presented in Section 8. Proofs of all key results can be found in the appendix.

5 Notation

Momentum Methods and Our Contributions

In this section we give a brief review of the relevant literature, and provide a summary of our contributions.

The baseline first-order method for minimizing a differentiable function $f$ is the gradient descent (GD) method,

where $\omega_{k}>0$ is a stepsize. For convex functions with $L$-Lipschitz gradient (function class ${\cal F}^{1,1}_{0,L}$), GD converges at at the rate of ${\cal O}(L/\epsilon)$. When, in addition, $f$ is $\mu$-strongly convex (function class ${\cal F}^{1,1}_{\mu,L}$), the rate is linear: ${\cal O}((L/\mu)\log(1/\epsilon))$ . To improve the convergence behavior of the method, Polyak proposed to modify GD by the introduction of a (heavy ball) momentum termArguably a much more popular, certainly theoretically much better understood alternative to Polyak’s momentum is the momentum introduced by Nesterov , leading to the famous accelerated gradient descent (AGD) method. This method converges nonassymptotically and globally; with optimal sublinear rate ${\cal O}(\sqrt{L/\epsilon})$ when applied to minimizing a smooth convex objective function (class ${\cal F}^{1,1}_{0,L}$), and with the optimal linear rate ${\cal O}(\sqrt{L/\mu}\log(1/\epsilon))$ when minimizing smooth strongly convex functions (class ${\cal F}^{1,1}_{\mu,L}$). Both Nesterov’s and Polyak’s update rules are known in the literature as “momentum” methods. In this paper, however, we focus exclusively on Polyak’s heavy ball momentum., $\beta(x_{k}-x_{k-1})$. This leads to the gradient descent method with momentum (mGD), popularly known as the heavy ball method:

More specifically, Polyak proved that with the correct choice of the stepsize parameters $\omega_{k}$ and momentum parameter $\beta$, a local accelerated linear convergence rate of ${\cal O}(\sqrt{L/\mu}\log(1/\epsilon))$ can be achieved in the case of twice continuously differentiable, $\mu$-strongly convex objective functions with $L$-Lipschitz gradient (function class $\mathcal{F}_{\mu,L}^{2,1}$). See the first line of Table 1.

Recently, Ghadimi et al. performed a global convergence analysis for the heavy ball method. In particular, the authors showed that for a certain combination of the stepsize and momentum parameter, the method converges sublinearly to the optimum when the objective function is convex and has Lipschitz gradient ($f\in\mathcal{F}_{0,L}^{1,1}$), and linearly when the function is also strongly convex ($f\in\mathcal{F}_{\mu,L}^{1,1}$). A particular, selection of the parameters $\omega$ and $\beta$ that gives the desired accelerated linear rate was not provided.

To the best of our knowledge, despite considerable amount of work on the on heavy ball method, there is still no global convergence analysis which would guarantee an accelerated linear rate for $f\in\mathcal{F}_{\mu,L}^{1,1}$. However, in the special case of a strongly convex quadratic, an elegant proof was recently proposed in . Using the notion of integral quadratic constraints from robust control theory, the authors proved that by choosing $\omega_{k}=\omega=4/(\sqrt{L}+\sqrt{\mu})^{2}$ and $\beta=(\sqrt{L/\mu}-1)^{2}/(\sqrt{L/\mu}+1)^{2}$, the heavy ball method enjoys a global asymptotic accelerated convergence rate of ${\cal O}(\sqrt{L/\mu}\log(1/\epsilon))$. The aforementioned results are summarized in the first part of Table 1.

Extensions of the heavy ball method have been recently proposed in the proximal setting , non-convex setting and for distributed optimization .

2 Stochastic heavy ball method

In contrast to the recent advances in our theoretical understanding of the (classical) heavy ball method, there has been less progress in understanding the convergence behavior of stochastic variants of the heavy ball method. The key method in this category is stochastic gradient descent with momentum (mSGD; aka: stochastic heavy ball method):

where $g_{k}$ is an unbiased estimator of the true gradient $\nabla f(x_{k})$. While mSGD is used extensively in practice, especially in deep learning , its convergence behavior is not very well understood.

In fact, we are aware of only two papers, both recent, which set out to study the complexity of mSGD: the work of Yang et al. , and the work of Gadat et al. . In the former paper, a unified convergence analysis for stochastic gradient methods with momentum (heavy ball and Nesterov’s momentum) was proposed; and an analysis for both convex and non convex functions was performed. For a general Lipschitz continuous convex objective function with bounded variance, a rate of $O(1/\sqrt{\epsilon})$ was proved. For this, the authors employed a decreasing stepsize strategy: $\omega_{k}=\omega_{0}/\sqrt{k+1}$, where $\omega_{0}$ is a positive constant. In , the authors first describe several almost sure convergence results in the case of general non-convex coercive functions, and then provided a complexity analysis for the case of quadratic strongly convex function. However, the established rate is slow. More precisely, for strongly convex quadratic and coercive functions, mSGD with diminishing stepsizes $\omega_{k}=\omega_{0}/k^{\beta}$ was shown to convergence as ${\cal O}(1/k^{\beta})$ when the momentum parameter is $\beta<1$, and with the rate $O(1/\log k)$ when $\beta=1$. The convergence rates established in both of these papers are sublinear. In particular, no insight is provided into whether the inclusion of the momentum term provides what is was aimed to provide: acceleration.

The above results are summarized in the second part of Table 1. From this perspective, our contribution lies in providing an in-depth analysis of mSGD (and, additionally, of SGD with stochastic momentum). Our contributions are discussed next.

3 Connection to incremental gradient methods

Assuming ${\cal D}$ is discrete distribution (i.e., we sample from $M$ matrices, ${\bf S}^{1},\dots,{\bf S}^{M}$, where ${\bf S}^{i}$ is chosen with probability $p_{i}>0$), we can write the stochastic optimization problem (1) in the finite-sum form

Choosing $x_{0}=x_{1}$, mSGD with fixed stepsize $\omega_{k}=\omega$ applied to (4) can be written in the form

where ${\bf S}_{t}={\bf S}^{i}$ with probability $p_{i}$. Problem (4) can be also solved using incremental average/aggregate gradient methods, such as the IAG method of Blatt et al. . These methods have a similar form to (5), with the main difference being in the way the past gradients are aggregated. While (5) uses a geometric weighting of the gradients, the incremental average gradient methods use a uniform/arithmetic weighting. The stochastic average gradient (SAG) method of Schmidt et al. can be also written in a similar form. Note that mSGD uses a geometric weighting of previous gradients, while the the incremental and stochastic average gradient methods use an arithmetic weighting. Incremental and incremental average gradient methods are widely studied algorithms for minimizing objective functions which can expressed as a sum of finite convex functions. For a review of key works on incremental methods and a detailed presentation of the connections with stochastic gradient descent, we refer the interested reader to the excellent survey of Bertsekas ; see also the work of Tseng .

In , an incremental average gradient method with momentum was proposed for minimizing strongly convex functions. It was proved that the method converges to the optimum with linear rate. The rate is always worse than that of the no-momentum variant. However, it was shown experimentally that in practice the method is faster, especially in problems with high condition number. In our setting, the objective function has a very specifc structure (1). It is not a finite sum problem as the distribution ${\cal D}$ could be continous; and we also do not assume strong convexity. Thus, the convergence analysis of can not be directly applied to our problem.

4 Summary of contributions

We now summarize the contributions of this paper.

New momentum methods. We study several classes of stochastic optimization algorithms (SGD, SN, SPP and SDSA) with momentum, which we call mSGD, mSN, mSPP and mSDSA, respectively (see the first and second columns of Table 2). We do this in a simplified setting with quadratic objectives where all of these algorithms are equivalent. These methods can be seen as solving three related optimization problems: the stochastic optimization problem (1), the best approximation problem (3) and its dual. To the best of our knowledge, momentum variants of SN, SPP and SDSA were not analyzed before.

Sublinear rate for Cesaro averages. We show that the Cesaro averages, $\hat{x}_{k}=\frac{1}{k}\sum_{t=0}^{k-1}x_{t}$, of all primal momentum methods enjoy a sublinear $O(1/k)$ rate (see line 3 of Table 3). This holds under weaker assumptions than those which lead to the linear convergence rate.

Primal-dual correspondence. We show that SGD, SN and SPP with momentum arise as affine images of SDSA with momentum (see Theorem 5). This extends the result of where this was shown for the no-momentum methods ($\beta=0$) and in the special case of the unit stepsize ($\omega=1$).

Stochastic momentum. We propose a new momentum strategy, which we call stochastic momentum. Stochastic momentum is a stochastic (coordinate-wise) approximation of the deterministic momentum, and hence is much less costly, which in some situations leads to computational savings in each iteration. On the other hand, the additional noise introduced this way increases the number of iterations needed for convergence. We analyze the SGD, SN and SPP methods with stochastic momentum, and prove linear convergence rates. We prove that in some settings the overall complexity of SGD with stochastic momentum is better than the overall complexity of SGD with momentum. For instance, this is the case if we consider the randomized Kaczmarz (RK) method as a special case of SGD, and if ${\bf A}$ is sparse.

Space for generalizations. We hope that the present work can serve as a starting point for the development of SN, SPP and SDSA methods with momentum for more general classes (beyond special quadratics) of convex and perhaps also nonconvex optimization problems. In such more general settings, however, the symmetry which implies equivalence of these algorithms will break, and hence a different analysis will be needed for each method.

5 No need for variance reduction

SGD is arguably one of the most popular algorithms in machine learning. Unfortunately, SGD suffers from slow convergence, which is due to the fact that the variance of the stochastic gradient as an estimator of the gradient does not naturally diminish. For this reason, SGD is typically used with a decreasing stepsize rule, which ensures that the variance converges to zero. However, this has an adverse effect on the convergence rate. For instance, SGD has a sublinear rate even if the function to be minimized is strongly convex. To overcome this problem, a new class of so-called variance-reduced methods was developed over the last 2-5 years, including SAG , SDCA , SVRG/S2GD , minibatch SVRG/S2GD , and SAGA .

Since we assume that the linear system (2) is feasible, it follows that the stochastic gradient vanishes at the optimal point (i.e., $\nabla f_{\mathbf{S}}(x_{*})=0$ for any ${\bf S}$). This suggests that additional variance reduction techniques are not necessary since the variance of the stochastic gradient drops to zero as we approach the optimal point $x_{*}$. In particular, in our context, SGD with fixed stepsize enjoys linear rate without any variance reduction strategy . Hence, in this paper we can bypass the development of variance reduction techniques, which allows us to focus on the momentum term.

Technical Preliminaries

A general framework for studying consistent linear systems via carefully designed stochastic reformulations was recently proposed by Richtárik and Takáč . In particular, given the consistent linear system (2), they provide four reformulations in the form of a stochastic optimization problem, stochastic linear system, stochastic fixed point problem and a stochastic intersection problem. These reformulations are equivalent in the sense that their solutions sets are identical. That is, the set of minimizers of the stochastic optimization problem is equal to the set of solutions of the stochastic linear system and so on. Under a certain assumption, for which the term exactness was coined in , the solution sets of these reformulations are equal to the solution set of the linear system.

where ${\bf H}$ is a random symmetric positive semidefinite matrix defined as ${\bf H}:={\bf S}({\bf S}^{\top}{\bf A}{\bf B}^{-1}{\bf A}^{\top}{\bf S})^{\dagger}{\bf S}^{\top}.$ By $\dagger$ we denote the Moore-Penrose pseudoinverse.

have the same spectrum. Matrix ${\bf B}$ is symmetric and positive semidefinite (with respect to the standard inner product). Let

be the eigenvalue decomposition of ${\bf W}$, where ${\bf U}=[u_{1},\dots,u_{n}]$ is an orthonormal matrix of eigenvectors, and $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{n}$ are the corresponding eigenvalues. Let $\lambda_{\min}^{+}$ be the smallest nonzero eigenvalue, and $\lambda_{\max}=\lambda_{n}$ be the largest eigenvalue. It was shown in that $0\leq\lambda_{i}\leq 1$ for all $i\in[n]$.

2 Three algorithms for solving the stochastic optimization problem

The authors of consider solving the stochastic optimization problem (1) via stochastic gradient descent (SGD)The gradient is computed with respect to the inner product $\langle{\bf B}x,y\rangle$.

where $\omega>0$ is a fixed stepsize and ${\bf S}_{k}$ is sampled afresh in each iteration from ${\cal D}$. Note that the gradient of $f_{\bf S}$ with respect to the ${\bf B}$ inner product is equal to

where ${\bf Z}:={\bf A}^{\top}{\bf H}{\bf A}$, and $x_{*}$ is any vector in ${\cal L}$.

They observe that, surprisingly, SGD is in this setting equivalent to several other methods; in particular, to the stochastic Newton methodIn this method we take the ${\bf B}$-pseudoinverse of the Hessian of $f_{{\bf S}_{k}}$ instead of the classical inverse, as the inverse does not exist. When ${\bf B}={\bf I}$, the ${\bf B}$ pseudoinverse specializes to the standard Moore-Penrose pseudoinverse.,

and to the stochastic proximal point methodIn this case, the equivalence only works for $0<\omega\leq 1$.

3 Stochastic fixed point problem

The stochastic fixed point problem considered in as one of the four stochastic reformulations has the form

a formula for ${\cal L}_{\bf S}$ is obtained by replacing ${\bf A}$ with ${\bf S}^{\top}{\bf A}$ everywhere.

The stochastic fixed point method (with relaxation parameter $\omega>0$) for solving (13) is defined by

4 Best approximation problem, its dual and SDSA

It was shown in that the above methods converge linearly to $x_{*}=\Pi^{{\bf B}}_{{\cal L}}(x_{0})$; the projection of the initial iterate onto the solution set of the linear system. Hence, besides solving problem (1), they solve the best approximation problem

The Fenchel dual of (16) is the (bounded) unconstrained concave quadratic maximization problem

Boundedness follows from consistency. It turns out that by varying ${\bf A},{\bf B}$ and $b$ (but keeping consistency of the linear system), the dual problem in fact captures all bounded unconstrained concave quadratic maximization problems.

The dual method—stochastic dual subspace ascent (SDSA)—has the form

where ${\bf S}_{k}$ is in each iteration sampled from ${\cal D}$, and $\lambda_{k}$ is chosen greedily, maximizing the dual objective $D$: $\lambda_{k}\in\arg\max_{\lambda}D(y_{k}+{\bf S}_{k}\lambda)$. Such a $\lambda$ might not be unique, however. SDSA is defined by picking the solution with the smallest (standard Euclidean) norm. This leads to the formula:

SDSA proceeds by moving in random subspaces spanned by the random columns of ${\bf S}_{k}$. In the special case when $\omega=1$ and $y_{0}=0$, Gower and Richtárik established the following relationship between the iterates $\{x_{k}\}$ produced by the primal methods (9), (11), (12), (15) (which are equivalent), and the dual method (19):

5 Other related work

Variants of the sketch-and-project methods have been recently proposed for solving several other problems. Xiang and Zhang show that the sketch-and-project framework is capable of expressing, as special cases, randomized variants of 16 classical algorithms for solving linear systems. Gower and Richtárik use similar ideas to develop of linearly convergent randomized iterative methods for computing/estimating the inverse and the pseudoinverse of a large matrix, respectively. A limited memory variant of the stochastic block BFGS method for solving the empirical risk minimization problem arising in machine learning was proposed by Gower et al. . Tu et al. utilize the sketch-and-project framework to show that breaking block locality can accelerate block Gauss-Seidel methods. In addition, they develop an accelerated variant of the method for a specific distribution ${\cal D}$. Loizou and Richtárik use the sketch-and-project method to solve the average consensus problem; and Hanzely et al. design new variants of sketch and project methods for the average consensus problem with privacy considerations (see Section 8.3 for more details regarding the average consensus problem).

Primal Methods with Momentum

where $\omega>0$ is a stepsize and $\beta\geq 0$ is a momentum parameter. Instead of marrying the momentum term with gradient descent, we can marry it with SGD. This leads to SGD with momentum (mSGD), also known as the stochastic heavy ball method:

Since SGD is equivalent to SN and SPP, this way we obtain momentum variants of the stochastic Newton (mSN) and stochastic proximal point (mSPP) methods. The method is formally described below:

To the best of our knowledge, momentum variants of SN and SPP were not considered in the literature before. Moreover, as far as we know, there are no momentum variants of even deterministic variants of (11), (12) and (15), such as incremental or batch Newton method, incremental or batch proximal point method and incremental or batch projection method; not even for a problem formulated differently.

In the rest of this section we state our convergence results for mSGD/mSN/mSPP.

satisfy $a_{1}+a_{2}<1$. Let $x_{*}=\Pi_{\mathcal{L}}^{\mathbf{B}}(x_{0})$. Then

where $q=\frac{a_{1}+\sqrt{a_{1}^{2}+4a_{2}}}{2}$ and $\delta=q-a_{1}$. Moreover, $a_{1}+a_{2}\leq q<1$.

In the above theorem we obtain a global linear rate. To the best of our knowledge, this is the first time that linear rate is established for a stochastic variant of the heavy ball method (mSGD) in any setting. All existing results are sublinear. These seem to be the first momentum variants of SN and SPP methods.

If we choose $\omega\in(0,2)$, then the condition $a_{1}+a_{2}<1$ is satisfied for all

If $\beta=0$, mSGD reduces to SGD analyzed in . In this special case, $q=1-\omega(2-\omega)\lambda_{\min}^{+}$, which is the rate established in . Hence, our result is more general.

Let $q(\beta)$ be the rate as a function of $\beta$. Note that since $\beta\geq 0$, we have

Clearly, the lower bound on $q$ is an increasing function of $\beta$. Also, for any $\beta$ the rate is always inferior to that of SGD ($\beta=0$). It is an open problem whether one can prove a strictly better rate for mSGD than for SGD.

Our next theorem states that $\Pi_{\cal L}^{\bf B}(x_{k})=x_{*}$ for all iterations $k$ of mSGD. This invariance is important, as it allows the algorithm to converge to $x_{*}$.

Note that in view of (6), $\nabla f_{{\bf S}}(x)={\bf B}^{-1}{\bf A}^{\top}{\bf H}({\bf A}x-b)\in{\rm Range}({\bf B}^{-1}{\bf A}^{\top})$. Since

and since $x_{0}=x_{1}$, it can shown by induction that $x_{k}\in x_{0}+{\rm Range}({\bf B}^{-1}{\bf A}^{\top})$ for all $k$. However, ${\rm Range}({\bf B}^{-1}{\bf A}^{\top})$ is the orthogonal complement to ${\rm Null}({\bf A})$ in the ${\bf B}$-inner product. Since ${\cal L}$ is parallel to ${\rm Null}({\bf A})$, vectors $x_{k}$ must have the same ${\bf B}$-projection onto ${\cal L}$ for all $k$: $\Pi^{\bf B}_{\cal L}(x_{0})=x_{*}$. ∎

2 Cesaro average: sublinear rate without exactness assumption

In this section we present the convergence analysis of the function values computed on the Cesaro average. Again our results are global in nature. To the best of our knowledge are the first results that show $O(1/k)$ convergence of the stochastic heavy ball method. Existing results apply in more general settings at the expense of slower rates. In particular, and get $O(1/\sqrt{k})$ and $O(1/k^{\beta})$ convergence, respectively. When $\beta=1$, gets $O(1/\log(k))$ rate.

Choose $x_{0}=x_{1}$ and let $\{x_{k}\}_{k=0}^{\infty}$ be the random iterates produced by mSGD/mSN/mSPP, where the momentum parameter $0\leq\beta<1$ and relaxation parameter (stepsize) $\omega>0$ satisfy $\omega+2\beta<2$. Let $x_{*}$ be any vector satisfying $f(x_{*})=0$. If we let $\hat{x}_{k}=\frac{1}{k}\sum_{t=1}^{k}x_{t}$, then

In the special case of $\beta=0$, the above theorem gives the rate

This is the convergence rate for Cesaro averges of the “basic method” (i.e., SGD) established in .

Our proof strategy is similar to in which the first global convergence analysis of the (deterministic) heavy ball method was presented. There it was shown that when the objective function has a Lipschitz continuous gradient, the Cesaro averages of the iterates converge to the optimum at a rate of $O(1/k)$. To the best of our knowledge, there are no results in the literature that prove the same rate of convergence in the stochastic case for any class of objective functions.

In the authors analyzed mSGD for general Lipshitz continuous convex objective functions (with bounded variance) and proved the sublinear rate $O(1/\sqrt{k})$. In , a complexity analysis is provided for the case of quadratic strongly convex smooth coercive functions. A sublinear convergence rate of $O(1/k^{\beta})$, where $\beta\in(0,1)$, was proved. In contrast to our results, where we assume fixed stepsize $\omega$, both papers analyze mSGD with diminishing stepsizes.

3 L​1𝐿1L1 convergence: accelerated linear rate

In this section we show that by a proper combination of the relaxation (stepsize) parameter $\omega$ and the momentum parameter $\beta$, mSGD/mSN/mSPP enjoy an accelerated linear convergence rate in mean.

Dual Methods with Momentum

In the previous sections we focused on methods for solving the stochastic optimization problem (1) and the best approximation problem (3). In this section we focus on the dual of the best approximation problem, and propose a momentum variant of SDSA, which we call mSDSA.

In our first result we show that the random iterates of the mSGD/mSN/mSPP methods arise as an affine image of mSDSA under the mapping $\phi$ defined in (18).

Let $x_{0}=x_{1}$ and let $\{x_{k}\}$ be the iterates of mSGD/mSN/mSPP. Let $y_{0}=y_{1}=0$, and let $\{y_{k}\}$ be the iterates of mSDSA. Assume that the methods use the same stepsize $\omega>0$, momentum parameter $\beta\geq 0$, and the same sequence of random matrices ${\bf S}_{k}$. Then

for all $k$. That is, the primal iterates arise as affine images of the dual iterates.

So, the sequence of vectors $\{\phi(y_{k})\}$ mSDSA satisfies the same recursion of degree as the sequence $\{x_{k}\}$ defined by mSGD. It remains to check that the first two elements of both recursions coincide. Indeed, since $y_{0}=y_{1}=0$ and $x_{0}=x_{1}$, we have $x_{0}=\phi(0)=\phi(y_{0})$, and $x_{1}=x_{0}=\phi(0)=\phi(y_{1})$. ∎

2 Convergence

We are now ready to state a linear convergence convergence result describing the behavior of mSDSA in terms of the dual function values $D(y_{k})$.

satisfy $a_{1}+a_{2}<1$. Let $x_{*}=\Pi_{\mathcal{L}}^{\mathbf{B}}(x_{0})$ and let $y_{*}$ be any dual optimal solution. Then

where $q=\frac{a_{1}+\sqrt{a_{1}^{2}+4a_{2}}}{2}$ and $\delta=q-a_{1}$. Moreover, $a_{1}+a_{2}\leq q<1$.

This follows by applying Theorem 1 together with Theorem 5 and the identity $\tfrac{1}{2}\|x_{k}-x_{0}\|^{2}_{\bf B}=D(y_{*})-D(y_{k})$. ∎

Methods with Stochastic Momentum

In this case, mSGD becomes the randomized Kaczmarz method with momentum (mRK), and the iteration (22) takes the explicit form

Note that the cost of one iteration of this method is ${\cal O}(\|{\bf A}_{j:}\|_{0}+n)$, where the cardinality term $\|{\bf A}_{j:}\|_{0}$ comes from the stochastic gradient part, and $n$ comes from the momentum part. When ${\bf A}$ is sparse, the second term will dominate. Similar considerations apply for many other (but clearly not all) distributions ${\cal D}$.

Hence, we replace the momentum term by an unbiased estimator, which allows us to cut the cost to ${\cal O}(\|{\bf A}_{j:}\|_{0})$.

We now propose a variant of the SGD/SN/SPP methods employing stochastic momentum (smSGD/smSN/smSPP). Since SGD, SN and SPP are equivalent, we will describe the development from the perspective of SGD. In particular, we propose the following method:

2 Convergence

In the next result we establish L2 linear convergence of smSGD/smSN/smSPP. For this we will require the matrix ${\bf B}$ to be equal to the identity matrix.

satisfy $a_{1}+a_{2}<1$. Let $x_{*}=\Pi_{\mathcal{L}}^{\mathbf{I}}(x_{0})$. Then

It is straightforward to see that if we choose $\omega\in(0,2)$, then the condition $a_{1}+a_{2}<1$ is satisfied for all $\beta$ belonging to the interval

The upper bound is similar to that for mSGD/mSN/mSPP; the only difference is an extra factor of $n$ next to the constant 16.

3 Momentum versus stochastic momentum

As indicated in the introduction, if we wish to compare mSGD to smSGD used with momentum parameter $\beta$, it makes sense to use momentum parameter $\beta n$ in smSGD. This is because the momentum term in smSGD will then be an unbiased estimator of the deterministic momentum term used in mSGD.

Let $q(\beta)$ be the convergence constant for mSGD with stepsize $\omega=1$ and an admissible momentum parameter $\beta\geq 0$. Further, let $\bar{a}_{1}(\beta),\bar{a}_{2}(\beta),\bar{q}(\beta)$ be the convergence constants for smSGD with stepsize $\omega=1$ and momentum parameter $\beta$. We have

Hence, the lower bound on the rate for smSGD is worse than the lower bound for mSGD.

The same conclusion holds for the convergence rates themselves. Indeed, note that since $\bar{a}_{1}(\beta n)-a_{1}(\beta)=2\beta^{2}(n-1)\geq 0$ and $\bar{a}_{2}(\beta n)-a_{2}(\beta)=2\beta^{2}(n-1)\geq 0$, we have

and hence the rate of mSGD is always better than that of smSGD.

However, the expected cost of a single iteration of mSGD may be significantly larger than that of smSGD. Indeed, let $g$ be the expected cost of evaluating a stochastic gradient. Then we need to compare ${\cal O}(g+n)$ (mSGD) against ${\cal O}(g)$ (smSGD). If $g\ll n$, then one iteration of smSGD is significantly cheaper than one iteration of mSGD. Let us now compare the total complexity to investigate the trade-off between the rate and cost of stochastic gradient evaluation. Ignoring constants, the total cost of the two methods (cost of a single iteration multiplied by the number of iterations) is:

and since $q(\beta)$ and $\bar{q}(\beta n)$ are continuous functions of $\beta$, then because $g+n>g$, for small enough $\beta$ we will have $C_{\text{mSGD}}(\beta)>C_{\text{smSGD}}(\beta n).$ In particular, the speedup of smSGD compared to mSGD for $\beta\approx 0$ will be close to

Thus, we have shown the following statement.

For small $\beta$, the total complexity of smSGD is approximately $1+n/g$ times smaller than the total complexity of mSGD, where $n$ is the number of columns of ${\bf A}$, and $g$ is the expected cost of evaluating a stochastic gradient $\nabla f_{{\bf S}}(x)$.

Special Cases: Randomized Kaczmarz with Momentum and Randomized Coordinate Descent with Momentum

The updates for smSGD can be derived by substituting the momentum term $\beta(x_{k}-x_{k-1})$ with its stochastic variant $n\beta e_{i_{k}}^{\top}(x_{k}-x_{k-1})e_{i_{k}}$. We do not aim to be comprehensive. For more details on the possible combinations of the parameters $\mathbf{S}$ and $\mathbf{B}$ we refer the interested reader to Section 3 of .

We now provide a discussion on mRCD (the method in the first row of Table 4). Let $\mathbf{B}=\mathbf{I}$ and let pick in each iteration the random matrix $\mathbf{S}=e_{i}$ with probability $p_{i}=\|\mathbf{A}_{i:}\|^{2}/\|\mathbf{A}\|_{F}^{2}$. In this setup the update rule of the mSGD simplifies to

The objective function takes the following form:

For $\beta=0$, this method reduces to the randomized Kaczmarz method with relaxation, first analyzed in . If we also have $\omega=1$, this is equivalent with the randomized Kaczmarz method of Strohmer and Vershynin . RK without momentum ($\beta=0$) and without relaxation ($\omega=1$) converges with iteration complexity of

For $\omega=1$ and $\beta=\left(1-\sqrt{0.99\lambda_{\min}^{+}}\right)^{2}=\left(1-\sqrt{\frac{0.99}{\|\mathbf{A}\|^{2}_{F}}\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})}\right)^{2}$, the iteration complexity of the mRK is:

For $\omega=\|\mathbf{A}\|^{2}_{F}/\lambda_{\max}(\mathbf{A}^{\top}\mathbf{A})$ and $\beta=\left(1-\sqrt{\frac{0.99\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})}{\lambda_{\max}(\mathbf{A}^{\top}\mathbf{A})}}\right)^{2}$ the iteration complexity becomes:

This is quadratic improvement on the previous best result (36).

The Kaczmarz method for solving consistent linear systems was originally introduced by Kaczmarz in 1937 . This classical method selects the rows to project onto in a cyclic manner. In practice, many different selection rules can be adopted. For non-random selection rules (cyclic, greedy, etc) we refer the interested reader to . In this work we are interested in randomized variants of the Kaczmarz method, first analyzed by Strohmer and Vershynin . In it was shown that RK converges with a linear convergence rate to the unique solution of a full-rank consistent linear system. This result sparked renewed interest in design of randomized methods for solving linear systems . All existing results on accelerated variants of RK use the Nesterov’s approach of acceleration . To the best of our knowledge, no convergence analysis of mRK exists in the literature (Polyak’s momentum). Our work fills this gap.

2 mRCD: randomized coordinate descent with momentum

We now provide a discussion on the mRCD method (the method in the second row of Table 4). If the matrix $\mathbf{A}$ is positive definite, then we can choose $\mathbf{B}=\mathbf{A}$ and $\mathbf{S}=e_{i}$ with probability $p_{i}=\frac{\mathbf{A}_{ii}}{{\rm Trace}(\mathbf{A})}$. It is easy to see that $\mathbf{W}=\frac{\mathbf{A}}{{\rm Trace}(\mathbf{A})}$. In this case, $\mathbf{W}$ is positive definite and as a result, $\lambda_{\min}^{+}(\mathbf{W})=\lambda_{\min}(\mathbf{W})$. Moreover, we have

For $\beta=0$ and $\omega=1$ the method is equivalent with randomized coordinate descent of Leventhal and Lewis , which was shown to converge with iteration complexity

In contrast, following Theorem 4, we can obtain the following $L_{1}$ iteration complexity results for mRCD:

For $\omega=1$ and $\beta=\left(1-\sqrt{\frac{0.99}{{\rm Trace}(\mathbf{A})}\lambda_{\min}(\mathbf{A})}\right)^{2}$, the iteration complexity is

For $\omega={\rm Trace}(\mathbf{A})/\lambda_{\max}(\mathbf{A})$ and $\beta=\left(1-\sqrt{\frac{0.99\lambda_{\min}(\mathbf{A})}{\lambda_{\max}(\mathbf{A})}}\right)^{2}$ the iteration complexity becomes

This is quadratic improvement on the previous best result (38).

It is known that if ${\bf A}$ is positive definite, the popular randomized Gauss-Seidel method can be interpreted as randomized coordinate descent (RCD). RCD methods were first analyzed by Lewis and Leventhal in the context of linear systems and least-squares problems , and later extended by several authors to more general settings, including smooth convex optimization , composite convex optimization , and parallel/subspace descent variants . These results were later further extended to handle arbitrary sampling distributions . Accelerated variants of RCD were studied in . For other non-randomized coordinate descent variants and their convergence analysis, we refer the reader to . To the best of our knowledge, mRCD and smRCD have never been analyzed before in any setting.

3 Visualizing the acceleration mechanism

We devote this section to the graphical illustration of the acceleration mechanism behind momentum. Our goal is to shed more light on how the proposed algorithm works in practice. For simplicity, we illustrate this by comparing RK and mRK.

The Projection: The projection step corresponds to the first part $x_{k}-\omega\nabla f_{\mathbf{S}_{k}}(x_{k})$ of the mRK update (22) and it means that the current iterate $x_{k}$ is projected onto a randomly chosen hyperplane $H_{i}$In the plots of Figure 1, the hyperplane of each update is chosen in an alternating fashion for illustration purposes. The value of the stepsize $\omega\in(0,2)$ defines whether the projection is exact or not. When $\omega=1$ (no relaxation) the projection is exact, that is the point $\Pi_{H_{i}}(x_{k})$ belongs in the hyperplane $H_{i}$. In Figure 1 all projections are exact.

Addition of the momentum term: The momentum term (right part of the update rule) $\beta(x_{k}-x_{k-1})$ forces the next iterate $x_{k+1}$ to be closer to the solution $x_{*}$ than the corresponding point $\Pi_{H_{i}}(x_{k})$. Note also that the vector $x_{k+1}-\Pi_{H_{i}}(x_{k})$ is always parallel to $x_{k}-x_{k-1}$ for all $k\geq 0$.

In the example of Figure 1, the performance of mRK is similar to the performance of RK until iterate $x_{3}$. After this point, the momentum parameter becomes more effective and the mRK method accelerates. This behavior appears also in our experiments in the next section where we work with matrices with many rows. There we can notice that the momentum parameter seems to become more effective after the first $m+1$ iterations.

Numerical Experiments

In comparing the methods with their momentum variants we use both the relative error measure $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}/\|x_{0}-x_{*}\|^{2}_{\mathbf{B}}$ and the function values $f(x_{k})$Remember that in our setting we have $f(x_{*})=0$ for the optimal solution $x_{*}$ of the best approximation problem; thus $f(x)-f(x_{*})=f(x)$. The function values $f(x_{k})$ refer to function (35) in the case of RK and to function (37) for the RCD. For block variants the objective function of problem (1) has also closed form expression but it can be very difficult to compute. In these cases one can instead evaluate the quantity $\|\mathbf{A}x-b\|^{2}_{\mathbf{B}}$.. In all implementations, except for the experiments on average consensus (Section 8.3), the starting point is chosen to be $x_{0}=0$. In the case of average consensus the starting point must be the vector with the initial private values of the nodes of the network. All the code for the experiments is written in the Julia programming language. For the horizontal axis we use either the number of iterations or the wall-clock time measured using the tic-toc Julia function.

This section is divided in three main experiments. In the first one we evaluate the performance of the mSGD method in the special cases of mRK and mRCD for solving both synthetic consistent Gaussian systems and consistent linear systems with real matrices. In the second experiment we computationally verify Theorem 8 (comparison between the mSGD and smSGD methods). In the last experiment building upon the recent results of we show how the addition of the momentum accelerates the pairwise randomized gossip (PRG) algorithm for solving the average consensus problem.

In this subsection we study the computational behavior of mRK and mRCD when they compared with their no momentum variants for both synthetic and real data.

The synthetic data for this comparison is generated as followsNote that in the first experiment we use Gaussian matrices which by construction are full rank matrices with probability 1 and as a result the consistent linear systems have unique solution. Thus, for any starting point $x_{0}$, the vector $z$ that used to create the linear system is the solution mSGD converges to. This is not true for general consistent linear systems, with no full-rank matrix. In this case, the solution $x_{*}=\Pi_{{\cal L}}^{\mathbf{B}}(x_{0})$ that mSGD converges to is not necessarily equal to $z$. For this reason, in the evaluation of the relative error measure $\|x_{k}-x_{*}\|^{2}_{\mathbf{B}}/\|x_{0}-x_{*}\|^{2}_{\mathbf{B}}$, one should be careful and use the value $x_{*}=x_{0}+\mathbf{A}^{\dagger}(b-\mathbf{A}x_{0})\overset{x_{0}=0}{=}\mathbf{A}^{\dagger}b$..

For each linear system we run mRK (Figure 2) and mRCD (Figure 3) for several values of momentum parameters $\beta$ and fixed stepsize $\omega=1$ and we plot the performance of the methods (average after 10 trials) for both the relative error measure and the function values. Note that for $\beta=0$ the methods are equivalent with their no-momentum variants RK and RCD respectively.

From Figures 2 and 3 it is clear that the addition of momentum term leads to an improvement in the performance of RK and RCD, respectively. More specifically, from the two figures we observe the following:

For the well conditioned linear systems ($1/\lambda_{\min}^{+}$ small) it is known that even the no-momentum variant converges rapidly to the optimal solution. In these cases the benefits of the addition of momentum are not obvious. The momentum term is beneficial for the case where the no-momentum variant ($\beta=0$) converges slowly, that is when $1/\lambda_{\min}^{+}$ is large (ill-conditioned linear systems).

For the case of fixed stepsize $\omega=1$, the problems with small condition number require smaller momentum parameter $\beta$ to have faster convergence. Note the first two rows of Figures 2 and 3, where $\beta=0.3$ or $\beta=0.4$, are good options.

We observe that both mRK and mRCD, with appropriately chosen momentum parameters $0<\beta\leq 0.5$, always converge faster than their no-momentum variants, RK and RCD, respectively. This is a smaller momentum parameter than $\beta\approx 0.9$ which is being used extensively with mSGD for training deep neural networks .

In a stochastic power iteration with momentum is proposed for principal component analysis (PCA). There it was demonstrated empirically that a naive application of momentum to the stochastic power iteration does not result in a faster method. To achieve faster convergence, the authors proposed mini-batch and variance-reduction techniques on top of the addition of momentum. In our setting, mere addition of the momentum term to SGD (same is true for special cases such as RK and RCD) leads to empirically faster methods.

1.2 Real Data

In Figure 4 the performance of all methods for both relative error measure $\|x_{k}-x_{*}\|^{2}/\|x_{*}\|^{2}_{\mathbf{B}}$ and function values $f(x_{k})$ is presented. Note again that $\beta=0$ represents the baseline RK method. The addition of momentum parameter is again often beneficial and leads to faster convergence. As an example, inspect the plots for the mushrooms dataset in Figure 4, where mRK with $\beta=0.5$ is much faster than the simple RK method in all presented plots, both in terms of iterations and time. In particular, the addition of a momentum parameter leads to visible speedup for the datasets mushrooms, splice, a9a and ionosphere. For these datasets the acceleration is obvious in all plots both in terms of relative error and function values. For the datasets australian, gisette and madelon the speedup is less obvious in the plots of the relative error, while for the plots of function values it is not present at all.

2 Comparison of momentum & stochastic momentum

In Theorem 8, the total complexities (number of operations needed to achieve a given accuracy) of mSGD and smSGD have been compared and it has been shown that for small momentum parameter $\beta$,

where $C_{\text{mSGD}}$ and $C_{\text{smSGD}}$ represent the total costs of the two methods. The goal of this experiment is to show that this relationship holds also in practice.

For this experiment we assume that the non-zeros of matrix $\mathbf{A}$ are not concentrated in certain rows but instead that each row has the same number of non-zero coordinates. We denote by $g$ the number the non-zero elements per row. Having this assumption it can be shown that for the RK method the cost of one projection is equal to $4g$ operations while the cost per iteration of the mRK and of the smRK are $4g+3n$ and $4g+1$ respectively. For more details about the cost per iteration of the general mSGD and smSGD check Table 6.

3 Faster method for average consensus

It was shown recently that several randomized methods for solving linear systems can be interpreted as randomized gossip algorithms for solving the AC problem when applied to a special system encoding the underlying network . As we have already explained both basic method and basic method with momentum (this paper) find the solution of the linear system that is closer to the starting point of the algorithms. That is, both methods converge linearly to $x_{*}=\Pi^{{\bf B}}_{{\cal L}}(x_{0})$; the projection of the initial iterate onto the solution set of the linear system and as a result (check Introduction) can be interpreted as methods for solving the best approximation problem (16). In the special case that

the starting point of the method are the initial values of the nodes $x_{0}=c$,

it is straightforward to see that the solution of the best approximation problem is a vector with all components equal to the consensus value $\bar{c}:=\tfrac{1}{n}\sum_{i}c_{i}$. Under this setting, the famous randomized pairwise gossip algorithm (randomly pick an edge $e\in E$ and replace the private values of its two nodes to their average) that was first proposed and analyzed in , is equivalent with the RK method without relaxation ($\omega=1$) .

In the gossip framework, the condition number of the linear system when RK is used has a simple structure and it depends on the characteristics of the network under study. More specifically, it depends on the number of the edges $m$ and on the Laplacian matrix of the networkMatrix $\mathbf{A}$ of the linear system is the incidence matrix of the graph and it is known that the Laplacian matrix is equal to $\mathbf{L}=\mathbf{A}^{\top}\mathbf{A}$, where $\|\mathbf{A}\|^{2}_{F}=2m$.:

where $\mathbf{L}=\mathbf{A}^{\top}\mathbf{A}$ is the Laplacian matrix of the network and the quantity $\lambda_{\min}^{+}(\mathbf{L})$ is the very well studied algebraic connectivity of the graph .

The convergence analysis in this paper holds for any consistent linear system $\mathbf{A}x=b$ without any assumption on the rank of the matrix $\mathbf{A}$. The lack of any assumption on the form of matrix $\mathbf{A}$ allows us to solve the homogeneous linear system $\mathbf{A}x=0$ where $\mathbf{A}$ is the incidence matrix of the network which by construction is rank deficient. More specifically, it can be shown that ${\rm rank}(\mathbf{A})=n-1$ . Note that many existing methods for solving linear systems make the assumption that the matrix $\mathbf{A}$ of the linear systems is full rank and as a result can not be used to solve the AC problem.

3.2 Numerical Setup

Our goal in this experiment is to show that the addition of the momentum term to the randomized pairwise gossip algorithm (RK in the gossip setting) can lead to faster gossip algorithms and as a result the nodes of the network will converge to the average consensus faster both in number of iterations and in time. We do not intend to analyze the distributed behavior of the method (this is on-going research work). In our implementations we use three of the most popular graph topologies in the literature of wireless sensor networks. These are the line graph, cycle graph and the random geometric graph $G(n,r)$. In practice, $G(n,r)$ consider ideal for modeling wireless sensor networks, because of their particular formulation. In the experiments the $2$-dimensional $G(n,r)$ is used which is formed by placing $n$ nodes uniformly at random in a unit square with edges only between nodes that have euclidean distance less than the given radius $r$. To preserve the connectivity of $G(n,r)$ a radius $r=r(n)=\log(n)/n$ is used . The AC problem is solved for the three aforementioned networks for both $n=100$ and $n=200$ number of nodes. We run mRK with several momentum parameters $\beta$ for 10 trials and we plot their average. Our results are available in Figures 6 and 7.

Note that the vector of the initial values of the nodes can be chosen arbitrarily, and the proposed algorithms will find the average of these values. In Figures 6 and 7 the initial value of each node is chosen independently at random from the uniform distribution in the interval $(0,1)$.

3.3 Experimental Results

By observing Figures 6 and 7, it is clear that the addition of the momentum term improves the performance of the popular pairwise randomized gossip (PRG) method . The choice $\beta=0.4$ as the momentum parameter improves the performance of the vanilla PRG for all networks under study and $\beta=0.5$ is a good choice for the cases of the cycle and line graph. Note that for networks such as the cycle and line graphs there are known closed form expressions for the algebraic connectivity . Thus, using equation (39), we can compute the exact values of the condition number $1/\lambda_{\min}^{+}$ for these networks. Interestingly, as we can see in Table 7 for $n=100$ and $n=200$ (number of nodes), the condition number $1/\lambda_{\min}^{+}$ appearing in the iteration complexity of our methods is not very large. This is in contrast with experimental observations from Section 8.1.1 where it was shown that the choice $\beta=0.5$ is good for very ill conditioned problems only ($1/\lambda_{\min}^{+}$ very large).

Appendix A Proof of Theorem 1

Fix $F_{1}=F_{0}\geq 0$ and let $\{F_{k}\}_{k\geq 0}$ be a sequence of nonnegative real numbers satisfying the relation

where $a_{2}\geq 0$, $a_{1}+a_{2}<1$ and at least one of the coefficients $a_{1},a_{2}$ is positive. Then the sequence satisfies the relation $F_{k+1}\leq q^{k}(1+\delta)F_{0}$ for all $k\geq 1,$ where $q=\frac{a_{1}+\sqrt{a_{1}^{2}+4a_{2}}}{2}$ and $\delta=q-a_{1}\geq 0$. Moreover,

with equality if and only if $a_{2}=0$ (in which case $q=a_{1}$ and $\delta=0$).

Choose any $\delta\geq 0$ satisfying $a_{2}\leq(a_{1}+\delta)\delta$. Adding $\delta F_{k}$ to both sides of (40), we get

We now claim that $\delta=\frac{-a_{1}+\sqrt{a_{1}^{2}+4a_{2}}}{2}$ satisfies the relations. Non-negativity of $\delta$ follows from $a_{2}\geq 0$, while the second relation follows from the fact that $\delta$ satisfies

Let us now argue that $0<q<1$. Nonnegativity of $q$ follows from nonnegativity of $a_{2}$. Clearly, as long as $a_{2}>0$, $q$ is positive. If $a_{2}=0$, then $a_{1}>0$ by assumption, which implies that $q$ is positive. The inequality $q<1$ follows directly from the assumption $a_{1}+a_{2}<1$. By unrolling the recurrence (42), we obtain $F_{k+1}\leq F_{k+1}+\delta F_{k}\leq q^{k}(F_{1}+\delta F_{0})=q^{k}(1+\delta)F_{0}.$

Finally, let us establish (42). Noting that $a_{1}=q-\delta$, and since in view of (43) we have $a_{2}=q\delta$, we conclude that $a_{1}+a_{2}=q+\delta(q-1)\leq q$, where the inequality follows from $q<1$. ∎

The following identities were established in . For completeness, we include different (and somewhat simpler) proofs here.

In view of (10), and since ${\bf Z}{\bf B}^{-1}{\bf Z}={\bf Z}$ (see ), we have

By taking expectations in the last identity with respect to the random matrix $\mathbf{S}$, we get $\langle\nabla f(x),x-x_{*}\rangle_{\mathbf{B}}=2f(x).$ ∎

Moreover, if exactness is satisfied, and we let $x_{*}=\Pi^{\mathbf{B}}_{{\cal L}}(x)$, we have

A.2 The Proof

We will now analyze the three expressions 1, 2, 3 separately. The first expression can be written as

We will now bound the second expression. First, we have

By substituting the bounds (51), (53), (54) into (50) we obtain

Now by first taking expectation with respect to ${\bf S}_{k}$, we obtain:

where in the second step we used the inequality $\langle\nabla f(x_{k}),x_{k-1}-x_{k}\rangle\leq f(x_{k-1})-f(x_{k})$ and the fact that $\omega\beta\geq 0$, which follows from the assumptions. We now apply inequalities (48) and (49), obtaining

It suffices to apply Lemma 9 to the relation (55). The conditions of the lemma are satisfied. Indeed, $a_{2}\geq 0$, and if $a_{2}=0$, then $\beta=0$ and hence $a_{1}=1-\omega(2-\omega)\lambda_{\min}^{+}>0$. The condition $a_{1}+a_{2}<1$ holds by assumption.

Appendix B Proof of Theorem 3

Let $p_{t}=\frac{\beta}{1-\beta}(x_{t}-x_{t-1})$ and $d_{t}=\|x_{t}+p_{t}-x_{*}\|_{{\bf B}}^{2}$. In view of (22), we can write

Taking expectation with respect to the random matrix $\mathbf{S}_{t}$ we obtain:

where the inequality follows from convexity of $f$. After rearranging the terms we get

where $\alpha=\frac{4\omega}{1-\beta}-\frac{2\omega^{2}}{(1-\beta)^{2}}>0$. Taking expectations again and using the tower property, we get

Finally, using Jensen’s inequality, we get

It remains to note that $\theta_{1}=\|x_{0}-x_{*}\|_{{\bf B}}^{2}+\frac{2\omega\beta}{(1-\beta)^{2}}f(x_{0}).$

Appendix C Proof of Theorem 4

In the proof of Theorem 4 the following two lemmas are used.

Consider the second degree linear homogeneous recurrence relation:

where M=\bigg{(}\sqrt{\frac{a_{1}^{2}}{4}+\frac{(-a_{1}^{2}-4a_{2})}{4}}\bigg{)}=\sqrt{-a_{2}} and $\theta$ is such that $a_{1}=2M\cos(\theta)$ and $\sqrt{-a_{1}^{2}-4a_{2}}=2M\sin(\theta)$.

We can now turn to the proof of Theorem 4. Plugging in the expression for the stochastic gradient, mSGD can be written in the form

Subtracting $x_{*}$ from both sides of (59), we get

Multiplying the last identity from the left by $\mathbf{B}^{1/2}$, we get

Taking expectations, conditioned on $x_{k}$ (that is, the expectation is with respect to $\mathbf{S}_{k}$):

Taking expectations again, and using the tower property, we get

which can be written in a coordinate-by-coordinate form as follows:

where $s_{k}^{i}$ indicates the $i$th coordinate of $s_{k}$.

We will now fix $i$ and analyze recursion (62) using Lemma 13. Note that (62) is a second degree linear homogeneous recurrence relation of the form (58) with $a_{1}=1+\beta-\omega\lambda_{i}$ and $a_{2}=-\beta$. Recall that $0\leq\lambda_{i}\leq 1$ for all $i$. Since we assume that $0<\omega\leq 1/\lambda_{\max}$, we know that $0\leq\omega\lambda_{i}\leq 1$ for all $i$. We now consider two cases:

Applying Theorem 2, we know that $x_{*}=\Pi_{\mathcal{L}}^{\mathbf{B}}(x_{0})=\Pi_{\mathcal{L}}^{\mathbf{B}}(x_{1})$. Using Lemma 12 twice, once for $x=x_{0}$ and then for $x=x_{1}$, we observe that $s_{0}^{i}=u_{i}^{\top}\mathbf{B}^{1/2}(x_{0}-x_{*})=0$ and $s_{1}^{i}=u_{i}^{\top}\mathbf{B}^{1/2}(x_{1}-x_{*})=0$. Finally, in view of (63) we conclude that

Since $0<\omega\lambda_{i}\leq 1$ and $\beta\geq 0$, we have $1+\beta-\omega\lambda_{i}\geq 0$ and hence

where the last inequality can be shown to holdThe lower bound on $\beta$ is tight. However, the upper bound is not. However, we do not care much about the regime of large $\beta$ as $\beta$ is the convergence rate, and hence is only interesting if smaller than 1. for $(1-\sqrt{\omega\lambda_{\min}^{+}})^{2}<\beta<1$. Applying Lemma 13 the following bound can be deduced

where $P_{i}$ is a constant depending on the initial conditions (we can simply choose $P_{i}=|C_{0}|+|C_{1}|$).

Now putting the two cases together, for all $k\geq 0$ we have

where $C=4\sum_{i:\lambda_{i}>0}P_{i}^{2}$.

Appendix D Proof of Theorem 7

The proof follows a similar pattern to that of Theorem 1. However, stochasticity in the momentum term introduces an additional layer of complexity, which we shall tackle by utilizing a more involved version of the tower property.

For simplicity, let $i=i_{k}$ and $r_{k}^{i}:=e_{i}^{\top}(x_{k}-x_{k-1})e_{i}$. First, we decompose

We shall use the tower property in the form

where $X$ is some random variable. We shall perform the three expectations in order, from the innermost to the outermost. Applying the inner expectation to the identity (66), we get

We will now analyze the three expressions 1, 2, 3 separately. The first expression is constant under the expectation, and hence we can write

We will now bound the second expression. Using the identity

By substituting the bounds (69), (72), (73) into (68) we obtain

We now take the middle expectation (see (67)) and apply it to inequality (75):

where in the second step we used the inequality $\langle\nabla f(x_{k}),x_{k-1}-x_{k}\rangle\leq f(x_{k-1})-f(x_{k})$ and the fact that $\omega\beta\geq 0$, which follows from the assumptions. We now apply inequalities (48) and (49), obtaining

It suffices to apply Lemma 9 to the relation (55). The conditions of the lemma are satisfied. Indeed, $a_{2}\geq 0$, and if $a_{2}=0$, then $\beta=0$ and hence $a_{1}=1-\omega(1-\omega)\lambda_{\min}^{+}>0$. The condition $a_{1}+a_{2}<1$ holds by assumption.

The convergence result in function values follows as a corollary by applying inequality (48) to (30).

References