Compelling ReLU Networks to Exhibit Exponentially Many Linear Regions at Initialization and During Training

## Abstract

Abstract A neural network with ReLU activations may be viewed as a composition of
piecewise linear functions. For such networks, the number of distinct linear regions
expressed over the input domain has the potential to scale exponentially with depth,
but it is not expected to do so when the initial parameters are chosen randomly.
Therefore, randomly initialized models are often unnecessarily large, even when approximating simple functions. To address this issue, we introduce a novel training strategy:
we first reparameterize the network weights in a manner that forces the network
to exhibit a number of linear regions exponential in depth. Training first on
our derived parameters provides an initial solution that can later be refined by
directly updating the underlying model weights. This approach allows us to learn
approximations of convex, one-dimensional functions that are several orders of
magnitude more accurate than their randomly initialized counterparts. We further demonstrate how to extend our approach to multidimensional and non-convex functions, with similar benefits observed.

## 1 Introduction

Beyond complementary advances in areas like hardware, storage, and networking, the success of neural networks is primarily due to their ability to efficiently capture and represent nonlinear functions (Gibou et al., ).
In a neural network, the goal of an activation function is to introduce nonlinearity between the network’s layers so that the network does not simplify to a single linear function.
The rectified linear unit (ReLU) has a unique interpretation in this regard. Since it either deactivates a neuron or acts as an identity, the resulting transformation on each individual input remains linear. However, each possible configuration of active and inactive neurons can produce a unique linear transformation over a particular region of input space. The number of these activation patterns and their corresponding linear regions provides a way to measure the expressivity of a ReLU network(^1^11The reader is referred to Chmielewski-Anders () for definitions of linear regions, activation patterns, and activation regions.) and can theoretically scale exponentially with the depth of the network
(Montufar et al., ; Serra et al., ).
Hence deep architectures may outperform shallow ones.

Surprisingly, though, a sophisticated theory of how to best encode functions into ReLU networks is lacking, and in practice, adding depth is often observed to help less than one might expect from this exponential intuition.
Lacking more advanced theory, practitioners typically use random parameter initialization and gradient descent, the drawbacks of which often lead to extremely inefficient solutions.
Hanin & Rolnick () show a rather disappointing bound pertaining to randomly initialized networks: they prove that the average number of linear regions formed upon initialization is entirely independent of the configuration of the neurons, so depth is not properly utilized.
They observed that gradient descent has a difficult time creating new activation regions and that their bounds approximately held after training.
As we will discuss later, the number of linear regions is not actually a model property that gradient descent can directly optimize.
Gradient descent is also prone to redundancy; Frankle & Carbin () show how around 95% of weights may ultimately be eliminated from a network without significantly degrading accuracy.

The present work aims to begin eliminating these inefficiencies, starting in a simple one-dimensional setting.
Drawing inspiration from existing theoretical ReLU constructions, our novel contributions include a special reparameterization of a ReLU network that forces it to maintain an exponential number of activation patterns over the input domain.
We then demonstrate a novel pretraining strategy, which trains these derived parameters before manipulating the underlying matrix weights.
This allows the network to discover solutions that are more accurate and unlikely to be found otherwise.
Our technique directly rectifies the issues raised by Hanin & Rolnick (). Firstly, we initialize immediately with an exponential number of linear regions, which would not be expected to happen otherwise. Secondly, we are not reliant on gradient descent to “discover” new regions; we only need to maintain the existing ones, which is already guaranteed by the reparameterization during the pretraining stage. Additionally, this pretraining step can avoid areas of the loss landscape where unassisted gradient descent might make short-sighted optimizations that eliminate activation regions.
Our results demonstrate that minimizing the reliance of network training on unassisted gradient descent can reliably produce error values orders of magnitude lower than a traditionally-trained network of equal size.
Although the preliminary results in this theoretical study pertain to relatively low-dimensional functions, the paper concludes with our views on extending these theoretical exponential benefits for ReLU networks to arbitrary smooth functions with arbitrary dimensionality, which would have significant practical utility.

## 2 Related Work

This work is primarily concerned with a novel training methodology, but it also possesses a significant approximation theoretic component.
In one dimension, where we initially focus, our work can be viewed as a first step towards generalizing the approximation to $x^{2}$ we review in this section. The reparameterization we employ modifies this method to become trainable in order to represent other convex one-dimensional functions, and then converts that result back into a matrix representation.

Figure: Figure 1: (Top to bottom) Composed triangle waves; using collections of the above function to approximate $x^{2}$; derivatives of the above approximations.
Refer to caption: x1.png

### 2.1 Related Work in Approximation

Infinitely wide neural networks are known to be universal function approximators, even with only one hidden layer (Hornik et al., ; Cybenko, ).
Infinitely deep networks of fixed width are universal approximators as well (Lu et al., ; Hanin, ).
In finite cases, the trade-off between width and depth is often studied in assessing a network’s ability to approximate (learn) a function.

Notably, there exist functions that can be represented with a sub-exponential number of neurons in a deep architecture, yet which require an exponential number of neurons in a wide and shallow architecture.
For example, Telgarsky () shows that deep neural networks with ReLU activations on a one-dimensional input are able to generate symmetric triangle waves with an exponential number of linear segments (shown in Figure 1 as the ReLU network $T(x)$).
This network functions as follows: each layer takes a one-dimensional input on $[0,1]$, and outputs a one-dimensional signal also on $[0,1]$.
The function they produce in isolation is a single symmetric triangle.
Together in a network, each layer inputs its output to the next, performing function composition.
Since each layer converts lines from 0 to 1 into triangles, it doubles the number of linear segments in its input signal, exponentially scaling with depth.

The same effect can be achieved with non-symmetric triangle waves (Huchette et al., ) (or any shape that sufficiently “folds” the input space (Montufar et al., )). Our reparameterization strategy focuses on non-symmetric triangle waves. The location on $(0,1)$ of a triangular function’s peak acts as an alternative parameterization of the function, instead of its weight matrix representation. Unlike the raw weights of a depth $d$ ReLU network, any values chosen for the triangle peak parameters will result in the creation of $2^{d}$ linear regions. Furthermore, this parameterization is trainable. By using the peak location to set the raw weights of a layer, the gradients can backpropagate through the raw weights to update the triangle peaks, effectively confining the network to a subspace of weights that generates many linear regions.

The dilated triangular waveforms produced in this manner are not particularly useful on their own. Their oscillations quickly become excessively rapid, and their derviatives do not exist everywhere (especially in the infinite depth limit). But these problems can be rectified by taking a sum over the layers of a network.
Yarotsky () and Liang & Srikant () demonstrate the ability to construct $y=x^{2}$ on $[0,1]$ with exponential accuracy using symmetric triangle waves.
To produce their approximation, one begins with $f_{0}(x)=x$, then computes $f_{1}(x)=f_{0}(x)-T(x)/4$, $f_{2}(x)=f_{1}(x)-T(T(x))/16$, $f_{3}(x)=f_{2}(x)-T(T(T(x)))/64$, and so forth, as pictured in Figure 1.
As these successive approximations are computed, Figure 1 plots their convergence to $x^{2}$, as well as the convergence of the derivative to $2x$. Our reparameterization generalizes this approximation to use non-symmetric triangle waves to approximate a wider class of convex, differentiable, one-dimensional functions.

The $x^{2}$ approximation is used by other theoretical works as a building block to guarantee exponential convergence rates in more complex systems.
One possible use case is to construct a multiplication gate.
Perekrestenko et al. () does so via the identity $(x+y)^{2}=x^{2}+y^{2}+2xy$.
The squared terms can all be moved to one side, expressing the product $xy$ as a linear combination of squared terms.
They then further assemble these multiplication gates into an interpolating polynomial, which can have an exponentially decreasing error when the interpolation points are chosen to be the Chebyshev nodes.
Polynomial interpolation does not scale well into high dimensions, so this and papers with similar approaches will usually come with restrictions that limit function complexity: Wang et al. () requires low input dimension, Montanelli et al. () uses band limiting, and Chen et al. () approximates low-dimensional manifolds. These works all make use of a fixed representation of $x^{2}$. If our networks were substituted in for the $x^{2}$ approximation, these works would provide theoretical guarantees about the capabilities of the resulting model. Even though the approximation rates will not scale well with input dimension, they serve as a bound that can be improved upon. In Section 5 we further elaborate on how to use our networks to represent higher-dimensional or non-convex functions.

Lastly, other works focus on showing how ReLU networks can encode and subsequently surpass traditional approximation methods (Lu et al., ; Daubechies et al., ).
Interestingly, certain fundamental themes from above like composition, triangles, or squaring are still present.
Another interesting comparison of the present work is to Ivanova & Kubat (), which uses decision trees as a means to initialize neural networks.
It is a sigmoid/classification analogy to this work, but rather than an attempting to improve neural networks with decision trees, it is an attempt to improve decision trees with neural networks. Finally, Park et al. () approaches the problem of linear region maximization from an information theory perspective and uses a loss penalty rather than a reparameterization to increase the number of linear regions.

### 2.2 Related Work in Initialization

Our work seeks to improve network initialization by making use of explicit theoretical constructs.
This stands in sharp contrast the current standard approach, which treats neurons homogeneously.
Two popular initialization methods implemented in PyTorch are the Kaiming (He et al., ) and Xavier initialization (Glorot & Bengio, ).
They use weight values that are sampled from distributions defined by the input and output dimension of each layer.
Aside from sub-optimal approximation power associated with random weights, a common issue is that an entire ReLU network can collapse into outputting a constant value at initialization. This is referred to as the dying ReLU phenomenon (Qi et al., ; Nag et al., ).
It occurs when the initial weights and biases cause every neuron in a particular layer to output a negative value.
The ReLU activation then sets the output of that layer to 0, blocking any gradient updates.
Worryingly, as depth goes to infinity, the dying ReLU phenomenon becomes increasingly likely (Lu et al., ).
Several papers propose solutions: Shin & Karniadakis () use a data-dependent initialization, while Singh & Sreejith () introduce an alternate weight distribution called RAAI that can reduce the likelihood of the issue and increase training speed.
We observed during our experiments that RAAI greatly reduces, but does not eliminate the likelihood of dying ReLU.
Our approach enforces a specific network structure that does not collapse in this manner.

## 3 Initialization and Pretraining Construction

We begin by discussing how to deliberately architect the weights of a 4-neuron-wide, depth $d$ ReLU network to induce a number of linear segments exponential in $d$.
As alluded to in Section 2.1, our initialization algorithm is simply to associate a nonsymmetric triangle function with each layer, and then to randomize the location of each triangle’s peak in $(0,1)$; any choice produces $2^{d}$ regions.
This yields a reparameterization of the ReLU network in terms of the triangle peak locations, $a_{i}$, rather than in terms of the network’s weights. The pretraining step of our algorithm then involves learning and updating these parameters.
This is done by setting individual neuron weights in accordance with the procedures described in this section, and then backpropagating the gradients to update the triangle peak parameters.
Pseudocode for our entire algorithm is provided in Appendix A.1.

To illustrate how we set weights, we first describe the mathematical functions that arise in our analysis.
Triangle functions are defined as

$$ $\displaystyle T_{i}(x)=\begin{cases}\frac{x}{a_{i}}&0\leq x\leq a_{i}\\ 1-\frac{x-a_{i}}{1-a_{i}}&a_{i}\leq x\leq 1\end{cases}$ $$

where $0<a_{i}<1$. This produces a triangular shape with a peak at $x=a_{i}$ and with both endpoints satisfying $y=0$. Each layer considered in isolation would compute these if directly fed the input signal. $T_{i}(x)$’s derivatives are the piecewise linear functions:

$$ $T^{\prime}_{i}(x)=\begin{cases}\frac{1}{a_{i}}&0<x<a_{i}\\ \frac{1}{1-a_{i}}&a_{i}<x<1\end{cases}$ (1) $$

In a deep network, the layers feed into each other, composing their respective triangle functions:

$$ $W_{i}(x)=\bigcirc_{j=0}^{i}T_{j}(x)=T_{i}(T_{i-1}(...T_{0}(x)))$ (2) $$

These triangle waves will have $2^{i}$ linear regions, doubling with each layer. The output of the network in pretraining will be a weighted sum over the triangle waves formed at each layer. Assuming the network to be infinitely deep, we have

$$ $F(x)=\sum_{i=0}^{\infty}s_{i}W_{i}(x)$ (3) $$

where $s_{i}$ are scaling coefficients on each of the composed triangular waveforms $W_{i}$.

Figure: Figure 2: On the right is a network representation of a triangle function. The middle shows that triangle function as a hidden layer of a network. The one-dimensional input and output of a triangle function is converted into shared weights. A full pretraining network is assembled on the left.
Refer to caption: x2.png

To encode these functions into the weights of a ReLU network, we begin with triangle functions (see the right subnetwork in Figure 2).
Its maximum output is $1$ at the peak location $a\in(0,1)$.
Neuron $t_{1}$ simply preserves the input signal.
Meanwhile, $t_{2}$ is negatively biased, deactivating it for inputs less than $a$.
Subtracting $t_{2}$ from $t_{1}$ changes the slope at the point where $t_{2}$ begins outputting a nonzero signal.
The weight $-1/(a-a^{2})=-(1/a+1/(1-a))$ is picked to completely negate $t_{1}$’s positive influence, and then produce a negative slope.
When these components are combined into a deep network, the individual triangles they form will be composed with those earlier in the network, but the two neurons will still behave analogously. Considering the output waveform of $t_{1}$ and $t_{2}$ neurons over the entire input domain $[0,1]$ and at arbitrary depth, we have that $t_{1}$ neurons are always active, outputting complete triangle waves. In contrast, $t_{2}$ neurons are deactivated for small inputs, so they output an alternating sequence of triangles and inactive regions (see Figure 3).

Figure: Figure 3: Each colored line shows the output signal of a neuron with respect to the input to the network. Colors match the corresponding neurons in Figure 2.
Refer to caption: x3.png

Naively using the output of one triangle generator (from Figure 2) as input to the next would form a $1\times 2\times 1\times 2\times 1...$ shape, but this is unnecessarily deep. We can replicate the one-dimensional function composition in the hidden layers on the left side of Figure 2 by using weight sharing instead. Any outgoing weight from $t_{1}$ or $t_{2}$ is shared; every neuron taking in a triangle wave as input does so by combining $t_{1}$ and $t_{2}$ in the same proportion. In this way, we can avoid having to use the extra intermediate neurons. This is evident in Figure 3 where the waves output by $t_{1}$ and $t_{2}$ neurons posses identical slopes in regions where both neurons are active.

There are two other neurons in the pretraining network’s layers, each with a specific role. The accumulation neuron (marked as “sum” in Figure 2) maintains a weighted sum of all previous triangle waves through each layer. This is similar to a residual connection (He et al., ), except that no other neurons take the sum as input (its output is convex, and therefore not rapidly oscillating, and detrimental to making more linear regions).
If the sum neuron were naively implemented, it would multiply the $t_{1}$ and $t_{2}$ weights by the sum coefficients.
Based on Theorem 4, these coefficients will be exponentially decaying, so learning these weights directly may cause conditioning issues.
Instead, all weights in each layer are multiplied by the ratios between successive scaling coefficients, so that the outputs of $t_{1}$ and $t_{2}$ neurons decay in amplitude in each layer.
A conventional bias will have no connections to prior layers, so it will be unable to adjust to the rescaling neuron outputs without having to learn exponentially small values.
Therefore, a fourth neuron is configured to output a constant signal, so the weight other neurons place on this output can replace their bias.
The bias neuron will connect to the previous layer’s bias neuron so that the constant signal can scale down gradually with each layer. We handle this scaling explicitly, rather than relying on regularizers like batch norm (Ioffe & Szegedy, ) or layer norm (Lei Ba et al., ), partially in order to preserve the mathematical transparency of our model, and partially because this solution arises naturally from our derivations.

### 3.1 Differentiable Model Output

Given the rapid oscillations of the triangle waves formed at each layer, the resulting output of the network will appear jagged or fractal with many choices of scaling parameters. This can lead to the network settling into solutions that are not predictive of the areas in between training points (or possibly not even interpolating the training data) due to rapid oscillation. Our main mathematical result addresses this issue by forcing differentiability of the network output. To achieve differentiability, it turns out that for this class of functions, the peak locations of the triangle waves uniquely determine the scales with which to sum them:

###### Theorem 3.1 .

$F^{\prime}(x)$ is defined on $[0,1]$ only if the scaling coefficients $s_{i}$ are selected based on the triangle peaks $a_{i}$ according to:

$$ $s_{i+1}=s_{i}(1-a_{i+1})a_{i+2}$ (4) $$

###### Theorem 3.1 .

## 4 One-Dimensional Experiments

The goal of deep learning is to train a network to approximate a (generally unknown) nonlinear function.
Accordingly, in this section, we implement our method to learn several convex one-dimensional curves.
These nonlinear functions are known—so that we can measure error—and we remark on more difficult functions in Sections 5 and 6.
The aim of these experiments is twofold: (1) we would like to determine how to learn the most effective function representations possible, and (2) to explore how the utilization of an increased number of linear regions can affect a network’s ability to capture underlying nonlinearity in its training data.
To demonstrate that our networks can learn function representations that better utilize depth, we benchmark against PyTorch’s (Paszke et al., ) default settings (nn.linear() uses Kaiming initialization), as well as the RAAI distribution from Singh & Sreejith (), and produce errors that are orders of magnitude lower than both.
For our experimental networks, we use a common set of initialization points where the triangle peaks and scaling parameters are chosen according to our main theorem (Equation 4).
We compare the effect of reparameterized pretraining against skipping immediately to standard gradient descent (from a reparameterized initialization as in Section 3).
We also compare training the scaling parameters freely during pretraining, instead of choosing them to achieve differentiability.
The intent of both comparisons is see if pretraining and further differentiability constraints facilitate smoother navigation of the loss landscape.
Lastly, we conduct a second round of tests to determine if pretrained networks display an enhanced predictive capacity on unseen data points, as might be expected if they can leverage greater nonlinearity in their outputs.

### 4.1 Experimental Setup

All models are trained using Adam (Kingma & Ba, ) as the optimizer with a learning rate of 0.001 for 1,000 epochs to ensure convergence (see Appendix for a brief analysis of different learning rates).
Each network is four neurons wide with five hidden layers, along with a one-dimensional input and output.
The loss function used is the mean squared error, and the average and minimum loss are recorded for 30 models of each type.
The networks unique to this paper share a common set of starting locations using the construction in Section 3, so that the effects of each training regimen are directly comparable.
As in related papers (Perekrestenko et al., ; Daubechies et al., ), we focus for the moment on one-dimensional examples, which are sufficient to demonstrate our proposed theory and methodology.
The four curves we approximate are $x^{3}$, $x^{11}$, $\tanh(3x)$, and a quarter period of a sine wave. The curves are chosen to capture a variety of convex one-dimensional functions.
To approximate the sine and the hyperbolic tangent, the triangle waves are added to the line $y=x$.
For the other approximations, the waves are subtracted.
This requires the first scaling factor to be $a_{0}*a_{1}$ instead of $(1-a_{0})*a_{1}$.
The first set of data is 500 evenly spaced points on the interval $[0,1]$ for each of the curves.
This is chosen to be very dense deliberately, to try to evoke the most accurate representations from each network.
We determine from these tests that pretraining with differentiability enforced produces the best results, so we compare it to standard networks in our second set of experiments.
We use a second set of data consisting of only 10 points, with a test set of 10 points spaced in between so as to be as far away from learned data as possible.
The goal of this set of experiments is to compare the predictive capacities of the networks on unseen data.

### 4.2 Numerical Results

**Table 1: Minimum and mean (30 samples) MSE error approximating $y=x^{3}$ and $x^{11}$.**
| Training Type | Min $x^{3}$ | Min $x^{11}$ | Mean $x^{3}$ | Mean $x^{11}$ |
| --- | --- | --- | --- | --- |
| Default Network | $2.11\times 10^{-5}$ | $2.19\times 10^{-5}$ | $7.20\times 10^{-2}$ | $2.82\times 10^{-2}$ |
| RAAI Distribution | $2.14\times 10^{-5}$ | $4.40\times 10^{-5}$ | $3.97\times 10^{-2}$ | $4.12\times 10^{-2}$ |
| Pretraining Skipped | $7.63\times 10^{-7}$ | $1.86\times 10^{-5}$ | $3.89\times 10^{-5}$ | $3.56\times 10^{-4}$ |
| Differentiability Not Enforced | $1.64\times 10^{-7}$ | $3.20\times 10^{-6}$ | $1.02\times 10^{-5}$ | $3.73\times 10^{-5}$ |
| Differentiability Enforced | $\mathbf{7.86\times 10^{-8}}$ | $\mathbf{8.86\times 10^{-7}}$ | $\mathbf{5.27\times 10^{-7}}$ | $\mathbf{7.87\times 10^{-6}}$ |

**Table 2: Minimum and mean (30 samples) MSE error approximating $y=\sin(x)$ and $y=\tanh(3x)$.**
| Training Type | Min $\sin(x)$ | Min $\tanh(3x)$ | Mean $\sin(x)$ | Mean $\tanh(3x)$ |
| --- | --- | --- | --- | --- |
| Default Network | $4.50\times 10^{-5}$ | $5.75\times 10^{-5}$ | $1.15\times 10^{-1}$ | $1.96\times 10^{-1}$ |
| RAAI Distribution | $3.59\times 10^{-5}$ | $1.09\times 10^{-5}$ | $3.63\times 10^{-2}$ | $2.31\times 10^{-2}$ |
| Pretraining Skipped | $1.96\times 10^{-7}$ | $1.07\times 10^{-6}$ | $1.93\times 10^{-5}$ | $8.38\times 10^{-5}$ |
| Differentiability Not Enforced | $\mathbf{4.41\times 10^{-8}}$ | $1.49\times 10^{-7}$ | $1.47\times 10^{-5}$ | $3.81\times 10^{-4}$ |
| Differentiability Enforced | $5.06\times 10^{-8}$ | $\mathbf{6.82\times 10^{-8}}$ | $\mathbf{2.21\times 10^{-7}}$ | $\mathbf{8.42\times 10^{-7}}$ |

**Table 3: Minimum errors on unseen points from training on sparse data.**
| Training Type | Min $x^{3}$ | Min $x^{11}$ | Min $\sin(x)$ | Min $\tanh(3x)$ |
| --- | --- | --- | --- | --- |
| Default Network | $2.41\times 10^{-4}$ | $2.14\times 10^{-3}$ | $2.27\times 10^{-5}$ | $1.60\times 10^{-4}$ |
| Differentiable Pretraining | $\mathbf{5.65\times 10^{-6}}$ | $\mathbf{6.53\times 10^{-4}}$ | $\mathbf{7.92\times 10^{-7}}$ | $\mathbf{5.09\times 10^{-6}}$ |

Our first set of results are shown in Tables 1 and 2, wherein we observe several important trends.
First, the worst performing networks are those that rely on purely randomized initializations.
Even the networks that forgo pretraining benefit from initializing with many activation regions.
When pretraining constraints are used, they are able to steer gradient descent to the best solutions, resulting in reductions in minimum error of three orders of magnitude over default networks.

Pretraining with differentiability enforced also closes the gap between the minimum and mean errors compared to other setups.
This indicates that these loss landscapes are indeed the most reliable to traverse. Enforcing differentiability during pretraining can impart a bias towards smoother solutions during subsequent unassisted gradient descent.

The last trend to observe is the poor average performance of default networks.
In a typical run of these experiments, around half of the default networks collapse.
RAAI is able to eliminate most, but not all of the dying ReLU instances due to its probabilistic nature, so it, too, has high mean error.

Figure: Figure 4: Comparison between standard initialization/gradient descent and pretraining with differentiablity enforced. Using more linear regions allows the curve to better predict the test points.
Refer to caption: x4.png

Our second set of results is shown in Table 3 and in Figure 4. Here the most important impact of utilizing exponentially many linear regions is demonstrated. Not only can more accurate representations of training data be learned, but more linear regions allow the network to better capture underlying nonlinearity to enhance its predictive power in regression tasks. This result is especially significant because it indicates that even in cases where there are fewer data points than linear regions, having the additional regions can still provide performance advantages.

### 4.3 Gradient Descent does not Directly Optimize Efficiency

Figure: Figure 5: Approximation produced with standard methodologies (shown before ReLU is applied). The neuron colors here are arbitrary; they do not correspond to Figure 2.
Refer to caption: x5.png

Figure 5 shows the interior of a default network. The layers here are shown before applying ReLU.
The default networks fail to make efficient use of ReLU to produce linear regions, even falling short of 1 bend per neuron, which can easily be attained by forming a linear spline (1 hidden layer) that interpolates some of the data points.
Examining the figure, the first two layers are wasted.
No neuron’s activation pattern crosses $y=0$, so ReLU is never used.
Layer 3 could be formed directly from the input signal.
Deeper in the network, more neurons remain either strictly positive or negative.
Those that intersect $y=0$ are monotonic, only able to introduce one bend at a time.
The core issue is that while more bends leads to better accuracy, networks that have few bends are not locally connected in parameter space to those that have many.
This is problematic since gradients can only carry information about the effects of infinitesimal parameter modifications.
If a bend exists, gradient descent can reposition it. But for a neuron that always outputs a strictly positive value (such as the red in layer 2), bends cannot be introduced by infinitesimal weight or bias adjustments.
Therefore, bend-related information will be absent from its gradients.
Gradient descent will only compel a network to bend by happenstance; indirectly related local factors must guide a neuron to begin outputting negative values.
Occasionally, these local incentives are totally absent, and the network outputs a bend-free line of best fit.

### 4.4 Effect of Reparameterized Initialization vs Pretraining

Figure: Figure 6: Pretraining (right) results in better structural retention in deeper layers.
Refer to caption: x6.png

Figure 6 compares the top performing reparameterized models for $x^{3}$ with and without pretraining.
We observe that without the guidance of the pretraining, gradient descent usually loses the triangle generating structure around layer 4 or 5, devolving into noisy patterns and resulting in higher errors.
Pretraining maintains structure at greater depths.
This behavior of gradient descent we observe is problematic since theoretical works often rely on specific constructions within networks to prove their results.
Gradient descent greedily abandons any such structure during its optimization in favor of models that can be worse in the long term.
A theoretical result that shows a certain representation exists in the set of neural networks will thus have a hard time actually learning it without a subsequent plan to control training.
This phenomenon highlights the need for theoretical work to derive more expressive training controls, which can give better guidance for the optimization of functions of greater dimension or complexity.

## 5 Extension to Non-Convex Functions and Higher Dimensions

Figure: Figure 7: Approximation of $y=x^{3}-x$ by difference of pretrained components, achieving a loss of $4\times 10^{-7}$. A standard $8\times 5$ network yielded a larger loss of $8\times 10^{-6}$.
Refer to caption: extracted/5905662/x3-x.png

The results presented so far are limited to one-dimensional convex functions, but these are not limitations of our method.
Our constructed networks can be applied to approximate higher-dimensional functions by using these networks as activation functions in a larger network (see Figure 9 in the Appendix for a schematic).
For instance, in two dimensions, we take two or more copies of our network and estimate a target function as a linear combination of the networks’ outputs.
Because the weights of the linear combination are randomized, each copy of the network will tend to learn different lower-dimensional aspects of the target function, whereby the combination can accurately capture the higher-dimensional function.

The same approach can address non-convex functions since many non-convex functions (including every continuous function with a bounded second derivative (Hessian)) can be expressed as the difference of two convex functions (see Zhang et al. (), which also describes how every ReLU network can be decomposed into a difference of two piecewise linear convex functions).
We do not claim that two copies of our network are sufficient for learning every pathological function, but preliminary results in this section are surprisingly promising (and of course, assembling more copies of our network into a deeper or wider larger network would give additional representational power).

We note that the approach in this section bears similarity to the recently popularized Kolmogorov-Arnold Networks (KAN) (Liu et al., ), which use spline functions as activations of a network where every weight is 1.
This is based on the Kolomogorov-Arnold representation theorem (Kolmogorov, ; Arnold, ; ), which gives a guarantee that every continuous multivariate function can be represented by a neural network-like structure with one-dimensional activation functions.
But by restricting the activation functions to splines, the KAN abandons its theoretical guarantees.
Using our networks as activations, however, maintains the theoretical ability to perform polynomial interpolation since they retain the ability to approximate $x^{2}$, and the multiplication gate discussed in Section 2.1 is based on a linear combination of squared terms.

Figures 7 and 8 show learning the non-convex function $y=x^{3}-x$ and the two-dimensional function $z=r^{3}$ on $[-1,1]$, each by using two of our convex function blocks.
Rather than directly using our networks as activation functions, we instead encode each of their layers as a block diagonal structure in one combined matrix (which has the same effect).
Two slight modifications are made to each of the component networks to better interact with ReLU.
First, two neurons are used for the sum so that a negative component may be maintained in the case of $x^{3}-x$, which would otherwise be zeroed out by ReLU.
Second, the triangular functions are generated upside-down, preventing ReLU from zeroing the function outside of the approximation region.

For both $x^{3}-x$ and $r^{3}$, we substantially outperform standard initialization and training procedures.
Especially striking is the increased number of linear regions in Figure 8.
Visible in Figure 7 is the gradient descent-induced destructuring discussed in Section 4.4.
This likely occurs due to the block diagonal structure of the larger weight matrices, which results in many zeros for the second pass of gradient descent to fill in unassisted.
This could be addressed either by limiting training only to the block diagonal entries, or by deriving a more advanced pretraining method that makes use of the full matrix. Note that since $x^{3}-x$ is twice differentiable we can write it as $(x^{3}-x+3x^{2})-3x^{2}$, a difference of functions convex on $[-1,1]$ obtained by finding a parabola based on trying to make the derivative of $x^{3}-x$ monotone. Note that this is not a unique decomposition.

Figure: Figure 8: Approximations of $z=\sqrt{x^{2}+y^{2}}^{3}$ using an $8\times 5$ regular network (left) and a union of two of our pretrained components (right). Losses are $1.5\times 10^{-4}$ and $3.5\times 10^{-6}$, respectively, demonstrating a nearly two orders of magnitude improvement using our techniques.
Refer to caption: extracted/5905662/lump.png

## 6 Concluding Remarks

This paper focused on exploiting the potential computational complexity advantages neural networks offer for the problem of efficiently learning nonlinear functions; in particular, compelling ReLU networks to approximate functions with exponential accuracy as network depth is linearly increased.
Our results showed improvements of one to several orders of magnitude in using our initialization and pretraining strategy to train ReLU networks to learn various nonlinear functions, including non-convex and multi-dimensional functions. This finding is particularly powerful since random initialization and gradient descent are not likely to produce an efficient solution on their own, even if it can be proven to exist in the set of sufficiently sized ReLU networks.
Although we proposed one preliminary strategy for extending our network for higher-dimensional functions,
this strategy has limitations, such as leaving many weights to be filled in by gradient descent (which it generally does not do effectively).
We anticipate continuing investigating strategies for improving the generalizability of our work, which we believe is critically important since it offers the possibility of an exponentially more accurate drop-in replacement of linear layers in any architecture.
We are hopeful that future works by our group and others will help illuminate a complete theory for harnessing the potential exponential power of depth in ReLU or even other networks.

## Appendix A Appendix

Figure: Figure 9: Extending our results into a larger network. ReLU should not be used in the layers between the blocks, as creating nonconvex functions this way will tend to make a lot of negative values.
Refer to caption: x8.png

### A.1 Initialization and Training Algorithm

The initialization step of our algorithm is to generate a vector $A=[a_{0},a_{1},...a_{n}]^{T}$, where each $a_{i}$ is randomized in $(0,1)$.
Given this $A$, the pretraining step of our algorithm sets the weights of the input ($I$), hidden ($H_{i}$, $1\leq i\leq n-1$), and output ($O$) layers of the network as follows:

$$ $\displaystyle I(x)=\begin{bmatrix}x\\ x\\ x\\ 0\end{bmatrix}+\begin{bmatrix}0\\ 0\\ -a_{0}\\ 1\end{bmatrix}$ $$

$$ $\displaystyle H_{i}(x)=\begin{bmatrix}1&S_{i}/a_{i}&-S_{i}/(a_{i}-a_{i}^{2})&0 \\ 0&S_{i}/a_{i}&-S_{i}/(a_{i}-a_{i}^{2})&0\\ 0&S_{i}/a_{i}&-S_{i}/(a_{i}-a_{i}^{2})&-S_{i}a_{i+1}\\ 0&0&0&S_{i}\end{bmatrix}\times ReLU(H_{i-1}(x))$ $$

$$ $\displaystyle O(x)=\begin{bmatrix}1&S_{n}/a_{n}&-S_{n}/(a_{n}-a_{n}^{2})&0\end {bmatrix}\times ReLU(H_{n-1}(x))$ $$

where $S_{i}$ can either be chosen independently, or chosen based on $A$. In the latter case, Equation 4 gives $S_{i}=s_{i}/s_{i-1}=(1-a_{i})a_{i+1}$. This assignment of $I$, $H$, and $O$ is used in each iteration of the pretraining algorithm.

Figure: Algorithm 1 Initialization and Pretraining

The network weights can then be initialized once more based on the learned vector $A$ so that the network weights can be updated via regular gradient descent. In our experiments, both phases of training (pretraining and gradient descent) ran for 1000 epochs.
Full pseudocode is listed in Algorithm 1.

### A.2 Derivations

For convenience we first restate the functions defined in the main body of the paper.

$$ $\displaystyle T_{i}(x)=\begin{cases}\frac{x}{a_{i}}&0\leq x\leq a_{i}\\ 1-\frac{x-a_{i}}{1-a_{i}}&a_{i}\leq x\leq 1\end{cases}$ $$

$$ $T^{\prime}_{i}(x)=\begin{cases}\frac{1}{a_{i}}&0<x<a_{i}\\ \frac{1}{1-a_{i}}&a_{i}<x<1\end{cases}$ $$

$$ $W_{i}(x)=\bigcirc_{j=0}^{i}T_{j}(x)=T_{i}(T_{i-1}(...T_{0}(x)))$ $$

$$ $F(x)=\sum_{i=0}^{\infty}s_{i}W_{i}(x)$ $$

Notationally, we will denote the sorted $x$-locations of the peaks and valleys of $W_{i}(x)$ by the lists $P_{i}=\{x:W_{i}(x)=1\}$ and $V_{i}=\{x:W_{i}(x)=0\}$. We will use the list $B_{i}$ to reference the locations of all non-differentiable points, which we refer to as bends. $B_{i}:=P_{i}\cup V_{i}$. $f_{i}(x)=\sum_{n=0}^{i-1}s_{n}W_{n}(x)$ will denote finite depth approximations up to but not including layer $i$. The error function $E_{i}(x)=\sum_{n=i}^{\infty}s_{n}W_{n}(x)=F(x)-f_{i}(x)$ will represent the error between the finite approximation and the infinite depth network. Although the derivative of an infinite sum of functions is not always equal to the sum of the individual functions derivatives, it is in this case, so notationally we use $F^{\prime}$ to describe the latter for brevity.

Figure: Figure 10: Triangle functions $T_{0}$ and $T_{1}$, and the triangle wave resulting from their composition. Note how $T_{1}$ is reflected in $W_{1}$.
Refer to caption: x9.png

The goal of this section will be to show how to select the $s_{i}$ based on $a_{i}$ in a manner where the derivative $F^{\prime}(x)$ is defined and on all of $[0,1]$. Fortunately, the left and right derivative limits $F^{\prime}_{+}(x)$ and $F^{\prime}_{-}(x)$ already exist, since each bend in each $W_{i}$ will use the slope of the line segment to its left or right respectively. The problem is that without choosing the scaling parameters appropriately, $F^{\prime}_{+}(x)$ and $F^{\prime}_{-}(x)$ will take different values so that $F^{\prime}(x)$ will be undefined on all bend points.

Figure 10 highlights some important properties about composing triangle functions. Peaks alternate with valleys. Peak locations in one layer become valleys in the next. Valleys in one layer remain valleys in all future layers since 0 is a fixed point of each $T_{i}$. To produce $W_{i}$, each line segment of $W_{i-1}$ becomes a dilated copy of $T_{i}$. Each triangle function has two distinct slopes $1/a_{i}$ and $-1/(1-a_{i})$ which are dilated by the chain rule during the composition. On negative slopes of $W_{i-1}$, the input to layer $i$ is reversed, so those copies of $T_{i}$ are reflected. Due to the reflection, the slopes of $W_{i}$ on each side of a peak or valley are proportional. Alternatively, one could consider that on each side of a peak in $W_{i-1}$, there is a neighborhood of points that are greater than $a_{i}$, and are composed with the same line segment of $T_{i}$ that has slope $-1/(1-a_{i})$. Either way, it’s important to note that the slopes on each side of a bend scale identically during each composition.

Before we begin reasoning about $F^{\prime}(x)$, it can simplify the analysis to only consider the derivative of the error function $E^{\prime}(x)$.

###### Lemma A.1 .

for $x\in P_{i}$, $F^{\prime}(x)$ is defined if and only if $E^{\prime}_{i}(x)$ is defined.

###### Proof.

All of $W_{n}(x)$ for $n<i$ are differentiable at $x\in P_{i}$ since $x$ will lie in the interior of a linear region of $W_{n}$. Therefore, $f^{\prime}_{i}(x)=\sum_{n=0}^{i-1}s_{n}W^{\prime}_{n}(x)$ exists at these points. Since $E^{\prime}_{i}+f^{\prime}_{i}=F^{\prime}$, $F^{\prime}(x)$ is defined if and only if $E^{\prime}_{i}(x)$ is defined.
∎

Here we compute the right derivative $(E_{i})^{\prime}_{+}(x)$ of the error function at a point $x$. The left derivative will only be different by a constant factor.

###### Lemma A.2 .

For all $x\in P_{i}$, $E_{+}^{\prime}(x)$ and $E_{-}^{\prime}(x)$ are proportional to

$$ $s_{i}-\frac{1}{1-a_{i+1}}\left(s_{i+1}+\sum_{n=i+2}^{\infty}s_{n}\prod_{k=i+2} ^{n}\frac{1}{a_{k}}\right).$ (5) $$

###### Proof.

Let $x_{n}$ be some point in $P_{i}$, and let $n$ be its index in any list it appears in. To calculate the value of $E_{+}^{\prime}(x)=\sum_{n=i}^{\infty}s_{n}(W_{n})^{\prime}_{+}(x)$, we will have to find the slope of the linear intervals to the immediate right of $x_{n}$ for all $W_{i}$. The terms in the summation Equation 5 mostly derive from the chain rule. We will use $R$ to represent $W_{i+}^{\prime}(x)$. The first term in the sum will be $Rs_{i}$. $W^{\prime}_{i+1}(x)=T^{\prime}_{i+1}(W_{i}(x))W^{\prime}_{i}(x)$ since the derivatives of composed functions will multiply. There are two different possible slope values of $T_{i+1}$ to multiply by, and the correct one to multiply is $-1/(1-a_{i+1})$ because $x_{n}$ is a peak of $W_{i}$, so $W_{i}(x)>a_{i+1}$ for $x\in(B_{i+1}[n-1],B_{i+1}[n+1])$. This gives $W^{\prime}_{i+1}(x)=-R\frac{s_{i+1}}{(1-a_{i+1})}$. Note that the second term has the opposite sign as the first.

For all remaining terms, since $x_{n}$ was in $P_{i}$, it is in $V_{j}$ for $j>i$. For $x\in(B_{j+1}[n-1],B_{j+1}[n+1])$, $W_{j}(x)<a_{j+1}$ and the chain rule applies the first slope $1/a_{j+1}$. Since this slope is positive, every term has the opposite sign as the first. Summing up all the terms with the coefficients $s_{i}$, and factoring out $R$ will yield the desired formula. Note that this same reasoning applies to the left side, the initial slope $L$ will just be different.
∎

###### Lemma A.3 .

If $E_{+}^{\prime}(x)=E_{-}^{\prime}(x)$, $E^{\prime}(x)$ must be equal to 0.

###### Proof.

Let $S$ represent Equation 5, and $R$ and $L$ be the constants of proportionality for the directional derivatives. If $E^{\prime}_{+}=E^{\prime}_{-}$, then $R_{x}S=L_{x}S$ for all $x\in P_{i}$. Since $W_{i}$ is comprised of alternating positive and negatively sloped line segments, $R_{x}$ and $L_{x}$ have opposite signs. The only way to satisfy the equation then is if $S=0$. Consequently, $E^{\prime}(x)=0$ for all $x\in P_{i}$.
∎

This lemma shows that to calculate the derivative at of $F(x)$ for any bend point $x$ (whenever it’s actually defined), one needs only to compute the derivative of the finite approximation $f_{i}$ (which excludes $W_{i}$). This will be useful later for proving other results. A visualization of this lemma is provided in Figure 11

###### Lemma A.4 .

For all $x\in P_{i}$:

$$ $F^{\prime}(x)=f^{\prime}_{i}(x)=\sum_{j=0}^{i-1}s_{j}W^{\prime}_{j}(x)$ (6) $$

###### Proof.

From the previous lemma we have $E^{\prime}(x)=0$ whenever the directional derivatives are equal. $F(x)=\sum_{j=0}^{i-1}s_{j}W_{j}(x)+E(x)$. The first $i-1$ terms are differentiable at the points $P_{i}$ since those points lie between the discontinuities in $B_{i-1}$. Therefore $F^{\prime}(x)$ is defined and can be calculated using the finite sum.
∎

We now prove our main theorem.
Intuitively, the theorem shows that there is a way to sum the triangular waveforms $W_{i}$ so that the resulting function approximation converges to a differentiable function, which, as mentioned at the start of this subsection, can aid in preventing overfitting when using exponentially many linear regions.
The idea of the proof is that much of the formula for $E^{\prime}(x)$ will be shared between two successive generations of peaks.
Once they are both valleys, they will behave the same, so the sizes of their remaining discontinuities will need to be proportional.

###### Theorem (3.1) .

$F^{\prime}(x)$ is defined on $[0,1]$ only if the scaling coefficients are selected based on $a_{i}$ according to:

$$ $s_{i+1}=s_{i}(1-a_{i+1})a_{i+2}$ $$

###### Proof.

Rewriting Equation 5 (which is equal to 0) for layers $i$ and $i+1$ in the following way:

$$ $s_{i}(1-a_{i+1})=s_{i+1}+\frac{1}{a_{i+2}}\left(s_{i+2}+\sum_{n=i+3}^{\infty}s _{n}\prod_{k=i+3}^{n}\frac{1}{a_{k}}\right)$ $$

$$ $s_{i+1}(1-a_{i+2})=s_{i+2}+\sum_{n=i+3}^{\infty}s_{n}\prod_{k=i+3}^{n}\frac{1} {a_{k}}$ $$

allows for a substitution to eliminate the infinite sum

$$ $s_{i}(1-a_{i+1})=s_{i+1}+\frac{1-a_{i+2}}{a_{i+2}}s_{i+1}$ $$

Collecting all the terms gives

$$ $s_{i+1}=\frac{s_{i}(1-a_{i+1})}{1+\frac{1-a_{i+2}}{a_{i+2}}}$ $$

which simplifies to the desired result.
∎

###### Lemma A.1 .

###### Proof.

###### Lemma A.2 .

###### Proof.

###### Lemma A.3 .

###### Proof.

###### Lemma A.4 .

###### Proof.

###### Theorem (3.1) .

###### Proof.

### A.3 Sufficiency for differentiability

We can show that in addition to being necessary for the existence of the derivative on the bend points, our choice of scaling is sufficient when $a_{i}$ are bounded away from 0 or 1.

###### Theorem A.5 .

$F^{\prime}(x)$ is well-defined if we can find $c>0$ such that $c\leq a_{i}\leq 1-c$ for all $i$ and choose all $s_{i}$ according to Equation 4.

###### Proof.

We begin by considering Equation 5 for layer $i$ (it equals 0 by Theorem A.3).

$$ $s_{i}=\frac{1}{1-a_{i+1}}\left(s_{i+1}+\sum_{n=i+2}^{\infty}s_{n}\prod_{k=i+2} ^{n}\frac{1}{a_{k}}\right)$ $$

We will prove our result by substituting Equation 4 into this formula, and then verifying that the resulting equation is valid. First we would like to rewrite each $s_{n}$ in terms of $s_{i}$. Equation 4 gives a recurrence relation. Converting it to an non-recursive representation we have:

$$ $s_{n}=s_{i}\left(\prod_{j=i+1}^{n}1-a_{j}\right)\left(\prod_{k=i+2}^{n+1}a_{k}\right)$ (7) $$

When we substitute this into Equation 5, three things happen: each term is divisible by $s_{i}$ so $s_{i}$ cancels out, every factor in the product except the last cancels, and $1-a_{i+1}$ cancels. This leaves

$$ $1=a_{i+2}+(1-a_{i+2})a_{i+3}+(1-a_{i+2})(1-a_{i+3})a_{i+4}+...=\sum_{n=i+2}^{ \infty}a_{n}\prod_{m=i+2}^{n-1}(1-a_{m})$ (8) $$

This equation has an interpretation that is important to the argument. 1 is the full size of the initial discontinuity at a point in $P_{i}$, and each term on the other side represents proportionally how much the discontinuity is closed for each triangle wave that is added. Every time a wave is added, it can be thought of as subtracting from each side the first term appearing on the right. We will now show that each term of the sum on the right accounts for a fraction (equal to $a_{i}$) of the remaining error. Inductively we can show:

$$ $1-\sum_{n=i+2}^{j}a_{n}\prod_{m=i+2}^{n-1}(1-a_{m})=\prod_{m=i+2}^{j}(1-a_{m})$ (9) $$

This means that as the first term appearing on the right is repeatedly subtracted, that term is always equal to $a_{n}$ times the left side. As a base case, we have $(1-a_{i+2})=(1-a_{i+2})$. Assuming the above equation holds for all previous values of $j$

$$ $\displaystyle 1-\sum_{n=i+2}^{j+1}a_{n}\prod_{m=i+2}^{n-1}(1-a_{m}))=1-\sum_{n =i+2}^{j}a_{n}\prod_{m=i+2}^{n-1}(1-a_{m}))-a_{j+1}\prod_{m=i+2}^{j}(1-a_{m}))=$ $$

using the inductive hypothesis to make the substitution

$$ $\displaystyle\prod_{m=i+2}^{j}(1-a_{m}))-a_{j+1}\prod_{m=i+2}^{j}(1-a_{m}))= \prod_{m=i+2}^{j+1}(1-a_{m}))$ $$

Since all $c<a_{i}<1-c$, the size of the discontinuity at the points $P_{i}$ is upper bounded by the exponentially decaying series $(1-c)^{n}$, which approaches zero.
∎

###### Theorem A.5 .

###### Proof.

### A.4 Error Decay

###### Lemma A.6 .

The ratio $s_{i+2}/s_{i}$ is at most $0.25$.

###### Proof.

by applying Equation 4 twice, we have

$$ $\displaystyle s_{i+2}=s_{i}(1-a_{i+1})(1-a_{i}+2)a_{i+2}a_{i+3}$ $$

To maximize $s_{i+2}$ we choose $a_{i+1}=0$ and $a_{i+3}=1$. The quantity $a_{i+2}-a_{i+2}^{2}$ is a parabola with a maximum of $0.25$ at $a_{i+2}=0.5$.
∎

Since each $W_{i}$ takes values between 0 and 1, its contribution to $F$ is bounded by $s_{i}$. Since the $s_{i}$ decay exponentially, one could construct a geometric series to bound the error of the approximation and arrive at an exponential rate of decay.

###### Lemma A.6 .

###### Proof.

### A.5 Second Derivatives

Here we show that any function represented by one of these networks that is not $y=x^{2}$ does not have a second derivative defined at the bend locations. To show this we will sample a discrete series of $\Delta y/\Delta x$ values from $F^{\prime}(x)$ and show that the limits of these series on the right and left are not the same (unless all $a_{i}=0.5$), which implies that $F^{\prime\prime}(x)$ does not exist (see Figure 11 below). First we will produce the series of $\Delta x$. Let $x$ be the location of a peak of $W_{i}$, and let $l_{n}$ and $r_{n}$ be its immediate neighbors in $B_{i+n}$.

###### Lemma A.7 .

If $c<a_{i}<1-c$ for all $i$, we have $\lim_{n\rightarrow\infty}r_{n}=\lim_{n\rightarrow\infty}l_{n}=x$. Furthermore, $r_{n},l_{n}\neq x$ for any finite $i$.

###### Proof.

Let $R$ and $L$ denote the magnitude of $W^{\prime}_{i}$ on the left and right of $x$. $x$ is a peak location of $W_{i}$, so the right side slope is negative and the left is positive. Solving for the location of $T_{i+1}(W_{i}(x))=1$ on each side will give $l_{1}=x-(1-a_{i+1})/L$ and $r_{1}=x+(1-a_{i+1})/R$.

On each subsequent iteration $i+n$ ($n\geq 2$), $x$ is a valley point and the $\Delta x$ intervals get multiplied by $a_{i+n}$. Since $x$ is a valley point the right slope is positive and the left is negative. The slope magnitudes are given by $\frac{1}{x-l_{n}}$ and $\frac{1}{r_{n}-x}$ since $W_{i+n}$ ranges from 0 to 1 over these spans. Solving for the new peaks again will give $l_{n+1}=x-a_{i+1}(x-l_{n})$ and $r_{n+1}=x+a_{i+1}(r_{n}-x)$. The resulting non-recursive formulas are:

$$ $x-l_{n}=\frac{1-a_{i+1}}{L}\prod_{m=2}^{n}a_{i+m}\text{ and }r_{n}-x=\frac{1-a _{i+1}}{R}\prod_{m=2}^{n}a_{i+m}$ (10) $$

The right hand sides will never be equal to zero with a finite number of terms since $a$ parameters are bounded away from 0 0 and $1$ by $c$.
∎

Figure: Figure 11: The derivatives of the first few stages of approximation. Notice that each time a constant segment“splits”, the two neighboring segments adjacent to the split converge back to the original value. The points in the series approaching $x$ from the right are marked. $a_{0}\neq 0.5$ and all other parameters are set to $0.5$, which will cause the left and right sets of points to lie on lines with different slopes.
Refer to caption: x10.png

Next we derive the values of $\Delta y$ to complete the proof.

###### Theorem A.8 .

$F(x)$ cannot be twice differentiable unless $F(x)=x^{2}$.

###### Proof.

The points $l_{n}$ and $r_{n}$ are all peak locations, Equation 6 gives their derivative values as $f^{\prime}_{i+n}(r_{n})$. In our argument for sufficiency, we reasoned about the sizes of the discontinuities in $f^{\prime}$ at $x$. Since $l_{n}$ and $r_{n}$ always lie on the linear intervals surrounding $x$ as $n\rightarrow\infty$, we can get the value of $f^{\prime}_{i}(x)-f^{\prime}_{i+n}(r_{n})$ using Equation 9 with the initial discontinuity size set to $R*s_{i}$ rather than 1. Focusing on the right hand side we get:

$$ $f^{\prime}_{i}(x)-f^{\prime}_{i+n}(r_{n})=R*s_{i}\prod_{m=2}^{n}(1-a_{i+m})$ $$

taking $\Delta y/\Delta x$ gives a series:

$$ $\frac{R^{2}s_{i}}{(1-a_{i+1})}\prod_{m=i+2}^{n}\frac{1-a_{m}}{a_{m}}$ $$

The issue which arises is that the derivation on the left is identical, except for a replacement of $R$ by $L$. The only way for these formulas to agree then is for $R=L$ which implies $a_{i}=1-a_{i}=0.5$. Since this argument applies at any layer, then all $a$ parameters must be $0.5$ (which approximates $y=x^{2}$.
∎

###### Lemma A.7 .

###### Proof.

###### Theorem A.8 .

###### Proof.

### A.6 Concavity

Each of the $f^{\prime}_{i}$ are composed of constant value segments, we will show that those values are monotonically decreasing (this can be seen in Figure11). This then allows us to show that $F$ is concave.

###### Lemma A.9 .

The function $F(x)$ is concave when all $s_{i}$ are chosen according to Equation 4.

###### Proof.

To establish this we will introduce the list $Y^{\prime}_{i}=[F^{\prime}(V_{i}[0]),f^{\prime}_{i}(V_{i}[n]),F^{\prime}(V_{i+
1}[2^{i}]]$ for $0\leq n\leq 2^{i}$, which tracks the values of $F^{\prime}$ at the $i$^th set of valley points. All but the first and last points will have been peaks at some point in their history, so Equation 6 gives the value of those derivatives as $f^{\prime}_{i}$.

We establish two inductive invariants. One is that the y-values in the list $Y_{i}$ remain sorted in descending order. The other is that $Y^{\prime}_{i}[n]\geq f^{\prime}_{i}(x)\geq Y^{\prime}_{i}[n+1]$ for $V_{i}[n]<x<V_{i}[n+1]$, indicating that the constant value segments of $f_{i}$ lie in between the limits in the list $Y_{i}$. Together, these two facts imply that each iteration of the approximation $f_{i}$ is concave. Which then can prove that their limit $F$ is also concave.

As the base case $f_{0}$ is a line with derivative 0, $V_{0}$ contains its two endpoints. $Y^{\prime}_{0}$ is positive for the left endpoint (negative for right) since on the far edges $F^{\prime}$ is a sum of a series of positive (or negative) slopes, Therefore both the points in $Y^{\prime}$ are in descending sorted order. The second part of the invariant is true since 0 is in between those values.

Consider an arbitrary interval $(V_{i}[n],V_{i}[n+1])$ of $f_{i}$, this entire interval is between two valley points, so $f^{\prime}_{i}$ (which hasn’t added $W_{i}$ yet) is some constant value, which we know from the second inductive hypothesis is in between $Y^{\prime}_{i}[n]$ and $Y^{\prime}_{i}[n+1]$. The point $x\in P_{i}\cap(V_{i}[n],V_{i}[n+1])$ will have $F^{\prime}(x)=f^{\prime}_{i}(x)$, and it will become a member of $V_{i+1}$. This means we will have $Y_{i+1}[2n]>Y_{i+1}[2n+1]>Y_{i+1}[2n+2]$, maintaining sorted order of $Y^{\prime}$.

Adding $s_{i}W_{i}$ takes $f_{i}$ to $f_{i}+1$ splitting each constant valued interval in two about the points $P_{i}$, increasing the left side, and decreasing the right side. Recalling from the derivation of Equation 5 all terms but the first in the sum have the same sign, so the limiting values in $Y^{\prime}_{i}$ are approached monotonically. Using the first inductive hypothesis, we have on the left interval $Y^{\prime}_{i}[n]=Y^{\prime}_{i+1}[2n]>f^{\prime}_{i+1}>f^{\prime}_{i}=Y^{
\prime}_{i+1}[2n+1]$ and on the right we have $f^{\prime}_{i}=Y^{\prime}_{i+1}[2n+1]>f_{i+1}>Y^{\prime}_{i}[n+1]=Y^{\prime}_{
i}[2n+2]$. And so each constant interval $f_{i+1}$ remains bounded by the limits in $Y^{\prime}_{i+1}$.

We will now show by contradiction that the limit $F$ of the sequence of concave $f_{i}$ is also concave. Assume that $F$ is non-concave. Then there exist points $a$, $b$, and $c$ such that $F(b)$ lies strictly below the line connecting the points $(a,F(a))$ and $(c,F(c))$. Lets say it’s below the line by an amount $\epsilon$. Since at each point $f_{i}$ converges to $F$, we can find $i_{a}$ such that $f_{ia}(a)-F(a)<\epsilon/2$, etc…, we take $i=\max(i_{a},i_{b},i_{c})$. Since $f_{i}(a)$ and $f_{i}(c)$ are no more than $\epsilon/2$ lower than their limiting values, the entire line connecting $(a,f_{i}(a))$ and $(c,f_{i}(c))$ is no more than $\epsilon/2$ lower than the line between $(a,F(a))$ and $(c,F(c))$. $f_{i}(b)$ is also no more than $\epsilon/2$ higher than $F(b)$, thus $f_{i}(b)$ must still lie below the line between $(a,f_{i}(a))$ and $(c,f_{i}(c))$, making $f_{i}$ non-concave and producing a contradiction.
∎

###### Lemma A.9 .

###### Proof.

### A.7 Learning Rates

Figure: Figure 12: Learning $y=x^{3}$ at various learning rates for 1000 epochs. The figure shows the log losses of the best model of 30 for networks trained according to the methods in this paper versus random initialization and gradient descent (“Default Network”). For high learning rates, neither learns because the steps are too large. For low learning rates (below $10^{-5}$), the steps are too small, and our networks are likely deriving an implicit advantage from being forcibly initialized to a convex function. Both methods are able to converge for learning rates in between $10^{-4}$ and $10^{-2}$; one could run for more epochs to see a similar advantage of our method for smaller learning rates, if desired.
Refer to caption: x11.png

All results in the main body of the paper used a constant learning rate of $10^{-3}$.
In this appendix, we considered an ablation study on the learning rate for the task of learning $y=x^{3}$. As seen in Figure 12, the learning rate we selected was approximately optimal for both our method as well as default network training.
We note that for this ablation study, a constant 1000 epochs were run, which explains why both methods perform worse as the learning rate becomes minuscule. At small learning rates, what is really measured is how close of a guess the initialization is to the target function. Our networks are preforming better here simply because they are always outputting convex functions. But this only accounts for losses on the order of $10^{-4}$, which indicates that in the more reasonable learning rate ranges, our pretraining is performing a meaningful function and enabling order-of-magnitude improvements over default network training.

## Potential Broader Impact

This paper presents work whose goal is to enable more efficient neural networks. While the present work is largely theoretical, future advances in this line of research could enable the use of much smaller networks in many practical applications, which could substantially mitigate the rapidly growing issue of energy usage in large learning systems.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: Please address the Weaknesses.
soundness: 2
strengths: -	The paper presents detailed introduction and explanation of the proposed method, including how to construct the initialisation and how to calculate the gradient. The paper is compressive, well-structed, and easy to follow. 
-	The paper present experiments for cases of one dimension and high-dimension non-convex problems.
-	I find the paper makes conceptual contributions of proposing a new initialisation strategy.
summary: This paper designs a novel training strategy: (1) reparameterise the network weights in to make it exhibit a number of linear regions exponential in depth; (2) train on the derived parameters for an initial solution; (3) refine the parameter by directly updating the underlying model weights. Experiments are given to support the method.
weaknesses: -	The experiments are not sufficient. The current experiments only cover quite shallow (three layers) ReLU neural networks on very simple tasks. It is unclear whether the results apply to complex scenarios, like deeper neural networks, transformer on fitting images, mining on text data, etc. Thus, the paper actually cannot help understand the success of deep learning.
- No comparison is given with other initialisation methods.
- The explanation of why this method works is not sufficient. This makes the method not convincing.
-	No theoretical results are provided. This is particularly severe given the experiments are insufficient.

The paper looks like in an early stage with insufficient validation. I suggest the authors to do more following on this paper.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: My review is as follows:

- I think this paper brings up important issues with standard neural network training practices (such as using relu for activation or gradient descent, etc).

- One thing that I think is potentially missing is the verification of the findings on a somewhat more realistic scenario. Could we expect the proposed method to outperform a standard neural network approach (e.g. a similar size relu network trained by SGD) when, say, predicting airline delays? Or, other more standard methods such as linear regression or decision tree?

- To make the results potentially more broad, I wonder if the proposed strategy could be somehow applied to classification (perhaps can test it on a simple dataset such as "two spirals" dataset). If that's not straightforward, I think it'd still help me understand the contributions better if the authors can comment on the challenges.

- Which of the findings of this paper could we expect to carry over to other non-linear activation functions such as sigmoid?
soundness: 3
strengths: Please see the Questions section.
summary: This paper focuses on the expressivity of ReLU networks and argues that the standard neural network training approaches lead to models that cannot utilize all of the linear regions that a ReLU network has the potential to exhibit. The paper contains an approach to overcome this issue.
weaknesses: Please see the Questions section.

---

confidence: 3
contribution: 3
initial_rating: 3
presentation: 3
questions: 1. How are ReLU activations guaranteed to generate symmetric triangle waves? There are many possible compositions of ReLU activations, but only a subset of them are symmetric triangle waves. If the network is reparameterized using only the triangle wave basis functions as proposed in the paper, will it lose some flexibility and expressivity as it is not possible for the reparameterized network to create other shapes or patterns within each layer? 

2. Could you provide some intuitions and/or theory regarding why pretraining helps maintain the triangle generating structure and avoid eliminating activation regions as the network gets deeper?

3. Does the proposed method improve convergence rate? Could you demonstrate it with some experiments and/or theory?
soundness: 3
strengths: 1. The proposed reparameterization of ReLU networks is novel and interesting. Based on the observation that ReLU networks can generate symmetric triangle waves, the proposed approach introduces an approach to directly reprent the ReLU neurons (rather than weights) as asymmetric triangle wave basis function with learnable locations of the peaks within [0, 1].
2. The proposed learning algorithm is simple and seems to be effective in training the reparameterized network for simple target functions as shown in the experiments.
3. The 1D demonstrations and experiments are nice, which illustrate how the proposed method learns useful patterns for fitting the target function.
summary: This paper proposes to reparameterize ReLU networks by parameterizing the peaks of the triangle wave basis functions generated by ReLU activations. This ensures that the number of linear regions grows exponentially with the depth of the network, which reduces the waste of representation capacities in randomly initialized ReLU networks. A learning algorithm is proposed to train the reparameterized ReLU network by first updating the derived parameters and then updating the actual weights underlying the model. The proposed method is empirically evaluated on 1D convex and 2D nonconvex target functions, which demonstrate its improved accuracy compared to randomly initialized networks.
weaknesses: 1. Theorem 3 seems to be an “only if” statement, and therefore setting $s_{i+1}$ according to Eq 4 is a necessary but not sufficient condition to guarantee differentiability of the reparameterized network.

2. It is unclear how useful the proposed method is in practice, since it is mostly evaluated on simple 1D convex target functions. As this is not a theory paper, the proposed method should at least be evaluated on some semi-real datasets, (e.g., some of the common UCI datasets).

3. Writing needs to be improved. The logic flow is a bit confusing. Also, several important concepts and building blocks are not well explained in the paper. For example, important definitions like definitions of linear regions, activation patterns, and activation regions should be stated in the paper or appendix.

4. The following closely related works which analyze compositions and/or reparameterization of ReLU activations are not discussed in the paper.

[1] K Eckle, J Schmidt-Hieber. A comparison of deep networks with relu activation function and linear spline-type methods. Neural Networks 2019.

[2] DM Elbrächter, J Berner, P Grohs. How degenerate is the parametrization of neural networks with the ReLU activation function? NeurIPS 2019.

[3] W Chen, H Ge. Neural characteristic activation analysis and geometric parameterization for ReLU networks. NeurIPS 2024.

[4] B Hanin, D Rolnick. Complexity of linear regions in deep networks. ICML 2019.

[5] M Raghu, B Poole, J Kleinberg, S Ganguli, J Sohl-Dickstein. On the expressive power of deep neural networks. ICML 2017.

[6] D Rolnick, K Kording. Reverse-engineering deep ReLU networks. ICML 2020.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: * How does the proposed reparameterization strategy perform in complex, high-dimensional tasks, and what are the challenges in scaling this method effectively to more realistic datasets?

* Have you considered testing this method on real-world datasets with more variability and noise? How robust is the technique in such scenarios, and are there any performance trade-offs when dealing with non-synthetic data?
soundness: 3
strengths: * The paper introduces a novel approach to reparameterize ReLU network weights, which forces the network to exhibit an exponential number of activation regions. This significantly enhances the expressivity of the network and addresses the inefficiencies of randomly initialized models, providing a more accurate and efficient approximation of nonlinear functions.

* The proposed pretraining strategy allows the network to initialize with exponentially more linear regions, thus reducing the reliance on gradient descent to discover new activation patterns. This results in faster convergence and much more accurate function approximations, as demonstrated through numerical experiments, showing orders of magnitude lower errors compared to standard initialization methods.
summary: This paper presents a novel strategy to improve the efficiency of ReLU neural networks. It focuses on overcoming the limitations of randomly initialized networks, which tend to be unnecessarily large and inefficient in approximating simple functions. The authors introduce a reparameterization of network weights that ensures an exponential number of activation patterns, thus maximizing the linear regions in the input space. Their approach includes a pretraining stage using derived parameters that enhances the expressivity of the network before standard gradient descent is applied. This method shows significant improvement in approximating both convex and non-convex functions, with better accuracy and efficiency compared to traditional networks. The paper's findings demonstrate that networks initialized with exponential linear regions can capture nonlinearity more effectively, leading to more accurate function approximations. It concludes with potential extensions to multidimensional and non-convex functions, positioning this strategy as a promising tool for more efficient deep learning models.
weaknesses: * While the paper demonstrates improvements in one-dimensional convex functions, the results for higher-dimensional functions and complex non-convex problems are not as thoroughly explored. The proposed method may face scalability challenges when extending to high-dimensional inputs, where the complexity of real-world tasks lies.

* The introduction of a pretraining step with specific reparameterization adds complexity to the network training pipeline. This may make the approach more difficult to implement or integrate into standard deep learning workflows, especially for practitioners looking for more straightforward techniques.

* The effectiveness of the method heavily relies on carefully derived theoretical constructs, such as triangle functions and their parameterization. While this works well in controlled scenarios, its practical robustness in more diverse and noisy real-world datasets is not fully tested or demonstrated.

Meta-review: This paper proposes a new training method that pretrains a ReLU network to a parameter having exponentially many linear regions before switching to naive gradient descent updates, which is demonstrated to have good approximation capabilities compared to standard initialization methods such as Kaiming initialization or RAAI. 

This submission has gone through a lot of activity during the discussion period. Since the initial submission only contained 1-d/2-d experiments, all reviewers expressed concerns whether the proposed algorithm will be useful in practical high-dimensional datasets. The authors updated the paper with results on a 8-d UCI dataset and a 2-d two-spiral dataset (as well as other edits), and also presented some preliminary experiments on CIFAR-10 and ImageNet datasets (albeit trained only for a few epochs).

After the revision and the rebuttal, many reviewers updated their scores and the current average numerical score sits at the decision boundary. After a careful reading of the paper, reviews, and the discussion, my final evaluation is that the paper is unfortunately placed slightly below the bar for acceptance at ICLR.

The paper proposes a new training algorithm applicable to any linear layers, with the guarantee-by-design (at least during the pretraining stage) that the network parameters stay in the region with many “linear pieces” in the input space. I believe this can be a nice idea that has a good potential to develop into a new and efficient training procedure for ReLU networks in general. However, whenever one challenges the “de facto standard”, one needs to be equipped with strong theoretical guarantees on the performance and/or extensive experiments demonstrating the superiority of the proposed method. In my honest opinion, the submission at the moment fails to offer neither of them.

On the **theory** side, I think the presented theoretical results do not successfully demonstrate the advantage of the proposed algorithm. The theory is mostly on the characteristics (differentiability, concavity etc) of the function $F$ constructed by the infinite depth network. The network is guaranteed to have exponentially many linear regions in the input space up to the pretraining stage, but there is no theoretical guarantee after this stage: the training dynamics, convergence rates, generalization bounds, etc in the naive GD update stage are not considered. One can only intuitively expect that starting at a pretrained parameter with exponentially many number of linear piece should be helpful for a better approximation in the final stage. 

In light of this, I must say that at least some of the authors’ claims regarding their theoretical guarantees are overly strong. For example,
- **“the theoretical results that explain why our method works are all proven”**: does not seem convincing from a learning theorist’s perspective.
-  **“theoretical guarantees are for convex functions and non-convex functions that are differences of convex functions”**: what kind of guarantees? In order to claim such a thing, we expect some universal approximation-type theorems. However, in the paper, I only see “any $F(x)$ with $s_i$ chosen according to Eq 3.1 is concave”, not something like “any convex function $g: [0,1] \to \mathbb{R}$ with $g(0) = 0$ can be approximated by $F$.”

On the **experiments** side, as all reviewers pointed out, in order to make significant impact, the proposed method must be tested on higher-dimensional real datasets. Although the authors put extensive efforts to address reviewers' concerns on the practical usefulness, I still wonder if the new results are significant. 

In the new Appendix A.9, it seems to me that the “advantages” of the proposed algorithm over the ordinary training can only be seen when we restrict training to block diagonal parameters. The authors assert this is for fair comparison, but it's hard for me to see why we should restrict ourselves in a low-performing regime (block diagonal variants have higher test losses for both methods). Also, in order to really ensure fairness, I think the authors should have also considered memory/FLOP footprints and wall-clock time consumption.

The authors also reported preliminary experimental results on CIFAR-10 and ImageNet. Although the initial results look promising and I understand the lack of time, as an AC, I cannot make premature conclusions based on plots that only display 2 or 6 epochs of training. Also, some important training details are missing here—whether dense or block diagonal version is used, how many epochs correspond to the pretraining stage (as “switch-off” happens at 500 epochs in Fig 14), why can CIFAR-10 losses be such small after only 6 epochs, etc. Such missing details prevent me from accurately assessing the validity and significance of the preliminary results.

Overall, I think the paper has clear merits but the results are not ripe enough to be presented at ICLR. At this time, I recommend rejection.

Decision: Reject