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(Input(shape=x_train[0].shape)) # input layer
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
dnn.add(Dense(units=1)) # output node
# Configure the DNN
loss_func = 'mean_absolute_error'
optimizer_algorithm = 'adam'
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.subplots_adjust(bottom=0.15, top=0.85, left=0.1, right=0.9, wspace=0.2)
plt.show()
