Year: 2021

Abstract:
Stochastic gradient descent is commonly employed as the most principled optimization algorithm for deep learning, and the dependence of the generalization error of neural networks on the given hyperparameters is crucial. However, the case in which the batch size and learning rate vary with time has not yet been analyzed, nor the dependence of them on the generalization error as a functional form for both the constant and time-varying cases has been expressed. In this study, we analyze the generalization bound for the time-varying case by applying PAC-Bayes and experimentally show that the theoretical functional form for the batch size and learning rate approximates the generalization error well for both cases. We also experimentally show that hyperparameter optimization based on the proposed model outperforms the existing libraries.

heading: 1 INTRODUCTION
text: While deep learning addresses a variety of practical tasks, the reason deep learning generalizes well is not fully understood. It is hoped that the performance of deep learning methods can be improved by understanding deep learning. Consequently, deep learning will become a cornerstone for the next generation of technology. Stochastic gradient descent (SGD) is a principled learning algorithm for deep learning, and because most deep learning algorithms are developed based on SGD, its theoretical analysis is crucial. He et al. (2019) modeled the generalization bounds when both the batch size and learning rate are constant in SGD; however, they did not model the generalization error. Furthermore, it is crucial to decrease the learning rate as the training progresses to make SGD converge (Robbins & Monro, 1951), which results in good performance (You et al., 2020), but He et al. (2019) did not consider the case in which these parameters are time-varying. We analyze the time-varying case in SGD and propose a generalization error model for both constant and time-varying cases. Based on these models, we also propose a novel hyperparameter search method. In this study, we utilized a PAC-Bayes-based generalization bound to generate a model of the generalization error. To validate the proposed model, we trained approximately 300 networks on both practical networks (VGG16 (Simonyan & Zisserman, 2015), Wide-ResNet 28-10 (Zagoruyko & Komodakis, 2016)) and data (CIFAR10, CIFAR100 (Krizhevsky et al., 2009)). We also showed that the proposed model can approximate the experimental results well. We also show that the proposed model can be used for hyperparameter search via Bayesian optimization and can yield better results, compared with the previous method for searching in the uniform distribution, and comparable results for searching in the logarithmic uniform distribution.
The contributions of the present study are summarized as follows:
• We analyzed the generalization bound using PAC-Bayes with the timevarying batch size and learning rate.
• We directly modeled the generalization error on the batch size and learning rate and proposed its functional form.
• We experimentally showed that it is possible to model the generalization error with high accuracy for practical data.
• Our experimental results demonstrated that the hyperparameter search method built on the proposed model can outperform current libraries when
searching over uniform distributions and has comparable performance when searching over logarithmic uniform distributions.

heading: 2 RELATED WORK
text: In this section, we present related work on SGD modeling, generalization bounds, prediction of model performance, and hyperparameter optimization.
Theoretical analysis of SGD: SGD is often treated as a continuous-time stochastic differential equation (SDE) by modeling the noise of the gradient and the loss function and applying the Euler– Maruyama approximation (Lin et al., 2018; Mandt et al., 2017; Li et al., 2017; Jastrzębski et al., 2017). This approximation can be solved analytically as an Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930). In this study, we also analyzed SGD using the same process as the previously mentioned studies. However, in most previous studies, the batch size and learning rate in SGD were treated as constants. Hence, in this study, we treat them as time-varying and seek an analytic solution.
Theoretical error bounds: There is a large body of research on theoretical explanations for the generalization capability of deep neural networks (Negrea et al., 2019; Arora et al., 2018; Cao & Gu, 2019; Deng et al., 2021). He et al. (2019), building on their analysis of SGD, provided a theoretical explanation of the effects of the batch size and learning rate on the generalization error. However, this research only handles the case of constant batch size and learning rate and does not model the generalization error in the function form. In this study, we derive a generalization bound based on an analytic solution for the time-varying case and derive a model of the generalization error.
Predicting model performance: Most studies are based on the results from training networks, assuming hyperparameters as inputs, and searching for the best parameters without training the networks. In neural architecture search, the performance of a model was predicted for a given model on a given dataset (Istrate et al., 2019), small data (Klein et al., 2017), and model size (Real et al., 2019), respectively. Several studies predict the performance of a model in a functional form by treating both the dataset and model sizes as input (Rosenfeld et al., 2020), or the pruning scale as input (Rosenfeld et al., 2021). In this study, we exploit the findings of the above studies to predict the performance of a model in a functional form, based on a theoretically derived generalization bound to the batch size and learning rate.
Hyperparameter optimization: Hyperparameter optimization algorithms are mainly categorized into two types: nonadaptive and adaptive. Nonadaptive algorithms search for the best hyperparameters without using any past search results and, thus, have the advantage of being easily parallelizable. Grid search and random search fall into this category, and many studies have shown that random search is more useful (Bergstra & Bengio, 2012; Li & Talwalkar, 2020). Adaptive algorithms search for the best hyperparameters by exploiting past search results and can obtain better hyperparameters with fewer searches. Adaptive algorithms include specific algorithms, such as the evolutionary algorithm (Young et al., 2015) and Bayesian optimization (Wu et al., 2019), and practical tools, such as Hyperopt (Bergstra et al., 2015) and Optuna (Akiba et al., 2019). These methods are flexible enough to be optimized for arbitrary functions. However, there is still room for more efficient search when searching only for specific hyperparameters. In this study, we aim to further improve the search efficiency of both the batch size and learning rate by employing a hyperparameter search based on the proposed generalization error model.
In the following sections, we will analyze the influence of the batch size and learning rate on the generalization error in SGD, based on assumptions about the distribution of the gradients and the shape of the risk function near minima, an approximation to SDE, and the PAC-Bayes method.

heading: 3 PRELIMINARIES
text: In this section, we describe the preliminary knowledge used to model the learning process of SGD.
Generalization error of stochastic algorithms. In general, machine learning algorithms are designed to select a hypothesis function Fθ from a hypothesis class {Fθ|θ ∈ Θ ⊂ Rd} such that the expected value of the risk function R computed from the loss function l is minimized. In a proba-
bilistic algorithm such as SGD, the parameter θ after learning is considered to lie on the distribution Q. The risk function for parameter θ is defined as follows:
R(θ) = E(X,Y )∼D [l(Fθ(X), Y )] . (1)
Therefore, R(Q) = Eθ∼Q [R(θ)]. However, one cannot compute the risk function, and in practice, the empirical risk R̂(θ) from the data is approximated as the risk functionR(θ).
R̂(θ) = 1 N N∑ i=1 l(Fθ(Xi), Yi), (2)
where (Xi, Yi) is the N training data consisting of the input data and label pairs. Thus, R̂(Q) = Eθ∼Q[R̂(θ)]. Stochastic gradient descent. Gradient descent is an algorithm that computes the gradient of the risk function defined in Eq.1 and updates the parameter θ based on the gradient. The gradient g(θ(t)) is defined as g(θ(t)) , ∇θ(t)R(θ(t)) = ∇θ(t)E(X,Y ) [ l(Fθ(t)(X), Y ) ] . We update the parameter θ by θ(t + 1) = θ(t) − ηg(θ(t)), where η is the learning rate. In practice, the true gradient of g(θ(t)) is approximated by the gradient of ĝS(θ(t)) computed from the empirical risk R̂(θ) defined in Eq.2, and the parameter θ is updated with ĝS(θ(t)). Thus, the gradient ĝS(θ(t)) is defined as ĝS(θ(t)) = ∇θ(t)R̂(θ(t)) = 1|S| ∑ i∈S ∇θ(t)l(Fθ(t)(Xi), Yi), and the parameter θ is updated by θ(t+1) = θ(t)−ηĝS(θ(t)). We use a mini-batch S to train the model in parallel, and |S| represents its size. Then, we treat li(θ) = l(Fθ(t)(Xi), Yi).
Analysis of the learning process. The loss function li(θ) and its mean R̂(θ) is an unbiased estimator of the risk functionR(θ), and the gradient of the loss function ∇θli(θ) and its mean ĝS(θ) is an unbiased estimator of the gradient g(θ) = ∇θR(θ). Hence, we have E[li(θ)] = E[R̂(θ)] = R(θ) and E[∇θli(θ)] = E[ĝS(θ)] = g(θ) = ∇θR(θ), where the expectations are computed from the data distribution (X,Y ). Here, as in Mandt et al. (2017), we assume that {∇li(θ)} is independently and identically sampled from a Gaussian distribution with mean g(θ) = ∇θR(θ). We can express ∇θli(θ) ∼ N (g(θ), C), where C is the covariance matrix, which is assumed to be a constant matrix for all θ. Then the gradient of ĝS(θ) computed for a mini-batch follows the following Gaussian distribution:
ĝS(θ) = 1 |S| ∑ i∈S ∇θli(θ) ∼ N ( g(θ), 1 |S| C ) . (3)
Because the covariance matrix is symmetric and (semi) positive definite, we suppose that C can be factorized as C = BB>. In an actual SGD algorithm, the parameters are iteratively updated using Eq.3 to minimize the risk functionR(θ):
∆θ(t) = θ(t+ 1)− θ(t) = −ηĝS(θ(t)) = −ηg(θ) + η√ |S| B∆W, ∆W ∼ N (0, I). (4)
Analytical solution of the stochastic differential equation. When the learning rate η is small, Eq.4 can be approximated by SDE, which is known as Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930). Here, we assume that the risk function is convex and second-order differentiable in the neighborhood of the minimum (He et al., 2019). In addition, when the risk function R(θ) is minimized at θ = θ∗, it is shifted so that it is minimized at θ = 0. Thus, the risk function R(θ) can be approximated by second-order statistic on the distribution Q as follows:
R(θ) = 1 2 θ>Aθ, (5)
whereA is the Hessian matrix around the minimum of the risk function and is assumed to be positive definite. This assumption was experimentally demonstrated by Li et al. (2018), and other studies have made similar assumptions (Jastrzębski et al., 2017; Poggio et al., 2017). Because E[R̂(θ)] = R(θ), we further assume that the empirical risk R̂(θ) is in the following quadratic form:
R̂(θ) = 1 2 (θ − θb)>Â(θ − θb) (6)
where θb is the bias term, Â is the Hessian matrix, and E[θb] = 0, E[Â] = A.
It is well known that Eq.4 has an analytic solution (Gardiner, 2004),
θ(t) = e−Atθ(0) +
√ η
|S| ∫ t 0 e−A(t−t ′)BdW (t′), (7)
where W (t′) is white noise and followsN (0, I) and Q is the distribution of θ in Eq.7 and Gaussian.

heading: 4 THEORETICAL GENERALIZATION BOUND
text: In this section, we derive the generalization bounds for both constant and time-varying cases.

heading: 4.1 CONSTANT CASE
text: Using Eq.7 and PAC-Bayes, we obtain a generalization bound for SGD as follows (see Appendix B.1 for the derivation).
Theorem 1 (extension of He et al. (2019), Theorem 2). We treatA as the Hessian of the risk function around the local minimum, C as the covariance matrix of the gradients calculated by single sample points, Q as the distribution of the output hypothesis function of SGD, Σ is the covariance matrix of the distribution Q, d is the dimension of the parameter θ (network size), and R is the search radius of the parameter θ.
When A and Σ commute, for any positive real δ ∈ (0, 1), with probability 1 − δ over a training sample set of size N, we have the following inequality for the distribution Q:
R(Q) ≤ R̂(Q)
+ √√√√ η2|S| tr(CA−1) + d log ( 2|S|η )− log (det (CA−1)) +R2 − d+ 2 log ( 1δ ) + 2 logN + 4 4N − 2 .
(8)
This generalization bound incorporates a new scale factor for the parameter search, R2, into Theorem 2 of He et al. (2019). Although similar results can be derived without assuming that A and Σ commute, we used the above assumption for simplicity. R represents the radius of the hypersphere in the parameter search, and ‖θ − θ0‖2 ≤ R2. We also show that this generalization bound also has a network initialization factor. A detailed description is provided in Appendix B.2.

heading: 4.2 TIME-VARYING CASE
text: It should be noted that, when solving Eq.4, the hyperparameters in SGD such as the batch size |S| and the learning rate η can be time-varying. While He et al. (2019) derived a generalization bound for the case in which these are constant, we further analyze the time-varying case. In the time-varying case, the solution takes the form of a hyperparameter term in Eq.7 entering the integral (Gardiner, 2004).
θ(t) = e−Atθ(0) + ∫ t 0 √ n(t′)e−A(t−t ′)BdW (t′), (9)
where n(t) = η(t)|S(t)| . Then, we treat the time-varying batch size and learning rate as a univariate function n(t). This allows us to treat changes in the learning rate (Robbins & Monro, 1951; You et al., 2020) and batch size (Smith et al., 2018) in a unified manner. In general, if n(t) is allowed to be an arbitrary function, this integral cannot be solved analytically. However, because the learning rate generally decreases and the batch size generally increases as the learning progresses, if we assume n(t) to be a monotonically decreasing function and apply the mean value theorem for integrals∫ t1 t0 n(t)F (t)dt = n(ξ) ∫ t1 t0 F (t)dt, ξ ∈ [t0, t1], then we obtain a generalization bound for the timevarying case as follows (see Appendix B.3 for the detailed derivation).
Theorem 2. We treat D = ∫ t
0 e−A(t−t ′)Ce−A(t−t ′)dt′. Using the same settings as theorem 1 and
when A and Σ commute and A and D commute, we have the following inequality for the timevarying case:
R(Q) ≤ R̂(Q)
+ √√√√ ninter2 tr(CA−1) + d log ( 2ninter )− log (det (CA−1)) +R2 − d+ 2 log ( 1δ ) + 2 logN + 4 4N − 2 ,
(10)
and ninter = (1− λ) min (n(t)) + λmax (n(t)), λ ∈ [0, 1]. (11)

heading: 5 FUNCTIONAL APPROXIMATION OF THE GENERALIZATION ERROR
text: Based on the theoretical analysis thus far, we consider approximating R(Q) in a functional form. From Eq.8 and Eq.10, we have already modeled the generalization bound (R(Q) − R̂(Q)). We can also design a empirical functional form based on He et al. (2019)’s experimental results on the effect of ratio of the batch size to learning rate on the generalization error. Thus, if we can model the training error R̂(Q), we can derive the functional form of the generalization error. We model the training error by substituting Eqs.7 and 9 into Eq.6. Subsequently, we show the specific functional form of the generalized error model for constant and time-varying cases.

heading: 5.1 CONSTANT CASE
text: First, we organize Eq.7 with respect to the batch size |S| and learning rate η and the learned parameter θ(t) is as follows:
θ(t) = a0 +
√ η
|S| a1, (12)
where a0,a1 are vectors that do not depend on |S|, η. By substituting this into Eq.6, the functional form of the training error R̂(Q) is as follows. The detailed derivation is in Appendix B.4.1.
R̂(Q) ' Eθ∼Q [ 1
2 (θ − θb)>Â(θ − θb)
] = E [ 1
2
( a′0 + √ η
|S| a1
)> Â ( a′0 + √ η
|S| a1
)]
= a0 + a1 η
|S| + a2
√ η
|S| , (13)
where a′0 = a0 − θb, and a0, a1, a2 are constants. Although we can model R(Q) directly by substituting Eq.7 into Eq.5, this cannot be treated as a practical functional form as it does not model the experimental increase of the generalization bound with respect to the decrease of η|S| (He et al., 2019). Thus, we address this problem by modeling the training error R̂(Q) and the generalization bound (R(Q) − R̂(Q)) separately. To model the generalization bound, we propose the following functional form based on the model derived in Eq.8 and He et al. (2019)’s experimental results. The detailed derivation is given in Appendix B.5.1.
R(Q)− R̂(Q) ' a3 + a4 √ |S| η
(14)
Here, a3, a4 are constants. By combining Eq.13 and Eq.14, we propose the following functional form as a model of the generalization errorR(Q).
R(Q) = ̃(S, η) = ( c0 η
|S| + c1
√ η
|S|
) + ( c2 √ |S| η ) + c3 (15)
where c0 ∼ c3 are constants. The first term models the training error R̂(Q), and the second term models the generalization bound (R(Q)−R̂(Q)). Thus, the first term denotes that the training error monotonically increases with respect to the hyperparameter η|S| , and the second term denotes that the generalization bound inversely decreases.

heading: 5.2 TIME-VARYING CASE
text: As discussed in Section 4.2, we approximate Eq.9 to a simple form similar to Eq.12, using an approximation of the mean value theorem for integrals. Based on the functional form obtained by substituting Eq.9 into Eq.6, the model of generalization bounds derived in Eq.10, and He et al. (2019)’s experimental results, we propose the following functional form as the model of the generalization error (S, η). A detailed derivation is given in Appendix B.4.2 and B.5.2.
̃(S, η) = ( c0n 2 1 + c1n1 ) + ( c2 √ 1
n2
) + c3 (16)
n1 = (1− λ1) min (√ n(t) ) + λ1 max (√ n(t) )
(17)
n2 = (1− λ2) min (n(t)) + λ2 max (n(t)) (18)
where c0 ∼ c3, λ1, λ2 are constants and λ1, λ2 ∈ [0, 1]. Note that, for n1 in the first term, the intermediate value is computed after taking the square root of n(t), and for n2 in the second term, the linear interpolation is computed.

heading: 5.3 PRACTICAL FUNCTIONAL FORM
text: In the discussion above, we constructed a rough model of the generalization error. However, the models in Eqs.15 and 16 have no upper bounds, and they don’t model the fact that the loss function
hits its head at the random guess level. Therefore, it is necessary to model the transition from the initial random guess level to the proposed model ̃(S, η). Thus, we propose the model ̂(S, η) using a soft-min function based on the LogSumExp function.
̂(S, η) = softmin(̃(S, η), 0; c4) = 1
c4
( − log ( e−c4̃(S,η) + e−c4 0 )) (19)
where c4 is a constant that controls the shape (sharpness) of the soft-min function (softmin(x, y) = − log (exp (−x) + exp (−y))). 0 is the value of the loss function at the random guess level, which is statistically determined and not explored. (e.g., for a balanced dataset in image classification, we can treat it as 0 = 1− (1/Nclasses).)

heading: 6 EXPERIMENT
text: To validate the proposed model ̂(S, η), we trained a practical network (VGG16 (Simonyan & Zisserman, 2015), Wide-ResNet 28-10 (Zagoruyko & Komodakis, 2016)) on the image classification tasks (CIFAR10, CIFAR100 (Krizhevsky et al., 2009)).

heading: 6.1 EXPERIMENTAL SETTINGS AND RESULT
text: In the constant case, the batch size and learning rate were fixed throughout the training process. To evaluate the generalization error, we trained 60 VGG models with a batch size |S| ∈ {16, 32, 64, 128, 256, 512} and a learning rate η ∈ {1/2i|i = 3, ..., 12} for CIFAR10, as well as 66 Wide-ResNet models with the same batch size |S| and learning rate η ∈ {1/2i|i = 2, ..., 12} for CIFAR100 over a total of 200 epochs. The experimental results for the constant case are shown in Figures.1a and 1b. We plotted the log scale of η|S| on the horizontal axis and the generalization error on the vertical axis. We found both a decreasing and an increasing property of the generalization error before and after η|S| ' 10
−4, which we discussed in Eq.15. Although He et al. (2019) claimed that increasing η|S| decreases the generalization error, the experimental results show that the generalization error increases when η|S| exceeds a certain value. This implies that, in addition to the generalization bound, the training error R̂(Q) must be modeled. In the time-varying case, only the learning rate is set to be time-varying, and we use an exponential step function. We chose ηinit and ηfinal and trained over 200 epochs, starting with a learning rate of ηinit, decaying with a decay rate of γ at 60, 120, and 160 epochs, and finally reached a learning rate ηfinal = γ3ηinit at 160 epochs.
In addition to the constant case, we trained a case in which the batch size |S| ∈ {16, 32, 64, 128, 256}, the learning rate of ηinit ∈ {1/2i|i = 3, ..., 4} and ηfinal ∈ {0.1/2i|i = 0, ..., 4} for CIFAR10, the same batch size |S|, the learning rate of ηinit ∈ {1/2i|i = 2, ..., 4} and the same ηfinal for CIFAR100, as well as a case in which the batch size |S| ∈ {128, 256}, the learning rate of ηinit ∈ {1/2i|i = 5, ..., 9} and ηfinal ∈ {0.1/2i|i = 5, ..., 9} for both datasets. In
total, we measured the generalization error for the 153 models for CIFAR10 and 184 models for CIFAR100. More implementation details can be found in Appendix C.1. The experimental results for the time-varying case are shown in Figures.1d and 1e. ηinit|S| and ηfinal |S| are plotted on the horizontal axis, whereas the generalization error (S, η) is plotted on the vertical axis. We also plot the results of projecting the generalization error onto ηinit|S| and ηfinal |S| , respectively. As in the constant case, we can observe that, when the value of n(t) = η(t)|S(t)| is not in the appropriate range (either quite large or quite small), the generalization error becomes large.

heading: 6.2 ERROR LANDSCAPE ESTIMATION
text: We estimate the obtained generalization error (S, η) using a parametric family, ̂(S, η;φ). Here, we treat the parameter φ as φ = {c0, c1, c2, c3, c4} for the constant case, and then {λ1, λ2} is added to the time-varying case. As in Rosenfeld et al. (2021), we search for φ such that the squared error of δ(S, η;φ) = ̂(S,η;φ)− (S,η) (S,η) is minimized. Implementation details are provided in Appendix C.2.
Figure.1c plots our results for the constant case. The mean µ and standard deviation σ of the relative error δ are |µ| < 2% and σ < 12%, respectively. Considering the small number of parameters used to estimate the generalization error (|φ| = 5), the proposed functional model can accurately estimate the complex generalization error landscape (60 and 66 models, respectively). In Figure.1f, we plot the results for the time-varying case. We achieved |µ| < 3% and σ < 17%, respectively. Although the relative error is larger than the constant case, we can still declare that the proposed functional model can estimate the generalization error sufficiently, considering that the generalization error of 153 or 184 models is modeled with a small number of parameters (|φ| = 7). In Appendix C.2.3, we provide a comparison between the proposed model, a model that analytically solves the time-varying case of the step function, and a model that uses only the final learning rate.
In addition, we verified the stability of the proposed model with respect to the number of networks. We randomly sampled a certain number of generalization errors obtained by training the networks, estimated the parameter φ from them only, and computed the relative error δ. We performed this procedure 100 times and computed µ and σ of the relative error. Details of the experimental setup are presented in Appendix C.3.1.
In Figure.2 we show the experimental results. Figures.2a and 2c show that, for the constant case of CIFAR10 and CIFAR100, we can accurately estimate the generalization error landscape by randomly training about 30 models. In contrast, Figures.2b and 2d show that, for the time-varying case, we need to train about 75 models for the estimation.

heading: 7 HYPERPARAMETER OPTIMIZATION
text: In this section, we show that the proposed model is useful for applications such as hyperparameter optimization. In this study, we used Bayesian optimization.
Bayesian optimization can compute a function y = f(x) that can be fitted to any sequence of data points {(xi, yi)}. Nogueira (2014) achieved this by using a Gaussian process with a Matern kernel, and many Bayesian optimization methods use generic kernel functions which can be fitted to any function. While this is practical for fitting unknown functions, ideally, the most accurate Bayesian optimization can be achieved by employing a kernel function constructed from a feature extraction function ψ(x) that can compute the unknown function with a linear combination of y = w>ψ(x). Therefore, for the constant case, we propose a kernel function k(xi, xj) for the input x = η|S| based on Eq.15 for hyperparameter optimization.
k(xi, xj) = ψ(xi) >ψ(xj) = √ xixj + xixj + 1 √ xixj + 1. (20)
Details are in Appendix C.4.1. In contrast, the time-varying case cannot be represented by a linear combination of y = w>ψ(x), because the product term of the parameters (c and λ) appears. Therefore, we employ a hyperparameter search by sequential model-based global optimization (SMBO) (Bergstra et al., 2011) with the proposed model in Eq.19 using Eq.16. For SMBO, we use the Levenberg–Marquardt algorithm (Moré, 1978) and Thompson Sampling (Thompson, 1933; Russo et al., 2017). More details are provided in Appendix C.4.2.
To validate the proposed method, we explored η|S| using the fitted generalization error model ̂(S, η;φ) as the training result of the network. The details are given in Appendix C.4.3. Figure.3 shows the experimental results. The horizontal axis plots the number of training networks, and the vertical axis plots the maximum test accuracy among the networks trained to date. A total of 100 experiments were conducted, and the 95% confidence interval was shaded. The methods used for comparison were Bayesian optimization using the Matern kernel and RBF kernel (Nogueira, 2014), Hyperopt (Bergstra et al., 2015), and Optuna (Akiba et al., 2019).
As shown in the Figure.3, in almost all cases, the proposed method achieves the highest accuracy with the least number of models in the search method with uniform distribution. This indicates that the proposed model represents a good feature extractor, and the proposed kernel function is useful for hyperparameter optimization. In addition, in the constant case, the proposed method achieved the highest accuracy with the same number of models as the comparison method, which searches for logarithmic uniform distributions. It indicates that theoretical model-based hyperparameter optimization can be applied practically. In practice, it is generally believed that it is better to search for hyperparameters S, η in the logarithmic uniform range (Rasmussen, 1997; Sundararajan & Keerthi, 2001). However, there is no theoretical foundation for this (this is confirmed by the fact that in Hyperopt, the results are not improved by exploring the logarithmic range). These results suggest that the proposed model-based hyperparameter optimization is both theoretically sound and practically useful when performing a hyperparameter search on the batch size and learning rate.

heading: 8 DISCUSSION AND CONCLUSION
text: In this paper, we presented a theoretical analysis and a model of the generalization error for the batch size and learning rate in SGD. We analyzed the case in which these parameters are constant during training and the more challenging and practical case when they vary with time. Through experiments for both cases, we show that the proposed generalization error model is practical, and we test the number of networks required to estimate the generalization error landscape (= stability of the generalization error model). In addition, we demonstrate that a hyperparameter search based on the proposed model outperform the conventional Bayesian optimization and hyperparameter optimization libraries, such as Hyperopt and Optuna. To the best of our knowledge, this is the first study that explores the batch size and learning rate based on the generalization error model. From the discussion, we conclude that this theoretical model is a practical and useful model for understanding the effects of the batch size and learning rate on the generalization error. We hope that this work will provide further theoretical insights into deep learning and new research directions for hyperparameter search based on theoretical models.

heading: A BACKGROUND

heading: A.1 PAC-BAYESIAN FRAMEWORK
text: The PAC-Bayesian framework originates from McAllester (1999a;b)’s work. In the PAC-Bayesian perspective, the hypothesis functions learned by probabilistic machine learning algorithms are drawn at random from the hypothesis classes based on some “rules.” The generalization ability of an algorithm is negatively correlated with the distance (often measured by the Kullback–Leibler (KL) divergence) between the output hypothesis distribution and the a priori distribution (typically a Gaussian or uniform distribution). The above results suggest that there is a tradeoff between minimizing the empirical risk and exploring further regions of the hypothesis space from the beginning (a priori).
Let P be the initial priori distribution in the parameter space Θ, and let Q be the distribution in Θ after learning. We can now define the expected value of the risk function under distribution Q as follows:
R(Q) = Eθ∼Q[R(θ)].
Similarly, the empirical risk function under distribution Q can be defined as follows:
R̂(Q) = Eθ∼Q[R̂(θ)].
Based on these, we obtain the classical result in that the expected value of the risk functionR(Q) is uniformly bounded by the expected value of the empirical risk function R̂(Q) and KL-divergence D(Q‖P ). Lemma 1 (See McAllester (1999a), Theorem1). For any positive real δ ∈ (0, 1) with a probability of at least 1− δ over a sample of size N , we have the following inequality for all distributions Q:
R(Q) ≤ R̂(Q) +
√ D(Q‖P ) + log ( 1δ ) + logN + 2
2N − 1 , (21)
where D(Q‖P ) is the KL divergence between the distributions Q and P and is defined as D(Q‖P ) = Eθ∼Q ( log Q(θ)
P (θ)
) . (22)

heading: B SUPPLEMENTAL PROOFS FOR THE PROPOSED METHOD
text: In this section, we describe the details of the proof, which we omitted in the proposed method.

heading: B.1 GENERALIZATION BOUND FOR THE CONSTANT CASE
text: To prove Theorem 1, we use an existing proof method that approximates the SGD update of the model to a stochastic differential equation (see, e.g., (Mandt et al., 2017)). With appropriate assumptions, the learning process of SGD can be approximated by the Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930). Because the Ornstein–Uhlenbeck process has an analytic solution, we can use the distribution Q of the analytic solution θ to derive the KL-divergence, which is used to derive a generalization error bound from Eq.21. Here, we modify and extend He et al. (2019)’s Theorems 1 and 2 to derive the generalization error bound for the constant case. Lemma 2 (cf. Gardiner (2004), pp.109-110). Under second-order differentiable assumption (Eq.5), if t is sufficiently large, the Ornstein–Uhlenbeck process (Eq.4)’s distribution
q(θ) = M exp { −1
2 θ>Σ−1θ
} , (23)
has the following property:
AΣ + ΣA = η
|S| C. (24)
This lemma is similar to that of Gardiner (2004).
Proof of Lemma 2. From the analytical solution of the Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930), the parameter θ can be expressed as follows:
θ(t) = e−Atθ(0) +
√ η
|S| ∫ t 0 e−A(t−t ′)BdW (t′), (25)
where W (t′) is a white noise and follows N (0, I). From Eq.23, we obtain Σ = Eθ∼Q [ θθ> ] . (26)
Thus, we obtain the following equation:
AΣ + ΣA = Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ η
|S| ∫ t 0 Ae−A(t−t ′)Ce−A(t−t ′)dt′
+ η
|S| ∫ t 0 e−A(t−t ′)Ce−A(t−t ′)Adt′
= Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ η
|S| ∫ t 0 d dt′ e−A(t−t ′)Ce−A(t−t ′)dt′
= Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ η |S| C − η |S| e−AtCe−At
' η |S| C, (27)
where A is the Hessian matrix and symmetric. The last transformation is guided by the assumption that A is a positive definite matrix and all its eigenvalues are positive and by the approximation of the term multiplied by e−At twice with zero, because t is sufficiently large.
The proof is completed.
Lemma 3 (extension of He et al. (2019), Appendix B.1). For any positive real δ ∈ (0, 1), with probability 1 − δ over a training sample set of size N, we have the following inequality for the distribution Q of the output hypothesis function of SGD:
R(Q) ≤ R̂(Q)+
√ η |S| tr(CA −1)− 2 log (det(Σ)) + 2R2 − 2d+ 4 log ( 1 δ ) + 4 logN + 8
8N − 4 , (28)
where A is the Hessian of the risk function around the local minimum, C is the covariance matrix of the gradients calculated by single sample points, Σ is the covariance matrix of the distribution Q, d is the dimension of the parameter θ(network size), and R is the search radius of the parameter θ.
Proof of Lemma 3. In the PAC-Bayesian framework ((Lemma 1), it is crucial to compute the KLdivergence between the parameter distribution of the output hypothesis function and the prior parameter distribution in the hypothesis space. The prior distribution can be interpreted as an initial parameter distribution, and in most cases, either a uniform distribution or a Gaussian distribution is used. In this study, we assume that the prior distribution is Gaussian N (−θ∗, I). θ∗ is a parameter for which the risk function R(θ) is minimized. The mean value −θ∗ of the prior distribution corresponds to the shift in θ in Eq.5, which actually eliminates prior knowledge of the posterior distribution. He et al. (2019) uses N (0, I) for the prior distribution. However, it is unnatural because the prior distribution has prior knowledge of the posterior distribution. Conversely, our setting
becomes more natural. The posterior distribution Q and prior distribution P are represented by p(θ) and q(θ) with respect to the parameter θ, given by the following equation:
p(θ) = 1√
(2π)d det(I) exp
( −1
2 (θ + θ∗)>I(θ + θ∗)
) (29)
q(θ) = 1√
(2π)d det(Σ) exp
( −1
2 θ>Σ−1θ
) . (30)
Here, d represents the dimension of the parameter θ, and Eq.30 computed the normalization term M in Eq.23. Thus, we obtain
log
( q(θ)
p(θ) ) = log (√ (2π)d det(I)
(2π)d det(Σ) exp
( −1
2 (θ + θ∗)>I(θ + θ∗)− 1 2 θ>Σ−1θ
))
= 1
2 log
( 1
det(Σ)
) + 1
2
( (θ + θ∗)>I(θ + θ∗)− θ>Σ−1θ ) . (31)
Substituting Eq.31 into Eq.22, we can compute the KL divergence between Q and P as follows. Here, we assume Θ = Rd.
D(Q‖P ) = Eθ∼Q ( Q(θ)
P (θ) ) =
∫ θ∈Θ log ( q(θ) p(θ) ) q(θ)dθ
= ∫ θ∈Θ [ 1 2 log ( 1 det(Σ) ) + 1 2 ( (θ + θ∗)>I(θ + θ∗)− θ>Σ−1θ )] q(θ)dθ
= 1
2 log
( 1
det(Σ)
) + 1
2 ∫ θ∈Θ θ>Iθq(θ)dθ + 1 2 ∫ θ∈Θ θ∗>Iθ∗q(θ)dθ + ∫ θ∈Θ θ∗>Iθq(θ)dθ
− 1 2 ∫ θ∈Θ θ>Σ−1θq(θ)dθ
= 1
2 log
( 1
det(Σ)
) + 1
2 Eθ∼N (0,Σ)
[ θ>Iθ ] + 1
2 Eθ∼N (0,Σ)
[ θ∗>Iθ∗ ] + Eθ∼N (0,Σ) [ θ∗>Iθ ] − 1
2 Eθ∼N (0,Σ)
[ θ>Σ−1θ ] = 1
2 log
( 1
det(Σ)
) + 1
2 tr(Σ) +
1 2 ‖θ∗‖2 − 1 2 tr(I)
≤ 1 2 log
( 1
det(Σ)
) + 1
2 tr(Σ) +
1 2 R2 − 1 2 tr(I). (32)
The last transformation is obtained because ‖θ∗‖2 ≤ R2 always holds when the search range of the parameter θ is in a hypersphere of radius R. From Eq.24, we obtain
AΣ + ΣA = η
|S| C.
Therefore, we obtain AΣA−1 + Σ = η
|S| CA−1. (33)
By computing the trace on both sides, the equation becomes tr ( AΣA−1 + Σ ) = tr ( η
|S| CA−1
) . (34)
The left-hand side (LHS) can be transformed as follows: LHS = tr ( AΣA−1 + Σ ) = tr ( AΣA−1 ) + tr (Σ)
= tr ( ΣA−1A ) + tr (Σ)
= tr (Σ) + tr (Σ)
= 2tr (Σ) . (35)
Thus, we can derive the following:
tr (Σ) = 1
2 tr
( η
|S| CA−1
) = 1
2
η |S| tr ( CA−1 ) . (36)
We can also compute the following: tr(I) = d, (37)
because d is the dimension of parameters θ and I ∈ Rd×d. Therefore, by substituting Eq.36 and Eq.37 into Eq.32, we obtain
D(Q‖P ) ≤ 1 4 η |S| tr(CA−1)− 1 2 log (det(Σ)) + 1 2 R2 − 1 2 d. (38)
Eq.38 is an upper bound on the distance between the distribution Q learned by SGD and the a priori distribution P . Therefore, by substituting Eq.38 into Eq.21, we can obtain a PAC-Bayesian generalization bound for SGD.
The proof is completed.
Lemma 4 (extension of He et al. (2019), Appendix B.2). Assuming that A and Σ commute, the KL divergence between the distribution Q of SGD and the prior distribution P satisfies the following inequality:
D(Q‖P ) ≤ η 4|S| tr(CA−1) + 1 2 d log ( 2|S| η ) − 1 2 log ( det(CA−1) ) + 1 2 R2 − 1 2 d. (39)
Lemma 4 provides the distance between the distribution Q obtained by SGD and the prior distribution of the parameters. The assumption that A and Σ commute implies that matrices A and Σ are simultaneously diagonalizable, and similar assumptions have been used in the work of Jastrzębski et al. (2017) and Liu et al. (2021). Although similar results can be derived without assuming that A and Σ commute, we used the above assumption for simplicity. Based on this assumption, we can compute a generalization bound for a constant case.
Proof of Lemma 4. Because A and Σ commute by assumption, the following equation can be obtained by transforming Eq.24:
AΣ + ΣA = η
|S| C
2ΣA = η
|S| C
Σ = η
2|S| CA−1. (40)
Thus,
det(Σ) = det
( η
2|S| CA−1
) = ( η
2|S|
)d det(CA−1), (41)
and we can compute
log (det(Σ)) = log
[( η
2|S|
)d det(CA−1) ]
= −d log (
2|S| η
) + log ( det(CA−1) ) . (42)
Therefore, Eq.39 can be obtained by substituting Eq.42 into Eq.38.
The proof is completed.
Thus, we can easily prove Theorem 1.
Proof of Theorem 1. It is obtained by substituting Eq.39 into Eq.21.
The proof is completed.

heading: B.2 ADDITIONAL TWO FACTORS IN THE PROPOSED GENERALIZATION BOUND
text: In addition to the three factors discussed in He et al. (2019)’s Theorem 2, the local geometry around minima, gradient fluctuations, and hyperparameters, the generalization bound derived in this study incorporates a factor for the range of parameter search and a factor for network initialization.
Parameter search space The norm of θ after learning appears as the third term R2 in Eq.32. Hence, the parameter search range and weight decay are naturally incorporated into the generalization bound.
Network initialization. Moreover, this generalization bound naturally incorporates the network initialization factor. Assuming that the prior distribution is N (−θ∗, λI), 0 < λ, we obtain the following:
p′(θ) = 1√
(2πλ)d det(I) exp
( − 1
2λ (θ + θ∗)>I(θ + θ∗)
) , 0 < λ. (43)
Eq.32 is as follows.
D(Q‖P ′) ≤ 1 2 log
( λd
det(Σ)
) + 1
2λ tr(Σ) +
1 2λ R2 − 1 2 tr(I). (44)
By differentiating the right-hand side (RHS) by λ, we can observe that RHS is minimized by RHS = d 2 log
( tr(Σ)+R2
d
) − 12 log (det (Σ)), when λ = tr(Σ)+R2 d . Assuming that d is sufficiently large and
tr(Σ) + R2 < d is valid (this property is also called overparameterization (Allen-Zhu et al., 2019; Brutzkus et al., 2018; Du et al., 2019)), we can theoretically guarantee that the generalization bound will be tighter when we change the covariance matrix I of the prior distribution to a diagonal matrix with even smaller eigenvalues, such as He’s initialization (He et al., 2015) or Xavier’s initialization (Glorot & Bengio, 2010). Therefore, the generalization bound proposed in this study, which is an extended version of He et al. (2019), also naturally incorporates factors related to the initialization of the network.

heading: B.3 GENERALIZATION BOUND FOR THE TIME-VARYING CASE
text: To prove Theorem 2, we only need to describe how the time-varying batch size and learning rate change Eq.24 in Lemma 2. To do so, we solve following Lemma 5: Lemma 5 (extension of Gardiner (2004), pp.115-116). Under second-order differentiable assumption (Eq.5), if t is sufficiently large and hyperparameter η|S| is time-varying, the Ornstein–Uhlenbeck process (Eq.4)’s distribution
q(θ) = M exp { −1
2 θ>Σ−1θ
} (45)
has the following property: AΣ + ΣA = n(ξ)C, ξ ∈ [0, t]. (46)
Here,
n(t) = η(t)
|S(t)| . (47)
Proof of Lemma 5. From the analytical solution of the Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930) in the time-varying case (Gardiner, 2004), the parameter θ can be expressed as follows:
θ(t) = e−Atθ(0) + ∫ t 0 √ n(t′)e−A(t−t ′)BdW (t′), (48)
where W (t′) is a white noise and follows N (0, I). From Eq.45, we obtain Σ = Eθ∼Q [ θθ> ] . (49)
Thus, we obtain the following equation:
AΣ + ΣA = Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ ∫ t 0 n(t′)Ae−A(t−t ′)Ce−A(t−t ′)dt′
+ ∫ t 0 n(t′)e−A(t−t ′)Ce−A(t−t ′)Adt′ (50)
whereA is the Hessian matrix and symmetric. We use the mean-value theorem of integrals to obtain the following equation: For some ξ1, ξ2 ∈ [0, t],∫ t
0
n(t′)Ae−A(t−t ′)Ce−A(t−t ′)dt′ = n(ξ1) ∫ t 0 Ae−A(t−t ′)Ce−A(t−t
′)dt′ (51)∫ t 0 n(t′)e−A(t−t ′)Ce−A(t−t ′)Adt′ = n(ξ2) ∫ t 0 e−A(t−t ′)Ce−A(t−t ′)Adt′ (52)
Here, we assume that the matrix D is
D = ∫ t 0 e−A(t−t ′)Ce−A(t−t ′)dt′ (53)
and that A and D commute. Then, because AD = DA,
ξ = ξ1 = ξ2 (54)
holds for Eq.51 and Eq.52. Thus, by substituting Eq.54 into Eqs.51 and 52, and then substituting the result into Eq.50, we derive
AΣ + ΣA = Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ n(ξ) ∫ t 0 Ae−A(t−t ′)Ce−A(t−t ′)dt′
+ n(ξ) ∫ t 0 e−A(t−t ′)Ce−A(t−t ′)Adt′
= Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA
+ n(ξ) ∫ t 0 d dt′ e−A(t−t ′)Ce−A(t−t ′)dt′
= Ae−At〈θ(0), θ>(0)〉e−At + e−At〈θ(0), θ>(0)〉e−AtA + n(ξ)C − n(ξ)e−AtCe−At
' n(ξ)C. (55)
The last transformation is guided by the assumption that A is a positive definite matrix and all its eigenvalues are positive and by the approximation of the term multiplied by e−At twice with zero, because t is sufficiently large.
The proof is completed.
Proof of Theorem 2. Using Eq.55, rather than Eq.27 in the proof of Theorem 1, leads to Eq.10. Furthermore, as exactly
min(n(t)) ≤ n(ξ) ≤ max(n(t)), (56)
by applying a certain λ ∈ [0, 1], we can express
n(ξ) = (1− λ) min(n(t)) + λmax(n(t)).
Therefore, Eq.11 was derived simultaneously.
The proof is completed.

heading: B.4 FUNCTIONAL MODEL OF THE TRAINING ERROR

heading: B.4.1 CONSTANT CASE
text: We derive a training error model for the constant case. We approximate the training error R̂(θ) by substituting Eq.7 into Eq.6 and then approximate R̂(Q) by computing the expected value in the parameter distribution Q. In the constant case, we can compute
R̂(θ) ' 1 2 (θ − θb)>Â(θ − θb) = 1 2
( e−Atθ(0)− θb + √ η
|S| ∫ t 0 e−A(t−t ′)BdW (t′) )> Â ( e−Atθ(0)− θb + √ η
|S| ∫ t 0 e−A(t−t ′)BdW (t′) ) = 1
2 (θ(0)>e−At − θ>b )Â(e−Atθ(0)− θb)
+ 1
2 (∫ t 0 e−A(t−t ′)BdW (t′) )> Â (∫ t 0 e−A(t−t ′)BdW (t′) ) η |S|
+ ((∫ t 0 e−A(t−t ′)BdW (t′) )> Â(e−Atθ(0)− θb) )√ η |S| . (57)
Therefore, by taking the expected value with the distribution Q, we obtain R̂(Q) = Eθ∼Q [ R̂(θ) ] ' Eθ∼Q [ 1
2 (θ − θb)>Â(θ − θb) ] = E [ 1
2 (θ(0)>e−At − θ>b )Â(e−Atθ(0)− θb) ] + E [ 1
2 (∫ t 0 e−A(t−t ′)BdW (t′) )> Â (∫ t 0 e−A(t−t ′)BdW (t′) )] η |S|
+ E [((∫ t 0 e−A(t−t ′)BdW (t′) )> Â(e−Atθ(0)− θb) )]√ η |S|
= a0 + a1 η
|S| + a2
√ η
|S| . (58)
This allows us to model R̂(Q), as shown in Eq.13.

heading: B.4.2 TIME-VARYING CASE
text: To compute the training error model for the time-varying case, we only need to consider how the time-varying batch size and learning rate change Eq.6. Thus, we use the mean value theorem of integrals in Eq.9 to obtain the following equation:
θ(t) = θ(0)e−At + ∫ t 0 √ n(t′)e−A(t−t ′)BdW (t′)
= θ(0)e−At + √ n(ξ) ∫ t 0 e−A(t−t ′)BdW (t′). (59)
Here, ξ ∈ [0, t] and as exactly min( √ n(t)) ≤ √ n(t) ≤ max( √ n(t)), (60)
for some λ ∈ [0, 1], we obtain√ n(ξ) = (1− λ) min( √ n(t)) + λmax( √ n(t)). (61)
Thus, we can compute the first term in Eq.16, which is the training error model for the time-varying case, by substituting Eq.61 into Eq.59 and then performing the same transformation as in Appendix B.4.1.

heading: B.5 FUNCTIONAL MODEL OF THE GENERALIZATION BOUND

heading: B.5.1 CONSTANT CASE
text: In this section, under appropriate assumptions, we show that Eq.8 monotonically decreases with respect to η|S| as a basis for deriving Eq.14 based on Eq.8. We recall He et al. (2019)’s proof to make our paper complete. We first make the following assumptions about dimension d of the parameter space Θ.
Assumption 1 (See He et al. (2019), Assumption 2). The network size is large enough:
d > tr(CA−1)η
2|S| , (62)
where d is the number of parameters, C is the magnitude of the individual gradient noise, A is the Hessian matrix around the global minima, η is the learning rate, and |S| is the batch size.
This assumption can be justified by the fact that network sizes of neural networks are often extremely large. This property is also called overparameterization (Allen-Zhu et al., 2019; Brutzkus et al., 2018; Du et al., 2019). We can obtain the following corollary by combining Eq.8 and Assumption 1.
Corollary 1 (See He et al. (2019), Appendix B.3). When all conditions of Theorem 1 and Assumption 1 hold, the generalization bound of the network is negatively correlated with the ratio of the learning rate to the batch size.
Proof of Corollary 1. We first define
I = η
2|S| tr(CA−1)+d log ( 2|S| η ) − log (det (CA−1))+R2−d+2 log (1 δ )+2 logN+4. (63)
Then, Eq.8 becomes
R(Q) ≤ R̂(Q) + √ I
4N − 2 . (64)
We calculate the derivative of I with respect to the ratio η|S| to check whether the generalization bound has a negative correlation with the ratio. For brevity, we set k = η|S| :
∂I ∂k = ∂ ∂k
[ k
2 tr(CA−1)− d log(2k)− log (det (CA−1)) +R2 − d+ 2 log (1 δ ) + 2 logN + 4 ] = tr(CA−1)
2 − d k . (65)
Therefore, when Assumption 1 holds, we have
d > tr(CA−1)η
2|S| = k 2 tr(CA−1). (66)
Thus, ∂I
∂k < 0. (67)
Then, I and the generalization bound have a negative correlation with the ratio of the learning rate to the batch size.
The proof is competed.

heading: B.5.2 TIME-VARYING CASE
text: In this section, the basis for deriving the second term in Eq.16 based on Eq.10, we will show that Eq.10 monotonically decreases with respect to n2 in Eq.18 under appropriate assumptions. First, we make the following assumptions about the dimension d of the parameter space Θ. Assumption 2 (cf. He et al. (2019), Assumption 2). The network size is large enough:
d > tr(CA−1)
2 n2, (68)
where d is the number of parameters, C is the magnitude of the individual gradient noise, A is the Hessian matrix around the global minima, and n2 is the intermediate value of η(t) |S(t)| in Eq.18, η(t) is the learning rate, and |S(t)| is the batch size.
This assumption can also be justified as Assumption 1 when the network sizes of neural networks are often extremely large. We can obtain the following corollary by combining Eq.10 and Assumption 2. Corollary 2 (extension of He et al. (2019), Appendix B.3). When all conditions of Theorem 2 and Assumption 2 hold, the generalization bound of the network is negatively correlated with the ratio of the learning rate to the batch size.
Proof of Corollary 2. Along the proof of Corollary 1, we use Eq.10, k = n2 and Assumption 2, rather than Eq.8, k = η|S| and Assumption 1. Then, we can obtain
∂I ∂k < 0. (69)
The proof is competed.
C IMPLEMENTATION DETAIL

heading: C.1 DATASETS AND MODELS

heading: C.1.1 DATASETS
text: We tested the proposed model on two popular image datasets. CIFAR10 and CIFAR100 (Krizhevsky et al., 2009): 60 K natural RGB images of 10 classes for CIFAR10 and 100 classes for CIFAR100 with a train/test split of 50K/10 K. For both datasets, we use PyTorch version1.

heading: C.1.2 MODELS
text: We tested the proposed model using two popular model networks. For CIFAR10, we used VGG16 (Simonyan & Zisserman, 2015), and for CIFAR100, we use WRN28-10 (Zagoruyko & Komodakis, 2016). For both model networks, we build on the code from the implementation of Kim (2018).

heading: C.1.3 TRAINING
text: We modified the implementation of Kim (2018). Both datasets are normalized by mean and standard deviation. In the main experiments, training was performed via SGD with a momentum of 0.9, and weight decay of 5e-4 for 200 epochs. We began training with a learning rate of ηinit, run for 200 epochs, and reduced by a multiplicative factor of γ after 60, 120, and 160 epochs to make ηfinal = γ 3ηinit.

heading: C.2 ERROR ESTIMATION EXPERIMENT

heading: C.2.1 EXPERIMENTAL DETAILS
text: As described in Section 6.2, we fit the parameters of φ to minimize δ(S, η;φ), using least squares regression. In doing so, we used scipy.optimize.curve_fit2. The optimal parameter φ∗ is given by
φ∗ = arg min φ ∑ S,η |δ(S, η;φ)|2 . (70)
To measure the performance of the proposed model, the mean µ and standard deviation σ of the relative error δ(S, η;φ∗) are computed based on the fitted parameter φ∗. To evaluate the stability of the model, we fit a parameter φ based on a randomly sampled generalization error and compute µ and σ from the relative error δ(S, η;φ). The values of µ and σ obtained after 100 repetitions are shaded within one standard deviation, as shown in Figure.2.

heading: C.2.2 FOUND PHI VALUES

heading: C.2.3 COMPARISON MODELS FOR THE TIME-VARYING CASE
text: For the time-varying case in this study, we used a step decay for the learning rate. Thus, rather than using the mean value theorem of integrals to make approximations, as in Eqs.55 and59, we can obtain a more accurate model by treating the integral value as a constant for each interval as follows:∫ t
0
n(t′)F (t′)dt′ = n(0) ∫ t0 0 F (t′)dt′ + n(t0) ∫ t t0 F (t′)dt′ = n(0)I0 + n(t0)I1. (71)
We call this the analytical solution model. We call the model using only the final learning rate as the stationary model and fit these two models to the generalization error. The comparison results can be found in Table.2.
Here, the analytical solution model uses step decay and can divide the integration interval into four, and then |φ| = 13. Meanwhile, the stationary model has the same number of parameters as the constant case, |φ| = 5. As shown in these results, the analytical solution model obtains the smallest relative error, and the proposed model has almost the same performance with a smaller number of parameters. Assuming that the relative error δ follows a normal distribution N (µ, σ) for the mean µ and standard deviation σ of the relative error δ in the fitted model, we can observe that Akaike’s information criterion (AIC) (Akaike, 1998) is minimized by the proposed model. Thus, the integral in the proposed model can be approximated sufficiently using the mean value theorem.

heading: C.3 STABILITY EXPERIMENT

heading: C.3.1 EXPERIMENTAL DETAILS
text: We use the same experimental setup as in Appendix C.2.1. However, there is instability in scipy.optimize.curve_fitwith respect to the initial value of the search parameterφ. Thus, we sampled an initial value of φ from a uniform distribution [0, 1] and fitted it several times (10 and 5 times for the constant and time-varying cases, respectively). Then, we used the parameter that best fit the sampled data for the stability experiment.

heading: C.4 HYPERPARAMETER OPTIMIZATION EXPERIMENT

heading: C.4.1 BAYESIAN OPTIMIZATION MODEL FOR THE CONSTANT CASE
text: In this section, we describe the implementation details of the original kernel function proposed in Eq.20. To construct the proposed kernel function, we implemented the following kernel function using sklearn.gaussian_process.kernels3:
from sklearn.gaussian_process.kernels import DotProduct, ConstantKernel
sigma_0_bounds = (1e-10, 1e10) equation:proposed_kernel_constant_case = (DotProduct(sigma_0_bounds=sigma_0_bounds) ** 0.5
+ DotProduct(sigma_0_bounds=sigma_0_bounds) + DotProduct(sigma_0_bounds=sigma_0_bounds) ** -0.5 + ConstantKernel())
sigma_0_bounds decides the range of the σ, which controls the inhomogenity of the kernel.

heading: C.4.2 ALGORITHM FOR THE TIME-VARYING CASE
text: The proposed model Eq.16 for the time-varying case cannot be represented by a linear combination of y = w>ψ(x), as the product term of the parameters of c, λ appears. Thus, we employed a hyperparameter search by sequential model-based global optimization (Bergstra et al., 2011) combining the Levenberg–Marquardt algorithm (LM) (Moré, 1978) and Thompson Sampling (Thompson, 1933; Russo et al., 2017), directly using the proposed model. The LM is used to solve non-linear least squares problems. It interpolates between the Gauss—Newton algorithm (GNA) and the GD method. It is more robust than the GNA, which implies that, in many cases, it finds a solution even if it starts very far from the final minimum. It can also be regarded as the GNA using a trust region approach. The distribution of the parameter φ is obtained by minimizing the relative error of δ(S, η;φ) = ̂(S,η;φ)− (S,η) (S,η) with respect to the parameter φ by the LM. The algorithm is as follows.
Algorithm 1 Hyperparameter search algorithm for the time-varying case H ← ∅ sample φ0 from U(0, 1) for t← 1 to T do
η |S| ∗ max , η|S| ∗ min ← arg min η |S|max , η|S|min δ(S, η;φt−1) Evaluate f( η|S| ∗ max , η|S| ∗ min
) . Expensive step H ← H∪ (( η|S| ∗ max , η|S| ∗ min ), f( η|S| ∗ max , η|S| ∗ min )) µ,Σ← LM(δ,H) sample φt from trunc-norm(µ, √ diag(Σ), blower, bupper)
end for returnH
Here, H represents the history of observations, and f( η|S| ∗ max , η|S| ∗ min
) represents the evaluation of the network for a hyperparameter ( η|S| ∗ max , η|S| ∗ min
). To keep the parameter φ in this range, Thompson sampling is performed for each parameter from the following truncated normal distribution using the lower bound blower and upper bound bupper.
f(x;µ, σ, blower, bupper) = 1
σ
ϕ(x−µσ )
Φ( bupper−µ σ )− Φ( blower−µ σ ) , (72)
where ϕ(x) is the probability density function of the standard normal distribution as follows:
ϕ(x) = 1√ 2π exp
( −1 2 x2 ) , (73)
3https://github.com/scikit-learn/scikit-learn
and Φ(x) is its cumulative distribution function:
Φ(x) = 1
2
( 1 + erf ( x√ 2 )) . (74)

heading: C.4.3 EXPERIMENTAL DETAILS
text: In this section, we present the implementation details of the hyperparameter optimization experiments for comparison. In both constant and time-varying cases, η|S| was allowed to explore in the range of [10−7, 10−2] for CIFAR10 and in the range of [10−6, 10−2] for CIFAR100. These were determined based on the experimental results for each dataset (Figures.1a and 1b). We computed the maximum test accuracy obtained over 30 rounds of hyperparameter search. By performing this for a total of 100 times, we computed a 95% confidence interval of the maximum accuracy. The proposed method uses a uniform distribution search, whereas the comparison method uses a uniform distribution search and log-uniform search. For comparison, we used Bayesian optimization4 (Nogueira, 2014) using the Matern kernel and RBF kernel, Hyperopt5 (Bergstra et al., 2015), and Optuna6 (Akiba et al., 2019).
4https://github.com/fmfn/BayesianOptimization 5https://github.com/hyperopt/hyperopt 6https://github.com/optuna/optuna

Review summaries:
The paper considers the functional relationship between (time-varying) learning rate & batch size and the generalization error. The contributions are summarized below: * Assuming quadratic risk function and the SDE limit, the authors derive a novel theoretical generalization bound when the (learning rate / batch size) ratio monotonically decreases. * Inspired by the generalization bound, the authors propose a heuristic functional form of the generalization error w.r.t. the (learning rate / batch size) ratio. * The authors conduct experiments on CIFAR10 and CIFAR100 datasets, and show that the heuristic generalization error model can fit the actual generalization curve relatively well. * The authors propose a novel kernel function for Bayesian optimization in hyperparameter optimization. The kernel function is inspired by the heuristic generalization error model. Empirical results show that the proposed kernel can match or outperform existing Bayes optimization baselines. 

The authors study the dependence of generalization error of trained neural networks on the batch size and the learning rate used in SGD. The primary motivation is to extend the previous works to the case where batch size and learning rate change over the trajectory. The authors present a functional form using PAC-based generalization bounds to model the desired dependence. They conduct extensive experimentation to show that the functional form approximates the generalization error well. In addition, the authors show that a hyperparameter search based on the proposed model outperforms existing hyperparameter optimization libraries like Hyperopt and Optuna.

This paper approximates the learning dynamics by an Ornstein-Uhlenbeck process and uses some PAC-Bayes results to obtain a functional form for the generalisation behaviour as a function of the ration of the learning rate to the batch size.  It obtains results for both constant and time-varying parameters.  It then uses these to develop a kernel function for hyper-parameter optimisation using Gaussian Processes modelling.

The authors study the learning rate over batch size dependence of generalization bounds and try to use them for hyperparameter tuning.  In section 4, the authors study generalization bounds for convex optimization with the motivation that this helps understand neural networks. In section 5, the authors use the bounds to propose a functional approximation for the generalization based as a function of the n=eta/|S| scaling of the fluctuations around the minimum, including for the possibility that n is time dependent. The generalization error is roughly a sum of integer powers of sqrt(n). In section 6, the authors the authors fit a set of points to the generalization formula they found. In section 7, the authors use this expression to look for optimal hyperparameters and compare them with other hyperparameter tuning methods.