回声状态网络

基本模型

输入层和循环层

对于一个输入信号$\boldsymbol{u}(n)\in\mathbb{R}^{N_u}$, 给定期望的目标输出信号$\boldsymbol{y}^{target}(n)\in\mathbb{R}^{N_y}$, 其中$n=1,2,\cdots,T$ 是离散的时间点, $T$是训练集中数据点的个数. 我们的目标是训练一个模型, 使得模型的输出$\boldsymbol{y}(n)\in\mathbb{R}^{N_y}$尽可能地匹配期望的目标输出信号$\boldsymbol{y}^{target}(n)\in\mathbb{R}^{N_y}$, 即最小化$\boldsymbol{y}$和$\boldsymbol{y}^{target}$之间的误差测度$E(\boldsymbol{y},\boldsymbol{y}^{target})$. 误差测度一般采用均方误差(MSE), 例如均方根误差(RMSE)

$$E(\boldsymbol{y},\boldsymbol{y}^{target})= \frac{1}{N_y} \sum_{i=1}^{N_y} \sqrt{\frac{1}{T}\sum_{n=1}^T(y_i(n)-y_i^{target}(n))^2} \\ = \frac{1}{N_y} \sum_{n=1}^{T} \sqrt{\frac{1}{T}\|\boldsymbol{y}_i-\boldsymbol{y}^{target}_i\|^2_2}$$

回声状态网络的结构和循环神经网络(RNN)类似, 它的循环层(储备池 , Reservoir)的更新式如下

$$\begin{gather} \tilde{\boldsymbol{x}}(n) = {\rm tanh}(\boldsymbol{W}^{in}\cdot[1;\boldsymbol{u}(n)]+\boldsymbol{W}\cdot\boldsymbol{x}(n-1)), \\ \boldsymbol{x}(n) = (1-\alpha)\cdot\boldsymbol{x}(n-1)+\alpha\cdot\tilde{\boldsymbol{x}}(n) \end{gather}$$

其中$\boldsymbol{x}(n)\in\mathbb{R}^{N_x}$是$n$时刻储备池神经元的激活向量(对应RNN中的循环层的状态); $\tilde{\boldsymbol{x}}(n)\in\mathbb{R}^{N_x}$ 是一个中间量, 用于求$\boldsymbol{x}(n)$的泄露整合(leaky integration); ${\rm tanh(\cdot)}$是采用的激活函数, 这里是按元素进行计算的(element-wise); $[\cdot;\cdot]$代表向量在竖直方向的连接, $[1;\boldsymbol{u}(n)]$这样表示是为了省略偏置项$\boldsymbol{b}^{in}$, 因为这样在矩阵$\boldsymbol{W}^{in}\in\mathbb{R}^{N_x\times(1+N_{u})}$第一列中已经隐含了偏置项; $\boldsymbol{W}^{in}\in\mathbb{R}^{N_x\times(1+N_{u})}$是输入系数矩阵, $\boldsymbol{W}\in\mathbb{R}^{N_x\times N_x}$是循环系数矩阵

这里采用泄露整合(leaky integration)来更新$\boldsymbol{x}(n)$, 其中参数$\alpha\in(0,1]$称为泄露率; 当$\alpha=1$时$\boldsymbol{x}(n)\equiv\tilde{\boldsymbol{x}}(n)$是模型的一种特殊情况, 此时模型不采用泄露整合.

输出层

输出层如下定义

$$\boldsymbol{y}(n) = \boldsymbol{W}^{out}\cdot[1;\boldsymbol{u}(n);\boldsymbol{x}(n)],$$

其中$\boldsymbol{y}(n)\in\mathbb{R}^{N_y}$是网络的输出, $\boldsymbol{W}^{out}\in\mathbb{R}^{N_y\times(1+N_u+N_x)}$是输出权重矩阵, $[\cdot;\cdot;\cdot]$还是代表向量在竖直方向的连接.

网络的整体结构

$$\begin{split} \boldsymbol{y}(n) &= \boldsymbol{W}^{out}\cdot[1;\boldsymbol{u}(n);\boldsymbol{x}(n)] \\\\ \boldsymbol{x}(n) &= {\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(n)]+\boldsymbol{W}\boldsymbol{x}(n-1)\right) \\ &= {\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(n)]+\boldsymbol{W} {\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(n-1)]+\boldsymbol{W}\boldsymbol{x}(n-2)\right)\right) \\ &\vdots \\ &= {\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(n)]+\boldsymbol{W} {\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(n-1)]+\cdots+\boldsymbol{W}{\rm tanh}\left( \boldsymbol{W}^{in} [1;\boldsymbol{u}(1)]+\boldsymbol{W}\boldsymbol{x}(0)\right)\right)\right) \end{split}$$

可以看到网络在$n$时刻的输出$\boldsymbol{y}(n)$是和前$n$个时刻的输入$\boldsymbol{u}(n),\boldsymbol{u}(n-1),\cdots,\boldsymbol{u}(1)$和初始状态$\boldsymbol{x}(0)$相关的.

储备池的生成

为了生成一个"好的"储备池, 首先我们要了解储备池的作用是什么.

储备池的作用

• 将输入向量非线性映射到高维空间
• 存储网络的状态(隐含了前面的所有输入数据)

各时刻的输入对结果的影响

记$n$时刻的输出$\boldsymbol{y}(n)$和目标$\boldsymbol{y}^{target}(n)$的误差为$E(n)$

$$E(n) = \|\boldsymbol{y}(n)-\boldsymbol{y}^{target}(n)\|_2^2$$

$$\frac{\partial E(n)}{\partial \boldsymbol{u}(n-k)}= \frac{\partial E(n)}{\partial \boldsymbol{x}(n)} \frac{\partial \boldsymbol{x}(n)}{\partial \boldsymbol{x}(n-1)} \frac{\partial \boldsymbol{x}(n-1)}{\partial \boldsymbol{x}(n-2)} \cdots \frac{\partial \boldsymbol{x}(n-k+1)}{\partial \boldsymbol{x}(n-k)} \frac{\partial \boldsymbol{x}(n-k)}{\partial \boldsymbol{u}(n-k)} \qquad k=0,1,2,\cdots,n-1$$

$$\begin{split} \frac{\partial E(n)}{\partial \boldsymbol{x}(n)} =\ \frac{\partial E(n)}{\partial \boldsymbol{y}(n)} \frac{\partial \boldsymbol{y}(n)}{\partial \boldsymbol{x}(n)} = 2\begin{bmatrix} y_1(n)-y^{target}_1(n),&\cdots&,y_{N_y}(n)-y^{target}_{N_y}(n) \end{bmatrix} \boldsymbol{W}^{out}_x \end{split}$$

$$\begin{split} \frac{\partial \boldsymbol{x}(n-i+1)}{\partial \boldsymbol{x}(n-i)} &= \frac{\partial \boldsymbol{x}(n-i+1)}{\partial \boldsymbol{\eta}(n-i+1)} \frac{\partial \boldsymbol{\eta}(n-i+1)}{\partial \boldsymbol{x}(n-i)} \\ &= \begin{bmatrix} \frac{\partial x_{1}(n-i+1)}{\partial \eta_{1}(n-i+1)}&0&\cdots&0\\ 0&\frac{\partial x_{2}(n-i+1)}{\partial \eta_{2}(n-i+1)}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\frac{\partial x_{N_x}(n-i+1)}{\partial \eta_{N_x}(n-i+1)} \end{bmatrix} \boldsymbol{W} \\\\ i&=1,2,\cdots,k \end{split}$$

$$\begin{split} \frac{\partial \boldsymbol{x}(n-k)}{\partial \boldsymbol{u}(n-k)} &= \frac{\partial \boldsymbol{x}(n-k)}{\partial \boldsymbol{\eta}(n-k)} \frac{\partial \boldsymbol{\eta}(n-k)}{\partial \boldsymbol{u}(n-k)} \\ &= \begin{bmatrix} \frac{\partial x_{1}(n-k)}{\partial \eta_{1}(n-k)}&0&\cdots&0\\ 0&\frac{\partial x_{2}(n-k)}{\partial \eta_{2}(n-k)}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\frac{\partial x_{N_x}(n-k)}{\partial \eta_{N_x}(n-k)} \end{bmatrix} \boldsymbol{W}^{out}_u \end{split}$$

这里采用分子布局进行表示

于是有

$$\frac{\partial E(n)}{\partial \boldsymbol{u}(n-k)}= \frac{\partial E(n)}{\partial \boldsymbol{x}(n)} \prod_{i=1}^{k} \left( \begin{bmatrix} \frac{\partial x_{1}(n-i+1)}{\partial \eta_{1}(n-i+1)}&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&\frac{\partial x_{N_x}(n-i+1)}{\partial \eta_{N_x}(n-i+1)} \end{bmatrix} \boldsymbol{W} \right) \frac{\partial \boldsymbol{x}(n-k)}{\partial \boldsymbol{u}(n-k)}$$

输出层的训练

岭回归

因为ESN的输出层是一个典型的线性前馈结构, 其表达式可以写成矩阵的形式如下

$$\boldsymbol{Y} = \boldsymbol{W}^{out}\boldsymbol{X},$$

其中$\boldsymbol{Y}\in\mathbb{R}^{N_y\times T}$是所有的$\boldsymbol{y}(n),n=1,\cdots,T$, $\boldsymbol{X}\in\mathbb{R}^{(1+N_u+N_x)\times T}$是所有的$[1;\boldsymbol{u}(n);\boldsymbol{x}(n)],n=1,\cdots,T$. 为书写方便起见, 这里只用一个$\boldsymbol{X}$来表示$[1;\boldsymbol{U};\boldsymbol{X}]$.

我们的目标是找到一个最优的输出权重矩阵$\boldsymbol{W}^{out}\in\mathbb{R}^{N_y\times(1+N_u+N_x)}$使得$\boldsymbol{Y}$和$\boldsymbol{Y}^{target}$之间的误差最小. 实际上这个问题等价于在求解一个超定线性方程组

$$\boldsymbol{Y}^{target} = \boldsymbol{W}^{out}\boldsymbol{X}$$

其中$\boldsymbol{Y}^{target}\in\mathbb{R}^{N_y\times T}$和$\boldsymbol{X}\in\mathbb{R}^{(1+N_u+N_x)\times T}$是已知的, 求$\boldsymbol{W}^{out}\in\mathbb{R}^{N_y\times(1+N_u+N_x)}$的最小二乘解. 这个方程组是超定的, 因为一般而言$T\gg 1+N_u+N_x$, 限定方程的个数要远大于待定参数的个数.

这里不是直接求最小二乘解, 我们给目标函数增加一个正则项

$$\begin{split} \boldsymbol{W}^{out} &= arg\min_{\boldsymbol{W}^{out}} \sum_{n=1}^T\sum_{i=1}^{N_y}\left((y_i(n)-y_i^{target}(n))^2+\beta(\boldsymbol{w}^{out}_{i,n})^2\right)\\ &= arg\min_{\boldsymbol{W}^{out}}\left( \|\boldsymbol{Y}-\boldsymbol{Y}^{target}\|_F^2 + \beta\|\boldsymbol{W}^{out}\|_F^2\right) \\ &= arg\min_{\boldsymbol{W}^{out}}\left( \|\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target}\|_F^2 + \beta\|\boldsymbol{W}^{out}\|_F^2\right) \end{split}$$

其中$\boldsymbol{w}_{i,n}^{out}$表示$\boldsymbol{W}^{out}$的第$i$行第$n$列的元素

$\|\cdot\|_F$代表Frobenius范数, 对于矩阵$A=[a_{ij}]_{m\times n}$ , 其F-范数定义为$\|A\|_F = \sqrt{tr(A^TA)}=\sqrt{\sum_{i=1}^m\sum_{j=1}^n a_{ij}^2}$.

对上面的式子关于$\boldsymbol{W}^{out}$求偏导, 由矩阵的求导公式有

$$\begin{split} \frac{\partial \|\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target}\|_F^2}{\partial \boldsymbol{W}^{out}} &= \frac{\partial {\rm tr}\left((\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})^T(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})\right)}{\partial \boldsymbol{W}^{out}} \\ &= \frac{\partial {\rm tr}\left((\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})^T(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})\right)}{\partial {\rm tr}\left(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target}\right)} \frac{\partial {\rm tr}\left(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target}\right)}{\partial \boldsymbol{W}^{out}} \\ &= (\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})\boldsymbol{X}^T \\ \notag \\ \frac{\partial \|\boldsymbol{W}^{out}\|_F^2}{\partial \boldsymbol{W}^{out}} &= \frac{\partial {\rm tr}\left((\boldsymbol{W}^{out})^T(\boldsymbol{W}^{out})\right)}{\partial \boldsymbol{W}^{out}} = \boldsymbol{W}^{out} \end{split}$$

令

$$(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})\boldsymbol{X}^T+\beta\boldsymbol{W}^{out} = \boldsymbol{0}$$

得到

$$\boldsymbol{W}^{out} = \boldsymbol{Y}^{target}\boldsymbol{X}^T\left (\boldsymbol{X}\boldsymbol{X}^T+\beta\boldsymbol{I}\right)^{-1}$$

其中当$\beta=0$时, $\boldsymbol{W}^{out}$对应的结果就是一般的最小二乘解.

在时间序列上的应用

Mackey-Glass时间序列

Mackey-Glass方程是由Michael Mackey和Leon Glass在1977年在Science上发表的论文生理控制系统的振荡和混沌(Oscillation and chaos in physiological control systems)中提出的, 该方程为非线性时滞微分方程, 它的一般形式如下

$$\frac{dx(t)}{dt} = \beta\frac{x(t-\tau)}{1+x^{n}(t-\tau)}-\gamma x(t), \qquad \beta,\gamma,n>0,$$

其中$\beta,\gamma,n$都是正的实参数, $\tau$表示相位差.

我们取参数$\beta=0.2,\gamma=0.1,n=10$, 得到方程的一种形式

$$\frac{dx(t)}{dt} = 0.2\frac{x(t-\tau)}{1+x^{10}(t-\tau)}-0.1 x(t)$$

在这个形式下, 当相位差$\tau\ge 17$时, 由方程构造的序列是无周期的, 混沌的并且不随时间收敛和也不随时间发散.

我们这这里使用欧拉方法对Mackey-Glass方程进行数值求解, 生成时间序列. 取相位差$\tau=17$, 初始值由余弦函数如下生成

$$\begin{equation*} x(0)=\cos(0),\quad x(1)=\cos(1)\quad ,\cdots,\quad x(17)=cos(17) \end{equation*}$$

递推公式如下

$$\begin{cases} x(n+1) = x(n) + h\left(0.2\frac{x(t-17)}{1+x^{10}(t-17)}-0.1 x(n)\right),& n \ge17\\ x(n) = \cos(n),&0\le n\le17 \end{cases}$$

其中$h$是时间步长, 这里取$h=1$

由此我们得到数据集, 下图是部分样本点的展示

下面是生成数据集的相图

可以看出数据不呈现任何周期性, 不收敛也不发散, 在一定的区间上呈现混沌状态

使用ESN预测Mackey-Glass时间序列

Python实现

• 导入数据集

trainLen = 2000 ## 设置训练集样本个数
testLen = 2000  ## 设置测试集样本个数
initLen = 100   ## 设置初始集样本个数(预热作用,不参与模型的训练)
## 导入数据集
data = np.loadtxt('MackeyGlass_t17.txt')

• 设置初始参数

inSize = outSize = 1 ## 输入和输出维度
resSize = 1000       ## 储备池的大小
a = 0.3              ## 泄露率

• 生成储备池

生成$[-0.5,0.5)$上均匀分布的随机矩阵

Win = np.random.rand(resSize,1+inSize) - 0.5 ## 输入矩阵
W = np.random.rand(resSize,resSize) - 0.5    ## 循环矩阵

• 计算谱半径

计算循环矩阵的谱半径并对其进行放缩

rhoW = max(abs(linalg.eig(W)[0])) ## 计算谱半径
W *= 1.25 / rhoW                  ## 对矩阵进行放缩

• 设置目标向量

使用下一时刻的数据作为上一时刻输入的目标值

$$\begin{equation*} \boldsymbol{y}^{target}(n) = \boldsymbol{u}(n+1) \end{equation*}$$

## 设置目标向量
Yt = data[None,initLen+1:trainLen+1] 

• 计算储备池的训练样本点的输出

按照下面的公式计算储备池的输出

$$\begin{gather*} \tilde{\boldsymbol{x}}(n) = {\rm tanh}(\boldsymbol{W}^{in}\cdot[1;\boldsymbol{u}(n)]+\boldsymbol{W}\cdot\boldsymbol{x}(n-1)), \\ \boldsymbol{x}(n) = (1-\alpha)\cdot\boldsymbol{x}(n-1)+\alpha\cdot\tilde{\boldsymbol{x}}(n) \end{gather*}$$

## 为储备池的输出向量分配内存
X = np.zeros((1+inSize+resSize,trainLen-initLen))

## 按照公式计算储备池的输出,并储存在变量X中 
x = np.zeros((resSize,1))
for t in range(trainLen):
    u = data[t]
    x = (1-a)*x + a*np.tanh( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) )
    if t >= initLen:
        X[:,t-initLen] = np.vstack((1,u,x))[:,0]

• 训练输出层的矩阵

根据公式计算输出矩阵的岭回归最小二乘解

$$\begin{equation*} \boldsymbol{W}^{out} = \boldsymbol{Y}^{target}\boldsymbol{X}^T\left (\boldsymbol{X}\boldsymbol{X}^T+\beta\boldsymbol{I}\right)^{-1} \end{equation*}$$

## 正则项系数
reg = 1e-8
## 计算Wout的岭回归最小二乘解
X_T = X.T
Wout = np.dot( np.dot(Yt,X_T), linalg.inv( np.dot(X,X_T) + \
    reg*np.eye(1+inSize+resSize) ) )

• 对训练好的模型进行评估

紧接着训练集, 将训练集最后一个数据的下一个数据作为网络的输入

u = data[trainLen]

将网络当前时刻的输出作为下一时刻的输入

Y = np.zeros((outSize,testLen))
## 循环计算data[trainLen]后2000个状态
for t in range(testLen):
    x = (1-a)*x + a*np.tanh( np.dot( Win, np.vstack((1,u)) ) + np.dot( W, x ) )
    y = np.dot( Wout, np.vstack((1,u,x)) )
    Y[:,t] = y
    ## 将网络当前时刻的输出作为下一时刻的输入
    u = y

计算前300个点的均方误差

errorLen = 300
mse = sum( np.square( data[trainLen+1:trainLen+errorLen+1] - 
    Y[0,0:errorLen] ) ) / errorLen
print('MSE = ' + str( mse ))

## 输出
MSE = 0.00108

对比预测序列和目标序列的差异

从上图可以看出模型对时间序列前300个点的预测都呈现良好的效果, 而随着时间的推移预测值逐渐偏离目标值.