16、非线性动态系统学习：算法、应用与扩展-优快云博客

本文链接：https://blog.youkuaiyun.com/spark7igniter/article/details/154627402

非线性动态系统学习：算法、应用与扩展

1. 算法效果评估

将算法恢复的结构与生成系统的已知结构进行比较，结果表明该算法能很好地恢复非线性动力学的形式。通过用与训练期间方差相似的高斯噪声激励已学习的模型，还可以生成“幻想”数据，这些生成的观测流在性质上与真实系统的时间序列非常相似。

为了量化拟合质量，我们比较了非线性模型下训练序列和新测试序列的对数似然与基本线性动态系统模型或静态模型（如因子分析）下的似然。结果显示，对于所有示例系统，非线性动态系统在训练数据和测试数据上的似然性都显著更优。不过，对于系统E，由于其存在强烈的滞后现象（锁模），线性动态系统比因子分析要好得多，因为前一个时间步的输出是当前输出的优秀预测器。

2. 天气数据应用实例

以澳大利亚墨尔本1981 - 1990年的每日最高和最低气温记录为例，使用一个具有两个内部状态变量、三个输出且无输入的模型进行训练。训练阶段的三个输出分别是每日最低和最高气温以及表示一年中时间（月份）的实值输出（范围为[0, 12]）。

模型在1500天的气温记录（略超过四个季节）上进行训练，然后在剩余的2150天进行测试。测试时，仅向模型提供每日最低和最高气温，尝试预测一年中的时间（月份）。具体步骤如下：
1. 使用EKS算法根据两个可用的观测流进行状态推断。
2. 完成状态推断后，使用模型学习到的输出函数预测一年中的时间。

由于春季和秋季气温与季节的静态关系不明确，这个预测问题本质上需要利用过去和/或未来的信息。训练后的预测结果表明，算法发现了隐藏状态与观测之间的关系，能够进行合理的预测。此外，模型还能根据推断的状态预测最低和最高气温。

通过在状态空间的多个点评估非线性输出函数，可以发现模型隐式地参数化了一年中时间与温度之间的关系。每个评估产生一个包含月份、最低气温和最高气温的三元组，将这些三元组相互绘制，可以展示模型发现的墨尔本季节性温度变化。

以下是不同模型对训练和测试数据分配的对数似然差异对比表格：
| 模型 | 训练数据 | 测试数据 |
| — | — | — |
| 因子分析（FA） | 较低 | 较低 |
| 线性动态系统（LDS） | 中等 | 中等 |
| 非线性动态系统模型（NLDS） | 较高 | 较高 |

3. 算法扩展

3.1 学习RBF的均值和宽度

可以放宽高斯径向基函数（RBF）具有固定中心和宽度的假设，但这会导致拟合算法更加复杂和缓慢。为了推导RBF中心$c_i$和宽度矩阵$S_i$的学习规则，需要考虑它们如何通过RBF核（6.14）影响成本函数（6.17）。对成本函数$c$的期望求导数，并交换期望和导数的顺序：
$\frac{\partial c}{\partial c_i} = \frac{\partial c}{\partial r_i} \frac{\partial r_i}{\partial c_i} = 2 \langle (yF - z)^> Q^{-1} h_i r_i(x) S^{-1}_i (x - c_i) \rangle$

由于$c_i$在该方程的多个地方非线性出现，因此无法以封闭形式求解$c_i$。不过，可以使用上述梯度移动中心$c_i$以降低成本，这相当于对$c_i$进行部分M步。方程（6.19）除了优化$y$和$Q$所需的一阶和二阶期望外，还需要计算三个三阶期望：$\langle r_i(x) r_k(x) x_l \rangle$，$\langle r_i(x) x_k x_l \rangle$和$\langle r_i(x) z_k x_l \rangle$。

类似地，对$S^{-1}_i$求成本的导数得到：
$\frac{\partial c}{\partial S^{-1}_i} = \frac{\partial c}{\partial r_i} \frac{\partial r_i}{\partial S^{-1}_i} = \langle [(yF - z)^> Q^{-1} h_i] r_i(x) [S_i - (x - c_i) (x - c_i)^>] \rangle$

此时还需要三个四阶期望：$\langle r_i(x) r_k(x) x_l x_m \rangle$，$\langle r_i(x) k x_l x_m \rangle$和$\langle r_i(x) z_k x_l x_m \rangle$。

这些额外的期望增加了算法的存储和计算时间，而通过小梯度步骤移动中心和宽度所带来的优势可能无法弥补这一成本。一种启发式方法是使用无监督技术（如高斯混合的EM算法）来放置中心和宽度，该方法仅考虑输入密度而不考虑输出非线性。或者，可以使用近似方法，例如$\langle r_i(x) \rangle \approx r_i(\langle x \rangle)$。

3.2 在线学习

之前介绍的算法是批量算法，即假设使用整个观测序列来估计模型参数。不过，可以相对容易地推导出该算法的在线版本，通过递归最小二乘法（RLS）算法在接收观测值时更新参数。RLS算法实际上是离散卡尔曼滤波器的一个特例。

对于标量$z$，在线算法如下：
$y_k = y_{k - 1} + (\langle zF \rangle_k - y_{k - 1} \langle FF^> \rangle_k) P_k$
$P_k = P_{k - 1} - \frac{P_{k - 1} \langle FF^> \rangle_k P_{k - 1}}{1 + \langle F^> P_{k - 1} F \rangle_k}$
$Q_k = Q_{k - 1} + \frac{1}{k} [\langle z^2 \rangle_k - y_k \langle Fz \rangle_k - Q_{k - 1}]$

初始化$y_0 = 0$，$Q_0 = I$，$P_0$非常大，经过几次迭代后，$y_k$的估计值将迅速收敛到最小二乘法得到的精确值。$Q$的估计值将收敛到正确值加上一个偏差，这是因为早期$Q$的估计是基于$y_k$而不是$\lim_{k \to \infty} y_k$的残差。$P_k$是$(\sum_{j = 1}^{k} \langle FF^> \rangle_j)^{-1}$的递归估计，通过矩阵求逆引理获得。

在线算法中的期望$\langle \cdot \rangle_k$是通过使用前一个参数$y_{k - 1}$运行一步扩展卡尔曼滤波器来计算的，而批量EM算法中的期望是通过使用当前参数估计对整个序列运行扩展卡尔曼平滑器来计算的。一般来说，除非时间序列是非平稳的，否则批量算法收敛后得到的参数估计对数据的建模效果会比在线算法更好。

有趣的是，这里描述的RLS在线算法的更新与双扩展卡尔曼滤波器方法用于系统识别的参数更新非常相似，这种相似性可以用于推导非平稳非线性动态系统的在线参数估计算法。

下面是在线学习算法的流程mermaid图：

graph LR
    A[接收观测值] --> B[计算期望]
    B --> C[更新y_k]
    C --> D[更新P_k]
    D --> E[更新Q_k]
    E --> A

3.3 处理非平稳性

为了处理非平稳时间序列，假设参数可以根据具有协方差$S_y$的高斯随机游走漂移：
$y_k = y_{k - 1} + E_k$，其中$E_k \sim N(0, S_y)$

同时，$z$变量与参数$y$和非线性核$F$的关系为：
$z_k = y_k F_k + w_k$，其中$w_k \sim N(0, Q)$

可以将其视为“状态变量”$y_k$具有时变“输出矩阵”$F_k$的观测模型。由于动力学和观测模型在$y$上都是线性的，且噪声是高斯的，因此可以应用以下卡尔曼滤波器递归计算漂移参数$y$的分布：
$\hat{y} k = \hat{y} {k - 1} + \frac{(\langle zF \rangle_k - \hat{y} {k - 1} \langle FF^> \rangle_k) P {k|k - 1}}{Q_{k - 1} + \langle F^> P_{k|k - 1} F \rangle_k}$
$P_{k|k - 1} = P_{k - 1} + S_y$
$P_k = P_{k|k - 1} - \frac{P_{k|k - 1} \langle FF^> \rangle_k P_{k|k - 1}}{Q_{k - 1} + \langle F^> P_{k|k - 1} F \rangle_k}$
$Q_k = Q_{k - 1} + \lambda (\langle z^2 \rangle_k - \hat{y} k \langle Fz \rangle_k - Q {k - 1})$

需要注意的是：
1. 这些方程描述的是一个普通的卡尔曼滤波器，但“输出”$z$和“输出矩阵”$F_k$都具有高斯分布的联合不确定性。
2. 通过在其重新估计方程中引入遗忘因子$\lambda$，假设输出噪声协方差可以漂移。

该在线算法可以看作是对参数和隐藏变量的贝叶斯后验的近似，通过两个独立的高斯分布（一个关于$(x; z)$，一个关于$y$）递归近似真实的后验分布。$y_k$的近似后验分布具有均值$\hat{y}_k$和协方差$P_k$。

3.4 使用贝叶斯方法进行模型选择和复杂度控制

与其他最大似然程序一样，本文描述的EM算法有过拟合数据集的风险，即可能在数据噪声中找到虚假模式，从而导致泛化能力差。在实现中，使用了一些岭回归（对$h_i$参数进行权重衰减正则化），在实践中似乎效果不错，但需要一些启发式方法来设置正则化参数。

最大似然方法面临的另一个相关问题是缺乏内置的模型选择程序。例如，选择状态空间$x$的维度和基函数的数量$I$时，更高的$x$维度和更多的基函数原则上总是会导致似然的更高最大值，这意味着更复杂的模型总是会比简单模型更受青睐，但这会导致过拟合。

贝叶斯方法提供了一个通用框架，可以一致地处理过拟合和模型选择问题。其关键思想是尽可能避免最大化，而是根据后验概率对可能的模型、结构和参数进行加权，并根据这个加权后验进行预测。

对于非线性动态系统，可以将参数$y$视为未知量。那么模型对时间$k + 1$的输出预测为：
$p(y_{k + 1} | u_{1:k + 1}, y_{1:k}) = \int dy p(y_{k + 1} | u_{k + 1}, y_{1:k}, u_{1:k}, y) p(y | y_{1:k}, u_{1:k})$
$= \int dy p(y | y_{1:k}, u_{1:k}) \int dx_{k + 1} p(y_{k + 1} | u_{k + 1}, x_{k + 1}, y) p(x_{k + 1} | u_{1:k + 1}, y_{1:k}, y)$

参数的后验分布可以通过贝叶斯规则递归获得：
$p(y | y_{1:k}, u_{1:k}) = \frac{p(y_k | u_{1:k}, y_{1:k - 1}, y) p(y | U_{1:k - 1}, y_{1:k - 1})}{p(y_k | u_{1:k}, y_{1:k - 1})}$

双扩展卡尔曼滤波器、联合扩展卡尔曼滤波器和非平稳在线算法都是这些贝叶斯递归的粗略近似。

可以使用贝叶斯规则根据模型结构的概率选择不同的模型结构：
$P(s_m | y_{1:k}, u_{1:k}) = \frac{p(y_{1:k} | u_{1:k}, s_m) P(s_m)}{\sum_n p(y_{1:k} | u_{1:k}, s_n) P(s_n)}$

获得所需积分的可处理近似有以下几种方法：
1. 马尔可夫链蒙特卡罗（MCMC）技术 ：用于对参数和隐藏变量进行采样。如果样本集中在参数和状态概率较高的区域，采样可以是计算高维积分的有效方法。例如，Gibbs采样已用于线性动态系统，而粒子滤波是一种有前途的非线性系统方法，其中样本（“粒子”）可用于表示每个时间步参数和隐藏状态的联合分布。
2. 自动相关性确定（ARD） ：对每个参数使用具有可调方差的零均值高斯先验。这些方差是控制模型参数先验分布的超参数，优化这些超参数并对参数进行积分会导致“无关”参数从模型中消除。例如，Tipping成功地将具有每个数据点中心的RBF网络的ARD用于非线性回归，并将其命名为“相关向量机”。
3. 变分方法 ：用于对模型结构后验概率进行下界估计。与EM可以被视为使用隐藏变量的分布对似然进行下界估计一样，变分贝叶斯方法使用隐藏变量和参数的分布对证据进行下界估计。变分贝叶斯方法已用于推断线性动态系统的结构，但将其推广到本文描述的非线性系统并不简单。

4. 模型的可识别性和表达能力

从上述实验可以看出，本文提出的算法能够为各种非线性时间序列学习良好的密度模型。该算法能够很好建模的非线性系统类别定义了其表达能力，而模型在原则上恢复特定非线性系统实际参数的能力则涉及模型的可识别性问题。这两个问题密切相关，因为它们都描述了实际非线性系统与模型参数设置之间的映射。

模型存在三种技术上导致不可识别的平凡退化情况，但无需过多关注：
1. 状态空间维度置换 ：可以对状态空间的维度进行置换，并相应地置换输出映射和动力学的定义域，从而获得完全等效的模型。
2. 状态变量线性变换 ：状态变量可以进行缩放或通过任何可逆线性映射进行变换，这种变换可以被输出和动力学函数吸收，产生具有相同输入 - 输出行为的模型。为了不失一般性，通常将状态演化噪声的协方差设置为单位矩阵，这既设定了状态空间的尺度，又在不降低模型表达能力的情况下禁止了某些状态变换。
3. 噪声吸收 ：将观测噪声与状态噪声视为不相关，且两者均值为零，因为这些可以在不损失一般性的情况下被吸收到$f$和$g$函数中。

然而，还存在其他更难以克服的不可识别形式。例如，如果$f$和$g$都是非线性的，那么（至少在无噪声情况下），对于状态的任何任意可逆变换，模型都可能不可识别。

以下是可识别性相关情况的表格总结：
| 情况 | 说明 | 是否影响表达能力 |
| — | — | — |
| 状态空间维度置换 | 可获得等效模型 | 否 |
| 状态变量线性变换 | 产生相同输入 - 输出行为模型 | 否 |
| 噪声吸收 | 可被函数吸收 | 否 |
| 其他非线性情况 | 可能导致不可识别 | 可能 |

非线性动态系统学习：算法、应用与扩展

5. 总结与展望

在前面的内容中，我们详细探讨了非线性动态系统学习的相关算法、应用以及各种扩展方法。下面对这些内容进行总结，并对未来的研究方向进行展望。

5.1 内容总结

算法效果评估 ：通过与已知结构对比以及对数似然比较，验证了算法在恢复非线性动力学形式和拟合数据方面的有效性，非线性动态系统在训练和测试数据上表现更优。
天气数据应用 ：以墨尔本气温数据为例，展示了算法在实际问题中的应用，能够学习到时间与温度的关系并进行合理预测。
算法扩展
- 学习RBF的均值和宽度 ：放宽RBF固定中心和宽度假设会增加算法复杂度，但可通过启发式方法或近似方法处理。
- 在线学习 ：基于RLS算法实现了在线版本，能在接收观测值时更新参数，但批量算法在平稳时间序列上建模效果更好。
- 处理非平稳性 ：采用高斯随机游走假设和卡尔曼滤波器处理非平稳时间序列，将在线算法视为贝叶斯后验的近似。
- 贝叶斯方法 ：用于模型选择和复杂度控制，避免过拟合，通过多种近似方法处理积分问题。
可识别性和表达能力 ：算法具有一定的表达能力，但存在一些导致模型不可识别的情况，部分情况不影响表达能力，而其他非线性情况可能带来挑战。

5.2 未来研究方向

虽然本文提出的算法在非线性时间序列建模方面取得了一定的成果，但仍有许多可以改进和拓展的地方。以下是一些值得深入研究的方向：
- 更高效的算法优化 ：在学习RBF的均值和宽度以及处理高阶期望时，算法的存储和计算时间增加。未来可以探索更高效的算法或近似方法，以减少计算成本，提高算法的实用性。
- 在线学习的改进 ：在线算法在处理非平稳时间序列时是对贝叶斯后验的近似，可进一步研究如何提高这种近似的精度，使其在非平稳环境下表现更优。
- 贝叶斯方法的完善 ：贝叶斯方法在模型选择和复杂度控制方面具有优势，但在处理高维积分和复杂模型时存在挑战。可以深入研究更有效的积分近似方法，如改进MCMC技术、优化ARD和变分方法，以更好地应用于非线性动态系统。
- 可识别性问题的解决 ：对于更复杂的不可识别情况，需要进一步研究有效的解决方案，以提高模型恢复实际参数的能力，增强模型的可靠性。

下面是未来研究方向的流程图：

graph LR
    A[算法优化] --> B[高效学习RBF]
    A --> C[降低计算成本]
    D[在线学习改进] --> E[提高近似精度]
    D --> F[适应非平稳环境]
    G[贝叶斯方法完善] --> H[改进积分近似]
    G --> I[优化现有方法]
    J[可识别性解决] --> K[克服复杂不可识别情况]
    J --> L[增强模型可靠性]