Neural networks [31], also known as artificial neural networks, are one of the most well-known algorithms used to recognize patterns and solve common problems in statistics, data science and machine learning. A neural network is comprised of an input layer, one or more hidden layers and an output layer, with several units called neurons in each layer. Neurons in the input layer bring the initial data into the network for further processing, while the result from this network is produced by neurons in the output layer. Each neuron in hidden layers connects to another via receiving input produced by neurons in the preceding layer and producing the output to neurons in the next layer. Neural network models, in general, are very flexible, which, on the other hand, brings challenges in terms of constructing the best network structures. See [29] for more details on neural networks.

In this day and age, problems, such as image classification, object detection, or natural language processing task, become increasingly complex. This calls for deep neural networks. A deep neural network refers to an artificial neural network with more than one hidden layer [2, 37]. As an example, consider a deep neural network with three hidden layers, which is shown in Figure 1. The number of neurons in the hidden layers is one of the common hyperparameters of interest. The number of neurons in the input layer is equal to the number of features in the data, while the number of neurons in the output layer depends on the type of problems we solve. For example, very often, there is one neuron in the output layer for a regression problem. If a multi-class classification problem is considered, one common choice is to use one neuron per a class.

Hyperparameters for neural networks are those variables whose values determine the network structure and/or the way a network is learned. The number of neurons in the hidden layers of Figure 1 is an example of a hyperparameter. Typically, a neural network requires us to determine values for a set of different hyperparameters including, but not limited to, learning rate, batch size, number of layers, dropout rate, number of epochs, number of training iterations, normalization, pooling, momentum and weight initialization. Those hyperparameter settings drastically affect the success and accuracy of the network. As a consequence, finding an appropriate hyperparameter combination is increasingly important.

Design of experiments is a systematic and efficient method that allows scientists and engineers to explore the relationship between multiple input variables and output variables. Traditional physical experiments favor various factorial designs for studying systems or processes [45]. But they are no longer proper when the systems or processes become complex, which calls for computer experiments and space-filling designs, a family of designs appropriate for computer experiments [10, 36].

Factorial Designs At the early stage of an investigation, due to the lack of enough prior knowledge, the experimenter conducts a screening experiment involving as many variables as possible and the primary goal is to identify some important factors from them using a first-order model. Two-level factorial designs [7, 33], thanks to their simple structure and nice statistical properties, are commonly used for this purpose. The experimenter often proceeds to the next stage by capturing the curvature in the response surface, which requires the consideration of three-level factorial designs or composite designs such as orthogonal array composite designs (OACDs) [48].

Space-Filling Designs A space-filling design refers to a design that scatters its design points uniformly over the whole design region. It is tremendously complicated to find a space-filling design that provides good coverage of the entire input space, especially in the high-dimensional input space. Instead, it is natural and reasonable to consider a design that enjoys the space-filling property in low-dimensional projections. This idea dates back to Latin hypercube designs (LHDs) proposed by [32]. Such designs are orthogonal arrays of strength one and enjoy the maximum space-filling property in all univariate projections – there is exactly one observation in any interval formed through dividing the range of an input variable into the same number of equally spaced intervals as the run size. The space-filling property of an LHD may be further evaluated by other optimality criteria such as orthogonality [5], a distance criterion [20] and a discrepancy criterion [11]. There has been extensive work along this line – see, for example, [12, 21, 27].

Strong orthogonal arrays (SOAs) introduced and studied by [17] are the other class of space-filling designs in the literature with a focus on low-dimensional projections. An SOA of strength t achieves the space-filling properties in all g ≤ t dimensions. Such an array does well as a comparable orthogonal array of strength t in all t dimensions but is more space-filling in all 1 ≤ g ≤ t − 1 dimensions than the latter. It can also be made to be an LHD through level expansion to achieve the maximum space-filling property in all one-dimensional projections. As a result, this SOA-based LHD enjoys the benefits of better space-filling properties than the comparable ordinary LHD in all 2 ≤ g ≤ t dimensions. More developments on SOAs can be found in [15, 18, 28, 38, 42].

The Kriging model originated in the areas of geosciences [23] and is now popular in computer experiments for building a surrogate model for complicated computer models. In a Kriging model, the response is referred to as a realization of a Gaussian process, and the predicted response at a target point can be represented by a weighted average of the responses at observed points. For more details on Kriging models, we refer to [8, 14, 22, 35, 44].

As most learning algorithms are stochastic, instead of using the Universal Kriging model for the deterministic case in computer experiments, we consider the Universal Kriging model with a noise term (2.1) y ( x ) = ∑ i = 1 m β i f i ( x ) + Z ( x ) + ϵ , where ∑ i = 1 m β i f i ( x ) is the mean function with f i being the ith basis function and β i being the ith coefficient, Z ( x ) is a second-order stationary random process with constant mean 0 and the covariance matrix given by Cov [ Z ( u ) , Z ( v ) ] = σ 2 ∏ i = 1 k R u i − v i , where u = ( u 1 , … , u k ) and v = ( v 1 , … , v k ) are two points (runs), R ( · ) is a spatial correlation function and σ 2 is the variance of the random process, and ϵ ∼ N ( 0 , τ 2 ) is independent of Z ( x ). In this study, we adopt the Ordinary Kriging model with a noise term, where the mean function in equation (2.1) is a constant. We choose the Matern correlation function with parameter ν = 5 / 2 in the study based on the conclusion drawn by [47] that the Matern correlation function provides a good balance between differentiability and smoothness. The explicit form is given by R ( h i ; θ i ) = 1 + 5 h i θ i + 5 h i 2 3 θ i 2 exp − 5 h i θ i , where h i = | u i − v i | and θ i is the unknown hyperparameter which can be estimated using the maximum likelihood estimation method.

MNIST Dataset The MNIST dataset (Modified National Institute of Standards and Technology database) is a subset of a larger dataset available at NIST, the National Institute of Standards and Technology, and has served as a canonical training dataset for many learning techniques and pattern recognition methods. It is a dataset of handwritten digits, first developed and released by [25], containing 70000 grayscale images of the 10 digits that are 28 × 28 pixels in width and height. Some examples are demonstrated in Figure 2. These images come with labels and thus are used for supervised learning tasks. MNIST is user-friendly, as images in this set have been split into a training set of 60,000 images, and a test set of 10,000 images. In this paper, we further consider a subset of the training set as a validation set. Each time, we randomly take out 15% of the full training set to serve as a validation set.

Hyperparameters We consider five common hyperparameters: learning rate, number of epochs, batch size, dropout rate and number of neurons in a hidden layer. These five hyperparameters are numerical factors: learning rate and dropout rate are continuous factors and number of epochs, batch size and number of neurons are discrete factors. The domain of each hyperparameter is given in Table 1.

Designs Various factorial designs and space-filling designs of 5 columns corresponding to five hyperparameters are examined. More specifically, we use three factorials: the 2 5 and 3 5 full factorials and an OACD of 34 runs combining the 2 V 5 − 1 factorial with a three-level 18-run orthogonal array. Two full factorials are generated using the R package DoE.base [16]. The 34-run OACD can be constructed using Table 2 in [48] and is given in the supplementary material. In addition to three factorials, this study includes a 243-run random LHD generated using the R package lhs [6] and several 32-run space-filling designs – an SOA of strength three, its corresponding LHD and four other types of LHDs: a random LHD, a maximin LHD, a maximum projection LHD and a uniform LHD, which are generated using R packages lhs, SLHD [1], MaxPro [21], and UniDOE [49], respectively. There are 32 levels in each column of these LHDs while the SOA of strength three has 8 levels and is listed in the supplementary material. An LHD based on this SOA can be obtained by expanding 8 levels to 32 levels following [17]. The two largest designs, 3 5 factorial and 243-run random LHD, serve as the benchmarks in comparison with designs of small run sizes.

The design matrices of the above designs in terms of natural units are provided in the supplementary material for reference, where the lowest and highest levels are given in Table 1 for each factor. We linearly interpolate other levels within the range. For discrete variables, we further round the levels to the nearest integers.

Data Collection All neural network implementations are built in TensorFlow: https://www.tensorflow.org. There is a stochastic component to the achieved accuracy, as images are shuffled for each model trained. The network based on the training dataset is further assessed on the test dataset of 10,000 images using test accuracy. The test accuracy is appended in the tables of design matrices with natural units in the supplementary material. We also record the cross-entropy loss of each hyperparameter combination in the column right before the test accuracy column in these tables for those interested parties.

We evaluate designs by considering their ability to collect informative data for building a statistical model that specifies the relationship between hyperparameters and test accuracy. In our comparison, the Kriging model is an appropriate consideration, as it can compensate for the effects of data clustering and give a better estimation of prediction error. More specifically, we use the Ordinary Kriging model with a noise term. As test accuracy has values with ( 0 , 1 ), we consider the Arcsine Transformation for test accuracy when building the model.

Essentially, the model is built using a set of data including hyperparameter combinations selected by a design we consider and their corresponding accuracy, which is then expected to be evaluated on the data comprising as many diverse hyperparameter combinations as possible over the whole hyperparameter region and their accuracy. This, however, is infeasible especially for the continuous hyperparameters. In this paper, we consider two extreme cases – generating test data from the 3 5 factorial and 243-run random LHD. The random LHD enjoys the maximum space-filling property in all one dimensions, while the 3 5 factorial covers the entire 5-dimensional input space in a uniform fashion. The root mean square error (RMSE) values for predicting the mean test accuracy on these two designs are listed in Table 2. Table 2 also includes results on the hybrid design combining the 3 5 factorial with 243-run random LHD.

When tested on the 3 5 factorial, from Table 2, it can be seen that a three-level design, 34-run OACD, outperforms all other small designs. This design performs as well as we expected because this 34-run OACD is a design comprising 34 runs taken from the 3 5 factorial. In other words, a subset of the test data is used to train the Kriging model in this case. Surprisingly, although the SOA is tested on a different type of design, the 3 5 factorial, it dramatically outclasses all other small designs except the OACD. Notably, both the OACD and SOA have smaller RMSE than the 243-run LHD when tested on the 3 5 factorial. Moreover, the model on the 3 5 factorial produces the smallest test RMSE, because the test RMSE is indeed the training RMSE.

Considering the test data generated from the 243-run random LHD, the results in Table 2 clearly show that all factorials are worse than those space-filling designs, and that, when the Kriging model is built on the 243-run random LHD, the test RMSE is the training RMSE and thus is the smallest. Though the maximin LHD produces as the smallest test RMSEs as the SOA in this case, its performance when tested on the 3 5 factorial is not as well as that of the SOA – the SOA performs exceedingly well compared to the other small space-filling designs.

Both the OACD and SOA perform excessively well when tested on the same type of design, the 3 5 factorial and 243-run random LHD, respectively. But their performances tested on the different types of designs considerably differ. The SOA is still very competitive on the 3 5 factorial, as it produces the smallest test RMSE among all comparable designs except the OACD. On the 243-run random LHD, however, the performance of OACD is the worst except for the performance of the 2 5 factorial. The SOA is a clear winner among all small designs we considered. The numerical results in the last column of Table 2 further support the above conclusion. In other words, for the setting we considered, this SOA allows us to efficiently collect the most informative data for building a Kriging model that specifies the relationship between hyperparameters and test accuracy. The SOA of strength three we use is indeed optimal in some sense – Theorem 4 of [40] implies that this SOA is optimal under the uniform projection criterion. Moreover, it enjoys better space-filling properties in all two dimensions than an ordinary SOA of strength three. Notably, the model built on the 2 5 design produces an extremely large test RMSE value when tested on the hybrid design, which implies that this two-level design is inadequate for building a Kriging model to capture the true surface of test accuracy.

Figure 3 and Figure 4 also sustain that the performance of this SOA is superb. Figure 3 displays the density plots of the observed test accuracy values for all designs. Each density plot in the first two rows of Figure 3 corresponds to a type of design that is used to generate the training data. The density curve of observed test accuracy values for the hybrid design combining the 3 5 factorial and 243-run random LHD is present in the last plot of Figure 3. It is bimodal with a higher peak within [ 0.8 , 1 ] and a lower peak of no more than 0.4. Clearly, only the SOA correctly captures this feature while all others fail to do so. We also provide the density curves of observed test accuracy values for the 3 5 factorial and 243-run random LHD, respectively, in the last row for reference. Figure 4 gives the histograms of the test errors which are the differences from predictions of test accuracy on the hybrid design to the observed values. The histograms of those LHDs are skewed to the left, which implies that the models based on these designs tend to overpredict the test accuracy on the 3 5 factorial, while models from the 2 5 factorial and OACD tend to underpredict the test accuracy on the 243-run LHD, as their histograms are skewed to the right. Only the SOA is superior because its histogram depicts an approximately normal distribution that is symmetric around 0.

The significant advantage in building a Kriging model of the SOA brings us to take a closer look at these small designs themselves. For all designs, each column is rescaled to [ − 1 , 1 ], and we then calculate the Euclidean distance from each design point to the center of the design and make a histogram with a density curve for the distances of each design. These plots are given in Figure 5. Histograms of the distances for the 3 5 factorial, 243-run random LHD and their hybrid design provided in the last row of Figure 5 serve as benchmarks for comparison. Only the SOA captures the feature of the hybrid design that the density plot has a light lower tail but a heavy upper tail, while all other designs fail to do so. Remarkably, though a space-filling design is expected to scatter its design points uniformly over the design region, these space-filling designs we considered seem not as space-filling as expected because Figure 5 shows that they do not provide points that are close to the centers or at corners of the regions.

In addition to using the Kriging model to determine the relationship between hyperparameters and test accuracy, we also consider the second-order model. Although the findings from this study are similar to the above, the second-order model is incapable of adequately capturing the test accuracy surface, as the test RMSE values are consistently larger than those obtained using the Kriging model for all cases listed in Table 2.