1. 项目概述：当概率分布“流动”起来
最近在琢磨一个挺有意思的方向，就是把Wasserstein距离这套几何工具，跟神经动力学模型给揉到一块儿去。这听起来可能有点抽象，但它的核心想法其实挺直观的：我们能不能把大脑里神经元群体的活动，或者更一般地，把一堆相互作用的“智能体”的状态变化，看作是在一个“概率分布空间”里的某种“流动”或“运动”？
传统上，我们描述一个动态系统，比如一组神经元，会用它们的膜电位、发放率这些变量，然后写出一组微分方程。但如果我们换个视角，把关注点从单个神经元转移到整个群体的“状态分布”上呢？比如，在某个时刻，有多少比例的神经元处于高兴奋状态，多少处于静息状态，这个分布本身是如何随时间演化的？Wasserstein几何，或者说最优传输理论，恰好为我们度量这种分布之间的“距离”和描述其“最短路径”（测地线）提供了完美的数学语言。而“随机测地投影”则是在这个框架下，处理那些支撑集有限（比如只有有限几种状态）的随机系统时，一个非常精巧且必要的技术工具。它要解决的，是如何让一个离散的、随机的动态过程，尽可能地“行走”在连续分布空间的测地线上。这不仅仅是理论上的美感，对于理解神经编码的稳健性、设计高效的机器学习算法（特别是生成模型和强化学习），乃至分析社会网络中意见的传播，都可能打开一扇新的窗户。
2. 核心思路拆解：从静态距离到动态流形
要理解这个项目，我们得先掰开揉碎几个核心概念，看看它们是怎么串联起来的。这不像调个模型参数那么简单，它涉及到底层数学框架的转换。
2.1 Wasserstein距离：不只是度量，更是“搬运工”的蓝图
首先得说清楚Wasserstein距离（也叫推土机距离）。它度量的是两个概率分布之间的差异。想象你有两堆土，分布P和分布Q，形状不一样。Wasserstein距离问的是：要把P这堆土重新塑形成Q，最小的“总搬运功”是多少？这里的“功”通常定义为土的质量乘以搬运距离的p次方（常用p=2）。这个“最优搬运方案”本身，就蕴含了从P到Q的一种“自然”变换路径。
在神经动力学的语境下，这两个“土堆”就是系统在两个不同时刻的状态概率分布。Wasserstein距离告诉我们，系统状态演化的“成本”有多大。更重要的是，最优传输计划给出了一个如何将初始时刻的“概率质量”分配到目标时刻的详细方案，这暗示了一种可能的、最经济的演化方式。这比单纯用KL散度这类只关心密度值差异、不关心支撑集几何结构的度量要物理得多，也更有解释力——神经元状态的改变总是需要“能量”或“资源”的，Wasserstein距离天然地刻画了这种成本。
2.2 神经动力学与概率演化：从微观规则到宏观流形
经典神经动力学模型，比如Wilson-Cowan模型或发放率模型，描述的是大量神经元群体的平均活动。它们通常给出的是均值（比如平均发放率）随时间变化的方程。但一个更丰富的描述是刻画整个发放率分布的演化。这引向了Fokker-Planck方程或平均场理论，它们描述了概率密度函数如何随微观神经元动力学（包含确定性漂移和随机噪声）而演变。
这里的核心思想是：这个概率分布的演化轨迹，可以被视为在一个无限维的流形（所有可能概率分布的集合）上的“一条路径”。Wasserstein几何为这个流形赋予了黎曼结构。在这个流形上，两点（两个分布）之间的最短路径就是Wasserstein测地线。那么，一个自然的问题就是：由底层微观动力学（神经方程）驱动的概率流，与这个Wasserstein流形上的测地线，有什么关系？如果系统是“自由”的（没有外部约束），它的演化会沿着测地线吗？通常不会，因为微观动力学有自己的驱动项。但我们可以问：如何控制或设计微观动力学，使得宏观分布演化尽可能接近测地线？这可能是能量效率最优的演化方式，对应着神经信息处理或计算中的某种优化原则。
2.3 有限支撑系统的挑战与随机测地投影的登场
前面说的概率分布，通常是连续分布。但在很多实际场景中，系统的状态是离散且有限的。例如：
- 一个神经元的状态可以简化为“静息”、“阈下兴奋”、“发放”等有限几种。
- 一个社会网络中的个体，其观点可能只有“赞成”、“反对”、“中立”几种。
- 一个强化学习智能体在离散状态空间中行动。
这时，系统的状态分布是一个离散概率向量，其支撑集就是那有限几个状态点。问题来了：在离散支撑集上，经典的Wasserstein几何（基于连续最优传输）需要调整。两个离散分布之间的Wasserstein距离计算（本质上是一个线性规划问题）和它们之间的测地线定义，都变得不同。
更棘手的是动力学部分。假设我们有一个离散状态的马尔可夫链（这正是描述许多随机神经动力学或群体决策过程的模型），它的转移概率矩阵定义了分布如何一步演化。这个演化过程是随机的，并且是在离散状态空间上跳跃。而Wasserstein测地线描述的是连续分布空间中的确定性光滑路径。如何将离散的、随机的微观演化，“投影”或“近似”到连续的、确定性的测地线路径上？
这就是“随机测地投影”要解决的问题。它的目标是为离散随机过程（如马尔可夫链）在每一步，找到一个转移概率核，使得由该核驱动的分布演化，在Wasserstein意义下，最接近那条理想的、连接当前分布和目标分布（可能是测地线上的下一点）的连续测地线。这是一种在离散约束下，对连续最优路径的最佳随机逼近。
注意：这里的“投影”不是简单的正交投影，而是在随机核的集合中，寻找一个使得下一步的期望分布与测地线上的目标分布之间的Wasserstein距离最小的那个核。这通常归结为一个带有线性约束的凸优化问题。
3. 核心环节实现：构建随机投影算子的数学框架
理论说得再多，不如看看具体怎么构造。我们一步步来拆解如何为一个有限状态空间的系统，实现基于Wasserstein几何的随机测地投影。
3.1 问题形式化：定义舞台与角色
假设我们有一个离散状态空间 $\mathcal{X} = {x_1, x_2, ..., x_N}$，例如N个可能的神经活动模式或认知状态。系统在时间 $t$ 的状态分布是一个概率向量 $\mu_t = (\mu_t^1, ..., \mu_t^N)$，满足 $\mu_t^i \geq 0$ 且 $\sum_i \mu_t^i = 1$。
目标：我们有一条在连续分布空间中预先计算好的Wasserstein测地线 $\rho_s$, $s \in [0,1]$，连接初始分布 $\rho_0$ 和最终分布 $\rho_1$。对于我们的离散系统当前分布 $\mu_t$，我们假设它近似等于 $\rho_{s_t}$（即测地线上的某一点）。我们希望系统在下一个（离散）时间步 $t+1$ 的分布 $\mu_{t+1}$，能够尽可能地接近测地线上的下一个点 $\rho_{s_t + \delta}$，其中 $\delta$ 是一个小的时间步长。
手段：我们通过设计一个随机转移矩阵$\Pi_t$ 来实现这一步演化。$\Pi_t$ 是一个 $N \times N$ 的矩阵，其中元素 $\pi_t^{ij}$ 表示从状态 $x_i$ 转移到状态 $x_j$ 的概率。那么，演化方程为： $$\mu_{t+1} = \mu_t \cdot \Pi_t$$ 即，$\mu_{t+1}^j = \sum_{i=1}^N \mu_t^i \pi_t^{ij}$。
核心优化问题：寻找转移矩阵 $\Pi_t$，使得：
- $\Pi_t$ 是合法的随机矩阵（行和为1，元素非负）。
- 由它生成的 $\mu_{t+1}$ 与目标分布 $\nu_t := \rho_{s_t + \delta}$ 之间的Wasserstein距离$W_2(\mu_{t+1}, \nu_t)$ 最小。
这里我们通常使用 $p=2$ 的Wasserstein距离（$W_2$），因为它导出的几何结构最友好。对于离散分布，$W_2^2(\mu, \nu)$ 的计算是一个线性规划问题，或者等价地，是求解一个最优耦合（联合分布）$\gamma$，使得其边缘分布分别为 $\mu$ 和 $\nu$，并最小化 $\sum_{i,j} \gamma^{ij} \cdot d(x_i, x_j)^2$，其中 $d(x_i, x_j)$ 是状态空间 $\mathcal{X}$ 上定义的距离（例如，汉明距离、欧氏距离等，取决于状态的具体含义）。
3.2 算法步骤：从理论到可计算
实际操作中，我们无法直接对无限维的连续分布 $\rho_s$ 操作。通常，我们会用一组离散的样本来近似表示连续分布，或者我们的目标分布 $\nu_t$ 本身就是一个设计目标（例如，我们希望系统流向某个特定的分布）。以下是实现随机测地投影的关键步骤：
步骤一：定义状态空间与度量首先，必须明确定义离散状态集 $\mathcal{X}$ 和其上的距离函数 $d(\cdot, \cdot)$。这个距离至关重要，它决定了Wasserstein几何的“形状”。例如，如果状态是二进制神经元发放模式，距离可以是汉明距离（不同比特的个数）。如果状态是某种嵌入空间中的点（如词向量），则可以使用欧氏距离。
步骤二：计算或给定目标演化路径我们需要一条目标路径。有两种主要方式：
- 解析/数值测地线：如果我们知道初始和目标分布（$\mu_0$ 和 $\mu_T$），并且它们是连续分布或可以近似为连续分布，我们可以先计算连接它们的 $W_2$ 测地线。对于一维空间或某些特殊分布（如高斯分布），测地线有解析解。更一般的情况，需要数值求解McCam的偏微分方程或使用Sinkhorn迭代等近似算法。得到连续路径后，再在离散时间点上采样得到 ${\nu_t}$。
- 预设目标序列：在控制或规划问题中，${\nu_t}$ 可能直接是我们希望系统跟随的分布序列，它本身就被认为是“理想”的测地路径。
步骤三：构建单步随机投影优化问题在每一个时间步 $t$，我们面临如下凸优化问题： $$ \begin{aligned} \min_{\Pi_t, \gamma} & \quad \sum_{i,j,k} \gamma^{jk} \cdot d(x_j, x_k)^2 \ \text{s.t.} & \quad \mu_{t+1}^k = \sum_j \gamma^{jk} = \sum_i \mu_t^i \pi_t^{ij} \quad \forall k \quad \text{(目标分布匹配/近似)} \ & \quad \sum_k \gamma^{jk} = \mu_{t+1}^j \quad \forall j \ & \quad \sum_j \pi_t^{ij} = 1, \quad \pi_t^{ij} \geq 0 \quad \forall i,j \ & \quad \gamma^{jk} \geq 0 \quad \forall j,k \end{aligned} $$ 这里 $\gamma$ 是 $\mu_{t+1}$ 和 $\nu_t$ 之间的耦合矩阵。第一个约束将 $\mu_{t+1}$ 与转移矩阵 $\Pi_t$ 联系起来，第二个和第三个约束是耦合矩阵的边缘分布约束。目标函数是 $W_2^2(\mu_{t+1}, \nu_t)$。
实操心得：这个优化问题变量较多（$\Pi_t$ 和 $\gamma$），约束复杂。一个常用的简化技巧是不严格要求 $\mu_{t+1} = \nu_t$，而是将目标函数改为 $W_2^2(\mu_{t+1}, \nu_t)$ 加上一个对 $\Pi_t$ 的正则项（例如，鼓励稀疏性以减少控制成本）。这样，问题可以分解或迭代求解：先固定 $\Pi_t$ 更新耦合 $\gamma$（一个标准的OT问题），再固定 $\gamma$ 更新 $\Pi_t$（一个带线性约束的二次规划或线性规划问题）。使用Python的
POT(Python Optimal Transport)库和CVXPY库可以相对方便地搭建和求解此类问题。
步骤四：迭代执行与系统演化求解得到 $\Pi_t^*$ 后，系统就按照这个转移矩阵随机演化一步，得到新的分布 $\mu_{t+1}$。然后，将 $t$ 更新为 $t+1$，目标分布更新为 $\nu_{t+1}$，重复步骤三，直到到达最终时间。
3.3 一个简化示例：三状态系统的路径跟踪
假设状态空间是 $\mathcal{X} = {A, B, C}$，我们定义距离矩阵（例如，$d(A,B)=1, d(A,C)=2, d(B,C)=1$）。初始分布 $\mu_0 = (0.8, 0.2, 0.0)$，我们希望系统在3个时间步内演化到目标分布 $\mu_{target} = (0.1, 0.3, 0.6)$。
首先，我们（假设性地）计算一条离散化的“测地路径”作为目标序列，例如：
- $\nu_0 = (0.8, 0.2, 0.0)$
- $\nu_1 = (0.5, 0.3, 0.2)$
- $\nu_2 = (0.2, 0.3, 0.5)$
- $\nu_3 = (0.1, 0.3, 0.6)$
在 $t=0$ 时，$\mu_0 = \nu_0$。我们需要求解 $\Pi_0$，使得 $\mu_1 = \mu_0 \Pi_0$ 尽可能接近 $\nu_1$。通过求解上述优化问题，我们可能得到一个如下的转移矩阵： $$ \Pi_0 = \begin{bmatrix} 0.7 & 0.3 & 0.0 \ 0.2 & 0.8 & 0.0 \ 0.0 & 0.0 & 1.0 \quad \text{(状态C无质量，转移任意定义)} \end{bmatrix} $$ 计算得 $\mu_1 = (0.62, 0.38, 0.0)$。这与目标 $\nu_1=(0.5,0.3,0.2)$ 有差距，因为一步之内从A和B转移质量到C的成本（距离）很高，优化器会权衡“紧跟目标”和“转移成本”。然后我们用 $\mu_1$ 和 $\nu_2$ 求解 $\Pi_1$，如此继续。
这个例子展示了即使目标路径很平滑，有限状态空间的随机演化也可能无法完美跟踪，但随机测地投影给出了在转移概率约束下的最优妥协方案。
4. 应用场景与价值：不止于理论优雅
这套框架的价值在于它提供了一个原则性的方法，来设计或分析有限状态随机系统的演化，使其宏观行为符合某种几何最优性。下面看几个潜在的应用方向。
4.1 神经编码与信息传输的优化
大脑在处理信息时，可能面临着能量约束下的效率优化问题。将神经元群体活动的分布演化建模为Wasserstein空间中的运动，测地线可能对应着信息传输“功耗”最小的路径。随机测地投影则可以解释，在神经元放电存在随机性（噪声）且状态离散（全或无动作电位）的约束下，神经网络是如何通过调整突触连接强度（即影响转移概率 $\Pi_t$）来逼近这种高效路径的。这为理解神经环路如何实现稳健、高效的信息编码与计算提供了新的理论工具。
4.2 生成模型与概率流的学习
在机器学习中，扩散模型和基于流的生成模型的核心思想，就是学习一个将简单分布（如高斯噪声）变换到数据分布的确定性或随机性过程。这个过程可以看作是在分布空间中的一条路径。Wasserstein测地线提供了某种“理想”的、距离最短的路径。对于离散数据（如文本、分类图像），标准的扩散过程需要特别设计。随机测地投影框架可以指导我们如何为离散状态空间设计前向噪声过程和反向生成过程，使得整个变换路径在Wasserstein意义下更高效，可能带来更快的采样速度或更好的生成质量。
4.3 多智能体系统与群体决策的协调
想象一组机器人或自动驾驶车辆，它们需要从一种队形分布变换到另一种队形分布。每个智能体的运动是局部的、带有随机性的。我们可以将整个群体的空间分布视为一个概率分布。目标分布序列 ${\nu_t}$ 就是期望的队形变换路径（一条测地线）。每个智能体根据自己当前的状态（位置），依据由随机测地投影求解出的转移概率 $\Pi_t$ 来决定下一步的移动策略。这样，在个体层面只有局部随机规则，但宏观上整个群体能协调一致地、近似最优地完成队形变换，并且对个体故障有一定容错性。
4.4 计算神经科学中的模型验证与拟合
在计算神经科学中，我们常建立复杂的随机网络模型来模拟实验数据。如何判断一个模型产生的群体活动动力学是“合理”或“高效”的？我们可以从真实神经数据中，估计出群体活动分布随时间演化的序列，并计算其近似测地线。然后，将模型产生的分布演化序列与这条测地线进行比较，计算它们之间的平均Wasserstein距离。这个距离可以作为一个新的、基于几何的模型拟合优度指标。同时，随机测地投影的理论可以启发我们，如何调整模型参数（如连接权重、噪声水平）来使模型产生的路径更接近测地线，从而可能发现神经系统遵循的优化原则。
5. 常见挑战与实战调优心得
在实际尝试实现和应用这个框架时，会遇到不少坑。这里分享一些从理论到代码过程中积累的经验。
5.1 计算复杂度：维度灾难与近似求解
最直接的挑战是计算。状态空间大小 $N$ 一旦增大，转移矩阵 $\Pi_t$ 的变量数就是 $O(N^2)$，耦合矩阵 $\gamma$ 也是 $O(N^2)$。求解包含这么多变量的优化问题，即使使用凸优化求解器，在 $N$ 超过几百后也会变得非常缓慢。
应对策略：
- 利用稀疏性：在许多应用中，状态转移通常只在“相似”状态之间发生。可以预先定义一个邻域关系，强制 $\Pi_t$ 在非邻域状态间的转移概率为零，从而大幅减少变量。这需要合理定义状态空间的拓扑结构。
- 参数化转移矩阵：不直接优化整个 $\Pi_t$ 矩阵，而是用一个参数化的函数（如一个小型神经网络）来生成它，优化函数的参数。这尤其适用于需要学习一个策略的强化学习场景。
- 使用熵正则化OT：标准的 $W_2$ 计算是线性规划，可以用熵正则化的Sinkhorn算法来近似，它能将问题转化为一系列迭代的矩阵缩放运算，计算效率高，且易于GPU加速。虽然这是对精确Wasserstein距离的近似，但在实践中效果很好，且带来一定的平滑性。
- 分层或粗粒化：对于超大规模状态空间，可以考虑先对状态进行聚类或粗粒化，在粗粒化空间进行规划，再将策略细化到原始空间。
5.2 距离矩阵的定义：决定几何的本质
$d(x_i, x_j)$ 的选择是根本性的。它定义了什么是“状态之间的差异”。在神经动力学中，如果状态是神经元发放模式，汉明距离是自然的选择。但如果状态是更抽象的（如认知状态），距离的定义可能依赖于任务或先验知识。一个糟糕的距离定义会导致测地线没有意义，投影也就失去了价值。
心得：距离矩阵最好能反映状态变化的“实际成本”或“难度”。有时可以从数据中学习这个距离度量。例如，在观察了系统大量的自然演化后，可以用状态之间的平均首次到达时间或转移概率的负对数来定义一个有效的“距离”。
5.3 目标路径的生成：测地线怎么算？
对于任意的两个离散分布，计算连接它们的 $W_2$ 测地线本身就是一个研究课题。在离散支撑集上，测地线通常不是唯一的，而且可能不是连续分布的路径，而是一系列离散分布的插值（McCann的“位移插值”）。
实操建议：
- 对于一维有序状态：情况最简单，测地线可以通过累积分布函数的逆的线性插值得到。
- 对于网格状状态空间（如图像）：可以使用基于熵正则化OT的Sinkhorn插值，这已经有一些成熟的库（如
POT）支持。 - 作为控制目标：在很多场景下，目标路径 ${\nu_t}$ 并非来自严格的测地线计算，而是人为设定的、平滑的分布变化序列。此时，随机测地投影框架的价值在于提供了一种在随机性约束下“跟踪”任何给定路径的最优方法。
5.4 随机性与确定性之间的权衡
随机测地投影产生的是一个随机策略（$\Pi_t$）。有时我们可能希望系统行为更确定。可以在优化目标中加入一个正则项，例如最小化转移矩阵的熵 $-\sum_{ij} \pi_t^{ij} \log \pi_t^{ij}$，这会鼓励选择概率更集中（即更确定）的转移。正则项的强度是一个超参数，用于控制“紧跟目标”和“减少随机性”之间的权衡。
5.5 与经典方法的对比与联系
需要意识到，这个方法不是万能的。对于某些问题，经典的控制理论方法（如线性二次型调节器LQR在均值上的应用）或基于策略梯度的强化学习可能更直接。随机测地投影框架的优势在于其直观的几何解释和直接在分布层面进行优化的特性。它特别适用于：
- 任务目标本身就是关于分布形态的（如“使群体分布均匀覆盖某个区域”）。
- 系统的性能或成本天然用分布之间的距离来衡量。
- 我们需要对系统的宏观涌现行为有一个基于几何的理解。
我个人在尝试将这套理论用于模拟小规模神经元网络的可塑性规则设计时，最大的体会是：它迫使你从“群体统计”的层面去思考动力学，而不是纠缠于单个神经元的参数。调试过程中，最关键也最耗时的部分往往是状态表示和距离定义——这本质上是在定义你所关心的问题的“语义空间”。一旦这个空间定义合理了，后续的优化和演化往往会展现出令人惊喜的、符合直觉的协调行为。另一个深刻的教训是，初始目标路径不要设得太“激进”，即相邻 $\nu_t$ 之间的Wasserstein距离不要太大，否则单步投影的误差会累积，导致系统最终严重偏离目标路径。这好比让一个行动迟缓的人去跟踪一个飞奔的运动员，每一步的“最优追赶”都力不从心，最终只会越落越远。合理的做法是让目标路径的“速度”与系统随机演化的内在时间尺度相匹配。