Self-Normalizing Neural Networks

Introduction

Deep Learning has set new records at different benchmarks and led to various commercial applications . Recurrent neural networks (RNNs) achieved new levels at speech and natural language processing, for example at the TIMIT benchmark or at language translation , and are already employed in mobile devices . RNNs have won handwriting recognition challenges (Chinese and Arabic handwriting) and Kaggle challenges, such as the “Grasp-and Lift EEG” competition. Their counterparts, convolutional neural networks (CNNs) excel at vision and video tasks. CNNs are on par with human dermatologists at the visual detection of skin cancer . The visual processing for self-driving cars is based on CNNs , as is the visual input to AlphaGo which has beaten one of the best human GO players . At vision challenges, CNNs are constantly winning, for example at the large ImageNet competition , but also almost all Kaggle vision challenges, such as the “Diabetic Retinopathy” and the “Right Whale” challenges .

However, looking at Kaggle challenges that are not related to vision or sequential tasks, gradient boosting, random forests, or support vector machines (SVMs) are winning most of the competitions. Deep Learning is notably absent, and for the few cases where FNNs won, they are shallow. For example, the HIGGS challenge, the Merck Molecular Activity challenge, and the Tox21 Data challenge were all won by FNNs with at most four hidden layers. Surprisingly, it is hard to find success stories with FNNs that have many hidden layers, though they would allow for different levels of abstract representations of the input .

To robustly train very deep CNNs, batch normalization evolved into a standard to normalize neuron activations to zero mean and unit variance . Layer normalization also ensures zero mean and unit variance, while weight normalization ensures zero mean and unit variance if in the previous layer the activations have zero mean and unit variance. However, training with normalization techniques is perturbed by stochastic gradient descent (SGD), stochastic regularization (like dropout), and the estimation of the normalization parameters. Both RNNs and CNNs can stabilize learning via weight sharing, therefore they are less prone to these perturbations. In contrast, FNNs trained with normalization techniques suffer from these perturbations and have high variance in the training error (see Figure 1). This high variance hinders learning and slows it down. Furthermore, strong regularization, such as dropout, is not possible as it would further increase the variance which in turn would lead to divergence of the learning process. We believe that this sensitivity to perturbations is the reason that FNNs are less successful than RNNs and CNNs.

Self-normalizing neural networks (SNNs) are robust to perturbations and do not have high variance in their training errors (see Figure 1). SNNs push neuron activations to zero mean and unit variance thereby leading to the same effect as batch normalization, which enables to robustly learn many layers. SNNs are based on scaled exponential linear units “SELUs” which induce self-normalizing properties like variance stabilization which in turn avoids exploding and vanishing gradients.

Self-normalizing Neural Networks (SNNs)

We consider the mapping $g$ that maps mean and variance of the activations from one layer to mean and variance of the activations in the next layer

A neural network is self-normalizing if it possesses a mapping $g:\Omega\mapsto\Omega$ for each activation $y$ that maps mean and variance from one layer to the next and has a stable and attracting fixed point depending on $(\omega,\tau)$ in $\Omega$. Furthermore, the mean and the variance remain in the domain $\Omega$, that is $g(\Omega)\subseteq\Omega$, where $\Omega=\{(\mu,\nu)\ |\ \mu\in[\mu_{\min},\mu_{\max}],\nu\in[\nu_{\min},\nu_{\max}]\}$. When iteratively applying the mapping $g$, each point within $\Omega$ converges to this fixed point.

Therefore, we consider activations of a neural network to be normalized, if both their mean and their variance across samples are within predefined intervals. If mean and variance of $\bm{x}$ are already within these intervals, then also mean and variance of $\bm{y}$ remain in these intervals, i.e., the normalization is transitive across layers. Within these intervals, the mean and variance both converge to a fixed point if the mapping $g$ is applied iteratively.

Therefore, SNNs keep normalization of activations when propagating them through layers of the network. The normalization effect is observed across layers of a network: in each layer the activations are getting closer to the fixed point. The normalization effect can also observed be for two fixed layers across learning steps: perturbations of lower layer activations or weights are damped in the higher layer by drawing the activations towards the fixed point. If for all $y$ in the higher layer, $\omega$ and $\tau$ of the corresponding weight vector are the same, then the fixed points are also the same. In this case we have a unique fixed point for all activations $y$. Otherwise, in the more general case, $\omega$ and $\tau$ differ for different $y$ but the mean activations are drawn into $[\mu_{\min},\mu_{\max}]$ and the variances are drawn into $[\nu_{\min},\nu_{\max}]$.

We aim at constructing self-normalizing neural networks by adjusting the properties of the function $g$. Only two design choices are available for the function $g$: (1) the activation function and (2) the initialization of the weights.

For the activation function, we propose “scaled exponential linear units” (SELUs) to render a FNN as self-normalizing. The SELU activation function is given by

SELUs allow to construct a mapping $g$ with properties that lead to SNNs. SNNs cannot be derived with (scaled) rectified linear units (ReLUs), sigmoid units, $\tanh$ units, and leaky ReLUs. The activation function is required to have (1) negative and positive values for controlling the mean, (2) saturation regions (derivatives approaching zero) to dampen the variance if it is too large in the lower layer, (3) a slope larger than one to increase the variance if it is too small in the lower layer, (4) a continuous curve. The latter ensures a fixed point, where variance damping is equalized by variance increasing. We met these properties of the activation function by multiplying the exponential linear unit (ELU) with $\lambda>1$ to ensure a slope larger than one for positive net inputs.

For the weight initialization, we propose $\omega=0$ and $\tau=1$ for all units in the higher layer. The next paragraphs will show the advantages of this initialization. Of course, during learning these assumptions on the weight vector will be violated. However, we can prove the self-normalizing property even for weight vectors that are not normalized, therefore, the self-normalizing property can be kept during learning and weight changes.

These integrals can be analytically computed and lead to following mappings of the moments:

The spectral norm of $\mathcal{J}(0,1)$ (its largest singular value) is $0.7877<1$. That means $g$ is a contraction mapping around the fixed point $(0,1)$ (the mapping is depicted in Figure 2). Therefore, $(0,1)$ is a stable fixed point of the mapping $g$.

A normalized weight vector $\bm{w}$ cannot be ensured during learning. For SELU parameters $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$, we show in the next theorem that if $(\omega,\tau)$ is close to $(0,1)$, then $g$ still has an attracting and stable fixed point that is close to $(0,1)$. Thus, in the general case there still exists a stable fixed point which, however, depends on $(\omega,\tau)$. If we restrict $(\mu,\nu,\omega,\tau)$ to certain intervals, then we can show that $(\mu,\nu)$ is mapped to the respective intervals. Next we present the central theorem of this paper, from which follows that SELU networks are self-normalizing under mild conditions on the weights.

We assume $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$. We restrict the range of the variables to the following intervals $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.95,1.1]$, that define the functions’ domain $\Omega$. For $\omega=0$ and $\tau=1$, the mapping Eq. (3) has the stable fixed point $(\mu,\nu)=(0,1)$, whereas for other $\omega$ and $\tau$ the mapping Eq. (3) has a stable and attracting fixed point depending on $(\omega,\tau)$ in the $(\mu,\nu)$-domain: $\mu\in[-0.03106,0.06773]$ and $\nu\in[0.80009,1.48617]$. All points within the $(\mu,\nu)$-domain converge when iteratively applying the mapping Eq. (3) to this fixed point.

Consequently, feed-forward neural networks with many units in each layer and with the SELU activation function are self-normalizing (see definition 1), which readily follows from Theorem 1. To give an intuition, the main property of SELUs is that they damp the variance for negative net inputs and increase the variance for positive net inputs. The variance damping is stronger if net inputs are further away from zero while the variance increase is stronger if net inputs are close to zero. Thus, for large variance of the activations in the lower layer the damping effect is dominant and the variance decreases in the higher layer. Vice versa, for small variance the variance increase is dominant and the variance increases in the higher layer.

However, we cannot guarantee that mean and variance remain in the domain $\Omega$. Therefore, we next treat the case where $(\mu,\nu)$ are outside $\Omega$. It is especially crucial to consider $\nu$ because this variable has much stronger influence than $\mu$. Mapping $\nu$ across layers to a high value corresponds to an exploding gradient, since the Jacobian of the activation of high layers with respect to activations in lower layers has large singular values. Analogously, mapping $\nu$ across layers to a low value corresponds to an vanishing gradient. Bounding the mapping of $\nu$ from above and below would avoid both exploding and vanishing gradients. Theorem 2 states that the variance of neuron activations of SNNs is bounded from above, and therefore ensures that SNNs learn robustly and do not suffer from exploding gradients.

The proof can be found in the Appendix Section A3. Thus, when mapped across many layers, the variance in the interval $$is mapped to a value below$3$. Consequently, all fixed points$(\mu,\nu)$of the mapping$g$(Eq. (3)) have$\nu<3$. Analogously, Theorem 3 states that the variance of neuron activations of SNNs is bounded from below, and therefore ensures that SNNs do not suffer from vanishing gradients.

In the derivative of the mapping (Eq. (3)), we used the central limit theorem (CLT) to approximate the network inputs $z=\sum_{i=1}^{n}w_{i}x_{i}$ with a normal distribution. We justified normality because network inputs represent a weighted sum of the inputs $x_{i}$, where for Deep Learning $n$ is typically large. The Berry-Esseen theorem states that the convergence rate to normality is $n^{-1/2}$ . In the classical version of the CLT, the random variables have to be independent and identically distributed, which typically does not hold for neural networks. However, the Lyapunov CLT does not require the variable to be identically distributed anymore. Furthermore, even under weak dependence, sums of random variables converge in distribution to a Gaussian distribution .

Experiments

We compare SNNs to other deep networks at different benchmarks. Hyperparameters such as number of layers (blocks), neurons per layer, learning rate, and dropout rate, are adjusted by grid-search for each dataset on a separate validation set (see Section A4). We compare the following FNN methods:

“MSRAinit”: FNNs without normalization and with ReLU activations and “Microsoft weight initialization” .

“BatchNorm”: FNNs with batch normalization .

“LayerNorm”: FNNs with layer normalization .

“WeightNorm”: FNNs with weight normalization .

“ResNet”: Residual networks adapted to FNNs using residual blocks with 2 or 3 layers with rectangular or diavolo shape.

The benchmark comprises 121 classification datasets from the UCI Machine Learning repository from diverse application areas, such as physics, geology, or biology. The size of the datasets ranges between $10$ and $130,000$ data points and the number of features from $4$ to $250$. In abovementioned work , there were methodological mistakes which we avoided here. Each compared FNN method was optimized with respect to its architecture and hyperparameters on a validation set that was then removed from the subsequent analysis. The selected hyperparameters served to evaluate the methods in terms of accuracy on the pre-defined test sets (details on the hyperparameter selection are given in Section A4). The accuracies are reported in the Table A11. We ranked the methods by their accuracy for each prediction task and compared their average ranks. SNNs significantly outperform all competing networks in pairwise comparisons (paired Wilcoxon test across datasets) as reported in Table 1 (left panel).

We further included 17 machine learning methods representing diverse method groups in the comparison and the grouped the data sets into “small” and “large” data sets (for details see Section A4). On 75 small datasets with less than 1000 data points, random forests and SVMs outperform SNNs and other FNNs. On 46 larger datasets with at least 1000 data points, SNNs show the highest performance followed by SVMs and random forests (see right panel of Table 1, for complete results see Tables A12 and A12). Overall, SNNs have outperformed state of the art machine learning methods on UCI datasets with more than 1,000 data points.

Typically, hyperparameter selection chose SNN architectures that were much deeper than the selected architectures of other FNNs, with an average depth of 10.8 layers, compared to average depths of 6.0 for BatchNorm, 3.8 WeightNorm, 7.0 LayerNorm, 5.9 Highway, and 7.1 for MSRAinit networks. For ResNet, the average number of blocks was 6.35. SNNs with many more than 4 layers often provide the best predictive accuracies across all neural networks.

The Tox21 challenge dataset comprises about 12,000 chemical compounds whose twelve toxic effects have to be predicted based on their chemical structure. We used the validation sets of the challenge winners for hyperparameter selection (see Section A4) and the challenge test set for performance comparison. We repeated the whole evaluation procedure 5 times to obtain error bars. The results in terms of average AUC are given in Table 2. In 2015, the challenge organized by the US NIH was won by an ensemble of shallow ReLU FNNs which achieved an AUC of 0.846 . Besides FNNs, this ensemble also contained random forests and SVMs. Single SNNs came close with an AUC of 0.845$\pm$0.003. The best performing SNNs have 8 layers, compared to the runner-ups ReLU networks with layer normalization with 2 and 3 layers. Also batchnorm and weightnorm networks, typically perform best with shallow networks of 2 to 4 layers (Table 2). The deeper the networks, the larger the difference in performance between SNNs and other methods (see columns 5–8 of Table 2). The best performing method is an SNN with 8 layers.

Since a decade, machine learning methods have been used to identify pulsars in radio wave signals . Recently, the High Time Resolution Universe Survey (HTRU2) dataset has been released with 1,639 real pulsars and 16,259 spurious signals. Currently, the highest AUC value of a 10-fold cross-validation is 0.976 which has been achieved by Naive Bayes classifiers followed by decision tree C4.5 with 0.949 and SVMs with 0.929. We used eight features constructed by the PulsarFeatureLab as used previously . We assessed the performance of FNNs using 10-fold nested cross-validation, where the hyperparameters were selected in the inner loop on a validation set (for details on the hyperparameter selection see Section A4). Table 3 reports the results in terms of AUC. SNNs outperform all other methods and have pushed the state-of-the-art to an AUC of $0.98$.

Conclusion

We have introduced self-normalizing neural networks for which we have proved that neuron activations are pushed towards zero mean and unit variance when propagated through the network. Additionally, for activations not close to unit variance, we have proved an upper and lower bound on the variance mapping. Consequently, SNNs do not face vanishing and exploding gradient problems. Therefore, SNNs work well for architectures with many layers, allowed us to introduce a novel regularization scheme, and learn very robustly. On 121 UCI benchmark datasets, SNNs have outperformed other FNNs with and without normalization techniques, such as batch, layer, and weight normalization, or specialized architectures, such as Highway or Residual networks. SNNs also yielded the best results on drug discovery and astronomy tasks. The best performing SNN architectures are typically very deep in contrast to other FNNs.

Acknowledgments

This work was supported by IWT research grant IWT150865 (Exaptation), H2020 project grant 671555 (ExCAPE), grant IWT135122 (ChemBioBridge), Zalando SE with Research Agreement 01/2016, Audi.JKU Deep Learning Center, Audi Electronic Venture GmbH, and the NVIDIA Corporation.

References

The references are provided in Section A7.

Appendix

This appendix is organized as follows: the first section sets the background, definitions, and formulations. The main theorems are presented in the next section. The following section is devoted to the proofs of these theorems. The next section reports additional results and details on the performed computational experiments, such as hyperparameter selection. The last section shows that our theoretical bounds can be confirmed by numerical methods as a sanity check.

The proof of theorem 1 is based on the Banach’s fixed point theorem for which we require (1) a contraction mapping, which is proved in Subsection A3.4.1 and (2) that the mapping stays within its domain, which is proved in Subsection A3.4.2 For part (1), the proof relies on the main Lemma 12, which is a computer-assisted proof, and can be found in Subsection A3.4.1. The validity of the computer-assisted proof is shown in Subsection A3.4.5 by error analysis and the precision of the functions’ implementation. The last Subsection A3.4.6 compiles various lemmata with intermediate results that support the proofs of the main lemmata and theorems.

A1 Background

For neural networks with scaled exponential linear units, the mean is of the activations in the next layer computed according to

and the second moment of the activations in the next layer is computed according to

The parameters $\alpha_{\rm 01}$ and $\lambda_{\rm 01}$ ensure

Figure 2 visualizes the mapping $g$ for $\omega=0$ and $\tau=1$ and $\alpha_{\rm 01}$ and $\lambda_{\rm 01}$ at few pre-selected points. It can be seen that $(0,1)$ is an attracting fixed point of the mapping $g$.

A2 Theorems

Theorem 1 shows that the mapping $g$ defined by Eq. (4) and Eq. (5) exhibits a stable and attracting fixed point close to zero mean and unit variance. Theorem 1 establishes the self-normalizing property of self-normalizing neural networks (SNNs). The stable and attracting fixed point leads to robust learning through many layers.

We assume $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$. We restrict the range of the variables to the domain $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.95,1.1]$. For $\omega=0$ and $\tau=1$, the mapping Eq. (4) and Eq. (5) has the stable fixed point $(\mu,\nu)=(0,1)$. For other $\omega$ and $\tau$ the mapping Eq. (4) and Eq. (5) has a stable and attracting fixed point depending on $(\omega,\tau)$ in the $(\mu,\nu)$-domain: $\mu\in[-0.03106,0.06773]$ and $\nu\in[0.80009,1.48617]$. All points within the $(\mu,\nu)$-domain converge when iteratively applying the mapping Eq. (4) and Eq. (5) to this fixed point.

A2.2 Theorem 2: Decreasing Variance from Above

The variance decreases in $$and all fixed points$(\mu,\nu)$of mapping Eq. (5) and Eq. (4) have$\nu<3$.

A2.3 Theorem 3: Increasing Variance from Below

We consider $\lambda=\lambda_{\rm 01}$, $\alpha=\alpha_{\rm 01}$ and the two domains $\Omega_{1}^{-}=\{(\mu,\omega,\nu,\tau)\ |\ -0.1\leqslant\mu\leqslant 0.1,-0.1\leqslant\omega\leqslant 0.1,0.05\leqslant\nu\leqslant 0.16,0.8\leqslant\tau\leqslant 1.25\}$ and $\Omega_{2}^{-}=\{(\mu,\omega,\nu,\tau)\ |\ -0.1\leqslant\mu\leqslant 0.1,-0.1\leqslant\omega\leqslant 0.1,0.05\leqslant\nu\leqslant 0.24,0.9\leqslant\tau\leqslant 1.25\}$.

A3 Proofs of the Theorems

We have to show that the mapping $g$ defined by Eq. (4) and Eq. (5) has a stable and attracting fixed point close to $(0,1)$. To proof this statement and Theorem 1, we apply the Banach fixed point theorem which requires (1) that $g$ is a contraction mapping and (2) that $g$ does not map outside the function’s domain, concretely:

Let $(X,d)$ be a non-empty complete metric space with a contraction mapping $f:X\to X$. Then $f$ has a unique fixed-point $x_{f}\in X$ with $f(x_{f})=x_{f}$. Every sequence $x_{n}=f(x_{n-1})$ with starting element $x_{0}\in X$ converges to the fixed point: $x_{n}\xrightarrow[n\to\infty]{\ }x_{f}$.

Contraction mappings are functions that map two points such that their distance is decreasing:

A function $f:X\to X$ on a metric space $X$ with distance $d$ is a contraction mapping, if there is a $0\leqslant\delta<1$, such that for all points $\bm{u}$ and $\bm{v}$ in $X$: $d(f(\bm{u}),f(\bm{v}))\leqslant\delta d(\bm{u},\bm{v})$.

To show that $g$ is a contraction mapping in $\Omega$ with distance $\|.\|_{2}$, we use the Mean Value Theorem for $u,v\in\Omega$

in which $M$ is an upper bound on the spectral norm the Jacobian $\mathcal{H}$ of $g$. The spectral norm is given by the largest singular value of the Jacobian of $g$. If the largest singular value of the Jacobian is smaller than 1, the mapping $g$ of the mean and variance to the mean and variance in the next layer is contracting. We show that the largest singular value is smaller than 1 by evaluating the function for the singular value $S(\mu,\omega,\nu,\tau,\lambda,\alpha)$ on a grid. Then we use the Mean Value Theorem to bound the deviation of the function $S$ between grid points. To this end, we have to bound the gradient of $S$ with respect to $(\mu,\omega,\nu,\tau)$. If all function values plus gradient times the deltas (differences between grid points and evaluated points) is still smaller than 1, then we have proofed that the function is below 1 (Lemma 12). To show that the mapping does not map outside the function’s domain, we derive bounds on the expressions for the mean and the variance (Lemma 13). Section A3.4.1 and Section A3.4.2 are concerned with the contraction mapping and the image of the function domain of $g$, respectively.

With the results that the largest singular value of the Jacobian is smaller than one (Lemma 12) and that the mapping stays in the domain $\Omega$ (Lemma 13), we can prove Theorem 1. We first recall Theorem 1:

We assume $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$. We restrict the range of the variables to the domain $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.95,1.1]$. For $\omega=0$ and $\tau=1$, the mapping Eq. (4) and Eq. (5) has the stable fixed point $(\mu,\nu)=(0,1)$. For other $\omega$ and $\tau$ the mapping Eq. (4) and Eq. (5) has a stable and attracting fixed point depending on $(\omega,\tau)$ in the $(\mu,\nu)$-domain: $\mu\in[-0.03106,0.06773]$ and $\nu\in[0.80009,1.48617]$. All points within the $(\mu,\nu)$-domain converge when iteratively applying the mapping Eq. (4) and Eq. (5) to this fixed point.

According to Lemma 12 the mapping $g$ (Eq. (4) and Eq. (5)) is a contraction mapping in the given domain, that is, it has a Lipschitz constant smaller than one. We showed that $(\mu,\nu)=(0,1)$ is a fixed point of the mapping for $(\omega,\tau)=(0,1)$.

The domain is compact (bounded and closed), therefore it is a complete metric space. We further have to make sure the mapping $g$ does not map outside its domain $\Omega$. According to Lemma 13, the mapping maps into the domain $\mu\in[-0.03106,0.06773]$ and $\nu\in[0.80009,1.48617]$.

Now we can apply the Banach fixed point theorem given in Theorem 4 from which the statement of the theorem follows. ∎

A3.2 Proof of Theorem 2

The variance decreases in $$and all fixed points$(\mu,\nu)$of mapping Eq. (5) and Eq. (4) have$\nu<3$.

Therefore we have according to Theorem 16:

We set $x=\nu\tau$ and $y=\mu\omega$ and obtain

The derivative to this sub-function with respect to $y$ is

The inequality follows from Lemma 24, which states that $ze^{z^{2}}\operatorname{erfc}(z)$ is monotonically increasing in $z$. Therefore the sub-function is increasing in $y$. The derivative to this sub-function with respect to $x$ is

The sub-function is increasing in $x$, since the derivative is larger than zero:

First inequality: We applied Lemma 22 two times.

Equalities factor out $\sqrt{2}\sqrt{x}$ and reformulate.

Second inequality part 2: we show that for $a=\frac{1}{10}\left(\sqrt{\frac{960+169\pi}{\pi}}-13\right)$ following holds: $\frac{8x}{\pi}-\left(a^{2}+2a(x+y)\right)\geqslant 0$. We have $\frac{\partial}{\partial x}\frac{8x}{\pi}-\left(a^{2}+2a(x+y)\right)=\frac{8}{\pi}-2a>0$ and $\frac{\partial}{\partial y}\frac{8x}{\pi}-\left(a^{2}+2a(x+y)\right)=-2a<0$. Therefore the minimum is at border for minimal $x$ and maximal $y$:

for $a=\frac{1}{10}\left(\sqrt{\frac{960+169\pi}{\pi}}-13\right)>0.878$.

Equalities only solve square root and factor out the resulting terms $(2(2x+y)+1)$ and $(2(x+y)+0.878)$.

We set $\alpha=\alpha_{01}$ and multiplied out. Thereafter we also factored out $x$ in the numerator. Finally a quadratic equations was solved.

The sub-function has its minimal value for minimal $x=\nu\tau=1.5\cdot 0.8=1.2$ and minimal $y=\mu\omega=-1\cdot 0.1=-0.1$. We further minimize the function

We now divide the $\mu$-domain into $-1\leqslant\mu\leqslant 0$ and $0\leqslant\mu\leqslant 1$. Analogously we divide the $\omega$-domain into $-0.1\leqslant\omega\leqslant 0$ and $0\leqslant\omega\leqslant 0.1$. In this domains $g$ is strictly monotonically.

For all domains $g$ is strictly monotonically decreasing in $\nu$ and strictly monotonically increasing in $\tau$. Note that we now consider the range $3\leqslant\nu\leqslant 16$. For the maximal value of $g$ we set $\nu=3$ (we set it to 3!) and $\tau=1.25$.

We consider now all combination of these domains:

$-1\leqslant\mu\leqslant 0$ and $-0.1\leqslant\omega\leqslant 0$:

$g$ is decreasing in $\mu$ and decreasing in $\omega$. We set $\mu=-1$ and $\omega=-0.1$.

$-1\leqslant\mu\leqslant 0$ and $0\leqslant\omega\leqslant 0.1$:

$g$ is increasing in $\mu$ and decreasing in $\omega$. We set $\mu=0$ and $\omega=0$.

$0\leqslant\mu\leqslant 1$ and $-0.1\leqslant\omega\leqslant 0$:

$g$ is decreasing in $\mu$ and increasing in $\omega$. We set $\mu=0$ and $\omega=0$.

$0\leqslant\mu\leqslant 1$ and $0\leqslant\omega\leqslant 0.1$:

$g$ is increasing in $\mu$ and increasing in $\omega$. We set $\mu=1$ and $\omega=0.1$.

Therefore the maximal value of $g$ is $-0.0180173$.

A3.3 Proof of Theorem 3

We consider $\lambda=\lambda_{\rm 01}$, $\alpha=\alpha_{\rm 01}$ and the two domains $\Omega_{1}^{-}=\{(\mu,\omega,\nu,\tau)\ |\ -0.1\leqslant\mu\leqslant 0.1,-0.1\leqslant\omega\leqslant 0.1,0.05\leqslant\nu\leqslant 0.16,0.8\leqslant\tau\leqslant 1.25\}$ and $\Omega_{2}^{-}=\{(\mu,\omega,\nu,\tau)\ |\ -0.1\leqslant\mu\leqslant 0.1,-0.1\leqslant\omega\leqslant 0.1,0.05\leqslant\nu\leqslant 0.24,0.9\leqslant\tau\leqslant 1.25\}$ .

The mean value theorem states that there exists a $t\in$ for which

Therefore we are interested to bound the derivative of the $\xi$-mapping Eq. (13) with respect to $\nu$:

The sub-term Eq. (322) enters the derivative Eq. (47) with a negative sign! According to Lemma 18, the minimal value of sub-term Eq. (322) is obtained by the largest largest $\nu$, by the smallest $\tau$, and the largest $y=\mu\omega=0.01$. Also the positive term $\operatorname{erfc}\left(\frac{\mu\omega}{\sqrt{2}\sqrt{\nu\tau}}\right)+2$ is multiplied by $\tau$, which is minimized by using the smallest $\tau$. Therefore we can use the smallest $\tau$ in whole formula Eq. (47) to lower bound it.

First we consider the domain $0.05\leqslant\nu\leqslant 0.16$ and $0.8\leqslant\tau\leqslant 1.25$. The factor consisting of the exponential in front of the brackets has its smallest value for $e^{-\frac{0.01\cdot 0.01}{2\cdot 0.05\cdot 0.8}}$. Since $\operatorname{erfc}$ is monotonically decreasing we inserted the smallest argument via $\operatorname{erfc}\left(-\frac{0.01}{\sqrt{2}\sqrt{0.05\cdot 0.8}}\right)$ in order to obtain the maximal negative contribution. Thus, applying Lemma 18, we obtain the lower bound on the derivative:

Next we consider the domain $0.05\leqslant\nu\leqslant 0.24$ and $0.9\leqslant\tau\leqslant 1.25$. The factor consisting of the exponential in front of the brackets has its smallest value for $e^{-\frac{0.01\cdot 0.01}{2\cdot 0.05\cdot 0.9}}$. Since $\operatorname{erfc}$ is monotonically decreasing we inserted the smallest argument via $\operatorname{erfc}\left(-\frac{0.01}{\sqrt{2}\sqrt{0.05\cdot 0.9}}\right)$ in order to obtain the maximal negative contribution.

Thus, applying Lemma 18, we obtain the lower bound on the derivative:

A3.4 Lemmata and Other Tools Required for the Proofs

In this section, we show that the largest singular value of the Jacobian of the mapping $g$ is smaller than one. Therefore, $g$ is a contraction mapping. This is even true in a larger domain than the original $\Omega$. We do not need to restrict $\tau\in[0.95,1.1]$, but we can extend to $\tau\in[0.8,1.25]$. The range of the other variables is unchanged such that we consider the following domain throughout this section: $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.8,1.25]$.

The definition of the entries of the Jacobian $\mathcal{J}$ is:

If the largest singular value of the Jacobian is smaller than 1, then the spectral norm of the Jacobian is smaller than 1. Then the mapping Eq. (4) and Eq. (5) of the mean and variance to the mean and variance in the next layer is contracting.

We show that the largest singular value is smaller than 1 by evaluating the function $S(\mu,\omega,\nu,\tau,\lambda,\alpha)$ on a grid. Then we use the Mean Value Theorem to bound the deviation of the function $S$ between grid points. Toward this end we have to bound the gradient of $S$ with respect to $(\mu,\omega,\nu,\tau)$. If all function values plus gradient times the deltas (differences between grid points and evaluated points) is still smaller than 1, then we have proofed that the function is below 1.

The singular values of the $2\times 2$ matrix

We used an explicit formula for the singular values . We now set ${\mathcal{H}}_{11}=a_{11},{\mathcal{H}}_{12}=a_{12},{\mathcal{H}}_{21}=a_{21},{\mathcal{H}}_{22}=a_{22}$ to obtain a formula for the largest singular value of the Jacobian depending on $(\mu,\omega,\nu,\tau,\lambda,\alpha)$. The formula for the largest singular value for the Jacobian is:

where $\mathcal{J}$ are defined in Eq. (62) and we left out the dependencies on $(\mu,\omega,\nu,\tau,\lambda,\alpha)$ in order to keep the notation uncluttered, e.g. we wrote ${\mathcal{J}}_{11}$ instead of ${\mathcal{J}}_{11}(\mu,\omega,\nu,\tau,\lambda,\alpha)$.

In order to bound the gradient of the singular value, we have to bound the derivatives of the Jacobian entries ${\mathcal{J}}_{11}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, ${\mathcal{J}}_{12}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, ${\mathcal{J}}_{21}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, and ${\mathcal{J}}_{22}(\mu,\omega,\nu,\tau,\lambda,\alpha)$ with respect to $\mu$, $\omega$, $\nu$, and $\tau$. The values $\lambda$ and $\alpha$ are fixed to $\lambda_{\rm 01}$ and $\alpha_{\rm 01}$. The 16 derivatives of the 4 Jacobian entries with respect to the 4 variables are:

The following bounds on the absolute values of the derivatives of the Jacobian entries ${\mathcal{J}}_{11}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, ${\mathcal{J}}_{12}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, ${\mathcal{J}}_{21}(\mu,\omega,\nu,\tau,\lambda,\alpha)$, and ${\mathcal{J}}_{22}(\mu,\omega,\nu,\tau,\lambda,\alpha)$ with respect to $\mu$, $\omega$, $\nu$, and $\tau$ hold:

where we used that (a) $J_{11}$ is strictly monotonically increasing in $\mu\omega$ and $|2-\operatorname{erfc}\left(\frac{0.01}{\sqrt{2}\sqrt{\nu\tau}}\right)|\leqslant 1.00584$ and (b) Lemma 47 that $|e^{\mu\omega+\frac{\nu\tau}{2}}\operatorname{erfc}\left(\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}\right)|\leqslant e^{0.01+\frac{0.64}{2}}\operatorname{erfc}\left(\frac{0.01+0.64}{\sqrt{2}\sqrt{0.64}}\right)=0.587622$ ∎

For the first term we have $0.434947\leqslant e^{\mu\omega+\frac{\nu\tau}{2}}\operatorname{erfc}\left(\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}\right)\leqslant 0.587622$ after Lemma 47 and for the second term $0.582677\leqslant\sqrt{\frac{2}{\pi\nu\tau}}e^{-\frac{\mu^{2}\omega^{2}}{2\nu\tau}}\leqslant 0.997356$, which can easily be seen by maximizing or minimizing the arguments of the exponential or the square root function. The first term scaled by $\alpha$ is $0.727780\leqslant\alpha e^{\mu\omega+\frac{\nu\tau}{2}}\operatorname{erfc}\left(\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}\right)\leqslant 0.983247$ and the second term scaled by $\alpha-1$ is $0.392294\leqslant(\alpha-1)\sqrt{\frac{2}{\pi\nu\tau}}e^{-\frac{\mu^{2}\omega^{2}}{2\nu\tau}}\leqslant 0.671484$. Therefore, the absolute difference between these terms is at most $0.983247-0.392294$ leading to the derived bound.

We assume $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$. We restrict the range of the variables to the domain $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.8,1.25]$.

We set $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$ and restrict the range of the variables to $\mu\in[\mu_{\rm min},\mu_{\rm max}]=[-0.1,0.1]$, $\omega\in[\omega_{\rm min},\omega_{\rm max}]=[-0.1,0.1]$, $\nu\in[\nu_{\rm min},\nu_{\rm max}]=[0.8,1.5]$, and $\tau\in[\tau_{\rm min},\tau_{\rm max}]=[0.8,1.25]$.

The absolute values of derivatives of the largest singular value $S(\mu,\omega,\nu,\tau,\lambda,\alpha)$ given in Eq. (71) with respect to $(\mu,\omega,\nu,\tau)$ are bounded as follows:

The Jacobian of our mapping Eq. (4) and Eq. (5) is defined as

from which follows using the bounds from Lemma 5:

Derivative of the singular value w.r.t. $\mu$:

where we used the results from the lemmata 5, 6, 7, and 9.

Derivative of the singular value w.r.t. $\omega$:

Derivative of the singular value w.r.t. $\nu$:

where we used the results from the lemmata 5, 6, 7, and 9.

Derivative of the singular value w.r.t. $\tau$:

We set $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$ and restrict the range of the variables to $\mu\in[\mu_{\rm min},\mu_{\rm max}]=[-0.1,0.1]$, $\omega\in[\omega_{\rm min},\omega_{\rm max}]=[-0.1,0.1]$, $\nu\in[\nu_{\rm min},\nu_{\rm max}]=[0.8,1.5]$, and $\tau\in[\tau_{\rm min},\tau_{\rm max}]=[0.8,1.25]$.

The distance of the singular value at $S(\mu,\omega,\nu,\tau,\lambda_{\rm 01},\alpha_{\rm 01})$ and that at $S(\mu+\Delta\mu,\omega+\Delta\omega,\nu+\Delta\nu,\tau+\Delta\tau,\lambda_{\rm 01},\alpha_{\rm 01})$ is bounded as follows:

The mean value theorem states that a $t\in$ exists for which

We now apply Lemma 10 which gives bounds on the derivatives, which immediately gives the statement of the lemma. ∎

We set $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$ and restrict the range of the variables to $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.8,1.25]$.

The the largest singular value of the Jacobian is smaller than 1:

Therefore the mapping Eq. (4) and Eq. (5) is a contraction mapping.

We set $\Delta\mu=0.0068097371$, $\Delta\omega=0.0008292885$, $\Delta\nu=0.0009580840$, and $\Delta\tau=0.0007323095$.

For a grid with grid length $\Delta\mu=0.0068097371$, $\Delta\omega=0.0008292885$, $\Delta\nu=0.0009580840$, and $\Delta\tau=0.0007323095$, we evaluated the function Eq. (71) for the largest singular value in the domain $\mu\in[-0.1,0.1]$, $\omega\in[-0.1,0.1]$, $\nu\in[0.8,1.5]$, and $\tau\in[0.8,1.25]$. We did this using a computer. According to Subsection A3.4.5 the precision if regarding error propagation and precision of the implemented functions is larger than $10^{-13}$. We performed the evaluation on different operating systems and different hardware architectures including CPUs and GPUs. In all cases the function Eq. (71) for the largest singular value of the Jacobian is bounded by $0.9912524171058772$.

A3.4.2 Lemmata for proofing Theorem 1 (part 2): Mapping within domain

We further have to investigate whether the the mapping Eq. (4) and Eq. (5) maps into a predefined domains.

The mapping Eq. (4) and Eq. (5) map for $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$ into the domain $\mu\in[-0.03106,0.06773]$ and $\nu\in[0.80009,1.48617]$ with $\omega\in[-0.1,0.1]$ and $\tau\in[0.95,1.1]$.

for all $\omega\in[-0.1,0.1]$ and $\tau\in[0.95,1.1]$.

A3.4.3 Lemmata for proofing Theorem 2: The variance is contracting

We consider the main sub-function of the derivate of second moment, $J22$ (Eq. (62)):

that depends on $\mu\omega$ and $\nu\tau$, therefore we set $x=\nu\tau$ and $y=\mu\omega$. Algebraic reformulations provide the formula in the following form:

For $\lambda=\lambda_{\rm 01}$ and $\alpha=\alpha_{\rm 01}$, we consider the domain $-1\leqslant\mu\leqslant 1$, $-0.1\leqslant\omega\leqslant 0.1$, $1.5\leqslant\nu\leqslant 16$, and, $0.8\leqslant\tau\leqslant 1.25$.

For $x$ and $y$ we obtain: $0.8\cdot 1.5=1.2\leqslant x\leqslant 20=1.25\cdot 16$ and $0.1\cdot(-1)=-0.1\leqslant y\leqslant 0.1=0.1\cdot 1$. In the following we assume to remain within this domain.

For $1.2\leqslant x\leqslant 20$ and $-0.1\leqslant y\leqslant 0.1$,

is smaller than zero, is strictly monotonically increasing in $x$, and strictly monotonically decreasing in $y$ for the minimal $x=12/10=1.2$.

The graph of the subfunction in the specified domain is displayed in Figure A3.

For $\lambda=\lambda_{\rm 01}$, $\alpha=\alpha_{\rm 01}$, $-1\leqslant\mu\leqslant 1$, $-0.1\leqslant\omega\leqslant 0.1$ $1.5\leqslant\nu\leqslant 16$, and $0.8\leqslant\tau\leqslant 1.25$, we first show that the derivative is positive and then upper bound it.

is negative. This expression multiplied by positive factors is subtracted in the derivative Eq. (118), therefore, the whole term is positive. The remaining term

of the derivative Eq. (118) is also positive according to Lemma 21. All factors outside the brackets in Eq. (118) are positive. Hence, the derivative Eq. (118) is positive.

First equality brings the expression into a shape where we can apply Lemma 15 for the the function Eq. (115).

First inequality: The overall factor $\tau$ is bounded by 1.25.

Second inequality: We apply Lemma 15. According to Lemma 15 the function Eq. (115) is negative. The largest contribution is to subtract the most negative value of the function Eq. (115), that is, the minimum of function Eq. (115). According to Lemma 15 the function Eq. (115) is strictly monotonically increasing in $x$ and strictly monotonically decreasing in $y$ for $x=1.2$. Therefore the function Eq. (115) has its minimum at minimal $x=\nu\tau=1.5\cdot 0.8=1.2$ and maximal $y=\mu\omega=1.0\cdot 0.1=0.1$. We insert these values into the expression.

Third inequality: We use for the whole expression the maximal factor $e^{-\frac{\mu^{2}\omega^{2}}{2\nu\tau}}<1$ by setting this factor to 1.

Fourth inequality: $\operatorname{erfc}$ is strictly monotonically decreasing. Therefore we maximize its argument to obtain the least value which is subtracted. We use the minimal $x=\nu\tau=1.5\cdot 0.8=1.2$ and the maximal $y=\mu\omega=1.0\cdot 0.1=0.1$.

Sixth inequality: evaluation of the terms.

is larger than zero when the term $\sqrt{\pi}\alpha e^{\frac{(\mu\omega+\nu\tau)^{2}}{2\nu\tau}}\operatorname{erfc}\left(\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}\right)-\frac{\sqrt{2}(\alpha-1)}{\sqrt{\nu\tau}}$ is larger than zero. This term obtains its minimal value at $\mu\omega=0.01$ and $\nu\tau=16\cdot 1.25$, which can easily be shown using the Abramowitz bounds (Lemma 22) and evaluates to $0.16$, therefore ${\mathcal{J}}_{12}>0$ in $\Omega^{+}$.

A3.4.4 Lemmata for proofing Theorem 3: The variance is expanding

We consider functions in $\mu\omega$ and $\nu\tau$, therefore we set $x=\mu\omega$ and $y=\nu\tau$.

For $\lambda=\lambda_{\rm 01}$ and $\alpha=\alpha_{\rm 01}$, we consider the domain $-0.1\leqslant\mu\leqslant 0.1$, $-0.1\leqslant\omega\leqslant 0.1$ $0.00875\leqslant\nu\leqslant 0.7$, and $0.8\leqslant\tau\leqslant 1.25$.

For $x$ and $y$ we obtain: $0.8\cdot 0.00875=0.007\leqslant x\leqslant 0.875=1.25\cdot 0.7$ and $0.1\cdot(-0.1)=-0.01\leqslant y\leqslant 0.01=0.1\cdot 0.1$. In the following we assume to be within this domain.

In this domain, we consider the main sub-function of the derivate of second moment in the next layer, $J22$ (Eq. (62)):

that depends on $\mu\omega$ and $\nu\tau$, therefore we set $x=\nu\tau$ and $y=\mu\omega$. Algebraic reformulations provide the formula in the following form:

For $0.007\leqslant x\leqslant 0.875$ and $-0.01\leqslant y\leqslant 0.01$, the function

smaller than zero, is strictly monotonically increasing in $x$ and strictly monotonically increasing in $y$ for the minimal $x=0.007=0.00875\cdot 0.8$, $x=0.56=0.7\cdot 0.8$, $x=0.128=0.16\cdot 0.8$, and $x=0.216=0.24\cdot 0.9$ (lower bound of $0.9$ on $\tau$).

For $\lambda=\lambda_{\rm 01}$, $\alpha=\alpha_{\rm 01}$ and the domain $-0.1\leqslant\mu\leqslant 0.1$, $-0.1\leqslant\omega\leqslant 0.1$, $0.00875\leqslant\nu\leqslant 0.7$, and $0.8\leqslant\tau\leqslant 1.25$. We are interested of the derivative of

The derivative of the equation above with respect to

$\tau$ is smaller than zero for maximal $\nu=0.7$, $\nu=0.16$, and $\nu=0.24$ (with $0.9\leqslant\tau$);

$y=\mu\omega$ is larger than zero for $\nu\tau=0.008750.8=0.007$, $\nu\tau=0.70.8=0.56$, $\nu\tau=0.160.8=0.128$, and $\nu\tau=0.24\cdot 0.9=0.216$.

A3.4.5 Computer-assisted proof details for main Lemma 12 in Section A3.4.1.

We investigate the error propagation for the singular value (Eq. (71)) if the function arguments $\mu,\omega,\nu,\tau$ suffer from numerical imprecisions up to $\epsilon$. To this end, we first derive error propagation rules based on the mean value theorem and then we apply these rules to the formula for the singular value.

For a real-valued function $f$ which is differentiable in the closed interval $[a,b]$, there exists $t\in$ with

It follows that for computation with error $\Delta x$, there exists a $t\in$ with

Therefore the increase of the norm of the error after applying function $f$ is bounded by the norm of the gradient $\left\|\nabla f(\bm{x}+t\Delta\bm{x})\right\|$.

We now compute for the functions, that we consider their gradient and its 2-norm:

$f(\bm{x})=x_{1}+x_{2}$ and $\nabla f(\bm{x})=(1,1)$, which gives $\left\|\nabla f(\bm{x})\right\|=\sqrt{2}$.

$f(\bm{x})=x_{1}-x_{2}$ and $\nabla f(\bm{x})=(1,-1)$, which gives $\left\|\nabla f(\bm{x})\right\|=\sqrt{2}$.

$f(\bm{x})=x_{1}x_{2}$ and $\nabla f(\bm{x})=(x_{2},x_{1})$, which gives $\left\|\nabla f(\bm{x})\right\|=\left\|\bm{x}\right\|$.

$f(\bm{x})=\frac{x_{1}}{x_{2}}$ and $\nabla f(\bm{x})=\left(\frac{1}{x_{2}},-\frac{x_{1}}{x_{2}^{2}}\right)$, which gives $\left\|\nabla f(\bm{x})\right\|=\frac{\left\|\bm{x}\right\|}{x_{2}^{2}}$.

$f(x)=\sqrt{x}$ and $f^{\prime}(x)=\frac{1}{2\sqrt{x}}$, which gives $\left|f^{\prime}(x)\right|=\frac{1}{2\sqrt{x}}$.

$f(x)=\exp(x)$ and $f^{\prime}(x)=\exp(x)$, which gives $\left|f^{\prime}(x)\right|=\exp(x)$.

If the values $\mu,\omega,\nu,\tau$ have a precision of $\epsilon$, the singular value (Eq. (71)) evaluated with the formulas given in Eq. (62) and Eq. (71) has a precision better than $292\epsilon$.

This means for a machine with a typical precision of $2^{-52}=2.220446\cdot 10^{-16}$, we have the rounding error $\epsilon\approx 10^{-16}$, the evaluation of the singular value (Eq. (71)) with the formulas given in Eq. (62) and Eq. (71) has a precision better than $10^{-13}>292\epsilon$.

We have the numerical precision $\epsilon$ of the parameters $\mu,\omega,\nu,\tau$, that we denote by $\Delta\mu,\Delta\omega,\Delta\nu,\Delta\tau$ together with our domain $\Omega$.

With the error propagation rules that we derived in Subsection A3.4.5, we can obtain bounds for the numerical errors on the following simple expressions:

Using these bounds on the simple expressions, we can now calculate bounds on the numerical errors of compound expressions:

Subsequently, we can use the above results to get bounds for the numerical errors on the Jacobian entries (Eq. (62)), applying the rules from Subsection A3.4.5 again:

We will show that our computations are correct up to 3 ulps. For our implementation in GNU C library and the hardware architectures that we used, the precision of all mathematical functions that we used is at least one ulp. The term “ulp” (acronym for “unit in the last place”) was coined by W. Kahan in 1960. It is the highest precision (up to some factor smaller 1), which can be achieved for the given hardware and floating point representation.

“Ulp$(x)$ is the gap between the two finite floating-point numbers nearest $x$, even if $x$ is one of them. (But ulp(NaN) is NaN.)”

“an ulp in $x$ is the distance between the two closest straddling floating point numbers $a$ and $b$, i.e. those with $a\leqslant x\leqslant b$ and $a\not=b$ assuming an unbounded exponent range.”

In the literature we find also slightly different definitions .

“IEEE-754 mandates four standard rounding modes:”

“Round-to-nearest: $r(x)$ is the floating-point value closest to $x$ with the usual distance; if two floating-point value are equally close to $x$, then $r(x)$ is the one whose least significant bit is equal to zero.”

“IEEE-754 standardises 5 operations: addition (which we shall note $\oplus$ in order to distinguish it from the operation over the reals), subtraction ($\ominus$), multiplication ($\otimes$), division ($\oslash$), and also square root.”

Consequently, the IEEE-754 standard guarantees that addition, subtraction, multiplication, division, and squared root is precise up to one ulp.

“With the Intel486 processor and Intel 387 math coprocessor, the worst- case, transcendental function error is typically $3$ or $3.5$ ulps, but is some- times as large as $4.5$ ulps.”

According to https://www.mirbsd.org/htman/i386/man3/exp.htm and http://man.openbsd.org/OpenBSD-current/man3/exp.3:

“exp$(x)$, log$(x)$, expm1$(x)$ and log1p$(x)$ are accurate to within an ulp”

which is the same for freebsd https://www.freebsd.org/cgi/man.cgi?query=exp&sektion=3&apropos=0&manpath=freebsd:

“The values of exp(0), expm1(0), exp2(integer), and pow(integer, integer) are exact provided that they are representable. Otherwise the error in these functions is generally below one ulp.”

The same holds for “FDLIBM” http://www.netlib.org/fdlibm/readme:

“FDLIBM is intended to provide a reasonably portable (see assumptions below), reference quality (below one ulp for major functions like sin,cos,exp,log) math library (libm.a).”

A3.4.6 Intermediate Lemmata and Proofs

Since we focus on the fixed point $(\mu,\nu)=(0,1)$, we assume for our whole analysis that $\alpha=\alpha_{\rm 01}$ and $\lambda=\lambda_{\rm 01}$. Furthermore, we restrict the range of the variables $\mu\in[\mu_{\rm min},\mu_{\rm max}]=[-0.1,0.1]$, $\omega\in[\omega_{\rm min},\omega_{\rm max}]=[-0.1,0.1]$, $\nu\in[\nu_{\rm min},\nu_{\rm max}]=[0.8,1.5]$, and $\tau\in[\tau_{\rm min},\tau_{\rm max}]=[0.8,1.25]$.

For bounding different partial derivatives we need properties of different functions. We will bound a the absolute value of a function by computing an upper bound on its maximum and a lower bound on its minimum. These bounds are computed by upper or lower bounding terms. The bounds get tighter if we can combine terms to a more complex function and bound this function. The following lemmata give some properties of functions that we will use in bounding complex functions.

Throughout this work, we use the error function $\operatorname{erf}(x):=\frac{1}{\sqrt{\pi}}\int_{-x}^{x}e^{-t^{2}}$ and the complementary error function $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$.

$\exp(x)$ is strictly monotonically increasing from at $-\infty$ to $\infty$ at $\infty$ and has positive curvature.

According to its definition $\operatorname{erfc}(x)$ is strictly monotonically decreasing from 2 at $-\infty$ to 0 at $\infty$.

Next we introduce a bound on $\operatorname{erfc}$:

The statement follows immediately from (page 298, formula 7.1.13). ∎

$e^{x^{2}}\operatorname{erfc}(x)$ is strictly monotonically decreasing for $x>0$ and has positive curvature (positive 2nd order derivative), that is, the decreasing slowes down.

A graph of the function is displayed in Figure A5.

The derivative of $e^{x^{2}}\operatorname{erfc}(x)$ is

Thus $e^{x^{2}}\operatorname{erfc}(x)$ is strictly monotonically decreasing for $x>0$.

The second order derivative of $e^{x^{2}}\operatorname{erfc}(x)$ is

Again using Lemma 22 (first inequality), we get

For the last inequality we added 1 in the numerator in the square root which is subtracted, that is, making a larger negative term in the numerator. ∎

The function $xe^{x^{2}}\operatorname{erfc}(x)$ has the sign of $x$ and is monotonically increasing to $\frac{1}{\sqrt{\pi}}$.

The derivative of $xe^{x^{2}}\operatorname{erfc}(x)$ is

We apply Lemma 22 to $x\operatorname{erfc}(x)e^{x^{2}}$ and divide the terms of the lemma by $x$, which gives

For $\lim_{x\to\infty}$ both the upper and the lower bound go to $\frac{1}{\sqrt{\pi}}$. ∎

$h_{11}(\mu,\omega)=\mu\omega$ is monotonically increasing in $\mu\omega$. It has minimal value $t_{11}=-0.01$ and maximal value $T_{11}=0.01$.

$h_{22}(\nu,\tau)=\nu\tau$ is monotonically increasing in $\nu\tau$ and is positive. It has minimal value $t_{22}=0.64$ and maximal value $T_{22}=1.875$.

𝜇𝜔𝜈𝜏2𝜈𝜏\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}). $h_{1}(\mu,\omega,\nu,\tau)=\frac{\mu\omega+\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}$ is larger than zero and increasing in both $\nu\tau$ and $\mu\omega$. It has minimal value $t_{1}=0.5568$ and maximal value $T_{1}=0.9734$.

The derivative of the function $\frac{\mu\omega+x}{\sqrt{2}\sqrt{x}}$ with respect to $x$ is

since $x>0.8\cdot 0.8$ and $\mu\omega<0.1\cdot 0.1$. ∎

𝜇𝜔2𝜈𝜏2𝜈𝜏\frac{\mu\omega+2\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}). $h_{2}(\mu,\omega,\nu,\tau)=\frac{\mu\omega+2\nu\tau}{\sqrt{2}\sqrt{\nu\tau}}$ is larger than zero and increasing in both $\nu\tau$ and $\mu\omega$. It has minimal value $t_{2}=1.1225$ and maximal value $T_{2}=1.9417$.

The derivative of the function $\frac{\mu\omega+2x}{\sqrt{2}\sqrt{x}}$ with respect to $x$ is

$h_{3}(\mu,\omega,\nu,\tau)=\frac{\mu\omega}{\sqrt{2}\sqrt{\nu\tau}}$ monotonically decreasing in $\nu\tau$ and monotonically increasing in $\mu\omega$. It has minimal value $t_{3}=-0.0088388$ and maximal value $T_{3}=0.0088388$.

$h_{4}(\mu,\omega,\nu,\tau)=\left(\frac{\mu\omega}{\sqrt{2}\sqrt{\nu\tau}}\right)^{2}$ has a minimum at 0 for $\mu=0$ or $\omega=0$ and has a maximum for the smallest $\nu\tau$ and largest $\left|\mu\omega\right|$ and is larger or equal to zero. It has minimal value $t_{4}=0$ and maximal value $T_{4}=0.000078126$.

$\frac{\sqrt{\frac{2}{\pi}}(\alpha-1)}{\sqrt{\nu\tau}}>0$ and decreasing in $\nu\tau$.

Statements follow directly from elementary functions square root and division. ∎

$2-\operatorname{erfc}\left(\frac{\mu\omega}{\sqrt{2}\sqrt{\nu\tau}}\right)>0$ and decreasing in $\nu\tau$ and increasing in $\mu\omega$.

Statements follow directly from Lemma 21 and $\operatorname{erfc}$. ∎

For $\lambda=\lambda_{\rm 01}$ and $\alpha=\alpha_{\rm 01}$, $\sqrt{\frac{2}{\pi}}\left(\frac{(\alpha-1)\mu\omega}{(\nu\tau)^{3/2}}-\frac{\alpha}{\sqrt{\nu\tau}}\right)<0$ and increasing in both $\nu\tau$ and $\mu\omega$.

We consider the function $\sqrt{\frac{2}{\pi}}\left(\frac{(\alpha-1)\mu\omega}{x^{3/2}}-\frac{\alpha}{\sqrt{x}}\right)$, which has the derivative with respect to $x$:

This derivative is larger than zero, since

The last inequality follows from $\alpha-\frac{3\cdot 0.1\cdot 0.1(\alpha-1)}{0.8\cdot 0.8}>0$ for $\alpha=\alpha_{\rm 01}$.

We next consider the function $\sqrt{\frac{2}{\pi}}\left(\frac{(\alpha-1)x}{(\nu\tau)^{3/2}}-\frac{\alpha}{\sqrt{\nu\tau}}\right)$, which has the derivative with respect to $x$:

which has as derivative with respect to $x$:

The derivative of the term $3(\alpha-1)\mu^{2}\omega^{2}-x(-\alpha+\alpha\mu\omega+1)-\alpha x^{2}$ with respect to $x$ is $-1+\alpha-\mu\omega\alpha-2\alpha x<0$, since $2\alpha x>1.6\alpha$. Therefore the term is maximized with the smallest value for $x$, which is $x=\nu\tau=0.8\cdot 0.8$. For $\mu\omega$ we use for each term the value which gives maximal contribution. We obtain an upper bound for the term: