The kernel-balanced equation for deep neural networks
Deep neural networks have shown many fruitful applications in this decade. A
network can get the generalized function through training with a finite
dataset. The degree of generalization is a realization of the proximity scale
in the data space. Specifically, the scale is not clear if the dataset is
complicated. Here we consider a network for the distribution estimation of the
dataset. We show the estimation is unstable and the instability depends on the
data density and training duration. We derive the kernel-balanced equation,
which gives a short phenomenological description of the solution. The equation
tells us the reason for the instability and the mechanism of the scale. The
network outputs a local average of the dataset as a prediction and the scale of
averaging is determined along the equation. The scale gradually decreases along
training and finally results in instability in our case.
Kenichi Nakazato
# The kernel-balanced equation for deep neural networks

###### Abstract

Deep neural networks have shown many fruitful applications in this decade. A network can get the generalized function through training with a finite dataset. The degree of generalization is a realization of the proximity scale in the data space. Specifically, the scale is not clear if the dataset is complicated. Here we consider a network for the distribution estimation of the dataset. We show the estimation is unstable and the instability depends on the data density and training duration. We derive the kernel-balanced equation, which gives a short phenomenological description of the solution. The equation tells us the reason for the instability and the mechanism of the scale. The network outputs a local average of the dataset as a prediction and the scale of averaging is determined along the equation. The scale gradually decreases along training and finally results in instability in our case.

## I Introduction

In the recent decade, data-driven modeling has been empowered with techniques from machine learning. Among them, deep neural networks are the most powerful ones with a large number of applications[1; 2; 3; 4; 5]. Despite the fruitful ones, we do not know much about the mechanism behind them[6; 7]. Specifically, the network can get generalized functions only with the finite dataset. In other words, the network can learn a generalized relation between input and output from the finite one. We can get predictions for unknown inputs but do not fully understand how it works.

A neural network, \(y=f(\mathbf{x},\mathbf{w})\), can be defined with an input, \(\mathbf{x}\), and output, \(y\). In training, we adjust the parameters, \(\mathbf{w}\), of the network so that the pre-defined relation, \(y_{i}=f(\mathbf{x}_{i})\), is satisfied as possible. As the pre-defined relation, we give a dataset, \(\{(\mathbf{x}_{i},y_{i})\}\), in advance. We usually update the parameters step by step along the gradient, \(\sum_{i}\nabla L(|y_{i}-f(\mathbf{x}_{i})|)\), of the minimizing function, \(L\). The name _neural network_ stems from the architectures of the function, \(f\), which is inspired by the brain network. In fact, one of the most famous architecture, convolutional neural networks, is originally a model of retinal structure[8; 9; 10; 11]. In this paper, we focus on them.

We can derive neural tangent kernels, NTKs, or training responses as a theoretical approach to understanding the generalization mechanism[12; 13; 14; 15; 16]. There, we can describe the training response, \(\Theta(\mathbf{x},\mathbf{x}_{i})\), which shows the influence on the output, \(f(\mathbf{x})\), by a training step, in the following equation,

\[f(\mathbf{x},\mathbf{w}+\mathbf{\delta}) \sim f+\frac{\partial f}{\partial\mathbf{w}}\cdot\mathbf{\delta} \tag{1}\] \[= f-\eta\frac{\partial f}{\partial\mathbf{w}}\cdot\frac{\partial f_{i }}{\partial\mathbf{w}}\frac{dL}{df_{i}}\] (2) \[\equiv f-\eta\Theta(\mathbf{x},\mathbf{x}_{i})\frac{dL}{df_{i}}, \tag{3}\]

where the model is trained with a single data, \((\mathbf{x}_{i},y_{i})\), with a minimizing target, \(L\), known as loss function and the parameter, \(\eta\), is a learning rate. We usually have a dataset with many data points, \(\{(\mathbf{x}_{i},y_{i})\}\), and train the model with that. In such a case we can write the equation with the sum of training responses,

\[\Delta f(\mathbf{x})\propto-\sum_{j}\Theta(\mathbf{x},\mathbf{x}_{j})\frac{dL}{df_{j}}. \tag{4}\]

Furthermore, in some cases, we can assume a simple ansatz for the training response with an aging effect. It can be expressed in the following,

\[\Theta(\mathbf{x},\mathbf{x}_{i})\propto t^{-\alpha}K(\mathbf{x},\mathbf{x}_{i}), \tag{5}\]

where the exponent, \(\alpha\), shows aging decay and the response kernel, \(K\), is a positive decreasing function of the distance, \(|\mathbf{x}-\mathbf{x}_{i}|\). The decreasing scale would be determined by the architecture, but we assume it is a similar one with an exponential curve, in this paper.

As a minimizing function, we consider it a more complicated problem than simple supervised training. In standard supervised training, we optimize a network so that the relation, \(y_{i}=f(\mathbf{x}_{i})\), is satisfied for any data point. On the contrary, we want to estimate the distribution of the dataset, \(\{(\mathbf{x}_{i},y_{i})\}\). To do that, we estimate the local
mean, \(\mu(\mathbf{x})\), and standard deviation, \(\sigma(\mathbf{x})\), for each input, \(\mathbf{x}\). As a network, we assume the following one,

\[\mu(\mathbf{x}) =g\circ f(\mathbf{x}) \tag{6}\] \[\sigma(\mathbf{x}) =h\circ f(\mathbf{x}). \tag{7}\]

We have a shared function, \(f\), and specific ones, \(g\) and \(h\), to estimate the mean and standard deviation, respectively. As an optimizing function, we minimize the following,

\[L\equiv-\log(\Pi_{i}\frac{1}{\sigma(\mathbf{x}_{i})\sqrt{2\pi}}\exp(-\frac{1}{2} (\frac{y_{i}-\mu(\mathbf{x}_{i})}{\sigma(\mathbf{x}_{i})})^{2})). \tag{8}\]

In other words, we want to fit with a Gaussian distribution. This type of problem is known as _uncertainty estimation_ in the field of machine learning[17; 18; 19]. We want to know both maximum likelihood estimation, \(\mu(\mathbf{x})\), and uncertainty of that, \(\sigma(\mathbf{x})\), at the same time. Our problem setting is a much simpler one than various applications.

In the case of the standard prediction, we usually minimize the distance, \(|\mu_{i}-y_{i}|\), and the optimal solution is the exact one, \(\mu_{i}=y_{i}\). However, we can have a different mean value because we simultaneously estimate standard deviation, \(\sigma\), in our network. As stated above, training with a data, \((\mathbf{x}_{i},y_{i})\), can influence on another prediction for a data point, \((\mathbf{x}_{j},y_{j})\). In sum, the predicted mean and standard deviation, \(\mu(\mathbf{x}_{i})\) and \(\sigma(\mathbf{x}_{i})\), can be a local estimation of the data distribution. Our main question is how the scale of estimation is determined. We assume two different inputs, \(\mathbf{x}_{a}\), and \(\mathbf{x}_{b}\), should have similar outputs, \(f(\mathbf{x}_{a})\sim f(\mathbf{x}_{b})\), depending on a distance between them, as the nature of prediction. However, we do not know the reality in the cases of deep neural networks. In other words, we want to know how and to what extent the estimation is generalized.

As a hypothesis, we can assume some specific scales of the structure of the dataset are reflected in the prediction. In other words, the scale may be the reflection of the semantic structure of the dataset. However, we use a randomly generated dataset rather than specific public ones, like the other studies in the statistical physics[20; 21]. The advantage of this way is that we can get a universal understanding of the nature of the neural networks independent of the dataset instance.

In the next section II, we introduce our model. There we describe our network and dataset. In addition, we note on training method as well. We show the training dynamics in the section III. We show the estimation is unstable. Furthermore, we introduce a phenomenological description, kernel-balanced equation, of the solution, which explains the instability. It gives answers, the scale and generalization, for our question. Finally, we show that the equation can be redisplayed with dynamics by training response.

## II Model

In general, deep neural networks consist of layers of nonlinear transformations, \(f_{i}\), and an input, \(\mathbf{x}\), and output, \(y\),

\[y=f_{n}\circ\cdots\circ f_{0}(\mathbf{x}). \tag{9}\]

In each layer, we often use combination of linear convolution, \(c_{jk}\), and non-linear activation function, \(R\),

\[\mathbf{h}_{k}=R(b_{k}+\sum_{j}c_{jk}(\mathbf{h}_{j})), \tag{10}\]

where j-th channel of input for a layer, \(\mathbf{h}_{j}\), is transformed into k-th channel of output, \(\mathbf{h}_{k}\). Hidden variables, \(\mathbf{h}_{i}\), have often multi-channels for more degree of freedom of the network. We call this as a convolution layer[8; 9; 10; 11].

Here, we focus on a simple convolutional network with \(n\) layers. We assume the input, \(\mathbf{x}\), is a 1-dimensional bit string with the size, \(b\). In other words, we assume the input, \(\mathbf{x}\), is a binary vector in this paper. Each convolution layer can be defined with the number of output channels, \(s_{c}\), and kernel size, \(s_{k}\), of the convolution. In our model, the number of channels, \(s_{c}\), in each convolution layer is fixed. Needless to say, the output of a mid-layer means the input to the next layer. In the final layer, we usually use a linear network,

\[y=R(\sum_{i}\mathbf{a}_{i}\cdot\mathbf{h}_{i}+b), \tag{11}\]

where the hidden variable, \(\mathbf{h}_{i}\), is the input for the last layer. Non-linear activation, \(R\), is applied after linear transformation with parameters, \(\mathbf{a}_{i}\) and \(b\). In numerical experiments, we use a setting, \(b=8\), \(s_{k}=3\), \(s_{c}=3\), with ELU as an activation function and SGD for the training algorithm with a learning rate, \(\eta=0.1\), without a learning momentum[22; 23; 24; 25].

Since we consider a network with two outputs, \(\mu\), and \(\sigma\), we have two linear networks, \(g\) and \(h\), after the convolution layers, \(f\), in eq. (6) and (7). We call the network as _variance network_, here. In addition, we also consider a simpler one with only one output, \(\mu\), for easier understanding, in eq. (6). We call that as _average network_.

As the dataset, we consider a random bit encoding[16]. The dataset, \(\{(\mathbf{x}_{i},y_{i})\}\), consists of pairs of an input, \(\mathbf{x}_{i}\), and output, \(y_{i}\). We randomly generate \(N\) pairs in advance and use them as a training dataset. We train a model
with the dataset and a loss function, eq. (8). We can focus on training dynamics itself independent of a specific dataset instance by testing randomly generated ones and analyzing its statistical features.

We also consider a simplified model for understanding training dynamics. As stated, training dynamics can be described simply in a simple equation, (3) or (5). The training response, \(\Theta(\mathbf{x},\mathbf{x}_{i})\), is known as a neural tangent kernel and can be constant during training in an infinity limit of network size[12; 13]. Even if the size is finite, it can be represented by a product of a time-dependent term, \(A(t)\), and almost constant kernel, \(K(\mathbf{x},\mathbf{x}_{i})\), like equation (5)[16]. In the case of an average network, we can write down simplified dynamics,

\[\frac{d\mu_{i}}{dt}\propto-\sum_{j}K(\mathbf{x}_{i},\mathbf{x}_{j})\frac{dL}{d\mu_{j}}, \tag{12}\]

where we ignore the time-dependent term in training response. In other words, we consider short-term training dynamics and call it as _response kernel dynamics_.

In many cases, the loss function, \(L\), evaluates the distance between the prediction, \(f_{i}\), and the answer, \(y_{i}\), e.g. mean squared error. In such a case, we can write it as follows,

\[\frac{df_{i}}{dt}\propto\sum_{j}K_{ij}(y_{j}-f_{j}), \tag{13}\]

where the response kernel, \(K_{ij}\), can be seen constant during training as assumption.

## III Results

### 1-point training

We start from the simplest one, where our dataset has only a pair, \((\mathbf{x}_{0},y_{0})\). We train a variance network with that. We call it _1-point training_, here. When the dataset consists of N-pairs, we call it _N-point training_.

Before showing results, we should confirm the form of the loss function, equation 8. Since the loss is the log-likelihood of the Gaussian, it can be rewritten with the summation, easily,

\[L\propto\Sigma_{i}(\sigma(\mathbf{x}_{i})+\frac{1}{2}(\frac{y_{i}-\mu(\mathbf{x}_{i})} {\sigma(\mathbf{x}_{i})})^{2}). \tag{14}\]

The first term, \(\Sigma_{i}\sigma(\mathbf{x}_{i})\), can be minimized by the minimal of the standard deviation, \(\sigma(\mathbf{x}_{i})=0\). However, the second term includes that in its denominator and can be divergent by the zero value. If the first term is minimized faster, the second term can change its value abruptly. In other words, we can see a numerical instability, there. On the contrary, if the scale of the standard deviation can be adjusted moderately, we do not see numerical instability. We can focus on the trajectory of the point, \((y_{i}-\mu(\mathbf{x}_{i}),\sigma(\mathbf{x}_{i}))\), to see the determination of the prediction scale.

We show the results of 1-point training in FIG. 1. In the figure, we show trajectories of training dynamics on the vector field along the loss gradient. All of them start from around the center, \(|\mu_{0}-y_{0}|\sim 0.5\) and \(\sigma_{0}\sim 0.5\), and converged into the optimal point, \(\mu_{0}\sim y_{0}\) and \(\sigma_{0}\sim 0\). As we can confirm, they move almost along the vector field. However, we cannot get to the optimal one, \((\mu_{0},\sigma_{0})=(y_{0},0)\), because that is numerically unstable. One of our network outputs, \(\sigma_{0}\), is in a denominator of the loss function, equation (8). In other words, our formulation of a variance network cannot have the optimal point as a solution for 1-point training. It is a reasonable result because we do not have a meaningful definition of standard deviation only with a data point.

### data density transition

Next, we consider N-point training with a variance network. We sample pairs of random encoding, \((\mathbf{x}_{i}),y_{i}\), with a size, \(N\), and train the network with it. Firstly we want to roughly grasp the feature of training dynamics. To do that we evaluate training results with variance. We can estimate it in two ways, mean squared error, \(V\equiv<(y_{i}-\mu_{i})^{2}>\), and another one, \(V^{*}\equiv<\sigma_{i}^{2}>\). In FIG. 2, we show those estimated variances with different sizes of datasets, \(N\). We can confirm the non-negligible difference between them in the cases with a little dataset, but it is negligible in the cases with a larger dataset. In FIG. 2, we show the difference after training of a fixed epoch, \(e_{mx}=2000\), but it can depend on the duration of the training epoch itself.

In FIG. 3, we show the dependence on the training epoch. In the figure, we plot the difference, \(|V-V^{*}|\), against the size of a dataset, \(N\), and the duration of the training epoch, \(e_{mx}\). As we can see, the difference tends to grow when the duration, \(e_{mx}\), is large. On the other hand, it tends to be reduced when the size, \(N\), is large. In other words, the difference can be negligible when the data density is large but it may grow after enough training.

As we already confirmed, 1-point training is numerically unstable. In addition, N-point training is also un
stable if the training dataset is sparse enough, FIG. 2. If we assume the difference, \(|V-V^{*}|\), stems from numerical instability, the results are reasonable. However, why the instability grows after longer training? We can expect that network outputs, \(\mu_{i}\) and \(\sigma_{i}\), depend not only on the local data point, \((\mathbf{x}_{i},y_{i})\), but also on other ones nearby it. We evaluate the output, \(\mu_{i}\), with a weighted average, \(\sum_{j}\exp(-\alpha|\mathbf{x}_{i}-\mathbf{x}_{j}|)y_{j}/N\). The parameter, \(\alpha\), means a spatial scale of the average. In FIG. 4, we show training dynamics with a dataset, \(N=20\). On the left, we show the dynamics of two variances. On the right, we show the growth of the scale, \(\alpha\). The scale is optimized so that the weighted average and the output, \(\mu_{i}\), should match with each other. The variances show no difference at first. But we see a significant difference between them in the end. At the same time, the scale, \(\alpha\), grows along the training. This suggests a spatial scale of prediction is reduced after enough training.

### kernel-balanced equation

We consider the simplified dynamics, (12) or (13), to understand the growth of scale, \(\alpha\). Firstly, we study the training dynamics of an average network with its simplified one,

\[\frac{d\mu_{i}}{dt}\propto-\sum_{j}K_{ij}(y_{j}-\mu_{j}). \tag{15}\]

Figure 1: Training dynamics with a single point data. A network is trained with a data point, \((\mathbf{x}_{0},y_{0})\), and the error, \(|\mu_{0}-y_{0}|\), and standard deviation, \(\sigma_{0}\), are shown on the map of loss function. In the figure, we show 4 results with different colors and initial conditions, but all trajectories finally end up with numerical instability. The vector field shows the gradient of the loss function. The color of the arrows shows the steepness of the gradient in the log scale. Here we used the setting, learning rate \(\eta=0.1\), the size of input \(b=8\), the number of layers 3. We used SGD and ELU as an optimizer and activation, respectively.

Figure 3: Numerical instability against sample size and training epochs, \(e_{mx}\). We show the difference between two predicted variances. The horizontal axis means the sample size of the training dataset. The vertical axis means training epochs. The difference is shown with color.
The matrix, \(K_{ij}\), consists of kernel distance terms, \(K(|\mathbf{x}_{i}-\mathbf{x}_{j}|)\), and we assume it can be expressed in an exponential form as already introduced,

\[K_{ij}\sim\exp(-\beta|\mathbf{x}_{i}-\mathbf{x}_{j}|). \tag{16}\]

We show eigenvectors and values with an ideally simple case, in FIG. 5. We constructed a matrix, \(K_{ij}=\exp(-|x_{i}-x_{j}|)\), with sorted 100 random values, \(0\leq x_{i}\leq 1\). In other words, we show the features of a response kernel with randomly distanced data points. As we can see, the eigenvalues decrease in a power-law manner. On the other hand, eigenvectors show Gabor wavelet-like forms[26; 27]. Major modes show broader waves than minor ones. This means that training dynamics reduces large scaled spatial error at first. Local error is reduced after enough training. These dynamics can be interpreted in our case as follows,

\[\frac{d\mu_{i}}{dt} \propto \sum_{j}K_{ij}(y_{j}-\mu_{j}) \tag{17}\] \[\sim \sum_{j}K_{ij}(y_{j}-\mu_{i})\] (18) \[\sim \sum_{k}\lambda_{k}\sum_{j}\tilde{K}_{ijk}(y_{j}-\mu_{i}), \tag{19}\]

where the values, \(\tilde{K}_{ijk}\sim\exp(-\beta_{k}|x_{i}-x_{j}|)\), show differently scaled kernel terms. If we can assume a relation, \(\mu_{j}\sim\mu_{i}\pm\delta\), around the point, \(\mathbf{x}_{i}\), we get the equation, 18. This means the local expectation, \(\mu_{i}\), should match the weighted average, \(K_{ij}y_{j}\). Finally, we rewrite it in a manner of multi-scale expansion, eq. 19, using the differently scaled kernels, \(\tilde{K}_{ijk}\), and the weight of each mode, \(\lambda_{k}\). Needless to say, we assume the weight should be larger for the global averaging term, which has a small parameter, \(\beta_{k}\). Thus, the weighted balance of kernel response,

\[\mu_{i}=\frac{\sum_{j}\tilde{K}_{ijk}y_{j}}{\sum_{j}\tilde{K}_{ijk}}, \tag{20}\]

is realized from more global modes to local ones through training.

In the same way, we can write down the kernel-balanced equation for a variance network in the following,

\[\frac{d\mu_{i}}{dt} \propto \sum_{j}\frac{K_{ij}}{2\sigma_{j}^{2}}(y_{j}-\mu_{j}) \tag{21}\] \[\sim \sum_{k}\lambda_{k}\sum_{j}\tilde{K}_{ijk}\frac{y_{j}-\mu_{i}}{2 \sigma_{j}^{2}}\] (22) \[\frac{d\sigma_{i}}{dt} \propto \sum_{j}\frac{K_{ij}}{\sigma_{j}^{3}}((y_{j}-\mu_{j})^{2}-\sigma _{j}^{2})\] (23) \[\sim \sum_{k}\lambda_{k}\sum_{j}\tilde{K}_{ijk}\frac{(y_{j}-\mu_{i})^ {2}-\sigma_{i}^{2}}{\sigma_{j}^{3}}. \tag{24}\]

We can notice the denominator, \(\sigma_{j}\), in the equation (22), as the difference from the previous one. In addition, we have one more equation on the other dynamics, in eq. (23) and (24). These equations suggest we can have

Figure 4: Training dynamics with a training dataset, \(N=20\). (a) Two predicted errors, \(<(y_{i}-\mu_{i})^{2}>\) and \(<\sigma_{i}^{2}>\), are shown. (b) We can approximate the predicted value, \(\mu_{i}\), as a weighted average, \(\sum_{j}\exp(-\alpha|\mathbf{x}_{i}-\mathbf{x}_{j}|)y_{j}/N\). Here we show the evolution of the prediction scale, \(\alpha\).

Figure 5: Eigenvalues and vectors of a response kernel, \(K_{ij}\). We calculated eigenvalues and vectors with a randomly distanced one within a range, \((0,1)\). We generated 100 random values, \(x_{i}\), and made a response kernel, \(\exp(-|x_{j}-x_{i}|)\), with those ones.
numerical instability again, because the point, \(\sigma_{i}=0\), is a local optimal as a result of training.

To confirm the numerically unstable dynamics, we show the training dynamics of N-point training, \(N=20\), in FIG. 6. All the trajectories, \((|y_{i}-\mu_{i}|,\sigma_{i})\), are plotted on the loss gradient, in the same manner as FIG. 1. At a first glance, we notice jumps in them. Needless to say, these suggest numerical instability.

We also show the training dynamics against the training epoch, in FIG. 7. We show some convergent trajectories, in (a), and unstable ones, in (b). All of the trajectories, \(\sigma_{i}\) and \(\mu_{i}\), are plotted, in (c) and (d). We can confirm that jumps happen at the timing when the convergent trajectories approach the optimal.

## IV Discussion

Machine learning is now one of the most powerful ways to construct a model for a given dataset. Among the techniques in the field, deep neural networks play central roles not only in business applications but also in scientific studies. If we train a model with a dataset, \(\{(\mathbf{x}_{i},y_{i})\}\), the trained model, \(f\), can predict outputs for any input, \(y=f(\mathbf{x})\), successfully even if the point, \(\mathbf{x}\), is not included in the dataset. This feature is known as generalization. We know generalization should be a reflection of proximity in the input space. In other words, two outputs, \(f(\mathbf{x}_{a})\) and \(f(\mathbf{x}_{b})\), should be similar if those inputs, \(\mathbf{x}_{a}\) and \(\mathbf{x}_{b}\), are similar with each other. However, the distribution of the dataset is often very complicated and the spatial scale of similarity can be complicated as well. In this paper, the structure of the prediction scale is focused on. We consider neural tangent kernel or training response because they can describe the spatial effect of a training step and be applied to training dynamics. We elucidate the determination mechanism of the prediction scale with simple convolution networks and simplified models for the dynamics.

As an estimation problem for the dataset distribution, we adopt a loss function for the fitting to Gaussian one. We minimize the loss and get an optimal fitting model by training with a dataset. Our model, a variance network, outputs an expectation and standard deviation, \(\mu(\mathbf{x})\) and \(\sigma(\mathbf{x})\). As the dataset, we simply set a randomly generated encoding from a 1D bit-string, \(\mathbf{x}_{i}\), to a binary output, \(y_{i}\). This problem set is a very simple one so that we can understand the dynamics of the prediction scale.

As shown in FIG. 1, the dynamics is numerically unstable in the simplest case, 1-point training, because we cannot determine a meaningful standard deviation. In a similar manner, N-point training suffers from numerical instability, as shown in FIG. 2 and FIG. 3. The numerical instability stems from the reduction of the prediction scale, shown in FIG. 4. The scale gradually reduces along training and finally two types of variances,

Figure 6: Training dynamics with multiple training points, \(N=20\). In the figure, we show all trajectories with different colors. The vector field shows the gradient of a loss function. The color of the arrows shows the steepness of the gradient in the log scale.

Figure 7: A demonstration of numerical instability. We show training dynamics with sample size, \(N=20\). (a) Trajectories of \(\sigma_{i}\) for the group, converged into the state, \(y_{i}^{*}-\mu_{i}=0,\sigma_{i}=0\). In the box, we show the pairs of input, \(x_{i}\), and output, \(y_{i}\). (b) Trajectories of \(\sigma_{i}\), showing first instability. (c) All trajectories of \(\sigma_{i}\). (d) All trajectories of \(\mu_{i}\).
and \(<\sigma_{i}^{2}>\), do not show consistency with each other.

To understand the dynamics of the prediction scale, response kernel dynamics for an average network are studied. The kernel matrix has Gabor wavelet-like eigenvectors with eigenvalues decaying along a power law. This suggests the network learns spatially larger scale patterns at first and local patterns later. These dynamics can be shown in a kernel-balanced equation, in eq. (20). An average network outputs a weighted mean value averaged over the answers, \(y_{i}\). The weight determines the prediction scale and it decreases along the training. In the case of variance network, the solution is not so straightforward, eq. (23) and (24), but the prediction scale decreases along training as well. In addition, the predictions suffer from numerical instability again. Once any prediction, \((\mu_{i},\sigma_{i})\), approaches to the optimal point, \(\sigma_{i}=0\), it destabilizes training dynamics of other predictions, \((\mu_{j},\sigma_{j})\).

The kernel-balanced equation suggests we can understand convolutional networks as a function that outputs the local average. The scale of averaging is not fixed and decreases along the training. In other words, the network can output the same values as the dataset after enough training, if the response kernel does not collapse during it. It is known that the kernel is constant in an ideal condition. Even if the network is not ideal and has a finite size, it can be seen as almost constant in some cases. We still need more studies on kernel stability but we believe the equation should be effective for a more variety of cases.

The problem we consider here is known as uncertainty estimation in the field of machine learning. In reality, we suffer from uncertainty for many reasons in practical usage. If we should assume some noises in the observation, we cannot regard the dataset as the ground truth anymore. Even if we can exclude such a noise in a way, the truths may not be deterministic. In addition, we have some more origins of it in the model and its training. The model can be redundant and have many solutions for the given dataset. We often use non-deterministic training algorithms and this can result in many solutions.

It is known that uncertainty can be divided into two classes, epistemic and aleatoric one[19]. What we consider here is the latter one. In such a case, we can model the dataset as a distribution at most. In this context, what we show here is the insufficiency of Gaussian modeling. Our formulation suggests Laplacian modeling is insufficient as well. On the contrary, a Bayesian approach, in which the parameter has its distribution, can be a solution. In fact, the numerical instability occurs at the optimal point, \(\sigma_{i}=0\), and can be overcome by replacement of it, \(\sigma_{i}\), with another probabilistic one, \(P(\sigma_{i})\). Actually, it is reported that t-distribution is an effective way of that[28]. Our formulation explains the reason for the effectiveness, therefore.

In our days, we often train very large models with large datasets. In such a case, the dataset usually distributes in a non-uniform manner. As pointed out, we often see model variability in such cases[33]. As we have shown with the kernel-balanced equation, the prediction can be different between training phases, especially for minor modes in the dataset. Since the training for minor modes is postponed to a much later stage depending on the initial condition, we can see model variability.

Our instability stems from the shape of the loss function, equation 8. As shown in the kernel-balanced equation, this sort of instability cannot be prevented in a straightforward way. However, when the model and the dataset are very large, training requires much longer epochs and we often need much larger computing time for a training epoch. In such a case, we may not see instability because of short training time, in reality. As a possibility, when we have drop-out layers in the network, we can have a non-zero variance in the prediction for the same input. If the standard deviation for the input, \(\mathbf{x}\), is not zero, we do not see the instability. In this meaning, the instability is not necessarily universal for real applications. However, the analysis with the kernel-balanced equation, especially for the average network, can be applicable to a wide range of applications, because of its simplicity. As the kernel-equation equations shows, training always starts from major modes and proceeds to minor modes later, we can carefully design the dataset density for more efficient training, though we still need more studies.

As such an application, we can apply it for variance estimation without a variance head for prediction of standard deviation, \(\sigma(\mathbf{x})\). Since the output after \(t\)-epochs, \(f_{t}(\mathbf{x})\), approximates a local average dependent on the prediction scale, the difference between two outputs, \((f_{t}(\mathbf{x})-f_{T}(\mathbf{x}))^{2}\), tells us the local variance in the limit, \(T\rightarrow\infty\). We do not know the quantitative precision of this formula yet, but it can be a convenient way for variance estimation.

We are getting more and more data through ubiquitous sensors and it is even available on the web. The field of data science sheds light on the complexity of the real world with them. Actually, astonishingly powerful applications emerge one after another with the aid of such huge datasets and machine powers[1; 2; 3; 4]. Machine learning algorithms are necessary not only in such practical applications but also in scientific studies in the real world's complexity[29; 30; 31]. However, we note that those algorithms still require further understanding. In reality, new technologies often focus on an innovative mathematical formulation and its implementation, but a dynamical
understanding must be necessary as well[16; 32]. We can design an efficient and even safe one with such an understanding. We believe the science of complexity can be an effective approach for the field.

## Acknowledgements

This work is motivated through discussions in Sense-Time Japan, HONDA, and an internship student WL.
