储备池计算（reservoir computing）是循环神经网络（recurrent neural network, RNN）的快速训练方法，优点在于可以避免梯度消失/爆炸的问题^[1]. 储备池计算的核心思想是使用固定、随机且稀疏连接的大规模储备池作为RNN的隐藏层，通过输入输出层之间的训练来学习解决特定的任务. 上述训练思想已在回声状态网络（echo state networks, ESN）中成功应用^[2-3]. ESN储备池内部的神经元存在循环连接（前一时刻与当前时刻的储备池之间存在连接权边），因此拥有记忆能力. ESN先随机生成输入连接权重和循环连接权重，再保持随机连接的循环网络不变以计算各时刻的内部状态，最后只需训练输出连接权重，即可得到可对网络内部状态进行分类/预测/回归的模型. 与传统的RNN相比，ESN具备更高的训练效率和便捷性，在各领域广受重视，尤其在时间序列预测任务上展现出独特优势^[4-5]. 关于ESN的研究仍存在瓶颈，如在储备池拓扑结构设计方面，传统储备池是随机网络，其不确定性导致性能波动较大. 如何让储备池具有自适应、自组织、自演化的复杂结构是新兴的研究课题.

自然界及人类社会中的大量复杂网络存在自适应、自组织、自演化的特点，例如生物神经网络^[6]、社交网络^[7]、交通网络^[8]等. 复杂网络^[9-10]由大量节点和连接构成，节点之间的连接关系呈现多样性和复杂性. 相关研究表明，具有无标度特性的复杂网络比随机网络更适合作为储备池的拓扑结构^[11]. 现有工作主要采用无标度网络的经典构造算法——BA（Barabasi-Albert）算法来生成储备池，该算法依赖于增长机制和偏好连接（preferential attachment, PA），具有幂律的度分布^[12]. 但BA算法也存在局限：一方面，该算法只能描述幂律指数γ=3的无标度网络；另一方面，网络在增长过程中，越旧的节点具有越高的度数，而新增的节点不可能出现度数大的情况. 基于BA算法的适应度模型由此诞生^[13]，该模型通过引入可调参数对偏好连接机制进行细化. Tsiotas^[14]提出适应度的时间动态概念，利用随网络增长而变化的拓扑适应度生成无标度网络，其中适应度由度、聚类系数、中间度、接近度和特征向量中心性控制. 为了更好解释社交网络的动态性，Behera等^[7]引入基于加权介数的模型，模型节点连接的标准是加权介数中心性而不是度数，介数中心性在网络增长的每一步都被计算，能更准确地描述各种真实社交网络. 还有方法通过修改偏好连接机制，使连接到定点的概率与该定点的度数非线性相关，这种机制被称为非线性偏好连接（nonlinear preferential attachment, NPA）^[15].

以上算法都是通过一个接一个添加节点的方式生成无标度网络，构建不灵活且速度慢. Wang等^[16-17]提出随机图构造算法，称为Chung-Lu（CL）算法，可以灵活简单地构建出复杂网络. Fasino等^[18]进一步证明通过适当选择参数，CL算法可以生成具有多种性质的无标度网络，如具有指定幂律指数的网络以及具有理想平均度和最大期望度的网络. 遗憾的是，尚未有研究工作将CL算法构造的无标度网络应用于ESN储备池设计. 本研究利用基于随机矩阵理论的CL算法生成具有理想平均度的无标度网络，并将其作为ESN储备池拓扑结构；针对CL算法在构造无标度网络过程中出现的度值偏差问题，引入随机剪枝和度值剪枝对网络连边进行剪枝.

传统的ESN包括输入层、大规模储备池和输出层，假设它们的神经元个数分别为I、N、O. ESN的特点是全连接的输入权重$ {\boldsymbol{W}}_{\text{in}}\text{}\in {\text{}\bf{R}}^{{N\times I}} $和稀疏连接的循环权重$ {{\boldsymbol{W}}}\text{}\in \text{}{\bf{R}}^{{N\times N}} $在随机生成后均保持不变，网络唯一要训练得到的权重是输出权重$ {\boldsymbol{W}}_{\text{out}}\text{}\in \text{}{\bf{R}}^{{O\times }\text{(}{N+I}\text{)}} $. ESN的训练过程包括3个步骤. 1）设有长度为l的训练输入输出对[u(t), y(t)]，1 ≤ t ≤ l，计算储备池内部状态$ \boldsymbol{x}\text{(}{t}\text{)}\text{}\in \text{}{\bf{R}}^{{N}} $：

(1)$ \boldsymbol{x}\left({t}\right)={f}\left({{W}}_{\text{in}}\boldsymbol{u}\left({t}\right)+\boldsymbol{Wx}\left({t}{-1}\right)+\boldsymbol{b}\right). $

式中：$ \boldsymbol{b}\text{}\in \text{}{\bf{R}}^{\text{N}} $为偏置向量，f (·)一般为非线性激活函数tanh，x(0)一般初始化为零向量或随机向量. 2）为了保证储备池具有“回声”属性，通常要舍弃输入序列的初始阶段. 假设初始阶段持续到s时刻，那么从第s+1时刻开始对内部状态进行收集. 对应的扩展状态矩阵X和输出矩阵Y分别为

(2)${\boldsymbol{ X}}=\left[\begin{array}{c}{[\boldsymbol{u}(s+1) ; \boldsymbol{x}(s+1)]^{\mathrm{T}}} \\\vdots \\{[\boldsymbol{u}(l) ; \boldsymbol{x}(l)]^{\mathrm{T}}}\end{array}\right] \in \mathbf{R}^{(l-s) \times(N+I)}， $

(3)$ {\boldsymbol{Y}}=\left[\begin{array}{c}{\boldsymbol{y}}(s+1)^{\mathrm{T}} \\\vdots \\{\boldsymbol{y}}(l)^{\mathrm{T}}\end{array}\right] \in \mathbf{R}^{(l-s) \times O} . $

式中：[·;·]为2个向量的上下拼接操作. 3）求解$ {\boldsymbol{W}}_{\text{out}} $，计算式为

(4)$ {\boldsymbol{W}}_{\text{out}}\text=\text{argmin}\,\,{||{{\boldsymbol{Y}}}-{{\boldsymbol{XW}}_{\text{out}}}^{\text{T}}||}_{2}^{2}\text{.} $

求解方法有伪逆法、岭回归法或弹性网络等^[19-21]. 在测试阶段，先根据测试数据的输入u(t)，利用式(1)计算内部状态x(t)，再利用训练得到的$ {{{\boldsymbol{W}}}}_{\text{out}} $计算测试数据的输出：

(5)$ \boldsymbol{y}(t)=f^{\text {out }}\left({\boldsymbol{W}}_{\text {out }}[\boldsymbol{u}(t) ; \boldsymbol{x}(t)]\right). $

式中：$ {{f}}^{\text{}\text{out}} $(·)为线性激活函数.

将ESN储备池抽象为无向加权图G = {V, E, A}；其中V为节点（神经元）的集合（$ \text{|}{V}\text{}\text{|}\text{}\text{=}\text{}{N} $）；E为相邻时刻神经元之间的连接边的集合（$ \text{|}{E}\text{|}\text{}\text{=}\text{}{M} $）；$ \boldsymbol{A}\text{}\in {\text{}\bf{R}}^{{N}\text{×}{N}} $为对称的邻接矩阵，矩阵元素表示边权重，范围在[0,1.0]. 若$ {{A}}_{{ij}} $= 0，则表示节点i和j不相连；若$ {{A}}_{{ij}} $> 0，则表示节点i和j相连，数值越接近1.0，节点之间的信息传递越强. $ {{d}}_{{i}} $为节点i的度是与该点相连的边的数量，$ \text{0}\text{}\leqslant{\text{}{d}}_{{i}}\text{}\leqslant\text{}{N}{-1} $；节点i的强度$s_i $是与它相连的边的权重之和，可由邻接矩阵A的行或列求和得到$ {{s}}_{{i}}\text{=}\displaystyle\sum {{A}}_{{ij}} $.

CL算法生成的图通常用G(w)表示，$ \boldsymbol{w}=[w_1, \cdots, w_i, \cdots, w_N ]^{\mathrm{T}} $为预先设定的理想度向量，$ {{w}}_{{i}} $为非负实数，表示节点i的理想度值，决定节点i与其他节点的连接程度. 对于节点i, j = 1,···, N,节点i和j存在边的概率为

(6)$ P_{i j}=\frac{w_i w_j}{\sigma}, \,\,\,\,\sigma=\sum_{k=1}^N w_k . $

与节点i相连的理想边数，即节点i的理想度值为

(7)$ w_i=\sum_{j=1}^N P_{i j}=\frac{w_i}{\sigma} \sum_{j=1}^N w_j . $

以上为生成N个节点的随机图$ {G}\text{}\in \text{}{G}\text{(}\boldsymbol{w}\text{)} $规则. 对于某些度向量，$ {{P}}_{{ij}} $可以大于1，为了保持式(6)的意义，通常对理想度向量w施加约束. 本研究采用定义1生成随机图.

定义1　当$ w_1 \geqslant \cdots \geqslant w_N \geqslant 0 $并且$ w_1^2 \leqslant \displaystyle{\sum}_{i=1}^N w_i $时，称$ \boldsymbol{w}\text{}\in {\text{}\bf{R}}^{{N}} $是满足CL算法生成基本条件的向量. B为所有满足条件的w的集合.

CL算法允许提前设定图G(w)的度分布，选择满足幂律度分布的参数w即可生成具有任意指数的幂律网络. 假设$ n_k \approx \alpha k^{-\gamma} $（α为常数）是有N个节点且$ \gamma > 1 $的无标度网络G的度分布，$ n_k $为度值为k的节点数. 那么度值大于或等于k的节点数N(k)可以近似表示为

(8)$ N(k)=\sum_{i=k}^{\infty} n_i \approx \alpha \int_k^{\infty} x^{-\gamma} \mathrm{d} x=\frac{\alpha}{\gamma-1} k^{1-\gamma} . $

根据定义1可知，G中节点的最大度值为$ {{w}}_{\text{1}} $，记作$ {{d}}_{\text{max}} $. 假设最大度值节点只有一个，即N($ {{d}}_{\text{max}} $)=1，结合式(8)可以得到

(9)$ d_{\max } \approx\left(\frac{\alpha}{\gamma-1}\right)^{\frac{1}{\gamma-1}}. $

同理可得G中节点的最小度值为$ {{w}}_{{N}} $，记作$ {{d}}_{\text{min}} $. 假设G中节点的度值都大于或等于$ {{d}}_{\text{min}} $，那么$N\left(d_{\min }\right)=N \approx \alpha d_{\min }^{1-\gamma} / (\gamma-1 )$，则

(10)$ d_{\min } \approx\left(\frac{\alpha}{N \times \gamma-N \times 1}\right)^{\frac{1}{\gamma-1}}. $

当$ {G}\text{}\in \text{}{G}\text{(}\boldsymbol{w}\text{)} $且w满足定义1时，节点1,2,···,i的理想度值大于或等于$ {{w}}_{{i}} $. 假设$ {N}\text{(}{{w}}_{{i}}\text{)}\text{}\text{≈}\text{}{i} $，由式(8)得到$ i \approx \alpha w_{i}^{1-\gamma} / (\gamma-1 )$，解得

(11)$ w_i \approx c_i^{-\frac{1}{\gamma-1}},\,\, c=\left(\frac{\gamma-1}{\alpha}\right)^{\frac{1}{\gamma-1}} . $

式(11)的假设存在一些问题. 对于节点i之后的其余节点，存在度值和节点i的度值相同的情况，因此对于$ {{w}}_{{i}} $的计算存在误差. 为此，引入参数$ {{i}}_{\text{0}} $来调节i，$ {{i}}_{\text{0}} $对于克服由式(11)引起的网络结构的限制很有用，该参数还会影响度值较小的节点的数量. 调整后公式为

(12)$ w_i \approx c\left(i_0+i\right)^{-\frac{1}{\gamma-1}}, \,\,i_0 > 1 . $

用式(12)生成的理想度向量w仍然存在问题，即$ \boldsymbol{w}\text{}\notin \text{}B $，不能确保w满足定义1中CL算法的生成基本条件. 这个条件可以通过选择适当的常数c来满足.

定义2　当$ \gamma $>1时，对于式(12)，当$0 \leqslant c \leqslant c_{\max }$时，使$ \boldsymbol{w}\in B $.

(13)$ c_{\max }=\left\{\begin{array}{ll}\dfrac{1}{p-1}, &1 < \gamma < 2 ; \\\lg \;(N+1), &\gamma=2 ; \\\dfrac{N^{1-p}-p}{1-p},& \gamma > 2 .\end{array}\right. $

式中：中间变量$ p=1 /(\gamma-1) $. 通过式(13)的约束，利用CL算法可以生成具有无标度性质的复杂网络. 为了使$ {{i}}_{\text{0}} $和c能够在适当的假设下，约定CL无标度网络的理想平均度值d和最大理想度m.

定义3　当$\gamma $> 2时，假设d = d(N)和m = m(N)是2个非递减函数，且$ 0 \leqslant d(N) \leqslant m(N) \leqslant N $，满足

(14)$ \lim _{N \rightarrow \infty} \frac{d(N)}{m(N)}=0 . $

存在常数$ \eta(0<\eta<1) $，使

(15)$ \eta N d(N) \geqslant m(N)^2. $

对于节点$i=1, \cdots, N; w_i=c\left(i_0+i\right)^{-p}; p=1 /(\gamma-1) $；当满足

(16)$ c=c(N)=(1-p) d(N) N^p， $

(17)$ i_0=i_0(N)=n\left(\frac{(1-p) d(N)}{m(N)}\right)^{\frac{1}{p}}-1. $

那么有以下结论：1）$ \boldsymbol{w}\in {B} $，符合CL算法的生成基本条件. 2）使用G(w)生成的图有指定$ \gamma $的度值分布，符合幂律性质，即对于任意一个理想度值$ k \geqslant w_k $，节点数$ n_k \approx \alpha k^{-\gamma} $. 3）使用G(w)生成的图的平均度值$ \operatorname{mean}(\boldsymbol{w})=\displaystyle{\sum}_{i=1}^N w_i / N $渐进于d(N)，即

(18)$ \lim _{N \rightarrow \infty} \frac{\operatorname{mean}(\boldsymbol{w})}{d(N)}=1. $

须注意式(15)必须满足，如果不满足该条件，则无法获得任意的度向量w. 如果d(N)的上限是常数，则该条件意味着最大理想度值m(N)的增长速度不得快于$ \sqrt{{N}} $. 定义3本质上是遵循幂律分布的w的最一般的形式. 考虑到现实中无标度网络中的平均度值都是大于1的，进一步增加约束使得CL算法满足该条件.

定义4　令$ \text{0}\text{}\text{ < }\text{}{p}\text{}\text{ < }\text{}\text{1} $，d为大于1的常数. 对于节点1,2,···, N；$ w_i=c\left(i_0+i\right)^{-p}; p=1 /(\gamma-1); \gamma>2 $；$ {{i}}_{\text{0}} $和c按式(16)和(17)选择. 为了生成满足定义1的度向量w，且图G的mean(w)渐进于d，最大理想度m须满足：

(19)$ m=\sqrt{\frac{d N}{2}}. $

基于以上4个定义，通过改进的CL算法可以生成满足无标度幂律性质的网络G(w). 如图1所示为改进CL算法生成的具有不同幂律指数的无标度网络的节点度分布. 该图显示了度值k与度值数量$ {{n}}_{{k}} $的关系，取N=10 000，d=150. 随着度值的增大，度值数量呈现指数级下降，形成尾部延伸，这与幂律分布特点相符. 随着幂律指数$ \gamma $增大，尾部延伸越来越长，进一步表明基于改进CL算法构造的复杂网络具有无标度性质.

本研究提出利用改进Chung-Lu算法生成具有理想平均度的无标度网络，并将其作为ESN储备池拓扑结构. 该储备池具有2个可调参数：幂律指数$ \gamma $和理想平均度值d. 幂律指数用于控制网络的度分布，较大的$ \gamma $表示存在高度连通的节点，较小的$ \gamma $表示节点度数分布更加均匀. 在现实网络中，$ \gamma $∈[2, 3]. 理想平均度值影响网络的真实平均度，当d较小时，节点之间的连接相对较少，可能导致网络稀疏；当d较大时，网络更加密集，节点之间的连接较多；但太大的d将导致网络缺少高度连通的中心节点，网络可能失去典型的无标度性质，影响网络的传播性能. 无标度储备池的生成步骤如下.

1）将ESN储备池抽象为无向加权图$ G = \{V, E, \boldsymbol{A}, \boldsymbol{w}\} $. $\boldsymbol{w}=\left[w_1, w_2, \cdots, w_N\right]^{\mathrm{T}} $为待求的理想度向量.

2）初始化幂律指数和理想平均度值，根据$ {p}{= 1/(}\gamma {-1)} $求得$ {p} $；根据式(16)和(19)初始化c和最大理想度m，根据c、m和(17)得到参数$ {{i}}_{\text{0}} $；由以上参数和式(12)初始化理想度向量$ \boldsymbol{w} $.

3）根据w计算储备池边的数量M. M=$\displaystyle {\sum}_{i=1}^N w_i / 2 $为储备池的理想边数，在生成过程中可能会产生重复边，为了抵消重复边，将储备池的理想边数设为M+e，$ e=\operatorname{mean}_2({\boldsymbol{w}})^2 / 2 $. $ \operatorname{mean}_2(\boldsymbol{w}) = \displaystyle{\sum}_{i=1}^N w_i^2 / \displaystyle{\sum}_{i=1}^N w_i $为w的二阶平均值. 这样生成的边数量不恒定.

4）计算理想度向量$ \boldsymbol{w} \in \mathbf{R}^N $的累加和，再把结果赋值给向量$ {\boldsymbol{w}}_{\text{sum}} $，并对$ {\boldsymbol{w}}_{\text{sum}} $进行归一化操作，归一化后的$ {\boldsymbol{w}}_{\text{sum}} $是递增的.

(20)$ \begin{split}& \boldsymbol{w}_{\text {sum}}=\left[w_{\text {sum}_0}, w_{\text {sum}_1}, \cdots, w_{\text {sum}_N}\right]^{\mathrm{T}} \in \mathbf{R}^{N+1}； \\& w_{\text {sum}_t}=w_{\text {sum}_{t-1}}+w_t, w_{\text {sum}_0}=0, \quad t\geqslant 1 .\end{split} $

5）生成储备池边的行列索引集合I和J. 先生成随机向量$ \boldsymbol{r}\text{}\in {\text{}\bf{R}}^{{M}\text{×1}} $，再根据随机值$ {{r}}_{{i}}\text{}\in \text{[0,}\text{}\text{1}\text{.0}\text{)} $和$ {\boldsymbol{w}}_{\text{sum}} $，将每个随机值$ {{r}}_{{i}} $分配到$ {\boldsymbol{w}}_{\text{sum}} $相应的区间，具体而言，当$ {{w}}_{{\text{sum}}_{{t}{-1}}}\text{}\text{ < }{\text{}{r}}_{{i}}\text{}\leqslant{\text{}{w}}_{{\text{sum}}_{{t}}} $，将时间步t加入行索引集合I中. 列索引集合J同理.

6）行索引集合和列索引集合的数量均为M个. 分别依次取出I和J中的元素，组成非零边的下标，最终得到具有M条边的稀疏矩阵A. 例如I={2,1,4}，J={3, 2, 6}，由于是无向图，那么A中$ {{A}}_{\text{23}} $和$ {{A}}_{\text{32}} $、$ {{A}}_{\text{12}} $和$ {{A}}_{\text{21}} $、$ {{A}}_{\text{46}} $和$ {{A}}_{\text{64}} $不为零.

7）在步骤3）中为了抵消重复边，额外增加了储备池的理想边数，会导致真实的平均度值与mean(w)有差距，且随着理想平均度值d的增大，差距会越来越大. 为了解决该问题，对稀疏矩阵A进行剪枝处理. 根据真实的平均度值与mean(w)的差值计算剪枝的总边数n，接着选择以下2种机制进行剪枝. 随机剪枝：随机选择稀疏矩阵A中的n条边，赋予$ {{A}}_{{ij}}\text{}\text{=}\text{}{{A}}_{{ji}}\text{}\text{=}\text{}\text{0} $. 度值剪枝：根据每个节点的度值计算每个节点的剪枝比例，按照比例得到每个节点须减掉的边数$ {{n}}_{{i}} $，再依次遍历度值大于mean(w)的节点，随机选择$ {{n}}_{{i}} $条边，赋予$ {{A}}_{{ij}\text{}}\text{=}\text{}{{A}}_{{ji}}\text{}\text{=}\text{}\text{0} $. 重复此步骤直到真实的平均度值与mean(w)的差值保持不变. 采取这2种剪枝机制的目的是模拟对网络的随机攻击或针对性攻击，以提升网络的鲁棒性. 在随机攻击模拟中，攻击是随机执行的，不考虑网络中各节点的重要性，这样的模拟可以反映当网络遭遇无差别的干扰时，网络稳定的抵抗能力. 在针对性攻击模拟中，攻击是有策略的，特别是通过识别并移除节点度大的节点，即被选中的是网络中作用最为关键的节点. 这种模拟能够解释网络在经受目标明确的攻击时的脆弱性. 通过这2种不同类型攻击的模拟，可以更全面地理解网络在面对各种挑战时的能力.

8）将矩阵A中的非零边赋予区间[−1, 1]中符合均匀分布的权重值.

基于上述步骤，引入剪枝机制的无标度ESN模型，训练过程见算法1.

算法1　无标度ESN模型的训练算法

输入：ESN模型参数、幂律指数$ \gamma $、理想平均度值d、比例因子$ \alpha_{\mathrm{W}} $、抛弃时刻数s

输出：训练完成的CL-ESN模型

1. 初始随机生成满足均匀分布条件的$ {{{\boldsymbol{W}}}}_{\text{in}} $.

2. 根据选定参数幂律指数$ \gamma $、理想平均度值d、$ p=1 /(\gamma-1) $，求得p.

3. 根据式(16)和(19)初始化c和最大理想度m，再根据c、m和式(17)得到参数$ {{i}}_{\text{0}} $. 最后通过式(12)初始化理想度向量w.

4. 根据w与M=$\displaystyle{\sum}_{i=1}^N w_i / 2 $计算储备池边的数量M. 将储备池的理想边数设为M+e，$e=\operatorname{mean}_2({\boldsymbol{w}})^2 / 2 $.

5. 根据式(20)计算理想度向量$ \boldsymbol{w}\text{}\in {\text{}\bf{R}}^{{N}} $的累加和，把结果赋值给向量$ {\boldsymbol{w}}_{\text{sum}} $并进行归一化操作.

6. 生成储备池边的行列索引集合I和J. 先生成随机向量$ \boldsymbol{r}\text{}\in {\text{}\bf{R}}^{{M}\text{×1}} $，再根据随机值$ {{r}}_{{i}}\text{}\in \text{[0,}\text{}\text{1.0)} $和$ {\boldsymbol{w}}_{\text{sum}} $，将每个随机值$ {{r}}_{{i}} $分配到$ {\boldsymbol{w}}_{\text{sum}} $相应的区$ {{w}}_{{\text{sum}}_{{t-}\text{1}}}\text{}{ < }{\text{}{r}}_{{i}}\text{}\leqslant{\text{}{w}}_{{\text{sum}}_{{t}}} $，将t加入行索引集合I中. 列索引集合J同理.

7. 对稀疏矩阵A，根据真实的平均度值与mean(w)的差值计算剪枝的总边数n，并进行剪枝.

8. 将A中的非零边赋予区间[−1,1]中符合均匀分布的权重值，生成具有无标度性质的初始储备池连接权重$ {\boldsymbol{W}}_{\text{0}}\text{}\text{=}\text{}\boldsymbol{A} $.

9. 通过比例因子$ \alpha_{\mathrm{W}}\left(0<\alpha_{\mathrm{W}}<1\right) $ 调整$ {\boldsymbol{W}}_{\text{0}} $，使其满足回声状态的必要条件，即$ \boldsymbol{W}=\left(\alpha_{\mathrm{W}} / \rho\left(\boldsymbol{W}_0\right)\right) \boldsymbol{W}_0 $，其中$ \rho\left({\boldsymbol{W}}_0\right) $为$ {\boldsymbol{W}}_{\text{0}} $的谱半径.

10. 设置储备池初始内部状态x(0) = 0，利用输入变量$ \boldsymbol{u}\text{(}{t}\text{)} $驱动储备池，根据式(1)计算内部状态$ \boldsymbol{x}\text{(}{t}\text{)} $.

11. 从s+1时刻开始收集储备池的内部状态，根据式(2)、(3)得到扩展状态矩阵X和输出矩阵Y.

12. 根据式(4)计算网络的输出权重矩阵$ {\boldsymbol{W}}_{\text{out}} $.

13. 训练结束.

为了测试所提算法构建的ESN储备池的有效性，采用混沌时间序列预测任务对模型性能进行测试，预测任务基于Mackey-Glass（MG）数据集^[22]和Multiple Superimposed Oscillator（MSO）数据集^[23]进行.

MG数据集是混沌时间序列预测模型的基准数据集，广泛用于测试神经网络对非线性混沌系统的辨识能力. 该数据集由时滞微分方程给出：

(21)$ \frac{\mathrm{d} u(t)}{\mathrm{d} t}=\frac{a u(t-\gamma)}{1+u^n(t-\gamma)}+b u(t) . $

当$ \gamma $>16.8时，MG数据集产生混沌吸引子，建模难度增加. 分别取$ \gamma $=16和$ \gamma $=17生成2个数据集（MG-16和MG-17），其他参数取值分别为n=10、a=0.2、b=−0.1，初始条件u(0)=1.2. 训练和测试样本数均为2 000. MG数据集的主要特点是具有长期记忆效应，即当前状态不仅受到前一时刻的状态影响，还受到更早时刻的状态影响. 此外，MG数据集的时间序列数据具有混沌性，即微小的扰动可能导致系统的行为发生不可预测的变化. 基于这些特点，MG数据集通常被用于测试非线性模型的性能，特别是在神经网络的时间序列预测方面.

某个正弦信号的预测任务相对简单，仅需有限的记忆容量就能预测网络的下一步输出. 不过预测不是整数倍频率叠加的正弦波是具有挑战性的任务，原因是正弦信号的波长可能非常长，导致预测困难. 为了解决这个问题，Fu等^[23]提出使用多重叠加振荡器（multiple superimposed oscillator, MSO）的定义：

(22)$ \begin{split} y(t_{\mathrm{mso}})& =\sin (0.2 t_{\mathrm{mso}})+\sin (0.311 t_{\mathrm{mso}})+\sin (0.42 t_{\mathrm{mso}})+ \\& \sin (0.51 t_{\mathrm{mso}})+\sin (0.63 t_{\mathrm{mso}})+\sin (0.74 t_{\mathrm{mso}}) . \\[-4pt]\end{split} $

式中：t_mso为MSO采样时间步. 在MSO数据集中，要求模型能够同时管理多个独立的内部状态. 传统的ESN模型中储备池神经元通常存在相互耦合现象，这与MSO数据集的要求相反. 为了更好处理MSO数据集，必须对ESN储备池的结构进行优化. 本节生成10 000个数据点，随后对MSO数据集进行标准化的预处理操作，以确保数据的一致性. 处理过后的数据如图2所示. 图中，A为振荡信号幅值.

选用均方误差MSE作为模型预测性能的评价指标，计算式为

(23)$ \mathrm{MSE}=\frac{1}{n} \sum_{i=1}^n\left(y_i-y_i^{\prime}\right)^2. $

式中：$ {{y}}_{{i}} $为真实值，$ {{y}}_{{i}}^{{'}} $为预测值. MSE越小，表示模型预测的结果越接近真实值，模型预测性能越好.

富人俱乐部系数RC用来评估网络中高度连通节点的紧密程度. 该参数通常基于节点的度来衡量节点的重要性，一般而言，富人俱乐部的定义主要关注节点的拓扑结构，不考虑权重：

(24)$ \mathrm{RC}=\frac{2 M_k}{N_k\left(N_k-1\right)} . $

式中：$ {{M}}_{{k}} $为连接度大于k的节点之间的连接数量，代表着富人俱乐部内部的连接数；$ {{N}}_{{k}} $为连接度大于k的节点的数量，代表富人俱乐部内部的节点数；$ {{N}}_{{k}} $−1是为了排除节点自身的连接. RC∈[0, 1.0]，当RC较高时，意味着度大于k的节点子群内部存在较高比例的连接，形成了紧密的子网络；RC较低时，节点子群内部缺乏连接，子网络更加松散随机.

平均绝对百分比误差MAPE的计算式为

(25)$ \text { MAPE }=\frac{100 {\text{%}}}{n} \sum_{i=1}^n\left|\frac{y_i-y_i^{\prime}}{y_i^{\prime}}\right| . $

MAPE越小，表示模型的准确性较高.

运行时间t_r是指从构造储备池开始到训练出输出连接矩阵的总时间.

基于前期工作基础，实验设置储备池N=500，CL算法参数幂律指数γ=2.3、2.6、2.9，理想平均度数的范围在[10, 200]，步长为10. 考虑到运行的随机性，每组参数迭代运行10次取均值.

如图3所示为引入度值剪枝或随机剪枝后CL算法的节点度分布. 观察可得到与图1相同的结论，证明引入度值剪枝或随机剪枝的CL算法仍具有无标度性质.

如图4所示为当$ \gamma $=2.6时，有无剪枝算法下理想平均度值d与度值k的关系. mean(w)为CL算法G(w)中的理想度向量的一阶均值，即$ \operatorname{mean}(\boldsymbol{w})= \displaystyle{\sum}_{i=1}^N w_i / N $，真实平均度值是储备池邻接矩阵A中非零值的个数与节点数N的比值. 观察图4(a)可以发现，随着d的增大，mean(w)和真实平均度值都在增大；虚线一直位于实线上方，意味着真实生成的网络平均度值一直大于mean(w)，而且两者之间的差距越来越大. 原因是在生成储备池过程中，为了抵消重复边，增加了储备池的理想边数，并且理想边数会随着mean(w)的增加而增加. 引入2种剪枝方法后，mean(w)与真实平均度值的差距不再增大，并且通过随机剪枝后的CL算法效果更好. 原因是度值剪枝按照比例确定每个节点的剪枝数，该过程会造成一定的误差；随机剪枝只涉及减法，造成的误差比度值剪枝小.

如图5所示为有无剪枝机制的算法关于不同幂律指数的理想平均度值d和富人俱乐部系数RC之间的关系，设置的富人俱乐部的阈值度为100. 由图可知，在3种幂律指数下，随着理想平均度值的增加，RC开始为0，然后整体呈上升趋势，但在局部出现一些波动. 这表明节点更倾向于与度值大的节点连接，使度值大的节点在网络中具有更重要的地位，形成不平等的连接，导致富者愈富的现象，即高度连接的节点越来越集中在少数节点上. 这种现象与无标度网络的特征是一致的. 另外，$ \gamma $越大，RC越大，原因是$ \gamma $增大，网络节点间形成了更多的边，导致RC更大.

在$ \gamma $=2.6的情况下，生成d=50、100、150、200符合无标度网络性质的理想度向量w，对比4种理想度向量下的储备池的稀疏度$ {{R}}_{\text{c}} $如表1所示，3种机制下基于CL算法构建的ESN的MSE、MAPE和t如表2~4所示，结果为10次运算的均值. 表中，下标0代表无剪枝机制，下标1代表度值剪枝机制，下标2代表随机剪枝机制. 由表1~3可知，1）在3种机制下，引入度值剪枝或随机剪枝的CL算法比无剪枝的CL算法具有更小的$ {{R}}_{\text{c}} $、MSE和MAPE. 这说明剪枝处理没有对网络的性能造成影响，反而有助于维持网络的稳定性和鲁棒性. 通过剪枝机制，网络在去除一些冗余连接的同时，仍保留关键连边，在面对各种攻击时展现出更强的抵抗能力. 2）具有随机剪枝的CL算法比度值剪枝的CL算法具有更小的$ {{R}}_{\text{c}} $、MSE和MAPE. 由于随机剪枝涉及非关键性节点，对网络结构的影响较小，这表明网络在一定程度上能够抵御随机攻击，并在失去部分边后仍然保持其核心功能. 由于度值剪枝涉及的是关键性节点，当这些节点受到攻击时，对网络结构的影响较大，从而对网络的性能造成影响. 由表4可知，1）在不同数据集下，3种机制CL算法的运行时间都随着d的增大而增大，说明d越大，网络越复杂. 结合表1可得，d增大，网络的稀疏度越大，网络边数增多，网络越复杂. 2）在不同数据集下，引入度值剪枝或随机剪枝的CL算法比无剪枝的CL算法的运行时间长，说明剪枝操作会引入额外的时间复杂度. 引入度值剪枝或随机剪枝后的算法，须对多的边数进行剪枝，步骤更加繁琐，造成运行时间延长，但最长不超过1.5 s. 3）在不同数据集下，引入随机剪枝的CL算法比引入度值剪枝的CL算法运行速度快. 原因是随机剪枝只需根据多的边数进行一次剪枝；度值剪枝按照比例确定每个节点的剪枝数，该过程会造成误差，须通过迭代来减小误差. 因此度值剪枝比随机剪枝的时间复杂度更高，运行时间更长，最大差异在1 s左右.

根据无标度性质的对比实验以及预测性能的对比实验可知，基于随机剪枝的CL算法优于无剪枝或度值剪枝的CL算法，为此选用所提算法与BA算法及其各种变体生成不同的储备池并得到相应的ESN，对比它们对混沌时间序列的预测精度，平均绝对百分比误差和运行时间. 对比实验算法为BA算法^[12]，MC^[24]，MCC^[25]，MCE^[26]，已经被证明可构建无标度网络. BA算法是在度值k控制下生成的适应算法；MC为在聚类系数控制下生成的适应算法；MCC为在紧密中心性（closeness centrality, CC）控制下生成的适应算法，紧密中心性用于衡量某节点与其他节点之间的紧密程度，表示该节点到其他节点的平均最短路径长度的倒数；MCE是在特征向量中心性（eigenvector centrality, CE）控制下生成的适应算法，特征向量中心性是基于节点在网络中的连接程度和邻居节点的中心性来计算的，简单来讲就是相邻的邻居节点越重要，该节点越重要. 如表5所示，设置参数$ {N}\text{=500}、{d}\text{=150}、 \gamma = 2.6 $，结果为运行10次的均值. 由表可知，1）随机剪枝的CL算法比对比算法的预测精度高. 说明本研究所设计的储备池结构在性能上优于对比算法构建的储备池，具有更好的数据拟合和处理能力. 2）4种对比算法的运行时间比所提算法运行时间长，其中BA算法运行时间较短，其他3种运行时间与所提算法差距很大. 主要原因是BA算法及其变体都依赖增长和偏好连接，每加入一个新的节点，都要与指定个数的旧节点进行连接；并且在选点过程中须计算图的k、聚类系数、CC或者CE，耗费时间成本很大. 该实验结果进一步说明，所提算法构建灵活且速度快.

脑科学的研究表明人脑可能具有内禀的网络特性，即无标度特性，它使得大脑内的神经元有枢纽和层级感. 为了模拟大脑网络的构造过程，本研究提出基于随机矩阵理论的Chung-Lu构造算法生成无标度网络，并将其作为ESN的储备池. 为了解决度值偏差，引入度值剪枝和随机剪枝2种机制删除冗余边，提升网络的鲁棒性. 度值剪枝有可能破坏网络中节点度大的节点，性能不如随机剪枝. 实验表明，基于剪枝机制的CL算法构建的网络，仍能保持无标度网络的幂律特性. 相比于BA算法及其变体，所提算法不再依赖于增长和偏好连接，因此构建灵活且速度快，并且在混沌时间序列预测方面具有更好的效果. 在未来的工作中，将着眼于利用不同类型的时间序列数据集，对提出的 ESN模型进行更全面的验证，以检验其在各种实际应用背景下的鲁棒性和有效性.