Examples
There is no need to write networks by hand. In torch a network with a three-dimensional input, 30 hidden relu neurons and a two-dimensional output is
import numpy as np
import torch
import matplotlib.pyplot as plt
torch.manual_seed(7)
net = torch.nn.Sequential(
torch.nn.Linear(3, 30),
torch.nn.ReLU(),
torch.nn.Linear(30, 2),
)
print(net[0].weight.shape, net[0].bias.shape)torch.Size([30, 3]) torch.Size([30])The weight matrix of the first layer has shape \(30 \times 3\), and there are 30 biases. The network computes what we expect it to compute: the matrix form of multilayer perceptrons.
x = torch.rand(3)
W1, b1, W2, b2 = net[0].weight, net[0].bias, net[2].weight, net[2].bias
print(net(x))
print(b2 + W2 @ torch.relu(b1 + W1 @ x))torch networks compute with 32-bit floating point numbers by default, and so must their inputs; torch.rand makes such numbers.
tensor([0.1888, 0.0847], grad_fn=<ViewBackward0>)
tensor([0.1888, 0.0847], grad_fn=<AddBackward0>)Networks can be arbitrarily deep and wide, with different activation functions in different layers, for example 20 inputs, hidden layers of 100, 50 and 100 neurons and 10 outputs.
deep = torch.nn.Sequential(
torch.nn.Linear(20, 100), torch.nn.ReLU(),
torch.nn.Linear(100, 50), torch.nn.Tanh(),
torch.nn.Linear(50, 100), torch.nn.Sigmoid(),
torch.nn.Linear(100, 10),
)
print(deep(torch.rand(20)).shape)torch.Size([10])Every hidden neuron adds one kink, and with more neurons the fit follows the data more closely. A relu neuron \(\max(0, wx + b)\) is a hinge like those of the splines, with its knot at \(-b/w\): in one dimension a relu network is a linear spline whose knots are found by gradient descent. The spline with fixed knots is fitted by least squares, which cannot get stuck, and up to 16 knots it has a test error as low as or lower than the network with as many neurons, with fewer parameters. With 40 knots on 50 points it overfits.
With \(p\) inputs a polynomial or a spline needs a column for every product of functions of the inputs: a polynomial of degree \(d\) with all interactions has \(\binom{p+d}{d}\) columns. A generalized additive model (GAM) keeps one spline per input and sums them, SplineTransformer followed by LinearRegression, so its columns grow only linearly with \(p\), but it cannot fit interactions. A relu neuron \(\max(0, w^\top x + b)\) is a hinge along a direction \(w\) that it learns, and its output does not change in the directions perpendicular to \(w\).
With one input all four fit well, and the polynomial and the GAM best. The function is not a sum of functions of single inputs, so from two inputs on the GAM misses most of it. kNN gets steadily worse, because with more inputs the nearest neighbours are farther away: with 50 inputs its test error, 0.44, is close to the variance of the function, 0.50. With 8 inputs the polynomial has 1287 columns for 2000 points and a test error of 0.10, and from 10 inputs on it has more columns than there are points. The network needs \(8(p + 2) + 1\) parameters, 417 for 50 inputs, and its test error grows slowly, to 0.035. The function varies along one direction only, and the network finds it.
Can a network predict the wind peak five hours ahead better than the linear regression of week 2? We standardize the input and the output, because the usual initializations of networks work best on standardized data, and fit a network with one hidden layer of 128 relu neurons by AdamW on batches of 128.
A network this large fits the training years better and better the longer it is trained, and the test years worse after a while. So we fit on 2015 and 2016, keep 2017 as a validation set, and stop at the epoch with the lowest validation error, as in early stopping. The test years 2019 and 2020 are used only once, at the end.
import copy
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
weather = pd.read_csv("https://go.epfl.ch/bio322-weather2015-2018.csv")
weather_test = pd.read_csv("https://go.epfl.ch/bio322-weather2019-2020.csv")
drop = ["LUZ_wind_peak", "time"]
X = weather.iloc[:-5].drop(columns=drop).values
y = weather["LUZ_wind_peak"].values[5:] # the wind peak five hours later
X_test = weather_test.iloc[:-5].drop(columns=drop).values
y_test = weather_test["LUZ_wind_peak"].values[5:]
year = weather["time"].values[:-5] // 1_000_000
fit, val = year < 2017, year == 2017
x_scaler = StandardScaler().fit(X[fit])
y_scaler = StandardScaler().fit(y[fit].reshape(-1, 1))
def to_tensor(A):
return torch.tensor(x_scaler.transform(A), dtype=torch.float32)
X_fit, X_val, X_te = to_tensor(X[fit]), to_tensor(X[val]), to_tensor(X_test)
y_fit = torch.tensor(y_scaler.transform(y[fit].reshape(-1, 1)), dtype=torch.float32)
def rmse(a, b):
return float(np.sqrt(np.mean((a - b) ** 2)))torch.manual_seed(31)
model = torch.nn.Sequential(
torch.nn.Linear(X_fit.shape[1], 128), torch.nn.ReLU(),
torch.nn.Linear(128, 1),
)
opt = torch.optim.AdamW(model.parameters())
loader = torch.utils.data.DataLoader(
torch.utils.data.TensorDataset(X_fit, y_fit), batch_size=128, shuffle=True)
def predict(A):
with torch.no_grad():
return y_scaler.inverse_transform(model(A).numpy())[:, 0]
curve, best = [], (np.inf, None, 0)
for epoch in range(1, 41):
for xb, yb in loader:
opt.zero_grad()
torch.mean((model(xb) - yb) ** 2).backward()
opt.step()
curve.append((rmse(predict(X_fit), y[fit]), rmse(predict(X_val), y[val])))
if curve[-1][1] < best[0]:
best = (curve[-1][1], copy.deepcopy(model.state_dict()), epoch)
model.load_state_dict(best[1])
linear = LinearRegression().fit(X[fit], y[fit])
print(f"stopped at epoch {best[2]}")
print(f"network: training rmse {rmse(predict(X_fit), y[fit]):.2f}, "
f"test rmse {rmse(predict(X_te), y_test):.2f}")
print(f"linear: training rmse {rmse(linear.predict(X[fit]), y[fit]):.2f}, "
f"test rmse {rmse(linear.predict(X_test), y_test):.2f}")stopped at epoch 5
network: training rmse 7.34, test rmse 8.72
linear: training rmse 7.84, test rmse 8.90curve = np.array(curve)
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(range(1, 41), curve[:, 0], label="training (2015-2016)")
ax.plot(range(1, 41), curve[:, 1], label="validation (2017)")
ax.axvline(best[2], color="#7a838b", ls="--", lw=1)
ax.set(xlabel="epoch", ylabel="rmse [km/h]")
ax.legend()
plt.show()
The network predicts the test years better than the linear regression, but a lot of the variability is not captured.
fig, ax = plt.subplots(figsize=(4.8, 4.2))
ax.plot(predict(X_te), y_test, ".", ms=2, alpha=0.3)
ax.plot([0, 80], [0, 80], color="#a8452f", lw=1.5)
ax.set(xlabel="predicted wind peak [km/h]", ylabel="measured wind peak [km/h]",
xlim=(0, 80), ylim=(0, 80))
plt.show()
If the model predicted the test data perfectly, all points would lie on the diagonal. They scatter widely around it, most likely because the irreducible error of this problem is large.
A network can output more than a mean. With two outputs it can give the mean \(m(x)\) and the variance \(s(x)\) of a normal distribution for \(Y\,|\,x\), both non-linear functions of the input. The loss is the negative log-likelihood of the normal distribution, without its constant, \[ \mathcal L = \frac1n\sum_{i=1}^n \left(\frac{(y_i - m(x_i))^2}{2s(x_i)} + \frac12\log s(x_i)\right), \] and the second output is passed through \(\mathrm{softplus}(z) = \log(1 + e^z)\), so that the variance is positive. We fit the logarithm of the wind peak five hours later against the pressure in Luzern, so that the fitted distribution of the wind peak itself is log-normal.
pressure = weather["LUZ_pressure"].values[:-5].reshape(-1, 1)
log_peak = np.log(weather["LUZ_wind_peak"].values[5:]).reshape(-1, 1)
px, py = StandardScaler().fit(pressure), StandardScaler().fit(log_peak)
Xp = torch.tensor(px.transform(pressure), dtype=torch.float32)
yp = torch.tensor(py.transform(log_peak), dtype=torch.float32)[:, 0]
torch.manual_seed(1)
net2 = torch.nn.Sequential(torch.nn.Linear(1, 50), torch.nn.Tanh(),
torch.nn.Linear(50, 2))
def normal_loss(out, y):
m, s = out[:, 0], torch.nn.functional.softplus(out[:, 1])
return torch.mean((m - y) ** 2 / (2 * s) + 0.5 * torch.log(s))
opt = torch.optim.AdamW(net2.parameters())
loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(Xp, yp),
batch_size=32, shuffle=True)
for epoch in range(20):
for xb, yb in loader:
opt.zero_grad()
normal_loss(net2(xb), yb).backward()
opt.step()grid = np.linspace(pressure.min(), pressure.max(), 200).reshape(-1, 1)
with torch.no_grad():
out = net2(torch.tensor(px.transform(grid), dtype=torch.float32)).numpy()
m = out[:, 0]
sd = np.sqrt(np.log1p(np.exp(out[:, 1])))
back = lambda z: np.exp(py.inverse_transform(z.reshape(-1, 1))[:, 0])
fig, ax = plt.subplots(figsize=(8.8, 3.6))
ax.plot(pressure[:, 0], np.exp(log_peak[:, 0]), "o", ms=4.5, mew=0.7, alpha=0.4,
mfc="none", mec="#7a838b", rasterized=True) # open grey circles, for both themes
ax.plot(grid, back(m), color="#a8452f", lw=2.5, label="mean")
ax.plot(grid, back(m + sd), color="#a8452f", lw=2, ls="--", label="± one standard deviation")
ax.plot(grid, back(m - sd), color="#a8452f", lw=2, ls="--")
ax.set(xlabel="LUZ_pressure [hPa]", ylabel="LUZ_wind_peak in 5 h [km/h]", ylim=(0, 115))
ax.legend()
plt.show()
The mean depends slightly non-linearly on the pressure, and the spread of the wind peak is smaller at high pressure than at low pressure, where there are few points and several of them have high wind peaks. The dashed lines are not symmetric around the mean, because the fitted distribution is log-normal.