Greedy Shallow Networks: An Approach for Constructing and Training Neural Networks

TL;DR

This paper introduces a greedy algorithm for constructing shallow neural networks with ReLU activation, using integral representations to efficiently select network parameters, potentially reducing reliance on backpropagation training.

Contribution

The paper proposes a novel greedy approach utilizing the ridgelet transform to efficiently build and initialize shallow neural networks, offering an alternative to traditional training methods.

Findings

The method effectively constructs neural networks with competitive performance.

It provides improved initializations that can sometimes replace full training.

Numerical experiments validate the approach's advantages over conventional techniques.

Abstract

We present a greedy-based approach to construct an efficient single hidden layer neural network with the ReLU activation that approximates a target function. In our approach we obtain a shallow network by utilizing a greedy algorithm with the prescribed dictionary provided by the available training data and a set of possible inner weights. To facilitate the greedy selection process we employ an integral representation of the network, based on the ridgelet transform, that significantly reduces the cardinality of the dictionary and hence promotes feasibility of the greedy selection. Our approach allows for the construction of efficient architectures which can be treated either as improved initializations to be used in place of random-based alternatives, or as fully-trained networks in certain cases, thus potentially nullifying the need for backpropagation training. Numerical experiments demonstrate the tenability of the proposed concept and its advantages compared to the conventional techniques for selecting architectures and initializations for neural networks.

Equations

$$g(x)=\sum _{n=1}^N{c}_n\sigma (a_{n}x+b_n),$$

$$\underset{c_1,\dots ,c_N\in \mathbb{R}}{\min }\Vert f(x_{tr})-\sum _{n=1}^N{c}_n\sigma (a_n{x}_{tr}+b_n)\Vert .$$

$$\mathcal{D}=\{\sigma ({\tilde{a}}_n{x}_{tr}+{\tilde{b}}_n){\}}_{n=1}^M,$$

$$g(x)=\sum _{n=1}^N{w}_n\sigma ({\bar{a}}_{n}x+{\bar{b}}_n)$$

$$\mathcal{D}=\{\sigma ({\bar{a}}_j{x}_{tr}+{\bar{b}}_j){\}}_{j\in \Lambda }\subset {\mathbb{R}}^{N_{tr}};$$

$$\underset{c_1,\dots ,c_N\in \mathbb{R}}{\min }\Vert f_{tr}-\sum _{n=1}^N{c}_n\sigma ({\bar{a}}_n{x}_{tr}+{\bar{b}}_n)\Vert ;$$

$$\text{Find elements}{g}_1,\dots ,g_N\in \mathcal{D}\text{and coefficients}{c}_1,\dots ,c_N\in \mathbb{R}$$

$$\text{such that}{\Vert f-\sum _{n=1}^N{c}_n{g}_n\Vert }_{\mathcal{X}}=\underset{\begin{array}{c}{g}_1,\dots ,g_N\in \mathcal{D}\\ c_1,\dots ,c_N\in \mathbb{R}\end{array}}{\min }{\Vert f-\sum _{n=1}^N{c}_n{g}_n\Vert }_{\mathcal{X}}.$$

$$g_{N}=\operatorname*{argmax}_{g\in\mathcal{D}}\big{|}\big{\langle}f_{N-1},\,g\big{\rangle}\big{|};$$

$$\text{find}{c}_1,\dots ,c_N\in \mathbb{R}\text{such that}\sum _{n=1}^N{c}_n{g}_n=\mathrm{proj}(f,\mathrm{span}\{g_1,\dots ,g_N\});$$

$$f_N=f-\mathrm{proj}(f,\mathrm{span}\{g_1,\dots ,g_N\})=f-\sum _{n=1}^N{c}_n{g}_n.$$

$$f(x)=\underset{{\mathbb{R}}^d\times \mathbb{R}}{\int }c(a,b)\sigma (a\cdot x+b)\text{d}a\text{d}b,$$

$${\mathcal{R}}_{\tau }f(a,b)=\underset{{\mathbb{R}}^d}{\int }f(x)\overline{\tau (a\cdot x+b)}\text{d}x,$$

$$(2\pi {)}^{d-1}\underset{-\infty }{\overset{\infty }{\int }}\frac{\hat{\sigma }(z)\overline{\hat{\tau }(z)}}{|z{|}^d}\text{d}z=1,$$

$$f(x)=\underset{{\mathbb{R}}^d\times \mathbb{R}}{\int }{\mathcal{R}}_{\tau }f(a,b)\sigma (a\cdot x+b)\text{d}a\text{d}b.$$

$$\tau (z)=-\frac{(\exp (-z^2/2){)}^{\prime \prime \prime \prime }}{2(2\pi {)}^{d-1/2}}=-\frac{z^4-6z^2+3}{2(2\pi {)}^{d-1/2}}{e}^{-z^2/2},$$

$$f(x)$$

$$=\underset{{\mathbb{S}}^d}{\int }{\mathcal{C}\mathcal{R}}_{\tau }f(\bar{a},\bar{b})\sigma (\bar{a}\cdot x+\bar{b})\text{d}(\bar{a},\bar{b}),$$

$${\mathcal{C}\mathcal{R}}_{\tau }f(\bar{a},\bar{b})=\underset{0}{\overset{\infty }{\int }}{\mathcal{R}}_{\tau }(r\bar{a},r\bar{b})r^{d+1}\text{d}r.$$

$$f(x)=\cos (2\pi x)\exp (x).$$

$$f(x)=\sin (2\pi x)\exp (-x^2)+\cos (17x)\exp (x).$$

$$f(x,y)=\sin (\pi x)\cos (\pi y)\exp (-(x^2+y^2)).$$

$$f(x,y)=\cos (5(x+y))\sin (3(x-y))\exp (-(x^2+y^2)).$$

$$f(\bar{x})=\sin (2\pi (x_1+x_2+x_3+x_4)).$$

Full text

We present a greedy-based approach to construct an efficient single hidden layer neural network with the ReLU activation that approximates a target function. In our approach we obtain a shallow network by utilizing a greedy algorithm with the prescribed dictionary provided by the available training data and a set of possible inner weights. To facilitate the greedy selection process we employ an integral representation of the network, based on the ridgelet transform, that significantly reduces the cardinality of the dictionary and hence promotes feasibility of the greedy selection. Our approach allows for the construction of efficient architectures which can be treated either as improved initializations to be used in place of random-based alternatives, or as fully-trained networks in certain cases, thus potentially nullifying the need for backpropagation training. Numerical experiments demonstrate the tenability of the proposed concept and its advantages compared to the conventional techniques for selecting architectures and initializations for neural networks. *

Keywords: Greedy algorithm, Integral representation, Shallow neural network, Architecture search, Ridgelet transform.

2000 Mathematics Subject Classification: 68T07, 41A46, 65D15.

2012 Computing Classification System:

• •

Mathematics of computing $\rightarrow$ Mathematical analysis $\rightarrow$ Functional analysis $\rightarrow$ Approximation

• •

Mathematics of computing $\rightarrow$ Mathematical analysis $\rightarrow$ Mathematical optimization $\rightarrow$ Continuous optimization $\rightarrow$ Nonconvex optimization

• •

Networks $\rightarrow$ Network architectures $\rightarrow$ Network design principles $\rightarrow$ Layering

11footnotetext: Computational and Applied Mathematics Group, Oak Ridge National Laboratory, Oak Ridge, TN, USA22footnotetext: Department of Mathematics, University of Tennessee at Knoxville, Knoxville, TN, USA

1 Introduction

This paper deals with the problem of designing and training a shallow neural network for a given set of training data. Such a problem comes up naturally in applications where neural networks are deployed and is the topic of active research, see, e.g., [34, 42, 43] and the references therein. For an overview of the area of neural network approximation we refer the reader to the books [18, 23].

Shallow neural networks are feed-forward neural networks with a single hidden layer. Such networks are known to be universal approximators, in the sense that any continuous target function can be approximated by a shallow network with any accuracy, see e.g. [14, 7]. However, finding an appropriate (or good) network architecture for approximating a given function is not always a straightforward task, and theoretical estimates rely on nonlinear approximation theory [10] and only exists in simple settings [27, 8, 33].

Generally one first has to select an appropriate architecture, which is commonly done according to a heuristic rather than the given data, see e.g. [17, 22]. Even though such an approach works, it often results in unnecessarily large networks where some nodes do not actually contribute to the approximation.

Next, a weight initialization scheme has to be specified. The conventional choices, e.g. [16, 21, 28] are based on sampling weights from an appropriate distribution to normalize the gradients during the backpropagation optimization. While such approaches are relevant for deep neural networks, they do not address the issues of shallow networks. As a result, randomly initialized neural networks often fail to achieve good approximation in practice, as explained in e.g., [35, 20].

Finally, the optimization algorithm and such hyperparameters as the learning and decay rates, batch size, etc. have to be selected, see, e.g., [9]. This choice is generally either based on a heuristic or is dictated by an extensive manual hyperparameter tuning, see e.g. [4].

Unfortunately, conventional approaches often lead the network to a bad local minimum cursed by the inherent non-convexity of the optimization problem, resulting in a low quality approximation, which is not stable and may converge slowly. In this paper we propose a greedy selection based method for constructing shallow network approximations that do not suffer the aforementioned drawbacks and offers a trade-off between the network complexity and approximation accuracy. We call networks obtained by this approach Greedy Shallow Networks (GSN).

The main contributions of our paper are the following:

• •

We present an alternative method for architecture selection and the initialization/training of shallow neural networks;

• •

We develop a novel approach for reducing the size of the dictionary by utilizing the ridgelet transform and the integral representation of neural networks;

• •

We derive the procedure for selecting the network architecture based on the available training data, rather than on an established heuristic.

1.1 Sample problem

We demonstrate our approach in the following simple setting. Let $f(x)$ be a one-dimensional real-valued target function we want to approximate with respect to a norm $\left\|\cdot\right\|$ by a neural network $g(x)$ with a single hidden layer of $N$ nodes, and an activation function $\sigma$, namely:

[TABLE]

where $a_{n},b_{n},c_{n}\in\mathbb{R}$. From the computational perspective it is beneficial to keep the number of nodes $N$ sufficiently small while still providing a good approximation to the target function $f(x)$. Hence one is required to select appropriate values for $N$ and inner and outer weights $\{(a_{n},b_{n},c_{n})\}_{n=1}^{N}$. This problem can first be simplified to only selecting the number of nodes $N$ and the corresponding inner weights $\{(a_{n},b_{n})\}_{n=1}^{N}$. Indeed, denote by $\{x_{tr},f_{tr}\}\in\mathbb{R}^{N_{tr}}$ the vectors of available training points and values respectively, then the optimal values of the outer weights $c_{1},\ldots,c_{N}$ can be obtained easily by solving the convex optimization problem

[TABLE]

In turn, we propose that the problem of selecting pairs $\{(a_{n},b_{n})\}_{n=1}^{N}$ with small $N$ can be attacked with an iterative greedy strategy, for instance, such as renown Matching Pursuit ([25]) or Orthogonal Matching Pursuit ([32]). Indeed, the classical greedy approximation task is as follows: for an element $f$, find a sparse linear approximation by the elements of a given set $\mathcal{D}$. The problem of constructing a network thus can be presented as follows: take sufficiently large number $M$ of possible values of the inner weights $\{(\tilde{a}_{n},\tilde{b}_{n})\}_{n=1}^{M}\subset\mathbb{R}^{2}$, and construct the dictionary

[TABLE]

which is then used to find an appropriate $m$-term combination $f(x)\approx\sum_{j=1}^{m}\gamma_{j}\,\sigma(\tilde{a}_{n_{j}}x+\tilde{b}_{n_{j}})$. However, in our context the problem of constructing such a linear combination is still hardly tangible as the domain of weights $a$ and $b$ is the whole real plane $\mathbb{R}^{2}$, which results in an unreasonably large dictionary $\mathcal{D}$. One way to make this approach tangible is to impose restrictions on the activation function. Specifically, in this paper we consider the ReLU activation $\sigma(z)=\max\{z,0\}$, which appears to be the standard choice of activation in most modern network architectures. Taking into account positive homogeneity of ReLU, we rewrite the network (1) by scaling the inner weights $\{(a_{n},b_{n})\}_{n=1}^{N}$ to be on the unit circle, i.e. $a,b\in\mathbb{R}\longrightarrow\bar{a},\bar{b}\in\mathbb{S}^{1}\subset\mathbb{R}^{2}$ such that

[TABLE]

with $w_{n}=c_{n}/\sqrt{a_{n}^{2}+b_{n}^{2}}$. The advantage of the above formulation is that each node is described by the single point $(\bar{a}_{n},\bar{b}_{n})$ on the circle $\mathbb{S}^{1}\subset\mathbb{R}^{2}$ rather than a point on the plane $\mathbb{R}^{2}$. Hence the dictionary (3) can now be obtained by only sampling the points from the circle rather than the whole plane. Next, we employ a greedy algorithm with the constructed dictionary $\mathcal{D}$ to find an $N$-term approximation to the target function $f$. Finally, once an appropriate approximation $f(x)\approx\sum_{n=1}^{N}\gamma_{n}\,\sigma(\bar{a}_{n}x+\bar{b}_{n})$ is constructed, the optimal output weights $c_{1},\ldots,c_{N}$ can be found via the convex minimization problem (2).

1.2 Proposed approach

In this section we summarize the proposed approach for a case of general multivariate target function $f:\mathbb{R}^{d}\to\mathbb{R}$. Denote by $f_{tr}\in\mathbb{R}^{N_{tr}}$ the vector of provided training values. Then the process of constructing the GSN consists of the following steps:

1. GSN 1.

Sample points $\{(\bar{a}_{j},\bar{b}_{j})\}_{j\in\Lambda}$ from the sphere $\mathbb{S}^{d}$ and construct the dictionary

[TABLE] 2. GSN 2.

Optionally, reduce the dictionary size, e.g. via the method proposed in Section 3.1; 3. GSN 3.

Use the Orthogonal Greedy Algorithm (defined in Section 2.1) to select the number of nodes $N$ and the corresponding inner weights $(\bar{a}_{n},\bar{b}_{n})$, $1\leq n\leq N$; 4. GSN 4.

Compute the outer weights $c_{n}$ by solving the convex optimization problem on the training data $\{x_{tr},f_{tr}\}$

[TABLE] 5. GSN 5.

Optionally, further train the network via a backpropagation algorithm using the constructed weights $\{(c_{n},\bar{a}_{n},\bar{b}_{n})\}_{n=1}^{N}$ as an initialization.

For convenience of the presentation in this paper we consider one-dimensional output functions $f:\mathbb{R}^{d}\to\mathbb{R}$; however our method is also applicable in a case of multi-dimensional output target function $f:\mathbb{R}^{d}\to\mathbb{R}^{t}$ with trivial modifications. Most notably, the convex optimization problem in the step GSN 4 would have to be solved for $c_{1},\ldots,c_{N}\in\mathbb{R}^{t}$.

1.3 Related work

Although greedy selection methods for neural network approximations has been studied in the past by various authors, e.g. [12, 2], they were primarily used as a theoretical tool rather than the practical approach for constructing neural networks.

Another approach that proposes the discretization of the inner weights space is the Random Vector Functional-Link (or RVFL) network, e.g. [30, 31], where the inner weights are generated randomly and the outer weights are found via a least square regularization, however such approach results in a large number of nodes and thus an unnecessarily complex networks.

The dictionary reduction step in our approach (defined in Section 3) relies on neural network integral representations, particularly on the ridgelet transform. Similar approach is considered in [36], where the authors suggest to use the values of the ridgelet transform to derive a probability distribution from which the inner weights are randomly selected. We replace the random selection step with a greedy selection step which, even though being more expensive as a selection method, results in a better approximation with a much fewer number of nodes. Moreover, the integral representation is completely optional in our approach, it is used only to decrease the dictionary size to accelerate the realization of the greedy algorithm.

Alternative methods for discretizing integral representations have been considered by various authors. In particular, in [3] the authors employ a greedy method to discretize the solution of a certain convex optimization problem, while the same problem is solved via the conditional gradient algorithm in [1]. A related Monte–Carlo discretization method for integral representations of radial basis function (RBF) neural networks is considered in [26].

2 Greedy selection

A problem of greedy approximation can be stated as follows: find the best $N$-term approximation for a given element $f$ of a Banach space $\mathcal{X}$ with respect to the dictionary $\mathcal{D}$ (a set of elements that is dense in $\mathcal{X}$), i.e.

[TABLE]

Greedy algorithms offer ways of approximately solving such a problem by iteratively expanding the support of the allowed linear combinations. Despite such an abstract setting, greedy algorithms are known to be successfully applied to many practical problems with various modifications, see e.g. [41, 2, 5, 13].

In this paper we are using the Orthogonal Greedy Algorithm (OGA), which is defined below. For the purpose of presentation we only state the algorithm in Euclidean space setting; however, like other greedy algorithms, it is applicable in a general setting of infinite-dimensional Banach spaces, thus allowing for an extension of our approach for other problem settings. The suitability of employing greedy algorithm for selecting the nodes is supported by theoretical research (including such results as convergence rate estimates and the provable near-best sparse approximation guarantees) in the approximation theory community, which is beyond the scope of this paper. For an in-depth overview of the field of greedy algorithms we refer the reader to the book [38].

2.1 Orthogonal Greedy Algorithm (OGA)

Let the set $\mathcal{D}$ be a dictionary in $\mathbb{R}^{N_{tr}}$, i.e. the elements of $\mathcal{D}$ are normalized and $\overline{\operatorname*{span}\mathcal{D}}=\mathbb{R}^{N_{tr}}$. Let $f\in\mathbb{R}^{N_{tr}}$ be an element to be approximated and let $\operatorname*{proj}(f,G)$ denote the orthogonal projection of $f$ onto the subspace $G\subset\mathbb{R}^{N_{tr}}$. Then the OGA constructs sparse approximation of $f$ with respect to $\mathcal{D}$ in the following way.

OGA: Set $f_{0}=f$ and on each iteration $N\geq 1$ perform the following steps:

Choose an element $g_{N}\in\mathcal{D}$ that maximizes the inner product with $f_{N-1}$

[TABLE] 2. 2.

Compute the orthogonal projection of $f$ onto $\operatorname*{span}\{g_{1},\ldots,g_{N}\}$

[TABLE] 3. 3.

Set the next remainder $f_{N}$ as

[TABLE]

This procedure can be continued until either the sufficient approximation accuracy (measured by $\left\|f_{N}\right\|$) is reached, or the maximal number of iterations is completed, hence allowing for the trade-off between the representation sparsity and the approximation accuracy.

We note that in a more general setting the performance of greedy algorithms suffers from the dictionary size, which is typically large in high dimensions. In particular, in the case of infinitely-dimensional spaces, step $1$ of the OGA might not even be feasible. Conventional ways to overcome this issue would be to consider a weak version of a greedy algorithm or to perform a restricted depth dictionary search, see e.g. [39, 40]. Additionally, the computations required in steps $1$ and $2$ can be relaxed by considering an approximate version of a greedy algorithm in which the computations are performed not precisely but with allowed inaccuracies of controllable magnitude, see e.g. [19, 15]. Such approaches are well studied and are known to simplify the realization of a greedy algorithm.

In the current manuscript we focus on the case of finitely-dimensional Hilbert spaces and consider the original formulation of the OGA with the exact calculations in steps $1$–$3$. The main challenge in our context is the large dimensionality of the dictionary $\mathcal{D}$, which we address by deriving a novel approach for domain clustering based on the ridgelet transform. Namely, in Section 3.1 we propose a method for reducing the dictionary size which is based on an integral representation of a neural network.

3 Integral representations

In this section we propose a method for the dictionary size reduction, based on an integral representation of the network, that restricts the sampling domain for the inner weights from the $d$-sphere $\mathbb{S}^{d}$ to its subset. Namely, for a given target function $f:\mathbb{R}^{d}\to\mathbb{R}$, we consider representations of the form

[TABLE]

called the neural network integral representation of $f$ with the kernel $c(a,b)$. Note that for a given function $f$ there can be infinitely many suitable kernels $c$ such that (4) holds. The idea of the dictionary reduction comes from the observation that the inner weights should be sampled only from the support of the kernel $c$ rather than the whole space.

To give a concrete example and to demonstrate the dictionary reduction property, we consider here the representation based on the ridgelet transform. In our numerical experiments we observe that the proposed technique allows for a dictionary reduction by approximately $50\%$ in case of one-dimensional target function $f$ (examples 4.1 and 4.2), and by $70\%$ in case of two-dimensional one (examples 4.3 and 4.4), without a noticeable affect on the quality of the constructed network.

We also note that there are more efficient ways to achieve search space reduction, such as those based on more fitting integral representations, see e.g. [29, 33]. However the overview of such methods is outside the scope of the current paper and will be discussed in future publications.

3.1 Dictionary reduction procedure

One of the integral representations (4) that we consider below originates from the harmonic analysis perspective to shallow neural networks, see, e.g. [6], and employs the ridgelet transform given by

[TABLE]

with some choice of function $\tau:\mathbb{R}\to\mathbb{C}$. It is known that under the admissibility condition

[TABLE]

where $\hat{\sigma},\hat{\tau}$ are the Fourier transforms of $\sigma$ and $\tau$, for a large class of target functions $f$ (see e.g. [37] for details) the following reconstruction formula holds

[TABLE]

In particular, the ReLU activation $\sigma(z)=\max\{z,0\}$ forms an admissible pair with the function

[TABLE]

which is the function that we use in our numerical experiments.

By switching to spherical coordinates and using the positive homogeneity of the ReLU, we can rewrite the reconstruction formula (6) in the following way

[TABLE]

where $\mathcal{CR}_{\tau}f$ is what we call the collapsed ridgelet transform, defined as

[TABLE]

Now the reduced dictionary $\mathcal{D}^{\prime}$ can by obtained by keeping only those vectors $\sigma(\bar{a}_{n}x_{tr}+\bar{b}_{n})$ in $\mathcal{D}$, whose weights $(\bar{a}_{n},\bar{b}_{n})$ satisfy $|\mathcal{CR}_{\tau}f(\bar{a}_{n},\bar{b}_{n})|>\epsilon$ with a problem dependent threshold $\epsilon>0$. The value of $\epsilon$ in this context controls the trade-off between the potential approximation accuracy and the computational complexity. Namely, larger values of $\epsilon$ promote the reduced cardinality of the dictionary $\mathcal{D}^{\prime}$, and thus simplify the greedy selection step, however at a possible cost of approximation accuracy.

The rationale of why thresholding the collapsed ridgelet transform results in efficient dictionary size reduction is supported by our theoretical findings which state that $\mathcal{CR}_{\tau}f$ is well-localized. For the purpose of this paper, however, we only support this claim with numerical experiments (Figures 3(b), 5(b), 7(a), and 9(a)).

A particular importance of the integral representation (7) is that, under the assumption that $\mathcal{CR}_{\tau}f\in L_{1}(\mathbb{S}^{d})$, it can be deduced that the vector of training values $f_{tr}$ is in the closure of the convex hull of the dictionary $\mathcal{D}=\{\sigma(\bar{a}\cdot x_{tr}+\bar{b})\}_{(\bar{a},\bar{b})\in\mathbb{S}^{d}}$. As follows from the greedy approximation theory [11], in that case the GSN initialization described in Subsection 1.2 approximates the training data $f_{tr}$ with the rate of $\mathcal{O}(N^{-1/2})$, where $N$ is the number of nodes (i.e. iterations of the greedy algorithm). Note that this is the worst-case estimate, in practice we generally observe exponential convergence rates.

4 Numerical experiments

All numerical experiments are implemented in Python 3.6. Neural networks are trained in Keras with the TensorFlow backend, and the rest of the computations are performed with the use of Numpy package. For reproducibility, the Numpy random seed is set to [math] for all presented numerical experiments. The code used to obtain the presented numerical results is publicly available at https://github.com/sukiboo/sgn.

In order to construct the dictionary $\mathcal{D}$ for a $d$-dimensional target function, we sample $10,000\,d$ points from the $d$-sphere $\mathbb{S}^{d}$. In case $d=1$ we sample the points uniformly from the interval $[-\pi,\pi)$. In case $d=2$ we employ the golden spiral method (shown in Figure 1) to sample the points from the sphere $\mathbb{S}^{2}$. For other values of $d$ we generate and normalize independent Gaussian vectors, which are uniformly distributed (in expectation) over $\mathbb{S}^{d}$.

For the settings $d=1,2$ we reduce the size of the dictionary $\mathcal{D}$ via the method described in Section 3. Namely, we compute the ridgelet transform $\mathcal{R}f(r,\phi)$ by discretizing the integral (5) on the provided training data, and the collapsed ridgelet transform $\mathcal{CR}f(\phi)$ by summing the ridgelet transform along the sampled directions $\phi\in\mathbb{S}^{d}$ using the formula (8). The values of $\mathcal{CR}f(\phi)$ are used to reduce the dictionary $\mathcal{D}$ by discarding the sampled points $\phi$ that correspond to the values of $|\mathcal{CR}f(\phi)|$ that are smaller than the $0.1\%$ of the maximal value of $|\mathcal{CR}f(\phi)|$.

The Orthogonal Greedy Algorithm is applied to the vector of training values with respect to the dictionary $\mathcal{D}$. The number of constructed nodes is chosen as the iteration of the greedy algorithm that minimizes the approximation error on the given validation set. That is, the resulting network architecture is the one that minimizes the validation error as shown in Figures 3(c), 5(c), 7(b), 9(b), and 10(a) for the respective examples.

For each target function we construct two network approximations and compare them. First, we obtain the GSN initialization as described in Subsection 1.2, and train it via the backpropagation. Next, we train the neural network with random initialization with ‘TruncatedNormal’ initial weights for both the hidden and the output layers. Both networks are trained with the ‘Adam optimizer’, see [24], with the initial learning rate $10^{-3}$ and the exponential decay rate of $4.6\times 10^{-4}$. The training is performed for $10,000$ epochs to ensure that the learning is completed; plots of the training losses are presented to demonstrate that every training process is completed and an appropriate hyperparameters are selected. Since it is known that networks with random initialization and small number of parameters are sensitive to the initial state and are prone to get stuck in a bad local minimum, in some examples we use the smaller batch sizes for randomly initialized networks than for the GSN ones; otherwise the choice of hyperparameters for the training of two networks is the same. Each randomly initialized network we train for $10$ times and report the best result the network achieved; networks with the GSN initialization are trained once since the initialization is deterministic.

4.1 Example 1

Consider the target function $f:[-1,1]\to\mathbb{R}$ given by

[TABLE]

We take $50$ training points, $15$ validation points, and $1,000$ test points. The batch size is set to be $50$ for the network with the GSN initialization and $1$ for the network with random initialization. The resulting approximations are presented in Figure 2.

The ridgelet and collapsed ridgelet transforms are shown in Figures 3(a) and 3(b) respectively. In this example the number of constructed nodes is $23$ (Figure 3(c)).

4.2 Example 2

Consider the target function $f:[-1,1]\to\mathbb{R}$ given by

[TABLE]

We take $100$ training points, $20$ validation points, and $1,000$ test points. The batch size is set to be $100$ for the network with the GSN initialization and $1$ for the network with random initialization. The resulting approximations are presented in Figure 4.

The ridgelet and collapsed ridgelet transforms are shown in Figures 5(a) and 5(b) respectively. In this example the number of constructed nodes is $40$ (Figure 5(c)).

Note that in this example networks with random initialization fails to achieve a reasonable approximation despite multiple attempts. This is likely caused by restricting the number of nodes and absence of additional regularization as explained in [20]. The network with the GSN initialization, on the other hand, under the same conditions recovers every oscillation of the target function.

4.3 Example 3

Consider the target function $f:[-1,1]^{2}\to\mathbb{R}$ given by

[TABLE]

We take $256$ training points, $50$ validation points, and $10,000$ test points. The batch size is set to be $256$ for the network with the GSN initialization and $3$ for the network with random initialization. The resulting approximations are presented in Figure 6.

The collapsed ridgelet transform is shown in Figure 7(a), the number of constructed nodes is $50$ (see Figure 7(b)).

4.4 Example 4

Consider the target function $f:[-1,1]^{2}\to\mathbb{R}$ given by

[TABLE]

We take $1,024$ training points, $200$ validation points, and $10,000$ test points. The batch size is set to be $1,024$ for the network with the GSN initialization and $11$ for the network with random initialization. The resulting approximations are presented in Figure 8.

The collapsed ridgelet transform is shown in Figure 9(a), the number of constructed nodes is $84$ (see Figure 9(b)).

4.5 Example 5

Consider the target function $f:[-1,1]^{4}\to\mathbb{R}$ given by

[TABLE]

We take $4,000$ training points, $400$ validation points, and $10,000$ test points. The number of constructed nodes is $109$ (see Figure 10(a)).

The approximation errors on the test set are the following: 3.20e-01 (GSN initialization), 6.99e-02 (network trained with the GSN initialization), 4.29e-01 (network trained with random initialization). Note that in this example the GSN initialization by itself provides lower accuracy than the network trained with random initialization, which is caused by a higher dimensionality of the training data and the iterative nature of the greedy algorithm. The network trained with the GSN initialization, however, still achieves significantly better accuracy when compared to a random initialization.

5 Discussion

In this paper we propose a novel way of initializing shallow neural networks with the ReLU activation via an iterative greedy strategy. GSN initialization provides an interpretable way of selecting an architecture by allowing for a measurable trade-off between the number of nodes and the approximation accuracy. Due to the nature of greedy selection, the GSN is computationally efficient in the sense that each constructed node is essential and contributes to the approximation of a target function.

Presented numerical results show that initializing a network with our approach stabilizes and improves the training of a shallow network. In particular, greedy initialization allows us to use the full training set as one batch and run the optimization algorithm for a much fewer number of epochs, which we did not find possible with the randomly initialized networks. Moreover, since, for a chosen set of parameters, the GSN initialization is deterministic, our approach is not affected by the issue of the initialization sensitivity. Additionally, we observe that in some cases the GSN initialization can serve as a fully-trained neural network, thus potentially removing the need of learning via the backpropagation.

We also note that the approximation accuracy in the presented numerical examples is by no means the best our method can achieve; the purpose of the numerical section is to compare the GSN initialization with the classical approach in terms of training stability and hyperparameter dependency.

While in this paper we only consider the ReLU activation, our method can be employed with any other positively homogeneous activation function. For a wider range of activations, the appropriate modifications are required, which is one of the topics of our future research projects. Similarly, the currently employed $\ell_{2}$-norm can be replaced with any other norm, which is in immediate interest for various applications. Finally, we mention that even though the current method cannot be directly extended for the case of deep neural networks, the work in this direction is our main priority.

Acknowledgment

This material is based upon work supported in part by: the U.S. Department of Energy, Office of Science, Early Career Research Program under award number ERKJ314; U.S. Department of Energy, Office of Advanced Scientific Computing Research under award numbers ERKJ331 and ERKJ345; the National Science Foundation, Division of Mathematical Sciences, Computational Mathematics program under contract number DMS1620280; and by the Laboratory Directed Research and Development program at the Oak Ridge National Laboratory, which is operated by UT-Battelle, LLC., for the U.S. Department of Energy under contract DE-AC05-00OR22725.