This package aims to implement a series of operations first described in Grohs et al. (2023), Petersen and Voigtlaender (2018), Grohs et al. (2022), Jentzen et al. (2023), and a broad extension to that framework of operations as seen in Rafi et al. (2024). Our main definitions will be from Rafi et al. (2024), but we will also delve deeper into the literature when necessary.
Our definition of neural networks will be ordered tuples of ordered pairs. A neural network is something like \(((W_1,b_1),(W_2,b_2), (W_3,b_3))\). Where each \(W_i\) is a weight matrix and each \(b_i\) is a bias vector. Create a neural network by supplying a numeric vector containing the number of neurons in each layer.
layer_architecture = c(4,5,6,5)
create_nn(layer_architecture)
#> [[1]]
#> [[1]]$W
#> [,1] [,2] [,3] [,4]
#> [1,] 0.7254202 0.4951218 -0.2429171 1.94309749
#> [2,] -0.1923639 -1.0268889 -0.3629262 -0.78896620
#> [3,] -0.2235433 -0.8877182 0.3747027 -0.99155173
#> [4,] 2.2619089 -0.6051185 1.0304409 -0.32990715
#> [5,] -0.3086433 -0.4449017 0.4229585 0.01492379
#>
#> [[1]]$b
#> [,1]
#> [1,] -0.7077897
#> [2,] -1.4321714
#> [3,] -1.0142027
#> [4,] -1.2299353
#> [5,] 0.3883264
#>
#>
#> [[2]]
#> [[2]]$W
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.5993240 0.1508196 0.7894957 0.5172492 0.7288341
#> [2,] 1.8225499 0.5283804 -0.2812573 -0.9892737 2.0574180
#> [3,] -0.0998676 -1.2374602 -0.4842872 -1.1660346 -0.8422486
#> [4,] 1.1164926 2.1424882 1.5219002 -0.6328137 0.1148362
#> [5,] -0.2622076 -0.1116592 -0.2900375 0.7329308 0.6834428
#> [6,] -1.8538367 0.7130021 -0.0697573 0.9475905 -1.4684740
#>
#> [[2]]$b
#> [,1]
#> [1,] -0.5961234
#> [2,] 1.7000471
#> [3,] 0.2384115
#> [4,] -0.2099502
#> [5,] -1.1450410
#> [6,] -0.3215038
#>
#>
#> [[3]]
#> [[3]]$W
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1.45857543 -0.2874806 -0.9408067 -0.2894312 -0.3512325 0.9282448
#> [2,] 1.15932595 0.7232562 -1.6152553 -0.5257204 -0.6376106 -1.0398877
#> [3,] 0.43857269 1.7720217 1.5960320 1.2630813 -0.4903045 0.2466736
#> [4,] 0.94583160 -0.0123493 -0.5306671 0.5577958 1.3712969 0.0293787
#> [5,] -0.08461945 0.4327464 -1.3102911 -0.4846275 -0.4935751 1.1897589
#>
#> [[3]]$b
#> [,1]
#> [1,] -0.6799475
#> [2,] -0.8055017
#> [3,] 1.0975811
#> [4,] -0.3931456
#> [5,] 0.1918469Note that each weight matrix and bias vector is populated from a standard normal distribution. Associated with each neural network will be a family of functions:
dep: showing the depth of a neural network:
param: showing the number of parameters (sum of the
total number of elements in each weight matrix and bias vector):
inn: the width of the input layer:
out: the width of the output layer:
hid: the number of hidden layers, always defined to
be one less than depth:
lay: a list specifying the width of each layer,
often will be called the layer architecture:
Instantiation refers to the act of applying an activation function between each layer and resulting in a continuous function. Three activation functions are included: ReLU, sigmoid, and hyperbolic tangent.
Because the current theoretical frameworks of Rafi et al. (2024), Jentzen et al. (2023), Grohs et al. (2023) only show theoretical results for ReLU activation, our Xpn, Sne, and Csn, among others will only show approximations under ReLU activations.
Instantiations must always be accompanied by the appropriate vector \(x\) of the same length as the input layer of the instantiated neural network. See examples:
create_nn(c(1, 3, 5, 6)) |> inst(ReLU, 8)
#> [,1]
#> [1,] 0.4570365
#> [2,] -0.6990080
#> [3,] -1.3157376
#> [4,] 0.7749714
#> [5,] 1.3164407
#> [6,] -0.5154481A numeric matrix is evaluated as a batch, with one input sample per
column. For repeated evaluation, realize_nn() validates the
network once and returns a reusable function.
nn <- create_nn(c(2, 4, 1))
batch <- matrix(c(1, 2, 3, 4), nrow = 2)
inst(nn, ReLU, batch)
#> [,1] [,2]
#> [1,] 0.8751035 0.8751035
f <- realize_nn(nn, ReLU)
f(batch)
#> [,1] [,2]
#> [1,] 0.8751035 0.8751035Instantiation will have a special symbol \(\mathfrak{I}\), and instantiation with ReLU will be denoted as \(\mathfrak{I}_{\mathsf{ReLU}}\).
Composition has the wonderful property that it works well with instantiation, i.e. instantiation of two composed neural networks are the same as the composition of the instantiated continuous functions, i.e.
\[ \mathfrak{I}_{\mathsf{ReLU}} \left( \nu_1 \bullet \nu_2\right)(x) = \mathfrak{I}_{\mathsf{ReLU}}\left( \nu_1\right) \circ \mathfrak{I}_{\mathsf{ReLU}}\left( \nu_2\right)(x) \]
Note: When composing, the output layer width of the innermost neural network must match the input layer width of the outer neural network.
c(1,5,6,3,3) |> create_nn() -> nu_1
c(3,4,6,3,2) |> create_nn() -> nu_2
nu_2 |> comp(nu_1)
#> [[1]]
#> [[1]]$W
#> [,1]
#> [1,] 2.3090371
#> [2,] 0.2389431
#> [3,] 0.4243359
#> [4,] 0.4688192
#> [5,] 0.3594160
#>
#> [[1]]$b
#> [,1]
#> [1,] 0.06481649
#> [2,] -0.14874294
#> [3,] 0.50669375
#> [4,] 0.26877941
#> [5,] 0.54888395
#>
#>
#> [[2]]
#> [[2]]$W
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.4504625 -0.85454553 0.4859195 -0.45763865 -0.8617718
#> [2,] -0.2065270 -0.06683044 0.9597480 -0.55238328 0.1700360
#> [3,] 0.8320086 0.03759043 -1.2904293 0.07695865 0.0323700
#> [4,] 0.6388684 1.10434640 0.3809041 -0.04362206 -1.9293359
#> [5,] 1.7297790 -1.00619833 1.5690024 -0.92953061 1.1665997
#> [6,] 0.9211430 0.41944605 1.2952370 -3.51554955 -0.2919620
#>
#> [[2]]$b
#> [,1]
#> [1,] -0.6438372
#> [2,] 1.5907587
#> [3,] 0.5518107
#> [4,] -1.4367229
#> [5,] 1.6571493
#> [6,] -0.8641815
#>
#>
#> [[3]]
#> [[3]]$W
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] -1.6524714 0.1642487 -0.7866409 1.478168 -0.9724772 -1.0221943
#> [2,] -0.7772266 -1.4475563 -0.1504985 -1.078204 0.6415636 -1.2042211
#> [3,] -0.1009975 -1.5028978 2.3356442 1.451244 0.6282726 -0.3406182
#>
#> [[3]]$b
#> [,1]
#> [1,] 0.76544544
#> [2,] 1.16369324
#> [3,] 0.04649973
#>
#>
#> [[4]]
#> [[4]]$W
#> [,1] [,2] [,3]
#> [1,] 0.5452922 0.32927408 0.4874979
#> [2,] -0.3209399 -0.09439638 0.2157386
#> [3,] -0.1866712 0.19523165 -0.6000300
#> [4,] -1.0134196 -0.31059285 -0.6481413
#>
#> [[4]]$b
#> [,1]
#> [1,] -2.6468603
#> [2,] 0.4068598
#> [3,] 1.8800760
#> [4,] 1.4005222
#>
#>
#> [[5]]
#> [[5]]$W
#> [,1] [,2] [,3] [,4]
#> [1,] 0.2901188 0.3535010 0.5005883 1.0264827
#> [2,] -0.4641124 1.0684200 -0.7910108 -0.7690487
#> [3,] -0.8284524 -1.2911631 -0.2212950 -1.0071846
#> [4,] -0.9667745 2.0256973 -1.1863459 0.7351982
#> [5,] 0.2284624 1.5531469 0.7945425 0.9380483
#> [6,] 0.2532401 0.8508438 -0.3855570 0.6117359
#>
#> [[5]]$b
#> [,1]
#> [1,] -0.90839368
#> [2,] 1.09902267
#> [3,] 0.64353422
#> [4,] -0.67048742
#> [5,] -0.07668242
#> [6,] -1.42032524
#>
#>
#> [[6]]
#> [[6]]$W
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] -0.5302223 2.09746241 0.4695967 0.5393372 1.2342237 -0.3123669
#> [2,] 0.3215652 -0.07544237 -2.4492050 -1.2693528 -0.6795681 -0.1651806
#> [3,] -2.0981813 0.45166960 1.2370330 0.1265308 -1.5072000 1.0882599
#>
#> [[6]]$b
#> [,1]
#> [1,] -1.4595610
#> [2,] 1.1000563
#> [3,] -0.1324809
#>
#>
#> [[7]]
#> [[7]]$W
#> [,1] [,2] [,3]
#> [1,] 0.3867486 -1.4995651 0.9513026
#> [2,] 1.4356521 0.7746126 1.2857973
#>
#> [[7]]$b
#> [,1]
#> [1,] 0.4969362
#> [2,] -0.4248948Given a neural network you may perform scalar left multiplication on a neural network. They instantiate in quite nice ways. Suppose if the neural network instantiated as \(\mathfrak{I}_{\mathsf{ReLU}}\left( \nu \right)(x)\), scalar left multiplication instantiates as:
\[ \mathfrak{I}_{\mathsf{ReLU}} \left( \lambda \triangleright \nu\right)(x) = \lambda \cdot \mathfrak{I}_{\mathsf{ReLU}}\left( \nu\right)(x) \]
Scalar left multiplication instantiates as:
\[ \mathfrak{I}_{\mathsf{ReLU}} \left(\nu\triangleleft \lambda\right)(x) = \mathfrak{I}_{\mathsf{ReLU}}\left(\nu\right)(\lambda \cdot x) \]
Here is the R code for this, compare the two outputs:
c(1,3,4,8,1) |> create_nn() -> nn
nn |> inst(ReLU,5)
#> [,1]
#> [1,] 12.69314
2 |> slm(nn) |> inst(ReLU, 5)
#> [,1]
#> [1,] 25.38628And now for scalar right multiplication, although the difference in output is not as obvious:
Neural networks may also be stacked on top each other. The instantiation of two stacked neural networks is the concatenation of the outputs of the instantiated networks individually. In other words, mathematically we may say that:
\[ \mathfrak{I}_{\mathsf{ReLU}}\left( \nu_1 \boxminus \nu_2\right)\left( \left[x \quad y\right]^\intercal\right) = \left[ \mathfrak{I}_{\mathsf{ReLU}}\left( \nu_1\right) \left( x\right)\quad \mathfrak{I}_{\mathsf{ReLU}}\left( \nu_2\right)\left( y\right)\right]^\intercal \]
To make this more concrete observe that:
c(3,4,6,3,7,1,3,4) |> create_nn() -> nn_1
c(2,6,4,5) |> create_nn() -> nn_2
(nn_1 |> stk(nn_2)) |> inst(ReLU, c(4,3,2,1,6))
#> [,1]
#> [1,] 0.9167291
#> [2,] -0.6516221
#> [3,] 3.3522944
#> [4,] 0.4923644
#> [5,] 7.0903153
#> [6,] 22.1888272
#> [7,] 6.0859709
#> [8,] -2.7866794
#> [9,] 2.6563157
print("Compare to:")
#> [1] "Compare to:"
nn_1 |> inst(ReLU, c(4,3,2))
#> [,1]
#> [1,] 0.9167291
#> [2,] -0.6516221
#> [3,] 3.3522944
#> [4,] 0.4923644
nn_2 |> inst(ReLU, c(1,6))
#> [,1]
#> [1,] 7.090315
#> [2,] 22.188827
#> [3,] 6.085971
#> [4,] -2.786679
#> [5,] 2.656316Note Stacking of unequal depth happens automatically
by something called “tunneling”. Use stk_many(list(...)) to
stack any number of networks in a single operation.
Neural networks may be added such that the continuous function created under ReLU instantiation is the sum of the individual instantiated functions, i.e.
\[ \mathfrak{I}_{\mathsf{ReLU}}\left( \nu \oplus \mu\right)\left( x\right) = \mathfrak{I}_{\mathsf{ReLU}}\left( \nu \right)(x) + \mathfrak{I}_{\mathsf{ReLU}}\left(\mu\right)\left( x\right) \]
The code works as follows:
c(1,5,3,2,1) |> create_nn() -> nu
c(1,5,4,9,1) |> create_nn() -> mu
nu |> inst(ReLU,4) -> x_1
mu |> inst(ReLU,4) -> x_2
x_1 + x_2 -> result_1
print("The sum of the instantiated neural network is:")
#> [1] "The sum of the instantiated neural network is:"
print(result_1)
#> [,1]
#> [1,] 2.34707
(nu |> nn_sum(mu)) |> inst(ReLU,4) -> result_2
print("The instantiation of their neural network sum")
#> [1] "The instantiation of their neural network sum"
print(result_2)
#> [,1]
#> [1,] 2.34707Use nn_sum_many(list(...)) to sum any number of
compatible networks in one operation.
Now that we have a basic operations for neural networks, we are able to go into more sophisticated functions. We start with some basic
We have neural networks for approximating squaring any real number. These only work with ReLU instantiates and most results in this field as well as most of the literature focuses on ReLU.
These must be supplied with two arguments, \(q\in (2,\infty)\) and \(\varepsilon \in (0,\infty)\). Accuracy increases the closer we are the \(2\) and \(\varepsilon\) respectively.
See the examples:
Similarly we may define the product neural network which approximates the product of two real numbers, given \(q \in (2,\infty)\), and \(\varepsilon \in (0,\infty)\). Accuracy increases the closer we are to \(2\) and \(\varepsilon\) respectively. These must be instantiated with a list of two real numbers:
Repeated applications of the \(\mathsf{Prd}^{q,\varepsilon}\). Will give us neural networks for approximating raising to a power.
These are called using the \(\mathsf{Pwr}^{q,\varepsilon}\) neural networks. These need three arguments \(q\in (2,\infty)\), \(\varepsilon \in (0,\infty)\), and \(n \in \mathbb{N}\), the power to which we are approximating. Power networks are constructed iteratively and reused, avoiding repeated recursive construction. Here is an example:
Pnm() implements a neural-network analogue of \(c_0 + c_1x + \cdots + c_nx^n\). Supply
coefficients in ascending power order. Power networks are constructed
once, padded to equal depth, scaled, and summed with one n-ary
operation.
We can do neural network sums, scalar left multiplication, and raising to a power. We thus have enough technology to create neural network polynomials, and more specifically finite power series approximations of common functions.
This is the neural network for approximating \(e^x\). These need three arguments \(q\in (2,\infty)\), \(\varepsilon \in (0,\infty)\), and \(n \in \mathbb{N}\), the power to which we will truncate the power series expansion. This may require significant computation time, depending on the power to which we are approximating. Here is an example of the code
Xpn(5,2.1,0.1) |> inst(ReLU, 2)
#> [,1]
#> [1,] 7.243428
print("Compare to:")
#> [1] "Compare to:"
exp(2)
#> [1] 7.389056By far the biggest improvement in accuracy will come from increasing the power series truncation limit. This will also contribute the most to computation time.
This is the neural network for approximating \(\cos(x)\). These need three arguments \(q\in (2,\infty)\), \(\varepsilon \in (0,\infty)\), and \(n \in \mathbb{N}\), the power to which we will truncate the power series expansion. This may require significant computation time, depending on the power to which we are approximating. Here is an example of the code
Csn(3,2.1,0.1) |> inst(ReLU, 0.4)
#> [,1]
#> [1,] 0.9212815
print("Compare to:")
#> [1] "Compare to:"
cos(0.4)
#> [1] 0.921061By far the biggest improvement in accuracy will come from increasing the power series truncation limit. This will also contribute the most to computation time.
This is the neural network for approximating \(\sin(x)\). These need three arguments \(q\in (2,\infty)\), \(\varepsilon \in (0,\infty)\), and \(n \in \mathbb{N}\), the power to which we will truncate the power series expansion. This may require significant computation time, depending on the power to which we are approximating. Here is an example of the code
Sne(3,2.1,0.1) |> inst(ReLU, 0.4)
#> [,1]
#> [1,] 0.3895701
print("Compare to:")
#> [1] "Compare to:"
sin(0.4)
#> [1] 0.3894183By far the biggest improvement in accuracy will come from increasing the power series truncation limit. This will also contribute the most to computation time.
A simple trapezoidal approximation can be done with neural networks given two sample points along the two legs of a trapezoid. These may be instantiated with any of the three activation functions implemented, ReLU, Sigmoid, or Tanh Observe:
h = 0.2
mesh = c(0,0+h)
samples = sin(mesh)
Trp(h) |> inst(ReLU, samples)
#> [,1]
#> [1,] 0.01986693
Trp(h) |> inst(Sigmoid, samples)
#> [,1]
#> [1,] 0.01986693
Trp(h) |> inst(Tanh, samples)
#> [,1]
#> [1,] 0.01986693We may extend this to what we call the “extended trapezoid”, which approximates the area under a curve \(f(x)\) using the trapezoidal rule once instantiated with function sample values from the full mesh. These require two parameters \(N \in \mathbb{N}\), representing the number of trapezoids that we want to split the area under \(f(x)\) into, and \(h\) the mesh width. Note we will always have one less trapezoid than the number of meshpoints.
These may be instantiated with any of the three activation functions implemented, ReLU, Sigmoid, or Tanh.
Observe:
Suppose you have a function \(f:\mathbb{R}^d \rightarrow \mathbb{R}\) with Lipschitz constant \(L\) with respect to the \(\ell^1\) norm.
Take sample points \(\{x_1,x_2,...,x_n\}\) within the domain. Then let \(f_1,f_2,...,f_n\) be a family of functions, define for all \(i \in \{1,2,...,n\}\) as: \[ f_i(x) = f(x_i) - L\lVert x-x_i\rVert_1 \]
These would create little “hills” with the tip of the hill being at whatever points on the function where the samples were take from. We may “saw off” the base of these hills giving us a sawtooth-like function that approximates our function \(f(x)\), i.e.
\[ \hat{f} (x) = \max_{i \in \{1,2,...,n\}} \left\{ f_i \left( x\right)\right\} \]
We will call this the “maximum convolution approximation”. The stated maximum-convolution construction and guarantee require ReLU activation.
In one dimension, vectors can still be supplied for backward compatibility:
seq(0, 2 * pi, length.out = 500) -> x
sin(x) -> y
1 -> L
MC(x, y, L) |> inst(ReLU, 2.5)
#> [,1]
#> [1,] 0.5970914
print("Compare to:")
#> [1] "Compare to:"
sin(2.5)
#> [1] 0.5984721For \(d>1\), put one sample point in each column of a matrix. The same convention is used when evaluating a batch.