Sampling weights of deep neural networks
We introduce a probability distribution, combined with an efficient sampling
algorithm, for weights and biases of fully-connected neural networks. In a
supervised learning context, no iterative optimization or gradient computations
of internal network parameters are needed to obtain a trained network. The
sampling is based on the idea of random feature models. However, instead of a
data-agnostic distribution, e.g., a normal distribution, we use both the input
and the output training data to sample shallow and deep networks. We prove that
sampled networks are universal approximators. For Barron functions, we show
that the $L^2$-approximation error of sampled shallow networks decreases with
the square root of the number of neurons. Our sampling scheme is invariant to
rigid body transformations and scaling of the input data, which implies many
popular pre-processing techniques are not required. In numerical experiments,
we demonstrate that sampled networks achieve accuracy comparable to iteratively
trained ones, but can be constructed orders of magnitude faster. Our test cases
involve a classification benchmark from OpenML, sampling of neural operators to
represent maps in function spaces, and transfer learning using well-known
architectures.
Erik Lien Bolager, Iryna Burak, Chinmay Datar, Qing Sun, Felix Dietrich
# Sampling weights of deep neural networks

###### Abstract

We introduce a probability distribution, combined with an efficient sampling algorithm, for weights and biases of fully-connected neural networks. In a supervised learning context, no iterative optimization or gradient computations of internal network parameters are needed to obtain a trained network. The sampling is based on the idea of random feature models. However, instead of a data-agnostic distribution, e.g., a normal distribution, we use both the input and the output training data of the supervised learning problem to sample both shallow and deep networks. We prove that the sampled networks we construct are universal approximators. We also show that our sampling scheme is invariant to rigid body transformations and scaling of the input data. This implies many popular pre-processing techniques are no longer required. For Barron functions, we show that the \(L^{2}\)-approximation error of sampled shallow networks decreases with the square root of the number of neurons. In numerical experiments, we demonstrate that sampled networks achieve comparable accuracy as iteratively trained ones, but can be constructed orders of magnitude faster. Our test cases involve a classification benchmark from OpenML, sampling of neural operators to represent maps in function spaces, and transfer learning using well-known architectures.

## 1 Introduction

Training deep neural networks involves finding all weights and biases. Typically, iterative, gradient-based methods are employed to solve this high-dimensional optimization problem. Randomly sampling all weights and biases before the last, linear layer circumvents this optimization. However, the drawback of this approach is that the high-dimensional probability distribution must be chosen. Random projection networks [46] or extreme learning machines [25] involve choosing the distribution completely problem- and data-agnostic, e.g., a normal distribution. In this work, we introduce a data-driven sampling scheme to construct weights and biases close to gradients of the target function. This idea provides a solution to three main challenges that have prevented randomly sampled networks to compete successfully against iterative optimization methods in the setting of supervised learning.

Figure 1: Random feature models choose weights in a data-agnostic way, compared to sampling them where it matters: at large gradients. The arrows illustrate where the network weights are placed.
**Deep neural networks** Random feature models and extreme learning machines are typically defined for networks with a single hidden layer. Our sampling scheme accounts for the high-dimensional ambient space that is introduced after this layer, and thus deep networks can be constructed efficiently.

**Approximation accuracy** Gradient-based, iterative approximation can find accurate solutions with a relatively small number of neurons. Randomly sampling weights using a data-agnostic distribution often requires thousands of neurons to compete. Our sampling scheme takes into account the given data points and function values, leading to accurate and width-efficient approximations.

**Interpretability** Sampling weights and biases completely randomly, i.e., without taking into account the given data, leads to networks that do not incorporate any information about the given problem. We employ and analyze a recently introduced weight construction technique [22] that uses the direction between pairs of data points to construct individual weights and biases. In addition, we propose a sampling distribution over these data pairs that leads to efficient use of weights; cf. Figure 1. The distribution also leads to invariance to orthogonal transformations and scaling of the input data, which makes many common pre-processing techniques redundant.

## 2 Related work

Sampled weights are closely related to random Fourier features, Li et al. [36] and Liu et al. [38] review and unify its theory and algorithms. Shallow [12] and deep networks [24] with normally distributed random weights have been analyzed regarding classification accuracy. Similar to finite element methods in scientific computing, earlier work placed activation function at equally spaced collocation points [42; 8]. Gallicchio and Scardapane [23] provide a review of deep random feature models, and discuss autoencoders and reservoir computing (resp. echo-state networks [29]) in this context. The latter are randomly sampled, recurrent networks to model dynamical systems [5]. Barron spaces are very useful to study approximation errors and convergence rates of neural networks [1; 16], as well as regularization techniques, esp. Tikohnov and Total Variation [35]. A lot of work [46; 15; 13; 48; 54] surrounds the approximation rate of \(\mathcal{O}(m^{-1/2})\) for neural networks with one hidden layer, originally proved by Barron [1]. The rate with respect to width \(m\), but not the constant, is independent of the input space dimension. This implies that neural networks can mitigate the curse of dimensionality, as opposed to many approximation methods with fixed, non-trained basis functions [11], including random feature models with data-agnostic probability distributions. Regarding networks with more than one hidden layer, E and Wojtowytsch [13; 14] discuss simple examples that are not Barron functions, i.e., cannot be represented by shallow networks. However, they can be represented by deep networks. The convergence rates of over-parameterized networks with one hidden layer is considered in [15], with a comparison to the Monte-Carlo approximation. Fiedler et al. [19] and Fornasier et al. [20] prove that for given, comparatively small networks with one hidden layer, all weights and biases can be recovered exactly by evaluating the network at specific points in the input space. The work of Spek et al. [50] showed a certain duality between weight spaces and data spaces, albeit in a purely theoretical setting. Recent work from Bollt [6] analyzes individual weights in networks by visualizing the placement of ReLU activation functions in space. Weights are also sampled in Bayesian neural networks [43; 4; 21; 51]. The goal is to learn a good posterior distribution and ultimately express uncertainty around both weights and the output of the network. These methods are computationally often on par with or worse than iterative optimization. Most relevant to our work is the weight construction method by Galaris et al. [22], who proposed to use pairs of data points to construct weights. Their primary goal was to randomly sample weights that nevertheless capture a potentially low-dimensional structure in the data. No further analysis was provided, and only a uniform distribution over the data pairs was used.

## 3 Mathematical framework

We introduce sampled networks, which are neural networks where for all hidden layers, each pair of weight and bias is completely determined by two points from the input space. This duality between weights and data has been shown theoretically [50], here, we provide an explicit relation. The weights are informed by the difference between the two points, and the bias is the inner product between the weight and one of the two points. After the weights and biases of the hidden layers are constructed, we must only solve an optimization problem for the coefficients of a linear layer, mapping the output
from the last hidden layer to the final output. To formalize the construction of sampled networks, we start by introducing some notation.

Let \(\mathcal{X}\subseteq\mathbb{R}^{D}\) be the input space, and \(\Phi\) a neural network with \(L\) hidden layers, parameters \(\{W_{l},b_{l}\}_{l=1}^{L+1}\), and activation function \(\phi:\mathbb{R}\rightarrow\mathbb{R}\). We write \(\Phi^{(l)}(x)=\phi(W_{l}\Phi^{(l-1)}(x)-b_{l})\) as the output of the \(l\)th layer, with \(\Phi^{(0)}(x)=x\) for all \(x\in\mathcal{X}\). The two activation functions we focus on in this work are the ReLU (rectified linear unit), \(\phi(x)=\max\{x,0\}\), and the tanh (hyperbolic tangent). We set \(N_{l}\) to be the number of neurons in the \(l\)th layer, with \(N_{0}=D\) and \(N_{L+1}\) as the output dimension. We write \(w_{l,i}\) for the \(i\)th row of \(W_{l}\) and \(b_{l,i}\) for the \(i\)th entry of \(b_{l}\). Building upon work of Galaris et al. [22], we now introduce sampled networks. We start with an arbitrary probability distribution over pairs of data points, and later refine it to include known function values.

**Definition 1**.: _Let \(\Phi\) be a neural network with \(L\) hidden layers. For \(l=1,\ldots,L\), let \((x_{0,i}^{(1)},x_{0,i}^{(2)})_{i=1}^{N_{l}}\) be pairs of points sampled over \(\mathcal{X}\times\mathcal{X}\). We say \(\Phi\) is a sampled network if the weights and biases of every layer \(l=1,2,\ldots,L\) and neurons \(i=1,2,\ldots,N_{l}\), are of the form_

\[w_{l,i}=s_{1}\frac{x_{l-1,i}^{(2)}-x_{l-1,i}^{(1)}}{\|x_{l-1,i}^{(2)}-x_{l-1, i}^{(1)}\|^{2}},\quad b_{l,i}=\langle w_{l,i},x_{l-1,i}^{(1)}\rangle+s_{2} \tag{1}\]

_where \(s_{1},s_{2}\in\mathbb{R}\) are constants, and \(x_{l-1,i}^{(j)}=\Phi^{(l-1)}(x_{0,i}^{(j)})\) for \(j=1,2\) and \(x_{l-1,i}^{(1)}\neq x_{l-1,i}^{(2)}\). The last set of weights and biases are \(W_{L+1},b_{L+1}=\arg\min\mathcal{L}(W_{L+1}\Phi^{(L)}(\cdot)-b_{L+1})\), where \(\mathcal{L}\) is a suitable loss._

The scalars are used to fix what values the points \(x^{(1)},x^{(2)}\) take on when the activation function is applied; cf. Figure 2. For ReLU, we set \(s_{1}=1\) and \(s_{2}=0\), to let \(x^{(1)}\) be mapped to zero and \(x^{(2)}\) to one. For tanh, we set \(s_{1}=2s_{2}\) and \(s_{2}=\ln(3)/2\), which implies \(x^{(1)}\) and \(x^{(2)}\) are mapped to \(\mp 1/2\) respectively, and the midpoint between the two points is assigned to zero. We choose the constants such that in a regression problem with ReLU, we linearly interpolate values between the two points. For classification, the tanh construction introduces a boundary if \(x^{(1)}\) belongs to a different class than \(x^{(2)}\). We will use this idea later to define a useful distribution over pairs of points (cf. Definition 2).

The space of functions that sampled networks can approximate is not immediately clear. First, we are only using points in the input space to construct both the weights and the biases, instead of letting them take on any value. Second, there is a dependence between the bias and the direction of the weight. Third, for deep networks, the sampling space changes after each layer. These apparent restrictions require investigation into which functions we can approximate. We assume that the input space in Theorem 1 and Theorem 2 extends into its ambient space \(\mathbb{R}^{D}\) as follows. Let \(\mathcal{X}^{\prime}\) be any compact space of \(\mathbb{R}^{D}\) with reach \(\tau>0\). We then set the input space \(\mathcal{X}\) to be the space of points including \(\mathcal{X}^{\prime}\) and those that are at most \(\epsilon_{I}\) away from \(X^{\prime}\), given the canonical distance of \(\mathbb{R}^{D}\), where \(0<\epsilon_{I}<\tau\). In Theorem 1, we also consider \(L\) layers to show that the construction of weights in deep networks does not destroy the space of functions that networks with one hidden layer can approximate, even though we alter the space of weights we can construct when \(L>1\).

**Theorem 1**.: _For any number of layers \(L\in\mathbb{N}_{>0}\), the space of sampled networks with \(L\) hidden layers, \(\Phi\colon\mathcal{X}\rightarrow\mathbb{R}^{N_{L+1}}\), with activation function ReLU is dense in \(C(\mathcal{X},\mathbb{R}^{N_{L+1}})\)._

_Sketch of proof:_ We show that the possible biases that sampled networks can produce are all we need inside a neuron, and the rest can be added in the last linear part and with additional neurons. We then show that we can approximate any neural network with one hidden layer with at most \(6\cdot N_{1}\) neurons

Figure 2: Placement of the point pairs \(x^{(1)},x^{(2)}\) for activation functions ReLU (left) and tanh (right).
-- which is not much, considering the cost of sampling versus backpropagation. We then show that we can construct weights so that we preserve the information of \(\mathcal{X}\) through the first \(L-1\) layers, and then we use the arbitrary width result applied to the last hidden layer. The full proof can be found in Appendix A.1. We also proved similar results for the case of tanh activation function with one hidden layer. The proof differs fundamentally because tanh is not positive homogeneous.

We now show existence of sampled networks for which the \(L^{2}\) approximation error of Barron functions is bounded (cf. Theorem 2). We later demonstrate that we can actually construct such networks (cf. Section 4.1). The Barron space [1, 16] is defined as

\[\mathcal{B}=\{f\colon f(x)=\int_{\Omega}w_{2}\phi(\langle w_{1},x\rangle-b)d \mu(b,w_{1},w_{2})\text{ and }\|f\|_{\mathcal{B}}<\infty\}\]

with \(\phi\) being the ReLU function, \(\Omega=\mathbb{R}\times\mathbb{R}^{D}\times\mathbb{R}\), and \(\mu\) being a probability measure over it. The Barron norm is defined as \(\|f\|_{\mathcal{B}}=\inf_{\mu}\max_{(b,w_{1},w_{2})\in\text{supp}(\mu)}\{|w_{ 2}|(\|w_{1}\|_{1}+|b|)\}\).

**Theorem 2**.: _Let \(f\in\mathcal{B}\) and \(\mathcal{X}=[0,1]^{D}\). For any \(N_{1}\in\mathbb{N}_{>0}\), \(\epsilon>0\), and an arbitrary probability measure \(\pi\), there exist sampled networks \(\Phi\) with one hidden layer, \(N_{1}\) neurons, and ReLU activation function, such that_

\[\|f-\Phi\|_{2}^{2}=\int_{\mathcal{X}}\lvert f(x)-\Phi(x)\rvert^{2}d\pi(x)< \frac{(3+\epsilon)\|f\|_{\mathcal{B}}^{2}}{N_{1}}.\]

_Sketch of proof:_ It quickly follows from the results of E et al. [16], which showed it for regular neural networks, and Theorem 1. The full proof can be found in Appendix A.2.

Up until now, we have been concerned with the space of sampling networks, but not with the distribution of the parameters. We found that the following distribution works well, putting emphasis on points that are close and differ a lot with respect to the output of the true function. The norms used to define the densities are arbitrary over their given space, denoted by the subscript.

**Definition 2**.: _Let \(\mathcal{X}\) be the input space, and let \(\mathcal{Y}=f(\mathcal{X})\) be given function values. Let \(P\) be a continuous probability distribution defined through its conditional distribution with density \(p\), for \(l=1,2,\ldots,L\) and \(x_{0}^{(1)},x_{0}^{(2)}\in\mathcal{X}\), setting \(x_{l-1}^{(1)},x_{l-1}^{(2)}=\Phi^{(l-1)}(x_{0}^{(1)}),\Phi^{(l-1)}(x_{0}^{(2)})\),_

\[p\left(x_{0}^{(1)},x_{0}^{(2)}\mid\{W_{j},b_{j}\}_{j=1}^{l-1}\right)\propto \begin{cases}\frac{\|f(x_{0}^{(2)})-f(x_{0}^{(1)})\|_{\mathcal{Y}}}{\|x_{l-1}^ {(2)}-\ x_{l-1}^{(1)}\|_{\mathcal{X}_{l-1}}},&x_{l-1}^{(1)}\neq x_{l-1}^{(2)} \\ 0,&\text{otherwise},\end{cases} \tag{2}\]

_where \(\Phi^{(l-1)}(\cdot)\) is the sampled network up to layer \(l-1\) with parameters \(\{W_{j},b_{j}\}_{j=1}^{l-1}\), and \(\mathcal{X}_{l-1}=\Phi^{(l-1)}(\mathcal{X})\)._

The conditional distribution \(P\) is defined over pairs of points, but as the network parameters are uniquely identified by these points through Definition 1, it is also a distribution over these parameters. Using this distribution also comes with the benefit that sampling and training are not affected by rigid body transformations (affine transformation with orthogonal matrix) and scaling, as long as the transformation is invariant with respect to the true function. If \(H\) is such a transformation, we say it is equivariant with respect to \(f\), if there exists a scalar and rigid body transformation \(H^{\prime}\) such that \(H^{\prime}f(x)=f(H(x))\) for all \(x\in\mathcal{X}\), and invariant if \(H^{\prime}\) is the identity function. We also assume that norms \(\|\cdot\|_{\mathcal{Y}}\) and \(\|\cdot\|_{\mathcal{X}_{0}}\) in Equation (2) are such that orthogonal matrix multiplication is an isometry.

**Theorem 3**.: _Let \(f\) be the target function and \(H\) be a scalar and rigid body transformation, equivariant with respect to \(f\). If we have two sampled networks, \(\Phi,\hat{\Phi}\), with the same number of hidden layers \(L\) and neurons \(N_{1},\ldots,N_{L}\), where \(\Phi\colon\mathcal{X}\to\mathbb{R}^{N_{L+1}}\) and \(\hat{\Phi}\colon H(\mathcal{X})\to\mathbb{R}^{N_{L+1}}\), then the following statements hold:_

1. _If_ \(\hat{x}_{0,i}^{(1)}=H(x_{0,i}^{(1)})\) _and_ \(\hat{x}_{0,i}^{(2)}=H(x_{0,i}^{(2)})\)_, for all_ \(i=1,2,\ldots,N_{1}\)_, then_ \(\Phi^{(1)}(x)=\hat{\Phi}^{(1)}(H(x))\)_, for all_ \(x\in\mathcal{X}\)_._
2. _If_ \(H\) _is invariant w.r.t._ \(f\)_, then for any parameters of_ \(\Phi\)_, there exists parameters of_ \(\hat{\Phi}\) _such that_ \(\Phi(x)=\hat{\Phi}(H(x))\) _for all_ \(x\in\mathcal{X}\)_, and vice versa._
_._
3. _The probability measure_ \(P\) _over the parameters is invariant under_ \(H\)_._

_Sketch of proof:_ Any neuron in the sampled network can be written as \(\phi(\langle s_{1}\frac{w}{\|w\|^{2}},x-x^{(1)}\rangle-s_{2})\). As we divide by the square of the norm of \(w\), the scalar in \(H\) cancels. There is a difference between two points in both inputs of \(\langle\cdot,\cdot\rangle\), which means the translation cancels. Orthogonal matrices cancel due to isometry. When \(H\) is invariant with respect to \(f\), the loss function is also unchanged and lead to the same output. Similar argument is made for \(P\), and the theorem follows. The complete proof can be found in Appendix A.3.

If the input space is embedded in a larger ambient space \(\mathbb{R}^{D^{\prime}}\), with \(D<D^{\prime}\), we sample from a subspace with dimension \(\tilde{D}=\dim(\overline{\operatorname{span}\{\mathcal{X}\}})\leq D^{\prime}\). All the results presented in this section still hold due to orthogonal decomposition. However, the standard approach of backpropagation and initialization allows the weights to take on any value in \(\mathbb{R}^{D^{\prime}}\), which implies a lot of redundancy when \(\tilde{D}\ll D^{\prime}\). The biases are also more relevant to the input space than when initialized to zero -- potentially avoiding the issues highlighted by Holzmuller and Steinwart [27]. For these reasons, we have named the proposed method Sampling Where It Matters (SWIM), which is summarized in Algorithm 1. For computational reasons, we consider a random subset of all possible pairs of training points when sampling weights and biases. The subset size can be changed depending on the size of the training set and the number of neurons.

```
Constant:\(\boldsymbol{\epsilon}\in\mathbb{R}_{>0}\), \(\varsigma\in\mathbb{N}_{>0}\), \(L\in\mathbb{N}_{>0}\), \(\{N_{l}\in\mathbb{N}_{>0}\}_{l=1}^{L+1}\), and \(s_{1},s_{2}\in\mathbb{R}\) Data:\(X=\{x_{i}\colon x_{i}\in\mathbb{R}^{D},i=1,2,\ldots,M\}\), \(Y=\{y_{i}\colon f(x_{i})=y_{i}\in\mathbb{R}^{N_{L+1}},i=1,2,\ldots,M\}\) \(\Phi^{(0)}(x)=x\); for\(l=1,2,\ldots,L\)do \(\tilde{M}\leftarrow\varsigma\cdot\lceil\frac{N}{M}\rceil\cdot M\) ; \(P^{(l)}\in\mathbb{R}^{M}\); \(\tilde{X}=\{(x_{i}^{(1)},x_{i}^{(2)})\colon\Phi^{(l-1)}(x_{i}^{(1)})\neq\Phi^{ (l-1)}(x_{i}^{(2)})\}_{i=1}^{\tilde{M}}\sim\text{Uniform}(X\times X)\); for\(i=1,2,\ldots,\tilde{M}\)do \(\tilde{x}_{i}^{(1)},\tilde{x}_{i}^{(2)}\leftarrow\Phi^{(l-1)}(x_{i}^{(1)}), \Phi^{(l-1)}(x_{i}^{(2)})\); \(\tilde{y}_{i}^{(1)},\tilde{y}_{i}^{(2)}=f(x_{i}^{(1)}),f(x_{i}^{(2)})\); \(P_{i}^{(l)}\leftarrow\frac{\|\tilde{y}_{i}^{(2)}-\tilde{y}_{i}^{(1)}\|_{ \mathcal{Y}}}{\max\{\|\tilde{x}_{i}^{(2)}-\tilde{x}_{i}^{(1)}\|_{\mathcal{X}_ {i-1}}\epsilon\}}\);  end for\(W_{l}\in\mathbb{R}^{N_{l},N_{l-1}},b_{l}\in\mathbb{R}^{N_{l}}\); for\(k=1,2,\ldots,N_{l}\)do  Sample \((x^{(1)},x^{(2)})\) from \(\tilde{X}\), with replacement and with probability proportional to \(P^{(l)}\); \(\tilde{x}^{(1)},\tilde{x}^{(2)}\leftarrow\Phi^{(l-1)}(x^{(1)}),\Phi^{(l-1)}(x ^{(2)})\); \(W_{l}^{(k;\cdot)}\gets s_{1}\frac{\tilde{x}^{(2)}-\tilde{x}^{(1)}}{\| \tilde{x}^{(2)}-\tilde{x}^{(1)}\|}\hat{x}_{l}^{(1)}\hat{x}_{l}^{(k)}\gets \langle W_{l}^{(k;\cdot)},\tilde{x}^{(1)}\rangle+s_{2}\);  end for\(\Phi^{(l)}(\cdot)\leftarrow\phi(W_{l}\,\Phi^{(l-1)}(\cdot)-b_{l})\);  end for\(W_{L+1},b_{L+1}\leftarrow\operatorname*{arg\,min}_{\mathcal{L}}(\Phi^{(L)} (X),Y)\); return\(\{W_{l},b_{l}\}_{l=1}^{L+1}\)
```

**Algorithm 1**The SWIM algorithm, for activation function \(\phi\), and norms on input, output of the hidden layers, and output space, \(\|\cdot\|_{\mathcal{X}_{0}}\),\(\|\cdot\|_{\mathcal{X}_{l}}\), and \(\|\cdot\|_{\mathcal{Y}}\) respectively. Also, \(\mathcal{L}\) is a loss function, which in our case is always \(L^{2}\) loss, and \(\operatorname*{arg\,min}_{\mathcal{L}}\mathcal{L}(\cdot,\cdot)\) becomes a linear optimization problem.

We end this section with a time and memory complexity analysis of Algorithm 1. In Table 1, we list runtime and memory usage for three increasingly strict assumptions. The main assumption is that the dimension of the output is less than or equal to the largest dimension of the hidden layers. This is true for the problems we consider, and the difference in runtime without this assumption is only reflected in the linear optimization part. The term \(\lceil N/M\rceil\) is required because only a subset of points are considered when sampling. For the full analysis, see Appendix F.
## 4 Numerical experiments

We now demonstrate the performance of Algorithm 1 on numerical examples. Our implementation is based on the numpy and scipy Python libraries, and we run all experiments on a machine with 32GB system RAM (256GB in Section 4.3 and Section 4.4) and a GeForce 4x RTX 3080 Turbo GPU with 10GB RAM. The Appendix contains detailed information on all experiments. In Section 4.1 we compare sampling to random Fourier feature models, on the approximation of a Barron function. In Section 4.2 we compare classification accuracy of sampled networks to iterative, gradient-based optimization in a classification benchmark with real datasets. In Section 4.3 we demonstrate that more specialized architectures can be sampled, by constructing deep neural architectures as PDE solution operators. In Section 4.4 we show how to use sampling of fully-connected layers for transfer learning. For the probability distribution over the pairs in Algorithm 1, we always choose the \(L^{\infty}\) norm for \(\|\cdot\|_{\mathcal{Y}}\) and for \(l=1,2,\ldots,L\), we choose the \(L^{2}\) norm for \(\|\cdot\|_{\mathcal{X}_{l-1}}\).

The code to reproduce the experiments from the paper can be found at

[https://gitlab.com/felix.dietrich/swimnetworks-paper](https://gitlab.com/felix.dietrich/swimnetworks-paper)

An up-to-date code base is maintained at

[https://gitlab.com/felix.dietrich/swimnetworks](https://gitlab.com/felix.dietrich/swimnetworks)

### Illustrative example: approximating a Barron function

We compare random Fourier features and our sampling procedure on a test function for neural networks [53]: \(f(x)=\sqrt{3/2}(\|x-a\|-\|x+a\|)\), with \(x\in\mathbb{R}^{D}\) and the constant vector \(a\in\mathbb{R}^{D}\) defined by \(a_{j}=2j/D-1\), \(j=1,2,\ldots,D\).

This function is special in that its Barron norm is equal to one for all input dimensions, and it can be represented exactly with an infinitely wide network with one hidden layer, ReLU activation, and weights uniformly distributed on a sphere of radius \(D^{1/2}\). We approximate the function using networks of up to three hidden layers. We compare relative \(L^{2}\) approximation errors over the data domain \(\mathcal{X}=[-1,1]^{2}\) for the three approximation techniques, with \(10,000\) data points sampled uniformly, separately for training and test sets. For a given approximation \(\hat{f}\), the error is defined by \(e^{2}_{\text{rel}}=\sum_{x\in\mathcal{X}}(f(x)-\hat{f}(x))^{2}/\sum_{x\in \mathcal{X}}f(x)^{2}\). For the random feature method, we use \(w\sim N(0,1),\ b\sim U(-\pi,\pi)\), as proposed in [46], and \(\phi=\sin\). For sampling, we also use \(\phi=\sin\) to obtain a fair comparison. Close to \(x=0\), the sine function is similar to \(\tanh\), for which we proved universal approximation. We also observed very similar accuracy results when trying the experiment with the \(\tanh\) function in this example. We iterate over the number of neurons \(m\) in each hidden layers, from \(m=64\) up to \(m=4096\). Figure 3 shows results in \(D=5\) and \(D=10\) (results for \(D=2,3,4\) are similar, and contained in the Appendix B). The sampled networks can be constructed approximately as fast as the random feature method (see Appendix B for the comparison). What we observe is that random features have comparable accuracy for a network with a single hidden layer, but very poor performance for deeper networks. This can be explained by the much larger ambient space dimension of the data after it is processed through the first hidden layer. Conversely, our sampling method can be used to obtain accurate results with more layers, even outperforming the case with one hidden layer. We also observe that the convergence rate at two hidden layers is faster than the theoretical rate in dimensions \(D<10\).

\begin{table}
\begin{tabular}{c c c}  & Runtime & Memory \\ \hline Assumption (I) & \(\mathcal{O}(L\cdot M(\min\{\lceil N/M\rceil,M\}+N^{2}))\) & \(\mathcal{O}(M\cdot\min\{\lceil N/M\rceil,M\}+LN^{2})\) \\ Assumption (II) & \(\mathcal{O}(L\cdot M(\lceil N/M\rceil+N^{2}))\) & \(\mathcal{O}(M\cdot\lceil N/M\rceil+LN^{2})\) \\ Assumption (III) & \(\mathcal{O}(M)\) & \(\mathcal{O}(M)\) \\ \end{tabular}
\end{table}
Table 1: Runtime and memory requirements for training sampled neural networks with the SWIM algorithm, where \(N=\max\{N_{0},N_{1},N_{2},\ldots,N_{L}\}\). Assumption (I) is that the output dimension is less than or equal to \(N\). Assumption (II) adds that \(N<M^{2}\), i.e., number of neurons and input dimension is less than the size of dataset squared. Assumption (III) requires a fixed architecture.
### Classification benchmark from OpenML

We use the "OpenML-CC18 Curated Classification benchmark" [3] to compare our sampling method with the Adam optimizer [31]. With both methods, we separately perform neural architecture search over networks with \(500\) neurons in all of their hidden layers, with the number of hidden layers from \(1\) to \(5\). We evaluate each architecture with 10-fold cross-validation. Pre-processing failed in 11 of the 72 datasets contained in the benchmark, and we restricted the maximum number of points to \(5000\) to speed up the cross-validation times. For the Adam optimizer, we employ early stopping with patience=3, monitoring the loss value. The regularization of the least squares problem in sampling was set to \(10^{-10}\). We pre-process the labels of the classification task with one-hot encoding. Figure 4 shows the benchmark results. On all benchmark tasks, sampling networks is faster than training them iteratively (on average, \(30\) times faster). The classification accuracy is comparable (cf. Figure 4, second and third plot). The ideal number of layers for each problem seems to be slightly higher for the Adam optimizer, but in several instances, sampling deeper networks seems to perform better.

### Deep neural operators

We sample several deep neural operators and compare their speed and accuracy to iterative gradient-based training of the same architectures. As a test problem, we consider the Burgers' equation on one-dimensional space, \(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu\frac{ \partial^{2}u}{\partial x^{2}},x\in(0,2\pi),t\in(0,1]\). We use periodic boundary conditions and set the viscosity parameter \(\nu=0.1\). The goal is to predict the solution at \(t=1\) from the initial condition at \(t=0\). Thus, the neural operators we construct will represent the map \(\mathcal{G}:u(x,0)\to u(x,1)\). We generate initial conditions by sampling five Fourier coefficients for the lowest frequencies and restoring the function values from these coefficients. Using a classical numerical solver, we generate 15000 pairs of \((u(x,0),u(x,1))\), and split them into the train (60%), validation (20%), and test sets (20%). Figure 5 shows samples from the generated dataset.

#### 4.3.1 Fully-connected network in signal space

The first baseline for the task is a fully-connected network (FCN) trained with tanh activation to predict the discretized solution from the discretized initial condition. We trained the classical version using the Adam optimizer and the mean squared error as a loss function. We also performed early

Figure 4: Fitting time, accuracy, and number of layers using weight sampling, compared to training with the Adam optimizer. The best architecture is chosen separately for each method and each problem, by evaluating 10-fold cross-validation error over \(1\)-\(5\) layers with 500 neurons each.

Figure 3: Relative \(L^{2}\) approximation error of a Barron function (test set), using random features and sampling, both with sine activation. Left: input dimension \(D=5\). Right: input dimension \(D=10\).
stopping based on the mean relative \(L^{2}\)-error on the validation set. **For sampling**, we use Algorithm 1 to construct a fully-connected network with tanh as the activation function.

#### 4.3.2 Fully-connected network in Fourier space

Similarly to Poli et al. [45], we train a fully-connected network in Fourier space. During training, we perform a Fourier transform on the initial condition and the solution, keeping only the lowest frequencies of the resulting vectors. We then train a standard FCN on the transformed data, as in the previous section. To compute the reported metrics, we apply an inverse Fourier transform to the network predictions and compare them to the original solutions in the dataset. **For sampling**, we split complex coefficients into real and imaginary parts before passing them to the sampled fully-connected layers. We also perform the same pre-processing step for the Adam training, to ensure a fair comparison between the methods.

#### 4.3.3 POD-DeepONet

The third architecture considered here is a variation of DeepONet, a deep operator network architecture [39]. The original DeepONet consists of two trainable components: the trunk net, which transforms the coordinates of an evaluation point, and the branch net, which transforms the function values on some grid. The outputs of these nets are then combined into the predictions of the whole network \(\mathcal{G}(u)(y)=\sum_{k=1}^{p}b_{k}(u)t_{k}(y)+b_{0},\) where \(u\) is a discretized input function; \(y\) is an evaluation point; \([t_{1},\dots,t_{p}]^{T}\in\mathbb{R}^{p}\) are the \(p\) outputs of the trunk net; \([b_{1},\dots,b_{p}]^{T}\in\mathbb{R}^{p}\) are the \(p\) outputs of the branch net; and \(b_{0}\) is a bias. DeepONet sets no restrictions on the architecture of the two nets, but often the fully-connected networks are used for one-dimensional input. POD-DeepONet proposed by Lu et al. [41] first assumes that evaluation points lie on the input grid. Then, it performs proper orthogonal decomposition (POD) of discretized solutions in the train data and uses its components instead of the trunk net to compute the outputs \(\mathcal{G}(u)(\xi_{j})=\sum_{k=1}^{p}b_{k}(u)\psi_{k}(\xi_{j})\ +\psi_{0}(\xi_{j}).\) Here \([\psi_{1}(\xi_{j}),\dots,\psi_{p}(\xi_{j})]\) are \(p\) precomputed POD components for a point \(\xi_{j}\), and \(\psi_{0}(\xi_{j})\) is the mean of discretized solutions evaluated at \(\xi_{j}\). Hence, only the branch net is trained in POD-DeepONet. We followed Lu et al. [41] and applied scaling of \(1/p\) to the network output when training the model with Adam. **For sampling**, we employ the orthogonality of the components and turn POD-DeepONet into a fully-connected network. Let \(\xi=[\xi_{1},\dots,\xi_{n}]\) be the grid used to discretize the input function \(u\) and evaluate the output function \(\mathcal{G}(u)\). Then we can compute the POD components of the train data \(\Psi(\xi)=[\psi_{1}(\xi),\dots,\psi_{p}(\xi)]\in\mathbb{R}^{n\times p}\) for this grid. If \(b(u)\in\mathbb{R}^{p}\) is the output vector of the trunk net, the POD-DeepONet transformation can be written in matrix form \(\mathcal{G}(u)(\xi)=\Psi b(u)+\psi_{0}(\xi).\) As \(\Psi^{T}\Psi=I_{p}\), we can express the output of the trunk net as \(b(u)=\Psi^{T}(\mathcal{G}(u)(\xi)-\psi_{0}(\xi)).\) Using this equation, we can transform the training data to sample a fully-connected network for \(b(u)\). Similarly to previous cases, we use tanh as the activation function for sampling.

#### 4.3.4 Fourier Neural Operator

The concept of a Fourier Neural Operator (FNO) was introduced by Li et al. [37] to solve PDEs by performing part of the training in Fourier space. An FNO consists of Fourier blocks, combining a linear operator in the Fourier space and a skip connection in signal space. As a first step, FNO lifts an input signal to a higher dimensional channel space. Let \(v_{t}\in\mathbb{R}^{n\times d_{v}}\) be an input to a Fourier block having \(d_{v}\) channels and discretized with \(n\) points. Then, the output \(v_{t+1}\) of the Fourier block is computed as \(v_{t+1}=\phi(F_{k}^{-1}(R\cdot F_{k}(v_{t}))+W\circ v_{t})\in\mathbb{R}^{n\times d _{v}}.\) Here, \(F_{k}\) is a discrete Fourier transform keeping only the \(k\) lowest frequencies and \(F_{k}^{-1}\) is the corresponding inverse transform; \(\phi\) is an activation function; \(\cdot\) is a spectral convolution; \(\circ\) is a \(1\times 1\) convolution with bias; and \(R\in\mathbb{C}^{k\times d_{v}\times d_{v}}\), \(W\in\mathbb{R}^{d_{v}\times d_{v}}\) are learnable parameters. FNO stacks several Fourier blocks and then projects the output signal to the target dimension. The projection operator, as well as the lifting operator, are usually parameterized with neural networks. **For sampling**, the construction of convolution kernels is not possible yet, so we cannot sample FNO directly. Instead, we use the general idea of FNO to construct a comparably accurate neural operator. Similarly to the original FNO, we append grid coordinates to the input signal before the lifting operation. Then, we draw the weights from a uniform distribution on [-1, 1] to compute the \(1\times 1\) lifting convolution. We also normalize the input data before the lifting. We first apply the Fourier transform to both input and target data, and then train a fully-connected network for each channel in Fourier space. We also use skip connections, as in the original FNO, by removing the input data from the lifted target function
during training, and then re-adding it before moving to the output of the block. After sampling and transforming the input data with the sampled networks in Fourier space, we apply the inverse Fourier transform. After the Fourier block(s), we sample a fully-connected network that maps the signal to the solution.

The results of the experiments in Figure 5 show that sampled models are comparable to the Adam-trained ones. The sampled FNO model does not directly follow the original FNO architecture, as we are able to only sample fully-connected layers. This shows the advantage of gradient-based methods: as of now, they are applicable to much broader use cases. These experiments showcase one of the main advantages of sampled networks: speed of training. We run sampling on the CPU; nevertheless, we see a significant speed-up compared to Adam training performed on GPU.

### Transfer learning

Training deep neural networks from scratch involves finding a suitable neural network architecture [17] and hyper-parameters [2; 55]. Transfer learning aims to improve performance on the target task by leveraging learned feature representations from the source task. This has been successful in image classification [30], multi-language text classification, and sentiment analysis [52; 7]. Here, we compare the performance of sampling with iterative training on an image classification transfer learning task. We choose the CIFAR-10 dataset [34], with 50000 training and 10000 test images. Each image has dimension \(32\times 32\times 3\) and must be classified into one of the ten classes. For this task, we consider ResNet50 [26], VGG19 [49], and Xception [9], all pre-trained on the ImageNet dataset [47]. We freeze the weights of all convolutional layers of these models and append two fully connected layers: one hidden layer and one output layer. We refer to these two layers as the classification head, which we sample with Algorithm 1 and compare the classification accuracy against iterative training with the Adam optimizer.

Figure 5: Left: Samples of initial conditions \(u_{0}\) and the corresponding solutions \(u_{1}\) for the Burgers’ equation. Right: Parameters of the best model for each architecture, the mean relative \(L^{2}\) error on the test set, and the time required to train this model. We average the metrics across three runs with different random seeds.

Figure 6: Left: Train and test accuracies with different widths for ResNet50 (averaged over 5 random seeds). Middle: Test accuracy with different models with and without fine-tuning (width = 2048). Right: Adam training and sampling times of the classification head (averaged over 5 experiments).
Figure 6 (left) shows that for a pre-trained ResNet50, the test accuracy using the sampling approach is higher than the Adam training approach for a width greater than 1024. We observe similar qualitative behavior for VGG19 and Xception (figures in Appendix E). Figure 6 (middle) shows that the sampling approach results in a higher test accuracy with all three pre-trained models. Furthermore, the deviation in test accuracy obtained with the sampling algorithm is very low, demonstrating that the sampling approach is more robust to changing random seeds than iterative training.

After fine-tuning the whole neural network with the Adam optimizer with a learning rate of \(10^{-5}\), the test accuracies of the sampled networks are close to the iterative approach. The highest test accuracy is obtained by sampling the classification head of the pre-trained Xception architecture, followed by fine-tuning. Thus, the sampled weights provide a good starting point for fine-tuning the entire model (a comparison of sampling and training approaches for the three models before and after fine-tuning for different widths is contained in Appendix E). Figure 6 (right) shows that the sampling approach can be up to 2 orders of magnitude faster than iterative training for lower widths, and is around 10 times faster for a width of 2048.

In summary, Algorithm 1 is much faster than iterative training, yields a higher test accuracy for certain widths before fine-tuning, and is more robust with respect to changing random seeds. The sampled weights also provide a good starting point for the fine-tuning of the entire model.

## 5 Broader Impact

Sampling weights through data pairs at large gradients of the target function offers improvements over random feature models. In terms of accuracy, networks with relatively large widths can even be competitive to iterative, gradient-based optimization. Constructing the weights through pairs of points also allows to sample deep architectures efficiently. Sampling networks offers a straightforward interpretation of the internal weights and biases, namely, which data pairs are important during sampling.

Given the recent critical discussions around faster advancements in artificial intelligence, publishing work that potentially speeds up the development (concretely, training speed) in this area by orders of magnitude may seem irresponsible. The solid mathematical underpinning of random feature models and, now, sampled networks, combined with much greater interpretability of the individual steps during network construction, should mitigate some of these concerns.

## 6 Conclusion

We present a data-driven sampling method for fully-connected neural networks that outperforms random feature models in terms of accuracy, and in many cases is competitive to iterative optimization. The time to obtain a trained network is orders of magnitude faster compared to iterative, gradient-based optimization. In addition, much fewer hyperparameters need to be optimized as well, mostly the regularization constant, as opposed to learning rate, number of training epochs, and type of optimizer.

Several open issues remain, we list the most pressing here. Many architectures like convolutional or transformer networks cannot be sampled with our method yet, and thus must still be trained with iterative methods. Implicit problems, such as the solution to PDE without any training data, are a challenge, as our distribution over the data pairs relies on known function values from a supervised learning setting. Iteratively refining an initially random guess may prove useful here. On the theory side, convergence rates for Algorithm 1 beyond the default Monte-Carlo estimate are not available yet, but are important for robust applications in engineering.

In the future, hyperparameter optimization, including neural architecture search, could benefit from improved training time of sampled networks. We already demonstrate benefits for transfer learning, which can now be exploited for other pre-trained models and tasks. Analyzing which data pairs are sampled during training may help to understand the datasets better. We did not show that our sampling distribution is optimal, so there is a possibility of even more efficient heuristics. Applications in scientific computing may benefit most from sampling networks, as data is often abundant, but accuracy and speed requirements are much higher than for classical tasks in machine learning.
## Acknowledgments and Disclosure of Funding

The authors are grateful for discussions with Erik M. Bollt and Constantinos Siettos. E.B., I.B., and F.D. are funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - project no. 468830823. C.D. is partially funded by the Institute for Advanced Study (IAS) at the Technical University of Munich. Q.S. is supported by the TUM Georg Nemetschek Institute - Artificial Intelligence for the Built World.
