这一篇文章叫做 数据驱动的坐标发现与方程发现算法。
想回答的问题很简单，“如何根据数据写方程”。

想想牛顿的处境，如何根据各种不同物体下落的数据，写出万有引力的数学公式的。这篇文章就是来做这件事的。当然，这篇论文并没有从牛顿视角，完全去思考牛顿所想。而是利用现有的深度学习技术 和 动力学重构的方法 。提出了一种框架，基于现有的计算机技术，去发现物质运动背后的物理规律。

这里直接给出这篇方法的核心思路图，我们后面会逐个讲解。

正文

Data-driven discovery of coordinates and governing equations
作者包括华盛顿大学的 Kathleen Champion 和 Steven L. Brunton。Brunton的SINDy方法，我们在 基于RNN进行动力学重构 的文章中也提到过。

这篇文章的主要创新点在于

1. 在数据降维这块，提出可以用神经网络的算法来进行数据编码，从而来实现降维。
   这在文中被称为自编码器，学习一个用于表示数据的坐标系，希望经过坐标系的变换，能够提取出更加明显和有效的动力学特征。
2. 结合了非线性动力学识别(SINDy)的方法，来进行参数的拟合。这在下文中也会介绍。

核心方法

首先，我们考虑一个 $n$ 维状态的动力学系统如公式(1)所示，其中 $\mathbf{x}(t)\in {\mathbb{R}}^n$ 。

$$\frac{d}{dt}\mathbf{x}(t)=\mathbf{f}(\mathbf{x}(t))\text{\hspace{1em}(1)}$$

我们的目标是

1. 寻求一组具有相关动态模型的约化坐标。即 $\mathbf{z}(t)=\phi (\mathbf{x}(t))\in {\mathbb{R}}^d(d\ll n)$，这个新坐标下的一个动力学维数远远小于原始系统。
2. 能够给出在这个简约坐标下对应的动力学方程
   $$\frac{d}{dt}\mathbf{z}(t)=\mathbf{g}(\mathbf{z}(t))\text{\hspace{1em}(2)}$$
3. 这个方法同时还提出对应的编码函数 $\mathbf{z}=\phi (\mathbf{x})$ 和 解码函数 $\mathbf{x}\approx \psi (\mathbf{z})$。实现在简化坐标和还原原始坐标之间的转化。

我们希望在对原数据进行 $\mathbf{z}=\phi (\mathbf{x})$ 的编码之后，能更方便我们使用公式(3)的形式来进行描述。其中，
基函数库 $\Theta (\mathbf{z})=\left[{\mathbf{\theta }}_1(\mathbf{z}),{\mathbf{\theta }}_2(\mathbf{z}),\dots ,{\mathbf{\theta }}_p(\mathbf{z})\right]$ 是由多项式或初等函数所组成。
稀疏系数向量 $\mathbf{\Xi }=\left[{\mathbf{\xi }}_1,\dots ,{\mathbf{\xi }}_d\right]$ 指的是其中的非零系数尽可能少，大部分 ${\mathbf{\xi }}_i=0$.

$$\frac{d}{dt}\mathbf{z}(t)=\mathbf{g}(\mathbf{z}(t))=\mathbf{\Theta }(\mathbf{z}(t))\mathbf{\Xi }\text{\hspace{1em}(3)}$$

(误差函数与编码器) 上式基函数库 $\Theta$ 是在训练之前由人类专家指定好的，而稀疏系数向量是在训练的过程中随之确定的。由于我们要求公式(3)是尽量成立的，再令 $\dot{\mathbf{z}}(t)={\nabla }_{\mathbf{x}}\phi (\mathbf{x}(t))\dot{\mathbf{x}}(t)$ ，可以得出我们误差函数 其实只是把公式(3)移项，让等式左右相减，尽量为0
$${\mathcal{L}}_{d\mathbf{z}/dt}={\Vert {\nabla }_{\mathbf{x}}\phi (\mathbf{x})\dot{\mathbf{x}}-\mathbf{\Theta }\left(\phi (\mathbf{x}{)}^T\right)\mathbf{\Xi }\Vert }_2^2\text{\hspace{1em}(4)}$$

(误差函数与解码器) 但公式(4)只与编码器 $\phi$ 有关，我们还希望得到正确解码器 $\psi$，要求解码后能够还原 $\dot{\mathbf{x}}$ 和 $\mathbf{x}$ 的时间序列。于是可以写出下面两个误差函数
$$\begin{array}{rl}{\mathcal{L}}_{d\mathbf{x}/dt}& ={\Vert \dot{\mathbf{x}}-\left({\nabla }_{\mathbf{z}}\psi (\phi (\mathbf{x}))\right)\left(\mathbf{\Theta }\left(\phi (\mathbf{x}{)}^T\right)\mathbf{\Xi }\right)\Vert }_2^2.\\ {\mathcal{L}}_{\text{recon }}& =\Vert \mathbf{x}-\psi (\phi (\mathbf{x})){\Vert }_2^2,\end{array}\text{\hspace{1em}(5,6)}$$

(误差函数与稀疏系数) 此外，我们还希望用于描述编码后系统，所使用的函数越少越好。于是我们对使用L1范数来进行约束。即 ${\mathcal{L}}_{\text{reg}}=||\mathbf{\Xi }|{|}_1$

备注：$\Vert x{\Vert }_0$ 表示非零元素的个数，是NP难题，一般用L1范数代替。
L1范数，$\Vert x{\Vert }_1={\sum }_{i=1}^n\left|x_i\right|$，一般用来进行稀疏优化。

(最终的误差函数) 现在，我们把上面误差函数汇总，得到了本文最终所使用的误差函数-公式(7)，其中 ${\lambda }_1,{\lambda }_2,{\lambda }_3$ 是超参数。

$${\mathcal{L}}_{\text{recon }}+{\lambda }_1{\mathcal{L}}_{d\mathbf{x}/dt}+{\lambda }_2{\mathcal{L}}_{d\mathbf{z}/dt}+{\lambda }_3{\mathcal{L}}_{\mathrm{reg}}\text{\hspace{1em}(7)}$$

除了使用 $L_1$ 正则化之外，作者还提出可以使用SINDy的顺序阈值方法来纳入训练过程。就是在训练的期间，固定的时间间隔内，低于某个阈值的系数会被强行固定为0，然后使用模型中剩余的项继续训练。

---

结果 Result

以Lorenz系统为例，原始系统为
$$\begin{array}{rl}& {\dot{z}}_1=\sigma \left(z_2-z_1\right)\\ & {\dot{z}}_2=z_1\left(\rho -z_3\right)-z_2\\ & {\dot{z}}_3=z_1{z}_2-\beta z_3\end{array}$$

假设我们采集的数据是由这三个状态变量的非线性组合形成的时间序列数据
$$\mathbf{x}(t)={\mathbf{u}}_1{z}_1(t)+{\mathbf{u}}_2{z}_2(t)+{\mathbf{u}}_3{z}_3(t)+{\mathbf{u}}_4{z}_1(t{)}^3+{\mathbf{u}}_5{z}_2(t{)}^3+{\mathbf{u}}_6{z}_3(t{)}^3.$$

我们看着这个组合的数据，去重构动力学方程。得到结果如下
$$\begin{array}{rl}& {\dot{z}}_1=-10.0z_1+10.0z_2\\ & {\dot{z}}_2=27.7z_1-0.9z_2-5.5z_1{z}_3\\ & {\dot{z}}_3=-2.7z_3+5.5z_1{z}_2\end{array}$$
可以验证还原出的系统与真实系统是基本等价的。

总结

回顾这篇文章开头的流程图，是否有加深理解呢？

相关链接

• SINDy源码
• 本文的论文原文
• Matlab - 非线性动力学/混沌系统/复杂性科学/系统科学常用工具
• Matlab - RK4的Lorenz仿真
• 论文泛读 - 基于RNN建模: Reconstructing computational system dynamics from neural data with recurrent neural