{'fraction_non_alphanumeric': 0.03607979184735473, 'fraction_numerical': 0.014277136566815665, 'mean_word_length': 4.280490382582963, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Stochastic Bundle Method for Interpolating Networks

We propose a novel method for training deep neural networks that are capable of interpolation, that is, driving the empirical loss to zero. At each iteration, our method constructs a stochastic approximation of the learning objective. The approximation, known as a bundle, is a pointwise maximum of linear functions. Our bundle contains a constant function that lower bounds the empirical loss. This enables us to compute an automatic adaptive learning rate, thereby providing an accurate solution. In addition, our bundle includes linear approximations computed at the current iterate and other linear estimates of the DNN parameters. The use of these additional approximations makes our method significantly more robust to its hyperparameters. Based on its desirable empirical properties, we term our method Bundle Optimisation for Robust and Accurate Training (BORAT). In order to operationalise BORAT, we design a novel algorithm for optimising the bundle approximation efficiently at each iteration. We establish the theoretical convergence of BORAT in both convex and non-convex settings. Using standard publicly available data sets, we provide a thorough comparison of BORAT to other single hyperparameter optimisation algorithms. Our experiments demonstrate BORAT matches the state-of-the-art generalisation performance for these methods and is the most robust.


Introduction
Training a deep neural network (DNN) is a challenging optimization problem: it involves minimizing the average of many high-dimensional non-convex loss functions. In practice, the main algorithms utilised are Stochastic Gradient Descent (SGD) (Robbins and Monro, 1951) and adaptive gradient methods such as AdaGrad (Duchi et al., 2011) or Adam (Kingma and Welling, 2014). It has been observed that SGD tends to provide better generalization performance than adaptive gradient methods (Wilson et al., 2017). However, the downside of SGD is that it requires the manual design of a learning-rate schedule. Many forms of schedule have been proposed in the literature, including piece wise constant (Huang et al., 2017), geometrically decreasing (Szegedy et al., 2015) and warm starts with cosine annealing (Loshchilov and Hutter, 2017). Examples of these schemes are plotted in Figure 1. Consequently, practitioners who wish to use SGD in a novel setting need to select a learning-rate schedule for their learning task. To that end, they first need to choose the parameterization of the schedule (e.g. picking one of the examples given above). Then they need to tune the corresponding parameters to get good empirical performance. This typically results in a cross-validation that searches over many critical and sensitive hyper-parameters. For example, a piece wise linear scheme requires a starting learning rate value, a decay factor and a list or metric to determine at which points in training to decay the learning rate. Due to the high dimensionality of this search space performing a grid search can mean training a large number of models. This number increases exponentially in combination with other hyperparameters such as regularisation and batch size. Thus, finding an SGD learning rate schedule that produces strong generalisation performance for a new task is time and computationally demanding, often requiring hundreds of GPU hours. In this work, we alleviate this issue by presenting a family of algorithms that achieve comparable generalisation performance to SGD with a highly refined learning rate schedule, while requiring far less tuning of hyperparameters. This in turn leads to a reduction in the time, cost and energy required when finding a highly accurate network for a new task.
In more detail, we present a novel family of proximal algorithms for the optimisation of DNNs that are capable of interpolation. A DNN is said to interpolate a data set if it has the ability to simultaneously drive the loss of all the training samples in a data set to their minimum value. Thus a lower bound of the objective function is known: for instance, it is close to zero for a model trained with the cross-entropy loss. Our algorithms approximate the loss within each proximal problem by a bundle, that is, a point-wise maximum over linear functions. By including the interpolation lower bound within this bundle, we obtain the following two modelling benefits: (i) our model more closely mimics the true loss function than existing baselines like SGD, and; (ii) the learning rate gets automatically re-scaled at each iteration of the algorithm, thereby providing accurate updates. By increasing the number of linear approximations in the bundle, the true loss is better modelled. As an upshot, we obtain additional stability to optimisation and task-specific hyperparameters. Based on its highly desirable empirical properties, we term these methods Bundle optimization for Robust and Accurate Training (BORAT).
All the variants of BORAT use a single learning rate hyperparameter that requires minimal tuning. In particular, the learning rate hyperparameter is kept constant throughout the training procedure, unlike the learning rate of SGD that needs to be decayed for good generalization. The BORAT family of algorithms obtain state-of-the-art empirical performance for single hyperparameter training of neural networks.

Contributions
• We design a family of adaptive algorithms that have a single hyperparameter that does not need any decaying schedule. In contrast, the related APROX (Asi and Duchi, 2019) and L4 (Rolinek and Martius, 2018) use respectively two and four hyperparameters for their learning-rate.
• For the deep learning setting we establish a link between stochastic optimisation methods with adaptive learning rates and proximal bundle methods.
• We provide convergence rates in various stochastic convex settings and for a class of non-convex problems.
• We derive a novel algorithm for solving small quadratic programs where the constraints define the probability simplex. This algorithm permits a parallel implementation allowing for efficient solution on modern hardware.
• We empirically demonstrate the increased stability to hyperparameters when increasing the bundle size. We show this on the CIFAR100 and Tiny ImageNet data sets.
• We achieve state-of-the-art results for learning a differentiable neural computer; training variants of residual networks on the SVHN and CIFAR data sets, and training a Bi-LSTM on the Stanford Natural Language Inference data set.
A preliminary version of this work appeared in the proceedings of ICML 2020 (Berrada et al., 2019b). While this previous work has considered bundles of size 2 resulting in the ALI-G algorithm detailed with in Section 3.3. This article significantly differs from the previous work by (i) considering bundles of size greater than 2; (ii) introducing an novel algorithm to compute the exact solution of each bundle; (iii) investigating the robustness towards hyperparameters; and (iv) showing how the use of large bundles permits easy application of BORAT to challenging non-smooth losses.

Preliminaries
Loss Function. We consider a supervised learning task where the model is parameterized by w ∈ R d . Usually, the objective function can be expressed as an expectation over z ∈ Z, a random variable indexing the samples of the training set: where each ℓ z is the loss function associated with the sample z. We assume that each ℓ z admits a known lower bound B. For simplicity we will often assume this lower bound is 0, which is the case for the large majority of loss functions used in machine learning. For instance, suppose that the model is a deep neural network with weights w performing classification. Then for each sample z, ℓ z (w) can represent the cross-entropy loss, which is always non-negative. Other non-negative loss functions include the structured or multi-class hinge loss, and the ℓ 1 or ℓ 2 loss functions for regression. Note that it is always possible to subtract a non-zero lower bound B from the loss function to define a new equivalent problem that satisfies the aforementioned assumption.
Interpolation. In this work we consider DNNs that can interpolate the data. Formally, we assume: where ǫ is a tolerance on the amount of error in the interpolation assumption. We will often want to make reference to the case when ǫ = 0. Following previous work (Ma et al., 2018b) we will refer to this setting as perfect interpolation. The interpolation property is satisfied in many practical cases, since modern neural networks are typically trained in an over-parameterized regime where the parameters of the model far exceed the size of the training data (Li et al., 2020). Additionally most modern DNN architectures can be easily increased in size and depth, allowing them to interpolate all but the largest data sets. Note, the data has to be labelled consistently for this to be possible. For instance it is impossible to interpolate a data set with two instances of the same image that have two different labels.
Regularisation. It is often desirable to encourage generalisation by the addition of a regularisation function φ(w) to the objective. Typical choices for φ include λ w 2 and λ w 1 where λ governs the strength of the regularisation. However, in this work we incorporate such regularisation as a constraint on the feasible domain: Ω = w ∈ R d : φ(w) ≤ r for some value of r. In the deep learning setting, this will allow us to assume that the objective function can be driven close to zero without unrealistic assumptions about the value of the regularisation term for the final set of parameters. Our framework can handle any constraint set Ω on which Euclidean projections are computationally efficient. This includes the feasible set induced by ℓ 2 regularization: Ω = w ∈ R d : w 2 2 ≤ r , for which the projection is given by a simple rescaling of w. Finally, note that if we do not wish to use any regularization, we define Ω = R d and the corresponding projection is the identity.
Problem Formulation. The learning task can be expressed as the problem (P) of finding a feasible vector of parameters w ⋆ ∈ Ω that minimizes f : We refer to the minimum value of f over Ω as f ⋆ : f ⋆ min w∈Ω f (w).
Proximal Perspective of Projected Stochastic Gradient Descent. In order to best introduce BORAT, we first detail the proximal interpretation of projected stochastic gradient descent (PSGD). The PSGD algorithm can be seen as solving a sequence of proximal problems. Within each proximal problem, a minimisation is performed over an approximate local model of the loss. This approximation is the first order Taylor's expansion of ℓ zt around the current iterate and a proximal term. At time step t, the PSGD proximal problem has the form: Here w t is the current iterate, z t is the index of the sample chosen and η t is the learning rate. For convex Ω, problem (3) can be solved in two steps: first solving the unconstrained problem and then using Euclidean projection onto Ω. Setting Ω = R d removes the need for projection and we recover SGD. When clear from the context we will use SGD to refer to both projected and un-projected variants. To solve the unconstrained problem one only needs to set the gradient to zero to recover the familiar closed form SGD update.

The BORAT Algorithm
In this section we detail the BORAT algorithm. We start by introducing BORAT's proximal problem that is exactly solved at each iteration. We explain its advantages and disadvantages in relation to the SGD proximal problem (3). Each BORAT proximal problem is best solved in the dual. Hence we next introduce the dual problem, which permits a far more efficient solution due to its low dimensionality. We then consider a special case of BORAT with minimal bundle size which we call Adaptive Learning-rates for Interpolation with Gradients (ALI-G). ALI-G permits a closed form solution and results in an automatically scaled gradient descent step. Specifically, it recovers a stochastic variant of the Polyak step size (Polyak, 1969), which offers competitive results in practice. This special case is used extensively in our experiments and in our analysis to establish the convergence rate of BORAT. Lastly, we detail our novel direct method for efficiently solving the general dual problem. This algorithm exploits the small size of the bundle to compute the exact optimum by solving a finite number of linear systems, removing the need for an inner iterative optimisation algorithm.

Primal Problem
For tackling problems of type (P) we identify two deficiencies in the proximal view of SGD (3). First, the approximation of the loss permits negative values even though the loss for (P) is defined to be non-negative. Second, the accuracy of a linear model quickly deteriorates for functions with high curvature away from the site of the approximation. Due to this crude model, the selection of η t for all t is critical for achieving good performance with SGD. BORAT aims to address these deficiencies by using a model composed of a pointwise maximum over N linear approximations to better model the loss. One of the linear approximations is chosen to be a constant function equal to the loss lower bound, that is, 0. By including this linear approximation, we address the first deficiency. The second deficiency is addressed by making use of the remaining N − 1 linear approximations. These extra approximations allow us to model variation over the parameter space and positive curvature of the loss. A model of this form in combination with a proximal term is commonly known as a bundle of size N . The main disadvantage of bundles is that they require multiple gradient evaluations to be performed and then held in memory. Hence we in this work only consider N ≤ 5, except where mentioned otherwise. Each linear approximation of the loss is constructed at a pointŵ n t using a different loss function ℓ z n t . The subscript t indicates the iteration number, and n indexes over the N linear approximations. With this notation we first introduce a bundle of size 1 as: If we setŵ 1 t to w t we recover the SGD proximal problem. Thus SGD effectively uses a bundle of size N = 1. We next introduce an expanded expression for a bundle of size N , before showing how to convert this into the compact form of max n∈[N ] a n⊤ (w − w t ) + b n . Within a bundle each linear approximation is formed around a different pointŵ n t . Hence in order to write each linear approximation in the aforementioned compact form we split each linear term in two. The first piece is a constant term, that does not depend on w, and is a multiple of the distance between the current iterate and the centre of each approximation w n t − w t . The second term depends on the distance w − w t , for all linear approximations. This gives the following expanded form for a bundle of size N as: For BORAT we always setŵ 1 t = w t . We additionally use the last linear approximation to enforce the lower bound on the loss. This is done by setting ∇ℓ zt ( We give details on how we selectŵ n t for n ∈ {2, ..., N − 1} in Section 3.5. Thus each bundle is composed of N − 1 linear approximations of the function, and the lower bound on the loss. These linear approximations of the loss need to be stored in memory during each step. Hence, in order to fit on a single GPU we only consider small bundle sizes in this work (N ≤ 5). For clarity we depict a 1D example for a bundle with N = 3 in Figure  2. Unlike SGD, the BORAT proximal problem (6) is not smooth and hence cannot be solved by simply setting the derivatives to zero. Instead we choose to solve each proximal problem in the dual. t and ℓ z 2 t of the loss function f are shown in red and blue (solid lines). The bundle is formed of the max of three linear approximations (dashed lines) and a proximal term. Two of these linear approximations are formed using the loss functions ℓ z 1 t and ℓ z 2 t , and the last enforces the known lower bound on the loss. Here the BORAT approximation gives a more accurate model than approximation used by ALI-G, which would only include the linearization atŵ 1 t and the lower bound B. In this simple example this improved accuracy allows for a larger step to be taken towards the minimum.

Dual Problem
The dual of (6) is a constrained concave quadratic maximisation over N dual variables α 1 , ..., α N , and can be concisely written as follows (see supplementary material for derivation): Here A t is a N × d matrix whose n th row is ∇ℓ zt (ŵ n t ). We define b t = [b 1 t , ..., b N t ] ⊤ , α = [α 1 , α 2 , . . . , α N ] ⊤ and ∆ N is a probability simplex over the N variables. The dual problem (7) has a number of features that make it more appealing for optimisation than the primal (6). First, the primal problem is defined over the parameters space w ∈ R d , where p is in the order thousands if not millions for modern DNNs. In contrast, the dual variables are of dimension N, where is N typically a small number > 10 due to the memory requirements. Second, the dual is smooth and hence allows for faster convergence with standard optimisation techniques. Furthermore, as will be seen shortly, we use the fact that the dual feasible region is a tractable probability simplex to design a customised algorithm for its solution. We detail this algorithm in this Section 3.6. Once we have found the dual solution α t , we recover the following update from the KKT conditions: The form of the update (8) deserves some attention. Since each row of A t is either the gradient of the loss ℓ z n t or a zero vector, and α t belongs to the probability simplex, the update step moves in the direction of a non-negative linear combination of negative gradients −∇ℓ z n t (ŵ n t ). Due to the definitions of ∇ℓ zt (ŵ N t ) = [0, ..., 0] ⊤ and b N t = 0, any weight given to α N reduces the magnitude of the resulting step. This has the effect that as the loss value gets close to zero BORAT automatically decreases the size of the step taken.

ALI-G (N=2)
We now consider a special case of BORAT with N = 2. Here the bundle is the point wise maximum over the linear approximation of the loss around the current point and the global lower bound B, which we assume is 0. This special case is worthy of extra attention for the following four reasons. First, this special case only requires one gradient evaluation per step and has a similar time complexity to SGD. Second, it admits a closed form solution. Third, we use this special case extensively in our analysis of the convergence rate of BORAT. Fourth, given the definitionsŵ 1 t = w t andŵ N t = 0, if we set N = 2 we recover an algorithm that simply scales the SGD learning rate. Specifically, it automatically scales down a maximal learning rate η by a factor α 1 ∈ [0, 1] to an appropriate value close to optimality. This is clear from the simplified version of Equation (8), which has the following form: Hence we will call this special case Adaptive Learning-rates for Interpolation with Gradients (ALI-G). For ALI-G the primal problem (6) simplifies to the following: Likewise, the dual problem (7) can be reduced to the following: The ALI-G dual is a maximisation over a constrained concave function in one dimension. Hence we can obtain the optimal point by projecting the unconstrained solution on to the feasible region. This results in the following closed form solution: This value of α 1 is then used in (9). The ALI-G update can be viewed as a stochastic analog of the Polyak step size (Polyak, 1969) with the addition of a maximal value η. Recall that, from the interpolation assumption, we have f ⋆ = 0. The ALI-G update is computationally cheap with the evaluation of ∇l zt (w t ) 2 being the only extra computation required over SGD. Hence when the interpolation assumption holds we suggest that ALI-G can be easily used in place of SGD.

The Advantage of Bundles with more than Two Pieces
While ALI-G has many favourable qualities, its local model of the loss is still crude. We next give two motivating examples to help demonstrate why using a more complex model of the loss often increases the robustness to η. Thus it may prove useful to use larger values of N in settings that are sensitive to step size η.
In the convex setting any function can be perfectly modelled by the point-wise maximum over an infinite number of linear approximations. While intractable, performing a minimisation over this model would recover the true optimum by definition. With this perfect model any value of η could be used. Setting η to a large enough value would recover the optimum in a single step. This example demonstrates that, at least asymptotically as the accuracy of the local model increases we can expect a reduced dependence on the correct scale of the step size.
w=− 3 5 w= 3 5 f : w → w 2 − |w| 3 Linearizations of f Figure 3: A simple example where the ALI-G step-size oscillates due to non-convexity. For this problem, ALI-G only converges ifs its maximal learning-rate η is less than 10. By contrast for the same example BORAT with N > 2 converges for all values of η. Additionally for η ≥ 10 it converges to the optimum in a single update.
Figure 3 provides a non-convex motivating example for use of larger values of N . Here we demonstrate a 1D symmetric function where ALI-G does not converge for large η. Instead it oscillates between the two values w = −3/5 and w = 3/5. However, if we were to use a bundle with a N ≥ 3 our model of the loss would include both the linear approximations at w = −3/5 and w = 3/5 simultaneously and hence when minimising over this model we converge to the optimum for any value of η. While this is a carefully constructed synthetic example it highlights why we would expect a more accurate model of the loss can help to reduce the dependence on the step size.

Selecting Additional Linear Approximations or the Bundle
When constructing bundles of size larger than two, we are faced with two design decisions regarding how to select additional linear approximations to add to the bundle. First, where in parameter space should we construct the additional linear approximationsŵ n t ? And second, should we use the same mini-batch of data when constructing the stochastic linear approximations, or should we sample a new batch to evaluate each linear approximation?
Selectingŵ n t . Ideally we would select the location of the linear approximationsŵ n t for n ∈ {2, ..., N −1} in order to maximise the progress made towards w ⋆ at each step. However, without the knowledge of w ⋆ a priori, due to the high dimensional and non-convex nature of problem (P) this is infeasible. Instead we make use of a a heuristic. Inspired by the work of previous bundle methods for convex problems (Smola et al., 2007, Asi andDuchi, 2019) we selectŵ n t using the following method: In other words we construct the bundle by recursively adding linear approximations centred at the current optimum. This method of selecting additional linear approximations is appealing as it requires no extra hyperparameters and helps refine the approximation in the region of parameter space that would be explored by the next update.
Re-sampling z for additional linear approximations. When constructing additional linear approximations we choose to re-sample z. Concretely, we use a new mini-batch of data to construct each stochastic linear approximation. While it is possible to construct all N − 1 non-zero linear approximations using the same batch of data we find this does not work well in practice. Indeed, such a method behaves similarly to taking multiple consecutive steps of SGD on the same mini-batch, which tends to produce poor optimisation performance.
Summary. We now summarise the bundle construction procedure for N > 2. We construct a bundle around w t by first using two linear approximations, one centred at w t and the second given by a known lower bound on the loss. We then sequentially add linear approximations until we have N . These extra linear approximations are constructed one at a time using new batches of data and centred around the point that is the current minimizer of the bundle. We note that each parameter update of BORAT uses N − 1 batches of data. Therefore BORAT updates the parameters N − 1 fewer times than SGD in an epoch (given the same batch-size). Once the bundle is fully constructed we update w t . At this stage we apply momentum (if enabled) and project back on the feasible set Ω. The construction of the bundle requires solving a minimization problem for each newly added piece when N ≥ 2. Therefore it is critical that such a problem gets solved very efficiently. We next detail how we do this by solving the corresponding dual problem for N ≥ 2.

Efficient Dual Algorithm to Compute N ≥ 2 Linear Pieces
In the general case (N > 2), the dual has more than 2 degrees of freedom and can not be written as a 1D minimisation. Thus the method derived for the case N = 2 is no longer applicable. This means it is not possible to obtain a simple closed form update. We must instead run an inner optimisation to solve (7) at each step. Many methods exist for maximising a concave quadratic objective over a simplex. Two algorithms particularly well suited to problems of the form (7) are (Frank and Wolfe, 1956) or Homotopy Methods (Bach et al., 2011). However, we propose a novel algorithm that exploits the fact that N is small to find the maximum directly. This method decomposes the problem of solving (7) into several sub-problems which provides two computational conveniences. First, the BORAT algorithm repeatedly searches for solutions to a bundle with only one newly added linear piece since the last search. As one might expect, this task shares a great number of sub-problems with the previous solution and allows for much of the computation to be reused. Second, our dual algorithm allows for a parallel implementation, which makes it very fast on the hardware commonly used for deep learning. To illustrate the efficiency of our dual method, the run time of the dual solution, that is finding α t once we have constructed (7), takes less than 5% of the time spent in the call of the optimiser. Note due to the large size of the the networks and the small size of N the majority of the call time is dominated computing A ⊤ t A t and A t α t . We now formally introduce our dual solution algorithm. Our method uses the observation that at the optimal solution in a simplex the partial derivatives will be equal for all nonzero dimensions. This observation can be formally stated as: In other words, the value of the partial derivative is shared among all coordinates of α * that are non-zero.
This proposition can be easily proved by contradiction. If the partial derivatives are not equal, then moving in the direction of the largest would result in an increase in function value. Likewise moving in the negative direction would produce a decrease. Hence the current point cannot be optimal. Please see Appendix B for a formal proof of Proposition 1.
In the following paragraphs we explain how this proposition can be used to break up the task of solving problem (7) into 2 N − 1 subproblems. Given a unique subset I of non-zero dimensions of α Each of these subproblems involves finding the unconstrained optimum and checking if this point lies within the simplex. We now give an example of a single subproblem. To simplify notation, let Q ηA ⊤ t A t . We note that: If we knew that α * had exclusively non-zero coordinates, then by applying Proposition (1) to the dual objective D we can recover a solution α * and by solving the following linear system: The first N rows of the system would enforce that α * satisfies the condition given by Proposition 1, and the last row of the linear system would ensure that the coordinates of α * sum to one. In the general case, we do not know which coordinates of α * are non-zero. However, since typical problems are in low-dimension N , we can enumerate all possibilities of subsets of non-zero coordinates for α * . We detail this further. We consider a non-empty subset I ⊆ [N ], for which we define: We then solve the corresponding linear subsystem: This φ (I) ∈ R |I| can then be lifted to R N by setting zeros at coordinates non-contained in I. Formally, we define ψ (I) ∈ R N such that: Therefore, given I ⊆ [N ], we can generate a candidate solution ψ (I) for problem (7) by solving a linear system in dimension |I|. In the following proposition, we establish that doing so for all possibilities of I guarantees to find the correct solution: Proposition 2 (Problem Equivalence) We define the set of feasible solutions reached by the different candidates ψ (I) : Then we have that: argmax In other words, we can find the optimal solution of (7) by enumerating the members of Ψ and picking the one with highest objective value.
This proposition is trivially true because Ψ is simply the intersection between (i) the original feasible set ∆ N and (ii) the set of vectors that satisfy the necessary condition of Optimality given by Proposition 1. This insight results in Algorithm 1 which returns a optimal solution to the dual problem. We characterise this claim in the following proposition.

Proposition 3 (Sets of Solutions)
Algorithm 1 returns a solution α * that satisfies α * ∈ argmax α∈∆ N D(α). This is true even when a the dual does not have a unique solution. Proof given in Appendix C.
Procedurally Algorithm 1 starts by computing Q = −ηA ⊤ t A t and b t to form a "master" system Qx = b. We consider each of the 2 N − 1 subsystems of Qx = b t defined by an element of the set set I (lines 1-3), where I represents the set of none-zero dimensions of each subsystem. For each subsystem we get a independent subproblem. To ensure the solution to each subproblem will satisfy N n=1 x n = 1 and all partial derivatives have equal value an extra row-column is introduced to each system (line 2). We then compute the point which satisfies the optimality conditions detailed in Proposition 1 for each subsystem, by solving for x. We then check if each of these points is feasible, that is, belongs to ∆ N by examining signs x n ≥ 0, ∀n ∈ {1, ..., n} (line 4). Note, we have by construction N n=1 x n = 1. Finally, we select α * as the feasible point with maximum dual value (lines 7-8). The optimal α * is then used to define the weight update (8). For example, if α * = [1, 0, ..., 0] ⊤ an SGD step w t+1 = w t − η∇ℓ zt (w t ) will be taken. The BORAT algorithm can be viewed as automatically picking the best step out of a maximum of (2 N − 1) possible options, where each option has a closed form solution. Although the computational complexity of this method is exponential in N , this algorithm is still very efficient in practice for two reasons. First, we only consider small N . Second, as sub problems can be solved independently, it permits a parallel implementation. With a fully parallel implementation the time complexity of this algorithm reduces to O N 3 . Empirically with such an implementation, for N ≤ 10 we observe approximately a linear relationship between N and time taken per training epoch. See Appendix F for a comparison of training epoch time between BORAT and SGD.

Computation considerations
While the method of adding linear approximations detailed in (13) requires running Algorithm 1 at each inner loop iteration, when adding an additional element α n if we keep track of the best dual value we need only compute the 2 N −1 new subproblems that include nonzero α n . Thus we only need to run Algorithm 1 once for each w t update, that is once per N − 1 batches.

Summary of Algorithm
The full BORAT method is outlined in Algorithm 2. The bundle is constructed in lines 4-6. The update is obtained in lines 7-9 using Algorithm 1. Finally, the updated parameters are projected to the feasible region Ω in line 9.
In the majority of our experiments, we accelerate BORAT with Nesterov momentum. We use Nesterov momentum as we find this helps produce strong generalise performance. The update step at line 9 of Algorithm 2 is then replaced by (i) a velocity update v t = µv t−1 − ηA t α and (ii) a parameter update w t+1 = Π Ω (w t + µv t ).

Justification and Analysis
The interpolation setting gives by definition, f ⋆ = 0. However, more subtly, it also allows the updates to rely on the stochastic estimate ℓ zt (w t ) rather than the exact but expensive f (w t ). Intuitively, this is possible because in the interpolation setting, we know the global minimum is achieved for each loss function ℓ zt (w t ) simultaneously. The following results formalise the convergence guarantee of BORAT in the stochastic setting. Note here we prove these result for BORAT with the minor modification, that is, all linear approximations are formed using the same mini-batch of data, ℓ z n t = ℓ zt for all n ∈ {2, ..., N − 1}. First, we consider the standard convex setting, where we additionally assume the interpolation assumption is satisfied and that each ℓ z is Lipschitz continuous. Next we consider an important class of nonconvex functions used for analysis in previous works related to interpolation (Vaswani et al., 2019a).

Theorem 1 (Convex and Lipschitz )
Let Ω be a convex set. We assume that for every z ∈ Z, ℓ z is convex and C-Lipschitz. Let w ⋆ be a solution of (P), and assume that we have perfect interpolation: ∀z ∈ Z, ℓ z (w ⋆ ) = 0. Then BORAT for N ≥ 2 applied to f satisfies: Hence BORAT recovers the same asymptotic rate as SGD without the need to reduce the learning rate η. In the Appendix D we show that for convex and β-smooth and the α-strongly convex settings BORAT recovers rates of O(1/T ) and O(exp(αT )) respectively. Note the earlier version of this work provides convergence results for ALI-G without perfect interpo-lation. We follow earlier work (Vaswani et al., 2019a) and provide a convergence rate for BORAT applied to non-convex functions that satisfy the Restricted Secant Inequality (RSI). A function is said to satisfy the RSI condition with constant µ over a the set Ω if the following holds: Theorem 2 (RSI ) We consider problems of type (P). We assume ℓ z satisfies the RSI with constant µ, smoothness constant β and perfect interpolation e.g. l z (w * ) = 0, ∀z ∈ Z. Then if set η ≤η = min{ 1 4β , 1 4µ , µ 2β 2 } then in the worst case we have: In this setting BORAT recovers the same asymptotic rate as SGD.

Related work
Bundle Methods. Bundle methods have been primarily proposed for the optimisation of non-smooth convex functions (Lemaréchal et al., 1995, Smola et al., 2007, Auslender, 2009). However, these works do not treat the stochastic case, as they consider small problems where the full gradient can be cheaply evaluated. To our knowledge our work is the first to propose a bundle method for the optimisation of stochastic non-convex problems.
Interpolation in Deep Learning. The interpolation property of DNNs was utilised by early efforts to analyse the convergence speed of SGD. These works demonstrate that SGD achieves the convergence rates of full-batch gradient descent in the interpolation setting (Ma et al., 2018a, Vaswani et al., 2019b, Zhou et al., 2019. Such works are complementary to ours in the sense that they provide a convergence analysis of an existing algorithm for deep learning. In a different line of work, Liu et al. (2019a) propose to exploit interpolation to prove convergence of a new acceleration method for deep learning. However, their experiments suggest that the method still requires the use of a hand-designed learning-rate schedule.
Adaptive Gradient Methods. Similarly to BORAT, most adaptive gradient methods also rely on tuning a single hyperparameter, thereby providing a more pragmatic alternative to SGD that needs a specification of the full learning-rate schedule. While the most popular ones are Adagrad (Duchi et al., 2011), RMSPROP (Tieleman and Hinton, 2012), Adam (Kingma and Welling, 2014) and AMSGrad (Reddi et al., 2018), there have been many other variants (Zeiler, 2012, Orabona and Pál, 2015, Défossez and Bach, 2017, Levy, 2017, Mukkamala and Hein, 2017, Zheng and Kwok, 2017, Bernstein et al., 2018, Chen and Gu, 2018, Shazeer and Stern, 2018, Zaheer et al., 2018, Chen et al., 2019, Loshchilov and Hutter, 2019, Luo et al., 2019. However, as pointed out in (Wilson et al., 2017), adaptive gradient methods tend to give poor generalization in supervised learning. In our experiments, the results provided by BORAT are significantly better than those obtained by the most popular adaptive gradient methods. Recently, Liu et al. (2019b) have proposed to "rectify" Adam with a learning-rate warm up, which partly bridges the gap in generalization performance between Adam and SGD. However, their method still requires a learning-rate schedule, and thus remains difficult to tune on new tasks.
Adaptive Learning-Rate Algorithms. Vaswani et al. (2019b) show that one can use line search in a stochastic setting for interpolating models while guaranteeing convergence. This work is complementary to ours, as it provides convergence results with weaker assumptions on the loss function, but is less practically useful as it requires up to four hyperparameters, instead of one for BORAT. Less closely related methods, included secondorder ones, adaptively compute the learning-rate without using the minimum (Schaul et al., 2013, Martens and Grosse, 2015, Tan et al., 2016, Zhang et al., 2017, Baydin et al., 2018, Wu et al., 2018, Li and Orabona, 2019, Henriques et al., 2019, but do not demonstrate competitive generalization performance against SGD with a well-tuned hand-designed schedule. L 4 Algorithm. The L 4 algorithm (Rolinek and Martius, 2018) also uses a modified version of the Polyak step-size. However, the L 4 algorithm computes an online estimate of f ⋆ rather than relying on a fixed value. This requires three hyperparameters, which are in practice sensitive to noise and crucial for empirical convergence of the method. In addition, L 4 does not come with convergence guarantees. In contrast, by utilizing the interpolation property and a single learning rate, our method is able to (i) provide reliable and accurate minimization with only a single hyperparameter, and (ii) offer guarantees of convergence in the stochastic convex setting.
Frank-Wolfe Methods. The proximal interpretation in Equation (10) allows us to draw additional parallels between ALI-G and existing methods. In particular, the formula of the learning-rate α 1 may remind the reader of the Frank-Wolfe algorithm (Frank and Wolfe, 1956) in some of its variants (Locatello et al., 2017), or other dual methods Jaggi, 2013, Shalev-Shwartz and. This is because such methods solve in closed form the dual of problem (10), and problems in the form of (10) naturally appear in dual coordinate ascent methods (Shalev-Shwartz and Zhang, 2016). When no regularization is used, ALI-G and Deep Frank-Wolfe (DFW) (Berrada et al., 2019a) are procedurally identical algorithms. This is because in such a setting, one iteration of DFW also amounts to solving (10) in closed-form -more generally, DFW is designed to train deep neural networks by solving proximal linear support vector machine problems approximately. However, we point out the two fundamental advantages of BORAT over DFW: (i) BO-RAT can handle arbitrary (lower-bounded) loss functions, while DFW can only use convex piece-wise linear loss functions; and (ii) as seen previously, BORAT provides convergence guarantees in the convex setting.
SGD with Polyak's Learning-Rate. Oberman and Prazeres (2019) extend the Polyak step-size to rely on a stochastic estimation of the gradient ∇ℓ zt (w t ) only, instead of the expensive deterministic gradient ∇f (w t ). However, they still require evaluating f (w t ), the objective function over the entire training data set, in order to compute its learning-rate, which makes the method impractical. In addition, since they do not do exploit the interpolation setting nor the fact that regularization can be expressed as a constraint, they also require the knowledge of the optimal objective function value f ⋆ . We also refer the interested reader to the recent analysis of Loizou et al. (2020), which provides a set of improved theoretical results.
aProx Algorithm. Most similar to this work (Asi and Duchi, 2019) have recently introduced the aProx algorithm, a family of proximal stochastic optimisation algorithms for convex problems. The aProx "truncated model" and the aProx "relatively accurate models" share aspects to ALI-G and BORAT respectively. However, there are four clear advantages of this work over (Asi and Duchi, 2019) in the interpolation setting, in particular for training neural networks. First, our work is the first to empirically demonstrate the applicability and usefulness of the algorithm on varied modern non-convex deep learning tasks -most of our experiments use several orders of magnitude more data and model parameters than the small-scale convex problems of (Asi and Duchi, 2019). Second, our analysis and insights allow us to make more aggressive choices of learning rate than (Asi and Duchi, 2019). Indeed, (Asi and Duchi, 2019) assume that the maximal learning-rate is exponentially decaying, even in the interpolating convex setting. In contrast, by avoiding the need for an exponential decay, the learning-rate of BORAT requires only one hyperparameters instead of two for aProx. Third, our analysis proves fast convergence in function space rather than iterate space. Fourth, unlike BORAT Asi and Duchi (2019) do not include the global lower bound in their aProx "relatively accurate models" and does not provide details on how to solve each proximal problem efficiently.

Experiments
We split the experimental results into two sections. The first section demonstrates the strong generalisation performance of ALI-G and BORAT on a wide range of tasks. Here we compare ALI-G and BORAT to other single hyperparameter optimisation algorithms. The second section shows for tasks where SGD and ALI-G are sensitive to the learning rate, using larger N increases both the robustness to the learning rate and task regularisation hyperparameter. For experiments we chosen to investigate BORAT with (N = 3) and (N = 5) and refer to the resulting algorithms as BORAT3 and BORAT5 respectively. For N ≥ 2 BORAT uses more than one batch of data for each update. In order to give a fair comparison we keep the number of passes through the data constant for all experiments. This has the effect that BORAT3 and BORAT5 respectively make a half and a quarter of the weight updates of SGD and ALIG.  (Paszke et al., 2017). We use our implementation of L 4 , which we unit-test against the official TensorFlow implementation (Abadi et al., 2015). we employ the official implementation of DFW ‡ and we re-use their code for the experiments on SNLI and CIFAR.

Effectiveness of ALI-G and BORAT
We empirically compare ALI-G and BORAT to the optimisation algorithms most commonly used in deep learning using standard loss functions. Where the hyperparameters of each algorithm are cross validated to select a best performing model. We consider a wide range of problems: training a wide residual network on SVHN, training a Bi-LSTM on the Stanford Natural Language Inference data set, training both wide residual networks and densely connected networks on the CIFAR data sets and lastly training a wide residual network on the Tiny Imagenet data set. For these problems, we demonstrate that ALI-G (N = 2) and BORAT for N ∈ {3, 5} obtains comparable performance to SGD with a hand-tuned learning rate schedule, and typically outperforms adaptive gradient methods. Finally, we empirically assess the performance of BORAT and its competitors in terms of training objective on CIFAR-100 and ImageNet, in order to demonstrate the scalability of BORAT to large-scale settings. Note that the tasks of training wide residual networks on SVHN and CIFAR-100 are part of the DeepOBS benchmark (Schneider et al., 2019), which aims at standardising baselines for deep learning optimisers. In particular, these tasks are among the most difficult ones of the benchmark because the SGD baseline benefits from a manual schedule for the learning rate where as ALI-G and BORAT uses a single fixed value. Despite this, our set of experiments demonstrate that ALI-G and BORAT obtains competitive performance in relation to SGD. In addition, our methods significantly outperforms adaptive gradient methods. All Experiments were performed on a single GPU (SVHN, SNLI, CIFAR) or on up to 4 GPUs (ImageNet). Method. For SGD, we use the manual schedule for the learning rate of Zagoruyko and Komodakis (2016). For L 4 Adam and L 4 Mom, we cross-validate the main learning-rate hyperparameter α to be in {0.0015, 0.015, 0.15} (0.15 is the value recommended by (Rolinek and Martius, 2018)). For other methods, the learning rate hyperparameter is tuned as a power of 10. The ℓ 2 regularization is cross-validated in {0.0001, 0.0005} for all methods but BORAT. For BORAT, the regularization is expressed as a constraint on the ℓ 2 -norm of the parameters, and its maximal value is set to 100. SGD, BORAT and BPGrad use a Nesterov momentum of 0.9. All methods use a dropout rate of 0.4 and a fixed budget of 160 epochs, following (Zagoruyko and Komodakis, 2016 does not converge on this task. However, note that L 4 Adam achieves accurate results. Even though SGD benefits from a hand-designed schedule, BORAT and other adaptive methods obtain comparable performance on this task.

Bi-LSTM on SNLI
Setting. We train a Bi-LSTM of 47M parameters on the Stanford Natural Language Inference (SNLI) data set (Bowman et al., 2015). The SNLI data set consists of 570k pairs of sentences, with each pair labeled as entailment, neutral or contradiction. This large scale data set is commonly used as a pre-training corpus for transfer learning to many other natural language tasks where labeled data is scarcer (Conneau et al., 2017) -much like ImageNet is used for pre-training in computer vision. We follow the protocol of (Berrada et al., 2019a) and we re-use their results for the baselines.
Method. For L 4 Adam and L 4 Mom, the main hyperparameter α is cross-validated in {0.015, 0.15} -compared to the recommended value of 0.15, this helped convergence and considerably improved performance. The SGD algorithm benefits from a hand-designed schedule, where the learning-rate is decreased by 5 when the validation accuracy does not improve. Other methods use adaptive learning-rates and do not require such a schedule. The value of the main hyperparameter η is cross-validated as a power of ten for the BORAT algorithm and for previously reported adaptive methods. Following the implementation by (Conneau et al., 2017), no ℓ 2 regularization is used. The algorithms are evaluated with the Cross-Entropy (CE) loss and the multi-class hinge loss (SVM), except for DFW which is designed for use with an SVM loss only. For all optimisation algorithms, the model is trained for 20 epochs, using a batch size of 64, following (Conneau et al., 2017).
Results.  In red, SGD benefits from a hand-designed schedule for its learning-rate. In black, adaptive methods have a single hyperparameter for their learning-rate. With an SVM loss, DFW and ALI-G are procedurally identical algorithms -but in contrast to DFW, ALI-G can also employ the CE loss. Methods in the format X * re-use results from (Berrada et al., 2019a). SGD † is the result from (Conneau et al., 2017).

Wide Residual Networks and Densely Connected Networks on CIFAR
Setting. We follow the methodology of our previous work (Berrada et al., 2019a). We test two architectures: a Wide Residual Network (WRN) 40-4 (Zagoruyko and Komodakis, 2016) and a bottleneck DenseNet (DN) 40-40 (Huang et al., 2017). We use 45k samples for training and 5k for validation. The images are centred and normalized per channel. We apply standard data augmentation with random horizontal flipping and random crops.
Method. We compare ALI-G and BORAT to other common single hyperparameter optimisers. Here we cross validate the hyperparameters in order to find a best performing network for each method. AMSGrad was selected in (Berrada et al., 2019a) because it was the best adaptive method on similar tasks, outperforming in particular Adam and Adagrad. In addition to these baselines, we also provide the performance of L 4 Adam, L 4 Mom, AdamW (Loshchilov and Hutter, 2019) and Yogi (Zaheer et al., 2018). We follow the cross validation scheme of (Berrada et al., 2019a) Table 3 details the results of the comparison of ALI-G and BORAT against other single hyperparameter optimisers for the CE loss only. In this setting, ALI-G and BORAT outperforms AMSGrad, AdamW, Yogi, L 4 Adam and L 4 Mom and constant step size SGD by a large margin. This is true for all values of N . The second best method is DFW which also has the restriction that it can only be used in conjunction with the hinge loss. ALI-G and BORAT produce the test accuracy achievable using SGD with the manual learning rate schedules from (Huang et al., 2017) and (Zagoruyko and Komodakis, 2016) for half of the model data set combinations considered here. For these tasks BORAT provides state of the art results without the need for the learning rate to be manually adapted through out training. For the remaining two combinations; training a DN on Cifar10 and training a WRN on Cifar100, BORAT lags in performance by approximately 0.2% and 2% respectively. With minor variation depending on which version of BORAT is used. Note SGD with a hand tuned learning rate schedule provides a reasonable upper limit of the generalisation performance achievable due to the amount of time that has been put into improving the schedule.  Table 3: Test accuracy of single hyperparameter optimisation methods. Each reported result is an average over three independent runs; the standard deviations and optimal hyperparameters are reported in Appendix F (the standard deviations are at most 0.4 for ALI-G and BORAT).

Wide Residual Networks on Tiny ImageNet
Setting. Tiny ImageNet contains 200 classes for training where each class has 500 images. The validation set contains 10,000 images. All images are RGB with 64x64 pixels. The test-set contains 10,000 images however ground truth labels are not freely available. We again use the Wide Residual Network (WRN) detailed in Section 6.1.1. The images are centred and normalized per channel. We apply standard data augmentation with random horizontal flipping and random crops.
Method. We investigate using SGD (with a constant step-size), ALI-G and BORAT to train a WRN on Tiny ImageNet. We use Adam (Kingma and Ba, 2015) as an indicator of what can be expected from a popular adaptive method applied to the same task. We make use of both the cross entropy (CE) and multi-class hinge (SVM) losses. The learning-rate hyperparameter η is cross validated in powers of 10, a batch-size of 128 and a training time of 200 epochs was used for all experiments. The ℓ 2 regularisation is cross validated in powers of 10 for ADAM. For constant step-size SGD, ALI-G and BORAT we make use of the constrained base regularisation detailed in Section 2, cross validating r ∈ {50, 100, 150, 200, 250}. Additionally all methods excluding Adam use a Nesterov momentum of 0.9.
Results. The best results for SGD, ALI-G, BORAT and Adam are shown in table 4. For both losses Adam performs worst, only achieving 2.4% validation accuracy when optimising the SVM loss. The best results are achieved by BORAT with N = 5 and N = 3 for the CE and SVM losses respectively. These results suggest the added gain of BORAT is more pronounced for larger more challenging data sets.  Table 4: Validation accuracy of single hyperparameter optimisation methods on the Tiny ImageNet data set for cross entropy and hinge loss.

Comparing Training Performance on CIFAR-100
In this section, we empirically assess the performance of ALI-G and its competitors in terms of training objective on CIFAR-100. In order to have comparable objective functions, the ℓ 2 regularization is deactivated. The learning-rate is selected as a power of ten for best final objective value, and the batch-size is set to its default value. For clarity, we only display the performance of SGD, Adam, Adagrad and BORAT (DFW does not support the cross-entropy loss). The L 4 methods diverge in this setting. Here SGD uses a constant learning-rate to demonstrate the need for adaptivity. Therefore all methods use one hyperparameter for their learning-rate. All methods use a fixed budget of 200 epochs for WRN-CIFAR-100 and 300 epochs for DN-CIFAR-100. As can be seen, ALI-G and BORAT provides better training performance than the baseline algorithms on both tasks. Here BORAT3 and BORAT5 are slightly slower than ALI-G due to the fewer number of parameter updates.

Training at Large Scale
We demonstrate the scalability of BORAT by training a ResNet18 (He et al., 2016) on the ImageNet data set. In order to satisfy the interpolation assumption, we employ a loss function tailored for top-5 classification (Lapin et al., 2016), and we do not use data augmentation. Our focus here is on the training objective and accuracy. ALI-G and BORAT uses the following training setup: a batch-size of 1024 split over 4 GPUs, a ℓ 2 maximal norm of 400 for w, a maximal learning-rate of 10 and no momentum. As can be seen in figure 5, ALI-G and BORAT reaches 99% top-5 accuracy in 12 epochs, and accurately minimises the objective function to 2 · 10 −4 within 90 epochs.

Robustness of BORAT
In this section we show that using additional linear approximations increases the robustness of BORAT to its learning rate η and the problem regularisation amount r. BORAT produces high accuracy models for a wider range of values than SGD and ALI-G. Hence on problems where SGD and ALI-G are sensitive to their learning rate BORAT with N ≥ 2 produces the best performance. These results demonstrate the advantage of using a larger bundle to model the loss in each proximal update. In order to illustrate this increased robustness we assess the stability of constant step size SGD, ALI-G, BORAT3 and BORAT5 to their hyperparameters. This is done by completing a grid search over η and r for a number of tasks while holding the batch size and epoch budget constant. This allows us to assess the range of values where these algorithms produce high accuracy models. We perform this grid search for six tasks split over the CIFAR100 and Tiny ImageNet data sets. We choose these data sets as they are challenging yet a model capable of interpolation can still fit on a single GPU. We compare against SGD with a constant learning rate for two reasons. First, SGD effectively uses a bundle of size 1, composed of only the linearization around the current iterate. Second, constant step size SGD has a single learning rate hyperparameter with the same scale as BORAT and also permits the constraint based regularisation described in Section 2.

Wide Residual Networks on CIFAR100 and Tiny ImageNet
Setting. For a description of the CIFAR100 and Tiny Imagenet data sets please refer to sections 6.1.3 and 6.1.4 respectively.
Method. We run four separate experiments on the CIFAR100 data set. The first two of these experiments examine the robustness to hyperparameters of SGD, ALI-G and BORAT when in combination with the cross entropy (CE) and multi-class hinge (SVM) losses. For the second two experiments on CIFAR100 we investigate the same algorithms in the presents of label noise. We limit our self to the CE loss for these two experiments. Label noise is applied by switching the label of the images in the training set with probability p to a random class label in the same super class. We repeat the experiments without label noise on the more challenging Tiny ImageNet data set, resulting in six different settings. For all algorithms we train a wide residual network (WRN) detailed in Section 6.1.3 and make use of a Nesterov momentum of 0.9 as we found its use produced superior generalisation performance. Step Size (η) Step Size (η) Colour represents test performance, where darker colours correspond to higher values. For both losses, increasing N allows for higher learning rates to be used while still producing convergent behaviour. Additionally for the CE loss and η ∈ {1.0, 10.0} larger N allows for a smaller r or greater levels of regularisation to be used. Consequently increasing N improves the over all robustness to hyperparameters, while not sacrificing generalisation performance.
Results. Figure 6 details the robustness of SGD, ALI-G and BORAT for the CE and SVM losses without label noise. In these two settings ALI-G exhibits similar robustness to SGD, however it outperforms SGD in terms of the performance of the best model trained. On the CE loss SGD and ALI-G are reasonably robust, providing good results for the majority of hyperparameter combinations with η ≤ 1. For the SVM loss all algorithms are sensitive to their learning rate η. SGD and ALI-G only produce an increase in accuracy for small values of η, hence no model achieves high accuracy in the 200 epoch budget. For both losses, increasing N improves the robustness and produces convergent behaviour for larger values of η and smaller values of r. This is particularly pronounced for the SVM loss. Here permitting larger η allows for a high accuracy network to be trained within the 200 epoch budget.
For both losses and η = 0.01 BORAT3 and BORAT5 slightly under perform compared to ALI-G and SGD. This is simply a consequence of these methods making significantly fewer updates. In spite of this the resultant effect of increasing the bundle size is a larger range of hyperparameters that produce good generalisation performance. The results from the label noise experiments are shown in Figure 7. For p = 0.1 the results closely mirror those where no label noise is used, shown in Figure 6, however here all models achieve roughly 0 − 5% worse test performance. When the noise is increased to p = 0.5 the test error increases drastically by 0 − 25%. In both cases p = 0.1 and p = 0.5, increasing N permits to obtain good results for a larger range of values of η and r. Interestingly, using N = 3, 5 results in the best performance when p = 0.5. This is somewhat expected since using a larger bundle size means more samples are used to calculate each parameter update and hence helps reduce the effect of the label noise. Step Size (η) Step Size (η) : Test accuracy of SGD, ALI-G, BORAT3 and BORAT5 robustness to hyperparameters when trained on CIFAR100 with noisy labels. Here we give results with two different levels of noise p = 0.1 (upper row) and p = 0.5 (lower row). For readability purposes, the colour encodes test accuracy, where darker colours correspond to higher values. Increasing N allows for higher learning rates and greater levels of regularisation to be used. Additionally, when the level of noise is high (p = 0.5, lower row), BORAT3 and BORAT5 significantly outperforms SGD (N = 1) and ALI-G (N = 2). Figure 8 details the robustness of SGD, ALI-G and BORAT for the Tiny Imagenet experiments. These results show BORAT offers improved robustness of on multiple data sets as we recover similar performance to CIFAR100. For the CE loss increasing N allows for slightly higher learning rates and greater levels of regularisation to be used for the CE loss. For the SVM loss BORAT produces models with reasonable accuracy with η ∈ {0.1, 1.0} opposed to SGD and ALI-G that produce poor results for all learning rates excluding η = 1.0. Consequently, when training with the SVM loss for 200 epochs, BORAT produces a drastically better models than SGD and ALI-G. This happens despite BORAT making significantly fewer updates of w t in the 200 epochs. Step Size (η) Step Size (η) Across all six experiments BORAT consistently produces high accuracy models for a larger range of hyperparameter combinations than SGD or ALI-G. This is particularly pronounced for the SVM loss experiments where the decrease sensitivity of BORAT makes all the difference between finding hyperparameters that result in a high accuracy model and not. Consequently, in comparison to its competitors, BORAT is systematically more robust to the choice of learning-rate and regularization hyper-parameter, and also offers better generalization more often than not. Therefore we particularly recommend using BORAT for tasks where hyper-parameter tuning is a highly time consuming task.

Discussion
In this work have introduced BORAT a bundle method designed for the optimisation of DNNs capable of interpolation. We detailed a special case of BORAT with a minimal bundle size that we name ALI-G. ALI-G produces a algorithm that automatically decays the steps size as the loss of the iterate approaches its minimal value. By only needing to find a good maximal learning rate that is automatically decayed we remove the labour and compute intensive task of finding a good learning rate schedule. Using standard publicly available data sets, we have shown that a ALI-G is highly effective in a number of settings.
When tuning a single hyperparameter for the learning-rate, ALI-G achieves state of the art performance while being simple to implement and having minor additional computational cost over SGD. For problems where ALI-G is sensitive to its learning rate its robustness and performance can be enhanced by using BORAT which makes use of a greater bundle size. BORAT produces high accuracy models for a larger range of hyperparameters and is particularly effective when using the multi-class hinge loss or in the presence of label noise. However, BORAT also achieves similar generalisation performance to ALI-G in all settings considered. Our results suggest that BORAT is presently the most robust single hyperparameter optimisation method for DNNs. It may be possible to modify BORAT so that it can make a parameter update after each gradient evaluation by cleverly selecting how linear approximations are added and removed from the bundle. If done correctly this may lead to a further increase in the speed of convergence. However we leave this to future work. When keeping the total number of gradient evaluations constant, a downside of using a larger bundle size for BORAT is that the number of parameter updates decreases. Despite this, in our experiments, BORAT is able to obtain good results within the same budget of passes through the data. Notably, this also means that BORAT has a time per epoch comparable to that of adaptive gradient methods. Another potential criticism of BORAT is that its memory footprint grows linearly with the bundle size. However in practice, our results show that using a bundle size of three, which corresponds to the same memory cost as ADAM, is often sufficient to obtain improved robustness. Finally, the applicability of BORAT can be limited by the required assumption of interpolation. While we argue that this interpolation property is satisfied in many interesting use cases, it may also not hold true for one or more of the following confounding factors: (i) limited size of the neural network, such as those used in embedded devices; (ii) large size of the training data set, which is becoming increasingly common in machine learning; and (iii) complexity of the loss function, such as in adversarial training. Furthermore, the concept of interpolation itself is ill-defined for unsupervised tasks. Thus, another interesting direction of future work would be to generalise BORAT to non-interpolating settings. This could be achieved by adapting an estimate of the minimum value of the objective function online whilst retaining the desirable quality of minimal hyperparameter tuning.

Appendix A. Dual Derivation
Lemma 1 The dual of the primal problem (6), can be written as follows: Where A t is defined as a d × N matrix, whose n th row is ∇ℓ z n t (ŵ n t ), b n t = [b 1 t , ..., b N t ] ⊤ and α = [α 1 , α 2 , . . . , α N ] ⊤ .
Proof : We start from the primal problem: The inner problem is now smooth inw and ζ. We write the KKT conditions: We define ∆ N {α ∈ R N : N n=1 α n = 1, α n ≥ 0, n = 1, ..., N } as probability simplex over the elements of α. Thus when we plug in the KKT conditions we obtain: Using the definitions of A t , b n t and α, we can recover the following compact form for the dual problem: Appendix B. Proof of Proposition 1 Proposition 1 Let F : R N → R be a concave function. Let us define α * = argmax α∈∆ F (α). Then there exists c ∈ R such that: ∀n ∈ [N ] such that α n * > 0, we have: In other words, the value of the partial derivative is shared among all coordinates of α * that are non-zero.
Proof : We start form the optimally condition for constrained problems: We consider the pointα. Without loss of generality we assumeα is equal to α * for all but two dimensions i and j. We letα i = 0 andα j = α * i + α * j . Hence we knowα ∈ ∆ as we have iα i = i α * i = 1 andα i ≥ 0, ∀i. Letting α =α in (27) give us: Rearranging we have: Hence for any α * i = 0, we have: Noticing, this results holds for any i and j gives us the following and proves the result.
Appendix C. Proof of Proposition 3 Proposition 3 Algorithm 1 returns a solution α * that satisfies α * ∈ argmax α∈∆ N D(α). This is true even when a the dual does not have a unique solution.
Proof : let x * ∈ argmax α∈∆N D(α) such that x * has a maximal number of zero coordinates.
x * exists as we know the solution set is empty. Let I be the set of non-zero coordinates of x * . We denote by S (I) the set of soultions to the linear system associated with I: S (I) is a polytope as the intersection intersection between the probability simplex and a linear sub-space. There for it admits and vertex representation: such that Conv(·) denotes the convex hull operation, and V (I) is a finite set. Since S (I) contains x * [I], S (I) is non-empty and neither is V (I) . let v be and element of V (I) . Note, that if v had one or more zero coordinates, it would be a solution after lifting to R N it would also be a optimal solution to the dual as good as x * while having more zero coordinates than x*, this is impossible by the definition of x * . Thus we can conclude v has exclusively non-zero coordinates.
Since v is a external point of S (I) , see section 2.1 of Bertsekas (2009) for a definition. It follows from proposition 2.1.4b of (Bertsekas, 2009) the columns ofQ are independent. Therefore the linear system admin a unique solution x * , which is found by Algorithm 1.

Lemma 2
Adding additional linear approximations to the bundle of BORAT can never result in a lower maximal dual value. Formally: Proof : Any vector of ∆ l can be lifted to ∆ m by appending m − l zeros to it, which does not impact the value of the objective function. The lifted set ∆ l is then a subset of ∆ m , hence the result (maximisation performed performed over a larger space).

Theorem 3 (Convex )
We assume that Ω is a convex set, and that for every z ∈ Z, ℓ z is convex and. Let w ⋆ be a solution of (P), and assume that we have perfect interpolation: ∀z ∈ Z, ℓ z (w ⋆ ) = 0. Then BORAT for N ≥ 2 applied to f satisfies: where D is defined in (7).
Proof : A consequence of Lemma 2 is that the convergence rate given by Theorem 5 improves as N increases.

D.1 Convex and C-Lipshits
Theorem 1 (Convex and Lipschitz ) We assume that Ω is a convex set, and that for every z ∈ Z, ℓ z is convex and C-Lipschitz. Let w ⋆ be a solution of (P), and assume that we have perfect interpolation: ∀z ∈ Z, ℓ z (w ⋆ ) = 0. Then BORAT for N ≥ 2 applied to f satisfies: Proof : We start from (50), hence we have: From Lemma 2 we additionally have that D 2 (α * ) ≤ D N >2 (α * ) hence consider the N = 2 as this provides the worse rate.
For N = 2 we have exactly two sub problems, and hence can write the dual in following compact form: We now introduce I T and J T as follows: Defining 1 to be the indicator function we can write: From our definition of I T for all t ∈ I T we have η g zt 2 ≥ ℓ zt (w t ) hence we can write: Summing over T : Using w T +1 − w ⋆ 2 ≥ 0, we obtain: Likewise, using the observation that ℓ z ≥ 0, we can write: We minimize the right hand-side overŵ, which gives: Thus by choosingŵ = w and re-ordering, we obtain the following result (a.k.a. the Polyak-Lojasiewicz inequality): However we have ℓ z (w) = 0 which concludes the proof.

Lemma 5
For any a, b ∈ R d , we have that: Proof : This is a simple application of the parallelogram law, but we give the proof here for completeness.

Lemma 6
Let z ∈ Z. Assume that ℓ z is α-strongly convex and achieves its (possibly constrained) minimum at w ⋆ ∈ Ω. Then we have: Proof : By definition of strong-convexity Bubeck (2015), we have: In addition, since w ⋆ minimizes ℓ z , then necessarily: Combining the two equations gives the desired result.
Taking expectations: Hence we recover the final rate: Theorem 6 (Strongly Convex -Small Eta) We assume that Ω is a convex set, and that for every z ∈ Z, ℓ z is α-strongly convex and β-smooth. Let w ⋆ be a solution of (P), and assume that we have prefect interpolation: ∀z ∈ Z, ℓ z (w ⋆ ) = 0. Then BORAT for N ≥ 2 and applied to f with a η ≤ 1 2β satisfies: Proof : We start from (56): We next use Lemma 3 write: Now we use our assumption on η to give: taking expectations: Using Lemma 6 We use a trivial induction over t and write: Given an arbitrary w ∈ R d , we now wish to relate the distance w − w ⋆ 2 to the function values f (w) − f (w ⋆ ). From smoothness, we have: w 2

E.0.2 subproblems
Algorithm 1 solves (2 N − 1) sub problem. Each of these linear systems corresponds to a loci of the simplex, defining the feasible set. We refer to each of these (2 N − 1) options as subproblems. However, each subproblem can also be interpreted as a sub-system of Qα = b, (see line 2 of Algorithm 1 for definitions). For (N = 3) the form of Qα = b is detailed in (121).
where c is defined in Theorem 1. Each subproblem is defined by setting a unique subset of the dual variables α n to zero, before solving for the remaining variables. For (N = 3) we have exactly seven subproblems, which we each give a unique label, see table 5. For clarity, we detail a few specific subproblems. The SGD subproblem corresponds to setting α 2 , α 3 = 0, and recovers the SGD update. Likewise, the ESGD step corresponds to setting α 1 , α 3 = 0 and recovers an update similar to the extra gradient method. Finally, by setting α 2 = 0 before solving the system we recover the a ALI-G like step, hence we label this subproblem as the ALI-G step. If a subproblem results in a α ∈ ∆ 3 we refer to that subproblem as feasible, conversely, if it results in a α / ∈ ∆ 3 we refer to that subproblem as infeasible. Algorithm 2 computes the dual value for all feasible subproblems and selects the largest. This subproblem's α is then used in (120) to update the parameters. The closed form solutions for α for each of the 7 subproblems are listed in table 5. We use a subscript to show that α belongs to a certain subproblem. For example α SGD = [1, 0, 0]. .
Due to spatial constraints we state these results without proof. The dual value for each subproblem is simply derived by inserting the closed form expression for α for each subproblem detailed in table 5 into (7). These dual values observe a tree like hierarchy, This hierarchy is a consequence of the subproblems maximising the dual over progressively larger spaces, ∆ n−1 ⊂ ∆ n .

E.0.4 Feasible Subproblems
To give a worst case rate on the convergence of BORAT with N = 3, we first prove for smooth functions and small η, only certain subproblems will be feasible. We start with a useful lemma, before proving this result.

Lemma 7
Let z ∈ Z. Assume that ℓ z is β-smooth. If we define w ′ = w − η∇ℓ z (w) and η ≤ 1 β then we have: Note that we do not assume that ℓ z is convex.

Lemma 8 (Feasible Subproblems)
If ℓ z n t has smoothness constant β and we set η ≤ 1 2β for the BORAT algorithm with (N = 3) detailed in Section E.0.1, the ALIG, EALIG and MAX3 subproblems will always be infeasible. Proof : We start by showing the ALIG step is infeasible. From Lemma 3 we have: .
Hence for any 1 2β ≥ η, η > ℓ z n t (wt) ||gz t || 2 cannot hold. We now use a similar argument to show that the EALIG subproblem is infeasible. For EALIG to be feasible we require α 3t EALIG > 0,

Lemma 10
Let z ∈ Z. Assume that ℓ z is β-smooth and non-negative on R d . Then we have: Note that we do not assume that ℓ z is convex.
Proof : Let w ∈ R d . By Lemma 3.4 of Bubeck (2015), we have: Therefore we can write: And since ∀ u, l z (u) ≥ 0, we have: We now choose u = − 1 β ∇l z (w), which yeilds: which gives the desired result.
In this section we derive the rate for each of the remaining feasible steps, that is, SGD, ESGD and M AX2.
Rearranging gives: Plugging in to the the definition of γ t , gives the desired result.
Taking expectations with respect to z t : Noting that w k does not depend on z t , and neither does w * due to the interpolation property: For the this step to be convergent we need the following condition to hold 2η 3 β 2 µ + η 2 β 2 − ηµ ≤ 0. However, for η ≤η = min{ 1 4β , 1 4µ , µ 2β 2 } we have: Hence, we recover the rate:

E.4 Worst Case Rate
It is clear by inspection that the worst case rate derived corresponds to the MAX2 subproblem. Hence in the worst case this step is taken for all t, and thus a trivial induction gives the result of Theorem 2.

F.1 Cifar Hyperparameters and Variance
Here we detail the hyperparameters and variance for the ALI-G and BORAT results reported in table 6. For other optimisation methods please refer to Appendix E of Berrada et al. (2019b).

F.2 Empirical Run Time
Here we detail the effect of increasing N on the run time of BORAT. Due to each update requiring N − 1 gradient evaluations BORAT with N ≥ 2 take significantly longer between updates than other methods. However, BORAT achieves good empirical convergence rates taking N − 1 fewer parameter updates than other methods, as shown in the results section.  reduces to O N 3 as discussed in Section 3.6. Additionally we need only run Algorithm 1 once every N − 1 batches so the per epoch time complexity is O N 2 . Of course in practice Algorithm 1 is only responsible for a fraction of the run time, where its contribution is determined by the relative size of the model and N . Table 7 show with a parallel implementation the effect of the this extra computation on the training epoch time isn't significant and the time complexity scales approximately linearly with N . Moreover, Table 7 shows for large scale learning problems, such as ImageNet the extra run time when increasing N is negligible.