Massively Parallel Exact Inference for Hawkes Processes

Abstract

Multivariate Hawkes processes are a widely used class of self-exciting point processes, but maximum likelihood estimation naively scales as $O(N^2)$ in the number of events. The canonical linear exponential Hawkes process admits a faster $O(N)$ recurrence, but prior work evaluates this recurrence sequentially, without exploiting parallelization on modern GPUs. We show that the Hawkes process intensity can be expressed as a product of sparse transition matrices admitting a linear-time associative multiply, enabling computation via a parallel prefix scan. This yields a simple yet massively parallelizable algorithm for maximum likelihood estimation of linear exponential Hawkes processes. Our method reduces the computational complexity to approximately $O(N/P)$ with $P$ parallel processors, and naturally yields a batching scheme to maintain constant memory usage, avoiding GPU memory constraints. Importantly, it computes the exact likelihood without any additional assumptions or approximations, preserving the simplicity and interpretability of the model. We demonstrate orders-of-magnitude speedups on simulated and real datasets, scaling to thousands of nodes and tens of millions of events, substantially beyond scales reported in prior work. We provide an open-source PyTorch library implementing our optimizations.

Figures (5)

• Figure 1: A visualization of the Blelloch scan upsweep and downsweep stages performed on an example input sequence $\{a_i\}_{i=1}^8$. Here, $p_{i:j}=\bigodot_{k=i}^j a_k$ denotes a partial aggregate.
• Figure 2: Scaling of our method across various sequence lengths, showing average epoch time (left) and peak memory usage (right). The gray region indicates where batching is enabled to maintain constant memory (with batch size $2^{19}=524288.0$). Dashed black lines indicate the maximum epoch time and memory usage across all runs, respectively. Diagonal gridlines with unit power (slope 1 on log-log plot) are shown for reference.
• Figure 3: Comparison between our method (black) and other implementations of the intensity computation. The red dotted lines indicate the 24-hour walltime limit (86.4 seconds per epoch for 1000 epochs) and the maximum available memory (80 GiB), respectively. Crosses indicate that the next run either ran out of memory or exceeded 24 hours.
• Figure 4: Interaction matrix for all crime types learned from Chicago crime data. Choropleth maps showing influence of particularly influential districts (District 7 and District 10).
• Figure 5: Interaction matrices for selected crime types.