基本模型

输入层和循环层

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

E(y,ytarget)=1Nyi=1Ny1Tn=1T(yi(n)yitarget(n))2=1Nyn=1T1Tyiyitarget22E(\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)的更新式如下

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

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

输出层

输出层如下定义

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

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

网络的整体结构

y(n)=Wout[1;u(n);x(n)]x(n)=tanh(Win[1;u(n)]+Wx(n1))=tanh(Win[1;u(n)]+Wtanh(Win[1;u(n1)]+Wx(n2)))=tanh(Win[1;u(n)]+Wtanh(Win[1;u(n1)]++Wtanh(Win[1;u(1)]+Wx(0))))\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}

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

储备池的生成

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

储备池的作用

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

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

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

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

E(n)u(nk)=E(n)x(n)x(n)x(n1)x(n1)x(n2)x(nk+1)x(nk)x(nk)u(nk)k=0,1,2,,n1\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

E(n)x(n)= E(n)y(n)y(n)x(n)=2[y1(n)y1target(n),,yNy(n)yNytarget(n)]Wxout\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}

x(ni+1)x(ni)=x(ni+1)η(ni+1)η(ni+1)x(ni)=[x1(ni+1)η1(ni+1)000x2(ni+1)η2(ni+1)000xNx(ni+1)ηNx(ni+1)]Wi=1,2,,k\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}

x(nk)u(nk)=x(nk)η(nk)η(nk)u(nk)=[x1(nk)η1(nk)000x2(nk)η2(nk)000xNx(nk)ηNx(nk)]Wuout\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}

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

于是有

E(n)u(nk)=E(n)x(n)i=1k([x1(ni+1)η1(ni+1)00xNx(ni+1)ηNx(ni+1)]W)x(nk)u(nk)\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的输出层是一个典型的线性前馈结构, 其表达式可以写成矩阵的形式如下

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

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

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

Ytarget=WoutX\boldsymbol{Y}^{target} = \boldsymbol{W}^{out}\boldsymbol{X}

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

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

Wout=argminWoutn=1Ti=1Ny((yi(n)yitarget(n))2+β(wi,nout)2)=argminWout(YYtargetF2+βWoutF2)=argminWout(WoutXYtargetF2+βWoutF2)\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}

其中wi,nout\boldsymbol{w}_{i,n}^{out}表示Wout\boldsymbol{W}^{out}的第ii行第nn列的元素

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

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

WoutXYtargetF2Wout=tr((WoutXYtarget)T(WoutXYtarget))Wout=tr((WoutXYtarget)T(WoutXYtarget))tr(WoutXYtarget)tr(WoutXYtarget)Wout=(WoutXYtarget)XTWoutF2Wout=tr((Wout)T(Wout))Wout=Wout\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}

(WoutXYtarget)XT+βWout=0(\boldsymbol{W}^{out}\boldsymbol{X}-\boldsymbol{Y}^{target})\boldsymbol{X}^T+\beta\boldsymbol{W}^{out} = \boldsymbol{0}

得到

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

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

在时间序列上的应用

Mackey-Glass时间序列

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

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

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

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

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

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

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

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

递推公式如下

{x(n+1)=x(n)+h(0.2x(t17)1+x10(t17)0.1x(n)),n17x(n)=cos(n),0n17\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}

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

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

Figure_10

下面是生成数据集的相图

Figure_4

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

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

Python实现

  • 导入数据集
1
2
3
4
5
trainLen = 2000 ## 设置训练集样本个数
testLen = 2000 ## 设置测试集样本个数
initLen = 100 ## 设置初始集样本个数(预热作用,不参与模型的训练)
## 导入数据集
data = np.loadtxt('MackeyGlass_t17.txt')
  • 设置初始参数
1
2
3
inSize = outSize = 1 ## 输入和输出维度
resSize = 1000 ## 储备池的大小
a = 0.3 ## 泄露率
  • 生成储备池

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

1
2
Win = np.random.rand(resSize,1+inSize) - 0.5 ## 输入矩阵
W = np.random.rand(resSize,resSize) - 0.5 ## 循环矩阵
  • 计算谱半径

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

1
2
rhoW = max(abs(linalg.eig(W)[0])) ## 计算谱半径
W *= 1.25 / rhoW ## 对矩阵进行放缩
  • 设置目标向量

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

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

1
2
## 设置目标向量
Yt = data[None,initLen+1:trainLen+1]
  • 计算储备池的训练样本点的输出

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

x~(n)=tanh(Win[1;u(n)]+Wx(n1)),x(n)=(1α)x(n1)+αx~(n)\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*}

1
2
3
4
5
6
7
8
9
10
## 为储备池的输出向量分配内存
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]
  • 训练输出层的矩阵

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

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

1
2
3
4
5
6
## 正则项系数
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) ) )
  • 对训练好的模型进行评估

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

1
u = data[trainLen]

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

1
2
3
4
5
6
7
8
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个点的均方误差

1
2
3
4
errorLen = 300
mse = sum( np.square( data[trainLen+1:trainLen+errorLen+1] -
Y[0,0:errorLen] ) ) / errorLen
print('MSE = ' + str( mse ))
1
2
## 输出
MSE = 0.00108

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

Figure_1

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