Curvature-Exploiting Acceleration of Elastic Net Computations

This paper introduces an efficient second-order method for solving the elastic net problem. Its key innovation is a computationally efficient technique for injecting curvature information in the optimization process which admits a strong theoretical performance guarantee. In particular, we show improved run time over popular first-order methods and quantify the speed-up in terms of statistical measures of the data matrix. The improved time complexity is the result of an extensive exploitation of the problem structure and a careful combination of second-order information, variance reduction techniques, and momentum acceleration. Beside theoretical speed-up, experimental results demonstrate great practical performance benefits of curvature information, especially for ill-conditioned data sets.


Introduction
Lasso, ridge and elastic net regression are fundamental problems in statistics and machine learning, with countless applications in science and engineering [33]. Elastic net regression amounts to solving the following convex optimization problem minimize x∈Ê d

2n
Ax for given data matrices A ∈ Ê n×d and b ∈ Ê n and regularization parameters γ 1 and γ 2 .
Setting γ 1 = 0 results in ridge regression, γ 2 = 0 yields lasso and letting γ 1 = γ 2 = 0 reduces the problem to the classical least-squares. Lasso promotes sparsity of the optimal solution, which sometimes helps to improve interpretability of the results. Adding the additional l 2 -regularizer helps to improve the performance when features are highly correlated [28,33].
The convergence rates of iterative methods for solving (1) are typically governed by the condition number of the Hessian matrix of the ridge loss, C + γ 2 I, where C = 1 n A ⊤ A is the sample correlation matrix. Real-world data sets often have few dominant features, while the other features are highly correlated with the stronger ones [13,28]. This translates to a rapidly decaying spectrum of C. In this paper, we demonstrate how this property can be exploited to reduce the effect of ill-conditioning and to design faster algorithms for solving the elastic net regression problem (1).

Related work
Over the past few years, there has been a great attention to developing efficient optimization algorithms for minimizing composite objective functions where f (x) = 1 n n i=1 f i (x) is a finite sum of smooth and convex component functions f i (x), and h (x) is a possibly non-smooth convex regularizer. In machine learning applications, the function f typically models the empirical data loss and the regularizer h is used to promote desired properties of a solution. For example, the elastic net objective can fit to this form with f i (x) = 1 2 a ⊤ i x − b i 2 + γ 2 x 2 2 , and h(x) = γ 1 x 1 .

First-Order methods
Standard deterministic first-order methods for solving (2), such as proximal gradient descent, enjoy linear convergence for strongly convex objective functions and are able to find an ǫapproximate solution in time O dnκ log 1 ǫ , where κ is the condition number of f . This runtime can be improved to O dn √ κ log 1 ǫ if it is combined with Nesterov acceleration [4,23]. However, the main drawback of these methods is that they need to access the whole data set in every iteration, which is too costly in many machine learning tasks.
For large-scale problems, methods based on stochastic gradients have become the standard choice for solving (2). Many linearly convergent proximal methods such as, SAGA [9] and Prox-SVRG [31], have been introduced and shown to outperform standard first-order methods under certain regularity assumptions. These methods improve the time complexity to O d (n +κ) log 1 ǫ , whereκ is a condition number satisfyingκ ≥ κ. When the component functions do not vary substantially in smoothness,κ ≈ κ, and this complexity is far better than those of deterministic methods above. By exploiting Nesterov momentum in different ways (see, e.g., [2,8,11,19]), one can improve the complexity to O d n + √ nκ log 1 ǫ , which is also optimal for this class of problems [30].

Second-order methods
Second-order methods are known to have superior performance compared to their first-order counterparts both in theory and practice, especially when the problem at hand is highly nonlinear and/or ill-conditioned. However, such methods often have very high computational cost per iteration. Recently, there has been an intense effort to develop algorithms which use second-order information with a more reasonable computational burden (see, e.g., [1,6,10,21,26,31,32] and references therein). Those methods use techniques such as random sketching, matrix sampling, and iterative estimation to construct an approximate Hessian matrix. Local and global convergence guarantees have been derived under various assumptions. Although many experimental results have shown excellent performance of those methods on many machine learning tasks, current second-order methods for finite-sum optimization tend to have much higher time-complexities than their first-order counterparts (see [32] for a detailed comparison).
Apart from having high time complexities, none of the methods cited above have any guarantees in the composite setting since their analyses hinge on differentiability of the objective function. Instead, one has to rely on methods that build on proximal Newton updates (see, e.g., [12,18,20,25]). However, these methods still inherit the high update and storage costs of conventional second-order methods or require elaborate tuning of several parameters and stopping criteria depending on a phase transition which occurs in the algorithm.

Ridge regression
For the smooth ridge regression problem, the authors in [13] have developed a preconditioning method based on linear sketching which, when coupled with SVRG, yields a significant speed-up over stochastic first-order methods. This is a rare second-order method that has a comparable or even better time complexity than stochastic first-order methods. More precisely, it has a guaranteed running time of O(d(n + κ H ) log 1 ǫ ), where κ H is a new condition number that can be dramatically smaller thanκ, especially when the spectrum of C decays rapidly. When d ≪ n, the authors in [29] combine sub-sampled Newton methods with the mini-batch SVRG to obtain some further improvements.

Contributions
Recently, the work [3] shows that under some mild algorithmic assumptions, and if the dimension is sufficiently large, the iteration complexity of second-order methods for smooth finite-sum problems composed of quadratics is no better than first-order methods. Therefore, it is natural to ask whether one can develop a second-order method for solving the elastic net problem which has improved practical performance but still enjoys a strong worst-case time complexity like the stochastic first-order methods do? It should be emphasized that due to the non-smooth objective, achieving this goal is much more challenging than for ridge regression. The preconditioning approach in [13] is not applicable, and the current theoretical results for second-order methods are not likely to offer the desired running time.
In this paper, we provide a positive answer to this question. Our main contribution is the design and analysis of a simple second-order method for the elastic net problem which has a strong theoretical time complexity and superior practical performance. The convergence bound adapts to the problem structure and is governed by the spectrum and a statistical

Algorithm
Time complexity 2nd-order yes measure of the data matrix. These quantities often yield significantly stronger time complexity guarantees for practical datasets than those of stochastic first-order methods (see Table 1). To achieve this, we first leverage recent advances in randomized low-rank approximation to generate a simple, one-shot approximation of the Hessian matrix. We then exploit the composite and finite-sum structure of the problem to develop a variant of the ProxSVRG method that builds upon Nesterov's momentum acceleration and inexact computations of scaled proximal operators, which may be of independent interest. We provide a simple convergence proof based on an explicit Lyapunov function, thus avoiding the use of sophisticated stochastic estimate sequences.

Preliminaries and Notation
Vectors are indicated by bold lower-case letters, and matrices are denoted by bold uppercase letters. We denote the dot product between x and y by x, y = x ⊤ y, and the Euclidean norm of x by x 2 = x, x . For a symmetric positive definite matrix H, x, y H = x ⊤ Hy is the H-inner product of two vectors x and y and x H = x, x H is the Mahalanobis norm of x. We denote by λ i (A) the ith largest eigenvalue of A. Finally, λ i denotes the ith largest eigenvalue of the correlation matrix C.
In the paper, we shall frequently use the notions of strong convexity and smoothness in the H-norm, introduced in the next two assumptions.
Assumption 1 implies that ∇f is L-Lipschitz: As a consequence, we have the following bound: Assumption 2. The function f (x) is µ-strongly convex w.r.t the H-norm, i.e, there exists a positive constant µ such that holds for all x, y ∈ Ê n and λ ∈ [0, 1].
Assumption 2 is equivalent to the requirement that holds for all x, y ∈ Ê n . We will use both of these definitions of strong convexity in our proofs.
At the core of our method is the concept of scaled proximal mappings, defined as follows: Definition 1 (Scaled Proximal Mapping). For a convex function h and a symmetric positive definite matrix H, the scaled proximal mapping of h at x is The scaled proximal mappings generalize the conventional ones: However, while many conventional prox-mappings admit analytical solutions, this is almost never the case for scaled proximal mappings. This makes it hard to extend efficient firstorder proximal methods to second-order ones. Fortunately, scaled proximal mappings do share some key properties with the conventional ones. We collect a few of them in the following result: ). The following properties hold: 1. prox H h (x) exists and is unique for x ∈ dom h. 2. Let ∂h (x) be the subdifferential of h at x, then Finally, in our algorithm, it will be enough to solve (3) approximately in the following sense: Definition 2 (Inexact subproblem solutions). We say that x + ∈ Ê d is an ǫ-optimal solution 3 Building Block 1: Randomized Low-Rank Approximation The computational cost of many Newton-type methods is dominated by the time required to compute the update direction d = H −1 g for some vector g ∈ Ê d and approximate Hessian H. A naive implementation using SVD would take O (nd 2 ) flops, which is prohibitive for large-scale data sets. A natural way to reduce this cost is to use truncated SVD. However, standard deterministic methods such as the power method and the Lanzcos method have run times that scale inversely with the gap between the eigenvalues of the input matrix. This gap can be arbitrarily small for practical data sets, thereby preventing us from obtaining the desired time complexity. In contrast, randomized sketching schemes usually admit gapfree run times [15]. However, unlike other methods, the block Lanczos method, detailed in Algorithm 1, admits both fast run times and strong guarantees on the errors between the true and the computed approximate singular vectors. This property turns out to be critical for deriving bounds on the condition number of the elastic net.
i be the SVD of A. Then, the following bounds hold with probability at least 9/10: Note that we only run Algorithm 1 once and ǫ ′ = 1/2 is sufficient in our work. Thus, the computational cost of this step is negligible, in theory and in practice.

Aproximating the Hessian
In this work, we consider the following approximate Hessian matrix of the ridge loss: Here, the first term is a natural rank r approximation of the true Hessian, while the second term is used to capture information in the subspace orthogonal to the column space of V r . The inverse of H in (6) admits the explicit expression so the evaluation of H −1 x has time complexity O (rd).

Bounding the Condition Number
We now turn our attention to studying how the approximate Hessian affects the relevant condition number of the elastic net problem. We first introduce a condition number that usually determines the iteration complexity of stochastic first-order methods under nonuniform sampling.
Definition 3. The average condition number of (1.1) is For the elastic net problem (1), the smooth part of the objective is the ridge loss .
Since we define smoothness and strong convexity of f i (x) in the H-norm, the relevant constants are For comparison, we also define the conventional condition numberκ, which characterizes the smoothness and strong convexity of f i (x) in the Euclidean norm. In this caseκ = and µ = λ d (C + γ 2 I) .
It will become apparent that κ H can be expressed in terms of a statistical measure of the ridge loss and that it may be significantly smaller thanκ. We start by introducing a statistical measure that has been widely used in the analysis of ridge regression (see, e.g., [16] and the references therein).

Definition 4 (Effective Dimension).
For a positive constant λ, the effective dimension of C is defined as The effective dimension generalizes the ordinary dimension and satisfies d λ ≤ d with equality if and only if λ = 0. It is typical that when C has a rapidly decaying spectrum, most of the λ i 's are dominated by λ, and hence d λ can be much smaller than d.
The following lemma bounds the eigenvalues of the matrix H −1/2 (C + γ 2 I) H −1/2 , which can be seen as the effective Hessian matrix.

Lemma 1 ([13]
). Invoking Algorithm 1 with data matrix 1 √ n A, target rank r, and target precision ǫ ′ = 1/2, it holds with probability at least 9/10 that Equipped with Lemma 1, we can now connect κ H with d λ using the following result.
Theorem 1. With probability at least 9/10, the following bound holds up to a multiplicative constant: Proof. See Appendix A.
compared toκ. If the spectrum of C decays rapidly, then the terms i>r λ i are negligible and the ratio is approximately i≤r λ i /(rλ r ). If the first eigenvalues are much larger than λ r , this ratio will be large. For example, for the australian data set [7], this ratio can be as large as 1.34 × 10 4 and 1.6 × 10 5 for r = 3 and r = 4, respectively. This indicates that it is possible to improve the iteration complexity of stochastic first-order methods if one can capitalize on the notions of strong convexity and smoothness w.r.t the H-norm in the optimization algorithm. Of course, this is only meaningful if there is an efficient way to inject curvature information into the optimization process without significantly increasing the computational cost. In the smooth case, i.e., γ 1 = 0, this task can be done by a preconditioning step [13]. However, this approach is not applicable for the elastic net, and we need to make use of another building block.

Building Block 2: Inexact Accelerated Scaled Proximal SVRG
In this section, we introduce an inexact scaled accelerated ProxSVRG algorithm for solving the generic finite-sum minimization problem in (2). We then characterize the convergence rate of the proposed algorithm.

Description of the Algorithm
To motivate our algorithm, we first recall the ProxSVRG method from [31]: For the sth outer iteration with the corresponding outer iteratex s , let x 0 =x s and for k = 0, 2, . . . , where i k is drawn randomly from {1, . . . , n} with probability p i k = L i k /(nL avg ). Since we are provided with an approximate Hessian matrix H, it is natural to use the following update: which can be seen as a proximal Newton step with the full gradient vector replaced by v k . Note that when h(·) is the ℓ 1 -penalty, ProxSVRG can evaluate (8) in time O(d), while evaluating (9) amounts to solving an optimization problem. It is thus is critical to keep the number of such evaluations small, which then translates into making a sufficient progress at each iteration. A natural way to achieve this goal is to reduce the variance of the noisy 10: end for 11:x s+1 ← x T 12: end for Output:x S gradient v k . This suggests to use large mini-batches, i.e., instead of using a single component function f i k , we use multiple ones to form: where B k ⊂ {1, . . . , n} is a set of indices with cardinality |B k | = b. It is easy to verify that v k is an unbiased estimate of ∇f (x k ). Notice that naively increasing the batch size makes the algorithm increasingly similar to its deterministic counterpart, hence inheriting a high-time complexity. This makes it hard to retain the runtime of ProxSVRG in the presence of 2nd-order information.
In the absence of second-order information and under the assumption that the proximal step is computed exactly, the work [24] introduced a method called AccProxSVRG that enjoys the same time complexity as ProxSVRG but allows for much larger mini-batch sizes. In fact, it can tolerate a mini-batch of size O √κ thanks to the use of Nesterov momentum. This indicates that an appropriate use of Nesterov momentum in our algorithm could allow for the larger mini-batches required to balance the computational cost of using scaled proximal mappings. The improved iteration complexity of the scaled proximal mappings will then give an overall acceleration in terms of wall-clock time. As discussed in [2], the momentum mechanism in AccProxSVRG fails to accelerate ProxSVRG unlessκ ≥ n 2 . In contrast, as we will see, our algorithm will be able to accelerate the convergence also in these scenarios. In summary, our algorithm is designed to run in an inner-outer fashion as ProxSVRG with large mini-batch sizes and Nesterov momentum to compensate for the increased computational cost of subproblems. The overall procedure is summarized in Algorithm 2.

Convergence Argument
In this subsection, we will show that as long as the errors in evaluating the scaled proximal mappings are controlled in an appropriate way, the iterates generated by the outer loop of Algorithm 2 converge linearly in expectation to the optimal solution. Recall that in Step 7 of Algorithm 2, we want to find an ǫ k -optimal solution in the sense of (5) to the following problem: The next lemma quantifies the progress made by one inner iteration of the algorithm. Our proof builds on a Lyapunov argument using a Lyapunov function on the form: Proof. See Appendix C.
Remark 1. Our proof is direct and based on natural Lyapunov functions, thereby avoiding the use of stochastic estimate subsequences as in [24] which is already very complicated even when the subproblems are solved exactly and H = I. We stress that the result in Lemma 2 also holds for smaller mini-batch sizes, namely b ∈ 1, . . . , O( L avg /µ) , provided that the step size η is reduced accordingly. In favor of a simple proof, we only report the large mini-batch result here.
Equipped with Lemma 2, we can now characterize the progress made by one outer iteration of Algorithm 2.
Proof. See Appendix D.
Remark 2. The theorem indicates that if the errors in solving the subproblems are controlled appropriately, the outer iterates generated by Algorithm 2 converge linearly in expectation to the optimal solution. Since V 0 depends on x ⋆ , it is difficult to provide a general closed-form expression for the target precisions ǫ k . However, we will show below that with a certain policy for selecting the initial point, it is sufficient to run the solver a constant number of iterations independently of x ⋆ . We stress that the results in this section are valid for minimizing general convex composite functions (2) and not limited to the elastic net problem.

Warm-Start
The overall time complexity of Algorithm 2 depends strongly on our ability to solve problem (10) in a reasonable computational time. If one naively starts the solver at a random point, it may take many iterations to meet the target precision. Thus, it is necessary to have a well-designed warm-start procedure for initializing the solver. Intuitively, the current iterate x k can be a reasonable starting point since the next iterate x k+1 should not be too far away from x k . However, in order to achieve a strong theoretical running time, we use a rather different scheme inspired by [19]. Let us first define the vector u k = y k − ηH −1 v k for k ∈ {0, 1, . . . , T − 1} and the function Then, the kth subproblem seeks for x k+1 such that where x ⋆ k+1 is the exact solution. We consider the initialization policy which can be seen as one step of the proximal gradient method applied to p (z, u k−1 ) starting at the current x k with step size γ.
The following proposition characterizes the difference in objective realized by z 0 and x ⋆ k+1 .
The proposition, together with (13), implies that it suffices to find x k+1 such that This is significant since one only needs to reduce the residual error by a constant factor independent of the target precision. Note also that the condition number κ sub is much smaller than κ ≈ λ 1 (H)/(λ d + γ 2 ), and computing the gradient of the smooth part of p (z, u) only takes time O (rd) instead of O (nd) as in the original problem. Those properties imply that the subproblems can be solved efficiently by iterative methods, where only a small (and known) constant number of iterations is needed. The next section develops the final details of convergence proof.

Global Time Complexity
We start with the time complexity of Algorithm 2. Let T (α) be the number of gradient evaluations that a subproblem solver takes to reduce the residual error by a factor α. Then, by Proposition 2, one can find an ǫ k -optimal solution to the kth subproblem by at most T (κ sub / (1 − ρ)) gradient evaluations, where one gradient evaluation is equivalent to O (d) flops. Consider the same setting of Theorem 2 and suppose that the subproblems are initialized by (14). Then, the time complexity of Algorithm 2 is given by: where the first summand is due to the full gradient evaluation at each outer loop; the second one comes from the fact that one needs O( √ κ H ) inner iterations, each of which uses a minibatch of size O( √ κ H ); and the third one is the result of O( √ κ H ) inner iterations, each of which solves a subproblem that needs T (κ sub / (1 − ρ)) gradient evaluations. We can now put things together and state our main result.
Proposition 3. Suppose that the approximate Hessian matrix H is given by (6) and that Algorithm 2 is invoked with f (x) = 1 2n Ax − b 2 2 + γ 2 2 x 2 2 and h (x) = x 1 . Assume further that the subproblems are solved by the accelerated proximal gradient descent method [4,23]. Our method can find an ǫ-optimal solution in time Proof. The task reduces to evaluating the term T (κ sub / (1 − ρ)) in (16). Recall that the iteration complexity of the accelerated proximal gradient descent method for minimizing the function F (x) = f (x) + h(x), where f is a smooth and strongly convex function and h is a possibly non-smooth convex regularizer, initiliazed at x 0 is given by where κ is the condition number. By invoking the about result with F (x) = p (x, u k ), k+1 , x 0 = z 0 , κ = κ sub , and ǫ is the right-hand side of (15), it follows that the number of iterations for each subproblem can be bounded by In addition, each iteration takes time O (rd) to compute the gradient implying the time complexity Finally, since κ H ≫ κ sub , the proof is complete.
We can easily recognize that this time complexity has the same form as the stochastic first-order methods discussed in Section 1.2.1 with the condition numberκ replaced by κ H . It has been shown in Theorem 1 that κ H can be much smaller thanκ, especially, when C has a rapidly decaying spectrum. We emphasize that the expression in (17) is available for free to us after having approximated the Hessian matrix. Hence no tuning is required to set ǫ k for solving the subproblems.

Experimental Results
In this section, we perform numerical experiments to verify the efficacy of the proposed method on real world data sets [7,14]. We compare our method with several well-known first-order methods: FISTA [4] with optimal step-size; Prox-SVRG [31] with epoch length 2n/b as suggested by the authors; Katyusha1 [2] with epoch length 2n/b, Katyusha momentum τ 2 = 0.5/b as suggested by the author; and our method with epoch length 2n/b. Since Katyusha1 can use a mini-batch of size √ n without slowing down the convergence, we set b = √ n for all methods except FISTA. Finally, to make a fair comparison, for each algorithm above, we tune only the step size, from the set η × {10 k , 2 × 10 k , 5 × 10 k |k ∈ {0, ±1, ±2}}, where η is the theoretical step size, and report the one having smallest objective value. Other hyper-parameters are set to their theory-predicted values. For Katyusha1, we also compute the largest eigenvalue of the Hessian matrix in order to set its parameter τ 1 . All methods are initialized at 0 and run for up to 100 epochs (passes through the full data set). For the subproblems in Algorithm 2, we just simply run FISTA with √ κ sub log κ sub iterations as discussed in the previous section, without any further tunning steps. The value of r is chosen as a small fraction of d so that the preprocessing time of Algorithm 1 is negligible. Note that the available spectrum of C after running Algorithm 1 also provides an insightful way to choose r. Figure 2 shows the suboptimality in objective versus the number of epochs for different algorithms solving the elastic net problem. We can see that our method systematically outperforms the others in all settings, and that there is a clear correspondence between the spectrum of C in Fig. 6    the first-order methods make almost no progress in the first 100 epochs, while our method can find a high-accuracy solution within tens of epochs, demonstrating a great benefit of second-order information. On the other hand, for a well-conditioned data set that does not exhibit high curvature such as real-sim, ProxSVRG is comparable to our method and even outperforms Katyusha1. This agrees with the theoretical time complexities summarized in Table 1. Finally, we demonstrate a hard instance for low-rank approximation methods via the cina0 data set, where C has a very large condition number and slowly decaying dominant eigenvalues. In this case, low-rank approximation methods with a small approximate rank r may not be able to capture sufficient curvature information. We can see that even in this case curvature information helps to avoid the stagnation experienced by FISTA and ProxSVRG.

Conclusions
We have proposed and analyzed a novel second-order method for solving the elastic net problem. By carefully exploiting the problem structure, we demonstrated that it is possible to deal with the non-smooth objective and to efficiently inject curvature information into the optimization process without (significantly) increasing the computational cost per iteration. The combination of second-order information, fast iterative solvers, and a well-designed warm-start procedure results in a significant improvement of the total runtime complexity over popular first-order methods. An interesting direction for future research would be to go beyond the quadratic loss. We believe that the techniques developed in this work can be extended to more general settings, especially when the smooth part of the objective function is self-concordant.
[17] Chonghai Hu, James T. Kwok, and Weike Pan. Accelerated gradient methods for stochastic optimization and online learning. In Advances in Neural Information Processing Systems, pages 781-789, Vancouver, Canada, Dec 2009.

A Proof of Theorem 1
First, by the definition of L i , it holds that where the last inequality follows since The term a ⊤ i H −1 a i can be bounded as It can be verified that where we have used Lemma 1 in the last step. It follows that which implies that On the other hand, it was shown in [13, Theorem 4] that Finally, recall that µ = λ d H −1/2 (C + γ 2 I) H −1/2 . By Lemma 1, we readily have µ ≥ c 1 γ 2 for some constant c 1 , which concludes the proof.

B Some Useful Auxiliary Results
To facilitate the analysis, we collect some useful inequalities regarding the Mahalanobis norm that are used in the subsequent proofs.

C Proof of Lemma 2
Before proving the lemma, we rewrite the sequences {x k , y k , z k } k≥0 generated by Algorithm 2 in the following recurrence form: We also recall the following well-known three-point identity: which holds for any symmetric matrix H ≻ 0 and any vectors a, b, c ∈ Ê d .
We are now ready to prove Lemma 2. From the definition of the Lyapunov function in (11), we have that where the equality (a) follows from (19), and (b) follows from(18c). Using the identity (19) again for the last term in (20), we obtain where the last equality follows from (18b). By noting that where we have used (18a) and the fact that With these observations, we are now ready to bound F (x k+1 ). In particular, by invoking Lemma 3 with y = y t , (25), it follows that Thus, combining (27)- (30) and using the fact that Invoking Lemma 4 with u k = (1 − τ ) x k + τ x ⋆ and then applying Lemma 5, we obtain Therefore, by taking the expectation on both sides of (31) and using (32), we obtain By choosing η = 1 Lavg , τ = µ 2Lavg , and b ≥ 60 Lavg µ , it is readily verified that the second and third terms on the right-hand side of (33) become nonpositive, which concludes the proof.

D Proof of Theorem 2
To begin with, let e t = E {F (x t ) − F (x ⋆ )} andẽ s = E {F (x s ) − F (x ⋆ )}, then by applying inequality (12) recursively, we obtain It can be verified that Thus, from (34) and the facts that τ = µ 2Lavg , ξ t H −1 ≤ 2ǫ t /L avg , it holds that where By the definition of V k , we have µ Multiplying both sides of the above inequality by (1 − τ ) −k yields H , then we can write the previous inequality as where we have separated the term α 0 u 0 from the sum and added the positive term α k u k in the last step. Note that α 0 = 2 can be further bounded by where c 1 = 1 + √ 2 1−τ . It is readily verified that {B k } k is an increasing sequence and that {S k 2 µ (c 1 V 0 + B k )} k is an increasing sequence satisfying S 0 > u 2 0 . Therefore, by invoking Lemma 7 with u t , α k , and S k , we have for any k ≥ 1 that where we have used √ a + b ≤ √ a + √ b for any a, b ≥ 0. Thus, for any t ∈ {0, 1, . . . , k}, it holds that Having upper bounds of E z t − x ⋆ H , we can now substitute them into (35) to get where the equality follows from the definition of α t , and in the last step, we have added a positive term α k in the first sum as well as added and subtracted α 0 in the second sum. Note that the function √ 1 + x − x is decreasing on the interval [0, ∞), it follows that We thus have where in the last step, we have applied the inequality (a + b) 2 ≤ (1 + β) a 2 + 1 + 1 β b 2 twice with β = 9 and β = 1, respectively. Using the definitions of B k and α t , the above inequality can be rewritten as With our choice of ǫ t , it follows that where the third equality follows since by definition, e k−1 ≤ E V k−1 ≤ δ k−1 , and sinceẽ s = e 0 ≤ V 0 , and the last inequality follows since ρ = 0.9τ . By applying the above inequality recursively, we obtain where c 3 = c 2 + 2/9. Therefore, when T ≥ 4 3ρ log (27c 3 ), it holds that It follows that where the first inequality follows from the triangle inequality and the last one follows from Property 1. Therefore, it holds that Denote the filtration F k by F k = {f 0 , f 1 , . . . , f k }, where f i denotes all the randomness incurring at time i for all i ∈ {1, . . . , k}. Note that the triple (x k , y k , z k ) depends on F k−1 , but not on f k , we thus have Taking the expectation on both sides of (43) completes the proof.
Lemma 5 (Bounding Variance). Assume that the indices in the mini-batch B k are sampled independently from {1, . . . , n} with probabilities P i = L i nLavg , then conditioned on x k and y k , it holds that Proof. We follow the original proofs in [24,31] with a few modifications to take the Mahalanobis norm into account. For any i ∈ {1, . . . , n}, consider the function φ i defined by Then, φ i (x ⋆ ) = min x φ i (x) since ∇φ i (x ⋆ ) = 0 and φ i is convex. It can be checked that ∇φ i is L i -Lipschitz, we thus have which is equivalent to Multiplying both sides of the previous inequality by Lavg L i and averaging from 1, . . . , n gives