Momentum and Stochastic Momentum for Stochastic Gradient, Newton, Proximal Point and Subspace Descent Methods
Nicolas Loizou Peter Richtárik
math.OC cs.LG cs.NA math.NA stat.ML




Momentum and Stochastic Momentum for Stochastic  Gradient, Newton, Proximal Point and Subspace Descent Methods
    Nicolas Loizou  School of Mathematics, The University of Edinburgh. —  E-mail: n.loizou@sms.ed.ac.uk
 Peter RichtárikCEMSE, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia —  School of Mathematics, The University of Edinburgh, United Kingdom — Moscow Institute of Physics and Technology (MIPT), Dolgoprudny, Moscow, Russia.— E-mail: peter.richtarik@ed.ac.uk
    December 22, 2017[A short version of this paper (5 pages) was posted on arXiv on 30 Oct 2017 <cit.>. The paper was accepted for presentation at the 2017 NIPS Optimization for Machine Learning workshop in a peer reviewed process. The accepted papers are listed on the website of the workshop, but are not published in any proceedings volume.]
==========================================================================================================================================================================================================================================================================================================================================================================================================

 

 
In this paper we study several classes of stochastic optimization algorithms enriched with  heavy ball momentum. Among the methods studied are: stochastic gradient descent, stochastic Newton, stochastic proximal point  and stochastic dual subspace ascent. This is the first time momentum variants of several of these methods are studied. We choose to perform our analysis in a setting in which all of the above methods are equivalent. We prove global nonassymptotic linear convergence rates for all methods and various measures of success, including primal function values, primal iterates (in L2 sense), and dual function values. We also show that the primal iterates converge at an accelerated linear rate in the L1 sense. This is the first time a linear rate is shown for the stochastic heavy ball method (i.e., stochastic gradient descent method with momentum). Under somewhat weaker conditions,  we establish a sublinear convergence rate for Cesaro averages of primal iterates.  Moreover, we propose a novel concept, which we call stochastic momentum, aimed at decreasing the cost of performing the momentum step. We prove linear convergence of several stochastic methods with stochastic momentum, and show that in some sparse data regimes and for sufficiently small momentum parameters, these methods enjoy better overall complexity than methods with deterministic momentum. Finally, we perform extensive numerical testing on artificial and real datasets, including data coming from average consensus problems. 




68Q25 · 68W20 · 68W40 · 65Y20 · 90C15 · 90C20 · 90C25 · 15A06 · 15B52 · 65F10 




§ 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 <cit.> 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 <cit.> 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.



 §.§ Setting


In this paper we study three closely related problems: 


  (i) stochastic optimization, 

  (ii) best approximation, and 

  (iii) (bounded) concave quadratic maximization.
 
These problems and the relationships between them are described in detail in Section <ref>. Here we only briefly outline some of the key relationships. By stochastic optimization we refer to the problem of the form

    min_x∈^n f(x)f_(x),

where the expectation is over random matrices  drawn from an arbitrary  distribution , and f_ is a stochastic convex quadratic function of a  least-squares type, depending on , and in addition on a matrix ∈^m× n, vector b∈^m,  and an n× n positive definite matrix  (see Section <ref> for full details). The problem is constructed in such a way that the set of minimizers of f is identical to the set of solutions of a given (consistent) linear system 
    x =b,

where ∈^m× n and b ∈^m. In this sense, (<ref>) can  be seen as the reformulation of the linear system (<ref>) into a stochastic optimization problem. Such reformulations  provide an explicit connection between the fields of  linear algebra and stochastic optimization, which may inspire future research by enabling the transfer of knowledge, techniques, and algorithms from one field to another. For instance, the randomized Kaczmarz method of Strohmer and Vershynin <cit.> for solving (<ref>) is equivalent to the stochastic gradient descent method applied to (<ref>), with  corresponding to a discrete distribution over unit coordinate vectors in ^m. However, the flexibility of being able to choose  arbitrarily allows for numerous generalizations of the randomized Kaczmarz method <cit.>. Likewise, provably faster variants of the randomized Kaczmarz method (for instance, by utilizing importance sampling) can be designed using the connection. 



 §.§ Three stochastic methods

Problem  (<ref>)  has several peculiar characteristics which are of key importance to this paper. For instance, the Hessian of f_ is a (random) projection matrix, which can be used to show that f_(x) = 12∇ f_(x)_^2 . Moreover, it follows that the Hessian of f has all eigenvalues bounded by 1, and so on. These characteristics can be used to show that several otherwise distinct  stochastic algorithms for solving the stochastic optimization problem (<ref>) are identical <cit.>. In particular, the following optimization methods for solving (<ref>) are identical[In addition, these three methods are identical to a stochastic fixed point method (with relaxation) for solving the fixed point problem x = Π__(x), where _ is the set of solutions of ^⊤ x = ^⊤ b, which is a  sketched version of the  linear system (<ref>), and can be seen as a stochastic approximation  of the set {x :  x = b}]: 


  * 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 ω>0.  The methods will be described in detail in Section <ref>; see also Table <ref> 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.



 §.§ Best approximation, duality and  stochastic dual subspace ascent
 It was shown in <cit.> in the ω=1 case and in <cit.> in the general ω>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 (<ref>). This naturally leads  to the best approximation problem, which is the problem of projecting[In 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 (<ref>):

    min_x∈^n12x-x_0_^2   subject to   x = b .



The dual of the best approximation problem is an unconstrained concave quadratic maximization problem <cit.>.  Consistency of x = b implies that the dual is bounded.  It follows from the results of <cit.> that for ω=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 (<ref>), known as


  * Method 4: stochastic dual subspace ascent (SDSA).

In this paper we show that this equivalence extends beyond the ω=1 case, specifically for 0<ω<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  randomly drawn in each iteration from . In this subspace, the method moves to the point which maximizes the dual objective, D(y). Since  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 <cit.>.



 §.§ Structure of the paper


The remainder of this work is organized as follows. In Section <ref> we summarize our contributions in the context of existing literature. In Section <ref> 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 <ref> we describe and analyze  primal  methods  with momentum (mSGD, mSN and mSPP), and in Section <ref> we describe and analyze the dual method with momentum (mSDSA). In Section <ref> we describe and analyze primal methods with stochastic momentum (smSGD, smSN and smSPP). Numerical experiments are presented in Section <ref>.  Proofs of all key results can be found in the appendix.




 

 §.§ Notation


The following notational conventions are used in this paper. Boldface upper-case letters denote matrices;  is the identity matrix. By  we denote the solution set of the linear system x=b. By _, where  is a random matrix, we denote the solution set of the sketched linear system ^⊤ x= ^⊤ b.  Throughout the paper,  is an n× n positive definite matrix giving rise to an inner product and norm on ^n. Unless stated otherwise, throughout the paper, x_* is the projection of x_0 onto  in the -norm: x_*=Π_^(x_0). We write [n]{1,2, … ,n}.




§ MOMENTUM METHODS AND OUR CONTRIBUTIONS
 


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




 §.§ Heavy ball method


The baseline first-order method for minimizing a differentiable function f is the gradient descent (GD) method, 
    x_k+1 = x_k - ω_k ∇ f(x_k),
 where ω_k>0 is a stepsize. For convex functions with L-Lipschitz gradient (function class ^1,1_0,L),  GD converges at at the rate of (L/ϵ). When, in addition, f is μ-strongly convex (function class ^1,1_μ,L), the rate is linear: ((L /μ) log(1/ϵ)) <cit.>. To improve the convergence behavior of the method, Polyak proposed to modify GD by the   introduction of a (heavy ball) momentum term[Arguably a much more popular, certainly theoretically much better understood alternative to Polyak's momentum is the momentum introduced by Nesterov <cit.>, leading to the famous accelerated gradient descent (AGD) method. This method converges nonassymptotically and globally; with optimal sublinear rate (√(L/ϵ))  <cit.> when  applied to minimizing a smooth convex objective function (class ^1,1_0,L), and with the optimal linear  rate (√(L/μ)log(1/ϵ)) when minimizing smooth strongly convex functions (class ^1,1_μ,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.], β(x_k-x_k-1). This leads to the gradient descent method with  momentum (mGD), popularly  known as the heavy ball method: 
    x_k+1 = x_k - ω_k ∇ f(x_k) + β (x_k-x_k-1).
  More specifically, Polyak proved that with the correct choice of the stepsize parameters ω_k and momentum parameter β, a local accelerated linear convergence rate of  (√(L/μ)log(1/ϵ))  can be achieved in the case of twice continuously differentiable, μ-strongly convex objective functions with L-Lipschitz gradient (function class ℱ_μ, L^2,1).  See the first line of Table <ref>. 

  
 
Recently, Ghadimi et al. <cit.> 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∈ℱ_0, L^1,1), and  linearly when the function is also strongly convex (f∈ℱ_μ, L^1,1). A particular, selection of the parameters ω and β 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∈ℱ_μ, L^1,1. However, in the special case of a strongly convex quadratic, an elegant proof was recently proposed in <cit.>. Using the notion of integral quadratic constraints from robust control theory, the authors proved that by choosing ω_k = ω=4 / (√(L)+√(μ))^2 and β=(√(L/μ)-1)^2/(√(L/μ)+1)^2, the heavy ball method enjoys a global asymptotic accelerated convergence rate of (√(L/μ)log(1/ϵ)). The aforementioned results are summarized in the first part of Table <ref>.


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



 §.§ 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):

    x_k+1 = x_k -ω_k g (x_k)  + β (x_k-x_k-1),

where g_k is an unbiased estimator of the true gradient ∇ f(x_k). While mSGD is used extensively in practice, especially in deep learning <cit.>, 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. <cit.>, and the work of Gadat et al. <cit.>. 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/√(ϵ)) was proved. For this, the authors employed a decreasing stepsize strategy: ω_k=ω_0/√(k+1), where ω_0 is a positive constant.  In <cit.>, 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 ω_k=ω_0 / k^β was shown to convergence as (1/k^β) when the momentum parameter is β <1, and with the rate O(1/log k) when β=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 <ref>. 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.
 

 


 §.§ Connection to incremental gradient methods


Assuming  is discrete distribution (i.e., we sample from M matrices, ^1,…,^M, where ^i is chosen with probability p_i>0), we can write the stochastic optimization problem (<ref>) in the finite-sum form

    min_x∈^n f(x) = ∑_i=1^M p_i f_^i(x).

Choosing x_0=x_1, mSGD with fixed stepsize ω_k=ω applied to (<ref>) can be written in the form

    x_k+1 =  x_k-ω∑_t=1^kβ^k-t∇ f__t (x_t)+β^k(x_1-x_0)
    = x_k-ω∑_t=1^kβ^k-t∇ f__t (x_t),
 
where _t = ^i with probability p_i. Problem (<ref>) can be also solved using incremental  average/aggregate gradient methods, such as the IAG method of Blatt et al. <cit.>. These methods have a similar form to (<ref>), with the main difference being in the way the past gradients are aggregated. While (<ref>) 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. <cit.> 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  <cit.>; see also the work of Tseng  <cit.>. 

In <cit.>, 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 (<ref>). It is not a finite sum problem as the distribution  could be continous; and we also do not assume strong convexity. Thus, the convergence analysis of <cit.> can not be directly applied to our problem.  




 §.§ 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 <ref>). 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 (<ref>), the best approximation problem (<ref>) and its dual.  To the best of our knowledge, momentum variants of SN, SPP and SDSA were not analyzed before. 




Linear rate. We prove several (global and non-asymptotic) linear convergence results for our primal momentum methods mSGD/mSN/mSPP. First, we establish a linear rate for the decay of x_k-x_*_^2 to zero (i.e., L2 convergence), for a range of stepsizes ω> 0 and momentum parameters β≥ 0. We show that the same rate holds for the decay of the expected function values f(x_k)-f(x_*) of (<ref>) to zero. Further, the same rate holds for  mSDSA, in particular, this is for the convergence of the dual objective to the optimum. For a summary of these results, and pointers to the relevat theorems, refer to lines 1, 2 and 6 of Table <ref>. Unfortunately, the theoretical rate for all our momentum methods is optimized for β = 0, and  gets worse as the momentum parameter increases. However, no prior linear rate for any of these methods with momentum are known. We give the first linear convergence rate for SGD with momentum (i.e., for the stochastic heavy ball method). 



Accelerated linear rate. We then study the decay of the larger quantity x_k-x_*_^2  to zero (i.e., L1 convergence).  In this case, we establish  an accelerated linear rate, which depends on the square root of the condition number (of the Hessian of f). This is a quadratic speedup when compared to  the no-momentum methods as these depend on the condition number. See lines 4 and 5 of Table <ref>. To the best of our knowledge, this is the first time an  accelerated rate is obtained for the stochastic heavy ball method (mSGD).  Note that there are no global non-asymptotic accelerated linear rates proved even in the non-stochastic setting (i.e., for the heavy ball method). Moreover, we are not aware of any accelerated linear convergence results for  the stochastic proximal point method.



Sublinear rate for Cesaro averages. We show that the Cesaro averages, x̂_k = 1/k∑_t=0^k-1x_t, of all primal momentum methods enjoy a sublinear O(1/k) rate (see line 3 of Table <ref>). 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 <ref>). This extends the result of <cit.> where this was shown for the no-momentum methods (β=0) and in the special case of the unit stepsize (ω=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  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.




 §.§ 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 <cit.>, SDCA <cit.>, SVRG/S2GD  <cit.>, minibatch SVRG/S2GD <cit.>, and SAGA <cit.>.
 
Since we assume that the linear system (<ref>) is feasible, it follows that the stochastic gradient vanishes at the optimal point (i.e., ∇ f_(x_*)=0 for any ). 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 <cit.>. 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áč <cit.>. In particular, given the consistent linear system (<ref>), 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 <cit.>, the solution sets  of these reformulations are equal to the solution set of the linear system.



 §.§ Stochastic optimization


 Stochasticity enters the reformulations via a user defined distribution  of matrices (all with m rows). In addition, the reformulations  utilize a positive definite matrix ∈^n× n as a parameter,  used to define an inner product in ^n via ⟨ x,z ⟩_⟨ x, z⟩ and the induced norm
x_ (x^⊤ x)^1/2.   In particular, the stochastic optimization reformulation (<ref>), i.e., 
min_x∈^n f(x) f_(x),
is defined by
setting  
    f_(x) 1/2 x - b_^2 = 1/2( x - b)^⊤ ( x - b),
 where  is a random symmetric positive semidefinite matrix defined as (^⊤^-1^⊤)^†^⊤. By † we denote the Moore-Penrose pseudoinverse.  



  
Hessian and its eigenvalues. Note that the Hessian[While the Hessian is not self-adjoint with respect to the standard inner product, it is self-adjoint with respect to the inner product ⟨ x, y⟩ which we use as the canonical inner product in ^n.] of f = f_ is given by
∇^2 f = ^-1,
where 

    ^⊤.

 Note that ∇^2 f and 
    ^-1/2^-1/2
 have the same spectrum. Matrix  is symmetric and positive semidefinite (with respect to the standard inner product). Let 

    = Λ^⊤ = ∑_i=1^n λ_i u_i u_i^⊤

be the eigenvalue decomposition of , where = [u_1,…,u_n] is an orthonormal matrix of eigenvectors, and λ_1≤λ_2 ≤⋯≤λ_n are the corresponding eigenvalues. Let λ_min^+ be the smallest nonzero eigenvalue, and λ_max  = λ_n be the largest eigenvalue. It was shown in <cit.> that 0≤λ_i ≤ 1 for all i∈ [n].




  
Exactness. Note that f_ is a convex quadratic, and that f_(x) = 0 whenever x∈{x :  x = b}. However, f_ can be zero also for points x outside of . Clearly, f(x) is nonnegative, and f(x)=0 for x∈. However, without further assumptions, the set of minimizers of f can be larger than . The exactness assumption mentioned above ensures that this does not happen. For necessary and sufficient conditions for exactness, we refer the reader to <cit.>. Here it suffices to remark that  a sufficient condition for exactness is to require  to be positive definite. This is easy to see by observing that

    f(x) = f_(x) =12 x - b^2_.






 §.§ Three algorithms for solving the stochastic optimization problem


The authors of <cit.> consider solving the stochastic optimization problem (<ref>)  via stochastic gradient descent (SGD)[The gradient is computed with respect to the inner product ⟨ x, y⟩.]

    x_k+1 = x_k - ω∇ f__k(x_k),

where ω>0 is a fixed stepsize and  _k is sampled afresh in each iteration from . Note that the gradient of f_ with respect to the  inner product is equal to

    ∇ f_(x) (<ref>)=^-1^⊤ ( x - b)  = ^-1^⊤ (x -x_*)  = ^-1 (x -x_*),

where
^⊤, and x_* is any vector in .



They observe that, surprisingly,  SGD is in this setting equivalent to several other methods; in particular, to the stochastic Newton method[In this method we take the -pseudoinverse of the Hessian of f__k instead of the classical inverse, as the inverse does not exist. When =, the  pseudoinverse specializes to the standard Moore-Penrose pseudoinverse.],

    x_k+1 = x_k - ω (∇^2 f__k(x_k))^†_∇ f__k(x_k),

and to the stochastic proximal point method[In this case, the equivalence only works for 0<ω≤ 1.]

    x_k+1 = min_x∈^n{ f__k(x) + 1-ω/2ωx-x_k_^2}.




 §.§ Stochastic fixed point problem


The stochastic fixed point problem considered in <cit.> as one of the four stochastic reformulations has the form

    x = Π^__(x),

where the expectation is taken with respect to ∼, and where Π^__(x) is the projection of x, in the  norm, onto  the sketched system _= {x ∈^n : ^⊤ x = ^⊤ b}. An explicit formula for the projection onto  is given by

    Π_^(x)min_x' ∈x'-x_ =x-^-1^⊤ (^-1 ^ ⊤ )^† ( x-b);

a formula for _ is obtained by replacing  with ^⊤ everywhere.


The stochastic fixed point method (with relaxation parameter ω>0) for solving 
(<ref>) is defined by 
    x_k+1 = ωΠ^___k(x_k) + (1-ω) x_k.








 §.§ Best approximation problem, its dual and SDSA


It was shown in <cit.> that the above methods converge linearly to x_*=Π^_(x_0); the projection of the initial iterate onto the solution set of the linear system. Hence, besides solving problem (<ref>), they solve the best approximation problem

    min_x∈^n P(x) 12x-x_0_^2   subject to   x = b .

The Fenchel dual of (<ref>) is
the (bounded) unconstrained concave quadratic maximization problem

    max_y∈^m D(y)  (b- x_0)^⊤ y - 12^⊤ y^2_^-1.

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



In the special case of unit stepsize, method (<ref>) was first proposed by Gower and Richtárik <cit.> under the name “sketch-and-project method”, motivated by the iteration structure which proceeds in two steps: i) replace the set {x∈^n  :  x = b} by its sketched variant __k, and then project the last iterate x_k onto __k. Analysis in <cit.> was done under the assumption that  be of full column rank. This assumption was lifted in <cit.>, and a duality theory for the method developed. In particular, for ω=1, the iterates {x_k} arise as  images of the iterates {y_k} produced by  a specific dual method for solving (<ref>) under the mapping ϕ:^m ↦^n given by

    ϕ (y)  x_0 + ^-1^⊤ y.

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

    y_k+1 = y_k + _k λ_k,

where _k is in each iteration sampled from , and λ_k is chosen greedily, maximizing the dual objective D: λ_k ∈max_λ D(y_k + _k λ). Such a λ might not be unique, however. SDSA is defined by picking the solution with the smallest (standard Euclidean) norm. This leads to the formula:

    λ_k =  (_k^⊤^-1^⊤_k )^†_k^⊤(b-(x_0 + ^-1^⊤ y_k) ).


SDSA proceeds by moving in random subspaces spanned by the random columns of _k. In the special case when ω=1 and  y_0=0,  Gower and Richtárik <cit.> established the following relationship between the iterates {x_k} produced by the primal methods (<ref>), (<ref>),  (<ref>),  (<ref>) (which are equivalent), and the dual method (<ref>):

    x_k = ϕ(y_k) (<ref>)=  x_0 + ^-1^⊤ y_k.




 §.§ Other related work
 Variants of the sketch-and-project methods have been recently proposed for solving several other problems. Xiang and Zhang <cit.> 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 <cit.> 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. <cit.>. Tu et al. <cit.>  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 .  Loizou and Richtárik <cit.> use the sketch-and-project method to solve the average consensus problem; and  Hanzely et al. <cit.> design new variants of sketch and project methods for the average consensus problem with  privacy considerations (see Section <ref> for more details regarding the average consensus problem).  




§ PRIMAL METHODS WITH MOMENTUM
 

 Applied to problem (<ref>), i.e., 
min_x∈^n f(x) = f_(x),
  the gradient descent method with momentum (also known as the  heavy ball method) of Polyak <cit.> takes the form

    x_k+1 = x_k - ω∇ f(x_k) + β(x_k - x_k-1),

where ω>0 is a stepsize and β≥ 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:

    x_k+1 = x_k - ω∇ f__k(x_k) + β(x_k-x_k-1).


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:





[!h][][b]0,9
mSGD / mSN / mSPP



Parameters: 
Distribution 𝒟 from which method samples matrices; positive definite matrix ∈^n× n; stepsize/relaxation parameter ω∈ the heavy ball/momentum parameter β.

Initialize: Choose initial points x_0,x_1∈^n



For k≥ 1 do


  * Draw a fresh _k ∼

  * Set 
    x_k+1=x_k -ω∇ f__k(x_k) + β(x_k - x_k-1)


Output: last iterate x_k





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 (<ref>), (<ref>) and (<ref>), 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. 



 §.§ L2 convergence and function values: linear rate

In this section we study L2 convergence of mSGD/mSN/mSPP; that is, we study the convergence of the quantity [x_k-x_*_^2] to zero. We show that for a range  of stepsize parameters ω > 0 and  momentum terms β≥ 0 the method enjoys global  linear convergence rate. To the best of our knowledge, these results are the first of their kind for the stochastic heavy ball method.
As a corollary of  L2 convergence, we  obtain convergence of the expected function values.









Choose x_0= x_1∈^n.  Assume exactness. Let {x_k}_k=0^∞ be the sequence of random iterates produced by mSGD/mSN/mSPP.  Assume 0< ω < 2 and β≥ 0 and that the expressions

    a_1  1+3β+2β^2 - (ω(2-ω) +ωβ)λ_min^+,     and    
    a_2 β +2β^2 + ωβλ_max

satisfy a_1+a_2<1. Let x_* = Π_ℒ^(x_0). Then 

    [x_k-x_*^2_] ≤ q^k (1+δ)  x_0-x_*^2_

and 

    [f(x_k)] ≤ q^k  λ_max/2 (1+δ) x_0-x_*^2_,

where  q=a_1+√(a_1^2+4a_2)/2 and δ=q-a_1. Moreover, a_1+a_2 ≤ q <1.

 See Appendix <ref>.




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 ω∈ (0,2), then the condition a_1+a_2<1 is satisfied for all   
    0≤β< 18( -4+ωλ_min^+-ωλ_max +√((4-ωλ_min^++ωλ_max)^2+16ω (2-ω) λ_min^+ )).




If β=0, mSGD reduces to  SGD analyzed in <cit.>. In this special case, q = 1-ω(2-ω)λ_min^+, which is the rate established in <cit.>. Hence, our result is more general. 


Let q(β) be the rate as a function of β. Note that since β≥ 0, we have

    q(β)    ≥    a_1 + a_2 
       =    1 + 4β + 4β^2 + ωβ(λ_max-λ_min^+) - ω(2-ω)λ_min^+ 
       ≥    1-ω(2-ω)λ_min^+ = q(0).

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








Our next theorem states that Π_^(x_k) = x_* for all iterations k of mSGD. This invariance is important, as it allows the algorithm to converge to x_*.


Let x_0=x_1∈^n be the starting points of the mSGD method and let {x_k} be the random iterates generated by mSGD. Then Π_^(x_k)=Π_^(x_0) for all k≥ 0.

 Note that in view of (<ref>), ∇ f_(x) = ^-1^⊤ ( x - b) ∈ Range(^-1^⊤). Since 

    x_k+1 = x_k -ω∇ f__k(x_k) + β(x_k - x_k-1),

and since x_0=x_1, it can shown by induction that x_k ∈ x_0 +  Range(^-1^⊤) for all k. However, Range(^-1^⊤) is the orthogonal complement to Null() in the -inner product. Since  is parallel to Null(), vectors x_k must have the same -projection onto  for all k: Π^_(x_0) = x_*.









 §.§ 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, <cit.> and <cit.> get O(1/√(k)) and O(1/k^β) convergence, respectively. When β=1, <cit.> gets O(1/log(k)) rate.



Choose x_0=x_1 and let {x_k}_k=0^∞ be the random iterates produced by mSGD/mSN/mSPP, where the momentum parameter 0≤β <1 and relaxation parameter (stepsize) ω > 0 satisfy ω + 2β <2. Let x_* be any vector satisfying f(x_*)=0. If we let x̂_k=1/k∑_t=1^kx_t, then

    [f(x̂_k)] ≤(1-β)^2x_0-x_*_^2 + 2ωβ f(x_0)/2ω(2-2β-ω) k.



 See Appendix <ref>.


In the special case of β=0, the above theorem gives the rate

    [f(x̂_k)] ≤x_0-x_*^2_/2ω (2 -ω) k .

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



Our proof strategy  is similar to  <cit.> 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 <cit.> the authors analyzed  mSGD for general Lipshitz continuous convex objective functions (with bounded variance) and proved the sublinear rate O(1/√(k)). In <cit.>, a complexity analysis is provided for the case of quadratic strongly convex smooth coercive functions. A  sublinear convergence  rate of O(1/k^β), where β∈ (0,1), was proved. In contrast to our results, where we assume fixed stepsize ω, both papers analyze mSGD with diminishing stepsizes. 






 §.§ L1 convergence: accelerated linear rate


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



Assume exactness. Let {x_k}_k=0^∞ be the sequence of random iterates produced by mSGD / mSN / mSPP, started with x_0, x_1 ∈^n satisfying the relation x_0-x_1 ∈ Range(^-1^ ⊤), with relaxation parameter (stepsize)  0<ω≤1/λ_max and momentum parameter  (1-√(ωλ_min^+))^2 < β <1. Let x_* = Π^_(x_0). Then there exists constant C >0 such that for all k≥0 we have 

    [x_k -x_*]_^2  ≤β^k C.




  (i) If we choose ω= 1 and β= (1- √(0.99 λ_min^+)) ^2 then
 [x_k -x_*]_^2  ≤β^k C
and the iteration complexity becomes
Õ(√(1/ λ_min^+)).

  (ii) If we choose ω= 1/λ_max and β= (1- √(0.99λ_min^+/λ_max))^2 then
[x_k -x_*]_^2  ≤β^k C
and the iteration complexity becomes
Õ(√(λ_max/ λ_min^+)).



 See Appendix <ref>.


Note that the convergence factor is precisely equal to the value of the momentum parameter β. Let x be any random vector in ^n with finite mean [x], and x_*∈^n is any reference vector (for instance, any solution of x = b). Then we have the identity (see, for instance <cit.>)

    x-x_*_^2 = x-x_*_^2 + x-[x]^2_.

This means that the quantity x-x_*_^2 appearing in our L2 convergence result (Theorem <ref>) is larger than 
x-x_*_^2 appearing in the L1 convergence result (Theorem <ref>), and hence harder to push to zero. As a corollary, L2 convergence implies L1 convergence. However, note that in Theorem <ref> we have established an accelerated rate.  A similar theorem, also obtaining an accelerated rate in the L1 sense, was established in <cit.> for an accelerated variant of SGD in the sense of Nesterov.




§ DUAL METHODS WITH MOMENTUM
 

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




[!h][][b]0,9
Stochastic Dual Subspace Ascent with Momentum (mSDSA)



Parameters: 
Distribution 𝒟 from which method samples matrices; positive definite matrix ∈^n× n; stepsize/relaxation parameter ω∈ the heavy ball/momentum parameter β. SDSA is obtained as a special case of mSDSA for β=0. 

Initialize: Choose initial points y_0 =  y_1 = 0 ∈^m



For k≥ 1 do


  * Draw a fresh _k ∼

  * Set λ_k = (_k^⊤^-1^⊤_k )^†_k^⊤(b-(x_0 + ^-1^⊤ y_k) )

  * Set y_k+1=y_k+ω_k λ_k +β(y_k-y_k-1)

Output: last iterate y_k







 §.§ Correspondence between primal and dual methods


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 ϕ defined in (<ref>).


   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 ω > 0, momentum parameter β≥ 0, and the same sequence of random matrices _k.  Then  

    x_k = ϕ(y_k) = x_0 + ^-1^⊤ y_k

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


First note that

    ∇ f__k(ϕ(y_k))    (<ref>)=   ^-1^⊤_k ( _k^⊤^-1^⊤_k)^†_k^⊤ (ϕ(y_k) - b)   =   - ^-1^⊤_k λ_k.

We now use this to show that

    ϕ(y_k+1)    (<ref>)=    x_0  + ^-1^⊤ y_k+1
       =    x_0+^-1^⊤[y_k+ω_k λ_k +β(y_k-y_k-1)]
       =   x_0+^-1^⊤ y_k_ϕ(y_k)+ω^-1^⊤_k λ_k_-∇ f__k(ϕ(y_k)) +β^-1^⊤ (y_k-y_k-1)
       =   ϕ(y_k)- ω∇ f__k(ϕ(y_k)) +β  (^-1^⊤ y_k-^-1^⊤ y_k-1)
       (<ref>)=   ϕ(y_k) -ω∇ f__k(ϕ(y_k)) +β  (ϕ(y_k) - ϕ(y_k-1)).

So, the sequence of vectors {ϕ(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 = ϕ(0) =  ϕ(y_0), and x_1 = x_0= ϕ(0) =  ϕ(y_1).





 §.§ 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).

 

Choose y_0= y_1∈^n.  Assume exactness. Let {y_k}_k=0^∞ be the sequence of random iterates produced by mSDSA.  Assume 0≤ω≤ 2 and β≥ 0 and that the expressions

    a_1  1+3β+2β^2 - (ω(2-ω) +ωβ)λ_min^+,     and    
    a_2 β +2β^2 + ωβλ_max

satisfy a_1+a_2<1. Let x_* = Π_ℒ^(x_0) and let y_* be any dual optimal solution. Then 

    [D(y_*)-D(y_k)] ≤ q^k (1+δ) [D(y_*)-D(y_0)]

where  q=a_1+√(a_1^2+4a_2)/2 and δ=q-a_1. Moreover, a_1+a_2 ≤ q <1.

 This follows by applying Theorem <ref> together with Theorem <ref> and the identity 12x_k-x_0^2_ = D(y_*) - D(y_k).



Note that for β=0, mSDSA simplifies to SDSA. Also recall that for unit stepsize (ω=1), SDSA was analyzed in  <cit.>. In the ω=1 and β=0 case, our result specializes to that established in  <cit.>. Following similar arguments to those in  <cit.>,  the same rate of convergence can be proved for the duality gap [P(x_k)-D(y_k)].





§ METHODS  WITH STOCHASTIC MOMENTUM
 

To motivate stochastic momentum, for simplicity fix =, and assume that _k is chosen as the jth random unit coordinate vector  of ^m with probability p_j>0. In this case, SGD (<ref>) reduces to the randomized Kaczmarz method for solving the linear system x = b, first analyzed for p_j∼_j:^2 by Strohmer and Vershynin <cit.>. 

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

    x_k+1 = x_k - ω_j: x_k -b_j/_j:^2_j:^ ⊤ + β(x_k-x_k-1).

Note that the cost of one iteration of this method is (_j:_0 + n), where the cardinality term _j:_0 comes from the stochastic gradient part, and n comes from the momentum part.
When  is sparse, the second term will dominate. Similar considerations apply for many other (but clearly not all) distributions . 

In such circumstances, we propose to replace the expensive-to-compute momentum term by a cheap-to-compute stochastic approximation thereof. In particular, we let i_k be chosen from [n] uniformly at random,  and replace x_k-x_k-1 with v_i_k e_i_k^⊤ (x_k-x_k-1)e_i_k^⊤,  where e_i_k∈^n is the i_kth unit basis vector in ^n, and β with nβ. Note that v_i_k can be computed in (1) time. Moreover, 
    _i_k[nβ v_i_k] =β( x_k-x_k-1).
 Hence,  we replace the momentum term by  an unbiased estimator, which allows us to cut the cost to (_j:_0). 





 §.§ Primal methods with stochastic momentum


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:

    x_k+1 = x_k - ω∇ f__k(x_k) + nβ e_i_k^⊤ (x_k-x_k-1)e_i_k.


The method is formalize below:




[!h][][b]0,9
smSGD/smSN/smSPP



Parameters: 
Distribution 𝒟 from which the method samples matrices;  stepsize/relaxation parameter ω∈ the heavy ball/momentum parameter β.

Initialize: Choose initial points x_1=x_0  ∈^n; set = ∈^n× n



For k≥ 1 do


  * Draw a fresh _k ∼

  * Pick i_k∈ [n] uniformly at random 

  * Set 
    x_k+1=x_k -ω∇ f__k(x_k) + β e_i_k^⊤ (x_k - x_k-1)e_i_k


Output: last iterate x_k










 §.§ Convergence


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

 Choose x_0= x_1∈^n.  Assume exactness. Let =.  Let {x_k}_k=0^∞ be the sequence of random iterates produced by smSGD/smSN/smSPP.  Assume 0< ω < 2 and β≥ 0 and that the expressions

    a_1  1+3βn+2β^2n - (ω(2-ω) +ωβn)λ_min^+,     and    
    a_2 1n(β +2β^2 + ωβλ_max)

satisfy a_1+a_2<1. Let x_* = Π_ℒ^(x_0). Then 

    [x_k+1-x_*^2] ≤ q^k (1+δ)  x_0-x_*^2

and 
[f(x_k)] ≤ q^k  λ_max/2 (1+δ) x_0-x_*^2,
where  qa_1+√(a_1^2+4a_2)/2 and δ q-a_1. Moreover, a_1+a_2 ≤ q < 1.


 See Appendix <ref>.


It is straightforward to see that if we choose ω∈ (0,2), then the condition a_1+a_2<1 is satisfied for all β belonging to the interval   
    0≤β < 18 (-4+ωλ_min^+-ωλ_max+√((4-ωλ_min^++ωλ_max)^2+16nω (2-ω)  λ_min^+)).

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.






 §.§ Momentum versus stochastic momentum




As indicated in the introduction, if we wish to compare mSGD to smSGD used with momentum parameter β, it makes sense to use momentum parameter β 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(β) be the convergence constant for mSGD with stepsize ω=1 and an admissible momentum parameter β≥ 0. Further, let a̅_1(β),a̅_2(β),q̅(β) be the convergence constants for smSGD with stepsize ω=1 and momentum parameter β.  We have

    q̅(β n) ≥a̅_1(β n) + a̅_2(β n)    (<ref>)=    1 + 4β + 4β^2n + β(λ_max-λ_min^+) - λ_min^+
       (<ref>)=    a_1(β) + a_2(β) + 4β^2 (n-1)
       ≥    a_1(β) + a_2(β).

 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 a̅_1(β n) - a_1(β) = 2β^2(n-1) ≥ 0 and a̅_2(β n) - a_2(β) = 2β^2 (n-1) ≥ 0, we have

    q̅(β n) = a̅_1(β n) + √(a̅^2_1(β n) +4 a̅_2(β n))/2≥a_1(β ) + √(a_1^2(β ) +4 a_2(β ))/2 = q(β),

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 (g+n) (mSGD) against (g) (smSGD). If g≪ 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:

    C_mSGD(β) g+n/1-q(β) = g+n/1-a_1(β) + √(a_1^2(β) + 4 a_2(β))/2,

and

    C_smSGD(β n) g/1-q̅(β n) = g/1-a̅_1(β n) + √(a̅_1^2(β n) + 4 a̅_2(β n))/2.


Since  
    q(0) = q̅(0 n),
 and since q(β) and q̅(β n) are continuous functions of β, then because g+n > g, for small enough β we will have 
C_mSGD(β) > C_smSGD(β n). In particular, the speedup of smSGD compared to mSGD for β≈ 0 will be close to

    C_mSGD(β)/C_smSGD(β n)≈lim_β' →_+ 0C_mSGD(β')/C_smSGD(β' n)(<ref>)+(<ref>)+(<ref>)=g+n/g = 1 + n/g.


Thus, we have shown the following statement.

  For small β, 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 , and g is the expected cost of evaluating a stochastic gradient ∇ f_(x).  






§ SPECIAL CASES: RANDOMIZED KACZMARZ WITH MOMENTUM AND RANDOMIZED COORDINATE DESCENT WITH MOMENTUM
 


In Table <ref> we specify several special instances of mSGD  by choosing distinct combinations of the parameters  and .  We use e_i to denote the ith unit coordinate vector in ^m, and _:C for the column submatrix of the m× m identity matrix indexed by (a random) set C. 






The updates for smSGD can be derived by substituting the momentum term β(x_k-x_k-1) with its stochastic variant nβ e_i_k^⊤ (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  and   we refer the interested reader to Section 3 of <cit.>. 

In the rest of this section we present in detail two special cases: the randomized Kaczmarz method with momentum (mRK) and the randomized coordinate descent method with momentum (mRCD). Further, we compare the L1 convergence rates (i.e., bounds on x_k-x_*_^2) obtained in this paper  with rates that can be inferred from known results for their no-momentum variants.



 §.§ mRK: randomized Kaczmarz with momentum
 We now provide a discussion on mRCD (the method in the first row of Table <ref>).  Let = and let pick in each iteration the random matrix =e_i with probability p_i=_i:^2 / _F^2. In this setup the update rule of the mSGD simplifies to 
    x_k+1=x_k -ω_i: x_k -b_i/_i:_2^2_i:^ ⊤ + β(x_k - x_k-1)

and 

    (<ref>)=   ^-1/2^⊤_∼[] ^-1/2  =  [ ^⊤] 
       =   ∑_i=1^mp_i_i:^⊤_i:/_i:^2_2  =  1/^2_F∑_i=1^m_i:^⊤_i:  = ^⊤/^2_F.


The objective function takes the following form:

    f(x)=_∼𝒟 [f_(x)]= ∑_i=1^m p_i f__i(x)= x -b^2_2/2^2_F.


For β=0, this method reduces to the randomized Kaczmarz method with relaxation, first analyzed in <cit.>. If we also have ω=1, this is equivalent with the randomized Kaczmarz method of Strohmer and Vershynin <cit.>.  RK without momentum (β=0)  and without relaxation (ω=1) 
converges  with iteration complexity <cit.> of 
    Õ(1/λ_min^+())=Õ(^2_F/λ_min^+(^⊤)).


In contrast, based on Theorem <ref> we have


  * For ω= 1 and β= (1- √(0.99 λ_min^+))^2=(1- √(0.99/^2_Fλ_min^+(^⊤)))^2, the iteration complexity of the mRK is: 
    Õ(√(^2_F/λ_min^+(^⊤))).


  * For ω= ^2_F/λ_max(^⊤) and β= (1- √(0.99 λ_min^+(^⊤)/λ_max(^⊤)))^2
the iteration complexity becomes:

    Õ(√(λ_max(^⊤)/λ_min^+(^⊤))).



 This is quadratic improvement on the  previous best result (<ref>).




  
Related Work. 
The Kaczmarz method for solving consistent linear systems was originally introduced by Kaczmarz in 1937 <cit.>. 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 <cit.>. In this work we are interested in randomized variants of the Kaczmarz method, first analyzed by Strohmer and Vershynin <cit.>. In <cit.> 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 <cit.>.  All existing results on accelerated variants of RK use the Nesterov's approach of acceleration <cit.>. To the best of our knowledge, no convergence analysis of mRK exists in the literature (Polyak's momentum). Our work fills this gap.




 §.§ mRCD: randomized coordinate descent with momentum
 

We now provide a discussion on the mRCD method (the method in the second row of Table <ref>). If the matrix  is positive definite, then we can choose =  and =e_i with probability p_i=_ii/ Trace(). It is easy to see that =/ Trace(). In this case,   is positive definite and as a result, λ_min^+()=λ_min(). Moreover, we have

    f(x)=_∼𝒟 [f_(x)]= ∑_i=1^m p_i f__i(x)= x -b^2_2/2 Trace().


For β=0 and ω=1 the method is equivalent with  randomized coordinate descent of Leventhal and Lewis <cit.>, which was shown to converge with iteration complexity  
    Previous  best result:    Õ( Trace()/λ_min()).


In contrast, following  Theorem <ref>, we can obtain the following  L_1 iteration complexity results for mRCD:
 
  
  * For ω= 1 and β= (1- √(0.99/ Trace()λ_min()))^2, the iteration complexity is 
    Õ(√( Trace()/λ_min())).


  * For  ω=  Trace()/λ_max() and β= (1- √(0.99 λ_min()/λ_max()))^2  the iteration complexity becomes 

    Õ(√(λ_max()/λ_min())).

 
 
 This is quadratic improvement on the previous best result (<ref>).
 
 


  
Related Work.  It is known that if  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 <cit.>, and later extended by several authors to more general settings, including smooth convex optimization <cit.>, composite convex optimization <cit.>, and parallel/subspace descent variants  <cit.>. These results were later further extended to handle  arbitrary sampling distributions  <cit.> . Accelerated variants of RCD were studied in <cit.>. For other non-randomized coordinate descent variants and their convergence analysis, we refer the reader to <cit.>. To the best of our knowledge, mRCD and smRCD have never  been analyzed before in any setting.



 §.§ 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.



In Figure <ref> we present in a simple ^2 illustration of the difference between the workings of RK and mRK. Our goal is to show graphically how the addition of momentum leads to acceleration.  Given iterate x_k, one can think of the update rule of the mRK (<ref>) in two steps:


  * The Projection: The projection step corresponds to the first part x_k - ω∇ f__k(x_k) of the  mRK update (<ref>) and it means that the current iterate x_k is projected onto a randomly chosen hyperplane H_i[In the plots of Figure  <ref>, the hyperplane of each update is chosen in an alternating fashion for illustration purposes]. The value of the stepsize ω∈ (0,2) defines whether the projection is exact or not. When ω=1 (no relaxation) the projection is exact, that is the point  Π_H_i(x_k) belongs in the hyperplane H_i.  In Figure <ref> all projections are exact.

  * Addition of the momentum term: The momentum term (right part of the update rule) β (x_k-x_k-1) forces the next iterate x_k+1 to be closer to the solution x_* than the corresponding point Π_H_i(x_k). Note also that the vector x_k+1-Π_H_i(x_k) is always parallel to x_k-x_k-1 for all k≥0.



In the example of Figure <ref>, 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 this section we study the computational behavior of the two proposed algorithms, mSGD and smSGD. In particular, we focus mostly on the evaluation of the performance of mSGD. To highlight the usefulness of smSGD, an empitical verification of Theorem <ref> is presented in subsection <ref>. As we have already mentioned, both mSGD and smSGD can be interpreted as sketch-and-project methods (with relaxation), and as a result a comprehensive array of well-known algorithms can be recovered as special cases by varying the main parameters of the methods (check Section <ref>). In our experiments we focus on the popular special cases of randomized Kaczmarz method (RK) and the randomized coordinate descent method (RCD) without relaxation (ω=1), and  show the practical benefits of the addition of the momentum term[The experiments were repeated with various values of the main parameters and initializations, and similar results were obtained in all cases.]. The choice of the stepsize ω=1 is not arbitrary. Recently, in <cit.> both relaxed RK and relaxed RCD were analyzed, and it was proved that the quantity x_k-x_*^2_ converges linearly to zero for ω∈ (0,2), and that the best convergence rate is obtained precisely for ω=1. Thus the comparison is with the best-in-theory  no-momentum variants. 


Note that, convergence analysis of the error x_k-x_*^2_ (L2 convergence) and of the expected function values f(x_k) in Theorem <ref> shows that mSGD enjoys global non-asymptotic linear convergence rate but not faster than the no-momentum method. 
The accelerated linear convergence rate has been obtained only in the weak sense (Theorem <ref>). Nevertheless, in practice as indicated from our experiments, mSGD is faster than its no momentum variant. Note also that in all of the presented experiments the momentum parameters β of the methods are chosen to be positive constants that do not depend on parameters that are not known to the users such as λ_min^+ and λ_max.

In comparing the methods with their momentum variants we use both the relative error measure x_k-x_*^2_ / x_0-x_*^2_ 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 (<ref>) in the case of RK and to function (<ref>) for the RCD.  For block variants the objective function of problem (<ref>) has also closed form expression but it can be very difficult to compute. In these cases one can instead evaluate the quantity x-b^2_.]. In all implementations, except for the experiments on average consensus (Section <ref>), 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 <ref> (comparison between the mSGD and smSGD methods). In the last experiment building upon the recent results of <cit.> we show how the addition of the momentum accelerates the pairwise randomized gossip (PRG) algorithm for solving the average consensus problem.





 §.§ Evaluation of mSGD


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.



  §.§.§ Synthetic Data


The synthetic data for this comparison is generated as follows[Note 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_*=Π_^(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_ / x_0-x_*^2_, one should be careful and use the value x_*=x_0+^† (b-  x_0)x_0=0=^† b.].

For mRK:
All elements of matrix ∈^m × n and of vector z ∈^n are chosen to be i.i.d 𝒩(0,1). Then the right hand side of the linear system is set to b= z. With this way the consistency of the linear system with matrix  and right hand side b is ensured. 

For mRCD:  A Gaussian matrix ∈^m × n is generated and then matrix = ^⊤∈^n × n is used in the linear system. The vector z∈^n is chosen to be i.i.d 𝒩(0,1) and again to ensure consistency of the linear system, the right hand side is set to b= z.

In particular for the evaluation of mRK we generate Gaussian matrices with m=300 rows and several columns while for the case of mRCD the matrix  is chosen to be Gaussian with m=500 rows and several columns[RCD converge to the optimal solution only in the case of positive definite matrices. For this reason = ^⊤∈^n × n is used which with probability 1 is a full rank matrix]. 
Linear systems of these forms were extensively studied <cit.> and it was shown that the quantity 1/λ_min^+(condition number) can be easily controlled. 


For each linear system we run mRK (Figure <ref>)  and mRCD (Figure <ref>) for several values of momentum parameters β and fixed stepsize ω=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 β=0 the methods are equivalent with their no-momentum variants RK and RCD respectively.

From Figures <ref> and <ref> 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/λ_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 (β=0) converges slowly, that is when  1/λ_min^+ is large (ill-conditioned linear systems). 

  * For the case of fixed stepsize ω=1, the problems with small condition number require smaller momentum parameter β to have faster convergence. Note the first two rows of Figures <ref> and <ref>, where β=0.3 or β =0.4, are good options. 

  * For large values of 1/λ_min^+, it seems that the choice of β=0.5 is the best. As an example for matrix ∈^300 × 280 in Figure <ref>, (where 1/λ_min^+=208,730), note that to reach relative error 10^-10, RK needs around 2 million iterations, while mRK with momentum parameter β=0.5 requires only half that many iterations. The acceleration is obvious also in terms of time where in 12 seconds the mRK with momentum parameter β=0.5 achieves relative error of the order 10^-9 and RK requires more than 25 seconds to obtain the same accuracy. 

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

  * In <cit.> 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.








  §.§.§ Real Data

In the following experiments we test the performance of mRK using real matrices (datasets) from the library of support vector machine problems LIBSVM <cit.>. Each dataset consists of a matrix ∈^m × n  (m features and n characteristics) and a vector of labels b ∈^m.  In our experiments we choose to use only the matrices of the datasets and ignore the label vector. As before, to ensure consistency of the linear system,  we choose a Gaussian vector z∈^n and the right hand side of the linear system is set to b= z. 
Similarly as in  the case of synthetic data, mRK is tested for several values of momentum parameters β and fixed stepsize ω=1. 

In Figure <ref> the performance of all methods for both relative error measure x_k-x_*^2 / x_*^2_ and function values f(x_k) is presented. Note again that β=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 <ref>, where mRK with β=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.








 §.§ Comparison of momentum & stochastic momentum


In Theorem <ref>, 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 β,

    C_β=C_mSGD(β)/C_smSGD(β n)≈  1 + n/g,
  where C_mSGD and C_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  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 <ref>.

As a first step a Gaussian matrix ∈^m × n is generated.
Then using this matrix several consistent linear systems are obtained as follows. Several values for g ∈ [1,n] are chosen and for each one of these a matrix _g ∈^m × n with the same elements as  but with n-g zero coordinates per row is produced. 
For every matrix _g, a Gaussian vector z_g ∈^n is drawn and to ensure consistency of the linear system, the right hand side is set to b_g=_g z.

We run both mSGD and smSGD with small momentum parameter β=0.0001 for solving the linear systems _g x=b_g for all selected values of g ∈ [1,n]. The starting point for each run is taken to be x_0=0 ∈^n.  The methods run until ϵ=x_k-x_*<0.001, where x_*=Π__g(x_0) and _g is the solution set of the linear system _g x=b_g. In each run the number of operations needed to achieve the accuracy ϵ have been counted. 
For each linear system the average after 10 trials of the value C_mSGD(β)/C_smSGD(β n) is computed. 



In Figure <ref> the actual ratio C_mSGD(β)/C_smSGD(β n) and the theoretical approximation 1 + n/g are plot and it is shown that they have similar behavior. Thus the theoretical prediction of Theorem <ref> is numerically confirmed.
In particular in the implementations we use the Gaussian matrices ∈^200 × 100 and ∈^1000 × 300.






 §.§ Faster method for average consensus





  §.§.§ Background

Average consensus (AC) is a fundamental problem in distributed computing and multi-agent systems <cit.>. Consider a connected undirected network =(,) with node set ={1,2,…,n} and edges , (||=m), where each node i ∈ owns a private value c_i ∈. The goal of the AC problem is each node of the network to compute the average of these private values, c̅1n∑_i c_i, via a protocol which allows communication between neighbours only. The problem comes up in many real world applications such as coordination of autonomous agents, estimation, rumour spreading in social networks, PageRank and distributed data fusion on ad-hoc networks and decentralized optimization. 

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 <cit.>. As we have already explained both basic method <cit.> 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_*=Π^_(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 (<ref>).
In the special case that


  * the linear system in the constraints of (<ref>) is the homogeneous linear system (x=0) with matrix ∈^ m × n being the incidence matrix of the undirected graph =(,), and 

  * 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 c̅1n∑_i c_i. Under this setting, the famous randomized pairwise gossip algorithm (randomly pick an edge e∈ E and replace the private values of its two nodes to their average) that was first proposed and analyzed in <cit.>, is equivalent with the RK method without relaxation (ω=1) <cit.>.

 


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 network[Matrix  of the linear system is the incidence matrix of the graph and it is known that the Laplacian matrix is equal to =^⊤, where ^2_F=2m.]:

    1/λ_min^+()(<ref>)=1/λ_min^+(^⊤/ ^2_F)^2_F=2m=2m/λ_min^+(^⊤)=2m/λ_min^+(),

where =^⊤ is the Laplacian matrix of the network and the quantity λ_min^+() is the very well studied algebraic connectivity of the graph <cit.>. 



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




  §.§.§ 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 <cit.>. 
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 β for 10 trials and we plot their average. Our results are available in Figures <ref> and <ref>.

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 <ref> and <ref> the initial value of each node is chosen independently at random from the uniform distribution in the interval (0,1).




  §.§.§ Experimental Results


By observing Figures <ref> and <ref>, it is clear that the addition of the momentum term improves the performance of the popular pairwise randomized gossip (PRG) method <cit.>.  The choice β=0.4 as the momentum parameter improves the performance of the vanilla PRG for all networks under study and β=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 <cit.>. Thus, using equation (<ref>), we can compute the exact values of the condition number 1/λ_min^+ for these networks. Interestingly, as we can see in Table <ref> for n=100 and n=200 (number of nodes), the condition number 1/λ_min^+ appearing in the iteration complexity of our methods is not very large. This is in contrast with experimental  observations from Section <ref> where it was shown that the choice β=0.5 is good  for very ill conditioned problems only (1/λ_min^+ very large). 



 
















§ PROOF OF THEOREM <REF>
 



 §.§ Lemmas


We start with a lemma.



Fix F_1=F_0≥ 0 and let {F_k}_k≥ 0 be a  sequence of nonnegative  real numbers satisfying the relation 

    F_k+1≤ a_1F_k +a_2 F_k-1,     ∀ k≥ 1,
 where  a_2 ≥ 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≤ q^k (1+ δ)  F_0 for all k≥ 1,
where q=a_1+√(a_1^2+4a_2)/2  and δ=q-a_1≥ 0. Moreover,

    q ≥ a_1 + a_2,

with equality if and only if a_2=0 (in which case q=a_1 and δ=0).



Choose any δ≥ 0 satisfying a_2 ≤ (a_1+δ)δ. Adding δ F_k  to both sides of (<ref>),  we get 

    F_k+1 + δ F_k  ≤   (a_1+δ)F_k + a_2 F_k-1
    ≤ (a_1+δ)(F_k + δ F_k-1) = q(F_k+δ F_k-1).



We now claim that   δ = -a_1+√(a_1^2+4a_2)/2 satisfies the relations. Non-negativity of δ  follows from a_2≥ 0, while the second   relation follows from the fact that δ satisfies 
    (a_1+δ)δ - a_2 = 0.
 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 (<ref>), we obtain
F_k+1≤ F_k+1 + δ F_k ≤ q^k (F_1+ δ F_0) = q^k(1+δ) F_0.

Finally, let us establish (<ref>). Noting that a_1 = q-δ, and since  in view of (<ref>) we have a_2=qδ, we conclude that a_1+a_2 = q + δ(q-1) ≤ q, where the inequality follows from q <1.



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


For all x ∈^n we have

    f_(x) = 1/2∇ f_(x)^2_.

Moreover, if x_*∈ (i.e., if x_* satisfies x_* =b), then for all x∈^n we have

    f_(x) = 1/2⟨∇ f_(x),x-x_* ⟩_,

and

    f(x) = 1/2⟨∇ f(x),x-x_* ⟩_.



 In view of (<ref>), and  since ^-1 = (see <cit.>), we have


    ∇ f_(x)^2_   (<ref>)=   ^-1 (x-x_*)^2_   =    (x-x_*)^⊤^-1 (x-x_*)    =    (x-x_*)^⊤ (x-x_*) 
       (<ref>)=   (x-x_*)^⊤^⊤ (x-x_*)    =    ( x-b)^⊤ (  x- b )   (<ref>)=   2f_(x).


Moreover,

    ⟨∇ f_(x),x-x_* ⟩_   (<ref>)=   ⟨^-1 (x-x_*),x-x_* ⟩_
        =    (x-x_*)^⊤^-1 (x-x_*)    =    2f_(x).

By taking expectations in the last identity with respect to the random matrix , we get
⟨∇ f(x),x-x_* ⟩_=2f(x).






For all x ∈^n and x_* ∈

    λ_min^+ f(x) ≤1/2∇ f(x) ^2_≤λ_max f(x)

and 

    f(x) ≤λ_max/2x-x_*^2_.

Moreover, if exactness is satisfied, and we let x_* =Π^_(x), we have 

    λ_min^+/2x-x_*^2_≤  f(x)  .








 §.§ The Proof

First, we decompose

    x_k+1-x_*^2_    =   x_k-ω∇ f__k(x_k)+β(x_k-x_k-1)-x_*^2_
        =   x_k-ω∇ f__k(x_k)-x_*^2__1
             +2⟨ x_k-ω∇ f__k(x_k)-x_*, β(x_k-x_k-1)  ⟩__2
             +β^2x_k-x_k-1^2__3.

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

    1   =   x_k-x_*^2_-2 ω⟨ x_k-x_*,∇ f__k(x_k) ⟩_ +ω^2 ∇ f__k(x_k)^2_
       (<ref>),(<ref>)=   x_k-x_*^2_-4ω f__k(x_k)+2ω^2 f__k(x_k) 
       =   x_k-x_*^2_-2ω(2-ω)f__k(x_k).

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

    2   = 2β⟨ x_k-x_*, x_k-x_k-1⟩_ +2ωβ⟨∇ f__k(x_k), x_k-1 - x_k  ⟩_
       = 2β⟨ x_k-x_*, x_k-x_*⟩_ + 2β⟨ x_k-x_*, x_*-x_k-1⟩_ +2ωβ⟨∇ f__k(x_k),x_k-1- x_k  ⟩_
       = 2βx_k-x_*^2_  + 2β⟨ x_k-x_*, x_*-x_k-1⟩_ +2ωβ⟨∇ f__k(x_k),x_k-1- x_k  ⟩_.

Using the fact that for arbitrary vectors a,b,c ∈^n we have the identity
2 ⟨ a-c,c-b ⟩_=a-b^2_-c-b^2_-a-c^2_,
we obtain

    2 ⟨ x_k-x_*, x_*-x_k-1⟩_=  x_k-x_k-1^2_- x_k-1-x_*^2_-x_k-x_*^2_.

Substituting this into (<ref>) gives

    2    = βx_k-x_*^2_+βx_k-x_k-1^2_-βx_k-1-x_*^2_ + 2ωβ⟨∇ f__k(x_k),x_k-1- x_k  ⟩_.


The third expression can be bound as

    3 =β^2(x_k-x_*)+(x_*-x_k-1)^2_≤  2β^2x_k-x_*^2_+2
    β^2x_k-1-x_*^2_.


By substituting the  bounds (<ref>), (<ref>), (<ref>) into  (<ref>) we obtain

    x_k+1-x_*^2_   ≤   x_k-x_*^2_-2ω(2-ω)f__k(x_k)
             + βx_k-x_*^2_+βx_k-x_k-1^2_-βx_k-1-x_*^2_
             + 2ωβ⟨∇ f__k(x_k),x_k-1- x_k  ⟩_  +  2β^2x_k-x_*^2_+2β^2x_k-1-x_*^2_
       ≤    (1+3β + 2β^2)x_k-x_*^2_+ (β + 2β^2 )x_k-1-x_*^2_-2ω(2-ω)f__k(x_k)
             + 2ωβ⟨∇ f__k(x_k),x_k-1- x_k  ⟩_.

Now by first taking expectation with respect to _k, we obtain:

    __k[x_k+1-x_*^2_]    ≤    (1+3β+2β^2)x_k-x_*^2_+ (β +2β^2)x_k-1-x_*^2_
             -2ω(2-ω)f(x_k) + 2ωβ⟨∇ f(x_k),x_k-1- x_k  ⟩_
       ≤    (1+3β+2β^2)x_k-x_*^2_+ (β +2β^2)x_k-1-x_*^2_
             -2ω(2-ω)f(x_k) + 2ωβ(f(x_k-1)-f(x_k))
        =     (1+3β+2β^2)x_k-x_*^2_+ (β +2β^2)x_k-1-x_*^2_
             - (2ω(2-ω) +2ωβ)f(x_k) + 2ωβ f(x_k-1).

where in the second step we used the inequality
 ⟨∇ f(x_k),x_k-1- x_k ⟩≤ f(x_k-1)-f(x_k) and the fact that ωβ≥ 0, which follows from the assumptions. We now apply  inequalities  (<ref>) and (<ref>), obtaining


    __k[x_k+1-x_*^2_]    ≤   (1+3β+2β^2 - (ω(2-ω) +ωβ)λ_min^+)_a_1x_k-x_*^2_
             + (β +2β^2 + ωβλ_max)_a_2x_k-1-x_*^2_.




By taking expectation again, and letting F_k[x_k-x_*^2_], we get  the relation

    F_k+1≤ a_1 F_k  + a_2 F_k-1 .


It suffices to apply  Lemma <ref> to the relation  (<ref>). The conditions of the lemma are satisfied. Indeed, a_2≥ 0, and if a_2=0, then β=0 and hence a_1=1-ω(2-ω)λ_min^+>0. The condition a_1+a_2<1 holds by assumption.


The convergence result in function values, [f(x_k)], follows as a corollary by applying inequality (<ref>) to  (<ref>).




§ PROOF OF THEOREM <REF>



Let p_t=β/1-β(x_t-x_t-1) 
and d_t = x_t + p_t -x_*_^2. In view of (<ref>),  we can write 

    x_t+1+p_t+1= x_t+p_t-ω/1-β∇ f__t(x_t),

and therefore

    d_t+1    =   x_t+p_t-ω/1-β∇ f__t(x_t) -x_* ^2_ 
        =    d_t -2 ω/1-β⟨ x_t+p_t-x_*,  ∇ f__t(x_t) ⟩_ + ω^2/(1-β)^2∇ f__t(x_t)^2_
        =    d_t -2ω/1-β⟨ x_t-x_*,  ∇ f__t(x_t) ⟩_ - 2 ωβ/(1-β)^2⟨ x_t-x_t-1,  ∇ f__t(x_t) ⟩_
             + ω^2/(1-β)^2∇ f__t(x_t)^2_.




Taking expectation with respect to the random matrix _t we obtain:

    __t[d_t+1]     =   __t[d_t] -2ω/1-β⟨ x_t-x_*,  ∇ f(x_t) ⟩_ - 2 ωβ/(1-β)^2⟨ x_t-x_t-1,  ∇ f(x_t) ⟩_
             + ω^2/(1-β)^2 2 f(x_t) 
       (<ref>)=   __t[d_t] -4ω/1-β  f(x_t) - 2 ωβ/(1-β)^2⟨ x_t-x_t-1,  ∇ f(x_t) ⟩_ + ω^2/(1-β)^2 2 f(x_t)
       ≤   __t[d_t] -4ω/1-β  f(x_t) - 2 ωβ/(1-β)^2   [f(x_t)-f(x_t-1)] + ω^2/(1-β)^2 2 f(x_t)
        =    __t[d_t] + [ -4ω/1-β - 2 ωβ/(1-β)^2 +2 ω^2/(1-β)^2] f(x_t)  +  2 ωβ/(1-β)^2 f(x_t-1),

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

    __t[d_t+1] +   2 ωβ/(1-β)^2 f(x_t) + α f(x_t)  ≤__t[d_t] + 2 ωβ/(1-β)^2 f(x_t-1),

where α =   4ω/1-β -2 ω^2/(1-β)^2 > 0. Taking expectations again and using the tower property, we get

    θ_t+1 + α[f(x_t)]  ≤θ_t,      t=1,2,…,

where θ_t = [d_t] + 2 ωβ/(1-β)^2[ f(x_t-1)]. By summing up (<ref>) for t=1,…, k we get

    ∑_t=1^k [f(x_t)] ≤θ_1-θ_k-1/α≤θ_1/α.

Finally, using Jensen's inequality, we get

    [f(x̂_k)] = [f(1/k∑_t=1^k x_t)] ≤[1/k∑_t=1^k f(x_t)] =  1/k∑_t=1^k [f(x_t)] (<ref>)≤θ_1/α k.


It remains to note that θ_1 = x_0-x_*_^2 + 2ωβ/(1-β)^2 f(x_0).





§ PROOF OF THEOREM <REF>
 

In the proof of Theorem <ref> the following two lemmas are used.



Assume exactness. Let x∈^n and x_* = Π_ℒ^(x). If λ_i=0, then u_i^⊤^1/2 (x-x_*)=0.





Consider the second degree linear homogeneous recurrence relation:

    r_k+1= a_1r_k+a_2 r_k-1

with initial conditions r_0,r_1 ∈.
Assume that the constant coefficients a_1 and a_2 satisfy the inequality a_1^2 +4a_2<0 (the roots of the characteristic equation t^2-a_1t-a_2=0 are imaginary). Then there are complex constants C_0 and C_1 (depending on the initial conditions r_0 and r_1) such that:

    r_k=2 M^k (C_0 cos( θ k) + C_1 sin(θ k))

where M= (√(a_1^2/4+(-a_1^2-4a_2)/4))=√(-a_2) and θ is such that a_1=2 M cos(θ) and √(-a_1^2-4a_2)=2 M sin(θ).



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

    x_k+1    =    x_k -ω∇ f__k(x_k) + β(x_k - x_k-1) 
       (<ref>)=    x_k- ω^-1_k(x_k-x_*) + β(x_k - x_k-1).


Subtracting  x_* from both sides of (<ref>), we get

    x_k+1-x_* 
        =     ( - ω^-1_k)(x_k-x_*) + β(x_k -x_* +x_* - x_k-1)
        =   ((1+β) - ω^-1_k)(x_k-x_*) - β(x_k-1-x_*).


Multiplying the last identity from the left by ^1/2, we get

    ^1/2 (x_k+1-x_*)
       =   ((1+β) - ω^-1/2_k ^-1/2) ^1/2(x_k -x_*) - β^1/2(x_k-1-x_*).


Taking expectations, conditioned on x_k (that is, the expectation is with respect to _k):

    ^1/2[x_k+1 -x_*  |  x_k]     =    ((1+β) - ω^-1/2[]  ^-1/2) ^1/2(x_k -x_*) - β^1/2(x_k-1-x_*) .


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

    ^1/2[x_k+1 -x_*]     =    ^1/2[[x_k+1 -x_*  |  x_k]]
       (<ref>)=   ((1+β) - ω^-1/2[] ^-1/2) ^1/2[x_k -x_*] - β^1/2[x_k-1-x_*].



Plugging the eigenvalue decomposition Λ^⊤ of the matrix =^-1/2[] ^-1/2 into the above, and multiplying both sides from the left by ^⊤, we obtain

    ^⊤^1/2[x_k+1 -x_*]  = ^⊤((1+β) - ωΛ^⊤)^1/2[x_k -x_*] - β^⊤^1/2[x_k-1-x_*].


Let us define s_k^⊤^1/2[x_k -x_*]  ∈^n. Then relation (<ref>) takes the form of the recursion

    s_k+1= [(1+β) - ωΛ ] s_k - β s_k-1,

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

    s_k+1^i= [(1+β) - ωλ_i ] s_k^i - β s_k-1^i     for all   i= 1,2,3,...,n,

where s_k^i indicates the ith coordinate of s_k.

We will now fix i and analyze recursion (<ref>) using Lemma <ref>. Note that (<ref>) is a second degree linear homogeneous recurrence relation of the form (<ref>) with a_1=1+β - ωλ_i and a_2=- β. Recall  that 0≤λ_i ≤1 for all i. Since we assume that 0< ω≤ 1/λ_max, we know that 0≤ωλ_i ≤ 1 for all i. We now consider two  cases:



  * λ_i =0.   

In this case, (<ref>) takes the form: 
    s_k+1^i=(1+β)s_k^i-β s_k-1^i .
 Applying Theorem <ref>, we know that x_*=Π_ℒ^(x_0)=Π_ℒ^(x_1). Using Lemma <ref> twice, once for x=x_0 and then for x=x_1, we observe that s_0^i=u_i^⊤^1/2 (x_0-x_*)=0 and s_1^i=u_i^⊤^1/2 (x_1-x_*)=0. Finally, in view of (<ref>) we conclude that 
    s_k^i=0   for all   k≥ 0 .





  * λ_i >0. 

Since 0<ωλ_i ≤ 1 and β≥ 0, we have 1+β - ωλ_i ≥ 0 and hence

    a_1^2+4a_2=(1+β -ωλ_i)^2-4β≤  (1+β -ωλ_min^+)^2-4β  < 0,

where the last inequality can be shown to hold[The lower bound on β is tight. However,  the upper bound is not. However, we do not care much about the regime of large β as β is the convergence rate, and hence is only interesting if smaller than 1.] for  (1-√(ωλ_min^+))^2 < β < 1. Applying Lemma <ref> the following bound can be deduced

    s_k^i     =    2(-a_2)^k/2 (C_0  cos(θ k) +C_1 sin(θ k))  ≤  2 β^k/2 P_i,

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≥0 we have

    [x_k -x_*]_^2   =   [x_k -x_*]^⊤[x_k -x_*]   =  [x_k -x_*] ^1/2^⊤^1/2[x_k -x_*]  
       =   ^⊤^1/2[x_k -x_*] _2^2   =  s_k^2  =  ∑_i=1^n (s_k^i)^2 
       =    ∑_i: λ_i=0 (s_k^i)^2  +  ∑_i:  λ_i >0 (s_k^i)^2  (<ref>)= ∑_i:  λ_i >0 (s_k^i)^2
       (<ref>)≤   ∑_i:  λ_i >0 4 β ^k P_i^2 
       =   β^k C,

where C=4∑_i:  λ_i >0  P_i^2.







§ PROOF OF THEOREM <REF> 
 

The proof follows a similar pattern to that of Theorem <ref>. 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^⊤(x_k-x_k-1)e_i. First, we decompose

    x_k+1-x_*^2 
        =   x_k-ω∇ f__k(x_k)+β r_k^i -x_*^2 
        =   x_k-ω∇ f__k(x_k)-x_*^2+2⟨ x_k-ω∇ f__k(x_k)-x_*, β r_k^i  ⟩ + β^2 r_k^i^2.

 We shall use the tower property in the form

    [[[ X   |  x_k,  _k]  |  x_k]] = [X],

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 (<ref>), we get

    [x_k+1-x_*^2  |  x_k, _k ] 
       =   [x_k-ω∇ f__k(x_k)-x_*^2 |  x_k, _k]_1
             +[2⟨ x_k-ω∇ f__k(x_k)-x_*, β r_k^i  ⟩ |  x_k, _k]_2
             +[β^2 r_k^i^2 |  x_k, _k]_3.


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

    1   =   x_k-ω∇ f__k(x_k)-x_*^2
       =   x_k-x_*^2-2 ω⟨ x_k-x_*,∇ f__k(x_k) ⟩ +ω^2 ∇ f__k(x_k)^2
       (<ref>) + (<ref>)=   x_k-x_*^2-4ω f__k(x_k)+2ω^2 f__k(x_k)
       =   x_k-x_*^2-2ω(2-ω)f__k(x_k).


We will now bound the second expression. Using the identity

    [r_k^i   |  x_k, _k] = _i [r_k^i] = ∑_i=1^n 1/nr_k^i = 1/n(x_k-x_k-1),

we can write

    2   =   [2⟨ x_k-ω∇ f__k(x_k)-x_*, β r_k^i  ⟩ |  x_k, _k]
       =    2⟨ x_k-ω∇ f__k(x_k)-x_*, β[r_k^i   |  x_k, _k] ⟩
       (<ref>)=    2⟨ x_k-ω∇ f__k(x_k)-x_*, βn (x_k - x_k-1)  ⟩
        =   
    2βn⟨ x_k-x_*, x_k-x_k-1⟩ +2ωβn⟨∇ f__k(x_k), x_k-1 - x_k  ⟩
       =    2βn⟨ x_k-x_*, x_k-x_*⟩ + 2βn⟨ x_k-x_*, x_*-x_k-1⟩ +2ωβn⟨∇ f__k(x_k),x_k-1- x_k  ⟩
       =    2βnx_k-x_*^2 + 2βn⟨ x_k-x_*, x_*-x_k-1⟩ +2ωβn⟨∇ f__k(x_k),x_k-1- x_k  ⟩.

Using the fact that for arbitrary vectors a,b,c ∈^n we have the identity
2 ⟨ a-c,c-b ⟩ =a-b^2-c-b^2-a-c^2,
we obtain

    2 ⟨ x_k-x_*, x_*-x_k-1⟩=  x_k-x_k-1^2- x_k-1-x_*^2-x_k-x_*^2.

Substituting this into (<ref>) gives

    2    = βnx_k-x_*^2+βnx_k-x_k-1^2-βnx_k-1-x_*^2 + 2ωβn⟨∇ f__k(x_k),x_k-1- x_k  ⟩.


The third expression can be bound as

    3   =   [β^2 r_k^i^2 |  x_k, _k] 
       =   β^2 _i[ r_k^i^2] 
       =   β^2 ∑_i=1^n 1n (x_k^i -x_k-1^i)^2 
       =   β^2nx_k-x_k-1^2 
       =   β^2n(x_k-x_*)+(x_*-x_k-1)^2 
       ≤   2β^2nx_k-x_*^2+ 2β^2nx_k-1-x_*^2.


By substituting the  bounds (<ref>), (<ref>), (<ref>) into  (<ref>) we obtain

    [x_k+1-x_*^2  |  x_k, _k ] 
       ≤   x_k-x_*^2-2ω(2-ω) f__k(x_k)
             + βnx_k-x_*^2+ βnx_k-x_k-1^2 -βnx_k-1-x_*^2 
             + 2ωβn⟨∇ f__k(x_k), x_k-1- x_k  ⟩  +  2β^2nx_k-x_*^2 
             + 2β^2nx_k-1-x_*^2 
       (<ref>)≤   (1+3βn + 2β^2n)x_k-x_*^2+ (βn + 2β^2n) x_k-1-x_*^2 
             - 2ω(2-ω)f__k(x_k) + 2ωβn⟨∇ f__k(x_k),x_k-1- x_k  ⟩.

We now take the middle expectation (see (<ref>)) and apply it to inequality (<ref>):

    [[x_k+1-x_*^2  |  x_k, _k ]  |  x_k]     ≤   (1+3βn + 2β^2n)x_k-x_*^2+ (βn + 2β^2n) x_k-1-x_*^2 
             -2ω(2-ω)f(x_k) + 2ωβn⟨∇ f(x_k),x_k-1- x_k  ⟩
       ≤   (1+3βn + 2β^2n)x_k-x_*^2+ (βn + 2β^2n) x_k-1-x_*^2 
             -2ω(2-ω)f(x_k) + 2ωβn(f(x_k-1)-f(x_k))
        =    (1+3βn + 2β^2n)x_k-x_*^2+ (βn + 2β^2n) x_k-1-x_*^2 
             - (2ω(2-ω) +2ωβn)f(x_k) + 2ωβn f(x_k-1).

where in the second step we used the inequality
 ⟨∇ f(x_k),x_k-1- x_k ⟩≤ f(x_k-1)-f(x_k) and the fact that ωβ≥ 0, which follows from the assumptions. We now apply  inequalities  (<ref>) and (<ref>), obtaining

    [[x_k+1-x_*^2  |  x_k, _k ]  |  x_k]   ≤   (1+3βn+2β^2n - (ω(2-ω) +ωβn)λ_min^+ )_a_1x_k-x_*^2 
             + 1n(β +2β^2 + ωβλ_max)_a_2x_k-1-x_*^2.






By taking expectation again (outermost expectation in the tower rule (<ref>)), and letting F_k[x_k-x_*^2_], we get  the relation

    F_k+1≤ a_1 F_k  + a_2 F_k-1 .


It suffices to apply  Lemma <ref> to the relation  (<ref>). The conditions of the lemma are satisfied. Indeed, a_2≥ 0, and if a_2=0, then β=0 and hence a_1=1-ω(1-ω)λ_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 (<ref>) to  (<ref>).



plain


