# Stochastic Learning under Random Reshuffling

In empirical risk optimization, it has been observed that stochastic gradient implementations that rely on random reshuffling of the data achieve better performance than implementations that rely on sampling the data uniformly. Recent works have pursued justifications for this behavior by examining the convergence rate of the learning process under diminishing step-sizes. This work focuses on the constant step-size case. In this case, convergence is guaranteed to a small neighborhood of the optimizer albeit at a linear rate. The analysis establishes analytically that random reshuffling outperforms uniform sampling by showing explicitly that iterates approach a smaller neighborhood of size O(μ^2) around the minimizer rather than O(μ). Furthermore, we derive an analytical expression for the steady-state mean-square-error performance of the algorithm, which helps clarify in greater detail the differences between sampling with and without replacement. We also explain the periodic behavior that is observed in random reshuffling implementations.

## Authors

• 11 publications
• 25 publications
• 24 publications
• 62 publications
• ### A Robust Variable Step Size Fractional Least Mean Square (RVSS-FLMS) Algorithm

In this paper, we propose an adaptive framework for the variable step si...
11/14/2017 ∙ by Shujaat Khan, et al. ∙ 0

• ### On the Influence of Momentum Acceleration on Online Learning

The article examines in some detail the convergence rate and mean-square...
03/14/2016 ∙ by Kun Yuan, et al. ∙ 0

• ### Finite-Time Error Bounds For Linear Stochastic Approximation and TD Learning

We consider the dynamics of a linear stochastic approximation algorithm ...
02/03/2019 ∙ by R. Srikant, et al. ∙ 0

• ### Constant Step Size Least-Mean-Square: Bias-Variance Trade-offs and Optimal Sampling Distributions

We consider the least-squares regression problem and provide a detailed ...
11/29/2014 ∙ by Alexandre Défossez, et al. ∙ 0

• ### Incremental Without Replacement Sampling in Nonconvex Optimization

Minibatch decomposition methods for empirical risk minimization are comm...
07/15/2020 ∙ by Edouard Pauwels, et al. ∙ 0

• ### On Uncensored Mean First-Passage-Time Performance Experiments with Multiwalk in R^p: a New Stochastic Optimization Algorithm

A rigorous empirical comparison of two stochastic solvers is important w...
12/06/2018 ∙ by Franc Brglez, et al. ∙ 0

• ### q-LMF: Quantum Calculus-based Least Mean Fourth Algorithm

Channel estimation is an essential part of modern communication systems ...
12/04/2018 ∙ by Alishba Sadiq, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Motivation

We consider minimizing an empirical risk function , which is defined as the sample average of loss values over a possibly large but finite training set:

 w⋆Δ=argminw∈RMJ(w)Δ=1NN∑n=1Q(w;xn) (1)

where the denotes the training data samples and the loss function is assumed convex and differentiable. We assume the empirical risk is strongly-convex which ensures that the minimizer, , is unique. Problems of the form (1

) are common in many areas of machine learning including linear regression, logistic regression and their regularized versions.

When the size of the dataset is large, it is impractical to solve (1) directly via traditional gradient descent by evaluating the full gradient at every iteration. One simple, yet powerful remedy is to employ the stochastic gradient method (SGD) [2, 3, 4, 5, 6, 7, 8, 9]. Rather than compute the full gradient over the entire data set, these algorithms pick one index at random at every iteration, and employ to approximate . Specifically, at iteration

, the update for estimating the minimizer is of the form

[10]:

 wi=wi−1−μ∇wQ(wi−1;xni), (2)

where

is the step-size parameter. Note that we are using boldface notation to refer to random variables. Traditionally, the index

is uniformly distributed over the discrete set

.

It has been noted in the literature [11, 12, 13, 14]

that incorporating random reshuffling into the gradient descent implementation helps achieve better performance. More broadly than in the case of the pure SGD algorithm, it has also been observed that applying random reshuffling in variance-reduction algorithms, like SVRG

[15], SAGA[16], can accelerate the convergence speed[17, 18, 19, 20]. The reshuffling technique has also been applied in distributed system to reduce the communication and computation cost[21].

In random reshuffling implementations, the data points are no longer picked independently and uniformly at random. Instead, the gradient descent algorithm is run multiple times over the data where each run is indexed by

and is referred to as an epoch. For each epoch, the original data is first reshuffled and then passed over in order. In this manner, the

-th sample of epoch can be viewed as , where the symbol represents a uniform random permutation of the indices. We can then express the random reshuffling algorithm for the th epoch in the following manner:

 wki=wki−1−μ∇wQ(wki−1;xσk(i)),i=1,…,N (3)

with the boundary condition:

 wk0=wk−1N (4)

In other words, the initial condition for epoch is the last iterate from epoch . The boldface notation for the symbols and in (3) emphasizes the random nature of these variables due to the randomness in the permutation operation. While the samples over one epoch are no longer picked independently from each other, the uniformity of the permutation function implies the following useful properties[19, 22, 23]:

 σk(i)≠ σk(j),1≤i≠j≤N (5) P[ σk(i)=n ]= 1N,1≤n≤N (6) P[σk(i+1)=n|σk(1:i)]= (7)

where represents the collection of permuted indices for the samples numbered through .

Several recent works [12, 13, 24] have pursued justifications for the enhanced behavior of random reshuffling implementations over independent sampling (with replacement). The work [13] examined the convergence rate of the learning process under diminishing step-sizes, i.e., , where is some positive constant. It analytically showed that, for strongly convex objective functions, the convergence rate under random reshuffling can be improved from in vanilla SGD[25] to . The incremental gradient methods[26, 27], which can be viewed as the deterministic version of random reshuffling, shares similar conclusions, i.e., random reshuffling helps accelerate the convergence rate from to under decaying step-sizes. Also, in the work [24], it establishes that random reshuffling will not degrade performance relative to the stochastic gradient descent implementation, provided the number of epochs is not too large.

In this work, we focus on a different setting than[12, 13, 24] involving random reshuffling under constant rather than decaying step-sizes. In this case, convergence is only guaranteed to a small neighborhood of the optimizer albeit at a linear rate. The analysis will establish analytically that random reshuffling outperforms independent sampling (with replacement) by showing that the mean-square-error of the iterate at the end of each run in the random reshuffling strategy will be in the order of . This is a significant improvement over the performance of traditional stochastic gradient descent, which is  [10]. Furthermore, we derive an analytical expression for the steady-state mean-square-error performance of the algorithm, which is exact for quadratic risks and provides a good approximation for general risks. This helps clarify in greater detail the differences between sampling with and without replacement We also explain the periodic behavior that is observed in random reshuffling implementations.

### I-a Overview of results

• Section II provides a stability proof, which shows that under constant step-size random reshuffling will converge into a small neighborhood around the minimizer. The radius of the neighborhood improves from under uniform sampling to under random reshuffling. — Theorem 1.

• Next, we examine more closely the value of the scaling constant in the factor by introducing a long term model and deriving an expression for its mean-square-deviation (MSD) performance — Theorem 2. The theorem reveals how the number of samples , step-size , and Hessian of the loss function impact performance.

• In Theorem 3 we provide an expression for an upper bound for the MSD performance at all points close to steady-state. The result of the theorem helps explain the periodic behavior that is observed in random reshuffling implementations.

• The mismatch between the original reshuffling and the long model is provided in Lemma 2.

• Inspired by quadratic risks, we simplify the MSD expressions in Theorems 2 and 3 by using the hyperbolic tanh functions — equations (76) and (80).

• In equations (81) – (85), we show that as the sample size increases, the established MSD expression in Theorem 2 will regress to the same expression as the uniform sampling case.

## Ii Analysis of the Stochastic Gradient Under Random Reshuffling

### Ii-a Properties of the Gradient Approximation

We start by examining the properties of the stochastic gradient under random reshuffling. One main source of difficulty that we shall encounter in the analysis of performance under random reshuffling is the fact that a single sample of the stochastic gradient is now a biased estimate of the true gradient and, moreover, it is no longer independent of past selections, . This is in contrast to implementations where samples are picked independently at every iteration. Indeed, note that conditioned on previously picked data and on the previous iterate, we have:

 E[∇wQσk(i)(wki−1)|wki−1,σk(1:i−1)] = 1N−i+1∑n∉σk(1:i−1)∇wQ(wki−1) ≠ ∇J(wk0) (8)

The difference (8) is generally nonzero in view of the definition (1). For the first iteration of every epoch however, it can be verified that the following holds:

 E[∇wQσk(i)(wk0)∣∣wk0](???)= 1NN∑n=1Q(wk0;xn) (???)= ∇J(wk0) (9)

since at the beginning of one epoch, no data has been selected yet. Perhaps surprisingly, we will be showing that the biased construction of the stochastic gradient estimate not only does not hurt the performance of the algorithm, but instead significantly improves it. In large part, the analysis will revolve around considering the accuracy of the gradient approximation over an entire epoch, rather than focusing on single samples at a time. Recall that by construction in random reshuffling, every sample is picked once and only once over one epoch. This means that the sample average (rather than the true mean) of the gradient noise process is zero since

 1NN∑i=1∇wQ(w;xσk(i))=∇J(w) (10)

for any and any reshuffling order . This property will become key in the analysis.

### Ii-B Convergence Analysis

We can now establish a key convergence and performance property for the random reshuffling algorithm, which provides solid analytical justification for its observed improved performance in practice.

To begin with, we assume that the risk function satisfies the following conditions, which are automatically satisfied by many learning problems of interest, such as mean-square-error or logistic regression analysis and their regularized versions — see, e.g.,

[28, 29, 30, 31, 32].

###### Assumption 1 (Condition on loss function).

It is assumed that is differentiable and has a -Lipschitz continuous gradient, i.e., for every and any :

 ∥∇wQ(w1;xn)−∇wQ(w2;xn)∥≤δn∥w1−w2∥ (11)

where . We also assume is -strongly convex:

 (∇wJ(w1)−∇wJ(w2))T(w1−w2) ≥ν∥w1−w2∥2 (12)

If we introduce , then each and are also -Lipschitz continuous.

The following theorem focuses on the convergence of the starting point of each epoch and establishes in (14) that it actually approaches a smaller neighborhood of size around . Afterwards, using this result, we also show that the same performance level holds for all iterates and not just for the starting points of the epochs.

To simplify the notation, we introduce the constant , which is the gradient noise variance at optimal point :

 KΔ=1NN∑n=1∥∇wQ(w⋆;xn)∥2 (13)
###### Theorem 1 (Stability of starting points).

Under assumption 1, the starting point of each run in (3), i.e., , satisfies

 limsupk→∞E∥wk0−w⋆∥2≤ 4μ2δ2N2ν2K=O(μ2) (14)

when the step-size is sufficiently small, namely, for 111 The proof in this theorem is based on the worst case scenario, which implies the inequalities hold for any realizations. Therefore, this proof is also applicable to the deterministic cyclic sampling case.. The convergence to steady-state regime occurs at an exponential rate, dictated by the parameter:

 αΔ=1−μνN/2 (15)
###### Proof.

See Appendix A

Having established the stability of the first point of every epoch, we can establish the stability of every point.

###### Corollary 1 (Full Stability).

Under assumption 1, it holds that

 limsupk→∞E∥wki−w⋆∥2= O(μ2) (16)

for all when the step-size is sufficiently small.

###### Proof.

See Appendix C

With the previous established Theorem 1, it is also easy to gain the convergence theorem under decaying step-sizes.

###### Corollary 2 (Convergence under decaying step-sizes).

Under assumption 1 and the decaying step-sizes is employed, the iterate converge to the minimizer exactly as with rate.

See Appendix D

## Iii Illustrating Behavior and Periodicity

In this section we illustrate the theoretical findings so far by numerical simulations. We consider the following logistic regression problem:

 minwJ(w)=1NN∑n=1Q(w;hn,γ(n)), (17)

where

is the feature vector,

is the scalar label, and

 Q(w;hn,γn)Δ=ρ∥w∥2+ln(1+exp(−γ(n)hTnw)). (18)

The constant is the regularization parameter. In the first simulation, we compare the performance of the standard stochastic gradient descent (SGD) algorithm (2) with replacement and the random reshuffling (RR) algorithm (3). We set and . Each

is generated from the normal distribution

, where is a diagonal matrix with each diagonal entry generated from the uniform distribution . To generate , we first generate an auxiliary random vector with each entry following . Next, we generate from a uniform distribution . If then is set as ; otherwise is set as . We select during all simulations. Figure 1 illustrates the mean-square-deviation (MSD) performance, i.e., , of the SGD and RR algorithms when . It is observed that the RR algorithm oscillates during the steady-state regime, and that the MSD at the is the best among all iterates during epoch . Furthermore, it is also observed that RR has better MSD performance than SGD. Similar observations also occur in Fig. 2, where . It is worth noting that the gap between SGD and RR is much larger in Fig. 2 than in Fig. 1.

Next, in the second simulation we verify the conclusion that the MSD for the starting point of each epoch for the random reshuffling algorithm, i.e., , can achieve instead of . We still consider the regularized logistic regression problem (17) and (18), and the same experimental setting. Recall that in Lemma 1, we proved that

 limsupk→∞E∥˜wk0∥2≤ O(μ2), (19)

which indicates that when is reduced a factor of 10, the MSD-performance should be improved by at least dB. We observe a decay of about 20dB per decade in Fig. 3 for a logistic regression problem with data points and 30dB per decade in Fig. 4 with .

## Iv Introducing a Long-Term Model

We proved in the earlier sections that the mean-square error under random reshuffling approaches a small neighborhood around the minimizer. Our objective now is to assess more accurately the size of the constant that multiplies in the result, and examine how this constant may depend on various parameters including the amount of data, , and the form of the loss function . To do that, we proceed in two steps. First, we introduce an auxiliary long-term model in (28) below and subsequently determine how far the performance of this model is from the original system described by (27) further ahead.

### Iv-a Error Dynamics

In order to quantify the performance of the random reshuffling implementation more accurately than the figure obtained earlier, we will need to impose a condition on the smoothness of the Hessian matrix of the risk function.

###### Assumption 2 (Hessian is Lipschitz continuous).

The risk function has a Lipschitz continuous Hessian matrix, i.e., there exists a constant , such that

 ∥∇2wJ(w1)−∇2wJ(w2)∥≤κ∥w1−w2∥ (20)

Under this assumption, the gradient vector, , can be expressed in Taylor expansion in the form[28, p. 378]:

 ∇wJ(w)=∇2wJ(w⋆)(w−w⋆)+ξ(w),∀w (21)

where the residual term satisfies:

 ∥ξ(w)∥≤κ2∥w−w⋆∥2 (22)

As such, we can rewrite algorithm (3) in the form:

 ˜wki= ˜wki−1+μ∇wJ(wki−1) +μ(∇wQ(wki−1;xσk(i))−∇wJ(wki−1)) = ˜wki−1−μ∇2wJ(w⋆)˜wki−1+μξ(wki−1) +μ(∇wQ(wki−1;xσk(i))−∇wJ(wki−1)) (23)

To ease the notation, we introduce the Hessian matrix and the gradient noise process:

 HΔ= ∇2wJ(w⋆) sσk(i)(wki−1)Δ= ∇wQ(wki−1;xσk(i))−∇wJ(wki−1) (24)

so that (23) is simplified as:

 ˜wki=(I−μH)˜wki−1+μξ(wki−1)+μsσk(i)(wki−1) (25)

Now property (9) motivates us to expand (25) into the following error recursion by adding and subtracting the same gradient noise term evaluated at :

 ˜wki= (I−μH)˜wki−1+μsσk(i)(wk0) +μ(sσk(i)(wki−1)−sσk(i)(wk0))noise mismatch+μξ(wki−1) (26)

Iterating (26) and using (4) we can establish the following useful relation, which we call upon in the sequel:

 ˜wk+10= (I−μH)N˜wk0+μN∑i=1(I−μH)N−isσk(i)(wk0) +μN∑i=1(I−μH)N−i(sσk(i)(wki−1)−sσk(i)(wk0)) +μN∑i=1(I−μH)N−iξ(wki−1) (27)

Note that recursion (27) relates to , which are the starting points of two successive epochs. In this way, we have now transformed recursion (3), which runs from one sample to another within the same epoch, into a relation that runs from one starting point to another over two successive epochs.

To proceed, we will ignore the last two terms in (27) and consider the following approximate model, which we shall refer to as a long-term model.

 ˜w′k+10=(I−μH)N˜w′k0+μN∑i=1(I−μH)N−isσk(i)(wk0)Δ=s′(wk0) (28)

Obviously, the state evolution will be different than (27) and is therefore denoted by the prime notation, . Observe, however, that in model (28) the gradient noise process is still being evaluated at the original state vector, , and not at the new state vector, .

### Iv-B Performance of the Long-Term Model across Epochs

Note that the gradient noise in (28) has the form of a weighted sum over one epoch. This noise clearly satisfies the property:

 E[s′(wk0)|wk0]= 0 (29)

We also know that satisfies the Markov property, i.e., it is independent of all previous and , where , conditioned on . To motivate the next lemma consider the following auxiliary setting.

Assume we have a collection of vectors in whose sum is zero. We define a random walk over these vectors in the following manner. At each time instant, we select a random vector uniformly and with replacement from this set and move from the current location along the vector to the next location. If we keep repeating this construction, we obtain behavior that is represented by the right plot in Fig. 5. Assume instead that we repeat the same experiment except that now we assume the data is first reshuffled and then vectors are selected uniformly without replacement. Because of the zero sum property, and because sampling is now performed without replacement, we find that in this second implementation we always return to the origin after selections. This situation is illustrated in the left plot of the same Fig. 5. The next lemma considers this scenario and provides useful expressions that allow us to estimate the expected location after or more (unitl ) movements. These results will be used in the sequel in our analysis of the performance of stochastic learning under RR.

###### Lemma 1.

Suppose we have a set of vectors with the constraint . Assume the elements of are randomly reshuffled and then selected uniformly without replacement. Let be any nonnegative constant, be any symmetric positive semi-definite matrix, and introduce

 RxΔ= 1NN∑i=1xixTi (30) Var(X)Δ= 1NN∑i=1∥xi∥2=\rm{\small Tr}(Rx) (31)

Define the following functions for any :

 f(n;X,β)Δ= E∥∥ ∥∥n∑j=1βn−jxσ(j)∥∥ ∥∥2 (32) F(n;X,B)Δ= E[n∑j=1Bn−jxσ(j)][n∑j=1xTσ(j)Bn−j] (33)

It then holds that

 f(n;X,β)= (∑n−1i=0β2i)N−(∑n−1i=0βi)2N−1Var(X) (34) F(n;X,B)= [n−1∑i=0BiRxBi]N−[n−1∑i=0Bi]Rx[n−1∑i=0Bi]N−1 (35)
###### Proof.

The proof is provided in Appendix E. ∎

We now return to the stochastic gradient implementation under random reshuffling. Recall from (10) that the stochastic gradient satisfies the zero sample mean property so that

 N∑i=1sσk(i)(w)=0 (36)

at any given point . Applying Lemma 1, we readily conclude that

 E[s′(wk0)s′(wk0)T|wk0] = N(∑N−1i=0(I−μH)iRks(I−μH)i)N−1 −[∑N−1i=0(I−μH)i]Rks[∑N−1i=0(I−μH)i]N−1 (37)

where

 RksΔ=1NN∑n=1sn(wk0)sn(wk0)T (38)

Similarly, we conclude for the gradient noise at the optimal :

 R′⋆sΔ= E[s′(w⋆)s′(w⋆)T] = N(∑N−1i=0(I−μH)iR⋆s(I−μH)i)N−1 −[∑N−1i=0(I−μH)i]R⋆s[∑N−1i=0(I−μH)i]N−1 (39)

where

 R⋆s= 1NN∑i=0∇Q(w⋆;xi)∇Q(w⋆;xi)T (40)
###### Theorem 2 (Performance of Long-term Model across Epochs).

Under assumptions 1 and 2, when the step size is sufficiently small, namely, for , the mean-square-deviation (MSD) of the long term model (28) is given by

 MSDltRRΔ= limsupk→∞∥w′k0−w⋆∥2 = μ2\rm{\small Tr}((I−(I−μH)2N)−1R′⋆s)+O(μ4) (41)

The convergence to steady-state regime occurs at an exponential rate, dictated by the parameter:

 αΔ=(1−μλmin(H))2N≈1−2μλmin(H)N (42)
###### Proof.

See Appendix F. ∎

The simulations in Fig. 6 show that the MSD expression (41) fits well the performance of the original random reshuffling algorithm. We will establish this fact analytically in the sequel. For now, the simulation is simply confirming that the performance of the long-term model is a good indication of the performance of the original stochastic gradient implementation under RR.

### Iv-C Performance of the Long-Term Model over Iterations

In the previous section we examined the performance of the long-term model at the starting points of successive epochs. In this section, we examine the performance of the same model at any iterate as time approaches . This analysis will help explain the oscillations that are observed in the learning curves in the simulations. First, similar to (39), we need to determine the covariance matrix for any . From Lemma 1, we immediately get that

 R′⋆s,iΔ= Es′i(w⋆)s′i(w⋆)T = N(∑i−1j=0(I−μH)jR⋆s(I−μH)j)N−1 −[∑i−1j=0(I−μH)j]R⋆s[∑i−1j=0(I−μH)j]N−1 (43)
###### Theorem 3 (Performance Upper-bound for Long- Term Model).

Under assumptions 1 and 2, when the step size satisfies , the upper-bound of mean-square-deviation (MSD) of the long term model (28) at all iterations is given by

 limk→∞E∥˜w′ki∥2 ≤ (1−μν)2iμ2\rm{\small Tr}((I−(I−μH)2N)−1R′⋆s) +(1−(1−μν)2i)μ2\rm{\small Tr% }((I−(I−μν)2i)−1R′⋆s,i) (44) Δ= ηiMSDltRR+(1−ηi)MSDltRR,i (45)
###### Proof.

See Appendix G.∎

We need to point out unlike that (41), expression (144) is an upper-bound rather than an actual performance expression. Still, this bound can help provide useful insights on the periodic behavior that is observed in the simulations. The expression (44) on the right-hand side is a convex combination of two performance measures as defined in (45), where the second term is always larger than the first term but approaching it as increases towards . This behavior will become clearer later in the context of an example and the hyperbolic representation in section V-B.

Before we continue, we would like to comment on the convergence curve under random reshuffling. Unlike the convergence curve under uniform sampling, we observe periodic fluctuations under random reshuffling in Figures 2 and 6. The main reason for this behavior is the fact that the gradient noise is no longer i.i.d. in steady-state. Specifically, the noise variance is now a function of the iterate and it assumes its lowest value at the beginning and end of every epoch. In lemma 1, we show that the variance of the random walk process resulting from random reshuffling at each iteration in Eq. (34). We plot the function for and in Fig. 7. Since the mean-square performance of the algorithm is related to the variance of the gradient noise, it is expected that this bell-shape behavior will be reflected in to the MSD curve as well, thus, resulting in better performance at the beginning and end of every epoch.

### Iv-D Mismatch Bound

Now we provide an upper bound on the mismatch between the long-term model (28) and the original algorithm (3).

###### Lemma 2 (Mismatch Bound).

After long enough iterations, i.e., , the difference between the long term model trajectory (28) and the original trajectory (3) is

 limsupk→∞ E∥˜w′k0−˜wk0∥2≤4μ2δ2N2ν2(N−1)K+O(μ3) (46)

Proof: See Appendix I.

## V Quadratic Risks and Hyperbolic Representation

Lastly, we consider an example involving a quadratic (least-squares) risk to show that, in this case, the long-term model provides the exact MSD for the original algorithm. The analysis will also provide some insights into expression (41). It also motivates a hyperbolic representation for the MSD, which helps provides some more insights into the MSD behavior.

Thus, consider the following quadratic risk function:

 minwJ(w)=12NN∑n=1∥Aw−xn∥2 (47)

where has full column rank. We have:

 ∇wJ(w)= ATAw−AT(1NN∑n=1xn)Δ=¯x (48) ∇wQ(w;xn)= ATAw−ATxn (49) ∇2J(wki)= ATA (50) sn(w)= AT(xn−¯x) (51)

Since the gradient noise is independent of , we have

 sn(wki)−sn(wk0)≡0 (52)

Moreover, since the risk is quadratic, it also holds that

 ξ(w)≡0 (53)

Therefore, the long-term model is exactly the same as the original algorithm. For this example, we can calculate the following quantities:

 w⋆= (ATA)−1AT¯x (54) R⋆s= AT1NN∑n=1(xn−¯x)(xn−¯x)TA=ATRxxA (55) Var(x)= 1NN∑n=1∥xn−¯x∥2 (56) I−μH= I−μATA (57)

In special case when the columns of are orthogonal and normalized, i.e., , we can simplify the MSD expression (41) by noting that

 R′⋆s = 1N−1(NN−1∑i=0(1−μ)2i−(N−1∑i=0(1−μ)i)2)ATRxxA = 1N−1⎛⎝N(1−(1−μ)2N)2μ−μ2−(1−(1−μ)N)2μ2⎞⎠ATRxxA (58)

and, hence,

 MSDRR= μ2\rm{\small Tr}((1−(1−μ)2N)−1R′⋆s) = μ2N−1(N2μ−μ2−(1−(1−μ)N)2μ2(1−(1−μ)2N))Var(x) = μ2N−1(N2μ−μ2−1−(1−μ)Nμ2(1+(1−μ)N))Var(x) (59)

In order to provide further insights on this MSD expression, we simplify it under a small assumption. We could introduce the Taylor series:

 (1−μ)N= 1−Nμ+O(N2μ2) (60)

However, this approximation can be bad if is large, which is not uncommon in big data. Instead, we appeal to:

 (1−μ)N=eNln(1−μ)=e−μN+O(μ2N)≈e−μN (61)

Notice it is instead of , and therefore (61) is a tighter approximation than (60) when is large. Based on this, we further approximate:

 1−(1−μ)N1+(1−μ)N≈tanh(μN/2) (62)

and arrive at the simplified expression:

 MSDRR≈ μN−1⎛⎜⎝N2−tanh(μN2)μ⎞⎟⎠Var(x) = μ2NN−1(1−2μNtanh(μN2))Var(x) (63)

For comparison purposes, we know that a simplified expression for MSD under uniform sampling has the following expression[10]:

 MSDus= μ2Var(x) (64)

Hence, the random reshuffling case has an extra multiplicative factor:

 mRRΔ=NN−1(1−2μNtanh(μN2)) (65)

We plot versus in the left plot of Fig. 8 where we ignore . Now it is clear from the figure that the smaller the step size or the smaller sample size are, the larger the improvement in performance is. In contrast, when goes to infinity, the term will converge to , i.e., the same performance as uniform sampling situation, which is consistent with the infinite-horizon case. Lastly, noting that

 R′⋆s,i = 1N−1(Ni−1∑j=0(1−μ)2j−(i−1∑j=0(1−μ)j)2)ATRxxA = 1N−1⎛⎝N(1−(1−μ)2i)2μ−μ2−(1−(1−μ)i)2μ2⎞⎠ATRxxA (66)

and using the approximation (61):

 R′⋆s,i≈NN−1(1−e−2μi2μ−(1−e−μi)2μ2N)ATRxxA (67)

Substituting in (44), we get for :

 limk→∞E∥˜w′ki∥2 ≈ e−2μiμ2NN−1(1−2μNtanh(μN2))Var(x) +(1−e−2μi)μ2NN−1(1−2μNtanh(μi2))Δ=mRR(i)Var(x) (68)

Since is monotonically increasing,