We have tested LeNet with several sizes (see Figure 9). In Figure 9a, we have plotted the training loss. Overall, ReLU and tanh perform well, along with some $\phi_{\theta}$ with small $\theta$ and $\varphi_{\delta,\omega}$, with $(\delta, \omega) = (0.99, 2)$. In Figure 9b, the results in terms of test accuracy are quite different: ReLU and tanh still achieve good accuracy, but the other functions achieving similar results on the training loss seem to be a bit behind.

So, in this standard setup, the functions we are proposing seem to make the neural network trainable and as expressive as with other activation functions, but with some overfitting. This is not surprising, since we are testing long-standing activation functions, ReLU and tanh, which have been selected both for their ability to make the neural network converge quickly with good generalization, against functions we have designed only according to their ability to propagate information. Therefore, according to these plots, taking into account generalization may be the missing piece of our study.

<figure><img src="image_18.png" /><figcaption>Figure 9: Training curves for LeNet with 3 different numbers of neurons per layer.</figcaption></figure>

**Multilayer perceptron.** We have trained a family of multilayer perceptrons on MNIST. They have a constant width $n_l \in \{3, 10\}$ and a depth $L \in \{3, 10, 30\}$. So, extreme cases such as a narrow and deep neural network ($n_l = 3, L = 30$) have been tested.

Two series of results are presented in Figure 10. The setups of Figure 10a and Figure 10a are identical, except for the initial random seed. In terms of training, the strength of our
Mufan Li, Mihai Nica, and Dan Roy. The neural covariance sde: Shaped infinite depth-and-width networks at initialization. *Advances in Neural Information Processing Systems*, 35:10795–10808, 2022.

Haitao Mao, Xu Chen, Qiang Fu, Lun Du, Shi Han, and Dongmei Zhang. Neuron campaign for initialization guided by information bottleneck theory. In *Proceedings of the 30th ACM International Conference on Information & Knowledge Management*, pp. 3328–3332, 2021.

James Martens, Andy Ballard, Guillaume Desjardins, Grzegorz Swirszcz, Valentin Dalibard, Jascha Sohl-Dickstein, and Samuel S Schoenholz. Rapid training of deep neural networks without skip connections or normalization layers using deep kernel shaping. *arXiv preprint arXiv:2110.01765*, 2021.

Alexander G de G Matthews, Jiri Hron, Mark Rowland, Richard E Turner, and Zoubin Ghahramani. Gaussian process behaviour in wide deep neural networks. In *International Conference on Learning Representations*, 2018.

Radford M Neal. Bayesian learning for neural networks, volume 118 of Lecture Notes in Statistics. Springer, New York, 1996.

Lorenzo Noci, Gregor Bachmann, Kevin Roth, Sebastian Nowozin, and Thomas Hofmann. Precise characterization of the prior predictive distribution of deep ReLU networks. *Advances in Neural Information Processing Systems*, 34:20851–20862, 2021.

Yann Ollivier. Online natural gradient as a Kalman filter. *Electronic Journal of Statistics*, 12(2):2930–2961, 2018.

Theodore Papamarkou, Maria Skoularidou, Konstantina Palla, Laurence Aitchison, Julyan Arbel, David Dunson, Maurizio Filippone, Vincent Fortuin, Philipp Hennig, Aliaksandr Hubin, Alexander Immer, Theofanis Karaletsos, Mohammad Emtiyaz Khan, Agustinus Kristiadi, Yingzhen Li, Jose Miguel Hernandez Lobato, Stephan Mandt, Christopher Nemeth, Michael A. Osborne, Tim G. J. Rudner, David Rügamer, Yee Whye Teh, Max Welling, Andrew Gordon Wilson, and Ruqi Zhang. Position Paper: Bayesian Deep Learning is Needed in the Age of Large-Scale AI. In *International Conference on Machine Learning*, pp. 39556–39586. PMLR, 2024.

Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In *International Conference on Machine Learning*, pp. 1310–1318. PMLR, 2013.

Stefano Peluchetti, Stefano Favaro, and Sandra Fortini. Stable behaviour of infinitely wide deep neural networks. In *International Conference on Artificial Intelligence and Statistics*, pp. 1137–1146. PMLR, 2020.

Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli. Exponential expressivity in deep neural networks through transient chaos. *Advances in Neural Information Processing Systems*, 29, 2016.

Hadi Pouransari, Zhucheng Tu, and Oncel Tuzel. Least squares binary quantization of neural networks. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops*, pp. 698–699, 2020.

Daniel A. Roberts, Sho Yaida, and Boris Hanin. *The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks*. Cambridge University Press, 2022.

Andrew M Saxe, Yamini Bansal, Joel Dapello, Madhu Advani, Artemy Kolchinsky, Brendan D Tracey, and David D Cox. On the information bottleneck theory of deep learning. *Journal of Statistical Mechanics: Theory and Experiment*, 2019(12):124020, 2019.

Samuel S Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein. Deep information propagation. In *International Conference on Learning Representations*, 2017.

Mariia Seleznova and Gitta Kutyniok. Analyzing finite neural networks: Can we trust neural tangent kernel theory? In *Mathematical and Scientific Machine Learning*, pp. 868–895. PMLR, 2022.
activation functions $\phi_{\theta}$ is more visible in the case of narrow neural networks ($n_l = 3$): in general, the loss decreases faster and attains better optima with $\phi = \phi_{\theta}$ than with $\phi = \tanh$ or ReLU.

We also notice that training a narrow and deep neural network with a ReLU activation function is challenging. This result is consistent with several observations we have made in Section 2.2: in narrow ReLU networks, the sequence of correlations $(C_{pq}^l)_{l}$ fails to converge to 1 (Fig. 1), and the pre-activations are far from being Gaussian (Fig. 2).

Finally, $\phi_{2.05}$ put aside, the results regarding the activation functions $\phi_{\theta}$ are consistent between the two runs. This is not the case with tanh and ReLU.

# 5 Discussion

**Generality of Constraints 1 to 4.** Provided that we want pre-activations that are Gaussian $N(0, 1)$, and we want the weights of each layer to be i.i.d. at initialization, the set of four constraints that we provide in Section 2 must hold. If one wants to relax these constraints, it becomes unavoidable to break symmetries or to solve harder problems. Let us consider two examples.

**Example 1:** instead of making all the pre-activations Gaussian, one may want to impose some other distribution. In this case, the problem to solve would be more difficult: we would have to decompose a non-Gaussian random variable $Z$ into a weighted sum $\frac{1}{\sqrt{n}} \sum_{j=1}^n W_j \phi(X_j)$. In this case, we cannot use Proposition 3, and we would have to make $W_j \phi(X_j)$ belong to a family of random variables stable by multiplication by a constant and by sum, for arbitrary $n$. If we aim for a non-Gaussian $Z$, this task is much harder. For a study of infinitely wide neural networks going beyond Gaussian pre-activations, see Peluchetti et al. (2020).

**Example 2:** instead of assuming that all the weights of a given layer are i.i.d., one may want to initialize them with different distributions or to introduce a dependence structure between them. If the goal remains to obtain Gaussian pre-activations, this kind of generalization should be feasible without changing drastically the constraints we are proposing.

**Hypothesis of independent pre-activations.** In the constraints, we have assumed that, for any layer, its inputs are independent. This assumption is discussed in Remark 7 and in Appendix B. According to the experimental results presented in Appendix B, we can build specific cases where the dependence between inputs breaks the Gaussianity of the outputted pre-activations, including with our pairs $(\phi_{\theta}, P_{\theta})$.

Also, when comparing the distribution of the pre-activations obtained with $(\phi_{\theta}, P_{\theta})$ in the case of independent inputs (Fig. 6b) and in the case of dependent inputs (Fig. 7c and 11), it is probable that the dependence between pre-activations tends to damage the Gaussianity after a certain number of layers, for pairs $(\phi_{\theta}, P_{\theta})$ with large $\theta$.
Ravid Shwartz-Ziv and Naftali Tishby. Opening the black box of deep neural networks via information. *arXiv preprint arXiv:1703.00810*, 2017.

Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Implicit neural representations with periodic activation functions. *Advances in Neural Information Processing Systems*, 33:7462–7473, 2020.

Nickolay Smirnov. Table for estimating the goodness of fit of empirical distributions. *The Annals of Mathematical Statistics*, 19(2):279–281, 1948.

P S Theocaris and A C Chrysakis. Numerical inversion of the Mellin transform. *IMA Journal of Applied Mathematics*, 20(1):73–83, 1977.

Naftali Tishby. The information bottleneck method. In *Proc. 37th Annual Allerton Conference on Communications, Control and Computing, 1999*, pp. 368–377, 1999.

Edward Charles Titchmarsh. *Introduction to the theory of Fourier integrals*. The Clarendon Press, Oxford, 2nd edition, 1937.

Dario Trevisan. Wide deep neural networks with Gaussian weights are very close to Gaussian processes. *arXiv preprint arXiv:2312.11737*, 2023.

Mariia Vladimirova, Jakob Verbeek, Pablo Mesejo, and Julyan Arbel. Understanding priors in Bayesian neural networks at the unit level. In *International Conference on Machine Learning*, pp. 6458–6467. PMLR, 2019.

Mariia Vladimirova, Julyan Arbel, and Stéphane Girard. Bayesian neural network unit priors and generalized Weibull-tail property. In *Asian Conference on Machine Learning*, pp. 1397–1412. PMLR, 2021.

Florian Wenzel, Kevin Roth, Bastiaan Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How good is the Bayes posterior in deep neural networks really? In *International Conference on Machine Learning*, pp. 10248–10259. PMLR, 2020.

Sho Yaida. Non-gaussian processes and neural networks at finite widths. In *Mathematical and Scientific Machine Learning*, pp. 165–192. PMLR, 2020.

Guodong Zhang, Aleksandar Botev, and James Martens. Deep learning without shortcuts: Shaping the kernel with tailored rectifiers. In *International Conference on Learning Representations*, 2022.
<figure><img src="image_19.png" /><figcaption>Figure 10: Training loss for a multilayer perceptron, narrow ($n_l \in \{3, 10\}$) and of various depths ($L \in \{3, 10, 30\}$).</figcaption></figure>
# A Activation function with infinite number of stable fixed points for V

## A.1 Proof that V admits an infinite number of fixed points when using $\phi = \varphi_{\delta,\omega}$

**Proposition 2.** For any $\delta \in (0, 1]$ and $\omega > 0$, let us pose the activation function $\phi = \varphi_{\delta,\omega}$. We consider the sequence $(v^l)_l$ defined by:

$$ \begin{aligned} \forall l \ge 0, \quad v^{l+1} &= \sigma_w^2 \int \varphi_{\delta,\omega} (\sqrt{v^l z})^2 \mathcal{D}z + \sigma_b^2, \\ v^0 &\in \mathbb{R}_+^*. \end{aligned} \quad (52) $$

Then there exists $\sigma_w > 0$, $\sigma_b \ge 0$, and a strictly increasing sequence of stable fixed points $(v_k^*)_{k \in \mathbb{Z}}$ of the recurrence Equation (52).

*Proof.* Let us define:

$$ \begin{aligned} \tilde{\mathcal{V}}(v) &:= \frac{1}{v} \mathcal{V}(v|\sigma_w = 1, \sigma_b = 0), \\ \text{so we have:} \quad \mathcal{V}(v|\sigma_w, \sigma_b) &= \sigma_w^2 v \tilde{\mathcal{V}}(v) + \sigma_b^2. \end{aligned} \quad (53) $$

In the following, we use a simplified notation: $\mathcal{V}(v) = \mathcal{V}(v|\sigma_w, \sigma_b)$.

Our goal is to find a sequence $(v_k^*)_{k \in \mathbb{Z}}$, $\sigma_w$ and $\sigma_b$ such that:

$$ \forall k \in \mathbb{Z}, \quad \mathcal{V}(v_k^*) = v_k^* \quad \text{and} \quad \mathcal{V}'(v_k^*) \in (-1, 1), $$

which would ensure that all $v_k^*$ are stable fixed points of $\mathcal{V}$. In order to understand how to build the sequence $(v_k^*)_{k \in \mathbb{Z}}$, let us consider $v > 0$. Let $\sigma_b^2 = 0$ and $\sigma_w^2 = 1/\tilde{\mathcal{V}}(v)$. So we have:

$$ \mathcal{V}(v) = v. $$

So, any $v > 0$ can possibly be a fixed point if we tune $\sigma_w$ accordingly. We just have to find a $v > 0$ such that: $\mathcal{V}'(v) \in (-1, 1)$, with:

$$ \mathcal{V}'(v) = \frac{1}{\tilde{\mathcal{V}}(v)} (\tilde{\mathcal{V}}(v) + v \tilde{\mathcal{V}}'(v)) = 1 + \tilde{\mathcal{V}}_e(\ln(v)), \quad (54) $$

where $\tilde{\mathcal{V}}_e: r \mapsto \ln(\tilde{\mathcal{V}}(\exp(r)))$. So, knowing that $\tilde{\mathcal{V}}_e$ is $\mathcal{C}^1$ and periodic, it is sufficient to prove that it is not constant to ensure that we can extract one $v_0^*$ such that $\mathcal{V}'(v_0^*) \in (-1, 1)$. Then, by periodicity of $\tilde{\mathcal{V}}_e$, we can build a sequence of stable fixed points $(v_k^*)_{k \in \mathbb{Z}}$.

We have:

$$ \begin{aligned} \tilde{\mathcal{V}}(v) &= \frac{1}{v} \int_{-\infty}^{\infty} \varphi_{\delta,\omega} (\sqrt{v}z)^2 \mathcal{D}z = \int_{-\infty}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln |\sqrt{v}z|) \right) \mathcal{D}z \\ &= 2 \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln(\sqrt{v}z)) \right) \mathcal{D}z = 2 \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin \left( \frac{\omega}{2} \ln v + \omega \ln z \right) \right) \mathcal{D}z. \end{aligned} $$

**Lemma 1.** $\tilde{\mathcal{V}}$ is not constant.

*Proof.*

$$ \begin{aligned} \tilde{\mathcal{V}}(v) &= \frac{1}{\sqrt{2\pi}} \frac{2}{v^{3/2}} \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln z) \right) \exp \left( -\frac{z^2}{2v} \right) \mathcal{D}z \\ &= \frac{1}{\sqrt{2\pi}} \frac{2}{v^{3/2}} \int_{0}^{\infty} \frac{z}{2\sqrt{z}} \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln \sqrt{z}) \right) \exp \left( -\frac{z}{2v} \right) \mathcal{D}z \\ &= \frac{1}{\sqrt{2\pi}} v^{-3/2} \mathcal{L}[\varphi_{\delta,\omega}(\sqrt{z})] \left( \frac{1}{2v} \right), \end{aligned} $$
**Other families of initialization distributions and activation functions.** Provided the constraints we have derived, we have made one choice to obtain our family of initialization distributions and activation functions: we have decided that the weights should be sampled from a symmetric Weibull distribution $W(\theta, 1)$. We have made this choice because Weibull distributions meet immediately Constraint 2, and we can modulate their Generalized Weibull Tail parameter easily (see Remark 9).

One may propose another family of initialization distributions, as long as it meets the constraints. However, such a family should be selected wisely: if one chooses a distribution with compact support, or GWT($\theta$) with $\theta \gg 2$, then the related activation function $\phi$ is likely to be almost linear. Intuitively, if we want $W\phi(X)$ to be $N(0, 1)$ with $X \sim N(0, 1)$ and very light-tailed $W$, then $\phi$ must “reproduce” the tail of its input $X$, so $\phi$ must be approximately linear around infinity.

**Preserve a characteristic during propagation or impose stable fixed points?** In previous works, both ideas have been used: Glorot & Bengio (2010) wanted to preserve the variance of the pre-activations, while Poole et al. (2016); Schoenholz et al. (2017) wanted to impose a specific stable fixed point for the correlation map $C$. According to the results we have obtained in Section 4.1, it is possible to preserve *approximately* the distribution of the pre-activations when passing through *one* layer. However, according to the results of Section 4.2, a drift can appear after several layers. So, when testing an initialization setup in the real world, with numerical errors and approximations, it is necessary to check the stable fixed points of the monitored characteristic. However, this is not easy to do in practice: without the Gaussian hypothesis, finding the possible limits of $(c_{ab}^l)_l$ can be difficult.

**Gaussian pre-activations and Neural Tangent Kernels (Jacot et al., 2018).** With our pair of activation functions and initialization procedure, we have been able to provide an exact Edge of Chaos, removing the infinite-width assumption. Since this infinite-width assumption is also fundamental in works on the Neural Tangent Kernels (NTKs), would it be possible to obtain the same kind of results within the NTK setup? We believe it is not. On one side, the infinite-width limit is used to end up with an NTK, which is *constant* during training, in order to provide exact equations of evolution of the trained neural network. On the other side, within our setup, we only ensure Gaussian pre-activations, and it does not imply that the NTK would be constant during training. Nevertheless, with Gaussian pre-activations, we might expect an improvement in the convergence rate of the neural network towards a stacked Gaussian process, as the widths of the layers tend to infinity. So, the NTK regime would be easier to attain.

**Taking into account the generalization performance.** In the EOC framework and ours, a common principle could be discussed and improved. The generalization performance is not taken into account at any step of the reasoning. As we have seen in the training experiments, the main difficulty encountered within our setup is the overfitting: as the training loss decreases without obstacle, the test loss is, at the end, worse than with ReLU
where $\mathcal{L}$ is the Laplace transform.

Let us suppose that $\tilde{\mathcal{V}}$ is constant: $\forall v > 0, \tilde{\mathcal{V}}(v) = c$. So, if we pose $v \leftarrow \frac{1}{2v}$, we have:

$$ \forall v > 0, \quad c = \frac{2}{\sqrt{\pi}} v^{3/2} \mathcal{L}[\varphi_{\delta,\omega}(\sqrt{z})](v), $$

that is:

$$ \mathcal{L}[\varphi_{\delta,\omega}(\sqrt{z})](v) = \frac{c\sqrt{\pi}}{2}v^{-3/2}. $$

Since the function $z \mapsto \varphi_{\delta,\omega}(\sqrt{z})$ is continuous on $\mathbb{R}^+$, then, almost everywhere (see Thm. 22.2, Billingsley, 1995):

$$ \varphi_{\delta,\omega}(\sqrt{z}) = \frac{c\sqrt{\pi}}{2}\mathcal{L}^{-1}[v^{-3/2}](z) = c\sqrt{z}, $$

which is impossible for $\delta \in (0, 1]$. Hence the result. $\square$

The function $\tilde{\mathcal{V}}_e : r \mapsto \ln(\tilde{\mathcal{V}}(\exp(r)))$ is continuous and $\frac{4\pi}{\omega}$-periodic:

$$ \tilde{\mathcal{V}}_e(r) = \ln \left[ 2 \int_0^\infty z^2 \exp \left( 2 \frac{\delta}{\omega} \sin \left( \frac{\omega}{2} r + \omega \ln z \right) \right) \mathcal{D}z \right], $$

so $\tilde{\mathcal{V}}_e$ is lower and upper bounded and reach its bounds (and, by Lemma 1, these bounds are different). We define:

$$ \tilde{\mathcal{V}}_e^+ = \max \tilde{\mathcal{V}}_e $$

$$ r_0^+ = \inf \{ r > 0 : \tilde{\mathcal{V}}_e(r) = \tilde{\mathcal{V}}_e^+ \}, $$

$$ \tilde{\mathcal{V}}_e^- = \min \tilde{\mathcal{V}}_e $$

$$ r_0^- = \inf \{ r > r_0^+ : \tilde{\mathcal{V}}_e(r) = \tilde{\mathcal{V}}_e^- \}. $$

By continuity, $\tilde{\mathcal{V}}_e(r_0^+) = \tilde{\mathcal{V}}_e^+$ and $\tilde{\mathcal{V}}_e(r_0^-) = \tilde{\mathcal{V}}_e^-$. Since $\tilde{\mathcal{V}}^+ > \tilde{\mathcal{V}}^-$ and $\tilde{\mathcal{V}}_e$ is $\mathcal{C}^1$, then there exists $r_0^* \in (r_0^+, r_0^-)$ such that $\tilde{\mathcal{V}}'_e(r_0^*) \in (-2, 0)$.

Since $\tilde{\mathcal{V}}_e$ is $\frac{4\pi}{\omega}$-periodic, we can define a sequence $(r_k^*)_{k \in \mathbb{Z}}$ such that:

$$ r_k^* := r_0^* + \frac{4k\pi}{\omega} $$

$$ \tilde{\mathcal{V}}_e(r_k^*) = \tilde{\mathcal{V}}_e(r_0^*) := \tilde{\mathcal{V}}_e^0 $$

$$ \tilde{\mathcal{V}}'_e(r_k^*) = \tilde{\mathcal{V}}'_e(r_0^*) \in (-2, 0). $$

So, by using Eqn. (53) and Eqn. (54) with $\sigma_b^2 = 0$ and $\sigma_w^2 = 1/\exp(\tilde{\mathcal{V}}_e^0)$:

$$ v_k^* := \exp(r_k), $$

$$ \mathcal{V}(v_k^*) = \frac{1}{\exp(\tilde{\mathcal{V}}_e^0)} v_k^* \tilde{\mathcal{V}}(v_k^*) = v_k^* $$

$$ \mathcal{V}'(v_k^*) = 1 + \tilde{\mathcal{V}}'_e(r_0^*) \in (-1, 1). $$

Thus, $(v_k^*)_{k \in \mathbb{Z}}$ is a sequence of stable fixed points of $\mathcal{V}(\cdot|\sigma_w, \sigma_b)$ for well-chosen $\sigma_w$ and $\sigma_b$. $\square$

## A.2 Practical computation of $\sigma_w^2$

We propose a practical method to ensure that $\mathcal{V}(\cdot|\sigma_w, \sigma_b = 0)$ admits an infinite number of stable fixed points when using activation function $\phi = \varphi_{\delta,\omega}$.
and tanh activation functions. To improve these generalization results, one may integrate
into the framework a separation between the training set and a validation set.

**Precise characterization of the pre-activations distributions $\mathcal{D}^l$.** Finally, finding a precise and usable characterization of the $\mathcal{D}^l$ remains an unsolved problem. In the EOC line of work, the problem has been simplified by using the Gaussian hypothesis. But, as shown in the present work, such a simplification is too rough and leads to inconsistent results. Nevertheless, our approach has its own drawbacks. Namely, we can ensure Gaussian pre-activations only when using specific initialization distributions and activation functions. So, we still miss a characterization of the distributions $\mathcal{D}^l$ which would apply to widely-used networks without harsh approximations, and be easy to use to achieve practical goals, such as finding an optimal initialization scheme. From this perspective, we hope that the problem representation of Section 2.4 will result in fruitful future research.

**Is it desirable to have Gaussian pre-activations?** We have provided several results that should help to answer this tough question, but it remains difficult to answer it definitively. One shall note a paradox regarding the ReLU activation. On one side, when ReLU is used in narrow and deep perceptrons, the pre-activations are far from being Gaussian, the sequence of correlations does not converge to 1, and training is difficult and unstable. On the other side, in LeNet-type networks with ReLU, training is easy, and the resulting networks generalize well. Besides, our activation functions $\phi_\theta$ perform quite differently depending on the setup: with $\phi_{2.05}$, LeNet can achieve good training losses, but the training of narrow and deep perceptrons may fail. We observe opposite results with $\phi_{10.00}$. Therefore, the strongest answer we can give is: with Gaussian pre-activations at initialization, a neural network is likely (but not sure) to be trainable, but it is impossible to predict its ability to generalize.

**Is the “Gaussian pre-activations hypothesis” a myth or a reality?** We have shown that, when using the ReLU activation function in several practical cases, it is largely a myth. But, when using tanh, the results depend on the neural network width: with a sufficient, but still reasonable, number of neurons per layer, it becomes a reality. In order to ensure that this reality remains tangible for any number of neurons per layer, we have established a set of constraints that the design of the neural networks must fulfill, and we have proposed a set of solutions fulfilling them. As a result, several of these solutions have all the prerequisite to become strong foundations of this Gaussian hypothesis, making it real in all tested cases.

## Acknowledgements

The project leading to this work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 834175) and from the French National Research Agency (ANR-21-JSTM-0001) in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02). We thank Thomas Dupic for proposing the trick involving the Laplace transform used in Appendix A.1.
In order to achieve this goal, we build $\sigma_w^2 = \sigma_\omega^2$ in the following way:

$$ 
\begin{align*}
\mathcal{V}_{\text{low}} &:= 2 \int_{0}^{\infty} z^2 \exp \left( -2 \frac{\delta}{\omega} \sin(\omega \ln(z)) \right) \mathcal{D}z \\
\mathcal{V}_{\text{upp}} &:= 2 \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln(z)) \right) \mathcal{D}z, \\
\sigma_{\omega}^2 &:= \left[ \frac{\mathcal{V}_{\text{low}} + \mathcal{V}_{\text{upp}}}{2} \right]^{-1}.
\end{align*}
$$

In practice, for $\delta = 0.99$, we obtain $\sigma_\omega$ for $\omega \in \{2, 3, 6\}$:

$$ \sigma_2 \approx 0.879, \quad \sigma_3 \approx 0.945, \quad \sigma_6 \approx 0.987. $$

## B Discussion about the independence of the pre-activations

In Proposition 3 and Constraint 1, we assume that, for any layer $l$, its inputs $(Z_j^l)_j$ are independent. This does not hold in full generality:

**Example 1.** Let $X \sim \mathcal{N}(0, 1)$ be some random input of a two-layer neural network. We perform the following operation:

$$ Z = \frac{1}{\sqrt{2}} [W_1^2 \phi(W_1^1 X) + W_2^2 \phi(W_2^1 X)], $$

where $(W_1^1, W_2^1, W_1^2, W_2^2)$ be i.i.d. random variables samples from some distribution P, and $\phi$ is some activation function.

Let P = $\mathcal{R}$, the Rademacher distribution, i.e., if $W \sim \mathcal{R}$, then $W = \pm 1$ with probability 1/2. Let $\phi = \text{Id}$. Then we have: $Y_1 := W_1^1 X \sim \mathcal{N}(0, 1)$ and $Y_2 := W_2^1 X \sim \mathcal{N}(0, 1)$. But they are not independent:

$$ W_1^2 W_1^1 X + W_2^2 W_2^1 X = (W'_1 + W'_2)X, $$

where $W_1^2 W_1^1$ and $W_2^2 W_2^1$ are two independent Rademacher random variables. So, $Z = 0$ with probability 1/2. So, $Z$ is not Gaussian.

In this example, we build a non-Gaussian random variable with a minimal neural network, in which we construct two dependent random variables. So, we should pay attention to this phenomenon when propagating the pre-activations in a neural network.

One should note that the structure of dependence of $W_1^2 W_1^1 X$ and $W_2^2 W_2^1 X$ does not involve their correlation (which is zero), and yet breaks the Gaussianity of their sum. So, in order to obtain a theoretical result about the distribution of $Z$, we should study finer aspects of the dependence structure.

But, on the practical side, we want to answer the question: to which extent does the relation of dependence between the pre-activations $(Z_j^l)_j$ affect the Gaussianity of $(Z_i^{l+1})_i$?

In order to answer this question, we propose a series of experiments on a two-layer neural network. We consider a vector of inputs $\mathbf{X} \in \mathbb{R}^{n_0}$, where the $(X_j)_{1 \le j \le n_0}$ are $\mathcal{N}(0, 1)$ and i.i.d. The scalar outputted by the network is:

$$ Z = \frac{1}{\sqrt{n_1}} \mathbf{W}^2 \phi \left( \frac{1}{\sqrt{n_0}} \mathbf{W}^1 \mathbf{X} \right), $$

where the weights $\mathbf{W}^1 \in \mathbb{R}^{n_1 \times n_0}$ and $\mathbf{W}^2 \in \mathbb{R}^{1 \times n_1}$ are i.i.d.

We test this setup with different initialization distributions P and activation functions $\phi$:

* usual ones: $\phi = \tanh$ or ReLU, $P = \mathcal{N}(0, \sigma_w^2)$, where $\sigma_w^2$ is such that the pair $(\sigma_w^2, \sigma_b^2 = 0.01)$ lies at the EOC;
We would like to thank an anonymous reviewer for pointing out an error in a previous version of the paper, and proposing an example which inspired Example 1 (see App. B).

## References

Sanjeev Arora, Nadav Cohen, Wei Hu, and Yuping Luo. Implicit regularization in deep matrix factorization. *Advances in Neural Information Processing Systems*, 32, 2019.

Patrick Billingsley. *Probability and Measure*. John Wiley & Sons, third edition, 1995.

Bruno De Finetti. La prévision: ses lois logiques, ses sources subjectives. In *Annales de l'Institut Henri Poincaré*, volume 7, pp. 1-68, 1937.

Benjamin Epstein. Some applications of the Mellin transform in statistics. *The Annals of Mathematical Statistics*, pp. 370-379, 1948.

Arthur Erdélyi, (Hans Heinrich) Wilhelm Magnus, Fritz Oberhettinger, Francesco Giacomo Tricomi, David Bertin, Watson B. Fulks, A. R. Harvey, Donald L. Thomsen, Jr., Maria A. Weber, E. L. Whitney, and Rosemarie Stampfel. *Higher transcendental functions*, volume I. Bateman Manuscript Project, 1953.

Vincent Fortuin, Adrià Garriga-Alonso, Sebastian W. Ober, Florian Wenzel, Gunnar Ratsch, Richard E Turner, Mark van der Wilk, and Laurence Aitchison. Bayesian neural network priors revisited. In *International Conference on Learning Representations*, 2022.

Bruno Gabutti and Laura Sacripante. Numerical inversion of the Mellin transform by accelerated series of Laguerre polynomials. *Journal of Computational and Applied Mathematics*, 34(2):191-200, 1991.

Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In *Proceedings of the 13th International Conference on Artificial Intelligence and Statistics*, pp. 249-256. JMLR Workshop and Conference Proceedings, 2010.

Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik. *Table of Integrals, Series, and Products*. Academic Press, eighth edition, 2014.

Alex Graves. Practical variational inference for neural networks. *Advances in Neural Information Processing Systems*, 24, 2011.

Soufiane Hayou, Arnaud Doucet, and Judith Rousseau. On the impact of the activation function on deep neural networks training. In *International Conference on Machine Learning*, pp. 2672-2680. PMLR, 2019.

Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In *Proceedings of the IEEE international conference on computer vision*, pp. 1026-1034, 2015.
• ours: $\phi = \phi_{\theta}$, $P = P_{\theta} = W(\theta, 1)$.

We study three cases:

• unfavorable case: $n_0 = 1$, various $n_1 \in [1, 10]$; the intermediary features $\mathbf{W}^1\mathbf{X}$ are weakly “mixed”, so it is credible that they lead to an output that is far from being Gaussian;

• favorable case: $n_1 = 2$, various $n_0 \in [1, 10]$;<fnref n="21" /> the intermediary features $\mathbf{W}^1\mathbf{X}$ are “mixed” with an increasing rate as $n_0$ increases;

• same-width case: $n_0 = n_1$.

<figure><img src="image_17.png" /><figcaption>Figure 11: Evolution of the distance of the standardized distribution of $Z$ to the $N(0, 1)$ according to the Kolmogorov-Smirnov statistic. Weight initialization is $W(\theta, 1)$ when using $\phi = \phi_{\theta}$ and is Gaussian according to the EOC when using $\phi = \tanh$ or ReLU. For each point, we have computed the KS-statistic over 200 000 samples of $(\mathbf{X}, \mathbf{W}^1, \mathbf{W}^2)$. In the first row, we have plotted the results for the **odd** activation functions $\phi_{\theta}^o$. In the second row, we have plotted the results for the **positive** activation functions $\phi_{\theta}^p$.</figcaption></figure>

According to Figure 11:

• unfavorable case (1st graph, $n_0 = 1$): $Z$ is far from being Gaussian for ReLU, tanh, and the $\phi_{\theta}^o$, but is much better for the $\phi_{\theta}^p$, especially with a small $\theta$;

• favorable case (2nd graph, $n_1 = 2$): $Z$ is closer to be Gaussian with our setups ($P_{\theta}, \phi_{\theta}^{o/p}$) than with $\phi = \tanh$ or ReLU, especially with larger $n_1$;

• same-width case (3rd graph, $n = n_0 = n_1$): the larger the width $n$, the closer $Z$ is to being Gaussian. For a fixed $n$, the distribution of $Z$ is closer to a Gaussian with smaller $\theta$. When we use ($P_{\theta}, \phi_{\theta}^{o/p}$), we are close to the performance of tanh or better.

<footnote n="21">²¹The case $n_1 = 1$ is trivial: there is no sum of dependent random variables in the second layer.</footnote>
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. *Journal of Machine Learning Research*, 2013.

Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. *Advances in Neural Information Processing Systems*, 31, 2018.

Diederik P Kingma and Jimmy L Ba. Adam: A method for stochastic optimization. In *International Conference on Learning Representations*, 2015.

Andrey Kolmogoroff. Confidence limits for an unknown distribution function. *The Annals of Mathematical Statistics*, 12(4):461–463, 1941.

Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324, 1998.

Haitao Mao, Xu Chen, Qiang Fu, Lun Du, Shi Han, and Dongmei Zhang. Neuron campaign for initialization guided by information bottleneck theory. In *Proceedings of the 30th ACM International Conference on Information & Knowledge Management*, pp. 3328–3332, 2021.

Alexander G de G Matthews, Jiri Hron, Mark Rowland, Richard E Turner, and Zoubin Ghahramani. Gaussian process behaviour in wide deep neural networks. In *International Conference on Learning Representations*, 2018.

Lorenzo Noci, Gregor Bachmann, Kevin Roth, Sebastian Nowozin, and Thomas Hofmann. Precise characterization of the prior predictive distribution of deep ReLU networks. *Advances in Neural Information Processing Systems*, 34:20851–20862, 2021.

Yann Ollivier. Online natural gradient as a Kalman filter. *Electronic Journal of Statistics*, 12(2):2930–2961, 2018.

Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In *International Conference on Machine Learning*, pp. 1310–1318. PMLR, 2013.

Stefano Peluchetti, Stefano Favaro, and Sandra Fortini. Stable behaviour of infinitely wide deep neural networks. In *International Conference on Artificial Intelligence and Statistics*, pp. 1137–1146. PMLR, 2020.

Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli. Exponential expressivity in deep neural networks through transient chaos. *Advances in Neural Information Processing Systems*, 29, 2016.

Hadi Pouransari, Zhucheng Tu, and Oncel Tuzel. Least squares binary quantization of neural networks. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops*, pp. 698–699, 2020.
The case $n_0 = n_1$ is the most realistic one: usually, the sizes of the layers of a neural network are of the same order of magnitude. In this case, our setup is better than or equivalent to the ones with tanh or ReLU.

**Mixing of inputs and interference phenomenon.** We see in Figure 11 that, in every graph, setups with small $\theta$ lead to better results than the ones with large $\theta$. This observation could be explained by Example 1. In this example, the inputs are weakly “mixed”: since $X_0$ is multiplied by a Rademacher random variable, it is possible to partially reconstruct $X_0$ after the first layer, and then build destructive interference (leading to an output $Z = 0$ half of the time).

So, if we want to avoid this “interference” phenomenon, we should use initialization distributions and activation functions such that every layer “mixes” strongly the inputs. Notably, initialization distributions should be far from being a combination of Dirac distributions. Typically, $P_{\theta=10}$ is close to the Rademacher distribution (see Fig. 5a).

Finally, this phenomenon can be greatly reduced by using our family $(\phi_{\theta}^P)_{\theta}$ of positive activation functions (see second row of Figure 11).

## C Constraints on the product of two random variables

**Proposition 4 (Density of a product of random variables at 0).** Let W, Y be two independent non-negative random variables and $Z = WY$. Let $f_W, f_Y, f_Z$ be their respective density. Assuming that $f_Y$ is continuous at 0 with $f_Y(0) > 0$, we have:

$$ \text{if } \lim_{w \to 0} \int_w^\infty \frac{f_W(t)}{t} dt = \infty, \quad \text{then} \quad \lim_{z \to 0} f_Z(z) = \infty. \qquad (55) $$

Moreover, if $f_Y$ is bounded:

$$ \text{if } \int_0^\infty \frac{f_W(t)}{t} dt < \infty, \quad \text{then} \quad f_Z(0) = f_Y(0) \int_0^\infty \frac{f_W(t)}{t} dt. \qquad (56) $$

*Proof.* Let $z, z_0 > 0$:

$$ \begin{aligned} f_Z(z) &= \int_0^\infty f_Y(t) \frac{1}{t} f_W \left( \frac{z}{t} \right) dt \\ &\ge \int_0^{z_0} f_Y(t) \frac{1}{t} f_W \left( \frac{z}{t} \right) dt \ge \inf_{[0, z_0]} f_Y \cdot \int_0^{z_0} \frac{1}{t} f_W \left( \frac{z}{t} \right) dt \ge \inf_{[0, z_0]} f_Y \cdot \int_{z/z_0}^\infty \frac{f_W(t)}{t} dt. \end{aligned} $$

Let us take $z_0 = \sqrt{z}$. We have:

$$ f_Z(z) \ge \inf_{[0, \sqrt{z})} f_Y \cdot \int_{\sqrt{z}}^\infty \frac{f_W(t)}{t} dt. $$

Then we take the limit $z \to 0$, hence:

* if $\int_0^\infty \frac{f_W(t)}{t} dt = \infty$, then: $\lim_{z \to 0} f_Z(z) = \infty$, which achieves (55);

* if $\int_0^\infty \frac{f_W(t)}{t} dt < \infty$, then: $f_Z(0) \ge f_Y(0) \int_0^\infty \frac{f_W(t)}{t} dt$, which achieves one half of (56);

Let us prove the second half of (56). Let $z, z_0 > 0$:

$$ \begin{aligned} f_Z(z) &= \int_0^\infty f_Y(t) \frac{1}{t} f_W \left( \frac{z}{t} \right) dt \\ &= \int_0^{z_0} f_Y(t) \frac{1}{t} f_W \left( \frac{z}{t} \right) dt + \int_{z_0}^\infty f_Y(t) \frac{1}{t} f_W \left( \frac{z}{t} \right) dt \\ &\le \sup_{[0, z_0]} f_Y \cdot \int_{z/z_0}^\infty \frac{f_W(t)}{t} dt + \int_1^\infty f_Y(z_0 t) \frac{1}{t} f_W \left( \frac{z}{z_0 t} \right) dt \end{aligned} $$
Andrew M Saxe, Yamini Bansal, Joel Dapello, Madhu Advani, Artemy Kolchinsky, Brendan D Tracey, and David D Cox. On the information bottleneck theory of deep learning. *Journal of Statistical Mechanics: Theory and Experiment*, 2019(12):124020, 2019.

Samuel S Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein. Deep information propagation. In *International Conference on Learning Representations*, 2017.

Ravid Shwartz-Ziv and Naftali Tishby. Opening the black box of deep neural networks via information. *arXiv preprint arXiv:1703.00810*, 2017. URL https://arxiv.org/pdf/1703.00810.pdf.

Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Implicit neural representations with periodic activation functions. *Advances in Neural Information Processing Systems*, 33:7462–7473, 2020.

Nickolay Smirnov. Table for estimating the goodness of fit of empirical distributions. *The Annals of Mathematical Statistics*, 19(2):279–281, 1948.

P S Theocaris and A C Chrysakis. Numerical inversion of the Mellin transform. *IMA Journal of Applied Mathematics*, 20(1):73–83, 1977.

Naftali Tishby. The information bottleneck method. In *Proc. 37th Annual Allerton Conference on Communications, Control and Computing*, 1999, pp. 368–377, 1999.

Edward Charles Titchmarsh. *Introduction to the theory of Fourier integrals*. The Clarendon Press, Oxford, 2nd edition, 1937.

Mariia Vladimirova, Jakob Verbeek, Pablo Mesejo, and Julyan Arbel. Understanding priors in Bayesian neural networks at the unit level. In *International Conference on Machine Learning*, pp. 6458–6467. PMLR, 2019.

Mariia Vladimirova, Julyan Arbel, and Stéphane Girard. Bayesian neural network unit priors and generalized Weibull-tail property. In *Asian Conference on Machine Learning*, pp. 1397–1412. PMLR, 2021.

Florian Wenzel, Kevin Roth, Bastiaan Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How good is the Bayes posterior in deep neural networks really? In *International Conference on Machine Learning*, pp. 10248–10259. PMLR, 2020.
Let $z_0 = \sqrt{z}$. We have:

$$f_Z(z) \leq \sup_{[0, \sqrt{z}]} f_Y \cdot \int_{\sqrt{z}}^{\infty} \frac{f_W(t)}{t} dt + \int_{1}^{\infty} f_Y(\sqrt{zt}) \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt,$$

where:

$$\int_{1}^{\infty} f_Y(\sqrt{zt}) \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt \leq \|f_Y\|_{\infty} \int_{1}^{\infty} \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt \leq \|f_Y\|_{\infty} \int_{0}^{\sqrt{z}} \frac{f_W(t)}{t} dt.$$

According to the hypotheses, we have, as $z \to 0$:

$$\sup_{[0, \sqrt{z}]} f_Y \to f_Y(0), \qquad \int_{\sqrt{z}}^{\infty} \frac{f_W(t)}{t} dt \to \int_{0}^{\infty} \frac{f_W(t)}{t} dt, \qquad \int_{0}^{\sqrt{z}} \frac{f_W(t)}{t} dt \to 0,$$

hence the result.

□

# D Activation functions with vertical tangent at 0

In the following lemma, we show that if we want the activation $Y$ to have a density that is 0 at 0, then the activation function $\phi$ should have a vertical tangent at 0. $G$ plays the role of pre-activation.

**Lemma 2.** Let $\phi$ be a function transforming a Gaussian random variable $G \sim N(0, 1)$ into a symmetrical random variable $Y$ with a density $f_Y$ such that $f_Y(0) = 0$. That is, $Y = \phi(G)$. Then $\phi$ has a vertical tangent at 0.

*Proof.* We have:

$$\phi(x) = F_Y^{-1}(F_G(x)),$$

where $F_G$ and $F_Y$ are the respective CDFs of $G$ and $Y$.

Thus:

$$\phi'(x) = F'_{G}(x) \frac{1}{F'_{Y}(F^{-1}_{Y}(F_{G}(x)))}$$

Therefore:

$$\phi'(0) = F'_{G}(0) \frac{1}{F'_{Y}(F^{-1}_{Y}(F_{G}(0)))} = \frac{1}{\sqrt{2\pi}} \frac{1}{F'_{Y}(F^{-1}_{Y}(1/2))} = \frac{1}{\sqrt{2\pi}} \frac{1}{F'_{Y}(0)} = \infty.$$

□

# E The Mellin transform

## E.1 Generalities

We assume that $G = WY \sim N(0, 1)$. Let us consider the random variables $|W|$, $|Y|$ and $|G| = |W| \cdot |Y|$. Let $f_{|W|}$, $f_{|Y|}$ and $f_{|G|}$ be their densities. Under integrability conditions, we can express the density $f_{|G|}$ of the product $|G| = |W||Y|$ with the product-convolution operator $*$:

$$f_{|G|}(z) = (f_{|W|} * f_{|Y|})(z), \quad \text{where} \quad (f_{|W|} * f_{|Y|})(z) = \int_{0}^{\infty} f_{|W|} \left(\frac{z}{t}\right) f_{|Y|}(t) \frac{1}{t} dt.$$

We can also express the CDF of $|G|$ this way:

$$F_{|G|}(z) = \int_{0}^{\infty} F_{|W|} \left(\frac{z}{t}\right) f_{|Y|}(t) dt. \qquad (57)$$
# A Activation function with infinite number of stable fixed points for V

## A.1 Proof that V admits an infinite number of fixed points when using $\phi = \varphi_{\delta,\omega}$

**Proposition 2.** For any $\delta \in (0, 1]$ and $\omega > 0$, let us pose the activation function $\phi = \varphi_{\delta,\omega}$. We consider the sequence $(v^l)_l$ defined by:

$$ \begin{array}{l} \forall l \ge 0, \quad v^{l+1} = \sigma_w^2 \int \varphi_{\delta,\omega} (\sqrt{v^l z})^2 \mathcal{D}z + \sigma_b^2, \qquad (20) \\ v^0 \in \mathbb{R}_*^+. \end{array} $$

Then there exists $\sigma_w > 0$, $\sigma_b \ge 0$, and a strictly increasing sequence of stable fixed points $(v_k^*)_{k \in \mathbb{Z}}$ of the recurrence Equation (20).

*Proof.* Let us define:

$$ \begin{array}{l} \tilde{\mathcal{V}}(v) := \frac{1}{v} \mathcal{V}(v|\sigma_w = 1, \sigma_b = 0), \\ \text{so we have:} \quad \mathcal{V}(v|\sigma_w, \sigma_b) = \sigma_w^2 v \tilde{\mathcal{V}}(v) + \sigma_b^2. \end{array} \qquad (21) $$

In the following, we use a simplified notation: $\mathcal{V}(v) = \mathcal{V}(v|\sigma_w, \sigma_b)$.

Our goal is to find a sequence $(v_k^*)_{k \in \mathbb{Z}}$, $\sigma_w$ and $\sigma_b$ such that:

$$ \forall k \in \mathbb{Z}, \quad \mathcal{V}(v_k^*) = v_k^* \quad \text{and} \quad \mathcal{V}'(v_k^*) \in (-1, 1), $$

which would ensure that all $v_k^*$ are stable fixed points of $\mathcal{V}$. In order to understand how to build the sequence $(v_k^*)_{k \in \mathbb{Z}}$, let us consider $v > 0$. Let $\sigma_b^2 = 0$ and $\sigma_w^2 = 1/\tilde{\mathcal{V}}(v)$. So we have:

$$ \mathcal{V}(v) = v. $$

So, any $v > 0$ can possibly be a fixed point if we tune $\sigma_w$ accordingly. We just have to find a $v > 0$ such that:

$$ \mathcal{V}'(v) \in (-1, 1). $$

with:

$$ \begin{align} \mathcal{V}'(v) &= \frac{1}{\tilde{\mathcal{V}}(v)} (\tilde{\mathcal{V}}(v) + v \tilde{\mathcal{V}}'(v)) \nonumber \\ &= 1 + \tilde{\mathcal{V}}_e'(\ln(v)), \tag{22} \end{align} $$

where $\tilde{\mathcal{V}}_e : r \mapsto \ln(\tilde{\mathcal{V}}(\exp(r)))$. So, knowing that $\tilde{\mathcal{V}}_e$ is $\mathcal{C}^1$ and periodic, it is sufficient to prove that it is not constant to ensure that we can extract one $v_0^*$ such that $\mathcal{V}'(v_0^*) \in (-1, 1)$. Then, by periodicity of $\tilde{\mathcal{V}}_e$, we can build a sequence of stable fixed points $(v_k^*)_{k \in \mathbb{Z}}$.
Then, we can use the following property of the Mellin transform $\mathcal{M}$:

$$ \mathcal{M}f_{|G|} = (\mathcal{M}f_{|W|}) \cdot (\mathcal{M}f_{|Y|}), \quad \text{where} \quad (\mathcal{M}f)(t) = \int_{0}^{\infty} x^{t-1}f(x) \, dx. $$

In short, $\mathcal{M}$ transforms a product-convolution into a product in the same manner as the Fourier transform $\mathcal{F}$ transforms a convolution into a product. We have then:

$$ f_{|Y|}(y) := \mathcal{M}^{-1} \left[ \frac{\mathcal{M}f_{|G|}}{\mathcal{M}f_{|W|}} \right](y). $$

Then, by symmetry, we can obtain $f_Y$ from $f_{|Y|}$. However, while $\mathcal{M}f_{|G|}$ and $\mathcal{M}f_{|W|}$ are easy to compute, the inverse Mellin transform $\mathcal{M}^{-1}$ seems to be analytically untractable in this case:

$$ \begin{gathered} (\mathcal{M}f_{|G|})(s) = \frac{2^{\frac{s}{2}-\frac{1}{2}}\Gamma(\frac{s}{2})}{\sqrt{\pi}}, \quad (\mathcal{M}f_{|W|})(s) = \Gamma\left(\frac{s-1}{\theta}+1\right), \\ \text{so:} \quad \frac{(\mathcal{M}f_{|G|})(s)}{(\mathcal{M}f_{|W|})(s)} = \frac{1}{\sqrt{\pi}} \frac{2^{\frac{s}{2}-\frac{1}{2}}\Gamma(\frac{s}{2})}{\Gamma\left(\frac{s-1}{\theta}+1\right)}. \end{gathered} $$

## E.2 Numerical inversion of the Mellin transform

**Computation of $f_{|Y|}$ by numerical inverse Mellin transform.** The Mellin transform of a function can be inverted by using Laguerre polynomials. Specifically, we use the method proposed by Theocaris & Chrysakis (1977) and slightly accelerated by the numerical procedure of Gabutti & Sacripante (1991):

$$ (\mathcal{M}^{-1}f)(z) = e^{-\frac{s}{2}} \sum_{k=0}^{\infty} c_{k+1} L_k \left(\frac{z}{2}\right), \quad \text{with } c_k := \sum_{n=1}^{k} \binom{k-1}{n-1} (-1)^{n-1} \frac{f(n)}{2^n \Gamma(n)}, \quad (58) $$

where the $(L_k)_k$ are the Laguerre polynomials (see Section 7.41, Gradshteyn & Ryzhik, 2014).

**Experiments.** A common way of computing the inverse Mellin transform consists of using Equation (58). Specifically, the sequence $(c_k)_k$ must be computed.

In order to compute the density $f_Y$ of $Q_\theta$ with $\theta = 2.05$, we have computed numerically $(c_k)_k$ for $k \in [1, 500]$. The results are plotted in Figure 12. We have tested three methods to compute the $c_k$:

* floating-point operations using 64 bits floats;

* floating-point operations using 128 bits floats;

* SymPy: make the whole computation using SymPy, a Python library of symbolic computation (very slow).

In all three cases, instabilities appear before the sequence $(c_k)_k$ has fully converged to 0. Moreover, the oscillations of $(c_k)_k$ around 0 have an increasing wavelength, which indicates that we may have to go far beyond $k = 500$ to get enough coefficients $(c_k)_k$ to reconstruct the wanted inverse Mellin transform.

In Figure 13, we have plotted two estimations of $f_{|Y|}$: the density obtained directly by using $(c_k)_{k \in [1,300]}$, and the density obtained by using a sequence $(c_k)_{k \in [1,20000]}$ where the values $(c_k)_{k \in [301,20000]}$ have been extrapolated from $(c_k)_{k \in [1,300]}$.$^{22}$

The resulting estimations of the density $f_{|Y|}$ take negative values and seem to be noised. So, more work is needed to obtain smooth and proper densities, especially if we want them to meet Constraints 3 and 4 (density at 0 and decay rate at $\infty$).

<footnote>$^{22}$The extrapolation has been performed by modeling the graph of $(c_k)_k$ as the product of a decreasing function and a cosine with decreasing frequency.</footnote>
We have:

$$
\begin{align*}
\tilde{\mathcal{V}}(v) &= \frac{1}{v} \int_{-\infty}^{\infty} \varphi_{\delta, \omega}(\sqrt{v}z)^2 \mathcal{D}z \\
&= \int_{-\infty}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln |\sqrt{v}z|) \right) \mathcal{D}z \\
&= 2 \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln(\sqrt{v}z)) \right) \mathcal{D}z \\
&= 2 \int_{0}^{\infty} z^2 \exp \left( 2 \frac{\delta}{\omega} \sin \left( \frac{\omega}{2} \ln v + \omega \ln z \right) \right) \mathcal{D}z.
\end{align*}
$$

**Lemma 1.** $\tilde{\mathcal{V}}$ is not constant.

Proof.

$$
\begin{align*}
\tilde{\mathcal{V}}(v) &= \frac{1}{\sqrt{2\pi}} \frac{2}{v^{3/2}} \int_0^\infty z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln z) \right) \exp \left( -\frac{z^2}{2v} \right) dz \\
&= \frac{1}{\sqrt{2\pi}} \frac{2}{v^{3/2}} \int_0^\infty \frac{z}{2\sqrt{z}} \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln \sqrt{z}) \right) \exp \left( -\frac{z}{2v} \right) dz \\
&= \frac{1}{\sqrt{2\pi}} v^{-3/2} \mathcal{L}[\varphi_{\delta,\omega}(\sqrt{z})] \left(\frac{1}{2v}\right),
\end{align*}
$$

where $\mathcal{L}$ is the Laplace transform.

Let us suppose that $\tilde{\mathcal{V}}$ is constant: $\forall v > 0, \tilde{\mathcal{V}}(v) = c$. So, if we pose $v \leftarrow \frac{1}{2v}$, we have:

$$
\forall v > 0, \quad c = \frac{2}{\sqrt{\pi}} v^{3/2} \mathcal{L}[\varphi_{\delta, \omega}(\sqrt{z})](v),
$$

that is:

$$
\mathcal{L}[\varphi_{\delta, \omega}(\sqrt{z})](v) = \frac{c\sqrt{\pi}}{2}v^{-3/2}.
$$

Since the function $z \mapsto \varphi_{\delta,\omega}(\sqrt{z})$ is continuous on $\mathbb{R}^+$, then, almost everywhere (see Thm. 22.2, Billingsley, 1995):

$$
\varphi_{\delta,\omega}(\sqrt{z}) = \frac{c\sqrt{\pi}}{2} \mathcal{L}^{-1}[v^{-3/2}](z) = c\sqrt{z},
$$

which is impossible for $\delta \in (0, 1]$. Hence the result.

<figure><img src="image_20.png" /></figure>

The function $\tilde{\mathcal{V}}_e : r \mapsto \ln(\tilde{\mathcal{V}}(\exp(r)))$ is continuous and $\frac{4\pi}{\omega}$-periodic:

$$
\tilde{\mathcal{V}}_e(r) = \ln \left[ 2 \int_0^\infty z^2 \exp \left( 2 \frac{\delta}{\omega} \sin \left( \frac{\omega}{2} r + \omega \ln z \right) \right) \mathcal{D}z \right],
$$
<figure><img src="image_18.png" /><figcaption>Figure 12: Evolution of various numerical computations of $c_k$ as $k$ grows.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 13: Numerical inverse Mellin transform with two different computations of $(c_k)_k$: direct computation of $(c_k)_{k \in [1,N]}$ with $N = 300$; extrapolation of $(c_k)_{k \in [301,20000]}$ from $(c_k)_{k \in [1,300]}$.</figcaption></figure>

**Conclusion.** We observe that this computation of the inverse Mellin transform has several intrinsic problems:

* the computation of the $c_k$ coefficients involves a sum of terms with alternating signs, which become larger (in absolute value) as $k$ grows and which are supposed to compensate such that $c_k \to 0$ as $k \to \infty$. Such a numerical computation, involving both large and small terms, makes the resulting $c_k$ very unstable as $k$ grows;

* when $\theta \approx 2$, the sequence $(c_k)_k$ tends extremely slowly to 0;

* if we approximate $M^{-1}f$ with the finite sum of the first $K$ terms of the series in Equation (58), we cannot guarantee the non-negativeness of the resulting function, which is meant to be a density.

So, this method is unpractical to compute the density of a distribution in our case.
so $\tilde{\mathcal{V}}_e$ is lower and upper bounded and reach its bounds (and, by Lemma 1, these bounds are different). We define:

$$ 
\begin{aligned}
\tilde{\mathcal{V}}_e^+ &= \max \tilde{\mathcal{V}}_e & r_0^+ &= \inf \{r > 0 : \tilde{\mathcal{V}}_e(r) = \tilde{\mathcal{V}}_e^+\}, \\
\tilde{\mathcal{V}}_e^- &= \min \tilde{\mathcal{V}}_e & r_0^- &= \inf \{r > r_0^+ : \tilde{\mathcal{V}}_e(r) = \tilde{\mathcal{V}}_e^-\}.
\end{aligned}
$$

By continuity, $\tilde{\mathcal{V}}_e(r_0^+) = \tilde{\mathcal{V}}_e^+$ and $\tilde{\mathcal{V}}_e(r_0^-) = \tilde{\mathcal{V}}_e^-$ . Since $\tilde{\mathcal{V}}^+ > \tilde{\mathcal{V}}^-$ and $\tilde{\mathcal{V}}_e$ is $\mathcal{C}^1$, then there exists $r_0^* \in (r_0^+, r_0^-)$ such that $\tilde{\mathcal{V}}'_e(r_0^*) \in (-2, 0)$.

Since $\tilde{\mathcal{V}}_e$ is $\frac{4\pi}{\omega}$-periodic, we can define a sequence $(r_k^*)_{k \in \mathbb{Z}}$ such that:

$$ 
\begin{aligned}
r_k^* &:= r_0^* + \frac{4k\pi}{\omega} \\
\tilde{\mathcal{V}}_e(r_k^*) &= \tilde{\mathcal{V}}_e(r_0^*) := \tilde{\mathcal{V}}_e^0 \\
\tilde{\mathcal{V}}'_e(r_k^*) &= \tilde{\mathcal{V}}'_e(r_0^*) \in (-2, 0).
\end{aligned}
$$

So, by using Eqn. (21) and Eqn. (22) with $\sigma_b^2 = 0$ and $\sigma_w^2 = 1/\exp(\tilde{\mathcal{V}}_e^0)$:

$$ 
\begin{aligned}
v_k^* &:= \exp(r_k), \\
\mathcal{V}(v_k^*) &= \frac{1}{\exp(\tilde{\mathcal{V}}_e^0)} v_k^* \tilde{\mathcal{V}}(v_k^*) = v_k^* \\
\mathcal{V}'(v_k^*) &= 1 + \tilde{\mathcal{V}}'_e(r_0^*) \in (-1, 1).
\end{aligned}
$$

Thus, $(v_k^*)_{k \in \mathbb{Z}}$ is a sequence of stable fixed points of $\mathcal{V}(\cdot|\sigma_w, \sigma_b)$ for well-chosen $\sigma_w$ and $\sigma_b$. $\square$

## A.2 Practical computation of $\sigma_w^2$

We propose a practical method to ensure that $\mathcal{V}(\cdot|\sigma_w, \sigma_b = 0)$ admits an infinite number of stable fixed points when using activation function $\phi = \varphi_{\delta,\omega}$.

In order to achieve this goal, we build $\sigma_w^2 = \sigma_\omega^2$ in the following way:

$$ 
\begin{aligned}
\mathcal{V}_{\text{low}} &:= 2 \int_0^\infty z^2 \exp \left( -2 \frac{\delta}{\omega} \sin(\omega \ln(z)) \right) \mathcal{D}z \\
\mathcal{V}_{\text{upp}} &:= 2 \int_0^\infty z^2 \exp \left( 2 \frac{\delta}{\omega} \sin(\omega \ln(z)) \right) \mathcal{D}z, \\
\sigma_\omega^2 &:= \left[ \frac{\mathcal{V}_{\text{low}} + \mathcal{V}_{\text{upp}}}{2} \right]^{-1}.
\end{aligned}
$$

In practice, for $\delta = 0.99$, we obtain $\sigma_\omega$ for $\omega \in \{2, 3, 6\}$:

$$ \sigma_2 \approx 0.879, \quad \sigma_3 \approx 0.945, \quad \sigma_6 \approx 0.987. $$
## F Obtaining $f_{|Y|}$: experimental details

We optimize the vector of parameters $\Lambda$ with respect to the following loss:

$$ \ell(\Lambda) := \|\hat{F}_{\Lambda} - F|_{G}\|_{\infty} \\ \hat{F}_{\Lambda}(z) := \int_{0}^{\infty} F|_{W|} \left(\frac{z}{t}\right) g_{\Lambda}(t) dt. $$

**Dataset.** We build the dataset $\mathcal{Z}$ of size $d$:

$$ \mathcal{Z} = \left\{ 0, z_{\max} \frac{1}{d-1}, z_{\max} \frac{2}{d-1}, \dots, z_{\max} \right\}. $$

In our setup, $d = 200$ and $z_{\max} = 5$.

**Computing the loss.** For each $z$ in $\mathcal{Z}$, we compute numerically $\hat{F}_{\Lambda}(z)$. Then, we are able to compute $\ell(\Lambda)$. We keep track of the computational graph with PyTorch, in order to backpropagate the gradient and train the parameters $\Lambda$ by gradient descent.

**Initialization of the parameters.** We initialize $\Lambda = (\alpha, \gamma, \lambda_1, \lambda_2)$ in the following way: $\alpha = 3; \gamma = 1; \lambda_1 = 1; \lambda_2 = 1$.

**Optimizer.** We use the Adam optimizer (Kingma & Ba, 2015) with the parameters: learning rate = 0.001; $\beta_1 = 0.9; \beta_2 = 0.999$; weight decay = 0. We train $\Lambda$ for 100 epochs.

**Learning rate scheduler.** We use a learning rate scheduler based on the reduction of the training loss. If the training loss does not decrease at least by a factor 0.01 for 20 epochs, then the learning rate is multiplied by a factor $1/\sqrt{10}$. After a modification of the learning rate, we wait at least 20 epochs before any modification.

**Scheduler for $\theta'$.** We recall that the definition of $g_{\Lambda}$ involves $\theta'$, defined by: $\frac{1}{\theta} + \frac{1}{\theta'} = \frac{1}{2}$. It is not a parameter to train. Empirically, we found that the following schedule improves the optimization process:

* from epoch 0 to epoch 49, $\theta'$ increases linearly from 2 to its theoretical value $(\frac{1}{2} - \frac{1}{\theta})^{-1}$;

* at the beginning of epoch 50, we reinitialize the optimizer and the learning rate scheduler;

* we finish the training normally, with $\theta' = (\frac{1}{2} - \frac{1}{\theta})^{-1}$.

## G The Kolmogorov-Smirnov test

**Description.** We describe here the Kolmogorov-Smirnov (KS) test (Kolmogoroff, 1941, Smirnov, 1948). Given a sequence $(Z'_1, \dots, Z'_s)$ of $s$ i.i.d. random variables sampled from $P'$:

1. we build the empirical CDF $F_s$ of this sample: $F_s(z) = \frac{1}{s} \sum_{k=1}^{s} 1_{Z'_k \le z}$;

2. we compare $F_s$ to the CDF $F_G$ of $G \sim \mathcal{N}(0, 1)$ by using the $\mathcal{L}^{\infty}$ norm: $D_s = \|F_s - F_G\|_{\infty}$, where $D_s$ is the “KS statistic”;

3. under the null hypothesis, i.e. $P' = \mathcal{N}(0, 1)$, we have: $\sqrt{s}D_s \xrightarrow{d} K$, where $K$ is the Kolmogorov distribution (Smirnov, 1948). We denote by $(K_\alpha)_\alpha$ the quantiles of $K$: $\mathbb{P}(K \le K_\alpha) = 1 - \alpha$, for all $\alpha \in [0, 1]$;

4. finally, we reject the null hypothesis at level $\alpha$ if: $\sqrt{s}D_s \le K_\alpha$.
## B Discussion about the independence of the pre-activations

In Proposition 3 and Constraint 1, we assume that, for any layer $l$, its inputs $(Z_j^l)_j$ are independent. In general, this is not true:

**Example 1.** Let $X \sim \mathcal{N}(0, 1)$ be some random input of a two-layer neural network. We perform the following operation:

$$Z = \frac{1}{\sqrt{2}} [W_1^2 \phi(W_1^1 X) + W_2^2 \phi(W_2^1 X)],$$

where $(W_1^1, W_2^1, W_1^2, W_2^2)$ be i.i.d. random variables samples from some distribution P, and $\phi$ is some activation function.

Let $P = \mathcal{R}$, the Rademacher distribution, i.e., if $W \sim \mathcal{R}$, then $W = \pm 1$ with probability 1/2.

Let $\phi = \text{Id}$. Then we have: $Y_1 := W_1^1 X \sim \mathcal{N}(0, 1)$ and $Y_2 := W_2^1 X \sim \mathcal{N}(0, 1)$. But they are not independent:

$$W_1^2 W_1^1 X + W_2^2 W_2^1 X = (W'_1 + W'_2)X,$$

where $W_1^2 W_1^1$ and $W_2^2 W_2^1$ are two independent Rademacher random variables. So, $Z = 0$ with probability 1/2. So, $Z$ is not Gaussian.

In this example, we build a non-Gaussian random variables with a minimal neural network, in which we construct two dependent random variables. So, we should pay attention to this phenomenon when propagating the pre-activations in a neural network.

One should note that the structure of dependence of $W_1^2 W_1^1 X$ and $W_2^2 W_2^1 X$ does not involve their correlation (which is zero), and yet breaks the Gaussianity of their sum. So, in order to obtain a theoretical result about the distribution of $Z$, we should study finer aspects of the dependence structure.

But, on the practical side, we want to answer the question: to which extent does the relation of dependence between the pre-activations $(Z_j^l)_j$ affect the Gaussianity of $(Z_i^{l+1})_i$?

In order to answer this question, we propose a series of experiments on a two-layer neural network. We consider a vector of inputs $X \in \mathbb{R}^{n_0}$, where the $(X_j)_{1 \le j \le n_0}$ are $\mathcal{N}(0, 1)$ and i.i.d. The scalar outputted by the network is:

$$Z = \frac{1}{\sqrt{n_1}} \mathbf{W}^2 \phi \left( \frac{1}{\sqrt{n_0}} \mathbf{W}^1 X \right),$$

where the weights $\mathbf{W}^1 \in \mathbb{R}^{n_1 \times n_0}$ and $\mathbf{W}^2 \in \mathbb{R}^{1 \times n_1}$ are i.i.d.

We test this setup with different initialization distributions P and activation functions $\phi$:

* usual ones: $\phi = \tanh$ or ReLU, $P = \mathcal{N}(0, \sigma_w^2)$, where $\sigma_w^2$ is such that the pair $(\sigma_w^2, \sigma_b^2 = 0.01)$ lies at the EOC;
**Limitations.** In the KS test presented above, the null hypothesis $H_0$ is $P' = \mathcal{N}(0, 1)$, which is exactly what we intend to demonstrate when using our family $\{(P_\theta, \phi_\theta^{o/p}) : \theta \in (2, \infty)\}$, while the alternative hypothesis $H_1$ is $P' \neq \mathcal{N}(0, 1)$. With this test design, we face a problem: it is impossible to claim that $H_0$ holds. More precisely, we have two possible outcomes: either $H_0$ is rejected, or it is not. Since $H_1$ is the complementary of $H_0$, a reject of $H_0$ means an accept of $H_1$. But the converse does not hold: if $H_0$ is not rejected, then it is impossible to conclude that $H_0$ is true: the sample size $s$ may simply be too small, or the KS test may be inadequate for our case.

So, to be able to conclude that the pre-activations are Gaussian, we should have built an alternative test, where the null hypothesis $H'_0$ is some misfit between $P'$ and $\mathcal{N}(0, 1)$. Or, at least, we should have computed the power of the current KS test.

However, our experimental results are sufficient to *compare* the quality of fit between $P'$ and $\mathcal{N}(0, 1)$ in the tested setups (see Figures 6 and 7). It remains clear that, with $\{(P_\theta, \phi_\theta) : \theta \in (2, \infty)\}$, the Gaussian pre-activations hypothesis is far more likely to hold than in setups involving tanh, ReLU and Gaussian weights.

# H Experiments

## H.1 Propagation of the correlations with MNIST

Within the setup of Section 3.2, we have plotted in Figure 14 the correlations propagated in a multilayer perceptron with inputs sampled for MNIST (10 samples per class, that is, 100 samples in total). The layers of the perceptron have $n_l = 10$ neurons each, and the activation function is $\phi = \text{ReLU}$.

We observe that, in this narrow NN case ($n_l = 10$), some irregularities appear as information propagates into the network: some classes seem to (de)correlate in an inconsistent way with the others. Specifically, while most of the classes tend to correlate exactly ($C^* \approx 1$), the third class (corresponding to the digit “2”) tends to a correlation $C^* \approx 0.7 \neq 1$. This result is contradictory to the EOC theory.

<figure><img src="image_20.png" /><figcaption>Figure 14: Propagation of correlations $c_{ab}^l$ in a multilayer perceptron with activation function $\phi = \text{ReLU}$ and inputs sampled from the MNIST dataset. The neural network is initialized at the EOC. Each plot displays a $10 \times 10$ matrix $C_{pq}^l$ whose entries are the average correlation between the pre-activations propagated by samples from classes $p, q \in \{0, \dots, 9\}$, at the input and right after layers $l \in \{10, 30, 50\}$. See also Figure 1 in Section 3.2 for results on CIFAR-10.</figcaption></figure>

## H.2 Propagation of the correlations with $\phi = \phi_\theta^o$

Within the setup of Section 3.2, we have plotted in Figure 15 the correlations propagated in a multilayer perceptron with $\phi = \phi_\theta$. The weights have been sampled from $\mathcal{W}(\theta, 1)$.
• ours: $\phi = \phi_{\theta}$, $P = P_{\theta} = W(\theta, 1)$.

We study three cases:

• unfavorable case: $n_0 = 1$, various $n_1 \in [1, 10]$; the intermediary features $\mathbf{W}^1\mathbf{X}$ are weakly “mixed”, so it is credible that they lead to an output that is far from being Gaussian;

• favorable case: $n_1 = 2$, various $n_0 \in [1, 10]$;<fnref n="20" /> the intermediary features $\mathbf{W}^1\mathbf{X}$ are “mixed” with an increasing rate as $n_0$ increases;

• same-width case: $n_0 = n_1$.

<figure><img src="image_21.png" /><figcaption>KS-statistic for various layer sizes in a 2-layer NN</figcaption></figure>

Figure 11: Evolution of the distance of the standardized distribution of $Z$ to the $N(0, 1)$ according to the Kolmogorov-Smirnov statistic. Weight initialization is $W(\theta, 1)$ when using $\phi = \phi_{\theta}$ and is Gaussian according to the EOC when using $\phi = \tanh$ or ReLU. For each point, we have computed the KS-statistic over 200 000 samples of $(\mathbf{X}, \mathbf{W}^1, \mathbf{W}^2)$.

According to Figure 11:

• unfavorable case (1st graph, $n_0 = 1$): $Z$ is far from being Gaussian, regardless of the activation function. One should note that $\tanh$ and $(P_{\theta}, \phi_{\theta})$ with small $\theta$ lead a distribution of $Z$ that is closer to the Gaussian than with ReLU and $(P_{\theta}, \phi_{\theta})$ with large $\theta$;

• favorable case (2nd graph, $n_1 = 2$): $Z$ is closer to be Gaussian with out setup $(P_{\theta}, \phi_{\theta})$ than with $\phi = \tanh$ or ReLU, especially with larger $n_1$;

• same-width case (3rd graph, $n = n_0 = n_1$): the larger the width $n$, the closer $Z$ is to being Gaussian. For a fixed $n$, the distribution of $Z$ is closer to a Gaussian with smaller $\theta$. When we use $(P_{\theta}, \phi_{\theta})$, we are close to the performance of $\tanh$ or better.

<footnote n="20">²⁰The case $n_1 = 1$ is trivial: there is no sum of dependent random variables in the second layer.</footnote>
In this setup, the results are consistent, and are consistent with the results for $\phi = \tanh$ (see Figure 1): the sequence $(C_{pq}^l)_l$ converges to 1, which was not the case for $\phi = \text{ReLU}$ and $n_l = 10$.

### H.3 Variance of the pre-activation when using $\phi = \phi_{\theta}^{0}$

As observed in Figure 6a, the product $W\phi_{\theta}(X)$ is above the KS threshold, corresponding to $s = 10^7$ samples and a $p$-value of 0.05. This is not necessarily the case in Figure 6b, where the samples are standardized. So, we suspect that the variance of $W\phi_{\theta}^{0}(X)$ is not exactly 1.

Therefore, we have reported in Table 2 the empirical standard deviation of the product $W\phi(X)$ (i.e., $Z'$ with $n = 1$), computed with $s = 10^7$ samples, where $X \sim \mathcal{N}(0, 1)$, $W \sim W(\theta, 1)$ if $\phi = \phi_{\theta}$ and $W \sim \mathcal{N}(0, 1)$ if $\phi = \tanh$ or $\text{ReLU}$. We observe that the standard deviation of $Z'$, which is expected to be 1 with $\phi = \phi_{\theta}^{0}$, is actually a bit different. This would largely explain the differences between Figure 6a and Figure 6b.

<table><caption>Table 2: Empirical standard deviation of $Z'$ with $n = 1$.</caption><thead><tr><th></th><th>2.05</th><th>2.5</th><th>3</th><th>4</th><th>5</th><th>7</th><th>10</th><th>tanh</th><th>ReLU</th></tr></thead><tbody><tr><td>θ</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>σ̄</td><td>1.003</td><td>0.994</td><td>0.989</td><td>0.984</td><td>0.981</td><td>0.979</td><td>0.978</td><td>0.628</td><td>0.707</td></tr></tbody></table>

### H.4 Propagation of the pre-activations

In this subsection, we show how the choice of the input data points affect the propagation of the pre-activations.

Let $\mathcal{D}^l$ be the distribution of the pre-activation $Z_1^l$ after layer $l$. In Figure 16, we have plotted the $\mathcal{L}^\infty$ distance between the CDFs of $\mathcal{D}^l$ and $\mathcal{N}(0, 1)$ for various input points, sampled from different classes to improve diversity between them. We have chosen the following setup: CIFAR-10 inputs, multilayer perceptron with 100 layers and 100 neurons per layer, Gaussian initialization at the EOC when using the activation function $\phi = \tanh$ or $\text{ReLU}$, and symmetric Weibull initialization $W(\theta, 1)$ when using $\phi = \phi_{\theta}^{0}$ or $\phi = \phi_{\theta}^{P}$.

This setup has been selected to illustrate clearly the variability of $\mathcal{D}^l$ when using various data points. This variability can be observed with 10 or 1000 neurons per layer, though it is less striking.

**Input normalization over the whole dataset.** We perform the usual normalization over the whole dataset:

$$ 
\begin{aligned}
\hat{x}_{a;ij} &:= \frac{x_{a;ij} - \mu_i}{\sigma_i}, \\
\mu_i &:= \frac{1}{Np_i} \sum_{\mathbf{x} \in \mathcal{D}} \sum_{j=1}^{p_i} x_{ij} \\
\sigma_i^2 &:= \frac{1}{Np_i - 1} \sum_{\mathbf{x} \in \mathcal{D}} \sum_{j=1}^{p_i} (x_{ij} - \mu_i)^2,
\end{aligned}
$$

where $\hat{x}$ is the normalized data point, $x_{a;ij}$ is the $j$-th component of the $i$-th channel of the input image $\mathbf{x}_a$, $p_i$ is the size of the $i$-th channel, and $N$ is the size of the dataset $\mathcal{D}$.

**Results.** In Figure 16, we observe that, when using the activation function $\phi = \tanh$, the shape of the curves is the same for all data points. This is not the case for $\text{ReLU}$: for the input point “bird”, the sequence $(\mathcal{D}^l)_l$ drifts away from $\mathcal{N}(0, 1)$ since the beginning, while for the input “car”, $(\mathcal{D}^l)_l$ first becomes closer to $\mathcal{N}(0, 1)$, then drifts away. It is also the case for $\phi_{\theta}$ with $\theta = 10$: for the input “deer”, the distance remains high, while it starts low and increases very slowly (input “dog”), or even decreases in the first place (input
The case $n_0 = n_1$ is the most realistic one: usually, the sizes of the layers of a neural network are of the same order of magnitude. In this case, our setup is better than or equivalent to the one with tanh.

**Remark 13** (Mixing of inputs and interference phenomenon). We see in Figure 11 that, in every graph, setups with small $\theta$ lead to better results than the ones with large $\theta$. This observation could be explained by Example 1. In this example, the inputs are weakly “mixed”: since $X_0$ is multiplied by a Rademacher random variable, it is possible to partially reconstruct $X_0$ after the first layer, and then build destructive interference (leading to an output $Z = 0$ half of the time).

So, if we want to avoid this “interference” phenomenon, we should use initialization distributions and activation functions such that every layer “mixes” strongly the inputs. Notably, initialization distributions should be far from being a combination of Dirac distributions. Typically, $P_{\theta=10}$ is close to the Rademacher distribution (see Fig. 5a).

**Remark 14** (Case $\phi = \tanh$ and $P = \mathcal{R}$). According to Remark 13, the relatively good performance of tanh should be explained by the choice of a Gaussian initialization $P$. However, the case $\phi = \tanh$ and $P = \mathcal{R}$, shown in Fig. 11 (dotted magenta line), contradicts partially this explanation. In this case, choosing a Rademacher initialization of the weights should lead to worse results than with a Gaussian initialization.

It is certainly true in the first graph ($n_0 = 1$), but it is obviously false in the third graph ($n = n_0 = n_1$) for larger $n$. Actually, the setup $\phi = \tanh$ and $P = \mathcal{R}$ is better than all the others in this graph.

So, the contradiction between Remarks 13 and 14 indicates that the dependence between the inputs remains to be investigated in depth. Specifically, we did not take into account the shape of the activation function in Remark 13. For instance, tanh-like functions make the network “forget” information about the input, which may improve the “mixing” process, and thus lead to better results

## C Constraints on the product of two random variables

**Proposition 4** (Density of a product of random variables at 0). Let $W, Y$ be two independent non-negative random variables and $Z = WY$. Let $f_W, f_Y, f_Z$ be their respective density. Assuming that $f_Y$ is continuous at 0 with $f_Y(0) > 0$, we have:

$$ \text{if } \lim_{w \to 0} \int_{w}^{\infty} \frac{f_W(t)}{t} dt = \infty, \quad \text{then} \quad \lim_{z \to 0} f_Z(z) = \infty. \qquad (23) $$

Moreover, if $f_Y$ is bounded:

$$ \text{if } \int_{0}^{\infty} \frac{f_W(t)}{t} dt < \infty, \quad \text{then} \quad f_Z(0) = f_Y(0) \int_{0}^{\infty} \frac{f_W(t)}{t} dt. \qquad (24) $$
<figure><img src="image_21.png" /><figcaption>Figure 15: Propagation of correlations $c_{ab}^l$ in a multilayer perceptron with activation function $\phi_{\theta}^o$ with $\theta \in \{2.05, 7.00\}$ and inputs sampled from the CIFAR-10 dataset. The weights are sampled from $W(\theta, 1)$ and the biases are zero. Each plot displays a $10 \times 10$ matrix $C_{pq}^l$ whose entries are the average correlation between the pre-activations propagated by samples from classes $p, q \in \{0, \dots, 9\}$, at the input and right after layers $l \in \{10, 30, 50\}$.</figcaption></figure>
Proof. Let $z, z_0 > 0$:

$$ 
\begin{aligned}
f_Z(z) &= \int_0^\infty f_Y(t) \frac{1}{t} f_W\left(\frac{z}{t}\right) dt \\
&\geq \int_0^{z_0} f_Y(t) \frac{1}{t} f_W\left(\frac{z}{t}\right) dt \\
&\geq \inf_{[0,z_0]} f_Y \cdot \int_0^{z_0} \frac{1}{t} f_W\left(\frac{z}{t}\right) dt \\
&\geq \inf_{[0,z_0]} f_Y \cdot \int_{z/z_0}^\infty \frac{f_W(t)}{t} dt.
\end{aligned}
$$

Let us take $z_0 = \sqrt{z}$. We have:

$$ f_Z(z) \geq \inf_{[0,\sqrt{z}]} f_Y \cdot \int_{\sqrt{z}}^{\infty} \frac{f_W(t)}{t} dt. $$

Then we take the limit $z \to 0$, hence:

• if $\int_0^\infty \frac{fw(t)}{t} dt = \infty$, then: $\lim_{z \to 0} f_Z(z) = \infty$, which achieves (23);

• if $\int_0^\infty \frac{fw(t)}{t} dt < \infty$, then: $f_Z(0) \geq f_Y(0) \int_0^\infty \frac{fw(t)}{t} dt$, which achieves one half of (24);

Let us prove the second half of (24). Let $z, z_0 > 0$:

$$ 
\begin{aligned}
f_Z(z) &= \int_0^\infty f_Y(t) \frac{1}{t} f_W\left(\frac{z}{t}\right) dt \\
&= \int_0^{z_0} f_Y(t) \frac{1}{t} f_W\left(\frac{z}{t}\right) dt + \int_{z_0}^\infty f_Y(t) \frac{1}{t} f_W\left(\frac{z}{t}\right) dt \\
&\leq \sup_{[0,z_0]} f_Y \cdot \int_{z/z_0}^\infty \frac{f_W(t)}{t} dt + \int_1^\infty f_Y(z_0 t) \frac{1}{t} f_W\left(\frac{z}{z_0 t}\right) dt
\end{aligned}
$$

Let $z_0 = \sqrt{z}$. We have:

$$ f_Z(z) \leq \sup_{[0,\sqrt{z}]} f_Y \cdot \int_{\sqrt{z}}^{\infty} \frac{f_W(t)}{t} dt + \int_1^{\infty} f_Y(\sqrt{z}t) \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt, $$

where:

$$ 
\begin{aligned}
\int_1^\infty f_Y(\sqrt{z}t) \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt &\leq \|f_Y\|_\infty \int_1^\infty \frac{1}{t} f_W\left(\frac{\sqrt{z}}{t}\right) dt \\
&\leq \|f_Y\|_\infty \int_0^{\sqrt{z}} \frac{f_W(t)}{t} dt.
\end{aligned}
$$
“truck”). But, in general, when the curves have converged after 100 layers, it seems that the limit is the same whatever the starting data point, and depends only on the choice of the activation function.

The most striking result lies in the difference between Figures 16 and 17. With $\phi = \phi_{\theta}^{P}$, the pre-activations tend always to be Gaussian, whatever $\theta$ or the input data point. This is not the case with ReLU, tanh and some $\phi_{\theta}^{P}$ with large $\theta$. However, $\phi_{\theta}^{P}$ remain “more Gaussian” than ReLU or tanh.

<figure><img src="image_22.png" /><figcaption>Figure 16: Activations functions tanh, ReLU and $\phi_{\theta}^{O}$.</figcaption></figure>

Propagation of the distribution $\mathcal{D}^l$ of the pre-activations across the layers $l \in [1, 100]$ when inputting various data points of CIFAR-10. Each curve represents the $\mathcal{L}^{\infty}$ distance between the CDF of $\mathcal{D}^l$ and the CDF of the Gaussian $N(0, 1)$. The dashed green line is the threshold of rejection of the Kolmogorov-Smirnov test with p-value 0.05: if a distribution $\mathcal{D}^l$ is represented by a point above this threshold, then the hypothesis $"\mathcal{D}^l = N(0, 1)"$ is rejected with p-value 0.05.
According to the hypotheses, we have, as $z \to 0$:

$$ \sup_{[0,\sqrt{z}]} f_Y \to f_Y(0) \\ \int_{\sqrt{z}}^{\infty} \frac{f_W(t)}{t} dt \to \int_{0}^{\infty} \frac{f_W(t)}{t} dt \\ \int_{0}^{\sqrt{z}} \frac{f_W(t)}{t} dt \to 0, $$

hence the result.

<figure><img src="image_22.png" /></figure>

## D Activation functions with vertical tangent at 0

In the following lemma, we show that if we want the activation $Y$ to have a density that is 0 at 0, then the activation function $\phi$ should have a vertical tangent at 0. $G$ plays the role of pre-activation.

**Lemma 2.** Let $\phi$ be a function transforming a Gaussian random variable $G \sim \mathcal{N}(0, 1)$ into a symmetrical random variable $Y$ with a density $f_Y$ such that $f_Y(0) = 0$. That is, $Y = \phi(G)$. Then $\phi$ has a vertical tangent at 0.

*Proof.* We have:

$$ \phi(x) = F_Y^{-1}(F_G(x)), $$

where $F_G$ and $F_Y$ are the respective CDFs of $G$ and $Y$.

Thus:

$$ \phi'(x) = F'_{G}(x) \frac{1}{F'_{Y}(F^{-1}_{Y}(F_{G}(x)))} $$

Therefore:

$$ \begin{align*} \phi'(0) &= F'_{G}(0) \frac{1}{F'_{Y}(F^{-1}_{Y}(F_{G}(0)))} \\ &= \frac{1}{\sqrt{2\pi} F'_{Y}(F^{-1}_{Y}(1/2))} \\ &= \frac{1}{\sqrt{2\pi} F'_{Y}(0)} \\ &= \infty. \end{align*} $$

<figure><img src="image_23.png" /></figure>
<figure><img src="image_23.png" /><figcaption>Figure 17: Activation functions $\phi_{\theta}^{P}$, same setup as in Figure 16.</figcaption></figure>

## H.5 Experimental details of the training procedure

**Training, validation, and test sets.** For MNIST and CIFAR-10, we split randomly the initial training set into two sets: the training set, which will be actually used to train the neural network, and the validation set, which will be used to stop training when the network begins to overfit.

The sizes of the different sets are as follows:

* MNIST: 50000 training samples; 10000 validation samples; 10000 test samples;

* CIFAR-10: 42000 training samples, 8000 validation samples; 10000 test samples.

The training sets are split into mini-batches with 200 samples each. No data augmentation is performed.

**Loss.** Given a classification task with $P$ classes, let $z^L \in \mathbb{R}^P$ be the pre-activation outputted by the last layer of the neural network. First, we perform a softmax operation:

$$y_p := \mathrm{softmax}(z_p^L) = \frac{\exp(z_p^L)}{\sum_{p'=1}^{P} \exp(z_{p'}^L)}.$$
## E Numerical inversion of the Mellin transform

A common way of computing the inverse Mellin transform consists of using Equation (18). Specifically, the sequence $(c_k)_k$ must be computed.

In order to compute the density $f_Y$ of $Q_\theta$ with $\theta = 2.05$, we have computed numerically $(c_k)_k$ for $k \in [1, 500]$. The results are plotted in Figure 12. We have tested three methods to compute the $c_k$:

* floating-point operations using 64 bits floats;

* floating-point operations using 128 bits floats;

* SymPy: make the whole computation using SymPy, a Python library of symbolic computation (very slow).

In all three cases, instabilities appear before the sequence $(c_k)_k$ has fully converged to 0. Moreover, the oscillations of $(c_k)_k$ around 0 have an increasing wavelength, which indicates that we may have to go far beyond $k = 500$ to get enough coefficients $(c_k)_k$ to reconstruct the wanted inverse Mellin transform.

In Figure 13, we have plotted two estimations of $f_{|Y|}$: the density obtained directly by using $(c_k)_{k \in [1,300]}$, and the density obtained by using a sequence $(c_k)_{k \in [1,20000]}$ where the values $(c_k)_{k \in [301,20000]}$ have been extrapolated from $(c_k)_{k \in [1,300]}$.<fnref n="21" />

The resulting estimations of the density $f_{|Y|}$ take negative values and seem to be noised. So, more work is needed to obtain smooth and proper densities, especially if we want them to meet Constraints 3 and 4 (density at 0 and decay rate at $\infty$).

## F Obtaining $f_{|Y|}$: experimental details

We optimize the vector of parameters $\Lambda$ with respect to the following loss:

$$ \ell(\Lambda) := \|\hat{F}_{\Lambda} - F_{|G|}\|_{\infty} \\ \hat{F}_{\Lambda}(z) := \int_{0}^{\infty} F_{|W|} \left(\frac{z}{t}\right) g_{\Lambda}(t) dt. $$

**Dataset.** We build the dataset $\mathcal{Z}$ of size $d$:

$$ \mathcal{Z} = \left\{ 0, z_{\max} \frac{1}{d-1}, z_{\max} \frac{2}{d-1}, \dots, z_{\max} \right\}. $$

In our setup, $d = 200$ and $z_{\max} = 5$.

<footnote n="21">²¹The extrapolation has been performed by modeling the graph of $(c_k)_k$ as the product of a decreasing function and a cosine with decreasing frequency.</footnote>
where the $(y_p)_p$ are the components of $\mathbf{y} \in \mathbb{R}^P$ and $(z_p^L)_p$ are the components of $\mathbf{z}^L$. Then, we compute the negative log-likelihood loss. For a target class $p \in \{1, \dots, P\}$, we pose:

$$ \ell(\mathbf{y}, p) := -\log(y_p). $$

**Optimizer.** We use the Adam optimizer (Kingma & Ba, 2015) with the parameters: learning rate = 0.001; $\beta_1$ = 0.9; $\beta_2$ = 0.999; weight decay = 0.

**Learning rate scheduler.** We use a learning rate scheduler based on the reduction of the training loss. If the training loss does not decrease at least by a factor 0.01 for 10 epochs, then the learning rate is multiplied by a factor $1/\sqrt{10}$.

**Early stopping.** We add an early stopping rule based on the reduction of the validation loss. If the validation loss does not decrease at least by a factor of 0.001 for 30 epochs, then we stop training.

## H.6 Additional training results: $\phi = \phi_{\theta}^P$

We provide here additional training experiments, with the same setup as in Section 5.4 but with our family of positive activation functions $(\phi_{\theta}^P)_{\theta}$.

**Results with CIFAR-10 + LeNet.** We have plotted in Figure 18 the evolution of the training loss and the test accuracy during training. The test loss is slightly better than with the family of odd activation functions $(\phi_{\theta}^P)_{\theta}$ (see Fig. 9b) and the training loss is slightly worse (see Fig. 9a). Overall, the plots are very similar with both types of activations functions.

<figure><img src="image_24.png" /><figcaption>Figure 18: Training curves for CIFAR-10 + LeNet with 3 different numbers of neurons per layer and $\phi = \phi_{\theta}^P$.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 12: Evolution of various numerical computations of $c_k$ as $k$ grows.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 13: Numerical inverse Mellin transform with two different computations of $(c_k)_k$: direct computation of $(c_k)_{k \in [1,N]}$ with $N = 300$; extrapolation of $(c_k)_{k \in [301,20000]}$ from $(c_k)_{k \in [1,300]}$.</figcaption></figure>

**Computing the loss.** For each $z$ in $\mathcal{Z}$, we compute numerically $\hat{F}_{\Lambda}(z)$. Then, we are able to compute $\ell(\Lambda)$. We keep track of the computational graph with PyTorch, in order to backpropagate the gradient and train the parameters $\Lambda$ by gradient descent.

**Initialization of the parameters.** We initialize $\Lambda = (\alpha, \gamma, \lambda_1, \lambda_2)$ in the following way:

$$ \alpha = 3, \quad \gamma = 1, \quad \lambda_1 = 1, \quad \lambda_2 = 1. $$

**Optimizer.** We use the Adam optimizer (Kingma & Ba, 2015) with the parameters:
**Results with MNIST + Multilayer perceptron.** We have plotted in Figure 19 the results of the training of a narrow multilayer perceptron, as in Figure 10. In this extreme setup, the positive activation functions $\phi_{\theta}^{P}$ perform worse than the odd ones $\phi_{\theta}^{O}$. In particular, many trainings failed with a narrow and deep network ($L = 30$).

<figure><img src="image_25.png" /><figcaption>Figure 19: Training loss for a multilayer perceptron, narrow ($n_l \in \{3, 10\}$) and of various depths ($L \in \{3, 10, 30\}$). Training curves averaged over 5 experiments.</figcaption></figure>
* learning rate: 0.001;

* $\beta_1 = 0.9$;

* $\beta_2 = 0.999$;

* weight decay: 0.

We train $\Lambda$ for 100 epochs.

**Learning rate scheduler.** We use a learning rate scheduler based on the reduction of the training loss. If the training loss does not decrease at least by a factor 0.01 for 20 epochs, then the learning rate is multiplied by a factor $1/\sqrt{10}$. After a modification of the learning rate, we wait at least 20 epochs before any modification.

**Scheduler for $\theta'$.** We recall that the definition of $g_\Lambda$ involves $\theta'$, defined by: $\frac{1}{\theta} + \frac{1}{\theta'} = \frac{1}{2}$. It is not a parameter to train. Empirically, we found that the following schedule improves the optimization process:

* from epoch 0 to epoch 49, $\theta'$ increases linearly from 2 to its theoretical value $(\frac{1}{2} - \frac{1}{\theta})^{-1}$;

* at the beginning of epoch 50, we reinitialize the optimizer and the learning rate scheduler;

* we finish the training normally, with $\theta' = (\frac{1}{2} - \frac{1}{\theta})^{-1}$.

# G Experiments

## G.1 Propagation of the correlations with MNIST

Within the setup of Section 2.2, we have plotted in Figure 14 the correlations propagated in a multilayer perceptron with inputs sampled for MNIST (10 samples per class, that is, 100 samples in total). The layers of the perceptron have $n_l = 10$ neurons each, and the activation function is $\phi = \text{ReLU}$.

We observe that, in this narrow NN case ($n_l = 10$), some irregularities appear as information propagates into the network: some classes seem to (de)correlate in an inconsistent way with the others. Specifically, while most of the classes tend to correlate exactly ($C^* \approx 1$), the third class (corresponding to the digit “2”) tends to a correlation $C^* \approx 0.7 \neq 1$. This result is contradictory to the EOC theory.

## G.2 Propagation of the correlations with $\phi = \phi_\theta$

Within the setup of Section 2.2, we have plotted in Figure 15 the correlations propagated in a multilayer perceptron with $\phi = \phi_\theta$. The weights have been sampled from $W(\theta, 1)$.
<figure><img src="image_26.png" /><figcaption>Figure 14: Propagation of correlations $c_{ab}^l$ in a multilayer perceptron with activation function $\phi = \text{ReLU}$ and inputs sampled from the MNIST dataset. The neural network is initialized at the EOC. Each plot displays a $10 \times 10$ matrix $C_{pq}^l$ whose entries are the average correlation between the pre-activations propagated by samples from classes $p, q \in \{0, \dots, 9\}$, at the input and right after layers $l \in \{10, 30, 50\}$. See also Figure 1 in Section 2.2 for results on CIFAR-10.</figcaption></figure>

In this setup, the results are consistent, and are consistent with the results for $\phi = \tanh$ (see Figure 1): the sequence $(C_{pq}^l)_l$ converges to 1, which was not the case for $\phi = \text{ReLU}$ and $n_l = 10$.

## G.3 Limits of the KS test

In the KS test presented above, the null hypothesis $H_0$ is $P' = \mathcal{N}(0, 1)$, which is exactly what we intend to demonstrate when using our family $\{(P_\theta, \phi_\theta) : \theta \in (2, \infty)\}$, while the alternative hypothesis $H_1$ is $P' \neq \mathcal{N}(0, 1)$. With this test design, we face a problem: it is impossible to claim that $H_0$ holds. More precisely, we have two possible outcomes: either $H_0$ is rejected, or it is not. Since $H_1$ is the complementary of $H_0$, a reject of $H_0$ means an accept of $H_1$. But the converse does not hold: if $H_0$ is not rejected, then it is impossible to conclude that $H_0$ is true: the sample size $s$ may simply be too small, or the KS test may be inadequate for our case.

So, to be able to conclude that the pre-activations are Gaussian, we should have built an alternative test, where the null hypothesis $H_0'$ is some misfit between $P'$ and $\mathcal{N}(0, 1)$. Or, at least, we should have computed the power of the current KS test.

However, our experimental results are sufficient to compare the quality of fit between $P'$ and $\mathcal{N}(0, 1)$ in the tested setups. It remains clear that, with $\{(P_\theta, \phi_\theta) : \theta \in (2, \infty)\}$, the Gaussian pre-activations hypothesis is far more likely to hold than in setups involving tanh, ReLU and Gaussian weights.
<figure><img src="image_27.png" /><figcaption>Figure 15: Propagation of correlations $c_{ab}^l$ in a multilayer perceptron with activation function $\phi_{\theta}$ with $\theta \in \{2.05, 7.00\}$ and inputs sampled from the CIFAR-10 dataset. The weights are sampled from $W(\theta, 1)$ and the biases are zero. Each plot displays a 10 × 10 matrix $C_{pq}^l$ whose entries are the average correlation between the pre-activations propagated by samples from classes $p, q \in \{0, \dots, 9\}$, at the input and right after layers $l \in \{10, 30, 50\}$.</figcaption></figure>
### G.4 Variance of the pre-activation when using $\phi = \phi_{\theta}$

As observed in Figure 6a, the product $W\phi_{\theta}(X)$ is above the KS threshold, corresponding to $s = 10^7$ samples and a $p$-value of 0.05. This is not necessarily the case in Figure 6b, where the samples are standardized. So, we suspect that the variance of $W\phi_{\theta}(X)$ is not exactly 1.

Therefore, we have reported in Table 2 the empirical standard deviation of the product $W\phi(X)$ (i.e., $Z'$ with $n = 1$), computed with $s = 10^7$ samples, where $X \sim \mathcal{N}(0, 1)$, $W \sim W(\theta, 1)$ if $\phi = \phi_{\theta}$ and $W \sim \mathcal{N}(0, 1)$ if $\phi = \tanh$ or ReLU. We observe that the standard deviation of $Z'$, which is expected to be 1 with $\phi = \phi_{\theta}$, is actually a bit different. This would largely explain the differences between Figure 6a and Figure 6b.

<table><caption>Table 2: Empirical standard deviation of $Z'$ with $n = 1$.</caption><thead><tr><th></th><th>2.05</th><th>2.5</th><th>3</th><th>4</th><th>5</th><th>7</th><th>10</th><th>tanh</th><th>ReLU</th></tr></thead><tbody><tr><td>σ̄</td><td>1.003</td><td>0.994</td><td>0.989</td><td>0.984</td><td>0.981</td><td>0.979</td><td>0.978</td><td>0.628</td><td>0.707</td></tr></tbody></table>

### G.5 Propagation of the pre-activations

In this subsection, we show how the choice of the input data points affect the propagation of the pre-activations.

Let $\mathcal{D}^l$ be the distribution of the pre-activation $Z_1^l$ after layer $l$. In Figure 16, we have plotted the $\mathcal{L}^\infty$ distance between the CDFs of $\mathcal{D}^l$ and $\mathcal{N}(0, 1)$ for various input points, sampled from different classes to improve diversity between them. We have chosen the following setup: CIFAR-10 inputs, multilayer perceptron with 100 layers and 100 neurons per layer, Gaussian initialization at the EOC when using the activation function $\phi = \tanh$ or ReLU, and symmetric Weibull initialization $W(\theta, 1)$ when using $\phi = \phi_{\theta}$.

This setup has been selected to illustrate clearly the variability of $\mathcal{D}^l$ when using various data points. This variability can be observed with 10 or 1000 neurons per layer, though it is less striking.

**With input normalization over the whole dataset.** In this setup, we perform the usual normalization over the whole dataset:

$$ 
\begin{align*}
\hat{x}_{a;ij} &:= \frac{x_{a;ij} - \mu_i}{\sigma_i}, \\
\mu_i &:= \frac{1}{N p_i} \sum_{\mathbf{x} \in \mathbb{D}} \sum_{j=1}^{p_i} x_{ij} \\
\sigma_i^2 &:= \frac{1}{N p_i - 1} \sum_{\mathbf{x} \in \mathbb{D}} \sum_{j=1}^{p_i} (x_{ij} - \mu_i)^2,
\end{align*}
$$
where $\hat{x}$ is the normalized data point, $x_{a;ij}$ is the j-th component of the i-th channel of the input image $\mathbf{x}_a$, $p_i$ is the size of the i-th channel, and $N$ is the size of the dataset $\mathcal{D}$.

We observe that, when using the activation function $\phi = \tanh$, the shape of the curves is the same for all data points. This is not the case for ReLU: for the input point “bird”, the sequence $(\mathcal{D}^l)_l$ drifts away from $\mathcal{N}(0, 1)$ since the beginning, while for the input “car”, $(\mathcal{D}^l)_l$ first becomes closer to $\mathcal{N}(0, 1)$, then drifts away. It is also the case for $\phi_\theta$ with $\theta = 10$: for the input “deer”, the distance remains high, while it starts low and increases very slowly (input “dog”), or even decreases in the first place (input “truck”).

But, in general, when the curves have converged after 100 layers, it seems that the limit is the same whatever the starting data point, and depends only on the choice of the activation function.

**With individual input normalization.** In this setup, we perform the individual normalization of the inputs. Without loss of generality, we define it for inputs that are order-1 tensors as follows:

$$ 
\begin{aligned}
\hat{x}_{a;i} &:= \frac{x_{a;i} - \mu}{\sigma}, \\
\mu &:= \frac{1}{p} \sum_{i=1}^{p} x_i \\
\sigma^2 &:= \frac{1}{p-1} \sum_{i=1}^{p} (x_i - \mu)^2,
\end{aligned}
$$

where $\hat{x}$ is the normalized data point, $x_{a;i}$ is the i-th component of $\mathbf{x}_a \in \mathbb{R}^p$.

According to Figure 17, the curves seem to be approximately identical for each data point, contrary to the setup with normalization over the whole dataset. And, naturally, the distance between $\mathcal{D}^1$ and $\mathcal{N}(0, 1)$ is close to 0 when using the activation functions $\phi_\theta$, which confirms that our method to build Gaussian pre-activations is efficient at least in the beginning.

**Conclusion.** In the specific case of CIFAR-10 images, it seems that the trajectory of the sequence $(\mathcal{D}^l)_l$ depends on the inputs, but only to the extent that they have different norms. As it can be seen in Figure 17, images from different classes which have been normalized individually lead to the same trajectories of distribution of pre-activations.

## G.6 Experimental details of the training procedure

**Training, validation, and test sets.** For MNIST and CIFAR-10, we split randomly the initial training set into two sets: the training set, which will be actually used to train the neural network, and the validation set, which will be used to stop training when the network begins to overfit.

The sizes of the different sets are as follows:

* MNIST: 50000 training samples; 10000 validation samples; 10000 test samples;
<figure><img src="image_28.png" /><figcaption>Figure 16: Setup with data points normalized over the whole dataset.</figcaption></figure>

Propagation of the distribution $\mathcal{D}^l$ of the pre-activations across the layers $l \in [1, 100]$ when inputting various data points of CIFAR-10. Each curve represents the $\mathcal{L}^\infty$ distance between the CDF of $\mathcal{D}^l$ and the CDF of the Gaussian $N(0, 1)$. The dashed green line is the threshold of rejection of the Kolmogorov-Smirnov test with $p$-value 0.05: if a distribution $\mathcal{D}^l$ is represented by a point above this threshold, then the hypothesis "$\mathcal{D}^l = N(0, 1)$" is rejected with $p$-value 0.05.

• CIFAR-10: 42000 training samples, 8000 validation samples; 10000 test samples.

The training sets are split into mini-batches with 200 samples each. No data augmentation is performed.
<figure><img src="image_29.png" /><figcaption>Figure 17: Setup with data points **normalized individually**. The rest of the setup is identical to the one used in Figure 16.</figcaption></figure>

**Loss.** Given a classification task with $P$ classes, let $\mathbf{z}^L \in \mathbb{R}^P$ be the pre-activation outputted by the last layer of the neural network. First, we perform a softmax operation:

$$y_p := \mathrm{softmax}(z_p^L) = \frac{\exp(z_p^L)}{\sum_{p'=1}^{P} \exp(z_{p'}^L)}.$$

where the $(y_p)_p$ are the components of $\mathbf{y} \in \mathbb{R}^P$ and $(z_p^L)_p$ are the components of $\mathbf{z}^L$. Then, we compute the negative log-likelihood loss. For a target class $p \in \{1, \dots, P\}$, we pose:

$$\ell(\mathbf{y}, p) := -\log(y_p).$$

**Optimizer.** We use the Adam optimizer (Kingma & Ba, 2015) with the parameters:

* learning rate: 0.001;
• $\beta_1 = 0.9$;

• $\beta_2 = 0.999$;

• weight decay: 0.

**Learning rate scheduler.** We use a learning rate scheduler based on the reduction of the training loss. If the training loss does not decrease at least by a factor 0.01 for 10 epochs, then the learning rate is multiplied by a factor $1/\sqrt{10}$.

**Early stopping.** We add an early stopping rule based on the reduction of the validation loss. If the validation loss does not decrease at least by a factor of 0.001 for 30 epochs, then we stop training.