HIGH-DIMENSIONAL WEIGHT-UPDATE HYPERPARAMETERS BY IMPLICIT DIFFERENTIATION
Year: 2022

Abstract:
Machine learning training methods depend plentifully and intricately on hyperparameters, motivating automated strategies for their optimisation. Many existing algorithms restart training for each new hyperparameter choice, at considerable computational cost. Some hypergradient-based one-pass methods exist, but these either cannot be applied to arbitrary optimiser hyperparameters (such as learning rates and momenta) or take several times longer to train than their base models. We extend these existing methods to develop an approximate hypergradient-based hyperparameter optimiser which is applicable to any continuous hyperparameter appearing in a differentiable model weight update, yet requires only one training episode, with no restarts. We also provide a motivating argument for convergence to the true hypergradient, and perform tractable gradient-based optimisation of independent learning rates for each model parameter. Our method performs competitively from varied random hyperparameter initialisations on several UCI datasets and Fashion-MNIST (using a one-layer MLP), Penn Treebank (using an LSTM) and CIFAR-10 (using a ResNet-18), in time only 2–3x greater than vanilla training.

heading: 1 INTRODUCTION
text: Many machine learning methods are governed by hyperparameters: quantities other than model parameters or weights which nonetheless influence training (e.g. optimiser settings, dropout probabilities and dataset configurations). As suitable hyperparameter selection is crucial to system performance (e.g. Kohavi & John (1995)), it is a pillar of efforts to automate machine learning (Hutter et al., 2018, Chapter 1), spawning several hyperparameter optimisation (HPO) algorithms (e.g. Bergstra & Bengio (2012); Snoek et al. (2012; 2015); Falkner et al. (2018)). However, HPO is computationally intensive and random search is an unexpectedly strong (but beatable; Turner et al. (2021)) baseline; beyond random or grid searches, HPO is relatively underused in research (Bouthillier & Varoquaux, 2020).
Recently, Lorraine et al. (2020) used gradient-based updates to adjust hyperparameters during training, displaying impressive optimisation performance and scalability to high-dimensional hyperparameters. Despite their computational efficiency (since updates occur before final training performance is known), Lorraine et al.’s algorithm only applies to hyperparameters on which the loss function depends explicitly (such as `2 regularisation), notably excluding optimiser hyperparameters.
Our work extends Lorraine et al.’s algorithm to support arbitrary continuous inputs to a differentiable weight update formula, including learning rates and momentum factors. We demonstrate our algorithm handles a range of hyperparameter initialisations and datasets, improving test loss after a single training episode (‘one pass’). Relaxing differentiation-through-optimisation (Domke, 2012) and hypergradient descent’s (Baydin et al., 2018) exactness allows us to improve computational and memory efficiency. Our scalable one-pass method improves performance from arbitrary hyperparameter initialisations, and could be augmented with a further search over those initialisations if desired.

heading: 2 WEIGHT-UPDATE HYPERPARAMETER TUNING
text: In this section, we develop our method. Expanded derivations and a summary of differences from Lorraine et al. (2020) are given in Appendix C.
exact hypergradients from the implicit function theorem (IFT) ( ), with our method’s approximate hypergradients ( ) superimposed. We target optimal validation loss ( ), adjusting weights w based on the training loss. Classical weight updates (for fixed hyperparameters λ) converge ( , ) to the best-response line w∗(λ) ( ); the IFT gives hyperparameter updates ( ) leading to a minimum of validation loss along w∗(λ). Our approximate hyperparameter updates ( ) differ in magnitude from these exact updates, but still give useful guidance.

heading: 2.1 IMPLICIT FUNCTION THEOREM IN BILEVEL OPTIMISATION
text: Consider some model with learnable parameters w, training loss LT , optimisation hyperparameters λ and validation loss LV . We use w∗,λ∗ to represent the optimal values of these quantities, found by solving the following bilevel optimisation problem:
(a) λ∗ = arg min λ LV (λ,w∗(λ)) , such that (b) w∗(λ) = arg min w LT (λ,w) . (1)
The optimal model parameters w∗ may vary with λ, making LV an implicit function of λ alone. We approach the outer optimisation (1a) similarly to Lorraine et al. (2020) and Majumder et al. (2019), using the hypergradient: the total derivative of LV with respect to the hyperparameters λ. One strategy for solving (1) is therefore to alternate between updating w for several steps using ∂LT∂w and updating λ using the hypergradient, as shown in Figure 1 (by and ).
Carefully distinguishing the total differential dλ and the partial differential ∂λ, we have dLV dλ = ∂LV ∂λ + ∂LV ∂w∗ ∂w∗ ∂λ . (2)
While derivatives of LV are easily computed for typical loss functions, the final derivative of the optimal model parameters (∂w ∗
∂λ ) presents some difficulty. Letting square brackets indicate the evaluation of the interior at the subscripted values, we may rewrite ∂w ∗
∂λ as follows:
Theorem 1 (Cauchy’s Implicit Function Theorem (IFT)) Suppose for some λ′ and w′ that[ ∂LT ∂w ] λ′,w′
= 0. If ∂LT∂w is a continuously differentiable function with invertible Jacobian, then there exists a function w∗(λ) over an open subset of hyperparameter space, such that λ′ lies in the open subset, [ ∂LT ∂w ] λ,w∗(λ) = 0 and
∂w∗
∂λ = − ( ∂2LT ∂w∂wT )−1 ∂2LT ∂w∂λT . (3)
w∗(λ) is called the best response of w to λ (Figure 1). While (3) suggests a route to computing ∂w ∗
∂λ , inverting a potentially high-dimensional Hessian in w is not computationally tractable.

heading: 2.2 APPROXIMATE BEST-RESPONSE DERIVATIVE
text: To develop and justify a computationally tractable approximation to (3), we mirror the strategy of Lorraine et al. (2020). Consider the broad class of weight optimisers with updates of the form
wi(λ) = wi−1(λ)− u(λ,wi−1(λ)) (4)
for some arbitrary differentiable function u, with i indexing each update iteration. We deviate here from the approach of Lorraine et al. (2020) by admitting general functions u(λ,w), rather than assuming the particular choice uSGD = η ∂LT∂w (see Appendix A.4 for more details). In particular, this allows λ to include optimiser hyperparameters. Differentiating (4) and unrolling the recursion gives[
∂wi ∂λ ] λ′ = − ∑
0≤j<i ( ∏ 0≤k<j [ I− ∂u ∂w ] λ′,wi−1−k )[ ∂u ∂λ ] λ′,wi−1−j , (5)
where j indexes our steps back through time from wi−1 to w0, and all w depend on the current hyperparameters λ′ — see Appendix C.2 for the full derivation. Now, we follow Lorraine et al. and assume w0, . . . ,wi−1 to be equal to wi. With this, we simplify the product in (5) to a jth power by evaluating all derivatives at wi. This result is then used to approximate ∂w ∗
∂λ by further assuming that wi ≈ w∗. These two approximations lead to[
∂wi ∂λ ] λ′ ≈ −  ∑ 0≤j<i ( I− ∂u ∂w )j ∂u ∂λ  λ′,wi(λ′) ≈ [ ∂w∗ ∂λ ] λ′ . (6)
We reinterpret i as a predefined look-back distance, trading off accuracy and computational efficiency.
The combination of these approximations implies wi = w∗ for all i, which is initially inaccurate, but which we would expect to become more correct as training proceeds. In mitigation, we perform several weight updates prior to each hyperparameter update. This means derivatives in earlier terms of the series of (6) (which are likely the largest, dominant terms) are evaluated at weights closer to w∗, therefore making the summation more accurate. In Section 4, we show that the approximations described here result in an algorithm that is both practical and effective.
Our approximate result (6) combines the general weight update of Majumder et al. (2019) with the overall approach and constant-weight assumption of Lorraine et al. (2020). The latter empirically show that an approximation similar to (6) leads to a directionally-accurate approximate hypergradient; we illustrate the approximate updates from our derivations in Figure 1 (by ).

heading: 2.3 CONVERGENCE TO BEST-RESPONSE DERIVATIVE
text: To justify the approximations in (6), note that the central part of that equation is a truncated Neumann series. Taking the limit i→∞ , when such a limit exists, results in the closed form[
∂w∗
∂λ ] λ′ ≈ − [( ∂u ∂w )−1 ∂u ∂λ ] λ′,w∗(λ′) . (7)
This is precisely the result of the IFT (Theorem 1) applied to u instead of ∂LT∂w ; that is, substituting the simple SGD update uSGD(λ,w) = η ∂LT∂w into (7) recovers (3) exactly. Thus, under certain conditions, our approximation (6) converges to the true best-response Jacobian in the limit of infinitely long look-back windows.

heading: 2.4 HYPERPARAMETER UPDATES
text: Substituting (6) into (2) yields a tractable approximation for the hypergradient dLVdλ , with which we can update hyperparameters by gradient descent. Our implementation in Algorithm 1 (Figure 4) closely parallels Lorraine et al.’s algorithm, invoking Jacobian-vector products (Pearlmutter, 1994) during gradient computation for memory efficiency via the grad_outputs argument, which also provides the repeated multiplication for the jth power in (6). Thus, we retain the O(|w| + |λ|) time and memory cost of Lorraine et al. (2020), where | · | denotes cardinality. The core loop to compute the summation in (6) comes from an algorithm of Liao et al. (2018). Note that Algorithm 1 approximates dLVdλ retrospectively by considering only the last weight update rather than any future weight updates.
Unlike differentiation-through-optimisation (Domke, 2012), Algorithm 1 crucially estimates hypergradients without reference to old model parameters, thanks to the approximate-hypergradient construction of Lorraine et al. (2020) and (6). We thus do not store network weights at multiple time
steps, so gradient-based HPO becomes possible on previously-intractable large-scale problems. In essence, we develop an approximation to online hypergradient descent (Baydin et al., 2018).
Optimiser hyperparameters generally do not affect the optimal weights, suggesting their hypergradients should be zero. However, in practice, w∗ is better reinterpreted as the approximately optimal weights obtained after a finite training episode. These certainly depend on the optimiser hyperparameters, which govern the convergence of w, thus justifying our use of the bilevel framework.
We emphasise training is not reset after each hyperparameter update — we simply continue training from where we left off, using the new hyperparameters. Consequently, Algorithm 1 avoids the time cost of multiple training restarts. While our locally greedy hyperparameter updates threaten a short-horizon bias (Wu et al., 2018), we still realise practical improvements in our experiments.

heading: 2.5 REINTERPRETATION OF ITERATIVE OPTIMISATION
text: Originally, we stated the IFT (3) in terms of minima of LT (zeros of ∂LT∂w ), and substituting uSGD = η ∂LT∂w into (7) recovers this form of (3). However, in general, (7) recovers the Theorem for zeros of u, which are not necessarily minima of the training loss. Despite this, our development can be compared to Lorraine et al. (2020) by expressing u as the derivative of an augmented ‘pseudo-loss’ function. Consider again the simple SGD update uSGD, which provides the weight update rule wi = wi−1 − η ∂LT∂w . By trivially defining a pseudo-loss L = ηLT , we may absorb η into a loss-like derivative, yielding wi = wi−1 − ∂L∂w . More generally, we may write L = ∫ u(λ,w) dw.
Expressing the update in this form suggests a reinterpretation of the role of optimiser hyperparameters. Conventionally, our visualisation of gradient descent has the learning rate control the size of steps over some undulating landscape. Instead, we propose fixing a unit step size, with the ‘learning rate’ scaling the landscape underneath. Similarly, we suppose a ‘momentum’ could, at every point, locally squash the loss surface in the negative-gradient direction and stretch it in the positive-gradient direction. In aggregate, these transformations straighten out optimisation trajectories and bring local optima closer to the current point. While more complex hyperparameters lack a clear visualisation in this framework, it nevertheless allows a broader class of hyperparameters to ‘directly alter the loss function’ instead of remaining completely independent, circumventing the problem with optimiser hyperparameters noted by Lorraine et al. (2020). Figure 2 and Appendix A.3 support this argument.

heading: 3 RELATED WORK
text: Kohavi & John (1995) first noted different problems respond optimally to different hyperparameters; Hutter et al. (2018, Chapter 1) summarise the resulting hyperparameter optimisation (HPO) field.
Black-box HPO treats the training process as atomic, selecting trial configurations by grid search (coarse or intractable at scale), random search (Bergstra & Bengio (2012); often more efficient) or population-based searches (mutating promising trials). Pure Bayesian Optimisation (Močkus et al., 1978; Snoek et al., 2012) guides the search with a predictive model; many works seek to exploit its sample efficiency (e.g. Swersky et al. (2014a); Snoek et al. (2015); Wang et al. (2016); Kandasamy et al. (2017); Lévesque et al. (2017); Perrone et al. (2019)). However, these methods require each proposed configuration to be fully trained, incurring considerable computational expense. Other techniques infer information during a training run — from learning curves (Provost et al., 1999; Swersky et al., 2014b; Domhan et al., 2015; Chandrashekaran & Lane, 2017; Klein et al., 2017), smaller surrogate problems (Petrak, 2000; van den Bosch, A. et al., 2004; Krueger et al., 2015; Sparks et al., 2015; Thornton et al., 2013; Sabharwal et al., 2016) or intelligent resource allocation (Jamieson & Talwalkar, 2016; Li et al., 2017; Falkner et al., 2018; Bertrand et al., 2017; Wang et al., 2018). Such techniques could be applied on top of our algorithm to improve performance.
In HPO with nested bilevel optimisation, hyperparameters are optimised conditioned on optimal weights, enabling updates during training, with a separate validation set mitigating overfitting risk; this relates to meta-learning (Franceschi et al., 2018). Innovations include differentiable unrolled stochastic gradient descent (SGD) updates (Domke, 2012; Maclaurin et al., 2015; Baydin et al., 2018; Shaban et al., 2019; Majumder et al., 2019), conjugate gradients or hypernetworks (Lorraine & Duvenaud, 2018; Lorraine et al., 2020; MacKay et al., 2019; Fu et al., 2017), solving one level while penalising suboptimality of the other (Mehra & Hamm, 2020), and deploying Cauchy’s implicit
function theorem (Larsen et al., 1996; Bengio, 2000; Luketina et al., 2016; Lorraine et al., 2020). Donini et al. (2020) extrapolate training to optimise arbitrary learning rate schedules, extending earlier work on gradient computation (Franceschi et al., 2017), but do not immediately accommodate other hyperparameters. While non-smooth procedures exist (Lopez-Ramos & Beferull-Lozano, 2021), most methods focus on smooth hyperparameters which augment the loss function (e.g. weight regularisation) and work with the augmented loss directly. Consequently, they cannot handle optimiser hyperparameters not represented in the loss function (e.g. learning rates and momentum factors; Lorraine et al. (2020)). Many of these methods compute hyperparameter updates locally (not over the entire training process), which may induce short-horizon bias (Wu et al., 2018), causing myopic convergence to local optima.
Modern developments include theoretical reformulations of bilevel optimisation to improve performance (Liu et al., 2020; Li et al., 2020a), optimising distinct hyperparameters for each model parameter (Lorraine et al., 2020; Jie et al., 2020), and computing forward-mode hypergradient averages using more exact techniques than we do (Micaelli & Storkey, 2020). Although these approaches increase computational efficiency and the range of tunable parameters, achieving both benefits at once remains challenging. Our algorithm accommodates a diverse range of differentiable hyperparameters, but retains the efficiency of existing approaches (specifically Lorraine et al. (2020)). Inevitably, we also inherit gradient-based methods’ difficulties with robustness and discrete hyperparameters.

heading: 4 EXPERIMENTS
text: Our empirical evaluation uses the hardware and software detailed in Appendix A.1, with code available at https://github.com/rmclarke/OptimisingWeightUpdateHyperparameters.
Throughout, we train models using SGD with weight decay and momentum. We uniformly sample initial learning rates, weight decays and momenta, using logarithmic and sigmoidal transforms (see Appendix B.3), applying each initialisation in the following eight settings:
Algorithm 1 Scalable One-Pass Optimisation of High-Dimensional Weight-Update Hyperparameters by Implicit Differentiation
8317
Random No HPO; hyperparameters held constant at their initial values, as in random search. Random (× LR) Random with an extra hyperparameter increasing or decreasing the learning rate by multiplication after every T = 10 update steps (emulating a hyperparameter update). Random (3-batched) Random reprocessed to retain only the best result of every three runs. Allows Random to exceed our methods’ computational budgets; imitates population training. Lorraine Lorraine et al. (2020)’s method optimising weight decay; other hyperparameters constant. Baydin Baydin et al. (2018)’s method optimising learning rate only; other hyperparameters constant. OursWD+LR Algorithm 1 updating weight decay and learning rate; other hyperparameters constant. OursWD+LR+M Optimising all hyperparameters (adding momentum to OursWD+LR) OursWD+HDLR+M OursWD+LR+M with independent learning rates for each model parameter Diff-through-Opt Optimising all hyperparameters using exact hypergradients (Domke, 2012)
Any unoptimised hyperparameters are fixed at their random initial values. Ideally, we seek resilience to poor initialisations, which realistically arise in unguided hyperparameter selection. Hyperparameters are tuned on the validation set; this is combined with the training set for our Random settings, so each algorithm observes the same data. We use the UCI/Kin8nm dataset split sizes of Gal & Ghahramani (2016) and standard 60%/20%/20% splits for training/validation/test datasets elsewhere, updating hyperparameters every T = 10 batches with look-back distance i = 5 steps (except Baydin, which has no such meta-hyperparameters, so updates hyperparameters at every batch (Baydin et al., 2018)).
Numerical data is normalised to zero-mean, unit-variance. For efficiency, the computational graph is detached after each hyperparameter update, so we never differentiate through hyperparameter updates — essentially, back-propagation and momentum histories are truncated. Our approximate hypergradient is passed to Adam (Kingma & Ba, 2015) with meta-learning rate κ = 0.05 and default β1 = 0.9, β2 = 0.999. While these meta-hyperparameters are not tuned, previous work indicates performance is progressively less sensitive to higher-order hyperparameters (Franceschi et al., 2017; 2018; Majumder et al., 2019). As some settings habitually chose unstable high learning rates, we clip these to [10−10, 1] throughout.
UCI Energy: Proof of Concept First, we broadly illustrate Algorithm 1 on UCI Energy, using a one-layer multi-layer perceptron (MLP) with 50 hidden units and ReLU (Glorot et al., 2011) activation functions, trained under OursWD+LR for 4 000 full-batch epochs (see Appendix B.1 for details). Figure 3 shows the evolution of learning rate and weight decay from a variety of initialisations, overlaid on the performance obtained when all hyperparameters are kept fixed during training. Notice the trajectories are attracted towards the region of lowest test loss, indicating our algorithm is capable of useful learning rate and weight decay adjustment.

heading: 4.1 UCI AND KIN8NM DATASETS: ESTABLISHING ROBUSTNESS
text: Next, we consider the same 50-hidden-unit MLP applied to seven standard UCI datasets (Boston Housing, Concrete, Energy, Naval, Power, Wine, Yacht) and Kin8nm, in a fashion analogous to Gal & Ghahramani (2016) (we do not consider dropout or Bayesian formulations). We train 200 hyperparameter initialisations for 4 000 full-batch epochs. Table 1 shows results for three datasets, with complete loss evolution and distribution plots in Appendix B.4. Some extreme hyperparameters caused numerical instability and NaN final losses, which our averages ignore. NaN results are not problematic: they indicate extremely poor initialisations, which should be easier for the user to rectify than merely mediocre hyperparameters. Error bars are based on 1 000 sets of bootstrap samples.
Given the sparse distribution of strong hyperparameter initialitions, random sampling unsurprisingly achieves generally poor averages. Random (× LR) applies a learning rate multiplier, uniformly chosen in [0.95, 1.01]; these limits allow extreme initial learning rates to revert to more typical values after 400 multiplications. This setting’s poor performance shows naïve learning rate schedules cannot match our algorithms’ average improvements. A stronger baseline is Random (3-batched), which harshly simulates the greater computational cost of our methods by retaining only the best of every three Random trials (according to validation loss). This setting comes close to, but cannot robustly beat, our methods — a claim reinforced by final test loss distributions (Figure 9, Appendix B.4).
Lorraine et al. (2020)’s algorithm is surprisingly indistinguishable from Random in our trials, though it must account for three random hyperparameters while varying only one (the weight decay). We surmise that learning rates and momenta are more important hyperparameters to select, and poor choices cannot be overcome by intelligent use of weight decay. Baydin et al. (2018)’s algorithm, however, comes closer to the performance of Random (3-batched), indicating more successful intervention in its sole optimisable hyperparameter (learning rate). Variance is also much lower than the preceding algorithms, suggesting greater stability. That said, despite concentrating on a more important hyperparameter, Baydin still suffers from being unable to control every hyperparameter.
Our scalar algorithms (OursWD+LR and OursWD+LR+M) appear generally more robust to these initialisations, with average losses beating Random, Lorraine and Baydin. Figures 8 and 9 (Appendix B.4) clearly distinguish these algorithms from the preceding: we achieve performant results over a wider space of initial hyperparameters. Given it considers more hyperparameters, OursWD+LR+M predictably outperforms OursWD+LR, although the difference is less prominent in certain datasets. Unlike pure HPO, the non-Random algorithms in Table 1 vary hyperparameters during training, combining aspects of HPO and schedule learning. They are thus more flexible than conventional, static-hyperparameter techniques, even in high dimensions (Lorraine et al., 2020).
Diff-through-Opt exactly differentiates the current loss with respect to the hyperparameters, over the same i = 5 look-back window and T = 10 update interval. As an exact version of OursWD+LR+M (though subject to the same short-horizon bias), it unsurprisingly matches the other algorithms (Figures 8 and 9, Appendix B.4). However, our scalar methods’ proximity to this exact baseline is reassuring given our much-reduced memory requirements and generally comparable error bars. In these experiments, lengthening Diff-through-Opt’s look-back horizon to all 4 000 training steps, and repeating those steps for 30 hyperparameter updates, did not improve its performance (Appendix B.7).
OursWD+HDLR+M theoretically mitigates short-horizon bias (Wu et al., 2018) by adapting appropriately to high- and low-curvature directions. While this substantially improves performance on some datasets, it is outperformed by scalar methods on others (Figures 8 and 9, Appendix B.4). We also see slightly reduced stability, with more trials diverging to NaN final losses. Intuitively, the risk of overfitting to the validation set is expected to increase as we introduce more hyperparameters, which could explain this behaviour — committing to large learning rates in a few selected directions, as this method often does, may prove harmful if the loss surface dynamics change suddenly. However, we leave detailed investigation of the high-dimensional dynamics to future work.

heading: 4.2 LARGE-SCALE DATASETS: PRACTICAL SCALABILITY
text: Fashion-MNIST: HPO in Multi-Layer Perceptrons We train the same single 50-unit hidden layer MLP on 10 epochs of Fashion-MNIST (Xiao et al., 2017), using 50-sample batches. Table 2 and Figure 12a (Appendix B.6) show average test set cross-entropies over 100 initialisations. Clearlyoutlying final losses (above 103) are set to NaN to stop them dominating our error bars.
Echoing Section 4.1, for arbitrary hyperparameter initialisations, our methods generally converge more robustly to lower losses, even when (as in Figure 12a, Appendix B.6) NaN solutions are included in our statistics. Importantly, we see mini-batches provide sufficient gradient information for our HPO task. Diff-through-Opt’s failure to beat our methods is surprising; we suppose noisy approximate gradients may regularise our algorithms, preventing them from seeking the short-horizon optimum so directly, thus mitigating short-horizon bias (see our sensitivity study in Section 4.3). Diff-through-Opt with long look-back horizons does not improve performance for equal computation (Appendix B.7). Finally, median and best loss evolution plots are shown in Figures 5b and 5c, the latter including results for a Bayesian Optimisation baseline. For more details, see Appendices B.5 and B.6.
Penn Treebank: HPO in Recurrent Networks Now, we draw inspiration from Lorraine et al. (2020)’s large-scale trials: a 2-layer, 650-unit LSTM (Hochreiter & Schmidhuber, 1997) with learnable embedding, trained on the standard Penn Treebank-3-subset benchmark dataset (Marcus et al., 1999) for 72 epochs. To focus our study, we omit the dropout, activation regularisation and predefined learning rate schedules used by Lorraine et al., though we retain training gradient clipping to a Euclidean norm of 0.25. Training considers length-70 subsequences of 40 parallel sequences, using 50 random hyperparameter initialisations.
Table 2 and Figure 12b (Appendix B.6) show final test perplexities. They reflect our intuition that adjusting progressively more hyperparameters reduces average test losses, continuing the trend we have seen thus far. Highly bimodal final loss distributions for some algorithms cause wide bootstrap-sampled error bars. Learning rate adjustments show particular gains: Our algorithms and Diff-through-Opt perform particularly well in less-optimal configurations.
CIFAR-10: HPO in Convolutional Networks Finally, to demonstrate scalability, we train a ResNet-18 (He et al., 2016) on CIFAR-10 (Krizhevsky, 2009) for 72 epochs. We use the vanilla normalised dataset (since unbiased data augmentation over both training and validation datasets exhausts our GPU memory) and optimise hyperparameters as before, using 100-image batches.
Table 2 and Figure 12c (Appendix B.6) show our results. Our gains are now more marginal, with this setting apparently robust to its initialisation, though the general ranking of algorithms remains similar, and we retain useful improvements over our baselines. However, our best final accuracies fall short of state-of-the-art, suggesting a more intricate and clever setting may yield further performance gains.

heading: 4.3 MISCELLANEOUS STUDIES
text: Experiment Runtimes For our larger-scale experiments, we illustrate comparative runtimes in Figure 14 (Appendix B.6.1, where they are discussed in detail). Note from our Fashion-MNIST sample (Figure 5a) that all HPO algorithms have a lower computational cost than naïvely training two typical fixed hyperparameter initialisations, despite achieving substantial HPO effect. This compares extremely favourably to HPO methods relying on repeated retraining. Additional experiments in Appendix B.8.2 combine OursWD+LR+M with Asynchronous Hyperband (Li et al., 2020b) or Population-Based Training (Jaderberg et al., 2017), showing improved and unchanged performance/time trade-offs, respectively.
UCI Energy: Sensitivity Study Finally, we consider a range of update intervals T and look-back distances i on UCI Energy, performing 400 hyperparameter updates from 100 random initialisations on each, using OursWD+LR+M. We plot the median final test losses for each choice of T and i to the right of Figure 4, and also the median performance with no HPO (Random), which we outperform in every case. Performance generally improves with larger T — likely because the total number of weight updates increases with T , since the number of hyperparameter updates is fixed at 400. While increasing i gives more terms in our approximating series (6), these extra terms become more and more biased: they are evaluated at the current weight values instead of at the progressively more different past weight values. We theorise this trade-off explains the more complex relationship between i and final performance. More details are given in Appendix B.2.

heading: 5 CONCLUSION AND FUTURE WORK
text: We have presented an algorithm for optimising continuous hyperparameters, specifically those appearing in a differentiable weight update, including optimiser hyperparameters. Our method requires only a single training pass, has a motivating true-hypergradient convergence argument, and demonstrates practical benefits at a range of experimental scales without greatly sacrificing training time. We also tackle traditionally-intractable per-parameter learning rate optimisation on non-trivial datasets. However, this setting surprisingly underperformed its scalar counterpart; further work is necessary to understand this result.
As future work, our myopic approach could be extended to longer horizons by incorporating the principles of recent work by Micaelli & Storkey (2020), which presents a promising research direction. We also depend inconveniently on meta-hyperparameters, which are not substantially tuned. Ultimately, we desire systems which are completely independent of human configuration, thus motivating investigations into the removal of these settings.
SOCIETAL IMPACT / ETHICS STATEMENT
Our work fits into the broad subfield of automatic machine learning (AutoML), which aims to use automation to strip away the tedious work that is necessary to implement practical ML systems. Our method focuses on automating the oft-labelled ‘black art’ of optimising hyperparameters. This contributes towards the democratisation of ML techniques, which we hope will improve accessibility to non-experts.
However, our method also has some associated risks. For one, developers of unethical machine learning applications (e.g. for mass surveillance, identity theft or automated weaponry) may use our techniques to improve their systems’ performance. The heavily metric-dominated nature of our field creates additional concerns — for instance, end-users of our method may not appreciate that optimising naïvely for training and validation loss alone may result in dangerously poorly-trained or unethical models if the chosen metric, developer’s intentions and moral good do not align.
More broadly, users may rely excessively on HPO techniques to optimise their models’ performance, which can lead to poor results and inaccurate comparisons if the HPO strategy is imperfect. Further, reducing the need to consider hyperparmeter tuning abstracts away an important component of how machine learning methods work in practice. Knowledge of this component may then become less accessible, inhibiting understanding and future research insights from the wider community.
Our datasets are drawn from standard benchmarks, and thus should not introduce new societal risks. Similarly, our work aims to decrease the computational burden of HPO, which should mitigate the environmental impact of ML training — an ever more important goal in our increasingly environmentally-conscious society.

heading: REPRODUCIBILITY STATEMENT
text: All datasets we use are publicly available; for the Penn Treebank dataset, we provide a link in Table 4 (Appendix A.2). What little data processing we perform is fully explained in the corresponding subsection of Section 4.
Our mathematical arguments are presented fully and with disclosure of all assumptions in Section C.
Source code for all our experiments is provided to reviewers, and is made available on GitHub (https: //github.com/rmclarke/OptimisingWeightUpdateHyperparameters). This source code contains a complete description of our experimental environment, configuration files and instructions on the reproduction of our experiments.

heading: ACKNOWLEDGEMENTS
text: We acknowledge computation provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).
We thank Austin Tripp for his feedback on drafts of this paper.
Ross Clarke acknowledges funding from the Engineering and Physical Sciences Research Council (project reference 2107369, grant EP/S515334/1).

heading: A NOTES

heading: A.1 EXPERIMENTAL CONFIGURATION
text: Our experiments are variously performed on one of three sets of hardware, as detailed in Table 3. While system configurations varied between experiments, this does not impair comparability. All Penn Treebank and CIFAR-10 experiments were performed on the CSD3 Cluster; otherwise, all experiments comparing runtimes were performed on the Local Cluster. We make use of GPU acceleration throughout, with the PyTorch (Paszke et al., 2019) and Higher (Grefenstette et al., 2019) libraries.

heading: A.2 DATASET LICENCES
text: Our datasets are all standard in the ML literature. For completeness, we outline the licences under which they are used in Table 4.

heading: A.3 EXPANDED INTUITIVE REINTERPRETATION OF ITERATIVE OPTIMISATION
text: In this section, we take the opportunity to give a more verbose elaboration on the ideas presented in Section 2.5, which seeks to address methodological concerns raised by Lorraine et al. (2020) and reinterpret iterative optimisation in light of our derived algorithm.
Lorraine et al. (2020) base their derivations on the training loss LT , which does not depend on optimiser hyperparameters such as the learning rate or momentum coefficient. Differentiating LT with respect to these hyperparameters thus cannot provide any useful information for computing their gradient-based updates. Similarly, the Implicit Function Theorem cannot be applied to motivate their approximation in this case, since LT is not an implicit function of these hyperparameters. Consequently, Lorraine et al. declare their method is methodologically restricted to operate on hyperparameters on which the training loss depends directly.
In contrast, we work with an arbitrary update u, which does then depend on the optimiser hyperparameters, allowing us to take the corresponding derivatives and apply the Implicit Function Theorem. However, the optimisation objective remains the same, so we should rightly be suspicious that the methodological issue raised by Lorraine et al. may still apply. To address these remaining concerns,
we look to recast our development in the form of Lorraine et al. (2020). We do this by defining a ‘pseudo-loss’ L such that our weight updates in (4) are equivalent to
w← w − [ ∂L ∂w ] λ,w . (8)
We must acknowledge that our optimisation objective is now for the update u to be zero (equivalent to a stationary point of L), rather than the training loss LT being at a minimum. But with that caveat, we may apply the machinery of Lorraine et al. directly: we have constructed a ‘loss function’ which depends directly on all the hyperparameters, so may legitimately deploy their strategy. By comparison with (4), we have that, for any particular problem,
L = ∫ u(λ,w) dw. (9)
Often, in simple cases, the mathematical construction in (9) has a geometric interpretation. Consider the elementary case of vanilla SGD with learning rate η, which gives the update rule uSGD = η ∂LT∂w . The required pseudo-loss is then L = ηLT , which suggests two equivalent geometric interpretations of the optimisation problem we are solving.
In the first (conventional) interpretation, we have a fixed loss surface LT whose minimum we seek. To optimise η in this problem, we imagine standing on the loss surface at our current point, then consider the effect of taking different-sized steps across it. Clearly, this interpretation decouples η from the underlying function LT , so suffers from Lorraine et al.’s methodological issue. In the second (our novel) interpretation, we seek the minimum of the loss surface L, which we may traverse only by taking steps of precisely the magnitude of the local gradient. Now, when we optimise η, we are rescaling the surface beneath our update step, and consequently also rescaling the update step size. In this way, we have coupled η and LT into one objective L, which varies with η and so adheres to the methodological restrictions noted by Lorraine et al.. The key idea is thus that we have reinterpreted the problem to allow optimiser hyperparameters to dynamically transform the objective surface, which legitimises our derivation.
When we consider SGD with momentum, the situation becomes more complicated, because there is not an immediately obvoius closed form for L. Suppose we are descending a steep-sided valley of the original loss surface LT . The momentum acts to ‘pull’ the update step in the direction it has accumulated in its buffer (down the axis of the valley), so it acts to augment gradients in the valley’s downhill direction, and diminishes the effect of gradients in other directions.
If we are working with the pseudo-loss L, we need to achieve the same effect, but we do not have direct control over the update step size, which is determined by the local gradient of L. So if the update step is aligned with the accumulated momentum, we must ‘squash’ the loss surface about our current point in this direction, such that our fixed update step lands at a point further down the valley in the original, unscaled space of LT . Accordingly, if the update step is misaligned with the accumulated momentum, this squashing in the direction of the momentum buffer must be combined with a ‘stretching’ in the update step direction, so that we make less progress in that direction of the original space. In aggregate, this is a complex and intricate transformation of the space, such that we bring the ‘correct’ new point (as specified by the optimiser) to lie underneath the update step given by − ∂L∂w . To summarise, the key principle we seek to convey in Section 2.5 is that of creating a new pseudo-loss surface L, which is a transformed version of LT , such that the update steps taken by a hyperparameterfree optimiser over the pseudo-loss L are equivalent to those taken by the true optimiser over the original training loss LT .
A.4 EXPANDED DISCUSSION OF UPDATE FUNCTION u(λ,w)
Throughout our derivation, we refer to an arbitrary weight update function u(λ,w). While we have not expressed this in its most general form for clarity, we outline here that this formulation supports a wide array of non-trivial weight update optimisation algorithms beyond SGD.
u(λ,w) is precisely the quantity which is added to the network weights w during an update step. For instance, in the PyTorch framework (Paszke et al., 2019), the pre-implemented optimisers first
compute some update tensor ∆w, then perform the final operation w ← w −∆w to update the weights. Thus, we have u(λ,w) = ∆w, and automatic differentiation trivially gives the necessary derivatives to apply our algorithm.
Some optimisers have internal state, for instance the momentum buffer of SGD with momentum and the two rolling average buffers of Adam. Although this dependency is not notated in our presentation of u, it does not invalidate our development — the only extra care required is to detach these buffers from the computational graph after each hyperparameter update. For PyTorch’s implementation of SGD with learning rate η, momentum coefficient µ and weight decay coefficient ξ, we thus have
v← µv + ∂LT ∂w + ξw umomentum = ηv , (10)
and need simply take care to detach v from the computational graph after each update. For Adam, we may read the definition of u straight from its original paper (Kingma & Ba, 2015); in their notation:
uAdam = α · m̂t√ v̂t + , (11)
where t indexes time and m̂t, v̂t are internal states, and we simply detach those states from the computational graph periodically. Those states depend directly on the additional states mt, vt, and ultimately on the gradient∇θft(θt−1) (which is [ ∂LT ∂w ] wi−1 in our notation).

heading: B ADDITIONAL RESULTS

heading: B.1 UCI ENERGY: PROOF OF CONCEPT
text: Our proof of concept graphic (Figure 3) is constructed using the same experimental configuration as Section 4.1. Training 500 random initialisations of the model, weight decay and learning rate for 4 000 network weight update steps on Random (without any HPO) yields source data for the heatmap forming the background of Figure 3. The heatmap itself is the predictive mean of a Gaussian process (Williams & Rasmussen, 1996) using the sum of an RBF and white noise kernel. We fit the Gaussian process in log (input and output) space using the defaults of the Scikit-learn implementation (Pedregosa et al., 2011), which automatically select hyperparameters for the kernels. Any NaN data points are replaced with a large final test loss (150) for the purpose of computing this regression. We then train a small number of additional initialisations using OursWD+LR and plot these trajectories. As always, the training set provides network weight updates, the validation set provides hyperparameter updates and our results plot final performance on the test set. While we do not expect the optimal trajectory during training to coincide with gradients of our background heatmap, the latter provides useful context.
Generally, given our limited training budget, larger learning rates improve final performance up to a point, then cause instability if taken to extremes, and the movement of our trajectories reflects this. Many trajectories select large weight decay coefficients before doubling back on themselves to smaller weight decays — suggestive of our algorithm using the weight decay to avoid converging to specific minima. Trajectories generally converge to the valley of good test performance at a learning rate of around 10−1 and weight decay below 10−2, further demonstrating appropriate behaviour. In short, we are able to make sensible and interpretable updates to the hyperparameters of this problem, which lends support to the applicability of our method.

heading: B.2 UCI ENERGY: HYPERGRADIENT STUDIES
text: Our sensitivity studies are based on a similar experimental configuration to Section 4.1: for each choice of update interval T and look-back distance i, we generate 100 random initialisations of weight decay, learning rate and momentum. Each initialisation is trained for 400 hyperparameter steps, so increasing T permits more network weight updates; averaging over the 100 final test losses yields each result cell in Figure 6 and Figure 7. As in Section 4.1, we consider the OursWD+LR+M and Diff-through-Opt algorithms, optimising hyperparameters using Adam with hyper-learning rate κ = 0.05 and truncating gradient histories at the previous hyperparameter update.
B.2.1 SENSITIVITY OF UPDATE INTERVAL AND LOOK-BACK DISTANCE
8317
(a) Random (×1000)
Figure 6 shows our final performance results in this setting: our algorithm’s strongest configurations (Figure 6b) compare favourably with results obtained using the exact hypergradients of Diff-throughOpt (Figure 6c). Noting that the conditions necessary to guarantee our algorithm’s convergence to the best-response derivative (specifically that the current network weights are optimal for the current hyperparameters) are almost certainly not met in practice, it is unsurprising that Diff-through-Opt achieves greater performance. This difference is much smaller, however, than that between both HPO algorithms and the Random base case; given our substantial efficiency benefits over Diff-through-Opt, we make an acceptable compromise between accuracy and scalability. Further, these results justify our choice of T = 10 and i = 5 to match the setting of Lorraine et al. (2020). That said, recall these figures enforce consistency of number of hyperparameter updates, not of total computational expense; there remains a case for using large numbers of less precise hyperparameter updates (small T and i), which may promote faster convergence in practice.
Our results illustrate the trade-off in selecting a look-back distance i. In principle, increasing i includes more terms in our Neumann series approximation to the best-response derivative (6), which should result in closer convergence to the true derivative. But our derivation of this expression replaces the sequence w0,w1, · · · ,wi−1 from (5) with wi at all time steps. Because the jth summation step of (5) involves a product of terms in wi−1 back to wi−1−j , the approximation of each summation step in (6) becomes less and less accurate as j increases, since we expect earlier model weights to become more different from wi. A larger i causes us to continue the summation up to larger j, so introduces progressively less accurate terms into our approximations. Thus, increasing i is only of value while this inaccuracy is counteracted by the improved convergence of our summation, giving the behaviour we see here.

heading: B.2.2 ACCURACY OF HYPERGRADIENTS
text: We use this same experimental setting to evaluate the accuracy of our approximate hypergradients, but disable learning rate clipping so as to allow the natural behaviour of each algorithm to be displayed. . Our approximate and exact settings are only guaranteed to be trying to compute the same hypergradient at the very first hyperparameter update, so we consider this time step only. Using those first hyperparameter updates from 50 random initialisations of the hyperparameters and model, we compute the absolute error of the approximate hypergradient (OursWD+LR+M) as a percentage of the exact hypergradient (Diff-through-Opt), then plot the means of these errors in Figure 7.
Notice our learning rate hypergradients do not exhibit any particularly strong trends over the metahyperparameters we consider, but our approximate weight decay hypergradients seem to benefit from a shorter look-back distance. Some of the latter average proportional errors are large; a closer examination of the error distributions reveals these cells have outlying results at high errors. We hypothesise that the generally very small weight decay coefficients present more opportunity for numerical computational error than the generally larger learning rates, in which case hypergradients for different hyperparameters may have markedly different accuracies. Indeed, the apparently quantised errors in the learning rate hypergradients indicate precision truncation is influencing the results of this study. However, our previous coarser studies suggest that, despite being approximate, our method does not suffer catastrophically from inaccurate hypergradients.

heading: B.3 EXPERIMENT DETAILS
text: Since the hyperparameters associated with SGD only take meaningful values on subsets of the real line, we apply gradient descent to transformed versions of the hyperparameters. For non-negative hyperparameters, such as learning rate and weight decay, we work with the base-10 logarithm of the actual value; for momentum, which lies in the unit interval [0, 1], we work with the inverse-sigmoid of the actual value. These transformed search spaces mean a naïve implementation of gradient descent will not update hyperparameters to nonsensical values.
Our hyperparameter initialisations are drawn uniformly from a broad range, intended to be a superset of typical values to permit analysis of the realistic case where we have little intuition for the optimal choices. Thus, our hyperparameter ranges include typical values for a variety of applications, while giving each algorithm the opportunity to briefly stretch into more extreme values should these be useful on occasion. Where we have applied a logarithmic transformation to a hyperparameter (learning rate and weight decay), we draw samples uniformly in logarithmic space; where inversesigmoid transforms are used (momentum), we draw samples in natural space to ensure adequate exploration of values near 0 and 1. Table 5 summarises our initialisation and optimisation strategy.
Throughout, we seek robustness to a range of different initialisations, which is indicated by average final performance across a variety of initial hyperparameter choices. We include metrics of best performance to show each algorithm is capable of reaching the same vicinity as the best possible case, but this is not the focus of our work. Intuitively, a more naïve algorithm (e.g. Random) could stumble across the best possible final model weights, in which case no other algorithm would be able to beat it; this would not necessarily mean we would prefer the Random algorithm in a general setting.

heading: B.4 UCI AND KIN8NM DATASETS
text: Results from all UCI datasets we considered and Kin8nm are presented in Figures 8 and 9. These show many of the features already discussed in Section 4.1. We account for experimental runs ending in NaN losses by appropriately scaling the maximum value of the CDFs in Figure 9 to be less than 1.0. For each dataset, one random configuration is selected, and the resulting learning rate evolutions plotted in Figure 10. Note that UCI Naval generated a large number of NaN runs, which explains the prominent failure of certain algorithms to reach the top of our CDF plots.

heading: B.5 UCI, KIN8NM AND FASHION-MNIST: BAYESIAN OPTIMISATION BASELINES
text: To contextualise our methods against the best possible performance, we also perform 32 repetitions of 30-sample Bayesian Optimisation (BO) on selected UCI datasets, Kin8nm and Fashion-MNIST, using the same procedure as Random to train each sampled configuration (but now with an separate validation set used to compare configurations), and record these in Tables 6 and 7. We use the Python BO implementation of Nogueira (2014), retaining the default configuration that the first five of our 30 initialisations be drawn randomly from the space. All runtimes are measured based on 8-way parallel execution per GPU. Although BO involves a much greater computational burden than the other methods we consider, it should more robustly estimate the optimal hyperparameters.
Our average performance falls short of the BO benchmark, but we come substantially closer than the Random and Lorraine baselines, and similarly as close as Diff-through-Opt. On Fashion-MNIST especially, our methods are very competitive, achieving performance in the vicinity of BO with 30 times less computation. Considering our best runs, we seem to suffer only a minor performance penalty (compared to the best hyperparameters we found) by adopting Algorithm 1.
As an additional comparison, we show in Figure 11 the evolution of incumbent best test loss over wall time during our previous training of 200 (UCI and Kin8nm)/100 (Fashion-MNIST) initialisations. Clearly we would expect BO’s greater computation to result in smaller final losses, which we observe, but our methods are capable of reaching performant results much faster than BO. This underscores the principal benefit of our algorithm: with only a single training episode and a relatively small runtime penalty, we achieve a meaningful HPO effect, which reaches the vicinity of the ideal performance found by much slower and more computationally onerous methods.

heading: B.6 LARGE-SCALE DATASETS
text: For Fashion-MNIST, Penn Treebank and CIFAR-10, Figure 12 shows the empirical distributions of final test performance, contextualising the values shown in Table 2. In addition, Figure 13 and Table 8 show our results in terms of error percentages for Fashion-MNIST and CIFAR-10.
B.6.1 RUNTIMES
For our larger-scale experiments, we illustrate comparative runtimes in Figure 14. More complex model architectures, with more learnable parameters, suffer a greater computational cost from taking second derivatives with respect to each model weight, which appears more prominently in our Penn Treebank and CIFAR-10 experiments. However, our extension of Lorraine et al. (2020)’s method to additional hyperparameters doesn’t dramatically increase runtimes: second derivative computations themselves are responsible for most added complexity. In all three cases, our scalar algorithms have runtimes within that of naïvely training two fixed hyperparameter initialisations (such as those considered by Random), comparing extremely favourably to HPO methods requiring repeated retraining. Naturally, the massive increase in the number of hyperparameters optimised by OursWD+HDLR+M causes a prominent jump in runtimes, but even this algorithm remains within a practical factor of the Random baseline.

heading: B.6.2 BATCH NORMALISATION AND FASHION-MNIST
text: As batch normalisation rescales the parameter gradients for each training step, it is capable of further influencing how final performance is governed by learning rates and their variations. To study this effect, we repeat the Fashion-MNIST experiments described in Section 4.2 using batch normalisation in our simple MLP model. Our results are plotted alongside the non-batch normalised CDFs in Figure 15, with averages summarised in Table 9.
Subjectively, performance does not seem to dramatically change when batch normalisation is invoked, with CDFs having broadly the same shape. Considering means and medians as our figures of merit (Table 9) suggests the same general pattern: generally, our algorithms benefit slightly from the absence of batch normalisation, and other algorithms benefit slightly from its presence, but there is not a unanimous trend. We conclude that, in the context of this work, batch normalisation does not substantially improve or hinder performance.

heading: B.7 EFFECT OF SHORT-HORIZON BIAS
text: Attempts to quantify the effect of short-horizon bias in our experiments are frustrated by the need to perform computationally-intensive longer-horizon experiments, which has limited the scope of evaluations we can perform. However, we are able to study the potential improvements of longer horizons in three experimental settings.

heading: B.7.1 TRUNCATED PROBLEM, FULL HORIZON
text: Firstly, we consider the four study datasets UCI Energy/Power, Kin8nm and Fashion-MNIST, shortening training so it only lasts for 200 (UCI and Kin8nm) / 1000 (Fashion-MNIST) network weight update steps. This shorter training episode means we cannot compare to our results in Section 4, but
does permit full-horizon differentiation-through-optimisation. Long Diff-through-Opt does exactly this, performing 30 (UCI and Kin8nm) / 50 (Fashion-MNIST) hyperparameter updates using the maximum look-back horizon and restarting training with the new hyperparameters after each update. In this sense, we use 30 or 50 times as much computation as in our other algorithms (here: Random, OursWD+LR+M and Diff-through-Opt). The latter update their hyperparameters at the same intervals as before, giving 20 (UCI and Kin8nm) / 100 (Fashion-MNIST) updates overall. We take 60 random hyperparameter initialisations and run each of these four methods on them.
Our results are shown in Figure 16 and Tables 10 and 11. They show a stark contrast between datasets — on UCI and Kin8nm trials, Long Diff-through-Opt manages, as expected, to achieve substantially improved final performance over its short-horizon competition, but on Fashion-MNIST, performance is much more similar between HPO algorithms. Wu et al. (2018) intuit short-horizon bias by describing the simple setting of a quadratic with different curvature magnitudes in each dimension, and we expect the simpler UCI/Kin8nm datasets to match this clean depiction of the optimisation space more closely than Fashion-MNIST. Thus, we speculate the space of non-trivial datasets, such as Fashion-MNIST, is sufficiently complicated that merely eliminating short-horizon bias is insufficient to prevent the optimiser from getting stuck in similarly-performant local optima.

heading: B.7.2 FASHION-MNIST, MEDIUM HORIZON
text: Next, we note our Fashion-MNIST experiments of Section 4 run for 10 000 network weight steps, which can accommodate a longer look-back horizon for the purposes of our study. We introduce a new algorithm, Medium Diff-through-Opt, which is identical to our previous Fashion-MNIST Diff-through-Opt configuration, except the look-back horizon is now 200 steps, and the number of hyperparameter updates reduced to 50 to keep the same computational burden. This setting explores the trade-off between hypergradient accuracy and number of hypergradients computed. We apply this new configuration to the same 100 iterations as before, plotting combined results in Figure 17 and Table 12.
Again, we see that Fashion-MNIST seems resilient to the short-horizon bias, with this mediumhorizon configuration unable to substantially improve on our existing results. Naturally, since Medium Diff-through-Opt performs fewer hyperparameter updates, its performance may suffer if the hyperparameters have not converged, which we expect to be the likely cause here. However, it is reassuring that our methods are capable of performing well despite their short-horizon foundations.

heading: B.7.3 UCI ENERGY, FULL HORIZON
text: Finally, we consider UCI Energy, whose size is more amenable to full-horizon studies. Our new setting here, Full Diff-through-Opt, extends the Diff-through-Opt configuration from Section 4 by using a full look-back horizon of 4 000 steps, restarting training after each hyperparameter update and performing 30 such updates. This permits a full-horizon study on the same problem as before, using a similar computational footprint to our Bayesian Optimisation trials of Appendix B.5 (that is, 30 times the cost of our original experiments). We consider 36 random initialisations for this algorithm, and combine the results with our original UCI Energy experiments, showing both in Figure 18 and Table 13.
Curiously, while Full Diff-through-Opt outclasses naïve Random approaches, it is very uncompetitive with its short-horizon counterparts, achieving much higher average final losses (the larger variance being, of course, expected due to the smaller number of repetitions). The principal difference between Full Diff-through-Opt and Diff-through-Opt, apart from the different look-back distances, is the number of hyperparameter updates performed, since these are far more expensive in the former case. While more hyperparameter updates would likely cause Full Diff-through-Opt’s performance to improve, the computational impact would be massive, with such a setting unlikely to be desired in practical applications. Again, we conclude that it is simplistic to claim short-horizon bias renders an algorithm ill-suited to HPO — in reality, it is one of many factors worthy of consideration in the selection of an approach. In this case, the larger number of updates we are able to perform by sacrificing horizon length leads to favourable performance improvements, even though this optimisation is technically ‘biased’ as a result.

heading: B.8 UCI ENERGY: FURTHER ANALYSIS
text: In light of our main results, we perform several additional experiments to investigate particular aspects of our algorithm. These are conducted on the UCI Energy dataset, using the same configuration as in Section 4.1 except where otherwise noted.
B.8.1 VALIDATION DATASET SIZE
Noting that results for OursWD+HDLR+M did not always outperform OursWD+LR+M as might have been expected, and considering the risk of overfitting to validation data, we consider the role the validation dataset might have to play on final test performance. The conventional training/validation/test dataset split sizes we use throughout this paper produce a training set which is substantially larger than validation or test. In OurWD+HDLR+M case, we have approximately the same number of hyperparameters as model weights, which is not the case in classical training. Thus, we posit that choosing a larger validation dataset may permit improved hyperparameter tuning, addressing this performance gap.
In Table 14 we show results over 200 random initialisations, with bootstrapped error bars computed as in Section 4.1. We fix the test set size at our default 10% of the available data (per Gal & Ghahramani (2016)), and vary the balance between training and validation data as shown.
Both the scalar and high-dimensional settings clearly benefit from a larger validation dataset, with a 67.5%/22.5%/10% split significantly improving performance. But OursWD+HDLR+M seems to exhibit less reliable, higher-variance behaviour at even larger validation splits, with OursWD+LR+M also suffering worse performance in these settings. The differences between best-case performance in these algorithms suggest there is yet more to be understood about OurWD+HDLR+M setting. However, it is clear that with larger numbers of hyperparameters, we must take care to balance the competing demands for training and validation data.

heading: B.8.2 ASYNCHRONOUS HYPERBAND AND POPULATION-BASED TRAINING
text: Finally, we consider two further non-gradient-based HPO strategies: Asynchronous Hyperband (ASHA; Li et al. (2020b)) and Population-Based Training (PBT; Jaderberg et al. (2017)). Both approaches exploit parallel computation to consider a population of hyperparameter selections, and the performance of each individual selection during training is used to focus computational effort on the most promising trials. Broadly speaking, ASHA does this by terminating poorly-performing trials early, while PBT locally perturbs trials or replaces them entirely with better-performing checkpoints. In this sense, ASHA and PBT are orthogonal to the type of algorithm we propose in this paper, which seeks to improve the performance obtained by single initialisations trained in isolation. However, our algorithm can readily be combined with ASHA or PBT: by training each proposed hyperparameter choice using our algorithm instead of fixing the hyperparameters as is conventional, we exploit our method’s benefits for each individual run, thus providing ASHA/PBT with a stronger population of
trials from which to make optimisation decisions. We demonstrate this synergy in the experiments of this section.
Throughout, we make use of the Ray Tune (Liaw et al., 2018) library’s implementation of these strategies, utilising the default settings unless noted, and consider the UCI Energy dataset using the same setup as in Section 4.1. For the Random setting, we do not combine the training and validation datasets, instead training with the former and reserving the latter for the evaluation of different initialisations by ASHA/PBT. For OurWD+LR+M setting, the validation set is used both for gradientbased hyperparameter updates during training and for configuration evaluation by ASHA/PBT. Otherwise, both settings are the same as in Section 4, and again we quote final losses on a held-out test set, with 1 000 full-dataset bootstrap samples used to compute standard errors. We also leave the hyperparameter search space unchanged from before (see Section B.3).
Asynchronous Hyperband For ASHA, we seek to match the 4 000 × 200 weight update steps performed in Section 4.1, with 4 000 steps still required for full training. Under the default settings of Ray Tune (no grace period, reduction factor 4, one bracket), we obtain computational budgets of (1, 4, 16, 64, 256, 1 024, 4 000) iterations, with the number of initialisations surviving to each budget step correspondingly decreasing by a factor of 4 each time. Non-parallelised Hyperband (Li et al., 2017) uses the same total computational effort for each budget step; applying this fact to the smallest computational budget gives that we must begin with 4 000×2007×1 ≈ 114 286 initialisations for computational comparability. While the asynchronous modifications of ASHA introduce additional stochasticity, we assume the expected workload to be the same as the non-parallelised case. We run eight workers in parallel to mimic the experimental setting of our earlier studies.
From the results in Table 15, we see our hypothesis largely confirmed. By virtue of curating the population of initialisations, rather than considering each one in isolation, ASHA trivially gives improved average performance over all earlier methods. However, combining ASHA with OursWD+LR+M yields a significantly more performant population of trials on average, which leads to low-valued, low-variance mean and median final test performance. In both ASHA settings, the best final losses in the resulting populations are towards the lower end of those found in our previous results, but do not beat the especially small loss found by Random (3-batched). From this, we suppose that the methods we consider are nearing the performance limits of the UCI Energy dataset on this model when trained in this way. But regardless, we clearly see that deploying our algorithm in a particular context tends to improve average-case performance, indicating a positive influence on training.
Population-Based Training For PBT, we again attempt to match the earlier computational burden by training 200 initialisations for 4 000 iterations using 8-way parallelism, acknowledging that PBT’s checkpoint restoration will cause most runs to perform more than 4 000 weight update steps. We set the perturbation interval to 100 steps, so each initialisation has 40 opportunities to be mutated or reset by PBT. In PBT + OursWD+LR+M, it is unclear how to accommodate PBT’s hyperparameter mutations, as our algorithm will have updated the hyperparameters PBT associates with each trial. Where we detect PBT has mutated the hyperparameters, we compute the ratio between the new and old values as seen by PBT, then apply that same ratio to mutate the true current value of the hyperparameters.
Our PBT results (Table 15) largely continue the pattern we have discussed. Combining our algorithm with PBT achieves substantial improvements in mean and median performance, giving the lowest values seen in any of our experiments on UCI Energy. PBT uses additional internal machinery to manage the population of hyperparameter initialisations, and we found this caused the total training time to considerably increase over the non-population-based methods of Section 4.1, so the performance improvements we see here are perhaps unsurprising (see also Figure 19 and our discussion of runtimes below). However, in situations where we have already committed to the additional cost of PBT, replacing the fixed-hyperparameter internal training routine with our algorithm clearly results in a superior population of models, improving the quality of the models available for downstream applications. Note that the quality of the search over hyperparameter initialisations affects our overall best loss found, so some variation in the Best column of Table 15 is to be expected.
Runtime Analysis Figure 19 shows a runtime analysis of these results on ASHA and PBT. We see clearly that the superior performance of OursWD+LR+M in isolated training of specific hyperparameter initialisations causes ASHA to reach low losses faster. In PBT, there is little difference in performance
vs. runtime when Random or OursWD+LR+M is used as the internal learning strategy. Subjectively, this is likely because the overhead of saving checkpoints and restarting from them dominated the overall cost of running PBT. In both cases, PBT is substantially slower than ASHA to reach the most performant losses. Note by comparison with Figure 11a that our methods alone are competitive with ASHA when we compare the evolution of incumbent best performance over time. Note also that ASHA and PBT do not have access to the intermediate test losses plotted in Figure 19, so the smallest test losses found are not necessarily achievable as an output of these methods.

heading: C COMPLETE DERIVATION
text: Here, we reprise the derivations of Section 2 in more verbose form. Our subsections mirror those used previously.
C.1 IMPLICIT FUNCTION THEOREM IN BILEVEL OPTIMISATION
Recall that we consider some model parameterised by parameters w, which we train to minimise a training loss LT . Our training algorithm is governed by hyperparameters λ, which we simultaneously train to minimise a validation loss LV . We target optimal w∗ and λ∗, defined by the bilevel optimisation problem
(a) λ∗ = arg min λ LV (λ,w∗(λ)) , such that (b) w∗(λ) = arg min w LT (λ,w) . (1 revisited)
(1b) simply reflects a standard model training procedure, with fixed λ, so we need only consider how to choose the λ∗ from (1a).
Suppose we seek a gradient-based hyperparameter optimisation method for its performance potential. We then require the total derivative of the validation objective LV with respect to the hyperparameters λ, which may be obtained by the chain rule:
dLV dλ = ∂LV ∂λ + ∂LV ∂w∗ ∂w∗ ∂λ . (2 revisited)
An exact expression for ∂w ∗
∂λ can be obtained by applying Cauchy’s Implicit Function Theorem (Theorem 1). Suppose that for some λ′ we have identified a locally optimal w′, and assume continuous differentiability of ∂LT∂w . Then, there exists an open subset of hyperparameter space containing λ′, over which we can define a function w∗(λ) relating locally optimal w for each λ by
∂w∗
∂λ = − ( ∂2LT ∂w∂wT )−1 ∂2LT ∂w∂λT . (3 revisited)
Of course, the inverse Hessian term in (3) becomes rapidly intractable for models at realistic scale.

heading: C.2 APPROXIMATE BEST-RESPONSE DERIVATIVE
text: To proceed, we require an approximate expression for ∂w ∗
∂λ . Suppose model training proceeds iteratively by the rule wi(λ) = wi−1(λ)− u(λ,wi−1(λ)) , (4 revisited) where i is the current iteration and u is an arbitrary function differentiable in both its arguments. Note that the trivial choice
u(λ,wi−1(λ)) = wi−1(λ)− f(λ,wi−1(λ)) (12) allows arbitrary differentiable weight updates f to be subsumed by this framework.
Differentiating with respect to λ at the current hyperparameters λ′, we have:[ ∂wi ∂λ ] λ′ = [ ∂wi−1 ∂λ ] λ′ − [ du dλ ] λ′,wi−1(λ′)
(13)
= [ ∂wi−1 ∂λ ] λ′ − [ ∂u ∂w ∂wi−1 ∂λ + ∂u ∂λ ] λ′,wi−1(λ′)
(14)
= [( I− ∂u
∂w ) ∂wi−1 ∂λ − ∂u ∂λ ] λ′,wi−1(λ′)
(15)
= − [ ∂u
∂λ ] λ′,wi−1(λ′) + [ I− ∂u ∂w ] λ′,wi−1(λ′) [( I− ∂u ∂w ) ∂wi−2 ∂λ − ∂u ∂λ ] λ′,wi−2(λ′)
(16)
= − ∑
0≤j<i  ∏ 0≤k<j [ I− ∂u ∂w ] λ′,wi−1−k(λ′) [∂u ∂λ ] λ′,wi−1−j(λ′)  (5 revisited)
To combine terms in (5), we must make approximations such that all bracketed expressions are evaluated at the same λ,w. For network weights w, we set w0 = wi = w∗ for all iterations i, effectively assuming we have been at the optimal w for some time; when we come to consider finite truncations of the summation, this approximation will be especially appealing. With this approximation, we may collapse the product and write[
∂wi ∂λ ] λ′ ≈ − ∑ 0≤j<i  ∏ 0≤k<j ( I− ∂u ∂w ) ∂u ∂λ  λ′,wi(λ′)
(17)
= −  ∑ 0≤j<i ( I− ∂u ∂w )j ∂u ∂λ  λ′,wi(λ′) . (6 revisited)
As i increases, our training algorithm should converge towards a locally optimal w∗, and the summation in (6) will accrue more terms. Instead, we reinterpret i as specifying the finite number of terms we consider as an approximation to the whole sum. With this redefined i, we arrive at our approximate best-response Jacobian[
∂w∗
∂λ ] λ′ ≈ −  ∑ 0≤j<i ( I− ∂u ∂w )j ∂u ∂λ  λ′,w∗(λ′) . (6 revisited)

heading: C.3 CONVERGENCE TO BEST-RESPONSE DERIVATIVE
text: To informally argue the convergence properties of our method, note that (6) is a truncated Neumann series: [
∂wi ∂λ ] λ′ ≈ −  ∑ 0≤j<i ( I− ∂u ∂w )j ∂u ∂λ  λ′,w∗(λ′) . (18)
Analogously to its scalar equivalent — the geometric series — this series converges as i→∞ if the multiplicative term (I− ∂u∂w ) is ‘contractive’ in some sense. For our purposes, it suffices to assume a Banach space and, for the operator norm ‖·‖, that ‖I− ∂u∂w‖ < 1. In this case, we can directly apply the closed-form limit expression for the Neumann series:
lim i→∞ [ ∂wi ∂λ ] λ′ ≈ − lim i→∞  ∑ 0≤j<i ( I− ∂u ∂w )j ∂u ∂λ  λ′,w∗(λ′)
(19)
= − [( I− ( I− ∂u
∂w
))−1 ∂u
∂λ ] λ′,w∗(λ′)
(20)
= −
[( ∂u
∂w
)−1 ∂u
∂λ ] λ′,w∗(λ′) . (7 revisited)
Note that (7) is exactly the result of the Implicit Function Theorem (3) applied to our more general weight update u. Concretely, suppose for some λ′′ and w′′ that u(λ′′,w′′) = 0; in words, that our training process does not update the network weights at this point in weight-hyperparameter space, and assume u is continuously differentiable with invertible Jacobian. Then, over some subset of hyperparameter space containing λ′′, we can say: w∗(λ) exists, we have u(λ,w∗(λ)) = 0 throughout the subspace and ∂w ∗
∂λ is given by (7) exactly.

heading: C.4 SUMMARY OF DIFFERENCES FROM LORRAINE ET AL. (2020)
text: In their derivation, which ours closely parallels, Lorraine et al. (2020) replace (4) with a hard-coded SGD update step; in our notation, they fix u(λ,w) = uSGD(λ,w) = η ∂LT∂w . They then isolate the
learning rate η, insisting all hyperparameters λ must be dependent variables of the training loss LT , before further differentiating the training loss in their derivation. Consequently, their algorithm cannot optimise optimiser hyperparameters such as learning rates and momentum factors, since these come into play only after the training loss has been computed and differentiated. While a careful choice of transformed loss function may allow certain optimiser hyperparameters to be handled by the framework of Lorraine et al. (2020) (e.g. setting L = ηLT to incorporate a learning rate), this manual rederivation is far from trivial in more complex cases (e.g. to incorporate momentum).
Instead, we define a general weight update step (4), without assuming any particular form. We may then differentiate the general update u directly, encompassing any hyperparameters appearing within it but outside the training loss, which Lorraine et al. (2020) would miss. Crucially, this includes the hyperparameters of most gradient-based optimisers, including SGD. Thus, we lift Lorraine et al.’s results from expressions about ∂LT∂w to expressions about u, then update hyperparameters inaccessible to the former. Concretely: we optimise weight decay coefficients, learning rates and momentum factors in settings where Lorraine et al. can only manage the weight decay alone.

Review summaries:
This paper tackles the hyperparameter optimization problem with a one-pass approach that alternatively optimizes over machine learning model parameters and hyperparameters. To optimize over a bilevel optimization problem, this paper generalizes Lorraine et al. (2019) by replacing the gradient update of parameters with a more general parametric form, thus allowing the hyperparameters to be updated alongside parameters by gradients coming from truncated Neumann series. 

The paper aims to accelerate HPO optimization by allowing for retraining when new hyperparameters are proposed, while being faster and more flexible that existing hypergradient-based one-pass methods. This results in a novel HPO one-pass technique that works with any continuous hyperparameter within a differentiable model weight update, such as learning rates. In a good range of experiments, the authors show that their proposal is faster than baselines.

This paper is concerned with tuning continuous hyperparameters of a neural network setup (model, loss function and optimizer parameters) by way of a cheap approximation to the hypergradient. It is closely related to earlier work (Lorraine, 2019), which however overlooked that the idea can be generalized to optimizer parameters (e.g., learning rate, momentum), which is effectively done here. In a range of well-done experiments, the method is shown to very convincingly outperform (Lorraine, 2019), and in some cases even the "gold standard" approach to differentiate through the whole learning trajectory. Code is provided.