Title: UNDERSTANDING AND EXPLOITING THE LOW-RANK STRUCTURE OF DEEP NETWORKS Anonymous

Abstract: Training methods for deep networks are primarily variants on stochastic gradient descent. Techniques that use (approximate) second-order information are rarely used because of the computational cost and noise associated with those approaches in deep learning contexts. However, in this paper, we show how feedforward deep networks exhibit a low-rank derivative structure. This low-rank structure makes it possible to use second-order information without needing approximations and without incurring a significantly greater computational cost than gradient descent. To demonstrate this capability, we implement Cubic Regularization (CR) on a feedforward deep network with stochastic gradient descent and two of its variants. There, we use CR to calculate learning rates on a per-iteration basis while training on the MNIST and CIFAR-10 datasets. CR proved particularly successful in escaping plateau regions of the objective function. We also found that this approach requires less problem-specific information (e.g. an optimal initial learning rate) than other first-order methods in order to perform well.

1. INTRODUCTION:

1.1 GRADIENT-BASED OPTIMIZATION AND DEEP LEARNING Gradient-based optimization methods use derivative information to determine intelligent search directions when minimizing a continuous objective function. The steepest descent method is the most basic of these optimization techniques, but it is known to converge very slowly in ill-conditioned systems. Even outside of these cases, it still only has a linear rate of convergence. Newton's method is a more sophisticated approach -one that uses second-order derivative information, which allows the optimizer to model the error surface more accurately and thus take more efficient update steps. When it converges, it does so quadratically, but Newton's method also has limitations of its own. Firstly, it does not scale well: it can be very expensive to calculate, store, and invert the objective function Hessian. Secondly, the method may fail if the Hessian is indefinite or singular.

A variety of methods have been developed to try and appropriate the strengths of each approach while avoiding their weaknesses. The conjugate gradient method, for example, uses only firstorder information but uses the history of past steps taken to produce a better convergence rate than steepest descent. Quasi-Newton methods, on the other hand, approximate the Hessian (or its inverse) using first-order information and may enforce positive-definiteness on its approximation. Other approaches like trust region methods use second-order information without requiring convexity. For further information about gradient-based optimization, see Nocedal & Wright (2006) .

Deep learning (DL) provides a set of problems that can be tackled with gradient-based optimization methods, but it has a number of unique features and challenges. Firstly, DL problems can be extremely large, and storing the Hessian, or even a full matrix approximation thereto, is not feasible for such problems. Secondly, DL problems are often highly nonconvex. Thirdly, training deep networks via mini-batch sampling results in a stochastic optimization problem. Even if the necessary expectations can be calculated (in an unbiased way), the variance associated with the batch sample calculations produces noise, and this noise can make it more difficult to perform the optimization. Finally, deep networks consist of the composition of analytic functions whose forms are known. As such, we can calculate derivative information analytically via back-propagation (i.e. the chain rule).

1.2. TRAINING METHODS FOR DEEP LEARNING:

These special characteristics of DL have motivated researchers to develop training methods specifically designed to overcome the challenges with training a deep neural network. One such approach is layer-wise pretraining (Bengio et al., 2007) , where pretraining a neural network layer-by-layer encourages the weights to initialize close to a optimal minimum. Transfer learning (Yosinski et al., 2014) works by a similar mechanism, relying on knowledge gained through previous tasks to encourage nice training on a novel task. Outside of pretraining, a class of optimization algorithms have been specifically designed for training deep networks. The Adam, Adagrad, and Adamax set of algorithms provide examples of using history-dependent learning rate adjustment (Kingma & Ba, 2014) . Similarly, Nesterov momentum provides a method for leveraging history dependence in stochastic gradient descent (Sutskever et al., 2013) . One could possibly argue that these methods implicitly leverage second order information via their history dependence, but the stochastic nature of mini-batching prevents this from becoming explicit. Some researchers have sought to use second-order information explicitly to improve the training process. Most of these methods have used an approximation to the Hessian. For example, the L-BFGS method can estimate the Hessian (or its inverse) in a way that is feasible with respect to memory requirements; however, the noise associated with the sampling techniques can either overwhelm the estimation or require special modifications to the L-BFGS method to prevent it from diverging (Byrd et al., 2016) . There have been two primary ways to deal with this: subsampling (Byrd et al., 2016; Moritz et al., 2016) and mini-batch reuse (Schraudolph et al., 2007; Mokhtari & Ribeiro, 2014) . Subsampling involves updating the Hessian approximation every L iterations rather than every iteration, as would normally be done. Mini-batch reuse consists of using the same minibatch on subsequent iterations when calculating the difference in gradients between those two iterations. These approximate second-order methods typically have a computational cost that is higher than, though on the same order of, gradient descent, and that cost can be further reduced by using a smaller mini-batch for the Hessian approximation calculations than for the gradient calculation (Byrd et al., 2011) . There is also the question of bias: it is possible to produce unbiased low-rank Hessian approximations (Martens et al., 2012) , but if the Hessian is indefinite, then quasi-Newton methods will prefer biased estimates -ones that are positive definite. Other work has foregone these kinds of Hessian approximations in favor of using finite differences (Martens, 2010) .

1.3. CONTRIBUTIONS:

In this paper, we prove, by construction, that the first and second derivatives of feedforward deep learning networks exhibit a low-rank, outer product structure. This structure allows us to use and manipulate second-order derivative information, without requiring approximation, in a computationally feasible way. As an application of this low-rank structure, we implement Cubic Regularization (CR) to exploit Hessian information in calculating learning rates while training a feedforward deep network. Finally, we show that calculating learning rates in this fashion can improve existing training methods' ability to exit plateau regions during the training process.

2. THE LOW-RANK STRUCTURE OF DEEP NETWORK DERIVATIVES:

Second-order derivatives are not widely used in DL, and where they are used, they are typically estimated. These derivatives can be calculated analytically, but this is not often done because of the scalability constraints described in Section 1.1. If we write out the first and second derivatives, though, we can see that they have a low-rank structure to them -an outer product structure, in fact. When a matrix has low rank (or less than full rank), it means that the information contained in that matrix (or the operations performed by that matrix) can be fully represented without needing to know every entry of that matrix. An outer product structure is a special case of this, where an mxn matrix A can be fully represented by two vectors A = uv T . We can then calculate, store, and use secondorder derivatives exactly in an efficient manner by only dealing with the components needed to represent the full Hessians rather than dealing with those Hessians themselves. Doing this involves some extra calculations, but the storage costs are comparable to those of gradient calculations.

In this section, we will illustrate the low-rank structure for a feedforward network, of arbitrary depth and layer widths, consisting of ReLUs in the hidden layers and a softmax at the output layer. A feedforward network with arbitrary activation functions has somewhat more complicated derivative formulae, but those derivatives still exhibit a low-rank structure. That structure also does not depend on the form of the objective function or whether a softmax is used, and it is present for convolutional and recurrent layers as well. The complete derivations for these cases are given in Appendix B.

In our calculations, we make extensive use of index notation with the summation convention (Ivancevic & Ivancevic, 2007) . In index notation, a scalar has no indices (v), a vector has one index (v as v i or v i ), a matrix has two (V as V ij , V i j , or V ij ), and so on. The summation convention holds that repeated indices in a given expression are summed over unless otherwise indicated. For example, a T b = i a i b i = a i b i . The pair of indices being summed over will often consist of a superscript and a subscript; this is a bookkeeping technique used in differential geometry, but in this context, the subscripting or superscripting of indices will not indicate covariance or contravariance. We have also adapted index notation slightly to suit the structure of deep networks better: indices placed in brackets (e.g. the k in v (k),j ) are not summed over, even if repeated, unless explicitly indicated by a summation sign. A tensor convention that we will use, however, is the Kronecker delta: δ ij , δ i j , or δ ij . The Kronecker delta is the identity matrix represented in index notation: it is 1 for i = j and 0 otherwise. The summation convention can sometimes be employed to simplify expressions containing Kronecker deltas. For example, δ j i v i = v j and δ ij V jk = V ik . Let us consider a generic feedforward network with ReLU activation functions in n hidden layers, a softmax at the output layer, and categorical cross-entropy as the objective function (defined in more detail in Appendix B. The first derivatives, on a per-sample basis, for this deep network are

EQUATION

EQUATION

where f is the per-sample objective function, v (k),j is the vector output of layer k, u i j is the matrix of weights in the softmax, and p j = u j i v (n),i (which is the vector quantity evaluated by the softmax). For the full derivation, and the definition of the matrix quantity η (n,k),m i , see Appendix B. In calculating these expressions, we have deliberately left ∂f ∂p j unevaluated. This keeps the expression relatively simple, and programs like TensorFlow (Abadi et al., 2015) can easily calculate this for us. Leaving it in this form also preserves the generality of the expression -there is no low-rank structure contained in ∂f ∂p j , and the low-rank structure of the network as a whole is therefore shown to be independent of the objective function and whether or not a softmax is used. In fact, as long as Equation 13 holds, any sufficiently smooth function of p j may be used in place of a softmax without disrupting the low-rank structure. The one quantity that needs to be stored here is η (n,k),j i for k = 1, 2, . . . , n -1; it will be needed in the second derivative calculations. Note, however, that this is roughly the same size as the gradient itself.

We can now see the low-rank structure: ∂f ∂u i j is the outer product (or tensor product) of the vectors ∂f ∂p i and v (n),j , and ∂f ∂w (k),i j is the outer product of ∂f ∂p l u l m η

(n,k),m i (which ends up being a rank-1 tensor) and v (k-1),j . The index notation makes the outer product structure clear. It is important to note that this low-rank structure only exists for each sample -a weighted sum of low-rank matrices is not necessarily (and generally, will not be) low rank. In other words, even if the gradient of f is low rank, the gradient of the expectation, F = E [f ], will not be, because the gradient of F is the weighted sum of the gradients of f . The second-order objective function derivatives are then

EQUATION

EQUATION

∂ 2 f ∂w (k),i j ∂w (r),s t = ∂ 2 f ∂p l ∂p q u l m η (n,k),m i v (k-1),j u q a η (n,r),a s v (r-1),t + ∂f ∂p l u l m ×    η (n,r),m s η (r-1,k),t i v (k-1),j r > k 0 r = k η (n,k),m i η (k-1,r),j s v (r-1),t r < k (5)

Calculating all of these second derivatives requires the repeated use of ∂ 2 f ∂p 2 . Evaluating that Hessian is straightforward given knowledge of the activation functions and objective used in the network, and storing it is also likely not an issue as long as the number of categories is small relative to the number of weights. For example, consider a small network with 10 categories and 1000 weights. In such a case, ∂ 2 f ∂p 2 would only contain 100 entries -the gradient would be 10 times larger. We now find that we have to store η (n,k),i j values in order to calculate the derivatives. In ∂ 2 f ∂w 2 , we also end up needing η (r,k),i j for r = n. In a network with n hidden layers, we would then have n(n-1) 2 of the η (r,k),i j matrices to store. For n = 10, this would be 45, for n = 20, this would be 190, and so on. This aspect of the calculations does not seem to scale well, but in practice, it is relatively simple to work around. It is still necessary to store η

(n,k),i j , k < n, but η (r,k),i j

, r < n, only actually shows up in one place, and thus it is possible to calculate each η (r,k),i j matrix, use it, and discard it without needing to store it for future calculations. The key thing to note about these second derivatives is that they retain a low-rank structure -they are now tensor products (or the sums of tensor products) of matrices and vectors. For example,

EQUATION

a is = ∂ 2 f ∂p i ∂p l u l m η (n,k),m s (7)

With these expressions, it would be relatively straightforward to extract the diagonal of the Hessian and store or manipulate it as a vector. The rank of the weighted sum of low rank components (as occurs with mini-batch sampling) is generally larger than the rank of the summed components, however. As such, manipulating the entire Hessian may not be as computationally feasible; this will depend on how large the mini-batch size is relative to the number of weights. The low rank properties that we highlight here for the Hessian exist on a per-sample basis, as they did for the gradient, and therefore, the computational savings provided by this approach will be most salient when calculating scalar or vector quantities on a sample-by-sample basis and then taking a weighted sum of the results. In principle, we could calculate third derivatives, but the formulae would likely become unwieldy, and they may require memory usage significantly greater than that involved in storing gradient information. Second derivatives should suffice for now, but of course if a use arose for third derivatives, calculating them would be a real option. Thus far, we have not included bias terms. Including bias terms as trainable weights would increase the overall size of the gradient (by adding additional variables), but it would not change the overall low-rank structure. Using the calculations provided in Appendix B, it would not be difficult to produce the appropriate derivations.

3. CUBIC REGULARIZATION IN DEEP LEARNING:

Cubic Regularization (CR) is a trust region method that uses a cubic model of the objective function:

f (x) ≈ m j (s j ) = f (x j ) + s T j ∂f ∂x + 1 2 s T j H j s j + 1 6 σ j s j 3 (8)

at the j-th iteration, where H j is the objective function Hessian and s j = xx j . The cubic term makes it possible to use information in the Hessian without requiring convexity, and the weight σ j on that cubic term can have its own update scheme (based on how well m (s j ) approximates f (Kohler & Lucchi, 2017)). Solving for an optimal s j value then involves finding the root of a univariate nonlinear equation (Nesterov & Polyak, 2006) . CR is not commonly used in deep learning; we have seen only one example of CR applied to machine learning (Kohler & Lucchi, 2017) and no examples with deep learning. This is likely the case because of two computationally expensive operations: calculating the Hessian and solving for s j . We can overcome the first by using the lowrank properties described above. The second is more challenging, but we can bypass it by using CR to calculate a step length (i.e. the learning rate) for a given search direction rather than calculating the search direction itself.

3.1. IMPLEMENTATION:

Our approach in this paper is to use CR as a metamethod -a technique that sits on top of existing training algorithms. The algorithm calculates a search direction, and then CR calculates a learning rate for that search direction. For a general iterative optimization process, this would look like x j+1 = x j + α j g j , where g j is the search direction (which need not be normalized), α j is the learning rate, and the subscript refers to the iteration. With the search direction fixed, m would then be a cubic function of α at each iteration. Solving ∂m ∂α = 0 as a quadratic equation in α then yields

α = -g T Hg ± (g T Hg) 2 -2 (g T ∇f ) σ g 3 σ g 3 (9)

If we assume that g T ∇f < 0 (i.e. g is a descent direction), then α is guaranteed to be real. Continuing under that assumption, of the two possible α values, we choose the one guaranteed to be positive. The sampling involved in mini-batch training means that there are a number of possible ways to get a final α j g j result. One option would be to calculate E [α j g j ]. This would involve calculating an α value with respect to the search direction produced by each sample point and then averaging the product αg over all of the sample points. Doing this should produce an unbiased estimate of α j g j , but in practice, we found that this approach resulted in a great deal sampling noise and thus was not effective. The second approach would calculate

E [α j ] × E [g j ].

To do this, we would calculate an α value with respect to the search direction produced by each sample point, as in the first option, calculate an average α value, and multiply the overall search direction by that average. This approach, too, suffered from excessive noise. In the interest of reducing noise and increasing simplicity, we chose a third option: once the step direction had been determined, we considered that fixed, took the average of g T Hg and g T ∇f over all of the sample points to produce m (α) and then solved for a single α j value. This approach was the most effective of the three.

3.2. COMPUTATIONAL RESULTS:

To test CR computationally, we created deep feedforward networks using ReLU activations in the hidden layers, softmax in the output layer, and categorical cross-entropy as the error function; we then trained them on the MNIST (LeCun et al., 1998) and CIFAR-10 (Krizhevsky & Hinton, 2009) data sets. This paper shows results from networks with 12 hidden layers, each 128 nodes wide.

For the purposes of this paper, we treat network training strictly as an optimization process, and thus we are not interested in network performance measures such as accuracy and validation error -the sole consideration is minimizing the error function presented to the network. As we consider that minimization progress, we will also focus on optimization iteration rather than wall clock time: the former indicates the behaviour of the algorithm itself, whereas the latter is strongly dependent upon implementation (which we do not want to address at this juncture). Overall computational cost per iteration matters, and we will discuss it, but it will not be our primary interest. Further implementation details are found in Appendix A.

Figure 1 shows an example of CR (applied on top of SGD). In this case, using CR provided little to no benefit. The average learning rate with CR was around 0.05 (a moving average with a period of 100 is shown in green on the learning rate plot both here and in the rest of the paper), which was close to our initial choice of learning rate. This suggests that 0.02 was a good choice of learning rate. Another reason the results were similar, though, is that the optimization process did not run into any plateaus. We would expect CR to provide the greatest benefit when the optimization gets stuck on a plateau -having information about the objective function curvature would enable the algorithm to increase the learning rate while on the plateau and then return it to a more typical value once it leaves the plateau. To test this, we deliberately initialized our weights so that they lay on a plateau: the objective function is very flat near the origin, and we found that setting the network weights to random values uniformly sampled between 0.1 and -0.1 was sufficient. Figure 2 shows the results of SGD with and without CR when stuck on a plateau. There, we see a hundred-fold increase in the learning rate while the optimization is on the plateau, but this rate drops rapidly as the optimization exits the plateau, and once it returns to a more normal descent, the learning rate also returns to an average of about 0.05 as before. The CR calculation enables the training process to recognize the flat space and take significantly larger steps as a result. Applying CR to SGD when training on CIFAR-10 (Figure 3 ) produced results similar to those seen on MNIST.

We then considered if this behaviour would hold true on other training algorithms: we employed CR with Adagrad (Duchi et al., 2011) and Adadelta (Zeiler, 2012)on MNIST. The results were similar. CR did not provide a meaningful difference when the algorithms performed well, but when those algorithms were stuck on plateaus, CR increased the learning rate and caused the algorithms to exit the plateau more quickly than they otherwise would have (Figures 4 and 5 ). The relative magnitudes of those increases were smaller than for SGD, but Adagrad and Adadelta already incorporate some adaptive learning rate behaviour, and good choices for the initial learning rate varied significantly from algorithm to algorithm. We also used a larger value for σ to account for the increased variability due to those algorithms' adaptive nature. The result with Adadelta showed some interesting learning rate changes: the learning rate calculated by CR dropped steadily as the algorithm exited the plateau, but it jumped again around iteration 1200 as it apparently found itself in a flat region of space. 

4. DISCUSSION:

We see this CR approach as an addition to, not a replacement for, existing training methods. It could potentially replace existing methods, but it does not have to in order to be used. Because of the low-rank structure of the Hessian, we can use CR to supplement existing optimizers that do not explicitly leverage second order information. The CR technique used here is most useful when the optimization is stuck on a plateau prior to convergence: CR makes it possible to determine whether the optimization has converged (perhaps to a local minimum) or is simply bogged down in a flat region. It may eventually be possible to calculate a search direction as well as a step length, which would likely be a significant advancement, but this would be a completely separate algorithm.

We found that applying CR to Adagrad and Adadelta provided the same kinds of improvements that applying CR to SGD did. However, using CR with Adam (Kingma & Ba, 2014) did not provide gains as it did with the other methods. Adam generally demonstrates a greater degree of adaptivity than Adagrad or Adadelta; in our experiments, we found that Adam was better than Adagrad or Adadelta in escaping the plateau region. We suspect that trying to overlay an additional calculated learning rate on top of the variable-specific learning rate produced by Adam may create interference in both sets of learning rate calculations. Analyzing each algorithm's update scheme in conjunction with the CR calculations could provide insight into the nature and extent of this interference, and provide ways to further improve both algorithms. In future work, though, it would not be difficult to adapt the CR approach to calculate layer-or variable-specific learning rates, and doing that could address this problem. Calculating a variable-specific learning rate would essentially involve rescaling each variable's step by the corresponding diagonal entry in the Hessian; calculating a layer-specific learning rate would involve rescaling the step of each variable in that layer by some measure of the block diagonal component of the Hessian corresponding to those variables. The calculations for variable-specific learning rates with CR are given in Appendix B.

There are two aspects of the computational cost to consider in evaluating the use of CR. The first aspect is storage cost. In this regard, the second-order calculations are relatively inexpensive (comparable to storing gradient information). The second aspect is the number of operations, and the second-order calculations circumvent the storage issue by increasing the number of operations. The number of matrix multiplications involved in calculating the components of Equation 9, for example, scales quadratically with the number of layers (see the derivations in Appendix B). Although the number of matrix multiplications will not change with an increase in width, the cost of naïve matrix multiplication scales cubically with matrix size. That being said, these calculations are parallelizable and as such, the effect of the computation cost will be implementation-dependent.

A significant distinction between CR and methods like SGD has to do with the degree of knowledge about the problem required prior to optimization. SGD requires an initial learning rate and (usually) a learning rate decay scheme; an optimal value for the former can be very problem-dependent and may be different for other algorithms when applied to the same problem. For CR, it is necessary to specify σ, but optimization performance is relatively insensitive to this -order of magnitude estimates seem to be sufficient -and varying σ has a stronger affect on the variability of the learning rate than it does on the magnitude (though it does affect both). If the space is very curved, the choice of σ matters little because the step size determination is dominated by the curvature, and if the space if flat, it bounds the step length. It is also possible to employ an adaptive approach for updating σ (Kohler & Lucchi, 2017), but we did not pursue that here. Essentially, using CR is roughly equivalent to using the optimal learning rate (for SGD).

5. CONCLUSIONS:

In this paper, we showed that feedforward networks exhibit a low-rank derivative structure. We demonstrate that this structure provides a way to represent the Hessian efficiently; we can exploit this structure to obtain higher-order derivative information at relatively low computational cost and without massive storage requirements. We then used second-order derivative information to implement CR in calculating a learning rate when supplied with a search direction. The CR method has a higher per-iteration cost than SGD, for example, but it is also highly parallelizable. When SGD converged well, CR showed comparable optimization performance (on a per-iteration basis), but the adaptive learning rate that CR provided proved to be capable of driving the optimization away from plateaus that SGD would stagnate on. The results were similar with Adagrad and Adadelta, though not with Adam. CR also required less problem-specific knowledge (such as an optimal initial learning rate) to perform well. At this point, we see it as a valuable technique that can be incorporated into existing methods, but there is room for further work on exploiting the low-rank derivative structure to enable CR to calculate search directions as well as step sizes. 

A EXPERIMENTAL PROCEDURE:

Starting at a point far from the origin resulted in extremely large derivative and curvature values (not to mention extremely large objective function values), and this could sometimes cause difficulties for the CR method. This was easy to solve by choosing an initialization point relatively near the origin; choosing an initialization relatively near the origin also provided a significantly better initial objective function value. We initialized the networks' weights to random values between an upper and lower bound: to induce plateau effects, we set, the bounds to ±0.1, otherwise, we set them to ±0.2.

All of the networks used a mini-batch size of 32 and were implemented in TensorFlow (Abadi et al., 2015) . The initial learning rate varied with network size; we chose learning rates that were large and reasonable but perhaps not optimal, and for optimization algorithms with other parameters governing the optimization, we used the default TensorFlow values for those parameters. For the learning rate decay, we used an exponential decay with a decay rate of 0.95 per 100 iterations. The σ value used is specified along with the initial learning rate for each network's results. This value was also not optimized but was instead set to a reasonable power of 10. Vector of inputs for a single sample v (k),j

Vector output of layer k w (k),j i

Matrix of weights for layer k A (•):

Activation function u j i

Matrix of output layer weights p j:

Vector of intermediate variables for the output layer ŷl

Vector of outputs for a single sample y l

Vector of labels for a single sample f Scalar objective function value for a single sample F Scalar objective function Equations 10-16 define a generic feedforward network with ReLU activation functions in the hidden layers, n hidden layers, a softmax at the output layer, and categorical cross-entropy as the objective function.

EQUATION

EQUATION

EQUATION

EQUATION

The relevant first derivatives for this deep network are

EQUATION

)

where there is no summation over j in Equation 18. We now define several intermediate quantities to simplify the derivation process:

EQUATION

EQUATION

where there is no summation over j in Equations 19 and 20. We can now complete our calculations of the first derivatives.

EQUATION

EQUATION

EQUATION

We then start our second derivative calculations by considering some intermediate quantities:

A (z) = 0 (29) ∂α (n,k),m q ∂w (r),s t = α (n,r),m a A w (r),a p v (r-1),p δ a s δ t b α (r-1,k),b q = η (n,r),m s α (r-1,k),t q (30) ∂γ (k),j s ∂ (•) = 0 (31) ∂η (n,k),m i ∂w (r),s t = ∂α (n,k),m q ∂w (r),s t γ (k),q i (32) ∂ 2 v (n),m ∂w (k),i j ∂w (r),s t =        ∂α (n,k),m q ∂w (r),s t γ (k),q i v (k-1),j r > k 0 r = k η (n,k),m i ∂v (k-1),j ∂w (l),s t r < k =    η (n,r),m s η (r-1,k),t i v (k-1),j r > k 0 r = k η (n,k),m i η (k-1,r),j s v (r-1),t r < k (33)

The second derivative of the ReLU vanishes, which simplifies the second derivative calculations significantly. Technically, the second derivative is undefined at the origin, but the singularity is removable, and thus we can define the second derivative to be 0 at the origin. We can then calculate the second-order objective function derivatives:

EQUATION

EQUATION

∂ 2 f ∂w (k),i j ∂w (r),s t = ∂ 2 f ∂p l ∂p q ∂p l ∂v (n),m ∂v (n),m ∂w (k),i j ∂p q ∂v (n),a ∂v (n),a ∂w (r),s t + ∂f ∂p l ∂p l ∂v (n),m ∂ 2 v (n),m ∂w (k),i j ∂w (r),s t = ∂ 2 f ∂p l ∂p q u l m η (n,k),m i v (k-1),j u q a η (n,r),a s v (r-1),t + ∂f ∂p l u l m ×    η (n,r),m s η (r-1,k),t i v (k-1),j r > k 0 r = k η (n,k),m i η (k-1,r),j s v (r-1),t r < k (36)

To use CR, we calculate α as

EQUATION

For a given iteration for the deep network described above (dropping the subscript j's so as not to interfere with the index notation), the quantities in this equation are

EQUATION

EQUATION

With these formulae in hand, and using the generalized binomial theorem ( √ 1 + ≈ 1 + /2 for 1), we can consider what happens to α in limiting cases:

EQUATION

The weight update scheme for a single learning rate at each iteration is

EQUATION

We could instead consider a weight-specific learning rate. If we assume that the baseline g vector for each variable is a unit step in the direction of that variable, and we ignore superscripts indicating the iteration number, the calculations for variable-specific learning rates α (u),i j for u i j and α (w,k),i j are as shown in Table 2 : 

α (u),i j ∂F ∂p i v (n),j ∂ 2 F ∂p i ∂p s v (n),j v (n),t 1 α (w,k),i j ∂F ∂p l u l m η (n,k),m i v (k-1),j ∂ 2 F ∂p l ∂p q u l m η (n,k),m i v (k-1),j u q a η (n,k),a s v (k-1),t 1 B.

2. CONVOLUTIONAL AND RECURRENT LAYERS:

Convolutional and recurrent layers preserve the low-rank derivative structure of the fully connected feedforward layers considered above, and we will show this in the following sections. Because we are only considering a single layer of each, we calculate the derivatives of the layer outputs with respect to the layer inputs -in a larger network, those derivatives will be necessary for calculating total derivatives via back-propagation.

B.2.1 CONVOLUTIONAL LAYER:

We can define a convolutional layer as

EQUATION

where x i j is the layer input, σ is the vertical stride, τ is the horizontal stride, A is the activation function, and v s t is the layer output. A convolutional structure can make the expressions somewhat complicated when expressed in index notation, but we can simplify matters by using the simplification z sl tk = x σs+l-1 τ t+k-1 . The layer definition is then

v s t = A z sl tk w l k (46)

The derivatives of the convolutional layer are with no summation over s and t in any of the expressions above. Using the simplification with z sl tk makes it significantly easier to see the low rank structure in these derivatives, but that structure is still noticeable without the simplification.

EQUATION

The conditional form of the expressions is more complicated, but it is also possible to see how the derivatives relate to w i j and submatrices of x i j .

B.2.2 RECURSIVE LAYER:

We can define our recursive layer as

EQUATION

where t indicates the number of times that the recursion has been looped through. If we inspect this carefully, we can actually see that this is almost identical to the hidden layers of the feedforward network: they are identical if we stipulate that the weights of the feedforward network are identical at each layer (i.e. w (k),i j = w i j ∀k) and if we treat the recursive loops like layers. This observation allows us to reuse some of our previous derivations. Primarily, we will use the fact that

∂ (•) ∂w i j = k ∂ (•) ∂w (k),m p ∂w (k),i j ∂w i j = k ∂ (•) ∂w (k),i j (60)

The first-order derivatives are then

∂v m (t) ∂w s r = k η (t,k),m s v r (t-1) (61) ∂v m (t) ∂x i = η (t,1),m s w s i (62)

If A is a ReLU, then the second derivatives are relatively simple

EQUATION

)

If A is not a ReLU, then we would use the results in the next section to calculate the second derivatives. Regardless of the exact form of A, though, we retain the low-rank structure as long as A is an entry-wise function of its arguments.

B.3 DEEP NETWORK WITH GENERAL ACTIVATION FUNCTIONS:

For a deep network with general entry-wise functions, the first derivatives are all identical to the derivations given in Section 2 save that A will be different. Before calculating second-order derivatives, though, we do some preliminary calculations that were not necessary before because A = 0 for ReLUs. First, we calculate the derivatives of γ (l),m r :

EQUATION

(l),m r ∂w (q),s t =    0 l < q λ (l),m rs v (l-1),t l = q λ (l),m rp w (l),p a η (l-1,q),a s v (q-1),t l > q (68)

where there is no summation over m in any of these equations. Next, we calculate the derivatives of α

(k,l),j m : α (k,l),j m = α (k,r),j b α (r,r-1),b a α (r-1,l),a m (69) α (r,r-1),b a = β (r),b a = γ (r),b s w (r),s a (70) ∂β (r),b a ∂w (q),s t =    0 r < q λ (r),b ps w (r),p a v (r-1),t + γ (r),b s δ t a r = q λ (r),b dp w (r),d a w (r),p m η (r-1,q),m s v (q-1),t r > q (71) ∂α (k,l),j m ∂w (q),s t = k r=l+1 α (k,r),j b ∂β (r),b a ∂w (q),s t α (r-1,l),a m =                          0 k < q k r=q+1 α (k,r),j b λ (r),b dp w (r),d a w (r),p c η (r-1,q),c s v (q-1),t α (r-1,l),a m +α (k,q),j b λ (q),b

ps w (q),p a v (q-1),t + γ (q),p a v (q-1),t + γ (k-1,q),a s v (q-1),t q < k (74) η (r-1,q),c s v (q-1),t η (r-1,k),a l +α (n,q),j b λ (q),b ps w (r),p a v (q-1),t η (q-1,k),a l + η

=                                     

=                                 

(n,q),j s η (q-1,k),t l q > k (k-1,q),a s v (q-1),t q < k (75)

∂ 2 v (n),j ∂w (k),l i ∂w (q),s t =      ∂η (n,k),j l ∂w (q),s t v (k-1),i q ≥ k ∂η (n,k),j l ∂w (q),s t v (k-1),i + η

(n,k),j l ∂v (k-1),i ∂w (q),s t q < k (76) η (r-1,q),c s v (q-1),t η (r-1,k),a l v (k-1),i +α (n,q),j b λ (q),b ps w (q),p a v (q-1),t η (q-1,k),a l v (k-1),i + η

=                                     

(n,q),j s η (q-1,k),t l v (k-1),i η (k-1,q),a s v (q-1),t v (k-1),i + η

(n,k),j l η (k-1,q),i s v (q-1),t q < k (77)

The second-order derivatives of the objective function are then

EQUATION

∂ 2 F ∂u i j ∂w ∂ 2 v (n),j ∂w (k),l i ∂w (q),s t = ∂ 2 F ∂p a ∂p b u a j η

(n,k),j l v (k-1),i u b c η (n,q),c s v (q-1),t + ∂F ∂p a u a j ∂ 2 v (n),j ∂w (k),l i ∂w (q),s t (80)

:



Figures and Tables:

Figure 1: Cubic Regularization (CR) applied to Stochastic Gradient Descent (SGD); initial learning rate = 0.01, σ = 100

Figure 3: Cubic Regularization (CR) applied to Stochastic Gradient Descent (SGD) on the CIFAR-10 Dataset; initial learning rate = 0.01, σ = 1000

Figure 4: Cubic Regularization (CR) applied to Adagrad; initial learning rate = 0.1, σ = 1000

σs + m -1 and j = τ t + p -

1,k),c s v (k-1),t η (r-1,k),a l +α (n,k),j m λ (k),m lp w (k),p a η

Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson. How transferable are features in deep neural networks? In Advances in neural information processing systems, pp. 3320-3328, 2014.

Table content: <table><tr><td>Matthew D Zeiler. Adadelta: an adaptive learning rate method. arXiv preprint arXiv:1212.5701,</td></tr><tr><td>2012.</td></tr></table>

Variable-Specific Learning Rate Calculations

Table content: <table><tr><td>Learning Rate</td><td>g t ∇f</td><td>g T Hg</td><td>g</td><td>3</td></tr></table>

Peer review (iclr2018):
 This paper proposes to set a global step size gradient-based optimization algorithms such as SGD and Adam using second order information. Instead of using second-order information to compute the update directly (as is done in e.g. Newton method), it is used to estimate the change of the objective function in a pre-computed direction. This is computationally much cheaper than full Newton because (a) the Hessian does not need to be inverted (b) vector-Hessian multiplication is only O(#parameters) for a single sample. There are many issues. ### runtime and computational issues ### Firstly, the paper does not clearly specify the algorithm it espouses. It states: "once the step direction had been determined, we considered that fixed, took the average of gT Hg and gT ∇f over all of the sample points to produce m (α) and then solved for a single αj value" You should present pseudo-code for this computation and not leave the reader to determine the detailed order of computation for himself. As it stands, it is not only difficult for the reader to infer these details, but also laborious to determine the computational cost per iteration on some network the reader might wish to apply your algorithm to. Since the paper discusses the computational cost of CR only in vague terms, you should at least provide pseudo-code. Specifically, consider equation (80) at the very end of the appendix and consider the very last term in that equation. It contains d^2v/dwdw. This is a "heavy" term containing the second derivative of the last hidden layer with respect to weights. You do not specify how you compute this term or quantities involving this term. In a ReLU network, this term is zero due to local linearity, but since you claim that your algorithm is applicable to general networks, this term needs to be analyzed further. While the precise algorithm you suggest is unclear, it*s purpose is also unclear. You only use the Hessian to compute the g^THg terms, i.e. for Hessian-vector multiplication. But it is well-known that Hessian-vector multiplication is "relatively cheap" in deep networks and this fact has been used for several algorithms, e.g. http://www.iro.umontreal.ca/~lisa/pointeurs/ECML2011_CAE.pdf and https://arxiv.org/pdf/1706.04859.pdf. How is your method for computing g^THg different and why is it superior? Also note that the low-rank structure of deep gradients is well-known and not a contribution of this paper. See e.g. https://www.usenix.org/system/files/conference/atc17/atc17-zhang.pdf ### Experiments ### The experiments are very weak. In a network where weights are initialized to sensible values, your algorithm is shown not to improve upon straight SGD. You only demonstrate superior results when the weights are badly initialized. However, there are a very large number of techniques already that avoid the "SGD on ReLU network with bad initial weights" problem. The most well-known are batch normalization, He initialization and Adam but there are many others. I don*t think it*s a stretch to consider that problem "solved". Your algorithm is not shown to address any other problems, but what*s worse is that it doesn*t even seem to address that problem well. While your learning curves are better than straight SGD, I suspect they are well below the respective curves for He init or batchnorm. In any case, you would need to compare your algorithm against these state-of-the-art methods if your goal is to overcome bad initializations. Also, in appendix A, you state that CR can*t even address weights that were initialized to values that are too large. You claim that your algorithm helps with "overcoming plateaus". While I have heard the claim that deep network optimization suffers from intermediate plateaus before, I have not seen a paper studying / demonstrating this behavior. I suggest you cite several papers that do this and then replicate the plateau situations that arose in those papers and show that CR overcomes them, instead of resorting to a platenau situation that is essentially artificially induced by intentionally bad hyperparameter choices. I do not understand why your initial learning rate for SGD in figures 2 and 3 (0.02 and 0.01 respectively) differ so much from the initial learning rate under CR. Aren*t you trying to show that CR can find the "correct" learning rate? Wouldn*t that suggest that initial learning rate for SGD should be comparable to the early learning rates chosen by CR? Wouldn*t that suggest you should start SGD with a learning rate of around 2 and 0.35 respectively? Since you are annealing the learning rate for SGD, it*s going to decline and get close to 0.02 / 0.01 anyway at some point. While this may not be as good as CR or indeed batchnorm or Adam, the blue constant curve you are showing does not seem to be a fair representation of what SGD can do. You say the minibatch size is 32. For MNIST, this means that 1 epoch is around 1500 iterations. That means your plots only show the first epoch of training. But MNIST does not converge in 1 epoch. You should show the error curve until convergence is reached. Same for CIFAR. "we are not interested in network performance measures such as accuracy and validation error" I strongly suspect your readers may be interested in those things. You should show validation classification error or at least training classification error in addition to cross-entropy error. "we will also focus on optimization iteration rather than wall clock time" Again, your readers care more about the latter. You need to show either error curves by clock time or the total time to convergence or supplement your iteration-based graphs with a detailed discussion of how long an iteration takes. The scope of the experiments is limited because only a single network architecture is considered, and it is not a state-of-the art architecture (no convolution, no normalization mechanism, no skip connections). You state that you ran experiments on Adam, Adadelta and Adagrad, but you do not show the Adam results. You say in the text that they were the least favorable for CR. This suggests that you omitted the detailed results because they were unfavorable to you. This is, of course, unacceptable! ### (Un)suitability of ReLU for second-order analysis ### You claim to use second-order information over the network to set the step size. Unfortuantely, ReLU networks do not have second-order information! They are locally linear. All their nonlinearity is contained in non-differentiable region boundaries. While this may lead to the Hessian being cheaper to compute, it means it is not representative of the actual behavior of the network. In fact, the only second-order information that is brought to bear in your experiments is the second-order information of the error function. I am not saying that this particular second-order information could not be useful, but you need to make a distinction in your paper between network second-order info and error function second-order info and make explicit that you only use the former in your experiments. As far as I know, most second-order papers use either tanh or a smoothed ReLU (such as the smoothed hinge used recently by Koh & Liang (https://arxiv.org/pdf/1703.04730.pdf)) for experiments to overcome the local linearity. ### The .sigma hyperparameter ### You claim that .sigma is not as important / hard to set as .alpha in SGD or Adam. You state: "We also found that this ap- proach requires less problem-specific information (e.g. an optimal initial learning rate) than other first-order methods in order to perform well." You have not provided sufficient evidence for this claim. You say that .sigma can be chosen by considering powers of 10. In many networks, choosing .alpha by considering powers of 10 is sufficient! Even if powers of 2 are considered for .alpha, this would reduce the search effort only by factor log_2(10). Also, what if the range of .sigma values that need to be considered is larger than the range of .alpha values? Then setting .sigma would take more effort. You do not give precise protocols how you set .sigma and how you set .alpha for non-CR algorithms. This should be clearly specified in Appendix A as it is central to your argument of easing hyperparameter search. ### Minor points ### - Your introduction could benefit from a few more citations - "The rank of the weighted sum of low rank components (as occurs with mini-batch sampling) is generally larger than the rank of the summed components, however." I don*t understand this. Every sum can be viewed as a weighted sum and vice versa. - Equation (8) could be motivated a bit better. I know it derives from Taylor*s theorem, but it might be good to discuss how Taylor*s theorem (and its assumptions) relate to deep networks. - why the name "cubic regularization"? shouldn*t it be something like "quadratic step size tuning"? . . . The reason I am giving a 2 instead of a 1 is because the core idea behind the algorithm given seems to me to have potential, but the execution is sorely lacking. A final suggestion: You advertise as one of your algorithms upsides that it uses exact Hessian information. Howwever, since you only care about the scale of the second-order term and not its direction, I suspect exact calculation is far from necessary and you could get away with very cheap approximations, using for example techniques such as mean field analysis (e.g. http://papers.nips.cc/paper/6322-exponential-expressivity-in-deep-neural-networks-through-transient-chaos.pdf).