A Representer Theorem for Hawkes Processes via Penalized Least Squares Minimization

AI summary

Overview

• Research area: Machine learning and statistics, specifically kernel methods (reproducing kernel Hilbert spaces) applied to point process models, in particular multivariate linear Hawkes processes.
• Technical level: Advanced. The paper relies on reproducing kernel Hilbert spaces, variational (functional) derivatives, Fredholm integral equations of the second kind, Mercer's theorem, and random Fourier features.
• Scope (1 sentence): The paper proves a new representer theorem for estimating the triggering kernels of linear multivariate Hawkes processes under penalized least squares, and shows that the resulting estimator is closed-form, requires no dual coefficient optimization, and scales far better than the prior state-of-the-art kernel method-based estimator.

What This Paper Is About

Hawkes processes model sequences of events where past events make future events more (or less) likely, and the "triggering kernel" is the function describing how strongly and how quickly one event type influences another. The authors ask how to estimate those kernels nonparametrically, without committing to a fixed parametric shape such as an exponential, while keeping the computation tractable for large event datasets. Their answer is a representer theorem for the non-discretized penalized least squares objective, in which the optimal kernel estimator is a linear combination of transformed "equivalent kernels" whose dual coefficients are all equal to one.

Key Contributions

1. A representer theorem for the non-discretized penalized least squares formulation of linear Hawkes processes. The authors state this is, to their knowledge, the first such result. Prior work (Bonnet & Sangnier, 2025) obtained a representer theorem only via a discretized approximation, which introduced additional optimization over dual coefficients.
2. Equivalent kernels defined by a system of simultaneous Fredholm integral equations, with dual coefficients analytically fixed to unity. In Theorem 1, the estimator of each triggering kernel is ĝ_ij(s) = Σ_{n∈N_i} α_n^{ij} h_j(s, t_n) − μ̂_i ∫_0^T h_j(s,t) dt with α_n^{ij} = 1, so no costly coefficient optimization is required.
3. Closed-form construction of the equivalent kernels via a degenerate kernel approximation. Proposition 3 solves the coupled integral equations in closed form using M feature maps, and Proposition 4 gives the optimal ĝ and μ̂ directly in terms of those feature maps, requiring only additive matrix operations and one matrix inversion.
4. A random Fourier feature implementation with fully analytic integrals. Using shift-invariant kernels and quasi-Monte Carlo random Fourier features (M = 100 in the experiments), all integrals are computed in closed form (Equations 16–17), in contrast to the Riemann/discretization approximation used by Bonnet & Sangnier (2025).

Main Findings

• Unit dual coefficients (Theorem 1): The optimal estimator of each triggering kernel is a linear combination of equivalent kernels evaluated at the observed event times, minus a baseline term, and every dual coefficient α_n^{ij} equals 1. The proof uses the path integral representation (Kim, 2021); an alternative derivation via Mercer's theorem is given in Appendix E.
• Closed-form baseline intensities (Proposition 2): Given the equivalent kernels, the baseline intensities {μ̂_i} have closed forms expressed via integrals of the equivalent kernels, involving |N_i|, the number of observed events on dimension i.
• Closed-form estimators (Propositions 3 and 4): With a degenerate kernel k(s,s') = Σ_{m=1}^M φ_m(s) φ_m(s'), the equivalent kernels are h_j(s,s') = φ(s)^T [(γ^{-1} I_{MU} + Ξ)^{-1} φ̃(s')] restricted to the appropriate M-length slice, where Ξ is a symmetric MU × MU block matrix built from event-pair integrals. Substituting these gives ĝ_ij(s) and μ̂_i in closed form.
• Analytic integrals: For shift-invariant kernels (which include Gaussian, Matérn, and Laplace kernels) and random Fourier features φ_m(s) = sqrt(2/M) cos(ω_m s + θ_m), the required integrals over the support window reduce to closed-form expressions involving Hadamard products and an unnormalized sinc function.
• Complexity advantage: Obtaining the estimators costs O(N̄² M² U² + M³ U³) time, where N̄ = max(|N_1|, ..., |N_U|), and O(M² U²) memory. The prior kernel-based method (Bonnet & Sangnier, 2025) costs O(N̄⁴ U² P) time, where P is the number of iterations for convergence, and O(N̄² U²) memory. The paper also notes that our method requires only a single matrix inversion and avoids tuning convergence criteria and learning rates.
• Why exactness matters: The number of feature maps M needed to approximate a one-dimensional RKHS kernel remains modest regardless of data size, whereas the number of discretization nodes needed for accurate integral evaluation grows with the data size.
• Mutually-exciting scenario results (Table 1, 10 trials, standard errors in brackets). For T = 2000 and average N̄ = 1318: Exp Δ² = 0.29 (0.01), cpu = 46.1 (23.3); Gau 0.20 (0.02), 3.42 (1.33); Ber 0.51 (0.16), 4.06 (2.03); Bonnet 0.21 (0.05), 135 (144); Ours 0.38 (0.15), 3.16 (2.42). For T = 3000, N̄ = 2055: Exp 0.29 (0.00), 74.7 (41.3); Gau 0.19 (0.02), 4.94 (3.00); Ber 0.69 (0.50), 6.61 (4.71); Bonnet 0.19 (0.04), 272 (298); Ours 0.27 (0.06), 4.76 (3.97). For T = 5000, N̄ = 4081: Exp 0.28 (0.00), 134 (64.9); Gau 0.18 (0.01), 12.2 (8.91); Ber 0.30 (0.06), 19.0 (15.4); Bonnet 0.15 (0.02), 1358 (1270); Ours 0.20 (0.06), 10.7 (7.71). For T = 7000, N̄ = 5380: Exp 0.28 (0.00), 180.4 (56.6); Gau 0.18 (0.02), 17.3 (8.00); Ber 0.27 (0.04), 29.0 (16.1); Bonnet 0.14 (0.04), 2070 (1210); Ours 0.16 (0.04), 13.1 (5.55). Bold entries in the table indicate performance not significantly different from the best under the Mann-Whitney U test (Holm, 1979) at p ≥ 0.01.
• Refractory scenario: Table 2 reports a refractory scenario over 10 trials with the same notation as Table 1, but the numeric values are not included in the provided paper content, so those results are not reported here.

Methodology in Plain English

The authors start from the linear multivariate Hawkes process intensity equation, in which each dimension's intensity equals a baseline plus a sum of triggering kernels applied to past events on all dimensions. Rather than maximizing the likelihood, they minimize a least squares contrast that has a quadratic form in the triggering kernels, plus a regularization penalty that constrains the kernels to lie in a reproducing kernel Hilbert space.

They then take the functional derivative of this penalized least squares objective with respect to each triggering kernel and set it to zero. Because the regularization penalty can be rewritten as a double integral involving the inverse kernel operator (a path integral representation), applying the integral operator to the stationarity equation turns it into a set of coupled Fredholm integral equations of the second kind. Crucially, the coefficients of those equations do not depend on the "target" dimension i, so the solution for all (i, j) pairs collapses to a single family of functions h_j, giving the representer form with all dual coefficients equal to one.

Because those integral equations cannot generally be solved directly, the authors approximate the kernel as a finite sum of M feature maps. Under that degenerate form, the coupled equations become a linear system in an MU × MU matrix Ξ, and the solution is written in terms of a single matrix inverse. For shift-invariant kernels they instantiate the feature maps with random Fourier features (quasi-Monte Carlo), which makes every integral appearing in Ξ and in the estimator closed form. In the experiments, each model is fitted on events in [0, 0.8T] and hyperparameters are selected by evaluating the negative log-likelihood (Gau, Ber, Bonnet) or the least squares loss (Ours) on [0.8T, T], searching over γ ∈ {0.1, 0.5, 1.0} and β ∈ {0.5, 1.0, 1.5}. Performance is measured by the integrated squared error Δ² between true and estimated triggering kernels over [0, A] with A = 5, and efficiency by CPU time. Baselines were run from the Python code of Bonnet & Sangnier (2025) under the MIT License; the proposed method used TensorFlow-2.10 on a MacBook Pro with a 12-core Apple M2 Max CPU and the GPU disabled.

Why This Matters

• Impact on research: The result shows that the least squares formulation for Hawkes processes has an unusually elegant structure, namely analytically fixed unit dual coefficients, that the more common likelihood-based or approximately-discretized formulations obscure. It connects the classical representer theorem machinery from kernel methods to the multivariate Hawkes setting without discretization, and gives a scalable alternative to iterative dual optimization.
• Finance: The paper lists finance (Bacry et al., 2015) among the application domains of multivariate Hawkes processes, where cross-excitation between assets, trades, or order arrivals is of interest.
• Neuroscience: The paper lists neuroscience (Gerhard et al., 2017), where mutually interacting spike trains from multiple neurons are a canonical multivariate Hawkes application.
• Social networks: The paper lists social networks (Zhou et al., 2013), where one user's activity can trigger activity by others.
• Seismology: The paper lists seismology (Ogata, 1988), the setting that originally motivated self-exciting point process models of earthquake aftershocks.
• Industry relevance: Because cost drops from O(N̄⁴ U² P) to O(N̄² M² U² + M³ U³) and only a single matrix inversion is needed, the estimator is described as lightweight, stable, and practical for large-scale event data, without convergence criteria or learning rate tuning. The broader class of nonparametric estimation problems the paper cites spans signal and image processing, system control, geostatistics, bioinformatics, and clinical studies.

Future Directions

1. Extending to non-linear Hawkes processes with link functions. The paper's model assumes the identity link, whereas Bonnet & Sangnier (2025) use a soft-plus link φ(x) = log(1 + e^{wx})/w with w = 100 so that intensity stays non-negative and both excitatory and inhibitory interactions can be captured. Whether the unit-dual-coefficient phenomenon survives a non-linear link is left open.
2. Evaluation on real event data. All reported empirical evaluations are on synthetic datasets, so behavior under model misspecification and real-world noise is not established by this work.
3. Practical selection rules for M, A, γ, and β. The experiments fix M = 100, A = 5, and search only over three values each for γ and β; the paper does not report a principled selection procedure or sensitivity analysis.
4. Comparison against a wider set of baselines and a full refractory-scenario analysis. Beyond the truncated Table 2, the paper does not report experiments against neural approaches such as RNN-based (Mei & Eisner, 2017) or Transformer-based (Zuo et al., 2020) event models, nor against the online RKHS regret-minimization approach of Yang et al. (2017), leaving these as natural comparisons.

Target Audience

Researchers and practitioners in machine learning and statistics who work with kernel methods, reproducing kernel Hilbert spaces, or point process models. It will be most useful to readers already comfortable with functional analysis and integral equations, and to applied scientists in finance, neuroscience, social network analysis, and seismology who need a scalable, nonparametric way to estimate interaction kernels from large event streams.

Authors’ abstract

The representer theorem is a cornerstone of kernel methods, which aim to estimate latent functions in reproducing kernel Hilbert spaces (RKHSs) in a nonparametric manner. Its significance lies in converting inherently infinite-dimensional optimization problems into finite-dimensional ones over dual coefficients, thereby enabling practical and computationally tractable algorithms. In this paper, we address the problem of estimating the latent triggering kernels--functions that encode the interaction structure between events--for linear multivariate Hawkes processes based on observed event sequences within an RKHS framework. We show that, under the principle of penalized least squares minimization, a novel form of representer theorem emerges: a family of transformed kernels can be defined via a system of simultaneous integral equations, and the optimal estimator of each triggering kernel is expressed as a linear combination of these transformed kernels evaluated at the data points. Remarkably, the dual coefficients are all analytically fixed to unity, obviating the need to solve a costly optimization problem to obtain the dual coefficients. This leads to a highly efficient estimator capable of handling large-scale data more effectively than conventional nonparametric approaches. Empirical evaluations on synthetic datasets reveal that the proposed method attains competitive predictive accuracy while substantially improving computational efficiency over existing state-of-the-art kernel method-based estimators.