Parameter Optimization of Leaky Integrator Echo State Network with Internal-point Penalty Function Method
-
摘要: 为了提升泄露积分型回声状态网（Leaky integrator echo state network，Leaky-ESN）的性能，提出利用罚函数内点法优化Leaky-ESN的全局参数，如泄漏率、内部连接权矩阵谱半径、输入比例因子等，这克服了通过反复试验法选取参数值而降低了Leaky-ESN模型的优越性和性能.Leaky-ESN的全局参数必须保障回声状态网满足回声状态特性，因此它们之间存在不等式约束条件.有学者提出利用随机梯度下降法来优化内部连接权矩阵谱半径、输入比例因子、泄露率三个全局参数，一定程度上提高了Leaky-ESN的逼近精度.然而，随机梯度下降法是解决无约束优化问题的基本算法，在利用随机梯度下降法优化参数时，没有考虑参数必须满足回声特性的约束条件（不等式约束条件），致使得到的参数值不是最优解.由于罚函数内点法可以求解具有不等式约束的最优化问题，应用范围广，收敛速度较快，具有很强的全局寻优能力.因此，本文提出利用罚函数内点法优化Leaky-ESN的全局参数，并以时间序列预测为例，检验优化后的Leaky-ESN的预测性能，仿真结果表明了本文提出方法的有效性.Abstract: To improve leaky integrator echo state network (Leaky-ESN) performance, internal-point penalty function (IPF) method is used to optimize the global parameters of Leaky-ESN, such as leakage rate, spectral radius of internal connection weight matrix, scaling of input, etc., which overcomes loss of superiority and performance of Leaky-ESN because of using trial and error method to select parameter values. The global parameters of Leaky-ESN have to guarantee the echo state network to meet the echo state property, thus inequality constraints exist between them. Some researchers put forward the method using the stochastic gradient descent (GD) to optimize leakage rate, spectral radius of internal connection weight matrix, and scaling of input, which can improve the approximation precision of the Leaky-ESN to some certain extent. However, the stochastic gradient descent method is a basic algorithm to solve unconstrained optimization problems. Without considering parameters which need satisfy the constraint conditions of the echo state property (inequality constraints) during using stochastic gradient descent method, the parameter value is not the optimal solution. Internal-point penalty function method can solve the optimized problem with inequality constraints, a wide scope of application, fast convergence speed, strong ability of global optimization. Therefore, in this paper, internal-point penalty function method is used to optimize the global parameters of Leaky-ESN, and time series prediction is selected as an example to examine the performance of the optimized Leaky-ESN. Simulation results show the effectiveness of the proposed approach.
-
回声状态网(Echo state network, ESN)是由Jaeger教授于2001年提出的一种新型递归神经网络模型[1-2], 不仅能够解决递归神经网络存在的问题, 而且利用其预测著名的Mackey-Glass混沌时间序列, 与以往报道的预测结果相比, 预测精度提高2400倍, 成为ESN研究的标志性成果[3].回声状态网在系统辨识[4]、动态模式识别[5]、滤波或者控制[6-8]、预测[9-10]等领域获得成功应用.与传统的递归神经网络相比, ESN具有如下显著特点: 1) ESN由随机稀疏连接的神经元组成储备池作为隐层.储备池的生成过程独立于回声状态网络的训练过程, 保证了ESN在训练过程中的稳定性. 2) ESN采用线性回归或最小二乘法训练储备池至输出层的权值, 使网络的训练过程得以简化, 保证了权值的全局最优性.然而, 除了输出权值外, ESN的全局参数还包括储备池处理单元个数、内部连接权矩阵谱半径、输入比例因子、输出反馈比例因子、泄露率等.这些全局参数之间存在复杂的非线性关系, 通常在训练之前根据经验选取, 因此使得ESN的逼近能力不能达到最优.文献[11]提出了一种ESN的改进模型, 称为泄露积分型回声状态网(Leaky integrator echo state network, Leaky-ESN).Leaky-ESN利用随机梯度下降法来优化内部连接权矩阵谱半径、输入比例因子、泄露率三个参数, 一定程度上提高了Leaky-ESN的预测精度.众所周知, 回声状态网需要满足回声特性.然而, 随机梯度下降法是解决无约束优化问题的基本算法, 在利用随机梯度下降法优化参数时, 没有考虑参数必须满足回声特性的约束条件(不等式约束条件).因此, 从某种意义上说, 得到的参数值不是最优解.罚函数内点法可以求解具有不等式约束的最优化问题, 应用范围广, 收敛速度较快, 具有很强的全局寻优能力.因此, 本文以Leaky-ESN为例, 提出利用罚函数内点法优化回声状态网的全局参数, 进一步提升回声状态网的预测性能.利用罚函数内点法实质上是将不等式约束条件下的优化问题转换为无约束的优化问题, 进而通过牛顿法来求解无约束优化问题, 获得回声状态网的全局最优参数.
1. 泄露积分型回声状态网
Leaky-ESN网络是ESN网络的一种改进模型, 能够学习慢变动态系统等.Leaky-ESN与ESN的拓扑结构相同, 其储备池是由泄露积分型神经元组成, 如图 1所示[11].在图 1中, 输入层有 $K$ 个输入节点; 储备池由 $N$ 个内部节点以及稀疏的节点连接权值构成; 输出层有 $L$ 个输出节点.图中实线表示了网络的必要连接, 而虚线表示了不同情况下可能存在的连接.假定网络的输入 $u(n)=[u_1(n), u_2(n), \cdots, u_K(n)]^{\rm T}$ , 储备池的神经单元状态为 $x(n)=[x_1(n), x_2(n), \cdots, x_N(n)]^{\rm T}$ , 网络的输出 $y(n)=[y_1(n), y_2(n), \cdots, y_L(n)]^{\rm T}$ .那么, Leaky-ESN网络的状态更新方程和输出方程为
$\begin{align} &x\left( n+1 \right)=\left( 1-a \right)x\left( n \right)+f\left( {{s}^{\text{in}}}{{W}^{\text{in}}}u\left( n+1 \right)+ \right. \\ &\quad \quad \quad \quad \left. \left( \rho W \right)x\left( n \right)+{{W}^{\text{fb}}}y\left( n \right) \right) \\ \end{align}$
(1) $y(n)=g\left( {{W}^{\text{out}}}\left[ x(n);u(n) \right] \right)$
(2) 其中, $n$ 是离散时刻, $a$ 是泄漏率, $a\in\left[{0, 1}\right]$ ; $\rho $ 是储备池连接权值 $W$ 的谱半径; $s^{\rm in}$ 是归一化之后输入单元的比例因子, 称为输入比例因子; $f$ 是储备池单元的激励函数(通常是tanh函数或者是S型函数); $g$ 表示输出激励函数(典型的是单位函数和tanh函数), 根据问题的不同 $g$ 可以取线性函数或者S型函数; $\left[; \right]$ 表示二个向量级联; $W^{\rm in}$ 、 $W$ 、 $W^{\rm out}$ 和 $W^{\rm fb}$ 分别表示输入、储备池、输出、输出反馈的连接权矩阵.在整个学习和训练过程中, $W$ , $W^{\rm in}$ 和 $W^{\rm fb}$ 在随机生成之后就固定不变, $W^{\rm out}$ 是通过网络训练计算得到.Leaky-ESN模型的训练过程就是计算输出权值 $W^{\rm out}$ 的过程, 就是在均方根误差的意义上的最小化.通常采用Tikhonov正则化(岭回归)的方法获得 $W^{\rm out}$ [12], 即:
${{W}^{\text{out}}}=Y{{X}^{\text{T}}}{{\left( X{{X}^{\text{T}}}+\theta I \right)}^{-1}}$
(3) 式中, $\theta $ 是一正则化系数.假设 $u_{\rm teach} ( n )$ 是训练样本的输入信号, $y_{\rm teach} ( n )$ 是训练样本的期望输出信号.把所有的 $u_{\rm teach} \left( n\right)$ 和 $x\left( n \right)$ 储存起来并用矩阵 $X\in {\bf R}^{(K+N)× T}$ 的行来表示, 对应的目标输出值 $y\left( n\right)$ 表示成 $Y\in {\bf R}^{L× T}$ 的行向量.为了评价Leaky-ESN模型训练效果, 通常采用标准均方根误差 $\left({NRMSE} \right)$ , 其表达式为:
$\text{NRMSE}\left( \text{y},{{\text{y}}_{\text{teach}}} \right)\text{=}\sqrt{\frac{\left\langle {{\left\| \text{y}\left( \text{n} \right)\text{-}{{\text{y}}_{\text{teach}}}\left( \text{n} \right) \right\|}^{\text{2}}} \right\rangle }{\left\langle {{\left\| {{\text{y}}_{\text{teach}}}\left( \text{n} \right)\text{-}\langle {{\text{y}}_{\text{teach}}}\left( \text{n} \right)\rangle \right\|}^{\text{2}}} \right\rangle }}$
(4) 式中, $\left\| \cdot \right\|$ 表示的是欧氏距离; $\left\langle \cdot\right\rangle$ 表示的是均值函数.
2. 罚函数内点法在Leaky-ESN全局参数优化中的应用
2.1 全局参数优化问题的数学模型
在Leaky-ESN的训练过程中, 本文除了优化 $W^{\rm out}$ , 还优化泄漏率 $a$ , 储备池连接权矩阵的谱半径 $\rho$ , 输入比例因子 $s^{\rm in}$ .因此, 需要建立这些全局参数优化问题的数学模型.由于这些参数必须保障Leaky-ESN具有回声状态特性, 它们之间存在不等式约束.因此, 具有不等式约束的优化问题表示为:
\begin{align}\min \quad\quad &\overline f \left( \zeta \right) \nonumber\\ {\rm s.t.}\quad\quad &\bar {g}_i (\zeta )\ge 0, \quad\quad i=1, 2, 3, \cdots, m \nonumber\\ &\zeta \in \overline X \label{eq5}\end{align}
(5) 式中 $\overline f $ , $\overline {g}_i $ 是定义在 ${\bf R}$ 上的函数, $\overline f \left( \zeta \right)$ 是目标函数, $\bar{g}_i \left( \zeta \right)\ge 0$ 是不等式约束条件; $\zeta =\left({\zeta _1, \zeta _2, \cdots, \zeta _{\overline n } }\right)$ 是一个 $\overline n $ 维向量; 集合约束 $\overline X $ 是 ${\bf R}$ 上的一个子集.其可行域记为
\begin{align}S=\left\{ {\zeta \vert \overline g _i \left( \zeta \right)\ge0, i=1, 2, 3, \cdots, m, \zeta \in \overline X } \right\} \nonumber\end{align}
罚函数法就是利用目标函数和约束函数构造具有惩罚性质的函数 $P\left( \zeta\right)=\overline P \left( {\overline f (\zeta ), \bar {g}_i (\zeta )}\right)$ , 使原约束优化问题转化为求 $P\left( \zeta\right)$ 最优解的无约束优化问题.罚函数内点法, 又称内罚函数法, 它是一类保持严格可行性的方法, 在求解无约束优化问题时, 严格要求迭代点在可行域的内部移动.当迭代点由可行域的内部接近可行域的边界时, 将有无穷大的障碍, 迫使迭代点返回可行域的内部[12].
2.2 目标函数和约束条件的建立
Leaky-ESN全局参数优化的目的是希望网络输出值与期望值之间尽可能的接近.因此, 在本文中建立了如下的目标函数:
\begin{align}L\left( n \right)=\frac{1}{2}\left\| {y_{\rm teach} \left( n\right)-y\left( n \right)} \right\|^2 \label{eq6}\end{align}
(6) 只有在满足回声状态特性的情况下, 回声状态网才是有效的, 能够具有很好的性能, 比如对时间序列的预测性能.回声状态特性就是随着迭代次数的增加, 回声状态网受到初始条件的影响以一定的速率越来越小直至消失, 对于输入驱动来说, 就是不受之前输入和内部状态的影响, 网络达到渐近稳定状态.递归神经网络的大多数应用是纯输入驱动情况(例如动态模式识别、时间序列预测、滤波或控制).对于这些应用, 回声状态网模型通常忽略输出反馈, 即 $W^{\rm fb}=0$ .
下面给出文献[11]的回声状态特性充分条件(定理1) 和必要条件(定理2).
定理1 [11]. 对于Leaky-ESN模型(1) 和(2), 如果满足如下条件:
1) $f$ 是tanh函数;
2) $g$ 是有界函数(比如tanh函数)或没有输出反馈( $W^{\rm fb}=0$ );
3) $\left| {1-(a-\sigma (W))} \right|＜1$ , $\sigma (\cdot)$ 表示矩阵的最大奇异值;
那么Leaky-ESN模型对所有输入都满足回声状态特性.
定理2 [11]. 对于Leaky-ESN模型(1) 和(2), 如果满足如下条件:
1) $f$ 是tanh函数;
2) 矩阵的谱半径 $\left| λ \right|_{\rm max} (\overline W )>1$ , 其中 $\overline W =\left( {\rho W} \right)+\left( {1-a} \right)I$ , ( $I$ 是单位矩阵);
那么Leaky-ESN模型不具有回声状态特性.
由文献[11]的定理1和定理2可知, 在Leaky-ESN不考虑输出反馈情况下(即 $W^{\rm fb}=0) $ , 如果矩阵 $\overline W =\left( {\rho W} \right)+\left( {1-a} \right)I$ 的谱半径 $\left|λ \right|_{\rm max} (\overline W)$ 大于1, 那么Leaky-ESN模型就不存在回声状态特性, 换句话说就是必须小于等于1才行. $\left| λ \right|_{\rm max} (\overline W )\le1$ 被称为有效谱半径, 能够保证对所有的任务都满足回声状态特性.因此根据 $\left| λ \right|_{\rm max} (\overline W )\le1$ 可得约束条件 $1-a+\rho \le 1$ .因此, Leaky-ESN全局参数优化问题的约束条件为
$a-\rho \text{ }\ge 0$
(7) $1-a>0$
(8) $\rho >0$
(9) 2.3 罚函数内点法的障碍函数建立
罚函数内点法需要通过构建障碍函数来解决约束优化问题.所谓障碍函数, 就是在可行域边界设置一道屏障, 以保证算法在进行优化时始终在可行域内.罚函数内点法的障碍函数形式如下:
\begin{equation}\label{eq10}\varphi \left( {\zeta, M^{\left( k \right)}} \right)=\overline f \left(\zeta \right)+M^{\left( k \right)}\sum\limits_{i=1}^m {G\left[{\bar {g}_i\left( \zeta \right)} \right]}\end{equation}
(10) 式中, $k=1, 2, \cdots, l$ 是内点法惩罚因子的递归次数; $\zeta $ 是集合约束; $m$ 是集合约束中不等式条件的个数; $M^{(k)}$ 是第 $k$ 步惩罚因子( $M^{(k)}>0) $ .惩罚函数 $G[\bar {g}_i(\zeta )]$ 通常具有如下两种形式:
$G\left[ {{{\bar{g}}}_{i}}(\zeta ) \right]=\frac{1}{{{{\bar{g}}}_{i}}(\zeta )}$
(11) $G\left[ {{{\bar{g}}}_{i}}(\zeta ) \right]=-\ln \left( {{{\bar{g}}}_{i}}(\zeta ) \right)$
(12) 本文采用倒数形式的惩罚项, 即式(11).因此, 构造的内点法障碍函数形式如下:
$\varphi \left( \zeta ,{{M}^{\left( k \right)}} \right)=\bar{f}\left( \zeta \right)+{{M}^{\left( k \right)}}\sum\limits_{i=1}^{m}{\frac{1}{{{{\bar{g}}}_{i}}(\zeta )}}$
(13) 式中 $\zeta=[a; \rho], m=3.$ 在优化过程中 $M^{(k)}$ 是一个严格递减序列, $M^{(0) }>M^{(1) }>\cdots>M^{(k)}$ .因此在实际操作过程中令
${{M}^{(k+1)}}=c{{M}^{(k)}}$
(14) 式中 $c$ 为缩减系数 $(0＜c＜1) $ .初始的惩罚因子 $M^{(0) }$ 对于整个障碍函数的性能起着很大的影响作用.如果太小, 那么它起的作用就会很小; 如果太大, 那么需要花费大量的时间来得到最优值, 因此采用如下的形式来计算:
$M^{\left( 0 \right)}=\left| {\frac{\overline f \left( {\zeta \left( 0\right)} \right)}{\sum\limits_{i=1}^m {\frac{1}{\bar {g}_i \left( {\zeta\left( 0 \right)} \right)}} }} \right|$
(15) 式中, $\overline f \left( {\zeta \left( 0 \right)}\right)$ 是目标函数初始值; $\zeta \left( 0 \right)$ 是参数的初始优化值.
通过式(13) 可知, 随着迭代次数 $k$ 的增加, 在障碍函数中, 惩罚因子越来越小, 因此障碍项也越来越小, 换句话说就是目标函数越来越接近最优值, 那么此时得到的解越来越接近最优解.
Leaky-ESN模型的惩罚项可以表示成:
$G[g_i \left( {a, \rho } \right)]=\left\{ {a, \rho \vert a-\rho \ge0, 1-a>0, \rho >0} \right\}$
(16) 因此, Leaky-ESN模型的障碍函数为:
\begin{align}\label{eq17} \varphi \left( {X, M^{\left( k \right)}}\right)=&\sum\limits_{n=1}^T {L\left( n \right)}+\nonumber\\&M^{\left( k \right)}\left( {\frac{1}{a-\rho}+\frac{1}{1-a}+\frac{1}{\rho }} \right)\end{align}
(17) 式中, $T$ 是批量计算中一个样本点所含的离散时刻数.
通过式(17), 我们把非线性不等式约束优化问题转换成了无约束优化问题.
2.4 罚函数内点法的优化过程及步骤
本文采用牛顿法来求解无约束优化问题的最优参数值.牛顿法的基本思想是利用二次近似多项式的极值点求法给出原函数的极值点的求法.其迭代公式如下[13]:
\begin{equation}\label{eq18}\zeta (n+1) =\zeta (n)-\left[{\nabla ^2\varphi (\zeta (n))}\right]^{-1}\nabla \varphi (\zeta (n))\end{equation}
(18) 根据式(18) 可知, 需要计算泄漏率 $a$ 和谱半径 $\rho$ 关于障碍函数的一阶偏导数和二阶偏导数.为了简化符号的表达, 我们用 $\varphi_1 $ , $\varphi _2 $ , $\varphi _{11} $ , $\varphi _{12} $ , $\varphi_{21} $ 和 $\varphi _{22} $ 来代替 $\frac{\partial \varphi }{\partial a}$ , $\frac{\partial \varphi }{\partial \rho }$ , $\frac{\partial\varphi ^2}{\partial a^2}$ , $\frac{\partial \varphi ^2}{\partial a\partial \rho }$ , $\frac{\partial \varphi ^2}{\partial \rho\partial a}$ 和 $\frac{\partial \varphi ^2}{\partial \rho^2}(\partial \varphi =\partial \varphi (\zeta, M^{(k)}))$ , 用 ${\rm0}_u =(0, \cdots, 0) ^{\rm T}$ 来表示 $\frac{\partial u(n)}{\partial q}(q\in \left\{ {a, \rho } \right\})$ .令 $\ell \left( n\right)={y_{\rm teach} \left( n \right)-y\left( n \right)}$ , 所有偏导数的推导如下所示:
\begin{align}\varphi _1 \left( n \right)=&-\ell ^{\rm T}\left( n\right)W^{\rm out}\left[{\frac{\partial x\left( n\right)}{\partial a};{\rm 0}_u }\right]-\nonumber\\&M^{(k)}\left( {\frac{1}{\left( {\rho -a} \right)^2}-\frac{1}{\left({a-1} \right)^2}} \right)\label{eq19}\end{align}
(19) \begin{align}\varphi _2 (n)=&-\ell ^{\rm T}(n)W^{\rm out}[\frac{\partial x(n)}{\partial \rho };{\rm 0}_u]+\nonumber\\&M^{(k)}\left( {\frac{1}{(\rho -a)^2}-\frac{1}{\rho^2}}\right)\label{eq20}\end{align}
(20) \begin{align}\varphi _{11} (n)=&\left( {W^{\rm out}\bigg[\frac{\partial x(n)}{\partial a};{\rm 0}_u \bigg]} \right)^2-\nonumber\\&\ell^{\rm T} (n)W^{\rm out}\bigg[\frac{\partial ^2x(n)}{\partial a^2};{\rm 0}_u \bigg]-\nonumber\\&M^{(k)}\left( {\frac{2}{(\rho -a)^3}+\frac{1}{(a-1) ^3}}\right)\label{eq21}\end{align}
(21) \begin{align}\varphi _{12} (n)\!=\!&\left( W^{\rm out}\bigg[\frac{\partial x(n)}{\partial \rho };{\rm 0}_u \bigg]\!\cdot\! W^{\rm out}\bigg[\frac{\partial x(n)}{\partial a};{\rm 0}_u \bigg] \right)-\nonumber\\&\ell ^{\rm T}(n)W^{\rm out}\bigg[\frac{\partial ^2x(n)}{\partial a\partial \rho };{\rm 0}_u \bigg]+M^{(k)} {\frac{3}{(\rho -a)^3}}\label{eq22}\end{align}
(22) \begin{align}\varphi _{21} (n)\!=\!&\left( {W^{\rm out}\bigg[\frac{\partial x(n)}{\partial a};{\rm 0}_u \bigg]\!\cdot\! W^{\rm out}\bigg[\frac{\partial x(n)}{\partial \rho };{\rm 0}_u \bigg]} \right)-\nonumber\\&\ell^{\rm T} (n)W^{\rm out}\bigg[\frac{\partial ^2x(n)}{\partial\rho\partial a};{\rm 0}_u \bigg]+M^{(k)}{\frac{3}{(\rho -a)^3}} \label{eq23}\end{align}
(23) \begin{align}\varphi _{22} (n)=&\left( {W^{\rm out}\bigg[\frac{\partial x(n)}{\partial \rho };{\rm 0}_u \bigg]} \right)^2-\nonumber\\&\ell^{\rm T} (n)W^{\rm out}\bigg[\frac{\partial ^2x(n)}{\partial\rho ^2};{\rm 0}_u \bigg]-\nonumber\\&M^{(k)}\left( {\frac{1}{(\rho -a)^3}-\frac{2}{\rho ^3}}\right)\label{eq24}\end{align}
(24) 通过式(19) $\sim$ (24) 可知, 我们还必须求得偏导数 $\frac{\partial x\left( n \right)}{\partial a}$ , $\frac{\partial x\left( n\right)}{\partial \rho }$ , $\frac{\partial ^2x\left( n\right)}{\partial a^2}$ , $\frac{\partial ^2x\left( n\right)}{\partial a\partial \rho }$ , $\frac{\partial ^2x\left( n\right)}{\partial \rho \partial a}$ 和 $\frac{\partial ^2x\left( n\right)}{\partial \rho ^2}$ .其中, 为了表示方便, 令 ${\overline \zeta}\left( n \right)=s^{\rm in}W^{\rm in}u\left( n \right)+\left( {\rho W} \right)x\left( n \right)$ , 可以得到如下偏导数:
\begin{align}\frac{\partial x(n)}{\partial a}=&(1-a)\frac{\partial x(n-1) }{\partial a}-x\left( {n-1} \right)+\nonumber\\&{f}'({\overline \zeta } (n))\cdot \left( {(\rho W)\frac{\partial x(n-1) }{\partial a}} \right)\label{eq25}\end{align}
(25) \begin{align}\frac{\partial x\left( n \right)}{\partial \rho }=&\left( {1-a}\right)\frac{\partial x\left( {n-1} \right)}{\partial \rho }+{f}'({\overline \zeta }(n))\cdot \nonumber\\&\left( {(\rho W)\frac{\partial x(n-1) }{\partial \rho }+Wx(n-1) }\right)\label{eq26}\end{align}
(26) \begin{align}\frac{\partial ^2x(n)}{\partial a^2}=&(1-a)\frac{\partial^2x(n-1) }{\partial a^2}-2\cdot \frac{\partial x(n-1) }{\partial a}+\nonumber\\&{f}"({\overline \zeta } (n))\cdot \left( {(\rho W)\frac{\partial x(n-1) }{\partial a}}\right)^2+\nonumber\\&{f}'({\overline \zeta } (n))\cdot \left( {\left( {\rho W}\right)\frac{\partial ^2x\left( {n-1} \right)}{\partial a^2}}\right)\label{eq27}\end{align}
(27) \begin{align}\frac{\partial ^2x(n)}{\partial a\partial \rho}=&(1-a)\frac{\partial ^2x(n-1) }{\partial a\partial \rho}-\nonumber\\&\frac{\partial x(n-1) }{\partial \rho}+{f}"({\overline \zeta } (n))\cdot\nonumber\\& \left( {(\rho W)\frac{\partial x(n-1) }{\partial \rho }+Wx(n-1) }\right)\cdot\nonumber\\&\left( {\left( {\rho W} \right)\frac{\partial x\left( {n-1}\right)}{\partial a}} \right)+{f}'({\overline \zeta } (n))\cdot \nonumber\\&\left( {W\frac{\partial x(n-1) }{\partial a}+\left( {\rho W}\right)\frac{\partial ^2x\left( {n-1} \right)}{\partial a\partial\rho }} \right)\label{eq28}\end{align}
(28) \begin{align}\frac{\partial ^2x(n)}{\partial \rho \partial a}=&(1-a)\frac{\partial ^2x(n-1) }{\partial \rho\partial a}-\nonumber\\&\frac{\partial x(n-1) }{\partial \rho}+{f}"({\overline \zeta } (n))\cdot \nonumber\\&\left( {(\rho W)\frac{\partial x(n-1) }{\partial \rho }+Wx(n-1) }\right)\cdot\nonumber\\&\left( {(\rho W)\frac{\partial x(n-1) }{\partial a}} \right)+{f}'({\overline \zeta } (n))\cdot \nonumber\\&\left( {W\frac{\partial x(n-1) }{\partial a}+(\rho W)\frac{\partial^2x(n-1) }{\partial \rho \partial a}} \right)\label{eq29}\end{align}
(29) \begin{align}\frac{\partial ^2x(n)}{\partial \rho ^2}=&(1-a)\frac{\partial^2x(n-1) }{\partial \rho ^2}+{f}"({\overline \zeta }(n))\cdot\nonumber\\& \left( {(\rho W)\frac{\partial x(n-1) }{\partial \rho }+Wx(n-1) } \right)^2+\nonumber\\& {f}'({\overline \zeta }(n))\cdot\bigg( W\frac{\partial x(n-1) }{\partial \rho}+\nonumber\\&(\rho W)\frac{\partial ^2x(n-1) }{\partial \rho^2}+W\frac{\partial x(n-1) }{\partial \rho } \bigg)\label{eq30}\end{align}
(30) 式(25) $\sim$ (30) 提供了一个由 $\frac{\partial x\left( {n-1}\right)}{\partial q}$ , $\frac{\partial ^2x\left( {n-1}\right)}{\partial q^2}$ 到 $\frac{\partial x\left( n\right)}{\partial q}$ , $\frac{\partial ^2x\left( n\right)}{\partial q^2}$ 的递归计算方法.这些递归方程的初始值可以是任意值(通常我们都设为零矩阵, 即 $\frac{\partial x\left( {n-1} \right)}{\partial q}={\rm 0}$ 和 $\frac{\partial^2x\left( {n-1} \right)}{\partial q^2}={\rm 0})$ .在递归方程中, $f$ 函数都是双曲正切函数(tanh函数), 因此我们有
\begin{equation}\label{eq31} {f}'\left( {\overline \zeta } \right)=\frac{4}{2+{\rm e}^{2\overline \zeta }+{\rm e}^{-2\overline \zeta }}\end{equation}
(31) \begin{equation}\label{eq32} {f}''\left( {\overline \zeta } \right)=\frac{8\left({{\rm e}^{-2\overline \zeta }-{\rm e}^{2\overline \zeta }}\right)}{\left( {2+{\rm e}^{2\overline \zeta }+{\rm e}^{-2\overline \zeta }} \right)^2}\end{equation}
(32) 下面给出利用罚函数内点法寻找Leaky-ESN全局参数(主要是泄漏率 $a$ 和谱半径 $\rho$ )的最优值的具体实现步骤.
步骤1. 对罚函数内点法初始化.初始参数 $\zeta(0) $ 的值需要选取在可行域内; 根据式(12) 可以计算得到 $M^{0}$ ; 缩减因子 $c$ 的选取; 收敛精度; 迭代步数 $k=0$ .
步骤2. 利用牛顿法求解无约束优化问题.首先要判断 $\left\| {\varphi _1(n)-\varphi _2 (n)} \right\|\le \varepsilon$ 是否成立, 如果满足收敛条件, 那么此时经过优化得到的参数值就是牛顿优化算法下的最优值, 若不然, 我们就利用牛顿法来继续优化参数, 即利用式(32) 优化参数 $\zeta$ 直到满足收敛条件.
\begin{align}\zeta (n)=&\zeta (n-1) -\left[{{\begin{array}{*{20}c} {\varphi _{11} (n-1) } \hfill&{\varphi _{12} (n-1) } \hfill \\ {\varphi _{21} (n-1) } \hfill&{\varphi _{22} (n-1) } \hfill \\\end{array} }} \right]^{-1}\!\!×\nonumber\\& \left[{{\begin{array}{*{20}c} {\varphi _1 (n-1) } \hfill \\ {\varphi _2 (n-1) } \hfill \\\end{array} }} \right]^{\rm T}\label{eq33}\end{align}
(33) 步骤3. 罚函数内点法寻找全局最优值.当步骤2收敛条件满足时, 此时的最优值 $\zeta(n)$ 可能只是相对的最优值, 并不一定是最终的全局最优值.因此, 我们还需要判断 $\left\|{\zeta (n)-\zeta (n-1) } \right\|\le \varepsilon$ 是否成立, 如果收敛条件成立, 那么此时的参数值就是要寻找的全局最优值, 否则就是相对最优值.因此, 我们将此时的相对最优值作为下次迭代的初始值, 令 $M^{(k+1) }=cM^{(k)}$ , $k=k+1$ , 然后重复上述两步, 直到满足所有的收敛条件为止.
步骤4. 现在通过前面3个步骤就可以得到全局最优参数了, 同时根据式(3) 和式(4) 得到输出权值 $W^{\rm out}$ 和训练误差.
3. 仿真与分析
为了评价罚函数内点法的优化效果, 本文以时间序列预测为例, 与利用文献[10]中的随机梯度下降法(Gradient descent, GD)优化Leaky-ESN模型参数进行比较和分析.程序运行总的步长是20000步, 500次为一个样本点( $T=500$ ), 总共40个样本点.当Leaky-ESN模型学习完成之后, 利用50个样本数据来评估优化后的模型的预测精度.为了保证数据的科学性与合理性, 训练误差NRMSE和测试NRMSE, 预测误差等都是经过50次随机测试之后所得出的平均值.
3.1 第一个时间序列
时间序列模型为 $u\left( n \right)=\sin \left( n \right)+\sin \left({0.51n} \right)+\sin \left( {0.1002n} \right)+\sin \left( {0.054n}\right)$ , 其期望输出 $d(n)=u(n-5) $ .Leaky-ESN模型的初始值设置为 $\left[{a\left( 0 \right);\rho \left(0 \right)} \right]=\zeta \left( 0 \right)=[0.43;0.2]$ , $s^{\rm in}=0.3$ , $W^{\rm fb}=0$ . $W^{\rm in}$ 是由具有均匀伪随机分布的[ $-1, 1$ ]元素构成的矩阵. $W$ 是由具有均匀伪随机分布的非零元素产生具有单位谱半径的随机稀疏矩阵, 其稀疏连接度为0.5.连接输入权值 $W^{\rm in}$ 和内部连接权值 $W$ 在整个训练过程中不变, 储备池元素个数 $N=20$ .罚函数内点法的惩罚因子 $M^{(0) }=6$ , 缩减因子 $c=0.7$ ; 收敛精度 $\varepsilon=10^{-6}$ .梯度下降法的初始值选取相同.
图 1是采用随机梯度下降法和罚函数内点法分别优化Leaky-ESN时的训练误差变化曲线.由图 1可以看出, 罚函数内点法和随机梯度下降法均能收敛, 利用内点法优化Leaky-ESN比随机梯度下降法具有更小的NRMSE.另外, 图 1表明, 在总迭代步数中, 运用罚函数内点法优化后的模型收敛的有效迭代步数(优化过程中趋于收敛真正需要的迭代步数)更少, 说明罚函数内点法在优化时具有更快的收敛速度.图 2和图 3是随机测试中某次时间序列预测的结果.从图 2和图 3中可以看出, 罚函数内点法优化的Leaky-ESN模型预测结果精度更高, 因此与它们相对应的误差与零值线(虚线)的距离更为接近.表 1是采用随机梯度下降法和罚函数内点法分别优化Leaky-ESN时的测试NRMSE (50次测试的统计平均值).由表 1可以看出, 利用罚函数内点法和随机梯度下降法优化后的Leaky-ESN模型均具有良好的预测性能, 其中罚函数内点法的测试误差明显要比随机梯度下降法小.图 4是第5000步数据的NRMSE的分布区间.其中, 柱形图纵坐标表示出现在每个区间总共的次数, 横坐标上数字1代表的误差出现的区间是 $[0, 0.05]$ (因为在此统计区间误差出现的次数为零, 所以人为地将第一个统计区间扩大, 避免在图中空置), 数字2代表的区间是 $[0.051, 0.071]$ , 数字2和数字3的区间间隔为 $0.02$ , 以此类推.从图中可知, 罚函数内点法优化模型误差在区间 $[0.072, 0.092]$ 出现次数最多达到 $23$ 次, 随机梯度下降法在区间 $[0.198, 0.218]$ 出现最多为 $25$ 次, 因此罚函数内点法优化模型的标准均方根误差出现在更小区间的次数更多, 进一步说明了优化后的模型具有更高的预测精度和更好的预测稳定性.
表 1 Leaky-ESN的测试NRMSETable 1 Testing NRMSE of Leaky-ESN方法 第一时间序列 第二时间序列 随机梯度下降法 0.21240 0.17710 罚函数内点法 0.08698 0.01397 3.2 第二个时间序列
Mackey-Glass混沌时间序列模型:
\begin{align}\label{eq34} x(n+1) =& x(n)+\nonumber\\& \Delta T\left(\frac{{\alpha x\left( {t-\dfrac{\tau }{\Delta T}}\right)}}{1+x\left( t-\dfrac{\tau }{\Delta T} \right)^{10}}+\gamma x( n ) \right)\end{align}
(34) 式中, $\alpha =0.2$ , $\tau =17$ , $\gamma=-0.1$ , $\Delta T=1/10$ . $\tau $ 表示延迟因子, 当 $\tau>16.8$ 时序列(34) 具有混沌特性. Leaky-ESN模型, 罚函数内点法初始值设置与第一个时间序列一样.通过仿真可知, 与随机梯度下降法相比, 用罚函数内点法优化的Leaky-ESN模型仍具有良好的逼近性能.图 5描述了随机梯度下降法和罚函数内点法的训练误差NRMSE.从图 5中可知, 罚函数内点法优化的模型具有更小的NRMSE.图 6和图 7是随机测试中某次时间序列预测的结果.从图 6和图 7中可以看出罚函数内点法优化的模型预测结果精度更高.由表 1也可以看出, 与随机梯度下降法相比, 利用罚函数内点法优化后的Leaky-ESN模型均具有良好的预测性能.图 8是第5000步数据的NRMSE的分布区间.其中, 柱形图纵坐标表示NRMSE出现在每个区间总共的次数, 数字1表示误差的统计区间 $\left[{0, 0.01}\right]$ (同样地, 为了避免过多的在图中空置, 因此我们就将第一个统计区间人为地放大一点), 数字2表示误差的统计区间 $\left[{0.011, 0.015}\right]$ , 数字2和数字3的区间间隔为0.004, 以此类推.罚函数内点法优化的模型误差主要集中在区间 $\left[{0.011, 0.015}\right]$ , 总共出现22次, 并且分布比较紧凑, 说明稳定性相对比较好; 而随机梯度下降法优化的模型主要集中在后面几个区间, 而且分布相对比较均匀, 说明模型训练时会有一定的随机性和不稳定性.
由以上两个仿真例子可以说明, 与随机梯度下降法优化的Leaky-ESN相比, 罚函数内点法优化的Leaky-ESN具有更好的逼近精度.一旦参数的初始值满足可行域要求, 罚函数内点法将始终使参数工作在可行域, 因此不用像随机梯度下降法那样去实时判断参数是否满足回声特性条件, 避免算法被强行中断重新开始, 这样不仅节约了整个优化过程的时间, 而且对可行域中的任意一点初值都具有很好的收敛性和稳定性.
4. 结论
为了保障Leaky-ESN模型的回声特性, Leaky-ESN的泄漏率和储备池连接权谱半径等全局参数之间满足不等式约束条件.通常需要合理地选择Leaky-ESN的全局参数值, 有效地提升Leaky-ESN模型的性能.因此, 本文利用罚函数内点法在参数的可行域内进行搜索, 充分利用了目标函数的二阶导数信息, 能够获取Leaky-ESN模型的最优全局参数值, 提高了Leaky-ESN模型的预测精度.另外, 在参数修正过程中增加了目标函数的二阶导数信息, 使得计算量增大, 但与随机梯度下降法相比, 罚函数内点法将始终使参数工作在可行域内, 不用实时判断参数是否满足回声特性条件, 避免算法被强行中断重新开始, 实际上节约了整个优化过程的时间, 而且对可行域中的任意一点初值都具有很好的收敛性和稳定性.
-
表 1 Leaky-ESN的测试NRMSE
Table 1 Testing NRMSE of Leaky-ESN
方法 第一时间序列 第二时间序列 随机梯度下降法 0.21240 0.17710 罚函数内点法 0.08698 0.01397 
下载: 导出CSV
-
[1] Jaeger H. Tutorial on Training Recurrent Neural Networks, Covering BPTT, RTRL, EKF, and the "Echo State Network" Approach. Technical Report GMD Report 159, German National Research Center for Information Technology, German, 2002. [2] Jaeger H. The "Echo State" Approach to Analysing and Training Recurrent Neural Networks. Technical Report GMD report 148, German National Research Center for Information Technology, German, 2001. [3] Jaeger H, Haas H. Harnessing nonlinearity:predicting chaotic systems and saving energy in wireless communication. Science, 2004, 304(5667):78-80 doi: 10.1126/science.1091277 [4] Lun S X, Wang S, Guo T T, Du C J. An Ⅰ-Ⅴ model based on time warp invariant echo state network for photovoltaic array with shaded solar cells. Solar Energy, 2014, 105:529-541 doi: 10.1016/j.solener.2014.04.023 [5] Skowronski M D, Harris J G. Noise-robust automatic speech recognition using a predictive echo state network. IEEE Transactions on Audio, Speech, and Language Processing, 2007, 15(5):1724-1730 doi: 10.1109/TASL.2007.896669 [6] Han S I, Lee J M. Precise positioning of nonsmooth dynamic systems using fuzzy wavelet echo state networks and dynamic surface sliding mode control. IEEE Transactions on Industrial Electronics, 2013, 60(11):5124-5136 doi: 10.1109/TIE.2012.2218560 [7] Li G Q, Niu P F, Zhang W P, Zhang Y. Control of discrete chaotic systems based on echo state network modeling with an adaptive noise canceler. Knowledge-Based Systems, 2012, 35(15):35-40 http://www.sciencedirect.com/science/article/pii/S0950705112001153 [8] Song R Z, Xiao W D, Sun C Y. A new self-learning optimal control laws for a class of discrete-time nonlinear systems based on ESN architecture. Science China Information Sciences, 2014, 57(6):Article No. 068202 doi: 10.1007%2Fs11432-013-4954-y.pdf [9] Lun S X, Yao X S, Qi H Y, Hu H F. A novel model of leaky integrator echo state network for time-series prediction. Neurocomputing, 2015, 159(1):58-66 http://linkinghub.elsevier.com/retrieve/pii/S0925231215001782 [10] Bianchi F M, Scardapane S, Uncini A, Rizzi A, Sadeghian A. Prediction of telephone calls load using echo state network with exogenous variables. Neural Networks, 2015, 71(C):204-213 https://www.researchgate.net/publication/281677244_Prediction_of_telephone_calls_load_using_Echo_State_Network_with_exogenous_variables [11] Jaeger H, Lukoševičius M, Popovici D, Siewert U. Optimization and applications of echo state networks with leaky-integrator neurons. Neural Networks, 2007, 20(3):335-352 doi: 10.1016/j.neunet.2007.04.016 [12] Lukoševičius M. A practical guide to applying echo state networks. Neural Networks:Tricks of the Trade (Second Edition). Berlin Heidelberg:Springer-Verlag, 2012. 659-686 [13] Nocedal J, Wright S J. Numerical Optimization (Second Edition). New York:Springer, 2006. 30-31 -
下载:
计量
- 文章访问数: 1923
- HTML全文浏览量: 142
- PDF下载量: 732
- 被引次数: 0
