Smooth Neural Point Processes via B-Splines
Michele Bellomo Politecnico di Milano Milan, Italy [email protected]
arXiv:2607.21098v1 [cs.LG] 23 Jul 2026
Alberto Dolara Politecnico di Milano Milan, Italy
Riccardo Ramaschi Politecnico di Milano Milan, Italy
Tomaso Aste University College London London, UK
Abstract Temporal point processes (TPPs) provide a general and flexible framework for modeling sequences of events in continuous time. Neural networks have been successfully employed to model TPPs in a highly expressive and data-driven way. Neural TPPs are typically trained via Maximum Likelihood Estimation (MLE) by minimizing the negative log-likelihood (NLL), which depends on both the conditional intensity function (CIF) and its integral over time, the compensator. Recent neural TPP approaches enable exact evaluation of the NLL without numerical integration. However, these methods typically model the compensator rather than the CIF directly, impose constraints on the neural network architecture, and are computationally expensive during training, as event contributions to the NLL are evaluated sequentially rather than in parallel. In this work, we propose a novel neural TPP model that directly parametrizes the CIF as a non-negative combination of B-spline basis functions, whose coefficients are predicted by a neural network. This formulation enables exact evaluation of the NLL, preserves full flexibility in the neural architecture, allows efficient parallelization during training, and naturally supports CIF smoothness regularization through the integrated squared second derivative. Experiments on both synthetic and real-world datasets show improved computational efficiency and predictive accuracy compared to the reference neural TPP baseline.
1
Introduction
Temporal point processes (TPPs) [1] provide a general and flexible framework for modeling sequences of events in continuous time. Event sequences are common in many disciplines, such as seismology [2], biology [3], epidemiology [4], social sciences [5], and finance [6, 7, 8]. A central object in the theory of TPPs is the conditional intensity function (CIF), which describes the instantaneous expected rate of events given the past history. Classical models such as Poisson [9] and Hawkes [10, 11] processes provide parametric forms for the CIF. These models are attractive due to their interpretability and tractability, but they may be too restrictive to capture complex temporal dependencies observed in real-world event streams. A straightforward extension of classical parametric models is to adopt a black-box approach and rely on neural networks to learn complex temporal dependencies directly from data. Despite its flexibility, this approach presents important computational challenges. Neural TPPs are typically trained via Maximum Likelihood Estimation (MLE) by minimizing the negative log-likelihood (NLL), which depends on both the CIF and its integral over time, the compensator. While evaluating the CIF is Preprint.
straightforward, computing the integral term can be difficult when the intensity is represented by a neural network. Some approaches [12, 13] address this issue by approximating the integral term via numerical quadrature. However, this solution is computationally expensive and may introduce instability during training. Omi et al. [14] proposed to use a network to model the compensator instead of the CIF, and then recover the CIF via automatic differentiation [15]. While this approach effectively avoids numerical integration, it requires enforcing monotonicity of the compensator through architectural constraints [16, 17]. These constraints limit the flexibility of the model and, importantly, prevent parallel evaluation, since the intensity of every time point is computed sequentially. As a result, training and inference are slow, making the approach infeasible in data-intensive settings. In this work, we propose a novel neural TPP model that directly parametrizes the CIF using a non-negative combination of B-spline basis functions, whose coefficients are predicted by a neural network. This formulation enables exact and efficient evaluation of the NLL, eliminating the need for numerical integration. At the same time, it preserves full flexibility in the choice of the neural architecture, as no structural constraints are required to ensure validity of the model. Importantly, the intensity of the time points can be evaluated in parallel over time, leading to significantly improved computational efficiency compared to existing approaches. Finally, the proposed representation naturally enables smoothness regularization of the CIF via penalization of the integrated squared second derivative, providing additional control over the learned dynamics and better generalization.
2
Background
2.1
Temporal Point Processes (TPPs)
Intuitively, we can think of a point process (Tn )n≥1 as an increasing sequence of random times 0 < T1 < T2 < . . . . These random variables represent the occurrence (or arrival) times of events. A point process can be associated with a counting process Nt , a stochastic process that counts the number of arrivals up to time t X Nt = 1{Tn ≤t} . n≥1
An important quantity for a point process is the conditional intensity function (CIF), which represents the instantaneous rate of occurrence of events at time t, given the history up to time t (t excluded): E(Nt+∆t − Nt | Ft ) . ∆t Ft is the filtration of the process, which represents the information available up to time t. In the following, we will use the simplified notation λ(t) := λ(t | Ft ), that omits the explicit dependence of the CIF on the filtration. λ(t | Ft ) = lim + ∆t→0
2.2
Maximum Likelihood Estimation (MLE)
The most general approach for fitting the parameters of a point process is Maximum Likelihood Estimation (MLE), which consists in maximizing the likelihood of the observed sequence of events under the assumed model. The general form of the likelihood for point processes on the time period [0, T ] is "N # T Y L= λ(Ti ) e−Λ(T ) i=1
where Λ(T ) =
RT 0
λ(s) ds is the compensator of the process.
For numerical reasons, it is usually preferable to maximize the log-likelihood ! NT X ℓ= ln(λ(Ti )) − Λ(T ) i=1
2
or, alternatively, to minimize the negative log-likelihood (NLL), i.e., the log-likelihood with the sign reversed. Often, it is more convenient to work with inter-event times τi = Ti − Ti−1 instead of absolute times Ti . In this case, the log-likelihood can be written as ℓ=
NT X
(ln(ϕ(τi )) − Φ(τi ))
i=1
where ϕ(τRi ) = λ(Ti ) is the CIF reparameterized in terms of the time elapsed since the last event, and τ Φ(τi ) = 0 i ϕ(s) ds is the relative integral. 2.3
Random Time Change Theorem
The Random Time Change Theorem states that, given an increasing sequence of time points {T1 , T2 , . . . } and a TPP with compensator Λ(·), the transformed sequence {Λ(T1 ), Λ(T2 ), . . . } is a realisation of a unit rate Poisson process if and only if the original sequence {T1 , T2 , . . . } is a realisation of the TPP defined by Λ(·). By the Random Time Change Theorem, when the observed event times of a TPP are transformed through its compensator, the resulting sequence follows a unit-rate Poisson process. A key property of a unit-rate Poisson process is that its inter-arrival times are independent and exponentially distributed with rate 1, i.e., Φ(τi ) = Λ(Ti ) − Λ(Ti−1 ) ∼ Exp(1), with Φ(τi ) being the cumulative integral of the CIF reparameterized in terms of the time elapsed since the last event. Since an Exp(1) distribution has median log(2), a median-based estimator τ̂i+1 for the waiting time until the next event can be obtained by finding τ̂i+1 such that Φ(τ̂i+1 ) = log(2) . 2.4
(1)
B-splines
Splines are piecewise polynomial functions defined on a partition of an interval [18]. A spline of degree p with knots {xj } is a function that is a polynomial of degree p on each subinterval [xj , xj+1 ) and has continuous derivatives up to order p − 1 at the knots. Spline spaces are dense in a wide class of function spaces, and therefore provide flexible approximations of target functions. In particular, by increasing the number of knots, spline functions can approximate arbitrarily well any target function under mild regularity conditions. A convenient way to represent splines is through a B-spline basis [19]. Any spline function s(x) can be expressed as K X s(x) = wk Bk (x) k=1
where wk ∈ R are coefficients and {Bk (x)}K k=1 are the non-negative B-spline basis functions. The B-spline basis has local support, meaning that each basis function is nonzero only on a limited number of adjacent knot intervals, leading to sparse and numerically stable representations. Moreover, since B-spline basis functions are non-negative, imposing non-negativity constraints on the coefficients wk ≥ 0 is sufficient to ensure that the resulting spline s(x) is non-negative for all x. Spline functions of degree p = 3, known as cubic splines, play a particularly important role. It can be shown that, among all functions defined on an interval [a, b] that interpolate a set of points {(xi , yi )}, cubic splines are the smoothest in the sense of minimizing the roughness Z b P(s) = (s′′ (x))2 dx. (2) a
In the context of nonparametric regression, cubic splines arise as the exact solution to the functional optimization problem n X 2 min (yi − s(xi )) + αP(s) (3) s
i=1
3
where α > 0 controls the strength of the regularization. The roughness penalty P(s) can be computed very efficiently. Indeed, it admits a quadratic form P(s) = w⊤ Rw,
(4)
where w is the vector of spline coefficients and R is a positive semidefinite matrix with entries Rb Rij = a Bi′′ (x)Bj′′ (x) dx that depends only on the chosen basis. Importantly, R can be computed once in advance, and its structure is typically sparse due to the local support of B-spline basis functions, enabling efficient evaluation of the penalty during optimization.
3
Methodology
3.1
Proposed model
We propose to model the CIF ϕ(τ ) on the inter-arrival times τi as a (cubic) B-spline function ϕ(τ ) =
K X
wk Bk (τ ).
k=1
The coefficients wk are predicted by a generic neural network that takes as input the history of past events and possibly additional contextual information. Since B-spline basis functions are nonnegative, non-negativity of the CIF is guaranteed by construction by imposing wk ≥ 0 through suitable output activations of the neural network, for example softplus. A key advantage of this parameterization is that the integral of the CIF, required for the NLL evaluation, can be computed in closed form as Z τ K X Φ(τ ) = ϕ(s) ds = wk Ik (τ ), 0
k=1
where Ik (τ ) are the integrated basis functions, which can be precomputed exactly before training. Knots {kj } can be distributed uniformly over the domain of inter-arrival times, or chosen according to empirical quantiles of the observed inter-arrival time distribution. To avoid oscillations in the learned CIF, one can either use a limited number of knots K, or increase K and add the roughness penalty P introduced in Equation 2 to the NLL in the loss function
−
NT X
(ln(ϕ(τi )) − Φ(τi )) + αP(ϕ).
(5)
i=1
While no optimality result like the regression problem in Equation 3 holds in the context of TPPs, penalizing the roughness of the CIF provides a principled way to control local fluctuations. Moreover, Equation 4 for the efficient computation of the penalization makes its inclusion during training essentially cost-free. 3.2
Multivariate extension
The proposed formulation naturally extends to multivariate TPPs with multiple event types m ∈ {1, . . . , M }, each associated with its own conditional intensity function λm (t). In this setting, each conditional intensity function is modeled as ϕm (τ ) =
Km X
wm,k Bm,k (τ ),
k=1 m where {Bm,k (τ )}K k=1 denotes the B-spline basis associated with event type m, and wm,k > 0 are the corresponding coefficients predicted by the neural network.
The extension is immediate, as it only requires increasing the dimensionality of the network output to PM match the total number m=1 Km of spline coefficients across all event types. 4
The roughness penalty naturally generalizes as a weighted sum over event types P =
M X
⊤ αm wm Rm wm ,
m=1 ⊤
where wm = (wm,1 , . . . , wm,Km ) and Rm is the precomputed roughness matrix associated with the basis of type m. The coefficients αm ≥ 0 control the strength of the smoothness regularization for each event type, allowing different levels of roughness penalization across the corresponding intensity functions. For simplicity, one can also choose to use a shared temporal basis across all event types, instead of defining a distinct spline basis with different knot placements and numbers for each type. 3.3
Goodness of fit
Similarly to [14], we evaluate the goodness of fit of the proposed model on the test data using the mean absolute error (MAE) computed on the median-based estimator τ̂i+1 for the next inter-arrival time described in Section 2.3 with respect to the true observed τi+1 . Our model allows to obtain the estimator τ̂i+1 in a very efficient way. In fact, instead of evaluating the CIF and the compensator in a unique future time as in [14], with a single forward pass our model outputs the spline coefficients from which the entire future trajectories of both the CIF and the compensator are obtained. This enables solving Equation 1 very efficiently, as no additional forward passes of the model are required during the iterations of the root-finding algorithm. Given the monotonicity of the compensator, a simple yet effective approach to solve Equation 1 is the bisection method. In case of a multivariate TPP, the median-based estimator τ̂i+1 for the next inter-arrival time is PM computed using the total CIF ϕ(τ ) = m=1 ϕm (τ ) and the corresponding total compensator PM Φ(τ ) = m=1 Φm (τ ). The estimated next event type can be determined by m̂i+1 = arg maxm∈{1,...,M } ϕm (τ̂i+1 ).
4
Experiments
4.1
Datasets
We test our model on the same synthetic and real datasets used in [14], whose model will be used as the baseline for comparison. Each synthetic dataset consists of a sequence of 100000 time points, generated by different TPP models. Every sequence is then split into train and test sets, using a standard 80%/20% division. The TPPs considered are: • a stationary Poisson process with λ(t) = 1; 2πt • a non-stationary Poisson process with λ(t) = 0.99 · sin( 2000 ) + 1; • a stationary Renewal process, in which the inter-event intervals τi are independent and identically distributed according to a log-normal probability distribution with mean µ = 1.0 and standard deviation σ = 6.0; • a non-stationary Renewal process, whose time points Ti are obtained by first generating independent and identically distributed inter-event times τi′ from a stationary Renewal process (Gamma distribution with mean µ = 1.0 and standard deviation σ = 0.5), and then R Ti rescaling the time by inverting the formula τi′ = Ti−1 r(s) ds, where r(t) is a non-negative 2πt trend function defined as r(t) = 0.99 sin 20000 + 1; P • a Self-correcting process, with λ(t) = exp(t − ti <t 1); n o P • a Hawkes process with single exponential kernel λ(t) = µ+ Ti <t α β exp − β(t − Ti ) and parameters µ = 0.2, α = 0.8, β = 1; 5
• a Hawkes process with multiple exponential kernel λ(t) = µ +
M XX
αj βj exp
n
o − βj (t − Ti )
Ti <t j=1
and parameters M = 2, µ = 0.2, α1 = α2 = 0.4, β1 = 0.4, β2 = 20. We additionally use the two publicly available real-world datasets employed in [14]. For both datasets, we apply the same preprocessing described in [14], and time is measured in hours. • Music dataset [20] contains last.fm users’ listening history in January 2009. 100 sequences from the 100 most active users are created, using the first 80% of events in the sequences for training and the last 20% for testing. • Meme dataset [21] contains popular phrases collected from several online resources. 50 sequences from the 50 most used phrases during August 2008 are created, using the first 40 sequences for training and the last 10 sequences for testing. 4.2
Model architecture and training procedure
We use the same architecture for all datasets: a minimal Transformer with embedding size 64, 4 attention heads, and a single vanilla decoder block as described in [22], with an inner feed-forward layer dimension of 64. The input embedding is constructed as the sum of two components: a sinusoidal positional encoding to represent the chronological order of past events, and a learned embedding of the inter-arrival times obtained through a linear dense layer. The CIF curve over future times is predicted by taking as input the inter-arrival times of the last 20 events. The CIF is parametrized as a cubic B-spline function with 20 internal knots distributed at quantiles of inter-arrival training times. The network outputs the positive coefficients of the B-splines using softplus activation in the last layer. The network is trained minimizing the penalized loss in Equation 5 using Adam optimizer with learning rate 0.001, batch size 256, and maximum training epochs 200. Early stopping is applied based on the validation set, obtained as 20% of the training data, with patience 10 epochs, the NLL without penalization as monitored validation loss, and restoring the best weights at the end of the training. The validation set is also used to select for each dataset the optimal value of the regularization parameter α, evaluated over a range of values from 10−3 to 10−10 , also including the case with no regularization α = 0. All other hyperparameters were chosen to be reasonable for the considered case studies and were not subject to any optimization. Our implementation relies on Keras [23] and KerasHub [24] libraries, and experiments were run on a T4 GPU provided by the Google Colab environment. 4.3
Computational performance evaluation
To ensure a fair comparison with the benchmark in [14], the predictive accuracy results are obtained under the same training setting, using the cut sequence of last inter-event times to predict only the CIF at the next event. However, our spline-based neural TPP model allows parallel evaluation of the CIF over multiple future times. In particular, given a sequence of events, our model is able to return the CIF evaluated at the time of each event in a single forward pass. This significantly reduces the training and testing time, and is especially beneficial when using a Transformer architecture. We compare the training computational cost of the CIF evaluation scheme used in [14] with the parallel multi-evaluation variant of our model. In the latter setting, the model processes longer sequences (in our experiments, of length 100) and directly outputs the CIF at each event time in a single forward pass. A warm-up period can be introduced, excluding the first n events from the loss computation, as they may not have sufficient conditioning history (in our experiments, n = 20). Finally, when using a Transformer architecture, it is also possible to limit the effective context by restricting attention to the most recent n past events through a custom attention mask (n = 20 in our experiments). 6
4.4
Results
Table 1 shows the results of the experiments described in Section 3.3 on the datasets introduced in Section 4.1. For each model, the best regularization parameter α is selected on the validation set, and then the MAE of the median-based estimator is computed on the test set and compared with the one reported by Omi et al. [14]. Table 1: Time prediction experiment results. Type
Dataset
MAE
Omi et al. MAE [14]
Synthetic
Stationary Poisson Non-stationary Poisson Stationary Renewal Non-stationary Renewal Self-correcting Hawkes 1 Hawkes 2
0.692 0.711 0.957 0.406 0.494 0.399 0.947
0.696 0.710 0.894 0.414 0.496 0.848 0.962
Real
Music Meme
0.183 0.135
0.783 0.811
Our model outperforms the baseline in 5 out of 7 synthetic datasets and in both real-world datasets. On the “Hawkes 1” dataset and on the real-world datasets, the improvement in MAE is substantial. In most experiments, the optimal smoothing parameter α was nonzero, suggesting that regularization was indeed beneficial. The effect of regularization on the learned CIF from non-stationary Renewal and Hawkes 2 processes is shown in Figures 1 and 2.
Figure 1: Regularization effect of different smoothing parameters α on the learned CIF from the non-stationary Renewal process.
Figure 2: Regularization effect of different smoothing parameters α on the learned CIF from the Hawkes 2 process. A higher value of the regularization parameter α reduces the fluctuations of the predicted CIF, as is evident in the non-stationary Renewal process in Figure 1. In both cases, values of α not far from 7
the optimal one do not lead to a significant degradation in performance, indicating that the model is reasonably robust to this parameter. Using the T4 GPU provided by the Google Colab environment, one epoch with the standard training setting takes on average 1750 ms on synthetic datasets, whereas with our parallel multi-evaluation approach described in Section 4.3, an epoch takes on average 150 ms, corresponding to a speed-up of approximately 12×.
5
Discussion
In this work, we introduced a novel neural TPP model that directly parametrizes the CIF using a nonnegative combination of B-spline basis functions, whose coefficients are predicted by a generic neural network. This method preserves exact evaluation of the NLL, but, unlike previous approaches such as [14], it allows full flexibility in the choice of the neural architecture and, thanks to the extensive parallelization, enables fast training on large datasets. Moreover, the smoothing regularization through the integrated squared second derivative prevents the oscillatory behavior typically associated with splines and leads to improved generalization. We evaluated the proposed approach on both synthetic and real-world datasets commonly used as benchmarks for neural TPPs. The results show that parametrizing the CIF with B-splines does not degrade predictive accuracy, and in fact leads to a significant improvement, in particular on the real-world datasets. In addition to improved predictive accuracy, the proposed model is computationally more efficient thanks to its flexibility, enabling the use of architectures such as Transformers and parallel processing across sequence events, unlike the sequential computation required by the recurrent neural network in [14]. This advantage becomes substantial when adopting the multi-evaluation approach introduced in Section 4.3, which allows the CIF to be evaluated for all events in the sequences in a single forward pass. This resulted in an empirical training speedup of approximately 12× on the synthetic datasets considered in our experiments.
6
Limitations and future work
The proposed model introduces an additional hyperparameter, namely the regularization coefficient α. While the model is significantly more efficient due to the use of a Transformer and the parallel multi-evaluation strategy, tuning α may partially offset these computational gains. In future work, we aim to investigate more efficient and theoretically grounded strategies for the automatic tuning of this parameter. The proposed methodology shows promising results, but it has so far been evaluated only on datasets of moderate size and using minimal architectures. In future work, we plan to assess the proposed method on in-depth real-world case studies involving large-scale datasets with multiple event types, where other neural TPP approaches may become infeasible and where the computational efficiency and scalability of the proposed framework can be more fully appreciated.
References [1] Daryl J Daley, David Vere-Jones, et al. An introduction to the theory of point processes: volume I: elementary theory and methods. Springer, 2003. [2] Yosihiko Ogata. Seismicity analysis through point-process modeling: A review. Seismicity patterns, their statistical significance and physical meaning, pages 471–507, 1999. [3] Wilson Truccolo, Uri T Eden, Matthew R Fellows, John P Donoghue, and Emery N Brown. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of neurophysiology, 93(2):1074–1089, 2005. [4] Peter J Diggle. Statistical analysis of spatial and spatio-temporal point patterns. CRC press, 2013. 8
[5] George O Mohler, Martin B Short, P Jeffrey Brantingham, Frederic Paik Schoenberg, and George E Tita. Self-exciting point process modeling of crime. Journal of the american statistical association, 106(493):100–108, 2011. [6] Emmanuel Bacry, Iacopo Mastromatteo, and Jean-François Muzy. Hawkes processes in finance. Market Microstructure and Liquidity, 1(01):1550005, 2015. [7] Alan G Hawkes. Hawkes processes and their applications to finance: a review. Quantitative Finance, 18(2):193–198, 2018. [8] Michele Bellomo, Francesco Grimaccia, et al. Hawkes processes in finance: An application in modeling transactions in eua futures. CONTEMPORARY MATHEMATICS, 6(5):5590–5611, 2025. [9] Michele Bellomo. Can we do without the hairdresser? a mathematical solution. The College Mathematics Journal, pages 1–6, 2025. [10] Alan G Hawkes. Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58(1):83–90, 1971. [11] Patrick J Laub, Young Lee, and Thomas Taimre. The elements of Hawkes processes. Springer, 2021. [12] How Jing and Alexander J Smola. Neural survival recommender. In Proceedings of the Tenth ACM International Conference on Web Search and Data Mining, pages 515–524, 2017. [13] Hongyuan Mei and Jason M Eisner. The neural hawkes process: A neurally self-modulating multivariate point process. Advances in neural information processing systems, 30, 2017. [14] Takahiro Omi, Kazuyuki Aihara, et al. Fully neural network based model for general temporal point processes. Advances in neural information processing systems, 32, 2019. [15] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey. Journal of machine learning research, 18(153):1–43, 2018. [16] Joseph Sill. Monotonic networks. Advances in neural information processing systems, 10, 1997. [17] Pawel Chilinski and Ricardo Silva. Neural likelihoods via cumulative distribution functions. In Conference on uncertainty in artificial intelligence, pages 420–429. PMLR, 2020. [18] Simon N Wood. Generalized additive models: an introduction with R. chapman and hall/CRC, 2017. [19] Peter J. Green and Bernard W. Silverman. Nonparametric Regression and Generalized Linear Models. Springer US, Boston, MA, 1994. [20] Oscar Celma. Last.fm dataset. http://ocelma.net/MusicRecommendationDataset/ lastfm-360K.html. [21] Jure Leskovec, Lars Backstrom, and Jon Kleinberg. Meme-tracking and the dynamics of the news cycle. In Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD ’09), pages 497–506. ACM, 2009. https: //snap.stanford.edu/data/memetracker9.html. [22] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017. [23] François Chollet et al. Keras. https://keras.io, 2015. [24] Matthew Watson, François Chollet, Divyashree Sreepathihalli, Samaneh Saadat, Ramesh Sampath, Gabriel Rasskin, Scott Zhu, Varun Singh, Luke Wood, Zhenyu Tan, Ian Stenbit, Chen Qian, Jonathan Bischof, et al. Kerashub. https://github.com/keras-team/keras-hub, 2024. 9