这一篇文章,主要是给非线性动力学,对微分方程模型参数拟合感兴趣的朋友写的。笼统的来说,这与混沌系统的预测有关;传统的机器学习的模式识别虽然也会谈论预测结果,但他们一般不会涉及连续的预测。这里我们考虑的是,连续的预测,而不仅是进行一步预测。比如生物神经元的发放序列,股票序列,天气序列等等。
需要大家已经具备的知识
- 了解什么是微分方程
- 了解最小二乘法
读完之后你会获得
文章目录
- 1 什么是动力学重构?有什么用?
- 如何直观理解
- 重构的数学表达
- 2 动力学重构论文调研
- 2017 - HOCC - High-order correlation computations
- HOCC-理论
- HOCC-例子:Lorenz系统
- 2022 Reconstruction of nonlinear flows from noisy time series
- 2023 Reconstructing equations of motion and networkstructures from noisy time series with hiddenvariables
- 附:相关链接
1 什么是动力学重构?有什么用?
动力学重构,在下面简称为重构。
如何直观理解
重构这个词,大家如果觉得陌生,那么参数拟合大家可能更好理解。只不过我们是对微分方程进行参数拟合。
比如,当我们手头有一个时间序列的数据(如下图这个FHN网络中,5个节点的膜电位序列),我们希望把这个时间序列背后的模型中的系数给拟合出来。
把系数拟合出来的好处有很多啦,比如:
(1)参数可能有物理意义上的可解释性,可以用来解释和理解系统,
(2)可以用来预测,这个系统之后会朝着怎样的方向演化,
(3)参数有了之后,其实模型就变成解析的了,非线性动力学有一套方法,比如寻找周期,分析临界态,分岔分析之类的方法,都能用上了。
重构的数学表达
考虑一个含有 N N N个状态变量的系统,它的状态变量可以表示为 x 1 , x 2 , . . . , x N x_1, x_2, ..., x_N x1,x2,...,xN,每一个状态变化规则,我们假设可以被一个连续性的函数所刻画,比如N=3的时候,我们可以写成
{ d x 1 d t = x ˙ 1 = f 1 ( x 1 , x 2 , x 3 , p ) d x 2 d t = x ˙ 2 = f 2 ( x 1 , x 2 , x 3 , p ) d x 2 d t = x ˙ 3 = f 3 ( x 1 , x 2 , x 4 , p ) \left\{\begin{array}{l} \frac{dx_1}{dt} = \dot{x}_1=f_1\left(x_1, x_2, x_3, p\right) \\ \frac{dx_2}{dt} = \dot{x}_2=f_2\left(x_1, x_2, x_3, p\right) \\ \frac{dx_2}{dt} = \dot{x}_3=f_3\left(x_1, x_2, x_4, p\right) \end{array}\right. ⎩ ⎨ ⎧dtdx1=x˙1=f1(x1,x2,x3,p)dtdx2=x˙2=f2(x1,x2,x3,p)dtdx2=x˙3=f3(x1,x2,x4,p)
可以看到,其中的每个变量都与其他变量有关, p p p 用来表示其中的待定参数。所谓的重构问题,是根据状态变量的 n n n个时间序列 x i ( 1 ) , x i ( 2 ) , . . . , x i ( n ) , i = 1 , 2 , . . . , N x_i(1),x_i(2),...,x_i(n), i=1,2,...,N xi(1),xi(2),...,xi(n),i=1,2,...,N,将其中参数还原出来。
如果有同学需要一个具体的微分方程组,也可以写出
v ˙ i = v i − v i 3 − w i − ∑ j = 1 N A i j ( v j − v i ) + Γ 1 , i w ˙ i = a i + b i + c i w i + Γ 2 , i (1) \begin{aligned} & \dot{v}_i=v_i-v_i^3-w_i-\sum_{j=1}^N A_{i j}\left(v_j-v_i\right)+\Gamma_{1, i} \\ & \dot{w}_i=a_i+b_i+c_i w_i+\Gamma_{2, i} \tag{1} \end{aligned} v˙i=vi−vi3−wi−j=1∑NAij(vj−vi)+Γ1,iw˙i=ai+bi+ciwi+Γ2,i(1)
这里的状态变量是 v i v_i vi和 w i w_i wi,里面的参数有 a i , b i , c i a_i,b_i,c_i ai,bi,ci,以及链接权重 A i j A_{ij} Aij;
事实上,这就是生成 【图1】所使用的微分方程模型。动力学重构就是基于【图1】的时间序列,将公式(1)中的参数都拟合出来。(文末会给出源码地址)
网络重构的概念
作为区分,我们这里补充介绍网络重构的概念。
类似于因果推断中,需要判断两个状态变量 x x x 和 y y 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 N 个状态的系统,其中有 n n n 个显状态, N − n N-n N−n 个隐状态。在每一个时刻都会面临一个高斯白噪声,数学表达如下:
x ˙ i ( t ) = f i ( x , y , α ) + Γ i ( t ) , i = 1 , 2 , ⋯ , n , y ˙ i ( t ) = f i ( x , y , α ) + Γ i ( t ) , i = n + 1 , n + 2 , ⋯ , N , (HOCC.1) \begin{aligned} & \dot{x}_i(t)=f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})+\Gamma_i(t), i=1,2, \cdots, n, \\ & \dot{y}_i(t)=f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})+\Gamma_i(t), i=n+1, n+2, \cdots, N, \tag{HOCC.1} \end{aligned} x˙i(t)=fi(x,y,α)+Γi(t),i=1,2,⋯,n,y˙i(t)=fi(x,y,α)+Γi(t),i=n+1,n+2,⋯,N,(HOCC.1)
对于这样一个系统,理论上我们可以计算 x x x 变量关于时间的二阶导数 x ¨ \ddot{x} x¨,因为其中只有 x , y \boldsymbol{x},\boldsymbol{y} x,y与时间有关,由链式法则可得
x ¨ i ( t ) = ∑ j = 1 n ∂ f i ( x , y , α ) ∂ x j x ˙ j + ∑ j = n + 1 N ∂ f i ( x , y , α ) ∂ y j y ˙ j + Γ i ′ ( t ) , (HOCC.2) \ddot{x}_i(t)=\sum_{j=1}^n \frac{\partial f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})}{\partial x_j} \dot{x}_j+\sum_{j=n+1}^N \frac{\partial f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})}{\partial y_j} \dot{y}_j+\Gamma_i^{\prime}(t), \tag{HOCC.2} x¨i(t)=j=1∑n∂xj∂fi(x,y,α)x˙j+j=n+1∑N∂yj∂fi(x,y,α)y˙j+Γi′(t),(HOCC.2)
为了简洁,做一个符号替换,定义2个符号。
ϕ
\phi
ϕ 是把(HOCC.2)右边搞简单,
Γ
\Gamma
Γ 是把噪声搞简单。
y
~
(
x
,
x
˙
,
α
)
\tilde{y}(x, \dot{x}, \alpha)
y~(x,x˙,α) 就是
y
(
x
,
x
˙
,
α
)
{y}(x, \dot{x}, \alpha)
y(x,x˙,α) 剔除噪声项后剩下的决定项。
ϕ
i
(
x
,
x
˙
,
α
)
=
[
∑
j
=
1
n
∂
f
i
(
x
,
y
,
α
)
∂
x
j
f
j
(
x
,
y
,
α
)
+
∑
j
=
n
+
1
N
∂
f
i
(
x
,
y
,
α
)
∂
y
j
f
j
(
x
,
y
,
α
)
]
∣
y
=
y
~
(
x
,
x
˙
,
α
)
\begin{aligned} \phi_i(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha})= & {\left[\sum_{j=1}^n \frac{\partial f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})}{\partial x_j} f_j(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})\right.} \left.+\sum_{j=n+1}^N \frac{\partial f_i(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})}{\partial y_j} f_j(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{\alpha})\right]_{\mid \boldsymbol{y}=\tilde{\boldsymbol{y}}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha})} \end{aligned}
ϕi(x,x˙,α)=[j=1∑n∂xj∂fi(x,y,α)fj(x,y,α)+j=n+1∑N∂yj∂fi(x,y,α)fj(x,y,α)]∣y=y~(x,x˙,α)
Γ i ′ ′ ( t ) = { Γ i ′ ( t ) + ∑ j = 1 N g i j ( x , x ˙ , α ) Γ j ( t ) i = 1 , 2 , ⋯ , N − n , Γ i ( t ) + ∑ j = 1 N g i j ( x , x ˙ , α ) Γ j ( t ) i = N − n + 1 , ⋯ , n , (HOCC.3) \Gamma_i^{\prime \prime}(t)=\left\{\begin{array}{c}\Gamma_i^{\prime}(t)+\sum_{j=1}^N g_{i j}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha}) \Gamma_j(t) \\ i=1,2, \cdots, N-n, \\ \Gamma_i(t)+\sum_{j=1}^N g_{i j}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha}) \Gamma_j(t) \\ i=N-n+1, \cdots, n,\end{array}\right. \tag{HOCC.3} Γi′′(t)=⎩ ⎨ ⎧Γi′(t)+∑j=1Ngij(x,x˙,α)Γj(t)i=1,2,⋯,N−n,Γi(t)+∑j=1Ngij(x,x˙,α)Γj(t)i=N−n+1,⋯,n,(HOCC.3)
于是,把上面两个公式带入(HOCC.2),得下面得形式
x ¨ i ( t ) = ϕ i ( x , x ˙ , α ) + Γ i ′ ′ ( t ) , i = 1 , 2 , ⋯ , N − n , x ˙ i ( t ) = f i ( x , y ( x , x ˙ , α ) , α ) + Γ i ′ ′ ( t ) , i = N − n + 1 , N − n + 2 , ⋯ , n , \begin{aligned} \ddot{x}_i(t)= & \phi_i(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha})+\Gamma_i^{\prime \prime}(t), i=1,2, \cdots, N-n, \\ \dot{x}_i(t)= & f_i(\boldsymbol{x}, \boldsymbol{y}(\boldsymbol{x}, \dot{\boldsymbol{x}}, \boldsymbol{\alpha}), \boldsymbol{\alpha})+\Gamma_i^{\prime \prime}(t), i=N-n+1, \\ & N-n+2, \cdots, n, \end{aligned} x¨i(t)=x˙i(t)=ϕi(x,x˙,α)+Γi′′(t),i=1,2,⋯,N−n,fi(x,y(x,x˙,α),α)+Γi′′(t),i=N−n+1,N−n+2,⋯,n,
下面作者这里做了一个关键操作,作者把 ϕ i \phi_i ϕi 直接完全用 x x x 来进行表示了,说是可以找到一个基向量 Z i Z_i Zi 来实现这样的效果。原文如下:
图2
ps 是否可以在 M i M_i Mi 取有限项的时候,将 ϕ i \phi_i ϕi 完美表示出来呢。为什么可以完美表示出来,如何表示,以及是否会存在例外?大家可以思考一下哦~
总之,通过这种基的选取,我们可以把方程全换成由
x
x
x 作为主导的形式,
x
¨
i
=
∑
μ
=
1
M
i
A
i
,
μ
Z
i
,
μ
(
x
,
x
˙
)
+
Γ
i
′
′
(
t
)
,
i
=
1
,
2
,
⋯
,
N
−
n
,
x
˙
i
=
∑
μ
=
1
M
i
A
i
,
μ
Z
i
,
μ
(
x
,
x
˙
)
+
Γ
i
′
′
(
t
)
,
i
=
N
−
n
+
1
,
⋯
,
n
.
\begin{aligned} & \ddot{x}_i=\sum_{\mu=1}^{M_i} A_{i, \mu} Z_{i, \mu}(\boldsymbol{x}, \dot{\boldsymbol{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}(\boldsymbol{x}, \dot{\boldsymbol{x}})+\Gamma_i^{\prime \prime}(t), i=N-n+1, \cdots, n . \end{aligned}
x¨i=μ=1∑MiAi,μZi,μ(x,x˙)+Γi′′(t),i=1,2,⋯,N−n,x˙i=μ=1∑MiAi,μZi,μ(x,x˙)+Γi′′(t),i=N−n+1,⋯,n.
现在来分析一下,上述方程中的
x
x
x 是已知的,选取的基的形式已知,所以,未知数只有参数矩阵
A
A
A。
如此简洁的线性形式,可以直接用线性最小二乘法来实现求解,且得到的是全局最优。
HOCC-例子:Lorenz系统
会不会觉得理论部分太抽象呢?我们下面以 Lorenz 系统为例,我们再做一个讲解。
首先,给出 Lorenz 系统的形式,
x
˙
=
α
1
y
+
α
2
x
+
Γ
1
(
t
)
,
y
˙
=
β
1
x
+
β
2
y
+
β
3
x
z
+
Γ
2
(
t
)
,
z
˙
=
γ
1
z
+
γ
2
x
y
+
Γ
3
(
t
)
.
\begin{aligned} & \dot{x}=\alpha_1 y+\alpha_2 x+\Gamma_1(t), \\ & \dot{y}=\beta_1 x+\beta_2 y+\beta_3 x z+\Gamma_2(t), \\ & \dot{z}=\gamma_1 z+\gamma_2 x y+\Gamma_3(t) . \end{aligned}
x˙=α1y+α2x+Γ1(t),y˙=β1x+β2y+β3xz+Γ2(t),z˙=γ1z+γ2xy+Γ3(t).
假设其中
y
y
y 是隐状态,先不管噪声;于是,当我们对
x
x
x 求二阶导,得
x
¨
=
α
1
y
˙
+
α
2
x
˙
:
可以将
y
˙
带入
=
α
1
(
β
1
x
+
β
2
y
+
β
3
x
z
)
+
α
2
x
˙
:
基于
α
1
y
=
x
˙
−
α
2
x
=
α
1
β
1
x
+
β
2
(
x
˙
−
α
2
x
)
+
α
1
β
3
x
z
+
α
2
⋅
x
˙
:
同类项合并
=
(
β
2
+
α
2
)
x
˙
+
(
α
1
β
1
−
α
2
β
2
)
x
+
α
1
β
3
x
z
(
HOCC-L
)
\begin{aligned} \ddot{x} & =\alpha_1 \dot{y}+\alpha_2 \dot{x} & \quad:可以将\dot{y}带入\\ & =\alpha_1\left(\beta_1 x+\beta_2 y+\beta_3 x z\right)+\alpha_2 \dot{x} & \quad:基于\alpha_1y=\dot{x}-\alpha_2 x \\ & = \alpha_1 \beta_1 x+\beta_2\left(\dot{x}-\alpha_2 x\right)+\alpha_1 \beta_3 x z+\alpha_2 \cdot \dot{x} & \quad: 同类项合并 \\ & =\left(\beta_2+\alpha_2\right) \dot{x}+\left(\alpha_1 \beta_1-\alpha_2 \beta_2\right) x +\alpha_1 \beta_3 x z &(\text{HOCC-L}) \end{aligned}
x¨=α1y˙+α2x˙=α1(β1x+β2y+β3xz)+α2x˙=α1β1x+β2(x˙−α2x)+α1β3xz+α2⋅x˙=(β2+α2)x˙+(α1β1−α2β2)x+α1β3xz:可以将y˙带入:基于α1y=x˙−α2x:同类项合并(HOCC-L)
类似的,我们可以把
z
z
z 也用
x
,
z
x,z
x,z 表示出来。如果如果我们把
x
1
x_1
x1 代表
x
x
x,
x
2
x_2
x2 代表
z
z
z,就可以得到
x
¨
1
=
(
β
2
+
α
2
)
x
˙
1
+
(
α
1
β
1
−
α
2
β
2
)
x
1
+
α
1
β
3
x
1
x
2
x
˙
2
=
γ
2
α
1
x
˙
1
x
1
−
γ
2
α
2
α
1
x
1
2
+
γ
1
x
2
\begin{aligned} \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{aligned}
x¨1=x˙2=(β2+α2)x˙1+(α1β1−α2β2)x1+α1β3x1x2α1γ2x˙1x1−α1γ2α2x12+γ1x2
此时,就实现了 显状态的演化只与显状态有关 的目的。后续使用最小二乘法求解即可。
大家都应该知道怎么求解吧(•_•),啰嗦一下:
就是把上面 x ¨ 1 , x ˙ 2 \ddot{x}_1, \dot{x}_2 x¨1,x˙2合并成一个矩阵 A A A,把他们右边所用到的项 x ˙ 1 , x 1 , x 1 x 2 , x 1 ˙ x 1 , x 1 2 , x 2 \dot{x}_1, x_1, x_1x_2, \dot{x_1}x_1, x_1^2, x_2 x˙1,x1,x1x2,x1˙x1,x12,x2设为另一个矩阵 B B B,然后参数设为矩阵 C C C;
于是方程可以被表达为 A = B C A=BC A=BC,
故参数矩阵: B = A C − 1 B=AC^{-1} B=AC−1,这在matlab中可以使用pinv
直接实现求解。
另外,很容易发现这个方法的一个“小”缺陷,
公式 (HOCC-L) 在进行代换时,用到了 α 1 y = x ˙ − α 2 x \alpha_1y=\dot{x}-\alpha_2 x α1y=x˙−α2x,隐状态竟然可以通过显状态直接表示,这么好运?
如果 x ˙ = α 1 x + α 2 y 2 + α 3 y z \dot{x}=\alpha_1 x+ \alpha_2 y^2+\alpha_3 yz x˙=α1x+α2y2+α3yz,并且 y , z y,z y,z 都是隐状态,由于无法完全将显状态分离变量,而且解的形式中仍然会有隐变量 z z z 的存在,那么构造完美的基函数将不那么容易 (if not impossible)
————
(对于这个问题,各位有没有解决办法呢?欢迎把回答写在评论区哦~)
留坑
2022 Reconstruction of nonlinear flows from noisy time series
2022北邮
2023 Reconstructing equations of motion and networkstructures from noisy time series with hiddenvariables
2023北邮
附:相关链接
代码 - Lorenz模型代码
代码 - FHN模型仿真代码
代码 - 2022,Wang代码
代码 - 2023,Yan代码
论文集