Massively Parallel Exact Inference for Hawkes Processes¶

会议: NeurIPS 2026
arXiv: 2604.01342
OpenReview: https://openreview.net/forum?id=6yoruSsZzc
代码: https://github.com/ahmrr/HawkesTorch
领域: 时间序列
关键词: 霍克斯过程, 并行前缀扫描, 精确极大似然估计, 点过程, GPU加速

一句话总结¶

将经典线性指数霍克斯过程的似然与梯度递推重构为增广状态下的稀疏转移矩阵乘法，结合 Blelloch 并行前缀扫描与分批状态传递，实现无近似、常数显存开销的百万至千万级事件极大似然并行精确推断。

研究背景与动机¶

多元霍克斯过程（Multivariate Hawkes Process, MHP）是一类经典的自激与互激时序点过程模型，广泛应用于高频量化金融、地震余震预测、流行病传播溯源、神经脉冲动力学以及社交媒体信息级联建模。在多元线性霍克斯模型中，各节点间的交叉互激矩阵天然对应有向信息流与 Granger 因果依存结构，具有极高且直观的可解释性。然而，极大似然估计（MLE）的计算复杂度严重制约了其向真实大规模场景的扩展：在包含 $N$ 个时序事件的序列中，朴素计算每个事件发生时的条件强度函数需要累加历史全部前序事件的激发衰减，全局时间复杂度高达 $\mathcal{O}(N^2)$。

经典线性指数衰减核由于具备马尔可夫无记忆性，强度函数允许展开为一维时序仿射递推，从而将时间复杂度压制在 $\mathcal{O}(N)$。但这种标准递推本质上是沿时间轴严格串行的，无法充分利用现代 GPU 上成千上万个计算核心的并行计算吞吐能力。为了缓解扩展瓶颈，以往研究通常引入低秩参数分解、随机优化、期望最大化（EM）或事件稀疏性截断等近似手段，但这不仅牺牲了精确对数似然的保真度，还往往削弱了模型本身体系清晰、物理含义明确的优势。

线性指数核的强度递推实质上具备一种特殊的稀疏仿射代数结构，只要经过恰当的代数增广，即可转化为满足结合律的线性矩阵乘积。核心 idea：通过状态增广将强度递推与参数梯度递推转化为具备结合律且单步仅需 $\mathcal{O}(M)$ 的稀疏转移矩阵乘积，直接借助 Blelloch 并行前缀扫描算法将事件维度推断并行化为 $\mathcal{O}(N/P)$，配合分批边界传递实现恒定显存的大规模精确求解。

方法详解¶

整体框架¶

算法的目标是在给定事件序列 $\mathcal{H}=\{(t_i, m_i)\}_{i=1}^N$（其中 $t_i$ 为递增时间戳，$m_i \in \{1,\dots,M\}$ 为节点索引）和 $K$ 个指数衰减核（衰减率 $\gamma_k$）的设定下，并行求解对数似然 $\mathcal{L}(\bm{\mu}, \mathcal{A}, \bm{\gamma})$ 及其关于基线强度 $\mu_p$、交互强度 $\alpha_{p,q,k}$ 以及衰减率 $\gamma_k$ 的精确梯度。

整个端到端流水线分为前向强度计算和反向解析梯度计算两个阶段。前向传播通过对增广状态进行并行前缀扫描，提取出每个事件发生时刻各核的状态向量 $\mathbf{R}_i^k$，计算条件强度并累加对数似然；反向传播避开 PyTorch Autograd 遍历前缀扫描计算图带来的巨大解释器开销，分别构造事件指示向量前缀扫描与时间加权状态前缀扫描，闭式解析求出全部参数的一阶导数。

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    In["输入事件序列 (t_i, m_i)<br/>N 个事件, M 个节点, K 个核"] --> FwdPre["状态增广与稀疏前缀扫描<br/>构建稀疏转移矩阵序列并执行 Blelloch 扫描"]
    FwdPre --> IntensityCalc["并行条件强度计算<br/>并行得到各事件时刻强度并累加似然"]
    IntensityCalc --> BwdPre["反向解析梯度并行前缀扫描<br/>双路前缀扫描计算参数导数所需的状态"]
    BwdPre --> BatchGrad["分批状态传递与恒定显存更新<br/>跨批传递边界状态, 累积全局梯度并更新参数"]

关键设计¶

1. 状态增广与稀疏前缀扫描：将串行仿射递推转化为高维线性矩阵链乘

标准线性指数核强度为 $\lambda_p(t) = \mu_p + \sum_{k=1}^K \gamma_k R_p^k(t)$，其中未加权激发状态满足仿射关系 $\mathbf{R}_i^k = e^{-\gamma_k \Delta t_i} \mathbf{R}_{i-1}^k + e^{-\gamma_k \Delta t_i} \bm{\alpha}_{:, m_{i-1}, k}$，其中 $\Delta t_i = t_i - t_{i-1}$。为消除常数偏置项带来的非线性，作者构造了 $(M+1)$ 维的增广状态向量 $\tilde{\mathbf{R}}_i^k = [\mathbf{R}_i^k; 1]$，并将单步转移算子表示为分块矩阵： $$ \mathbf{M}i^k = \begin{pmatrix} e^{-\gamma_k \Delta t_i} I_M & e^{-\gamma_k \Delta t_i} \bm{\alpha} $$ 由于所有转移矩阵均具备对角标量与单列偏置的特殊稀疏结构，两个转移矩阵相乘具有闭合保持性：}, k} \ \mathbf{0}^\top & 1 \end{pmatrix$(s_A I_M, \mathbf{v}_A) \times (s_B I_M, \mathbf{v}_B) = (s_A s_B I_M, s_A \mathbf{v}_B + \mathbf{v}_A)$。该矩阵乘法满足结合律，且计算与显存复杂度从朴素稠密矩阵的 $\mathcal{O}(M^3)$ 和 $\mathcal{O}(M^2)$ 严谨降至线性的 $\mathcal{O}(M)$。因此，前向状态序列直接化简为前缀矩阵连乘 $\tilde{\mathbf{R}}_i^k = (\prod_{j=1}^i \mathbf{M}_j^k) \tilde{\mathbf{R}}_0^k$，进而可无缝调用 Blelloch 并行前缀扫描算法，将事件维度计算耗时从串行 $\mathcal{O}(N)$ 压低至 $\mathcal{O}(N/P + \log N)$（当处理器核心数 $P$ 充足时逼近 $\mathcal{O}(N/P)$）。

2. 并行条件强度计算：数值稳定的事件间隔差分与似然评估

直接展开历史累加式若将公共因子 $e^{-\gamma_k t_i}$ 提出，会导致在长序列推断中 $e^{\gamma_k t_j}$ 迅速溢出（Overflow），同时 $e^{-\gamma_k t_i}$ 产生下溢（Underflow）。本方法直接以相邻事件时间差 $\Delta t_i$ 构建单步衰减因子 $e^{-\gamma_k \Delta t_i} \in (0, 1]$，并在并行前缀扫描中执行浮点结合律运算，在根源上彻底规避了长时间序列下的数值溢出与浮点下溢隐患。在前缀扫描完成后，各事件的瞬时条件强度 $\lambda_{m_i}(t_i)$ 即可在 GPU 线程级并行求和： $$ \lambda_{m_i}(t_i) = \mu_{m_i} + \sum_{k=1}^K \gamma_k R_{m_i}^k(t_i) $$ 对数似然的前半部分 $\sum_{i=1}^N \log \lambda_{m_i}(t_i)$ 从而能在微秒级完成全序列归约。

3. 反向解析梯度并行前缀扫描：避免动态计算图开销的双路状态推导

若直接利用 PyTorch Autograd 对 Blelloch 树形扫描图（包含约 $2N$ 个节点与 $\lceil \log_2 N \rceil$ 树层）执行自动求导，Python 解释器在逐层串行遍历计算图节点时会造成严重停顿，且无法发挥同层算子的 GPU 并行优势。作者通过解析微分推导，证明参数梯度所依赖的充分统计量同样遵循稀疏仿射递推。具体引入两个反向辅助状态：节点激发响应状态 $\mathbf{K}_i^k$ 与时间加权状态 $\mathbf{L}_i^k$，其增广转移矩阵形式与前向高度同构： $$ \mathbf{K}i^k = e^{-\gamma_k \Delta t_i} \mathbf{K}}^k + e^{-\gamma_k \Delta t_i} \mathbf{e{m}}, \quad \mathbf{Li^k = e^{-\gamma_k \Delta t_i} \mathbf{L}}^k + t_{i-1} e^{-\gamma_k \Delta t_i} \bm{\alpha{:, m $$ 其中 }, k$\mathbf{e}_m$ 为仅在第 $m$ 个分量为 1 的独热向量。通过执行两路与前向算子相同的并行前缀扫描，梯度算子 $\nabla \mu_p$、$\nabla \alpha_{p,q,k}$ 以及 $\nabla \gamma_k$ 均可在常数次张量广播与点乘下精确闭式求解，彻底省去了自动微分图的构建与遍历。

4. 分批状态传递与恒定显存更新：突破单卡显存墙的大规模分块流水线

当节点规模 $M \gtrsim 10^3$ 且事件数 $N \gtrsim 10^6$ 时，若直接在 GPU 显存中保存整条序列的前缀中间状态，将引发高达 $\mathcal{O}(KMN)$ 的显存占用爆炸。基于结合律性质，长序列可被划分为若干大小为 $N_b$ 的连续事件块 $\mathcal{I}_b$。由于第 $b$ 个分块内任意事件的状态仅依赖于前一块边界处的末尾状态 $\tilde{\mathbf{R}}_{(b-1)N_b}^k$，算法可以在批内局部执行前缀扫描，再以单步仿射变换融入前驱边界向量： $$ \tilde{\mathbf{R}}i^k = \left[\prod}^i \mathbf{Mj^k\right] \tilde{\mathbf{R}}_b $$ 每个 Batch 在前向完成后立即执行反向梯度累加并释放批内中间激活值，使得单卡显存消耗被严格锁定在与总事件数 }^k, \quad \forall i \in \mathcal{I$N$ 无关的 $\mathcal{O}(K M N_b)$ 恒定常数级别，同时该状态转移协议天然适配跨多卡甚至分布式集群流水线并行。

损失函数 / 训练策略¶

模型采用带截断 $\ell_1$ 稀疏惩罚的负对数似然损失函数进行优化，以鼓励节点间交叉交互的稀疏性，同时避免对固有的强自激效应过度收缩： $$ \mathrm{NLL}(\bm{\mu}, \mathcal{A}) = -\frac{1}{N}\left[\sum_{i=1}^N \log \lambda_{m_i}(t_i) - \sum_{m=1}^M \int_0^T \lambda_m(t) dt\right] + \lambda_1 \sum_{p,q,k} \alpha_{p,q,k} \mathbf{1}{\alpha_{p,q,k} < h} $$ 其中 $\lambda_1 \ge 0$ 为稀疏惩罚系数，$h > 0$ 为铰链截断阈值。优化器统一采用 Adam 进行参数全更新，学习率设为 0.05。

实验关键数据¶

主实验¶

在真实社交媒体网络 MemeTracker 信息级联数据集（涵盖 10 个具有代表性的事件子集，最大包含超 3,200 万个真实事件、约 2,500 个站点节点）上，对比现有利用事件稀疏性与低秩近似的基线方法 Lazy MHP 的单轮平均训练耗时。

子集名称 (Subset)	节点数 ($M$)	事件数 ($N$)	对数似然 (Log-Lik.)	本文耗时 (Ours)	Lazy MHP 耗时	加速比 (Gain)
Bail Out	835	64,138	3.612	0.06 s	0.47 s	7.8x
Miami Heat	1,173	133,451	1.680	0.21 s	1.60 s	7.6x
Amy Winehouse	1,561	226,247	2.043	0.34 s	3.30 s	9.7x
Arab Spring	1,377	400,199	2.796	0.46 s	4.60 s	10.0x
Greece	1,925	1,077,125	2.280	1.24 s	—	可行扩展
Libya	2,003	2,486,030	2.679	2.97 s	—	可行扩展
Crisis	2,266	5,085,759	2.495	6.94 s	—	可行扩展
Leader	2,344	10,604,129	2.674	14.91 s	—	可行扩展
North Korea...	2,451	16,408,957	2.747	23.73 s	—	可行扩展
Prince William...	2,494	32,687,780	2.827	48.33 s	—	可行扩展

在涵盖 22 个警区的芝加哥犯罪时空事件数据集（2001-2025 年超 840 万条精确记录）上，拟合多核 Hawkes 模型的耗时如下： - 全部犯罪类型联合序列（$N = 7,834,628$）：平均单 Epoch 耗时仅 1.012 s（Log-Likelihood: 1.012）； - 盗窃类（Theft, $N = 2,071,626$）：平均单 Epoch 耗时 0.296 s； - 攻击类（Assault, $N = 563,535$）：平均单 Epoch 耗时 0.082 s。

消融实验¶

在合成数据集上，针对不同似然及梯度实现架构进行效率与显存消融测试（固定最大运行时间 24 小时、显存上限 80 GiB）。同时在合成尺度自由网络（$N=1,000,000$）上检验精确 MLE 的参数恢复误差。

实现配置 / 架构	时间复杂度 (N 维)	峰值显存占用特征	最大可支持事件规模	运行瓶颈与终止原因
Naive 暴力实现	$\mathcal{O}(N^2)$	随序列平方急剧暴增	$N < 10^5$	严重显存溢出 (OOM) 且运算严重超限
Sequential (GPU 串行递推)	$\mathcal{O}(N)$	常数级别（低显存）	极低规模截断	单核计算效率低下，单步耗时极高超时
Sequential (CPU 64核串行)	$\mathcal{O}(N)$	常数级别（低显存）	约 $N \approx 5 \times 10^5$	递推本质串行，多核无法分摊，耗时超限
Prefix Scan + PyTorch Autograd	$\mathcal{O}(N/P + \log N)$	$\mathcal{O}(N)$ 计算图开销	$N \approx 10^6$	计算图庞大且求导解释器开销大，显存溢出
Full Model (Prefix Scan + Manual Grad + Batching)	$\mathcal{O}(N/P)$	常数显存 (固定 $N_b$)	$N > 6.7 \times 10^7$	线性稳定扩展，单轮仅需数十秒至数百秒

在参数恢复实验中，当每个交互参数的平均事件样本数充分（$N/M^2 \ge 16$）时，基线率 $\bm{\mu}$ 和交互矩阵 $\mathcal{A}$ 的相对 $\ell_2$ 误差均严谨收敛至 5% 以下；衰减率 $\gamma$ 在所有设置下相对误差均稳定低于 2%。

关键发现¶

运算加速随着事件规模增长更加显著：在短序列（6.4 万事件）下加速比约为 7.8 倍，而当事件数增加到 40 万时加速比提升至 10 倍，在大规模场景下由于充分填满 GPU 核心流水线，相对串行递推实现了跨量级的断层式提速。
显存控制机制表现稳健：启用分批（Batch Size 取 $2^{19} = 524,288$）后，处理序列长度从 $2^{10}$ 翻倍暴增至 $2^{26}$（约 6,700 万事件）过程中，峰值显存始终压制在 47.31 GiB，完全满足单张 A100（80GB）的运行要求。
自动微分成为不可忽视的效率暗礁：Blelloch 算法的多级归约图在 Autograd 中形成了极其庞大的 DAG，导致前向与反向遍历产生巨大的 CPU 调度时延与显存占用，手动推导解析梯度是实现极限加速的关键支撑。

亮点与洞察¶

仿射递推矩阵化代数构造型巧：洞察到指数衰减 Hawkes 过程状态方程在齐次增广下的封闭性，将非线性时间序列递推还原为结合律线性算子，为点过程领域引入前缀扫描铺平了道路。
零近似保真度与正交兼容性：与以往依靠截断、变分推断或低秩投影的加速思路截然不同，该方法完全保留了原汁原味的无损极大似然目标，能够与现有的低秩参数化、先验稀疏化正交叠加以进一步提升性能。
广阔的工程迁移价值：此类基于结合律的「增广矩阵 + Blelloch 前缀扫描 + 边界状态跨批接力」范式，可直接借鉴推广至连续时间状态空间模型（SSM）、神经常微分方程（Neural ODEs）以及线性自回归注意力等时序长程依赖计算中。

局限与展望¶

硬件架构依赖性强：由于核心收益完全依赖高并发处理单元，在单核 CPU 或弱算力边缘终端上，前缀扫描的算法常数反而不及轻量级串行递推。
核函数形式受限：当前状态增广与结合律仅严格成立于指数衰减核（Exponential Kernels），对于幂律核（Power-law / Omori-Utsu 核）或非紧凑无记忆核，无法直接写为低秩分块对角矩阵。
未来改进可尝试探索高阶核函数或多项式基函数的指数混合逼近（Sum-of-Exponentials），将并行前缀扫描机制推广至更为复杂的非参数化与非线性 Hawkes 过程推断中。

评分¶

新颖性: ⭐⭐⭐⭐☆ 首次发掘线性指数 Hawkes 过程状态方程的稀疏增广矩阵结合律，成功将并行前缀扫描引入点过程精确推断。
实验充分度: ⭐⭐⭐⭐⭐ 涵盖 6700 万级合成事件、3200 万级社交多峰级联以及 840 万级芝加哥时空犯罪实测，对比了暴力、CPU/GPU串行及 Autograd 完备消融。
写作质量: ⭐⭐⭐⭐⭐ 理论推导紧凑严谨，增广转移算子构造与梯度闭式解析清晰，代码与工程落地指引详实。
价值: ⭐⭐⭐⭐⭐ 彻底打破了霍克斯过程数十年来受困于大规模事件难以精确拟合的算力魔咒，开源 PyTorch 库具备极高的复用价值。

实现配置 / 架构	时间复杂度 (N 维)	峰值显存占用特征	最大可支持事件规模	运行瓶颈与终止原因
Naive 暴力实现	\(\mathcal{O}(N^2)\)	随序列平方急剧暴增	\(N < 10^5\)	严重显存溢出 (OOM) 且运算严重超限
Sequential (GPU 串行递推)	\(\mathcal{O}(N)\)	常数级别（低显存）	极低规模截断	单核计算效率低下，单步耗时极高超时
Sequential (CPU 64核串行)	\(\mathcal{O}(N)\)	常数级别（低显存）	约 \(N \approx 5 \times 10^5\)	递推本质串行，多核无法分摊，耗时超限
Prefix Scan + PyTorch Autograd	\(\mathcal{O}(N/P + \log N)\)	\(\mathcal{O}(N)\) 计算图开销	\(N \approx 10^6\)	计算图庞大且求导解释器开销大，显存溢出
Full Model (Prefix Scan + Manual Grad + Batching)	\(\mathcal{O}(N/P)\)	常数显存 (固定 \(N_b\))	\(N > 6.7 \times 10^7\)	线性稳定扩展，单轮仅需数十秒至数百秒