extfieldsofstudy:
Mathematics
Computer Science

text:
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.


s2fieldsofstudy:
Computer Science
Business

text:
Research on the Model and Characteristics of Underground Magnetic Induction Communication Channel

The traditional electromagnetic wave wireless communication in the underground environment has the problem of unstable channel path loss, large antenna size, high path loss, etc. To address these issues, the channel models of magnetic induction communication and magnetic induction waveguide communication based on quasi-static field coupling are proposed, and the characteristics of magnetic field strength, path loss, bandwidth, and channel capacity are analyzed in detail. The results show that the magnetic induction communication system channel is stable, compared with the ordinary induction communication, and the path loss of magnetic induction waveguide communication is reduced a lot, even in the case of high noise and transmission distance increased by more than 20 times. But the bandwidths of the two ways are small and similar. The path loss and bandwidth decide the system capacity, and system capacity is also affected by the number of turns, working frequency, coil resistance, and size.

respectively, and r is the distance between them. If the sinusoidal current is fed into the transmitting coil, that is, I = I 0 · e −jwt , the receiving coil generates sinusoidal current due to the interaction between two coupling coils, thus realizing information transmission. As shown in Figs. 1(b) and (c), two equivalent circuit models are established [8][9][10][11][12][13].
In model 1 in Fig. 1(b) [8]: where R t and R r are resistances of transmitting and receiving coils, respectively; L t and L r are selfinductances of the two; Z t and Z r are self-impedances of the two; Z rt is the influence of receiving coil on transmitting coil; Z tr is the influence of transmitting coil on receiving coil; M is mutual inductance between the two; Z l is the load impedance; U s is the transmitter voltage; U r is the inductive voltage of receiving coil. The transmitting power P t is the power consumed by the main circuit, and the receiving power P r is the power consumed by the load impedance Z l . To maximize the receiving power, the load impedance is designed as the complex conjugate of the output impedance of the secondary circuit. When a low resistance loop with high frequency and multiple turns is adopted, the receiving power P r is [8]: where R 0 is the resistance per unit length of the wire, μ the permeability of the earth, and N r the number of turns of the transmitting loop. In model 2 in Fig. 1(c), C t and C r are capacitances of transmitting and receiving coils, respectively; R l is the load impedance; R s is the power resistance. The receiving power P r is [11][12][13]: where Q t = ωL t /R t and Q r = ωL r /R r are the quality factors of transmitting and receiving coils, respectively; η t = R s /(R s + R t ) and η r = R l /(R l + R r ) are the efficiencies of transmitting and receiving coils, respectively. If R s R t and R l R r , the transmitting coil radius a t is far less than the communication distance r, then the received power P r can be expressed as [10]:

Magnetic Induction Waveguide Communication Channel Model
As shown in Fig. 2, the magnetic induction waveguide is modeled as a multi-stage transformer, and the distance r between the two coils is much greater than the radius a of the coils, so just considering the mutual inductance between adjacent coils [9]. By adding capacitance to each coil, the resonant coil can be formed by designing the appropriate capacitance value, which can effectively transmit the signal. If the sending and receiving coils and relay coil all use the same design parameters such as number of turns, impedance, matching capacitance, and self-induction, then the power P r of the receiving coil is [8,9]: where Z (i−1)i is the reflection impedance generated by the (i − 1) coil on the i coil, L the self-inductance of coils, R the impedance of coils, C the matching capacitance, N the winding turns, n the total number of coils, and {b i , i = 0, 1, 2, 3, . . . , n − 1} are the polynomial coefficients.

CHANNEL CHARACTERISTIC RESEARCH AND SIMULATION ANALYSIS.
This paper uses MATLAB to simulate and analyze the channel characteristics of magnetic induction and magnetic induction waveguide communication system. Assuming that all the coil parameters are the same, set radius of 0.1 m, coil number of turns N of 10, per unit length of wire coil impedance R 0 of 0.02 Ω/m, frequency f of 10 MHZ, and current I of 1 A. The observation point P is located at the altitude h 0 and the radial coordinate ρ, and the magnetic field generated by the buried circular coil is [14]:

Magnetic Field Intensity
where S is the coil area, and P 0 and Q 0 are transmission losses caused by attenuation of the conductive ground [15].
where D = ρ/h, Z = z/h, and H = (μωσ) 1/2 h. H is the ratio of source depth to skin depth. J 0 and J 1 are Bessel functions of the first kind. The observation point Q is located at depth h (z = −h), the radial coordinate ρ, where the magnetic field generated by the circular coil on the ground is [14]: where where A = a/h, and Z = h 0 /h. To observe the difference between the vertical component of the magnetic field generated by the uplink loop and the downlink loop, both loops measure the vertical component of the magnetic field at the same distance (10 m) from the transmitting loop. Fig. 4 is obtained from Equations (7) and (10), and Fig. 4 shows the magnitude of the vertical component of the magnetic field generated by the underground loop relative to the surface loop when the frequency is 1-10 MHz under the condition of low and high earth conductivities.
As can be seen from Fig. 4, under the same transmission distance, the vertical component of the magnetic field generated by the underground loop is larger. For the earth with low conductivity, the difference between the two is small, but for the earth with high conductivity, the difference is large.

Path Loss and Link Budget
The transmitting power of the magnetic induction communication system consists of the induction power consumed by the receiving coil and the power consumed by the coil resistance, since the magnetic induction communication is realized through the near-field coupling. The coil resistance is small, and the radiated power is small, so the limited transmitted energy will not be wasted on the radiation in the surrounding space, and most of the power is transferred to the receiving coil.
According to Equations (2) and (4), the path losses of magnetic induction model 1 and model 2 respectively are: The path loss of electromagnetic wave in soil is [16]: L EM (r) = −10 lg P r P t 6.4 + 20 lg r + 20 lg β + 8.69αr (15) where . The result of MATLAB calculation from Equations (13)-(15) is shown in Fig. 5. Fig. 5 shows that the two magnetic induction models have similar path losses, and the path loss of magnetic induction communication is approximately a lg function of distance r, while the path loss of electromagnetic wave system mainly comes from material absorption loss, so it is an approximately linear function of distance r. The conductivity has a great influence on the propagation of the electromagnetic wave. When the conductivity is large, the path loss of the electromagnetic wave is obviously greater than that of magnetic induction communication. Even in low conductivity environments, when the transmission distance is long enough, the magnetic induction system can have less path loss than the electromagnetic wave.  As can be seen from Fig. 6, with the increase of working frequency, the transmission distance of the magnetic induction system increases, while the transmission distance of electromagnetic wave decreases, and electromagnetic wave works more efficiently at lower frequencies, while magnetic induction works better at higher frequencies. As the number of turns increases, the transmission distance of magnetic induction model 1 increases, while the transmission distance of model 2 decreases. Therefore, model 1 is suitable for the case of multiple turns, while model 2 is more suitable for the case of low turns.
We can know from Equation (5) that for the magnetic induction waveguide system, path loss can usually be reduced by reducing the ratio of the coil distance to the coil radius, increasing the operating frequency or the number of turns of the coil, and reducing the wire resistance. However, reducing the ratio between the coil distance and coil radius will increase the difficulty of equipment deployment, and due to the constraints of matching capacitance, the working frequency and the number of turns of the coil cannot be increased without limit. Reducing the wire resistance will increase the cost, and too low resistance will cause the signal fluctuation in the band; therefore, these parameters should be balanced to achieve the best performance of the system. According to Equation (5), the path loss of the magnetic induction waveguide system can be obtained, as shown in Fig. 7. Figure 7(a) shows that compared with the common magnetic induction system, magnetic induction waveguide path loss can be greatly reduced and can be further reduced by reducing the distance and resistance of the relay coil. Fig. 7(b) shows the path losses of the magnetically induced waveguide system with a different number of repeating coils. It can be found that when the transmission distance is relatively close (< 40 m), increasing the number of repeating coils will have no obvious effect and cause greater path loss. When the transmission distance is large, increasing the number of repeating coils will have a good effect on reducing the path loss.

Bit Error Rate
Bit error rate mainly depends on the three factors of path loss, noise, and system modulation mode.
Assuming that the noise power is 83 dBm in the case of high noise, 103 dBm in the case of low noise, the transmitting power P t is set as 10 dBm, and 2PSK is adopted in the modulation mode. Figures 8(a) and (b) respectively show the bit error rate of magnetic induction system model 1, model 2, and electromagnetic wave system under low noise and high noise conditions. The bit error rate of model 1 is similar to that of model 2, and the bit error rate of model 1 is smaller. In the case of low noise, the transmission distance of the magnetic induction system is greater than that of the electromagnetic wave system regardless of the conductivity. In the case of high noise, the transmission distance of the magnetic induction system is between the electromagnetic wave systems with low conductivity and high conductivity.  Figure 9 shows the relationship between bit error rate and transmission distance of the magnetic induction waveguide system with different coil numbers, resistance, and spacing under high noise. It can be seen from the figure that even in the case of high noise, the transmission range of the magnetic induction waveguide system is more than 20 times higher than that of the magnetic induction system and electromagnetic wave system. By analyzing the bit error rate, the transmission distance can be improved by reducing the distance and resistance of the loop or increasing the number of the loops.

Bandwidth and System Capacity
In a magnetic induction communication system, only a center frequency can realize secondary loop output impedance of the load matching, and any deviation from the central frequency will cause power reflection and increase path loss. Similarly, in a magnetically induced waveguide system, only one central frequency can make all coils reach the resonant state. Any deviation from the central frequency will make the resonant state disappear, and the load does not match the system. Therefore, it is necessary to analyze the bandwidth of the magnetic induction system and magnetic induction waveguide system. Figure 10(a) shows the frequency response of the magnetic induction communication at the transmission distance of 5 m. It can be seen from the figure that the 3 dB bandwidth of the system is about 2 kHz, which is much smaller than the electromagnetic wave system. Fig. 10(b) shows the frequency response of the magnetic induction waveguide system at the transmission distance of 100 m at different relay distances and different line resistances. The results show that the 3 dB bandwidth of the magnetic induction waveguide system is similar to that of the common magnetic induction system. The lower line resistance can reduce the path loss, but it will increase the frequency response fluctuation in the band. The reduction of the relay distance can increase the bandwidth, but it also means that more repeating coils need to be deployed, which increases the cost and difficulty of deployment.
According to Shannon's theorem, channel capacity can be expressed as: where B is the channel bandwidth; 3 dB bandwidth is selected as the channel bandwidth in this paper; P r is the received power; and N 0 is the total noise and interference power. The capacity of the magnetic induction system is [12]: where Q is the quality factor of the coil, and f 0 is the center operating frequency. If the transmitting power is set at 10 dBm, the noise power set at −105 dBm, and the received power given by Equations (2) and (4), then the channel capacities of model 1 and model 2 of the magnetic induction system are shown in Fig. 11. Figure 11(a) shows that with the increase of transmission distance, channel capacity is reduced, that the capacities of the two models are similar, and that model 1 is larger. Fig. 11(b) shows the influence of wire resistance and frequency on channel capacity. The results show that frequency has less influence on channel capacity, while wire resistance has more influence on capacity. Channel capacity can be improved by increasing wire resistance.  In the magnetic induction waveguide system, the path loss at frequency f 0 + 0.5B is twice as much as the path loss at the center frequency f 0 , which can be obtained by formula (5): When the relay coil distance is large, the influence of the highest order variable in the polynomial is the largest, and the bandwidth B is much smaller than the center frequency f 0 . Therefore, Eq. (18) is approximate: Then the bandwidth and capacity of the magnetic induction waveguide system are expressed as [17] where received power P r is given by Equation (5); the transmission power is set to 10 dBm; the noise power is set to −105 dBm. Fig. 12 shows the channel capacity of the magnetic induction waveguide with different coil resistances R and coil spacings r. It can be found that reducing the coil resistance can reduce the path loss, but it will also reduce the bandwidth; increasing the spacing between coils will increase the bandwidth, but also increase the path loss. Therefore, there is an optimal coil resistance and coil spacing to maximize the channel capacity. Figure 13 shows that the relationship between channel capacity and transmission distance of the magnetic induction waveguide system under different coil radii and operating frequencies. Different from  the traditional wireless channel, the signal-to-noise ratio and bandwidth decrease significantly with the increase of transmission distance, so the channel capacity decreases rapidly. Increasing the coil radius can reduce the bandwidth and greatly reduce path loss, especially in the long-distance transmission. Therefore, increasing the coil radius will reduce the channel capacity, but the attenuation speed is slow. Increasing the working frequency can increase mutual inductance and improve the channel capacity to some extent.

CONCLUSION
The Conclusion is based on the numerical comparisons of the EM wave and MI communication system results, illustrated in the figures.
(1) When the transmission distance is the same, the vertical component of the magnetic field generated by the underground loop is larger than that generated by the ground loop. For the earth with low conductivity, the difference between the two is small, but for the earth with high conductivity, the difference is large.
(2) The path losses of the two magnetic induction models are similar. When the conductivity is large, the path loss of the electromagnetic wave is obviously greater than that of magnetic induction communication. Even in a low conductivity environment, when the transmission distance is large enough, the magnetic induction system has a path loss smaller than the electromagnetic wave, and the electromagnetic wave has a higher working efficiency at lower frequencies, while the magnetic induction system is more suitable for working at higher frequencies.
(3) The magnetic induction waveguide technology can greatly reduce the path loss, and the path loss can be further reduced by reducing the distance and resistance of the relay coil. When the transmission distance is relatively short, increasing the number of relay coils will have no obvious effect and cause greater path loss. When the transmission distance is relatively large, increasing the number of relay coils will have a good effect on reducing the path loss.
(4) The bandwidth of the magnetic induction communication system is similar to that of the magnetic induction waveguide system, and both of them are small. The bandwidth of the magnetic induction waveguide system can be increased by reducing the relay distance.
(5) The channel capacities of the two models of magnetic induction system are similar, and the frequency has less influence on the capacity, while the line resistance has more influence on the capacity. The channel capacity can be increased by increasing the line resistance. For the magnetic induction waveguide system, there is an optimal coil resistance and coil spacing to maximize the channel capacity.