动力学重构/微分方程参数拟合 - 基于模型

这一篇文章，主要是给非线性动力学，对微分方程模型参数拟合感兴趣的朋友写的。笼统的来说，这与混沌系统的预测有关；传统的机器学习的模式识别虽然也会谈论预测结果，但他们一般不会涉及连续的预测。这里我们考虑的是，连续的预测，而不仅是进行一步预测。比如生物神经元的发放序列，股票序列，天气序列等等。

需要大家已经具备的知识

1. 了解什么是微分方程
2. 了解最小二乘法

读完之后你会获得

1 什么是动力学重构？有什么用？

动力学重构，在下面简称为重构。

如何直观理解

重构这个词，大家如果觉得陌生，那么参数拟合大家可能更好理解。只不过我们是对微分方程进行参数拟合。

比如，当我们手头有一个时间序列的数据（如下图这个FHN网络中，5个节点的膜电位序列），我们希望把这个时间序列背后的模型中的系数给拟合出来。

把系数拟合出来的好处有很多啦，比如：

(1)参数可能有物理意义上的可解释性，可以用来解释和理解系统，
(2)可以用来预测，这个系统之后会朝着怎样的方向演化，
(3)参数有了之后，其实模型就变成解析的了，非线性动力学有一套方法，比如寻找周期，分析临界态，分岔分析之类的方法，都能用上了。

图1

重构的数学表达

考虑一个含有$N$个状态变量的系统，它的状态变量可以表示为 $x_1,x_2,...,x_N$，每一个状态变化规则，我们假设可以被一个连续性的函数所刻画，比如N=3的时候，我们可以写成

$$\{\begin{array}{l}\frac{d{x}_1}{dt}={\dot{x}}_1=f_1\left(x_1,x_2,x_3,p\right)\\ \frac{d{x}_2}{dt}={\dot{x}}_2=f_2\left(x_1,x_2,x_3,p\right)\\ \frac{d{x}_2}{dt}={\dot{x}}_3=f_3\left(x_1,x_2,x_4,p\right)\end{array}$$

可以看到，其中的每个变量都与其他变量有关，$p$ 用来表示其中的待定参数。所谓的重构问题，是根据状态变量的$n$个时间序列 $x_i(1),x_i(2),...,x_i(n),i=1,2,...,N$，将其中参数还原出来。

如果有同学需要一个具体的微分方程组，也可以写出
$$\begin{array}{rl}& {\dot{v}}_i=v_i-v_i^3-w_i-\sum _{j=1}^N{A}_{ij}\left(v_j-v_i\right)+{\Gamma }_{1,i}\\ & {\dot{w}}_i=a_i+b_i+c_i{w}_i+{\Gamma }_{2,i}\end{array}\text{\hspace{1em}(1)}$$
这里的状态变量是 $v_i$和 $w_i$，里面的参数有$a_i,b_i,c_i$，以及链接权重$A_{ij}$；
事实上，这就是生成 【图1】所使用的微分方程模型。动力学重构就是基于【图1】的时间序列，将公式(1)中的参数都拟合出来。

（文末会给出源码地址）

网络重构的概念
作为区分，我们这里补充介绍网络重构的概念。
类似于因果推断中，需要判断两个状态变量 $x$ 和 $y$ 是否存在因果相关性。而网络重构的为了判断在一个网络中（每个节点可以视为一个状态变量），两个节点的是否存在连边。
目前上海交通大学周栋焯和相关团队所创作的：Causal connectivity measures for pulse-output network reconstruction: Analysis and applications，介绍了网络重构领域的一些基本方法：时延相关系数、时延互信息、格兰杰因果关系和传输熵。他们的工作开源在github，论文公开地址。

2 动力学重构论文调研

动力学重构是一个发展很快的领域，与网络重构（因果推断）不同，我们不仅要判定两个变量是有因果关系的，而且我们还需要用具体的数学模型将他们的公式给刻画出来。

如果没有公式，给出公式，这称为方程重构。如果有公式，根据数据计算将公式中的系数计算出来，这被称为参数重构。他们都属于动力学重构的范畴。

2017 - HOCC - High-order correlation computations

2017年，北京师范大学团队Yang Chen提出了基于High-order correlation（高阶导数相关性）来在含有隐节点和噪声的情况下进行重构。简述一下思想

HOCC-理论

定义一个 $N$ 个状态的系统，其中有 $n$ 个显状态（用 $x_i$ 表示），$N-n$ 个隐状态（用 $y_i$ 表示）。在每一个时刻都会面临一个高斯白噪声，数学表达如下：

$$\begin{array}{rl}& {\dot{x}}_i(t)=f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })+{\Gamma }_i(t),i=1,2,\cdots ,n,\\ & {\dot{y}}_i(t)=f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })+{\Gamma }_i(t),i=n+1,n+2,\cdots ,N,\end{array}\text{\hspace{1em}(HOCC.1)}$$

补充显状态和隐状态的含义
比如我们知道对于Lorenz系统而言，其方程为：
$$\begin{array}{rl}& \dot{x}(t)=a(y(t)-x(t)),\\ & \dot{y}(t)=rx(t)-y(t)-x(t)z(t),\\ & \dot{z}(t)=xy(t)-bz\end{array}$$
我们其实可以认为是有3个状态变量 $x,y,z$ 在相互影响，而这种影响关系——比如：
$\dot{x}$ 表达式中除了自身之外，还含有 $y$。
$\dot{y}$ 除了自身之外含 $x,z$，
$\dot{z}$ 除了自身之外含 $x,y$，将这样的关系，通过网络图来表示如下：

现在来考虑隐变量的含义，假设我们有一台机器，只能够测量 $x,z$ 状态的时间序列，而无法测量 $y$ 的时间序列。这时，就称这个问题为：在 $y$ 为隐状态的情况下进行动力学重构的问题。

对于这样一个系统，理论上我们可以计算 $x$ 变量关于时间的二阶导数 $\ddot{x}$，因为其中只有$\mathbf{x},\mathbf{y}$与时间有关，由链式法则可得

$${\ddot{x}}_i(t)=\sum _{j=1}^n\frac{\partial f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })}{\partial x_j}{\dot{x}}_j+\sum _{j=n+1}^N\frac{\partial f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })}{\partial y_j}{\dot{y}}_j+{\Gamma }_i^{\prime }(t),\text{\hspace{1em}(HOCC.2)}$$

为了简洁，做一个符号替换，定义2个符号。$\varphi$ 是把(HOCC.2)右边搞简单，$\Gamma$ 是把噪声搞简单。$\tilde{y}(x,\dot{x},\alpha )$ 就是 $y(x,\dot{x},\alpha )$ 剔除噪声项后剩下的决定项。
$$\begin{array}{rl}{\varphi }_i(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha })=& [\sum _{j=1}^n\frac{\partial f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })}{\partial x_j}{f}_j(\mathbf{x},\mathbf{y},\mathbf{\alpha }){+\sum _{j=n+1}^N\frac{\partial f_i(\mathbf{x},\mathbf{y},\mathbf{\alpha })}{\partial y_j}{f}_j(\mathbf{x},\mathbf{y},\mathbf{\alpha })]}_{\mid \mathbf{y}=\tilde{\mathbf{y}}(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha })}\end{array}$$

$${\Gamma }_i^{\prime \prime }(t)=\{\begin{array}{c}{\Gamma }_i^{\prime }(t)+{\sum }_{j=1}^N{g}_{ij}(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha }){\Gamma }_j(t)\\ i=1,2,\cdots ,N-n,\\ {\Gamma }_i(t)+{\sum }_{j=1}^N{g}_{ij}(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha }){\Gamma }_j(t)\\ i=N-n+1,\cdots ,n,\end{array}\text{\hspace{1em}(HOCC.3)}$$

于是，把上面两个公式带入(HOCC.2)，得下面得形式

$$\begin{array}{rl}{\ddot{x}}_i(t)=& {\varphi }_i(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha })+{\Gamma }_i^{\prime \prime }(t),i=1,2,\cdots ,N-n,\\ {\dot{x}}_i(t)=& f_i(\mathbf{x},\mathbf{y}(\mathbf{x},\dot{\mathbf{x}},\mathbf{\alpha }),\mathbf{\alpha })+{\Gamma }_i^{\prime \prime }(t),i=N-n+1,\\ & N-n+2,\cdots ,n,\end{array}$$

下面作者这里做了一个关键操作，作者把 ${\varphi }_i$ 直接完全用 $x$ 来进行表示了，说是可以找到一个基向量 $Z_i$ 来实现这样的效果。原文如下：

图2

ps 是否可以在 $M_i$ 取有限项的时候，将 ${\varphi }_i$ 完美表示出来呢。为什么可以完美表示出来，如何表示，以及是否会存在例外？大家可以思考一下哦~

总之，通过这种基的选取，我们可以把方程全换成由 $x$ 作为主导的形式，
$$\begin{array}{rl}& {\ddot{x}}_i=\sum _{\mu =1}^{M_i}{A}_{i,\mu }{Z}_{i,\mu }(\mathbf{x},\dot{\mathbf{x}})+{\Gamma }_i^{\prime \prime }(t),i=1,2,\cdots ,N-n,\\ & {\dot{x}}_i=\sum _{\mu =1}^{M_i}{A}_{i,\mu }{Z}_{i,\mu }(\mathbf{x},\dot{\mathbf{x}})+{\Gamma }_i^{\prime \prime }(t),i=N-n+1,\cdots ,n.\end{array}$$

现在来分析一下，上述方程中的 $x$ 是已知的，选取的基的形式已知，所以，未知数只有参数矩阵 $A$。
如此简洁的线性形式，可以直接用线性最小二乘法来实现求解，且得到的是全局最优。

HOCC-例子：Lorenz系统

会不会觉得理论部分太抽象呢？我们下面以 Lorenz 系统为例，我们再做一个讲解。

首先，给出 Lorenz 系统的形式，
$$\begin{array}{rl}& \dot{x}={\alpha }_{1}y+{\alpha }_{2}x+{\Gamma }_1(t),\\ & \dot{y}={\beta }_{1}x+{\beta }_{2}y+{\beta }_{3}xz+{\Gamma }_2(t),\\ & \dot{z}={\gamma }_{1}z+{\gamma }_{2}xy+{\Gamma }_3(t).\end{array}$$

假设其中 $y$ 是隐状态，先不管噪声；于是，当我们对 $x$ 求二阶导，得
$$\begin{array}{rlr}\ddot{x}& ={\alpha }_1\dot{y}+{\alpha }_2\dot{x}& :可以将\dot{y}带入\\ & ={\alpha }_1\left({\beta }_{1}x+{\beta }_{2}y+{\beta }_{3}xz\right)+{\alpha }_2\dot{x}& :基于{\alpha }_{1}y=\dot{x}-{\alpha }_{2}x\\ & ={\alpha }_1{\beta }_{1}x+{\beta }_2\left(\dot{x}-{\alpha }_{2}x\right)+{\alpha }_1{\beta }_{3}xz+{\alpha }_2\cdot \dot{x}& :同类项合并\\ & =\left({\beta }_2+{\alpha }_2\right)\dot{x}+\left({\alpha }_1{\beta }_1-{\alpha }_2{\beta }_2\right)x+{\alpha }_1{\beta }_{3}xz& (\text{HOCC-L})\end{array}$$

类似的，我们可以把 $z$ 也用 $x,z$ 表示出来。如果如果我们把 $x_1$ 代表 $x$，$x_2$ 代表 $z$，就可以得到
$$\begin{array}{rl}{\ddot{x}}_1=& \left({\beta }_2+{\alpha }_2\right){\dot{x}}_1+\left({\alpha }_1{\beta }_1-{\alpha }_2{\beta }_2\right)x_1\\ & +{\alpha }_1{\beta }_3{x}_1{x}_2\\ {\dot{x}}_2=& \frac{{\gamma }_2}{{\alpha }_1}{\dot{x}}_1{x}_1-\frac{{\gamma }_2{\alpha }_2}{{\alpha }_1}{x}_1^2+{\gamma }_1{x}_2\end{array}$$

此时，就实现了 显状态的演化只与显状态有关 的目的。后续使用最小二乘法求解即可。

大家都应该知道怎么求解吧(•_•)，啰嗦一下：
就是把上面${\ddot{x}}_1,{\dot{x}}_2$合并成一个矩阵$A$，把他们右边所用到的项${\dot{x}}_1,x_1,x_1{x}_2,\dot{x_1}{x}_1,x_1^2,x_2$设为另一个矩阵$B$，然后参数设为矩阵$C$；
于是方程可以被表达为 $A=BC$，
故参数矩阵：$B=A{C}^{-1}$，这在matlab中可以使用pinv直接实现求解。

另外，很容易发现这个方法的一个“小”缺陷，
公式 (HOCC-L) 在进行代换时，用到了 ${\alpha }_{1}y=\dot{x}-{\alpha }_{2}x$，隐状态竟然可以通过显状态直接表示，这么好运？
如果 $\dot{x}={\alpha }_{1}x+{\alpha }_2{y}^2+{\alpha }_{3}yz$，并且 $y,z$ 都是隐状态，由于无法完全将显状态分离变量，而且解的形式中仍然会有隐变量 $z$ 的存在，那么构造完美的基函数将不那么容易 (if not impossible)
————
（对于这个问题，各位有没有解决办法呢？欢迎把回答写在评论区哦~）

---

留坑/草稿

2022 Reconstruction of nonlinear flows from noisy time series

2022北邮, Juanjuan Wang, Yueheng Lan

2023 Reconstructing equations of motion and network structures from noisy time series with hidden variables - [Yan,2023]

2023 北邮，Zishuo Yan, Yueheng Lan

这篇论文是在含有隐变量和噪声的情况下去进行重构，我们先看看，它在含有隐变量的 Kuramoto系统 中的表现。作者从一个32大小的网络中，选取了4个节点作为隐节点。并用论文方法来对4个隐节点的自然频率进行推断。一些基本设定

• 网络大小：32，

• 隐节点：4个，并且两两不为邻居。原文：In the following, the hidden nodes are randomly selected from the network and no pair are neighbours，但我看图(A)的上面两个绿色节点就好像是邻居？（也许图有问题，但就算隐节点是邻居，也是可以推断出来的。图14就是这样的实验）

• 连接矩阵中的元素 $a_{ij}$：1或0；（这里连接的结构应当是已知的？

• 耦合强度: $K=1$

• 自然频率$\omega$：选取范围是从0到20的均匀分布
  
  

• 从图11来看，连接强度的误差为0.03，还不错。但自然频率的话，好像是会有细微波动。

• 图12显示，随着节点增加精度会下降；随着隐节点比例的增多，精度会下降。

• 图14是用与隐状态连接的节点，来推断隐节点的结构的拟合结果，基本上都推断出来了。

• 图15是自然频率的范围的不同，对于推断精度的影响。发现如果网络节点的自然频率有非常大的差别，会有利于推断精度的提高。

• 图16(F) 展示了随机自然频率范围 $\mu$ 的变大，序参量(刻画同步性的量)开始降低，精度开始提高，并且GP关联维数也开始提高。

附：相关链接

Kuramoto系统简介

GP关联维数与计算源码

Flowus博客

代码 - Lorenz模型代码

代码 - FHN模型仿真代码

代码 - 2022，Wang代码

代码 - 2023，Yan代码