Dying ReLU and Initialization: Theory and Numerical Examples

TL;DR

This paper provides a theoretical analysis of the dying ReLU problem in deep neural networks, proves its inevitability with certain initializations, and introduces a new asymmetric initialization method that effectively prevents ReLU neurons from dying, supported by numerical evidence.

Contribution

The paper offers a rigorous proof of the dying ReLU phenomenon in deep networks and proposes a theoretically grounded asymmetric initialization to mitigate this issue.

Findings

Deep ReLU networks will almost surely have dying neurons as depth increases.

Symmetric initialization methods suffer from the dying ReLU problem.

The proposed asymmetric initialization effectively prevents neurons from dying, as demonstrated by numerical examples.

Abstract

The dying ReLU refers to the problem when ReLU neurons become inactive and only output 0 for any input. There are many empirical and heuristic explanations of why ReLU neurons die. However, little is known about its theoretical analysis. In this paper, we rigorously prove that a deep ReLU network will eventually die in probability as the depth goes to infinite. Several methods have been proposed to alleviate the dying ReLU. Perhaps, one of the simplest treatments is to modify the initialization procedure. One common way of initializing weights and biases uses symmetric probability distributions, which suffers from the dying ReLU. We thus propose a new initialization procedure, namely, a randomized asymmetric initialization. We prove that the new initialization can effectively prevent the dying ReLU. All parameters required for the new initialization are theoretically designed. Numerical examples are provided to demonstrate the effectiveness of the new initialization procedure.

Figures

Equations

$${\mathcal{N}}^{\ell }(\text{x})$$

$$\varphi (\text{x})=\text{ReLU}(\text{x}):=\left(\max \{x_1,0\},\cdots ,\max \{x_{\text{fan-in}},0\}\right),\text{where}\text{x}=(x_1,\cdots ,x_{\text{fan-in}}).$$

$$\mathcal{L}({\mathbf{\theta }}_{\vec{\mathbf{N}}},{\mathcal{T}}_m)=\frac{1}{m}\sum _{(\text{x},y)\in {\mathcal{T}}_m}\ell ({\mathcal{N}}^L(\text{x};{\mathbf{\theta }}_{\vec{\mathbf{N}}}),y),$$

$${\mathfrak{J}}_{{\vec{\mathbf{N}}}_L}:=\{\omega \in \Omega |{\mathcal{N}}^L(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_L}(\omega ))\text{ is a constant function in }\mathcal{D}\}.$$

$$P\left({\mathbf{W}}_j^{\ell }\in {\mathbb{R}}_{-}^{N_{\ell -1}},{\mathbf{b}}_j^{\ell }<0\right)\ge p>0,\forall 1\le j\le N_{\ell },\forall 1\le \ell \le L,$$

$$\begin{array}{r}\underset{L\to \infty }{\lim }P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=1,\end{array}$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})\le 1-\prod _{\ell =1}^{L-1}(1-(1/2{)}^{N_{\ell }}),$$

$$\begin{array}{r}\underset{L\to \infty }{\lim }P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=1,\underset{N\to \infty }{\lim }P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=0.\end{array}$$

$$P(⟨{\mathbf{W}}_j^{\ell },{\mathbf{v}}_1⟩>0,⟨{\mathbf{W}}_j^{\ell },{\mathbf{v}}_2⟩<0|{\mathbf{v}}_1,{\mathbf{v}}_2)\le \frac{1}{4},0\ne {\mathbf{v}}_1,{\mathbf{v}}_2\in {\mathbb{R}}_{+}^{N_{\ell -1}},\forall 1\le j\le N_{\ell },$$

$$p_{\text{low}}(L,N)\le P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})\le 1-\prod _{\ell =1}^{L-1}(1-(1/2{)}^N),$$

$$p_{\text{low}}(L,N)=1-a_1^{L-2}+\frac{(1-2^{-N+1})(1-2^{-N})}{1+(N-1)2^{-N}}(-a_1^{L-2}+a_2^{L-2}).$$

$${\mathcal{N}}^L(\text{x};{\mathbf{\theta }}^{\ast })={\mathbf{c}}^{\ast }=\underset{\mathbf{c}\in {\mathbb{R}}^{N_L}}{\mathrm{argmin}}\mathbb{E}\left[\ell (\mathbf{c},f(\text{x})))\right],$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})\le 1-(1-2^{-N}{)}^{L-1}\le 1-(1-(L-1)2^{-N})\le L{2}^{-N}=\delta .$$

$${\mathbf{v}}_{-k}={\left(v_1,\cdots ,v_{k-1},v_{k+1},\cdots ,v_{n+1}\right)}^T\in {\mathbb{R}}^n.$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})\le 1-\prod _{\ell =1}^{L-1}\left(1-{\left(1/2-{\gamma }_{\ell }\right)}^{N_{\ell }}\right),$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})\le 1-\prod _{\ell =1}^{L-1}\left(1-(1/\gamma {)}^N\right).$$

$$\begin{array}{r}\underset{L\to \infty }{\lim }P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=1,\underset{N\to \infty }{\lim }P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=0.\end{array}$$

$$\frac{{\mathcal{A}}_{low,\ell }}{2}\mathbb{E}[q^{\ell }(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell }})]+{\sigma }_{b,\ell }^2\le \mathbb{E}[q^{\ell +1}(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell +1}})]\le \frac{{\mathcal{A}}_{upp,\ell }}{2}\mathbb{E}[q^{\ell }(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell }})]+{\sigma }_{b,\ell }^2,$$

$${\mathcal{A}}_{upp,\ell }=\frac{{\sigma }_{b,\ell +1}^2}{{\sigma }_{b,\ell }^2}\left(\frac{N_{\ell -1}{\sigma }_w^2+2\sqrt{2/\pi }{N}_{\ell }{\mu }_{\ell ,1}^{\prime }{\sigma }_w+2N_{\ell }{\mu }_{\ell ,2}^{\prime }}{N_{\ell -1}+1}\right).$$

$$\frac{{\mathcal{A}}_{low,\ell }}{2}\mathbb{E}[q^{\ell }(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell }})]+{\sigma }_{b,\ell }^2\le \mathbb{E}[q^{\ell +1}(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell +1}})]\le \frac{{\mathcal{A}}_{upp,\ell }}{2}\mathbb{E}[q^{\ell }(\text{x};{\mathbf{\theta }}_{{\vec{\mathbf{N}}}_{\ell }})]+{\sigma }_{b,\ell }^2,$$

$${\sigma }_w=\sqrt{2}\left(-\frac{{\mu }_1^{\prime }}{\sqrt{\pi }}+\sqrt{\frac{{\mu }_1^{\prime 2}}{\pi }+1-{\mu }_2^{\prime }}\right),$$

$$\begin{array}{r}{f}_1(x)=|x|,f_2(x)=x\sin (5x),f_3(x)=1_{\{x>0\}}(x)+0.2\sin (5x).\end{array}$$

$$f_4(x_1,x_2)=\left[\begin{array}{c}|x_1+x_2|\\ |x_1-x_2|\end{array}\right].$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=P(\exists \ell \in \{1,\dots ,L-1\}\text{ such that }\varphi ({\mathcal{N}}^{\ell }(\text{x}))=0\forall \text{x}\in \mathcal{D}),$$

$$({\mathcal{N}}^{L-1}(\text{x}){)}_j={\mathbf{W}}_j^{L-1}\varphi ({\mathcal{N}}^{L-2}(\text{x}))+{\mathbf{b}}_j^{L-1}={\mathbf{W}}_j^{L-1}\varphi ({\mathcal{N}}^{L-2}(0))+{\mathbf{b}}_j^{L-1}=({\mathcal{N}}^{L-1}(0){)}_j>0.$$

$$G_{L-1}:=\{{\mathbf{W}}_j^{L-1}\varphi ({\mathcal{N}}^{L-2}(\text{x}))={\mathbf{W}}_j^{L-1}\varphi ({\mathcal{N}}^{L-2}(0)),\forall \text{x}\in \mathcal{D}\},$$

$$R_{L-2}:=\{\varphi ({\mathcal{N}}^{L-2}(\mathbf{x}))=\varphi ({\mathcal{N}}^{L-2}(0)),\forall \text{x}\in \mathcal{D}\}.$$

$$P(G_{L-1})=P(G_{L-1}|R_{L-2})P(R_{L-2})+P(G_{L-1}|R_{L-2}^c)P(R_{L-2}^c)=P(R_{L-2}).$$

$$P({\mathfrak{J}}_{{\vec{\mathbf{N}}}_L})=P(\exists \ell \in \{1,\dots ,L-1\}\text{ such that }\varphi ({\mathcal{N}}^{\ell }(\text{x})))=0\forall \text{x}\in \mathcal{D}).$$

$$\begin{array}{rl}{A}_{\ell }& :={\mathfrak{J}}_{{\vec{\mathbf{N}}}_{\ell }}\stackrel{a.s.}{=}\{\exists j\in \{1,\dots ,\ell -1\}\text{ such that }\varphi ({\mathcal{N}}^j(\text{x}))=0\forall \text{x}\in \mathcal{D}\},\end{array}$$

$$P(A_{\ell })\ge P(C_{\ell -1,0}).$$

Full text

\ams

60J05, 62M45, 68U99

1 Introduction

The rectified linear unit (ReLU), $\max\{x,0\}$, is one of the most successful and widely-used activation functions in deep learning [30, 40, 35]. The success of ReLU is based on its superior training performance [16, 47] over other activation functions such as the logistic sigmoid and the hyperbolic tangent [15, 31]. The ReLU has been used in various applications including image classification [29, 49], natural language processes [33], speech recognition [23], and game intelligence [44], to name a few.

The use of gradient-based optimization is inevitable in training deep neural networks. It has been widely known that the deeper a neural network is, the harder it is to train [46, 9]. A fundamental difficulty in training deep neural networks is the vanishing and exploding gradient problem [39, 17, 6]. The dying ReLU is a kind of vanishing gradient, which refers to a problem when ReLU neurons become inactive and only output 0 for any input. It has been known as one of the obstacles in training deep feed-forward ReLU neural networks [50, 1]. To overcome this problem, a number of methods have been proposed. Broadly speaking, these can be categorized into three general approaches. One approach modifies the network architectures. This includes but not limited to the changes in the number of layers, the number of neurons, network connections, and activation functions. In particular, many activation functions have been proposed to replace the ReLU [33, 20, 7, 28]. However, the performance of other activation functions varies on different tasks and data sets [40] and it typically requires a parameter to be turned. Thus, the ReLU remains one of the popular activation functions due to its simplicity and reliability. Another approach introduces additional training steps. This includes several normalization techniques [24, 42, 51, 4, 53] and dropout [45]. One of the most successful normalization techniques is the batch normalization [24]. It is a technique that inserts layers into the deep neural network that transform the output for the batch to be zero mean unit variance. However, batch normalization increases by 30% the computational overhead to each iteration [34]. The third approach modifies only weights and biases initialization procedure without changing any network architectures or introducing additional training steps [31, 15, 20, 43, 34]. The third approach is the topic of our work presented in this paper.

The intriguing ability of gradient-based optimization is perhaps one of the major contributors to the success of deep learning. Training deep neural networks using gradient-based optimization fall into the noncovex nonsmooth optimization. Since a gradient-based method is either a first- or a second-order method, and once converged, the optimizer is either a local minimum or a saddle point. The authors of [12] proved that the existence of local minima poses a serious problem in training neural networks. Many researchers have been putting immense efforts to mathematically understand the gradient method and its ability to solve nonconvex nonsmooth problems. Under various assumptions, especially on the landscape, many results claim that the gradient method can find a global minimum, can escape saddle points, and can avoid spurious local minima [32, 2, 13, 14, 57, 52, 55, 8, 10, 9, 25]. However, these assumptions do not always hold and are provably false for deep neural networks [41, 26, 3]. This further limits our understanding on what contributes to the success of the deep neural networks. Often, theoretical conditions are impossible to be met in practice.

Where to start the optimization process plays a critical role in training and has a significant effect on the trained result [36]. This paper focuses on a particular kind of bad local minima due to a bad initialization. Such a bad local minimum causes the dying ReLU. Specifically, we consider the worst case of dying ReLU, where the entire network dies, i.e., the network becomes a constant function. We refer this as the dying ReLU neural networks (NNs). This phenomenon could be well illustrated by a simple example. Suppose $f(x)=|x|$ is a target function we want to approximate using a ReLU network. Since $|x|=\text{ReLU}(x)+\text{ReLU}(-x)$, a 2-layer ReLU network of width 2 can exactly represent $|x|$. However, when we train a deep ReLU network, we frequently observe that the network is collapsed. This trained result is shown in Fig. 1. Our 1,000 independent simulations show that there is a high probability (more than 90%) for the deep ReLU network to collapse to a constant function. In this example, we employ a 10-layer ReLU network of width 2 which should perfectly recover $f(x)=|x|$.

Almost all common initialization schemes in training deep neural networks use symmetric probability distributions around 0. For example, zero mean uniform distributions and zero mean normal distributions were proposed and used in [31, 15, 20]. We show that when weights and biases are initialized from symmetric probability distributions around 0, the dying ReLU NNs occurs in probability as the number of depth goes to infinite. To the best of our knowledge, this is the first theoretical work on the dying ReLU. This result partially explains why training extremely deep networks is challenging. Furthermore, it says that the dying ReLU is inevitable as long as the network is deep enough. Also, our result implies that it is the network architecture that decides whether an initialization procedure is good or bad. Our analysis reveals that a specific network architecture can avoid the dying ReLU NNs with high probability. That is, for any $\delta>0$, when a symmetric initialization is used and $L=\Omega(\log_{2}(N/\delta))$ is satisfied where $L$ is the number of depth and $N$ is the number of width at each layer, with probability $1-\delta$, the dying ReLU NNs will not happen. This explains why wide networks are popularly used in practice. However, the use of deep and wide networks requires one to train a gigantic number of parameters and this could be computationally very expensive and time consuming.

On the other hands, deep and narrow networks play fundamental roles both in theory and in practice. For example, in the work of [19], it was shown that deep and narrow ReLU networks are inevitable in constructing finite element’s basis functions. This represents an application of deep ReLU networks in finite element method for solving PDEs. Also, many theoretical works [38, 54, 18] on the expressivity of ReLU networks rely heavily on deep and narrow networks in approximating polynomials through sparse concatenations. Despite the importance of deep and narrow networks, due to aforementioned reasons, they have not been successfully employed in practice. Thereby, we propose an initialization scheme, namely, a randomized asymmetric initialization (RAI), to directly overcome the dying ReLU. On the theoretical side, we show that RAI has a smaller upper bound of the probability of the dying ReLU NNs than standard symmetric initialization such as the He initialization [20]. All hyper-parameters for RAI are theoretically driven to avoid the exploding gradient problem. This is done by the second moment analysis [20, 17, 39].

The rest of the paper is organized as follows. After setting up notation and terminology in Section 2, we present the main theoretical results in Section 3. In Section 4, upon introducing a randomized asymmetric initialization, we discuss its theoretical properties. Numerical examples are provided in Section 5, before the conclusion in Section 6.

2 Mathematical setup

Let $\mathcal{N}^{L}:\mathbb{R}^{d_{\text{in}}}\to\mathbb{R}^{d_{\text{out}}}$ be a feed-forward neural network with $L$ layers and $N_{\ell}$ neurons in the $\ell$-th layer ($N_{0}=d_{\text{in}}$, $N_{L}=d_{\text{out}}$). Let us denote the weight matrix and bias vector in the $\ell$-th layer by $\bm{W}^{\ell}\in\mathbf{R}^{N_{\ell}\times N_{\ell-1}}$ and $\mathbf{b}^{\ell}\in\mathbb{R}^{N_{\ell}}$, respectively. Given an activation function $\phi$ which is applied element-wise, the feed-forward neural network is recursively defined as follows: $\mathcal{N}^{1}(\textbf{x})=\bm{W}^{1}\textbf{x}+\bm{b}^{1}$ and

[TABLE]

Here $\mathcal{N}^{L}$ is called a $L$-layer neural network or a $(L-1)$-hidden layer neural network. In this paper, the rectified linear unit (ReLU) activation function is employed, i.e.,

[TABLE]

Let $\vec{\bm{N}}=(N_{0},N_{1},\cdots,N_{L})$ be a vector representing the network architecture. Let $\bm{\theta}_{\vec{\bm{N}}}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ be the set of all weight matrices and bias vectors. Let $\mathcal{T}_{m}=\{\textbf{x}_{i},y_{i}\}_{1\leq i\leq m}$ be the set of $m$ training data and let $\mathcal{D}=\{\textbf{x}_{i}\}_{i=1}^{m}$ be the training input data. We assume that $\mathcal{D}\subset B_{r}(0)$ for some $r>0$. Given $\mathcal{T}_{m}$, in order to train $\bm{\theta}_{\vec{\bm{N}}}$, we consider the standard loss function $\mathcal{L}(\bm{\theta}_{\vec{\bm{N}}},\mathcal{T}_{m})$:

[TABLE]

where $\ell:\mathbb{R}^{d_{\text{out}}}\times\mathbb{R}^{d_{\text{out}}}\to\mathbb{R}$ is a loss criterion. In training neural networks, the gradient-based optimization is typically employed to minimize the loss $\mathcal{L}$. The first step for training would be to initialize weight matrices and bias vectors. Typically, they are initialized according to certain probability distributions. For example, uniform distributions around 0 or zero-mean normal distributions are common choices.

In this paper, we focus on the worst case of dying ReLU, where the entire network dies, i.e., the network becomes a constant function. We refer this as the dying ReLU neural network. We then define two phases: (1) a network is dead before training, and (2) a network is dead after training. The phase 1 implies the phase 2, but not vice versa. When the phase 1 happens, we say the network is born dead (BD).

3 Dying probability bounds and asymptotic properties

We investigate the probability of the dying ReLU neural networks at random initialization. Let $(\Omega,\mathcal{F},P)$ be a probability space. For a network architecture $\vec{\bm{N}}_{L}=(d_{\text{in}},N_{1},\cdots,N_{L-1},d_{\text{out}})$, let $\bm{\theta}_{\vec{\bm{N}}_{L}}:\Omega\to\mathbb{R}^{|\vec{\bm{N}}_{L}|}$ be a vector-valued random variable, where $|\vec{\bm{N}}_{L}|:=\sum_{\ell=1}^{L}N_{\ell}(N_{\ell-1}+1)$ is the total number of network parameters. Let us consider the following BD event:

[TABLE]

Then, $P(\mathfrak{J}_{\vec{\bm{N}}_{L}})$ is the BD probability (BDP). The following theorem shows that a deep ReLU network will eventually be BD in probability as the number of depth $L$ goes to infinity.

Theorem 3.1**.**

For $\vec{\bm{N}}_{L}=(d_{\text{in}},N,\cdots,N,d_{\text{out}})\in\mathbb{N}^{L+1}$, let $\bm{\theta}_{\vec{\bm{N}}_{L}}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ and $\mathcal{N}^{L}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{L}})$ be a ReLU neural network with $L$ layers, each having $N$ neurons. Suppose that all weights and biases are randomly initialized from probability distributions, which satisfy

[TABLE]

for some constant $p>0$, where $\bm{W}^{\ell}_{j}$ is the $j$-th row of the $\ell$-th layer weight matrix and $\bm{b}^{\ell}_{j}$ is the $j$-th component of the $\ell$-th layer bias vector. Then

[TABLE]

where $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ is defined in (3).

Proof 3.2**.**

The proof can be found in Appendix A.

We remark that Equation 4 is a very mild condition and it can be satisfied in many cases. For example, when symmetric probability distributions around 0 are employed, the condition is met with $p=2^{-N-1}$. Theorem 3.1 implies that the fully connected ReLU network will be dead at the initialization as long as the network is deep enough. This explains theoretically why training a very deep network is hard.

Theorem 3.1 shows that the ReLU network asymptotically will be dead. Thus, we are now concerned with the convergence behavior of the probability of NNs being BD. Since almost all common initialization procedures use symmetric probability distributions around 0, we derive an upper bound of the born dead probability (BDP) for symmetric initialization.

Theorem 3.3**.**

*For $\vec{\bm{N}}_{L}=(d_{\text{in}},N_{1},\cdots,N_{L-1},d_{\text{out}})\in\mathbb{N}^{L+1}$, let $\bm{\theta}_{\vec{\bm{N}}_{L}}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ and $\mathcal{N}^{L}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{L}})$ be a ReLU neural network with $L$ layers, each having $N_{1},\cdots,N_{L-1}$ neurons. Suppose that all weights are independently initialized from symmetric probability distributions around 0 and all biases are either drawn from a symmetric distribution or set to zero. Then *

[TABLE]

where $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ is defined in (3). Furthermore, if $\vec{\bm{N}}_{L}=(d_{\text{in}},N,\cdots,N,d_{\text{out}})\in\mathbb{N}^{L+1}$,

[TABLE]

Proof 3.4**.**

The proof can be found in Appendix B.

Theorem 3.3 provides an upper bound of the BDP. It shows that at a fixed depth $L$, the network will not be BD in probability as the number of width $N$ goes to infinite. In order to understand how this probability behaves with respect to the number of width and depth, a lower bound is needed. We thus provide a lower bound of the BDP of ReLU NNs at $d_{\text{in}}=1$.

Theorem 3.5**.**

For a network architecture $\vec{\bm{N}}_{L}=(1,N,\cdots,N,d_{\text{out}})\in\mathbb{N}^{L+1}$, let $\bm{\theta}_{\vec{\bm{N}}_{L}}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ and $\mathcal{N}^{L}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{L}})$ be a ReLU neural network with $L$ layers, each having $N$ neurons. Suppose that all weights are independently initialized from continuous symmetric probability distributions around 0, which satisfies

[TABLE]

where $\bm{W}_{j}^{\ell}$ is the $j$-th row of the $\ell$-th layer weight matrix, and all biases are initialized to zero. Then

[TABLE]

where $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ is defined in (3), $a_{1}=1-(1/2)^{N}$, $a_{2}=1-(1/2)^{N-1}-(N-1)(1/4)^{N}$, and

[TABLE]

Proof 3.6**.**

The proof can be found in Appendix C.

Theorem 3.5 reveals that the BDP behavior depends on the network architecture. In Fig. 2, we plot the BDP with respect to increasing the number of layers at varying width from $N=2$ to $N=5$. A bias-free ReLU feed-forward NN with $d_{\text{in}}=1$ is employed with weights randomly initialized from symmetric distributions. The results of one million independent simulations are used to calculate each probability estimation. Numerical estimations are shown as symbols. The upper and lower bounds from Theorem 3.5 are also plotted with dash and dash-dot lines, respectively. We see that when the NN gets narrower, the probability of NN being BD grows faster as the depth increases. Also, at a fixed width $N$, the BDP grows as the number of layer increases. This is expected by Theorems 3.1 and 3.5.

Once the network is BD, we have no hope to train the network successfully. Here we provide a formal statement of the consequence of the network being BD.

Theorem 3.7**.**

For any $\vec{\bm{N}}_{L}$, suppose $\bm{\theta}_{\vec{\bm{N}}_{L}}$ is initialized to be in $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ defined in (3). Then, for any loss function $\mathcal{L}$, and for any gradient based method, the ReLU network is optimized to be a constant function, which minimizes the loss.

Proof 3.8**.**

The proof can be found in Appendix D.

Theorem 3.7 implies that no matter what gradient-based optimizers are employed including stochastic gradient desecent (SGD), SGD-Nesterov [48], AdaGrad [11], AdaDelta [56], RMSProp [22], Adam [27], BFGS [37], L-BFGS [5], the network is trained to be a constant function which minimizes the loss.

If the online-learning or the stochastic gradient method is employed, where the training data are independently drawn from a probability distribution $P_{\mathcal{D}}$, the optimized network is

[TABLE]

where the expectation $\mathbb{E}$ is taken with respect to $\textbf{x}\sim P_{\mathcal{D}}$. For example, if $L^{2}$-loss is employed, i.e., $\ell(\mathcal{N}^{L}(\textbf{x}),f(\textbf{x}))=(\mathcal{N}^{L}(\textbf{x})-f(\textbf{x}))^{2}$, the resulting network is $\mathbb{E}[f(\textbf{x})]$. If $L^{1}$ loss is employed, i.e., $\ell(\mathcal{N}^{L}(\textbf{x}),f(\textbf{x}))=|\mathcal{N}^{L}(\textbf{x})-f(\textbf{x}))|$, the resulting network is the median of $f(\textbf{x})$ with respect to $\textbf{x}\sim P_{\mathcal{D}}$. Note that the mean absolute error (MAE) and the mean squared error (MSE) used in practice are discrete versions of $L^{1}$ and $L^{2}$ loss, respectively, if the size of minibatch is large.

When we design a neural network, we want the BDP to be small, say, less than 1% or 10%. Then, the upper bound (Equation 5) of Theorem 3.3 can be used for designing a specific network architecture, which has a small probability of NNs being born dead.

Corollary 3.9**.**

*For fixed depth $L$, let $\vec{\bm{N}}_{L}=(d_{\text{in}},N,\cdots,N,d_{\text{out}})\in\mathbb{N}^{L+1}$. For $\delta>0$, suppose the width $N$ is $N=\log_{2}L/\delta$. Then, $P(\mathfrak{J}_{\vec{\bm{N}}_{L}}^{c})>1-\delta$, where $\mathfrak{J}_{\vec{\bm{N}}_{L}}^{c}$ is the complement of $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ defined in (3). *

Proof 3.10**.**

This readily follows from

[TABLE]

As a practical guide, we constructed a diagram shown in Fig. 3 that includes both theoretical predictions and our numerical tests. We see that as the number of layers increases, the numerical tests match closer to the theoretical results. It is clear from the diagram that a 10-layer NN of width 10 has a probability of dying less than 1% whereas a 10-layer NN of width 5 has a probability of dying greater than 10%; for width of three the probability is about 60%. Note that the growth rate of the maximum number of layers is exponential which is expected by Corollary 3.9.

Our results explain why deep and narrow networks have not been popularly employed in practice and reveal a limitation of training extremely deep networks. We remark that modern practical networks happen to be sufficiently wide enough to have small BDPs. We believe that the born dead network problem is a major impediment for one to explore a potential use of deep and narrow networks.

4 Randomized asymmetric initialization

The so-called ‘He initialization’ [20] is perhaps one of the most popular initialization schemes in deep learning community, especially when the ReLU activation function is concerned. The effectiveness of the He initialization has been shown in many machine learning applications. The He initialization uses mean zero normal distributions. Thus, as we discussed earlier, it suffers from the dying ReLU. We thus propose a new initialization procedure, namely, a randomized asymmetric initialization. The motivation is in twofolds. One is to mimics the He initialization so that the new scheme can produce similar generalization performance. The other is to alleviate the problem of dying ReLU neural networks. We hope that the proposed initialization opens a new way to explore the potential use of deep and narrow networks in practice.

For ease of discussion, we introduce some notation. For any vector $\bm{v}\in\mathbb{R}^{n+1}$ and $k\in\{1,\cdots,n+1\}$, we define

[TABLE]

In order to train a $L$-layer neural network, we need to initialize $\bm{\theta}_{L}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$. At each layer, let $\textbf{V}^{\ell}=[\bm{W}^{\ell},\bm{b}^{\ell}]\in\mathbb{R}^{N_{\ell}\times(N_{\ell-1}+1)}$. We denote the $j$-th row of $\textbf{V}^{\ell}$ by $\textbf{V}^{\ell}_{j}=[\bm{W}^{\ell}_{j},\bm{b}^{\ell}_{j}]\in\mathbb{R}^{N_{\ell-1}+1}$, $j=1,\cdots,N_{\ell}$ where $N_{0}=d_{\text{in}}$ and $N_{L}=d_{\text{out}}$.

4.1 Proposed initialization

We propose to initialize $\textbf{V}^{\ell}$ as follows. Let $\text{P}_{\ell}$ be a probability distribution defined on $[0,M_{\ell}]$ for some $M_{\ell}>0$ or $[0,\infty)$. Note that $\text{P}_{\ell}$ is asymmetric around 0. At the first layer of $\ell=1$, we employ the ‘He initialization’ [20], i.e., $\textbf{W}^{1}_{ij}\sim N(0,2/d_{\text{in}})$ and $\textbf{b}^{1}=\textbf{0}$. For $\ell\geq 2$, and each $1\leq j\leq N_{\ell}$, we initialize $\bm{V}_{j}^{\ell}$ as follows:

Randomly choose $k^{\ell}_{j}$ in $\{1,2,\cdots,N_{\ell-1},N_{\ell-1}+1\}$. 2. 2.

Initialize $(\textbf{V}^{\ell}_{j})_{-k^{\ell}_{j}}\sim\mathcal{N}(0,\sigma_{\ell}^{2}\bm{I})$ and $(\textbf{V}^{\ell}_{j})_{k^{\ell}_{j}}\sim\text{P}_{\ell}$.

Since an index is randomly chosen at each $\ell$ and $j$ and a positive number is randomly drawn from an asymmetric probability distribution around 0, we name this new initialization a randomized asymmetric initialization. Only for the first layer, the He initialization is employed. This is because since an input could have a negative value, if the weight which corresponds to the negative input were to be initialized from $\text{P}_{\ell}$, this could cause the dying ReLU. We note that the new initialization requires us to choose $\sigma_{\ell}^{2}$ and $\text{P}_{\ell}$. In Subsection 4.2, these will be theoretically determined. One could choose multiple indices in the step 1 of the new initialization. However, for simplicity, we constraint ourselves to a single index case.

We first show that this new initialization procedure results in a smaller upper bound of the BDP.

Theorem 4.1**.**

For $\vec{\bm{N}}_{L}=(d_{\text{in}},N_{1},\cdots,N_{L-1},d_{\text{out}})\in\mathbb{N}^{L+1}$, let $\bm{\theta}_{\vec{\bm{N}}_{L}}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ and $\mathcal{N}^{L}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{L}})$ be a ReLU neural network with $L$ layers, each having $N_{1},\cdots,N_{L-1}$ neurons. Suppose $\bm{\theta}_{\vec{\bm{N}}_{L}}$ is initialized by the randomized asymmetric initialization (RAI). Then,

[TABLE]

where $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ is defined in (3), $\gamma_{1}=0$ and $\gamma_{j}$’s are some constants in $(0,0.5]$, which depend on $\{N_{\ell}\}_{\ell=1}^{L-1}$ and the training input data $\mathcal{D}$.

Proof 4.2**.**

The proof can be found in Appendix E.

When a symmetric initialization is employed, $\gamma_{j}=0$ for all $1\leq j<N_{L}$, which results in Equation 5 of Theorem 3.3. Although the new initialization has a smaller upper bound compared to those by symmetric initialization, as Theorem 3.1 suggests, it also asymptotically suffers from the dying ReLU.

Corollary 4.3**.**

Assuming the same conditions in Theorem 4.1, and $N_{\ell}=N$ for all $\ell$. Let $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ be the event defined in (3). Then, there exists $2<\gamma$, which depends on $N,L$ and the training input data $\mathcal{D}$, such that

[TABLE]

*For fixed depth $L$ and $\delta>0$, if the width $N$ is $N=\log_{\gamma}L/\delta$, we have $P(\mathfrak{J}_{\vec{\bm{N}}_{L}}^{c})>1-\delta$, where $\mathfrak{J}_{\vec{\bm{N}}_{L}}^{c}$ is the complement of $\mathfrak{J}_{\vec{\bm{N}}_{L}}$. Furthermore, *

[TABLE]

Proof 4.4**.**

The proof is readily followed from Theorem 3.1, 4.1 and Corollary 3.9.

4.2 Second moment analysis

The proposed randomized asymmetric initialization (RAI) described in Subsection 4.1 requires us to determine $\sigma_{\ell}^{2}$ and $\text{P}_{\ell}$. Similar to the He initialization ([21]), we aim to properly choose initialization parameters from the length map analysis. Following the work of [39], we present the analysis of a single input propagation through the deep ReLU network. To be more precise, for a network architecture $\vec{\bm{N}}_{\ell}=(d_{\text{in}},N_{1},\cdots,N_{\ell-1},N_{\ell})\in\mathbb{N}^{\ell+1}$, we track the expectation of the normalized squared length of the input vector at each layer, $\mathbb{E}[q^{\ell}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{\ell}})]$, where $q^{\ell}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{\ell}})=\frac{\|\mathcal{N}^{\ell}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{\ell}})\|^{2}}{N_{\ell}}$. The expectation $\mathbb{E}$ is taken with respect to all weights and biases $\bm{\theta}_{\vec{\bm{N}}_{\ell}}$.

Theorem 4.5**.**

For $\vec{\bm{N}}_{\ell}=(d_{\text{in}},N_{1},\cdots,N_{\ell-1},N_{\ell})\in\mathbb{N}^{\ell+1}$, let $\bm{\theta}_{\vec{\bm{N}}_{\ell}}=\{\bm{W}^{j},\bm{b}^{j}\}_{1\leq j\leq\ell}$ be initialized by the RAI with the following $\sigma_{\ell}^{2}$ and $\text{P}_{\ell}$. Let $\sigma_{\ell}^{2}:=\frac{\sigma_{w}^{2}}{N_{\ell-1}}$ for some $\sigma_{w}^{2}$ and let $\text{P}_{\ell}$ be a probability distribution whose support is $[0,M_{\ell}]\subset\mathbb{R}^{+}$ such that $\mu_{\ell,i}^{\prime}=E[X_{\ell}^{i}]<\infty$, for $i=1,2$, where $X_{\ell}\sim\text{P}_{\ell}$ and $\mu_{\ell,1}^{\prime}\geq M_{\ell}/2$. Then for any input $\textbf{x}\in\mathbb{R}^{d_{in}}$, we have

[TABLE]

where $\sigma_{b,\ell}^{2}=\frac{\mu_{\ell+1,2}^{\prime}+\sigma_{\ell+1}^{2}N_{\ell+1}}{N_{\ell+1}+1}$, $\mathcal{A}_{low,\ell}=\frac{\sigma_{b,\ell+1}^{2}}{\sigma_{b,\ell}^{2}}\left(\frac{N_{\ell}\mu_{\ell,2}^{\prime}+N_{\ell-1}\sigma_{w}^{2}}{N_{\ell-1}+1}\right)$, and

[TABLE]

Proof 4.6**.**

The proof can be found in Appendix F.

Corollary 4.7**.**

Under the same conditions of Theorem 4.5, if $N_{\ell}=N$, $M_{\ell}=M$, $E[X_{\ell}^{i}]=\mu_{i}^{\prime}$, $i=1,2$ for all $\ell$, and $\mu_{1}^{\prime}\geq M/2$, we have

[TABLE]

*where $\sigma_{b,\ell}^{2}=\frac{\mu_{2}^{\prime}+\sigma_{w}^{2}}{N+1}$, $\mathcal{A}_{low,\ell}=\frac{N(\mu_{2}^{\prime}+\sigma_{w}^{2})}{N+1}$, and $\mathcal{A}_{upp,\ell}=\frac{N(\sigma_{w}^{2}+2\sqrt{2/\pi}\mu_{1}^{\prime}\sigma_{w}+2\mu_{2}^{\prime})}{N+1}$. *

Since $\sigma_{b,\ell}^{2}>0$, $\lim_{\ell\to\infty}\mathbb{E}[q^{\ell}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{\ell}})]$ cannot be zero. In order for $\lim_{\ell\to\infty}\mathbb{E}[q^{\ell}(\textbf{x};\bm{\theta}_{\vec{\bm{N}}_{\ell}})]<\infty$, the initialization parameters $(\sigma_{w}^{2},\mu_{\ell,1}^{\prime},\mu_{\ell,2}^{\prime})$ have to be chosen to satisfy $\mathcal{A}_{upp,\ell}<2$. Assuming $\mu_{2}^{\prime}<1$, if $\sigma_{w}$ is chosen to be

[TABLE]

we have $\frac{\mathcal{A}_{upp,\ell}}{2}=\frac{N}{N+1}$ which satisfies the condition.

4.3 Comparison against other initialization procedures

In Fig. 4, we demonstrate the probability that the network is BD by the proposed randomized asymmetric initialization (RAI) method. Here we employ $\text{P}=\text{Beta}(2,1)$ and $\sigma_{w}=-\frac{2\sqrt{2}}{3\sqrt{\pi}}+\sqrt{1+\frac{8}{9\pi}}\approx 0.6007$ from Equation 7 (see Appendix G for the Python code). To compare against other procedures, we present the results by the He initialization [20]. We also present the results of existing asymmetric initialization procedures; the orthogonal [43] and the layer-sequential unit-variance (LSUV) [34] initializations. The LSUV is the orthogonal initialization combined with rescaling of weights such that the output of each layer has unit variance. Because weight rescaling cannot make the output escape from the negative part of ReLU, it is sufficient to consider the orthogonal initialization. We see that the BDPs by the orthogonal initialization are very close to and a little lower than those by the He initialization. This implies that the orthogonal initialization cannot prevent the dying ReLU network. Furthermore, we show the results by the He initialization with positive constant bias of 10 and 100. Naively speaking, having a big positive bias will help in preventing dying ReLU neurons, as the input of each layer is pushed to be positive, although this might cause the exploding gradient problem in training. We see that the BDPs by the He with bias 10 and 100 are lower than those by the He with bias 0 and the orthogonal initialization. However, it is clearly observed that our proposed initialization (RAI) drastically drops the BDPs compared to all others. This is implied by Theorem 4.1.

As a practical guide, we constructed a diagram shown in Fig. 5 to demonstrate the 1% safe operation region by the RAI. It is clear from the diagram that the RAI drastically expands its 1% safe region compared to those by the He initialization. Note that the growth rate of the maximum number of layers is exponential which is expected by Corollary 4.3.

5 Numerical examples

We demonstrate the effectiveness of the proposed randomized asymmetric initialization (RAI) in training deep ReLU networks.

Test functions include one- and two-dimensional functions of different regularities. The following test functions are employed as unknown target functions. For one dimensional cases,

[TABLE]

For two dimensional case,

[TABLE]

We employ the network architecture having the width of $d_{\text{in}}+d_{\text{out}}$ at all layers. Here $d_{\text{in}}$ and $d_{\text{out}}$ are the dimensions of the input and output, respectively. It was shown in [18] that the minimum number of width required for the universal approximation is less than or equal to $d_{\text{in}}+d_{\text{out}}$. We thus choose this specific network architecture. as it theoretically guarantees to approximate any continuous function. In all numerical examples, we employ one of the most popular first-order gradient-based optimization, Adam [27] with the default parameters. The minibatch size is chosen to be either 64 or 128. The standard $L_{2}$-loss function $\ell(\mathcal{N}^{L}(\textbf{x},\bm{\theta}),(\textbf{x},y))=(\mathcal{N}^{L}(\textbf{x};\bm{\theta})-y)^{2}$ is used on 3,000 training data. The training data are randomly uniformly drawn from $[-\sqrt{3},\sqrt{3}]^{d_{\text{in}}}$. Without changing any setups described above, we present the approximation results based on different initialization procedures. The results by our proposed randomized asymmetric initialization are referred to ‘Rand. Asymmetric’ or ‘RAI’. Specifically, we use $\text{P}=\text{Beta}(2,1)$ with $\sigma_{w}$ defined in Equation 7. To compare against other methods, we also show the results by the He initialization [20]. We present the ensemble of 1,000 independent training simulations.

In one dimensional examples, we employ a 10-layer ReLU network of width 2. It follows from Fig. 4 that we expect to observe at least 88% training results by the symmetric initialization and 22% training results by the RAI are collapsed. In the two dimensional example, we employ a 20-layer ReLU network of width 4. According to Fig. 4, we expect to see at least 63% training results by the symmetric initialization and 3.7% training results by the RAI are collapsed.

Fig. 6 shows all training outcomes of our 1,000 simulations for approximating $f_{1}(x)$ and its corresponding empirical probabilities by different initialization schemes. For this specific test function, we observe only 3 trained results shown in A, B, C. In fact, $f_{1}(x)$ can be represented exactly by a 2-layer ReLU network with width 2, $f_{1}(x)=|x|=\text{ReLU}(x)+\text{ReLU}(-x)$. It can clearly be seen that the He initialization results in the collapse with probability more than 90%. However, this probability is drastically reduced to 40% by the RAI. These probabilities are different from the probability that the network is BD. This implies that even though the network wasn’t BD, there are cases that after training, the network dies. In this example, 5.6% and 18.3% of results by the symmetric and our method, respectively, are not dead at the initialization, however, they are ended up with collapsing after training. The 37.3% of training results by the RAI perfectly recover the target function $f_{1}(x)$, however, only 2.2% of results by the He initialization achieve this success. Also, 22.4% of the RAI and 4.2% of the He initialization produce the half-trained results which correspond to Fig. 6 (B). We remark that the only difference in training is the initialization. This implies that our new initialization does not only prevent the dying ReLU network but also improves the quality of the training in this case.

The approximation results for $f_{2}(x)$ are shown in Fig. 7. Note that $f_{2}$ is a $C^{\infty}$ function. It can be seen that 91.9% of training results by the symmetric initialization and 29.2% of training results by the RAI are collapsed which correspond to Fig. 7 (A). This indicates that the RAI can effectively alleviate the dying ReLU. In this example, 3.9% and 7.2% of results by the symmetric and our method, respectively, are not dead at the initialization, however, they are ended up with collapsing after training. Except for the collapse, other training outcomes are not easy to be classified. Fig. 7 (B,C,D) shows three training results among many others. We observe that the behavior and result of training are not easily predictable in general. However, we consistently observe partially collapsed results after training. Such partial collapses are also observed in Fig. 7 (B,C,D). We believe that this requires more attention and postpone the study of this partial collapse to future work.

Similar behavior is observed for approximating a discontinuous function $f_{3}(x)$. The approximation results for $f_{3}(x)$ and its corresponding empirical probabilities are shown in Fig. 8. We see that 93.8% of training results by the He initialization and 32.6% of training results by the RAI are collapsed which correspond to Fig. 8 (A). In this example, the RAI drops the probability of collapsing by 60.3 percentage point. Again, this implies that the RAI can effectively avoid the dying ReLU, especially when deep and narrow ReLU networks are employed. Fig. 8 (B,C,D) shows three trained results among many others. Again, we observe partially collapsed results.

Next we show the approximation result for a multi-dimensional inputs and outputs function $f_{4}(\textbf{x})$ defined in Equation 9. We observe similar behavior. Fig. 9 shows some of approximation results for $f_{4}$ and its corresponding probabilities. Among 1,000 independent simulations, the collapsed results are obtained by the He initialization with 76.8% probability and by the RAI with 9.6% probability. From Fig. 4, we expect to observe at least 63% and 3.7% of results by the symmetric and the RAI to be collapsed. Thus, in this example, 13.8% and 5.9% of results by the symmetric and our method, respectively, are not dead at the initialization, however, they are ended up with Fig. 9 (A) after training. This indicates that the RAI can also effectively overcome the dying ReLU in multi-dimensional inputs and outputs tasks.

As a last example, we demonstrate the performance of the RAI on the MNIST dataset. For the training, we employ the cross-entropy loss and the mini-batch size of 100. The networks are trained using Adam [27] with its default values. In Fig. 10, the convergence of the test accuracy is shown with respect to the number of epochs. On the left and right, we employ the ReLU network of depth 2 and width 1024 and of depth 50 and with 10, respectively. We see that when the shallow and wide network is employed, the RAI and the He initialization show similar generalization performance (test accuracy). However, when the deep and narrow network is employed, the RAI performs better than the He initialization. This indicates that the proposed RAI not only reduces the BDP of deep networks, but also has good generalization performance.

6 Conclusion

In this paper, we establish, to the best of our knowledge, the first theoretical analysis on the dying ReLU. By focusing on the worst case of dying ReLU, we define ‘the dying ReLU network’ which refers to the problem when the ReLU network is dead. We categorize the dying process into two phases. One phase is the event where the ReLU network is initialized to be a constant function. We refer to this event as ‘the network is born dead’. The other phase is the event where the ReLU network is collapsed after training. Certainly, the first phase implies the second, but not vice versa. We show that the probability that the network is born dead goes to 1 as the depth goes infinite. Also, we provide an upper and a lower bound of the dying probability in $d_{\text{in}}=1$ when the standard symmetric initialization is used.

Furthermore, in order to overcome the dying ReLU networks, we propose a new initialization procedure, namely, a randomized asymmetric initialization (RAI). We show that the RAI has a smaller upper bound of the probability of NNs being born dead. By establishing the expected length map relation (second moment analysis), all parameters needed for the new method are theoretically designed. Numerical examples are provided to demonstrate the performance of our method. We observe that the RAI does not only overcome the dying ReLU but also improves the training and generalization performance.

Acknowledgements

This work received support by the DARPA EQUiPS grant N66001-15-2-4055, the AFOSR grant FA9550-17-1-0013, and the DARPA AIRA grant HR00111990025. The research of the third author was partially supported by the NSF of China 11771083 and the NSF of Fujian 2017J01556, 2016J01013.

Appendix A Proof of Theorem 3.1

The proof starts with the following lemma.

Lemma A.1**.**

Let $\mathcal{N}^{L}(\textbf{x})$ be a $L$-layer ReLU neural network with $N_{\ell}$ neurons at the $\ell$-th layer. Suppose all weights are randomly independently generated from probability distributions satisfying $P(\bm{W}_{j}^{\ell}\bm{z}=\bm{0})=0$ for any nonzero vector $\bm{z}\in\mathbb{R}^{N_{\ell-1}}$ and any $j$-th row of $\bm{W}^{\ell}$. Then

[TABLE]

where $\mathcal{D}\subset B_{r}(\bm{0})=\{\textbf{x}\in\mathbb{R}^{d_{\text{in}}}|\|x\|<r\}$ for any $r>0$ and $\mathfrak{J}_{\vec{\bm{N}}_{L}}$ is the event defined in (3).

Proof A.2**.**

Suppose $\mathcal{N}^{L}(\textbf{x})=\mathcal{N}^{L}(\bm{0})$ for all $\textbf{x}\in\mathcal{D}\subset B_{r}(\bm{0})$. Then $\phi(\mathcal{N}^{L-1}(\textbf{x}))=\phi(\mathcal{N}^{L-1}(\bm{0}))$ for all $\textbf{x}\in\mathcal{D}$. If $\phi(\mathcal{N}^{L-1}(\textbf{x}))=\phi(\mathcal{N}^{L-1}(\bm{0}))=\bm{0}$, we are done as $\ell=L-1$. If it is not the case, there exists $j$ in $\{1,\cdots,N_{L-1}\}$ such that for all $\textbf{x}\in\mathcal{D}$,

[TABLE]

*Thus we have $\bm{W}^{L-1}_{j}\left(\phi(\mathcal{N}^{L-2}(\textbf{x}))-\phi(\mathcal{N}^{L-2}(\bm{0}))\right)=0$ for all $\textbf{x}\in\mathcal{D}$. Let consider the following events: *

[TABLE]

Note that $P(G_{L-1}|R_{L-2})=1$. Also, since $P(\bm{W}^{\ell}\bm{z})=0$ for any nonzero vector $\bm{z}$, we have $P(G_{L-1}|R_{L-2}^{c})=0$. Therefore,

[TABLE]

Thus we can focus on $\phi(\mathcal{N}^{L-2}(\textbf{x}))=\phi(\mathcal{N}^{L-2}(\bm{0})),\forall\textbf{x}\in\mathcal{D}$. If $\phi(\mathcal{N}^{L-2}(\textbf{x}))=\phi(\mathcal{N}^{L-2}(\bm{0}))=\bm{0}$, we are done as $\ell=L-2$. If it is not the case, it follows from the similar procedure that $\phi(\mathcal{N}^{L-3}(\textbf{x}))=\phi(\mathcal{N}^{L-3}(\bm{0}))$ in $\mathcal{D}$. By repeating these, we conclude that

[TABLE]

Proof A.3**.**

Let $\mathcal{D}\subset B_{r}(0)\subset\mathbb{R}^{d_{\text{in}}}$ be a training domain where $r$ is any positive real number. We consider a probability space $(\Omega,\mathcal{F},P)$ where all random weight matrices and bias vectors are defined on. For every $\ell\geq 1$, let $\mathcal{F}_{\ell}$ be a sub-$\sigma$-algebra of $\mathcal{F}$ generated by $\{\bm{W}^{j},\bm{b}^{j}\}_{1\leq j\leq\ell}$. Since $\mathcal{F}_{k}\subset\mathcal{F}_{\ell}$ for $k\leq\ell$, $(\mathcal{F}_{\ell})$ is a filtration. Let us define the events of our interest $\{A_{\ell}\}_{2\leq\ell}$ where

[TABLE]

where the second equality is from Lemma A.1. Note that $A_{\ell}$ is measurable in $\mathcal{F}_{\ell-1}$. Here $\{\bm{b}^{\ell}\}_{1\leq\ell}$ could be either [math] or random vectors. Since $\mathcal{N}^{1}(\textbf{x})=\bm{W}^{1}\textbf{x}+\bm{b}^{1}$, $P(A_{1})=0$. To calculate $P(A_{\ell})$ for $\ell\geq 2$, let consider another event $C_{\ell,k}$ on which exactly $(N_{\ell}-k)$-components of $\phi(\mathcal{N}^{\ell}(\textbf{x}))$ are zero on $\mathcal{D}$. For notational completeness, we set $C_{1,k}=\emptyset$ for $0\leq k<N_{1}$ and $C_{1,N_{1}}=\Omega$. Then since $C_{\ell-1,0}\subset A_{\ell}$, we have

[TABLE]

We want to show $\lim_{\ell\to\infty}P(C_{\ell,0})=1$. Since $\{C_{\ell-1,k}\}_{0\leq k\leq N_{\ell-1}}$ is a partition of $\Omega$, by the total law of probability, we have

[TABLE]

where $P(C_{\ell,0}|C_{\ell-1,0})=1$, and $P(C_{\ell,s}|C_{\ell-1,0})=0,\forall s\geq 1$. Since $\bm{W}_{ij}^{\ell}$ and $\bm{b}^{\ell}_{j}$ are independently initialized, we have

[TABLE]

where the second and the third inequalities hold from the assumption. Here $p_{k}>0$ does not depend on $\ell$. If $\bm{W}_{ij}^{\ell},\bm{b}_{j}^{\ell}$ are randomly initialized from symmetric distributions around 0,

[TABLE]

Let define a transition matrix $V_{\ell}$ of size $(N_{\ell-1}+1)\times(N_{\ell}+1)$ such that the $(i+1,j+1)$-component is defined to be

[TABLE]

Then given

[TABLE]

we have

[TABLE]

Suppose $N_{\ell}=N$ for all $\ell\geq 1$. Note that the first row of $V_{\ell}$ is $[1,0,\cdots,0]$ for all $\ell\geq 2$. Thus we have the following strictly increasing sequence $\{a_{\ell}\}_{\ell=1}^{\infty}$:

[TABLE]

Since $a_{\ell}\leq 1$, it converges, say, $\lim_{\ell\to\infty}a_{\ell}=a\leq 1$. Suppose $a\neq 1$, i.e., $a-1<0$ and let $p_{*}=\min_{1\leq k\leq N}p_{k}$. Then since $a_{k+1}=(\pi_{k}V_{k+1})_{1}$, we have

[TABLE]

Thus

[TABLE]

By taking limit on the both sides, we have

[TABLE]

which leads a contradiction. Therefore, $a=\lim_{\ell\to\infty}P(C_{\ell,0})=1$. It then follows from Equation 11 that

[TABLE]

which completes the proof.

Appendix B Proof of Theorem 3.3

Proof B.1**.**

Based on Lemma A.1, let consider

[TABLE]

Then if $\textbf{x}\in\mathcal{D}$, $\tilde{A}_{\ell,\textbf{x}}^{c}\subset A_{\ell}^{c}$. Thus it suffices to compute $P(\tilde{A}_{\ell,\textbf{x}}^{c})$ as

[TABLE]

For $\textbf{x}\neq\bm{0}$, let consider

[TABLE]

Note that $\bigcup_{1\leq j<\ell}U_{j,\textbf{x}}=\tilde{A}_{\ell,\textbf{x}}$ and $\tilde{A}_{\ell,\textbf{x}}^{c}=\bigcap_{1\leq j<\ell}U_{j,\textbf{x}}^{c}$. Since $P(\tilde{A}_{j,\textbf{x}}^{c}|\tilde{A}_{j-1,\textbf{x}})=0$ for all $j$, we have

[TABLE]

Note that $P(\tilde{A}_{1,\textbf{x}}^{c})=1$. Also, note that since the rows of $\bm{W}^{\ell}$ and $\bm{b}_{j}^{\ell}$ are independent,

[TABLE]

*Since the weight and biases are randomly drawn from symmetry distribution around 0 and

$P\left(\bm{W}^{j-1}_{s}\phi(\mathcal{N}^{j-2}(\textbf{x}))+\bm{b}^{j-1}_{s}=0|\tilde{A}_{j-1,\textbf{x}}^{c}\right)=0$, we obtain*

[TABLE]

Therefore, $P(\tilde{A}_{j,\textbf{x}}|\tilde{A}_{j-1,\textbf{x}}^{c})=2^{-N_{j-1}}$ and thus, $P(\tilde{A}_{j,\textbf{x}}^{c}|\tilde{A}_{j-1,\textbf{x}}^{c})=1-2^{-N_{j-1}}$. It then follows from Equation 12 that

[TABLE]

which completes the proof.

Appendix C Proof of Theorem 3.5

Proof C.1**.**

We now assume $d_{\text{in}}=1$, $N_{\ell}=N$ and without loss of generality let $\mathcal{D}\subset[-r,r]$ be a training domain for any $r>0$. Also we assume that all weights are initialized from continuous symmetric probability distributions around 0. And the biases are set to zeros.

Since $d_{\text{in}}=1$, for each $\ell$, there exist non-negative vectors $\bm{v}_{\pm}\in\mathbb{R}^{N_{\ell}}_{+}$ such that

[TABLE]

Let $B_{\ell,0}$ be the event where $\phi(\mathcal{N}^{\ell}(x))=0$, and let $B_{\ell,1}$ be the event where

[TABLE]

for some $\bm{v}_{\pm}$, $\bm{v}\in\mathbb{R}_{+}^{N_{\ell}}$ and $b>0$. Let $B_{\ell,2}$ be the event where

[TABLE]

for some linearly independent vectors $\bm{v}_{\pm}$. Then it can be checked that $P(B_{\ell+1,k}|B_{\ell,s})=0$ for all $2\geq k>s\geq 0$. Thus, it suffices to consider $P(B_{\ell+1,k}|B_{\ell,s})$ where $0\leq k\leq s\leq 2$. At $\ell=1$, since $\phi(\mathcal{N}^{1}(\textbf{x}))=\phi(\bm{W}^{1}\textbf{x})$ and $\mathcal{D}\subset[-r,r]$, we have

[TABLE]

For $\ell>1$, it can be checked that $P(B_{\ell,0}|B_{\ell-1,0})=1$, $P(B_{\ell,0}|B_{\ell-1,1})=2^{-N}$, and thus $P(B_{\ell,1}|B_{\ell-1,1})=1-2^{-N}$. For $P(B_{\ell,0}|B_{\ell-1,2})$ and $P(B_{\ell,1}|B_{\ell-1,2})$, we observe the followings. In $B_{\ell,2}$, we have

[TABLE]

Since $\bm{v}_{\pm}$ are nonzero vectors in $\mathbb{R}^{N_{\ell}}_{+}$ and thus satisfy $\langle\bm{v}_{+},\bm{v}_{-}\rangle\geq 0$, by the assumption, we obtain

[TABLE]

*Thus we have *

[TABLE]

where the expectation is taken under $p_{\bm{v}_{\pm}}$. For $P(B_{\ell,1}|B_{\ell-1,2})$, there are only three ways to move from $B_{\ell-1,2}$ to $B_{\ell,1}$. That is

[TABLE]

Thus $P(B_{\ell,1}|B_{\ell-1,2})=P(B_{\ell,1}^{(A)}|B_{\ell-1,2})+P(B_{\ell,1}^{(B)}|B_{\ell-1,2})+P(B_{\ell,1}^{(C)}|B_{\ell-1,2}).$ Note that due to the symmetry, $P(B_{\ell,1}^{(A)}|B_{\ell-1,2})=P(B_{\ell,1}^{(B)}|B_{\ell-1,2})$. Thus we have

[TABLE]

Note that

[TABLE]

Since $P(A_{\ell+1})=P(B_{\ell,0})$ where $A_{\ell}$ is defined in Equation 10, we aim to estimate $P(B_{\ell,0})$. Let $V_{\ell}$ be the transition matrix whose $(i+1,j+1)$-component is $P(B_{\ell,j}|B_{\ell-1,i})$. Then $P(B_{\ell,0})=\pi_{1}V_{2}\cdots V_{\ell}$ where $(\pi_{1})_{j}=P(B_{1,j-1})$. By letting

[TABLE]

we obtain

[TABLE]

where $D_{\ell}=\text{diag}(V_{\ell})$.

To find a lower bound of $P(B_{\ell,0})$, we consider the following transition matrix $\mathcal{P}$ of size $3\times 3$ which is defined to be

[TABLE]

It can be checked that

[TABLE]

where $\gamma=\frac{\mathcal{P}_{32}}{\mathcal{P}_{22}-\mathcal{P}_{33}}=\frac{2^{-N+1}+(N-2)4^{-N}}{2^{-N}+(N-1)4^{-N}}=1+\frac{2^{-N}-4^{-N}}{2^{-N}+(N-1)4^{-N}}$ and $D=\text{diag}(\mathcal{P})$. Thus we have

[TABLE]

Similarly, we obtain

[TABLE]

where $g(x)=1-2^{-N+1}-(N-1)x$, $p_{\ell}=(V_{\ell})_{31}=P(B_{\ell,0}|B_{\ell-1,2})$, $\gamma_{\ell,i}=\gamma_{i}$ for $1\leq i\leq\ell$, $\gamma_{\ell,\ell+1}=1$,

[TABLE]

We want to show that

[TABLE]

where

[TABLE]

Let denote $\bar{\gamma_{\ell,i}}:=\gamma_{\ell,i}-1$ and $g_{i}:=g(p_{i})$. It then suffices to show that

[TABLE]

Note that $4^{-N}=p_{1}\leq p_{j}<2^{-N}$, $\mathcal{P}_{22}>g(p_{j})$, and thus $\mathcal{P}_{22}^{\ell}>(g(p_{1}))^{\ell}\geq\prod_{j=1}^{\ell}g(p_{j})$. Also,

[TABLE]

Thus we have

[TABLE]

Therefore,

[TABLE]

which implies that

[TABLE]

Furthermore, it can be checked that

[TABLE]

Appendix D Proof of Theorem 3.7

Proof D.1**.**

Since $\mathcal{N}^{L}(\mathbf{x})$ is a constant function, it follows from Lemma A.1 that with probability 1, there exists $\ell$ such that $\phi(\mathcal{N}^{\ell}(\textbf{x}))\equiv\mathbf{0}$. Then the gradients of the loss function with respect to the weights and biases in the $1,\cdots,\ell$-th layers vanish. Hence, the weights and biases in layers $1,\dots,\ell$ will not change when a gradient based optimizer is employed. This implies that $\mathcal{N}^{L}(\textbf{x})$ always remains to be a constant function as $\phi(\mathcal{N}^{\ell}(\textbf{x}))\equiv\mathbf{0}$. Furthermore, the gradient method changes only the weights and biases in layers $\ell+1,\dots,L$. Therefore, the ReLU NN can only be optimized to a constant function, which minimizes the loss function $\mathcal{L}$.

Appendix E Proof of Theorem 4.1

Lemma E.1**.**

Let $\bm{v}\in\mathbb{R}^{n+1}$ be a vector such that $\bm{v}_{k}\sim\text{P}$ and $\bm{v}_{-k}\sim N(0,\sigma^{2}\bm{I}_{n})$ where $\bm{v}_{-k}$ is defined in Equation 6. For any nonzero vector $\textbf{x}\in\mathbb{R}^{n+1}$ whose $k$-th element is positive (i.e., $\textbf{x}_{k}>0$), let

[TABLE]

Then

[TABLE]

where $F_{\text{P}}(t)$ is the cdf of P.

Proof E.2**.**

We first recall some properties of the normal distribution. Let $Y_{1},\cdots,Y_{n}$ be i.i.d. random variables from $N(0,\sigma^{2})$ and let $\bm{Y}_{n}=(Y_{1},\cdots,Y_{n})$. Then for any vector $\bm{a}\in\mathbb{R}^{n}$,

[TABLE]

Suppose $\bm{v}$ is a random vector generated in the way described in Subsection 4.1. Then for any $\bm{x}\in\mathbb{R}^{n+1}$,

[TABLE]

where $\tilde{\sigma}^{2}=\|\tilde{\bm{x}}_{-k}\|^{2}\sigma^{2}$ and $Z\sim N(0,\tilde{\sigma}^{2})$. If $\bm{x}_{k}>0$ and $\|{\textbf{x}}_{-k}\|^{2}>0$, we have

[TABLE]

Therefore, it suffices to compute $P(\bm{v}_{k}>Z)$. Let $f_{Z}(z)$ be the pdf of $Z$ and $f_{\text{P}}(x)$ be the pdf of $\bm{v}_{k}\sim\text{P}$. Then

[TABLE]

*which completes the proof. *

Proof E.3**.**

For each $\ell$ and $s$, let $k_{s}^{\ell}$ be randomly uniformly chosen in $\{1,\cdots,N_{\ell-1}+1\}$. Let $\bm{V}^{\ell}_{s}=[\bm{W}^{\ell}_{s},\bm{b}^{\ell}_{s}]$ and $\textbf{n}^{\ell}(\textbf{x})=\left[\phi(\mathcal{N}^{\ell}(\textbf{x})),1\right]$. Recall that for a vector $\textbf{v}\in\mathbb{R}^{N+1}$,

[TABLE]

To emphasize the dependency of $k_{s}^{\ell}$, we denote $\bm{V}^{\ell}_{s}$ whose $k$-th component is generated from P by $\bm{V}^{\ell}_{s}(k)$.

The proof can be started from Equation 13. It suffices to compute

[TABLE]

Note that for each $s$,

[TABLE]

Also, from Lemma E.1, we have

[TABLE]

where

[TABLE]

Note that if $\textbf{n}^{\ell-1}_{k}(\textbf{x})=0$, the above probability is simply $1/2$. Also, since the event $\tilde{A}_{j-1,\textbf{x}}^{c}$ is given, we know that $\phi(\mathcal{N}^{\ell}(\textbf{x}))\neq 0$. Thus,

[TABLE]

We thus denote

[TABLE]

where $\delta_{\ell-1,\textbf{x}}$ is $\mathcal{F}_{\ell-2}$ measurable whose value lies in $(0,0.5]$ and it depends on x. It then follows from Equation 12 that

[TABLE]

where ${P}$ is the probability distribution with respect to $\{\bm{W}^{j},\bm{b}^{j}\}_{1\leq j\leq\ell}$. Note that since the weight matrix and the bias vector of the first layer is initialized from a symmetric distribution, $\delta_{1,\textbf{x}}=0$. Let $\delta_{1,\textbf{x}}^{*}=0$. By the mean value theorem, there exists some numbers $\delta_{2,\textbf{x}}^{*},\cdots,\delta_{\ell-1,\textbf{x}}^{*}\in(0,0.5]$ such that

[TABLE]

Let $\bm{\delta}^{*}_{\textbf{x}}=(\delta_{1,\textbf{x}}^{*},\cdots,\delta_{\ell-1,\textbf{x}}^{*})$. By setting

[TABLE]

the proof is completed.

Appendix F Proof of Theorem 4.5

Let $X_{\ell}\sim\text{P}_{\ell}$ whose pdf is denoted by $f_{P_{\ell}}(x)$. Suppose $0<\mu_{\ell,i}^{\prime}=E[X_{\ell}^{i}]<\infty$ for $i=1,2$. We then define three probability distribution functions:

[TABLE]

Also we denote its corresponding cdfs by $F^{(i)}_{P_{\ell}}(t)=\int_{0}^{t}f^{(i)}_{P_{\ell}}(x)dx$. For notational convenience, let us assume that $\text{P}_{\ell}=\text{P}$ and $M_{\ell}=M$ for all $\ell$ and define

[TABLE]

We denote the pdf of normal distributions by $f_{Y}^{\ell,k}(y)$ where

[TABLE]

where $1\leq k\leq N_{\ell-1}+1$. Recall that for a vector $\textbf{v}\in\mathbb{R}^{N+1}$,

[TABLE]

For any probability distribution P defined on $[0,M]$ whose pdf is denoted by $f_{P}(x)$, let

[TABLE]

Let $Y_{\ell,k}\sim N(0,\tilde{\sigma}_{\ell,k}^{2})$ and $X\sim\text{P}$. Then

[TABLE]

Therefore, $0\leq\mathcal{J}^{\ell,k}_{\text{P}}\leq 0.5$.

The proof of Theorem 4.5 starts with the following lemma.

Lemma F.1**.**

Suppose $Y\sim N(0,\tilde{\sigma}_{\ell,k}^{2})$, $X\sim\text{P}_{\ell}$ and $\mu_{\ell,1}^{\prime}\geq M/2$. Then

[TABLE]

Proof F.2**.**

Let $I_{\ell,k}$ be the quantity of our interest:

[TABLE]

Without loss of generality, we can assume $M=1$. This is because

[TABLE]

where $u=x/M$ and $s=y/M$. Here $MU\sim\text{P}_{\ell}$ and $S\sim N(0,\tilde{\sigma}^{2}_{\ell,k}/M^{2})$.

Let us first consider the inner integral. Let $G(y)=\int_{y}^{1}(x^{2}-2xy+y^{2})f_{P_{\ell}}(x)dx$. It then can be written as follow:

[TABLE]

Note that since $0\leq y\leq 1$, and thus $y^{2}\leq y$, we have

[TABLE]

Since $2\mu_{\ell,1}^{\prime}\geq 1$ and

[TABLE]

we have $G(y)\leq\mu_{\ell,2}^{\prime}\int_{y}^{1}f^{(2)}_{P_{\ell}}(x)dx$. Therefore,

[TABLE]

*where $\mathcal{J}^{\ell,k}_{\text{P}}$ is defined in Equation 16. Since $0\leq\mathcal{J}^{\ell,k}_{\text{P}}\leq 0.5$, we conlcude that $I_{\ell,k}\leq\mu_{\ell,2}^{\prime}/2$. *

Proof F.3**.**

It follows from

[TABLE]

that we have

[TABLE]

Since $k_{j}^{\ell}$ is randomly uniformly selected, by taking the expectation with respect to $k_{j}^{\ell}$, we have

[TABLE]

Since the above is independent of $j$, and $q^{\ell}(\textbf{x})=\|\mathcal{N}^{\ell}(\textbf{x})\|^{2}/N_{\ell}$,

[TABLE]

At a fixed $\ell$ and for each $t=1,\cdots,N_{\ell}$, let $\textbf{v}^{\ell,t}=[\textbf{W}^{\ell}_{t},\textbf{b}^{\ell}_{t}]\in\mathbb{R}^{N_{\ell-1}+1}$ be a random vector such that $\textbf{v}_{-k_{t}}^{\ell,t}\sim N(0,\sigma_{\ell}^{2}\bm{I}_{N_{\ell-1}})$ and $\textbf{v}_{k_{t}}^{\ell,t}\sim\text{P}_{\ell}$. Then $\mathcal{N}^{\ell}_{t}(\textbf{x})=\langle\textbf{v}^{\ell,t},\textbf{n}^{\ell}(\textbf{x})\rangle$. Given $\mathcal{N}^{\ell-1}(\textbf{x})$, assuming $(\textbf{n}^{\ell-1}_{k_{t}}(\textbf{x}))>0$, one can view $\mathcal{N}^{\ell}_{t}(\textbf{x})$ as

[TABLE]

Thus $\phi(\mathcal{N}^{\ell}_{t}(\textbf{x}))^{2}=\left(\textbf{n}^{\ell-1}_{k_{t}}(\textbf{x})\right)^{2}\phi\left(\sigma_{\ell,k_{t}}Z+X_{\ell}\right)^{2}$ and we obtain

[TABLE]

By multiplying $(\textbf{n}^{\ell-1}_{k_{t}}(\textbf{x}))^{2}$ in the both sides, we have

[TABLE]

Since $k_{t}$ is randomly selected, by taking expectation w.r.t. $k_{j}$, we have

[TABLE]

By the Cauchy-Schwarz inequality, the third term in the right hand side of Equation 18 can be bounded by

[TABLE]

where $\sigma_{w}=\sigma_{\ell}\sqrt{N_{\ell-1}}$. Let $I_{\ell,k_{t}}=\int_{0}^{M}\int_{y}^{M}(x-y)^{2}dF_{P}(x)f_{Y}^{\ell,k_{t}}(y)dy$. It then follows from Lemma F.1 that $I_{\ell,k}\leq\mu_{\ell,2}^{\prime}/2$ for any $k$. Thus the second term in the right hand side of Equation 18 can be bounded by

[TABLE]

Therefore, we obtain

[TABLE]

where

[TABLE]

Since $\mathcal{C}_{\ell}$ is independent of $t$,

[TABLE]

It then follows from Equation 17,

[TABLE]

that

[TABLE]

Thus we have

[TABLE]

where

[TABLE]

By taking expectation with respect to $\mathcal{N}^{\ell-1}(\textbf{x})$, we obtain

[TABLE]

The lower bound can be obtained by dropping the second and the third terms in Equation 18. Thus

[TABLE]

where

[TABLE]

Appendix G Randomized asymmetric initialization

⬇

1import numpy as np

2

3

4def RAI(fan_in, fan_out):

5 """Randomized asymmetric initializer.

6 It draws samples using RAI where fan_in is the number of input units in the weight tensor and fan_out is the number of output units in the weight tensor.

7 """

8 V = np.random.randn(fan_out, fan_in + 1) * 0.6007 / fan_in ** 0.5

9 for j in range(fan_out):

10 k = np.random.randint(0, high=fan_in + 1)

11 V[j, k] = np.random.beta(2, 1)

12 W = V[:, :-1].T

13 b = V[:, -1]

14 return W.astype(np.float32), b.astype(np.float32)