Year: 2019

Abstract:
We investigate the use of ellipsoidal trust region constraints for second-order optimization of neural networks. This approach can be seen as a higher-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 show 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. 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-ofthe-art first-order methods in terms of backpropagations, but further advances in hardware are needed to render Newton methods competitive in terms of time.

heading: 1 INTRODUCTION
text: We consider finite-sum optimization problems of the form
min w∈Rd
[ L(w) :=
n∑ i=1 `(f(w,xi,yi))
] , (1)
which typically arise in neural network training, e.g. for empirical risk minimization over a set of data points (xi,yi) ∈ Rin × Rout, i = 1, . . . , n. Here, ` : Rout × Rout → R+ is a convex loss function and f : Rin × Rd → Rout represents the neural network mapping parameterized by the concatenation of layers w ∈ Rd, which is non-convex due to its multiplicative nature and potentially non-linear activation functions. We assume that L is lower bounded and twice differentiable, i.e. L ∈ C2(Rd,R) and consider finding a first- and second-order stationary point w̄ for which ‖∇L(w̄)‖ ≤ g and λmin ( ∇2L(w̄) ) ≥ − H . Non-convex optimization problems are ubiquitous in machine learning. Among the most prominent examples are present-day deep neural networks, that have achieved outstanding results on core tasks such as collaborative filtering (Wang et al., 2015), sentence classification (Kim, 2014) and image classification (Krizhevsky et al., 2012).
In the era of big data and deep neural networks, stochastic gradient descent (SGD) is one of the most widely used training algorithms (Bottou, 2010). What makes SGD so attractive is its simplicity and per-iteration cost that are 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 (Nesterov, 2013; Shalev-Shwartz et al., 2017) and thus adaptive gradient methods that employ dynamic, coordinate-wise learning rates based on past gradients—including Adagrad (Duchi et al., 2011), RMSprop (Tieleman & Hinton, 2012) and Adam (Kingma & Ba, 2014)—have become a popular alternative, often providing significant speed-ups over SGD. Yet, there exist no theoretical proofs that these methods are faster than gradient descent (Li & Orabona, 2018).
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 (Conn et al., 2000). While second-order algorithms have a long-standing history even in the realm of neural network training (Hagan & Menhaj, 1994; Becker et al., 1988), they were mostly considered as too computationally and memory
expensive for practical applications. Yet, the seminal work of Martens (2010) 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 (Pearlmutter, 1994; Schraudolph, 2002). Among the class of regularized Newton methods, trust region (Conn et al., 2000) and cubic regularization algorithms (Cartis et al., 2011) are the most principled approaches in the sense that they yield the strongest convergence guarantees. Recently, stochastic extensions have emerged (Xu et al., 2017b; Yao et al., 2018; Kohler & Lucchi, 2017; Gratton et al., 2017), 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 (Xu et al., 2017a; Yao et al., 2018). 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. Furthermore, we empirically find that many neural network problems exhibit diagonally dominated Hessian matrices which suggests the effectivity of diagonal preconditioning. (Section 3) • We investigate the first application of ellipsoidal TR methods for deep learning. In Theorem 1
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 2). • Finally, we provide an experimental benchmark across different real-world datasets and
architectures. We also compare against adaptive gradient methods and show results in terms of backprogations, epochs, and wall-clock time; a comparison we were not able to find in the literature. (Section 5)
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.

heading: 2 RELATED WORK
text: First-order methods The prototypical method for optimizing Eq. (1) is SGD (Robbins & Monro, 1951). While the practical success of SGD in non-convex optimization is unquestioned, the theoretical foundation of this phenomenon is still rather limited. Recent findings suggest the ability of this method to escape saddle points and reach local minima in polynomial time for general non-convex problems, but they either need to artificially add noise to the iterates (Ge et al., 2015; Lee et al., 2016) or make an assumption on the inherent noise of vanilla SGD (Daneshmand et al., 2018). For neural network training, a recent line of research proclaims the effectiveness of SGD, but the results usually come at the cost of fairly strong assumptions such as heavy overparametrization and Gaussian inputs (Du et al., 2017; Brutzkus & Globerson, 2017; Li & Yuan, 2017; Du & Lee, 2018; Allen-Zhu et al., 2018). Adaptive gradient methods (Duchi et al., 2011; Tieleman & Hinton, 2012; Kingma & Ba, 2014) build on the intuition that larger learning rates for smaller gradient components and smaller learning rates for larger gradient components balance their respective influences and thereby make the methods behave as if they were optimizing a more isotropic surface. Such approaches have first been suggested for neural networks by LeCun et al. (2012). Recently, convergence guarantees for such methods are starting to appear (Ward et al., 2018; Li & Orabona, 2018). However, these are not superior to the O( −2g ) worst-case complexity of standard gradient descent (Cartis et al., 2012b).
Regularized Newton methods The most principled class of regularized Newton methods are trust region (TR) and adaptive cubic regularization algorithms (ARC) (Conn et al., 2000; Cartis et al., 2011), 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 O( −2g ), ARC only takes at most O( −3/2 g ). However, simple modifications to the TR framework allow these methods to obtain the same accelerated rate (Curtis et al., 2017). Both methods take at most O( −3H ) iterations to find an H approximate second-order stationary point (Cartis et al., 2012a). These rates are optimal for second-order Lipschitz continuous functions (Carmon et al., 2017; Cartis et al., 2012a) and they can be retained even when only sub-sampled gradient and Hessian information is used (Kohler & Lucchi, 2017; Yao et al., 2018; Xu et al., 2017b; Blanchet et al., 2016; Liu et al., 2018; Cartis & Scheinberg, 2017). Furthermore, the involved Hessian information can be computed solely based on Hessian-vector products, which are implementable efficiently for neural networks (Pearlmutter, 1994). This makes these methods particularly attractive for deep learning, but the empirical evidence of their applicability is so far very limited. We are only aware of the works of Liu et al. (2018) and Xu et al. (2017a), which report promising first results but these 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 framework1 (LeCun et al., 2012; Martens, 2010; Martens & Grosse, 2015). These methods have been termed hessian-free since only access to GGN-vector products is required. As the GGN matrix is always positive semidefinite, they cannot leverage negative curvature to escape saddles and hence, there exist no second-order convergence guarantees. Furthermore, there are cases in neural network training where the Hessian is better conditioned than the GGN matrix (Mizutani & Dreyfus, 2008). Nevertheless, the above works report promising preliminary results, most notably Grosse & Martens (2016) report 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 (Ma et al., 2019). Moreover, Xu et al. (2017a) reports many cases where TR and GGN algorithms perform similarly.
This line of work is to be seen as complementary to our approach since it is straight forward to replace the Hessian in the TR framework with the GGN matrix. Furthermore, the preconditioners used in Martens (2010) and Chapelle & Erhan (2011), namely diagonal estimates of the empirical Fisher and Fisher matrix, respectively, can directly be used as matrix norms in our ellipsoidal TR framework.

heading: 3 AN ALTERNATIVE VIEW ON ADAPTIVE GRADIENT METHODS
text: Adaptively preconditioned gradient methods update iterates as wt+1 = wt − ηtA−1/2t gt, where gt is a stochastic estimate of ∇L(wt) and At is a positive definite symmetric pre-conditioning matrix. In Adagrad, Aada,t is the un-centered second moment matrix of the past gradients computed as
Aada,t := GtG ᵀ t + I, (2)
where > 0, I is the d× d identity matrix and Gt = [g1,g2, . . . ,gt]. 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
Arms,t := ( (1− β)Gt diag(βt, . . . , β0)Gᵀt ) + I, (3)
where β ∈ (0, 1). In order to save computational efforts, the diagonal versions diag(Aada) and diag(Arms) 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. In that way, the optimization methods can account for gradients of potentially different scales arising from e.g. different layers of the networks.

heading: 3.1 ADAPTIVE PRECONDITIONING AS ELLIPSOIDAL TRUST REGION
text: Starting from the fact that adaptive methods employ coordinate-wise stepsizes, one can take a principled view of these methods. Namely, their update steps arise from minimizing a first-order
1This algorithm is a simplified TR method, initially tailored for non-linear least squares problems (Nocedal & Wright, 2006)
Taylor model of the function L within an ellipsoidal search space around the current iterate wt, where the diameter of the ellipsoid along a particular coordinate is implicitly given by ηt and ‖gt‖A−1t . Correspondingly, vanilla (S)GD optimizes the same first-order model within a spherical constraint. Fig. 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.2 Since the models are linear, the constrained minimizer is always found on the boundary.
It is well known that GD struggles to progress towards the minimizer of quadratics along lowcurvature directions (see e.g., Goh (2017)). While this effect is negligible for well-conditioned objectives (Fig. 1, left), it leads to drastically 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 w∗ along the low curvature direction which accelerates convergence.
Theorem 1 (Preconditioned gradient methods as TR). A preconditioned gradient step
wt+1 −wt = st := −ηtA−1t gt (4)
with stepsize ηt > 0, symmetric positive definite preconditioner At ∈ Rd×d and gt 6= 0 minimizes a first-order model around wt ∈ Rd in an ellipsoid given by At in the sense that
st := arg min s∈Rd
[ m1t (s) = L(wt) + sᵀgt ] , s.t. ‖s‖At ≤ ηt‖gt‖A−1t . (5)
Corollary 1 (Rmsprop). The step srms,t := −ηtA−1/2rms,tgt minimizes a first-order Taylor model around wt in an ellipsoid given by A 1/2 rms,t (Eq. 3) in the sense that
srms,t := arg min s∈Rd
[ m1t (s) = L(wt) + sᵀgt ] , s.t. ‖s‖
A 1/2 rms,t ≤ ηt‖gt‖A−1/2rms,t . (6)
Equivalent results can be established for Adam using gadam,t := (1−β) ∑t k=0 β
t−kgt as well as for Adagrad by replacing the matrix Aada into the constraint in Eq. (6). 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.
2For illustrative purposes, we only plot every other trust region.

heading: 3.2 DIAGONAL VERSUS FULL PRECONDITIONING
text: A closer look at Fig. 1 reveals that the first two problems come with level sets that are perfectly axis-aligned, which makes these objectives particularly attractive for diagonal preconditioning. For comparison, on the right of Fig. 1, we report another quadratic problem instance, where the Hessian is no longer zero on the off-diagonals. As can be seen, the interaction between coordinates introduces a tilt in the level sets and reduces the superiority of diagonal Adagrad over plain GD. However, using the full preconditioner Aada 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 Agarwal et al. (2018) for a recent exception).
Interestingly, early results by Becker et al. (1988) on the curvature structure of neural nets report a strong diagonal dominance of the Hessian matrix ∇2L(w). This suggests that the loss surface is indeed somewhat axis-aligned. However, the reported numbers are only for tiny feed-forward networks of at most 256 parameters. Therefore, we generalize these findings in the following to larger networks. Furthermore, we contrast the diagonal dominance of real Hessian matrices to the expected behavior of random Wigner matrices. Of course, true Hessians do not have i.i.d. entries but the symmetry of Wigner matrices suggests that this baseline is not completely off. Furthermore, we also compare Hessians of Ordinary Least Squares problems (OLS) with random inputs. For this purpose, let δA define the ratio of diagonal to overall mass of a matrix A, i.e. δA := ∑ i |Ai,i|∑
i ∑ j |Ai,j |
as in (Becker et al., 1988).
Proposition 1 (Diagonal share of Wigner matrix). For random Gaussian3 Wigner matrix W formed as
Wi,j = Wj,i := { ∼ N (0, σ21), i < j ∼ N (0, σ22), i = j,
(7)
where ∼ stands for i.i.d. draws (Wigner, 1993), the diagonal mass of the expected absolute matrix amounts to
δE[|W|] = 1
1 + (d− 1)σ2σ1 . (8)
Thus, if we suppose the Hessian at any given point w were a random Wigner matrix we would expect the share of diagonal mass to fall with O(1/d) as the network grows in size. A similar result can be derived for the large n limit in the case of OLS Hessians.
Proposition 2 (Diagonal share of OLS Hessian). Let X ∈ Rd×n and assume each xi,j is generated i.i.d. with zero-mean finite second moment σ2 > 0. Then the share of diagonal mass of the expected matrix E [|Hols|] amounts to
δE[|Hols|] n→∞→
√ n
√ n+ (d− 1) √ 2 π . (9)
Again, we expect the diagonal mass to fall in d and interestingly, empirical simulations suggest that the result of Proposition 2 holds already in small n settings (see Fig.D.2) and it is likely that finite n results can be derived when adding assumptions such as Gaussian data. As can be seen in Fig. 2 below, even for a practical batch size of n = 32 the diagonal mass δH of neural networks stays above both benchmarks when its input is random as well as with real-world data at random initialization, during training and after convergence.
3The argument naturally extends to any distribution with positive expected absolute values.
These findings are in line with Becker et al. (1988) and suggest that full matrix preconditioning is most probably not worth the additional computational cost for neural networks. Consequently, we use diagonal preconditioning for both first- and second-order methods in all of our experiments in Section 5. Further theoretical elaboration of these findings present an interesting direction of future research.

heading: 4 SECOND-ORDER TRUST REGION METHODS
text: Cubic regularization (Nesterov & Polyak, 2006; Cartis et al., 2011) 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 L around the current iterate wt.4 In particular, in iteration t the update step of the trust region algorithm is computed as
min s∈Rd
[ mt(s) := L(wt) + gᵀt s + 1
2 sᵀBts
] , s.t. ‖s‖At ≤ ∆t (10)
where ∆t > 0 and gt and Bt are either ∇L(wt) and ∇2L(wt) or suitable approximations. The matrix At induces the shape of the constraint set. So far, the common choice for neural networks is At := I, ∀t which gives rise to spherical trust regions (Xu et al., 2017a; Liu et al., 2018). By solving the constrained problem (10), TR methods overcome the problem that pure Newton steps may be ascending, attracted by saddles or not even computable. See Appendix B for more details.

heading: 4.1 CONVERGENCE OF ELLIPSOIDAL TRUST REGION METHODS
text: 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 ‖w‖At := (wᵀAtw) 1/2 induced by symmetric positive definite matrices At are called uniformly equivalent, if ∃µ ≥ 1 such that 1
µ ‖w‖At ≤ ‖w‖2 ≤ µ‖w‖At , ∀w ∈ Rd,∀t = 1, 2, . . . . (11)
We now establish a result which shows that the RMSProp ellipsoid is indeed uniformly equivalent.
4In the following we only treat TR methods, but we would like to emphasize that the use of matrix induced norms can directly be transferred to the cubic regularization framework.
Lemma 1 (Uniform equivalence). Suppose ‖gt‖2 ≤ L2H for all wt ∈ Rd, t = 1, 2, . . . Then there always exists > 0 such that the proposed preconditioning matrices Arms,t (Eq. 3) are uniformly equivalent, i.e. Def. 1 holds. The same holds for the diagonal variant.
Consequently, the ellipsoids Arms,t can directly be applied to any convergent TR framework without losing convergence guarantees (Conn et al. (2000), Theorem 6.6.8).5 Interestingly, this result cannot be established for Aada,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 (Tieleman & Hinton, 2012) and Adadelta (Zeiler, 2012).
Why ellipsoids? There are many sources for ill-conditioning in neural networks such as un-centered and correlated inputs (LeCun et al., 2012), saturated hidden units, and different weight scales in different layers (Van Der Smagt & Hirzinger, 1998). While the quadratic term of model (10) 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 Chapter 7.7 in Conn et al. (2000)). The ellipsoids we propose are such that they allow for longer steps along coordinates that have seen small gradient components in past and vice versa. Thereby the TR shape is adaptively adjusted to fit the current region of the non-convex loss landscape. This procedure is not only effective when the iterates are in an ill-conditioned neighborhood of a minimizer (Figure 1), but it also helps to escape elongated plateaus (see autoencoder in Section 5).

heading: 4.2 A STOCHASTIC TR FRAMEWORK FOR NEURAL NETWORK TRAINING
text: 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. In particular, we adapt the framework of Yao et al. (2018); Xu et al. (2017b) to the case of iteration-dependent norm constraints (Algorithm 1). We prove that the O ( max { −2g −1 H , −3 H }) worst-case iteration complexity is retained, which is optimal if the function is second-order smooth (Cartis et al., 2012a). Towards this end, we assume that the derivative estimates are sufficiently accurate in the following sense.
Assumption 1 (Sufficiently accurate derivatives). The approximations of the gradient and Hessian at step t satisfy
‖gt −∇L(wt)‖ ≤ δg and ‖Bt −∇2L(wt)‖ ≤ δH ,
where δg ≤ (1−η) g4 and δH ≤ min { (1−η)v H 2 , 1 } , 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 (Xu et al., 2017b; Kohler & Lucchi, 2017; Tripuraneni et al., 2017). Following these references, we assume access to the function value in each iteration for our theoretical analysis but we note that convergence can be retained even for fully stochastic trust region methods (Gratton et al., 2017; Chen et al., 2018; Blanchet et al., 2016) and indeed our experiments in Section 5 use sub-sampled function values due to memory constraints.
5Note that the assumption of bounded batch gradients, i.e. smooth objectives, is common in the analysis of stochastic algorithms (Allen-Zhu, 2017; Defazio et al., 2014; Schmidt et al., 2017; Duchi et al., 2011).
Algorithm 1 Stochastic Ellipsoidal Trust Region Method 1: Input: w0 ∈ Rd, γ > 1, 1 > η > 0, ∆0 > 0 2: for t = 0, 1, . . . , until convergence do 3: Compute approximations gt and Bt. If ‖gt‖ ≤ g , set gt := 0. 4: Set At := Arms,t or At := diag (Arms,t) (see Eq. (3)). 5: Obtain st by solving mt(st) approximately. 6: Compute ratio of function over model decrease
ρt = L(wt)− L(wt + st) mt(0)−mt(st)
(12)
7: Set
∆t+1 = { γ∆t if ρS,t > η ∆t/γ if ρS,t < η , wt+1 = { wt + st if ρt ≥ η (successful) wt otherwise (unsuccessful)
8: end for
Theorem 2 (Convergence rate of Algorithm 1). Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 and 2 hold. Then Algorithm 1 finds an O( g, H) first- and second-order stationary point in at most O ( max { −2g −1 H , −3 H }) iterations.

heading: 5 EXPERIMENTS
text: Trust region methods 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. As can be seen in Figure 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 TR method. Interestingly, this makes an actual difference in the image reconstruction quality of autoencoders (see Figure 12). We thus draw the clear conclusion that the ellipsoidal trust region constraints we propose are to be preferred over their spherical counterpart when training neural networks. Both the experimental and architectural details are provided in Appendix C.
Benchmark with SGD To put the previous results into context, we also benchmark several stateof-the-art gradient methods. We fix their sample size to 32 (as advocated e.g. in Masters & Luschi (2018)) but grid search the stepsize since it is the ratio of these two quantities that effectively determines the level of stochasticity (Jastrzebski et al., 2017). As the TR methods have a larger batch size6 of 128–512, we report results both in terms of number of backpropagations and epochs for a fair comparison. A close look at Figure 4 and 10 (Appendix) indicates that the ellipsoidal TR methods can be slightly superior in terms of backprops but at best manage to keep pace with first-order methods in terms of epochs. Furthermore, the limit points of both first- and second-order methods yield the same order of loss in most experiments. When taking gradient norms into account (plot omitted), we indeed find no spurious local minima and only the autoencoders give rise to saddle points.

heading: 6 CONCLUSION
text: 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. Bekas et al. (2007)) 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., Osawa et al. (2018); Dünner et al. (2018)) will be needed before Newton-type methods can replace (stochastic) gradient methods in deep learning.

heading: Appendix A: Proofs

heading: A NOTATION
text: 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 ‖ · ‖ we denote an arbitrary norm. For a symmetric positive definite matrix A we introduce the compact notation ‖w‖A = (wᵀAw)1/2, where w ∈ Rd.

heading: B EQUIVALENCE OF PRECONDITIONED GRADIENT DESCENT AND FIRST-ORDER TRUST REGION METHODS
text: Theorem 3 (Theorem 1 restated). A preconditioned gradient step
wt+1 −wt = st := −ηtA−1t gt (13)
with stepsize ηt > 0, symmetric positive definite preconditioner At ∈ Rd×d and gt 6= 0 minimizes a first-order local model around wt ∈ R in an ellipsoid given by At in the sense that
st := arg min s∈Rd
[ m1t (s) = L(wt) + sᵀgt ] , s.t. ‖s‖At ≤ ηt‖gt‖A−1t . (14)
Proof. We start the proof by noting that the optimization problem in Eq. (14) is convex. For ηt > 0 the constraint satisfies the Slater condition since 0 is a strictly feasible point. As a result, any KKT point is a feasible minimizer and vice versa.
Let L(s, λ) denote the Lagrange dual of Eq. (5) L(s, λ) := L(wt) + sᵀgt + λ ( ‖s‖A − ηt‖gt‖A−1t ) . (15)
Any point s is a KKT point if and only if the following system of equations is satisfied
∇sL(s, λ) = gt + λ
‖s‖At Ats = 0 (16) λ ( ‖s‖At − ηt‖gt‖A−1t ) = 0. (17)
‖s‖At − ηt‖gt‖A−1t ≤ 0 (18) λ ≥ 0. (19)
For st as given in Eq. (4) we have that ‖st‖At = √ η2t gt(A −1 t ) ᵀ AtA −1 t gt = ηt √ gtA −1 t gt = ηt‖gt‖A−1t . (20)
and thus 17 and 18 hold with equality such that any λ ≥ 0 is feasible. Furthermore,
∇sL(st, λ) = ∇f(wt) + λ
‖st‖At Atst
(4) = gt − ηt
λ
ηt‖gt‖A−1t AtA
−1 t gt = gt −
λ
‖gt‖A−1t gt
(21) is zero for λ = ‖gt‖A−1 ≥ 0. As a result, st 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 methods7 on an ill-conditioned quadratic problem. While the method 1st TR
7Essentially Algorithm 1 with mt based on a first order Taylor expansion, i.e. m1t (s) as in Eq. (14).
optimizes a linear model within a ball in each iteration, 1st TRada optimizes the same model over the ellipsoid given by the Adagrad matrix Aada. The results in Figure 5 show that the methods behave very similar to their constant stepsize analogues.

heading: C CONVERGENCE OF ELLIPSOIDAL TR METHODS
text: 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 µ ≥ 1 such that
1 µ ‖w‖At ≤ ‖w‖2 ≤ µ‖w‖At , ∀t,∀w ∈ Rd. (22)
Towards this end, Conn et al. (2000) identify the following sufficient condition on the basis of which we will prove that our proposed ellipsoid Arms is indeed uniformly equivalent under some mild assumptions. Lemma 2 (Theorem 6.7.1 in Conn et al. (2000)). Suppose that there exists a constant ζ ≥ 1 such that
1 ζ ≤ σmin (At) ≤ σmax (At) ≤ ζ ∀t, (23)
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 (Conn et al., 2000)). 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 Yao et al. (2018) (see Theorem 2 below).
First, we show that the proposed Arms,t ellipsoid satisfies Definition 1.
Lemma 3 (Lemma 1). Suppose ‖gt‖2 ≤ L2H for all wt ∈ Rd, t = 1, 2, . . . Then there always exists > 0 such that the proposed preconditioning matrices Arms,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 Gt := [g1,g2, . . . ,gt].
We consider Arms which is built up as follows8
8This is a generalization of the diagonal variant proposed by Tieleman & Hinton (2012), which precondition the gradient step by an elementwise division with the square-root of the following estimate gt = (1− β)gt−1 + β∇L(wt)2.
Arms,t := (1− β)Gdiag(βt, βt−1, . . . , β0)︸ ︷︷ ︸ :=D Gᵀ + I. (24) From the construction of Arms,t it directly follows that for any unit length vector u ∈ Rd \ {0}, ‖u‖2 = 1 we have
uᵀ ((1− β)GDGᵀ + I)u =(1− β)uᵀGD1/2(D1/2)ᵀGᵀu + ‖u‖22 =(1− β) ( (D1/2)ᵀGᵀu )ᵀ ( (D1/2)ᵀGᵀu ) + ‖u‖22
≥ > 0,
(25)
which proves the lower bound for ζ = 1/ . Now, let us consider the upper end of the spectrum of Arms,t. Towards this end, recall the geometric series expansion
t∑ i=0 βt−i = t∑ i=0 βi = 1− βt+1 1− β (26)
and the fact that GG> is a sum of exponentially weighted rank-one positive semi-definite matrices of the form gig ᵀ i . Thus
λmax(gig ᵀ i ) = Tr(gig ᵀ i ) = ‖∇gi‖ 2 ≤ L2H ,
where the latter inequality holds per assumption for any sample size |S|. Combining these facts we get that
uᵀ ((1− β)GDGᵀ + I)u =(1− β)uᵀGDGᵀu + ‖u‖22
=(1− β) t∑ i=0 βt−1uᵀgig ᵀ i u + ‖u‖ 2 2
≤(1− β) t∑ i=0 βt−iL2H‖u‖22 + ‖u‖22
=(1− βt+1)L2H + .
(27)
As a result we have that
≤ λmin (Arms,t) ≤ λmax (Arms,t) ≤ ( 1− βt+1 ) L2H + (28)
Finally, to achieve uniform equivalence we need the r.h.s. of (28) to be bounded by 1/ . This gives rise to a quadratic equation in , namely
2 + ( 1− βt+1 ) L2H − 1 ≤ 0 (29)
which holds for any t and any β ∈ (0, 1) as long as
0 ≤ ≤ 1 2
( √ L4H + 4− L 2 H). (30)
Such an 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 = 10−8 is valid for all L2H < 9.9 · 107. All of the above arguments naturally extend to the diagonal preconditioner diag(Arms).
Second, we note that it is no necessary to compute the update step by minimizing Eq. (10) to global optimality. Instead, it suffices to do better than the Cauchy- and Eigenpoint simultaneously (Conn et al., 2000; Yao et al., 2018). We here adapt this assumption for the case of iteration dependent norms (compare (Conn et al., 2000) Chapter 6). restate this assumption here
Assumption 2 (Approximate model minimization). Each update step st yields at least as much model decrease as the Cauchy- and Eigenpoint simultaneously, i.e.
mt(st) ≤ mt(sCt ) and mt(st) ≤ mt(sEt ), (31) where
sCt := arg min 0≤α≤∆t mt(−α gt ‖gt‖t ) and sEt := arg min |α|≤∆t mt(αut), (32)
where ut is an approximation to the corresponding negative curvature direction, i.e., for some 0 < ν < 1, uᵀtHtut ≤ ν ( ‖ut‖ ‖ut‖t )2 λmin(Bt) and ‖ut‖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 mt(st) ≤ mt(sEt ). Thus more elaborate Krylov descent methods such as Lanczos method might have to be employed for second-order criticality (See also Appendix B.2 and Conn et al. (2000) Chapter 7).
We now restate two results from Conn et al. (2000) that precisely quantify the model decrease guaranteed by Assumption 2.
Lemma 4 (Model decrease: Cauchy Point (Theorem 6.3.1. in Conn et al. (2000))). Suppose that sCt is computed as in Eq. (32). Then
mt(0)−mt(sCt ) ≥ 1
2 ‖gt‖min{ ‖gt‖ 1 + ‖Bt‖ ,∆t ‖gt‖ ‖gt‖t }. (33)
Lemma 5 (Model decrease: Eigenpoint (Theorem 6.6.1 in Conn et al. (2000))). Suppose that λmin(Bt) < 0 and sEt is computed as in Eq. (32). Then
mt(0)−mt(sEt ) ≥ − 1
2 νλmin(Bt) ( ‖ut‖ ‖ut‖t )2 ∆2t (34)
We are now ready to prove the final convergence results. Towards this end, we closely follow the line of arguments developed in Yao et al. (2018). First, we restate the following lemma which holds independent of the trust region constraint choice.
Lemma 6 (Yao et al. (2018)). Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 hold. Then
F (xt + st)− F (xt)−mt(xt) ≤ sᵀt (∇F (xt)− gt) + 1
2 δH‖st‖2 +
1 2 LH‖st‖3. (35)
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 g and H .
Lemma 7 (Eventually successful iteration - ‖gt‖ ≥ g). Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 and 2 hold and suppose that ‖gt‖ ≥ g as well as
δg < 1− η 4µ2 g, ∆t ≤ min  µ g1 + Lg , √ (1− η) g 12LH 1 µ4 , (1− η) g 3µ4  , (36) then the step st is successful.
Proof. First, by Assumption 2, Lemma 4, ‖gt‖ ≥ g and Lemma 1, we have
−mt(st) ≥ 1
2 ‖gt‖min{ ‖gt‖ 1 + ‖Bt‖ ,∆t ‖gt‖ ‖gt‖t }
≥ 1 2 ‖gt‖min{ g 1 + ‖Bt‖ ,∆t ‖gt‖ ‖gt‖t }
≥ 1 2 ‖gt‖min{ g 1 + ‖Bt‖ , ∆t µ }
= 1
2 g ∆t µ ,
(37)
where the last equality uses the above assumed upper bound on ∆t of Eq. (36). Using this result together with Lemma 6 and the fact that ‖st‖2 ≤ µ‖st‖t ≤ µ∆t due to Lemma 1, we find
1− ρt = L(wt + st)− L(wt)−mt(st)
−mt(st) (38)
≤ δg∆tµ+
1 2δH∆ 2 tµ 2 + 12LH∆ 3 tµ 3
1 2 g ∆t µ
(39)
= 2 δg g µ2 + δh g ∆tµ 3 + LH g ∆2tµ 4 (40)
≤ 1− η 2 + ( δH g ∆t + LH g ∆2t ) µ4, (41)
where the last inequality makes use of the upper bound assumed on δg. Now, we re-use the result of Lemma 10 in Yao et al. (2018), which states that ( δH g ∆t + LH g ∆2t ) ≤ 1−η2 for
∆t ≤ min {√
(1−η) g 12LH , (1−η) g 3
} to conclude that ( δH g ∆t + LH g ∆2t ) µ4 ≤ 1−η2 for our assumed
bound on ∆t in Eq. (36). As a result, Eq. (38) yields
1− ρt ≤ 1− η,
which implies that the iteration t is successful.
Lemma 8 (Eventually successful iteration - λmin(Bt) ≤ − H ). Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 and 2 hold and suppose that ‖gt‖ < g and λmin(Bt) < − H . If
δH < 1− η
2 ν H ,∆t ≤ (1− η) 2µ ν H LH
(42)
then iteration t is successful.
Proof. First, recall Eq. (35) and note that, since both st and −st are viable search directions, we can assume sᵀt∇F (wt) ≤ 0 w.l.o.g.. Then
L(wt + st)− L(wt)−mt(wt) ≤ 1
2 δH‖st‖2 +
1 2 LH‖st‖3
Therefore, recalling Eq. (34) as well as the fact that ‖ut‖2‖ut‖t ≤ µ and ‖st‖2 ≤ µ‖st‖t ≤ µ∆t due to Lemma 1
1− ρt = L(wt + st)− L(wt)−mt(st)
−mt(st)
≤ 1 2δH‖st‖ 2 + 12LH‖st‖ 3
ν 2 |λmin(Bt)|∆ 2 tµ 2
≤ 1 2δH‖st‖ 2 + 12LH‖st‖ 3
ν 2 H∆ 2 tµ 2
≤ 1 2δH∆ 2 tµ 2 + 12LH∆ 3 tµ 3
ν 2 H∆ 2 tµ 2
= δH ν H + LH∆tµ
ν H < 1− η,
(43)
where the last second inequality is due to the conditions in Eq. (42). Therefore, ρt ≥ η and the iteration is successful.
Together, these two results allow us to establish a lower bound on the trust region radius ∆t. Lemma 9. Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 and 2 hold. Suppose
δg < 1− η
4 g, δH < min{ 1− η 2 ν H , 1}.
then for Algorithm 1 we have
∆t ≥ 1
γ min
{ gµ
1 + Lg , √ (1− η) g 12LHµ8 , (1− η) g 3µ4 , (1− η) 2µ ν H LH } , ∀t = 1, 2, . . . (44)
Proof. The proof follows directly from ∆t ≥ ∆t−1/γ as well as the fact that any step is successful as soon as ∆t falls below min {
gµ 1+Lg
, √
(1−η) g 12LHµ8 , (1−η) g 3µ4 , (1−η) 2µ ν H LH
} due to Lemma 7 and 8.
Lemma 10 (Number of successful iterations). Under the same setting as Lemma 9, the number of successful iterations taken by Algorithm 1 is upper bounded by
|Tsucc| ≤ L(x0)− L(x∗) C H min{ 2g, 2H} ,
where C := ηmin {C1, C2}, C1 := 12 min { 1 1+Lg , Cg } , C2 := νµ 2 2 min { C2g , C 2 H } ,Cg :=
min
{ gµ 1+Lg , √ (1−η) g 12LHµ8 , (1−η) g 3µ4 } , CH := (1−η) 2µ ν H LH
Proof. Suppose Algorithm 1 does not terminate at iteration t. Then either ‖gt‖ ≥ g or λmin(B) ≤ − h. If ‖gt‖ ≥ g , according to (33) and Lemma 1, we have
−mt(st) ≥ 1
2 ‖gt‖min{ ‖gt‖ 1 + ‖Ht‖ ,∆t 1 µ }
≥ 1 2 g min{ g 1 + Lg , Cg g, CH H}
≥ C1 g min{ g, H}.
Similarly, in the second case λmin(Bt) ≤ − h, from Lemma 1 and 5 we have
−mt(st) ≥ 1
2 ν|λmin(Bt)|∆2tµ2 ≥ C2 H min{ 2g, 2H}.
Let Tsucc denote the number of successful iterations. Since L(w) is monotonically decreasing, we have
L(w0)− L(w∗) ≥ ∑ t=0 L(wt)− L(wt+1)
≥ ∑ t∈Tsucc L(wt)− L(wt+1)
≥ ∑ t∈Tsucc −mt(st)η
≥ ∑ t∈Tsucc C H min{ 2g, 2H}
≥ |Tsucc|C H min{ 2g, 2H}, which proves the assertion.
We are now ready to prove the final result. Particularly, given the lower bound on ∆t established in Lemma 9 we find an upper bound on the number of un-successful iterations, which combined with the result of Lemma 10 on the number of successful iterations yields the total iteration complexity of Algorithm 1.
Theorem 4 (Theorem 2 restated). Assume that L(w) is second-order smooth with Lipschitz constants Lg and LH . Furthermore, let Assumption 1 and 2 hold. Then Algorithm 1 finds an O( g, H) first- and second-order stationary point in at mostO ( max { −2g −1 H , −3 H }) iterations.
Proof. The result follows by combining the lemmas 9 and 10 as in Theorem 1 of Xu et al. (2017a).

heading: D DIAGONAL DOMINANCE IN NEURAL NETWORKS
text: In the following, we make statements about the diagonal share of random matrices. As E[ 1x ] 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 W formed as
Wi,j = Wj,i := { ∼ N (0, σ21), i < j ∼ N (0, σ22), i = j,
(45)
where ∼ stands for i.i.d. draws (Wigner, 1993), the diagonal mass of the expected absolute matrix amounts to
δE[|W|] = 1
1 + (d− 1)σ2σ1 . (46)
Proof.
δE,W = ∑d k=1 E [|Wk,k|]∑d
k=1 ∑d l=1 E [|Wk,l|] = dE [|W1,1|] dE [|W1,1|] + d(d− 1)E [|W1,2|]
= dσ1 √ 2/π
dσ1 √ 2/π + d(d− 1)σ2 √ 2/π =
1
1 + d(d−1)σ2
√ 2/π
dσ1 √ 2/π
= 1
1 + (d− 1)σ2σ1
(47)
which simplifies to 1d if the diagonal and off-diagonal elements come from the same Gaussian distribution (σ1 = σ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 Hols := XᵀX, where X ∈ Rd×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 D.2) and it is likely that finite n results can be derived when adding assumptions such as Gaussian data.
Proposition 4 (Proposition 2 restated). Let X ∈ Rn×d and assume each xi,k is generated i.i.d. with zero-mean finite second moment σ2 > 0. Then the share of diagonal mass of the expected matrix E [|Hols|] amounts to
δE[|Hols|] n→∞→
√ n
√ n+ (d− 1) √ 2 π
(48)
Proof.
δE[|Hols|] = ∑d k=1 E [|(Hols)k,k|]∑d
k=1 ∑d l=1 E [|(Hols)k,l|] =
∑d k=1 E [ | ∑n i=1 x 2 i,k| ]
∑d k=1 ∑d l=1 E [| ∑n i=1 xi,kxi,l|]
= d ∑n i=1 E [ x2i,1 ]
d ∑n i=1 E [ x2i,1 ] + d(d− 1)E [| ∑n i=1 xi,1xi,2|]
(49)
Where we used the fact that all xi,k are i.i.d. variables. Per assumption, we have E [ x2i,1 ]
= σ2,∀i. Furthermore, the products xi,1xi,2 are i.i.d with expectation 0 and variance σ4. By the central limit theorem
ZN = 1√ n n∑ i=1 xi,1xi,2 → Z
in law, with Z ∼ N (0, σ4), since E [ Z2n ] = 1nE [ ( ∑n i=1 xi,1xi,2) 2 ] = 1n ∑n i=1 E[x2i,1]E[x2i,2] = σ4
due to the independence assumption. Then E ( |ZN | 1|ZN |≥R ) ≤ E(|ZN |2/R) ≤ σ4/R. This implies that
E (|ZN |)→ E(|Z|) = √ 2
π σ2
As a result, we have that in the limit of large n
δE[|Hols|] n→∞→ dnσ
2 dnσ2 + d(d− 1) √ n √
2 πσ
2 =
1
1 + (d−1) √ 2 π√
n
=
√ n
√ n+ (d− 1) √ 2 π
(50)

heading: Appendix B: Background on second-order optimization

heading: A NEWTON’S METHOD
text: 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
wt+1 = wt −∇2L(wt)−1∇L(wt), (51)
which is equivalent to optimizing the local quadratic model
mN (wt) := L(wt) +∇L(wt)ᵀs + 1
2 sᵀ∇2L(wt)s (52)
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 Nocedal & Wright (2006) Theorem 3.5), whereas gradient descent at best achieves linear local convergence (Nesterov, 2013).
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. 52) 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(nd2 + d3) and are thus prohibitively high for applications in large dimensional problems.

heading: B TRUST REGION METHODS
text: B.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 mt but constrain the subproblem in such a way that the stepsize is restricted to stay within a certain radius ∆t within which the model is trusted to be sufficiently adequate
min s∈Rd
mt(s) = L(wt) +∇L(wt)ᵀs + 1
2 sᵀ∇2L(wt)s, s.t. ‖s‖ ≤ ∆t. (53)
Hence, contrary to line-search methods this approach finds the step st and its length ‖st‖ simultaneously by optimizing (53). Subsequently the actual decrease L(wt) − L(wt + st) is compared to the predicted decrease mt(0) − mt(st) and the step is only accepted if the ratio ρ := L(wt) − L(wt + st)/(mt(0) −mt(st)) exceeds some predefined success threshold η1 > 0. Furthermore, the trust region radius is decreased whenever ρ falls below η1 and it is increased whenever ρ exceeds the "very successful" threshold η20. 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 derivatives9 – and adapts the effective length along which the model is trusted accordingly. See Conn et al. (2000) for more details.
As a consequence, the plain Newton step sN,t = − ( ∇2Lt )−1∇Lt 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 (Conn et al., 2000).
B.2 SUBPROBLEM SOLVER
Interestingly, there is no need to optimize Eq. (53) to global optimality to retain the remarkable global convergence properties of TR algorithms. Instead, it suffices to do better than the Cauchy- and Eigenpoint10 simultaneously. One popular approach is to minimize mt(s) in nested Krylov subspaces. These subspaces naturally include the gradient direction as well as increasingly accurate estimates of the leading eigendirection
span{gt,Btgt,B2tgt, . . . ,B j tgt} (54)
until (for example) the stopping criterion
‖∇mt(sj)‖ ≤ ‖∇L(wt)‖min{κK , ‖∇L(wt)‖θ}, κK < 1, θ ≥ 0 (55) 9Note that the second-order Taylor models assume constant curvature.
10which are the model minimizers along the gradient and the eigendirection associated with its smallest eigenvalue, respectively.
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 (Conn et al., 2000). We here employ the preconditionied Steihaug-Toint CG method (Steihaug, 1983) in order to cope with possible boundary solutions of (53) 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.

heading: C DAMPED (GAUSS-)NEWTON METHODS
text: 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 λ > 0 multiple of the identity matrix to the second-order term in the model, which leads to the update step
min s∈Rd
mt(s) = L(wt) +∇L(wt)ᵀs + 1
2 sᵀ(∇2L(wt) + λI)s
= L(wt) +∇L(wt)ᵀs + 1
2 sᵀ∇2L(wt)s + λ‖s‖2.
(56)
This can also be solved hessian-free by conjugate gradients (or other Krylov subspace methods). The penalty parameter λ is acting inversely to the trust region radius ∆ 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 (Nocedal & Wright, 2006) but they have been proposed for neural network training already early on (Hagan & Menhaj, 1994).
Many algorithms in the existing literature replace the use of ∇2L(wt) in (56) with the Generalized Gauss Newton matrix (Martens, 2010; Chapelle & Erhan, 2011) or an approximation of the latter (Martens & Grosse, 2015). This matrix constitutes the first part of the well-known Gauss-Newton decomposition
∇2L(·) = 1 n n∑ i=1
`′′(fi(·))∇fi(·)∇fi(·)ᵀ︸ ︷︷ ︸ :=AGGN
+ 1
n n∑ i=1 `′(fi(·))∇2fi(·), (57)
where l′ and l′′ are the first and second derivative of l : Rout → R+ assuming that out = 1 (binary classification and regression task) for simplicity here.
It is interesting to note that the GGN matrix AGGN of neural networks is equivalent to the Fisher matrix used in natural gradient descent (Amari, 1998) 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 (Pascanu & Bengio, 2013). As can be seen in (57) the matrix AGGN 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 57). 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 `′(fi(·)) terms in 57 go to zero for all i ∈ {1, . . . , 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 (Bottou et al., 2018).Finally, Dauphin et al. (2014) propose the use of the absolute Hessian instead of the GGN matrix in a framework similar to 56. This method has been termed saddle-free Newton even though its manifold of attraction to a given saddle is non-empty11.
11It 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 (Lee et al., 2016).
C.1 COMPARISON TO TRUST REGION
Contrary to TR methods, the Levenberg-Marquardt methods never take plain Newton steps since the regularization is always on (λ > 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 (Mizutani & Dreyfus, 2008) (also see Figure 9). While some scholars believe that positive-definiteness is a desirable feature (Martens, 2010; Chapelle & Erhan, 2011), 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 λ outside a saddle since they guarantee O(|λ|3) progress, whereas a gradient descent step yields at least ‖∇f(w)‖2 progress (both under certain stepsize conditions) and one cannot conclude a-priori which one is better in general (Curtis & Robinson, 2017; Alain et al., 2018). Despite these theoretical considerations, many methods based on GGN matrices have been applied to neural network training (see Martens (2014) and references therein) and particularly the hessian-free implementations of (Martens, 2010; Chapelle & Erhan, 2011) can be implemented very cheaply (Schraudolph, 2002).

heading: D USING HESSIAN INFORMATION IN NEURAL NETWORKS
text: 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., Osawa et al. (2018); Dünner et al. (2018)) will be needed before Newton-type methods can replace (stochastic) gradient methods in deep learning.

heading: Appendix C: Experiment details

heading: A EXPERIMENTAL RESULTS
text: A.1 ELLIPSOIDAL TRUST REGION VS. FIRST-ORDER OPTIMIZER

heading: B FURTHER EXPERIMENT DETAILS
text: B.1 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 FashionMNIST12. While Fashion-MNIST consists of greyscale 28 × 28 images, CIFAR-10 are colored images of size 32× 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′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′880 parameters.
The Fashion-MNIST autoencoder has the same architecture as the one used in Hinton & Salakhutdinov (2006); Xu et al. (2017a); Martens (2010); Martens & Grosse (2015). 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′833′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 (Paszke et al., 2017) library.
B.2 RECONSTRUCTED IMAGES FROM AUTOENCODERS
12Both datasets were accessed from https://www.tensorflow.org/api_docs/python/tf/keras/datasets

Review summaries:
In this paper, the authors investigate the use of ellipsoidal trust region constraints for second order optimization. The authors first show that adaptive gradient methods can be viewed as first-order trust region methods with ellipsoid constraints. The authors then show that the preconditioning matrix of RMSProp and Adam can be used as norm inducing matrices for second order trust region methods. This ellipsoidal trust region method is empirically compared with first order gradient methods and spherical second order trust region methods.