Second-Order Optimization for Non-Convex Machine Learning: An Empirical Study

Introduction

The large-scale nature of many modern ML problems poses computational challenges which have rendered many classical optimization methods (developed for scientific computing and other related areas) inefficient or inapplicable. In this light, first-order optimization methods, such as SGD and its variants, have been the workhorse for ML applications due to their simplicity and versatility. On the other hand, since they consider only first-order gradient information, these methods come with well known deficiencies. These include relatively-slow convergence and sensitivity to hyper-parameter settings such as learning rate. Furthermore, without dedicated “baby-sitting” of these methods, they can stagnate at high training loss and have difficulty in escaping saddle points and flat regions. These deficiencies are particularly problematic in highly non-convex ML problems such as those that arise in neural network applications. By incorporating second-order information, i.e., curvature, second-order optimization methods hold the promise to solve these well-known deficiencies. When implemented naïvely, however, second-order methods are clearly not computationally competitive with first order alternatives. This, in turn, has unfortunately lead to the conventional wisdom that second-order methods are not appropriate for large-scale ML applications.

Within the scientific computing community, more so than the ML community, it is well-known that not only stochastic Newton-type methods in general, and Gauss-Newton in particular, can be made scalable , but more importantly, and unlike first-order methods, they are also very resilient to a variety of adversarial effects . Perhaps the most well known example is their resilience to ill-conditioning. A subtle, yet potentially more severe, example is that the success of most first-order methods is tightly intertwined with fine-tunning (often many) hyper-parameters, most importantly the step size . It is rare that these methods exhibit acceptable performance on first try, and it often takes many trials and errors before one can see reasonable results. In fact, the “true training time”, which almost always includes the time it takes to appropriately tune these parameters, can be frustratingly long. This sort of brute force hyper-parameter tuning is naturally computationally as well as financially expensive. In contrast, second-order optimization algorithms involve much less parameter tuning, and are less sensitive to specific choices of their hyper-parameters . A third example that has received attention due to the popularity of non-convex deep learning problems has to do with avoiding (possibly degenerate) saddle points and finding a local minimum. While some first-order algorithms can guarantee convergence to an approximate second-order critical point where gradient is very small and Hessian is almost positive semi-definite., e.g., , the vast majority of them lack such performance guarantees. Instead, their convergence can, at best, only be ensured to first-order critical points, which include saddle points. It has been argued that converging to saddle points can be undesirable for obtaining low generalization errors . In addition, important cases have been demonstrated where SGD stagnates at such high training error points . In contrast, employing the curvature information in the form of Hessian, in addition to the advantages mentioned above, can help with achieving convergence to second-order criticality.

Here, we aim to provide an empirical evaluation of variants of Newton-type methods for non-convex ML problems, and study whether they can address, in a computationally competitive manner, the main aforementioned challenges associated with first order methods, i.e., relatively-slow convergence, sensitivity to hyper-parameter settings, stagnation, and entrapment near saddle points. To do so, we exploit recent theoretical developments in , and focus on variants of trust-region and (adaptive) cubic regularization algorithms, in which the Hessian is suitably approximated. More specifically, in the context of several non-convex machine learning applications, we study the empirical performance of sub-sampled versions of these algorithms and set out to paint a more complete picture of their practical impact. In the process, we highlight the above-mentioned shortcomings of first-order methods and, in their light, assess the effectiveness of employing curvature information as remedy.

Main Questions

To accomplish our objective, we set out to answer the following questions, all of which are designed to shed light on various practical aspects of this paper’s underlying thesis.

(Computational Efficiency) Can these sub-sampled Newton-type methods be computationally efficient enough to be competitive with hand-tuned SGD with momentum?

(Robustness to Hyper-parameters) Does the performance of such Newton-type methods exhibit robustness to hyper-parameter tuning?

(Escaping Saddle Point) Does employing Hessian information help with avoiding saddle points and converging to lower training errors in highly non-convex problems?

(Generalization Performance) Can second-order methods be beneficial for obtaining low generalization error in machine learning problems?

(Benefits of Sub-sampling) How does sub-sampling in general, and various sub-sampling schemes in particular, affect the performance of TR and ARC methods?

(Comparison Among Second-Order Methods) How do sub-sampled TR and ARC compare with other second-order methods like (sub-sampled) Gauss-Newton (GN) and limited-memory BFGS (L-BFGS) methods?

Methodology

To minimize the effect of many confounding factors involved in almost all empirical evaluations, rather than chasing to beat the state-of-the-art performances on benchmark tasks, we focus our attention on two simple, yet illustrative, classes of non-convex ML problems, i.e., multi-layer perceptron (MLP) networks and non-linear least squares (NLS). More specifically, to address Q.1, Q.2, Q.3, and Q.4, which are concerned with comparisons among first and second-order methods, we consider training (deep) MLP networks with no additional bells and whistles; and to address Q.5 and Q.6, which involve comparisons among various second-order methods, we consider a simpler NLS problem. These questions, references to their theoretical studies, and sections of this paper involving their empirical treatment with relevant figures, are all gathered in Table 1.

In Section 3.2, we address Q.1–Q.4. In particular, in the context of image classification (Section 3.2.1) and deep auto-encoder (Section 3.2.2), we study the efficiency of sub-sampled TR method, which incorporates inexactness in both Hessian and the sub-problem solver, as compared with hand-tuned SGD with momentum. We answer item Q.1 by measuring the convergence speed over total computational cost. We then treat Q.2 by demonstrating the resiliency/sensitivity of various algorithms with respect to their main hyper-parameter. We do this through multiple simulations of all these examples with several choices of the main hyper-parameter for each algorithm. The main hyper-parameter is understood as the one for which, in practice, there is no “typical” value. For SGD with momentum, the learning rate is considered as the main parameter, since the momentum parameter is typically set to $\approx 0.9$. For trust region, the main hyper-parameter is the initial trust region, as there are typical values for other parameters of the algorithm.For Q.3, we consider various initialization schemes, including those that are close to high-level saddle points. We study the behavior of SGD-based methods near such regions and evaluate whether an appropriate use of curvature can indeed help with escaping such undesirable saddle points and making continued progress towards areas with lower training error. To address Q.4, we present test performances for all our experiments and assess the generalization errors obtained by all the methods considered.

In Section 3.3, we turn our attention to Q.5–Q.6, where we consider the NLS problem arising from binary classification task with least-squares loss. On several real datasets, we then demonstrate the effect of various sub-sampling strategies on the performance of sub-sampled TR and ARC methods, as compared with other second-order algorithms.

Background

In this section, we give a brief review of the general formulation of the optimization problems considered in our study as well as the Newton-type algorithms in question.

Following many machine learning applications, we consider the “finite-sum” optimization problem

2 Main Algorithms: Sub-Sampled TR and ARC

Arguably, line-search is the most straightforward approach for globalization of many Newton-type algorithms. However, near saddle points where the gradient magnitude can be small, traditional line search methods can be very ineffective and in fact produce iterates that can get stuck at a saddle point . Trust region and cubic regularization methods are two elegant globalization alternatives that, specially recently, have attracted much attention.

In large-scale non-convex settings, where the application of exact Hessian can be computationally infeasible, recently, theoretically studied the variants of TR and ARC algorithms in which the Hessian is suitably approximated. The details of the resulting TR and ARC algorithms are give in Algorithms 1 and 2, respectively. Iterations of these algorithms involve the following sub-problems:

Hessian Sub-Sampling We consider (P1) in large-scale regime where $n,d\gg 1$. In such settings, the mere evaluations of the Hessian and the gradient increase linearly in $n$. As studied in , given a sampling distribution $\{p_{i}\}_{i=1}^{n}$ over the set of indices $\{1,2,\cdots,n\}$, the sub-sampled Hession has the form

where $\mathcal{S}\subset\{1,2,\cdots,n\}$ is a random sample collection. When $\left|\mathcal{S}\right|\ll n$, sub-sampling can offer significant computational savings; see for examples of studies in convex settings. Recently theoretically showed that randomized sub-sampling can also be seamlessly extended to non-convex settings.

further showed that, in certain settings, one can construct more “informative” distributions over the indices in $\{1,2,\ldots,n\}$, as opposed to oblivious uniform sampling. Indeed, it is typically advantageous to bias the probability distribution towards picking indices corresponding to those $f_{i}$’s which are more relevant, in certain sense, in forming the Hessian. One such setting where this is possible is the finite-sum optimization of the form,

Inexact Sub-problem Solver In Algorithms 1 and 2, it is imperative that the sub-problems (1a) and (1b), respectively, are solved only approximately. Indeed, in large-scale problems, where the exact solution of sub-problems is a computational bottleneck, this relaxation is crucial; see for precise definitions as well as ways to obtain the approximate solution of the sub-problems (1a) and (1b).

Numerical Experiments

We are now ready to empirically evaluate the performance of the Newton-type methods considered in this paper (and studied theoretically in ) in several settings. In particular, we study the answers to Q.1–Q.6 posed at the outset in Section 1 in the context of two simple, yet illustrative, classes of non-convex optimization problems, i.e. (deep) MLP networks (Section 3.2) and NLS (Section 3.3). See the supplementary materials for more experiments.

In all of our experiments, we plot various quantities vs. total number of propagations , which is equivalent to measuring the number of oracle calls of function, gradient and Hessian-vector product. This is so since comparing algorithms in terms of “wall-clock” time can be highly affected by their particular implementation details as well as system specifications. In contrast, counting the number of propagations, as an implementation and system independent unit of complexity, is most appropriate and fair. Specifically, in neural nets, for a given data at the input layer, evaluation of network’s output layer involves one forward propagation. Performing one additional backward propagation gives the corresponding gradient. Each of the Hessian-vector products, required to solve the respective sub-problems of the second-order methods, is equivalent to two gradient evaluations, i.e., compared to the gradient, it involves one additional forward and backward propagations . Combining the batch size in each iteration, we summarize the number of propagations per iteration for each algorithm in Table 2.

In all of the following experiments, to approximately solve the sub-problems (1a) and (1b), respectively, we use CG-Steihaug method and the generalized Lanczos method with a maximum of 250 Lanczos iterations.

2 Hessian-Free Optimization for MLPs

In this section, we set out to empirically study the main questions Q.1–Q.4 of Section 1. For this, we consider two simple but non-trivial MLP netwroks under various settings, namely (a) 1-hidden layer neural networks for image classification, and (b) deep auto-encoder. For the following experiments, the empirical performance of the following algorithms, in light of Questions Q.1–Q.4, are compared We excluded ARC variants for reasons alluded to in Section 4.:

TR Uniform: Algorithm 1 with uniform sub-sampling,

GN Uniform: sub-sampled variants of Gauss-Newton method with modifications introduced in , and

SGD with Momentum: mini-batch SGD with momentum term and fixed step-size.

For Algorithm 1, we set $\eta_{1}=10^{-4}$, $\eta_{2}=0.8$, $\gamma_{1}=1.2$, and $\gamma_{2}=2$ (these are some values typical used in literature). The sample size used for Hessian sub-sampling is set to $5\%$ of the total training set, i.e., $0.05n$. For momentum SGD, the mini-batch size is chosen as $0.05n$ and the momentum parameter is set as $0.9$ (which is typically set in literature). For GN method, we use the same parameter settings as in .

2.1 One-Hidden Layer Neural Network

Here we consider a one-hidden layer neural network for the task of image classification using cifar10 data set. The hidden layer size is $512$, amounting to $d=1,578,506$. Initializations are done by setting $\mathbf{x}_{0}$ to a normalized vector drawn from standard normal distribution and all-zeros vector. Figure 1 shows training loss, training error and test error for all methods.

In light of Q.1–Q.4, we make the following observations:

(Computational Efficiency) In Figures 1(a)(b)(c), all algorithms start from a normalized random vector. From Figure 1(a), we observe that in terms of training loss, sub-sampled TR algorithm, though very competitive, can be slightly slower than SGD as long as SGD’s step size is appropriately fine-tuned. However, this does not necessarily translate to better test error; see Figure 1(c). In particular, while all TR runs achieve similar test errors, SGD’s generalization performance does not mirror its training behavior and appears highly “chaotic”, i.e., the step-size that achieved the fastest training (red dashed line), generalizes very poorly (compare red and blue dashed lines).

(Robustness to Hyper-parameters) From Figures 1(a)(b)(c), we notice the dependence of SGD’s performance on the choice of its learning-rate. Well-tuned step size can give both fast convergence of training process and good generalization performance. Too small a step size, however, can lead to slow convergence, while too large a step size, can cause SGD divergence or poor generalization performance. In contrast, Algorithm 1 with drastically different initial trust-region radii exhibit comparable performances; similar phenomenon are seen in Figures 1(d)(e)(f).

(Escaping Saddle Point) In Figure 1(d)(e)(f), all algorithms start from the origin. We can clearly see that SGD and GN easily get trapped at/near saddle points and/or flat regions (it is easy to check the gradients there are extremely small) and can barely make any progress. In contrast, sub-sampled TR, which effectively utilizes the Hessian, seamlessly escapes these regions and makes continued progress.

(Generalization Performances) In all initialization schemes, as shown in Figures 1, sub-sampled TR obtains competitive, if not better, generalization performance.

2.2 Deep Auto-Encoder

Here, we consider the deep auto-encoder problem and use the same model architectures as well as loss functions as in . The dataset and network architectures are given in Table 3. The experiments in Figures 2 and 3 are each done with initialization to a vector drawn from standard normal distribution as well as the all-zeros vector.

Although deep auto-encoder networks here are much more complex than 1-hidden layer network of Section 3.2.1, but we can still make very similar observations as they relate to Q.1–Q.4. In particular, on both datasets, Algorithm 1 converges comparably as fast, or faster than, other methods in terms of number of propagations. More importantly, it exhibits great robustness to its hyper-parameter, the initial trust-region radius, in contrast to SGD with momentum, which is heavily dependent on the choice of step size. Since these are rather complex networks, the optimization landscape is riddled with saddle and/or very flat regions. In this light, SGD and GN algorithms can both, rather easily, get trapped at/near these regions, which is quite contrary to the behavior of Algorithm 1 (Figures 2(b) and 3(b)). In terms of test error, Algorithm 1 is competitive to other methods on some experiments (Figures 2(a) and 3(a)) and clearly outperforms on others (Figures 2(b) and 3(b)).

We further examine two additional initialization strategies: normalized random initial point (Figure 4(a)) as well as a random vector scaled by $0.25$ (Figure 4(b)). Similar observations as earlier can also be made here. In particular, unlike Section 3.2.1, for this problem normalized random initial point (Figure 4(a)) seems to paint a different picture, i.e., SGD with momentum as well as GN all get trapped at high training levels while Algorithm 1 makes continued progress. It is worth mentioning that among all the initialization schemes we considered here, only under the particular scaled random initialization SGD can obtain desirable performance. This further demonstrates the versatility of Algorithm 1.

3 Non-Linear Least Squares

Table 4 summarizes the real data sets used for the experiments of this section. All datasets are from LIBSVM library . The following algorithms are compared (exact Hessian refers to $\mathbf{H}_{t}=\nabla^{2}F(\mathbf{x}_{t})$):

TR Full/Uniform/Non-Uniform: Algorithm 1 with full or uniform/non-uniform estimation of Hessian, respectively,

ARC Full/Uniform/Non-Uniform: Algorithm 2 with full or uniform/non-uniform estimation of Hessian, respectively,

GN Full/Uniform/Non-Uniform: GN with full or uniform/non-uniform estimation of Hessian, respectively,

LBFGS-100: the standard L-BFGS method with history size $100$ and using line-search.

For both Algorithms 1 and 2, we set $\eta_{1},\eta_{2},\gamma_{1}$ and $\gamma_{2}$, the same as Section 3.2. The sampling ratios, i.e., $|\mathcal{S}|/n$, for uniform and non-uniform sampling are set to $1\%$ and $0.1\%$, respectively. For all datasets, we set $\Delta_{0}=10$ for Algorithm 1 and $\sigma_{0}=10^{-4}$ for Algorithm 2. For GN, we do not regularize Hessian as in . Figure 5 gathers all comparison results of this section.

(Benefits of Sub-Sampling) Overall, both Algorithms 1 and 2 compare well with the classical TR and ARC methods. In fact, sub-sampling can, at times, help increase the efficiency, e.g., TR variants for covtype2. However, with too small a sample, the performance can hurt; ARC Full vs. ARC Non-Uniform and ARC Uniform for covtype2 and mnist2, respectively. The benefits of non-uniform sampling over uniform alternative are far more pronounced in the performance of Algorithm 2 than Algorithm 1. This can be attributed mainly to their respective sub-problem solvers in terms of total number of performed Hessian-vector products . In particular, CG-Steihaug used for the sub-problem (1a) of Algorithm 1 typically terminates in a handful of iterations whereas the generalized Lanczos method for solving the sub-problem (1b) of Algorithm 2 usually exhausts the allotted 250 iterations.

(Comparison Among Second-Order Methods) One can observe the consistent poor performances of L-BFGS-100 (green dotted lines) methods on all datasets, in particular with all $1$’s initialization vector. This is rather expected as contrary to popular belief, BFGS is not quite a “full-fledged” second-order method. Indeed, BFGS merely employs first-order information, i.e. gradients, to approximate the curvature, and starting from all $1$’s vector, L-BFGD cannot capture enough curvature information to navigate its way out of this region effectively. Gauss-Newton (dash lines), which has been specifically designed to effectively solve NLS problems, performs very well with random initialization. Starting from all $1$’s vector, however, where the gradient is very small, GN performs poorly. This is also expected because GN, similar in spirit to BFGS, does not fully utilize the Hessian information. In particular, in exchange for obtaining a positive definite approximation matrix, GN completely ignores the information from negative curvature, which is critical for allowing to escape from regions with small gradient. We can also observe that ARC is consistently no better than TR. This is an empirical evidence that the optimal worst-case complexity of ARC , though theoretically highly interesting, might be hard to observe in many practical settings.

Conclusion

In this paper, we aimed at painting a more complete picture of the practical advantages of Newton-type algorithms in general, and sub-sampled variants of trust-region and adaptive cubic regularization methods in particular, as compared with first-order alternatives. In the context of (deep) multi-layer perceptron networks as well as non-linear least squares, two simple, yet illustrative, non-convex machine learning applications, by making the following observations, we empirically attempted to make a case for the application of such second-order methods for machine learning.

(Computational Efficiency) The randomized sub-sampling approaches described here, and studied in detail in , can effectively make Newton-type methods computationally efficient enough to be competitive with popular first-order methods, widely used in machine learning, e.g., SGD with momentum. This is indeed due to the amortized combination of (1) low per-iteration cost offered by randomized sampling, and (2) small number of overall iterations due to the application of curvature information.

(Robustness to Hyper-parameters) In contrast to first-order algorithms whose performance is greatly affected by the choice of hyper-parameters, most notably step-size, the performance of the proposed Newton-type methods exhibit great robustness to such parameter tuning.

(Escaping Saddle Point) A greatly beneficial advantage of employing Hessian information is that it allows for such Newton-type algorithms, unlike many first-order alternatives, to seamlessly escape regions near saddle points.

(Generalization Performance) Second-order methods prove beneficial for the downstream machine learning objective of obtaining good generalization error. In particular, one can obtain very good levels of prediction accuracy only after a few iterations of such methods. This is highly beneficial in, say, distributed settings where the communication across the network is the main computational bottleneck.

(Benefits of Sub-sampling) On several real datasets, we validated the effectiveness of sub-sampled Newton-type methods, as compared with classical versions, in speeding up computations. Also, the advantages of non-uniform sampling over the oblivious uniform alternative was verified.

(Comparison Among Second-Order Methods) There are clear advantages in using (sub-sampled) TR and ARC methods over other second-order alternatives, e.g., L-BFGS and GN, in terms of effective exploitation of curvature.

For the examples of Section 3.2, despite the best of our efforts, we were unable to obtain the expected performance of Algorithm 2 using a variety of implementations, e.g., our own hand-written code as well as some existing packages such as GALAHAD . We believe that this is tightly connected to the choice of the sub-problem solver in all these implementations. Since we were unable to pinpoint the source of the problem, we did not include Algorithm 2 in examples of Section 3.2.

Finally, we acknowledge that, although we presented various experiments, the empirical study of methods such as the ones considered here, takes more than a single “proof-of-concept” paper, and the results presented here should be viewed as merely a glimpse into their various properties.

Acknowledgment

We would like to acknowledge ARO, DARPA, Cray, and NSF for providing partial support of this work. FR gratefully acknowledges the support of the Australian Research Council through a Discovery Early Career Researcher Award (DE180100923). We also sincerely thank Profs. Dominique Orban, Nicholas I.M. Gould, and Coralia Cartis for kindly helping us with the code for adaptive cubic regularization as well as setting up GALAHAD package. We also greatly appreciate Dr. Felix Lenders’ help with the installation of the trlib package. We would like to thank Dr. Amir Gholaminejad for valuable comments on our empirical evaluations and suggestions for improving them.

References

Appendix A Image Classification with Cifar10

Here, we evaluate the sensitivity of various algorithms, when initialized with different random seeds. We do this by consider similar set up as in section 3.2.1. We present the 10 different runs of Algorithm 1 and SGD with momentum with fixed the configuration. Overall, we observed that algorithms did not show much sensitivity with respect to random seed.

Appendix B Non-Linear Least Squares

In this section, we provide more results on different datasets as the experiment in section 3.3. We use the exact same setup as before and the datasets we considered here are gathered in table 5. In this set of experiments, we consider the sub-sampling ratio $1\%$ and $5\%$ of the training data for all sub-sampling methods.