Deep graph kernel point processes over networks

In plain words

Many systems produce timestamped events that happen at places in a network: a traffic sensor counts a jam, a patient's lab value goes out of range, a theft is logged near a pub. A point process model answers the question "how likely is the next event, given what came before?" The classic Hawkes model says each past event adds a bump of extra risk that fades with time. The shape of that bump is called the influence kernel.

This paper keeps that picture but makes the bump graph-aware. The effect of an earlier event depends on how long ago it happened and on which node it happened at, and the nodes are linked by a graph. A graph neural network describes that dependence, so the model can pick up effects that spread several steps along the links, effects that raise the risk, and effects that lower it. When the links are not known, an attention layer can learn which ones matter. Training includes a penalty that pushes its predicted risk to stay positive. The authors report better predictions of when and where the next event happens than several earlier neural models.

TL;DR

GraDK writes the Hawkes influence kernel as a low-rank sum over pairs of temporal basis functions and graph filter bases §3.2 §3.3 eq. 4. Both kinds of basis are small neural networks or learned matrices, and two graph filters are run: L3Net when the graph is known and GAT when it must be learned §3.3 app. C. Training uses the negative log-likelihood or a least-squares loss, plus a log-barrier penalty that keeps the intensity positive on a grid §3.4. On the paper's own tables the method gives the lowest Time MAE and Type KLD on most data sets, but not always the best log-likelihood (SAHP-G is higher on Sepsis and Theft) tab. 1 tab. 3. All results are author-reported.

Problem

Discrete events arrive at nodes of a graph, and the goal is to predict future event times and types and to recover how events influence each other. Earlier deep kernel point processes handle time but not graph structure, and they work in Euclidean space §1 §2. Graph point processes that model all pairwise influence need parameters that grow with the square of the node count; the authors argue that a low-rank graph-filter kernel needs only on the order of R·C·|V| parameters, where C is the average neighbourhood size §3.3 app. C. Classical multivariate Hawkes models use scalar node-to-node coefficients, which the authors say cannot capture time-varying or multi-hop influence §2.

Method

The model is a marked point process on a graph G = (V, E). Each event has a time t and a node v. The edges define where influence may flow; their direction and strength are learned §3.1.

1. Intensity. The conditional intensity at node v is a base rate plus a sum over past events of an influence kernel k(t', t, v', v) eq. 2. The kernel may be negative, which allows inhibition §2 §3.3.
2. Kernel decomposition. The kernel is written as a sum over D basis pairs, each with a weight σ_d, a temporal basis g_d and a graph basis h_d eq. 3. The authors note that a generalisation of Mercer's theorem for asymmetric kernels justifies this form §3.2.
3. Temporal basis. Each g_d is a sum over L pairs of functions: ψ_l(t') for the past event time and φ_l(t − t') for the lag, combined with weights β_dl. The displacement t − t' rather than t is used because it gives a low-rank approximation to general kernels §3.2. ψ_l and φ_l are independent fully connected networks of the same architecture §3.2.
4. Graph basis. Each h_d is a weighted sum of R localized graph filter bases B_r(v', v) §3.3. The paper allows algebraic or spectral filters and names ChebNet, L3Net, GAT and the GPS graph transformer as options §2 app. C. The experiments use L3Net, where B_r reaches the neighbours up to order o_r, and GAT, where B_r is a learned attention-like matrix app. B app. C.
5. Combined kernel. Multiplying out the two expansions gives a kernel with coefficients α_rl = Σ_d σ_d β_dl γ_dr eq. 4. The nominal rank D disappears; only the temporal rank L and the graph rank R need to be chosen §3.3.
6. Efficient evaluation. The paper evaluates ψ_l only at event times, evaluates φ_l on a uniform grid with linear interpolation, cuts the kernel at a maximum lag τ_max, and computes the integral term from cumulative sums on the same grid. The cost is O(n + |U_t|) for n events and a grid U_t app. E.
7. Losses. The negative log-likelihood (NLL) and a least-squares (LS) loss are both defined §3.4 app. D. Both integrate the intensity over time and sum over nodes.
8. Positivity. Because the kernel may be negative, a log-barrier term on a dense grid of times and nodes keeps the intensity above a bound b. The objective is the loss plus a barrier weighted by 1/w §3.4 app. E.
9. Training. Algorithm 1 is stochastic gradient descent over mini-batches. After each epoch the barrier weight w is multiplied by a constant a > 1, and the bound b is set from the smallest intensity seen on the grid alg. 1 app. F.
10. Rank selection. L and R are chosen by validation log-likelihood, with splits made at trajectory or subsequence level, and by inspecting the size of the learned α_rl §3.3 app. G.

Key equations (as in the paper; equation numbers are from the main text):

\lambda(t, v) = \mu + \sum_{(t', v') \in \mathcal{H}_t} k(t', t, v', v) \quad \text{[eq. 2]}

k(t', t, v', v) = \sum_{d=1}^{D} \sigma_d\, g_d(t', t - t')\, h_d(v', v) \quad \text{[eq. 3]}

k(t', t, v', v) = \sum_{r=1}^{R} \sum_{l=1}^{L} \alpha_{rl}\, \psi_l(t')\, \varphi_l(t - t')\, B_r(v', v) \quad \text{[eq. 4]}

The NLL and LS objectives, the barrier and the grid-based forms are given without equation numbers in §3.4 and in app. D and app. E.

# Graph kernel intensity lambda(t, v)  [§3.2] [§3.3] [eq. 2-4]
# Input:  history: list of (t_j, v_j) with t_j < t; query node v
# Output: scalar intensity; positivity is enforced by the barrier [§3.4]

def intensity(t, v, history, mu, alpha, psi, phi, B, tau_max):
    # Start from the base rate mu  [eq. 2]
    total = mu
    for t_prev, v_prev in history:
        # Lags beyond tau_max contribute nothing  [app. E]
        lag = t - t_prev
        if lag > tau_max:
            continue
        for r in range(R):
            for l in range(L):
                # Temporal part psi_l(t') * phi_l(t - t'), times the
                # graph part B_r(v', v), weighted by alpha[r, l]  [eq. 4]
                temporal = psi[l](t_prev) * phi[l](lag)
                total += alpha[r, l] * temporal * B[r][v_prev, v]
    return total

# Objective for one sequence: NLL plus log-barrier  [§3.4] [app. D] [app. E]
# Input:  events [(t_i, v_i)] on [0, T]; grid_t (U_t); grid_b (U_bar, t)
# Output: scalar loss; w is the barrier weight, b the lower bound

def objective(events, grid_t, grid_b, b, w):
    # Integral term: intensity summed over nodes and integrated in time
    integral = 0.0
    for t in grid_t:
        for v in nodes:
            # history_before(t) returns the events with time below t
            integral += intensity(t, v, history_before(t)) * dt
    # Log term at each observed event  [§3.4]
    log_sum = 0.0
    for t_i, v_i in events:
        log_sum += log(intensity(t_i, v_i, history_before(t_i)))
    # Barrier over the penalty grid, averaged  [§3.4]
    barrier = 0.0
    for t in grid_b:
        for v in nodes:
            barrier += log(intensity(t, v, history_before(t)) - b)
    barrier = -barrier / (len(grid_b) * len(nodes))
    return (integral - log_sum) + barrier / w

# Epoch update of barrier bound b and weight, as the code does it
# [alg. 1] [code: src/utils.py#L345-L355]
# The code stores t = 1/w, so the penalty enters as log-barrier / t
# Input:  lam_min (smallest grid intensity this epoch), t, b
# Output: updated t, b

def update_barrier(lam_min, t, b):
    # Shipped configs set b_bar = 0.0125 and b_upp = -0.0125, so this
    # test is lam_min < 0  [code: config/train_16-node-graph_mle.yaml#L22-L23]
    if lam_min < b_bar + b_upp:
        # Violation: restart with a strong barrier, bound set from lam_min
        t = t_init
        b = lam_min - b_bar
    else:
        # No violation: relax the barrier geometrically, capped at t_upp
        t = min(t * t_mul, t_upp)
        b = b_upp
    return t, b

Hyperparameters

Values stated in the paper for GraDK's own runs, with the code value where the paper is silent.

Baseline training settings, as stated in the paper app. H (learning rate / batch size):

Further baseline settings app. H: RMTPP embedding dimension 32 (chosen from 32, 64, 128); FullyNN embedding 64, two hidden layers of width 64, input embedding 10 (and an extra mark network the authors added); DNSK temporal and mark basis rank 1 and 3; THP-S and SAHP-G use their own code defaults.

Results

Author-reported. The synthetic setups are described in §4.1 and app. H; the real data sets are the Atlanta traffic counts (5 sensors), Sepsis (13 event types), California wildfire (25 clusters) and Valencia theft (52 pubs) §4.2. Each cell gives the mean and, in parentheses, the standard error over three runs. Time MAE is the error in predicted event frequency and Type KLD the divergence of predicted event types, both from 100 generated sequences §4.1. Hardware is not stated.

Synthetic data, main comparison tab. 1

Synthetic data, other sets app. H, tab. H3

Real data, Sepsis and Theft tab. 3

Real data, traffic and wildfire app. H, tab. H4

Ablation: spatio-temporal baselines on wildfire and theft (Euclidean distances) tab. 2

Displacement ablation: kernel with and without t − t' (16-node and wildfire, L3Net + NLL) app. H.5, tab. H5

Training cost (seconds per epoch, app. E, tab. E1; hardware not stated)

Rank search app. G, tabs. G1 and G2: testing ℓ on 16-node, 50-node, wildfire and theft for L in {1, 2} and R in {1, 2, 3}. The chosen cells are L = 1 with R = 3 (16-node, wildfire, theft) and L = 1 with R = 2 (50-node). Full values are in the appendix tables; the wildfire row for L = 1 is −3.900, −3.654, −3.625 for R = 1, 2, 3 app. G, tab. G1.

Other points the paper reports: NLL and LS give similar runtimes §4.1; LS is described as slightly more numerically stable §4.1; the learned wildfire and theft graphs are in fig. 3; the Sepsis dynamic kernel is in fig. 5; traffic and Sepsis graphs are in app. H. On theft the learned kernel shows multi-hop structure where DNSK does not fig. 4. Table 2 includes MHP as a classical benchmark, and GraDK uses 411 parameters against 4742 for DNSK-stpp tab. 2.

Paper vs code

Official repository xieyao14/graph-kernel at 0f9de2ac… (MIT; copyright Yao Xie). The paper links it in app. H. Code notes below cite files in that repository.

• Optimiser and learning rate. The paper uses Adam with learning rate 1e-2 and batch 32 app. H. The code uses Adadelta with the config learning rate utils.py. The shipped 16-node config sets lr 0.1 and batch 64 for 500 epochs train_16-node-graph_mle.yaml; the wildfire config sets lr 0.1, batch 32 and 1000 epochs train_wildfire_mle.yaml. Algorithm 1 in the paper describes plain SGD app. F, so the paper itself gives two different update rules.
• Extra regulariser not in the paper. The code adds a penalty on the time basis at long lags: the squared basis values weighted by a term that is zero below half of tau_max and rises to 1 at tau_max model.py, also added to the integral term model.py. It is added to the training loss with weight lam_reg, which is 100 for 16-node train_16-node-graph_mle.yaml and 1e4 for wildfire train_wildfire_mle.yaml utils.py. Neither §3.4 nor app. D contains it.
• Batch normalisation of the loss. The paper's NLL and LS are averaged over M sequences, with a 1/M factor §3.4 app. D. The code sums the log-likelihood and the integral over the whole mini-batch model.py model.py. The barrier is a mean over its grid utils.py. The relative weight of the penalty therefore depends on batch size.
• Log clamp. The code computes log(max(λ, 1e-5)) model.py. The paper's formulas take the plain log §3.4.
• LS normalisation. The code divides the LS objective by the window length T model.py. The paper's ℓ_LS has no such factor §3.4.
• Barrier schedule and bound. The paper sets b_temp to the smallest grid intensity minus a small ε, and takes the epoch's b from it, so b is tuned toward zero as learning progresses alg. 1 §3.4. The code keeps b at b_upp = −0.0125 unless a violation happens; then it sets b = lam_min − b_bar and resets t to t_init utils.py. Here ε is the constant b_bar = 0.0125 utils.py. A violation is lam_min < b_bar + b_upp, which is lam_min < 0 with the shipped values utils.py. The weight schedule is t = 1/w, multiplied by 1.3 per epoch and capped at 10, which is the paper's "multiply w by a" with a cap alg. 1 utils.py.
• Algorithm sign. Algorithm 1 writes the objective as −ℓ(θ) + (1/w)p, while §3.4 defines ℓ_NLL as the negative log-likelihood alg. 1 §3.4. Read literally, the two disagree on sign.
• GAT basis. The paper defines GAT with learned weights W and a on node features, with attention coefficients computed from the features of the two nodes app. C (item 3). The code does not compute attention from features. It learns a full n × n matrix of logits, applies a softmax over rows, masks it and normalises the columns to sum to one model.py. The mask is all ones at initialisation, so every pair is trainable model.py. The parameter count is O(R|V|²), not the O(RC|V|) the paper gives for masked attention app. C.
• L3Net order. The paper says the r-th filter has order r, for r = 1 to R app. G. The code uses orders 0 to R − 1, so the first filter is the identity train_synthetic_data_L3.py train_real_data_L3.py.
• Filter initialisation. Filters start as Gaussian noise with standard deviation sqrt(1 / (R² · mean degree)), masked to the k-hop neighbourhood model.py. The paper does not state an initialisation for B_r.
• Basis networks. The code uses softplus with beta = 100 model.py, Xavier initialisation with gain 0.5 and bias 0.01 model.py, and feeds the time inputs divided by T model.py. The paper gives none of these.
• Kernel coefficients. The code initialises α uniformly in [−1e-2, 1e-2] (init_std = 1e-2) and applies no softplus to them model.py train_synthetic_data_L3.py. The paper does not state an initialisation.
• Base rate. μ is a fixed constant, not learned model.py. The paper lists μ in θ §3.4 but gives no value for its experiments.
• Data split. The paper uses 80/20 splits app. H. The 16-node config uses 500 sequences, split 400 and 100 train_16-node-graph_mle.yaml; the paper describes 1,000 sequences app. H. The evaluation script tests on sequences 900 to 999 evaluation.py. Three different test sets are implied.
• Time MAE. The evaluation computes the mean event count per unit time for the generated sequences minus the same quantity for the test set, with no absolute value and no per-sequence error evaluation.py. The paper describes a mean absolute error §4.1.
• Grid sizes. The paper says only that the grids are dense or of order n §3.4 app. E. The synthetic script uses 50 to 100 for the evaluation and integration grids and 100 for the penalty grid train_synthetic_data_L3.py; the real-data script uses 100 and 200 train_real_data_L3.py.
• Data in the repository. The repository includes the 16-node synthetic set and the wildfire set only . The other seven data sets are not in it. The theft data is described as proprietary app. H.

ChebNet is implemented in the code model.py but is not used in the experiments app. C. The GPS graph transformer is not in the repository. The baseline models (RMTPP, FullyNN, DNSK, THP-S, SAHP-G) are not in it either, so their numbers could not be checked against code.

Pitfalls

• Unstated constants. The paper gives no value for τ_max, μ, the barrier constants (w, a, b_0, ε), the penalty and integration grids, the epoch count, the seed or the initialisation §3.4 app. E app. F. Implementing the paper as written will not reproduce the code's numbers without these. The code values are listed above.
• Positivity is only enforced on a grid. The barrier is evaluated on a finite grid of times and nodes §3.4. Between grid points the intensity can go negative, so a finer grid or a stricter bound is needed in practice. The paper notes that negative intensities are common early in training app. F.
• Log of a small intensity. The NLL takes the log of the intensity at each event. The paper says this can cause numerical trouble near zero §3.4. The code clamps at 1e-5 model.py, which biases the likelihood when intensities are tiny.
• Batch scaling. Because the code sums losses over the batch and the paper averages them, the barrier and regulariser weights are not directly comparable across batch sizes. Use the code's convention, or rescale.
• Cost. The paper's runtime claim depends on grid sizes: O(n + |U_t|), with the grid chosen around n app. E. The timing table is per epoch on unstated hardware app. E. The grid choice trades accuracy against cost app. E.
• Rank search versus final ranks. The rank search covers R = 1 to 3 and L = 1 to 2 app. G, but the final models use R = 4 on three data sets app. H. The paper does not say how these were chosen.
• Rank selection by coefficient size. The paper describes choosing L and R from the size of α_rl, which it calls a heuristic §3.3 app. G. The validation split is at trajectory level, and the paper states that adjacent subsequences may still be dependent app. G.
• Likelihood is not the target. On Sepsis and Theft, SAHP-G has the highest testing ℓ tab. 3. The authors argue that predictive error matters more here app. H.4. If the goal is the likelihood, the baseline may be the better choice.
• Wildfire sequence count. The paper gives 2,428 events and an average length of 436 per one-year sequence app. H, which implies about 5.6 sequences. The wildfire config uses 5 train_wildfire_mle.yaml.
• Negative influence needs the barrier. Inhibition is represented by negative kernel values. Without the barrier the intensity can become negative and the log term becomes undefined §3.4.
• Displacement matters. Using t instead of t − t' (the non-displacement kernel) gives worse Time MAE and Testing ℓ on both sets tested app. H.5, so the parametrisation should not be simplified.
• Results are author-reported. The paper reports three runs per cell with standard errors tab. 1. THP-S and SAHP-G use the default parameters of their own code implementations app. H.2.

Lineage

• Builds on: Hawkes process §2; the deep non-stationary kernel of Dong et al. (2022), which uses the log-barrier and the displacement parametrisation §2 §3.2 §3.4; L3Net (Cheng et al., 2020), GAT, ChebNet and the GPS transformer as graph filter options §2 app. C. None of these is a card in this wiki yet.
• Compared against: RMTPP, FullyNN, DNSK-mtpp, THP-S, SAHP-G (graph-aware attention), MHP, ETAS and DNSK-stpp §4 app. H.
• Superseded by: none known.