Sinusodial ESN example from Herbert Jäger

  • Apply RC on Sine Functions
  • Input function: $u(n)_{teach} = sin(n/5)$. Target: $y(n)_{teach} = sin^7(n/5)/2$

Information on notation

Recipe

  1. Let the network run for n = 0 to nmax = 300, starting from a zero network state
  2. Dismiss an initial transient of 100 steps after which the effects of the initial state have died out (state forgetting property)
  3. Collect the network states x(n) for n = nmin = 101, . . . , nmax = 300
  4. Compute the weights $\textbf{w}_i$ offline from these collected states, such that the mse becomes minimal.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge,LinearRegression
from scipy import linalg,stats # for eigenvalue and distribution for sparse connectivity
import sys
sys.path.append("./../")
from utils.timer import *
np.random.seed(42)

Create Data

In [2]:
data = np.array([np.sin(n/5) for n in range(601)])
print("Data size: ",data.size)
#1-D time series
inSize = outSize = 1

Data size:  601

Generate the reservoir

100 sparsely connected units. Achieved by using a distribution of values 0, +0.4 and −0.4 with probabilities 0.95, 0.025, 0.025. This leads to a sparse connectivity of 5 %.

The input weights are set to values of +1, -1 with equal probability.

Spectral Radius

We need its absolute value to be smaller than 1 to have echo states since we have the zero input sin(0) for n = 0.

In [3]:
resSize = 100

#Creating input-reservoir and reservoir-reservoir connections

Win = np.random.uniform(size=(resSize,inSize+1))<0.5
Win = np.where(Win==0, -1, Win)
W = 0.1 * stats.rv_discrete(name='sparse', values=([0,4,-4], [0.95,0.025,0.025])).rvs(size=(resSize,resSize),random_state=42) 

print('Computing spectral radius...')
rhoW = abs(linalg.eig(W)[0]).max()
print(f'Done: {rhoW}')

Computing spectral radius...
Done: 0.8945709957166343

Restricted case

$x(n+1) = \sigma(\textbf{W}^{in}\cdot u(n+1)_{train} + \textbf{Wx}(n) )$

Training

We want to minimize via linear regression $\textbf{W}^{out} = argmin_\textbf{W} ||y(n) - \text{tanh}(\textbf{W} \cdot x(n))||^2$, where $n \in \mathbb{N}_{train}$ the time indices, $u$ the inputs and $x$ the reservoir states, tanh as the activation on the output layer. Making use of echo state poperty of ESNs, we initialize the reservoir states for 100 timesteps. After that we will have echo states given by echo functions.

In [4]:
trainLen = 301
valLen = data.size - trainLen
initLen = 100 #We initialize the reservoir states for 100 timesteps
sigmoid = lambda k: 1 / (1 + np.exp(-k))
#Array shapes: (column vector, time)

y_train = 0.5*data[initLen:trainLen]**7
X = np.zeros((resSize,trainLen-initLen)) 
x = np.random.rand(resSize,1)  #reservoir states random init
for t in range(1,trainLen):
    u = data[t]
    x = sigmoid( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) )
    if t >= initLen:
        X[:,t-initLen] = x.ravel()

Linear Regression

In [5]:
regr = LinearRegression(fit_intercept=False)
with timer("Linear Regression"):
    regr.fit(X.T,np.arctanh(y_train.T))
    Wout = regr.coef_
print(f"MSE train: {np.square(np.dot(Wout,X)-np.arctanh(y_train)).mean()}")

Linear Regression: 0.00546574592590332 seconds
MSE train: 3.412469853669858e-09

Plot

In [6]:
plt.figure(figsize=(30,2),dpi=200)
incr = 1
plt.plot(range(0,trainLen-initLen,incr),y_train[::incr],'blue',label='Training Data')
plt.plot(range(0,trainLen-initLen,incr),np.dot(Wout,X).ravel()[::incr],'r.',alpha = 0.8,label='Fit');
plt.legend(fontsize=15,loc="upper right");
plt.xlim(0,trainLen-initLen);

Validation

In [7]:
y_pred = np.zeros((outSize,valLen))
X_val = np.zeros((resSize,valLen))
y_val = 0.5*data[trainLen:]**7

x = X[:,-1,None]
with timer("RC Forecast time"):
    for t in range(valLen):
        u = data[trainLen+t]
        x = sigmoid( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) )
        X_val[:,t] = x.ravel()
        y_pred[:,t] = np.tanh( np.dot( Wout, x ) )

mse = np.square(y_val - y_pred.ravel()).mean()
print('MSE validation: ' + str( mse ))

RC Forecast time: 0.012334108352661133 seconds
MSE validation: 3.5206412194234626e-09

Plot

In [8]:
plt.figure(figsize=(30,4),dpi=200)
incr = 1
plt.plot(range(0,valLen,incr),y_val[::incr],'blue',label='Ground Truth')
plt.plot(range(0,valLen,incr),y_pred.ravel(),'r.',alpha = 0.8,label='Prediction');
plt.legend(fontsize=15,loc="lower left");
plt.xlim(0,valLen);

Reservoir Unit States Plot

In [9]:
plt.figure(figsize=(20,5),dpi=200)
randints = np.random.randint(resSize,size=5)
randints.sort()
plt.plot(X_val[randints,:].T)
plt.legend(labels=[f"State {i+1}" for i in randints]);

General case

$x(n+1) = \sigma(\textbf{W}^{in}\cdot u(n+1)_{train} + \textbf{Wx}(n) + \textbf{W}^{back}y(n)_{teach})$

Training

In [10]:
Wback = np.random.uniform(size=(resSize,inSize+1))<0.5
Wback = np.where(Wback==0, -1, Wback)
y_train = 0.5*data[:trainLen]**7

X = np.zeros((resSize,trainLen-initLen)) 
x = np.random.rand(resSize,1)  #reservoir states random init x(0)
for t in range(1,trainLen):
    u = data[t]
    x = sigmoid( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) + np.dot(Wback, np.vstack((1,y_train[t-1]))))
    if t >= initLen:
        X[:,t-initLen] = x.ravel()
        
y_train = 0.5*data[initLen:trainLen]**7

Linear Regression

In [11]:
regr = LinearRegression(fit_intercept=False)
with timer("Linear Regression"):
    regr.fit(X.T,np.arctanh(y_train.T))
    Wout = regr.coef_
print(f"MSE train: {np.square(np.dot(Wout,X)-np.arctanh(y_train)).mean()}")

Linear Regression: 0.004116058349609375 seconds
MSE train: 6.567920152568097e-18

Plot

In [12]:
plt.figure(figsize=(30,2),dpi=200)
incr = 1
plt.plot(range(0,trainLen-initLen,incr),y_train[::incr],'blue',label='Training Data')
plt.plot(range(0,trainLen-initLen,incr),np.dot(Wout,X).ravel()[::incr],'r.',alpha = 0.8,label='Fit');
plt.legend(fontsize=15,loc="upper right");
plt.xlim(0,trainLen-initLen);

Validation

In [13]:
y_pred = np.zeros((outSize,valLen))
X_val = np.zeros((resSize,valLen))
y_val = 0.5*data[trainLen-1:]**7

x = X[:,-1,None]
with timer("RC Forecast time"):
    for t in range(valLen):
        u = data[trainLen+t]
        x = sigmoid( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) + np.dot(Wback, np.vstack((1,y_val[t]))))
        X_val[:,t] = x.ravel()
        y_pred[:,t] = np.tanh( np.dot( Wout, x ) )

y_val = 0.5*data[trainLen:]**7

mse = np.square(y_val - y_pred.ravel()).mean()
print('MSE validation: ' + str( mse ))

RC Forecast time: 0.014831066131591797 seconds
MSE validation: 6.323126578792129e-18

Plot

In [14]:
plt.figure(figsize=(30,4),dpi=200)
incr = 1
plt.plot(range(0,valLen,incr),y_val[::incr],'blue',label='Ground Truth')
plt.plot(range(0,valLen,incr),y_pred.ravel(),'r.',alpha = 0.8,label='Prediction');
plt.legend(fontsize=15,loc="lower left");
plt.xlim(0,valLen);

Reservoir Unit States Plot

In [15]:
plt.figure(figsize=(20,5),dpi=200)
randints = np.random.randint(resSize,size=4)
randints.sort()
plt.plot(X_val[randints,:].T)
plt.legend(labels=[f"State {i+1}" for i in randints]);

Comparison with RNN

Importing Tensorflow

In [16]:
import tensorflow as tf
from tensorflow.keras import Model, Sequential
from tensorflow.keras.layers import Dense, Dropout, BatchNormalization, Input, Flatten, LSTM
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau, Callback, TensorBoard

Create the RNN

In [17]:
model = Sequential()
model.add(LSTM(20, return_sequences=True, input_shape=(1,1)))
model.add(BatchNormalization())
model.add(LSTM(20, return_sequences=False))
model.add(BatchNormalization())
model.add(Dense(10,'relu'))
model.add(BatchNormalization())
model.add(Dense(outSize))
model.compile(optimizer=Adam(lr=0.001), loss='mse')
model.summary()

Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
lstm (LSTM)                  (None, 1, 20)             1760      
_________________________________________________________________
batch_normalization (BatchNo (None, 1, 20)             80        
_________________________________________________________________
lstm_1 (LSTM)                (None, 20)                3280      
_________________________________________________________________
batch_normalization_1 (Batch (None, 20)                80        
_________________________________________________________________
dense (Dense)                (None, 10)                210       
_________________________________________________________________
batch_normalization_2 (Batch (None, 10)                40        
_________________________________________________________________
dense_1 (Dense)              (None, 1)                 11        
=================================================================
Total params: 5,461
Trainable params: 5,361
Non-trainable params: 100
_________________________________________________________________

Training

Optimization Problem: $argmin_{\theta} ||\text{RNN}^\theta(u(n)) - 0.5*u^7(n+1)||,n \in \mathbb{N}_{train}$

In [18]:
class CustomLogger(Callback):
    def __init__(self, epochs_tot):
        super().__init__()
        self.epochs=epochs_tot
        self.threshold = 1000

    def on_epoch_begin(self, epoch, logs=None):
        self.time_begin = time.time()
    def on_epoch_end(self,epoch,logs=None):
        duration = time.time() - self.time_begin
        if (epoch+1) % self.threshold == 0 or epoch==0:
            print(f"The average {self.model.loss} in epoch {epoch+1} is {logs['loss']:1.8f}. Duration: {duration} seconds")
            
epochs = 10000
batch_size = 301
reduceLR = ReduceLROnPlateau(monitor='loss', factor=0.5, patience=epochs//10, min_lr=1e-7, verbose=1)
earlyStopping = EarlyStopping(monitor='loss', patience=epochs//2)
logger = CustomLogger(epochs)
with timer("RNN Training time"):
    history = model.fit(data[:trainLen,None,None],0.5*data[:trainLen,None,None]**7,batch_size=batch_size,epochs=epochs, verbose=0, callbacks=[earlyStopping,logger])

mse in epoch 1 is 0.63467389. Duration: 2.1376068592071533 seconds
mse in epoch 1000 is 0.00014003. Duration: 0.002429962158203125 seconds
mse in epoch 2000 is 0.00007305. Duration: 0.0025076866149902344 seconds
mse in epoch 3000 is 0.00005781. Duration: 0.0024487972259521484 seconds
mse in epoch 4000 is 0.00004566. Duration: 0.0024340152740478516 seconds
mse in epoch 5000 is 0.00003229. Duration: 0.0025577545166015625 seconds
mse in epoch 6000 is 0.00001658. Duration: 0.0024099349975585938 seconds
mse in epoch 7000 is 0.00001212. Duration: 0.0024518966674804688 seconds
mse in epoch 8000 is 0.00001044. Duration: 0.0025589466094970703 seconds
mse in epoch 9000 is 0.00001231. Duration: 0.0025331974029541016 seconds
mse in epoch 10000 is 0.00002528. Duration: 0.0024619102478027344 seconds
RNN Training time: 28.195341110229492 seconds

Plot

In [19]:
plt.figure(figsize=(30,2),dpi=200)
incr = 1
plt.plot(range(0,trainLen,incr),0.5*data[:trainLen][::incr]**7,'blue',label='Ground Truth')
plt.plot(range(0,trainLen,incr),history.model.predict(data[:trainLen,None,None])[::incr],'r.',alpha = 0.8,label='Prediction');
plt.legend(fontsize=15,loc="lower left");
plt.xlim(0,trainLen);

Validation

In [20]:
with timer("RNN Forecast time"):
    RNN_val = history.model.predict(data[trainLen:,None,None])
RNN_mse = np.square(RNN_val-y_val).mean()
print('MSE validation: ', RNN_mse)

RNN Forecast time: 0.03289008140563965 seconds
MSE validation:  0.09876100335921453

Plot

In [21]:
plt.figure(figsize=(30,4),dpi=200)
incr = 1
plt.plot(range(0,valLen,incr),y_val[::incr],'blue',label='Ground Truth')
plt.plot(range(0,valLen,incr),RNN_val[::incr],'r.',alpha = 0.8,label='Prediction');
plt.legend(fontsize=15,loc ="lower left");
plt.xlim(0,valLen);