Resampling Methods in Neural Networks: From Point to Interval Application to Coordinate Transformation

Arising from the field of machine learning, Artificial Neural Networks, also referred to as neural networks, have been used to model the nonlinear relationship between two or more variables and to recognize complex patterns in the data. They provide powerful data analysis tools to handle various types of problems in many fields of engineering sciences, such as:

One of the most important tasks in developing a neural network is related to partitioning the dataset at hand. Most neural network designers randomly and uniformly split the known dataset into three categories: training, validation, and test. The training set is often used to estimate the unknown parameters of net (e.g., weights in a regression feed-forward neural network). Ultimately, the will have been well learned over its training set. However, a low error in the training set does not tell us much about the quality of the network. What would happen if that trained network was applied to another input dataset? This question is highly relevant. To answer that question, one should use the validation set. The validation set helps one to make decisions on how to train the network and how to stop it to avoid overfitting. However, until now, the total training process of the neural network depends not only on the training set, but also the validation set. It means that any choice regarding training and validation splitting affects the performance of the network. In other words, without the test set we cannot conclude anything about that choice and which choice is good enough in terms of network generalization. In machine learning language, therefore, test sets are unseen data that play an important role in the evaluation of the generalized performance of the network.

This idea of splitting data in a random and uniform way and using the test set to measure the predictive accuracy and stability of the selected trained network model is known as Hold-out cross-validation. There is, however, significant limitations of using this method, as pointed out by Ziggah et al. (2019b) as follows: (1) the results are based on uncontrolled chosen split sets; (2) improper splitting of the data set could have an adverse effect on the model performance; and (3) the results are dependent on the existence of a large dataset making it unsuitable to be applied in data-insufficient situations (i.e., sparse dataset). To overcome these problems, K-fold cross-validation has been recommended as an alternative method (Burman 1989; Reitermanová 2010). In this method, the data is separated into K disjoint subsets (K-fold) of approximately equal size, so that each subset will be in a test set just once. For instance, for each K subset of cross-validation the union of K-1 folds is used as a training set for model development while the remaining part is used as the test set for model validation. Therefore, with this method, the division of data is less of a concern than the case of Hold-out cross-validation, because the entire dataset is employed as both training and testing sets.

The development of a neural network always relies on a certain degree of randomness. In general, the randomness arises in two different parts in neural network modeling: (1) choosing the network architecture, initialization, and learning procedure; and (2) choosing where to split the data (training, validation, and testing sets). There is usually no best choice for this randomness (Pan 1998). However, currently, we have a multitude of neural network models and the continuous development of new architectures to meet different purposes. Consequently, the choice about the network itself is less of a concern than data splitting. Many researchers usually draw many realizations for a given split data set to decide among the different networks which one best suits a certain problem. Although the K-fold presents advantages over the Hold-out cross-validation, it is still not clear on how to choose the K disjoint subsets and in some cases, the subsets are not of equal sizes, which does not guarantee a balanced version of cross-validation. Therefore, splitting data methods are still challenging.

Now the most important question is how to maximize the prediction result since it varies according to the randomness both over data splitting and over the network itself including architecture, initialization, and learning procedures. In light of this issue, we provide a framework that can capture the randomness related to the network itself. To accomplish this, we introduce the resampling methods and apply them in the neural networks context applied to surveying engineering.

The resampling methods replace some theoretical assumptions and derivations by repeatedly resampling the original data and making inferences from the re-sampled data (Efron 1982; Wood 2005). Here, we list the reasons for using resampling methods in real time:

One of the resampling methods which has been proposed for geodetic data processing is known as Jackknife. This method was developed by Quenouille (1949). Jackknife was initially proposed to estimate the bias of an estimator by iteratively deleting only one single observation each time from the sample, and then recomputing the statistic or estimator on the rest of the sample data. This Jackknife method usually refers to the classic Delete-1 Jackknife estimator. Based on the availability of powerful computers, several researchers have worked to improve the resampling methods, as shown in Tukey (1958), Miller (1974), Efron (1992), McIntosh (2016), Wang et al. (2022), and Wang and Yu (2020, 2021).

In the context of neural networks, the Delete-1 Jackknife method is a balanced version of cross-validation, in which one observation is used to validate the performance of the neural network trained by the remaining ($\mathrm{n}-1$) observations. In machine learning language, the Delete-1 Jackknife method is commonly referred to as leave-one-out procedure. The other balanced version of cross-validation is known as the Delete-d Jackknife method, or grouped Jackknife. In that case, instead of validating one observation at a time, we validate the subset of size $d$ and train the model with the remaining ($n-d$) observations at a time. The subset size d is selected from all of the observations without replacement. For instance, if we had a set of 10 observations, the application of the Delete-1 Jackknife method would provide 10 cross-validations, with each observation predicted only once. For this same example, by now taking the Delete-d Jackknife method for $\mathrm{d}=2$, we would have $C_{10}^2=45$ independent cross-validation subsets, with each observation predicted/estimated nine times. Different from K-fold cross-validation, Jackknife considers all possible combinations in splitting the dataset without replacement.

In the context of adjustment theory, Jackknife resampling methods have recently been combined with the weighted total least squares to improve both parameter and precision estimation, as can be seen in the works from Wang et al. (2022, 2020) and Wang and Yu (2020, 2021). Although in other branches of science there are applications of resampling methods in neural networks (see e.g., Pan 1998), to the best of our knowledge, there have been no investigations aimed at resampling methods in neural networks with applications for geodesy and surveying engineering.

In this contribution, we adopt the idea of resampling methods in neural networks. First, we extended the Delete-1 Jackknife method by drawing many trials in order to capture the randomness related to the network itself (architecture, initialization, and learning procedure). We called this method Delete-1 Jackknife Trials. Similar to the Delete-1 Jackknife Trials, Hold-out is also extended to its Hold-Out Trials version by randomly drawing many realizations, so that each prediction point is repeated a certain number of times. Further, we have also applied the Delete-d Jackknife. The user knows exactly the number of times each point prediction will be replicated for the cases where Delete-d Jackknife and Delete-1 Jackknife Trials are in play. On the other hand, in the case of the Hold-Out Trials method, the user chooses the percentage of division of the sample for each set (training, validation and test) and the number of realizations (trials). Consequently, the advantages of the Delete-d Jackknife, Delete-1 Jackknife Trials, and Hold-Out Trials methods are that they allow us to describe the probability distribution for the predictor rather than simply providing a point prediction, which is what is traditionally used in neural networks. This gives us an opportunity to answer the following question: How can we rely only on a single neural network?. It is important to say that the resampling-based cross-validation methods used here are slightly different from the model selection of the linear or nonlinear traditional adjustment methods. In our case, we care less about the best selection of the parameters. Instead, our main concern is with respect to the quality and reliability of the neural network predictor.

The paper is organized as follows. First, a short overview is given on the supervised back-propagation neural network. Next, the three resampling methods are presented in the framework of neural networks, namely Delete-d Jackknife, Delete-1 Jackknife Trials, and Hold-Out Trials. Here, we introduce the resampling methods into the regression neural network model. Beyond the point prediction, it is often desirable and also important to know the confidence associated with the neural network predictor, called interval prediction. Indeed, the uncertainty and bias related to the point prediction are also quantified. The results here are restricted to the coordinate transformation problem without loss of generality. Nonetheless, those resampling methods can be applied to any problem that can be formulated as a neural network model. In fact, they can be extended to any predictive models. Finally, we provide a comparison of the pros and cons of the proposed methods.

A neural network may often be considered a nonlinear or non-parametric regression model (Haykin 1999)which defines the relationship between the known $n$ vector of the output variables $\mathit{y}$ and the known $n$ vector of the input variables $\mathit{x}$. The term $\mathbf{ϵ}$ can be interpreted as a random error vector of size $n$ and the function $\mathcal{N}$ depends on some unknown parameter $\mathit{\theta }$ which may have infinite dimensions.

Typically, the unknown parameters $\theta$ are represented by weights and biases in a supervised feed-forward neural network. A feed-forward neural network is an artificial neural network wherein connections between the neurons do not form a cycle. In that case, the network is said to be fully connected in the sense that every node (also called neuron) in each layer of the network is connected to every node in the adjacent forward layer. For example, it can be formulated aswhere ${\mathit{y}}_{\mathit{o}}$ is the output vector with $n$ observations, ${\mathcal{F}}_o$ is the activation function for the neuron $o$ in the output layer (green node), $w_{k,o}$ (green lines) are the unknown weights between the neuron $o$ in the output layer and each neuron $k$ in the hidden layer (blue nodes), ${\mathcal{F}}_k$ is the activation function for each neuron $k$ in the hidden layer, $w_{i,k}$ (blue lines) are the unknown weights relative to the neurons $k$ in the hidden layer and each observation $i$ of the input data in the input-layer (red node), $b_k$ is the unknown bias for each neuron $k$ in hidden layer, and $b_0$ denotes the unknown bias for the output layer $o$. Note that there is a weight $w_{i,k}$ for each observation of the input data. In other words, the larger the input sample, the more weights there are in that neural network. This is a simple case of considering only an input pattern and an output pattern. Obviously, the model in [Eq. (2)] can be extended to the case where we have more than one hidden layer and several inputs and outputs.

The activation functions in the hidden layer (${\mathcal{F}}_k$) are often either logistic-sigmoid (logsig) or hyperbolic tangent-sigmoid ($\tanh$)where $v$ is any input variable. In the case of [Eq. (2)], the $v$ is the result of the linear combination of the observations of the input data, weights, and bias (i.e., the second summation term). The output units usually have a purely linear transfer function of the type ${\mathcal{F}}_o(v)=v$.

Usually, the back-propagation learning algorithm is employed in conjunction with a training algorithm to determine the weights and biases of a neural network. Since the transfer function involved in the neural networks is not linear [see Eq. (3)], a nonlinear least squares method that iteratively minimizes the mean squared errors between the desired outputs $y_{h,n}$ and the outputs response ${\widehat{y}}_{h,n}$ from the network is required (note: the index $h$ represents a given output attribute and $n$ an nth element of that attribute). In each epoch (iteration), therefore, the weights and biases are updated in a supervised manner by the following minimizationwhere $\epsilon$ represents the average squared error over the sample size $N$ for $O$ outputs.

Here, we use the Levenberg–Marquardt algorithm (LM), also known as the damped least squares method. The LM combines two minimization methods: the gradient-descent method (Newton-Raphson) and the Gauss-Newton method (Levenberg 1944). In the gradient-descent method, the mean squared errors ($\epsilon$) is reduced by updating the parameters in the steepest-descent direction. In the Gauss-Newton method, $\epsilon$ is reduced by assuming the least squares function is locally quadratic and finding the minimum of the quadratic form. The LM method acts more like a gradient-descent method when the parameters are far from their optimal value and acts more like the Gauss-Newton method when the parameters are close to their optimal value (Gavin 2020).

As mentioned, the available data are often split into three categories: training set, validation set, and test set. In the back-propagation algorithm, the first step in constructing a neural network is to adjust the weights and biases by using the training set. The performance of the model is often evaluated on the validation set at the end of each iteration (epoch). The aim of the validation set is to improve the generalization of the network learning. When the network starts to over-fit the training set, the performance on the validation set typically begins to decrease. When the validation error increases for a specified number of epochs the training is stopped, and the weights and biases at the minimum of the validation error are returned. This method for avoiding overfitting is called early stopping (Reitermanová 2010). Finally, the test set is used to evaluate the performance of the trained neural network, which theoretically guarantees an unbiased evaluation.

In practice, however, the test set does not guarantee an unbiased evaluation and prediction. If the data is a random sample, an estimate $\widehat{\mathit{\theta }}$ in [Eq. (1)] differs from the true value of $\mathit{\theta }$ because of the sampling error. Consequently, if this network is used to determine a prediction ${\widehat{\mathit{y}}}_{\mathit{i}}=\mathcal{N}({\mathit{x}}_{\mathit{i}},\widehat{\mathit{\theta }})$ for some new input vector ${\mathit{x}}_{\mathit{i}}$, the prediction $\mathcal{N}({\mathit{x}}_{\mathit{i}},\widehat{\mathit{\theta }})$ will be different from the true prediction. Furthermore, as mentioned before, there is a randomness related to both the data splitting and randomness related to the network itself including architecture (number of hidden layers and the structure of the connections), weights and biases, initialization, and learning procedure. In this paper, resampling approaches are used to approximate a sampling distribution of the prediction and estimate parameters of that distribution like its mean value, variance, percentiles, etc.

Here we use resampling methods in neural networks and apply them to coordinate transformation. The coordinate transformation is a classic problem that always arises in various fields of engineering sciences, such as to convert coordinates from one geodetic reference system into another one (Ziggah et al. 2019a), to transform of 3-D point clouds from terrestrial laser scanning for deformation monitoring (Wujanz et al. 2018), to design and optimize geodetic networks (Teunissen 1985) and so on.

As a truly practical application, we present a problem of spatial transformation. In spatial transformation, each point of an input coordinate system is mapped to an output coordinate system. In this study, the input system corresponds to the image’s pixel indices (numCols, numRows), where numCols and numRows are the number of columns and rows in the image, respectively. We can specify locations in terms of their columns ($x$) and rows ($y$) and, therefore, it is a left-handed coordinate system. On the other hand, the output system is related to the grid body-fixed coordinate system, which in this case is a grid designed on white paper (Fig. 2).

The origin of the output system is defined by the intersection of the first vertical line on the left and the first horizontal line at the bottom, which has $(0.0,0.0)$ coordinates. The horizontal and vertical axes are parallel to the grid width and height, respectively. It is therefore a right-handed coordinate system. The output system does not allow coordinates with negative values.

In general, we have the fixed coordinates in the output system $(X_i,Y_i)$ and the observed coordinates in the input system $(x_i,y_i)$, where $i$ denotes the point number. Since the grid dimension is $420\mathrm{mm}\times 297\mathrm{mm}$ and each point is offset by 42 mm in the X-direction and 29.7 mm in the Y-direction, we have a total of 121 points. All of these points were also selected in the image system. The image used here was captured by a cell phone camera at a distance of $\approx 1.5m$ from the target (grid). The dimensions of this image are $817\times 680\mathrm{pixels}$. The image was acquired with a camera oriented obliquely.

The choice related to the architecture, initialization, and learning procedure for the coordinate transformation problem was fixed here. The experiments were designed by considering one hidden layer composed of ten neurons (Fig. 3). The input data are the coordinates from the image system $(x_i,y_i)$ and the output data are the coordinates from the grid body-fixed coordinate system $(X_i,Y_i)$.

The training algorithm used here is the Levenberg–Marquardt back-propagation algorithm (LM). The Nguyen-Widrow algorithm was used to initialize the weights and biases (Nguyen and Widrow 1990). Levenberg’s contribution provided a “damped version” of the minimization problem for [Eq. (4)] (Gavin 2020)where $\mathit{J}$ is the Jacobian matrix containing first derivatives of $\epsilon$ with respect to the weights and biases; $\mathit{I}$ is the identity matrix, providing increment ${\mathit{h}}_{\mathit{lm}}$ to the weights and biases estimated parameters $\widehat{\mathit{\theta }}$. Note that $\mathit{y}$ are the desired outputs and $f(\widehat{\mathit{\theta }})$ are the output’s response from the network.

The non-negative damping factor ($\lambda$) in Eq. (10) is adjusted at each iteration. If the reduction of the mean squared errors $\epsilon$ in equation Eq. (4) is quick, a smaller value can be used, bringing the algorithm closer to the Gauss-Newton algorithm, but if an iteration gives insufficient reductions to the error, $\lambda$ can be increased, which would provide a step closer to the gradient-descent direction. Note that the gradient of $\epsilon$ with respect to $\widehat{\mathit{\theta }}$ is $\mathit{g}=-2{[{\mathit{J}}^T(\mathit{y}-f(\widehat{\mathit{\theta }})]}^T$. This means that for large values of damping factor $\lambda$, the step will be taken approximately in the direction opposite to the gradient. In other words, if any iteration occurs to result in a worse approximation, i.e., $\epsilon (\widehat{\mathit{\theta }}+{\mathit{h}}_{\mathit{lm}})>\epsilon (\widehat{\mathit{\theta }})$, then $\lambda$ is increased. On the other hand, as the solution improves, $\lambda$ is decreased and the solution tends to go to the local minimum and, therefore, LM approaches the Gauss-Newton method (Madsen et al. 2004; Marquardt 1963; Gavin 2020).

The initial value of the damping factor adopted here was $\lambda =0.001$, with a decrease factor of 0.1 and increase factor of 10. The performance of the model was evaluated on the validation set at the end of each iteration (epoch). The aim of the validation set was to improve the generalization of the LM. When the network starts to over-fit the data, the performance on the validation set typically begins to decrease. When the validation error increases for a specified number of epochs the training is stopped, and the weights and biases at the minimum of the validation error are returned. This method for avoiding overfitting is called early stopping. The performance function is also given by the mean squared error $\epsilon$ in Eq. (4). We used the command-line functionality of the Matlab’s Neural Network toolbox (R2019b) for training and validating the neural network. Thus, the learning process was configured to stop when any of these conditions were met:

The inputs and outputs values were normalized according to the z-score as follows:where $z$ represents the sample of the input or output of the neural network, $\overline{x}$ is the mean of the sample, and $S$ is the standard deviation of the sample. In this case, the normalized function gives values between $[-1,1]$.

From Fig. 5, we observe that both the magnitude and the direction of the biases occurred in different ways at each point. In other words, the biases had a random behavior. The biases were less than 1 mm for all resampling methods. Considering all points, the biases averaged 0.31 mm ($\pm 0.17$) for all resampling methods. In terms of the evaluation of the resampling methods on the out-of-sample test set, Jack-1 was also slightly better than the other two methods, as seen in Table 2. From Fig. 6, again, the biases had a random behavior, and in that case, both biases and the RMSE exhibited values less than 1 mm.

Table 1 and Fig. 5 show that, in terms of RMSE, Jack-1 is slightly better [0.50 mm ($\pm 0.21$)] than Jack-d [0.56 mm ($\pm 0.29$)] and H-out [0.55 mm ($\pm 0.29$)]. This is also true for the results of the application of neural-network-based predictors derived from the resampling methods on the out-of-sample test set in Table 2. This can be explained by the fact that the training set is almost equal to the number of observations $n$. In other words, choosing $d=1$ in Eq. (6) should give the most accurate assessment. To support this analysis, we also apply resampling methods with different configurations and their results are displayed in Table 3.

Table 3 shows that the H-out method for $d=1$ and trials of $t=\mathrm{9,900}$ provides similar results of Jack-1 for $t=99$ in Table 1. It is important to note that there are no advantages to increasing both the number of trials and the number of the validation subsets $d$. Note that the RMSE worsened when increasing $t$ even to $d=1$. This is due to the fact that there are more chances of overfitting or underfitting when increasing the trials $t$. Consequently, in using the resampling process, there are certain chances that one or more predictions will deviate significantly from their true value, as can be seen by the maximum RMSE in Table 3.

Resampling methods can be used with powerful computers to construct inferential procedures in modern statistical data analysis. They replace theoretical derivations required in applying traditional methods by repeatedly resampling the original data and making inferences from the re-samples. In this contribution, we introduce three resampling methods in neural networks, namely Delete-d Jackknife (Jack-d), Delete-1 Jackknife Trials (Jack-1), and Hold-Out Trials (H-out).

The Jack-d method consists of separating the dataset into a subset of size $d$ and training the remaining each of the $(n-d)$ observations one at a time. The subset size d is selected from all of the observations without replacement. The Jack-1 method is an extension of the classical Delete-1 Jackknife by taking $d=1$ and training each of the remaining $n-1$ observations one at a time and by drawing many realizations (trials). Similar to Jack-1, Hold-out is also extended to its H-out version by randomly drawing many realizations, so that each prediction point is repeated a certain number of times. The user knows exactly the number of times each point prediction will be replicated for the cases where Jack-d and Jack-1 are in play. On the other hand, in the case of the H-out method, the user chooses the percentage of division of the sample for each set (training, validation and test) and the number of realizations (trials). An advantage of these resampling techniques is that they allow us to determine the approximate probability distribution for the predictor rather than simply providing a point prediction. A consequence of this is that in addition to all points passing through the test set, they are often tested using different training processes. By nature, neural networks are non-deterministic due to random initialization of the weights, biases, and different optimization techniques. Therefore, the training of a neural network is different for each trial due to this randomness associated with the neural network.

A neural-network-based on resampling methods was applied to a problem of coordinate transformations. The results suggest that, in general, the accuracy of the resampling methods is based on the choice of the size of the test set and the number of realizations (trials). The smaller the test group, the better the accuracy of the neural network model, since the training set is almost equal to the dimension of the dataset. There are no advantages to increasing both the number of trials and the number of the test subsets $d$. Wang and Yu (2020) showed that the optimal $d$ is obtained by approaching its value to $\mathrm{n}/2$ for the case of parameter estimation by the weighted total least squares. In terms of a neural network for prediction, however, our results show that the optimal choice is driven in favor of a smaller $d$. However, in higher dimensional problems, both Jack-1 and Jack-d can lead to an increase in computational time for neural network learning. In that case, H-out with a larger test set should be more appropriate. Here, we are concerned less about the best selection of the neural network parameters. Instead, our key point is with respect to describing the quality of the neural network predictor. In future works, we will address efforts in the application of resampling methods to find the best network parameters.

Furthermore, due to the randomness associated with the neural network, there are certain chances that one or more predictions by the resampling process will significantly deviate from their true value. These deviations are understood to be outliers. Specifically, the Levenberg–Marquardt method used here is very sensitive to the initial network weighs. In addition, it does not consider outliers in the data. In fact, some realizations may have been overfitted. Since we have several predictions for each variable of output of the neural network, it is possible to remove those outliers. Therefore, outlier filtering methods will be studied and discussed in future works.

In the case of having independent and identically distributed variables (IIDs), the choice of Jack-1 should give the most accurate evaluation, since the training set is almost equal to the size of the dataset. The problem is that independence is not always guaranteed. For instance, the independence assumption is invalid for time series data. Therefore, Jack-d would be more suitable in this case. Actually, if the variables are not IIDs with respect to the underlying distribution, then an even more suitable way of resampling is to use the idea of H-out. In future studies, this issue will be addressed in detail.

Some or all data, models, or code generated used during the study are available in a repository online in accordance with funder data retention policies (Rofatto and Matsuoka 2022).