Nys-Newton: Nystr\"om-Approximated Curvature for Stochastic Optimization

Second-order optimization methods are among the most widely used optimization approaches for convex optimization problems, and have recently been used to optimize non-convex optimization problems such as deep learning models. The widely used second-order optimization methods such as quasi-Newton methods generally provide curvature information by approximating the Hessian using the secant equation. However, the secant equation becomes insipid in approximating the Newton step owing to its use of the first-order derivatives. In this study, we propose an approximate Newton sketch-based stochastic optimization algorithm for large-scale empirical risk minimization. Specifically, we compute a partial column Hessian of size ($d\times m$) with $m\ll d$ randomly selected variables, then use the \emph{Nystr\"om method} to better approximate the full Hessian matrix. To further reduce the computational complexity per iteration, we directly compute the update step ($\Delta\boldsymbol{w}$) without computing and storing the full Hessian or its inverse. We then integrate our approximated Hessian with stochastic gradient descent and stochastic variance-reduced gradient methods. The results of numerical experiments on both convex and non-convex functions show that the proposed approach was able to obtain a better approximation of Newton\textquotesingle s method, exhibiting performance competitive with that of state-of-the-art first-order and stochastic quasi-Newton methods. Furthermore, we provide a theoretical convergence analysis for convex functions.


Introduction
The problem of the optimization of various function is among the most critical and popular topics in machine learning and mathematical optimization. Let {(x i , y i )} n i=1 be given n training samples, where x i ∈ R d and y i ∈ {−1, 1}, and f be an objective function defined as follows.
where f i (w) = f (w; x i , y i ) : R d → R is the loss function. Well-known convex loss functions include logistic loss f i (w) = log(1 + exp (−y i x i w)).
To optimize (1), first-order optimization methods such as stochastic gradient descent (SGD) [36], AdaGrad [13], stochastic variance-reduced gradient (SVRG) [21], SAGA [11], Adam [24], and the stochastic recursive gradient algorithm (SARAH), [33], possibly augmented with momentum, are preferred for large-scale optimization problems owing to their more affordable computational costs, which are linear in dimensions per epoch O(nd). However, the convergence of the first-order methods is notably slow, and they are sensitive to hyperparameter choices and ineffective for ill-conditioned problems.
In contrast, Newton's method does not depend on the parameters of specific problems and requires only minimal hyperparameter tuning for self-concordant functions, such as 2 -regularized logistic regression. However, Newton's method involves a computational complexity of Ω(nd 2 + d 2.37 ) [1] per iteration and thus is not suitable for large-scale settings. To reduce this computational complexity, the sub-sampled Newton's method and random projection (or sketching) are commonly used to reduce the dimensionality of the problem and solve it in a lower-dimensional subspace. The sub-sampled Newton method performs well for large-scale but relatively low-dimensional problems by computing the Hessian matrix on a relatively small sample. However, the approximated Hessian can deviate significantly from the full Hessian in high-dimensional problems. Randomized algorithms [26,35] estimate the Hessian in Newton's method using a random embedding matrix S ∈ R m×n H S (w) := (∇ 2 f (w) 1 2 ) S S∇ 2 f (w) 1 2 . Specifically, their approximation used the square root of the generalized Gauss-Newton (GGN) matrix as a low-rank approximation instead of deriving it from actual curvature information, whereas S is a random projection matrix of size (m × n). Thus, the approximation may deviate significantly from the actual Hessian. Other sketch matrix-based approximations to reduce the computational complexity of the Newton step calculation have been proposed with deterministic updates which are therefore unsuitable for large-scale stochastic settings where n is significantly larger.
The limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm [30] is a widely used stochastic quasi-Newton method. More specifically, it estimates the Hessian inverse using the past difference of gradients and updates. The online BFGS (oBFGS) [37] method is a stochastic version of regularized BFGS and L-BFGS with gradient descent. [25] proposed two variants of a stochastic quasi-Newton method incorporating a variance-reduced gradient. The first variant used a sub-sampled Hessian with singular value thresholding, which is numerically weaker. The second variant used the LBFGS method to approximate the Hessian inverse. Because these quasi-Newton methods approximate the Hessian inverse using first-order gradients, they are often unable to provide significant curvature information. The stochastic quasi-Newton method (SQN) [6] used the Hessian vector product computed on a subset of each mini-batch instead of approximating the Hessian inverse from the difference between the current and previous gradients, as in LBFGS. SVRG-SQN [32] also incorporated variance-reduced gradients. Hence, both SQN and SVRG-SQN perform well with low-dimensional datasets, whereas their computational time cost is drastically increased for high-dimensional datasets.
In this study, we overcome the limitations of the random projection/sketch-based and quasi-Newton methods to further enhance the curvature information used to approximate the Newton step through the k-rank approximation of the Hessian matrix. In contrast to quasi-Newton methods that estimate the curvature information using first-order gradients only, we use Nyström approximation on the partial Hessian matrix (d × m) constructed for the m d randomly selected columns only to estimate the Hessian matrix. Thus, our approximation is also more suitable for high-dimensional cases. Note that the proposed method can be represented in terms of the sketch method; more details can be found in Section 2.2. We used our approximation with stochastic gradient descent (SGD) and stochastic variance-reduced gradient methods (SVRG), and theoretically prove their convergence for convex functions. In addition, we evaluated the performance of the proposed approach on several large-scale datasets with a wide range of dimensions for both convex and nonconvex functions. The experimental results show that the proposed methods consistently performed better than or competitively with existing methods, whereas the behavior of existing methods depended significantly on the problem. Contribution: The contribution of this study are summarized as follows.
• We propose Nyström-approximated Newton sketch-based methods to perform stochastic optimization.
• We show theoretically that the proposed optimization technique can achieve linear and linear-quadratic convergence in a stochastic convex optimization setup for SVRG and SGD, respectively.
• We empirically show that the performance of the proposed method compared favorably with that of existing methods in training deep learning models.

Proposed method
In this section, we propose Nyström stochastic gradient descent (Nyström-SGD) and stochastic variance-reduced gradient (Nyström-SVRG) methods. The Nyström method is a widely used kernel technique. However, surprisingly, Nyström method for second-order stochastic optimization has not been well studied.

Nyström SGD and SVRG methods
In the proposed method, to incorporate a better approximation of the curvature information, we use a low-rank approximation of the Hessian matrix. For example, the Hessian matrix of a binary logistic regression is positive semi-definite. However, if we use the low-rank approximation of the Hessian matrix, the approximated Hessian cannot be invertible. Moreover, the computation of the approximated Hessian is computationally expensive if computed in each step. Thus, we propose the following regularized variant as a Newton-like update.
where v t denotes the appropriate stochastic gradient, N τ is the rank k-approximated Hessian matrix at epoch τ , η ≥ 0 is a hyperparameter, and ρ > 0 is a regularization parameter. In this study, we denote B τ = (N τ + ρI) −1 . As (N τ + ρI) is positive definite, its matrix inverse can be computed.
The key challenge is to estimate a low-rank Hessian matrix for each step. To this end, we employed the Nyström method and computed the approximated Hessian at each epoch τ .
Definition 1 (Nyström approximation). Let H ∈ R d×d be a symmetric positive semi-definite matrix. Then, the k-rank (k d) approximation N k of the matrix H is given by where Z = CU k Σ −1/2 k ∈ R d×k , and C ∈ R d×m is a matrix consisting of m columns (m d) of H, and M k is the best k-rank approximation of M , which is formed by the intersection between those m columns of H and the corresponding m rows of H, and M † k is the pseudo-inverse of M k , and the rank of M is k ≤ m. Note that the number of columns m is a hyperparameter. To obtain the best k rank approximation, it can be computed using the singular value decomposition (SVD) of M k as M k = U k Σ k U k , where U k ∈ R m×k are singular vectors and Σ k ∈ R k×k are singular values. The pseudo-inverse can be computed as At each epoch, we uniformly sample m columns as Ω ⊂ {1, · · · , d}, where m d, unlike all columns, as in [6,43]. Then, we compute a partial Hessian matrix C ∈ R d×m as the Jacobian of ∇f (w) with respect to w Ω . Then, by using Nyström method in (3), the k-rank approximate Hessian can be computed as To efficiently compute the inverse in (2), we compute B τ v t−1 in (2) using the inversion lemma as where Here, (I + 1 ρ Z τ Z τ ) ∈ R k×k , and its inverse can be computed much more quickly than the inverse of (N τ + ρI) directly. Note that estimating the Hessian matrix via a low-rank approximation is a common technique, e.g. [10]. The key contribution here is the use of the Nyström approximation for the Hessian approximation.

Nyström approximation as Newton-sketch
In this section, we show that Nyström approximation can be interpreted as a Newton sketch-based method [35,26]. This sketching representation of the Nyström approximation is a key to analyze the convergence of the Nyström-SGD method. See Section 3 for a theoretical analysis of Nyström-SGD.
Let H = ∇ 2 f (w) be a Hessian of f (w) of the form H = X X, where X is an n × d matrix. It is always possible to assume that H = X X because H is a symmetric positive semi-definite (SPSD). The embedding W ∈ R d×m can be constructed as follows.
Moreover, C := HW , which is the sampled column matrix of the true Hessian, and M := W HW , which is the intersection matrix in (3). We define XW = " V . Then, we obtain, The right-hand side of (6) is similar to the Newton sketch with two differences, in that 1) X is replaced with the square root of the GGN matrix, and 2) the natural orthogonal matrix " U k in proposed method is replaced by a randomized embedding matrix S , which is expected to be orthogonal in principle, whereas in practice, it may deviate significantly.
If we let X = ∇ 2 f (w) 1/2 then, our approximation is of the form of More generally, the approximation above can be written in the form of an embedding matrix as follows. Let G = ρI d , and Then, by defining the embedding matrixS = ï S m×n 0 m×d 0 d×n I d ò and partial Hessian which is identical to the (7) and hence B τ = H −1 τ,S is non-singular, where H τ,S is the Nyström approximation H S at iteration τ .

Relation to 2 regularization:
The Hessian of 2 -regularized function, i.e., f (w) + λ 2 w 2 can be given as H = H f + λI, where λ ≥ 0. In the Nyström approximation for the regularized objective function, we approximate the H by the regularized Nyström approximation by adding an independent regularizer ρ: That is, one can treat ZZ T + (λ + ρ)I as the approximated Hessian of f (w) + (λ+ρ) 2 w 2 2 , which makes approximation invertible even in the case where λ = 0. Usually, ρ is obtained by trust region techniques which also depends on user defined constant. Therefore, we treat ρ as a hyperparameter and we can obtain the best ρ using grid search.

Theoretical Analysis
We show the theoretical analysis of Nyström-SGD and Nyström-SVRG methods for convex functions. Note that we provide a theoretical analysis for the objective function f given in (1), which does not include the 2 regularization.

Assumptions and Lemmas
In this section, we present a convergence analysis, which involves the following assumptions. Assumption 1. We assume that each f i is convex and twice continuously differentiable. Assumption 2. There exists two positive constants µ and Λ such that The lower bound always holds in the regularized case. Assumption 3. We assume that the gradient of each f i is Λ-Lipschitz continuous, i.e., Under this assumption, it is clear that ∇f is also Λ-Lipschitz continuous.
Theorem 1 provides the error bound for the Nyström approximation with respect to the best k-rank approximation.
Next, we connect the Nyström error bound (Theorem 1) with the upper bound of the Nyström approximation. Lemma 1. Let N = ZZ satisfy Theorem 1: Then, Nyström approximation matrix N is bounded, as follows: Then, we provide the bound of the regularized Nyström approximation with the help of the Assumption 1 and 2.
The following lemma provides the bounds to the inverse approximation B = (N + ρI) −1 using the above lemma.
The following lemma provides the upper bound of the variance of the variance-reduced gradient: Then, the expectation with respect to B is bounded; This bound follows from the [32, Lemma 6] which closely follows from [21,Theorem 1]. At this point, we apply Lemma 5, which states a result of a strongly convex function.
Lemma 5. Suppose that f is continuously differentiable and strongly convex with parameter µ. Let w * be a unique minimizer of f . Subsequently, for any w ∈ R, we have

Convergence analysis of Nyström SGD
Here, we show the convergence of Algorithm 1.
Assumption 4. We assume that the Hessian of the objective function is Lipschitz continuous with Lipschitz constant L H , such that Because the approximation (7) is similar to that of the Newton sketch, we follow a similar analysis for the Newton-like iteration for the Nyström SGD algorithm. We consider the Newton-like iteration: as Nyström-SGD is a specific case of Newton-like iteration. Note that Algorithm 1 is a special case of (11) where B τ is constant for iterations.
Definition 2 (Gaussian width). For a d-dimensional compact set C, Gaussian width is given by where Let L is constraint set for the given problem f (w * ) and K be a tangent cone at w * , such that Note that in the case of unconstrained optimization problem, L = R d and therefore K = R d . Sketch dimension m depends on the Gaussian width. Since, our approximation is similar to Newton-sketch, we take the same lower bound for m, whereˆ is user defined tolerance and c is constant. The next theorem shows linear-quadratic convergence for the Nyström SGD algorithms.
, and H s be an approximation given in . Consider the Nyström sketch update (11) initialized with w 0 such that w 0 − w * ≤ γ 8L H , and a sketch dimension m satisfy the lower bound Then with probability at least 1 − c 1 T e −c2m , the Eucliden error satisfies the bound where T is the total iteration, W is the Gaussian width and K is the tangent cone at w * .
Since the Nyström approximation can be interpreted as a Newton sketch (8), we can obtain Theorem 2 by using [35,Theorem 3.1]. Note that we can use the same lower bound on m as our embedding matrix S has similar properties as of [35].

Convergence analysis of Nyström-SVRG
The next theorem states the main result for the convergence of Algorithm 2. We follow an analysis similar to that of [32].
Theorem 3. Suppose Assumptions 1 -3, and Lemma 1 hold, and let w * be a unique minimizer of the objective function (1).

Closeness to Newton's method
Let H be the Hessian of the 2 -regularized objective function where λ is 2 -regularizer. Then, the inverse of Hessian of is given by H −1 w = (∇ 2 f (w) + λI) −1 . Let the approximated Nyström at w be given by (Z w Z w + λI) −1 . The distance of the regularized inverse matrix is then given as where  (20), where v ny and v ne are the search direction of the Nyström-SGD and Newton's method, respectively. This implies that superlinear convergence is quite difficult in the case of Nyström approximation. This problem remains as an interesting future work.

Related work
Newton methods: The evolution of the approximation of the Hessian or its inverse began with DFP [9,15], Broyden [4], SR1 [7], and well-known BFGS methods [5,14,17,39]. The BFGS method uses secant equations to approximate the Hessian inverse, and requires O(d 2 ) memory to store the approximate Hessian. [30] proposed a limited-memory BFGS (L-BFGS) to solve large-scale unconstrained optimization problems. They used L-BFGS to compute an approximation of the Hessian inverse by using O(md) memory instead of O(d 2 ), where m d. To address the increased demand for optimization in machine and deep learning, various first-order stochastic optimization methods [21,13,11,24,33] have been proposed, but their sensitivity to the choice of hyperparameters remains a concern. Owing to the low sensitivity of the choice of hyperparameters, stochastic quasi-Newton methods have gained popularity in recent years. A stochastic quasi-Newton method by [6] proposed an approximated sub-sampled Hessian using the Hessian vector product instead of using the difference of gradients in the secant equation. [32] proposed a stochastic L-BFGS using a similar technique along with mini-batch variance-reduced gradients. [25] proposed a strategy to compute the sub-sample Hessian followed by the singular value thresholding. As a generalization of quasi-Newton methods, [38] proposed the use of robust symmetric multi-secant updates to approximate the Hessian.
Recently, [43] approximated the Hessian inverse B by decomposing it into two parts such that B = B 1 + B 2 , where B 1 is either a sub-sampled Hessian, generalized Gauss-Newton matrix, Fisher information matrix (FIM), or any other approximation and B 2 is the L-BFGS approximation. Observe that using the sub-sampled Hessian can efficiently compute B 1 for low-dimensional problems. However, for high-dimensional data where d is significantly large, computing the subsample Hessian or FIM B 1 can be computationally expensive because it requires O(nd 2 ) time and O(d 2 ) memory. Additionally, computing the sub-sample Hessian or FIM at each iteration makes it more computationally intensive. AdaHessian [44] is an alternative Hessian approximation method, which estimates the diagonal of Hessian matrix with Hutchinson's method. In contrast, our proposed methods approximates the full Hessian matrix by using the Nyström method. [40] proposed the Nyström logistic regression algorithm, where the Nyström method is used to approximate the Hessian of the regularized logistic regression. Thus, it can be regarded as a variant of Nyström-SGD. However, [40] only considered the regularized logistic regression, in which the Hessian can be explicitly obtained, with deterministic optimization. In contrast, we propose the Nyström method for a general Hessian matrix for stochastic optimization and show its efficacy for deep learning models. Moreover, we elucidate the theoretical properties of the Nyström-SGD algorithm for convex functions (See Appendix B for details). Natural gradient methods: The Newton method is closely related the natural gradient descent [2] methods, which can be used for probabilistic models. More specifically, in probabilistic models, because the negative Hessian is equivalent to FIM, we can update the Newton-like algorithm using FIM. However, FIM computations are expensive. To handle this computational problem, several structural approximations of FIM have been developed [18,31,42,29,28,22,10]. However, these approximation methods are based on FIM and focus more on the computational aspects. In contrast, in this study, we propose the use of the Nyström approximation of the Hessian matrix and elucidate its theoretical properties for Newton-like algorithms. Newton-sketch method: Because the cost of Hessian computation and its inverse is not affordable, some significant work has been performed on sketch-based methods. [26] and [35] proposed the use of Newton-sketch-based methods. They have used various embedding techniques such as a subsampled randomized Hadamard transform (SRHT), sub-Gaussian, randomized orthogonal systems (ROS), and sparse Johnson-Lindenstrauss transform (SJLT). It should be noted that the approximation of the Hessian matrix and its accuracy highly depend on the embedding matrix. One of the drawbacks of the sketching method is that it needs to compute the square root of the Hessian matrix although the size of the embedding matrix is m × n, which makes it impractical for both large-scale and high-dimensionality data. In contrast, the Nyström approximation only needs to compute the partial Hessian matrix. It is important to note that Adaptive Newton sketch [26] increases the m size eventually, whereas the proposed method does not change the value of m through out the algorithm.

Experiments
In this section, we validate our proposed method for both convex and non-convex setups.

Real-world experiments (Convex setup)
We evaluated the performance of the proposed methods for logistic regression and 2 -svm on several benchmark datasets (See Table 2 in the supplementary material). Experimental Setup: We compared the proposed methods with existing state-of-the-art first and second order optimization methods, namely, Adam, SVRG-LBFGS [25], oBFGS [37], SVRG-SQN [32], and SQN [6]. We compare various variants of the proposed methods, namely, Nyström SGD (Nys-SGD) and Nyström SVRG (Nys-SVRG) for optimization error and classification accuracy on both training and test data. Nys-SGD rho=0.01 Nys-SGD rho=0.1 Figure 3: Effect of the ρ on the test accuracy for ResNet18 on imagenet dataset. x-axis is the number of epochs.
We demonstrate the performance of the proposed and existing methods on the 2 -regularized logistic regression problem with the regularizer λ ∈ {10 −2 , 10 −3 , 10 −4 }. For each λ, all methods were tuned for other hyperparameters such as η ∈ {10 0 , 10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 }, batch size, epochs, etc. The parameter m = 50 and ρ ∈ {10 0 , 10 −1 , 10 −2 , 10 −3 } were specific to the proposed variants only. The memory used in the quasi-Newton method was set to 20, which is a commonly used value [25,6]. We report the optimization error on the training set (Opt. Error) and testing set (test error) with respect to epochs and training CPU time cost per epoch. The best-performing model was selected based on the minimum optimization error on the training set and presented its corresponding validation error. We implemented the existing and proposed methods in MATLAB using the SGDLibrary [23]. We computed the results on an Intel(R) Xeon(R) CPU E5-2698v4 @ 2.20GHz with 40 cores running MATLAB R2019a.  Results: Figure 1 show the comparisons of the proposed methods with Adam and existing quasi-Newton methods. We present the numerical results per epoch as well CPU time. It may be observed from the results that the proposed methods performed better than or very competitively with the existing methods consistently on all the datasets with the varying number of dimensions and samples. Additionally, the performance of the proposed methods was almost stable when increasing the regularization parameters, whereas that of the existing methods changed significantly. The existing methods such as SVRG-LBFGS and OBFGS sometimes performed very competitively to the proposed methods. However, their performance was not very stable as they could fluctuate significantly even after reaching very close to the optimal point. It may be observed from the results that the proposed methods generated robust curvature information using only 50 columns of Hessian in the Nyström approximation, compared to the existing quasi-Newton methods using traditional secant equation. As we took the initial point w 0 to be the solution of the least squares problem, the first-order method Adam was unable to significantly reduce the optimization cost. It may be observed that existing stochastic methods require long CPU times for problems with higher dimensionality. More precisely, for the real-sim dataset, SQN and SVRG-SQN require a considerable period of time to reach the optimal point. (See additional results for logistic regression and 2 -svm in appendix)

Real-world experiment (Non-convex setup)
We also evaluated the performance of the Nyström SGD on the well-known deep models on Imagenet dataset. Experimental Setup: We compared our method with the first-order methods SGD and the well-known approximate second-order method KFAC on ResNet152 [19] and EfficientNet [41] models. For Nyström SGD, we used ρ = 0.1 and fixed the m = log 2 |w|, where |w| is the number of parameters in the respective model. We used a batch size of 128. We used a random sample of size of min{6400, n × 0.01} to compute the partial Hessian C for Nyström SGD. The update frequency used to re-estimate the preconditioner in KFAC and its variants is set to 200, as used in their experiments. The ImageNet results were computed on a Quadro RTX 8000 GPU. Results: Figure 2 present the results of the ResNet152 and EfficientNet on the ImageNet dataset. The proposed method outperformed both the SGD and KFAC for both the models in terms of training loss as well as test accuracy, showing the better optimization and generalization ability of the trained models. Table 1 shows the computational time comparison of methods. The per update computational time of the Nyström SGD on ResNet152 is 1.703 seconds which is slightly slower than SGD and KFAC. To further speed up the Nyström SGD is an interesting future work. Figure 3 shows the effect of the ρ parameter for ResNet18. As can be seen, the ρ parameter affects the model performance. We found setting ρ = 0.1 performs well in practice.

Conclusion
In this paper, we propose to approximate the Hessian matrix using Nyström method and to use the approximated Hessian for stochastic optimization. In terms of optimal cost, the proposed Nyström based methods performed better than the first-order methods Adam, SGD, and SVRG, and performed better than or comparably with the stochastic quasi-Newton methods SVRG-LBFGS, oBFGS, SVRG-SQN, and SQN. In terms of computational time, our methods performed compares favorably with the stochastic quasi-Newton methods. Additionally, experimental results on several publicly available benchmark datasets show that the proposed method exhibited much more stable and generalized behavior in comparison to both firstand second-order methods, which can be attributed to the better preservation of the curvature information through Nyström approximation of the Hessian matrix.

A.3 Proof of Lemma 3
Proof. From the Lemma 2, it is easy to see that,

A.4 Proof of Lemma 5
Proof. From the strong convexity of f , where the last inequality holds from the minimizer z = w − 1 µ ∇f (w) of quadratic model: By setting z = w * in (22) gives

A.5 Proof of Theorem 3
Proof. Using Assumption 3, we get Taking expectation on f , where the last inequality follows from the Lemma 5 and Lemma 4 for the second and the third term, respectively. Now, taking sum over t = 1, . . . , , we get Rearranging above inequality gives The second inequality follows from the fact that f (w * ) ≤ f (w ). From the assumption of η < µ∆/2Λ 2 δ 2 , which implies where, Choosing and η to satisfy (17), it follows that the α < 1.      To show the quality of the approximate curvature information used in the proposed methods, we present a comparison of the proposed methods with first-order methods (SGD and SVRG) used as the base for our algorithms in Figure 7 on w8a and adult datasets. The difference between SGD and Nys-SGD on both datasets clearly shows the benefits of using the approximate curvature information and the quality of the approximation. methods, Figure 7 also presents a performance comparison with the recently proposed state-of-the-art S4QN method [43]. For S4QN, we used two variants with sub-sampled Newton (S4QN-N) and Fisher information matrix (S4QN-F) as the base matrix and set r 1 = 0.1 and r 2 = 10. The time taken in the initial epoch by S4QN is significantly higher, as it computes Hessian on a sub-sample at each update. However, the use of a growing gradient batch size reduces the frequency of updates per epoch to speed up subsequent epochs, but decreases the quality of curvature information. Hence, no significant improvement may be observed in the optimization error. However, the proposed Nys-SGD and Nys-SVRG outperformed both variants of the S4QN, indicating that the proposed methods use a better approximation of the curvature information.

I Quality of Hessian approximation
In order to show the difference between the adaptive Newton sketch [26] and the Nyström approximation, we compute their Hessian approximations for logistic regression model on w8a dataset with same w. The adaptive Newton sketch uses the concept of effective dimension. d λ = trace(∇ 2 f (w)(∇ 2 f (w) + λI) −1 ).
First, we computed the effective dimensions for various λ = [10 −5 , 10 −4 10 −3 , 10 −2 , 10 −1 , 10 0 ]. Next, we compute the adaptive Newton sketch where S ∈ R m×n and ∇ 2 f (w) 1/2 ∈ R n×d . Since, their method is deterministic, to have a fair comparison, we compare the approximation only. We computed the adaptive Newton sketch H S with m = eff(d λ ) and H S1 with m = 10 eff(d λ ) using the randomized orthogonal matrix. Then, we computed the Nyström N = ZZ (without using λ and ρ) as in (3)  , where H = ∇ 2 f (w) is computed without using λ. Figure (11) shows the approximations and its difference with original Hessian.
It may be easily observed that the effective dimension can fluctuate with a larger λ and becomes insipid to approximate the Hessian of ill-condition or high-dimensional problems. Moreover, Nyström approximation showed that small column matrix with m = 10 is sufficient to approximate the Hessian with lowest error.