Linos Blog
Go back

从时间序列数据中发现混沌

来自Colorado大学的神秘算法

今天看到一份只有两页的讲义。它给出了一套从时间序列计算李雅普诺夫指数的算法：找一个邻居，跟着它跑，跑远了再换一个，然后把几段距离变化加起来处理一下，就得到了对李指数的估计。看完后我的第一反应是，这个算法步骤过于简单而显得有些不真实。

李雅普诺夫指数是通用的衡量混沌是否出现的黄金标准，当指数为正时，预示着动力系统存在混沌。从定义来讲，李指数说的是两条无限接近的轨迹怎样分开。在时间序列分析中，我们手里明明只有一串观测数据，算法到底从哪里弄来了第二条轨迹？而且它还会在中途把这条轨迹扔掉，再换一条新的。如此东拼西凑出来的东西，为什么能代表对初始轨迹的扰动？

要回答这些个问题，就得把讲义展开分析分析。从什么叫“一个状态”开始，一路经过动力系统、局部线性、离散采样和延迟嵌入，最后回到讲义中那张示意图。

运动由多个参数来加以描述

先从一个最简单的摆锤说起。

假设某一刻摆锤正好处于最低点。只知道这个位置，我们无法判断它下一刻会去哪里。它可能正在向左冲，也可能正在向右冲，还可能已经没有力气，准备停在这里。位置相同，未来却完全不同，是因为我们还不知道它的速度。 位置和速度放在一起，才足以描述这个摆锤当下的处境。动力系统里把这种能够决定后续演化的完整信息叫作“状态”。至于究竟要多少个数字才能写出一个状态，则取决于我们面对的系统。摆锤至少需要位置和速度；三维空间里飞行的刚体还需要更多位置、方向和转动信息。

摆锤的位置、速度与状态空间

图一：单独的位置会把两个不同状态压到一起；加上速度以后，它们才在状态空间中分开。

把系统所有可能的状态放在一起，就得到一张列出它所有可能处境的地图，也就是状态空间（也称相空间）。系统随时间变化，相当于一个点在这张地图上不断移动。这个点走过的路线，就是轨迹。

连续运动通常写成 dxdt=f(x)\frac{dx}{dt} = f(x)。系统走到状态 xx 时，规则 ff 告诉它接下来朝哪个方向走。离散系统则可以写成 xn+1=F(xn)x_{n+1} = F(x_n)，也就是拿现在的状态交给规则 FF，直接得到下一步的状态。这里要注意一个重点：xx 是个矢量，它可以由多个参数组成，所以才能表征一个系统的状态。 比如在我们举的摆锤这个例子中，x=(v,θ)x = (v, \theta)。

无论是连续的还是离散的微分方程，它们表征的都是“状态加规则决定未来”这个结构。只要这个结构存在，我们就可以沿着一条轨迹向前走，也可以比较两条起点略有差异的轨迹后来会变成什么样。

混沌真正混乱在什么地方

日常说一个东西“很混乱”，通常是说我们看不出规律。但动力系统所说的混沌反而可以由一套完全确定的规则产生。起点确定以后，后面每一步都照章办事，过程中不会跑出来个爱因斯坦拿着骰子临时决定路线。

那混沌它“混”在哪呢？究其根本是起点不可能被无限精确地知道。

取两个非常接近的初始状态，让它们遵守同一套规则。如果系统处在混沌状态，它们的距离会在一段时间内近似这样增长：

∥δx(t)∥≈∥δx0∥eλt\lVert \delta x(t) \rVert \approx \lVert \delta x_0 \rVert e^{\lambda t}

δx0\delta x_0 是最初那点小差异，δx(t)\delta x(t) 是经过时间 tt 后的差异。指数 λ\lambda 决定这点误差被放大的速度。比如 λ\lambda 为正，误差就会随时间指数增长；我们多知道几位小数，只是把预测精度丢失的时刻往后推一些，无法彻底解决问题。

因此可以这样定义单条轨迹的李雅普诺夫指数：

λ=lim⁡T→∞1Tln⁡(∥δx(T)∥∥δx0∥)\lambda = \lim_{T \to \infty} \frac{1}{T} \ln\left(\frac{\lVert \delta x(T) \rVert}{\lVert \delta x_0 \rVert}\right)

说白了，先看一小段初始误差在很长时间里总共被放大了多少倍，再取对数，把这个倍数摊回每一个单位时间。λ\lambda 小于零，邻近状态平均会靠拢；等于零，长期既没有指数收拢，也没有指数发散；大于零，则至少有一个方向会不断把微小差异拉开。

一个系统可以同时有许多方向，也就可能有一组不同的李雅普诺夫指数。我们往往只关心其中最大的一个。只要最初的扰动在最快伸长方向上有一点分量，反复演化以后，其他增长较慢的分量就会相形见绌，整个扰动向量看上去越来越像那个最快的方向。这就是讲义算法最终能得到最大指数的原因。

我们看到的只有一帧一帧的世界

真实运动可以连续发生（至少目前人类物理学是这么假定的），测量设备却只能隔一段时间记录一次。温度计每秒报一个数，相机每秒拍若干帧，我们能拿到的都是 x0,x1,x2,…x_0, x_1, x_2, \ldots 这样一格一格的数据。

假设两次采样之间相隔 Δt\Delta t，那么连续系统从一帧走到下一帧的完整演化可以记成 ΦΔt\Phi^{\Delta t}。给它换一个短名字 FF，就有：

F=ΦΔt,xn+1=F(xn)F = \Phi^{\Delta t}, \qquad x_{n+1} = F(x_n)

这一步常被含混地叫作“离散化”，但这里做的只是每隔 Δt\Delta t 给连续运动拍一张照片。

从这种一帧一帧的视角出发，一个微小的扰动从第零帧走到第一帧，会被这一段运动拉伸和旋转一次；从第一帧走到第二帧，又会被下一段运动处理一次。走过很多帧，总效果就是很多个局部变化连续相乘。然后因为对数能把这种难处理的连乘拆成逐段相加，我们在外面套一个对数。嗯，恭喜你再次发明了离散李雅普诺夫指数。

两条轨迹跑远以后，必须把它们拉回来

按照定义，我们似乎只要找两个非常接近的点，让它们一直往前走，最后量一下距离就行了。实际上，这个最直观的办法几乎一定会在靠后的迭代步骤中出问题。

混沌吸引子通常是有界的。可以把它想成一团不断被揉动的面团。动力学一边把面团拉长，一边又把它折回来。然而不在这个吸引子管辖范围内的面团，往往会很快地铺满整个宇宙。也就是说，有界性要求，两点在局部被拉开以后，迟早会跑到吸引子两个相隔较远的位置，但此时它们之间的距离最多只能和整个吸引子差不多大，不会继续保持指数增长。

更重要的是，李雅普诺夫指数讨论的是“无限接近”时的增长。两点跑远以后，距离开始反映吸引子的弯曲、折叠和整体尺寸，局部线性近似已经失效。拿这一段饱和距离继续平均，只会把结果越算越接近零。

数值算法的解决办法很像拿短尺量一条长路。尺子用到头以后，先记下这一段从 LiL_i 伸长到了 Li′L_i'，然后把扰动缩回一个很小的长度，继续量下一段。几段增长合起来时：

总伸长倍数=∏iLi′Li\text{总伸长倍数} = \prod_i \frac{L_i'}{L_i}

取对数以后，便有：

总对数伸长=∑iln⁡(Li′Li)\text{总对数伸长} = \sum_i \ln\left(\frac{L_i'}{L_i}\right)

这种反复缩短、逐段累计的思路早已用于已知方程的李雅普诺夫指数计算，Benettin 等人的方法就是经典来源。方程已知时，程序可以直接把扰动向量缩短到原来的长度。但讲义中算法面对的是观测数据，无法随手创造一个新状态，所以它只能从现有轨迹中另找一个足够近的点来接班。

一串数据里怎么凭空多出一条轨迹

现在回到开头的问题。实验只运行了一次，也没有人偷偷复制一台设备，再把初始状态拧动十亿分之一。第二条轨迹究竟从哪里来？

首先得意识到，我们测到的那一列数字往往连第一条完整轨迹都不是。传感器可能只记录摆锤的位置，却没有记录速度；也可能只测到一个化学反应里某种物质的浓度，而真实系统中还有许多变量同时变化。单独一个 sns_n 只是完整状态投在某根坐标轴上的影子。

解决办法是把同一个传感器最近几帧的记忆绑在一起：

Xn=(sn,sn−τ,…,sn−(m−1)τ)X_n = \left(s_n, s_{n-\tau}, \ldots, s_{n-(m-1)\tau}\right)

这里的 τ\tau 是取样之间的延迟，mm 是一共保留多少帧。摆锤现在处在什么位置，再加上它稍早以前处在什么位置，就已经透露了它大概朝哪里运动。保留更多合适的历史，便可能把原来挤在同一个观测值上的不同状态重新分开。

Takens 的延迟嵌入定理给这件事提供了理论基础。在一般的平滑观测、确定性系统和足够高的嵌入维数下，延迟坐标能在不把不同状态粘在一起的前提下重构原来的吸引子。后来 Sauer、Yorke 和 Casdagli 的 Embedology又把这套思想扩展到了更贴近混沌吸引子的情形。

重构出来的图形通常会被拉长、压扁和扭曲，但这种平滑换坐标只会把局部距离乘上有限倍数。一个有限倍数取对数后，再除以越来越长的观察时间，贡献会趋近于零。因此，在理想条件下，长期的李雅普诺夫指数能够被保留下来。

用了这种重构手段，第二条轨迹就出现了。混沌轨迹会在吸引子上长期游走，并在不同时间回到相近区域。我们把主轨迹在 t0t_0 时刻的状态记作 X(t0)X(t_0)，再去数据的其他时间寻找一个非常接近它的状态 Z(t0)Z(t_0)。从这两个时刻分别向后读取同样多的数据，就像让两个近邻初始状态同时向前演化。

这里必须排除时间上紧挨着 t0t_0 的点。昨天的我当然和今天的我很近，可这只说明时间是连续的，并没有提供一条独立的邻近轨迹。实际计算通常会设置一个时间隔离窗口，也常被叫作 Theiler window，只在窗口之外寻找邻居。

Wolf 算法其实是一场陪跑接力

到了这里，我们终于可以理解讲义上的步骤。

主轨迹 XX 是我们一路跟到底的运动轨迹，ZiZ_i 是第 ii 位陪跑轨迹。第 ii 次陪跑开始时，两者距离是 LiL_i；跑到局部线性近似快要失效时，距离变成 Li′L_i'。阈值 ε\varepsilon 决定一个陪跑轨迹什么时候应该下场，角度 θ\theta 则限制下一段陪跑轨迹和上一段方向之间可以差多远。

算法实际上是下面八步：

  1. 用延迟坐标把时间序列展开成状态空间轨迹。
  2. 在数据较前的位置选择主轨迹起点，给后面留下足够长的路。
  3. 排除时间上过近的点，寻找最近邻状态作为第一段陪跑轨迹起点。
  4. 让主轨迹和陪跑轨迹同步向前移动，并持续计算距离。
  5. 距离超过 ε\varepsilon 时，记录这一段的 ln⁡(Li′/Li)\ln(L_i'/L_i)。
  6. 在主轨迹当前位置寻找新的近邻，并让新方向尽量贴近上一段的分离方向。
  7. 一路接力到数据末尾；找不到合格邻居时就停止并报告。
  8. 用所有对数伸长之和除以实际经过的总时间。

最终估计式是：

λ^=∑iln⁡(Li′/Li)∑iΔTi\hat{\lambda} = \frac{\sum_i \ln\left(L_i'/L_i\right)}{\sum_i \Delta T_i}

讲义使用了 log⁡2\log_2。这会把结果写成每单位时间经历多少次“距离翻倍”；使用自然对数 ln⁡\ln，则得到文献中更常见的指数单位。底数改变不会导致量纲改变，也不影响系统发散还是收敛。

最容易被略过的是方向条件。扰动经过一段时间以后，会越来越靠近最快伸长的方向。若换陪跑轨迹时只顾距离近，随手挑一段从完全不同方向出发的陪跑轨迹，算法就会反复丢掉这段已经形成的方向。讲义第二页画出的角度圆锥，目的正是让新陪跑轨迹“保留”一些上一段的方向记忆。距离和方向往往无法同时达到完美，只能在足够近与足够同向之间取一个实际可用的平衡。

讲义算法示意图

图二：讲义中的算法示意图

这场接力跑为什么能算出最大指数

两条轨迹足够接近时，它们的差可以近似成主轨迹旁边的一根切向扰动。动力学中经常提到的局部线性，就是在这段很小的尺度上，真实动力学对两条轨迹差值的作用，可以用一个线性变换近似。于是观测两条近邻轨迹怎样分开，就近似于观测一根无穷小扰动怎样增长。

延迟嵌入又保证，在理想条件下，重构空间是原吸引子的平滑改写。它会改变我们看到的形状，却不会改变长期指数增长率。这样一来，从重构轨迹中找到的近邻，才有资格充当原系统里的近邻状态。

每次距离达到 ε\varepsilon 就换轨迹，是为了把测量留在局部线性范围。换陪跑轨迹时保留大致方向，是为了一段又一段地继续追踪同一个最快伸长方向。这个连接当然不可能像已知方程时的向量缩放一样精确，所以 Wolf 算法得到的是一个估计值，但数学上它最终逼近于目标。

最后，各段伸长倍数原本应该相乘。对数把乘法变成求和，除以总时间以后，恰好回到“长期平均对数增长率”这个定义。只要初始扰动没有极其巧合地完全避开最快伸长方向，长期演化就会让最快分量占据主导，因此结果落在最大的李雅普诺夫指数上。

Wolf 算法没有真的回到过去改动初始条件。它在一条足够长的轨迹中找到另一个曾经出现过的相近状态，用这个状态借来一根扰动向量。等这根向量长得太大，就记下来，再从数据里借一根方向相近的新向量。只要每一段都够短、方向接续得够好，整场接力测到的就是那根无限小扰动长期以来的增长速度。

算法的前提条件

到了这一步，我们可以说算法在数学想法上有来历，却还不能看到程序吐出一个正数就宣布发现了混沌。我们得再检查下数据和算法是否满足了算法成立的前提。可能导致最终结果不准确的地方如下：

实际判断时应当改变 τ,m,ε,θ\tau, m, \varepsilon, \theta 和时间隔离窗口，观察结果是否在一段合理区域内保持稳定。数据长度增加以后，估计值也应该逐步收敛。若参数稍微一动，答案就从正数跳到负数，那结果数字显然就没法表征混沌与否。

更深入的分析，可以去看看Abarbanel 等人的综述，他们专门提醒了轨迹追踪法对数据量、邻居方向和噪声的敏感；Bradley的文章也强调，有限精度、有限数据和噪声会限制从时间序列恢复动力学性质的能力。另外，随机过程同样可能在有限尺度上制造看似漂亮的发散曲线，所以一个孤立的正指数估计数值只能成为证据的一部分。

参考资料

Share this post:
Share this post via email
Previous Post
评‘广东女孩诬告亲父猥亵’
Next Post
和Rex讨论新保守主义