Neural networks and Reservoir Computing¶

Introduction¶

In this lecture we will cover the following topics:

Windowed approaches and nonlinear models for time series forecasting.
Neural networks and the Multi-Layer Perceptron.
A brief overview of Recurrent Neural Networks.
The ESN, a randomized RNN from the family of Reservoir Computing.
An introduction to dimensionality reduction with Principal Component Analysis.
Examples of forecasting with MLP and ESN on real-world electricity data.

In [39]:

import sys

# Install dependencies if the notebook is running in Colab
if 'google.colab' in sys.modules:
    !pip install -U -qq reservoir_computing

Windowed approaches for prediction¶

Limitation of linear models¶

Linear models are widely used for their simplicity and interpretiveness.
However, they assume a constant relationship between input and target variables.
As such, they cannot capture these nonlinearities and interactions effectively.

Consider the following examples:
1. Compute the position $y$ of a robotic arm given the positions $x_1, x_2, \dots$ of the joints.
2. Compute a trajectory $\dot{y}$ given the angular velocities $\omega_1, \omega_2, \dots$ of the joints.

A linear model will fail in the second task.

Multi-Layer Perceptron¶

A Multi-Layer Perceptron (MLP) is a simple Neural Network that consists of at least three layers of nodes:
- An input layer.
- One or more hidden layers.
- An output layer.
The nodes in the a hidden layers apply a nonlinear activation function $\sigma$ .
MLP are generally trained with a supervised learing technique called backpropagation.

The MLP's layers are fully connected, meaning each node in one layer connects with a certain weight to every node in the following layer.
The first layer receives the input.
The output of each layer is the input for the next layer until the final layer produces the output of the MLP.
The weights are stored in matrices $W$ and are the trainable parameters of the model.
The MLP can learn complex mappings from inputs to outputs to perform a wide range of data modeling and prediction tasks.

Example: a simple MLP¶

Let's consider a concrete example of an MLP with one input layer, two hidden layers (with $h_1$ and $h_2$ units respectively), and one output layer.
The input vectors are of size $d_\text{in}$ , and the output vectors are of size $d_\text{out}$ .

Input layer

The size of the input layer corresponds to the size of the input vectors, $d_\text{in}$ .
This layer simply passes the input to the first hidden layer without applying any transformation.

First hidden layer

This layer has $h_1$ nodes.
Each node connects to all the $d_{\text{in}}$ inputs.
It applies a weight matrix $W_\text{in} \in \mathbb{R}^{d_\text{in} \times h_1}$ that transforms the input from vectors in $\mathbb{R}^{d_\text{in}}$ to vectors in $\mathbb{R}^{h_1}$ .
Then, it applies a nonlinear activation function $\sigma$ , and outputs $h_1$ intermediate results.

Second hidden layer

This layer has $h_2$ nodes.
It takes the $h_1$ outputs from the first hidden layer.
It applies another set of weights and a nonlinear activation function to produce $h_2$ intermediate results.

Output Layer

The output layer has $d_\text{out}$ units.
It takes the $h_2$ outputs from the second hidden layer, applies weights and possibly a different activation function that depends on the task (usually, softmax for classification or no activation for regression).
The final output vector is of size $d_\text{out}$ .

MLP for time series forecasting¶

The MLP can be used as a windowed technique to perform time series forecasting.
The sequence of past values $x(t-P), \dots, x(t-1), x(t)$ will represent our input vector ( $d_\text{in} = P$ ).
The future value $x(t+\tau)$ will be the output that the model will learn to predict ( $d_\text{out} = 1$ ).
We can also train the MLP to predict a sequence of $d_\text{out} = H$ future predictions, e.g., $x(t+\tau), x(t+\tau+1), \dots, x(t+\tau+H)$ .

Key for the training of the MLP is the creation of input-output pairs $\{\mathbf{x}_i, \mathbf{y}_i\}$ used for training.
Given the input $\mathbf{x}_i$ the model must learn to predict the output $\mathbf{y}_i$ .
This is done by adapting the model weights to minimize the discrepancy between the prediction $\mathbf{\hat y}_i$ and the desired output $\mathbf{y}_i$ .
The weights are modified according to the gradient taken with respect to a loss function $\mathcal{L}(\mathbf{\hat y}_i, \mathbf{y}_i)$ such as the MSE, e.g., $W \leftarrow W + \delta \frac{\partial \mathcal{L}}{\partial W}$ .
$\delta$ is a small constant defining the gradient step (learning rate).

Next, we define the function to create input-output pairs.
For example, we create as input windows of size $P=12$ by specifying window_size=12.
Similarly, we create as ouput windows of size forecast_horizon=12.
If we are interested in predicting only one sample, e.g., $\tau=12$ and $H=1$ , we simply take the last element of the output window.

In [6]:

def create_windows(data, window_size, forecast_horizon):
    X, y = [], []
    for i in range(len(data) - window_size - forecast_horizon + 1):
        X.append(data[i:(i + window_size)])
        y.append(data[i + window_size:i + window_size + forecast_horizon])
    return np.array(X), np.array(y)

# Define window size and forecast horizon
window_size = 12
forecast_horizon = 12

# Create input-output pairs
X_train, y_train = create_windows(tr, window_size, forecast_horizon)
X_test, y_test = create_windows(te, window_size, forecast_horizon)

Next, we define the MLP.
We specify the size of the hidden layers as $h_1=16$ and $h_2=8$ . Using larger values increases the model capacity, but can lead to overfit.
As the activation function we use a ReLU.
As the algorithms to compute the gradients and update the values of the parameters $W$ we use Adam.
Finally, we train the model for 1000 iterations.

In [10]:

# Define and train the neural network
mlp = MLPRegressor(hidden_layer_sizes=(16,8), activation='relu', solver='adam', max_iter=1000)
mlp.fit(X_train_scaled, y_train)

Out[10]:

MLPRegressor(hidden_layer_sizes=(16, 8), max_iter=1000)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Is the prediction good or not?
To get an idea, we can compare the performance against a simple baseline.
In this case, we are making a prediction at a forecast horizon $\tau=24$ equal to the main seasonality of the time series.
The most natural baseline is to use the values of the previous day as the prediction for the next day.

In [13]:

# Compute the mse between the original time series and the time series shifhted by 24 time steps
mse = mean_squared_error(y_test[24:], y_test[:-24])
print(f"Mean Squared Error: {mse:.2f}")

plt.figure(figsize=(7, 3))
plt.plot(y_test[24:])
plt.plot(y_test[:-24])
plt.grid()
plt.show()

Mean Squared Error: 27.14

Any reasonable forecasting method must beat this baseline.

Windowed approaches drawbacks¶

What is the optimal value of $P$ ?
If $P$ is fixed, how can we deal with different temporal dependencies?
Let's consider the example of predicting the next word in a sentence.

RC is a familily of randomized RNNs, popularized in machine learning by Echo State Networks (ESNs).
RC and ESNs are terms ofen used interchangeably.
There are two main differences that separate an ESN from an RNN:

The output weights $W_o$ is the only part of the ESN that is trained.
Optimizing $W_o$ can be done with a simple linear regression algorithm.
The workflow is the following.

Generate a sequence of reservoir states $\mathbf{H} = \{\mathbf{h}(0), \mathbf{h}(1), \dots, \mathbf{h}(T)\}$ .
- This is done by applying the state update equation
ht=σ(Wix(t)+Whh(t−1))
- for each time step $t=1,2, \dots, T$ of the input sequence.
- The nonlinearity $\sigma$ is usually an hyperbolic tangent ( $\texttt{tanh}$ ).

Apply a linear regression algorithm to compute a linear mapping $g(\cdot)$ between the reservoir states and the desired output sequence $\mathbf{y} = \{y(1), y(2), \dots, y(T) \}$ :

$\mathbf{y} = g(\mathbf{H})$

The function $g(\cdot)$ is called readout.
In a forecasting setting, $y(t)$ correspond to a future value of the input, e.g., $y(t) = x(t+\tau)$ .

Readout¶

The readout $g(\cdot)$ is usually implemented through a linear regression model, e.g., Ridge Regression.
In this case $g(\cdot)$ corresponds to a weight matrix $\mathbf{W}_o$ .
However, any other regression model can be used to implement the readout, including an MLP.

Why the Reservoir approach works¶

The Reservoir extracts a rich pool dynamical features from the time series.
These are embedded into the high-dimensional Reservoir state $\mathbf{h}(t)$ .
Contrary to a fixed window, $\mathbf{h}(t)$ maintains a memory of all the previous inputs, back to the origin of the series $x(0)$ .
Some of the features are relevant for the task at hand, while others are not.
The task of the readout is to select those features relevant for the task.

Spectral radius¶

Cleary, the contribution of more recent inputs must have a stronger influence on the current state.
This means that the Reservoir should gradually forget its past states.
This property, called echo state property, gives the name to the Echo State Network.
It ensures to not model noise, to forget sporadic shocks, and the initial state $\mathbf{h}(0)$ that is uninformative.
In control theory, this translate in having a dynamic that is contractive.
A dynamic is contractive if two initially different states eventually converge.

Recall the state update equation

$\mathbf{h}(t) = \sigma \left(\mathbf{W}_i x(t)+\mathbf{W}_h \mathbf{h}(t-1)\right)$

The spectral radius is the largest eigenvalue of the state transtion matrix $W_h$ .
We can set the spectral radius by computing the largest eigenvalue $\lambda_\text{max}$ of $W_h$ and then letting $W_h = \rho \frac{W_h}{\lambda_\text{max}}$ .

A rule of thumb is to set $\rho$ just below 1.
However, to achieve good performance it is often necessary to fine-tune $\rho$ to values that can be lower or even higher than 1.

Let's check this with a practical example.
We consider two intializations, $\mathbf{h}_1(0) = [0,0,\dots,0]$ and $\mathbf{h}_2(0) = [1,1,\dots, 1]$ .
For this example, we set the input $x$ to be always zero to not let it affect the evolution of the Reservoir state.

In [14]:

# Initial states
initial_state_0 = np.zeros((1, 100), dtype=float) 
initial_state_1 = np.ones((1, 100), dtype=float)

x = np.zeros((1, 100, 1)) # Zero input, it does not contribute to the state
rhos = [.3, 0.99, 1.3]    # We will use three different spectral radii

In [16]:

plot_states_evolution(x, rhos, initial_state_0, initial_state_1)

Reservoir units are usually equipped with a $\texttt{tanh}$ activation.
A small value $\omega_\text{in}$ maps the Reservoir inputs towards the centre of the $\texttt{tanh}$ where is more linear.
Therefore, a small $\omega_\text{in}$ reduces the amount of nonlinearity in the Reservoir.
On the other hand, a large $\omega_\text{in}$ translates into a more nonlinear beahvior as the $\texttt{tanh}$ is closer to saturation.

Like for the spectral radius, $\omega_\text{in}$ should be tuned carefully.
A good starting value is usually around 0.1.

Another hypeparameter is the number of Reservoir units $N_h$ .
A larger value can give better performance at the cost of higher computation time.
A good starting point is usually $N_h=300$ , to be increased until there is no more gain in performance.

There are also other hyperparameters in the ESN, such as the sparsity of the Readout and an optional noise to inject in the state update equation.
These are usually less critical than $\rho$ and $\omega_\text{in}$ and can be left to their default value in most cases.
Tuning hyperparameters in randomized architectures such as the ESN is much more important than in trainable neural networks, since there is no training that can compensate for poorly initialized models.

Then, we compute the sequence of the Reservoir states states_tr and states_te associated with the training and test data, respecitvely.
Since the initial states of the Reservoir mostly depend on the initialization, e.g., $\mathbf{h}(0) = [0,0,\dots,0]$ , we want to get rid of the initial transient phase.
This can be done by setting n_drop that specifies how many of the initial states we drop.

In [21]:

n_drop=10
states_tr = res.get_states(Xtr[None,:,:], n_drop=n_drop, bidir=False)
states_te = res.get_states(Xte[None,:,:], n_drop=n_drop, bidir=False)
print(f"states_tr shape: {states_tr.shape}\nstates_te shape: {states_te.shape}")

states_tr shape: (1, 2666, 900)
states_te shape: (1, 266, 900)

Fit a linear readout¶

We are now ready to train the readout to predict the desired output given the sequence of Reservoir states.
We start by using a linear readout implemented by a Ridge regressor.
After the readout is trained, we use it to compute the predictions $\hat{Y}_\text{te}$ of the test data.

In [22]:

# Fit the ridge regression model
ridge = Ridge(alpha=1.0) 
time_start = time.time()
ridge.fit(states_tr[0], Ytr[n_drop:,:])
print(f"Training time: {time.time()-time_start:.4f}s")

# Compute the predictions
time_start = time.time()
Yhat = ridge.predict(states_te[0])
print(f"Test time: {time.time()-time_start:.4f}s")

Training time: 0.0344s
Test time: 0.0003s

Our data is a sequence of length $T$ of Reservoir states, each one of size $N_h$ .
They can be arranged in a matrix $\mathbf{H} \in \mathbb{R}^{T \times N_h}$ .
To illustrate the procedure, let's first consider a toy example with only $N_h = 3$ features.
In addition, we will create some structure in the data, by dividing the samples in 4 groups/clusters.
This will help us to see how PCA preserves the structure in the data.

In [25]:

# Generate 4 clusters of points in 3 dimensions
T = 300 # number of samples
N_h = 3 # number of features
H, clust_id = make_blobs(n_samples=T, n_features=N_h, 
                         centers=4, cluster_std=1.5, random_state=1)

Reducing the data dimensionality with PCA invole the following steps:

Standardization

Scale the data so that each feature has a mean of 0 and a standard deviation of 1.
This is important because PCA is affected by scale.

Covariance matrix computation:

Calculate the empirical covariance matrix $\mathbf{H}^T\mathbf{H}$ .
The matrix shows how changes in one variable are associated with changes in another variable.

Eigenvalues and eigenvectors computation

The eigenvectors of the covariance matrix represent the directions of maximum variance.
In the context of PCA, the eigenvectors are the principal components.
The eigenvalues indicate the variance explained by each principal component.

Sorting eigenvectors:

The eigenvectors are sorted by decreasing eigenvalues.
The top- $k$ eigenvectors are selected, where $k$ is the number of dimensions we want to keep.
In our case, we keep $k=2$ dimensions.

Projection onto the new feature space:

The original data are projected onto the selected principal components.
In our case, the 3-dimensional data are projected onto the plane spanned by the first two principale components.
The projected data are the reduced-dimensional data.

Exercises¶

Exercise 1¶

In the MLP example, change the window_size ( $P$ ) from 12 to 24.
Comment on the results and explain why the performance improve.

Exercise 2¶

In the MLP, change forecast_horizon from 24 to 36.
Do the same change also in the ESN.
With a forecast horizon of size 24 we used as a baseline the MSE between the time series and its shifted version by 24 lags, mse = mean_squared_error(y_test[24:], y_test[:-24]). Is this a good baseline also when the forecast horizon changes? What could be a better baseline?

Exercise 3¶

Contrarily from model-based approach such as ARIMA, the MLP and the ESN do not strctly require the data to be stationary. However, by removing trend and seasonality the MLP and the ESN can focus only on predicting the residuals.

Remove trend and seasonality and train both MLP and ESN to train the residuals.
In testing, remember to put back the trend and seasonality when computing the predictions.
Compare the results with those obtained by training the MLP and the ESN on the original data.

Exercise 4¶

For this exercise, we will consider a new dataset consisting of hourly water temperatures from Gulf of Mexico near Key West, Florida. The hourly temperatures are provided from October 3, 2016 to October 3, 2017 The dataset consists of a data frame with 6572 observations on the following 3 variables.

DateTime Date and time of reading (format mm/dd/yyyy h:00).
WaterTemp: Water temperature (in degrees Fahrenheit).
t: Time index (1 to 673).

In [38]:

temps = sm.datasets.get_rdataset("KeyWestWater", "Stat2Data").data
print("Raw data:\n------------------\n", temps.head())
plt.figure(figsize=(14,4))
plt.plot(temps["WaterTemp"])
plt.grid()
plt.show()

Raw data:
------------------
          DateTime  WaterTemp  t
0  10/3/2016 0:00       86.2  1
1  10/3/2016 1:00       86.2  2
2  10/3/2016 2:00       86.2  3
3  10/3/2016 3:00       86.2  4
4  10/3/2016 4:00       86.0  5

Apply the MLP on the temperature data.
Perform a small grid search to try different values for the following hyperparameters:
- hidden_layer_sizes (e.g., try more layers [16, 16, 8], layers with more units [32, 16], etc...).
- learning_rate_init.
- activation (try at least tanh and relu).
Apply the ESN on the temperature data.
Perform a small grid search to try different values for the following hyperparameters:
- spectral_radius.
- input_scaling.
- n_internal_units.
- connectivity.

Report the configuration of MLP and ESN that achieves the best performance on the test data.

Neural networks and Reservoir Computing¶

Introduction¶

Windowed approaches for prediction¶

Limitation of linear models¶

Example: forecasting electricity demand¶

Advantages of Nonlinear Models¶

Neural networks¶

Multi-Layer Perceptron¶

Example: a simple MLP¶

MLP for time series forecasting¶

Python example¶

Windowed approaches drawbacks¶

Recurrent Neural Networks¶

Challenges and shortcomings in RNNs¶

Reservoir Computing¶

Readout¶

Why the Reservoir approach works¶

Reservoir configuration¶

Spectral radius¶

Input scaling¶

Python example: electricity load forecasting with ESN¶

Load the data¶

Define the Reservoir and compute the states¶

Fit a linear readout¶

Reducing the dimensionality of the Reservoir states¶

Principal Component Analysis (PCA)¶