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

最新推荐文章于 2024-07-02 01:31:01 发布

comli_cn

最新推荐文章于 2024-07-02 01:31:01 发布

阅读量7.8k

点赞数 21

分类专栏：机器学习相关

本文链接：https://blog.csdn.net/comli_cn/article/details/109394553

版权

机器学习相关专栏收录该内容

11 篇文章 5 订阅

订阅专栏

1. ESN的任务

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

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

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

2. ESN的结构和训练步骤

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

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

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

整个计算过程如下所示：

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

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

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

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

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

在这里插入图片描述

3. ESN的预测步骤

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

$u(t)=W_{RO}*r(t)$

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

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

···

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

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

在这里插入图片描述

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

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

每个输入 $u (t)$ 都会和 $N / M$ 个库中的节点相连，因为输入个数时M，库中有N个节点，
即： $u\in R^{M*1}$ ， $r\in R^{N*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()