Multivariate Neural Hawkes Process

• The Neural Hawkes Process is a continuous-time recurrent model that replaces fixed additive triggering kernels with a dynamic LSTM architecture for multityped event streams.
• It integrates gating mechanisms to model both excitation and inhibition, enabling nonlinear and order-sensitive dependencies between events.
• Empirical studies show the model improves prediction accuracy and missing data handling, with extensions to graph and spatio-temporal domains for practical applications.

Searching arXiv for recent and foundational papers on Neural Hawkes processes and related graph/spatio-temporal extensions. {"query":"\"Neural Hawkes\" multivariate Hawkes process arXiv", "max_results": 10} {"query":"\"Graph Hawkes Neural Network\" temporal knowledge graphs arXiv", "max_results": 10} {"query":"\"Multivariate Spatio-Temporal Neural Hawkes Processes\" arXiv", "max_results": 10} Multivariate Neural Hawkes Process denotes a neurally self-modulating multivariate point process for sequences of typed events in continuous time. In the formulation introduced by Mei and Eisner, an observed stream $\{(t_i,k_i)\}_{i=1}^n$, with event types $k_i\in\{1,\dots,K\}$, is modeled through mark-specific conditional intensities, but the classical additive Hawkes kernel is replaced by a continuous-time recurrent architecture whose hidden state jumps at event times and decays analytically between them. This construction preserves the generative semantics of multivariate point processes while allowing future intensities to depend on past events through nonlinear, order-sensitive, and potentially inhibitory dynamics rather than only through positive additive excitation (Mei et al., 2016).

1. Formal setting and relation to classical Hawkes models

A multivariate point process on $[0,T]$ is specified by typed events $(t_i,k_i)$ and mark-specific intensities

$$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$

with total intensity $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$. For an observed sequence, the log-likelihood is

$$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$

The event term rewards high intensity at observed event times; the integral term accounts for the probability of no events over infinitesimal intervals (Mei et al., 2016).

The classical multivariate Hawkes process uses linear additive excitation. With exponential triggering kernels,

$$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$

or, in the more general exponential form discussed in the Neural Hawkes literature,

$$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$

This formulation enforces positive, additive, exponentially decaying influence from past events. A decomposable self-modulating variant relaxes the sign constraints on $\mu_k$ and $k_i\in\{1,\dots,K\}$0, then restores nonnegativity with a learned scaled softplus,

$k_i\in\{1,\dots,K\}$1

where $k_i\in\{1,\dots,K\}$2 remains additive in the past. That extension allows inhibition and inertia but still retains decomposability in the pre-activation intensity (Mei et al., 2016).

The multivariate Neural Hawkes Process departs from both classical Hawkes and decomposable self-modulation by replacing the additive decomposition itself with a continuous-time recurrent state. The shift is conceptual as well as parametric: past events no longer act through a fixed sum of kernels, but through a latent dynamical system that conditions future event probabilities on the entire event history.

2. Continuous-time LSTM parameterization

The defining mechanism of the Neural Hawkes Process is a continuous-time LSTM (CT-LSTM) with memory cells $k_i\in\{1,\dots,K\}$3 and hidden state $k_i\in\{1,\dots,K\}$4. At each event time $k_i\in\{1,\dots,K\}$5, the model reads the event type $k_i\in\{1,\dots,K\}$6, encoded in the original formulation as a one-hot vector, and updates gates analogous to those of a standard LSTM: input, forget, candidate, output, barred input, barred forget, and positive decay-rate gates. The output of these gates sets a new cell value $k_i\in\{1,\dots,K\}$7, a new steady-state target $k_i\in\{1,\dots,K\}$8, an output gate $k_i\in\{1,\dots,K\}$9, and positive decay rates $[0,T]$0 (Mei et al., 2016).

Between events, the state evolves analytically. On the interval $[0,T]$1,

$[0,T]$2

and

$[0,T]$3

Thus the cell decays exponentially toward a learned steady state, but the decay is dimension-specific and controlled by event-dependent gates. On each open interval, $[0,T]$4, $[0,T]$5, and $[0,T]$6 are piecewise constant, while $[0,T]$7 and $[0,T]$8 have closed-form trajectories (Mei et al., 2016).

Type-specific intensities are generated from the hidden state through another scaled softplus: $[0,T]$9 This positivity transform ensures strictly positive intensities while allowing signs and magnitudes of $(t_i,k_i)$0 to produce either excitation-like or inhibition-like effects before the nonlinearity. Positive decay rates are likewise enforced by a softplus or other positive transfer, guaranteeing stable exponential decay between events (Mei et al., 2016).

This architecture supports several behaviors absent from linear Hawkes models. Because different hidden dimensions can decay at different rates toward different steady states, $(t_i,k_i)$1 may evolve non-monotonically between events. Delayed responses can arise when a drifting cell trajectory crosses the steep region of $(t_i,k_i)$2, and order-sensitive interactions can emerge through the recurrent state update rather than through pairwise kernel superposition. In this sense, the model is neurally self-modulating rather than merely nonlinear at the output layer.

3. Likelihood, optimization, prediction, and simulation

Training uses the standard multivariate point-process log-likelihood, but the integral term

$(t_i,k_i)$3

must now be evaluated for a state trajectory defined by recurrent jumps and continuous decay. Although $(t_i,k_i)$4 is piecewise analytic, the original implementation estimates the integral by Monte Carlo rather than by analytic or interval-wise quadrature. Uniform samples $(t_i,k_i)$5 are drawn, $(t_i,k_i)$6 is evaluated, and

$(t_i,k_i)$7

This yields an unbiased estimator of the integral and of its gradient by exchanging gradient and expectation. In practice, $(t_i,k_i)$8 is chosen proportional to the number of observed events $(t_i,k_i)$9: $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$0 during training and $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$1 for evaluation when feasible (Mei et al., 2016).

Backpropagation proceeds through both the algebraic gate updates at event times and the closed-form exponential decays between events. Because the between-event dynamics are analytic, no ODE solver is required. For a model with $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$2 parameters, one sequence with $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$3 observed events and $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$4 Monte Carlo samples costs $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$5, with $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$6 in practice. Optimization in the original experiments used Adam with default settings and early stopping on a development set; unregularized log-likelihood was used because $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$7 regularization did not help in pilot runs (Mei et al., 2016).

Inference for the next event preserves the generative point-process semantics. Conditional on the history before event $$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$8, the next arrival time has density

$$\lambda_k(t\mid H_t)=\lim_{\Delta t\to 0}\frac{\Pr\!\big(N_k(t+\Delta t)-N_k(t)=1\mid H_t\big)}{\Delta t},$$9

and, given that the next event occurs at time $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$0, its type is distributed as

$\lambda(t)=\sum_{k=1}^K \lambda_k(t)$1

Minimum Bayes risk prediction under squared loss uses the conditional expectation of $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$2, while type prediction integrates the time-dependent type posterior over the arrival-time density. The paper evaluates these quantities by Monte Carlo sampling of $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$3 from $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$4, with shared time samples across types used for variance reduction (Mei et al., 2016).

Simulation uses Ogata’s thinning. A constant upper bound is constructed for the aggregate intensity on each inter-event interval by bounding the hidden-state contributions to each $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$5; proposed times are drawn from a homogeneous Poisson process at the bound rate, accepted with probability $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$6, and then assigned a type according to $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$7. Because hidden states and cells are bounded on each interval, such constant bounds can be built analytically (Mei et al., 2016).

4. Expressivity, missing data, and empirical behavior

The empirical motivation for the multivariate Neural Hawkes Process is that many event streams exhibit more than monotone mutual excitation. The CT-LSTM parameterization supports excitation and inhibition, non-additive and order-dependent effects, long-range temporal interactions, delayed responses, and non-monotonic intensity evolution between events. These properties arise from the recurrent state dynamics rather than from a fixed library of triggering kernels (Mei et al., 2016).

In synthetic experiments, the model fits streams generated by Hawkes, decomposable self-modulating, and neural self-modulating processes, and it is reported to be essential when the data-generating process itself is neural. On neural-generated intensities, it predicts true intensities accurately, with approximately $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$8 MSE versus $\lambda(t)=\sum_{k=1}^K \lambda_k(t)$9 for Hawkes. On retweet cascades, it consistently outperforms Hawkes and the decomposable self-modulating model at all training sizes. On MemeTrack, where $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$0 and sequences are short, the shared hidden state yields dramatic gains: cross-entropy per event is reported around $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$1 for Neural Hawkes versus $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$2 for the decomposable self-modulating model, and the hidden-state parameter sharing reduces parameters relative to Hawkes from $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$3 to $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$4 (Mei et al., 2016).

On prediction tasks from finance, ICU diagnoses, and Stack Overflow badges, the model achieves lower type error on most splits, while time RMSE is mixed and has no consistent winner. This pattern is consistent with the architecture’s emphasis on richer dependence of event type on detailed history. It does not imply uniformly better time forecasting under all evaluation protocols, but it does indicate that type prediction benefits substantially from continuous-time latent state dynamics (Mei et al., 2016).

A further empirical result concerns missing data. When subsets of event types are censored from Hawkes-generated data and all models are trained as if the observed streams were complete, the Neural Hawkes Process achieves better likelihood than Hawkes on every censored dataset. The explanation given in the original study is that the hidden state can encode a sufficient statistic for predicting future observations and can implicitly approximate the Bayesian posterior over unobserved history. This suggests a practical advantage in partially observed multivariate streams, where linear additivity is too rigid to represent “explaining away” patterns among latent causes (Mei et al., 2016).

5. Graph and spatio-temporal generalizations

A central limitation of the original multivariate formulation is that it requires enumerating event types. In graph-structured event streams, where an event may be a directed, labeled triple $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$5, the number of possible types scales as $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$6, which makes a separate intensity per type infeasible. The Graph Hawkes Neural Network (GHNN) addresses this by preserving the neurally self-modulating idea while replacing per-type parameters with shared entity and relation embeddings and a continuous-time LSTM over relevant graph history (Han et al., 2020).

For object prediction in a temporal knowledge graph, GHNN aggregates concurrent neighbors at each past timestamp by mean pooling,

$$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$7

feeds the concatenated input $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$8 into a cLSTM, and defines the candidate-object hidden state by gating the decayed memory with $$\log p(\{(t_i,k_i)\}_{i=1}^n) = \sum_{i=1}^n \log \lambda_{k_i}(t_i) - \int_0^T \sum_{k=1}^K \lambda_k(t)\,dt.$$9. The resulting triple intensity is

$$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$0

with $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$1 a scaled softplus. Directionality is handled by separate histories for outgoing and incoming edges, and relation labels enter both the cLSTM input and the intensity parameterization (Han et al., 2020).

Training in GHNN combines two link-prediction cross-entropies, for $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$2 and $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$3, with a time-to-event MSE term based on the expected next occurrence time of a specific triple. Integrals are approximated numerically via the trapezoidal rule, Adam with weight decay is used for optimization, and histories are truncated to maximum length $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$4 for efficiency. On ICEWS14, GHNN reports filtered link prediction MRR $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$5, Hits@1 $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$6, Hits@3 $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$7, Hits@10 $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$8, and time-prediction MAE $$\lambda_k(t)=\mu_k+\sum_{i:t_i<t}\alpha_{k,k_i}\,\beta e^{-\beta(t-t_i)},$$9 days; on GDELT, it reports MRR $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$0, Hits@1 $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$1, Hits@3 $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$2, Hits@10 $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$3, and MAE $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$4 hours. The paper also emphasizes a timestamp-aware filtering metric that removes only corrupted events occurring at the same timestamp as the query (Han et al., 2020).

A different extension incorporates spatial coordinates directly into the latent state. The Multivariate Spatio-Temporal Neural Hawkes Process (MSTNHP) models events $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$5 and defines

$$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$6

where the cell decays jointly in time and distance: $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$7 The time-space integral in the likelihood is approximated by Monte Carlo sampling over $$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$8, again avoiding an ODE solver because the between-event dynamics remain closed form. In simulation studies, the model is reported to recover sensible temporal intensity shapes and spatial intensity maps, whereas the temporal-only counterpart fails when the data are inherently spatio-temporal despite similar likelihood values. In a terrorism application to Pakistan, the fitted intensities show TTP with the highest overall intensity for most of the year, localized BRA surges, and meaningful date-specific spatial reversals between groups (Chukwuemeka et al., 27 Feb 2026).

Together, these extensions show that the multivariate Neural Hawkes construction is not limited to a flat catalog of event types. It can be reparameterized over graph structure, spatial coordinates, and other covariates so long as the latent state remains a continuous-time causal summary of the past and intensities remain positive point-process hazards.

6. Stability theory, high-dimensional analysis, and limitations

The original Neural Hawkes work is primarily algorithmic, but several theoretical results for generalized multivariate Hawkes processes bear directly on neural parameterizations. In the high-dimensional nonlinear Hawkes framework,

$$\lambda_k(t)=\mu_k+\sum_{h:t_h<t}\alpha_{k_h,k} e^{-\delta_{k_h,k}(t-t_h)}.$$9

the key stability condition is a Lipschitz-spectral bound. If each $\mu_k$0 is $\mu_k$1-Lipschitz and

$\mu_k$2

then existence of a stationary solution follows from Brémaud–Massoulié stability theory. The same analysis develops a thinning-based coupling construction that works beyond mutual excitation, accommodates mixed-sign kernels, and yields concentration inequalities for second-order statistics and smoothed cross-covariances in high dimensions (Chen et al., 2017).

This theory is not a CT-LSTM analysis, but it is directly relevant to neural Hawkes parameterizations in which the nonnegative link is Lipschitz and past events enter through learned filters. The paper explicitly notes that softplus is $\mu_k$3-Lipschitz and logistic is at most $\mu_k$4-Lipschitz, so stability can be promoted by controlling kernel $\mu_k$5 norms and the Lipschitz constants of neural links. This suggests one principled route toward stability-regularized training for neural point-process models, although the theorem itself is stated for the generalized convolutional form rather than for hidden-state dynamics (Chen et al., 2017).

Related graph-limit theory further clarifies how heterogeneity alters multivariate Hawkes behavior on networks. For inhomogeneous random graphs, the macroscopic limit intensity solves a nonlinear Volterra equation,

$\mu_k$6

and, in the linear case, the large-time threshold is governed by the spectral radius $\mu_k$7 of the graphon-induced operator $\mu_k$8, with stationary fields becoming position-dependent when the connectivity kernel $\mu_k$9 is nonconstant. A plausible implication is that graph-aware neural variants should inherit analogous sensitivity to structured heterogeneity even when their finite-dimensional parameterization differs from graphon-based Hawkes limits (Agathe-Nerine, 2021).

Several limitations recur across the literature. In the base Neural Hawkes model, the integral term is estimated by Monte Carlo, so gradients inherit sampling noise; interpretability is weaker than in pairwise-kernel Hawkes models; and the original paper identifies immediate or batch events, richer drift processes, hybrid decomposable-plus-neural designs, sparse interaction priors, and transformer-style alternatives as natural extensions (Mei et al., 2016). In GHNN, survival summation over all entities is $k_i\in\{1,\dots,K\}$00, histories are truncated, and concurrent events within a slice are aggregated by mean pooling under a conditional-independence assumption, so more expressive slice-level dependencies remain open (Han et al., 2020). In MSTNHP, the spatial decay is isotropic, time and space enter additively in the exponent, no explicit boundary correction is used, and integration cost scales with the number of Monte Carlo samples and the study-region area (Chukwuemeka et al., 27 Feb 2026).

Taken together, these results place the multivariate Neural Hawkes Process at the intersection of point-process likelihood theory, continuous-time recurrent modeling, and structured event forecasting. Its distinctive contribution is to retain the hazard-based semantics of multivariate Hawkes models while replacing fixed triggering kernels with a continuous-time latent dynamical system that can be shared, extended, and reparameterized across increasingly structured event domains.