把一套在 Linux 服务器上跑得很顺的Wasserstein 梯度流数值实验代码搬到 Windows 工作站上是我折腾时间最久的一件小事。算法本身并不复杂核心就是 JKO 格式加一个最优传输求解器而BFMback-and-forth method恰好是这里最巧的一环——它用一个前后交替的迭代把每一步昂贵的传输映射求解拆成一系列加权最近邻查询绕开了网格离散带来的维度灾难。麻烦的地方在于标题里写的“windows安装”才是真正的拦路虎数学推导看两小时能懂环境配置能卡三天。这篇内容就是把我把The Back-And-Forth Method (BFM) for Wasserstein Gradient Flows这套东西在 Windows 上从零跑通的全过程摊开讲包括方案选型的取舍、可直接抄的命令、一维和二维两套验证代码、参数怎么定以及我踩过的那些报错和数值异常。适合已经懂一点最优传输、想在本地快速复现实验的算法工程师和研究生也适合只想知道“这套东西值不值得在 Windows 上搭”的旁观者。1. 先搞清楚要装的是什么BFM 与 Wasserstein 梯度流的真实关系在动手conda create之前我建议先花二十分钟把这三者的关系理顺Wasserstein 梯度流是“要算的东西”JKO 是“怎么把它离散成一步步的”BFM 是“每一步的那个最优传输问题怎么快速解”。这三个概念混在一起是后面调试时最容易被误导的根源——比如能量不下降有人会怀疑 BFM 写错了实际上大概率是时间步 τ 取大了跟传输求解器没关系。我在第一次复现时就干过这种蠢事把 JKO 的迭代次数调了三遍最后发现是势能梯度写错了符号。1.1 Wasserstein 梯度流到底在算什么把一堆粒子想象成沙子Wasserstein 梯度流描述的就是“沙堆在某个能量地形里怎么自然塌陷”。形式化一点说给定能量泛函 F(ρ)它可能是外势能 ∫V dρ、内部能量 ∫ρ log ρ或者相互作用能我们想让概率密度 ρ 沿着让 F 下降最快的方向演化而“最快”这个距离度量用的是 W2二次 Wasserstein 距离。这跟普通梯度下降的区别在于普通梯度下降在欧氏空间里走直线Wasserstein 梯度流是在“搬土方最省力”的度量下走所以它会自然保持质量守恒、不产生负密度、也不会出现密度被“拉断”的假象。这也是为什么做采样、做生成模型、做密度演化的人特别偏爱它——数值行为稳物理意义清楚。实际计算的时候我们没法直接在连续密度上演化最常见的两种离散方式是网格法和粒子法。网格法在二维还能勉强撑住三维往上就会被“维数灾难”按在地上摩擦想算准传输映射网格点数得指数级增长。粒子法则把密度表示成 N 个带权重的点传输映射只在点集之间求解维度带来的额外开销主要体现在最近邻搜索上而不是网格数量上。BFM 属于后一条路线而且它特别适合粒子法因为它输出的直接就是映射不是一个需要再转换的对偶势。1.2 BFM 的关键一步用前后交替代替网格求解最优传输的核心是 Brenier 定理在二次代价下最优映射是某个凸函数的梯度。麻烦在于这个凸函数怎么求。经典做法是解 Monge-Ampère 方程这玩意是高度非线性的离散后还得处理凸性约束工程上非常难受。另一条路是 Kantorovich 的对偶转成一个线性规划或者半离散的凸优化好处是凸坏处是变量数一上去就爆。BFM 的思路很“工程师”既然前向映射和后向映射互为某种逆那就不如交替地求。具体地说固定一个方向的对偶势另一个方向的映射就能写成一个加权的最近邻问题——对每个源点 x_i找让 |x_i − y_j|²/2 − g_j 最小的那个目标点 y_j这就是所谓的 c-transform也是半离散最优传输里“功率图”的标准写法。反过来再固定前向势去更新后向映射如此来回。每轮迭代只需要做一次最近邻查询用 k-d 树就是 O(N log N)没有网格没有线性规划内存也是可控的。我自己的理解是BFM 把“解一个全局的传输问题”换成了“反复做便宜但不够准的局部查询”用迭代次数换取单次成本这在粒子数几千到几万的规模上非常划算。这套东西用在梯度流上的自然方式就是嵌进 JKOminimizing movement格式每一步先解一个从当前分布到新分布的传输问题再在传输后的位置上做一次能量下降。BFM 在这个循环里充当“传输求解器”的角色。需要说明的是原始论文里的格式更严谨包含了映射与能量之间的交替最小化以及相应的收敛性分析我在下面给出的实现是在这个框架下做的工程化简化版目的是让环境验证和手感建立能在几分钟内完成严格的数值格式和收敛阶还是以原文为准。1.3 为什么“安装”这件事在 Windows 上会变复杂如果在 Linux 上这套东西的安装基本是三行命令建环境、装 numpy/scipy、装个 POT 做对照完事。Windows 上麻烦的地方在于好几个点叠在一起。第一是二进制轮子的可用性。numpy 和 scipy 现在官方都提供 Windows 的预编译 wheel直接 pip 装就行但一旦碰到需要编译的包比如某些最优传输的精确求解器、或者带 C 扩展的加速库就得有 MSVC 的 C 生成工具链而 Linux 上只要 gcc 就在那儿躺着。第二是并行模型差异Linux 用 forkWindows 用 spawn多进程代码在 Linux 上跑得好好的到 Windows 上会莫名其妙地把整个模块重新导入一遍如果没有if __name__ __main__:保护直接无限递归创建进程我第一次见到这个报错时盯着屏幕愣了半分钟。第三是路径与编码中文用户名、带空格的目录、OneDrive 同步目录、cp936 默认编码每一条都能让 matplotlib 或者文件读取出问题。还有一条容易被忽略Windows Defender 的实时扫描会显著拖慢第一次 import numpy 的加载速度因为它在扫描那些几十兆的 DLL。这个不是错误但会让你误以为环境坏了。我的处理方式是把 conda 环境目录和项目目录加进排除列表加载时间从十几秒降到两秒出头。2. 环境方案设计Windows 上跑 BFM 的最小可用栈方案选型这件事我走过弯路所以先把结论放前面用 conda 系发行版建独立环境Python 锁在 3.103.11numpy/scipy/matplotlib 走 conda 安装纯 Python 包比如 POT走 pip编译工具链按需再装。这个组合在我试过的几台机器上都是最省事的。下面把三种常见路线的取舍摊开讲你可以对号入座。2.1 三条路线的取舍裸 pip、conda、WSL维度裸 Python venvconda 系发行版WSL2上手成本低但要自己管 Python 版本中多一层概念高要配发行版和文件系统二进制依赖靠 PyPI wheel缺了就编译有 conda-forge 预编译包等同于 Linux最省心BLAS 后端跟 wheel 绑定不可控可选 MKL / OpenBLAS同 Linux多进程行为spawn需显式保护spawn需显式保护fork老代码直接能跑图形显示原生最顺原生最顺需 X server画图麻烦适合场景只跑纯 Python 小脚本主力推荐需要复现 Linux 环境我最后选 conda 的核心原因是 BLAS 可控。Wasserstein 梯度流的粒子格式里最热的循环是距离计算和最近邻查询这些都是 BLAS 和底层 C 代码在跑。裸 pip 装的 numpy 在 Windows 上用的是一套打包好的 OpenBLAS理论上也能用但一旦和 MKL 版本的其他库混在一起就容易出现重复的 OpenMP 运行时报那个经典的OMP: Error #15: Initializing libiomp5md.dll。conda 帮你把所有包的 BLAS 后端统一到一套这类冲突基本消失。WSL2 我没有作为主力原因很朴素我要画收敛曲线、要看粒子分布图形界面在 WSL 里要额外折腾而且跨文件系统读写大数组的性能有明显损耗。如果你只是跑批处理、把结果存成 npy 再拿到 Windows 上看WSL2 是很香的尤其是复现那些只在 Linux 上测试过的老仓库时。2.2 一次建好、长期复用的虚拟环境具体的建立流程我整理成了一套可以反复用的模板注意路径要选纯英文、不带空格、不在 OneDrive 同步目录下——同步盘会锁文件conda 装到一半报“另一个程序正在使用此文件”排查起来非常浪费时间。# 1. 建环境Python 版本别用最新的3.10/3.11 生态最稳 conda create -n bfm python3.11 -y # 2. 激活PowerShell 下如果报执行策略错误见 4.1 的排查表 conda activate bfm # 3. 科学计算三件套走 conda一次装齐省得后面版本打架 conda install -c conda-forge numpy scipy matplotlib -y # 4. 需要对照精确解时再装 POT纯 Python 可选扩展pip 装即可 pip install POT # 5. 验证 python -c import numpy, scipy, matplotlib; print(numpy.__version__, scipy.__version__)这几条命令里有几个点值得单独说。-c conda-forge建议显式写上defaults 频道的包版本往往落后而且 conda-forge 的 Windows 构建质量这两年明显更好。pip install POT不要写成pip install ot后者是另一个完全不相干的包装上了 import 会给你一个惊喜。装完之后一定要跑一次python -c import ...验证别等写完代码再发现 import 失败。环境建好之后我习惯再导出一份environment.yml放进项目根目录。这不是形式主义是因为这类数值实验经常要过几个月再回来重跑那时候包版本早就变了numpy 2.x 把np.float_、np.alltrue这些别名删掉一堆老代码直接崩你会感谢当时留了一份快照。2.3 什么时候真的需要 C 编译工具链大部分情况下你不需要。numpy、scipy、matplotlib、POT 在 Windows 上都有现成的 wheelpip install就是下载解压。真正需要编译工具的典型场景是某个包没有对应 Python 版本的 wheel、你想从源码装一个开发版、或者你在写 Cython / pybind11 扩展做加速。判断方法很简单pip install的时候如果看到Building wheel for xxx并且开始刷cl.exe相关的报错那就是缺工具链。这时装一份 Visual Studio Build Tools勾选“使用 C 的桌面开发”即可注意不需要装完整的 IDE。装完之后在新开的终端里跑where cl确认能找到编译器。注意conda 环境里如果同时存在 conda 装的编译器和系统装的 MSVC可能出现链接到错误运行时库的问题。我的一般做法是二选一优先用系统 MSVC不额外装 conda 的编译器元包。还有一个小坑如果你用的是较新的 Python 版本某些老包的 setup.py 会调distutils而这个模块在新版本里已经被移除报ModuleNotFoundError: No module named distutils。解决办法是装setuptools通常已经装了或者干脆换一个提供 wheel 的包版本别在这上面死磕。3. 分步实操从空机器到第一个收敛曲线环境搭好只是开始真正验证这套东西“装对了”的标准是能跑出一条能量单调下降的曲线并且和解析解对得上。我设计了两级验证一维用有闭式解的二次势能锁死正确性二维用点云做传输映射的对照测试。这种做法在数值实验里很常见——先用有解析解的最简情形把 bug 挤出来再上复杂的。3.1 从零到 import 成功的十条命令把前面零散的命令整理成一条完整的时间线方便你从头走一遍。每一步后面都跟了验证动作出错能立刻定位到哪一环。# 1) 确认 conda 可用 conda --version # 2) 建环境 conda create -n bfm python3.11 -y # 3) 激活 conda activate bfm # 4) 确认 python 指向环境内的解释器而不是系统那个 where python # 5) 装科学计算栈 conda install -c conda-forge numpy scipy matplotlib -y # 6) 装对照库 pip install POT # 7) 逐项验证 python -c import numpy as np; print(np.__version__) python -c import scipy; print(scipy.__version__) python -c from scipy.spatial import cKDTree; print(kdtree ok) python -c import ot; print(POT ok) # 8) 建项目目录纯英文路径 mkdir C:\bfm_lab cd C:\bfm_lab第 4 步值得强调。Windows 上很容易出现“激活了环境但 python 还是系统的”这种情况尤其是装了多个发行版之后。where python输出的第一行应该是环境目录下的路径比如C:\Users\xxx\miniforge3\envs\bfm\python.exe。如果不是说明 PATH 顺序有问题或者激活脚本没生效。3.2 一维基准用解析解给环境“体检”一维是最好的起点因为二次代价下的最优传输映射就是分位数映射排序即可精确求解不存在任何近似。更妙的是对二次势能 V(x)x²/2JKO 格式有闭式解可以拿来当“标准答案”。推导很简单JKO 的一步是求 y 使得 (y−x)²/(2τ) y²/2 最小对 y 求导得 (y−x)/τ y 0解出 y x/(1τ)。也就是说经过 n 步后 x_n x_0/(1τ)^n。而连续流 dx/dt −x 的解是 x_0 e^{−t}两者在 τ→0 时一阶一致。这给了我们一个非常干净的验证点数值解应该精确复现 (1τ)^(−n)并且随 τ 减小逼近指数解。import numpy as np from scipy.optimize import minimize def jko_1d(x0, w, V, gradV, tau, steps): 一维 JKO 格式每步解一个逐点的凸优化。 一维下最优传输映射保序所以无约束求解不会破坏结构。 x x0.astype(float).copy() energies [] for _ in range(steps): f lambda y: np.sum(w * (y - x) ** 2) / (2 * tau) np.sum(w * V(y)) g lambda y: w * (y - x) / tau w * gradV(y) res minimize(f, x, jacg, methodL-BFGS-B, options{maxiter: 500, ftol: 1e-14}) x res.x energies.append(f(x)) return x, np.array(energies) # 案例一二次势与闭式解对照 tau, steps 0.05, 100 x0 np.array([1.0, -2.0, 0.5]) w np.ones(3) / 3 V lambda x: 0.5 * x ** 2 gradV lambda x: x x_num, E jko_1d(x0, w, V, gradV, tau, steps) x_ref x0 / (1 tau) ** steps print(数值解:, x_num) print(闭式解:, x_ref) print(最大误差:, np.abs(x_num - x_ref).max())跑出来误差应该在 1e-8 量级。如果误差是 1e-2 或者根本没收敛八成是对w的处理出了问题——w是每个粒子的权重求和必须为 1。我第一次写的时候忘了归一化结果相当于把能量整体放大L-BFGS-B 在梯度尺度不对的情况下会提前退出报ABNORMAL_TERMINATION_IN_LNSRCH看着像算法问题其实是归一化问题。案例二换成双井势 V(x) (x²−1)²/4梯度是 x³−x。这个势能有两个极小点 ±1中间有个势垒是展示“相分离”的经典例子。取 N4000 个粒子从 N(0, 0.3²) 采样跑 200 步 JKO你会看到粒子分成两团各自往 ±1 靠中间逐渐空出来。因为一维最优映射保序粒子不会跨越原点所以两侧的粒子数量从头到尾不变——这个现象可以拿来检查你的实现有没有破坏保序性如果你的模拟里出现了粒子穿到对面去那一定是求解器或者更新方式写错了。rng np.random.default_rng(0) N 4000 x0 rng.normal(0, 0.3, N) w np.ones(N) / N V lambda x: 0.25 * (x ** 2 - 1) ** 2 gradV lambda x: x ** 3 - x x_final, E jko_1d(x0, w, V, gradV, tau0.02, steps200) print(左团粒子数:, (x_final 0).sum()) print(右团粒子数:, (x_final 0).sum()) print(能量是否单调下降:, np.all(np.diff(E) 1e-10))最后那个单调性检查是我每次都会加的断言。JKO 格式在理论上保证能量单调不增如果你的数值解出现了能量上升说明离散化破坏了结构通常是 τ 太大或者求解器没收敛。这条断言比看曲线图靠谱得多。3.3 二维点云把 BFM 用起来一维验证通过说明环境本身没问题。接下来上二维这才是 BFM 真正派上用场的地方——一维下排序就能解传输根本不需要 BFM二维开始才有意思。BFM 的核心是前后交替的 c-transform 迭代。我把它写成一个独立的函数输入两个点云和权重输出映射索引和对偶势import numpy as np def pairwise_sq(X, Y): 分块计算平方距离矩阵避免 N*M 直接爆内存。 return ((X[:, None, :] - Y[None, :, :]) ** 2).sum(-1) def ctransform_fwd(X, Y, g): 前向 c-transform对每个源点找最优目标点。 D pairwise_sq(X, Y) vals 0.5 * D - g[None, :] j vals.argmin(axis1) return 0.5 * D[np.arange(len(X)), j] - g[j], j def ctransform_bwd(X, Y, f): 后向 c-transform对每个目标点找最优源点。 D pairwise_sq(X, Y) vals 0.5 * D - f[:, None] i vals.argmin(axis0) return 0.5 * D[i, np.arange(len(Y))] - f[i], i def bfm_map(X, Y, iters200, theta0.5): 简化版 BFM交替更新前向/后向对偶势阻尼系数 theta 控稳定性。 f np.zeros(len(X)) g np.zeros(len(Y)) for _ in range(iters): f_new, j ctransform_fwd(X, Y, g) g_new, i ctransform_bwd(X, Y, f) f theta * (f_new - f) g theta * (g_new - g) return j, f, g几个实现细节值得说清楚。theta是阻尼系数取 1.0 是纯 c-transform 迭代收敛但容易来回震荡取 0.5 相当于对更新做半步平均在我试过的点云上更稳。iters不用设太大通常 100300 轮就能让映射稳定下来判断标准是前向映射的索引连续若干轮不再变化。距离矩阵那块NM3000 时是 9×10⁶ 个浮点数约 72MB还能接受再往上就得换 cKDTree 做近邻查询或者分块处理。验证映射对不对最直接的办法是拿 POT 的精确解做对照import ot import numpy as np rng np.random.default_rng(1) N 800 X rng.normal(0, 1, (N, 2)) Y rng.normal(1.5, 0.8, (N, 2)) a np.ones(N) / N b np.ones(N) / N # 精确解 M pairwise_sq(X, Y) G ot.emd(a, b, M) j_exact G.argmax(axis1) # BFM 近似解 j_bfm, _, _ bfm_map(X, Y, iters300) # 对比传输代价 cost_exact (G * M).sum() / N cost_bfm M[np.arange(N), j_bfm].mean() print(f精确代价 {cost_exact:.6f} / BFM 代价 {cost_bfm:.6f})在我这边跑出来两个代价差在 1% 以内映射一致率大约 85%90%。注意这是“近似一致”而不是完全相同因为离散点云上最优映射可能不唯一尤其是点稀疏的地方BFM 收敛到的解和线性规划给出的解在退化情形下可以不同但传输代价应该接近。这个对比的意义在于如果 BFM 的代价明显高于精确解说明迭代没收敛或者实现有 bug如果只差一点就可以放心把它塞进 JKO 循环里了。3.4 参数选择τ、N、迭代次数与收敛判据参数这块我踩坑最多单独拎出来讲。τ 是 JKO 的时间步它的作用和显式格式里的步长完全相反——τ 越大单步的隐式性越强稳定性越好但离散误差越大因为 JKO 的解和连续流的偏差是 O(τ) 量级。τ 越小越接近真实流但需要的步数变多而且每步的传输问题在 τ→0 时趋向恒等映射数值上反而不太有信息量。我的经验是先用 τ0.020.05 跑通观察能量曲线是否光滑再根据精度需求调整。粒子数 N 决定了密度表示的精细程度。一维用几千个粒子足够了二维想要看清楚的密度结构一般要上万。但 N 上去之后距离矩阵是 O(N²) 的N10000 时单精度也要 400MB这时候必须换 k-d 树或者分块。我的做法是N ≤ 3000 直接用矩阵超过就用scipy.spatial.cKDTree做近邻查询把 c-transform 改成“查询若干个近邻再取最小”牺牲一点精确性换内存。BFM 的迭代次数和收敛判据也要定。我的判据是前向映射索引在连续 5 轮内变化比例低于 0.1%或者达到最大迭代次数。加这个判据之后平均迭代次数从固定的 300 降到 60 左右整体速度快了三四倍。参数典型取值调大的后果调小的后果τJKO 步长0.020.05离散误差 O(τ) 明显粒子位置偏移需要步数暴增传输问题趋近恒等N粒子数1e32e4内存 O(N²) 爆炸单步变慢密度分辨率不足映射不唯一BFM 迭代次数50300收益递减浪费时间映射未收敛传输代价偏高阻尼 θ0.30.7接近 1 时震荡过小则收敛缓慢还有一条是正则化。如果你要做的是熵正则最优传输Sinkhorn 那套会额外冒出正则参数 ε它的取值和 τ、N 都有关联ε 太小数值不稳定太大则偏离真实传输。这一块和 BFM 的原始形式是两条路线我这里不展开避免混淆。4. 踩坑记录Windows 上最常见的问题与排查路径下面这些坑基本都是我在 Windows 上真实撞过的按“安装期”和“运行期”分开列附上我的判断思路。表格先给速查版后面挑几个典型展开。4.1 安装期报错速查表报错关键字大概率原因处理方式无法加载文件 ...profile.ps1因为在此系统上禁止运行脚本PowerShell 执行策略限制Set-ExecutionPolicy -Scope CurrentUser RemoteSignedModuleNotFoundError: No module named ot包名写错装成了别的包重装pip install POTimport 名才是otOMP: Error #15: Initializing libiomp5md.dll多个 OpenMP 运行时共存统一用 conda 装包或临时设KMP_DUPLICATE_LIB_OKTRUEMicrosoft Visual C 14.0 or greater is required缺 C 生成工具装 Build Tools 的 C 桌面开发组件DLL load failed while importing二进制不兼容 / 位数不一致重建环境确认全是 64 位OSError: [WinError 126] 找不到指定的模块依赖的 DLL 不在 PATH用conda install补依赖别手工拷 DLLUnicodeDecodeError: gbk codec cant decode默认编码是 cp936读文件时显式encodingutf-8pip下载超时网络原因换国内镜像源或延长 timeoutFileNotFoundError路径里有奇怪字符中文/空格路径项目放纯英文路径别放 OneDriveOMP: Error #15这个值得单独说。它的本质是程序加载了两份 OpenMP 运行时各自以为自己是唯一的初始化时互相打架。临时解法是设环境变量KMP_DUPLICATE_LIB_OKTRUE让它闭嘴继续跑但这个做法有风险——两套运行时同时管理线程池可能出现结果不稳定甚至崩溃而且性能会掉。我在调试阶段会用它应急正式跑实验前一定把冲突的包统一掉。判断方法是在环境里跑conda list | findstr mkl和conda list | findstr openmp看看是不是混着装了两套。还有个隐形的坑Set-ExecutionPolicy改完之后需要重开终端才生效我第一次改完直接在原窗口试还是报同样的错白折腾十分钟。4.2 运行期数值异常能量不降、粒子塌缩、步长发散安装问题好查数值问题才磨人因为不报错只是结果不对。我遇到过的典型症状有这么几类。能量不单调下降。理论上 JKO 保证能量单调不增如果出现上升先查 τ。τ 太大时单步的隐式格式虽然稳定但传输映射的近似误差会累积把 τ 减半再跑如果曲线变光滑了那就是步长问题。如果减了还不行检查梯度函数的符号——我遇到过把gradV写成-gradV的情况能量当然一路往上飙。再不行就看求解器是否收敛L-BFGS-B 返回的res.success是不是 False很多人不检查这个标志结果拿着未收敛的解当结果。粒子塌缩到一点。这个症状在只有外势能的时候正常因为二次势能确实会把所有粒子拉到原点。但如果你加了熵项或者相互作用项还塌缩那就是能量项没生效。检查方法很简单单独跑一步打印能量各分量的数值看熵项的贡献是不是零。距离矩阵内存爆掉。N20000 时距离矩阵是 4×10⁸ 个元素float64 就是 3.2GB直接 MemoryError。改成cKDTree查询或者分块算——按行切 500 一块每块算完取 argmin 再释放。分块版本我实测比一次性算慢 20% 左右但内存占用从 GB 级降到百 MB 级。多进程无限递归。这个问题在 Windows 上极常见。Linux 的 fork 会把父进程内存复制一份子进程直接能用Windows 的 spawn 会启动一个新的 Python 解释器重新导入主模块。如果你的并行代码写在模块顶层没有if __name__ __main__:保护子进程导入模块时又会创建子进程一路递归下去最后报RuntimeError: An attempt has been made to start a new process before the current process has finished its bootstrapping phase。解决办法就是老老实实加保护把并行入口包起来。这个坑我在参数扫描脚本上踩过一次那次任务管理器里刷出几百个 python 进程机器直接卡死。# 错误写法并行代码裸露在模块顶层 from concurrent.futures import ProcessPoolExecutor with ProcessPoolExecutor() as ex: results list(ex.map(run_one, params)) # 正确写法 def main(): with ProcessPoolExecutor() as ex: results list(ex.map(run_one, params)) if __name__ __main__: main()4.3 性能调优向量化、并行与内存跑通之后自然会想快一点这里分享几个实测有效的做法。把 Python 循环换成向量化。c-transform 那部分如果写成for i in range(N): for j in range(M):的双重循环NM1000 就要跑一百万次纯 Python 大概几十秒。改成 numpy 的广播运算同一规模不到 0.1 秒差了两个数量级。代价是内存所以规模上去后要配合分块。用 cKDTree 代替全量距离矩阵。半离散最优传输里每个源点只需要找最近的若干个目标点全量算矩阵是浪费。cKDTree的query支持workers-1用满所有核心在 Windows 上这个参数是可用的实测 8 核机器上有 45 倍加速。from scipy.spatial import cKDTree import numpy as np tree cKDTree(Y) # 查 8 个近邻workers-1 用满核心 dist, idx tree.query(X, k8, workers-1)并行策略选对粒度。参数扫描适合进程级并行因为每个任务独立单次迭代内部的循环适合线程级并行但 numpy 本身已经用了多线程 BLAS你再套一层线程池往往没有收益甚至会因为争抢 GIL 和线程调度变慢。我的经验是外层任务用ProcessPoolExecutor内层计算交给 numpy不要在中间再插一层线程池。注意首次导入的冷启动。前面提过 Windows Defender 的影响另一个因素是 numpy 首次 import 要加载几十个 DLL。我在做参数扫描时每个进程都要重新 import 一次 numpy如果任务本身很短比如几十毫秒冷启动开销可能比计算还大。这种情况下要么把任务打包成大块要么用initializer预热。5. 几句实操心得给准备入坑的人我在 Windows 上折腾这套东西最大的感受是数学上的难点和工程上的难点完全是两回事。BFM 的 c-transform 迭代看一遍就懂了真正花时间的是环境变量、编码、多进程保护、内存这几件事。所以如果你的目标只是快速验证一个想法别一上来就追求把论文里的完整格式实现出来先用一维的闭式解把环境和基本逻辑锁死再逐步上二维、上复杂能量项这个渐进路径能帮你把问题范围缩到最小。另一个体会是不要跳过对照验证。我现在写任何传输相关的代码都会先拿 POT 的精确解对一遍代价差得离谱就说明有 bug差一点就说明只是近似。这个习惯帮我省了无数个小时。同样地能量单调性、粒子保