回声状态网络(ESN)教程

基础概念

回声状态网络(Echo State Network)提出于2001年，曾经是研究的热点，但近年来随着RNN，LSTM与其它一些变种的网络的出现，现在研究比较少了，但是其在时间序列预测上还有着很不错的应用。传统的MLP网络的隐层是一层层的全连接的神经元，而ESN引入了一个储备池计算模式来替代原始的隐层，这个储备池是什么呢？先来看下下图：
回声状态网络

网络结构依次是输入层，储备池和输出层，所谓的储备池就是中间的部分。这个储备池的特点是: (1)储备池中神经元的连接状态是随机的，即神经元之间是否建立连接并不是我们人工确定的；(2)储备池中的连接权重是固定的，不像传统的MLP网络使用梯度下降进行权重的更新。这样做的好处是：(1)大大降低了训练的计算量；(2)一定程度上避免了梯度下降的优化算法中出现的局部极小情况；(3)此外，在很多问题上确实有着不错的建模能力。ESN的基本思想就是由储备池生成一个随输入不断变化的复杂动态空间，当这个状态空间足够复杂时，就可以利用这些内部状态, 线性地组合处所需要的对应输出~(实际上就是传统的MLP拟合的能力)。那么，到底是如何训练这个网络的呢？不要着急，我们先来对ESN 做一些数学符号的定义。

数学定义

假设这个回声状态网络有 $N$ 个中间节点，即储备池中的神经元个数为 $N$ ，输入层和输出层神经元个数都为 $D$ 。用 $u(t)\in \mathbb{R}^D,x(t)\in \mathbb{R}^N, f(t)\in \mathbb{R}^D$ 分别表示t时刻的输入，网络状态(储备池的状态)和输出; $V\in \mathbb{R}^{N\times D},R\in \mathbb{R}^{N\times N}, W\in \mathbb{R}^{D\times N}$ 分别表示输入权重、中间权值及输出权重矩阵， $tanh(\cdot )$ 为激活函数。注意，有些文献会把 $V$ 定义为 $W_{in}$ ， $W$
定义为 $W_{out}$ ， $R$ 定义为 $W$ 。
则储备池的状态更新方式以及网络的输出为：

\begin{matrix} (122) & x (t) = t a n h (R x (t - 1) + V u (t)) \end{matrix}

$\begin{equation} x(t)=tanh(Rx(t-1)+Vu(t)) \end{equation}$

\begin{matrix} (123) & f (t) = W x (t) \end{matrix}

$\begin{equation} f(t)=Wx(t) \end{equation}$

构造过程

首先我们来看下ESN的构造过程：初始化，训练以及使用(测试)，如下图所示。

回声状态网络的构造过程。

首先进行初始化操作，在这步里，我们先确定储备池的大小，即神经元的个数。与传统的MLP相同，节点数越多，拟合能力越强。由于ESN仅仅通过调整输出权值(即图1中的 $W$ ，我们最终就是利用储备池的状态信息来确定这个 $W$ )来线性拟合输出结果，所以一般ESN需要远大于一般神经网络的节点规模。
接下来是随机生成连接矩阵，这个矩阵表示了哪些神经元之间是有连接的，以及连接的方向和权值，实际上就是有向图的矩阵表示~(即图1中的 $R$ )。接下来的缩放矩阵实际上就是归一化的操作，有些时候我们会直接使用一个缩放因子，使用该缩放因子乘以原始的随机生成的矩阵，相比于使用特征值来缩放更加快速，但是很大程度上也丧失了精度。为什么要进行这个操作呢？原因和我们在初始化一些神经网络的权值时是类似的，对于这些神经网络，我们通常会将权重初始化为0-1之间(或着-1至1)，有两个原因：~(1)受激活函数的影响，比如图3所示的 $sigmoid$ 和 $tanh$ 激活函数，在0与1之间区分度比较大，但是大于1之后，激活值变化不大；~(2)我们对激活函数求导，可以看出在大于1时，图像比较平坦，其导数接近于0，由此在计算梯度时会导致梯度过小，无法顺利实现权重的更新。对于ESN我们不使用梯度来更新权值，主要是第一点的问题。最后随机生成输入权值 $V$ 和输出权值 $W$ 。
这些参数会影响到网络短期记忆时间的长短。输入权值越小而内部矩阵的谱半径越接近1，网络短期记忆时间越长。但是，增强记忆能力的同时，这种操作也造成了网络对“快速变化”系统的建模能力下降。

sigmoid(左)和tanh(右)

第二步进行训练。值得注意的是这个“空转”过程，实际上就是初始化储备池的状态。为什么要进行这个操作呢？因为储备池的内部连接是随机的，最开始的输入序列得到储备池状态的噪声会比较大，所以会先使用一些数据来初始化储备池的状态，从而降低噪声的影响。对于使用线性回归确定输出权值，在下一部分我们来一起进行推导。

数学推导

我们来对前面所说的输出权值进行求解。问题如下：这里假定对输出权重 $W$ 做了2范数正则，正则化系数为 $\lambda$ 。记网络状态矩阵为 $X$ ，输出序列矩阵
为 $Y$ ，请写出输出权重 $W$ 的计算公式。

首先我们要优化的目标是：

\begin{matrix} (124) & m i n {‖ W X - Y ‖}_{2}^{2} + λ {‖ W ‖}_{2}^{2} \end{matrix}

$\begin{equation} min\left \| WX-Y \right \|_2^2 + \lambda \left \| W \right \|_2^2 \end{equation}$
求解非常简单，只需要对其求导，令其导数为0，解出

W

$W$ 即可，实际上就是最小二乘法求解，可得：

\begin{matrix} (125) & W = Y X^{T} (X X^{T} + λ I)^{- 1} \end{matrix}

$\begin{equation} W=YX^T(XX^T+\lambda I)^{-1} \end{equation}$
进一步说，就是岭回归(带2范数惩罚项的最小二乘回归)。

考虑一个等价的问题：考虑概率情况， $y(t)=f(t)+\varepsilon$ ，其中 $\varepsilon \sim \boldsymbol{N}(0,\beta^{-1}I)$ ，
且对于 $W_i$ 的分布有 $W_i \sim \boldsymbol{N}(0,\alpha^{-1}I)$ 。证明最大化后验 $p(W_i|X,Y_i)$ 得到的输出权重与第1问中最小二乘法得到的输出权重等价。

最大化后验，即：

\begin{matrix} (126) & \begin{aligned} m a x P (W | X, Y) = m a x \frac{P (W, X, Y)}{P (X, Y)} \\ = m a x \frac{P (X | W, Y) P (W)}{P (X, Y)} \end{aligned} \end{matrix}

$\begin{equation} \begin{split} max\,\, P(W|X, Y)= max\,\, \frac{P(W,X,Y)}{P(X,Y)} \\ =max\,\, \frac{P(X|W,Y)P(W)}{P(X,Y)} \end{split} \end{equation}$
由于分母与

W

$W$ 无关，所以可以省略，由此上式等价于：

\begin{matrix} (127) & m a x P (X | W, Y) P (W) \end{matrix}

$\begin{equation} max\,\, P(X|W,Y)P(W) \end{equation}$
假设训练样本数为

L

$L$ ，则上式等价于：

\begin{matrix} (128) & m a x \prod_{i = 1}^{L} P (X_{i} | W_{i}, Y) P (W_{i}) \end{matrix}

$\begin{equation} max\,\, \prod_{i=1}^{L}P(X_i|W_i,Y)P(W_i) \end{equation}$
加上对数处理后，上式可变为：

\begin{matrix} (129) & m a x \sum_{i = 1}^{L} (l o g P (X_{i} | W_{i}, Y_{i}) + l o g P (W_{i})) \end{matrix}

$\begin{equation} max\,\, \sum_{i=1}^{L}(logP(X_i|W_i,Y_i)+logP(W_i)) \end{equation}$
由题目中的条件可知：

Y_{i} - W_{i} X_{i} \sim N (0, β^{- 1} I)

$Y_i-W_iX_i \sim \boldsymbol{N}(0,\beta^{-1}I)$ ，所以上式
可写作：

\begin{matrix} (130) & \begin{aligned} m a x \sum_{i = 1}^{L} (l o g \frac{1}{\sqrt{2 π β^{- 1} I}} e x p (- \frac{(Y_{i} - W_{i} X_{i})^{2}}{2 β^{- 1} I}) + \\ l o g \frac{1}{\sqrt{2 π α^{- 1} I}} e x p (- \frac{W_{i}^{2}}{2 α^{- 1} I})) \end{aligned} \end{matrix}

$\begin{equation} \begin{split} max\,\,\sum_{i=1}^{L}(log\frac{1}{\sqrt{2\pi \beta^{-1}I}}exp(-\frac{(Y_i-W_iX_i)^2}{2\beta^{-1}I })+\\ log\frac{1}{\sqrt{2\pi \alpha^{-1}I}}exp(-\frac{W_i^2}{2\alpha^{-1}I })) \end{split} \end{equation}$
展开上式，即等价于：

\begin{matrix} (131) & \begin{aligned} m a x \sum_{i = 1}^{L} (- \frac{(Y_{i} - W_{i} X_{i})^{2}}{2 β^{- 1} I} - \frac{W_{i}^{2}}{2 α^{- 1} I}) \end{aligned} \end{matrix}

$\begin{equation} \begin{split} max\,\,\sum_{i=1}^{L}(-\frac{(Y_i-W_iX_i)^2}{2\beta^{-1}I }-\frac{W_i^2}{2\alpha^{-1}I }) \end{split} \end{equation}$
将负号去掉，则上式变为：

\begin{matrix} (132) & \begin{aligned} m i n (\sum_{i = 1}^{L} \frac{(Y_{i} - W_{i} X_{i})^{2}}{2 β^{- 1} I} + \sum_{i = 1}^{L} \frac{W_{i}^{2}}{2 α^{- 1} I}) \end{aligned} \end{matrix}

$\begin{equation} \begin{split} min\,\,(\sum_{i=1}^{L}\frac{(Y_i-W_iX_i)^2}{2\beta^{-1}I }+\sum_{i=1}^{L}\frac{W_i^2}{2\alpha^{-1}I }) \end{split} \end{equation}$
使用矩阵的形式来表示，可得：

\begin{matrix} (133) & \begin{aligned} m i n \frac{{‖ W X - Y ‖}_{2}^{2}}{2 β^{- 1} I} + \frac{{‖ W ‖}_{2}^{2}}{2 α^{- 1} I} \end{aligned} \end{matrix}

$\begin{equation} \begin{split} min\,\,\frac{\left \| WX-Y \right \|^2_2}{2\beta^{-1}I }+\frac{\left \| W \right \|^2_2}{2\alpha^{-1}I } \end{split} \end{equation}$
进一步地等价于：

\begin{matrix} (134) & \begin{aligned} m i n {‖ W X - Y ‖}_{2}^{2} + \frac{α}{β} {‖ W ‖}_{2}^{2} \end{aligned} \end{matrix}

$\begin{equation} \begin{split} min\left \| WX-Y \right \|_2^2 + \frac{\alpha }{\beta } \left \| W \right \|_2^2 \end{split} \end{equation}$
我们令

λ = \frac{α}{β}

$\lambda=\frac{\alpha }{\beta }$ ，即可得到最终的式子：

\begin{matrix} (135) & \begin{aligned} m i n {‖ W X - Y ‖}_{2}^{2} + λ {‖ W ‖}_{2}^{2} \end{aligned} \end{matrix}

$\begin{equation} \begin{split} min\left \| WX-Y \right \|_2^2 + \lambda \left \| W \right \|_2^2 \end{split} \end{equation}$
哈，可以看到与公式1是相同的，由此得证。

Matlab代码实现

这里给出Matlab代码的实现，代码来源于http://minds.jacobs-university.de/mantas：

% A minimalistic Echo State Networks demo with Mackey-Glass (delay 17) data 
% in "plain" Matlab.
% by Mantas Lukosevicius 2012
% http://minds.jacobs-university.de/mantas

% load the data
trainLen = 3000;
testLen = 1000;
initLen = 100;

data = load('MackeyGlass-t17.txt');

% plot some of it
% figure(10);
% plot(data(1:1000));
% title('A sample of data');

% generate the ESN reservoir
inSize = 1; outSize = 1;
resSize = 1000;
a = 0.3; % leaking rate

%rand( 'seed', 42 );
Win = (rand(resSize,1+inSize)-0.5) .* 1;
W = rand(resSize,resSize)-0.5;
% Option 1 - direct scaling (quick&dirty, reservoir-specific):
% W = W .* 0.13;
% Option 2 - normalizing and setting spectral radius (correct, slower):
disp 'Computing spectral radius...';
opt.disp = 0;
rhoW = abs(eigs(W,1,'LM',opt));
disp 'done.'
W = W .* ( 1.25 /rhoW);

% allocated memory for the design (collected states) matrix
X = zeros(1+inSize+resSize,trainLen-initLen);
% set the corresponding target matrix directly
Yt = data(initLen+2:trainLen+1)';

% run the reservoir with the data and collect X
x = zeros(resSize,1);
for t = 1:trainLen
    u = data(t);
    x = (1-a)*x + a*tanh( Win*[1;u] + W*x );
    if t > initLen
        X(:,t-initLen) = [1;u;x];
    end
end

% train the output
reg = 1e-8;  % regularization coefficient
X_T = X';
% Wout = Yt*X_T * inv(X*X_T + reg*eye(1+inSize+resSize));
Wout = Yt*X_T / (X*X_T + reg*eye(1+inSize+resSize));
% Wout = Yt*pinv(X);

% run the trained ESN in a generative mode. no need to initialize here, 
% because x is initialized with training data and we continue from there.
Y = zeros(outSize,testLen);
u = data(trainLen+1);
for t = 1:testLen 
    x = (1-a)*x + a*tanh( Win*[1;u] + W*x );
    y = Wout*[1;u;x];
    Y(:,t) = y;
    % generative mode:
    u = y;
    % this would be a predictive mode:
    %u = data(trainLen+t+1);
end

errorLen = 1000;
mse = sum((data(trainLen+2:trainLen+errorLen+1)'-Y(1,1:errorLen)).^2)./errorLen;
disp( ['MSE = ', num2str( mse )] );

% plot some signals
figure(1);
plot( data(trainLen+2:trainLen+testLen+1), 'color', [0,0.75,0] );
hold on;
plot( Y', 'b' );
hold off;
axis tight;
title('Target and generated signals y(n) starting at n=0');
legend('Target signal', 'Free-running predicted signal');

figure(2);
plot( X(1:20,1:200)' );
title('Some reservoir activations x(n)');

figure(3);
bar( Wout' )
title('Output weights W^{out}');

后记

确定性跳跃循环状态网络（CRJ）是ESN的一个变种，如图所示：

确定性跳跃循环状态网络

在下一节我们会基于前面的ESN的代码一起来实现一个CRJ。