Title: EXPONENTIALLY DECAYING FLOWS FOR OPTIMIZA-TION IN DEEP LEARNING

Abstract: The field of deep learning has been craving for an optimization method that shows outstanding property for both optimization and generalization. We propose a method for mathematical optimization based on flows along geodesics, that is, the shortest paths between two points, with respect to the Riemannian metric induced by a non-linear function. In our method, the flows refer to Exponentially Decaying Flows (EDF), as they can be designed to converge on the local solutions exponentially. In this paper, we conduct experiments to show its high performance on optimization benchmarks (i.e., convergence properties), as well as its potential for producing good machine learning benchmarks (i.e., generalization properties).

1. INTRODUCTION:

Due to recent progress in the field of machine learning, it becomes more and more important to develop and sophisticate methods of solving hard optimization problems. At the same time, in this field, such methods are additionally required to elicit decent generalization performance from statistical models. An efficient method of mathematical optimization, however, does not always produce sufficient generalization properties, since these are involved with two distinct mathematical problems; The former is to find one of the solutions which minimize a given (possibly non-convex) objective function, and the latter is to adjust parameters so that a statistical estimator achieves its best result. To address such a hard issue, we introduce a new mathematical perspective on optimization, and develop a method for machine learning based on this perspective. We then empirically show its rapid convergence rate and high compatibility with deep learning techniques, as well as good statistical properties.

In this field, many optimization methods have been proposed and modified so that they fit specific problems or models. One of the current standard methods is the gradient descent method. The method tends to converge slowly in general optimization problems. However, with various specific techniques, such as mini-batch training and batch normalization (Ioffe & Szegedy (2015) ), it has been found to be efficient for state-of-the-art purposes in the field of deep learning. Another class of methods that are now becoming popular and standard is adaptive methods, such as AdaGrad (Duchi et al. (2011) ) and Adam (Kingma & Ba (2015) ). Compared to the gradient descent method, these methods have been shown to improve convergence rates with almost the same computational cost as the gradient descent method, but are reported to result in poor statistical outcomes in some cases of machine learning (Wilson et al. (2017) ).

Other class of methods that have been thoroughly studied in the theory of mathematical optimization is second-order methods, such as the Newton method and the Gauss-Newton method. These methods possess great convergence properties, and in particular, have a potential to overcome plateau's problems (Dauphin et al. (2014) ). Furthermore, when it comes to applications in stochastic settings, the method based on the Gauss-Newton Matrix (or Fisher information Matrix) is shown to asymptotically attain the best statistical result, which is called Fisher efficiency (see Amari (1998) ). Despite these attractive characteristics, the methods have not yet been spotlighted in the field of machine learning due to several severe drawbacks; They suffer from high computational cost in general and their useful properties are no longer guaranteed in practical settings (see Section 12 in Martens (2014) ). One of the continuously developing second-order methods in this field, K-FAC ( Ba et al. (2017) , Grosse & Martens (2016) ), successfully produced high convergence rate empirically with relatively low computational cost. However, it still requires much effort to become compatible with some deep learning techniques. In addition, it is unclear whether the method has advantages in generalization performance.

In our approach, by introducing a Riemannian metric induced by non-linear functions, we constitute dynamical systems which describe motions along the shortest route from arbitrary initial points to the zeros of non-linear functions on the corresponding Riemannian manifold, that is, geodesic with respect to the Riemannian metric. One of the remarkable characteristics of our approach is that it enables us to flexibly design flows of such dynamical systems to control convergence rates. The results for the flows are then applicable to mathematical optimization problems, in particular, with deep neural network (DNN) models. In this paper, after providing mathematical ground of our methods, we experimentally demonstrate their performance in various aspects, from convergence rates to statistical properties.

2. GEODESIC FLOWS FOR NON-LINEAR EQUATIONS:

We start by establishing some essential properties of dynamics which are effective for the analysis of non-linear equations. Let F : R N → R N be a smooth function and J be the Jacobian with variable w ∈ R N , that is, J = ∂F/∂w. In this section, we deal with the well-posed case that there exists a connected closed subset Ω ⊂ R N , where J is regular and the equation has a unique solution ξ. Therefore, the positive matrix G = J T J induces a Riemannian metric g on Ω and (Ω, g) becomes a Riemannian manifold under some appropriate conditions. Let us then consider time evolution of variable w on this manifold. Our main purpose in this section is to study the characteristics of dynamical systems in which w(t) moves on geodesics between any point in Ω and ξ with respect to the metric g.

Let L be a Lagrangian given by

EQUATION

with v = dw/dt (also written as ẇ). The Euler-Lagrange equation for L is then expressed as

dp dt -∇ w L = 0 (2)

with momentum vector p = Gv. If the boundary condition at two points in Ω, w(t 0 ) = w 0 , w(t 1 ) = w 1 , is imposed on (2), a geodesic between w 0 and w 1 is obtained as the solution. In contrast, if we give an appropriate initial condition, w describes the motion along the geodesic from w 0 to w 1 . In fact, the following statement holds; Theorem 2.1. Let w 0 be an arbitrary point in Ω and (w(t), p(t)) be a solution of equation ( 2) with the following initial condition;

EQUATION

Then w satisfies F w(t) = (1 -t)F w 0 , (4) for t ∈ [0, 1]. In particular, w(t) passes through the point which is a solution of non-linear equation

F (w) = 0 at t = 1, that is, ξ = w(1).

We briefly describe the outline of the proof for the statement above. Throughout this paper, we regard functions of w as those of t in the trivial way; for instance, F (t) = F (w(t)). Note that p can be expressed as p = J T dF/dt. Then using the Beltrami identity for (2) with the initial condition above leads to the equation

EQUATION

where F 0 = F (0). Thus, a closed form expression is obtained as

EQUATION

which gives F (1) = 0 as asserted in Theorem 2.1. Now we take a different expression that the coefficient (1 -t) in ( 6) is replaced by a different monotonically decreasing function, that is,

EQUATION

where ρ denotes a monotonically decreasing smooth function from (t 0 , t 1 ) onto (0, 1). Then, we give the following differential equation whose solution is of the closed form (7);

EQUATION

A motion described by this differential equation differs from the one that is described by (2), but these two motions are along the same geodesic. Theorem 2.2. Let w 0 ∈ Ω be an arbitrary point. The differential equation

EQUATION

with an initial condition w 0 = w(t 0 ) has a unique solution that satisfies F (w(t 1 )) = 0. Furthermore, the orbit under flow f defined by f (w 0 , t) = w(t) coincides with that of the geodesic equation ( 2).

Note that since equation ( 9) is equivalent to equation ( 8), the orbit is invariant under coordinate transformations.

With respect to the choice of ρ, the end point t 1 can be set as ∞ under some appropriate conditions and Theorem 2.2 still holds in the sense that

EQUATION

In particular, if we set ρ(t) = e -t , then χ(t) = -1 and F can be represented as F (t) = e -t F 0 , so that the convergence rate of F is exponential. Definition 2.3. Let w be a solution of the differential equation

EQUATION

with an initial condition w 0 = w(t 0 ). The flow f (w 0 , t) = w(t) is called an exponentially decaying flow (EDF) of non-linear equation F (w) = 0.

For the end of this section, we present a connection between EDF and the classical Newton method.

If we apply the Euler method to the differential equation ( 9) with step size t, the corresponding iteration step can be written as

EQUATION

which recovers the Newton method with step size η i = t • χ(τ i ).

3. APPLICATION TO MATHEMATICAL OPTIMIZATIONS:

Consider smooth functions ϕ : R N → R M and L : R M → R. In this section, we develop a method based on EDF for a mathematical optimization problem of the form

EQUATION

In this section, the area Ω ⊂ R N is supposed to be a compact subset where there exists no stationary point except for a minima.

Let F : R N → R M denote derivative ∂L/∂ϕ. An example of such problems is least squares one in which L is given by

EQUATION

)

In this case, F = ϕ. In particular, if M = N and a minimal value of loss function L = L • ϕ is zero, the optimization problem is equivalent to solving non-linear equation F (w) = 0.

For optimization problems, we set up a standard equation that the gradient of loss function is zero, that is, ∇ L(w) = 0. Note that the gradient can be expressed as ∇ L(w) = J T ϕ F (w) with Jacobian J ϕ = ∂ϕ/∂w. Applying Theorem 2.2, we obtain the differential equation

EQUATION

)

where H is the Hessian of loss function L with respect to w. Since second order differentiation is typically computationally heavy, we seek other equations which possess almost the same property as (15) especially in terms of asymptotic behavior.

Let us decompose momentum H ẇ as

EQUATION

where J F denotes the Jacobian of F , and G is a symmetric matrix defined by

EQUATION

)

with Hessian matrix H L of L with respect to ϕ. We then consider the following equation instead of (15);

EQUATION

)

This equation no longer describes the motion along the geodesic related to (15) in general. However, if M = N and J ϕ is invertible, then w moves on another geodesic with respect to a different metric G F = J T F J F . In addition, if ρ(t) = e -t , F converges exponentially, which implies that ∇ L = J ϕ F also converges exponentially in (18) as well as in (15). In general cases, if a condition that d dt

1 2 F 2 = χ J F G -1 J T ϕ F, F ≤ χ F 2 (19)

is satisfied, then F converges to 0. This shows that in the neighborhood of solution ξ of equation ∇ L = 0, the momentum G ẇ sufficiently approximates H ẇ by ( 16). Definition 3.1. The flow given by

EQUATION

is referred to EDF of type H and written as EDF-H. Similarly, the flow given by

G dw dt = -J T ϕ F (21)

is referred to EDF of type G and written as EDF-G.

4. MODIFICATION SCHEMES FOR EDF-BASED METHODS:

Like second order methods, in EDF-based methods, matrix inversion has to be carried out, which requires expensive computational cost, particularly in large systems. Moreover, we often encounter rank deficient matrices in optimization problems. To deal with rank deficiency, in general, we need pseudo-inverse matrices which are more computationally heavy. Therefore, instead of the inverse of matrix A = G, H, for fixed v ∈ R M , we consider a projection which maps r = arg min x∈R M Ax -v to a vector P k (A, v) in the k-th order Krylov subspace

K k (A, v) = span{v, Av, A 2 v, . . . , A k-1 v} such that r = P ∞ (A, v).

One of the basic methods to construct such a projection for indefinite symmetric systems is the minimum residual method (Paige & Saunders (1975)), which requires only the matrix multiplication. Therefore, for numerical computations, we use the following differential equation approximated by the k-th order Krylov subspace;

EQUATION

k is a hyperparameter that interpolates between the gradient method and the method based on (15) or (18). In fact, in the case that k = 1, equation ( 22) has the form

EQUATION

which reduces to a kind of gradient methods, but the coefficient c conveys information about G or H unlike the standard gradient methods.

Next, similar to the Levenberg-Marquardt algorithm, we modify equation ( 18) by adding a damping factor to G in order that the method becomes more stable. So, we set

EQUATION

)

where λ is a continuous positive function of t. Then, we take the same approximation as ( 22) with A = G+λI. The damping factor λ plays several important roles in solving practical problems. First, it has solutions well-behaved near singularities. Even in the case that G is not invertible, equation ( 24) can be defined and solved unlike (18). If we choose λ such that it approaches to 0 as rapidly as the gradient in ( 18), the asymptotic behavior of ( 24) is almost the same as that of (18). Particularly, in the case that χ = -1, we set λ = a J T ϕ F b with a, b > 0, so that the convergence rate of (24) stays exponential. Second, the damping factor makes the method compatible with stochastic approaches such as mini-batch training, in deep learning models. Since the orbit of ( 24) is affected by the gradient due to the damping factor λ, the method based on this equation could take advantage of stochastic approaches as other gradient methods do. (For implementation of the algorithm, see Appendix A.) Finally, to accelerate the EDF-based methods, it is sometimes effective to change equation ( 18) into a second-order differential equation, particularly in the case in which the approximation method with k is applied. Specifically, we take the equation

EQUATION

)

where κ is a real-valued function of t. The convergence properties of the methods are controlled by the following differential equation;

EQUATION

There are two ways of determining κ. The first one is to set κ = α with constant α (W1), which leads to a similar scheme to the momentum gradient decent method. The other one is to set κ(t) = αt -1 (W2). In this setting, the equation can be discretized in a similar way as described in Su et al. (2014) , which is analogous to the Nesterov's acceleration scheme.

5. OPTIMIZATION PROBLEMS IN THE FIELD OF DEEP LEARNING:

In this section, we present how optimization problems are set up in the field of deep learning. Let x = {x j } n-1 j=0 and y = {y j } n-1 j=0 denote training datasets of input and output, respectively, where n is the size of dataset. Let d x and d y be dimensions of input data and output data, respectively. We write ϕ nn for neural networks with parameters w ∈ R N , and define ϕ by the direct sum of vectors {ϕ j } given by ϕ j (w) = ϕ nn (x j , w), that is, ϕ = ⊕ϕ j . Note that M = n × d y in this case. Then finding a minima of a given loss function is proposed as a standard optimization problem to train networks.

For the image classification tasks, there are two typical cases of setting loss functions. In the first case, the loss is set as ( 14). As already mentioned, in this case, F = ϕ and H L = I. In the second case, the loss is given by cross entropy with softmax function, that is,

EQUATION

where θ denotes the softmax function and θ j = θ(ϕ j ). In this case, F is expressed by the direct sum such that F = ⊕F j with

EQUATION

where s j denotes a sum of all elements in vector y j for each j. Note that if each y j is given as a probability vector, then s j = 1. Moreover, H L is expressed as H L = ⊕H j with

EQUATION

where ⊗ denotes the outer product. In both cases, the loss functions take the minimum value 0 if and only if F = 0. 

6. EXPERIMENTS:

In this study, we conducted following three groups of experiments; First, we examined optimization performance of the EDF-based methods on a data-fitting problem with a simple convolutional neural network model. Second, we tested both optimization and generalization performance in standard settings of classification tasks in which we employed residual network models (Resnet He et al. ( 2016)) and CIFAR-10/100 datasets. Third, we incorporated some techniques of data augmentation and regularization in the training into our experiment as we would like to measure effectiveness of our methods in more practical settings.

The primary purpose of the experiments in this paper is not to pursue the state-of-the-art performance, but to explore how the statistical results pertain to those of optimization when the loss functions are minimized. Therefore, we tuned hyperparameters of each method in accordance with the optimization benchmark.

It should be noted that the conjecture concerning non-convex optimization problems in the field of deep learning (Baldi & Hornik (1989) , Choromanska et al. (2015) ) is still an open problem (studied for linear cases in Baldi & Lu (2012) , Kawaguchi (2016) ). Hence, for experiments in this paper, we do not discuss whether each optimization method actually reaches to a global solution.

6.1. OPTIMIZATION PERFORMANCE ON A STANDARD DATA-FITTING PROBLEM:

We evaluated convergence performance of EDF-based methods (type G and type H) on a data-fitting problem of CIFAR-10, that is, full-batch training in the context of deep learning. The model we employed in these experiments was the convolutional neural network that consisted of two convolution filters with rectified linear units and max pooling layers, followed by two fully connected layers. For EDF-based methods, the step size was fixed to 1.0, and no damping factors were used. In addition, the acceleration technique derived from W2 of (25) was adapted, since W2 achieved better performance than W1. A similar experiment with a different type of second-order methods was conducted in Sohl-Dickstein et al. (2014) .

First, we examined change in convergence performance of EDF-G depending on hyperparameter k, the order of Krylov subspace. The results are illustrated in the left-hand side of Figure 1 . The presented trend is consistent with the theoretical fact that k interpolates between the gradient method (for small k) and the method based on dynamics (18) (for large k). In other words, when k is small, the method is similar to a gradient method, which converges slow, but as k becomes larger, the method leads to better approximation of the inverse matrix, which gives a rapidly decaying flow.

Next, we compared EDF-G (k = 1, 30) with EDF-H and other standard optimizers in deep learning: gradient descent methods with Polyaks momentum scheme (Momentum) and Nesterov's acceleration scheme (NAG), and Adam. The step sizes for Momentum, NAG, and Adam were fixed to 0.01, 0.001, and 0.001, respectively. The right-hand side of Figure 1 shows that both EDF-based methods made the loss decrease more rapidly than other standard optimizers. Even when EDF-based methods were reduced to the gradient methods at k = 1, EDF-G outperformed standard optimizers.

The figure also presents the difference between EDF-G and EDF-H. While their convergence rates around extremes were almost the same, the overall convergence rates of EDF-G were better than that of EDF-H.

As has been found, second-order-methods on full-batch training converge to the solution within a few iterations. However, with such a training, the quality of generalization becomes worse compared to the time when stochastic approach with a mini-batch of small size is taken. Like other secondorder-methods, EDF also suffers from a negative impact of full-batch training on generalization performance, though setting the hyperparameter k small makes EDF be compatible with stochastic approaches and take their advantages (see Appendix B).

6.2. PERFORMANCES ON CLASSIFICATION TASKS IN DEEP LEARNING:

In the following experiments, we compared both optimization and generalization performance of EDF-G with other methods, momentum stochastic gradient decent method (MSGD) and Adam on classification tasks for CIFAR-10/100. The experiments were conducted, employing residual network models with batch normalization (Resnet-56 for CIFAR-10 and Resnet-110 for CIFAR-100), and working with a mini-batch of size 250.

For CIFAR-100, during pre-investigation of the dataset, we found that 9 pairs of data, each of which has the same image but different labels, contaminated the neural network and might have had a non-negligible influence on both optimization and generalization performance (See Appendix C). Therefore, in our experiments on CIFAR-100, we excluded them from the training dataset and used the rest of the data (size = 49982).

At the end of each epoch, we computed training loss in the following manner. First, the loss was defined as (27) where n is the total number of data. Second, to calculate the statistics for batch normalization in the loss, as described in Ioffe & Szegedy (2015) , we adopted so-called "inference mode" in which moving averages of mean and variance over mini-batches are used for batch-normalization.

For EDF, we tested the case k = 1 and k = 2, with the damping factor set as a = b = 1. The momentum coefficient was fixed to α = 0.9 for the W1 acceleration (see 25).

For each of the tasks, we ran tests on 10 different learning rates; for EDF between 0.1 and 2.0, for Momentum SGD between 0.1 and 10.0, for Adam between 0.0001 and 0.1. We then chose the one with the best convergence performance for each optimizer. For EDF with W2, to obtain stability in convergence, the learning rate was set to 0.2 times that of the initial values at the end of the 20-th epoch. Such a change in learning rate in the middle of optimization did not bring advantages to either optimization or generalization performances for other methods including EDF with W1 (see Appendix D).

The results are presented in Figures 2 and 3 . As shown in the figures, with respect to the optimization performance, EDF reached an optimal solution with smaller error at higher convergence rate than Adam and Momentum SGD, even when k = 1, 2. Moreover, EDF overall attained better generalization performance than other methods.

6.3. PERFORMANCES WITH TECHNIQUES IN A MORE PRACTICAL SETTING:

For classification tasks, generalization performance often improves by adopting several techniques, such as data augmentation and regularization. In this group of experiments, employing the data augmentation based on random horizontal flip and shift and the L 2 regularization, we conducted comparisons between EDF and other methods that are similar to those in the previous sections.

When adopting these techniques, rate decay scheme has as large impact on optimization performance as learning rate itself. Because effectiveness of each rate decay scheme much depends on optimization methods, it is almost impossible to find the best scheme that is applicable to all the methods. For this reason, for MSGD and Adam, as initial learning rates, we chose one of the standard rates for MSGD and Adam, 0.1 and 0.001, respectively, and at the ends of the 100-th, 140-th, 180-th epochs, reset the rates to 0.1 times the rates at the moment. For EDF, we ran tests on two different rates 0.5 and 0.75, the optimal rates found in the experiments of Section 6.2, and reset to 0.2 times those of the initial values, only at the end of the 100-th epoch, in order to demonstrate performance of the EDF more clearly. Among the results obtained with EDF, we chose the ones with the best optimization performance.

The result of the comparison is presented in Figures 4 and 5 . As can be seen, we found a condition in which EDF achieved better performance than other methods on optimization while achieving sufficient levels of generalization.

7. CONCLUSION:

Obtaining good statistical results from limited available data is a critical goal in machine learning. To reach this goal, while developing an effective model is an essential approach, eliciting the best performance from the fixed model through optimization is important as well. In our study, to examine the performance of our optimization methods, Exponentially Decaying Flows (EDF) based methods, we explored their generalization properties pertaining to results of optimization. Our experiments showed that EDF-based methods are more likely to achieve optimal solutions which generalize the test data well than other standard optimizers are. Therefore, EDF-based methods are considered to be optimization methods that have a high potential in their application to various tasks, and thus, are worthwhile to be sophisticated through future studies. Figure 6 shows the results of experiments using EDF for simple examples that compare a full-batch training with stochastic trainings. In this example, the convolutional network similar to that used in Section 6.1 was employed on MNIST. The curve labeled as "EDF F" depicts the result of Full-batch training per step, and those labeled as "EDF S" illustrate the results of stochastic trainings per epoch with a mini-batch of size 500.

Let us divide the dataset {x, y} = {x j , y j } n-1 j=0 into distinct subsets of size p such that {x, y} = l i=1

x (i) , y (i) , x (i) , y (i) = {x j , y j } pi-1 j=p(i-1) , pl = n.

(30)

Then, gradients and Jacobian-vector-multiplications are calculated as

EQUATION

where ϕ (i) and F (i) are the subcomponents of ϕ and F , respectively, corresponding to the decomposition above. Thus, for a fixed k, if we set p as the same size as a mini-batch, then the computational cost 

:



Figures and Tables:

Figure 1: Convergence curves for various optimizers

Figure 2: Training loss and test error for Resnet-56 on CIFAR-10

Figure 4: Training loss and test error for Resnet-56 on CIFAR-10 with data augmentation

Figure 6: Comparison between full-batch and stochastic trainings on MNIST

Peer review (iclr2019):
 I found the paper poorly written, with many typos and incorrect formulations. It contains several sophisticated mathematical concepts (probably from differentiable geometry, dynamical systems and differential equations) that I believe a majority of ICLR audience would not be so familiar with; and these notions are not defined nor explained (e.g Riemannian metrics/manifold, Krylov subspaces). Section 2 presents some theoretical results and formal definition of exponentially decaying flows but I could not see the intuition or the goal of these results. In Section 3, for the main optimization problem (13), it is assumed that "there exists no stationary point except for a minima." and so, I understand that there is only "a minima" or at least that all the local minima are assumed to have the same value, which not realistic for deep learning. Comments and questions Section 1: After speaking of gradient descent methods, the authors say "Another class of methods... is adaptive methods such as AdaGrad and Adam": These are all based on stochastic gradient descent algorithm, so this distinction is surprising. Section 2: "Jacobian with variable w": does not make sense (although I can guess the author mean w.r.t.) "J is regular and the equation has a unique solution": what is regular here? No equation has been mentioned before "...becomes a Riemannian manifold under some appropriate conditions": which conditions? Section 3: "Applying Theorem 2.2, we obtain the differential equation": I can*t see how applying Theorem 2.2 leads to a differential equation. Section 4: "we consider a projection which maps r to a vector P_k(A, v) in the k-th order Krylov subspace such that r = P_.infty(A, v)": I don*t understand. "Particularly, in the case that χ = −1, we set λ = akJ φ T F k b with a, b > 0, so that the convergence rate of (24) stays exponential": why is the case χ = −1 relevant to point out? "Finally, to accelerate the EDF-based methods, it is sometimes effective to change equation (18) into a second-order differential equation": is it an empirical observation? or is there an explanation? Section 5: Again, so many notations (not all conventional: e.g. .theta denotes the softmax function), just to present a standard loss for a deep learning model. Section 6: "As has been found, second-order-methods on full-batch training converge to the solution within a few iterations.": Reference? The method is compared to Nesterov accelerated gradient (NAG) for the data-fitting problem (fig 1) but not on the other tasks (figs 2 to 5). In the end, it was still unclear to me what the training consists of with the introduced exponentially decaying flows. I understand that a differential equation was formulated but then how is it solved iteratively? What is the cost of solving it? In Section 6, different algorithms are compared in terms of number of steps. What about the computational cost?