# What is the collection of functions that a given finite neural network can approximate with ease?

To my understanding, one version of the universal approximation theorem runs as follows: Let $$\Phi$$ be the family of (trained) feedforward neural networks of bounded width, arbitrary depth, and mild constraints on the activation functions. Let $$F$$ be the space of continuous functions $$f: X \rightarrow \mathbb{R}$$ where $$X$$ is a compact subset of $$\mathbb{R}^n$$ such as $$[0,1]^n$$. Then for any $$f \in F$$ and $$\varepsilon > 0$$, there exists an $$\hat{f} \in \Phi$$ such that for all $$x \in X$$, $$\left|\hat{f}(x) - f(x)\right| < \varepsilon.$$

At first glance, this might seem to suggest that for any given $$f \in F$$ we could select an untrained finite-node neural network of sufficient size and train it by the usual methods to closely approximate $$f$$.

And indeed, this works well for many $$f$$. Suppose that we construct $$m$$ samples of the form $$x_k \in [0,1]^n$$ and $$y_k \in \mathbb{R}$$, where $$y_k = f(x_k)$$ and $$m = 3000$$. Additionally, suppose that $$n = 30$$, and $$x_k$$ is drawn independently and uniformly at random from $$[0,1]^n$$. Then it is straightforward to train a deep neural network (DNN) with 4 hidden layers and $$n$$ nodes per hidden layer to accurately estimate any of the following functions: $$\begin{array}{} f(x_k) = x_{k,1}x_{k,2} + x_{k,3} - x_{k,4}/2, \\[0mm] f(x_k) = x_{k,14}\log(x_{k,1}+3) + x_{k,3} - x_{k,9}^{1/3}/2, \\[0mm] f(x_k) = \text{sgn}(x_{k,1} - 1/2) |x_{k,1} - 1/2|^{1/5}, \\[0mm] f(x_k) = (x_{k,1}^2 + x_{k,2} - 11)^2 + (x_{k,1} + x_{k,2}^2 - 7)^2, \quad\text{or} \\[0mm] f(x_k) = x_{k,1}^{\Large x_{k,2}^{\LARGE x_{k,3}^{\LARGE x_{k,4}^{\LARGE x_{k,5}}}}}. \end{array}$$

For instance, below are the results of the DNN training and estimation for the first function $$f(x_k)$$ above. The figure was generated with the Python script at the end of this question.

However, for various other $$f$$, training a DNN to accurately predict $$f(x_k)$$ from $$x_k$$ under the above conditions appears far from easy. For example, if again $$n = 30$$, $$m = 3000$$, and $$x_k$$ is selected randomly from $$[0,1]^n$$ as before, it appears much harder to train a neural network to estimate the value of any of the following functions:

$$\begin{array}{} f(x_k) = \sum_{j=1}^n (x_{k,j}-1/2)^8, \\[0mm] f(x_k) = \sum_{j=1}^n \exp(-20x_{k,j}), \\[0mm] f(x_k) = \sum_{j=1}^n \max(1/10-x_{k,j}, 0), \\[0mm] f(x_k) = \min_{1 \leq j \leq n} x_{k,j}, \;\text{or} \\[0mm] f(x_k) = \sin(8\pi x_{k,1}). \\[0mm] \end{array}$$

For the first $$f(x_k)$$ above, for instance, the results of the DNN training and estimation are shown below.

It may be worth noting that the first four of the above five examples can perhaps be viewed as variants of the same example — they each might be regarded as loosely similar to a summation over either the largest, smallest, or largest and smallest elements of $$x_k$$. The fifth example, $$f(x_k) = \sin(8\pi x_{k,1})$$, appears to be a different type of example and exhibits an abrupt transition from nearly random estimation ($$r^2 \approx 0$$) to accurate estimation ($$r^2 \approx 1$$) as the coefficient $$8\pi$$ is decreased toward zero.

DNN training also fails for an even simpler $$f(x_k)$$: Choose $$y_k$$ uniformly at random from $$[0,1]$$ and define one of the variables of $$x_k$$ as a superposition of another variable of $$x_k$$ and $$y_k$$ multiplied by a sufficiently small constant. For instance, let $$x_{k,n} = x_{k,1} + \varepsilon y_k$$ where $$\varepsilon = 0.05$$. In this case, $$x_{k,n}$$ is no longer drawn independently and uniformly from $$[0,1]$$, but $$y_k$$ remains a deterministic and continuous function of $$x_k$$. Like the sinusoidal $$f(x_k)$$ above, this example exhibits an abrupt transition; as $$\varepsilon$$ is increased from $$0.05$$ to $$0.1$$, $$r^2$$ grows from near $$0$$ to near $$1$$.

(Critically, not all functions fall neatly into one of the above two camps. For instance, the described DNN predicts $$f(x_k) = \prod_{j=1}^n x_{k,j}^{1/n}$$ with an $$r^2$$ of roughly $$0.5$$.)

Taken together, these examples and observations leave me with several closely related questions:

• Are there other interesting examples of continuous, closed-form functions that neural networks struggle to approximate — especially examples that cannot be viewed as variants of the above.
• Intuitively, why do finite-node neural networks appear to have great difficulty approximating some continuous, closed-form functions and no difficulty approximating others?
• Are there any known theorems or other results that demarcate the hard-to-approximate functions from those that are easy to approximate?

(Answers to these questions might also apply to this question, which seems related.)

Python script

Below is a self-contained Python script that can be used to confirm the above observations by commenting in and out the functions in the for loop.

import os
os.environ['TF_ENABLE_ONEDNN_OPTS'] = '0'
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '1'
import numpy as np
from time import perf_counter
from numpy.random import random
from sklearn.preprocessing import StandardScaler
from keras.models import Sequential
from keras.layers import Input, Dense
import matplotlib.pyplot as plt

# ----------------------------------------------------------------------------------------------------
# Define constants and settings

n = 30  # number of elements in x_k
m = 3000  # number of samples used to train DNN
kc = 200  # Only samples with an index k > kc are used to train DNN.  The remaining samples are used to test DNN.
max_epochs = 300  # number of epochs to use when training DNN
pc0 = perf_counter()  # used to measure runtime

# ----------------------------------------------------------------------------------------------------
# Construct the input variables
# (The functions below were primarily designed for the case that n = 30 and m = 3000.)

# Construct x and y
x = []
y = []
for k in range(m):
xk = np.array(random(n))

# Functions for which DNN estimates are accurate (r_sqr of about 1)
# yk = xk[0]*xk[1] + xk[2] - xk[3]/2
# yk = xk[1]*np.log(xk[0]+3) + xk[2] - xk[3]**(1/3)/2
# yk = np.sign(xk[0] - 1/2) * abs(xk[0] - 1/2)**(1/5)
# yk = (xk[0]**2 + xk[1] - 11)**2 + (xk[0] + xk[1]**2 - 7)**2
# yk = xk[0]**xk[1]**xk[2]**xk[3]**xk[4]

# Functions for which DNN estimates are nearly random (r_sqr near 0)
yk = sum((xk-1/2)**8)
# yk = sum(np.exp(-20*xk))
# yk = sum([max(1/10-xk[j], 0) for j in range(n)])
# yk = min(xk)
# yk = np.sin(8*np.pi*xk[0])

# An additional example for which DNN estimates are nearly random
# xkr = np.array(random(n-1))
# yk = random()
# xkc = np.array([xkr[0] + 0.05*yk])
# xk = np.concatenate((xkr, xkc))

# A function for which DNN estimates are semi-accurate
# yk = np.prod(xk)**(1/n)  # r_sqr of about 0.5

x.append(xk)
y.append(yk)

x = np.array(x)
y = np.array([[y[k]] for k in range(m)])

# Standardize and split x and y
x = StandardScaler().fit_transform(x)
y = StandardScaler().fit_transform(y)
x_test, x_train = x[:kc], x[kc:]
y_test, y_train = y[:kc], y[kc:]

print('---')
print('Input statistics')
print('  Number of variables in each xk: ' + str(x.shape[1]))
print('  Number of samples in x: ' + str(x.shape[0]))
print('  Number of variables in each yk: ' + str(y.shape[1]))
print('  Number of samples in y: ' + str(y.shape[0]))
xk = [x[k][0] for k in range(m)]
print('  Min of first variable of x: ' + str(round(min(xk), 4)))
print('  Max of first variable of x: ' + str(round(max(xk), 4)))
print('  Mean of first variable of x: ' + str(round(np.mean(xk), 4)))
yk = [y[k][0] for k in range(m)]
print('  Min of y: ' + str(round(min(yk), 4)))
print('  Max of y: ' + str(round(max(yk), 4)))
print('  Mean of y: ' + str(round(np.mean(yk), 4)))
print('  Number of samples in x_train: ' + str(x_train.shape[0]))
print('  Number of samples in x_test: ' + str(x_test.shape[0]))
print('  Number of samples in y_train: ' + str(y_train.shape[0]))
print('  Number of samples in y_test: ' + str(y_test.shape[0]))

# ----------------------------------------------------------------------------------------------------
# Train and test the DNN

# Construct, configure, and train the DNN
dnn = Sequential()
dnn.add(Dense(units=n, activation='relu'))  # first hidden layer
dnn.add(Dense(units=n, activation='relu'))  # second hidden layer
dnn.add(Dense(units=n, activation='relu'))  # third hidden layer
dnn.add(Dense(units=n, activation='relu'))  # fourth hidden layer

# Configure the DNN
loss_func = 'mean_absolute_error'
dnn.compile(loss=loss_func, optimizer=optimizer_algorithm)

# Train the DNN
history = dnn.fit(x_train, y_train, validation_data=(x_test,y_test), verbose=0, epochs=max_epochs)

# Use the DNN to estimate the values in y_test
y_test_pred = dnn.predict(x_test, verbose=0)

# ----------------------------------------------------------------------------------------------------
# Analyze and plot the DNN estimates

# Compute r (Pearson correlation coefficient) and r_sqr
mean_y_test = np.mean(y_test)
mean_y_test_pred = np.mean(y_test_pred)
r = sum((y_test - mean_y_test) * (y_test_pred - mean_y_test_pred))[0] / \
np.sqrt(sum((y_test - mean_y_test)**2)[0] * sum((y_test_pred - mean_y_test_pred)**2)[0])
r_sqr = r**2

print('---')
print('Output statistics')
print('  Final loss value: ' + str(round(history.history['loss'][-1], 4)))
print('  Final val_loss value: ' + str(round(history.history['val_loss'][-1], 4)))
print('  mean_y_test = ' +  str(round(mean_y_test, 4)))
print('  r = ' + str(round(r, 4)))
print('  r_sqr = ' + str(round(r_sqr, 4)))

pc1 = perf_counter()  # We stop measuring runtime at this point.
print('---')
print('Runtime statistics')
print('  ' + str(round(pc1-pc0, 1)) + ' sec total runtime.')
print('---')

# Make figure
plt.figure(num=1, figsize=(11,5))
mngr = plt.get_current_fig_manager()
mngr.canvas.manager.window.wm_geometry("+%d+%d" % (100, 150))

# Plot loss and val_loss versus epoch
plt.subplot(1, 2, 1)
ymax = 1.05*max(max(history.history['val_loss'][1:]), max(history.history['loss'][1:]))
plt.plot(history.history['val_loss'], label='val_loss', color=[0.9, 0.5, 0.1])
plt.plot(history.history['loss'], label='loss', color=[0.1, 0.4, 0.8])
plt.gca().set_title('loss function: ' + loss_func + ', optimizer: ' + optimizer_algorithm + '\n', fontsize=11, pad=0)
plt.xlabel('Epoch', fontsize=10)
plt.ylabel('Loss histories', fontsize=10)
plt.legend(frameon=False, loc='best')
plt.gca().set_box_aspect(1)
plt.ylim([0, ymax])

# Plot scatter of estimates versus true values
plt.subplot(1, 2, 2)
ymin = min(np.concatenate([np.array([-0.05]), min(y_test), min(y_test_pred)]))
ymax = max(np.concatenate([np.array([1]), max(y_test), max(y_test_pred)]))
lims = [ymin, ymax]
plt.plot(lims, lims, linewidth=0.75, color=[0.75, 0.75, 0.75])
plt.plot(lims, [0, 0], linewidth=0.75, color=[0.75, 0.75, 0.75])
plt.plot([0, 0], lims, linewidth=0.75, color=[0.75, 0.75, 0.75])
plt.scatter(y_test, y_test_pred, s=1, edgecolors='none')
plt.gca().set_title('n = ' + str(n) +', m = ' + str(m) + ', r_sqr = ' + str(round(r_sqr, 4)) + '\n', fontsize=11, pad=0)
plt.xlabel('True values', fontsize=10)
plt.ylabel('DNN estimates', fontsize=10)
plt.gca().set_box_aspect(1)
plt.xlim(lims)
plt.ylim(lims)

plt.show()


I think the numbers of effective parameters needed to approximate these functions are very different. With a fixed sample size $$m=3000$$, there could be insufficient information to estimate too many effective parameters. This seems to be, at least in part, a demonstration of the curse of dimensionality.
Take the $$\sin(8\pi x_{k,1})$$ example, with $$8\pi$$ replaced by a small value, a linear or quadratic polynomial might be a good enough approximation; but with $$8\pi$$ bending the curve back and forth, a lot more effective parameters are needed to encode where and how the curve should bend, not to say the noise added from other components of $$x_k$$.
As a simple test, $$n$$ might be reduced to $$2$$ to see if the performance improves.
• Thanks for these comments and ideas. I think they provide useful clarification and move us closer to an intuitive explanation. If we decrease $n$ to $5$ (so just 126 parameters) and increase the $-20$ inside of the exponential to $-40$, we obtain a similar outcome as described above: the first set of five functions is still much more accurately estimated with a DNN than the second set. I believe this leaves open the original question of whether we can determine from analysis (rather than training) whether a function $f \in F$ can be accurately approximated by a given finite neural network. Commented Jul 17 at 3:44
• Based on your answer, it seems like one path to approach this would be to try to compute from $f$ an estimate of or bound on the number of "effective parameters" -- or other measurement of neural network size -- needed to accurately approximate $f$. Are any results of this sort known, even for just a neural network of a single hidden layer? Commented Jul 17 at 3:44
• @SapereAude I agree with you that some functions are harder to approximate, when we don't know the functional form. Did you still see $R^2\approx 0$ after bringing the dimensions down? If it improved, then dimensionality would at least a partial reason. I have to leave it for other site users to recommend relevant references. Commented Jul 17 at 17:35
• For the second set of five functions, decreasing $n$ to $5$ increases $R^2$ very roughly from $0.1$ to $0.6$ -- a big jump! So the curse of dimensionality is a partial reason for the quantitative difference between the $R^2$ values of the two sets of five functions. However, the first five functions have an $R^2$ of $~0.99$ for $n =5, 30$, so we could also say the curse of dimensionality explains very little about the qualitative discrepancy in $R^2$ between the two sets of functions. To explain this, some $f \in F$ must have a property which makes them easier for a DNN to approximate. Commented Jul 19 at 9:45