回声状态网络（ESN）的公式推导及代码实现

1. ESN的任务

给定一段信号：$u(0)，u(1)，\cdot \cdot \cdot ，u(N_t-1)$
和目标值：$v(1)，v(2)，\cdot \cdot \cdot ，v(N_t)$
学习一个黑箱模型M使得我们可以预测$v(N_t+1)，v(N_t+2)，\cdot \cdot \cdot$

优势：与传统的递归神经网络相比，ESN最大的优势是简化了网络的训练过程，解决了传统递归神经网络结构难以确定、训练算法过于复杂的问题，同时也可以克服递归网络存在的记忆消减等问题（ESN的训练方法与传统的递归神经网络有本质不同）。

ESN解决问题的思想：使用大规模随机稀疏网络（储备池）作为信息处理媒介，将输入信号从低维的输入空间映射到高维的状态空间，在高维的状态空间采用线性回归方法对网络的部分连接权进行训练，而其他连接权随机产生，并在网络训练过程中保持不变。这种思想在Steil关于传统递归神经网络的经典算法（Atiya-Parlos）的研究中也得到了验证：递归神经网络输出连接权改变迅速，而内部连接权则以高度耦合的方式缓慢改变。也就是说，如果递归神经网络内部连接权选择合适，在对网络进行训练时可以忽略内部连接权的改变。

2. ESN的结构和训练步骤

ESN的结构如上图所示，其中：
中间的大圆圈叫做“储备池”，它有以下特点：
（1）包含数目较多的神经元（与经典神经网络相比）；
（2）神经元之间的连接关系随机产生；
（3）神经元之间的连接具有稀疏性。

$u\in R^{M\ast 1}$，$W_{IR}\in R^{N\ast M}$，$W_{res}\in R^{N\ast N}$，$r\in R^{N\ast 1}$，$W_{RO}\in R^{L\ast N}$，$v\in R^{L\ast 1}$

上面的两个参数矩阵$W_{IR}\in R^{N\ast M}$和$W_{res}\in R^{N\ast N}$都是事先给定的数值，在训练的过程中只需要计算$W_{RO}\in R^{L\ast N}$即可。

整个计算过程如下所示：

（1）从输入到储备池(reservoir)的运算：$W_{IR}\ast u(t)$

（2）储备池中$r(t)$的更新：$r(t+\Delta t)=f[W_{res}\ast r(t)+W_{IR}\ast u(t)]$

（3）从储备池到输出：$W_{RO}\ast r(t)$

（4）损失函数：$L={\sum }_{t=d+1}^{N_t}|v(t)-W_{RO}\ast r(t){|}^2$

（6）使损失函数最小化的推导过程：

3. ESN的预测步骤

当$W_{RO}$确定之后，库的输出为：

$u(t)=W_{RO}\ast r(t)$

$r(t+\Delta t)=f[W_{res}\ast r(t)+W_{IR}\ast u(t)]$

$u(t+\Delta t)=W_{RO}\ast r(t+\Delta t)$

···

热启动方式：使用训练步中最后一个阶段的库状态作为预测中的$r(t)$；

冷启动方式：使用一个新的数据作为库的初始值

预测的时候不会再给单独的输入了，而是会将输出作为输入进行递推计算。

一般$W_{IR}$各元素会初始化为$[-\alpha ,\alpha ]$之间的均匀分布。

每个输入$u(t)$都会和$N/M$个库中的节点相连，因为输入个数时M，库中有N个节点，
即：$u\in R^{M\ast 1}$，$r\in R^{N\ast 1}$

$W_{res}$通常是一个大型，稀疏，有向或无向的随机网络，平均度为k，谱半径$\rho (W_{res})$是$W_{res}$最大的特征值。

$W_{res}$会初始化为一个稀疏矩阵。

4. ESN的代码

import pickle import numpy as np import matplotlib.pyplot as plt class ESN(): def __init__(self, data, N=1000, rho=1, sparsity=3, T_train=2000, T_predict=1000, T_discard=200, eta=1e-4, seed=2050): self.data = data self.N = N # reservoir size 库的大小 self.rho = rho # spectral radius 谱半径 self.sparsity = sparsity # average degree 平均度 sparsity：稀疏性 self.T_train = T_train # training steps self.T_predict = T_predict # prediction steps self.T_discard = T_discard # discard first T_discard steps discard：丢弃 self.eta = eta # regularization constant 正则化常数 self.seed = seed # random seed def initialize(self): """ 对连接权矩阵W_IR和W_res进行初始化 其中W_IR(N*1)是从输入到库的连接权矩阵，W_res(N*N)是从库到输出的连接权矩阵 """ if self.seed > 0: np.random.seed(self.seed) # 生成形状为N * 1的，元素为[-1, 1]之间的随机值的矩阵 self.W_IR = np.random.rand(self.N, 1) * 2 - 1 # [-1, 1] uniform # 生成形状为N * N的，元素为[0, 1]之间的随机值的矩阵 W_res = np.random.rand(self.N, self.N) # 将W_res中大于self.sparsity / self.N的元素置0 W_res[W_res > self.sparsity / self.N] = 0\ # np.linalg.eigvals(W_res)求出W_res的特征值，W_res矩阵除以自身模最大的特征值的模 W_res /= np.max(np.abs(np.linalg.eigvals(W_res))) # 在乘以谱半径 W_res *= self.rho # set spectral radius = rho self.W_res = W_res def train(self): u = self.data[:, :self.T_train] # traning data T_train = 2000 assert u.shape == (1, self.T_train) r = np.zeros((self.N, self.T_train + 1)) # initialize reservoir state r(N*(T_train + 1)) for t in range(self.T_train): # @是Python3.5之后加入的矩阵乘法运算符 r[:, t+1] = np.tanh(self.W_res @ r[:, t] + self.W_IR @ u[:, t]) # disgard first T_discard steps r丢弃前T_discard步变成r_p self.r_p = r[:, self.T_discard+1:] # length=T_train-T_discard v = self.data[:, self.T_discard+1:self.T_train+1] # target self.W_RO = v @ self.r_p.T @ np.linalg.pinv( self.r_p @ self.r_p.T + self.eta * np.identity(self.N)) train_error = np.sum((self.W_RO @ self.r_p - v)  2) print('Training error: %.4g' % train_error) def predict(self): u_pred = np.zeros((1, self.T_predict)) # u_pred是形状为(1, self.T_predict)的全零矩阵 r_pred = np.zeros((self.N, self.T_predict)) # r_pred是形状为(N, self.T_predict)的全零矩阵 r_pred[:, 0] = self.r_p[:, -1] # warm start 热启动 for step in range(self.T_predict - 1): u_pred[:, step] = self.W_RO @ r_pred[:, step] r_pred[:, step + 1] = np.tanh(self.W_res @ r_pred[:, step] + self.W_IR @ u_pred[:, step]) u_pred[:, -1] = self.W_RO @ r_pred[:, -1] self.pred = u_pred def plot_predict(self): ground_truth = self.data[:, self.T_train: self.T_train + self.T_predict] plt.figure(figsize=(12, 4)) plt.plot(self.pred.T, 'r', label='predict', alpha=0.6) plt.plot(ground_truth.T, 'b', label='True', alpha=0.6) plt.show() def calc_error(self): ground_truth = self.data[:, self.T_train: self.T_train + self.T_predict] rmse_list = [] for step in range(1, self.T_predict+1): error = np.sqrt( np.mean((self.pred[:, :step] - ground_truth[:, :step])  2)) rmse_list.append(error) return rmse_list if __name__ == "__main__": # http://minds.jacobs-university.de/mantas/code data = np.load('mackey_glass_t17.npy') # data.shape = (10000,) data = np.reshape(data, (1, data.shape[0])) # data.shape = (1, 10000) print(data.shape) esn = ESN(data) esn.initialize() esn.train() esn.predict() esn.plot_predict() 

参考