Adaptive norms for deep learning with regularized Newton methods

TL;DR

This paper introduces an adaptive norm approach for second-order optimization in neural networks, showing theoretical convergence guarantees and empirical advantages over spherical constraints, with potential for future hardware improvements.

Contribution

It establishes a connection between adaptive gradient methods and trust region methods, providing convergence proofs and demonstrating empirical benefits of ellipsoidal constraints.

Findings

Ellipsoidal constraints outperform spherical ones in backpropagation count and loss value.

Preconditioning matrices in Adam and RMSProp meet convergence conditions for trust region methods.

Newton methods show promise but need hardware advances for competitive computational time.

Abstract

We investigate the use of regularized Newton methods with adaptive norms for optimizing neural networks. This approach can be seen as a second-order counterpart of adaptive gradient methods, which we here show to be interpretable as first-order trust region methods with ellipsoidal constraints. In particular, we prove that the preconditioning matrix used in RMSProp and Adam satisfies the necessary conditions for provable convergence of second-order trust region methods with standard worst-case complexities on general non-convex objectives. Furthermore, we run experiments across different neural architectures and datasets to find that the ellipsoidal constraints constantly outperform their spherical counterpart both in terms of number of backpropagations and asymptotic loss value. Finally, we find comparable performance to state-of-the-art first-order methods in terms of backpropagations, but further advances in hardware are needed to render Newton methods competitive in terms of computational time.

Tables

Table 1. Table 1: Default parameters

	$\left|{𝒮}_0\right|$	${\Delta }_0$	${\Delta }_{\text{max}}$	${\eta }_1$	${\eta }_2$	${\gamma }_1$	${\gamma }_2$	${\kappa }_K$ (krylov tol.)
TR${}_{\text{uni}}$	512	${10}^{-4}$	10	${10}^{-4}$	0.95	1.1	1.5	0.1
TR${}_{\text{ada}}$	512	${10}^{-4}$	10	${10}^{-4}$	0.95	1.1	1.5	0.1
TR${}_{\text{rms}}$	512	${10}^{-4}$	10	${10}^{-4}$	0.95	1.1	1.75	0.1

Equations

$$\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }\left[\mathcal{L}(\mathbf{w}):=\sum _{i=1}^n\ell (f(\mathbf{w},{\mathbf{x}}_i,{\mathbf{y}}_i))\right],$$

$${\mathbf{A}}_{ada,t}:={\mathbf{G}}_t{\mathbf{G}}_t^{⊺}+ϵ\mathbf{I},$$

$${\mathbf{A}}_{rms,t}:=\left((1-\beta ){\mathbf{G}}_t\mathrm{diag}({\beta }^t,\dots ,{\beta }^0){\mathbf{G}}_t^{⊺}\right)+ϵ\mathbf{I},$$

$${\mathbf{w}}_{t+1}-{\mathbf{w}}_t={\mathbf{s}}_t:=-{\eta }_t{\mathbf{A}}_t^{-1}{\mathbf{g}}_t$$

$${\mathbf{s}}_t:=\arg \underset{\mathbf{s}\in {\mathbb{R}}^d}{\min }\left[m_t^1(\mathbf{s})=\mathcal{L}({\mathbf{w}}_t)+{\mathbf{s}}^{⊺}{\mathbf{g}}_t\right],\text{s.t.}\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_t}\le {\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_t^{-1}}.$$

$${\mathbf{s}}_{rms,t}:=\arg \underset{\mathbf{s}\in {\mathbb{R}}^d}{\min }\left[m_t^1(\mathbf{s})=\mathcal{L}({\mathbf{w}}_t)+{\mathbf{s}}^{⊺}{\mathbf{g}}_t\right],\text{s.t.}\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_{rms,t}^{1/2}}\le {\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_{rms,t}^{-1/2}}.$$

$$\underset{\mathbf{s}\in {\mathbb{R}}^d}{\min }\left[m_t(\mathbf{s}):=\mathcal{L}({\mathbf{w}}_t)+{\mathbf{g}}_t^{⊺}\mathbf{s}+\frac{1}{2}{\mathbf{s}}^{⊺}{\mathbf{B}}_t\mathbf{s}\right],\text{s.t.}\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_t}\le {\Delta }_t,$$

$$\frac{1}{\mu }\Vert \mathbf{w}{\Vert }_{{\mathbf{A}}_t}\le \Vert \mathbf{w}{\Vert }_2\le \mu \Vert \mathbf{w}{\Vert }_{{\mathbf{A}}_t}.$$

$$\Vert {\mathbf{g}}_t-\nabla \mathcal{L}({\mathbf{w}}_t)\Vert \le {\delta }_g\text{ and }\Vert {\mathbf{B}}_t-{\nabla }^2\mathcal{L}({\mathbf{w}}_t)\Vert \le {\delta }_H,$$

$${\Delta }_{t+1}$$

$${\mathbf{w}}_{t+1}-{\mathbf{w}}_t={\mathbf{s}}_t:=-{\eta }_t{\mathbf{A}}_t^{-1}{\mathbf{g}}_t$$

$${\mathbf{s}}_t:=$$

$$L(\mathbf{s},\lambda ):=\mathcal{L}({\mathbf{w}}_t)+{\mathbf{s}}^{⊺}{\mathbf{g}}_t+\lambda \left(\Vert \mathbf{s}{\Vert }_{\mathbf{A}}-{\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_t^{-1}}\right).$$

$${\nabla }_{\mathbf{s}}L(\mathbf{s},\lambda )={\mathbf{g}}_t+\frac{\lambda }{\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_t}}{\mathbf{A}}_t\mathbf{s}=0$$

$$\lambda \left(\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_t}-{\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_t^{-1}}\right)=0.$$

$$\Vert \mathbf{s}{\Vert }_{{\mathbf{A}}_t}-{\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_t^{-1}}\le 0$$

$$\lambda \ge 0.$$

$$\Vert {\mathbf{s}}_t{\Vert }_{{\mathbf{A}}_t}=\sqrt{{\eta }_t^2{\mathbf{g}}_t{({\mathbf{A}}_t^{-1})}^{⊺}{\mathbf{A}}_t{\mathbf{A}}_t^{-1}{\mathbf{g}}_t}={\eta }_t\sqrt{{\mathbf{g}}_t{\mathbf{A}}_t^{-1}{\mathbf{g}}_t}={\eta }_t\Vert {\mathbf{g}}_t{\Vert }_{{\mathbf{A}}_t^{-1}}.$$

$${\nabla }_{\mathbf{s}}L({\mathbf{s}}_t,\lambda )$$

$$\frac{1}{\mu }\Vert \mathbf{w}{\Vert }_{{\mathbf{A}}_t}\le \Vert \mathbf{w}{\Vert }_2\le \mu \Vert \mathbf{w}{\Vert }_{{\mathbf{A}}_t},\forall t,\forall \mathbf{w}\in {\mathbb{R}}^d.$$

$$\frac{1}{\zeta }\le {\sigma }_{\min }\left({\mathbf{A}}_t\right)\le {\sigma }_{\max }\left({\mathbf{A}}_t\right)\le \zeta \forall t,$$

$${\mathbf{A}}_{rms,t}:=\left((1-\beta )\mathbf{G}\underset{:=\mathbf{D}}{\underbrace{\mathrm{diag}({\beta }^t,{\beta }^{t-1},\dots ,{\beta }^0)}}{\mathbf{G}}^{⊺}\right)+ϵ\mathbf{I}.$$

$${\mathbf{u}}^{⊺}\left((1-\beta )\mathbf{G}\mathbf{D}{\mathbf{G}}^{⊺}+ϵ\mathbf{I}\right)\mathbf{u}$$

$$=$$

$$\ge$$

$$\sum _{i=0}^t{\beta }^{t-i}=\sum _{i=0}^t{\beta }^i=\frac{1-{\beta }^{t+1}}{1-\beta }$$

$${\lambda }_{max}({\mathbf{g}}_i{\mathbf{g}}_i^{⊺})=\mathrm{Tr}({\mathbf{g}}_i{\mathbf{g}}_i^{⊺})=\Vert \nabla {\mathbf{g}}_i{\Vert }^2\le L_H^2,$$

$${\mathbf{u}}^{⊺}\left((1-\beta )\mathbf{G}\mathbf{D}{\mathbf{G}}^{⊺}+ϵ\mathbf{I}\right)\mathbf{u}$$

$$=$$

$$\le$$

$$=$$

$$ϵ\le {\lambda }_{min}\left({\mathbf{A}}_{rms,t}\right)\le {\lambda }_{max}\left({\mathbf{A}}_{rms,t}\right)\le \left(1-{\beta }^{t+1}\right)L_H^2+ϵ$$

$${ϵ}^2+\left(1-{\beta }^{t+1}\right)L_H^2ϵ-1\le 0$$

Full text

1 Introduction

We consider finite-sum optimization problems of the form

[TABLE]

which typically arise in neural network training, e.g. for empirical risk minimization over a set of data points $({\bf x}_{i},{\bf y}_{i})\in\mathbb{R}^{in}\times\mathbb{R}^{out},i=1,\ldots,n$. Here, $\ell:\mathbb{R}^{out}\times\mathbb{R}^{out}\rightarrow\mathbb{R}^{+}$ is a convex loss function and $f:\mathbb{R}^{in}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{out}$ represents the neural network mapping parameterized by the concatenation of the weight layers ${\bf w}\in\mathbb{R}^{d}$, which is non-convex due to its multiplicative nature and potentially non-linear activation functions. We assume that $\mathcal{L}$ is lower bounded and twice differentiable, i.e. $\mathcal{L}\in C^{2}(\mathbb{R}^{d},\mathbb{R})$ and consider finding a first- and second-order stationary point $\bar{{\bf w}}$ for which $\|\nabla\mathcal{L}(\bar{{\bf w}})\|\leq\epsilon_{g}$ and $\lambda_{\min}\left(\nabla^{2}\mathcal{L}(\bar{{\bf w}})\right)\geq-\epsilon_{H}$.

In the era of deep neural networks, stochastic gradient descent (SGD) is one of the most widely used training algorithms [9]. 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 [49, 59] and thus adaptive gradient methods that employ dynamic, coordinate-wise learning rates based on past gradients—including Adagrad [27], RMSprop [61] and Adam [36]—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 [19]. While second-order algorithms have a long-standing history even in the realm of neural network training [33, 6], they were mostly considered as too computationally and memory expensive for practical applications. Yet, the seminal work of [44] 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 [55, 58]. Among the class of regularized Newton methods, trust region [19] and cubic regularization algorithms [14] are the most principled approaches in the sense that they yield the strongest convergence guarantees. Recently, stochastic extensions have emerged  [67, 68, 37, 31], 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 [66, 68]. 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.

2 Related work

First-order methods

The prototypical method for optimizing Eq. (1) is SGD [56]. 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 [29, 39] or make an assumption on the inherent noise of SGD [22]. 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 [25, 11, 41, 26, 4]. Adaptive gradient methods [27, 61, 36] 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 [38] and convergence guarantees are starting to appear [64, 40]. However, these are not superior to the $\mathcal{O}(\epsilon_{g}^{-2})$ worst-case complexity of standard gradient descent [15].

Regularized Newton methods

The most principled class of regularized Newton methods are trust region (TR) and adaptive cubic regularization algorithms (ARC) [19, 14], 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 [21]. Both methods take at most $\mathcal{O}(\epsilon_{H}^{-3})$ iterations to find an $\epsilon_{H}$ approximate second-order stationary point [13]. These rates are optimal for second-order Lipschitz continuous functions [12, 13] and they can be retained even when only sub-sampled gradient and Hessian information is used [37, 68, 67, 8, 42, 16]. Furthermore, the involved Hessian information can be computed solely based on Hessian-vector products, which are implementable efficiently for neural networks [55]. 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 [42] and [66], which report promising first results but are by no means fully encompassing.

Gauss-Newton methods

An interesting line of research proposes to replace the Hessian by (approximations of) the generalized-Gauss-Newton matrix (GGN) within a Levenberg-Marquardt framework111This algorithm is a simplified TR method, initially tailored for non-linear least squares problems [51] [38, 44, 46]. 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 [48]. Nevertheless, the above works report promising preliminary results, most notably [32] 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 [43]. Moreover, [66] 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 [44] and [17], namely diagonal estimates of the empirical Fisher and Fisher matrix, respectively, can directly be used as matrix norms in our ellipsoidal TR framework.

3 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

[TABLE]

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

[TABLE]

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.

3.1 Adaptive preconditioning as ellipsoidal Trust Region

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.222We 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., [30]). 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.

Lemma 1** (Preconditioned gradient methods as TR).**

A preconditioned gradient step

${\bf w}_{t+1}-{\bf w}_{t}={\bf s}_{t}:=-\eta_{t}{\bf A}_{t}^{-1}{\bf g}_{t}$

(4)

with stepsize $\eta_{t}>0$, symmetric positive definite preconditioner ${\bf A}_{t}\in\mathbb{R}^{d\times d}$ and ${\bf g}_{t}\not=0$ minimizes a first-order model around ${\bf w}_{t}\in\mathbb{R}^{d}$ in an ellipsoid given by ${\bf A}_{t}$ in the sense that

$\displaystyle{\bf s}_{t}:=\arg\min_{{\bf s}\in\mathbb{R}^{d}}\left[m_{t}^{1}({\bf s})=\mathcal{L}({\bf w}_{t})+{\bf s}^{\intercal}{\bf g}_{t}\right],\>\text{s.t.}\quad\|{\bf s}\|_{{\bf A}_{t}}\leq\eta_{t}\|{\bf g}_{t}\|_{{\bf A}_{t}^{-1}}.$

(5)

Corollary 1** (Rmsprop).**

The step ${\bf s}_{rms,t}:=-\eta_{t}{\bf A}_{rms,t}^{-1/2}{\bf g}_{t}$ minimizes a first-order Taylor model around ${\bf w}_{t}$ in an ellipsoid given by ${\bf A}_{rms,t}^{1/2}$ (Eq. 3) in the sense that

[TABLE]

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.

3.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 [1] for an exception). Interestingly, early results by [6] 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.333Of 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 [6].

Proposition 1** (Diagonal share of Wigner matrix).**

For a random Gaussian444The argument naturally extends to any distribution with positive expected absolute values. Wigner matrix ${\bf W}$ (see Eq. (42)) the diagonal mass of the expected absolute matrix amounts to: $\delta_{\mathbb{E}\left[|{\bf W}|\right]}=\frac{1}{1+(d-1)\frac{\sigma_{2}}{\sigma_{1}}}.$

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.

Proposition 2** (Diagonal share of OLS Hessian).**

Let ${\bf X}\in\mathbb{R}^{d\times n}$ and assume each ${\bf x}_{i,j}$ is generated i.i.d. with zero-mean and finite second moment $\sigma^{2}>0$. Then the share of diagonal mass of the expected matrix $\mathbb{E}\left[|{\bf H}_{\text{ols}}|\right]$ amounts to: $\delta_{\mathbb{E}\left[|{\bf H}_{\text{ols}}|\right]}\overset{n\rightarrow\infty}{\rightarrow}\frac{\sqrt{n}}{\sqrt{n}+(d-1)\sqrt{\frac{2}{\pi}}}.$

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 [6] 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.

4 Second-order Trust Region Methods

Cubic regularization [50, 14] 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}$.555In 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

[TABLE]

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 [66, 42]. 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.

Why ellipsoids?

There are many sources for ill-conditioning in neural networks such as un-centered and correlated inputs [38], saturated hidden units, and different weight scales in different layers [63]. 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 [19]). 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.

4.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.

Definition 1** (Uniformly equivalent norms).**

The norms $\|{\bf w}\|_{{\bf A}_{t}}:=\left({\bf w}^{\intercal}{\bf A}_{t}{\bf w}\right)^{1/2}$ induced by symmetric positive definite matrices ${\bf A}_{t}$ are called uniformly equivalent, if $\exists\mu\geq 1$ such that $\forall{\bf w}\in\mathbb{R}^{d},\forall t=1,2,\ldots$

[TABLE]

We now establish a result which shows that the RMSProp ellipsoid is indeed uniformly equivalent.

Lemma 2** (Uniform equivalence).**

Suppose $\|{\bf g}_{t}\|^{2}\leq L_{H}^{2}$ for all ${\bf w}_{t}\in\mathbb{R}^{d},$ $t=1,2,\ldots$ Then there always exists $\epsilon>0$ such that the proposed preconditioning matrices ${\bf A}_{rms,t}$ (Eq. 3) are uniformly equivalent, i.e. Def. 1 holds. The same holds for the diagonal variant.

Consequently, the ellipsoids ${\bf A}_{rms,t}$ can directly be applied to any convergent TR framework without losing the guarantee of convergence ([19], Theorem 6.6.8).666Note that the assumption of bounded batch gradients, i.e. smooth objectives, is common in the analysis of stochastic algorithms [3, 24, 57, 27]. 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 [61] and Adadelta [69].

4.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.

Assumption 1** (Sufficiently accurate derivatives).**

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

[TABLE]

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 [67, 37, 62]. 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 [31, 18, 8] and indeed our experiments in Section 5 use sub-sampled function values due to memory constraints. Secondly, we adapt the framework of [68, 67], which allows for cheap inexact subproblem minimization, to the case of iteration-dependent constraint norms (Alg. 1).

Assumption 2** (Approximate model minimization).**

*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 [68].

Theorem 1** (Convergence rate of Algorithm 1).**

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.

5 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 [60] 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.

6 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. [7]) 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.,  [52, 28]) 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 [32] 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.

Appendix A: Proofs

Appendix A Notation

Throughout this work, scalars are denoted by regular lower case letters, vectors by bold lower case letters and matrices as well as tensors by bold upper case letters. By $\|\cdot\|$ we denote an arbitrary norm. For a symmetric positive definite matrix ${\bf A}$ we introduce the compact notation $\|{\bf w}\|_{\bf A}=\left({\bf w}^{\intercal}{\bf A}{\bf w}\right)^{1/2}$, where ${\bf w}\in\mathbb{R}^{d}$.

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

Theorem 2** (Theorem 1 restated).**

A preconditioned gradient step

${\bf w}_{t+1}-{\bf w}_{t}={\bf s}_{t}:=-\eta_{t}{\bf A}_{t}^{-1}{\bf g}_{t}$

(9)

with stepsize $\eta_{t}>0$, symmetric positive definite preconditioner ${\bf A}_{t}\in\mathbb{R}^{d\times d}$ and ${\bf g}_{t}\not=0$ minimizes a first-order local model around ${\bf w}_{t}\in\mathbb{R}$ in an ellipsoid given by ${\bf A}_{t}$ in the sense that

$\displaystyle{\bf s}_{t}:=$ $\displaystyle\arg\min_{{\bf s}\in\mathbb{R}^{d}}\left[m_{t}^{1}({\bf s})=\mathcal{L}({\bf w}_{t})+{\bf s}^{\intercal}{\bf g}_{t}\right],\quad\text{s.t.}\;\|{\bf s}\|_{{\bf A}_{t}}\leq\eta_{t}\|{\bf g}_{t}\|_{{\bf A}_{t}^{-1}}.$

(10)

Proof.

We start the proof by noting that the optimization problem in Eq. (10) is convex. For $\eta_{t}>0$ the constraint satisfies the Slater condition since [math] is a strictly feasible point. As a result, any KKT point is a feasible minimizer and vice versa.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 methods777Essentially 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

C.1 Proof sketch

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. [68]). 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.

Definition 2** (Definition 1 restated).**

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

[TABLE]

Towards this end, [19] 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.

Lemma 3** (Theorem 6.7.1 in [19]).**

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

[TABLE]

then Definition 1 holds.

Having uniformly equivalent norms is sufficient to prove convergence of ellipsoidal TR methods (se AN.1 and Theorem 6.6.8 in [19]). 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 [68] (see Theorem 1 below).

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

Lemma 4** (Lemma 2 restated).**

Suppose $\|{\bf g}_{t}\|^{2}\leq L_{H}^{2}$ for all ${\bf w}_{t}\in\mathbb{R}^{d},$ $t=1,2,\ldots$ Then there always exists $\epsilon>0$ such that the proposed preconditioning matrices ${\bf A}_{rms,t}$ (Eq. 3) are uniformly equivalent, i.e. Def. 1 holds. The same holds for the diagonal variant.

Proof.

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 follows888This is a generalization of the diagonal variant proposed by [61], 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}$.

[TABLE]

From the construction of ${\bf A}_{rms,t}$ it directly follows that for any unit length vector ${\bf u}\in\mathbb{R}^{d}\setminus\{0\},\|{\bf u}\|_{2}=1$ we have

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

As a result we have that

[TABLE]

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

[TABLE]

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

[TABLE]

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 [19, 68]. We here adapt this assumption for the case of iteration dependent norms (compare [19] Chapter 6). restate this assumption here

Assumption 3** (Approximate model minimization).**

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

[TABLE]

where

[TABLE]

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 [19] Chapter 7).

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

Lemma 5** (Model decrease: Cauchy Point (Theorem 6.3.1. in [19])).**

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

[TABLE]

Lemma 6** (Model decrease: Eigenpoint (Theorem 6.6.1 in [19])).**

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

[TABLE]

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

Lemma 7** ([68]).**

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

[TABLE]

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}$ .

Lemma 8** (Eventually successful iteration - $\|{\bf g}_{t}\|\geq\epsilon_{g}$).**

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

[TABLE]

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

Proof.

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

[TABLE]

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

[TABLE]

where the last inequality makes use of the upper bound assumed on $\delta_{g}$. Now, we re-use the result of Lemma 10 in [68], 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

[TABLE]

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

Lemma 9** (Eventually successful iteration - $\lambda_{min}({\bf B}_{t})\leq-\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}\|<\epsilon_{g}$ and $\lambda_{min}({\bf B}_{t})<-\epsilon_{H}$. If

[TABLE]

then iteration $t$ is successful.

Proof.

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

[TABLE]

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

[TABLE]

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}$.

Lemma 10**.**

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

[TABLE]

then for Algorithm 1 we have

[TABLE]

Proof.

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. ∎

Lemma 11** (Number of successful iterations).**

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

[TABLE]

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}}$

Proof.

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

[TABLE]

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

[TABLE]

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

[TABLE]

which proves the assertion. ∎

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.

Theorem 3** (Theorem 1 restated).**

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.

Proof.

The result follows by combining the lemmas 10 and 11 as in Theorem 1 of [66]. 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

[TABLE]

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

[TABLE]

∎

Appendix D Diagonal Dominance in Neural Networks

In the following, we make statements about the diagonal share of random matrices. As $\mathbb{E}[\frac{1}{x}]$ might not exist for a random variable $x$, we cannot compute the expectation of the diagonal share but rather of for computing the diagonal share of the expectation of the random matrix in absolute terms. Note that this notion is still meaningful, as the average of many non-diagonally dominated matrices with positive entries cannot become diagonally dominated.

D.1 Proof of Proposition 1

Proposition 3** (Proposition 1 restated).**

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

[TABLE]

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

[TABLE]

Proof.

[TABLE]

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

When considering regression tasks, a direct competitor to neural network models is the classical Ordinary Least Squares (OLS) regression, which minimizes a quadratic loss over a linear model. In this case the Hessian simply amounts to the input-covariance matrix ${\bf H}_{\text{ols}}:={\bf X}^{\intercal}{\bf X}$, where ${\bf X}\in\mathbb{R}^{d\times n}$. We here show that the diagonal share of the expected matrix itself also decays in $d$, when $n$ grows to infinity. However, empirical simulations suggest the validity of this result even for much smaller values of $n$ (see Figure 6) and it is likely that finite $n$ results can be derived when adding assumptions such as Gaussian data.

Proposition 4** (Proposition 2 restated).**

Let ${\bf X}\in\mathbb{R}^{n\times d}$ and assume each ${\bf x}_{i,k}$ is generated i.i.d. with zero-mean finite second moment $\sigma^{2}>0$. Then the share of diagonal mass of the expected matrix $\mathbb{E}\left[|{\bf H}_{\text{ols}}|\right]$ amounts to

[TABLE]

Proof.

[TABLE]

Where we used the fact that all ${\bf x}_{i,k}$ are i.i.d. variables. Per assumption, we have $\mathbb{E}\left[{\bf x}_{i,1}^{2}\right]=\sigma^{2},\forall i$. Furthermore, the products ${\bf x}_{i,1}{\bf x}_{i,2}$ are i.i.d with expectation 0 and variance $\sigma^{4}$. By the central limit theorem

[TABLE]

in law, with $Z\sim{\cal N}(0,\sigma^{4})$, since $\mathbb{E}\left[Z_{n}^{2}\right]=\frac{1}{n}\mathbb{E}\left[\left(\sum_{i=1}^{n}{\bf x}_{i,1}{\bf x}_{i,2}\right)^{2}\right]=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[{\bf x}_{i,1}^{2}]\mathbb{E}[{\bf x}_{i,2}^{2}]=\sigma^{4}$ due to the independence assumption. Then $E\left(\;|Z_{N}|\;1_{|Z_{N}|\geq R}\right)\ \leq\ E(|Z_{N}|^{2}/R)\leq\sigma^{4}/R$. This implies that

[TABLE]

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

[TABLE]

∎

Appendix B: Background on second-order optimization

Appendix E Newton’s Method

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

[TABLE]

which is equivalent to optimizing the local quadratic model

[TABLE]

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 [51] Theorem 3.5), whereas gradient descent at best achieves linear local convergence [49].

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

F.1 Outer iterations

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

[TABLE]

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 derivatives999Note that the second-order Taylor models assume constant curvature. – and adapts the effective length along which the model is trusted accordingly. See [19] 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 [19].

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 Eigenpoint101010which 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

[TABLE]

until (for example) the stopping criterion

[TABLE]

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 [19]. We here employ the preconditionied Steihaug-Toint CG method [60] 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

[TABLE]

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 [51] but they have been proposed for neural network training already early on [33].

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 [44, 17] or an approximation of the latter [46]. This matrix constitutes the first part of the well-known Gauss-Newton decomposition

[TABLE]

where $l^{\prime}$ and $l^{\prime\prime}$ are the first and second derivative of $l:\mathbb{R}^{out}\rightarrow\mathbb{R}^{+}$ assuming that $out=1$ (binary classification and regression task) for simplicity here.

It is interesting to note that the GGN matrix ${\bf A}_{GGN}$ of neural networks is equivalent to the Fisher matrix used in natural gradient descent [5] in many cases like linear activation function and squared error as well as sigmoid and cross-entropy or softmax and negative log-likelihood for which the extended Gauss-Newton is defined [53]. As can be seen in (54) the matrix ${\bf A}_{GGN}$ is positive semidefinite (and low rank if $n<d$). As a result, there exist no second-order convergence guarantees for such methods on general non-convex problems. On the other end of the spectrum, the GGN also drops possibly positive terms from the Hessian (see 54). Hence it is not guaranteed to be an upper bound on the latter in the PSD sense. Essentially, GGN approximations assume that the network is piece-wise linear and thus the GGN and Hessian matrices only coincide in the case of linear and ReLU activations or non-curved loss functions. For any other activation the GGN matrix may approximate the Hessian only asymptotically and if the $\ell^{\prime}(f_{i}(\cdot))$ terms in 54 go to zero for all $i\in\{1,\ldots,n\}$. In non-linear least squares such problems are called zero-residual problems and GN methods can be shown to have quadratic local convergence there. In any other case the convergence rate does not exceed the linear local convergence bound of gradient descent. In practice however there are cases where deep neural nets do show negative curvature in the neighborhood of a minimizer [10].Finally, [23] propose the use of the absolute Hessian instead of the GGN matrix in a framework similar to 53. This method has been termed saddle-free Newton even though its manifold of attraction to a given saddle is non-empty111111It is the same as that for GD, which renders the method unable to escape e.g. when initialized right on a saddle point. To be fair, the manifold of attraction for GD constitutes a measure zero set [39]..

G.1 Comparison to trust region

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 [48] (also see Figure 8). While some scholars believe that positive-definiteness is a desirable feature [44, 17], 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 [20, 2]. Despite these theoretical considerations, many methods based on GGN matrices have been applied to neural network training (see [45] and references therein) and particularly the hessian-free implementations of [44, 17] can be implemented very cheaply [58].

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.,  [52, 28]) will be needed before Newton-type methods can replace (stochastic) gradient methods in deep learning.

Appendix C: Experiment details

Appendix I Ellipsoidal Trust Region vs. First-order Optimizers

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 [47]) but grid search the stepsize, since it is the ratio of these two quantities that effectively determines the level of stochasticity [35]. The TR methods use a batch size of 128 for the ResNet architecture and 512 otherwise121212We 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

Parameters

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.

Datasets

We use two real-world datasets for image classification, namely CIFAR-10 and Fashion-MNIST131313Both 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.

Network architectures

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 [34, 66, 44, 46]. 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 [54] library.

Appendix K Reconstructed Images from Autoencoders