NIS+

NIS+（Neural Information Squeezer Plus）是一个机器学习框架，可学习宏观动力学，并量化因果涌现的程度。该框架通过最大化有效信息，得到一个宏观动态模型。研究者通过在模拟数据和实际数据上的实验，证明了该框架的有效性。且该框架在不同测试环境里表现出更强的泛化能力。

问题与背景

在自然界和人类社会中，存在着许多由无数相互作用的元素构成的复杂系统，如气候系统、生态系统、鸟群、蚁群、细胞和大脑等。这些系统展现出丰富的非线性动力学行为，如果我们仅关注微观的尺度，会发现它们的行为非常复杂且难以预测。当我们从更宏观的尺度观察这些系统时，可以用更加简洁的规律来解释和预测这些系统，这便是复杂系统独有的涌现现象。

涌现现象与因果关系紧密相连。一方面，涌现现象是系统内部各元素间复杂非线性相互作用的结果；另一方面，这些涌现的特性又会对系统中的个体产生影响。因果涌现这个概念最早由Erik Hoel正式提出并定义，即因果涌现描述了宏观层面相对于微观层面在因果效应上的增强，这揭示了宏观与微观状态之间的差异和联系。

因果涌现的概念不仅将因果推理与复杂系统的涌现特性相结合，而且为解决一系列哲学问题提供了一个定量化的视角。

目前，关于如何定义因果涌现，有四个主要代表，分别是：①Hoel等基于粗粒化方法的因果涌现理论、②Rosas等基于信息分解的因果涌现理论、③张江等人基于奇异值分解的因果涌现理论、④Barnett等的动力学解耦。

具体内容可以参考词条因果涌现。

因果涌现识别

现实中，当我们要应用涌现的量化框架时，会出现微观状态的数据难以获取、粗粒化策略不好确定的问题，所以，我们需要找到从可观测数据中直接识别因果涌现的方法。

识别出复杂系统中的因果涌现，一方面可以节约预测成本，增强预测效果，尤其是泛化效果；一方面也可以诱导或者预防涌现的发生。在实际应用中，一些微观数据难以捕捉，我们往往只能收集到观测数据，并且，无法得到系统的真实动力学。因此，从可观测数据中识别系统中因果涌现的发生是一个关键的问题。近年来，基于神经网络的机器学习方法取得了突破性进展，借助此方法，以数据驱动的方式自主发现复杂系统的因果关系甚至动力学成为可能。此外，机器学习和神经网络还可以帮助我们找到更好的粗粒化策略、从数据中发现宏观层面的因果关系。

下面介绍两种因果涌现的识别方法：①基于信息分解的因果涌现识别近似方法、②神经信息压缩器方法系列（Neural Information Squeezer，简称NIS）。

基于信息分解的因果涌现识别

Rosas等学者通过信息分解框架给出了和Hoel等人不同的对因果涌现的新定义，并基于此识别量化因果涌现。但是信息分解框架中定义的信息原子难以计算，所以作者推导出只需要计算互信息的近似公式，提出了判定因果涌现发生的充分条件，即[math]\displaystyle{ \Psi_{t, t+1}(V) }[/math]，具体公式如下：

[math]\displaystyle{ \Psi_{t, t+1}(V):=I\left(V_t ; V_{t+1}\right)-\sum_j I\left(X_t^j ; V_{t+1}\right) }[/math]

式中，[math]\displaystyle{ V }[/math]是给定的宏观状态，[math]\displaystyle{ V_t ; V_{t+1} }[/math]代表两个连续时间的宏观状态变量，[math]\displaystyle{ X_t^j }[/math]表示第 j 维t时刻的微观变量。

当[math]\displaystyle{ \Psi_{t, t+1}(V) \gt 0 }[/math]，系统发生因果涌现。但当[math]\displaystyle{ \mathrm{\Psi}\lt 0 }[/math]，我们不能确定系统是否发生因果涌现。

需要指出的是，Hoel的定义和识别方法基于Judea Pearl因果，而此方法是基于格兰杰因果，利用机器学习框架，计算互信息的组合，没有引入do干预。

该方法因为是格兰杰因果，所以计算比较方便，不需要找到底层的动力学机制。且对系统的动力学没有马尔可夫性的假设和要求。但是也存在一些缺点：

1）该方法只是基于互信息计算，且得到的仅仅是发生因果涌现的近似的充分条件；

2）该方法需要手动给定粗粒化策略和宏观变量，而不同的选择会对结果造成显著影响；

3）高维系统中，[math]\displaystyle{ \Psi }[/math]作为近似条件，误差非常大，很容易得到负值，从而无法判断是否有因果涌现发生。

为了能够自动找到最佳的粗粒化策略，这套理论框架也发展出了相应的机器学习方法。Kaplanis等人基于机器学习的方法学习宏观态[math]\displaystyle{ V }[/math]以及最大化[math]\displaystyle{ \mathrm{\Psi} }[/math]：使用神经网络来学习将微观输入粗粒化成宏观输出，同时使用两个神经网络来分别学习互信息的计算,最后通过最大化两者之间的差(即[math]\displaystyle{ \mathrm{\Psi} }[/math])来优化学习。

NIS系列

张江等学者尝试基于神经网络和数据驱动提出了一种方法，能从时间序列数据中识别系统中的因果涌现，并且自动提取有效的粗粒化策略和宏观动力学，即神经信息压缩器（Neural Information Squeezer，简称NIS）。

模型由编码器（encoder）、动力学学习器（[math]\displaystyle{ f }[/math]）以及解码器（decoder）三个部分构成，编码器和解码器主要由可逆神经网络（Invertible Neural Network，简称INN）构建，动力学学习器由多层感知机（Multilayer Perceptron，简称MLP）构建。此模型框架可以看成是一个神经信息压缩器，将包含噪音的微观态压缩成宏观态，丢弃无用的信息，从而使得宏观动力学的因果性更强。NIS方法的模型框架如右图所示。

图中，模型输入是微观状态[math]\displaystyle{ X_t\ (X_t^1,X_t^2,…,X_t^p ) }[/math]，[math]\displaystyle{ p }[/math]表示输入数据的维数。[math]\displaystyle{ ϕ }[/math]是粗粒化函数（编码器），将输入的[math]\displaystyle{ p }[/math]维数据映射到[math]\displaystyle{ q }[/math]维数据上，得到宏观变量[math]\displaystyle{ Y_t }[/math]，此步丢失[math]\displaystyle{ p-q }[/math]维信息。[math]\displaystyle{ f }[/math]是动力学学习器，在宏观层面上学习有效的马尔可夫动力学。[math]\displaystyle{ \hat{y}_{t+1} }[/math]是预测的t+1时刻的宏观状态。由于此时数据是[math]\displaystyle{ q }[/math]维的，为了使用反粗粒化函数[math]\displaystyle{ ϕ^† }[/math]（解码器），需要用[math]\displaystyle{ p-q }[/math]维高斯随机向量填充数据。宏观变量经过反粗粒化函数之后可以得到预测的微观变量[math]\displaystyle{ \hat{x}_{t+1} }[/math]。而[math]\displaystyle{ x_t+1 }[/math]和[math]\displaystyle{ \hat{x}_{t+1} }[/math]之间的差值即为损失。

该方法的目标函数是希望在保证微观状态预测误差很小的条件下最大化有效信息，保证预测误差约束足够小是为了避免平凡解（粗粒化丢失太多信息，以致宏观动力学无效）的出现。最终希望得到有效的粗粒化维度[math]\displaystyle{ q }[/math]、粗粒化策略[math]\displaystyle{ \mathrm{\phi}_q }[/math]和宏观动力学[math]\displaystyle{ {\hat{f}}_{\mathrm{\phi}_q} }[/math]。

具体内容请参考NIS。

机器学习领域的分布外泛化问题

目前，机器学习需要在训练数据和测试数据是独立同分布（Independent and Identically Distributed，简称IID）的基础之上，这是一种理想化的假设。现实应用中，几乎不能满足IID假设，所以模型的性能都会有不同程度的下降。

以一个狗识别的图像判别任务为例。训练图像的背景一般是在草地上，可以看见整个狗。经过大量数据训练后，若给模型一个在草地上的狗的图像，那模型的正确率可以在95%以上；若给一张在水泥地上，且遮挡了狗部分身体的图像，那么模型的正确率可能是一半。

机器学习经过训练学习得到的很大程度上是变量之间的相关关系，这种相关关系可能来自样本选择偏差（Sample Selection Bias）。当考虑了其他变量后，原本看似有相关性的变量可能实际上并无关联，是一种虚假的相关关系。在上面的狗识别的例子中，机器学习学到的可能是根据草地和狗鼻子这两个因素来判断是不是狗，但是有因果关系的是狗鼻子，草地只是由于样本的分布不均匀，这就导致了模型不能有效识别出在水里、沙漠里的狗。简单来说，草地和标签之间是强相关、弱因果，狗鼻子和标签之间是强相关、强因果。

而在动力学系统中，模型可能对初始条件非常敏感，即使是微小的初始条件变化也可能导致系统行为的显著差异。如果训练数据的初始条件分布与测试数据不同，模型可能无法准确预测测试数据下的系统行为。我们希望的是，模型可以透过初始条件的差异，抓住系统动力学的本质。甚至有时候，我们期待模型可以泛化到某些动力学参数不同的系统测试数据上，依旧有不错的表现。在实际应用中，我们希望学习到的是因果关系（Causation），结果更加稳定且具有可解释性。

分布外泛化（Out-of-Distribution Generalization，简称OOD Generalization）是指模型在面对与训练数据分布不同的新数据时，仍然能够保持较好的性能，涉及到模型对未知或未见情况的适应能力，削弱样本选择偏差的情况。此问题通常表现为模型在训练集上运行良好，但在面对分布外的各种数据时，模型性能有所下降。造成这个问题的原因可能是过度拟合训练数据，导致泛化能力差；也可能是由于训练数据无法覆盖所有可能的输入情况，特别是对于罕见或极端的情况。

为了增强分布外泛化能力，学者可以生成多样化的数据，模拟不同的测试环境，还可以通过域适应技术，使模型可以适应不同的测试数据分布，还有不变性学习、元学习等方法。

神经信息压缩机（NIS）介绍

概述

假设复杂动态系统的行为数据是时间序列 [math]\displaystyle{ \{x_t\} }[/math]，时间步长[math]\displaystyle{ t = 1,2，…， T }[/math]，维数是p，它们构成了可观测的微观状态，我们假设不存在未观测变量。一个粗粒化策略（编码器） [math]\displaystyle{ ϕ: R_p → R_q }[/math] ，其中 [math]\displaystyle{ q ≤ p }[/math]， [math]\displaystyle{ q }[/math]是宏观状态的维度，作为超参数给定；一个相应的反粗粒化策略（解码器）[math]\displaystyle{ ϕ^\dagger: R_q → R_p }[/math]，以及一个宏观层面的马尔可夫动力学（动力学学习器）[math]\displaystyle{ f_q }[/math] ，使得[math]\displaystyle{ f_q }[/math]的有效信息（[math]\displaystyle{ \mathcal{J} }[/math]）值在通过[math]\displaystyle{ ϕ }[/math]、[math]\displaystyle{ f_q }[/math]和[math]\displaystyle{ ϕ^\dagger }[/math]预测出的[math]\displaystyle{ x_{t+1} }[/math]与[math]\displaystyle{ x_{t+1} }[/math]的实际数据的差距最小的约束下最大化。[math]\displaystyle{ \epsilon }[/math]是给定的常数。它们的关系用方程表示为：

[math]\displaystyle{ \begin{aligned}&\max_{\phi,f_{q},\phi^{+}}\mathcal{J}(f_{q}),\\&s.t.\begin{cases}\parallel\hat{x}_{t+1}-x_{t+1}\parallel\lt \epsilon,\\\hat{x}_{t+1}=\phi^{\dagger}\left(f_{q}(\phi(x_{t})\bigr)\right).\end{cases}\end{aligned} }[/math]

(1)

为了数学性质和可解释性，以及降低模型参数量，采用了可逆神经网络。作者将编码过程分解为了两个步骤：

①双射映射，[math]\displaystyle{ ψ: R_p → R_p }[/math]，此步无信息丢失，由可逆神经网络实现。

②投影运算，[math]\displaystyle{ \chi_q }[/math]，此步将输入的[math]\displaystyle{ p }[/math]维数据映射到[math]\displaystyle{ q }[/math]维数据上，得到宏观变量[math]\displaystyle{ Y_t }[/math]，此步丢失[math]\displaystyle{ p-q }[/math]维信息。

[math]\displaystyle{ y = \phi(x) = \chi_q (\psi(x)) }[/math]

解码器也分为两个步骤：

①用[math]\displaystyle{ p-q }[/math]维高斯随机向量填充数据得到[math]\displaystyle{ p }[/math]维数据。式中，ξ为p−q维的标准高斯随机向量，⊕表示向量拼接操作。

②使用反粗粒化函数得到预测的微观变量[math]\displaystyle{ \hat{x}_{t+1} }[/math]。由于使用可逆神经网络，此步和编码器中的双射映射共享参数。

[math]\displaystyle{ \phi^{\dagger}(x)=\psi_{\omega}^{-1}(x\oplus\xi) }[/math]

在编码器一步，虽然有丢失信息，但是通过数学证明可以得知，当神经网络训练收敛时，数据中相邻宏观状态之间的互信息倾向于接近两个相邻时间步的微观态互信息，丢失的信息与预测几乎无关，所有有效信息都汇聚到动力学学习器的信息瓶颈之中。且在此框架下，信息瓶颈的应用更加明显、更具可解释性。

然而由于该目标函数是一个泛函优化问题，往往很难优化。为了解决这个问题，将优化过程分为两个阶段。

第一个阶段：在给定宏观尺度[math]\displaystyle{ q }[/math]的情况下，最小化预测误差[math]\displaystyle{ \left\|\phi_q^{\dagger}(Y(t+1))-X_{t+1}\right\| }[/math]。可以采用随机梯度下降技术，通过训练神经网络解决。在此步，我们可以得到依赖于q的有效粗粒化策略和宏观动力学。

第二个阶段：由于参数 [math]\displaystyle{ q }[/math] 是一维的，且 [math]\displaystyle{ 0 \lt q \lt p }[/math]，简单迭代找到有效信息最大的宏观尺度。此步将复杂的泛函问题转化成一维空间中线性搜索问题。

具体请阅读词条NIS。

神经网络上的EI计算

如果神经网络的输入[math]\displaystyle{ X=(x_1,x_2,\cdot\cdot\cdot,x_n)\in [-L,L]^n }[/math]，则表示X定义在一个大小为L的超立方体上，其中L是一个非常大的整数。输出为[math]\displaystyle{ Y=(y_1,y_2,\cdot\cdot\cdot,y_m) }[/math]， [math]\displaystyle{ Y=\mu(X) }[/math]。这里µ是神经网络实现的确定性映射:[math]\displaystyle{ \mu: \mathcal{R}^n\rightarrow \mathcal{R}^m }[/math]，它在X处的雅可比矩阵是[math]\displaystyle{ \partial_{X'} \mu(X)\equiv \left\{\frac{\partial \mu_i(X')}{\partial X'_j}\left|_{X'=X}\right.\right\}_{nm} }[/math]。如果神经网络可以看作是给定X条件下的高斯分布，则神经网络的有效信息(EI)可以用以下方法计算:

[math]\displaystyle{ EI_L(μ)=I(do(X\sim U([-L,L]^{n};Y) ≈ -\frac{m+m\ln(2\pi)+\sum_{i=1}^m\sigma_i^2/m}{2}+n\ln(2L)+\mathbf{E}_{X\sim U([-L,L]^n} \left(\ln |\det(\partial_{X'} f(X))|\right) }[/math]

(2)

式中，[math]\displaystyle{ \Sigma=diag(\sigma_1^2,\sigma_2^2,\cdot\cdot\cdot,\sigma_m^2) }[/math] 是协方差矩阵， [math]\displaystyle{ \sigma_i }[/math]是输出[math]\displaystyle{ y_i }[/math]的标准差，可由[math]\displaystyle{ y_i }[/math]的均方误差估计得到，[math]\displaystyle{ U([-L,L]^n) }[/math]为[math]\displaystyle{ [-L,L]^n }[/math]上的均匀分布，[math]\displaystyle{ | · | }[/math]为绝对值，det为行列式。如果[math]\displaystyle{ \det(\partial_{X'} \mu(X))\equiv 0 }[/math]对于所有X，则令EI≈0。

但是，此公式不能直接应用于实际情况，因为它会随着输入n或输出m的维数增加而增加。解决这一问题的方法是通过除以输入维数来定义维度平均有效信息(dEI)，记为：

[math]\displaystyle{ \mathcal{J}_L = \frac{EI_L(\mu)}{n} }[/math]

当输入输出数相同时(m = n)，则：

[math]\displaystyle{ \mathcal{J}_L(\mu)=-\frac{1+\ln(2\pi)+\sum_{i=1}^n\sigma_i^2/n}{2}+\ln(2L)+\frac{1}{n}\mathbf{E}_{X\sim U([-L,L]^n} \left(\ln |\det(\partial_{X'} f(X))|\right) }[/math]

不过，此式子依赖于参数L，输入数据的域范围。为了消除L的影响，作者计算维度平均CE。对于维度q的宏观动力学[math]\displaystyle{ f_M }[/math]和维度p的微观动力学[math]\displaystyle{ f_m }[/math]，我们将维度平均CE定义为:

[math]\displaystyle{ \Delta\mathcal{J}_L(f_M,f_m)=\mathcal{J}_L(f_M)-\mathcal{J}_L(f_m)=\frac{EI_L(f_M)}{q}-\frac{EI_L(f_m)}{p} }[/math]

如果[math]\displaystyle{ f_M }[/math]和[math]\displaystyle{ f_m }[/math]由q维的[math]\displaystyle{ µ_M }[/math]和p维的[math]\displaystyle{ µ_m }[/math]的神经网络参数化，则：

[math]\displaystyle{ \Delta\mathcal{J}=\left(\frac{1}{q}\mathbf{E}_{X_M}\ln|\det \partial_{X_M}\mu_M|-\frac{1}{p}\mathbf{E}_{X_m}\ln|\det \partial_{X_m}\mu_m|\right)-\left(\frac{1}{q}\sum_{i=1}^{q}\ln\sigma_{i,M}^2-\frac{1}{p}\sum_{i=1}^{p}\ln\sigma_{i,m}^2\right) }[/math]

式中，其中[math]\displaystyle{ \sigma_{i,M} }[/math]和[math]\displaystyle{ \sigma_{i,m} }[/math]分别是[math]\displaystyle{ \mu_M }[/math]和[math]\displaystyle{ \mu_m }[/math]在第i维上的标准差。

此时，完全消除了输入或输出维度以及参数L的影响，使其成为一个更可靠的指标。从数据中识别因果涌现的结果只受到雅可比矩阵的方差的相对值和行列式的对数值的影响。

NIS的缺陷

此模型在弹簧振子、马尔科夫链、布尔网络上都进行了实验，实验效果理想，表明此模型的正确性、有效性。

但是此模型还有一些缺点和尚未改进的地方：

1）此框架目前只能在小数据集上使用，难以拓展到更大规模的复杂系统上。

2）虽然我们可以说明什么是粗粒度函数，并将其清晰地分解为信息转换和信息丢弃两个部分，但变量的分组方法隐式编码在可逆神经网络中，宏观动力学学习器的机制也没有完全透明。框架仍然缺乏可解释性，

3）模型需要基于马尔科夫性的假设。

4）此方法的数学形式是一个泛函问题，无法直接进行优化，虽然可以通过两个阶段得到结果，但是NIS没有真正地最大化有效信息。在NIS+中，学者将通过计算变分下界解决泛函问题。

具体请阅读词条NIS。

强化版神经信息压缩机（NIS+）

概述

为了最大化式1中定义的EI，我们将NIS的框架扩展为NIS+。在NIS+中，我们首先使用互信息和变分不等式的公式将互信息的最大化问题转化为机器学习问题，其次，使用[math]\displaystyle{ y_{t+1}=\phi(x_{t+1}) }[/math]来预测[math]\displaystyle{ y_{t} }[/math]，从而保证互信息最大化。最后，利用样本重加权技术来解决均匀分布干预的挑战，从而优化EI。所有这些技术组成了增强版神经信息压缩机(NIS+)。在此框架中，输入可观测的数据，输出是因果涌现的程度、宏观动力学、涌现斑图以及粗粒化策略。

数学推导

式1中，数学形式是一个泛函问题，无法直接进行优化，学者将通过计算并优化变分下界来解决泛函优化问题。同时，在NIS+框架中，学者使用了编码器将p维的输入数据进行粗粒化，得到q维的宏观数据，下面编码器的通用逼近定理将证明编码器的可以近似任意复杂的粗粒化函数。

在此章节，我们将使用大写字母来表示相应的随机变量。

例如，[math]\displaystyle{ X_{t} }[/math]表示时间t的微观状态[math]\displaystyle{ x_{t} }[/math]的随机变量，[math]\displaystyle{ Y_t+1 }[/math]表示时间t+1的宏观状态[math]\displaystyle{ y_t+1 }[/math]对应的随机变量。对于任意随机变量X， [math]\displaystyle{ \tilde{X} }[/math]表示干预后的同一随机变量X。[math]\displaystyle{ \hat{X} }[/math]表示神经网络对X的预测。

宏观EI的变分下界

原始的有约束的目标优化公式如式1所示。

在此方程中[math]\displaystyle{ \hat{X}_{t+1}=\psi_{\omega}^{-1}(\hat{Y}_{t+1}\bigoplus \xi) }[/math]，其中[math]\displaystyle{ \psi_{\omega}^{-1} }[/math]是可逆映射，根据引理1和引理2以及互信息的性质，我们可以得到：

[math]\displaystyle{ I(Y_t,\hat{Y}_{t+1})=I(Y_t,\hat{X}_{t+1})=H(Y_t)-H(Y_t|\hat{X}_{t+1}) }[/math]

令[math]\displaystyle{ U_q }[/math]为宏观变量的均匀分布，那么[math]\displaystyle{ H(\tilde{Y}_t)=H(U_q) }[/math]。故：

[math]\displaystyle{ \mathcal{J}(f_{\theta,q})=I(\tilde{Y}_{t},\hat{X}_{t+1})=H(U_q)-H(\tilde{Y}_t|\tilde{X}_{t+1}) }[/math]

由于[math]\displaystyle{ H(U_q) }[/math]为常数，所以[math]\displaystyle{ \mathcal{J}(f_{q}) }[/math]的优化可以转化为条件熵[math]\displaystyle{ H(\tilde{Y}_t|\tilde{X}_{t+1}) }[/math]的最小化优化问题。根据引理3可以得知：

[math]\displaystyle{ H(\tilde{Y}_t|\tilde{X}_{t+1}) \le -\iint \tilde{p}(\boldsymbol{y}_t, \boldsymbol{x}_{t+1})\ln g(\boldsymbol{y}_t|\boldsymbol{x}_{t+1}) \mathrm{d}\boldsymbol{y}_t \mathrm{d}\boldsymbol{x}_{t+1} = -\iint \tilde{p}(\boldsymbol{y}_{t})\tilde{p}(\boldsymbol{x}_{t+1}|\boldsymbol{y}_t)\ln g(\boldsymbol{y}_t|\boldsymbol{x}_{t+1}) \mathrm{d}\boldsymbol{y}_t \mathrm{d}\boldsymbol{x}_{t+1} }[/math]

式中，[math]\displaystyle{ \tilde{p} }[/math]表示随机变量[math]\displaystyle{ Y_t }[/math]被干预情况下的概率分布函数。

作者使用神经网络来拟合分布[math]\displaystyle{ g(y_t|x_t+1) }[/math]，根据引理3，[math]\displaystyle{ g(y_t|x_t+1) }[/math]可以是任何分布，在这里，假设[math]\displaystyle{ g(y_t|x_t+1) }[/math]为正态分布，即[math]\displaystyle{ g(\boldsymbol{y}_t|\boldsymbol{x}_{t+1})\sim N(\mu,\Sigma) }[/math]，其中[math]\displaystyle{ \mu =g_{\theta'}(\phi(\boldsymbol{x}_{t+1})) }[/math]，[math]\displaystyle{ \Sigma=diag(\sigma_1, \sigma_2,\cdot\cdot\cdot,\sigma_q) }[/math]是常数对角矩阵，进一步，假设[math]\displaystyle{ \sigma_i }[/math]是有界的，则[math]\displaystyle{ \sigma_i\in[\sigma_m,\sigma_M] }[/math]，其中[math]\displaystyle{ \sigma_m }[/math]和[math]\displaystyle{ \sigma_M }[/math]分别是MSE的最小值和最大值。则[math]\displaystyle{ g(y_t|x_t+1) }[/math]的对数概率密度函数为：

[math]\displaystyle{ \ln g(\boldsymbol{y}_t|\boldsymbol{x}_{t+1})\approx \ln \frac{1}{(2\pi)^{\frac{m}{2}}|\Sigma|^\frac{1}{2}} e^{-\frac{(\boldsymbol{y}_t-g_{\theta'}(\phi(\boldsymbol{x}_{t+1})))^2}{2|\Sigma|}} = -\frac{(\boldsymbol{y}_t-g_{\theta'}(\phi(\boldsymbol{x}_{t+1})))^2}{2|\Sigma|}+\ln \frac{1}{(2\pi)^{\frac{m}{2}}|\Sigma|^\frac{1}{2}} ≥ -\frac{(\boldsymbol{y}_t-g_{\theta'}(\phi(\boldsymbol{x}_{t+1})))^2}{2|\Sigma|}+\ln \frac{1}{(2\pi)^{\frac{m}{2}}|\Sigma|_{max}^\frac{1}{2}} }[/math]

如果训练足够充分，那么[math]\displaystyle{ \tilde{p}(\boldsymbol{x}_{t+1}|\boldsymbol{y}_t)\approx p(\boldsymbol{x}_{t+1}|\boldsymbol{y}_t) }[/math]，故：

[math]\displaystyle{ H(\tilde{Y}_t|\tilde{X}_{t+1}) \le \iint \tilde{p}(\boldsymbol{y}_{t})p(\boldsymbol{x}_{t+1}|\boldsymbol{y}_t)\left[\frac{(\phi (\boldsymbol{x}_t)-g_{\theta'}(\phi(\boldsymbol{x}_{t+1})))^2}{2|\Sigma|_{min}}-\ln \frac{1}{(2\pi)^{\frac{m}{2}}|\Sigma|_{max}^\frac{1}{2}}\right] \mathrm{d}\boldsymbol{y}_t \mathrm{d}\boldsymbol{x}_{t+1} }[/math]

式中，[math]\displaystyle{ |\Sigma|_{min}=\sigma_{min}^q }[/math]，[math]\displaystyle{ |\Sigma|_{max}=\sigma_{max}^q }[/math]。

由于[math]\displaystyle{ \tilde{p}(\boldsymbol{y}_{t})p(\boldsymbol{x}_{t+1}|\boldsymbol{y}_t)=\frac{\tilde{p}(\boldsymbol{y}_{t})}{p(\boldsymbol{y}_{t})}p(\boldsymbol{x}_{t+1},\boldsymbol{y}_t) }[/math]，其中[math]\displaystyle{ \tilde{p}(\boldsymbol{y}_{t}) }[/math]是目标分布，[math]\displaystyle{ p(\boldsymbol{y}_{t}) }[/math]是数据的原始分布。我们定义逆概率权重[math]\displaystyle{ w(\boldsymbol{x}_t) \equiv \frac{\tilde{p}(\boldsymbol{y}_{t})}{p(\boldsymbol{y}_{t})} }[/math]，[math]\displaystyle{ z=\frac{(\phi (\boldsymbol{x}_t)-g_{\theta'}(\phi(\boldsymbol{x}_{t+1})))^2}{2|\Sigma|_{min}}-\ln \frac{1}{(2\pi)^{\frac{m}{2}}|\Sigma|_{max}^\frac{1}{2}} }[/math]。且由于我们使用离散样本[math]\displaystyle{ {x_t} }[/math]来训练神经网络，所以我们可以使用样本均值作为期望的近似估计。因此，[math]\displaystyle{ H(\tilde{Y}_t|\tilde{X}_{t+1}) }[/math]的变分上界可写成:

[math]\displaystyle{ H(\tilde{Y}_t|\tilde{X}_{t+1}) \le \iint w(\boldsymbol{x}_t)p(\boldsymbol{x}_{t+1},\boldsymbol{y}_t)z \mathrm{d}\boldsymbol{y}_t \mathrm{d}\boldsymbol{x}_{t+1} \approx \frac{1}{T}\sum_{i=0}^{T-1}w(\boldsymbol{x}_t)z }[/math]

进而，我们可以得到原目标函数的变分下界：

[math]\displaystyle{ \mathcal{J}(f_{\theta,q}) \ge H(U_q)-\frac{1}{T}\sum_{i=0}^{T-1}w(\boldsymbol{x}_t)z }[/math]

原优化问题（式1）便转换成：

[math]\displaystyle{ \min_{\omega,\theta,\theta'} \sum_{i=0}^{T-1}w(\boldsymbol{x}_t)|\phi(\boldsymbol{x}_t)-g_{\theta'}(\phi(\boldsymbol{x}_{t+1}))|^2 }[/math]

[math]\displaystyle{ s.t. || \hat{\boldsymbol{x}}_{t+1}-\boldsymbol{x}_{t+1} || \lt \epsilon }[/math]

式中，[math]\displaystyle{ \omega }[/math]，[math]\displaystyle{ \theta }[/math]， [math]\displaystyle{ \theta' }[/math]分别表示NIS+框架中三个神经网络[math]\displaystyle{ \psi }[/math]， [math]\displaystyle{ f_{\theta} }[/math]， [math]\displaystyle{ g_{\theta'} }[/math]的参数。

然后构造拉格朗日函数：

[math]\displaystyle{ L(\omega,\theta,\theta',\lambda)=\sum_{i=0}^{T-1}w(\boldsymbol{x}_t)|\phi(\boldsymbol{x}_t)-g_{\theta'}(\phi(\boldsymbol{x}_{t+1}))|^2+\lambda|| \phi^\dagger(\boldsymbol{y}_{t+1})-\boldsymbol{x}_{t+1} || }[/math]

优化目标（式1）便转化为：

[math]\displaystyle{ \begin{aligned}&\min_{f,g,\phi,\phi\dagger}\sum_{t=1}^{T-1}w(\boldsymbol{x}_t)\parallel\boldsymbol{y}_t-g(\boldsymbol{y}_{t+1})\parallel+\lambda\parallel\hat{x}_{t+1}-\boldsymbol{x}_{t+1}\parallel,\\&s.t.\begin{cases}y_{t}=\phi(x_{t}),\\\hat{y}_{t+1}=f(y_t),\\\hat{x}_{t+1}=\phi^{\dagger}\left(f_{q}(\phi(x_{t})\bigr)\right),\\y_{t+1}=\phi(x_{t+1}).\end{cases}\end{aligned} }[/math]

(3)

式中，[math]\displaystyle{ x_{t} }[/math]、[math]\displaystyle{ x_{t+1} }[/math]表示可观测的微观数据，[math]\displaystyle{ y_{t} }[/math]、[math]\displaystyle{ y_{t+1} }[/math]表示经过粗粒化函数[math]\displaystyle{ \phi }[/math]得到的宏观数据，[math]\displaystyle{ \hat{y}_{t+1} }[/math]表示[math]\displaystyle{ y_{t} }[/math]经过宏观动力学[math]\displaystyle{ f: R_q → R_q }[/math]得到的预测[math]\displaystyle{ t+1 }[/math]时刻宏观变量值，[math]\displaystyle{ \hat{x}_{t+1} }[/math]表示[math]\displaystyle{ \hat{y}_{t+1} }[/math]经过反粗粒化函数[math]\displaystyle{ \phi^{\dagger} }[/math]得到的预测的[math]\displaystyle{ t+1 }[/math]时刻微观变量值，[math]\displaystyle{ g: R_q → R_q }[/math]表示反宏观动力学函数，可以根据[math]\displaystyle{ t+1 }[/math]时刻宏观变量值[math]\displaystyle{ y_{t+1} }[/math]推出预测的[math]\displaystyle{ t }[/math]时刻的宏观变量值[math]\displaystyle{ \hat{y}_{t} }[/math]，[math]\displaystyle{ λ }[/math]作为拉格朗日乘子，在实验框架内被认为是一个可调的超参数。

引理1——双射映射不影响互信息：

对于任意给定的连续随机变量X和Z，如果存在一个双射(一对一)映射f与另一个随机变量Y，使得对于任意[math]\displaystyle{ x\in Dom (X) }[/math]存在一个[math]\displaystyle{ y=f (x)\in Dom (Y) }[/math]，反之亦然，其中[math]\displaystyle{ Dom (X) }[/math]表示变量X的域，则X与Z之间的互信息等于Y与Z之间的互信息，即：

[math]\displaystyle{ I (X;Z)=I (Y;Z) }[/math]

引理2——连续自变量不影响互信息：

若[math]\displaystyle{ X\in Dom (X) }[/math]与[math]\displaystyle{ Y\in Dom (Y) }[/math]构成一条马尔可夫链[math]\displaystyle{ X\rightarrow Y }[/math]，且[math]\displaystyle{ Z\in Dom (Z) }[/math]是一个与X、Y均无关的随机变量，则：

[math]\displaystyle{ I (X;Y)=I (X;Y\bigoplus Z) }[/math]

引理3——条件熵的变分上界：

给定一个条件熵[math]\displaystyle{ H(\boldsymbol{y}|\boldsymbol{x}) }[/math]，其中[math]\displaystyle{ \boldsymbol{x}\in \mathcal{R}^s }[/math]，[math]\displaystyle{ \boldsymbol{y}\in \mathcal{R}^q }[/math]，则该条件熵存在一个变分上界：

[math]\displaystyle{ H(Y|X)\le -\iint p(\boldsymbol{y}, \boldsymbol{x})\ln g(\boldsymbol{y}|\boldsymbol{x}) \mathrm{d}\boldsymbol{y} \mathrm{d}\boldsymbol{x} }[/math]

其中[math]\displaystyle{ H(Y|X) \in R^q × R^s }[/math]是任意分布。

编码器的通用逼近定理

首先，我们扩展基本编码器的定义，引入一个新的运算[math]\displaystyle{ \eta_{p,s}: \mathcal{R}^p\rightarrow \mathcal{R}^s }[/math]，它表示原始变量的自复制。

[math]\displaystyle{ \eta_{p,s}(\boldsymbol{x})=\boldsymbol{x}\bigoplus \boldsymbol{x}_{s-p} }[/math]

向量[math]\displaystyle{ \boldsymbol{x}_{s-p} }[/math]是[math]\displaystyle{ s-p }[/math]维，其中每个维都是x中特定维的重复。例如，若[math]\displaystyle{ \boldsymbol{x}=(0.1,0.2,0.3) }[/math]，则[math]\displaystyle{ \eta_{2,5}(\boldsymbol{x})=(0.1,0.2,0.3,0.1,0.2) }[/math]。

学者们曾提出一般前馈神经网络和可逆神经网络的通用逼近定理，将其作为桥梁，可以证明任何前馈神经网络都可以用一系列双射映射(ψ)、投影(χ)和向量扩展(η)过程来模拟。对向量展开进行扩展后的基本编码器可表示为：

[math]\displaystyle{ \phi= Proj_q \circ\psi_{s} \circ \eta_{p,s}\circ \psi_{p} }[/math]

式中，函数[math]\displaystyle{ \psi_s: \mathcal{R}^s\rightarrow \mathcal{R}^s }[/math]和[math]\displaystyle{ \psi_p: \mathcal{R}^p\rightarrow \mathcal{R}^p }[/math]表示两个可逆映射。保留的最终维数q可能大于初始维数p。φ为降维算子。

根据通用逼近定理，对于定义在[math]\displaystyle{ K\times \mathcal{R}^p }[/math]上的任意函数[math]\displaystyle{ f }[/math]，其中[math]\displaystyle{ K\in \mathcal{R}^p }[/math]是紧集，且[math]\displaystyle{ p\gt q\in \mathcal{Z^+} }[/math]，则存在整数[math]\displaystyle{ s }[/math]和[math]\displaystyle{ W\in\mathcal{R}^{s\times p}, W'\in\mathcal{R}^{q\times s}, b\in\mathcal{R}^{s} }[/math]，使得:

[math]\displaystyle{ W'\cdot \sigma(W+b)\simeq f }[/math]

式中，[math]\displaystyle{ \sigma(\boldsymbol{x})=1/(1+\exp(-\boldsymbol{x})) }[/math]是向量上的sigmoid函数。

根据引理4，[math]\displaystyle{ +b }[/math]和[math]\displaystyle{ \sigma(\cdot) }[/math]都是可逆算子，因此，存在可逆神经网络[math]\displaystyle{ \psi_{q},\psi_{s}',\psi_{s},\psi_{p} }[/math]和两个整数[math]\displaystyle{ s_1,s_2 }[/math]（矩阵[math]\displaystyle{ W' }[/math]和[math]\displaystyle{ W }[/math]的秩），使得：

[math]\displaystyle{ (\psi_{q}\circ\eta_{s_2,q}\circ\chi_{s,s_2}\circ\psi_{s}')\circ(\psi_{s}\circ\eta_{s_1,s}\circ\chi_{p,s_1}\circ\psi_{p})\simeq W'\cdot\sigma(W\cdot+b) }[/math]

式中，[math]\displaystyle{ \psi_{s}\circ\eta_{s_1,s}\circ\chi_{p,s_1}\circ\psi_{p} }[/math]近似(模拟)函数[math]\displaystyle{ \sigma(W\cdot+b) }[/math]，[math]\displaystyle{ \psi_{q}\circ\eta_{s_2,q}\circ\chi_{s,s_2}\circ\psi_{s}' }[/math]近似(模拟)函数[math]\displaystyle{ W'\cdot }[/math]。

因此，如果令[math]\displaystyle{ \phi_{p,s,q}=(\psi_{q}\circ\eta_{s_2,q}\circ\chi_{s,s_2}\circ\psi_{s}')\circ(\psi_{s}\circ\eta_{s_1,s}\circ\chi_{p,s_1}\circ\psi_{p}) }[/math]，那么[math]\displaystyle{ \phi_{p,s,q}\simeq f }[/math]。

在实际应用中，虽然基本编码器和扩展版本不包括展开运算符，但我们总是在输入向量为编码器输入之前展开它。因此，有理由相信此定理仍然适用于堆叠编码器。

综上可知，编码器通用逼近定理：

对于任何连续函数[math]\displaystyle{ f }[/math]，定义在[math]\displaystyle{ K\times \mathcal{R}^p }[/math]，[math]\displaystyle{ K\in \mathcal{R}^p }[/math]是一个紧集，[math]\displaystyle{ p\gt q\in \mathcal{Z^+} }[/math]，存在整数[math]\displaystyle{ s }[/math]和扩展堆编码器[math]\displaystyle{ \phi_{p,s,q}: \mathcal{R}^p\rightarrow \mathcal{R}^q }[/math]（有[math]\displaystyle{ s }[/math]隐藏层）和扩展操作[math]\displaystyle{ \eta_{p,s} }[/math]，使得：

[math]\displaystyle{ \phi_{p,s,q}\simeq f }[/math]

此后，扩展堆叠编码器具有通用逼近性质，这意味着它可以近似(模拟)任何定义在[math]\displaystyle{ \mathcal{R}^p\times \mathcal{R}^q }[/math]粗粒化函数。

引理4——编码器的信息瓶颈：

对于任意向量[math]\displaystyle{ X\in \mathcal{R}^p }[/math]和矩阵[math]\displaystyle{ W\in \mathcal{R}^{s\times p} }[/math]，其中[math]\displaystyle{ s,p\in \mathcal{N} }[/math]，存在一个整数[math]\displaystyle{ s_1\leq \min(s,p) }[/math]和两个编码器的基本单位：[math]\displaystyle{ \psi_{s}\circ\eta_{s_1,s} }[/math]和[math]\displaystyle{ \chi_{p,s_1}\circ \psi_{p} }[/math]使得：

[math]\displaystyle{ W\cdot X\simeq(\psi_{s}\circ\eta_{s_1,s})\circ(\chi_{p,s_1}\circ \psi_{p})(X) }[/math]

式中，[math]\displaystyle{ \simeq }[/math]表示近似或模拟。

机器学习算法

样本重加权

输入数据的分布会对结果产生一定的影响，为了减少此影响，学者需要对样本进行重加权，也就是通过为数据中的每个单元分配适当的权重来改变输入数据的分布，解决选择偏差问题。在因果机器学习、因果推断领域、因果特征学习（Causal feature learning）和稳定学习（Stable learning）中发挥着重要作用。可以使结果减少偏见的影响，提高模型在未知环境下的泛化能力。

为了使用逆概率加权技术，我们需要估计样本的概率分布。KDE(Kernel Density Estimation)是一种常用的估计方法，它可以有效地消除离群值对整体概率分布估计的影响。

对于一个样本[math]\displaystyle{ (x_1, x_2, ..., x_n) }[/math]，我们有如下的核密度估计：

[math]\displaystyle{ \hat{f}_h(x)=\frac{1}{nh}\sum_{i=1}^n K(\frac{x-x_i}{h}) }[/math]

式中，n表示样本数，超参数h表示带宽（根据数据的粗略范围确定），文章中，作者一般设为0.05。K是核函数，具体来说是标准正态密度函数。

得到估计函数后，对每个样本点单独求值，得到每个样本点对应的概率值[math]\displaystyle{ p(y_t) }[/math]。将目标分布(均匀分布)的概率[math]\displaystyle{ \tilde{p}(y_t) }[/math]除以我们估计的概率值，得到每个样本点对应的逆概率权值[math]\displaystyle{ w(\boldsymbol{x}_t) }[/math]。为了覆盖所有的样本点，需要用参数L来限制均匀分布的范围，以保证边长为2L的正方形可以覆盖所有维度的所有样本点。

此时获得的权重与样本量有关，因此需要对它们进行归一化。然而，这将导致与样本对应的权重值非常小。因此，有必要将它们乘以样本数量以将权重值放大回其正常范围。然后，将该权值乘以每个样本的训练损失，增强稀疏区域的训练，达到逆概率加权的目的。由于每次迭代时编码器参数都会发生变化，导致[math]\displaystyle{ y_t }[/math]的分布发生变化，因此作者将每1000次迭代使用KDE重新估计整个样本的概率分布。

反向动力学

正向动力学训练是最小化预测误差[math]\displaystyle{ \left\|\phi_q^{\dagger}(Y(t+1))-X_{t+1}\right\| }[/math]，保证动力学预测未来的准确性，但是EI作为一种特殊的互信息，不仅与确定性有关，还与简并性有关，我们需要在提高动力学学习器的确定性的同时，提高它的非简并性。因此，学者在NIS的框架基础之上，加入了反向动力学，用以反向预测，即输入[math]\displaystyle{ y_{t+1} }[/math]，通过动力学学习器[math]\displaystyle{ g }[/math]之后，得到宏观量的反向预测值[math]\displaystyle{ \hat{y}_{t} }[/math]，使[math]\displaystyle{ y_{t+1} }[/math]和[math]\displaystyle{ \hat{y}_{t} }[/math]之间的误差值最小化。通过训练反向动力学学习器[math]\displaystyle{ g }[/math]，我们可以影响编码器，进而影响隐空间中的数据分布，从而使得动力学学习器[math]\displaystyle{ f }[/math]可以学到一个简并性低的动力学。

分阶段训练

在实践中，为了减轻机器学习的压力，可以先让模型学到比较靠谱的编码策略，然后通过反向动力学将编码策略调整至最优解。NIS+的训练过程分为两个阶段。

第一阶段：只训练前向神经网络，最小化预测误差[math]\displaystyle{ w(\boldsymbol{x}_t)\parallel\boldsymbol{y}_t-g(\boldsymbol{y}_{t+1})\parallel }[/math]。

第二阶段，只训练反向神经网络，最小化预测误差[math]\displaystyle{ \parallel\hat{x}_{t+1}-\boldsymbol{x}_{t+1}\parallel }[/math]，本质上是用于训练[math]\displaystyle{ \phi_{ω} }[/math]。

在整个模型的训练过程中，总的目标（式3）是最小化第一阶段的误差和第二阶段误差的拉格朗日函数，在不同的实验中，[math]\displaystyle{ λ }[/math]取不同的值。

面对大规模复杂系统的拓展

在实际应用中，如果系统不是小规模系统，而是类似元胞自动机的大规模的复杂系统，我们需要对此框架进行拓展，将编码器（解码器）进行组合，从而减轻模型训练的压力和难度。

首先，在处理高维复杂系统时，一次丢弃多个维度会给训练神经网络带来很大的挑战。我们可以将一系列基本编码器堆叠（串联）在一起并逐渐丢弃维度，降低训练难度。如右图a所示。

此外，一些大规模复杂系统具有先验知识，我们可以根据先验知识对微观维度进行分组，分组之后，对每一组都用编码器进行编码，相当于把编码器进行了并联。并行编码器之间共享参数，故神经网络依然可以高效、准确地获得粗粒化规则。最后，将从所有编码器获得的宏观变量连接成一个向量，以导出总体的宏观变量。这种并行结构也可以与卷积神经网络等架构结合起来。如右图b所示。

为了提高搜索最优尺度的效率，我们可以利用堆叠编码器获得隐藏空间的多个尺度，同时训练多个不同尺度的动力学学习器（相当于搜索不同q的宏观动力学），从而避免重新训练编码器，提高模型效率。如右图c所示。

将NIS+的编码器替换为堆叠编码器与并行编码器的任意组合时，式3的最优化函数依然适用（引理5、引理6）。

引理5——堆叠编码器不影响互信息：

若[math]\displaystyle{ X\in Dom(X) }[/math]和[math]\displaystyle{ Y\in Dom(Y) }[/math]形成一个马尔可夫链[math]\displaystyle{ X\rightarrow Y }[/math]， [math]\displaystyle{ \phi_L }[/math]和[math]\displaystyle{ \phi^\dagger_L }[/math]分别表示L层堆叠的编码器和解码器，则：

[math]\displaystyle{ I(X;Y)=I(X;\Phi^\dagger_L(Y)) }[/math]

引理6——并行编码器不影响互信息：

若[math]\displaystyle{ X\in Dom(X) }[/math]和[math]\displaystyle{ Y\in Dom(Y) }[/math]形成一个马尔可夫链[math]\displaystyle{ X\rightarrow Y }[/math]， [math]\displaystyle{ \phi_T }[/math]和[math]\displaystyle{ \phi^\dagger_T }[/math]分别表示由T个普通编码器或解码器组成的并行编码器和解码器，则：

[math]\displaystyle{ I(X;Y)=I(X;\Phi^\dagger_T(Y)) }[/math]

NIS+ 框架有效性的验证（数值实验）

SIR模型实验

在这个实验中，SIR（易感、感染、恢复或死亡）模型是一个简单的宏观动态系统，而微观变量则是通过在宏观变量中引入噪声产生的。主要目的是评估我们的模型是否能够有效地去除噪声、发现有意义的宏观动态、识别CE，并展示超越训练数据集分布的泛化能力。

(a) SIR模型的相空间，以及具有相同感染率和恢复率或死亡率的四个示例轨迹。同时显示用于训练的完整数据集(蓝色区域)和部分数据集(虚线区域)，分别由63,000和42,000个均匀分布的数据点组成。(b)曲线描述了不同模型的维平均有效信息(J)随训练时间的变化。线表示平均值，而带宽表示五次重复实验的标准差。(c)比较了SIR动力学的矢量场、NIS+的习得宏观动力学和习得编码器的雅可比矩阵变换后的宏观动力学。每个箭头代表一个方向，以及在该坐标点的动力学导数的大小。(d)将整体与局部进行比较，以评估在部分数据集(缺失42,000个数据点)或完整数据集上训练的不同模型的多步预测误差。这些模型包括NIS+、NIS、前馈神经网络(NN)、具有逆概率加权和反向动力学学习技术的前馈神经网络(NN+)、变分自编码器(VAE)及其重加权和反向动力学版本(VAE+)。(e).微观动力学(J (fm))和宏观动力学(J (fm))的CE(∆J)和EI的变化绘制为观测噪声变化的标准差σ。所有这些指标都是跨维度平均的。按照Rosas对CE的定义和计算方法，黄线表示Rosas ' Ψ的变化。垂直线表示归一化MAE的阈值为0.3。当σ大于阈值时，则不符合式1的误差约束，结果不可靠。(f)将SIR动力学的矢量场、NIS的习得宏观动力学和NIS的编码器雅可比矩阵变换后的宏观动力学与(c)进行比较。

比率导致观察到CE的发生。然而，当Rosas的Ψ < 0时，我们不能做出明确的判断，因为Ψ只能为CE提供充分条件。这两个指标在σ = 10−2处达到峰值，这与我们模拟中使用的时间步长(dt = 0.01)的大小相对应，反映了微观状态的变化水平。另一方面，如果噪声过大，有限的观测数据将使NIS+难以从数据中准确识别正确的宏观动力学。因此，CE的∆J的程度降至零。虽然NIS+判定σ > 10时不存在CE，但这一结果并不可靠，因为σ = 10−2后的归一化预测误差已经超过了所选阈值0.3(垂直虚线和虚线)。

因此，这些实验表明，通过最大化EI和学习独立的因果机制，NIS+可以有效地忽略数据中的噪声，准确地学习真实宏观动态，并泛化到未见数据。此外，NIS+在量化CE方面表现出优越的性能。

Boids 模型实验

Boids模型是一个著名的模拟鸟类集体行为的多代理模型，在这个实验中测试了NIS+ 在具有内在和外在噪声的不同环境中识别突发集体行为和CE 量化的能力。为了提高训练粗粒度策略的可解释性，将尝试在学习到的宏观状态和微观状态之间给出明确的对应关系。在300×300画布上使用N = 16个boids进行模拟以生成训练数据。为了评估NIS+发现有意义的宏观状态的能力，我们将所有boids分为两组，并人为地修改boids模型，为每组引入不同的恒定转向力。这种修改确保了两组具有不同转弯角度的独立轨迹。每个物体在每个时间步长的微观状态包括它们的水平和垂直位置，以及它们的二维速度。所有生物的微观状态形成一个4N维实数向量，作为训练NIS+的输入。

NIS+ boid

上图为NIS+学习Boids模型的集体群集行为的实验结果。(a)和(e)给出了各种条件下星体轨迹的实际和预测数据。具体来说，他们给出了在两个分离组和随机偏转角度条件下的多步(50步)预测的比较结果。它们的本征噪声水平α分别为0.001和0.4。(b)展示了多步预测的平均绝对误差(MAE)的上升，因为半径r(表示(a)中物体的初始位置范围)超出了训练数据的限制。(c)用不同的q超参数(代表不同宏观状态的尺度)描述了NIS+的训练时代中维度平均因果出现(∆J)的变化趋势。(d)为显著性图，直观地描述了每个宏观维度与每个体的空间坐标之间的关联。我们用橙色点突出了每个宏观状态维度中最重要的对应微观状态，这些微观状态是使用应用于模型的积分梯度(IG)方法确定的。横轴表示16个物体在微观状态下的x和y坐标，纵轴表示8个宏观维度。淡蓝色的虚线区分了不同个体体的坐标，而钢蓝色的实线分隔了实体组。

(f)和(g)表示不同噪声水平下∆J和归一化MAE的变化，(f)表示外部噪声(观测，加入微观状态)(δmax)， (g)表示本征噪声(α，通过修改Boids模型的动力学规则加入)。在(f)和(g)中，水平线表示违反式1中误差约束的阈值。当归一化MAE大于阈值0.3时，违反约束，结果不可靠。

NIS+可以通过最大化EI来学习最优宏观动态和粗粒度策略。这种最大化增强了它对超出训练数据范围的情况的泛化能力。学习到的宏观状态有效地识别了平均群体行为，并且可以使用IG方法将其归因于个体位置。此外，CE的程度随外在噪声的增加而增加，而随内在噪声的减少而减少。这一观察结果表明，通过粗粒化可以消除外在噪声，而内在噪声则不能。

生命游戏模型实验

本文以Conway’s Game of Life作为实验对象，其中每个细胞有两种状态作为二维状态输入:alive(1)或dead(0)，每个细胞受到八个相邻细胞的影响。生命游戏的进化只受输入状态及其更新规则的影响，其中生命游戏有四条进化规则，分别对应细胞繁殖和死亡等。Game of Life的更新规则如下表所示:Game of Life的训练样本生成过程如下:首先初始化状态xt。当考虑两个步骤的时间粗粒度时，根据更新规则生成状态xt+1、xt+2和xt+3的后续三个步骤，并将其输入到机器学习模型中。两个输入状态分别为xt和xt+1，微动力学输出分别为xt+1和xt+2。由于使用了时空粗粒度，宏观动力学将输出一个宏观状态，并将其解码为微观状态xt+2和xt+3。此过程重复多次(50,000个样本)，并生成图10d中用于训练的数据。而在其他实验中，我们生成了50万个样本。

然后，我们测试了在滑翔机模式上捕获动态模式的能力，其中模型是基于两个滑翔机模式进行训练的。该模型具有良好的预测效果，结果如图11所示。

此外，关于其他模型参数的更详细信息请参见表2。

大脑fMRI 时间序列数据模型实验

本实验在830名受试者的真实大脑fMRI时间序列数据上测试了NIS+模型。fMRI扫描数据是在受试者观看同一段电影时收集的。因此，在相似的自然刺激下，被试的相似经历是预期的，这对应于具有不同初始条件的相同动力学的时间序列。通过Schaefer atlas方法对原始数据进行预处理，将每个受试者的时间序列维数从大约140,000(不同受试者的维数不同)降至100，使NIS+能够运行并获得更清晰的结果。然后，选取前800个时间序列数据进行训练，剩下的30个时间序列进行测试。将结果与另一个fMRI数据集进行了比较，其中50名受试者处于静息状态。

上图为学习结果、因果涌现程度以及NIS+和NIS对大脑fMRI数据的归因分析。(a)在测试数据集上不同尺度(q)下，多步预测的平均误差随预测步长增加而增加。(b)比较不同模型和不同数据集(包括观影fMRI(视觉fMRI)和静息fMRI)的CE测量值(平均维度，∆J)。柱状图为10次重复实验的平均结果，误差柱状图为标准差。(c)给出了各子网络在Schaefer Atlas下的平均属性，采用积分梯度(IG)分析方法在编码器上以q = 1的尺度计算。误差条表示标准误差。(d)显示观影(视觉)fMRI数据的归因图。

为了证明NIS+对微观状态的预测能力，图(a)展示了在不同的超参数q下，归一化后的MAE随着微动力学对测试数据的预测步骤的变化。显然，当q = 27和q = 1时，NIS+的预测效果更好。具体来说，随着预测步骤的增加，与q = 1的曲线相比，q = 27的曲线显示出较慢的增长率。这表明选择超参数q为27可能比选择1更合适。然而，图(b)显示了不同的结果。当比较不同超参数q(绿色条)的CE(∆J)程度时，当q = 1时观察到最高的∆J。反之，当q = 27时，∆J值为负。这表明，当q = 27时，预测结果的改善可能归因于过拟合。因此，就∆J而言，q = 1优于q的其他值。尽管观察到q = 1时∆J的标准差较大，但NIS框架(红色条形图)支持了这一发现。此外，我们还将CE的结果与静息数据进行了比较，观察到NIS(深蓝色条)和NIS+(黄色条)在q = 7处达到峰值，这正好是Schaefer atalas中子系统的数量。因此，我们可以得出结论，当受试者看电影时，不同大脑区域的活动可以在每个时间步上用一个实数来表示。

为了研究NIS+粗粒度如何将输入数据转化为单维宏观状态，我们还利用IG方法来识别微观状态的最重要维度。结果如图 (c)和(d)所示。我们观察到视觉(VIS)子网络表现出最高的归因(图(c))。这些视觉子网络代表了受试者在观看电影片段时使用的功能系统。此外，我们可以在大脑地图上更详细地看到活动区域(图(d))，其中较深的颜色表示更大程度上归因于单一宏观状态。因此，NIS+识别出的具有相似最深颜色的区域，对应于大脑的深度视觉处理区域，可能代表了大脑在积极观看电影时的“协同核心”。这些区域的数字神经元可能协作并共同起作用。然而，通过将微观状态和宏观状态之间的相互信息分解为协同信息、冗余信息和唯一信息，还可以进一步证实和量化这一结论。 综上所述，NIS+展示了其学习和粗化来自大脑的复杂fMRI信号的能力，允许使用单个宏观状态模拟复杂的动态。这一发现揭示了电影观看过程中复杂的神经活动可以被一维宏观状态编码，主要集中在视觉(VIS)子网络内的区域。

总结与展望

NIS+框架的显著特征是，它专注于在保持有效性约束的同时最大化学习宏观动力学的有效信息(EI)。这使得学习到的涌现宏观动力学能够识别尽可能独立于输入数据分布的不变因果机制。这个特性不仅使NIS+能够识别不同环境中的数据CE，而且还增强了它对不同于训练数据的环境进行泛化的能力。因此，NIS+扩展了Hoel的CE理论，使其既适用于离散和连续动力系统，也适用于实际数据。

上述实验表明，通过最大化EI, NIS+在多步预测和模式捕获等任务中优于其他机器学习模型，即使在训练过程中没有遇到的环境中也是如此。因此，NIS+能够在潜在空间中获得更稳健的宏观动力学。NIS+比Rosas的Ψ指标更能合理地量化CE。NIS+可以有效地学习复杂系统的粗粒度化过程，并使用积分梯度(IG)方法可以将宏观状态和微观状态之间的关系可视化，从而可以识别微观状态中最重要的变量。在大脑实验中，发现存储在唯一宏观状态下的最关键信息与视觉子网络的微观变量表现出很强的相关性。这一发现突出了NIS+识别与有意义的生物学解释相一致的适当的粗粒度策略的能力。

NIS+在气候系统、集体行为、流体动力学、大脑活动和交通流等真实复杂系统的数据驱动建模方面具有各种应用潜力。通过学习更稳健的宏观动力学，这些系统的预测能力可以得到增强。除了动态学习之外，EI最大化方法可以在其他问题领域进一步扩展，如图像分类和自然语言理解。这可以促使诸如图像识别、对象检测、语言理解和机器翻译等任务的进步。

因果涌现与因果表征学习(causal representation learning, CRL)之间的关系。NIS+框架可以看作是CRL的一种形式，其中宏观动态作为因果机制，粗粒度策略作为表示。因此，可以应用CRL中使用的技术来查找数据中的CE。相反，可以在CRL中利用CE和粗粒度的概念来增强模型的可解释性。

NIS+的另一个有趣的优点是它对涌现理论的潜在贡献，它调和了关于涌现是一个客观概念还是一个依赖于观察者的认知概念的争论。通过设计一个最大化EI的机器，可以提取出客观的涌现特征和动态。

NIS+也存在一些局限性，首先，NIS+需要大量的训练数据来学习宏观动态和粗粒度策略，这在很多现实情况下可能是不可行的。如果培训不充分，可能导致CE的错误识别。因此，有必要结合其他数值方法，如Rosas的ΦID，以做出准确的判断。其次，神经网络的可解释性，特别是对于宏观动态学习者来说，仍然是一个难点。最后，当前的框架主要是为马尔可夫动力学设计的，而许多真实的复杂系统表现出长期记忆或涉及不可观察的变量。扩展NIS+框架以适应非马尔可夫动力学是未来研究的一个重要领域，解决这些限制并探索这些改进途径可以让NIS+应用于更广泛的复杂系统的数据分析。

其他

Rosas等学者通过将过量熵（excess entropy）分解成不重叠的部分来识别量化因果涌现，且提出了判定因果涌现发生的充分条件。

为了简化问题且不失通用性，以有两个输入变量（X₁、X₂）和一个输出变量（Y）的系统为例，目标变量和联合源变量的互信息[math]\displaystyle{ I(X_1,X_2; Y) }[/math]可以通过部分信息分解（Partial Information Decomposition，简称PID）分解成三种信息，分别是冗余信息（Redundant information，简称Red）、独特信息（Unique information，简称Un）、协同信息（Synergistic information，简称Syn），具体公式如下：

[math]\displaystyle{ I(X_1,X_2; Y) = Red(X_1,X_2; Y) + Un(X_1; Y |X_2) + Un(X_2; Y |X_1) + Syn(X_1,X_2; Y) }[/math]

不过，PID只适用于单个目标变量的情景，无法有效的应用于跨时间步长的系统。故，Rosas等学者提出了整合信息分解（Integrated Information Decomposition，简称ΦID）。当且仅当协同信息大于0的时候，系统有能力发生因果涌现。

[math]\displaystyle{ Syn(X_{t};X_{t+1}) \gt 0 }[/math]

在给定宏观状态[math]\displaystyle{ V }[/math]的情况下，[math]\displaystyle{ Syn(X_{t};X_{t+1}) ≥ Un(V_t;X_{t+1}|X_t) }[/math]。所以，如果宏观变量（V_t）所持有的关于微观变量独特信息大于0，则出现因果涌现。

[math]\displaystyle{ Syn(X_{t};X_{t+1}) ≥ Un(V_t;X_{t+1}|X_t) \gt 0 }[/math]

通过互信息的相关计算公式，可以得知：

[math]\displaystyle{ Un(V_t;X_{t+1}|X_t) ≥ I\left(V_t ; V_{t+1}\right)-\sum_j I\left(X_t^j ; V_{t+1}\right) + Red(V_t, V_{t+1};X_t) }[/math]

式中，[math]\displaystyle{ X_t^j }[/math]表示第 j 维t时刻的微观变量，[math]\displaystyle{ V_t ; V_{t+1} }[/math]代表两个连续时间的宏观状态变量。

由于[math]\displaystyle{ Red(V_t, V_{t+1};X_t) }[/math]为非负数，所以可以提出一个充分非必要条件[math]\displaystyle{ \Psi_{t, t+1}(V) }[/math]，用于测量两个时间步宏观变量的互信息减去每个t时刻微观变量和t+1时刻宏观变量的互信息。

当[math]\displaystyle{ \Psi_{t, t+1}(V) \gt 0 }[/math]，系统发生因果涌现。但当[math]\displaystyle{ \mathrm{\Psi}\lt 0 }[/math]，我们不能确定系统是否发生因果涌现。 [math]\displaystyle{ \Psi_{t, t+1}(V):=I\left(V_t ; V_{t+1}\right)-\sum_j I\left(X_t^j ; V_{t+1}\right) }[/math]

该方法避开讨论粗粒化策略，但是也存在很多缺点：

1）该方法只是基于互信息计算没有考虑因果，且得到的仅仅是发生因果涌现的近似的充分条件；

2）该方法需要预设宏观变量，且对宏观变量的不同选择会对结果造成显著影响；

3）当系统具有大量冗余信息或具有许多变量时，该方法的计算复杂度仍然很高。

Kaplanis等人基于机器学习的方法学习宏观态[math]\displaystyle{ V }[/math]以及最大化[math]\displaystyle{ \mathrm{\Psi} }[/math]：使用[math]\displaystyle{ f_{\theta} }[/math]神经网络来学习将微观输入[math]\displaystyle{ X_t }[/math]粗粒化成宏观输出[math]\displaystyle{ V_t }[/math]，同时使用神经网络[math]\displaystyle{ g_{\phi} }[/math]和[math]\displaystyle{ h_{\xi} }[/math]来分别学习[math]\displaystyle{ I(V_t;V_{t+1}) }[/math]和[math]\displaystyle{ \sum_i(I(V_{t+1};X_{t}^i)) }[/math]两者互信息的计算,最后通过最大化两者之间的差(即[math]\displaystyle{ \mathrm{\Psi} }[/math])来优化学习。 （暂定）

机器学习经过训练学习得到的很大程度上是变量之间的相关关系，相关关系的来源可以分成三种：

1）因果关系（Causation）：当且仅当T发生改变时，Y发生改变。比如在夏天的时候，冰淇淋的销量会上升。

2）混杂因素（Confounding）：X导致T和Y都上升，但是模型没有学到背后的X，只学出了T的上升导致Y上升，是一种虚假的相关关系。比如当人的年龄增加，很多人会吸烟和发福，而模型可能学出了是因为人吸烟，所以人发福。当给模型输入一组年龄一样，但是吸烟情况不一样的数据时，模型的有效性变会下降。

3）样本选择偏差（Sample Selection Bias）：当考虑了其他变量后，原本看似有相关性的变量可能实际上并无关联，是一种虚假的相关关系。在上面的狗识别的例子中，机器学习学到的可能是根据草地和狗鼻子这两个因素来判断是不是狗，但是有因果关系的是狗鼻子，草地只是由于样本的分布不均匀，这就导致了模型不能有效识别出在水里、沙漠里的狗。简单来说，草地和标签之间是强相关、弱因果，狗鼻子和标签之间是强相关、强因果。

根据以上公式，对于给定的维度q，在NIS中的不等式约束下的最大化问题方程1即可转化为无约束的损失函数最小化问题，这也是所定义的有约束目标函数的下界，即：

[math]\displaystyle{ \min_{f,g,\phi,\phi\dagger}\sum_{t=1}^{T-1}w(\boldsymbol{x}_t)||\boldsymbol{y}_t-g(\boldsymbol{y}_{t+1})||+\lambda||\hat{x}_{t+1}-\boldsymbol{x}_{t+1}|| }[/math]

[math]\displaystyle{ y_{t}=\phi(x_{t}) }[/math]

[math]\displaystyle{ y_{t+1}=\phi(x_{t+1}) }[/math]

[math]\displaystyle{ \hat{x}_{t+1}=\phi^{\dagger}\left(f_{q}(\phi(x_{t})\bigr)\right) }[/math]

[math]\displaystyle{ \hat{y}_{t+1}=f(y_t) }[/math]

式中，[math]\displaystyle{ x_{t} }[/math]、[math]\displaystyle{ x_{t+1} }[/math]表示可观测的微观数据，[math]\displaystyle{ y_{t} }[/math]、[math]\displaystyle{ y_{t+1} }[/math]表示经过粗粒化函数[math]\displaystyle{ \phi }[/math]得到的宏观数据，[math]\displaystyle{ \hat{y}_{t+1} }[/math]表示[math]\displaystyle{ y_{t} }[/math]经过宏观动力学[math]\displaystyle{ f: R_q → R_q }[/math]得到的预测[math]\displaystyle{ t+1 }[/math]时刻宏观变量值，[math]\displaystyle{ \hat{x}_{t+1} }[/math]表示[math]\displaystyle{ \hat{y}_{t+1} }[/math]经过反粗粒化函数[math]\displaystyle{ \phi^{\dagger} }[/math]得到的预测的[math]\displaystyle{ t+1 }[/math]时刻微观变量值，[math]\displaystyle{ g: R_q → R_q }[/math]表示反宏观动力学函数，可以根据[math]\displaystyle{ t+1 }[/math]时刻宏观变量值[math]\displaystyle{ y_{t+1} }[/math]推出预测的[math]\displaystyle{ t }[/math]时刻的宏观变量值[math]\displaystyle{ \hat{y}_{t} }[/math]，[math]\displaystyle{ λ }[/math]作为拉格朗日乘子，在实验框架内被认为是一个可调的超参数。[math]\displaystyle{ ω(y_t) }[/math]表示的逆概率权重具有以下定义：

[math]\displaystyle{ w(\boldsymbol{x}_t)=\frac{\tilde{p}(\boldsymbol{y}_t)}{p(\boldsymbol{y}_t)}=\frac{\tilde{p}(\phi(\boldsymbol{x}_t))}{p(\phi(\boldsymbol{x}_t))} }[/math]

其中，[math]\displaystyle{ \tilde{p} }[/math]为[math]\displaystyle{ do(y_t \sim U_q) }[/math]干预后宏观状态[math]\displaystyle{ y_t }[/math]的新分布，[math]\displaystyle{ p }[/math]为数据的自然分布。在实际运用中，使用核密度估计(Kernel Density Estimation，简称KDE)来近似[math]\displaystyle{ p(y_t) }[/math]。假设干预后的分布[math]\displaystyle{ \tilde{p}(y_t) }[/math]是均匀分布的。因此，权重[math]\displaystyle{ ω }[/math]由原始分布与修改后的分布之比决定

[math]\displaystyle{ \min_{f,g,\phi,\phi\dagger}\sum_{t=1}^{T-1}w(\boldsymbol{x}_t)||\boldsymbol{y}_t-g(\boldsymbol{y}_{t+1})||+\lambda||\hat{x}_{t+1}-\boldsymbol{x}_{t+1}|| }[/math]