Year: 2021

Abstract:
Stochastic gradient descent and other first-order variants, such as Adam and AdaGrad, are commonly used in the field of deep learning due to their computational efficiency and low-storage memory requirements. However, these methods do not exploit curvature information. Consequently, iterates can converge to saddle points and poor local minima. To avoid this, directions of negative curvature can be utilized, which requires computing the second-derivative matrix. In Deep Neural Networks (DNNs), the number of variables (n) can be of the order of tens of millions, making the Hessian impractical to store (O(n)) and to invert (O(n)). Alternatively, quasi-Newton methods compute Hessian approximations that do not have the same computational requirements. Quasi-Newton methods re-use previously computed iterates and gradients to compute a low-rank structured update. The most widely used quasi-Newton update is the L-BFGS, which guarantees a positive semi-definite Hessian approximation, making it suitable in a line search setting. However, the loss function in DNNs are non-convex, where the Hessian is potentially non-positive definite. In this paper, we propose using a Limited-Memory Symmetric Rank-1 quasi-Newton approach which allows for indefinite Hessian approximations, enabling directions of negative curvature to be exploited. Furthermore, we use a modified Adaptive Regularized Cubics approach, which generates a sequence of cubic subproblems that have closed-form solutions. We investigate the performance of our proposed method on autoencoders and feed-forward neural network models and compare our approach to state-of-the-art first-order adaptive stochastic methods as well as L-BFGS.

heading: 1 INTRODUCTION
text: Most deep learning problems involve minimization of the empirical risk of estimation
min Θ f(x; Θ), (1)
where Θ ∈ Rn is the set of weights and f is some scalar-valued loss function. To solve (1), various optimization approaches have been implemented, which we describe below. Throughout this paper, we write f(Θ) and f(x; Θ) interchangeably.
Gradient and adaptive gradient methods are widely used for training deep neural networks (DNN) for their computational efficiency. The most common approach is Stochastic Gradient Descent (SGD) which, despite its simplicity, performs well over a wide range of applications. However, in a sparse training data setting, SGD performs poorly due to limited training speed (Luo et al. (2019)). To address this problem, adaptive methods such as AdaGrad (Duchi et al. (2011)), AdaDelta (Zeiler (2012)), RMSProp (Hinton et al. (2012)) and Adam (Kingma & Ba (2014)) have been proposed. These methods take the root mean square of the past gradients to influence the current step. Amongst all of these adaptive methods, Adam is arguably the most widely used in a deep learning setting due to it rapid training speed.
Newton’s method has the potential to exploit curvature information from the second-order derivative (Hessian) matrix (see e.g., Gould et al. (2000)). Generally, the iterates are defined by Θk+1 = Θk − αk∇2f(Θk)−1∇f(Θk), where αk > 0 is a steplength defined by a linesearch criterion (Nocedal & Wright (2006)). In a DNN setting, we know that the number of parameters (n) of the network can be of the order of millions. Thus storing the Hessian which takes O(n2)
memory, becomes impractical. In addition, the inversion of the Hessian matrix, which takes O(n3) operations, is also impractical. Even though Newton’s method achieves convergence in fewer steps, the method becomes computationally intractable to use on large-scale DNNs.
Quasi-Newton methods are alternatives to Newton methods. They compute Hessian approximations, Bk+1, that satisfy the secant condition given by yk = Bk+1sk, where sk = Θk+1 −Θk and yk = ∇f(Θk+1) − ∇f(Θk). The most commonly used quasi-Newton method, including in the realm of deep learning, is the limited-memory BFGS update, or L-BFGS (see e.g., Liu & Nocedal (1989)), where the Hessian approximation is given by
Bk+1 = Bk + yky
> k
y>k sk − BkskskB
> k
s>k Bksk . (2)
The generic L-BFGS quasi-Newton update scheme is described in Algorithm 1, and numerous variants of L-BFGS exist (see Goldfarb et al. (2020); Moritz et al. (2016); Gower et al. (2016)). One advantage of using an L-BFGS update is that the Hessian approximation can be guaranteed to be definite, which is highly suitable in line-search settings because the update sk is guaranteed to be a descent direction, meaning there is some step length along this direction that results in a decrease in the objective function (see Nocedal & Wright (2006), Algorithm 6.1). However, because the L-BFGS update is positive definite, it does not readily detect directions of negative curvature for avoiding saddle points. In contrast, the Symmetric-Rank One (SR1) quasi-Newton update is not guarateed to be positive definite and can result in ascent directions for line-search methods. However, in trust-region settings where indefinite Hessian approximations are an advantage because they capture directions of negative curvature, the limited-memory SR1 (L-SR1) has been shown to outperform L-BFGS in DNNs for classification (see Erway et al. (2020)). We discuss this in more detail in Section 2 but in the context of Adaptive Regularization using Cubics.
Algorithm 1 L-BFGS Quasi-Newton Method with Line Search Require: Initial weights Θ0, batch size d, learning rate α, dataset D, loss function f(Θ).
for k = 0, 1, 2, . . . do Sample mini-batch of size d : Dk ⊆ D Perform the forward backward pass over the current mini-batch Compute the limited memory approximation Bk using (2) Compute step sk = αB−1k ∇Θf(Θk), where α is the line-search step length end for

heading: 2 L-SR1 ADAPTIVE REGULARIZATION USING CUBICS METHOD
text: We begin by discussing the L-SR1 update and the adaptive regularizion using cubics methods for large-scale optimization.
Unlike the BFGS update (2), which is a rank-two update, the SR1 update is a rank-one update, which is given by
Bk+1 = Bk + 1
s>k (yk −Bksk) (yk −Bksk)(yk −Bksk)> (3)
(see Khalfan et al. (1993)). As previously mentioned, Bk+1 in (3) is not guaranteed to be definite. However, it can be shown that the SR1 matrices can converge to the true Hessian (see Conn et al. (1991) for details). We note that the pair (sk,yk) is accepted only when |s>k (yk − Bksk)| > ε‖yk −Bksk‖22, for some constant ε > 0 (see Nocedal & Wright (2006), Sec. 6.2, for details). The SR1 update can be defined recursively as
Bk+1 = B0 + k∑ j=0
1
s>j (yj −Bjsj) (yj −Bjsj)(yj −Bjsj)>. (4)
In limited-memory SR1 (L-SR1) settings, only the last m n pairs of (sj ,yj) are stored and used. If Sk+1 = [ s0 s1 · · · sk ] and Yk+1 = [ y0 y1 · · · yk ], then Bk+1 admits a compact
representation of the form
Bk+1 = B0 + [ Ψk+1 ][ Mk+1 ] [ Ψ>k+1 ] , (5)
where
Ψk+1 = Yk+1−B0Sk+1 and Mk+1 = (Dk+1+Lk+1+L>k+1−S>k+1B0Sk+1)−1, (6)
where Lk+1 is the strictly lower triangular part, Uk+1 is the strictly upper triangular part, and Dk+1 is the diagonal part of S>k+1Yk+1 = Lk+1 + Dk+1 + Uk+1 (see Byrd et al. (1994) for further details).
Because of the compact representation of Bk+1, its partial eigendecomposition can be computed (see Burdakov et al. (2017)). In particular, if we compute the QR decomposition of Ψk+1, then we can write Bk+1 = B0 = U‖Λ̂k+1U>‖ , where U‖ ∈ R n×(k+1) has orthonormal columns and Λ̂ ∈ R(k+1)×(k+1) is a diagonal matrix. If B0 = δkI (see e.g., Lemma 2.4 in Erway et al. (2020)), where 0 < δk < δmax is some scalar and I is the identity matrix, then we obtain the eigendecomposition Bk+1 = Uk+1Λk+1U>k+1, where Uk+1 = [U‖ U⊥], with U⊥ ∈ Rn×(n−(k+1)) and U>k+1Uk+1 = I. Here, (Λk+1)i = δk + λ̂i for i ≤ k+ 1, where λ̂i is the ith diagonal in Λ̂k+1, and (Λ)i = δk for i > k + 1.
Since the SR1 Hessian approximation can be indefinite, some safeguard must be implemented to ensure that the resulting search direction sk is a descent direction. One such safeguard is to use a “regularization” term.
The Adaptive Regularization using Cubics (ARCs) method (Griewank (1981); Cartis et al. (2011)) can be viewed as an alternative to line-search and trust-region methods. At each iteration, an approximate global minimizer of a local (cubic) model,
min s∈Rn
mk(s) ≡ g>k s + 1
2 s>Bks + µk 3 (Φk(s)) 3, (7)
is determined, where gk = ∇f(Θk), µk > 0 is a regularization parameter, and Φk is a function (norm) that regularizes s. Typically, the Euclidean norm is used. In this work, we propose an alternative “shape-changing” norm that allows us to solve each subproblem (7) exactly. This shape-changing norm was proposed in Burdakov et al. (2017), and it is based on the partial eigendecomposition of Bk. Specifically, if Bk = UkΛkU>k is the eigendecomposition of Bk, then we can define the norm ‖s‖Uk def = ‖U>k s‖3. Applying a change of basis with s̄ = U>k s and ḡk = U>k gk, we can redefine the cubic subproblem as
min s̄∈Rn
m̄k(s) = ḡ > k s̄ +
1 2 s̄>Λks̄ + µk 3 ‖s̄‖33 . (8)
With this change of basis, we can find a closed-form solution of (8) easily. The proposed Adaptive Regularization using Cubics with L-SR1 (ARCSLSR1) algorithm is given in Algorithm 2.

heading: 2.1 CONTRIBUTIONS
text: The main contributions of this paper are as follows: 1. L-SR1 quasi-Newton methods: The most commonly used quasi-Newton approach is the L-BFGS method. In this work, we use the LSR1 update to better model potentially indefinite Hessians of the non-convex loss function. 2. Adaptive Regularization using Cubics (ARCs): Given that the quasi-Newton approximation is allowed to be indefinite, we use an Adaptive Regularized using Cubics approach to safeguard each search direction. 3. Shape-changing regularizer: We use a shape-changing norm to define the cubic regularization term, which allows us to compute the closed form solution to the cubic subproblem (7). 4. Computational complexity: Let m be the number of previous iterates and gradients stored in memory. The proposed LSR1 ARC approach is comparable to L-BFGS in terms of storage and compute complexity (see Table 1).
Algorithm 2 Limited-Memory Symmetric Rank-1 Adaptive Regularization using Cubics 1: Given: Θ0, γ2 ≥ γ1, 1 > η2 ≥ η1 > 0, and σ0 > 0 2: for k = 0, 1, 2 . . . do 3: Obtain Sk = [ s0 · · · sk ], Yk = [ y0 · · · yk ] 4: Solve the generalized eigenproblem S>k Yku = Λ̂S > k Sku and let δk = min{λ̂i}
5: Compute Ψk = Yk − δkSk 6: Perform QR decomposition of Ψ = QR 7: Compute the eigendecomposition RMR> = PΛP> 8: Assign U‖ = QP and U>‖ = P >Q> 9: Define C‖ = diag(c1, . . . , cm), where ci = 2 λi+ √ λ2i +4µ|ḡi| and ḡ‖ = U>‖ g
10: Compute α∗ = 2 δk+ √ δ2k+4µ‖g⊥‖ where g⊥ = g −U‖ḡ‖

heading: 11: Compute step s∗ = −α∗g + U‖(α∗Im −C‖)U>‖
text: 12: Compute m(s∗) and ρk = (f(Θk)− f(Θk+1))/m(s∗) 13: Set
Θk+1 = { Θk + sk, if ρk ≥ η1, Θk, otherwise and µk+1 =  0.5µk if ρk > η2, 0.5µk(1 + γ1) if η1 ≤ ρk ≤ η2, 0.5µk(γ1 + γ2) otherwise
14: end for

heading: 2.2 IMPLEMENTATION
text: Because full gradient computation is very expensive to perform, we impement a stochastic version of the proposed ARCs-LSR1 method. In particular, we use the batch gradient approximation
g̃k ≡ 1 |Bk| ∑ i∈Bk ∇fi(Θk).
In defining the SR1 matrix, we use the quasi-Newton pairs (sk, ỹk), where ỹk = g̃k+1 − g̃k (see e.g., Erway et al. (2020)).

heading: 3 CONVERGENCE ANALYSIS
text: In this section, we prove convergence properties of the proposed method (ARCs-LSR1 in Algorithm 2). The following theoretical guarantees follow the ideas from Cartis et al. (2011) and Benson & Shanno (2018).
First, we make the following mild assumptions:
A1. The loss function f(Θ) is continuously differentiable, i.e., f ∈ C1(Rn). A2. The loss function f(Θ) is bounded below.
Next, we prove that the matrix Bk in (4) is bounded.
Lemma 1 The SR1 matrix Bk+1 in (4) satsifies
‖Bk+1‖F ≤ κB for all k ≥ 1 for some κB > 0.
Proof: Using the limited-memory SR1 update with memory parameter m in (4), we have
‖Bk+1‖F ≤ ‖B0‖F + k∑
j=k−m+1
‖(yj −Bjsj)(yj −Bjsj)>‖F |s>j (yj −Bjsj)| .
Using a property of the Frobenius norm, namely, for real matrices A, ‖A‖2F = trace(AA >), we have that ‖(yj − Bjsj)(yj − Bjsj)>‖F = ‖yj − Bjsj‖22. Since the pair (sj ,yj) is accepted only when |s>j (yj −Bjsj)| > ε‖yj −Bjsj‖22, for some constant ε > 0, and B0 = δkI for some 0 < δk < δmax, we have
‖Bk+1‖F ≤ δmax + m
ε ≡ κB .
Given the bound on ‖Bk+1‖F , we obtain the following result, which is similar to Theorem 2.5 in Cartis et al. (2011).
Theorem 1 Under Assumptions A1 and A2, if Lemma 1 holds, then
lim inf k→∞ ‖gk‖ = 0.
Finally, we consider the following assumption, which can be satisfied when the gradient, g(Θ), is Lipschitz continuous on Θ.
A3. If {Θti} and {Θli} are subsequences of {Θk}, then ‖gti−gli‖ → 0 whenever ‖Θti−Θli‖ → 0 as i→∞. If we further make Assumption A3, we have the following stronger result (which is based on Corollary 2.6 in Cartis et al. (2011)):
Corollary 1 Under Assumptions A1, A2, and A3, if Lemma 1 holds, then
lim k→∞ ‖gk‖ = 0.

heading: 4 EXPERIMENTS
text: To empirically compare the efficiency of the method against popular optimization methods like SGD, ADAGRAD, ADAM, RMSProp and L-BFGS, we focus on two broad deep learning problems: image classification and image reconstruction. We choose these tasks due to their broad importance and availability of reproducible model architectures. We run each experiments on an average of 5 times with a random initialization in each experiment. The number of parameters, convolutional layers and fully connected layers are mentioned in Table 3.
Dataset: We measure the classification performance of each optimization method on 4 image datasets: MNIST (LeCun et al. (2010)), FashionMNIST (Xiao et al. (2017)), IRIS (Dua & Graff (2017)) and CIFAR10 (Krizhevsky et al.). We have provided a comprehensive view of the experiments in Table 2
Hyperparameter tuning: We empirically fine-tune the hyperparameters and select the best for each update scheme. We have made a comprehensive list of all the learning rates for the gradient and adaptive gradient based algorithms in Table 4 in the Appendix. The additional parameters are defined as follows:
• ADAM: We apply an perturbation of 1.0 × 10−6. β0 and β1 are chosen to be 0.9 and 0.999, respectively. • ADAGRAD: The initial accumulator value is set to 0. The perturbation is set to 1.0× 10−10. • SGD: We use a momentum of 0.9. • RMSPROP: We set α = 0.99. The perturbation is set 1.0× 10−8. • L-BFGS: The table 4 in Appendix A presents the initial learning rate for the stochastic step in
L-BFGS. We set the default learning rate to 1.0. We choose a history size m of 10 and max iterations to 10. The tolerance on function value/parameter change is set to 1.0 × 10−9 and the first-order optimality condition for termination is defined as 1.0× 10−9.
• ARC-LSR1: We choose the same parameters as L-BFGS.
Network architecture: For each problem, we define the model architecture in Table 3 in the appendix. We define the process of the forward and backward pass of a DNN in Algorithm 3 in the appendix.
Testbed and software: All experiments were conducted using open-source software PyTorch (Paszke et al. (2019)), SciPy (Virtanen et al. (2020)) and NumPy (Harris et al. (2020)). We use an Intel Core i7-8700 CPU with a clock rate of 3.20 GHz and an NVIDIA RTX 2080 Ti graphics card.

heading: 5 RESULTS
text: We have divided the sections into two categories: classification and image reconstruction. We present both the training results and the testing results for all methods.

heading: 5.1 CLASSIFICATION RESULTS
text: For each classification problem, we define the network architecture, the corresponding hyperparameters (other than the learning rate) for each optimization scheme.
IRIS: Since this dataset is relatively small, we assume a small network for our deep-learning model. The model is described in 3. We set the history size for the proposed approach and L-BFGS to 10 and the number of iterations to 10. Figure 1 shows the comparative performance of all the methods. Note that our proposed method (ARCLSR1) achieves the highest classification accuracy in the fewest number of epochs.
MNIST: We trained the network for 20 epochs with a batch size of 256 images each. We keep the same history size and number of iterations as the IRIS dataset for L-BFGS and the proposed ARCLSR1 approach. For training, it can be seen in Figure 2 that nearly all methods achieve optimal training accuracy. However, closely inspecting the testing curve, we notice that the proposed approach achieves higher accuracy than all the existing methods.
FMNIST: We train the network for 20 epochs with a batch size of 256 images. We keep the history size the same as the IRIS and MNIST experiments for the proposed approach and L-BFGS. For this method, the proposed ARCLSR1 approach is comparable to L-BFGS but outperforms the adaptive methods (see Figure 3).
CIFAR10: We use the same parameters presented in Table 4 in the previous section for the adaptive methods. For ARCLSR1 and L-BFGS, we have a history size of 100 with a maximum number of iterations of 100 and a batch size of 1024. Figure 4(a) represents the training loss (cross-entropy
loss). Figure 4(b) represents the testing accuracy, i.e., number of sample correctly predicted in the testing set. To demonstrate the efficacy of the proposed method on larger networks, additional experimentation on the ResNet50 architecture can be found in the appendix (Figure 8).

heading: 5.2 IMAGE RECONSTRUCTION RESULTS
text: The image reconstruction problem involves feeding a feedforward convolutional autoencoder model (with randomly initialized weights) a batch of the dataset. It follows the same deep learning convention as mentioned in Algorithm 3 in Appendix A. The loss function is defined between the reconstructed image and the original image.
MNIST: An image x ∈ Rn is fed to the network, compressed into a latent space z ∈ Rl, where l n, and reconstructed back to its original image size x̄ ∈ Rn. We compute the mean-squared loss error between the reconstruction and the true image. The weights are initialized randomly. Each experiment has been conducted 5 times and we considered a batch size of 256 images each with 50 epochs. The results for the image reconstruction can be seen in Figure 5.
One can notice that the initial descent provided by the proposed approach provides a significant decrease in the objective function. To understand better, we provide the details of the results during the initial epoch (Figure 9(a)) and the final epoch (Figure 9(b)). We notice that the ARCLSR1 method has minimized efficiently in the first half of the first epoch. This is empirical evidence that the method converges to the minimizer in fewer steps in comparison to the adaptive methods. In Figure 5 (b), we notice that all the adaptive methods eventually converge to the same point. For training response results on the F-MNIST dataset, see Section B in the appendix.

heading: 5.3 TIME COMPLEXITY ANALYSIS
text: We understand that the proposed approach performs competitively against all existing methods. We now analyze the time-constraints of each method. We choose to clock the computationally demanding algorithm here - CIFAR10 classification. We chose a maximum iterations of 100 with a history size of 100 for L-BFGS and the ARCs LSR1, with a batch size of 1024 images. Figure 7 plots the time required by each of the methods to reach non-overtrained minima with a batch size of 1024 images. As can be seen, the proposed approach is able to reach the desired minima in much less time than the rest of the algorithms. L-BFGS finds it hard to converge due to a very noisy loss function and a small batch size, thus causing the algorithm to break. Ozyildirim & Kiran (2020) argue that a large batch size is required for quasi-Newton methods to perform well. However, the ARCLSR1 method performs well with a small batch size as well.

heading: 6 CONCLUSION
text: In this paper, we proposed a novel quasi-Newton approach in a modified adaptive regularized cubic setting. We were able to empirically and theoretically show how an L-SR1 quasi-Newton approximation in an ARCs setting was able to perform either better or comparably to most of the state of the art optimization schemes. Even though the approach has yielded exceptional results, we need to test the method’s efficacy when the network size and dataset size is large and when availability of data is sparse.

Review summaries:
The paper describes a limited memory quasi-Newton method based on SR1 updating using a variant on the adaptive regularized cubic (ARC) approach to globalization.  The algorithm is applied to training deep neural networks for image classification and autoencoding, and compared to L-BFGS and various SGD variants.  The authors claim the main contributions to be the different optimizer ingredients (L-SR1, ARC safeguarding, and the particular shape-changing regularizer) as well as computational complexity similar to L-BFGS.

This paper investigates the application of a certain Quasi-Newton algorithm, the Limited-Memory Symmetric Rank-1 (L-SR1) algorithm, in deep learning problems. The benefit of this technique over similar more widely investigated methods that use a positive definite approximation of the Hessian, such as stochastic L-BFGS, is the fact that the L-SR1 approximation is not guaranteed to be definite, and thus has the potential to have a more accurate approximation of the true Hessian. Since in this case line-search methods may return ascent directions, authors propose a specific form of Adaptive Regularization using Cubics (ARC) as an alternative. Numerical simulations are provided comparing the performance of the proposed algorithm to SGD, adaptive methods such as Adam and a naive L-BFGS implementation.