Adaptive norms for deep learning with regularized Newton methods

Jonas Kohler, Leonard Adolphs, Aurelien Lucchi

Introduction

We consider finite-sum optimization problems of the form

In the era of deep neural networks, stochastic gradient descent (SGD) is one of the most widely used training algorithms . What makes SGD so attractive is its simplicity and per-iteration cost that is independent of the size of the training set ( $n$ ) and scale linearly in the dimensionality ( $d$ ). However, gradient descent is known to be inadequate to optimize functions that are ill-conditioned and thus adaptive gradient methods that employ dynamic, coordinate-wise learning rates based on past gradients—including Adagrad , RMSprop and Adam —have become a popular alternative, often providing significant speed-ups over SGD.

From a theoretical perspective, Newton methods provide stronger convergence guarantees by appropriately transforming the gradient in ill-conditioned regions according to second-order derivatives. It is precisely this Hessian information that allows regularized Newton methods to enjoy superlinear local convergence as well as to provably escape saddle points . While second-order algorithms have a long-standing history even in the realm of neural network training , they were mostly considered as too computationally and memory expensive for practical applications. Yet, the seminal work of renewed interest for their use in deep learning by proposing efficient Hessian-free methods that only access second-order information via matrix-vector products which can be computed at the cost of an additional backpropagation . Among the class of regularized Newton methods, trust region and cubic regularization algorithms are the most principled approaches in the sense that they yield the strongest convergence guarantees. Recently, stochastic extensions have emerged , which suggest their applicability for deep learning.

We here propose a simple modification to make TR methods even more suitable for neural network training. Particularly, we build upon the following alternative view on adaptive gradient methods:

While gradient descent can be interpreted as a spherically constrained first-order TR method, preconditioned gradient methods—such as Adagrad—can be seen as first-order TR methods with ellipsoidal trust region constraint.

This observation is particularly interesting since spherical constraints are blind to the underlying geometry of the problem, but ellipsoids can adapt to local landscape characteristics, thereby allowing for more suitable steps in regions that are ill-conditioned. We will leverage this analogy and investigate the use of the Adagrad and RMSProp preconditioning matrices as ellipsoidal trust region shapes within a stochastic second-order TR algorithm . Since no ellipsoid fits all objective functions, our main contribution lies in the identification of adequate matrix-induced constraints that lead to provable convergence and significant practical speed-ups for the specific case of deep learning. On the whole, our contribution is threefold:

We provide a new perspective on adaptive gradient methods that contributes to a better understanding of their inner-workings.

We investigate the first application of ellipsoidal TR methods for deep learning. We show that the RMSProp matrix can directly be applied as constraint inducing norm in second-order TR algorithms while preserving all convergence guarantees (Theorem 1).

Finally, we provide an experimental benchmark across different real-world datasets and architectures (Section 5). We compare second-order methods also to adaptive gradient methods and show results in terms of backpropagations, epochs, and wall-clock time; a comparison we were not able to find in the literature.

Our main empirical results demonstrate that ellipsoidal constraints prove to be a very effective modification of the trust region method in the sense that they constantly outperform the spherical TR method, both in terms of number of backprogations and asymptotic loss value on a variety of tasks.

Related work

The prototypical method for optimizing Eq. (1) is SGD . The practical success of SGD in non-convex optimization is unquestioned and theoretical explanations of this phenomenon are starting to appear. Recent findings suggest the ability of this method to escape saddle points and reach local minima in polynomial time, but they either need to artificially add noise to the iterates or make an assumption on the inherent noise of SGD . For neural networks, a recent line of research proclaims the effectiveness of SGD, but the results come at the cost of strong assumptions such as heavy over-parametrization and Gaussian inputs . Adaptive gradient methods build on the intuition that larger (smaller) learning rates for smaller (larger) gradient components balance their respective influences and thereby the methods behave as if optimizing a more isotropic surface. Such approaches have first been suggested for neural nets by and convergence guarantees are starting to appear . However, these are not superior to the $\mathcal{O}(\epsilon_{g}^{-2})$ worst-case complexity of standard gradient descent .

The most principled class of regularized Newton methods are trust region (TR) and adaptive cubic regularization algorithms (ARC) , which repeatedly optimize a local Taylor model of the objective while making sure that the step does not travel too far such that the model stays accurate. While the former finds first-order stationary points within $\mathcal{O}(\epsilon_{g}^{-2})$ , ARC only takes at most $\mathcal{O}(\epsilon_{g}^{-3/2})$ . However, simple modifications to the TR framework allow these methods to obtain the same accelerated rate . Both methods take at most $\mathcal{O}(\epsilon_{H}^{-3})$ iterations to find an $\epsilon_{H}$ approximate second-order stationary point . These rates are optimal for second-order Lipschitz continuous functions and they can be retained even when only sub-sampled gradient and Hessian information is used . Furthermore, the involved Hessian information can be computed solely based on Hessian-vector products, which are implementable efficiently for neural networks . This makes these methods particularly attractive for deep learning, but the empirical evidence of their applicability is rather limited. We are only aware of the works of and , which report promising first results but are by no means fully encompassing.

An interesting line of research proposes to replace the Hessian by (approximations of) the generalized-Gauss-Newton matrix (GGN) within a Levenberg-Marquardt frameworkThis algorithm is a simplified TR method, initially tailored for non-linear least squares problems . As the GGN matrix is always positive semidefinite, these methods cannot leverage negative curvature to escape saddles and hence, there exist no second-order convergence guarantees. Furthermore, there are cases in neural networks where the Hessian is better conditioned than the GGN matrix . Nevertheless, the above works report promising preliminary results, most notably find that K-FAC can be faster than SGD on a small convnet. On the other hand, recent findings report performance at best comparable to SGD on the much larger ResNet architecture . Moreover, reports many cases where TR and GGN algorithms perform similarly. This line of work can be seen as complementary to our approach since it is straightforward to replace the Hessian in the TR framework with the GGN matrix. Furthermore, the preconditioners used in and , namely diagonal estimates of the empirical Fisher and Fisher matrix, respectively, can directly be used as matrix norms in our ellipsoidal TR framework.

An alternative view on adaptive gradient methods

Adaptively preconditioned gradient methods update iterates as ${\bf w}_{t+1}={\bf w}_{t}-\eta_{t}{\bf A}_{t}^{-1/2}{\bf g}_{t},$ where ${\bf g}_{t}$ is a stochastic estimate of $\nabla\mathcal{L}({\bf w}_{t})$ and ${\bf A}_{t}$ is a positive definite symmetric pre-conditioning matrix. In Adagrad, ${\bf A}_{ada,t}$ is the un-centered second moment matrix of the past gradients computed as

where $\epsilon>0$ , $\mathbf{I}$ is the $d\times d$ identity matrix and ${\bf G}_{t}=[{\bf g}_{1},{\bf g}_{2},\ldots,{\bf g}_{t}]$ . Building up on the intuition that past gradients might become obsolete in quickly changing non-convex landscapes, RMSprop (and Adam) introduce an exponential weight decay leading to the preconditioning matrix

where $\beta\in(0,1)$ . In order to save computational efforts, the diagonal versions $\operatorname{diag}({\bf A}_{ada})$ and $\operatorname{diag}({\bf A}_{rms})$ are more commonly applied in practice, which in turn gives rise to coordinate-wise adaptive stepsizes that are enlarged (reduced) in coordinates that have seen past gradient components with a smaller (larger) magnitude.

Starting from the fact that adaptive methods employ coordinate-wise stepsizes, one can take a principled view on these methods. Namely, their update steps arise from minimizing a first-order Taylor model of the function $\mathcal{L}$ within an ellipsoidal search space around the current iterate ${\bf w}_{t}$ , where the diameter of the ellipsoid along a particular coordinate is implicitly given by $\eta_{t}$ and $\|{\bf g}_{t}\|_{{\bf A}_{t}^{-1}}$ . Correspondingly, vanilla (S)GD optimizes the same first-order model within a spherical constraint. Figure 1 (top) illustrates this effect by showing not only the iterates of GD and Adagrad but also the implicit trust regions within which the local models were optimized at each step.We only plot every other trust region. Since the models are linear, the minimizer is always on the boundary.

It is well known that GD struggles to progress towards the minimizer of quadratics along low-curvature directions (see e.g., ). While this effect is negligible for well-conditioned objectives (Fig. 1, left), it leads to a drastic slow-down when the problem is ill-conditioned (Fig. 1, center). Particularly, once the method has reached the bottom of the valley, it struggles to make progress along the horizontal axis. Here is precisely where the advantage of adaptive stepsize methods comes into play. As illustrated by the dashed lines, Adagrad’s search space is damped along the direction of high curvature (vertical axis) and elongated along the low curvature direction (horizontal axis). This allows the method to move further horizontally early on to enter the valley with a smaller distance to the optimizer ${\bf w}^{*}$ along the low curvature direction which accelerates convergence.

Let us now formally establish the result that allows us to re-interpret adaptive gradient methods from the trust region perspective introduced above.

Equivalent results can be established for Adam using ${\bf g}_{adam,t}:=(1-\beta)\sum_{k=0}^{t}\beta^{t-k}{\bf g}_{t}$ as well as for Adagrad by replacing the matrix ${\bf A}_{ada}$ into the constraint in Eq. (LABEL:eq:rms_step_tr). Of course, the update procedure in Eq. (5) is merely a reinterpretation of the original preconditioned update, and thus the employed trust region radii are defined implicitly by the current gradient and stepsize.

2 Diagonal versus full preconditioning

A closer look at Figure 1 reveals that the first two problems are perfectly axis-aligned, which makes these objectives particularly attractive for diagonal preconditioning. For comparison, we report another quadratic instance, where the Hessian is no longer zero on the off-diagonals (Fig. 1, right). As can be seen, this introduces a tilt in the level sets and reduces the superiority of diagonal Adagrad over plain GD. However, using the full preconditioner ${\bf A}_{ada}$ re-establishes the original speed up. Yet, non-diagonal preconditioning comes at the cost of taking the inverse square root of a large matrix, which is why this approach has been relatively unexplored (see for an exception). Interestingly, early results by on the curvature of neural nets report a strong diagonal dominance of the Hessian matrix $\nabla^{2}\mathcal{L}({\bf w})$ . However, the reported numbers are only for tiny networks of at most 256 parameters. We here take a first step towards generalizing these findings to modern day networks. Furthermore, we contrast the diagonal dominance of real Hessians to the expected behavior of random Wigner matrices.Of course, Hessians do not have i.i.d. entries but the symmetry of Wigner matrices suggests that this baseline is not completely off. For further evidence, we also compare Hessians of Ordinary Least Squares (OLS) problems with random inputs. For this purpose, let $\delta_{\bf A}$ define the ratio of diagonal to overall mass of a matrix ${\bf A}$ , i.e. $\delta_{\bf A}:=\frac{\sum_{i}|{\bf A}_{i,i}|}{\sum_{i}\sum_{j}|{\bf A}_{i,j}|}$ as in .

Thus, if we suppose the Hessian at any given point ${\bf w}$ were a random Wigner matrix we would expect the share of diagonal mass to fall with $\mathcal{O}(1/d)$ as the network grows in size. In the following, we derive a similar result for the large $n$ limit in the case of OLS Hessians.

Empirical simulations suggest that this result holds already in small $n$ settings (see Figure 6) and finite $n$ results can be likely derived under assumptions such as Gaussian data. As can be seen in Figure 2 below, even for a practical batch size of $n=32$ the diagonal mass $\delta_{\bf H}$ of neural networks stays above both benchmarks for random inputs as well as with real-world data.

These results are in line with and suggest that full matrix preconditioning might indeed not be worth the additional computational cost. Consequently, we use diagonal preconditioning for both first- and second-order methods in all of our experiments in Section 5. Further theoretical elaborations of these findings present an interesting direction of future research.

Second-order Trust Region Methods

Cubic regularization and trust region methods belong to the family of globalized Newton methods. Both frameworks compute parameter updates by optimizing regularized (former) or constrained (latter) second-order Taylor models of the objective $\mathcal{L}$ around the current iterate ${\bf w}_{t}$ .In the following we only treat TR methods, but we emphasize that the use of matrix induced norms can directly be transferred to the cubic regularization framework. In particular, in iteration $t$ the update step of the trust region algorithm is computed as

where $\Delta_{t}>0$ and ${\bf g}_{t}$ and ${\bf B}_{t}$ are either $\nabla\mathcal{L}({\bf w}_{t})$ and $\nabla^{2}\mathcal{L}({\bf w}_{t})$ or suitable approximations. The matrix ${\bf A}_{t}$ induces the shape of the constraint set. So far, the common choice for neural networks is ${\bf A}_{t}:=\mathbf{I},\;\forall t$ which gives rise to spherical trust regions . By solving the constrained problem (7), TR methods overcome the problem that pure Newton steps may be ascending, attracted by saddles or not even computable. Please see Appendix B for more details.

There are many sources for ill-conditioning in neural networks such as un-centered and correlated inputs , saturated hidden units, and different weight scales in different layers . While the quadratic term of model (7) accounts for such ill-conditioning to some extent, the spherical constraint is completely blind towards the loss surface. Thus, it is advisable to instead measure distances in norms that reflect the underlying geometry (see Chap. 7.7 in ). The ellipsoids we propose are such that they allow for longer steps along coordinates that have seen small gradient components in the past and vice versa. Thereby the TR shape is adaptively adjusted to fit the current region of the loss landscape. This is not only effective when the iterates are in an ill-conditioned neighborhood of a minimizer (Fig. 1), but it also helps to escape elongated plateaus (see autoencoder in Sec. 5). Contrary to adaptive first-order methods, the diameter ( $\Delta_{t}$ ) is updated directly depending on whether or not the local Taylor model is an adequate approximation at the current point.

1 Convergence of ellipsoidal Trust Region methods

Inspired by the success of adaptive gradient methods, we investigate the use of their preconditioning matrices as norm inducing matrices for second-order TR methods. The crucial condition for convergence is that the applied norms are not degenerate during the entire minimization process in the sense that the ellipsoids do not flatten out (or blow up) completely along any given direction. The following definition formalizes this intuition.

Consequently, the ellipsoids ${\bf A}_{rms,t}$ can directly be applied to any convergent TR framework without losing the guarantee of convergence (, Theorem 6.6.8).Note that the assumption of bounded batch gradients, i.e. smooth objectives, is common in the analysis of stochastic algorithms . In Theorem 1 we extend this result by showing the (to the best of our knowledge) first convergence rate for ellipsoidal TR methods. Interestingly, similar results cannot be established for ${\bf A}_{ada,t}$ , which reflects the widely known vanishing stepsize problem that arises since squared gradients are continuously added to the preconditioning matrix. At least partially, this effect inspired the development of RMSprop and Adadelta .

2 A stochastic ellipsoidal TR framework for neural network training

Since neural network training often constitutes a large-scale learning problem in which the number of datapoints $n$ is high, we here opt for a stochastic TR framework in order to circumvent memory issues and reduce the computational complexity. To obtain convergence without computing full derivative information, we first need to assume sufficiently accurate gradient and Hessian estimates.

The approximations of the gradient and Hessian at step $t$ satisfy

where $\delta_{g}\leq\frac{(1-\eta)\epsilon_{g}}{4}$ and $\delta_{H}\leq\min\left\{\frac{(1-\eta)v\epsilon_{H}}{2},1\right\}$ , for some $0<v<1$ .

For finite-sum objectives such as Eq. (1), the above condition can be met by random sub-sampling due to classical concentration results for sums of random variables . Following these references, we assume access to the full function value in each iteration for our theoretical analysis but we note that convergence can be retained even for fully stochastic trust region methods and indeed our experiments in Section 5 use sub-sampled function values due to memory constraints. Secondly, we adapt the framework of , which allows for cheap inexact subproblem minimization, to the case of iteration-dependent constraint norms (Alg. 1).

Each update step ${\bf s}_{t}$ yields at least as much model decrease as the Cauchy- and Eigenpoint simultaneously, i.e. $m_{t}({\bf s}_{t})\leq m_{t}({\bf s}_{t}^{C})$ and $m_{t}({\bf s}_{t})\leq m_{t}({\bf s}_{t}^{E}),$ where ${\bf s}_{t}^{C}$ and ${\bf s}_{t}^{E}$ are defined in Eq.(28).

Finally, given that the adaptive norms induced by ${\bf A}_{rms,t}$ satisfy uniform equivalence as shown in Lemma 1, the following Theorem establishes an ${\cal O}\left(\text{max}\left\{\epsilon_{g}^{-2}\epsilon_{H}^{-1},\epsilon_{H}^{-3}\right\}\right)$ worst-case iteration complexity which effectively matches the one of .

Assume that $\mathcal{L({\bf w})}$ is second-order smooth with Lipschitz constants $L_{g}$ and $L_{H}$ . Furthermore, let Assumption 1 and 2 hold. Then Algorithm 1 finds an $\mathcal{O}(\epsilon_{g},\epsilon_{H})$ first- and second-order stationary point in at most $\mathcal{O}\left(\max\left\{\epsilon_{g}^{-2}\epsilon_{H}^{-1},\epsilon_{H}^{-3}\right\}\right)$ iterations.

The proof of this statement is a straight-forward adaption of the proof for spherical constraints, taking into account that the guaranteed model decrease changes when the computed step $s_{t}$ lies outside the Trust Region. Due to the uniform equivalence established in 1, the altered diameter of the trust region along that direction and hence the change factor is always strictly positive and finite.

Experiments

To validate our claim that ellipsoidal TR methods yield improved performance over spherical ones, we run a set of experiments on two image datasets and three types of network architectures. All methods run on (almost) the same hyperparameters across all experiments (see Table 1 in Appendix B) and employ the preconditionied Steihaug-Toint CG method to solve the subproblems (Eq. 50) with the classical stopping criterion given in Eq. (52).

As depicted in Fig. 3, the ellipsoidal TR methods consistently outperform their spherical counterpart in the sense that they reach full training accuracy substantially faster on all problems. Moreover, their limit points are in all cases lower than those of the uniform method. Interestingly, this makes an actual difference in the image reconstruction quality of autoencoders (see Figure 12), where the spherically constrained TR method struggles to escape a saddle. We thus draw the clear conclusion that the ellipsoidal constraints we propose are to be preferred over spherical ones when training neural nets with second-order methods. More experimental and architectural details are provided in App. C.

To put the previous results into context, we also benchmark several state-of-the-art gradient methods. For a fair comparison, we report results in terms of number of backpropagations, epochs and time. All figures can be found in App. C. Our findings are mixed: For small nets such as the MLPs the TR method with RMSProp ellipsoids is superior in all metrics, even when benchmarked in terms of time. However, while Fig. 9 indicates that ellipsoidal TR methods are slightly superior in terms of backpropagations even for ResNets and Autoencoders, a close look at the Figures 10 and 11 (App. C) reveals that they at best manage to keep pace with first-order methods in terms of epochs and are inferior in time. Furthermore, only the autoencoders give rise to saddles, which adaptive gradient methods escape faster than vanilla SGD, similarly to the case for second-order methods in Fig. 3.

Conclusion

We investigated the use of ellipsoidal trust region constraints for neural networks. We have shown that the RMSProp matrix satisfies the necessary conditions for convergence and our experimental results demonstrate that ellipsoidal TR methods outperform their spherical counterparts significantly. We thus consider the development of further ellipsoids that can potentially adapt even better to the loss landscape such as e.g. (block-) diagonal hessian approximations (e.g. ) or approximations of higher order derivatives as an interesting direction of future research.

Yet, the gradient method benchmark indicates that the value of Hessian information for neural network training is limited for mainly three reasons: 1) second-order methods rarely yield better limit points, which suggests that saddles and spurious local minima are not a major obstacle; 2) gradient methods can run on smaller batch sizes which is beneficial in terms of epoch and when memory is limited; 3) The per-iteration time complexity is noticeably lower for first-order methods (Figure 11). These observations suggest that advances in hardware and distributed second-order algorithms (e.g., ) will be needed before Newton-type methods can replace gradient methods in deep learning.

As a side note, we reported elevated levels of diagonal dominance in neural network architectures, which may partially explain the success of diagonal preconditioning in first-order method. Further empirical and theoretical investigations of this phenomenon with a particular focus on layer-wise dependencies constitute an interesting direction of future research as for example algorithms such as K-FAC seem to achieve good results with block-diagonal preconditioning.

Broader impact

We consider our work fundamental research with no specific application other than training neural networks in general. Hence a broader impact discussion is not applicable.

References

Appendix A: Proofs

Appendix B Equivalence of Preconditioned Gradient Descent and first-order Trust Region Methods

Let $L({\bf s},\lambda)$ denote the Lagrange dual of Eq. (5)

Any point ${\bf s}$ is a KKT point if and only if the following system of equations is satisfied

For ${\bf s}_{t}$ as given in Eq. (4) we have that

and thus 13 and 14 hold with equality such that any $\lambda\geq 0$ is feasible. Furthermore,

is zero for $\lambda=\|{\bf g}_{t}\|_{{\bf A}^{-1}}\geq 0$ . As a result, ${\bf s}_{t}$ is a KKT point of the convex problem 5 which proves the assertion.

To illustrate this theoretical result we run gradient descent and Adagrad as well as the two corresponding first-order TR methodsEssentially Algorithm 1 with $m_{t}$ based on a first order Taylor expansion, i.e. $m_{t}^{1}({\bf s})$ as in Eq. (10). on an ill-conditioned quadratic problem. While the method 1st TR optimizes a linear model within a ball in each iteration, 1st TR ${}_{\text{ada}}$ optimizes the same model over the ellipsoid given by the Adagrad matrix ${\bf A}_{ada}$ . The results in Figure 4 show that the methods behave very similarly to their constant stepsize analogues.

Appendix C Convergence of ellipsoidal TR methods

At a high level, the proof can be divided into two steps: 1) establish that each update decreases the model value and 2) relate the model decrease to the function decrease, therefore proving that the function decreases.

Based on Assumption 2 the proof first relates the model decrease in each iteration to the gradient norm $\|{\bf g}_{t}\|$ and the magnitude of the smallest eigenvalue $|\lambda_{\min}({\mathcal{B}}_{t})|$ as well as $\Delta_{t}$ . In the case of interior solutions ( $\|{\bf s}_{t}\|_{t}<\Delta$ ), nothing changes compared to spherical Trust Region methods. When the computed step $s_{t}$ lies outside the Trust Region, however, the guaranteed model decrease changes by a constant factor, which accounts for the altered diameter of the trust region along that direction. Due to the uniform equivalence established in 1 this factor is always strictly positive and finite.

More specifically, the first step of the proof relies on Assumption 2 in order to relate the model decrease at each iteration $t$ to three quantities of interest: i) the gradient norm $\|{\bf g}_{t}\|$ , ii) the magnitude of the smallest eigenvalue $|\lambda_{\min}({\mathcal{B}}_{t})|$ , and iii) the trust region radius $\Delta_{t}$ . In the case of interior solutions ( $\|{\bf s}_{t}\|_{t}<\Delta$ ), the model decrease is shown as in the spherical Trust Region methods. When the computed step $s_{t}$ lies outside the Trust Region, however, the guaranteed model decrease changes by a constant factor, which accounts for the altered diameter of the trust region along that direction. Due to the uniform equivalence established in Lemma 2 this factor is always strictly positive and finite.

From here on, the proof proceeds in a standard fashion (see e.g. ). That is, a lower bound on $\Delta_{t}$ is established which (i) upper bounds the number of unsuccessful steps and (ii) lower bounds the guaranteed model decrease introduced above, which in turn allows to bound the overall number of successful steps as a fraction of the initial suboptimality in $\mathcal{L}$ . Assumption 1 together with the smoothness assumptions on $\mathcal{L}$ allow to finally relate the progress of each successful step to the actual function decrease. Finally, since the function decreases and because it is lower bounded, there is a finite number of steps which we can upper bound.

C.2 Proof

In order to prove convergence results for ellipsoidal Trust Region methods one must ensure that the applied norms are coherent during the complete minimization process in the sense that the ellipsoids do not flatten out (or blow up) completely along any given direction. This intuition is formalized in Assumption 1 which we restate here for the sake of clarity.

There exists a constant $\mu\geq 1$ such that

Towards this end, identify the following sufficient condition on the basis of which we will prove that our proposed ellipsoid ${\bf A}_{rms}$ is indeed uniformly equivalent under some mild assumptions.

Suppose that there exists a constant $\zeta\geq 1$ such that

Having uniformly equivalent norms is sufficient to prove convergence of ellipsoidal TR methods (se AN.1 and Theorem 6.6.8 in ). However, it is so far unknown how the ellipsoidal constraints influence the convergence rate itself. We here prove that the specific ellipsoidal TR method presented in Algorithm 1 preserves the rate of its spherically-constrained counterpart proposed in (see Theorem 1 below).

First, we show that the proposed ${\bf A}_{rms,t}$ ellipsoid satisfies Definition 1.

The basic building block of our ellipsoid matrix consists of the current and past stochastic gradients ${\bf G}_{t}:=[{\bf g}_{1},{\bf g}_{2},\ldots,{\bf g}_{t}].$

We consider ${\bf A}_{rms}$ which is built up as followsThis is a generalization of the diagonal variant proposed by , which preconditions the gradient step by an elementwise division with the square-root of the following estimate $g_{t}=(1-\beta)g_{t-1}+\beta\nabla\mathcal{L}({\bf w}_{t})^{2}$ .

which proves the lower bound for $\zeta=1/\epsilon$ . Now, let us consider the upper end of the spectrum of ${\bf A}_{rms,t}$ . Towards this end, recall the geometric series expansion

and the fact that ${\bf G}{\bf G}^{\top}$ is a sum of exponentially weighted rank-one positive semi-definite matrices of the form ${\bf g}_{i}{\bf g}_{i}^{\intercal}$ . Thus

where the latter inequality holds per assumption for any sample size $|S|$ . Combining these facts we get that

Finally, to achieve uniform equivalence we need the r.h.s. of (24) to be bounded by $1/\epsilon$ . This gives rise to a quadratic equation in $\epsilon$ , namely

which holds for any $t$ and any $\beta\in(0,1)$ as long as

Such an $\epsilon$ always exists but one needs to choose smaller and smaller values as the upper bound on the gradient norm grows. For example, the usual value $\epsilon=10^{-8}$ is valid for all $L_{H}^{2}<9.9\cdot 10^{7}$ . All of the above arguments naturally extend to the diagonal preconditioner $\operatorname{diag}({\bf A}_{rms})$ .

Second, we note that it is no necessary to compute the update step by minimizing Eq. (7) to global optimality. Instead, it suffices to do better than the Cauchy- and Eigenpoint simultaneously . We here adapt this assumption for the case of iteration dependent norms (compare Chapter 6). restate this assumption here

[A.2 restated] Each update step ${\bf s}_{t}$ yields at least as much model decrease as the Cauchy- and Eigenpoint simultaneously, i.e.

where ${\bf u}_{t}$ is an approximation to the corresponding negative curvature direction, i.e., for some $0<\nu<1$ , ${\bf u}_{t}^{\intercal}{\cal H}_{t}{\bf u}_{t}\leq\nu\left(\frac{\|{\bf u}_{t}\|}{\|{\bf u}_{t}\|_{t}}\right)^{2}\lambda_{min}({\mathcal{B}}_{t})\;\text{ and }\;\|{\bf u}_{t}\|_{t}=1.$

In practice, improving upon the Cauchy point is easily satisfied by any Krylov subspace method such as Conjugate Gradients, which ensures convergence to first order critical points. However, while the Steihaug-Toint CG solver can exploit negative curvature, it does not explicitly search for the most curved eigendirection and hence fails to guarantee $m_{t}({\bf s}_{t})\leq m_{t}({\bf s}_{t}^{E})$ . Thus more elaborate Krylov descent methods such as Lanczos method might have to be employed for second-order criticality (See also Appendix F.2 and Chapter 7).

We now restate two results from that precisely quantify the model decrease guaranteed by Assumption 2.

Suppose that ${\bf s}_{t}^{C}$ is computed as in Eq. (28). Then

Suppose that $\lambda_{min}({\mathcal{B}}_{t})<0$ and ${\bf s}_{t}^{E}$ is computed as in Eq. (28). Then

We are now ready to prove the final convergence results. Towards this end, we closely follow the line of arguments developed in . First, we restate the following lemma which holds independent of the trust region constraint choice.

Assume that $\mathcal{L({\bf w})}$ is second-order smooth with Lipschitz constants $L_{g}$ and $L_{H}$ . Furthermore, let Assumption 1 hold. Then

Second, we show that any iterate of Algorithm 1 is eventually successful as long as either the gradient norm or the smallest eigenvalue are above (below) the critical values $\epsilon_{g}$ and $\epsilon_{H}$ .

Assume that $\mathcal{L({\bf w})}$ is second-order smooth with Lipschitz constants $L_{g}$ and $L_{H}$ . Furthermore, let Assumption 1 and 2 hold and suppose that $\|{\bf g}_{t}\|\geq\epsilon_{g}$ as well as

then the step ${\bf s}_{t}$ is successful.

First, by Assumption 2, Lemma 5, $\|{\bf g}_{t}\|\geq\epsilon_{g}$ and Lemma 2, we have

where the last equality uses the above assumed upper bound on $\Delta_{t}$ of Eq. (32). Using this result together with Lemma 7 and the fact that $\|{\bf s}_{t}\|_{2}\leq\mu\|{\bf s}_{t}\|_{t}\leq\mu\Delta_{t}$ due to Lemma 2, we find

where the last inequality makes use of the upper bound assumed on $\delta_{g}$ . Now, we re-use the result of Lemma 10 in , which states that $\left(\frac{\delta_{H}}{\epsilon_{g}}\Delta_{t}+\frac{L_{H}}{\epsilon_{g}}\Delta_{t}^{2}\right)\leq\frac{1-\eta}{2}$ for $\Delta_{t}\leq\min\left\{\sqrt{\frac{(1-\eta)\epsilon_{g}}{12L_{H}}},\frac{(1-\eta)\epsilon_{g}}{3}\right\}$ to conclude that $\left(\frac{\delta_{H}}{\epsilon_{g}}\Delta_{t}+\frac{L_{H}}{\epsilon_{g}}\Delta_{t}^{2}\right)\mu^{4}\leq\frac{1-\eta}{2}$ for our assumed bound on $\Delta_{t}$ in Eq. (32). As a result, Eq. (34) yields

which implies that the iteration $t$ is successful. ∎

First, recall Eq. (31) and note that, since both ${\bf s}_{t}$ and $-{\bf s}_{t}$ are viable search directions, we can assume ${\bf s}_{t}^{\intercal}\nabla F({\bf w}_{t})\leq 0$ w.l.o.g.. Then

Therefore, recalling Eq. (30) as well as the fact that $\frac{\|{\bf u}_{t}\|_{2}}{\|{\bf u}_{t}\|_{t}}\leq\mu$ and $\|{\bf s}_{t}\|_{2}\leq\mu\|{\bf s}_{t}\|_{t}\leq\mu\Delta_{t}$ due to Lemma 2

where the last second inequality is due to the conditions in Eq. (38). Therefore, $\rho_{t}\geq\eta$ and the iteration is successful. ∎

Together, these two results allow us to establish a lower bound on the trust region radius $\Delta_{t}$ .

Assume that $\mathcal{L({\bf w})}$ is second-order smooth with Lipschitz constants $L_{g}$ and $L_{H}$ . Furthermore, let Assumption 1 and 2 hold. Suppose

The proof follows directly from $\Delta_{t}\geq\Delta_{t-1}/\gamma$ as well as the fact that any step is successful as soon as $\Delta_{t}$ falls below $\min\left\{\frac{\epsilon_{g}\mu}{1+L_{g}},\sqrt{\frac{(1-\eta)\epsilon_{g}}{12L_{H}\mu^{8}}},\frac{(1-\eta)\epsilon_{g}}{3\mu^{4}},\frac{(1-\eta)}{2\mu}\frac{\nu\epsilon_{H}}{L_{H}}\right\}$ due to Lemma 8 and 9. ∎

Under the same setting as Lemma 10, the number of successful iterations taken by Algorithm 1 is upper bounded by

where $C:=\eta\min\left\{C_{1},C_{2}\right\}$ , $C_{1}:=\frac{1}{2}\min\left\{\frac{1}{1+L_{g}},C_{g}\right\}$ , $C_{2}:=\frac{\nu\mu^{2}}{2}\min\left\{C_{g}^{2},C_{H}^{2}\right\}$ , $C_{g}:=\min\left\{\frac{\epsilon_{g}\mu}{1+L_{g}},\sqrt{\frac{(1-\eta)\epsilon_{g}}{12L_{H}\mu^{8}}},\frac{(1-\eta)\epsilon_{g}}{3\mu^{4}}\right\}$ , $C_{H}:=\frac{(1-\eta)}{2\mu}\frac{\nu\epsilon_{H}}{L_{H}}$

Suppose Algorithm 1 does not terminate at iteration $t$ . Then either $\|{\bf g}_{t}\|\geq\epsilon_{g}$ or $\lambda_{min}({\bf B})\leq-\epsilon_{h}$ . If $\|{\bf g}_{t}\|\geq\epsilon_{g}$ , according to (29) and Lemma 2, we have

Similarly, in the second case $\lambda_{min}({\bf B}_{t})\leq-\epsilon_{h}$ , from Lemma 2 and 6 we have

Let $T_{\text{succ}}$ denote the number of successful iterations. Since $\mathcal{L}({\bf w})$ is monotonically decreasing, we have

We are now ready to prove the final result. Particularly, given the lower bound on $\Delta_{t}$ established in Lemma 10 we find an upper bound on the number of un-successful iterations, which combined with the result of Lemma 11 on the number of successful iterations yields the total iteration complexity of Algorithm 1.

The result follows by combining the lemmas 10 and 11 as in Theorem 1 of . Specifically, suppose that Algorithm 1 terminates at iteration T. Then the total number of iterations $T=T_{\text{succ}}+T_{\text{unsucc}}$ and $\Delta_{T}=\Delta_{0}\cdot\gamma^{T_{\text{succ}}-T_{\text{unsucc}}}$ . From Lemma 10 we have $\Delta_{T}\geq\frac{1}{\gamma}\min\left\{\frac{\epsilon_{g}\mu}{1+L_{g}},\sqrt{\frac{(1-\eta)\epsilon_{g}}{12L_{H}\mu^{8}}},\frac{(1-\eta)\epsilon_{g}}{3\mu^{4}},\frac{(1-\eta)}{2\mu}\frac{\nu\epsilon_{H}}{L_{H}}\right\}:=\Delta_{inf}$ . Hence, $(T_{\text{succ}}-T_{\text{unsucc}})\log(\gamma)\geq\log(\Delta_{inf}/\Delta_{0})$ , which implies

Finally, combining Eq. 41 with the upper bound on successful steps from Lemma 11 yields

Appendix D Diagonal Dominance in Neural Networks

For random Gaussian Wigner matrix ${\bf W}$ formed as

where $\sim$ stands for i.i.d. draws , the diagonal mass of the expected absolute matrix amounts to

which simplifies to $\frac{1}{d}$ if the diagonal and off-diagonal elements come from the same Gaussian distribution ( $\sigma_{1}=\sigma_{2}$ ). ∎

For the sake of simplicity we only consider Gaussian Wigner matrices but the above argument naturally extends to any distribution with positive expected absolute values, i.e. we only exclude the Dirac delta function as probability density.

D.2 OLS Baseline

As a result, we have that in the limit of large $n$

Appendix B: Background on second-order optimization

The canonical second-order method is Newton’s methods. This algorithm uses the inverse Hessian as a scaling matrix and thus has updates of the form

which is equivalent to optimizing the local quadratic model

to first-order stationarity. Using curvature information to rescale the steepest descent direction gives Newton’s method the useful property of being linearly scale invariant. This gives rise to a problem independent local convergence rate that is super-linear and even quadratic in the case of Lipschitz continuous Hessians (see Theorem 3.5), whereas gradient descent at best achieves linear local convergence .

However, there are certain drawbacks associated with applying classical Newton’s method. First of all, the Hessian matrix may be singular and thus not invertible. Secondly, even if it is invertible the local quadratic model (Eq. 49) that is minimized in each NM iteration may simply be an inadequate approximation of the true objective. As a result, the Newton step is not necessarily a descent step. It may hence approximate arbitrary critical points (including local maxima) or even diverge. Finally, the cost of forming and inverting the Hessian sum up to $O(nd^{2}+d^{3})$ and are thus prohibitively high for applications in large dimensional problems.

Appendix F Trust Region Methods

Trust region methods are among the most principled approaches to overcome the above mentioned issues. These methods also construct a quadratic model $m_{t}$ but constrain the subproblem in such a way that the stepsize is restricted to stay within a certain radius $\Delta_{t}$ within which the model is trusted to be sufficiently adequate

Hence, contrary to line-search methods this approach finds the step ${\bf s}_{t}$ and its length ${}\left\|{\bf s}_{t}\right\|$ simultaneously by optimizing (50). Subsequently the actual decrease $\mathcal{L}({\bf w}_{t})-\mathcal{L}({\bf w}_{t}+{\bf s}_{t})$ is compared to the predicted decrease $m_{t}(0)-m_{t}({\bf s}_{t})$ and the step is only accepted if the ratio $\rho:=\mathcal{L}({\bf w}_{t})-\mathcal{L}({\bf w}_{t}+{\bf s}_{t})/(m_{t}(0)-m_{t}({\bf s}_{t}))$ exceeds some predefined success threshold $\eta_{1}>0$ . Furthermore, the trust region radius is decreased whenever $\rho$ falls below $\eta_{1}$ and it is increased whenever $\rho$ exceeds the ”very successful” threshold $\eta_{2}0$ . Thereby, the algorithm adaptively measures the accuracy of the second-order Taylor model – which may change drastically over the parameter space depending on the behaviour of the higher-order derivativesNote that the second-order Taylor models assume constant curvature. – and adapts the effective length along which the model is trusted accordingly. See for more details.

As a consequence, the plain Newton step ${\bf s}_{N,t}=-\left(\nabla^{2}\mathcal{L}_{t}\right)^{-1}\nabla\mathcal{L}_{t}$ is only taken if it lies within the trust region radius and yields a certain amount of decrease in the objective value. Since many functions look somehow quadratic close to a minimizer the radius can be shown to grow asymptotically under mild assumptions such that eventually full Newton steps are taken in every iteration which retains the local quadratic convergence rate .

F.2 Subproblem solver

Interestingly, there is no need to optimize Eq. (50) to global optimality to retain the remarkable global convergence properties of TR algorithms. Instead, it suffices to do better than the Cauchy- and Eigenpointwhich are the model minimizers along the gradient and the eigendirection associated with its smallest eigenvalue, respectively. simultaneously. One popular approach is to minimize $m_{t}({\bf s})$ in nested Krylov subspaces. These subspaces naturally include the gradient direction as well as increasingly accurate estimates of the leading eigendirection

until (for example) the stopping criterion

is met, which requires increased accuracy as the underlying trust region algorithm approaches criticality. Conjugate gradients and Lanczos method are two iterative routines that implicitly build up a conjugate and orthogonal basis for such a Krylov space respectively and they converge linearly on quadratic objectives with a square-root dependency on the condition number of the Hessian . We here employ the preconditionied Steihaug-Toint CG method in order to cope with possible boundary solutions of (50) but similar techniques exist for the Lanczos solver as well for which we also provide code. As preconditioning matrix for CG we use the same matrix as for the ellipsoidal constraint.

Appendix G Damped (Gauss-)Newton methods

An alternative approach to actively constraining the region within which the model is trusted is to instead penalize the step norm in each iteration in a Lagrangian manner. This is done by so-called damped Newton methods that add a $\lambda>0$ multiple of the identity matrix to the second-order term in the model, which leads to the update step

This can also be solved hessian-free by conjugate gradients (or other Krylov subspace methods). The penalty parameter $\lambda$ is acting inversely to the trust region radius $\Delta$ and it is often updated accordingly. Such algorithms are commonly known as Levenberg-Marquardt algorithms and they were originally tailored towards solving non-linear least squares problems but they have been proposed for neural network training already early on .

Many algorithms in the existing literature replace the use of $\nabla^{2}\mathcal{L}({\bf w}_{t})$ in (53) with the Generalized Gauss Newton matrix or an approximation of the latter . This matrix constitutes the first part of the well-known Gauss-Newton decomposition

Contrary to TR methods, the Levenberg-Marquardt methods never take plain Newton steps since the regularization is always on ( $\lambda>0$ ). Furthermore, if a positive-definite Hessian approximation like the Generalized Gauss Newton matrix is used, this algorithm is not capable of exploiting negative curvature and there are cases in neural network training where the Hessian is much better conditioned than the Gauss-Newton matrix (also see Figure 8). While some scholars believe that positive-definiteness is a desirable feature , we want to point out that following negative curvature directions is necessarily needed to escape saddle points and it can also be meaningful to follow directions of negative eigenvalue $\lambda$ outside a saddle since they guarantee $\mathcal{O}(|\lambda|^{3})$ progress, whereas a gradient descent step yields at least $\|\nabla f({\bf w})\|^{2}$ progress (both under certain stepsize conditions) and one cannot conclude a-priori which one is better in general . Despite these theoretical considerations, many methods based on GGN matrices have been applied to neural network training (see and references therein) and particularly the hessian-free implementations of can be implemented very cheaply .

Appendix H Using Hessian information in Neural Networks

While many theoretical arguments suggest the superiority of regularized Newton methods over gradient based algorithms, several practical considerations cast doubt on this theoretical superiority when it comes to neural network training. Answers to the following questions are particularly unclear: Are saddles even an issue in deep learning? Is superlinear local convergence a desirable feature in machine learning applications (test error)? Are second-order methods more ”vulnerable” to sub-sampling noise? Do worst-case iteration complexities even matter in real-world settings? As a result, the value of Hessian information in neural network training is somewhat unclear a-priori and so far a conclusive empirical study is still missing.

Our empirical findings indicate that the net value of Hessian information for neural network training is indeed somewhat limited for mainly three reasons: 1) second-order methods rarely yield better limit points, which suggests that saddles and spurious local minima are not a major obstacle; 2) gradient methods can indeed run on smaller batch sizes which is beneficial in terms of epoch and when memory is limited; 3) The per-iteration time complexity is noticeably lower for first-order methods. In summary, these observations suggest that advances in hardware and distributed second-order algorithms (e.g., ) will be needed before Newton-type methods can replace (stochastic) gradient methods in deep learning.

Appendix C: Experiment details

To put the previous results into context, we also benchmark several state-of-the-art gradient methods. We fix their sample size to 32 (as advocated e.g. in ) but grid search the stepsize, since it is the ratio of these two quantities that effectively determines the level of stochasticity . The TR methods use a batch size of 128 for the ResNet architecture and 512 otherwiseWe observed weaker performance when running with smaller batches, presumably because second-order methods are likely to ”overfit” noise in small batches in any given iteration as they extract more information of each batch per step by computing curvature.. For a fair comparison, we thus report results in terms of number of backpropagations, epochs and time . The findings are mixed: For small nets such as the MLPs the TR method with RMSProp ellipsoids is superior in all metrics, even when benchmarked in terms of time. However, while Fig. 9 indicates that ellipsoidal TR methods are slightly superior in terms of backpropagations even for bigger nets (ResNets and Autoencoders), a close look at the Figures 10 and 11 (App. C) reveals that they at best manage to keep pace with first-order methods in terms of epochs and are inferior in time. Furthermore, only the autoencoders give rise to a saddle point, which adaptive gradient methods escape faster than vanilla SGD, just like it was the case for second-order methods (see Fig. 3).

Appendix J Default parameters, architectures and datasets

Table 1 reports the default parameters we consider. Only for the larger ResNet18 on CIFAR-10, we adapted the batch size to $128$ due to memory constraints.

We use two real-world datasets for image classification, namely CIFAR-10 and Fashion-MNISTBoth datasets were accessed from https://www.tensorflow.org/api_docs/python/tf/keras/datasets. While Fashion-MNIST consists of greyscale $28\times 28$ images, CIFAR-10 are colored images of size $32\times 32$ . Both datasets have a fixed training-test split consisting of 60,000 and 10,000 images, respectively.

The MLP architectures are simple. For MNIST and Fashion-MNIST we use a $784-128-10$ network with tanh activations and a cross entropy loss. The networks has $101^{\prime}770$ parameters. For the CIFAR-10 MLP we use a $3072-128-128-10$ architecture also with tanh activations and cross entropy loss. This network has $410^{\prime}880$ parameters.

The Fashion-MNIST autoencoder has the same architecture as the one used in . The encoder structure is $784-1000-500-250-30$ and the decoder is mirrored. Sigmoid activations are used in all but the central layer. The reconstructed images are fed pixelwise into a binary cross entropy loss. The network has a total of $2^{\prime}833^{\prime}000$ parameters. The CIFAR-10 autoencoder is taken from the implementation of https://github.com/jellycsc/PyTorch-CIFAR-10-autoencoder.

For the ResNet18, we used the implementation from torchvision for CIFAR-10 as well as a modification of it for Fashion-MNIST that adapts the first convolution to account for the single input channel.

In all of our experiments each method was run on one Tesla P100 GPU using the PyTorch library.