import numpy as np
from scipy.integrate import quad
def f(x):
return np.exp(-2 * (x - 1)**2) + 0.6 * np.exp(-4 * (x + 1.5)**2)
def a(n):
return quad(lambda x: f(x) * np.cos(n * x), -np.pi, np.pi)[0] / np.pi
def b(n):
return quad(lambda x: f(x) * np.sin(n * x), -np.pi, np.pi)[0] / np.pi
def fourier(x, N):
y = a(0) / 2 * np.ones_like(x)
for n in range(1, N + 1):
y += a(n) * np.cos(n * x) + b(n) * np.sin(n * x)
return yInterpretations of what a neural network is
machine learning
python
A Fourier series is a one-hidden-layer neural network, and every neural network is a learned nonlinear feature transform followed by a linear model.
There are many ways to think about what a neural network is. In this series we look at a few of them. Each one is a little different, and each one makes some property of neural networks easier to understand.
Part 1: a basis set expansion
Let’s start with something familiar: the Fourier series. On the interval \([-\pi, \pi]\) we can write a function as
\[ f(x) \approx \frac{a_0}{2} + \sum_{n=1}^{N} \left[ a_n \cos(n x) + b_n \sin(n x) \right] \]
This is a basis set expansion: \(f\) is approximated by a weighted sum of fixed functions \(\{1, \cos x, \sin x, \cos 2x, \sin 2x, \ldots\}\). The only unknowns are the coefficients \(a_n\) and \(b_n\).
Computing the coefficients
In general, if we expand \(f\) in some basis \(\phi_k(x)\),
\[ f(x) \approx \sum_{k} c_k \phi_k(x), \]
we find the coefficients by minimizing the squared error \(\int (f - \sum_k c_k \phi_k)^2 dx\). Setting the derivative with respect to each \(c_j\) to zero gives the normal equations
\[ \sum_k \langle \phi_j, \phi_k \rangle \, c_k = \langle \phi_j, f \rangle \qquad \text{or} \qquad \mathbf{G}\mathbf{c} = \mathbf{b}, \]
where \(\langle g, h \rangle = \int_{-\pi}^{\pi} g(x) h(x) dx\) is the inner product, \(G_{jk} = \langle \phi_j, \phi_k \rangle\) is the Gram matrix, and \(b_j = \langle \phi_j, f \rangle\). For an arbitrary basis this is a linear regression problem: every coefficient depends on every other one through \(\mathbf{G}\).
The sines and cosines are special because they are orthogonal on \([-\pi, \pi]\):
\[ \begin{aligned} \int_{-\pi}^{\pi} \cos(m x) \cos(n x) \, dx &= \pi \, \delta_{mn} \quad (m, n \ge 1) \\ \int_{-\pi}^{\pi} \sin(m x) \sin(n x) \, dx &= \pi \, \delta_{mn} \\ \int_{-\pi}^{\pi} \sin(m x) \cos(n x) \, dx &= 0 \end{aligned} \]
That makes \(\mathbf{G}\) diagonal, and the normal equations decouple into one independent equation per coefficient, \(c_k = \langle \phi_k, f \rangle / \langle \phi_k, \phi_k \rangle\). For the Fourier series that is
\[ a_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos(n x) \, dx, \qquad b_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin(n x) \, dx . \]
No regression is required. Each coefficient is a projection of \(f\) onto one basis function, and it is computed with an integral. There is a second benefit hiding in that formula: \(a_n\) does not depend on \(N\). If we add more terms to the series, the coefficients we already have do not change. We just compute the new ones.
An example
Here is a lumpy function that is not obviously made of sines and cosines.
Top: f(x) and its Fourier partial sum. Bottom: the coefficients, each one computed from a single integral. Move the slider and the faded bars join the sum without any of the bars changing height.
With eight terms the series is indistinguishable from the function (the largest error is about 0.0015). It converges this quickly partly because \(f\) almost vanishes at \(\pm\pi\), so its periodic extension is smooth. A function that does not match up at the ends would show Gibbs ringing near them and converge much more slowly. Every coefficient came from a single integral, and adding a term just adds a bar to the chart; it never changes the ones already there.
Adding terms does not change the coefficients
To make the second point concrete, let’s fit the same function by least squares on a grid of points, first with 3 harmonics and then with 6. We do the same with a polynomial basis \(\{1, x, x^2, \ldots\}\), which is not orthogonal.
xg = np.linspace(-np.pi, np.pi, 2001)[:-1] # one period, evenly spaced
yg = f(xg)
def fourier_basis(x, N):
cols = [np.ones_like(x)]
for n in range(1, N + 1):
cols += [np.cos(n * x), np.sin(n * x)]
return np.column_stack(cols)
def poly_basis(x, N):
return np.column_stack([x**k for k in range(N + 1)])
for name, basis, Ns in [("Fourier", fourier_basis, (3, 6)),
("polynomial", poly_basis, (3, 6))]:
print(name)
for N in Ns:
c = np.linalg.lstsq(basis(xg, N), yg, rcond=None)[0]
print(f" N = {N}: first four coefficients", np.round(c[:4], 4))Fourier
N = 3: first four coefficients [ 0.2841 0.2015 0.1376 -0.2312]
N = 6: first four coefficients [ 0.2841 0.2015 0.1376 -0.2312]
polynomial
N = 3: first four coefficients [ 0.4682 0.1336 -0.056 -0.0188]
N = 6: first four coefficients [ 0.2641 0.3645 0.3232 -0.1281]The first four Fourier coefficients are the same whether we keep 3 harmonics or 6 (and they are the projections \(a_0/2\), \(a_1\), \(b_1\), \(a_2\) from before). The first four polynomial coefficients all change when we add \(x^4, x^5, x^6\), because those functions overlap with the ones already in the expansion. With a non-orthogonal basis, every time you add a function you have to solve for all the coefficients again.
One caveat: the grid matters. The sines and cosines are orthogonal as integrals, and they stay orthogonal as sums only on a uniform grid that covers exactly one period. That is why the last point is dropped with [:-1], since \(-\pi\) and \(\pi\) are the same point on the circle. On that grid the largest off-diagonal entry of \(\mathbf{H}^T\mathbf{H}\) is below \(10^{-12}\). On a non-uniform or non-periodic grid the discrete Fourier basis is not orthogonal, and the least squares coefficients would shift when you add terms too.
The Fourier series as a network
Now let’s draw the Fourier series the way we usually draw a neural network. There is one input, \(x\). It feeds a layer of “hidden units”, each of which computes \(\cos(n x)\) or \(\sin(n x)\). The outputs of those units are multiplied by the coefficients and summed to give \(f(x)\).
A truncated Fourier series drawn as a neural network with one hidden layer. The input weights are fixed integer frequencies; the output weights are the Fourier coefficients, drawn with a width proportional to their size.
It looks a lot like a neural network, and it really is one. Since \(\sin(nx) = \cos(nx - \pi/2)\), every hidden unit computes
\[ h_k(x) = \cos(w_k x + \beta_k), \]
with weight \(w_k = n\) and bias \(\beta_k \in \{0, -\pi/2\}\). The output is a linear combination \(f(x) = \sum_k c_k h_k(x)\). That is exactly a one-hidden-layer network with a cosine activation function. It is a special one in two ways:
- The first-layer weights and biases are fixed. We chose them ahead of time (integer frequencies, two phases) rather than learning them.
- Because of that choice the hidden units are orthogonal, so the output weights can be computed one at a time by projection, and they never change when we add more hidden units.
A typical neural network uses something like \(\tanh\) or ReLU instead of \(\cos\). The logistic sigmoid, and later \(\tanh\), came from the picture of a neuron firing: each is a smooth step. ReLU is a more modern choice, made because it makes deep networks easier to optimize. We also give up the integer frequencies of the Fourier basis. This gives us even more flexibility, but it comes at a cost: the hidden units are no longer orthogonal, and the first-layer weights now have to be fitted instead of computed. We have to find all the weights together by regression, and when we add more hidden units the old weights change.
There is a middle ground that makes the point even more clearly. In random-feature models and extreme learning machines, the first-layer weights are drawn at random and then frozen, and only the output weights are fitted. It is the same picture as the Fourier series, with a random basis in place of an orthogonal one. Because the basis is not orthogonal, the output weights come from least squares instead of projection.
Part 2: a feature generator for linear regression
Here is another way to read the same picture. A neural network is a nonlinear transformation of the input into a new, usually higher dimensional, vector space, followed by linear regression in that space. The hidden layers generate features. The last layer is plain old linear regression on those features.
We have actually already done both of these in Part 1. Let’s look at the Fourier network again, one step at a time, for a single input \(x\).
Step 1: expand the dimension. The first layer maps the scalar \(x\) to a vector with one element per hidden unit. For \(N = 3\) that is a vector in \(\mathbb{R}^7\):
\[ \mathbf{z} = \mathbf{w} x + \boldsymbol{\beta}, \qquad \mathbf{w} = \begin{bmatrix} 0 \\ 1 \\ 1 \\ 2 \\ 2 \\ 3 \\ 3 \end{bmatrix}, \quad \boldsymbol{\beta} = \begin{bmatrix} 0 \\ 0 \\ -\pi/2 \\ 0 \\ -\pi/2 \\ 0 \\ -\pi/2 \end{bmatrix}. \]
This step is linear. It does not add any information; it just makes copies of \(x\) that are stretched and shifted by different amounts. (The constant unit is \(\cos(0 \cdot x + 0) = 1\).)
Step 2: transform nonlinearly. Apply the activation function to each element:
\[ \mathbf{h} = \sigma(\mathbf{z}) = \cos(\mathbf{z}) = \begin{bmatrix} 1 & \cos x & \sin x & \cos 2x & \sin 2x & \cos 3x & \sin 3x \end{bmatrix}^T. \]
Now the seven elements are genuinely different functions of \(x\). These are the features.
Step 3: linear regression. The output is linear in the features,
\[ f(x) = \mathbf{c} \cdot \mathbf{h}(x). \]
If we have \(n\) data points, we stack the feature vectors as rows of a matrix \(\mathbf{H}\) (\(n \times 7\)), and finding \(\mathbf{c}\) is the linear least squares problem \(\min_\mathbf{c} \| \mathbf{H}\mathbf{c} - \mathbf{y} \|^2\), with solution \(\mathbf{c} = (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{y}\). That is the problem np.linalg.lstsq solved with fourier_basis in Part 1. It uses an SVD rather than forming that inverse, but the answer is the same. In the Fourier case \(\mathbf{H}^T\mathbf{H}\) happens to be diagonal, which is why the integrals were enough.
Here are the three steps written as a network layer, with matrices instead of loops.
w = np.array([0, 1, 1, 2, 2, 3, 3])
beta = np.array([0, 0, -np.pi / 2, 0, -np.pi / 2, 0, -np.pi / 2])
Z = xg[:, None] * w + beta # step 1: (n, 1) -> (n, 7)
H = np.cos(Z) # step 2: nonlinear activation
c = np.linalg.lstsq(H, yg, rcond=None)[0] # step 3: linear regression
print(H.shape)
print(np.round(c, 4))(2000, 7)
[ 0.2841 0.2015 0.1376 -0.2312 0.2014 -0.1485 0.1126]These are the same coefficients as before.
The output is linear in the new space
The point of the transformation is that a function which is nonlinear in \(x\) becomes linear in \(\mathbf{h}\). It is the same trick as fitting \(y = a + b/T\): the model is nonlinear in \(T\), but it is linear in the feature \(1/T\), so we compute that feature and fit a line. (Linearizing the Arrhenius equation, \(\ln k = \ln k_0 - \frac{E_a}{R}\frac{1}{T}\), is a close cousin. There we transform the output to \(\ln k\) as well, while a network only transforms the input.)
We can see this geometrically with only one harmonic, so the feature space is small enough to draw. Take \(y = 0.5 + 0.8\cos x - 0.6 \sin x\). As a function of \(x\) it is a wave. In the feature space \((h_1, h_2) = (\cos x, \sin x)\), the input interval is mapped onto a circle, and the data lie exactly on a flat plane \(y = c_0 + c_1 h_1 + c_2 h_2\) above that circle.
Left: the data are nonlinear in x. Right: the same data in the feature space (cos x, sin x) lie on a plane, so finding the coefficients is linear regression. Hover a point to find it in both panels.
The grey circle on the floor is where the inputs land in feature space. The points sit on the plane, and fitting a plane is linear regression.
Deep networks
So far there is one hidden layer. What changes when we stack more of them? Let’s write a network with two hidden layers using matrices, for a batch of \(n\) inputs \(\mathbf{X}\) (\(n \times d\)):
\[ \begin{aligned} \mathbf{H}_1 &= \sigma(\mathbf{X}\mathbf{W}_1 + \mathbf{b}_1) && (n \times d) \cdot (d \times m_1) \to (n \times m_1) \\ \mathbf{H}_2 &= \sigma(\mathbf{H}_1\mathbf{W}_2 + \mathbf{b}_2) && (n \times m_1) \cdot (m_1 \times m_2) \to (n \times m_2) \\ \mathbf{y} &= \mathbf{H}_2 \mathbf{c} + c_0 && (n \times m_2) \cdot (m_2 \times 1) \to (n \times 1) \end{aligned} \]
Each hidden layer is a linear map to a new dimension followed by an elementwise nonlinearity, just like Step 1 and Step 2 above. Stacking them is function composition. Call the layer maps \(\mathbf{g}_1(\mathbf{x}) = \sigma(\mathbf{x}\mathbf{W}_1 + \mathbf{b}_1)\) and \(\mathbf{g}_2(\mathbf{u}) = \sigma(\mathbf{u}\mathbf{W}_2 + \mathbf{b}_2)\). The features at the end are
\[ \mathbf{h}_2(\mathbf{x}) = \sigma\big(\sigma(\mathbf{x}\mathbf{W}_1 + \mathbf{b}_1)\mathbf{W}_2 + \mathbf{b}_2\big) = (\mathbf{g}_2 \circ \mathbf{g}_1)(\mathbf{x}), \]
and the output is still linear in them,
\[ f(\mathbf{x}) = \mathbf{c} \cdot \mathbf{h}_2(\mathbf{x}) + c_0 . \]
However many layers there are, the last one is a linear model on the features, and with a squared-error loss, fitting it is linear regression. The columns of \(\mathbf{H}_2\) are \(m_2\) basis functions, each one a composition of the layers before it. The only difference from the Fourier series is that these basis functions are not chosen in advance; their shapes are set by \(\mathbf{W}_1, \mathbf{b}_1, \mathbf{W}_2, \mathbf{b}_2\), which are fitted along with \(\mathbf{c}\).
A network with two hidden layers in graph form (top) and matrix form (bottom). The hidden layers transform n inputs into n rows of 3 features, the columns of H₂. The output layer is linear regression on them. (Three features keep the drawing small; the trained example below uses six.)
Checking that the last layer is linear regression
Let’s train a small network on the function from Part 1 and look inside it. This one has two hidden layers of six \(\tanh\) units. To keep it to numpy and scipy, we pack all the weights into one vector and minimize the mean squared error with BFGS.
from scipy.optimize import minimize
sizes = [1, 6, 6] # input, hidden layer 1, hidden layer 2
shapes = []
for m, k in zip(sizes[:-1], sizes[1:]):
shapes += [(m, k), (k,)] # W, b for each hidden layer
shapes += [(sizes[-1],), ()] # c, c0 for the output layer
def unpack(p):
out, i = [], 0
for s in shapes:
k = int(np.prod(s))
out.append(p[i:i + k].reshape(s))
i += k
return out
def features(p, X):
"The hidden layers: H = tanh(... tanh(X W1 + b1) ... W + b)."
*hidden, c, c0 = unpack(p)
H = X
for W, b in zip(hidden[::2], hidden[1::2]):
H = np.tanh(H @ W + b)
return H
def nn(p, X):
*_, c, c0 = unpack(p)
return features(p, X) @ c + c0
xt = np.linspace(-np.pi, np.pi, 200)
Xt, yt = xt[:, None], f(xt)
def loss(p):
return np.mean((nn(p, Xt) - yt)**2)
p0 = np.random.default_rng(42).normal(size=sum(int(np.prod(s)) for s in shapes))
sol = minimize(loss, p0, method="BFGS", options={"maxiter": 5000})
print(f"{len(p0)} parameters, MSE = {sol.fun:.2e}")61 parameters, MSE = 4.81e-06Now take the trained hidden layers, compute the features \(\mathbf{H}_2\) for the training inputs, and do ordinary linear regression of \(y\) on them. If the last layer is linear regression, we should get back the output weights the network learned.
H2 = features(sol.x, Xt)
A = np.column_stack([H2, np.ones(len(xt))])
c_lstsq = np.linalg.lstsq(A, yt, rcond=None)[0]
*_, c_nn, c0_nn = unpack(sol.x)
print("network:", np.round(np.append(c_nn, c0_nn), 3))
print("lstsq: ", np.round(c_lstsq, 3))
print(f"cond(A) = {np.linalg.cond(A):.0f}")network: [-1.044 0.962 -1.522 1.387 -0.859 0.774 -0.994]
lstsq: [-1.044 0.962 -1.523 1.387 -0.859 0.776 -0.994]
cond(A) = 159If you run this yourself, expect different numbers. The path BFGS takes depends on the numpy, scipy and BLAS versions, so the trained weights will differ. What should reproduce is that the two rows agree.
They agree to within the tolerance the optimizer stopped at, and that is not a coincidence. Write \(\mathbf{A} = [\mathbf{H}_2 \;\; \mathbf{1}]\) and \(\boldsymbol{\theta} = (\mathbf{c}, c_0)\). With the hidden weights held fixed, the loss \(\frac{1}{n}\|\mathbf{A}\boldsymbol{\theta} - \mathbf{y}\|^2\) is a convex quadratic in \(\boldsymbol{\theta}\). Its gradient, \(\frac{2}{n}\mathbf{A}^T(\mathbf{A}\boldsymbol{\theta} - \mathbf{y})\), vanishes at a stationary point of the full training problem, and that condition is exactly the normal equations. As long as \(\mathbf{A}\) has full rank (its condition number is modest here), the network’s output weights are the least squares solution on its own features.
Finally, here are the six basis functions the network built, the columns of \(\mathbf{H}_2\). Compare them to the sines and cosines in Part 1: these are not orthogonal and they are not pretty, but a weighted sum of them reproduces \(f(x)\).
Left: the columns of H₂, the basis functions learned by the hidden layers. Switch to cₖ·hₖ to see them scaled by their output weights. Right: their linear combination. Toggle basis functions off to see that the fit needs all of them.
So a neural network does the same two things as the Fourier series: it transforms the input into a new vector space with a nonlinear function, and it does linear regression in that space. The Fourier series fixes the transformation in advance and gets the regression for free. A neural network learns the transformation, and pays for it with a nonlinear optimization over all the weights at once.
You can see from the learned basis functions that there is a lot of cancellation: each value of the output is a sum of positive and negative contributions from the hidden units. That is a common property of learned features. They are not orthogonal, and they are rarely interpretable on their own; individually they mean little, and only their linear combination reproduces the output.
We would never ask “What did the Fourier series learn about our data?” The Fourier basis is fixed, and all the data-dependence is in the coefficients. A neural network is different in exactly one place: the features are learned too, so that is where most of what it “learned” lives. That is also why transfer learning works, since features learned on one problem are often a good basis for another. So these viewpoints do not make a network fully transparent, but they do make its structure transparent. It is a learned nonlinear feature transform followed by a linear model. The structure is simple; the learned features usually are not.
We used a simple dense network with a squared-error loss here, which makes the last layer linear regression. The same reading applies to other heads with a small change. A classifier with a softmax output is doing multinomial logistic regression on its features, so in general the last layer is a generalized linear model. It also applies to convolutional networks, recurrent networks, transformers, graph neural networks and many other architectures. What differs is how the features are generated.