Numerical influence of ReLU’(0) on backpropagation
Year: 2022

Abstract:
In theory, the choice of ReLU′(0) in [0, 1] for a neural network has a negligible influence both on backpropagation and training. Yet, in the real world, 32 bits default precision combined with the size of deep learning problems makes it a hyperparameter of training methods. We investigate the importance of the value of ReLU′(0) for several precision levels (16, 32, 64 bits), on various networks (fully connected, VGG, ResNet) and datasets (MNIST, CIFAR10, SVHN, ImageNet). We observe considerable variations of backpropagation outputs which occur around half of the time in 32 bits precision. The effect disappears with double precision, while it is systematic at 16 bits. For vanilla SGD training, the choice ReLU′(0) = 0 seems to be the most efficient. For our experiments on ImageNet the gain in test accuracy over ReLU′(0) = 1 was more than 10 points (two runs). We also evidence that reconditioning approaches as batch-norm or ADAM tend to buffer the influence of ReLU′(0)’s value. Overall, the message we convey is that algorithmic differentiation of nonsmooth problems potentially hides parameters that could be tuned advantageously.

heading: 1 Introduction
text: Nonsmooth algorithmic differentiation: The training phase of neural networks relies on firstorder methods such as Stochastic Gradient Descent (SGD) [14, 9] and crucially on algorithmic differentiation [15]. The fast “differentiator” used in practice to compute mini-batch descent directions is the backpropagation algorithm [29, 4]. Although designed initially for differentiable problems, it is applied indifferently to smooth or nonsmooth networks. In the nonsmooth case this requires surrogate derivatives at the non regularity points. We focus on the famous ReLU := max(0, ·), for which a value for s = ReLU′(0) has to be chosen. A priori, any value in [0, 1] bears a variational sense as it corresponds to a subgradient [27]. Yet in most libraries s = 0 is chosen; it is the case for TensorFlow [2], PyTorch [26] or Jax [10]. Why this choice? What would be the impact of a different value of s? How this interacts with other training strategies? We will use the notation backprops to denote backpropagation implemented with ReLU′(0) = s for any given real number s.1
1Definition in Section 2. This can be coded explicitly or cheaply emulated by backprop0. Indeed, considering fs : x 7→ ReLU(x) + s(ReLU(−x)− ReLU(x) + x) = ReLU(x), we have backprop0[fs](0) = s.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
What does backpropagation compute? A popular thinking is that the impact of s = ReLU′(0) should be limited as it concerns the value at a single point. This happens to be exact in theory but for surprisingly complex reasons related to Whitney stratification (see Section 2 and references therein):
— For a given neural network architecture, backprops outputs a gradient almost everywhere, independently of the choice of s ∈ R, see [7] and [5] for a detailed treatment of ReLU networks. — Under proper randomization of the initialization and the step-size of SGD, with probability 1 the value of s has no impact during training with backprops as a gradient surrogate, see [8, 6].
In particular, the set of network parameters such that backprop0 6= backprop1 is negligible and the vast majority of SGD sequences produced by training neural networks are not impacted by changing the value of s = ReLU′(0). These results should in principle close the question about the role of ReLU′(0). Yet, this is not what we observe in practice with the default settings of usual libraries.
A surprising experiment and the bifurcation zone: An empirical observation on MNIST triggered our investigations: consider a fully connected ReLU network, and let (θk,0)k∈N and (θk,1)k∈N be two training weights sequences obtained using SGD, with the same random initialization and the same mini-batch sequence but choosing respectively s = 0 and s = 1 in PyTorch [26]. As depicted in Figure 1 (left), the two sequences differ. A closer look shows that the sudden divergence is related to what we call the bifurcation zone, i.e., the set S0,1 of weights such that backprop0 6= backprop1. As mentioned previously this contradicts theory which predicts that the bifurcation zone is met with null probability during training. This contradiction is due to the precision of floating point operations and, to a lesser extent, to the size of deep learning problems. Indeed, rounding schemes used for inexact arithmetics over the reals (which set to zero all values below a certain threshold), may “thicken” negligible sets. This is precisely what happens in our experiments (Figure 1).
The role of numerical precision: Contrary to numerical linear algebra libraries such as numpy, which operates by default under 64 bits precision, the default choice in PyTorch is 32 bits precision (as in TensorFlow or Jax). We thus modulated machine precision to evaluate the importance of the bifurcation zone in Section 3. In 32 bits, we observed that empirically the zone occupies from about 10% to 90% of the network weight space. It becomes invisible in 64 bits precision even for quite large architectures, while, in 16 bits, it systematically fills up the whole space. Although numerical precision is the primary cause of the apparition of the zone, we identify other factors such as network size, sample size. Let us mention that low precision neural network training is a topic of active research [33, 20, 11, 16], see also [28] for an overview. Our investigations are complementary as we focus on the interplay between nonsmoothness and numerical precision.
Impact on learning: The next natural question is to measure the impact of the choice of s = ReLU′(0) in machine learning terms. In Section 4, we conduct extensive experiments combining different architectures (fully connected, VGG, ResNet), datasets (MNIST, SVHN, CIFAR10,
ImageNet) and other learning factors (Adam optimizer, batch normalization, dropout). In 32 bits numerical precision (default in PyTorch or Tensorflow), we consistently observe an effect of choosing s 6= 0. We observe a significant decrease in terms of test accuracy as |s| increases; this can be explained by chaotic oscillatory behaviors induced during training. In some cases gradients even explode and learning cannot be achieved. The sensitivity to this effect highly depends on the problem at hand, in particular, on the network structure and the dataset. On the other hand the choice s = 0 seems to be the most stable. We also observe that both batch normalization [21] and—to a lesser degree—the Adam optimizer [22] considerably mitigate this effect. All our experiments are done using PyTorch [26]; we provide the code to generate all figures presented in this manuscript.
One important message is that, even if the default choice s = 0 seems to be the most stable, our experiments show a counter-intuitive phenomenon that illustrates the interplay between numerical precision and nonsmoothness, and calls for caution when learning nonsmooth networks.
Outline of the paper: In Section 2 we recall elements of nonsmooth algorithmic differentiation which are key to understand the mathematics underlying our experiments. Most results were published in [7, 8]; we provide more detailed pointers to this literature in Appendix A.1. In Section 3 we describe investigations of the bifurcation zone and factors influencing its importance using fully connected networks on the MNIST dataset. Neural network training experiments are detailed in Section 4 with further experimental details and additional experiments reported in Appendix D.

heading: 2 On the mathematics of backpropagation for ReLU networks
text: This section recalls recent advances on the mathematical properties of backpropagation, with in particular the almost sure absence of impact of ReLU′(0) on the learning phase (assuming exact arithmetic over the reals). The main mathematical tools are conservative fields developed in [7]; we provide a simplified overview which is applicable to a wide class of neural networks.

heading: 2.1 Empirical risk minimization and backpropagation
text: Given a training set {(xi, yi)}i=1...N , the supervised training of a neural network f consists in minimizing the empirical risk:
min θ∈RP
J(θ) := 1
N N∑ i=1 `(f(xi, θ), yi) (1)
where θ ∈ RP are the network’s weight parameters and ` is a loss function. The problem can be rewritten abstractly, for each i = 1, . . . , N and θ ∈ RP , `(f(xi, θ), yi) = li(θ) where the function li : RP → R is a composition of the form
li = gi,M ◦ gi,M−1 ◦ . . . ◦ gi,1 (2) where for each j = 1, . . . ,M , the function gi,j is locally Lipschitz with appropriate input and output dimensions. A concrete example of what the functions gi,j look like is given in Appendix A.2 in the special case of fully connected ReLU networks. Furthermore, we associate with each gi,j a generalized Jacobian Ji,j which is such that Ji,j(w) = Jacgi,j (w) whenever gi,j is differentiable at w and Jac denotes the usual Jacobian. The value of Ji,j at the nondifferentiability loci of gi,j can be arbitrary for the moment. The backpropagation algorithm is an automatized implementation of the rules of differential calculus: for each i = 1, . . . , N , we have
backprop li(θ) = Ji,M (gi,M−1 ◦ . . . ◦ gi,1(θ))× Ji,M−1 (gi,M−2 ◦ . . . ◦ gi,1(θ))× . . .× Ji,1(θ). (3)
Famous autograd libraries such as PyTorch [26] or TensorFlow [1] implement dictionaries of functions g with their corresponding generalized Jacobians J , as well as efficient numerical implementation of the quantities defined in (3).

heading: 2.2 ReLU networks training
text: Our main example is based on the function ReLU: R→ R defined by ReLU(x) = max{x, 0}. It is differentiable save at the origin and satisfies ReLU′(x) = 0 for x < 0 and ReLU′(x) = 1 for x > 0.
The value of the derivative at x = 0 could be arbitrary in [0, 1] as we have ∂ReLU(0) = [0, 1], where ∂ denotes the subgradient from convex analysis. Let us insist on the fact that any value within [0, 1] has a variational meaning. For example PyTorch and TensorFlow use ReLU′(0) = 0.
Throughout the paper, and following the lines of [8], we say that a function g : Rp → Rq is elementary log-exp if it can be described by a finite compositional expression involving the basic operations +,−,×, / as well as the exponential and logarithm functions, inside their domains of definition. Examples include the logistic loss log(1 + e−x) on R, the multivariate Gaussian density (2π)−K/2 exp(− ∑K k=1 x 2 k/2) on RK , and the softmax function ( exi/ ∑K k=1 e xk ) 1≤i≤K on R
K . The expressions x3/x2 and exp(−1/x2) do not fit this definition because evaluation at x = 0 cannot be defined by the formula. Roughly speaking a computer evaluating an elementary log-exp expression should not output any NaN error for any input.
Definition 1 (ReLU network training). Assume that in (1), the function ` is elementary log-exp, the network f has an arbitrary structure and the functions gi,j in (2) are either elementary log-exp or consist in applying the ReLU function to some coordinates of their input. We then call the problem in (1) a ReLU network training problem. Furthermore, for any s ∈ R, we denote by backprops the quantity defined in (3) when ReLU′(0) = s for all ReLU functions involved in the composition (2).
Other nonsmooth activation functions: The ReLU operation actually allows to express many other types of nonsmoothness such as absolute values, maxima, minima, quantiles (med for median) or soft-thresholding (st). For any x, y, z ∈ R, we indeed have |x| = ReLU(2x) − x, 2 max(x, y) = (|x − y| + x + y), min(x, y) = −max(−x,−y), med(x, y, z) = min(max(x, y),max(x, z),max(y, z)), st(x) = ReLU(x− 1)− ReLU(−x− 1). Definition 1 is thus much more general than it may seem since it allows, for example, to express convolutional neural networks with max-pooling layers such as VGG or ResNet architectures which correspond to the models considered in the experimental section (although we do not re-program pooling using ReLU). The following theorem is due to [7], with an alternative version in [8].
Theorem 1 (Backprop returns a gradient a.e.). Consider a ReLU network training problem (1) as in Definition 1 and T ≥ 1. Define S ⊂ RP as the complement of the set{ θ ∈ RP , li differentiable at θ, backpropsli(θ) = ∇li(θ), ∀i ∈ {1, . . . , N}, s ∈ [−T, T ] } .
Then S is contained in a finite union of embedded differentiable manifolds of dimension at most

heading: P − 1 (and in particular has Lebesgue measure zero).
text: Although this theorem looks natural, this is a nontrivial result about the backpropagation algorithm that led to the introduction of conservative fields [7]. It implies that all choices for s = ReLU′(0) in [0, 1] = ∂ReLU(0) are essentially equivalent modulo a negligible set S. Perhaps more surprisingly, s can be chosen arbitrarily in R without breaking this essential property of backprop. The set S is called the bifurcation zone throughout the manuscript. For ReLU network training problems, the bifurcation zone is a Lebesgue zero set and is actually contained locally in a finite union of smooth objects of dimension strictly smaller than the ambient dimension. This geometric result reveals a surprising link between backpropagation and Whitney stratifications, as described in [7, 8]. In any case the bifurcation zone is completely negligible. Note that the same result holds if we allow each different call to the ReLU function to use different values for ReLU′(0).

heading: 2.3 ReLU network training with SGD
text: Let (Bk)k∈N denote a sequence of mini-batches with sizes |Bk| ⊂ {1, . . . , N} for all k and αk > 0 the learning rate. Given initial weights θ0 ∈ RP and any parameter s ∈ R, the SGD training procedure of f consists in applying the recursion
θk+1,s = θk,s − γ αk |Bk| ∑ n∈Bk backprops[`(f(xn, θk,s), yn)], (4)
where γ > 0 is a step-size parameter. Note that we explicitly wrote the dependency of the sequence in s = ReLU′(0). According to Theorem 1 if the initialization θ0 is chosen randomly, say, uniformly in a ball, a hypercube, or with iid Gaussian entries, then with probability 1, θ0 does not fall in the bifurcation zone S. Intuitively, since S is negligible, the odds that one of the iterates produced by the
algorithm fall on S are very low. As a consequence, varying s in the recursion (4) does not modify the sequence. This rationale is actually true for almost all values of γ. This provides a rigorous statement of the idea that the choice of ReLU′(0) “does not affect” neural network training. The following result is based on arguments developed in [8], see also [6] in a probabilistic context. Theorem 2 (ReLU′(0) does not impact training with probability 1). Consider a ReLU network training problem as in Definition 1. Let (Bk)k∈N be a sequence of mini-batches with |Bk| ⊂ {1, . . . , N} for all k and αk > 0 the associated learning rate parameter. Choose θ0 uniformly at random in a hypercube and γ uniformly in a bounded interval I ⊂ R+. Let s ∈ R, set θ0,s = θ0, and consider the recursion given in (4). Then, with probability one, for all k ∈ N, θk,s = θk,0.

heading: 3 Surprising experiments on a simple feedforward network

heading: 3.1 ReLU′(0) has an impact
text: Even though the ReLU activation function is non-differentiable at 0, autograd libraries such as PyTorch [26] or TensorFlow [1] implement its derivative with ReLU′(0) = 0. What happens if one chooses ReLU′(0) = 1? The popular answer to this question is that it should have no effect. Theorems 1 and 2 provide a formal justification which is far from being trivial.
A 32 bits MNIST experiment We ran a simple experiment to confirm this answer. We initialized two fully connected neural networks f0 and f1 of size 784×2000×128×10 with the same weights θ0,0 = θ0,1 ∈ RP which are chosen at random. Using the MNIST dataset [24], we trained f0 and f1 with the same sequence of mini-batches (Bk)k∈N (minibatch size 128), using the recursion in (4) for s = 0 and s = 1 and with a fixed αk = 1, and γ chosen uniformly at random in [0.01, 0.012]. At each iteration k, we computed the sum ‖θk,0 − θk,1‖1 of the absolute differences between the coordinates of θk,0 and θk,1. As a sanity check, we actually computed θk,0 a second time, denoting this by θ̃k,0, using a third network to control for sources of divergence in our implementation. Results are reported in Figure 1. The experiment was run using PyTorch [26] on a CPU.
ReLU’(0) has an impact First we observe no difference between θk,0 and θ̃k,0, which shows that we have controlled all possible sources of divergence in PyTorch. Second, while no differences between θ0,0 and θ0,1 is expected (Theorem 2), we observe a sudden deviation of ‖θk,0 − θk,1‖1 at iteration 65 which then increases in a smooth way. The deviation is sudden as the norm is exactly zero before iteration 65 and jumps above one after. Therefore this cannot be explained by an accumulation of small rounding errors throughout iterations, as this would produce a smooth divergence starting at the first iteration. So this suggests that there is a special event at iteration 65.
The center part of Figure 1 displays the minimal absolute value of neurons of the first hidden layer evaluated on the current mini-batch, before the application of ReLU. It turns out that at iteration 65, this minimal value is exactly 0, resulting in a drop in the center of Figure 1. This means that the divergence is actually due to an iterate of the sequence falling exactly on the bifurcation zone. According to Theorem 2, this event is so exceptional that it should never been seen.
Practice and Theory: Numerical precision vs Genericity This contradiction can be solved as follows: the minimal absolute value in Figure 1 oscillates between 10−6 and 10−8 which is roughly the value of machine precision in 32 bits float arithmetic. This machine precision value is of the order 10−16 in 64 bits floating arithmetic which is orders of magnitude smaller than the typical value represented in Figure 1. And indeed, performing the same experiment in 64 bits precision, the divergence of ‖θk,0 − θk,1‖1 disappears and the algorithm can actually be run for many epochs without any divergence between the two sequences. This is represented in Figure 7 in Appendix B. We also report similar results using ReLU6 [19] in place of ReLU on a similar network.

heading: 3.2 Relative volume of the bifurcation zone and relative gradient variation
text: The previous experiment suggests that mini-batch SGD algorithm (4) crossed the bifurcation zone:
S01 = {θ ∈ RP : ∃i ∈ {1, . . . , N}, backprop0[li](θ) 6= backprop1[li](θ)} ⊂ S. (5) This unlikely phenomenon is due to finite numerical precision arithmetic which thickens the negligible set S01 in proportion to the machine precision threshold. This is illustrated on the right of
Figure 1, which represents the bifurcation zone at iteration 65 in a randomly generated hyperplane (uniformly among 2 dimensional hyperplanes) centered around the weight parameters which generated the bifurcation in the previous experiment (evaluation on the same mini-batch). The white area corresponds to some entries being exactly zero, i.e., below the 32 bits machine precision threshold, before application of ReLU. On the other hand, in 64 bits precision, the same representation is much smoother and does not contain exact zeros (see Figure 7 in Appendix B). This confirms that the level of floating point arithmetic precision explains the observed divergence. We now estimate the relative volume of this bifurcation zone by Monte Carlo sampling (see Appendix C for details). All experiments are performed using PyTorch [26] on GPU.
Experimental procedure – weight randomization: We randomly generate a set of parameters {θj}Mj=1, with M = 1000, for a fully connected network architecture f composed of L hidden layers. Given two consecutive layers, respectively composed of m and m′ neurons, the weights of the corresponding affine transform are drawn independently, uniformly in [−α, α] where α = √ 6/ √ m. This is the default weight initialization scheme in PyTorch (Kaiming-Uniform [17]). Given this sample of parameters, iterating on the whole MNIST dataset, we approximate the proportion of θj for which backprop0(li)(θj) 6= backprop1(li)(θj) for some i, for different networks and under different conditions (see Appendix C for details).
Impact of the floating-point precision: Using a fixed architecture of three hidden layers of 256 neurons each, we empirically measured the relative volume of S01 using the above experimental procedure, varying the numerical precision. Table 1 reports the obtained estimates. As shown in Table 1, line 1, at 16 bits floating-point precision, all drawn weights {θi}Mi=1 fall within S01. In sharp contrast, when using a 64 bits precision, none of the sampled weights belong to S01. This proportion is 40% in 32 bits precision. For the rare impacted mini-batches (Table 1 line 2), the relative change in norm is above a factor 20, higher in 16 bits precision (Table 1 line 3). These results confirm that the floating-point arithmetic precision is key in explaining the impact of ReLU′(0) during backpropagation, impacting both frequency and magnitude of the differences.
Impact of sample and mini-batch size: Given a training set {(xn, yn)}n=1...N and a random variable θ ∈ RP with distribution Pθ, we estimate the probability that θ ∈ S01 ⊂ S. Intuitively, this probability should increase with sample size N . We perform this estimation for a fixed architecture of 4 hidden layers composed of 256 neurons each while varying the sample size. Results are reported in Figure 2. For both the 16 and the 32 bits floating-point precisions, our estimation indeed increases with the sample size while we do not find any sampled weights in S01 in 64 bits precision. We also found that the influence of mini-batch size is not significative.
Impact of network size: To evaluate the impact of the network size, we carried out a similar Monte Carlo experiment, varying the depth and width of the network. Firstly, we fixed the number of hidden layers to 3. Following the same experimental procedure, we empirically estimated the relative volume of S01, varying the width of the hidden layers. The results, reported in Figure 2, show that increasing the number of neurons by layer increases the probability to find a random θ ∼ Pθ in S01 for both 16 and 32 floating-point precision. In 64 bits, even with the largest width tested (1024 neurons), no sampled weight parameter is found in S01. Similarly we repeated the experiment varying the network depth and fixing, this time, the width of the layers to 256 neurons. Anew, the
results, reported in Figure 2, show that increasing the network depth increases the probability that θ ∼ Pθ belongs to S01 for both 16 and 32 bits floating-point precision while this probability is zero in 64 bits precision. This shows that the size of the network, both in terms of number of layers and size of layers is positively related to the effect of the choice of s in backpropagation. On the other hand, the fact that neither the network depth, width, or the number of samples impact the 64 bits case suggests that, within our framework, numerical precision is the primary factor of deviation.

heading: 4 Consequences for learning

heading: 4.1 Benchmarks and implementation
text: Datasets and networks We further investigate the effect of the phenomenon described in Section 3 in terms of learning using the CIFAR10 dataset [23] and the VGG11 architecture [31]. To confirm our findings in alternative settings, we also use the MNIST [24], SVHN [25] and ImageNet [12] datasets, fully connected networks (3 hidden layers of size 2048), and the ResNet18 and ResNet50 architectures [18]. Additional details on the different architectures and datasets are found in Appendix D.1. By default, unless stated otherwise, we use the SGD optimizer. We also investigated the effect of batch normalization [21], dropout [32], the Adam optimizer [22] as well as numerical precision. All the experiments in this section are run in PyTorch [26] on GPU. For each training experiment presented in this section (except ImageNet experiments), we use the optuna library [3] to tune learning rates for each experimental condition; see also Appendix D.2.

heading: 4.2 Effect on training and test errors
text: We first consider training a VGG11 architecture on CIFAR10 using the SGD optimizer. For different values of ReLU′(0), we train this network ten times with random initializations under 32 bits arithmetic precision. The results are depicted in Figure 3. Without batch normalization, varying the value of ReLU′(0) beyond a magnitude of 0.1 has a detrimental effect on the test accuracy, resulting in a concave shaped curve with maximum around ReLU′(0) = 0. On the other hand, the decrease of the training loss with the number of epochs suggests that choosing ReLU′(0) 6= 0 induces jiggling behaviors with possible sudden jumps during training. Note that the choice ReLU′(0) = 0 leads to a smooth decrease and that the magnitude of the jumps for other values is related to the magnitude of the chosen value. This is consistent with the interpretation that changing the value of ReLU′(0) has an excitatory effect on training. We observed qualitatively similar behaviors for a fully connected network on MNIST and a ResNet18 on CIFAR10 (Appendix D). Sensitivity to the magnitude of ReLU′(0) depends on the network architecture: our fully connected network on MNIST is less sensitive to this value than VGG11 and ResNet18. The latter shows a very high sensitivity since for values above 0.2, training becomes very unstable and almost impossible.
These experiments complement the preliminary results obtained in Section 3. In particular choosing different values for ReLU′(0) has an effect, it induces a chaotic behavior during training which affects test accuracy. The default choice ReLU′(0) = 0 seems to provide the best performances.
To conclude, we conducted four training experiments for ResNet50 on the ImageNet dataset [12] using the SGD optimizer. These were conducted with fixed learning rate, contrary to results reported above. We observe that switching from ReLU′(0) = 0 to 1 results in a massive drop from around 75% to 63% or 55% for two runs.

heading: 4.3 Mitigating factors: numerical precision, batch-normalization and Adam
text: We analyze below the combined effects of the variations of ReLU′(0) with numerical precision or classical reconditioning methods: Adam, batch-normalization and dropout.
Batch-normalization: As represented in Figure 3, batch normalization [21] not only allows to attain higher test accuracy, but it also completely filters out the effect of the choice of the value of ReLU′(0), resulting in a flat shaped curve for test accuracy. This is consistent with what we observed on the training loss (Figure 12 in Appendix D.3) for which different choices of ReLU′(0) lead to indistinguishable training loss curves. This experiment suggests that batch normalization has a significant impact in reducing the effect of the choice of ReLU′(0). This observation was confirmed with a very similar behavior on the MNIST dataset with a fully connected network (Figure 9 in Appendix D.2). We could observe a similar effect on CIFAR 10 using a ResNet18 architecture (see Appendix D.5), however in this case the value of ReLU′(0) still has a significative impact on test error, the ResNet18 architecture being much more sensitive.
Using the Adam optimizer: The Adam optimizer [22] is among the most popular algorithms for neural network training; it combines adaptive step-size strategies with momentum. Adaptive stepsize acts as diagonal preconditioners for gradient steps [22, 13] and therefore can be seen as having an effect on the loss landscape of neural network training problems. We trained a VGG11 network on both CIFAR 10 and SVHN using the Adam optimizer. The results are presented in Figure 4. We observe a qualitatively similar behavior as in the experiments of Section 4.2 but a much lower sensitivity to the magnitude of ReLU′(0). In other words, the use of the Adam optimizer mitigates the effect of this choice, both in terms of test errors and by buffering the magnitude of the sometimes chaotic effect induced on training loss optimization (Figure 13 in Appendix D.3).
Increasing numerical precision: As shown in Appendix D.3, using 64 bits floating precision on VGG11 with CIFAR10 cancels out the effect of ReLU′(0) = 1, in coherence with Section 3. More specifically ReLU′(0) = 1 in 64 bits precision obtains similar performances as ReLU′(0) = 0 in 32 bits precision. Furthermore, the numerical precision has barely any effect when ReLU′(0) = 0. Finally, we remark that in 16 bits with ReLU′(0) = 1, training is extremely unstable so that we were not able to train the network in this setting.
Combination with dropout: Dropout [32] is another algorithmic process to regularize deep neural networks. We investigated its effect when combined with varying choices of ReLU′(0). We used a fully connected neural network trained on the MNIST dataset with different values of dropout probability. The results are reported in Figure 11 in Appendix D.2. We did not observe any joint effect of dropout probability and magnitude of ReLU′(0) in this context.

heading: 4.4 Back to the bifurcation zone: more neural nets and the effects of batch-norm
text: The training experiments are complemented by a similar Monte Carlo estimation of the relative volume of the bifurcation zone as performed in Section 3 (same experimental setting). To avoid random outputs we force the GPU to compute convolutions deterministically. Examples of results are given in Figure 5. Additional results on fully connected network with MNIST and ResNet18 on CIFAR10 are shown in Section D. We consistently observe a high probability of finding an example on the bifurcation zone for large sample sizes. Several comments are in order.
Numerical precision: Numerical precision is the major factor in the thickening of the bifurcation zone. In comparison to 32 bits experiments, 16 bits precision dramatically increases its relative importance. We also considered 64 bits precision on SVHN, a rather large dataset. Due to the computational cost, we only drew 40 random weights and observed no bifurcation on any of the terms of the loss whatsoever. This is consistent with the experiments conducted in Section 3 and suggests that, within our framework, 64 bit precision is the main mitigating factor for our observations.
Batch normalization: In all our relative volume estimation experiments, we observe that batch normalization has a very significant effect on the proportion of examples found in the bifurcation zone. In 32 bits precision, the relative size of this zone increases with the addition of batch normalization, similar observations were made in all experiments presented in Appendix D. This is a counter-intuitive behavior as we have observed that batch normalization increases test accuracy and mitigates the effect of ReLU′(0). Similarly in 16 bits precision, the addition of batch normalization seems to actually decrease the size of the bifurcation zone. Batch normalization does not result in the same qualitative effect depending on arithmetic precision. These observations open many questions which will be the topic of future research.
4.5 Total number of ReLU calls during training
We consider the MNIST dataset with a fully connected ReLU network, varying the number of layers in 1, . . . , 6 and neuron per layers in 16, 64, 256, 512. For 16 and 32 bits precisions, we perform 100
epochs and proportion of ReLU(0) over all ReLU calls during training. Figure 6 suggests that, for a given precision, the number of times the bifurcation zone is met is roughly proportional to the total number of ReLU calls during training, independently of the architecture used (number of layers, neurons per layers). This suggests that the total number of ReLU calls during training is an important factor in understanding the bifurcation phenomenon. Broader investigation of this conjecture will be a matter of future work.

heading: 5 Conclusions and future work
text: The starting point of our work was to determine if the choice of the value s = ReLU′(0) affects neural network training. Theory tells that this choice should have negligible effect. Performing a simple learning experiment, we discovered that this is false in the real world and the first purpose of this paper is to account for this empirical evidence. This contradiction between theory and practice is due to finite floating point arithmetic precision while idealized networks are analyzed theoretically within a model of exact arithmetic on the field of real numbers. Owing to the size of deep learning problems, rounding errors due to numerical precision occur at a relatively high frequency, and virtually all the time for large architectures and datasets under 32 bit arithmetic precision (the default choice for TensorFlow and PyTorch libraries).
Our second goal was to investigate the impact of the choice of s = ReLU′(0) in machine learning terms. In 32 bits precision it has an effect on test accuracy which seems to be the result of inducing a chaotic behavior in the course of empirical risk minimization. This was observed consistently in all our experiments. However we could not identify a systematic quantitative description of this effect; it highly depends on the dataset at hand, the network structure as well as other learning parameters such as the presence of batch normalization and the use of different optimizers. Our experiments illustrate this diversity. We observe an interesting robustness property of batch normalization and the Adam optimizer, as well as a high sensitivity to the network structure.
Overall, the goal of this work is to draw attention to an overlooked factor in machine learning and neural networks: nonsmoothness. The ReLU activation is probably the most widely used nonlinearity in this context, yet its nondifferentiability is mostly ignored. We highlight the fact that the default choice ReLU′(0) = 0 seems to be the most robust, while different choices could potentially lead to instabilities. For a general nonsmooth nonlinearity, it is not clear a priori which choice would be the most robust, if any, and our investigation underlines the potential importance of this question. Our research opens new directions regarding the impact of numerical precision on neural network training, its interplay with nonsmoothness and its combined effect with other learning factors, such as network architecture, batch normalization or optimizers. The main idea is that mathematically negligible factors are not necessarily computationally negligible.
Acknowledgments and Disclosure of Funding
The authors thank anonymous referees for constructive suggestions which greatly improved the paper. The authors acknowledge the support of the DEEL project, the AI Interdisciplinary Institute ANITI funding, through the French “Investing for the Future – PIA3” program under the Grant agreement ANR-19-PI3A-0004, Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbers FA9550-19-1-7026, FA9550-18-1-0226, and ANR MaSDOL 19-CE23-0017-01. J. Bolte also acknowledges the support of ANR Chess, grant ANR-17-EURE0010 and TSE-P.

Review summaries:
The paper demonstrates empirically that on a variety of architectures and datasets, the derivative of ReLU at 0 is evaluated a significant number of times during training in 16 and 32 bit precision. These events are frequent enough such that defining ReLU'(0) to be anything other than 0 (default value in common DNN frameworks) hurts generalization, and sometimes training as well.  The authors perform a variety of investigations around this effect, such as the influence of training set size, network width, depth, batch size, BatchNorm, Dropout. 

This paper analyzes the impact of the choice of ReLU'[0] on the fidelity of neural networks. Even though one would expect that the value of ReLU'[0] should not have any noticeable impact, the authors notice the divergence in learned parametric values that depend on this choice. They also notice that this effect is more pronounced at lower precision (due to rounding effects that are dominant at fp16 and fp32 vs. fp64). They also study this effect for neural network sizes, topologies (including networks with and without batch normalization) and datasets.

The ReLU function is not differentiable at 0 which has a mathematical impact on ReLU neural networks. Nevertheless they are optimized using gradient based methods. The folk wisdom is that singular points only occur with zero probability and therefore the issue is irrelevant. This paper challenges this folk wisdom and identifies that the choice of the derivative at 0 does have a tangible impact on the optimization and classification performance.

This paper studies a long-neglected problem: the gradient of the non-smooth activation function in training a neural network. Specifically, the authors study the numerical influence of the ReLU derivative at zero on gradient descent. First, a theory is presented to illustrates the trivialness of ReLU'(0), then, the authors give counter-theory experiments to show that ReLU'(0) can affect training results.  In the last two sections, authors demonstrate that low-precision training, SGD optimizer, no batch norm training can magnify the impact of ReLU'(0).