Year: 2021

Abstract:
We present a novel adaptive optimization algorithm for large-scale machine learning problems. Equipped with a low-cost estimate of local curvature and Lipschitz smoothness, our method dynamically adapts the search direction and step-size. The search direction contains gradient information preconditioned by a well-scaled diagonal preconditioning matrix that captures the local curvature information. Our methodology does not require the tedious task of learning rate tuning, as the learning rate is updated automatically without adding an extra hyperparameter. We provide convergence guarantees on a comprehensive collection of optimization problems, including convex, strongly convex, and nonconvex problems, in both deterministic and stochastic regimes. We also conduct an extensive empirical evaluation on standard machine learning problems, justifying our algorithm’s versatility and demonstrating its strong performance compared to other start-of-theart first-order and second-order methods.

heading: 1 INTRODUCTION
text: This paper presents an algorithm for solving empirical risk minimization problems of the form:
minw∈RdF (w) := 1 n ∑n i=1f(w;x i, yi) = 1n ∑n i=1fi(w), (1)
where w is the model parameter/weight vector, {(xi, yi)}ni=1 are the training samples, and fi : Rd → R is the loss function. Usually, the number of training samples, n, and dimension, d, are large, and the loss function F is potentially nonconvex, making this minimization problem difficult to solve.
In the past decades, significant effort has been devoted to developing optimization algorithms for machine learning. Due to easy implementation and low per-iteration cost, (stochastic) first-order methods (Robbins & Monro, 1951; Duchi et al., 2011; Schmidt et al., 2017; Johnson & Zhang, 2013; Nguyen et al., 2017; 2019; Kingma & Ba, 2014; Jahani et al., 2021a; Recht et al., 2011) have become prevalent approaches for many machine learning applications. However, these methods have several drawbacks: (i) they are highly sensitive to the choices of hyperparameters, especially learning rate; (ii) they suffer from ill-conditioning that often arises in large-scale machine learning; and (iii) they offer limited opportunities in distributed computing environments since these methods usually spend more time on “communication” instead of the true “computation.” The main reasons for the aforementioned issues come from the fact that first-order methods only use the gradient information for their updates.
On the other hand, going beyond first-order methods, Newton-type and quasi-Newton methods (Nocedal & Wright, 2006; Dennis & Moré, 1977; Fletcher, 1987) are considered to be a strong family of optimizers due to their judicious use of the curvature information in order to scale the gradient. By exploiting the curvature information of the objective function, these methods mitigate many of the issues inherent in first-order methods. In the deterministic regime, it is known that these methods are relatively insensitive to the choices of the hyperparameters, and they handle ill-conditioned problems with a fast convergence rate. Clearly, this does not come for free, and these methods can have memory requirements up to O(d2) with computational complexity up to O(d3) (e.g., with a naive use of the Newton method). There are, of course, efficient ways to solve the Newton system with significantly lower costs (e.g., see Nocedal & Wright (2006)). Moreover, quasi-Newton methods require lower memory and computational complexities than Newton-type methods. Recently, there has been shifted attention towards stochastic second-order (Roosta-Khorasani & Mahoney, 2018;
Byrd et al., 2011; Martens, 2010; Jahani et al., 2020a; Xu et al., 2017; Roosta et al., 2018; Yao et al., 2018) and quasi-Newton methods (Curtis, 2016; Berahas et al., 2016; Mokhtari & Ribeiro, 2015; Jahani et al., 2021b; Berahas et al., 2019; Jahani et al., 2020b) in order to approximately capture the local curvature information. These methods have shown good results for several machine learning tasks (Xu et al., 2020; Berahas et al., 2020; Yao et al., 2019). In some cases, however, due to the noise in the Hessian approximation, their performance is still on par with the first-order variants. One avenue for reducing the computational and memory requirements for capturing curvature information is to consider just the diagonal of the Hessian. Since the Hessian diagonal can be represented as a vector, it is affordable to store its moving average, which is useful for reducing the impact of noise in the stochastic regime. To exemplify this, AdaHessian Algorithm (Yao et al., 2020) uses Hutchinson’s method (Bekas et al., 2007) to approximate the Hessian diagonal, and it uses a second moment of the Hessian diagonal approximation for preconditioning the gradient. AdaHessian achieves impressive results on a wide range of state-of-the-art tasks. However, its preconditioning matrix approximates the Hessian diagonal only very approximately, suggesting that improvements are possible if one can better approximate the Hessian diagonal.
1 0 0 0 0 S a m p le s
In this paper, we propose the dOubly Adaptive Scaled algorIthm for machine learning using Second-order information (OASIS). OASIS approximates the Hessian diagonal in an efficient way, providing an estimate whose scale much more closely approximates the scale of the true Hessian diagonal (see Figure 1). Due to this improved scaling, the search direction in OASIS contains gradient information, in which the components are well-scaled by the novel preconditioning matrix. Therefore, every gradient component in each dimension is adaptively scaled based on the approximated curvature for that dimension. For this reason, there is no need to tune the learning rate, as it would be updated automatically based on a local approximation of the Lipschitz smoothness parameter (see Figure 2). The well-scaled preconditioning matrix coupled with the adaptive learning rate results in a fully adaptive step for updating the parameters. Here, we provide a brief summary of our main contributions:
• Novel Optimization Algorithm. We propose OASIS as a fully adaptive method that preconditions the gradient information by a well-scaled Hessian diagonal approximation. The gradient component in each dimension is adaptively scaled by the corresponding curvature approximation. • Adaptive Learning Rate. Our methodology does not require us to tune the learning rate, as it is updated automatically via an adaptive rule. The rule approximates the Lipschitz smoothness parameter, and it updates the learning rate accordingly. • Comprehensive Theoretical Analysis. We derive convergence guarantees for OASIS with respect to different settings of learning rates, namely the case with adaptive learning rate for convex and strongly convex cases. We also provide the convergence guarantees with respect to fixed learning rate and line search for both strongly convex and nonconvex settings. • Competitive Numerical Results. We investigate the empirical performance of OASIS on a variety of standard machine learning tasks, including logistic regression, nonlinear least squares problems, and image classification. Our proposed method consistently shows competitive or superior performance in comparison to many first- and second-order state-of-the-art methods.
Notation. By considering the positive definite matrix D, we define the weighted Euclidean norm of vector x ∈ Rd with ‖x‖2D = xTDx. Its corresponding dual norm is shown as ‖ · ‖∗D. The operator is used as a component-wise product between two vectors. Given a vector v, we represent the corresponding diagonal matrix of v with diag(v).

heading: 2 RELATED WORK
text: In this paper, we analyze algorithms with the generic iterate updates:
wk+1 = wk − ηkD̂k−1mk, (2)
where D̂k is the preconditioning matrix, mk is either gk (the true gradient or the gradient approximation) or the first moment of the gradient with momentum parameter β1 or the bias corrected first moment of the gradient, and ηk is the learning rate. The simple interpretation is that, in order to update the iterates, the vector mk would be rotated and scaled by the inverse of preconditioning matrix D̂k, and the transformed information would be considered as the search direction. Due to limited space, here we consider only some of the related studies with a diagonal preconditioner. For more general preconditioning, see Nocedal & Wright (2006). Clearly, one of the benefits of a well-defined diagonal preconditioner is the easy calculation of its inverse. There are many optimization algorithms that follow the update in (2). A well-known method is stochastic gradient descent (SGD).The idea behind SGD is simple yet effective: the preconditioning matrix is set to be D̂k = Id, for all k ≥ 0. There are variants of SGD with and without momentum. The advantage of using momentum is to smooth the gradient (approximation) over the past iterations, and it can be useful in the noisy settings. In order to converge to the stationary point(s), the learning rate in SGD needs to decay. Therefore, there are many important hyperparameters that need to be tuned, e.g., learning rate, learning-rate decay, batch size, and momentum. Among all of them, tuning the learning rate is particularly important and cumbersome since the learning rate in SGD is considered to be the same for all dimensions. To address this issue, one idea is to use an adaptive diagonal preconditioning matrix, where its elements are based on the local information of the iterates. One of the initial methods with a non-identity preconditioning matrix is Adagrad (Duchi et al., 2011). In Adagrad, the momentum parameter is set to be zero (mk = gk), and the preconditioning matrix is defined as:
D̂k = diag( √∑k i=1gk gk). (3)
As is clear from the preconditioning matrix D̂k in (3), every gradient component is scaled with the accumulated information of all the past squared gradients. It is advantageous in the sense that every component is scaled adaptively. However, a significant drawback of D̂k in (3) has to do with the progressive increase of its elements, which leads to rapid decrease of the learning rate. To prevent Adagrad’s aggressive, monotonically decreasing learning rate, several approaches, including Adadelta (Zeiler, 2012) and RMSProp (Tieleman & Hinton, 2012), have been developed. Specifically, in RMSProp, the momentum parameter β1 is zero (or mk = gk) and the preconditioning matrix is as follows:
D̂k = √ β2D̂k−12 + (1− β2)diag ( gk gk ) , (4)
where β2 is the momentum parameter used in the preconditioning matrix. As we can see from the difference between the preconditioning matrices in (3) and (4), in RMSProp an exponentially decaying average of squared gradients is used, which prevents rapid increase of preconditioning components in (3). Another approach for computing the adaptive scaling for each parameter is Adam (Kingma & Ba, 2014). Besides storing an exponentially decaying average of past squared gradients like Adadelta and RMSprop, Adam also keeps first moment estimate of gradient, similar to SGD with momentum. In Adam, the bias-corrected first and second moment estimates, i.e., mk and D̂k in (2), are as follows:
mk = 1−β1 1−βk1
∑k i=1β k−i 1 gi, D̂k = √ 1−β2 1−βk2 ∑k i=1β k−i 2 diag(gi gi). (5)
There have been many other first-order methods with adaptive scaling (Loshchilov & Hutter, 2017; Chaudhari et al., 2019; Loshchilov & Hutter, 2016; Shazeer & Stern, 2018). The methods described so far have only used the information of the gradient for preconditioning mk in (2). The main difference of second-order methods is to employ higher order information for scaling and rotating the mk in (2). To be precise, besides the gradient information, the (approximated) curvature information of the objective function is also used. As a textbook example, in Newton’s method D̂k = ∇2F (wk) and mk = gk with ηk = 1. Diagonal Approximation. Recently, using methods from randomized numerical linear algebra, the AdaHessian method was developed (Yao et al., 2020). AdaHessian approximates the diagonal of the Hessian, and it uses the second moment of the diagonal Hessian approximation as the preconditioner D̂k in (2). In AdaHessian, Hutchinson’s method1 is used to approximate the Hessian diagonal as follows: Dk ≈ diag(E[zk ∇2F (wk)zk]), (6) where zk is a random vector with Rademacher distribution. Needless to say,2 the oracle ∇2F (wk)zk, or Hessian-vector product, can be efficiently calculated; in particular, for AdaHessian, it is computed with two back-propagation rounds without constructing the Hessian explicitly. In a nutshell, the first momentum for AdaHessian is the same as (5), and its second order momentum is:
D̂k = √
1−β2 1−βk2
∑k i=1β k−i 2 D 2 i . (7)
The intuition behind AdaHessian is to have a larger step size for the dimensions with shallow loss surfaces and smaller step size for the dimensions with sharp loss surfaces. The results provided by AdaHessian show its strength by using curvature information, in comparison to other adaptive first-order methods, for a range of state-of-the-art problems in computer vision, natural language processing, and recommendation systems (Yao et al., 2020). However, even for AdaHessian, the preconditioning matrix D̂k in (7) does not approximate the scale of the actual diagonal of the Hessian particularly well (see Figure 1). One might hope that a better-scaled preconditioner would enable better use of curvature information. This is one of the main focuses of this study. Adaptive Learning Rate. In all of the methods discussed previously, the learning rate ηk in (2) is still a hyperparameter which needs to be manually tuned, and it is a critical and sensitive hyperparameter. It is also necessary to tune the learning rate in methods that use approximation of curvature information (such as quasi-Newton methods like BFGS/LBFGS, and methods using diagonal Hessian approximation like AdaHessian). The studies (Loizou et al., 2020; Vaswani et al., 2019; Chandra et al., 2019; Baydin et al., 2017; Malitsky & Mishchenko, 2020) have tackled the issue regarding tuning learning rate, and have developed methodologies with adaptive learning rate, ηk, for first-order methods. Specifically, the work (Malitsky & Mishchenko, 2020) finds the learning rate by approximating the Lipschitz smoothness parameter in an affordable way without adding a tunable hyperparameter which is used for GD-type methods (with identity norm). Extending the latter approach to the weighted-Euclidean norm is not straightforward. In the next section, we describe how we can extend the work (Malitsky & Mishchenko, 2020) for the case with weighted Euclidean norm. This is another main focus of this study. In fact, while we focus on AdaHessian, any method with a positive-definite preconditioning matrix and bounded eigenvalues can benefit from our approach.

heading: 3 OASIS
text: In this section, we present our proposed methodology. First, we focus on the deterministic regime, and then we describe the stochastic variant of our method.
1For a general symmetric matrix A, E[z Az] equals the diagonal of A (Bekas et al., 2007). 2Actually, it needs to be said: many within the machine learning community still maintain the incorrect belief
that extracting second order information “requires inverting a matrix.” It does not.
3.1 DETERMINISTIC OASIS
Similar to the methods described in the previous section, our OASIS Algorithm generates iterates
Algorithm 1 OASIS Input:w0, η0,D0, θ0 = +∞
1: w1 = w0 − η0D̂0−1∇F (w0) 2: for k = 1, 2, . . . do 3: Form Dk via (8) and D̂k via (9) 4: Update ηk based on (10) 5: Set wk+1 = wk − ηkD̂k−1∇F (wk) 6: Set θk = ηkηk−1 7: end for according to (2). Motivated by AdaHessian, and by the fact that the loss surface curvature is different across different dimensions, we use the curvature information for preconditioning the gradient. We now describe how the preconditioning matrix D̂k can be adaptively updated at each iteration as well as how to update the learning rate ηk automatically for performing the step. To capture the curvature information, we also use Hutchinson’s method and update the diagonal approximation as follows: Dk = β2Dk−1 + (1− β2) diag(vk), where vk := zk ∇2F (wk)zk. (8) Before we proceed, we make a few more comments about the Hessian diagonal Dk in (8). As is clear from (8), a decaying exponential average of Hessian diagonal is used, which can be very useful in the noisy settings for smoothing out the Hessian noise over iterations. Moreover, it approximates the scale of the Hessian diagonal with a satisfactory precision, unlike AdaHessian Algorithm (see Figure 1). More importantly, the modification is simple yet very effective. Similar simple and efficient modification happens in the evolution of adaptive first-order methods (see Section 2). Further, motivated by (Paternain et al., 2019; Jahani et al., 2021b), in order to find a well-defined preconditioning matrix D̂k, we truncate the elements of Dk by a positive truncation value α. To be more precise:
(D̂k)i,i = max{|Dk|i,i, α}, ∀i ∈ [d]. (9) The goal of the truncation described above is to have a well-defined preconditioning matrix that results in a descent search direction; note the parameter α is equivalent to in Adam and AdaHessian). Next, we discuss the adaptive strategy for updating the learning rate ηk in (2). By extending the adaptive rule in (Malitsky & Mishchenko, 2020) and by defining θk := ηkηk−1 , ∀k ≥ 1, our learning rate needs to satisfy the inequalities: i) η2k ≤ (1 + θk−1)η2k−1, and ii) ηk ≤ ‖wk−wk−1‖D̂k
2‖∇F (wk)−∇F (wk−1)‖∗ D̂k
.
(These inequalities come from the theoretical results.) Thus, the learning rate can be adaptively calculated as follows: ηk = min{ √
1 + θk−1ηk−1, ‖wk − wk−1‖D̂k/(2‖∇F (wk)−∇F (wk−1)‖ ∗ D̂k )}. (10) As is clear from (10), it is only required to store the previous iterate with its corresponding gradient and the previous learning rate (a scalar) in order to update the learning rate. Moreover, the learning rate in (10) is controlled by the gradient and curvature information. As we will see later in the theoretical results, ηk ≥ α2L where L is the Lipschitz smoothness parameter of the loss function. It is noteworthy to highlight that due to usage of weighted-Euclidean norms the required analysis for the cases with adaptive learning rate is non-trivial (see Section 4). Also, by setting β2 = 1, α = 1, and D0 = Id, our OASIS algorithm covers the algorithm in (Malitsky & Mishchenko, 2020).
3.2 STOCHASTIC OASIS
In every iteration of OASIS, as presented in the previous section, it is required to evaluate the gradient and Hessian-vector product on the whole training dataset. However, these computations are prohibitive in the large-scale setting, i.e., when n and d are large. To address this issue, we present a stochastic variant of OASIS (see Appendix B for details) that only considers a small mini-batch of training data in each iteration. The Stochastic OASIS chooses sets Ik, Jk ⊂ [n] randomly and independently, and the new iterate is computed as: wk+1 = wk − ηkD̂k−1∇FIk(wk), where ∇FIk(wk) = 1|Ik| ∑ i∈Ik ∇Fi(wk) and D̂k is the truncated variant of Dk = β2Dk−1 + (1 − β2) diag(zk ∇2FJk(wk)zk) with ∇2FJk(wk) = 1|Jk| ∑ j∈Jk ∇ 2Fj(wk).
3.3 WARMSTARTING AND COMPLEXITY OF OASIS
Both deterministic and stochastic variants of our methodology share the need to obtain an initial estimate of D0. The importance of this is illustrated by the rule (10), which regulates the choice of
ηk, which is highly dependent on Dk. In order to have a better approximation of D0, we propose to sample some predefined number of Hutchinson’s estimates before the training process.
The main overhead of our methodology is the Hessian-vector product used in Hutchinson’s method for approximating the Hessian diagonal. With the current advanced hardware and packages, this computation is not a bottleneck anymore. To be more specific, the Hessian-vector product can be efficiently calculated by two rounds of back-propagation. We also present how to calculate the Hessian-vector product efficiently for various well-known machine learning tasks in Appendix B.

heading: 4 THEORETICAL ANALYSIS
text: In this section, we present our theoretical results for OASIS. We show convergence guarantees for different settings of learning rates, i.e., (i) adaptive learning rate, (ii) fixed learning rate, and (iii) with line search. Before the main theorems are presented, we state the following assumptions and lemmas that are used throughout this section. For brevity, we present the theoretical results using line search in Appendix A. The proofs and other auxiliary lemmas are in Appendix A.
Assumption 4.1. (Convex). The function F is convex. Assumption 4.2. (L−smooth). The gradients of F are L−Lipschitz continuous for all w ∈ Rd, i.e., ∃L > 0 such that ∀w, w′ ∈ Rd,F (w) ≤ F (w′) + 〈∇F (w′), w − w′〉+ L2 ‖w − w
′‖2. Assumption 4.3. The function F is twice continuously differentiable. Assumption 4.4. (µ−strongly convex). The function F is µ−strongly convex, i.e., there exists a constant µ > 0 such that ∀w, w′ ∈ Rd, F (w) ≥ F (w′) + 〈∇F (w′), w − w′〉+ µ2 ‖w − w ′‖2.
Here is a lemma regarding the bounds for Hutchinson’s approximation and the diagonal differences.
Lemma 4.5. (Bound on change of Dk). Suppose that Assumption 4.2 holds, then i) |(vk)i| ≤ Γ ≤√ dL, where vk = zk ∇2F (wk)zk; ii) ∃ δ ≤ 2(1− β2)Γ such that ∀k : ‖Dk+1 −Dk‖∞ ≤ δ.

heading: 4.1 ADAPTIVE LEARNING RATE
text: Here, we present theoretical convergence results for the case with adaptive learning rate using (10).
Theorem 4.6. Suppose that Assumptions 4.1, 4.2 and 4.3 hold. Let {wk} be the iterates generated by Algorithm OASIS, then we have: F (ŵk) − F ∗ ≤ LCk + 2L(1 − β2)Γ Qk k , where C = 2‖w1−w∗‖2D̂0+‖w1−w0‖ 2 D̂0
2 +2η1θ1(F (w0)−F (w ∗)), andQk = ∑k i=1( 2ηiθi+α 2α ‖wi−1−wi‖
2 + L2ηiθi+α
α ‖wi − w ∗‖2).
Remark 4.7. The following remarks are made regarding Theorem 4.6:
1. If β2 = 1, α = 1, and D0 = Id, the results in Theorem 4.6 completely match with (Malitsky & Mishchenko, 2020). 2. By considering an extra assumption as in (Reddi et al., 2019; Duchi et al., 2011) regarding bounded iterates, i.e., ‖wk − w∗‖2 ≤ B ∀k ≥ 0, one can easily show the convergence of OASIS to the neighborhood of optimal solution(s).
The following lemma provides the bounds for the adaptive learning rate for smooth and stronglyconvex loss functions. The next theorem provides the linear convergence rate for the latter setting.
Lemma 4.8. Suppose that Assumptions 4.2, 4.3, and 4.4 hold, then ηk ∈ [ α2L , Γ 2µ ].
Theorem 4.9. Suppose that Assumptions 4.2, 4.3, and 4.4 hold, and let w∗ be the unique solution for (1). Let {wk} be the iterates generated by Algorithm 1. Then, for all k ≥ 0 and β2 ≥ max{1− α 4µ4
4L2Γ2(α2µ2+LΓ2) , 1− α3µ3 4LΓ(2α2µ2+L3Γ2)} we have: Ψ k+1 ≤ (1− α 2 2Γ2κ2 )Ψ k, where Ψk+1 =
‖wk+1 − w∗‖2D̂k + 1 2‖wk+1 − wk‖ 2 D̂k + 2ηk(1 + θk)(F (wk)− F (w∗)).

heading: 4.2 FIXED LEARNING RATE
text: Here, we provide theoretical results for fixed learning rate for deterministic and stochastic regimes.
Remark 4.10. For any k ≥ 0, we have αI D̂k ΓI where 0 < α ≤ Γ.

heading: 4.2.1 DETERMINISTIC REGIME
text: Strongly Convex. The following theorem provides the linear convergence for the smooth and strongly-convex loss functions with fixed learning rate. Theorem 4.11. Suppose that Assumptions 4.2, 4.3, and 4.4 hold, and let F ∗ = F (w∗), where w∗ is the unique minimizer. Let {wk} be the iterates generated by Algorithm 1, where 0 < ηk = η ≤ α2
LΓ ,
and w0 is a starting point. Then, for all k ≥ 0 we have: F (wk)− F ∗ ≤ (1− ηµΓ ) k[F (w0)− F ∗].
Nonconvex. The following theorem provides the convergence to the stationary points for the nonconvex setting with fixed learning rate.
Assumption 4.12. The function F (.) is bounded below by a scalar F̂ . Theorem 4.13. Suppose that Assumptions 4.2, 4.3, and 4.12 hold. Let {wk} be the iterates generated by Algorithm , where 0 < ηk = η ≤ α2
LΓ , and w0 is a starting point. Then, for all T > 1 we have:
1 T ∑T k=1‖∇F (wk)‖2 ≤ 2Γ[F (w0)−F̂ ] ηT T→∞−−−−→ 0. (11)

heading: 4.2.2 STOCHASTIC REGIME
text: Here, we use EIk [.] to denote conditional expectation given wk, and E[.] to denote the full expectation over the full history. The following standard assumptions as in (Bollapragada et al., 2019; Berahas et al., 2016) are considered for this section. Assumption 4.14. There exist a constant γ such that EI [‖∇FI(w)−∇F (w)‖2] ≤ γ2. Assumption 4.15. There exist a constant σ2 <∞ such that EI [‖∇FI(w∗)‖2] ≤ σ2. Assumption 4.16. ∇FI(w) is an unbiased estimator of the gradient, i.e., EI [∇FI(w)] = ∇F (w).
Strongly Convex. The following theorem presents the convergence to the neighborhood of the optimal solution for the smooth and strongly-convex case in the stochastic setting. Theorem 4.17. Suppose that Assumptions 4.2, 4.3, 4.4, 4.15 and 4.16 hold. Let {wk} be the iterates generated by Algorithm 1 with ηk =η ∈ (0, α 2µ ΓL2 ), then, for all k ≥ 0,
E[F (wk)− F ∗] ≤ (1− c)k (F (w0)− F ∗) + η 2Lσ2
cα2 , (12)
where c = 2ηµΓ − 2η2L2 α2 ∈ (0, 1). Moreover, if ηk =η ∈ (0, α2µ 2ΓL2 ) then E[F (wk)− F ∗] ≤ ( 1− ηµΓ )k (F (w0)− F ∗) + ηΓLσ 2 α2µ . (13)
Nonconvex. The next theorem provides the convergence to the stationary point for the nonconvex loss functions in the stochastic regime. Theorem 4.18. Suppose that Assumptions 4.2, 4.3, 4.12, 4.14 and 4.16 hold. Let {wk} be the iterates generated by Algorithm 1, where 0 < ηk = η ≤ η 2
LΓ , and w0 is the starting point. Then, for all k ≥ 0, E [
1 T ∑T−1 k=0 ‖∇F (wk)‖2 ] ≤ 2Γ[F (w0)−F̂ ]ηT + ηΓγ2L α2 T→∞−−−−→ ηΓγ 2L α2 .
The previous two theorems provide the convergence to the neighborhood of the stationary points. One can easily use either a variance reduced gradient approximation or a decaying learning rate strategy to show the convergence to the stationary points (in expectation). OASIS’s analyses are similar to those of limited-memory quasi-Newton approaches which depends on λmax and λmin (largest and smallest eigenvalues) of the preconditioning matrix (while in (S)GD λmax= λmin = 1). These methods in theory are not better than GD-type methods. In practice, however, they have shown their strength.

heading: 5 EMPIRICAL RESULTS
text: In this section, we present empirical results for several machine learning problems to show that our OASIS methodology outperforms state-of-the-art first- and second-order methods in both deterministic and stochastic regimes. We considered: (1) deterministic `2-regularized logistic regression (strongly convex); (2) deterministic nonlinear least squares (nonconvex), and we report results on 2 standard machine learning datasets rcv1 and ijcnn13; and (3) image classification tasks on
3https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/
0 100 200 300 400 Number of Effective Passes
0.3 0.4 0.5 0.6 0.7 0.8 0.9 Te st A cc ur ac y
ijcnn1
300 350 4000.913
0.914
0.915
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
10 11 10 9 10 7 10 5 10 3 10 1 F( w ) F *
ijcnn1
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
0.3 0.4 0.5 0.6 0.7 0.8 0.9 Te st A cc ur ac y
rcv1
300 350 4000.944
0.945
0.946
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
10 13 10 11 10 9 10 7 10 5 10 3 10 1 101 F( w ) F *
rcv1 AdGD AdaHessian OASIS
Figure 3: Comparison of optimality gap and Test Accuracy for different algorithms on Logistic Regression Problems.
0 100 200 300 400 Number of Effective Passes
0.3 0.4 0.5 0.6 0.7 0.8 0.9 Te st A cc ur ac y
ijcnn1
300 350 4000.9025
0.9075
0.9125
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
100
2 × 100
F( w
)
ijcnn1
300 350 4000.985
0.990
0.995
1.000
1.005
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
0.5
0.6
0.7
0.8
0.9
Te st
A cc
ur ac
y
rcv1
AdGD AdaHessian OASIS
0 100 200 300 400 Number of Effective Passes
100
6 × 10 1
F( w
)
rcv1 AdGD AdaHessian OASIS
Figure 4: Comparison of objective function (F (w)) and Test Accuracy for different algorithms on Non-linear Least Square Problems.
MNIST, CIFAR10, and CIFAR100 datasets on standard network structures.In the interest of space, we report only a subset of the results in this section. The rest can be found in Appendix C.
To be clear, we compared the empirical performance of OASIS with algorithms with diagonal preconditioners. In the deterministic regime, we compared the performance of OASIS with AdGD (Malitsky & Mishchenko, 2020) and AdaHessian (Yao et al., 2020). Further, for the stochastic regime, we provide experiments comparing SGD (Robbins & Monro, 1951), Adam (Kingma & Ba, 2014), AdamW (Loshchilov & Hutter, 2017), and AdaHessian. For the logistic regression problems, the regularization parameter was chosen from the set λ ∈ { 110n , 1 n , 10 n }. It is worth highlighting that we ran each method for each of the following experiments from 10 different random initial points. Moreover, we separately tuned the hyperparameters for each algorithm, if needed. See Appendix C for details. The proposed OASIS is robust with respect to different choices of hyperparameters, and it has a narrow spectrum of changes (see Appendix C).
Logistic Regression. We considered `2-regularized logistic regression problems, F (w) = 1 n ∑n i=1 log(1 + e −yixTi w) + λ2 ‖w‖ 2. Figure 3 shows the performance of the methods in terms of optimality gap and test accuracy versus number of effective passes (number of gradient and Hessian-vector evaluations). As is clear, the performance of OASIS (with adaptive learning rate and without any hyperparameter tuning) is on par or better than that of the other methods.
Non-linear Least Square. We considered non-linear least squares problems (described in Xu et al. (2020)): F (w) = 1n ∑n i=1(yi− 1/(1 + e−x T i w))2. Figure 4 shows that our OASIS Algorithm always outperforms the other methods in terms of training loss function and test accuracy. Moreover, the behaviour of OASIS is robust with respect to the different initial points.
Image Classification. We illustrate the performance of OASIS on standard bench-marking neural network training tasks: MNIST, CIFAR10, and CIFAR100. The results for MNIST and the details of the problems are given in Appendix C. We present the results regarding CIFAR10 and CIFAR100. CIFAR10. We use standard ResNet-20 and ResNet-32 (He et al., 2015) architectures for comparing the performance of OASIS with SGD, Adam, AdamW and AdaHessian4. Specifically, we report 3 variants of OASIS: (i) adaptive learning rate, (ii) fixed learning rate (without first moment) and (iii) fixed learning rate with gradient momentum tagged with “Adaptive LR,” “Fixed LR,” and “Momentum,” respectively. For settings with fixed learning rate, no warmstarting is used in order to obtain an initial D0 approximation. For the case with adaptive learning case, we used the warmstarting strategy to approximate the initial Hessian diagonal. More details regarding the exact parameter values and hyperparameter search can be found in the Appendix C. The results on CIFAR10 are shown in the Figure 5 (the left and middle columns) and Table 1. As is clear, the simplest variant of OASIS with fixed learning rate achieves significantly better results, as compared to Adam, and a performance comparable to SGD. For the variant with an added momentum, we get similar accuracy as AdaHessian, while getting better or the same loss values, highlighting the advantage of using different preconditioning schema. Another important observation is that OASISAdaptive LR, without too much tuning efforts, has better performance than Adam with sensitive
4Note that we follow the same Experiment Setup as in (Yao et al., 2020), and the codes for other algorithms and structures are brought from https://github.com/amirgholami/adahessian.
hyperparameters. All in all, the performance of OASIS variants is on par or better than the other state-of-the-art methods especially SGD. As we see from Figure 5 (left and middle columns), the lack of momentum produces a slow, noisy training curve in the initial stages of training, while OASIS with momentum works better than the other two variants in the early stages. All three variants of OASIS get satisfactory results in the end of training. More results and discussion regarding CIFAR10 dataset are in Appendix C.
Table 1: Results of ResNet-20/32 on CIFAR10
Setting ResNet-20 ResNet-32
SGD 92.02 ± 0.22 92.85 ± 0.12 Adam 90.46 ± 0.22 91.30 ± 0.15 AdamW 91.99 ± 0.17 92.58 ± 0.25 AdaHessian 92.03 ± 0.10 92.71 ± 0.26 OASIS- Adaptive LR 91.20 ± 0.20 92.61 ± 0.22 OASIS- Fixed LR 91.96 ± 0.21 93.01 ± 0.09 OASIS- Momentum 92.01 ± 0.19 92.77 ± 0.18
Table 2: Results of ResNet-18 on CIFAR100.
Setting ResNet-18
SGD 76.57 ± 0.24 Adam 73.40 ± 0.31 AdamW 72.51 ± 0.76 AdaHessian 75.71 ± 0.47 OASIS- Adaptive LR 76.93 ± 0.22 OASIS- Fixed LR 76.28 ± 0.21 OASIS- Momentum 76.89 ± 0.34
CIFAR-100. We use the hyperparameter settings obtained by training on CIFAR10 on ResNet-20/32 to train CIFAR100 on ResNet-18 network structure.5 We similarly compare the performance of our method and its variants with SGD, Adam, AdamW and AdaHessian. The results are shown in Figure 5 (right column) and Table 2. In this setting, without any hyperparameter tuning, fully adaptive version of our algorithm immediately produces results surpassing other state-of-the-art methods especially SGD.

heading: 6 FINAL REMARKS
text: This paper presents a fully adaptive optimization algorithm for empirical risk minimization. The search direction uses the gradient information, well-scaled with a novel Hessian diagonal approximation, which itself can be calculated and stored efficiently. In addition, we do not need to tune the learning rate, which instead is automatically updated based on a low-cost approximation of the Lipschitz smoothness parameter. We provide comprehensive theoretical results covering standard optimization settings, including convex, strongly convex and nonconvex; and our empirical results highlight the efficiency of OASIS in large-scale machine learning problems. Future research avenues include: (1) deriving the theoretical results for stochastic regime with adaptive learning rate; (2) employing variance reduction schemes in order to reduce further the noise in the gradient and Hessian diagonal estimates; and (3) providing a more extensive empirical investigation on other demanding machine learning problems such as those from natural language processing and recommendation system (such as those from the original AdaHessian paper (Yao et al., 2020)).
5https://github.com/uoguelph-mlrg/Cutout.

heading: ETHICS STATEMENT
text: This work presents a new algorithm for training machine learning models. We do not foresee any ethical concerns. All datasets used in this work are from the public domain and are commonly used benchmarks in ML papers.

heading: REPRODUCIBILITY STATEMENT
text: We uploaded all the codes used to make all the experiments presented in this paper. We have used random seeds to ensure that one can start optimizing the ML models from the same initial starting point as was used in the experiments. We have used only datasets that are in the public domain, and one can download them from the following website https://www.csie.ntu.edu.tw/ ˜cjlin/libsvmtools/datasets/. After acceptance, we will include a link to the GitHub repository where we will host the source codes.

heading: A THEORETICAL RESULTS AND PROOFS

heading: A.1 ASSUMPTIONS
text: Assumption 4.1. (Convex). The function F is convex, i.e., ∀w, w′ ∈ Rd, F (w) ≥ F (w′) + 〈∇F (w′), w − w′〉. (14)
Assumption 4.2. (L−smooth). The gradients of F are L−Lipschitz continuous for all w ∈ Rd, i.e., there exists a constant L > 0 such that ∀w, w′ ∈ Rd,
F (w) ≤ F (w′) + 〈∇F (w′), w − w′〉+ L 2 ‖w − w′‖2. (15)
Assumption 4.3. The function F is twice continuously differentiable. Assumption 4.4. (µ−strongly convex). The function F is µ−strongly convex, i.e., there exists a constant µ > 0 such that ∀w, w′ ∈ Rd,
F (w) ≥ F (w′) + 〈∇F (w′), w − w′〉+ µ 2 ‖w − w′‖2. (16)
Assumption 4.12. The function F (.) is bounded below by a scalar F̂ . Assumption 4.14. There exist a constant γ such that EI [‖∇FI(w)−∇F (w)‖2] ≤ γ2. Assumption 4.16. ∇FI(w) is an unbiased estimator of the gradient, i.e., EI [∇FI(w)] = ∇F (w), where the samples I are drawn independently.

heading: A.2 PROOF OF LEMMA 4.5
text: Lemma 4.5. (Bound on change ofDk). Suppose that Assumption 4.2 holds, i.e., ∀w : ∇2f(w) LI , then
1. |(vk)i| ≤ Γ ≤ √ dL, where vk = zk ∇2F (wk)zk.
2. there ∃δ ≤ 2(1− β2)Γ such that ‖Dk+1 −Dk‖∞ ≤ δ, ∀k.
Proof. By Assumption 4.2, we have that ‖∇2F (w)‖2 ≤ L and hence
‖vk‖∞ ≤ ‖∇2F (w)‖∞ ≤ √ d‖∇2F (w)‖2 ≤ √ dL,
which finishes the proof of case 1.
Now, from (8) we can derive
Dk+1 −Dk (8) = (β2 − 1)Dk + (1− β2) zk ∇2F (wk)zk
and hence
‖Dk+1 −Dk‖∞ = (1− β2)‖Dk − zk ∇2F (wk)zk‖∞ ≤ (1− β2)‖Dk − vk‖∞ ≤ (1− β2) (‖Dk‖∞ + ‖vk‖∞) ≤ 2(1− β2)Γ.

heading: A.3 LEMMA REGARDING SMOOTHNESS WITH WEIGHTED NORM
text: In the following, we present a lemma for smoothness with weighted norm. Theorem 2.1.5 in (Nesterov et al., 2018) provides the same analysis for any norm, and the following lemma can be seen as a special case of Theorem 2.1.5 (Nesterov et al., 2018) with respect to the weighted Euclidean norm. In the following, we provide the proof for completeness. Lemma A.1. (Smoothness with norm D) Suppose that Assumptions 4.1 and 4.2 hold, and D 0, then we have:
‖∇F (x)−∇F (y)‖∗D ≤ L̃‖x− y‖D, (17)
where L̃ = L
λmin(D) .
Proof. By the equality L̃ = L
λmin(D) , we conclude that L̃D LI which results in:
F (y) ≤ F (x) + 〈∇F (x), y − x〉+ L̃ 2 ‖x− y‖2D. (18)
Extending proof of Theorem 2.1.5. in (Nesterov, 2013) we define φ(y) = F (y)− 〈∇F (x), y〉. Then, clearly, x ∈ arg minφ(y) and (18) is still valid for φ(y). Therefore
φ(x) ≤ φ(y − D −1
L̃ ∇φ(y)),
φ(x) (18) ≤ φ(y) + 〈∇φ(y),−D −1
L̃ ∇φ(y)〉+ L̃ 2 ‖ − D
−1
L̃ ∇φ(y)‖2D,
F (x)− 〈∇F (x), x〉 ≤ F (y)− 〈∇F (x), y〉 − 1 L̃ 〈∇φ(y), D−1∇φ(y)〉+ 1 2L̃ ‖D−1∇φ(y)‖2D,
F (x) ≤ F (y) + 〈∇F (x), x− y〉 − 1 2L̃ (‖∇φ(y)‖∗D)2, F (x) ≤ F (y) + 〈∇F (x), x− y〉 − 1 2L̃ (‖∇F (y)−∇F (x)‖∗D)2.
Thus F (x) + 〈∇F (x), y − x〉+ 1
2L̃ (‖∇F (y)−∇F (x)‖∗D)2 ≤ F (y). (19)
Adding (19) with itself with x swapped with y we obtain
1 L̃ (‖∇F (y)−∇F (x)‖∗D)2 ≤ 〈∇F (y)−∇F (x), y − x〉
= 〈D−1(∇F (y)−∇F (x)), D(y − x)〉 ≤ ‖∇F (y)−∇F (x))‖∗D‖y − x‖D,
which implies that
‖∇F (x)−∇F (y)‖∗D ≤ L̃‖x− y‖D. (20)

heading: A.4 PROOF OF THEOREM 4.6
text: Lemma A.2. Let f : Rd → R be a convex function, and x∗ is one of the optimal solutions for (1). Then, for the sequence of {wk} generated by Algorithm 1 we have:
‖wk+1 − w∗‖2D̂k+ 1 2 ‖wk − wk+1‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))+
2(1− β2)Γ (
( ηkθk α + 1 2 )‖wk−1 − wk‖2 + ( L2ηkθk α
+ 1)‖wk − w∗‖2 ) . (21)
Proof. We extend the proof in (Malitsky & Mishchenko, 2020). We have
‖wk+1 − w∗‖2D̂k = ‖wk+1 − wk + wk − w ∗‖2 D̂k
= ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2〈wk+1 − wk, D̂k(wk − w∗)〉 = ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2ηk〈∇F (wk), w∗ − wk〉 ≤ ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2ηk(F (w ∗)− F (wk)) (22) = ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k − 2ηk(F (wk)− F (w∗)),
where the third equality comes from the OASIS’s step and the inequality follows from convexity of F (w). Now, let’s focus on ‖wk+1 − wk‖2D̂k . We have
‖wk+1 − wk‖2D̂k =2‖wk+1 − wk‖ 2 D̂k − ‖wk+1 − wk‖2D̂k
=2〈−ηkD̂k −1 ∇F (wk), D̂k(wk+1 − wk)〉 − ‖wk+1 − wk‖2D̂k =− 2ηk〈∇F (wk), wk+1 − wk〉 − ‖wk+1 − wk‖2D̂k =− 2ηk〈∇F (wk)−∇F (wk−1), wk+1 − wk〉 − 2ηk〈∇F (wk−1), wk+1 − wk〉− ‖wk+1 − wk‖2D̂k =2ηk〈∇F (wk)−∇F (wk−1), wk − wk+1〉+ 2ηk〈∇F (wk−1), wk − wk+1〉− ‖wk+1 − wk‖2D̂k .
Now,
2ηk〈∇F (wk)−∇F (wk−1), wk − wk+1〉 ≤ 2ηk‖∇F (wk)−∇F (wk−1)‖∗D̂k‖wk − wk+1‖D̂k (10)
≤ ‖wk − wk−1‖D̂k‖wk − wk+1‖D̂k
≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k ,
where the first inequality comes from Cauchy-Schwarz and the third one follows Young’s inequality. Further,
〈∇F (wk−1), wk − wk+1〉 = 1
ηk−1 〈D̂k−1(wk−1 − wk), wk − wk+1〉
= 1
ηk−1 〈D̂k−1(wk−1 − wk), ηkD̂k−1∇F (wk)〉
= ηk ηk−1 (wk−1 − wk)T D̂k−1D̂k−1∇F (wk) = ηk ηk−1 (wk−1 − wk)T∇F (wk)+
ηk ηk−1
(wk−1 − wk)T ( D̂k−1D̂k −1 − I ) ∇F (wk), (23)
where the first two qualities are due the OASIS’s update step. The second term in the above equality can be upperbounded as follows:
(wk−1 − wk)T ( D̂k−1D̂k −1 − I ) ∇F (wk) ≤ (1− β2) 2Γ
α · ‖wk−1 − wk‖ · ‖∇F (wk)‖, (24)
where multiplier on the left is obtained via the bound on ‖ · ‖∞ norm of the diagonal matrix D̂k−1D̂k −1 − I:
‖D̂k−1D̂k−1 − I‖∞ = max i ∣∣∣(D̂k−1D̂k−1 − I) i ∣∣∣ = max i ∣∣(D̂k−1 − D̂k)i∣∣∣∣(D̂k−1)i∣∣ ≤ (1− β2)2Γα , which also represents a bound on the operator norm of the same matrix difference. Next we can use the Young’s inequality and∇F (w∗) = 0 to get
(1− β2) 2Γ
α · ‖wk−1 − wk‖ · ‖∇F (wk)‖ ≤ (1− β2)
Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2 ) . (25)
Therefore, we have
〈∇F (wk−1), wk − wk+1〉 ≤ ηk ηk−1 (wk−1 − wk)T∇F (wk)+
(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2).
Also,
‖wk+1 − wk‖2D̂k ≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k + 2ηkθk(wk−1 − wk) T∇F (wk)+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)− ‖wk+1 − wk‖2D̂k
≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k + 2ηkθk(F (wk−1)− F (wk))+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)− ‖wk+1 − wk‖2D̂k .
Finally, we have
‖wk+1 − w∗‖2D̂k ≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k + 2ηkθk(F (wk−1)− F (wk))+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)− ‖wk+1 − wk‖2D̂k
+ ‖wk − w∗‖2D̂k − 2ηk(F (wk)− F (w ∗)).
By simplifying the above inequality, we have:
‖wk+1 − w∗‖2D̂k+ 1 2 ‖wk − wk+1‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ ‖wk − w∗‖2D̂k + 1 2 ‖wk − wk−1‖2D̂k + 2ηkθk(F (wk−1)− F (w ∗))+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)
= ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)+
‖wk − w∗‖2D̂k−D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−D̂k−1
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))+
2ηkθk(1− β2) Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2)+
2(1− β2)Γ(‖wk − w∗‖2 + 1
2 ‖wk − wk−1‖2)
= ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))+
2(1− β2)Γ (ηkθk
α ‖wk−1 − wk‖2 + L2ηkθk α ‖wk − w∗‖2 + ‖wk − w∗‖2+
1 2 ‖wk − wk−1‖2
) .
Theorem 4.6. Suppose that Assumptions 4.1, 4.2 and 4.3 hold. Let {wk} be the iterates generated by Algorithm 1, then we have:
F (ŵk)− F ∗ ≤ LC
k + 2L(1− β2)Γ Qk k ,
where
C = ‖w1 − w∗‖2D̂0 + 1 2 ‖w1 − w0‖2D̂0 + 2η1θ1(F (w0)− F (w ∗))
Qk = k∑ i=1 ( ( ηiθi α + 1 2 )‖wi−1 − wi‖2 + ( L2ηiθi α + 1)‖wi − w∗‖2 ) .
Proof. By telescoping inequality (21) in Lemma A.2 we have:
‖wk+1 − w∗‖2D̂k+ 1 2 ‖wk − wk+1‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
+ 2 k−1∑ i=1 [ηi(1 + θi)− ηi+1θi+1](F (wk)− F (w∗))
≤ ‖w1 − w∗‖2D̂0 + 1 2 ‖w1 − w0‖2D̂0 + 2η1θ1(F (w0)− F (w ∗))︸ ︷︷ ︸ C +
2(1− β2)Γ k∑ i=1 ( ( ηiθi α + 1 2 )‖wi−1 − wi‖2 + ( L2ηiθi α + 1)‖wi − w∗‖2 )
︸ ︷︷ ︸ Qk
.
Moreover, by the rule for adaptive learning rule we know ηi(1 + θi)− ηi+1θi+1 ≥ 0, ∀i. Therefore, we have
2ηk(1 + θk)(F (wk)− F (w∗))+2 k−1∑ i=1 [ηi(1 + θi)− ηi+1θi+1](F (wk)− F (w∗))
≤ C + 2(1− β2)ΓQk.
By setting ŵ = ηk(1 + θk)wk +
∑k−1 i=1 (ηi(1 + θi)− ηi+1θi+1)wi
Sk , where Sk := ηk(1 + θk) +∑k−1
i=1 (ηi(1 + θi)− ηi+1θi+1), and by using Jensens’s inequality, we have:
F (ŵk)− F ∗ ≤ C
2Sk + (1− β2)Γ Qk Sk .
By the fact that ηk ≥ 1
2L thus
1 Sk ≤ 2L k , we have
F (ŵk)− F ∗ ≤ LC
k + 2L(1− β2)Γ Qk k .

heading: A.5 PROOF OF LEMMA 4.8
text: Lemma 4.8. Suppose that Assumptions 4.2, 4.3 and 4.4 hold, then ηk ∈ [ α
2L ,
Γ
2µ
] .
Proof. By (17) and the point that L̃ = L
λmin(D̂k) , we conclude that:
‖wk − wk−1‖D̂k 2‖∇F (wk)−∇F (wk−1)‖∗D̂k ≥ 1 2L̃ = λmin(D̂k) 2L ≥ α 2L .
Now, in order to find the upperbound for ‖wk−wk−1‖D̂k
2‖∇F (wk)−∇F (wk−1)‖∗ D̂k
, we use the following inequality
which comes from Assumption 4.4:
F (wk) ≥ F (wk−1) + 〈∇F (wk−1), (wk − wk−1)〉+ µ
2 ‖w − wk−1‖2. (26)
By setting µ̃ = µ
λmax(D̂k) , we conclude that µ̃D̂k µI which results in:
F (wk) ≥ F (wk−1) + 〈∇F (wk−1), (wk − wk−1)〉+ µ̃
2 ‖w − wk−1‖2D̂k . (27)
The above inequality results in:
µ̃‖w − wk−1‖2D̂k ≤ 〈∇F (wk)−∇F (wk−1), (wk − wk−1)〉
= 〈D̂k−1(∇F (wk)−∇F (wk−1)), D̂k(wk − wk−1)〉 ≤ ‖∇F (wk)−∇F (wk−1))‖∗D̂k‖wk − wk−1‖D̂k . (28)
Therefore, we obtain that
‖wk − wk−1‖D̂k 2‖∇F (wk)−∇F (wk−1)‖∗D̂k (28) ≤ 1 2µ̃ = λmax(D̂k) 2µ ≤ Γ 2µ ,
and therefore, by the update rule for ηk, we conclude that ηk ∈ [ α
2L ,
Γ
2µ
] .

heading: A.6 PROOF OF THEOREM 4.9
text: Lemma A.3. . Suppose Assumptions 4.2, 4.3 and 4.4 hold and let w∗ be the unique solution of (1). Then for (wk) generated by Algorithm 1 we have:
‖wk+1 − w∗‖2D̂k + 1 2 ‖wk+1 − wk‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))
+ ( (1− β2)Γ ( 1 +
2θkηk α
) − µηkθk ) ‖wk − wk−1‖2
+ ( (1− β2)Γ ( 2 +
2L2θkηk α
) − µηk ) ‖wk − w∗‖2.
Proof. By the update rule in Algorithm 1 we have:
‖wk+1 − w∗‖2D̂k = ‖wk+1 − wk + wk − w ∗‖2 D̂k
= ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2〈wk+1 − wk, D̂k(wk − w∗)〉 = ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2ηk〈∇F (wk), w∗ − wk〉 ≤ ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2ηk(F (w ∗)− F (wk)) (29) = ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k − 2ηk(F (wk)− F (w∗)),
where the inequality follows from convexity of F (w). By strong convexity of F (.) we have:
F (w∗) ≥ F (wk) + 〈∇F (wk), w∗ − wk〉+ µ
2 ‖wk − w∗‖2. (30)
In the lights of the strongly convex inequality we can change (29) as follows
‖wk+1 − w∗‖2D̂k (29),(30) ≤ ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k + 2ηk(F (w ∗)− F (wk)− µ 2 ‖wk − w∗‖2)
= ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k − 2ηk(F (wk)− F (w∗))− µηk‖wk − w∗‖2.
(31)
Now, let’s focus on ‖wk+1 − wk‖2D̂k . We have
‖wk+1 − wk‖2D̂k =2‖wk+1 − wk‖ 2 D̂k − ‖wk+1 − wk‖2D̂k
=2〈−ηkD̂k −1 ∇F (wk), D̂k(wk+1 − wk)〉 − ‖wk+1 − wk‖2D̂k =− 2ηk〈∇F (wk), wk+1 − wk〉 − ‖wk+1 − wk‖2D̂k =− 2ηk〈∇F (wk)−∇F (wk−1), wk+1 − wk〉 − 2ηk〈∇F (wk−1), wk+1 − wk〉− ‖wk+1 − wk‖2D̂k =2ηk〈∇F (wk)−∇F (wk−1), wk − wk+1〉+ 2ηk〈∇F (wk−1), wk − wk+1〉− ‖wk+1 − wk‖2D̂k . (32)
Now, 2ηk〈∇F (wk)−∇F (wk−1), wk − wk+1〉 ≤ 2ηk‖∇F (wk)−∇F (wk−1)‖∗D̂k‖wk − wk+1‖D̂k
≤ ‖wk − wk−1‖D̂k‖wk − wk+1‖D̂k ≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k , (33)
where the first inequality is due to Cauchy–Schwarz inequality, and the second inequality comes from the choice of learning rate such that ηk ≤ ‖wk−wk−1‖D̂k
2‖∇F (wk)−∇F (wk−1)‖∗ D̂k
. By plugging (33) into (32) we
obtain:
‖wk+1 − wk‖2D̂k (32) ≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k + 2ηk〈∇F (wk−1), wk − wk+1〉−
‖wk+1 − wk‖2D̂k . (34)
Now, we can summarize that
‖wk+1 − w∗‖2D̂k (31) ≤ ‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k − 2ηk(F (wk)− F (w∗))− µηk‖wk − w∗‖2
(34)
≤ 1 2 ‖wk − wk−1‖2D̂k + 1 2 ‖wk − wk+1‖2D̂k + 2ηk〈∇F (wk−1), wk − wk+1〉−
‖wk+1 − wk‖2D̂k + ‖wk − w ∗‖2 D̂k − 2ηk(F (wk)− F (w∗))− µηk‖wk − w∗‖2
= 1
2 ‖wk − wk−1‖2D̂k −
1 2 ‖wk − wk+1‖2D̂k + 2ηk〈∇F (wk−1), wk − wk+1〉+
‖wk − w∗‖2D̂k − 2ηk(F (wk)− F (w ∗))− µηk‖wk − w∗‖2
= 1
2 ‖wk − wk−1‖2D̂k−1 + ‖wk − w ∗‖2 D̂k−1 − 1 2 ‖wk − wk+1‖2D̂k+
2ηk〈∇F (wk−1), wk − wk+1〉 − 2ηk(F (wk)− F (w∗))+ 1 2 ‖wk − wk−1‖2D̂k−D̂k−1 + ‖wk − w ∗‖2 D̂k−D̂k−1 − µηk‖wk − w∗‖2. (35)

heading: Next, let us bound 〈∇F (wk−1), wk − wk+1〉.
text: By the OASIS’s update rule, wk+1 = wk − ηkD̂k−1∇F (wk), we have,
〈∇F (wk−1), wk − wk+1〉 = 1
ηk−1 〈D̂k−1(wk−1 − wk), wk − wk+1〉
= 1
ηk−1 〈D̂k−1(wk−1 − wk), ηkD̂k−1∇F (wk)〉
= ηk ηk−1 (wk−1 − wk)T D̂k−1D̂k−1∇F (wk) = ηk ηk−1 (wk−1 − wk)T∇F (wk)+
ηk ηk−1
(wk−1 − wk)T ( D̂k−1D̂k −1 − I ) ∇F (wk). (36)
The second term in the above equality can be bounded from above as follows:
(wk−1 − wk)T ( D̂k−1D̂k −1 − I ) ∇F (wk) ≤ (1− β2) 2Γ
α · ‖wk−1 − wk‖ · ‖∇F (wk)‖, (37)
where multiplier on the left is obtained via the bound on ‖ · ‖∞ norm of the diagonal matrix D̂k−1D̂k −1 − I:
‖D̂k−1D̂k−1 − I‖∞ = max i ∣∣∣(D̂k−1D̂k−1 − I) i ∣∣∣ = max i ∣∣(D̂k−1 − D̂k)i∣∣∣∣(D̂k−1)i∣∣ ≤ (1− β2)2Γα , which also represents a bound on the operator norm of the same matrix difference. Next we can use the Young’s inequality and∇F (w∗) = 0 to get
(1− β2) 2Γ
α · ‖wk−1 − wk‖ · ‖∇F (wk)‖ ≤ (1− β2)
Γ
α
( ‖wk−1 − wk‖2 + L2‖wk − w∗‖2 ) . (38)
By the inequalities (72), (36), (37) and (38), and the definition θk = ηk ηk−1 we have:
‖wk+1 − w∗‖2D̂k + 1 2 ‖wk+1 − wk‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))
+
( 1
2 ‖wk − wk−1‖2D̂k−D̂k−1 + (1− β2) 2Γθkηk α ‖wk − wk−1‖2 − µηkθk‖wk − wk−1‖2 ) + ( ‖wk − w∗‖2D̂k−D̂k−1 + (1− β2) 2ΓL2θkηk α ‖wk − w∗‖2 − µηk‖wk − w∗‖2 )
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))
+ ( (1− β2)Γ ( 1 +
2θkηk α
) − µηkθk ) ‖wk − wk−1‖2
+ ( (1− β2)Γ ( 2 +
2L2θkηk α
) − µηk ) ‖wk − w∗‖2.
Theorem 4.9. Suppose that Assumptions 4.2, 4.4, and 4.3 hold and let w∗ be the unique solution for (1). Let {wk} be the iterates generated by Algorithm 1. Then, for all k ≥ 0 we have: If β2 ≥ max{1− α4µ4
4L2Γ2(α2µ2 + LΓ2) , 1− α
3µ3
4LΓ(2α2µ2 + L3Γ2) }
Ψk+1 ≤ (1− α 2
2Γ2κ2 )Ψk, (39)
where
Ψk+1 = ‖wk+1 − w∗‖2D̂k + 1 2 ‖wk+1 − wk‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗)).
Proof. By Lemma A.3, we have:
‖wk+1 − w∗‖2D̂k + 1 2 ‖wk+1 − wk‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ ‖wk − w∗‖2D̂k−1 + 1 2 ‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗))
+ ( (1− β2)Γ ( 1 +
2θkηk α
) − µηkθk ) ‖wk − wk−1‖2
+ ( (1− β2)Γ ( 2 +
2L2θkηk α
) − µηk ) ‖wk − w∗‖2.
Lemma 4.8 gives us ηk ∈ [ α
2L ,
Γ
2µ ] , so we can choose large enough β2 ∈ (0, 1) such that(
(1− β2)Γ (
1 + 2θkηk α
) − µηkθk ) ≤ 0,(
(1− β2)Γ ( 2 + 2L2θkηk
α
) − µηk ) ≤ 0.
In other words, we have the following bound for β2:
β2 ≥ max{1− α4µ4 2L2Γ2(α2µ2 + LΓ2) , 1− α 3µ3 2LΓ(2α2µ2 + L3Γ2) }. (40)
However, we can push it one step further, by requiring a stricter inequality to hold, to get a recursion which can allows us to derive a linear convergence rate as follows
( (1− β2)Γ ( 1 +
2θkηk α
) − µηkθk ) ≤ −1
2 µηkθk,( (1− β2)Γ ( 2 + 2L2θkηk
α
) − µηk ) ≤ −1
2 µηk,
which requires a correspondingly stronger bound on β2:
β2 ≥ max{1− α4µ4 4L2Γ2(α2µ2 + LΓ2) , 1− α 3µ3 4LΓ(2α2µ2 + L3Γ2) }. (41)
Combining this with the condition for ηk and definition for θk, we obtain
( (1− β2)Γ ( 1 +
2θkηk α
) − µηkθk ) ‖wk − wk−1‖2 ≤ − 1
2 µηkθk‖wk − wk−1‖2
≤ − α 2
4Γ2κ2 ‖wk − wk−1‖2D̂k−1 ,
and ( (1− β2)Γ ( 2 +
2L2θkηk α
) − µηk ) ‖wk − w∗‖2 ≤ − 1
2 µηk‖wk − w∗‖2
≤ − α 4Γκ ‖wk − w∗‖2D̂k−1 .
where κ = Lµ . Using it to supplement the main statement of the lemma, we get
‖wk+1 − w∗‖2D̂k + 1 2 ‖wk+1 − wk‖2D̂k + 2ηk(1 + θk)(F (wk)− F (w ∗))
≤ (1− α 4Γκ )‖wk − w∗‖2D̂k−1 + 1 2 (1− α
2
2Γ2κ2 )‖wk − wk−1‖2D̂k−1 + 2ηkθk(F (wk−1)− F (w ∗)).
Mirroring the result in (Malitsky & Mishchenko, 2020) we get a contraction in all terms, since for function value differences we also can show
ηkθk ηk(1 + θk) = 1− ηk ηk(1 + θk) ≤ 1− α 2Γκ .

heading: A.7 PROOF OF THEOREM 4.11
text: Theorem 4.11. Suppose that Assumptions 4.2, 4.4, and 4.3 hold, and let F ∗ = F (w∗) where w∗ is the unique minimizer. Let {wk} be the iterates generated by Algorithm 1, where 0 < ηk = η ≤ α2
LΓ ,
and w0 is a starting point. Then, for all k ≥ 0 we have:
F (wk)− F ∗ ≤ (1− ηµ
Γ )k[F (w0)− F ∗]. (42)
Proof. By smoothness of F (.) we have:
F (wk+1) = F (wk − ηD̂−1k ∇f(wk))
≤ F (wk) +∇F (wk)T (−ηD̂−1k ∇f(wk)) + L 2 ‖ηD̂−1k ∇f(wk)‖ 2
≤ F (wk)− η
Γ ‖∇F (wk)‖2 +
η2L 2α2 ‖∇F (wk)‖2
= F (wk)− η( 1 Γ − ηL 2α2 )‖∇F (wk)‖2 (43) ≤ F (wk)− η 1
2Γ ‖∇F (wk)‖2, (44)
where the first inequality comes from Assumption 4.2, and the second inequality is due to Remark 4.10, and finally, the last inequality is by the choice of η. Since F (.) is strongly convex, we have 2µ(F (wk)− F ∗) ≤ ‖∇F (wk)‖2, and therefore,
F (wk+1) ≤ F (wk)− ηµ
Γ (F (wk)− F ∗).
Also, we have: F (wk+1)− F ∗ ≤ (1− ηµ
Γ )(F (wk)− F ∗).

heading: A.8 PROOF OF THEOREM 4.13
text: Theorem 4.13. Suppose that Assumptions 4.3, 4.12 and 4.2 hold. Let {wk} be the iterates generated by Algorithm 1, where 0 < ηk = η ≤ α2
LΓ , and w0 is a starting point. Then, for all T > 1 we have:
1
T T∑ k=1 ‖∇F (wk)‖2 ≤ 2Γ[F (w0)− F̂ ] ηT T→∞−−−−→ 0. (45)
Proof. By starting from (44), we have:
F (wk+1 ≤ F (wk)− η
2Γ ‖∇F (wk)‖2. (46)
By summing both sides of the above inequality from k = 0 to T − 1 we have: T−1∑ k=0 (F (wk+1)− F (wk)) ≤ − T−1∑ k=0 η 2Γ ‖∇F (wk)‖2.
By simplifying the left-hand-side of the above inequality, we have T−1∑ k=0 [F (wk+1)− F (wk)] = F (wT )− F (w0) ≥ F̂ − F (w0),
where the inequality is due to F̂ ≤ F (wT ) (Assumption 4.12). Using the above, we have T−1∑ k=0 ‖∇F (wk)‖2 ≤ 2Γ[F (w0)− F̂ ] η . (47)

heading: A.9 PROOF OF THEOREM 4.17
text: Lemma A.4 (see Lemma 1 (Nguyen et al., 2018) or Lemma 2.4 (Gower et al., 2019)). Assume that ∀i the function Fi(w) is convex and L-smooth. Then, ∀w ∈ Rd the following hold:
EI [‖∇FI(w)‖2] ≤ 4L(F (w)− F (w∗)) + 2σ2. (48)
Theorem A.5. Suppose that Assumptions 4.2, 4.3, 4.4, 4.15 and 4.16 hold. Let parameters α > 0, η > 0, β2 ∈ (0, 1] are chosen such that α < 2ΓLµ , η ≤ α 2L , β2 > 1 − ηµα 2Γ(Γ−ηµ) > 0. Let {wk} be the iterates generated by Algorithm 1, then, for all k ≥ 0,
E[‖wk − w∗‖2D̂k ] ≤ (1− c) k‖r0‖2D̂0 + (1 + 2(1−β2)Γ α ) 2σ2η2 αc , (49)
where c = ηµα−(Γ−ηµ)2(1−β2)ΓΓα ∈ (0, 1). Remark A.6. If we choose β2 := 1− ηµα4Γ(Γ−ηµ) > 0 then
E[‖wk − w∗‖2D̂k ] ≤ ( 1− ηµ4Γ )k ‖w0 − w∗‖2D̂0 + 4 2Γ−ηµ(Γ−ηµ)α Γµσ2η. Proof. In order to bound ‖rk+1‖2D̂k+1 by ‖rk+1‖ 2 D̂k we will use that
0 ≺ D̂k+1 = D̂k + D̂k+1 − D̂k D̂k + ‖D̂k+1 − D̂k‖∞I D̂k + ‖D̂k+1 − D̂k‖∞I D̂k + ‖Dk+1 −Dk‖∞I D̂k + 2(1− β2)ΓI
D̂k + 2(1− β2)Γ 1αD̂k (1 + 2(1−β2)Γ α )D̂k. (50)
Then
‖rk+1‖2D̂k+1 (50) ≤ (1 + 2(1−β2)Γα )‖rk+1‖ 2 D̂k = (1 + 2(1−β2)Γα )‖rk − ηD̂k −1∇FIk(wk)‖2D̂k
= (1 + 2(1−β2)Γα )‖rk‖ 2 D̂k − 2η(1 + 2(1−β2)Γα )〈rk,∇FIk(wk)〉 + (1 + 2(1−β2)Γα )η 2 ( ‖∇FIk(wk)‖∗D̂k )2 ≤ (1 + 2(1−β2)Γα )‖rk‖ 2 D̂k − 2η(1 + 2(1−β2)Γα )〈rk,∇FIk(wk)〉
+ (1 + 2(1−β2)Γα ) η2
α ‖∇FIk(wk)‖2. (51)
Now, taking an expectation with respect to Ik conditioned on the past, we obtain
E[‖rk+1‖2D̂k+1 ] (51) ≤ (1 + 2(1−β2)Γα )‖rk‖ 2 D̂k + 2η(1 + 2(1−β2)Γα ) ( F ∗ − F (wk)− µ 2 ‖wk − w∗‖2 ) + (1 + 2(1−β2)Γα ) η2
α
( 4L(F (wk)− F ∗) + 2σ2 ) ≤ (1 + 2(1−β2)Γα )‖rk‖ 2 D̂k + 2η(1 + 2(1−β2)Γα ) ( F ∗ − F (wk)− µ 2Γ ‖rk‖2D̂k
) + (1 + 2(1−β2)Γα ) η2
α
( 4L(F (wk)− F ∗) + 2σ2 ) = (1 + 2(1−β2)Γα ) ( 1− 2η µ
2Γ
) ‖rk‖2D̂k
+ (1 + 2(1−β2)Γα ) ( 4Lη 2 α − 2η ) ((F (wk)− F ∗)) + (1 + 2(1−β2)Γα ) η2 α 2σ 2
≤ (1 + 2(1−β2)Γα ) ( 1− η µ Γ ) ‖rk‖2D̂k + (1 + 2(1−β2)Γ α ) η2 α 2σ 2, (52)
where we used the fact that 2L ηα − 1 ≤ 0→ η ≤ α 2L . In order to achieve a convergence we need to have (1 + 2(1−β2)Γα ) ( 1− µη
Γ
) =: 1− c < 1.
We have
(1 + 2(1−β2)Γα ) ( 1− ηµ Γ ) = 1− ηµα− (Γ− ηµ)2(1− β2)Γ
Γα︸ ︷︷ ︸ c
Now, we need to choose β2 such that
1 > β2 > 1− ηµα
2Γ(Γ− ηµ) > 0.
With such a choice we can conclude that
E[‖rk‖2D̂k ] (52) ≤ (1− c)‖rk−1‖2D̂k−1 + (1 + 2(1−β2)Γ α ) η2 α 2σ 2,
≤ (1− c)k‖r0‖2D̂0 + k−1∑ i=0 (1− c)i(1 + 2(1−β2)Γα ) η2 α 2σ 2,
≤ (1− c)k‖r0‖2D̂0 + ∞∑ i=0 (1− c)i(1 + 2(1−β2)Γα ) η2 α 2σ 2,
≤ (1− c)k‖r0‖2D̂0 + (1 + 2(1−β2)Γ α ) 2σ2η2 αc .
Theorem 4.17. Suppose that Assumptions 4.2, 4.3, 4.4, 4.15 and 4.16 hold. Let {wk} be the iterates generated by Algorithm 1 with ηk =η ∈ (0, α 2µ ΓL2 ), then, for all k ≥ 0,
E[F (wk)− F ∗] ≤ (1− c)k (F (w0)− F ∗) + η2Lσ2
cα2 , (53)
where c = 2ηµΓ − 2η2L2 α2 ∈ (0, 1). Moreover, if ηk =η ∈ (0, α2µ 2ΓL2 ) then
E[F (wk)− F ∗] ≤ (
1− ηµ Γ
)k (F (w0)− F ∗) + ηΓLσ2
α2µ . (54)
Proof. First, we will upper-bound the F (wk+1) using smoothness of F (w)
F (wk+1) = F (wk − ηD̂k−1∇FIk(wk))
≤ F (wk) +∇F (wk)T (−ηD̂k−1∇FIk(wk)) + L 2 ‖ηD̂k−1∇FIk(wk)‖2 ≤ F (wk)− η∇F (wk)T D̂k−1∇FIk(wk) + η2L
2α2 ‖∇FIk(wk)‖2,
where the first inequality is because of Assumptions 4.2 and 4.4, and the second inequality is due to Remark 4.10. By taking the expectation over the sample Ik, we have
EIk [F (wk+1)] ≤ F (wk)− ηEIk [∇F (wk)T D̂k−1∇FIk(wk)] + η2L
2α2 EIk [‖∇FIk(wk)‖2]
= F (wk)− η∇F (wk)T D̂k−1∇F (wk) + η2L
2α2 EIk [‖∇FIk(wk)‖2]
≤ F (wk)− η
Γ ‖∇F (wk)‖2 +
η2L 2α2 ( 4L(F (wk)− F ∗) + 2σ2 ) ≤ F (wk) + η
Γ 2µ (F ∗ − F (wk)) +
η2L 2α2 ( 4L(F (wk)− F ∗) + 2σ2 ) (55)
Subtracting F ∗ from both sides, we obtain
EIk [F (wk+1)− F ∗] (55) ≤ (
1− η Γ 2µ+ η2L 2α2 4L
) (F (wk)− F ∗) + η2Lσ2
α2
≤ 1− (2ηµΓ − 2η2L2α2 )
︸ ︷︷ ︸ c
 (F (wk)− F ∗) + η2Lσ2α2 . (56)
By taking the total expectation over all batches I0, I1, I2,... and all history starting with w0, we have
E[F (wk)− F ∗] (56) ≤ (1− c)k (F (w0)− F ∗) + k−1∑ i=0 (1− c)i η 2Lσ2 α2
≤ (1− c)k (F (w0)− F ∗) + ∞∑ i=0 (1− c)i η 2Lσ2 α2
= (1− c)k (F (w0)− F ∗) + η2Lσ2
cα2
and the (53) is obtained. Now, if η ≤ α 2µ
2ΓL2 then
c = 2ηµ Γ − 2η
2L2
α2 ≥ 2η
( µ
Γ − α
2µ 2ΓL2 L2 α2
) = 2η (µ Γ − µ 2Γ ) = ηµ Γ
and (54) follows.

heading: A.10 PROOF OF THEOREM 4.18
text: Theorem 4.18. Suppose that Assumptions 4.3, 4.12, 4.2, 4.14 and 4.16 hold. Let {wk} be the iterates generated by Algorithm 1, where 0 < ηk = η ≤ η2
LΓ , and w0 is the starting point. Then, for all
k ≥ 0,
E
[ 1
T T−1∑ k=0
‖∇F (wk)‖2 ] ≤ 2Γ[F (w0)− F̂ ]
ηT + ηΓγ2L α2 T→∞−−−−→ ηΓγ 2L α2 .
Proof. We have that F (wk+1) = F (wk − ηD̂k−1∇FIk(wk))
≤ F (wk) +∇F (wk)T (−ηD̂k−1∇FIk(wk)) + L 2 ‖ηD̂k−1∇FIk(wk)‖2 ≤ F (wk)− η∇F (wk)T D̂k−1∇FIk(wk) + η2L
2α2 ‖∇FIk(wk)‖2,
where the first inequality is because of Assumptions 4.2 and 4.4, and the second inequality is due to Remark 4.10. By taking the expectation over the sample Ik, we have
EIk [F (wk+1)] ≤ F (wk)− ηEIk [∇F (wk)T D̂k−1∇FIk(wk)] + η2L
2α2 EIk [‖∇FIk(wk)‖2]
= F (wk)− η∇F (wk)T D̂k−1∇F (wk) + η2L
2α2 EIk [‖∇FIk(wk)‖2] ≤ F (wk)− η ( 1
Γ − ηL 2α2
) ‖∇F (wk)‖2 + η2γ2L
2α2 (57)
≤ F (wk)− η
2Γ ‖∇F (wk)‖2 +
η2γ2L
2α2 , (58)
where the second inequality is due to Remark 4.10 and Assumption 4.14, and the third inequality is due to the choice of the step length.
By inequality (58) and taking the total expectation over all batches I0, I1, I2,... and all history starting with w0, we have
E[F (wk+1)− F (wk)] ≤ − η
2Γ E[‖∇F (wk)‖2] +
η2γ2L
2α2 .
By summing both sides of the above inequality from k = 0 to T − 1, T−1∑ k=0 E[F (wk+1)− F (wk)] ≤ − η 2Γ T−1∑ k=0 E[‖∇F (wk)‖2] + η2γ2LT 2α2
= − η 2Γ E [ T−1∑ k=0 ‖∇F (wk)‖2 ] + η2γ2LT 2α2 .
By simplifying the left-hand-side of the above inequality, we have T−1∑ k=0 E [F (wk+1)− F (wk)] = E[F (wT )]− F (w0) ≥ F̂ − F (w0),
where the inequality is due to F̂ ≤ F (wT ) (Assumption 4.12). Using the above, we have
E [ T−1∑ k=0 ‖∇F (wk)‖2 ] ≤ 2Γ[F (w0)− F̂ ] η + ηΓγ2LT α2 .

heading: A.11 PROOFS WITH LINE SEARCH
text: Given the current iterate wk, the steplength is chosen to satisfy the following sufficient decrease condition
F (wk + ηkpk) ≤ F (wk)− c1ηk∇F (wk)T D̂k−1∇F (wk), (59)
where pk = −D̂k−1∇F (wk) and c1 ∈ (0, 1). The mechanism works as follows. Given an initial steplength (say ηk = 1), the function is evaluated at the trial point wk + αkpk and condition (59) is checked. If the trial point satisfies (59), then the step is accepted. If the trial point does not satisfy (59), the steplength is reduced (e.g., ηk = τηk for τ ∈ (0, 1)). This process is repeated until a steplength that satisfies (59) is found.
Algorithm 2 Backtracking Armijo Linesearch (Nocedal & Wright, 2006) Input: wk, pk 1: Select ηinitial, c1 ∈ (0, 1) and τ ∈ (0, 1) 2: η0 = ηinitial, j = 0 3: while F (wk + ηkpk) > F (wk) + c1ηk∇F (wk)T pk do 4: Set ηj+1 = τηj
5: Set j = j + 1 6: end while
Output: ηk = ηj
By following from the study (Berahas et al., 2019), we have the following theorems.

heading: A.11.1 DETERMINITIC REGIME - STRONGLY CONVEX
text: Theorem A.7. Suppose that Assumptions 4.3 and 4.2 and 4.4 hold. Let {wk} be the iterates generated by Algorithm 1, where ηk is the maximum value in {τ−j : j = 0, 1, . . . } satisfying (59) with 0 < c1 < 1, and w0 is the starting point. Then for all k ≥ 0,
F (wk)− F ? ≤ ( 1− 4µα 2c1(1− c1)τ
Γ2L
)k [F (w0)− F ?] .
Proof. Starting with (43) we have F (wk − ηkD̂k−1∇F (wk)) ≤ F (wk)− ηk ( 1
Γ − ηk
L
2α2
) ‖∇F (wk)‖2. (60)
From the Armijo backtracking condition (59), we have
F (wk − ηkD̂k−1∇F (wk)) ≤ F (wk)− c1ηk∇F (wk)T D̂k−1∇F (wk)
≤ F (wk)− c1ηk
Γ ‖∇F (wk)‖2. (61)
Looking at (60) and (61), it is clear that the Armijo condition is satisfied when
ηk ≤ 2α2(1− c1)
ΓL . (62)
Thus, any ηk that satisfies (62) is guaranteed to satisfy the Armijo condition (59). Since we find ηk using a constant backtracking factor of τ < 1, we have that
ηk ≥ 2α2(1− c1)τ
ΓL . (63)
Therefore, from (60) and by (62) and (63) we have F (wk+1) ≤ F (wk)− ηk ( 1
Γ − ηkL 2α2
) ‖∇F (wk)‖2
≤ F (wk)− ηkc1
Γ ‖∇F (wk)‖2
≤ F (wk)− 2α2c1(1− c1)τ
Γ2L ‖∇F (wk)‖2. (64)
By strong convexity, we have 2µ(F (w)− F ?) ≤ ‖∇F (w)‖2, and thus
F (wk+1) ≤ F (wk)− 4µα2c1(1− c1)τ
Γ2L (F (w)− F ?). (65)
Subtracting F ? from both sides, and applying (65) recursively yields the desired result.

heading: A.11.2 DETERMINISTIC REGIME - NONCONVEX
text: Theorem A.8. Suppose that Assumptions 4.3, 4.12 and 4.2 hold. Let {wk} be the iterates generated by Algorithm 1, where ηk is the maximum value in {τ−j : j = 0, 1, . . . } satisfying (59) with 0 < c1 < 1, and where w0 is the starting point. Then,
lim k→∞
‖∇F (wk)‖ = 0, (66)
and, moreover, for any T > 1,
1
T T−1∑ k=0 ‖∇F (wk)‖2 ≤ Γ2L[F (w0)− F̂ ] 2α2c1(1− c1)τT τ→∞−−−−→ 0.
Proof. We start with (64) F (wk+1) ≤ F (wk)−
2α2c1(1− c1)τ Γ2L ‖∇F (wk)‖2.
Summing both sides of the above inequality from k = 0 to T − 1, T−1∑ k=0 (F (wk+1)− F (wk)) ≤ − T−1∑ k=0 2α2c1(1− c1)τ Γ2L ‖∇F (wk)‖2.
The left-hand-side of the above inequality is a telescopic sum and thus, T−1∑ k=0 [F (wk+1)− F (wk)] = F (wT )− F (w0) ≥ F̂ − F (w0),
where the inequality is due to F̂ ≤ F (wT ) (Assumption 4.12). Using the above, we have T−1∑ k=0 ‖∇F (wk)‖2 ≤ Γ2L[F (w0)− F̂ ] 2α2c1(1− c1)τ . (67)
Taking limits we obtain,
lim τ→∞ τ−1∑ k=0 ‖∇F (wk)‖2 <∞,
which implies (66). Dividing (67) by T we conclude
1
T T−1∑ k=0 ‖∇F (wk)‖2 ≤ Γ2L[F (w0)− F̂ ] 2α2c1(1− c1)τT .

heading: A.12 PROOFS WITH DECAYING STEP-SIZE
text: In this section, we provide the analysis with decaying stepsize for both strongly convex and nonconvex cases in the stochastic regime. The proofs for the following theorems follows from the proofs in Theorems 4.17 and 4.18, and Theorems 4.7 and 4.9 in Bottou et al. (2018).

heading: A.12.1 STOCHASTIC REGIME – STRONGLY CONVEX
text: Theorem A.9. Suppose that Assumptions 4.2, 4.3, 4.4, 4.15 and 4.16 hold., and let F ? = F (w?), where w? is the minimizer of F . Let {wk} be the iterates generated by Algorithm 1, where ηk is a sequence of stepsize such that, for all k ≥ 0,
ηk = φ
ζ + k for some φ >
Γ µ and ζ > 0 such that η0 ≤
α2µ
2ΓL2 (68)
Then, for all k ≥ 0, the expected optimality gap satisfies E[F (wk)− F ?] ≤ ν
ζ + k (69)
where
ν := max { φ2Lσ2
( φµ
Γ − 1)α2
, ζ(F (w0)− F ?) }
(70)
Proof. By taking total expectation from (56), and replacing η with ηk we have:
E[F (wk+1)− F ∗] ≤ 1− (2ηkµΓ − 2η2kL2α2 )
︸ ︷︷ ︸ c E[F (wk)− F ∗] + η2kLσ2α2 . (71) By the assumption η0 ≤ α 2µ 2ΓL2 , and the way that ηk is designed, i.e., in a decaying way, we have ηk < η0 ≤ α 2µ
2ΓL2 , therefore
c = 2ηkµ Γ − 2η
2 kL 2
α2 ≥ 2ηk
( µ
Γ − α
2µ 2ΓL2 L2 α2
) = 2ηk (µ Γ − µ 2Γ ) = ηkµ Γ .
Thus, we have:
E[F (wk+1)− F ∗] ≤ (
1− ηkµ Γ
) E[F (wk)− F ∗] + η2kLσ 2
α2 . (72)
We do the rest of the analysis of (69) with induction. First, for k = 0, (69) holds by the definition of ν. Next, assuming (69) holds for k ≥ 0, by (72) we have:
E[F (wk+1)− F ?] ≤ 1− φµΓ k̂  ν k̂ + φ2Lσ2 k̂2α2
=  k̂ − φµΓ k̂2  ν + φ2Lσ2 k̂2α2
= ( k̂ − 1 k̂2 ) ν−  φµΓ − 1 k̂2  ν + φ2Lσ2 k̂2α2︸ ︷︷ ︸
nonpositive by definition of ν
≤ ν k̂ + 1
(73)

heading: A.12.2 STOCHASTIC REGIME – NONCONVEX
text: Theorem A.10. Suppose that Assumptions 4.3, 4.12, 4.2, 4.14 and 4.16 hold. Let {wk} be the iterates generated by Algorithm 1, where ηk is a steplength sequence satisfying
∞∑ k=0 ηk =∞ and ∞∑ k=0 η2k <∞. (74)
and w0 is the starting point. Then, with AT := ∑T−1 k=0 ηk,
lim T→∞
E[ ∑T−1 k=0 ηk‖∇F (wk)‖2] <∞, (75)
and therefore E[ 1
AT
∑T−1 k=0 ηk‖∇F (wk)‖2] T→∞−−−−→ 0. (76)
Proof. The steplength sequence {ηk} goes to zero due to the second condition in (74). meaning that w.l.o.g., we may assume that ηk ≤ 2α2
2ΓL for all k ≥ 0. By starting from and taking total expectation
from (57), we have
E[F (wk+1)]− E[F (wk)] ≤ −ηk ( 1
Γ − ηkL 2α2
) E[‖∇F (wk)‖2] + η2kγ 2L
2α2
≤ − ηk 2Γ E[‖∇F (wk)‖2] + η2kγ 2L 2α2 .
Summing both sides of this inequality for k ∈ {0, . . . , T − 1} gives:
F̂ − E[F (w0)] ≤ E[F (wT )]− E[F (w0)] ≤ − 1
2Γ
∑T−1 k=0 ηkE[‖∇F (wk)‖2] + γ2L 2α2 ∑T−1 k=0 η 2 k.
(77)
Multiplying by 2Γ, and rearranging we get∑T−1 k=0 ηkE[‖∇F (wk)‖2] ≤ 2Γ(E[F (w0)]− F̂ ) + Γγ2L α2 ∑T−1 k=0 η 2 k. (78)
The second condition in (74) guarantees that the right-hand side of the above inequality converges to a finite limit when T increases, therefore, it implies (75). By the first condition in (74), we conclude that AT →∞ as T →∞, which results in (76).

heading: A.13 CONVERGENCE BOUNDS INTERPRETATION
text: In this section, we provide the interpretation behind the convergence bounds in our theoretical results.
• In Theorem 4.6 there is the term (1 − β2)Γ Qk k
, which essentially highlights that the parameter β2 controls the size of the neighborhood of optimal solution(s). In other words, if β2 is very close to 1 (which is the case for our algorithm), the convergence neighborhood would be very small. Meaning that, with the assumption regarding bounded iterates, our algorithm converges to a very small neighborhood of optimal solution as β2 is very close to 1 (while in SGD the step-size controls the size of this neighborhood and accordingly, the step-size needs to go zero in order to converge to the small neighborhood of the optimal solution(s)). Note this convergence guarantee is for smooth and convex case with adaptive learning rate.
• Theorem 4.9 provides a linear convergence rate for the smooth and strongly convex case with adaptive learning rate. The analysis behind this theorem is based on the Lyapunov energy that decreases linearly.
• Remark 4.10 guarantees that the OASIS preconditioning matrix is positive definite and its eigenvalues are uniformly bounded above and below.
• Theorems 4.11 and 4.13 provides linear and sublinear rate to the stationary point(s) for deterministinc strongly convex and nonconvex cases, respectively. The learning rate is fixed in both mentioned cases, and the convergence rates are similar to those of limited-memory quasi-Newton approaches (e.g., L-BFGS).
• Theorems 4.17 and 4.18 provides linear and sublinear rate to the neighborhood of stationary point(s) (in expectation) for deterministinc strongly convex and nonconvex cases , respectively. Please note that we provided another convergence guarantee for the stochastic strongly convex case (please see Theorem A.5). Again, as the deterministic cases, the provided convergence guarantees completely match with those of the classical quasi-Newton methods (such as L-BFGS); the learning rate is fixed in the latter cases. Please note that these methods (such as OASIS and classical limited-memory quasi-Newton methods) in theory are not better than GD-type methods. In practice, however, they have shown their strength.
• Theorems A.9 and A.10 provides convergence guarantees for the stochastic setting with decaying learning rate. The mentioned theorems provides the convergence guarantees to the stationary point(s) (in expectation) for the stochastic strongly convex and nonconvex cases.
• Theorems A.7 and A.8 provides linear and sublinear convergence rates to the stationary point(s) for the deterministic strongly convex and nonconvex cases with utilizing linesearch.

heading: B ADDITIONAL ALGORITHM DETAILS

heading: B.1 RELATED WORK
text: As was mentioned in Section 2, we follow the generic iterate update:
wk+1 = wk − ηkD̂k−1mk.
Table 3 summarizes the methodologies discussed in Section 2.
B.2 STOCHASTIC OASIS
Here we describe a stochastic variant of OASIS in more detail. As mentioned in Section 3, to estimate gradient and Hessian diagonal, the choices of sets Ik, Jk ⊂ [n], are independent and correspond to only a fraction of data. This change results in Algorithm 3.
Algorithm 3 Stochastic OASIS Input: w0, η0, Ik, D0, θ0 = +∞, β2, α
1: w1 = w0 − η0D̂0−1∇F (w0)Ik 2: for k = 1, 2, . . . do 3: Calculate Dk = β2Dk−1 + (1− β2) diag(zk ∇2FJk(wk)zk) 4: Calculate D̂k by setting (D̂k)i,i = max{|Dk|i,i, α}, ∀i ∈ [d] 5: Update ηk = min{ √ 1 + θk−1ηk−1,
‖wk−wk−1‖D̂k 2‖∇FIk (wk)−∇FIk (wk−1)‖
∗ D̂k
}
6: Set wk+1 = wk − ηkD̂k−1∇FIk(wk) 7: Set θk =
ηk ηk−1
8: end for
Additionally, in order to compare the performance of our preconditioner schema independent of the adaptive learning-rate rule, we also consider variants of OASIS with fixed η. We explore two modifications, denoted as “Fixed LR” and “Momentum” in Section 5. “Fixed LR” is obtained from Algorithm 3 by simply having a fixed scalar η for all iterations, which results in Algorithm 4. “Momentum” is obtained from “Fixed LR” by considering a simple form of first-order momentum with a parameter β1, and this results in Algorithm 5. In Section C.5, we show that OASIS is robust with respect to the different choices of learning rate, obtaining a narrow spectrum of changes.
Algorithm 4 OASIS- Fixed LR Input: w0, η, Ik, D0, β2, α
1: w1 = w0 − ηD̂0−1∇FIk(w0) 2: for k = 1, 2, . . . do 3: Calculate Dk = β2Dk−1 + (1 − β2) diag(zk ∇2FJk(wk)zk) 4: Calculate D̂k by setting (D̂k)i,i = max{|Dk|i,i, α}, ∀i ∈
[d]
5: Set wk+1 = wk − ηD̂k−1∇FIk(wk) 6: end for
Algorithm 5 OASIS- Momentum Input: w0, η, Ik, D0, β1, β2 α
1: Set m0 = ∇FIk(w0) 2: w1 = w0 − ηD̂0−1m0 3: for k = 1, 2, . . . do 4: Calculate Dk = β2Dk−1 + (1 − β2) diag(zk ∇2FJk(wk)zk) 5: Calculate D̂k by setting (D̂k)i,i = max{|Dk|i,i, α}, ∀i ∈ [d] 6: Calculate mk = β1mk−1 + (1− β1)∇FIk(wk) 7: Set wk+1 = wk − ηD̂k−1mk 8: end for
In our experiments, we used a biased version of the algorithm with Ik = Jk. One of the benefits of this choice is to compute the Hessian-vector product efficiently. By reusing the computed gradient, the overhead of computing gradients with respect to different samples is significantly reduced (see Section B.3).
The final remark is related to a strategy to obtain D0, which is required at the start of OASIS Algorithm. One option is to do warmstarting, i.e., spend some time before training in order to sample some number of Hutchinson’s estimates. The second option is to use bias corrected rule for Dk
similar to the one used in Adam and AdaHessian
Dk = β2Dk−1 + (1− β2)diag(zk ∇2FJk(wk)zk),
Dcork = Dk
1− βk+12 ,
which allows to obtain D0 by defining D−1 to be a zero diagonal.

heading: B.3 EFFICIENT HESSIAN-VECTOR COMPUTATION
text: In the OASIS Algorithm, similar to AdaHessian (Yao et al., 2020) methodology, the calculation of the Hessian-vector product in Hutchinson’s method is the main overhead. In this section, we present how this product can be calculated efficiently. First, lets focus on two popular machine learning problems: (i) logistic regression; and (ii) non-linear least squares problems. Then, we show the efficient calculation of this product in deep learning problems.
Logistic Regression. One can note that we can show the `2-regularized logistic regression as:
F (w) = 1
n
( 1n ∗ log ( 1 + e−Y X Tw )) + λ
2 ‖w‖2, (79)
where 1n is the vector of ones with size 1× n, ∗ is the standard multiplication operator between two matrices, and X is the feature matrix and Y is the label matrix. Further, the Hessian-vector product for logistic regression problems can be calculated as follows (for any vector v ∈ Rd):
∇2F (w) ∗ v = 1 n  X T ∗ Y Y e−Y XT w( 1 + e−Y XT w )2  X ∗ v︸ ︷︷ ︸
1©︸ ︷︷ ︸ 2©︸ ︷︷ ︸
3©
 + λv. (80)
The above calculation shows the order of computations in order to compute the Hessian-vector product without constructing the Hessian explicitly.
Non-linear Least Squares. The non-linear least squares problems can be written as following:
F (w) = 1
2n ‖Y − φ(XTw)‖2. (81)
The Hessian-vector product for the above problem can be efficiently calculated as:
∇2F (w) ∗ v =
1
n  −XT ∗ [ φ(X T w) (1− φ(XTw)) (Y − 2(1 + Y ) φ(XTw) + 3φ(XTw) φ(XTw)) ] X ∗ v︸ ︷︷ ︸ 1©︸ ︷︷ ︸ 2©︸ ︷︷ ︸
3©
 . (82)
Deep Learning. In general, the Hessian-vector product can be efficiently calculated by:
∇2F (w) ∗ v = ∂ 2F (w) ∂w∂w ∗ v = ∂ ∂w ( ∂F (w)T ∂w v). (83)
As is clear from (83), two rounds of back-propagation are needed. In fact, the round regarding the gradient evaluation is already calculated in the corresponding iteration; and thus, one extra round of back-propagation is needed in order to evaluate the Hessian-vector product. This means that the extra cost for the above calculation is almost equivalent to one gradient evaluation.

heading: C DETAILS OF EXPERIMENTS

heading: C.1 TABLE OF ALGORITHMS
text: Table 4 summarizes the algorithms implemented in Section 5.
To display the optimality gap for logistic regression problems, we used Trust Region (TR) Newton Conjugate Gradient (Newton-CG) method (Nocedal & Wright, 2006) to find a w such that ‖∇F (w)‖2 < 10−19. Hereafter, we denote F (w)− F (w∗) the optimality gap, where we refer w∗ to the solution found by TR Newton-CG.

heading: C.2 PROBLEM DETAILS
text: Some metrics for the image classification problems are given in Table 5.
The number of parameters for ResNet architectures is particularly important, showcasing that OASIS is able to operate in very high-dimensional problems, contrary to the widespread belief that methods using second-order information are fundamentally limited by the dimensionality of the parameter space.

heading: C.3 BINARY CLASSIFICATION
text: In Section 5 of the main paper, we studied the empirical performance of OASIS on binary classification problems in a deterministic setting, and compared it with AdGD and AdaHessian. Here, we present the extended details on the experiment.

heading: C.3.1 PROBLEM AND DATA
text: As a common practice for the empirical research of the optimization algorithms, LIBSVM datasets6 are chosen for the exercise. Specifically, we chose 5 popular binary class datasets, ijcnn1, rcv1, news20, covtype and real-sim. Table 6 summarizes the basic statistics of the datasets.
6Datasets are available at https://www.csie.ntu.edu.tw/˜cjlin/libsvmtools/ datasets/.
Let (xi, yi) be a training sample indexed by i ∈ [n] := {1, 2, ..., n}, where xi ∈ Rd is a feature vector and yi ∈ {−1,+1} is a label. The loss functions are defined in the forms
fi(w) = log(1 + e −yixTi w) +
λ 2 ‖w‖2, (84)
fi(w) = ( yi −
1
1 + e−x T i w
)2 , (85)
where (84) is a regularized logistic regression of a particular choice of λ > 0 (and we used λ ∈ {0.1/n, 1/n, 10/n} in the experiment), and hence a strongly convex function; and (85) is a non-linear least square loss, which is apparently non-convex.
The problem we aimed to solve is then defined in the form
min w∈Rd
{ F (w) := 1
n n∑ i=1 fi(w)
} , (86)
and we denote w∗ the global optimizer of (86) for logistic regression.

heading: C.3.2 CONFIGURATION OF ALGORITHM
text: To better evaluate the performance of the algorithms, we configured OASIS, AdGD and AdaHessian with different choices of parameters.
• OASIS was configured with different values of β2 and α, where β2 can be any values in {0.95, 0.99, 0.995, 0.999} and α can be any value in the set {10−3, 10−5, 10−7} (see Section C.5 where OASIS has a narrow spectrum of changes with respect to different values of α and β2). In addition, we adopted a warmstarting approach to evaluate the diagonal of Hessian at the starting point w0. • AdGD was configured with 12 different values of the initial learning rate, i.e., η0 ∈ {10−11, 10−10, ..., 0.1, 1.0}. • AdaHessian was configured with different values of the fixed learning rate: for logistic regression, we used 12 values between 0.1 and 5; for non-linear least square, we used 0.01, 0.05, 0.1, 0.5, 1.0 and 2.0.
To take into account the randomness of the performance, we used 10 distinct random seeds to initialize w0 for each algorithm, dataset and problem.

heading: C.3.3 EXTENDED EXPERIMENTAL RESULTS
text: This section presents the extended results on our numerical experiments. For the whole experiments in this paper, we ran each method 10 times starting from different initial points.
For logistic regression, we show the evolution of the optimality gap in Figure 6 and the ending gap in Figure 7; the evolution of the testing accuracy and the maximum achieved ones are shown in Figure 8 and Figure 9 respectively.
For non-linear least square, the evolution of the objective and its ending values are shown in Figure 10; the evolution of the testing accuracy along with the maximum achived ones are shown in Figure 11.
C.4 IMAGE CLASSIFICATION
In the following sections, we provide the results on standard bench-marking neural network training tasks: CIFAR10, CIFAR100, and MNIST.
C.4.1 CIFAR10
In our experiments, we compared the performance of the algorithm described in Table 4 in 2 settings - with and without weight decay. The procedure for choosing parameter values differs slightly for these two cases, so they are presented separately. Implementations of ResNet20/32 architectures are taken from AdaHessian repository.7
Analogously to (Malitsky & Mishchenko, 2020), we also modify terms in the update formula for ηk
ηk = min{ √ 1 + θk−1ηk−1, ‖wk−wk−1‖D̂k
2‖∇F (wk)−∇F (wk−1)‖∗ D̂k
}.
Specifically, we incorporate parameter γ into √
1 + γθk−1 and use more optimistic bound of 1/Lk instead of 1/2Lk, which results in a slightly modified rule
ηk = min{ √ 1 + γθk−1ηk−1, ‖wk−wk−1‖D̂k
‖∇F (wk)−∇F (wk−1)‖∗ D̂k
}.
In order to find the best value for γ, we include it into hyperparameter tuning procedure, ranging the values in the set {1, 0.1, 0.05, 0.02, 0.01} (see Section C.5 where OASIS has a narrow spectrum of changes with respect to different values of γ).
Moreover, we used a learning-rate decaying strategy as considered in (Yao et al., 2020) to have the same and consistent settings. In the aforementioned setting, ηk would be decreased by a multiplier on some specific epochs common for every algorithms (the epochs 80 and 120).
7https://github.com/amirgholami/adahessian
With weight decay: It is important to mention that due to variance in stochastic setting, the learning rate needs to decay (regardless of the learning rate rule is adaptive or fixed) to get convergence. One can apply adaptive learning rate rule without decaying by using variance-reduced methods (future work). Experimental setup is generally taken to be similar to the one used in (Yao et al., 2020). Parameter values for SGD, Adam, AdamW and AdaHessian are taken from the same source, meaning that learning rate is set to be 0.1/0.001/0.01/0.15 and, where they are used, β1 and β2 are set to be 0.9 and 0.999. For OASIS, due to a significantly different form of the preconditioner aggregation, we conducted a small scale search for the best value of β2, choosing from the set {0.999, 0.99, 0.98, 0.97, 0.96, 0.95}, which produced values of 0.99/0.95/0.98 for “Fixed LR,” “Momentum” and adaptive variants correspondingly for ResNet20 and 0.99/0.99/0.98 for “Fixed LR,” “Momentum” and adaptive variants correspondingly for ResNet32, while β1 for momentum variant is still set to 0.9. Additionally, in the experiments with fixed learning rate, it is tuned for each architecture, producing values 0.1/0.1 for ResNet20 and 0.1/0.1 for ResNet32 for “Fixed LR” and “Momentum” variants of OASIS correspondingly. For the adaptive variant, γ is set to 0.1/1.0 for ResNet20/32. α is set to 0.1 for all OASIS variants. We train all methods for 160 epochs and for all methods. With fixed learning rate we employ identical scheduling, reducing step size by a factor of 10 at epochs 80 and 120. For the adaptive variant scheduler analogue is implemented, multiplying step size by ρ at epochs 80 and 120, where ρ is set to be 0.1/0.5 for ResNet20/32. Weight decay value for all optimizers with fixed learning rate is 0.0005 and decoupling for OASIS is done similarly to AdamW and AdaHessian. For adaptive OASIS weight decay is set to 0.001 without decoupling. Batch size for all optimizers is 256.
Without weight decay: In this setting we tune learning rate for SGD, Adam, AdamW and AdaHessian, obtaining the values 0.15/0.005/0.005/0.25 and 0.125/0.01/0.01/0.25 for ResNet20/32 accordingly. Where relevant, β1, β2 are taken to be 0.9, 0.999. For OASIS we still try to tune β2, which produces values of 0.999 for ResNet20 and adaptive case of ResNet32 and 0.99 for “Fixed LR,” “Momentum” in the case of ResNet32; for the momentum variant β1 = 0.9. Learning rates are 0.025/0.05 for ResNet20 and 0.025/0.1 for ResNet32 for “Fixed LR” and “Momentum” variants of OASIS correspondingly. For adaptive variant γ is set to 0.01/0.01 for ResNet20/32. α is set to 0.1 for adaptive OASIS variant, while hyperoptimization showed, that for “Fixed LR” and “Momentum” a value of 0.01 can be used. We train all methods for 160 epochs and for all methods with fixed learning rate we employ identical scheduling, reducing step size by a factor of 10 at epochs 80 and 120. For the adaptive variant scheduler analogue is implemented, multiplying step size by ρ at epochs 80 and 120, where ρ is set to be 0.1/0.1 for ResNet20/32. Batch size for all optimizers is 256.
Finally, all of the results presented are based on 10 runs with different seeds, where parameters are chosen amongst the best runs produced at the preceding tuning phase.
All available results can be seen in Table 7 and Figures 12 and 13. We ran our experiments on an NVIDIA V100 GPU.
C.4.2 CIFAR100
Analogously to CIFAR10, we do experiments in settings with and without weight decay. In both cases we took best learning rates from the corresponding CIFAR10 experiment with ResNet20 architecture without any additional tuning, which results in 0.1 for all. β2 is set to 0.99 in the case with weight decay and to 0.999 in the case without it. ResNet18 architecture implementation is taken from a github repository.8 We train all methods for 200 epochs. For all optimizers learning rate is decreased (or effectively decreased in case of fully adaptive OASIS) by a factor of 5 at epochs 60, 120 and 160.
All of the results presented are based off of 10 runs with different seeds. All available results can be seen in Table 8 and Figure 14. We ran our experiments on an NVIDIA V100 GPU.

heading: C.4.3 MNIST
text: Similar to the previous results, we ran the methods with 10 different random seeds. We considered Net DNN which has 2 convolutional layers and 2 fully connected layers with ReLu non-linearity, mentioned in the Table 5. In order to tune the hyperparameters for other algorithms, we considered the set of learning rates {100, 10−1, 10−2, 10−3}. For OASIS, we used the set of {10−1, 10−2} and the
8https://github.com/uoguelph-mlrg/Cutout
set for truncation parameter α ∈ {10−1, 10−2}. As is clear from the results shown in Figure 15 and Table 9, OASIS with momentum has the best performance in terms of training (lowest loss function), and and it is comparable with the best test accuracy for the cases with and without weight decay. It is worth mentioning that OASIS with adaptive learning rate got satisfactory results with lower number of parameters, which require tuning, which is vividly important, especially in comparison to first-order methods that are sensitive to the choice of learning rate. We ran our experiments on a Tesla K80 GPU.
C.5 SENSITIVITY ANALYSIS OF OASIS
It is worth noting that we designed and analyzed OASIS for the deterministic setting with adaptive learning rate. It is safe to state that no sensitive tuning is required in the deterministic setting (the learning rate is updated adaptively, and the performance of OASIS is completely robust even if the rest of the hyperparameters are not perfectly hand-tuned); see Figures 16 and 17.
We also provide sensitivity analysis for the image classification tasks. These problems are in stochastic setting, and we need to tune some of the hyperparameters in OASIS. As is clear from the following figures, OASIS is robust with respect to different settings of hyperparameters, and the spectrum of changes is narrow enough; implying that even if the hyperparameters aren’t properly tuned, we still get acceptable results that are comparable to other state-of-the-art first- and second-order approaches with less tuning efforts. Figure 18 shows that OASIS is robust for different values of learning rate (left figure), while SGD is completely sensitive with respect to learning rate choices. One possible reason for this is because OASIS uses well-scaled preconditioning, which scales each gradient component with regard to the local curvature at that dimension, whereas SGD treats all components equally.
Furthermore, in order to use an adaptive learning rate in the stochastic setting, an extra hyperparameter, γ, is used in OASIS (similar to (Malitsky & Mishchenko, 2020)). Figure 19 shows that OASIS is also robust with respect to different values of γ (unlike the study in (Malitsky & Mishchenko, 2020)). The aforementioned figures are for CIFAR10 dataset on the ResNet20 architecture; the same behaviour is observed for the other network architectures.
Finally, we show here that the performance of OASIS is also robust with respect to different values of β2. Figure 20 shows the robustness of OASIS across a range of β2 values. This figure is for CIFAR100 dataset on the ResNet18 architecture.

Review summaries:
The paper designs and analyses an algorithm for minimizing a separable function. It provides deterministic as well as stochastic versions, which either fully compute the gradient or sample it. The algorithm estimates the diagonal elements of the Hessian via Hessian-vector products. The algorithm makes use of this information for finding a better search direction and the step length, eliminating the need for a line search. It provides convergence guarantees and a number of experiments on classical ML problems as well as on deep nets.

The paper presents a novel adaptive second order method: ''OASIS'', for large scale optimization problems (convex and nonconvex). The search direction is obtained by preconditioning the gradient information with a matrix obtained by approximating the Hessian diagonal matrix (via Hutchinson's method with a momentum term). The learning rate is updated adaptively by approximating the Lipschitz smoothness parameter. On the theoretical front, convergence analysis is provided for the adaptive learning rate case, for the convex and strongly convex setups. Similar analysis is also provided for the fixed learning rate (LR) case, for the strongly convex and nonconvex settings. Finally, extensive empirical results are provided which show that the proposed method achieves comparable, and sometimes better results than other existing methods.

This work proposes OASIS, a second-order method which approximates the diagonal of Hessian matrix and uses the information to rescale the gradient vector. The main difference between OASIS and the existing method AdaHessian (Yao et al., 2020) is on the ways they approximate the diagonal of Hessian, that AdaHessian uses a formula similar to Adam and OASIS uses an exponential average. Moreover, OASIS also incorporates the adaptive stepsize in (Mishchenko & Malitsky, 2020). The authors established the convergence guarantees of OASIS under various settings including convex, strongly convex and nonconvex cases, using various learning rate schedulers such as the adaptive learning rate, fixed learning rate and line search. Empirical results on various machine learning tasks are provided to evaluate OASIS.