Wan et al.: Efficient Temporal Point Processes via Monotone Alternating Splines

1 Introduction↩︎

Event sequences are prevalent across various domains, such as criminology [1], [2], social networks [3], [4], finance [5], [6], seismology [7], [8], and neuroscience [9], [10]. These sequences are typically characterized by irregular timestamps in continuous time, which sets them apart from conventional regular time series data. Understanding the dynamics within these sequences provides valuable insights for forecasting future events. Temporal point processes (TPPs) offer a statistical framework for modeling such event sequences.

The key to characterizing a TPP is the next timestamp’s distribution given its history. One of the most common approaches to parameterize a TPP is to specify a conditional intensity function (CIF). Numerous methods [11]–[13] leverage deep learning approaches to model the CIF and recover the model parameters via maximum likelihood estimation (MLE). However, CIF-based paradigms suffer from a persistent computational bottleneck. Evaluating the closed-form log-likelihood of a TPP sequence requires integrating the CIF over a time window. If the CIF is based on complicated neural networks, the integral must be approximated by numerical integration, such as Monte Carlo [14]. This introduces cumulative errors and severely limits computational efficiency during both training and inference.

To entirely bypass integration, an equivalent alternative is to directly model the Cumulative Conditional Intensity Function (CCIF) [15]. By definition, the CIF is the exact derivative of the CCIF. Consequently, the required CIF can be recovered via automatic differentiation [16]. CCIF-based parameterization provides an exact, integral-free paradigm that eliminates integration errors and accelerates likelihood evaluations.

Despite its structural elegance, the current CCIF parameterization is severely constrained. Because the underlying CIF is strictly non-negative, a valid CCIF must mathematically be a continuous, monotonically increasing function. Current literature uniformly enforces this constraint using Monotone Neural Networks (MNNs) [15], [17]. Typically, MNNs are implemented using Multi-layer Perceptrons (MLPs) with strictly positive weights and non-decreasing activation functions. However, in this paper, we identify that the MNN architecture inevitably encounters three structural deadlocks: (1) The Convexity Restriction: Using purely convex activations (e.g., ReLU, Softplus) restricts the CCIF to be strictly convex, fundamentally preventing the predicted CIF from decaying after excitation. (2) The Saturation Limit: Even when bounded, non-convex activations (e.g., Sigmoid or Tanh) are used instead of convex ones in the MNN, the network faces a different issue: the predicted CCIF will stop growing and eventually stagnate. Consequently, the network becomes incapable of modeling long-horizon point processes where events persistently occur over time. (3) Violations of CCIF Requirements: Finally, to forcibly bypass the convexity of MNNs, some methods apply non-convex temporal transformations in advance, such as \(t\rightarrow\log(t)\). However, such transformations actually violate the fundamental mathematical requirements of the CCIF. Additionally, they heuristically force the intensity function to decay over time, which severely limits the model’s flexibility.

To resolve these structural deadlocks, we propose a new framework–Monotone Alternating Splines (MAS)–to model CCIF. MAS replaces the rigid MNN with piecewise monotone splines, which offer greater generality, flexibility, and efficiency, and enjoy strong theoretical guarantees for both fitting and generalization. Specifically, MAS separates the estimated CCIF into two distinct sub-components:

Firstly, the MAS interpolation explicitly divides the CCIF into multiple fine-grained intervals and leverages monotone splines to accurately fit each segment. The parameters of each interpolation segment are determined by a sequence encoder (e.g., a Transformer), and then the monotone splines are connected with \(C^1\) continuity. Through this fine-grained segmented fitting, MAS can accurately fit complex TPP sequences. Theoretically, it possesses an approximation capability comparable to existing parameterization methods, alleviating the MNN convexity restriction.

Secondly, because standard piecewise monotone splines cannot maintain global monotonicity, the MAS extrapolation component is introduced. Specifically, outside the interpolation intervals, MAS employs a simple monotonically increasing function (e.g., a linear or linear-plus-exponential function) that smoothly connects with the tail of the internal interpolation part. The extrapolation guarantees global monotonicity over the unbounded time horizon. Furthermore, it allows MAS to recover several classical TPP baselines, thereby securing robust generalization for distant future predictions.

Beyond the architecture, we establish a comprehensive theoretical foundation for the MAS framework. First, we analyze the generalization error of MAS, explicitly decomposing it into the interpolation error, the extrapolation error, and the complexity error. Based on this established generalization bound, we theoretically validate how the extrapolation component enhances the generalization capability of MAS. Furthermore, by comparing the approximation error bounds between MAS and MNNs, we theoretically prove that MAS possesses superior fitting and generalization capabilities compared to MNNs. Finally, we provide theoretical guidance on how to configure optimal MAS hyperparameters.

Our contribution can be summarized as follows:

• We expose the fundamental limitations of existing MNNs for CCIF modeling. We demonstrate that MNN architecture inherently traps the model in three structural deadlocks, critically restricting model expressivity for complex TPP dynamics.

• We propose Monotone Alternating Splines. By introducing a highly expressive interpolation component and a globally stable extrapolation component, MAS satisfies all mathematical requirements for CCIF modeling, delivering higher flexibility and computational efficiency.

• We establish a comprehensive theoretical foundation for MAS. By decomposing and bounding MAS’s generalization error, we prove its significant superiority over MNNs. Moreover, we explicitly provide a theoretical guideline for MAS hyperparameter selection. Extensive experiments across diverse synthetic and real-world event datasets validate the effectiveness of MAS.

The remainder of this paper is organized as follows. 2 reviews related work. 3 formulates the preliminaries of CCIF and its required conditions. Next, we analyze the limitations of MNNs in CCIF modeling in 4. 5 details the mathematical formulation of the MAS framework. In 6, we detail the theoretical analyses comparing structural approximation bounds and formalizing generalization guidelines. 7 presents the experimental evaluations. Finally, 8 concludes the paper.

2 Related Work↩︎

This section reviews related work on TPPs and monotone splines.

2.1 Temporal Point Processes↩︎

Temporal Point Processes provide the mathematical framework for modeling asynchronous event sequences in continuous time [18]. Classical statistical methods typically parameterize the CIF with predefined functional forms under specific assumptions. For example, the homogeneous Poisson process models independent event occurrences [19]; the Hawkes process captures self-exciting phenomena, where past events increase the probability of future arrivals [5], [20]; and the Self-correcting process models stress accumulation that decreases after an event [21]. However, these classical models, whether optimized via Maximum Likelihood Estimation (MLE) or Bayesian frameworks, still struggle to capture the highly nonlinear, multimodal temporal dependencies common in modern large-scale datasets.

2.2 CIF-based Neural TPPs↩︎

To improve structural expressiveness, deep learning has shifted TPP modeling toward data-driven neural methods. Early approaches used recurrent architectures to parameterize the exact intensity function. Specifically, RMTPP [11] and NHP [12] map historical event embeddings into a dynamic continuous-time CIF. Later methods leveraged self-attention to capture long-range interactions, leading to autoregressive models such as THP [13] and SAHP [22]. In parallel, continuous-time state formulations based on Ordinary Differential Equations [23], [24] and Stochastic Differential Equations [25] were proposed to flexibly model intensity trajectories between events. Despite improved predictive performance, these CIF-based models share a key computational bottleneck: exact Negative Log-Likelihood (NLL) evaluation requires integrating the CIF over continuous time. Since neural CIFs generally lack closed-form integrals, their integration must rely on Monte Carlo sampling [14] or numerical quadrature [26], which introduces approximation error and slows training.

2.3 Alternative TPP Parameterizations and CCIF Modeling↩︎

To bypass the CIF bottleneck, recent work has explored alternative parameterizations and estimation strategies. Some methods avoid likelihood integration entirely by using Normalizing Flows [27], or Score Matching frameworks [28], [29]. Meanwhile, other approaches model intermediate variables equivalent to the CIF to avoid integration. For example, some methods model the conditional distribution of the next event time using mixture distributions [30], [31], while others directly model the CCIF, which is the integral of the CIF, and then obtain the exact CIF via automatic differentiation. Existing CCIF-based models [15], [17], [32] typically adopt MNNs [33] to enforce the positive-intensity prerequisite. Specifically, our method focuses on improving current CCIF approaches.

2.4 Monotone Splines in Continuous Modeling↩︎

Monotone splines have a long history in statistical numerical analysis as tools for modeling monotone data [34]. Early work established methods for constructing smooth, monotonically increasing trajectories using piecewise basis functions [35]. Later methods introduced Rational Linear Splines (RLS) [36], Rational Quadratic Splines (RQS) [37], and Rational Cubic Splines (RCS) [38]. These rational monotone splines guarantee strict monotonicity and high-order smoothness, and they also support stable analytical inverses. Moreover, each segmented interval can be controlled independently by localized parameters without affecting global monotonicity, which provides high flexibility [39]. In this work, we identify the structural deadlocks of previous CCIF modeling methods and leverage monotone splines to introduce MAS as a more expressive alternative with stronger approximation guarantees.

3 Preliminaries↩︎

In this section, we provide the mathematical background on TPPs, CCIF, and how to structurally model the CCIF sequence following certain events.

3.1 Temporal Point Processes↩︎

A TPP is a probabilistic model that describes sequences of events occurring within a continuous time window \([0, T]\). A realization of a multivariate TPP is represented as an event sequence \(S = \{(t_n, k_n)\}_{n=1}^{N}\), where \(N\) is a random variable denoting the total number of events. Here, \(0 < t_1 < \dots < t_N < T\) are the event timestamps in continuous time, and \(k_n \in \{1,\dots,K\}\) is the mark (or event type) associated with the \(n\)-th event.

TPPs are primarily characterized by the CIF \(\lambda_k^*(t \mid \mathcal{H}_{t^-})\), which defines the instantaneous occurrence rate of a type-\(k\) event given the historical data \(\mathcal{H}_{t^-} = \{(t_n, k_n) : t_n < t\}\): \[\lambda_k^*(t) = \lim_{\delta_t \to 0} \frac{p(\text{type-}k\text{ event} \in [t,t+\delta_t] \mid \mathcal{H}_{t^-})}{\delta_t},\] where the asterisk indicates conditioning on the history up to but not including time \(t\). Given the CIF, the log-likelihood of a multivariate TPP for an observed sequence \(S\) is: \[\log p(S) = \sum_{n=1}^{N} \log\lambda_{k_n}^*(t_n) -\int_0^T \lambda^*(t) dt, \label{likelihood}\tag{1}\] where \(\lambda^*(t)=\sum_{k=1}^K\lambda^*_k(t)\) is the total intensity across all event types. Model parameters are estimated by maximizing the log-likelihood in 1 . Additionally, given the history \(\mathcal{H}_{t_n}\), the probability density function of the next event timestamp \(t_{n+1}\) is: \[p(t_{n+1} \mid \mathcal{H}_{t_n}) = \lambda^*(t_{n+1}) \exp\left(-\int_{t_{n}}^{t_{n+1}} \lambda^*(t)dt\right).\] We can estimate the expected next event timestamp \(\widehat{t}_{n+1}\) and the most likely mark \(\widehat{k}_{n+1}\) using the following estimators: \[\widehat{t}_{n+1}=\int_{t_n}^\infty tp(t \mid \mathcal{H}_{t_n})dt, \;\;\widehat{k}_{n+1}=\arg\max_k\frac{\lambda_{k}^*(\widehat{t}_{n+1})}{\lambda^*(\widehat{t}_{n+1})}. \label{predicting32t95n431}\tag{2}\]

3.2 CCIF Parameterization↩︎

A major drawback of standard CIF-based modeling is that the integral \(\int_0^T \lambda^*(t) dt\) in 1 generally lacks a closed-form solution when parameterized by complex neural networks, thereby requiring computationally inefficient numerical integration. An alternative, mathematically equivalent formulation is to model the CCIF directly: \(\Lambda_k^*(t) = \int_0^t \lambda_k^*(\tau) d\tau\). Under this framework, 1 can be rewritten as: \[\log p(S) = \sum_{n=1}^{N} \log\frac{d}{dt}\Lambda_{k_n}^*(t_n^-) - \Lambda^*(T), \label{CCIF-likelihood}\tag{3}\] where \(\Lambda^*(t) = \sum_{k=1}^K \Lambda_{k}^*(t)\), and \(\frac{d}{dt}\Lambda_{k_n}^*(t_n^-)\) represents the left-hand derivative of \(\Lambda_k^*(t)\) evaluated at \((t_n, k_n)\). Incorporating the CCIF allows MLE to bypass integral computations entirely. The required instantaneous intensity is easily and precisely recovered via automatic differentiation.

3.3 Modeling CCIF after Certain Events↩︎

Following each observed event \(t_1,t_2,\cdots,t_N\), the continuous global function \(\Lambda^*(t)\) can be decomposed and be expressed as: \[\Lambda^*(t) = \sum_{n=0}^{N} {I}_{(t_n,t_{n+1}]}\Lambda^{*(n)}(t), \label{eq:32decomposes}\tag{4}\] where \(t_0=0\), \(t_{N+1}=T\), and \({I}_{(t_n,t_{n+1}]}\) is an indicator function taking the value 1 if \(t \in (t_n, t_{n+1}]\) and 0 otherwise. Here, \(\Lambda^{*(n)}(t) = \int_{t_n}^t \lambda^*(\tau) d\tau + \Lambda^{*(n-1)}(t_n)\) represents the CCIF following the \(n\)-th event, with the base case \(\Lambda^{*(-1)}(0) = 0\). We emphasize that the natural support domain of \(\Lambda^{*(n)}(t)\) is \((t_n, \infty)\) because, during inference, the subsequent event \(t_{n+1}\) is unknown and could theoretically occur indefinitely far in the future. During training, \(t_{n+1}\) is observed, allowing us to cleanly truncate \(\Lambda^{*(n)}(t)\) to \((t_n, t_{n+1}]\) via the indicator function, subsequently concatenating \(\Lambda^{*(n)}(t)\) and \(\Lambda^{*(n+1)}(t)\) to reconstruct the exact global CCIF \(\Lambda^*(t)\).

4 Failures of Monotone Neural Networks↩︎

In this section, we investigate the principles of using MNNs for CCIF modeling. We show that positive-weighted MLPs restrict MNN expressiveness and prevent them from modeling complex TPP dynamics.

4.1 MNNs for TPP Modeling↩︎

We first explain the MNN-based parameterization of CCIF. Specifically, an \(L\)-layer MNN \(\Lambda^*(t)\) can be expressed as: \[\begin{align} &\boldsymbol{h}^{(1)} = \sigma\!\left(\boldsymbol{W}^{(1)} t + \boldsymbol{b}^{(1)}\right), \\ &\boldsymbol{h}^{(l)} = \sigma\!\left(\boldsymbol{W}^{(l)} \boldsymbol{h}^{(l-1)} + \boldsymbol{b}^{(l)}\right), \quad l=2,\dots,L-1, \\ &\Lambda^*(t) = \boldsymbol{W}^{(L)} \boldsymbol{h}^{(L-1)} + \boldsymbol{b}^{(L)}, \end{align}\] where \(t\) is the elapsed time, \(\sigma(\cdot)\) is the element-wise activation function, \(\boldsymbol{W}^{(l)}\) and \(\boldsymbol{b}^{(l)}\) are the connection weights and biases, respectively. The corresponding CIF is derived as: \[\label{equ:32MNN32CIF} \lambda^*(t) = \frac{\partial \Lambda^*(t)}{\partial t} = \boldsymbol{W}^{(L)} \left( \prod_{l=L-1}^{1} \mathrm{diag}\!\left[\sigma'(\boldsymbol{z}^{(l)})\right] \boldsymbol{W}^{(l)} \right).\tag{5}\]

The core motivation behind using MNNs for modeling CCIF is grounded in the following theorem.

While 1 ensures MNN’s monotonicity, it simultaneously traps the model in three structural deadlocks, as discussed in the following subsection.

4.2 Three Structural Deadlocks of MNNs↩︎

To construct a practical MNN, designers must make choices regarding the non-decreasing activation function \(\sigma(\cdot)\) in 1. However, we demonstrate that these choices force the MNN into one of three structural deadlocks: the Convexity Restriction, the Saturation Limit, or the Violations of CCIF Requirements.

4.2.1 Convexity Restriction↩︎

First, when equipped with standard convex activation functions, the resulting MNN becomes strictly convex, preventing the model from capturing any decaying intensity patterns. While unbounded and convex activations like ReLU or Softplus are favored by modern deep learning paradigms, these common activation functions conflict with the MNN architecture in 1.

2 reveals a severe side-effect: positive linear combinations of convex functions preserve \(\Lambda^*(t)\)’s convexity strictly. If common convex activation functions are applied in MNNs, the magnitude of the CIF \(\lambda^*(t)\) can only monotonically increase or remain completely flat. This severely limits the model’s expressive power. In real-world dynamics, such as an exponentially decaying Hawkes process, the CIF spikes after an event occurrence and then smoothly decays (i.e., \(\lambda^{*\prime}(t) < 0\)). An MNN bounded by 2 inherently fails to fit such a decreasing CIF, regardless of its depth or width.

4.2.2 Saturation Limit↩︎

Second, to bypass the aforementioned convexity restriction, an intuitive approach is to replace ReLU with S-shaped activations (e.g., Tanh or Sigmoid), which offer concave regions (\(\sigma'' < 0\)), thus permitting the derivative \(\lambda^*(t)\) to decline. Unfortunately, this exposes the MNN to the following saturation limit.

As is implied in 3, non-convex activation functions force the CCIF \(\Lambda^*(t)\) into an enforced upper limit \(M\), precluding the mechanism from modeling persistent event streams. For any ordinary sequence demonstrating a baseline rate \(\mu > 0\), the mathematically correct total expected events must approach infinity over time (\(\lim_{t \to \infty} \Lambda^*(t) = \infty\)). Consequently, Tanh- or Sigmoid-based MNNs fundamentally fail to represent a stationary long-term trend; the learned intensity \(\lambda^*(t)\) will invariably exhaust itself and force the predicted sequence to permanently freeze.

4.2.3 Violations of CCIF Requirements↩︎

Since using either convex or bounded, non-convex activation functions leads to the shortcomings of MNN-based \(\Lambda^*(t)\), the remaining option is to apply a non-convex transformation \(t\rightarrow F(t)\) to time \(t\) in 5 in advance, yielding: \[\Lambda^*(t) = \mathrm{MNN}(F (t)),\] where \(\mathrm{MNN}(\cdot)\) is a positive weighted MLP. For example, FullyNN [15] leverages a log-transformation \(F(t) = \log(t)\). Because \(\log t\) is intrinsically concave, this pre-processing step avoids the Convexity Restriction without requiring the \(\mathrm{MNN}(\cdot)\) itself to possess non-convex activations.

Nevertheless, this mathematical tactic still introduces inevitable structural flaws. It violates the boundary condition of CCIF and restricts the expressive power of the model:

• Violation of the Initial Physical Condition. Pre-transforming \(t\) permanently introduces an unavoidable intercept at \(t=0\), causing inaccurate NLL optimization. By definition, the cumulative integration at exact time zero must be zero: \(\Lambda^*(0) = \int_0^0 \lambda^*(\tau) d\tau = 0\). However, evaluating the deformed CCIF sequence at \(t \to 0^+\) yields a limit breakdown: \[\Lambda^*(0) = \lim_{t \to 0^{+}} \mathrm{MNN}(\log t) = \mathrm{MNN}(-\infty)\neq 0.\] Because MNN is a nested composition of affine layers and non-negative activations, its value at \(-\infty\) is unpredictable and can only be approximated during initialization. This inevitably leaves a persistent positive intercept \(\Lambda^*(0) = C \neq 0\). Such a non-zero starting point directly interferes with the exactness of the likelihood computation in 4 .

• Enforced Heuristic Decay Constraints. Additionally, deforming time \(t\) restricts the intrinsic trend of the derived CIF, severely degrading the network’s capability to model diverse TPP behaviors. According to the chain rule, CIF \(\lambda^*(t)\) is tied to the derivative of \(\log(t)\): \[\lambda^*(t) = \frac{\partial \mathrm{MNN}(\log t)}{\partial (\log t)} \cdot \frac{1}{t}=O(t^{-1}).\] This mathematical entanglement dictates an \(O(t^{-1})\) deterministic decay pattern acting as a global multiplier. While this \(1/t\) allows the MNN to roughly approximate TPPs with decaying CIFs, it contradicts other basic TPP behaviors. For instance, in a Self-Correcting process, the CIF accumulates and increases over time between events. The inherent \(1/t\) penalty suppresses this upward trend, making the model poorly suited for ascending dynamics.

Consequently, MNN-based CCIF methods face three fundamental deadlocks. As shown in 1, continuous-time MNNs inevitably suffer from expressive limitations and boundary violations, regardless of activation design or time deformation. To overcome this bottleneck, MAS replaces MNNs with piecewise monotone splines for explicit CCIF parameterization, yielding a more general and flexible framework.

5 Methodology↩︎

In this section, we introduce a novel CCIF modeling method: MAS, and detail its implementation. For complex TPPs, the exact CCIF after event \(t_n\), \(\Lambda^{*(n)}(t)\), can be highly intricate across the interval \([t_n, \infty)\). Directly fitting this complete trajectory using a single elementary function or a rigid MNN is generally suboptimal.

Consequently, a natural idea is to divide the entire CCIF \(\Lambda^{*(n)}(t)\) into smaller segments that are easier to fit. Inspired by this idea, we propose a flexible parameterization framework termed MAS to model the CCIF. MAS consists of two parts: an interpolation component to fit the recent CCIF dynamics right after an event, and an extrapolation component to model the long-term CCIF tail. Formally, MAS models the CCIF after the \(n\)-th event as: \[\label{eq:MAS95for95CCIF} \Lambda_\theta^{*(n)}(t)= \underbrace{\sum_{m=1}^M {I}_{(w^{(n)}_{m-1},w^{(n)}_m]}(t) f^{*(n)}_{m,\theta}(t)}_{\text{MAS interpolation}} + \underbrace{{I}_{(w^{(n)}_{M},\infty)}(t) g^{*(n)}_\theta(t)}_{\text{MAS extrapolation}},\tag{6}\] where the interval \((t_n, \infty)\) is partitioned into \(M+1\) contiguous segments. The first \(M\) segments form the interior \((w^{(n)}_0, w^{(n)}_1], \dots, (w^{(n)}_{M-1}, w^{(n)}_M]\), with the first knot starting exactly at the event time \(w^{(n)}_0 = t_n\). The \(m\)-th piecewise function \(f^{*(n)}_{m,\theta}(t)\) models the MAS interpolation within its respective interval, while the function \(g^{*(n)}_\theta(t)\) represents the MAS extrapolation over \((w^{(n)}_M, \infty)\). The vector \(\theta\) encapsulates the underlying neural network parameters driving the model.

5.1 MAS Interpolation↩︎

The MAS interpolation directly determines the near-term variation trend of the CCIF. By segmenting \(\Lambda_\theta^{*(n)}(t)\), the MAS interpolation provides high flexibility. However, because \(\Lambda_\theta^{*(n)}(t)\) is the integral of a positive CIF, a valid MAS interpolation must satisfy the following requirements:

(1) Global Monotonicity: Since the intensity \(\lambda^*(t) \ge 0\), \(\Lambda_\theta^{*(n)}(t)\) must be non-decreasing. For any \(t_n < t_1 \le t_2 < \infty\), it must hold that \(\Lambda_\theta^{*(n)}(t_1) \le \Lambda_\theta^{*(n)}(t_2)\).

(2) \(C^1\)-Continuity Between Timestamps: For any continuous interval strictly between two events \(t_n < t < t_{n+1}\), \(\Lambda^*(t)\) must be continuous and possess a continuous derivative \(\Lambda^{*\prime}(t)\).

(3) \(C^0\)-Continuity At Timestamps: At each timestamp \(t_n\), the CCIF is continuous (\(\Lambda^{*}(t_n^-) = \Lambda^{*}(t_n^+)\)), but its derivative is permitted to be discontinuous (\(\Lambda^{*\prime}(t_n^-) \neq \Lambda^{*\prime}(t_n^+)\)).

The second and third conditions ensure that the CIF remains continuous between any two consecutive timestamps but may exhibit jumps at the timestamps, which is a core characteristic of many history-dependent TPPs, such as the Hawkes process and the Self-Correcting process. To ensure the generated \(\Lambda_\theta^{*(n)}(t)\) perfectly satisfies monotonicity and \(C^1\)-continuity within \((t_n, t_{n+1})\), the adjacent functions in MAS must be seamlessly connected at the partition knots (i.e., \(w^{(n)}_1, \dots, w^{(n)}_M\)) in both function values and first derivatives: \[\label{eq:conditions95for95MAS} {!}{ \begin{align} f^{*(n)}_{m,\theta}(w^{(n)}_{m}) = f^{*(n)}_{m+1,\theta}(w^{(n)}_{m}),& f^{*(n)\prime}_{m,\theta}(w^{(n)}_{m}) = f^{*(n)\prime}_{m+1,\theta}(w^{(n)}_{m}), \\ f^{*(n)}_{M,\theta}(w^{(n)}_{M}) = g_{\theta}^{*(n)}(w^{(n)}_{M}),& f^{*(n)\prime}_{M,\theta}(w^{(n)}_{M}) = g_{\theta}^{*(n)\prime}(w^{(n)}_{M}). \end{align}}\tag{7}\]

Piecewise monotone splines perfectly satisfy these requirements, because their analytical forms are directly parameterized by the endpoint values and derivatives. Taking the Rational Quadratic Spline (RQS) as an example, given an interval \([w_{m-1}, w_{m}]\), boundary values \(\{y_{m-1},y_m\}\), and boundary derivatives \(\{\delta_{m-1},\delta_{m}\}\), the \(m\)-th piecewise RQS is: \[\label{eq:RQS95example} f_{m}(t) = y_{m-1} + \frac{\left(y_{m} - y_{m-1}\right) \left[s_m \tau^2 + \delta_{m-1} \tau (1 - \tau)\right]}{s_{m} + \left(\delta_{m} + \delta_{m-1} - 2s_{m}\right) \tau (1 - \tau)},\tag{8}\] where \(s_m=\frac{y_{m}-y_{m-1}}{w_m-w_{m-1}}\) and \(\tau=\frac{t-w_{m-1}}{w_m-w_{m-1}} \in [0, 1]\). By substituting \(\tau = 0\) (i.e., \(t = w_{m-1}\)) and \(\tau = 1\) (i.e., \(t = w_{m}\)) into Eq. 8 and its corresponding derivative, we obtain: \[\label{eq:RQS95property} \begin{align} f_{m}(w_{m-1}) &= y_{m-1}, \quad &f_{m}(w_{m}) &= y_{m}, \\ f'_{m}(w_{m-1}) &= \delta_{m-1}, \quad &f'_{m}(w_{m}) &= \delta_{m}. \end{align}\tag{9}\] This implies that an RQS spline piece is uniquely determined by its boundary constants \((y_{m-1}, y_m)\) and derivatives \((\delta_{m-1}, \delta_{m})\). Therefore, as long as adjacent splines share the same structural parameters \(\{(w_m, y_m, \delta_m)\}_{m=0}^M\), they are connected with exact \(C^1\)-continuity. Consequently, all existing monotone splines can act as the basis for MAS interpolation. In Appendix B, we discuss several specific variations and their implementation.

Parameterization of MAS Interpolation: Next, we discuss how the MAS interpolation is implemented. In practice, the defining parameters of MAS—the interpolation knots \(\{w^{(n)}_m\}_{m=1}^M\), cumulative function values \(\{y^{(n)}_m\}_{m=1}^M\), and derivatives \(\{\delta^{(n)}_m\}_{m=1}^M\)—are all unknown and must be derived from historical data using autoregressive or recurrent encoders. Specifically, at timestamp \(t_n\), we employ an encoder (e.g., a Transformer [13], parameterized by \(\theta_1\)) to compress the event history \(\mathcal{H}_{t_n}\) into an embedding \(\mathbf{h}_n\). Then, an MLP (parameterized by \(\theta_2\)) followed by a Softplus activation function maps \(\mathbf{h}_n\) to strictly positive step intervals (\(\Delta w_m, \Delta y_m, \delta_m\)). The knots are obtained cumulatively: \[\begin{align} w^{(n)}_m = w^{(n)}_{m-1} + \Delta w_m,\;\; y^{(n)}_m = y^{(n)}_{m-1} + \Delta y_m. \end{align}\] To ensure the seamless concatenation of consecutive CCIF sequences at event boundaries (the \(C^0\)-continuity requirement), we set \(w^{(n)}_0 = t_n\) and \(y^{(n)}_0 = \Lambda^{*(n-1)}(t_n)\). For adjacent interpolation segments, the network only predicts a single common derivative \(\delta^{(n)}_{m}\) at knot \(w^{(n)}_m\), natively forcing \(f^{*(n)}_{m,\theta}\) and \(f^{*(n)}_{m+1,\theta}\) to share the identical gradient.

Fitting Ability of MAS Interpolation: By introducing piecewise monotone splines, the fitting capability of MAS is significantly enhanced. Unlike MNNs, the theoretically superior expressiveness of MAS can be rigorously supported. Specifically, shorter interpolation intervals allow the MAS interpolation to approximate the ground-truth \(\Lambda^{*(n)}(t)\) more accurately. Formally, we obtain the following theorem.

4 states that MAS fits any ground-truth CIF and CCIF provided that the interpolation intervals are sufficiently small, which implies its fitting performance is comparable to any other TPP parameterization.

5.2 MAS Extrapolation↩︎

Next, we discuss another necessary component, MAS extrapolation, and how it is implemented. 4 assumes that the MAS interpolation support can extend to infinity, i.e., \(w^{(n)}_{M} \to \infty\). However, in practice, this is not possible. Therefore, an additional monotonic extrapolation function \(g^{*(n)}_\theta(t)\) must be used to model \(\Lambda^{*(n)}(t)\) over \((w^{(n)}_{M},\infty)\). It is worth noting that most monotone splines (without any additional tail design) are based on polynomials, which do not guarantee monotonicity outside the interpolation interval. Taking the RQS in 8 as an example, when \(s^{(n)}_m\) is smaller than \(\delta^{(n)}_{m-1}\) and \(t\) is sufficiently large, \(f^{(n)}_{m}(t)\) becomes monotonically decreasing and can even take negative values. See 3 for an illustration. Thus, an alternative parameterization is required for \(g^{*(n)}_\theta(t)\) in 6 to ensure the global monotonicity of \(\Lambda^{*(n)}(t)\) over \((w^{(n)}_{M},\infty)\).

In implementation, \(g^{*(n)}_\theta(t)\) can be a linear function \(g^{*(n)}_\theta(t) = a t + b\) with \(a,b > 0\), or a linear-plus-exponential function \(g^{*(n)}_\theta(t) = c + b t - \exp(-a t)\) with \(a,b,c > 0\). Similar to the MAS interpolation, the parameters \(a\), \(b\), and \(c\) are predicted from the history embedding \(\mathbf{h}_n\) by the MLP. Notably, when connecting the interpolation and the extrapolation at the final knot \(w^{(n)}_M\), we must guarantee 7 . To ensure this smoothness easily, we first determine the extrapolation function \(g^{*(n)}_\theta(t)\) from the network. Then, we strictly assign the final knot target as \(y^{(n)}_M = g^{*(n)}_{\theta}(w^{(n)}_M)\) and the final derivative target as \(\delta^{(n)}_M = g^{*(n)\prime}_{\theta}(w^{(n)}_M)\). According to the properties of monotone splines in 9 , substituting these targets shapes the final interpolation piece \(f_M\), satisfying the strict \(C^1\)-continuity of \(\Lambda^{*(n)}(t)\).

Equivalence to Classical TPPs: As an important hyperparameter in 6 , a larger \(w^{(n)}_{M}\) increases the contribution of piecewise monotone spline interpolation, enhancing the model’s fitting performance. Conversely, a smaller \(w^{(n)}_{M}\) amplifies the influence of the extrapolation component, improving the model’s generalization performance (further discussed in 6). Notably, when \(w^{(n)}_{M} = t_n\), MAS reduces to classical TPP models.

Similar to the two corollaries above, other classical TPPs can also be represented within the MAS framework. By incorporating the piecewise monotone spline interpolation component, MAS extends classical TPPs and significantly enhances the model’s flexibility.

5.3 Efficiency of MAS↩︎

MAS also enjoys superior efficiency compared to MNN-based methods. During the model training process using MLE, MAS does not require automatic differentiation to compute time derivatives at each timestamp in 3 , since all monotone splines have analytical derivatives. While automatic differentiation is necessary for MNN-based methods to recover the intensity from the CCIF, MAS avoids this derivative-computation step and further improves computational efficiency. In Appendix B, we provide a derivation of the analytical derivative for RQS.

5.4 MAS for Multivariate TPPs↩︎

MAS can be extended to multivariate TPPs. If the multivariate TPPs consist of \(K\) different event types, we can use \(K\) MAS models to represent \(K\) CCIFs. Given the event history \(\{(t_1, k_1), \dots, (t_n, k_n)\}\), a Transformer (parameterized by \(\theta_1\)) is employed to encode the historical information into a history embedding \(\mathbf{h}_n\). Then, \(\mathbf{h}_n\) is mapped by \(K\) different MLPs followed by Softplus (parameterized by \(\theta_{2,1}, \dots, \theta_{2,K}\)) to \(K\) sets of distinct variables for \(K\) MAS models. Consequently, we obtain \(\Lambda_1^{*}(t), \dots, \Lambda_K^{*}(t)\). The overall model parameters are given by \(\theta = \{\theta_1, \theta_{2,1}, \dots, \theta_{2,K}\}\). Subsequently, we train the model using 3 and utilize 2 to predict the next timestamp and event type.

6 Generalization Analysis↩︎

This section explores the theoretical generalization guarantees of Monotone Alternating Splines (MAS). We systematically outline our analysis in three aspects: first, we establish the generalization bound of MAS; second, we formally prove its theoretical superiority over conventional MNNs in overcoming inherent approximation barriers; and finally, we provide principled guidelines for MAS hyperparameter selection based on the established bounds.

6.1 Generalization Bound of MAS↩︎

To assess the model’s reliability on unseen data, we evaluate the gap between the empirical negative log-likelihood and its expectation. Suppose the training set consists of \(Z\) i.i.d. point process sequences \(S_1, \dots, S_Z\) on the time window \([0,T]\). The \(n\)-th timestamp of the \(z\)-th sequence is denoted as \(t_{zn}\). For simplicity, we consider the univariate case, where the max sequence length is \(N_0\). Without loss of generality, we assume the interpolation interval is a constant \(\Delta>0\). The empirical loss is defined as the averaged negative log-likelihood over all sequences: \(\mathcal{\widehat{L}}:=-\frac{1}{Z}\sum_{z=1}^{Z}\frac{1}{T}\log p(S_z)\). We train the model by minimizing \(\mathcal{\widehat{L}}\) and are interested in its expected value over the data distribution, denoted as \(\mathcal{{L}}:=\mathbb{E}[\mathcal{\widehat{L}}]\). We can first obtain the following lemma.

Then, we can derive the error bound between \(\mathcal{L}\) and \(\widehat{\mathcal{L}}\).

5 offers the bound between the empirical negative log-likelihood and its expectation, implying minimizing the empirical loss \(\widehat{\mathcal{L}}\) helps reduce the expected loss \(\mathbb{E}[\widehat{\mathcal{L}}]\). Taking a step forward, we combine 4 and 5 to derive the upper bound for MAS’s generalized error.

In Theorem 6, excluding the probability error, the core estimation errors are decomposed into three interpretable parts. First, the interpolation error \(C_0\left(\frac{\Delta}{C_4}+\frac{\Delta^2}{2}\right)\) decreases as the interpolation interval \(\Delta\) shrinks, as a smaller \(\Delta\) leads to a more accurate fit to the ground-truth CCIF. Second, the extrapolation error \(\frac{B}{T\cdot L}\) decreases as the total interpolation length \(L\) increases. As the contribution of the interpolation part becomes larger, the influence of the extrapolation part diminishes. Finally, the complexity error \(\frac{C_2}{\sqrt{Z}}\frac{\sqrt{L}}{\sqrt{\Delta}}R(L,\Delta)\) increases as the number of parameters in MAS increases, implying higher model complexity.

6.2 Theoretical Superiority over MNNs↩︎

As highlighted in Section 4, due to the inevitable convexity restriction, convex MNNs fundamentally fail to capture TPPs with decaying CIFs. We formally quantify this limitation as an irreducible error bound.

7 reveals that, for any sufficiently fine-grained MAS configuration satisfying \(\Delta < \frac{\gamma}{2 C_0}\), MAS achieves a strictly smaller approximation error than MNNs. The error gap mathematically proves that MAS breaks the inherent approximation deadlock of continuous MNNs.

6.3 Hyperparameter Selection for MAS↩︎

Beyond theoretical comparisons, 6 also offers guidance on selecting hyperparameters. Given a fixed sample size \(Z\), two important hyperparameters, \(\Delta\) and \(L\), can be chosen by minimizing the bound in ?? , as is presented in the following corollary.

According to Corollary 3, hyperparameters \(L\) and \(\Delta\) depend on two dataset-specific quantities: (1) the expected inter-event time \(\mathbb{E}[t_i]\), and (2) the degree of temporal variability in the TPP intensity, captured by its Lipschitz constant \(l_1\). As \(\mathbb{E}[t_i]\) grows, a longer interpolation horizon, i.e., a larger \(L\), is required to track changes in the CCIF accurately. Moreover, as the intensity fluctuation, reflected by \(l_1\), becomes sharper, MAS should provide a finer-grained approximation, which requires a smaller \(\Delta\). A larger sample size \(Z\) also supports a finer interpolation resolution, while a larger per-knot freedom \(q\) or longer maximum sequence length \(N_0\) requires a coarser interval to control complexity. Furthermore, Corollary 3 provides a simple plug-in strategy for selecting \(L\) and \(\Delta\): For synthetic datasets where the functional form of the intensity is known, \(l_1\), \(\lambda_{\min}\), and \(\mathbb{E}[t_i]\) can be estimated directly from samples, yielding an approximate lower bound for \(L\) and an approximate upper bound for \(\Delta\). For general real-world datasets, inspired by Scott’s rule [40], we suggest assuming that the data follow a classic TPP, e.g., exponential-kernel Hawkes process, and estimating \(L\) and \(\Delta\) accordingly. Given \(L\) and \(\Delta\), the number of interpolation knots \(p\) can also be calculated as \(p=L/\Delta\).

It is worth noting that, although 6 shows that an excessively large \(L\) and an overly small \(\Delta\) may lead to overfitting and thus increased generalization error, such overfitting can be effectively mitigated by standard training techniques (e.g., dropout [41] or SGD-based optimization [42]). In practice, we observe that this overfitting effect is mild. Therefore, we recommend using a moderately large \(L\) and a relatively fine \(\Delta\), consistent with the bounds provided in Corollary 3.

7 Experiments↩︎

This section presents the performance of MAS compared to other baseline models. We also evaluate the efficiency of MAS relative to other CCIF- or CIF-based methods. An ablation study is conducted at the end of the experiment.

7.1 Baselines, Evaluation Metrics and Experiment Setups↩︎

Baselines: For experiment baselines, we systematically categorize thirteen representative models into four paradigms:

• Parametric Models: Inhomogeneous Poisson Process (IPP) [18] and the standard Hawkes process [20].

• CIF Models: Approaches that directly parameterize the CIF, including RMTPP [11], THP [13], and SAHP [22].

• MNN-based CCIF Models: Approaches that leverage continuous MNNs to model the CCIF, like FullyNN [15], EMTPP [17], and CuFun [32].

• Integral/Intensity-Free and other state-of-the-art Models: Alternative flow-based or score-matching-based models, including IFTPP [30], TriTPP [27] and WSM [29]. We also compare other state-of-the-art models like HYPRO [43] and DTPP [31].

Some baselines, such as IPP, Hawkes, FullyNN, and TriTPP, only consider the univariate case in their original implementations. Therefore, our experiments are divided into two parts. First, we compare the performance of MAS with all baselines on univariate synthetic data. Then, we evaluate MAS on multivariate real-world datasets, comparing it with multivariate baselines, including RMTPP, THP, SAHP, EMTPP, IFTPP, WSM, CuFun, HYPRO, and DTPP.

Metrics: We use three evaluation metrics: negative log-likelihood (NLL), root-mean-square error (RMSE) of predicted timestamps, and accuracy (ACC) of predicted event types. The log-likelihood, predicted timestamps, and event types are computed based on 1 2 .

Experimental Setup: To ensure the performance gains solely originate from our proposed spline-based CCIF formulation rather than the capacity of the history encoder, we rigorously align the encoder architecture. Specifically, MAS, FullyNN, THP, EMTPP, and WSM are all uniformly equipped with an identically configured Transformer network (1 layer, 16 attention heads, 64-dimensional model size) to extract the history embedding. Furthermore, since continuous MNNs (e.g., FullyNN, EMTPP) are highly sensitive to large numerical inputs and prone to encountering convergence failures during optimization, we applied a uniform timestamp scaling strategy (e.g., dividing inter-event times by 10 or 100 depending on the absolute scale of the datasets). All baselines were retrained and evaluated under the exact same scaled configurations, identical learning rate (\(0.001\)), fixed batch size (\(64\)), and identical \(60\% - 20\% - 20\%\) train-validation-test splits. For MAS, we incorporate 10 knots for the interpolation component, using RQS as the interpolation monotone spline and a linear-plus-exponential function \(g(t)= c + bt - \exp(-at)\) as the extrapolation function. The interpolation support \(L\) is set to a constant. All models are trained with the Adam optimizer, conducted on an RTX 4090 with 24GB of memory. Details refer to Appendix C.

7.2 Synthetic Univariate Data↩︎

Datasets: In the synthetic univariate experiment, we consider five datasets constructed by FullyNN [15]: two Hawkes processes (1,2) with different intensities, two renewal processes 1,2 (stationary/non-stationary), and a self-correcting process. Details refer to Appendix C.

Results: ¿tbl:tab:hawkes95results? presents the performance comparison between all baselines and MAS. MAS achieved the best performance across most metrics on the majority of datasets. This validates our theory that, by introducing piecewise monotone spline as the interpolation component, MAS is capable of fitting various ground-truth CCIFs derived from different classical TPP. In contrast, other traditional models (e.g., IPP) or deep Hawkes processes (e.g., THP) become less effective in characterizing other TPP. Meanwhile, MAS also outperforms models that parametrize CCIF using MNNs in terms of fitting ability. 4 shows MAS’s predicted CIF and CCIF for two datasets, implying MAS successfully captures the intensity patterns of different TPP. A detailed comparison between MAS and FullyNN is provided in Appendix C, further validating the superiority of MAS over classic CCIF-based methods.

7.3 Real Multivariate Data↩︎

Datasets: We consider four real multivariate datasets: Taxi [44], Taobao [45], Retweet [46] and Earthquake [45]. Detailed information is in Appendix C.

Results: ¿tbl:tab:performance95comparison? presents the performance comparison between all baselines and MAS. On real-world datasets, MAS once again achieved the best results across multiple metrics. This further demonstrates the advantages of MAS in terms of flexibility and generalization, indicating that in multivariate scenarios, MAS successfully captures the interactions among multiple variables. The success of MAS can be attributed to two factors: first, the inclusion of more interpolation knots along the timeline, which makes MAS more flexible than baseline models; second, the introduction of the extrapolation component, which ensures that MAS correctly simulates the CCIF while maintaining strong generalization ability. Therefore, our model achieves greater flexibility compared to baseline models without compromising generalization.

7.4 Efficiency Test↩︎

We further compare the computational efficiency of MAS against THP, RMTPP, and two CCIF-based models, EMTPP and FullyNN. To ensure a fair comparison, we roughly align the parameter counts of all models on the Hawkes1 dataset and record the actual training time over 50 epochs under identical settings. As shown in 4 (e), MAS runs significantly faster than all four baselines, demonstrating that it is highly competitive in computational efficiency alongside its modeling accuracy. The details are provided in Appendix C. For EMTPP and FullyNN, automatic differentiation invariably increases the per-step computational overhead and exacerbates backpropagation costs over long sequences and high-dimensional parameter spaces. In contrast, the explicitly decoupled piecewise interpolation and extrapolation components of MAS enable direct intensity evaluations within each time segment, eliminating the need for gradient flows or numerical approximations. Ultimately, MAS delivers significantly lower training times and superior computational efficacy.

7.5 Long Horizon Prediction↩︎

We further evaluate the long-horizon prediction capabilities of MAS on the univariate Hawkes1 and multivariate TAXI datasets, comparing its performance against THP, EMTPP, and CuFun in forecasting the next 3, 5, and 10 events. As shown in ¿tbl:tab:long95horizon?, MAS not only maintains high accuracy in short-term predictions but also consistently outperforms the baselines as the forecasting horizon extends. The sustained superiority of MAS in multi-step scenarios indicates that it does not merely overfit localized intensity fluctuations immediately following an event; rather, it effectively leverages the extrapolation component to secure stable long-term trends, thereby significantly mitigating cumulative error propagation. For sequences exhibiting strong self-exciting dynamics, such as Hawkes1, MAS precisely captures both the rapid intensity decay and subsequent re-excitations. Meanwhile, on multivariate datasets like TAXI, MAS seamlessly integrates shared history representations with highly expressive piecewise monotone splines, effectively modeling complex mutual interactions across different event types and consequently enhancing the joint forecasting accuracy for multiple future events.

Long-horizon results. The best is in bold.

7.6 Ablation Study↩︎

Parameter Sensitivity: We conduct ablation studies to verify the contribution of each component in MAS. We consider the following factors: (1) the type of interpolation monotone spline (RQS/RLS), (2) the type of extrapolation function (linear/linear-plus-exponential), (3) the interpolation support length (\(L=4/6/10\)), (4) the number of interpolation knots (\(p=5/10/50/100\)). The experiments are performed on the Taxi and Earthquake datasets, using NLL and ACC as metrics. We set RQS interpolation, linear-plus-exponential extrapolation, \(L = 6\), and \(p = 10\) as the default choices.

As shown in ¿tbl:tab:ablation95study?, as the number of interpolation knots \(p\) increases and the interpolation support length \(L\) extends, the model’s performance on both datasets first improves and then declines. This validates our theoretical analysis: the increase in interpolation support length and the number of interpolation knots leads to a rise in the complexity error term in 6. This overfitting, however, is not severe. Furthermore, MAS is insensitive to different interpolation and extrapolation functions. We speculate that under the default settings (L=6, p=10), MAS already fits the CCIFs at sufficiently fine granularity, making different functions yield negligible differences.

8 Conclusion↩︎

We propose a novel approach, MAS, for modeling the CCIF of TPPs. Compared to MNNs, MAS offers greater generality, flexibility, and efficiency, and enjoys strong theoretical guarantees for both fitting and generalization. We prove that the interpolation component of MAS enhances fitting performance, while the extrapolation component of MAS enhances generalization performance. Extensive experiments demonstrate that MAS outperforms other CCIF- and CIF-based modeling methods on most univariate and multivariate benchmarks across multiple evaluation metrics.

9 Notations↩︎

All notations used in the paper are listed in the 1.

10 Monotone Splines↩︎

This appendix section introduces some common monotone splines and discusses their parameterizations.

10.1 Rational Quadratic Splines↩︎

A piece of RQS on the interval \([w_{m-1},w_m]\) can be written as: \[f_{m}(t)= y_{m-1} + \frac{\left(y_{m} - y_{m-1}\right) \left(s_m \tau^2 + \delta_{m-1} \tau (1 - \tau)\right)}{s_{m} + \left(\delta_{m} + \delta_{m-1} - 2s_{m}\right) \tau (1 - \tau)}, \label{RQS95another32form}\tag{10}\] where \(s_m=\frac{y_{m}-y_{m-1}}{w_m-w_{m-1}}\) and \(\tau=\frac{t-w_{m-1}}{w_m-w_{m-1}}\in[0,1]\). When \(t=w_{m-1}\), i.e., \(\tau=0\), \(f_{m}(t)=y_{m-1}\), \(f^{\prime}_{m}(t)=\delta_{m-1}\). When \(t=w_{m}\), i.e., \(\tau=1\), \(f_{m}(t)=y_{m}\), \(f^{\prime}_{m}(t)=\delta_{m}\). This implies that once \(w_{m-1}, w_{m}, \delta_{m-1},\delta_{m},y_{m-1},y_{m}\) are specified, the RQS piece is uniquely determined.

The derivative with respect to \(\tau\) can be written as: \[\begin{align} f_m'(\tau) &= (y_m - y_{m-1}) \frac{ \bigl[2s_m\tau + \delta_{m-1}(1-2\tau)\bigr] \; }{ J_m }\\ &-(y_m - y_{m-1})\frac{\; \bigl[s_m\tau^2 + \delta_{m-1}\,\tau(1-\tau)\bigr] V_m(1-2\tau)}{J_m^2}, \end{align}\] where \(V_m=\delta_m+\delta_{m-1}-2s_m\), \(J_m=\bigl[s_m + V_m\,\tau(1-\tau)\bigr]\). The time derivative is obtained by multiplying the expression above by \(\frac{d\tau}{dt}=\frac{1}{w_m-w_{m-1}}\).

10.2 Rational Cubic Splines↩︎

A piece of RCS on the interval \([w_m,w_{m+1}]\) can be written as: \[\begin{align} f_m(t) &= \frac{(1 - \tau)^3 v_m y_m + \tau(1 - \tau)^2 \left[ (2u_m v_m + v_m) y_m + v_m h_m d_m \right] }{(1 - \tau)^2 v_m + 2 u_m v_m \tau(1 - \tau) + \tau^2 u_m}\notag\\ &+\frac{\tau^2(1 - \tau)\left[ (2u_m v_m + u_m) y_{m+1} - u_m h_m \delta_{m+1} \right] + \tau^3 u_m y_{m+1}}{(1 - \tau)^2 v_m + 2 u_m v_m \tau(1 - \tau) + \tau^2 u_m}. \end{align}\] where \(\tau=\frac{t-w_m}{w_{m+1}-w_m}\in[0,1]\) and \(h_m=w_{m+1}-w_m\), \(f_{m}(w_m)=y_m\), \(f_{m}(w_{m+1})=y_{m+1}\), \(f'_{m}(w_m)=\delta_m\), \(f'_{m}(w_{m+1})=\delta_{m+1}\). \(u_m\) and \(v_m\) are two shape parameters, \(u_m,v_m>0\).

10.3 Monotone Cubic Splines↩︎

On \([w_m, w_{m+1}]\), let \(h_m = w_{m+1}-w_m,\;\tau = \frac{ t- w_m}{h_m}\), then a piece of monotone cubic spline can be expressed as: \[\begin{align} f_m(t) &= y_m\,h_{00}(\tau) + h_m\,\delta_m\,h_{10}(\tau) \\ &+ y_{m+1}\,h_{01}(\tau) + h_m\,\delta_{m+1}\,h_{11}(\tau), \end{align}\] where \[\begin{align} h_{00}(\tau) = 2\tau^3 - 3\tau^2 + 1,\;& h_{10}(\tau) = \tau^3 - 2\tau^2 + \tau,\\ h_{01}(\tau) = -2\tau^3 + 3\tau^2,\;& h_{11}(\tau) = \tau^3 - \tau^2. \end{align}\]

Again, \(f_{m}(w_m)=y_m\), \(f_{m}(w_{m+1})=y_{m+1}\), \(f'_{m}(w_m)=\delta_m\), \(f'_{m}(w_{m+1})=\delta_{m+1}\).

10.4 Rational Linear Splines↩︎

On the interval \([w_{m-1},w_{m}]\), a piece of RLS can be written as:

\[\label{eq:2461} f(t) \;=\; a_m \;+\; \frac{c_m\,(t - w_m)}{1 + d_m\,(t - w_m)}, \quad t \in [w_{m-1},\,w_m].\tag{11}\]

\(f(t)\) is monotone increasing as long as \(c_m>0\) and the denominator remains positive on the interval.

11 Experiments↩︎

11.1 Datasets↩︎

We consider the following five univariate datasets:

(1) Hawkes Process (1, 2): We consider two exponential-decay Hawkes processes with CIF: \(\lambda^*(t)=0.2+\sum_{t_n<t}\sum_{i=1}^M\alpha_i\exp(-\beta_i(t-t_n))\). We set \(M=1,\alpha_1=0.8,\beta_1=1\) for the first dataset and \(M=2,\alpha_1=0.4,\beta_1=\beta_2=1,\alpha_2=0.4\) and \(\beta_2=20\) for the second dataset.

(2) Self-correcting Process: The CIF is given as \(\lambda^*(t)=\exp(t-\sum_{t_n<t}1)\).

(3) Renewal Process (Stationary / Non-stationary): For the stationary case, timestamp intervals are independent, log-normal distributed with a mean \(1.0\) and standard deviation \(0.5\). For the non-stationary case, the trend function is set to \(r(t)=0.99\sin(2\pi t/20000)+1\).

For multivariate datasets, their information is listed as follows:

(1) Taxi [44]: This dataset tracks the time-stamped taxi pick-up and drop-off events across the New York City. The event types are categorized into \(K=10\) types and there are \(2000\) sequences with an average sequence length of \(40\).

(2) Taobao [45]: This dataset contains time-stamped user click behaviors on Taobao shopping pages. The event types are categorized into \(K=20\) types and there are \(4800\) sequences with an average sequence length of \(150\).

(3) Retweet [46]: This dataset contains time-stamped user retweet event sequences. The events are categorized into \(K = 3\) types. There are \(9000\) sequences with an average sequence length of \(90\).

(4) Earthquake [45]: This dataset contains time-stamped earthquake event sequences. The events are categorized into \(K = 7\) types. There are \(3000\) sequences with an average sequence length of \(20\).

11.2 Experimental Setup↩︎

Because MNN is sensitive to the scale of its outputs, we applied a timestamp transformation to part of our data. Following the approach of [29], we divided every timestamp in the retweet dataset by 100 and divided every timestamp in the earthquake dataset by 5. In addition, we divided the timestamps of each of the five synthetic datasets by 10. This ensures that all models converge within 100 epochs.

All models were trained on a single NVIDIA RTX\(\cdot\)​4090\(\cdot\)D (24 GB). For fairness, MAS, FullyNN, EMTPP, THP, and WSM all encode historical information using Transformers of the same size. The Transformer hyperparameters are as follows: number of attention heads = 16, number of layers = 1, model dimension = 64, inner feed‑forward dimension = 8, key dimension = 16, value dimension = 16, and dropout rate = 0.1. For MAS, when using a Rational Quadratic Spline (RQS) as the interpolation function, we set the lower bound for both the interpolation interval width and length to 0.01, and the minimum derivative at interpolation points to 0.01.

We fixed the batch size at 64 and the learning rate at 0.001 for all models. Each model was trained for 100 epochs on the train set and converged within 100 epochs. We then evaluated their performance on the test set.

When comparing the accuracy of MAS with that of other models, we observed that—even for the same model on the same dataset—the reported performance can vary dramatically across different studies. For common baseline models such as SAHP, THP, and RMTPP, we referenced their best-reported accuracy from the literature (for example, [25], [45], and [29]) and our experiments, and compared those to our MAS. For instance, on the Taobao dataset, THP achieved an accuracy of 0.467 according to [25], 0.531 according to [45], and 0.594 in [29]. all below MAS’s performance of 0.600.

For datasets like Taobao, where the raw time intervals are very large, we first scaled them down by a factor of 1/100 to ensure rapid convergence of models such as EMTPP. However, this scaling can cause the model’s NLL to change substantially, making it difficult to compare against past work (since they may have used different scaling factors). Therefore, we retrained all baseline models on the same scaled dataset—keeping the data, training and testing strategy, and number of epochs identical—and then reported their NLL on the same test set.

11.3 Supplementary Results↩︎

11.3.1 Cross Validation↩︎

To rigorously demonstrate that MAS indeed outperforms other baselines, we conduct a cross-validation study comparing MAS with a classical CCIF-based method, FullyNN. We evaluate FullyNN under different network depths, widths, activation functions, and training data sizes to ensure that MAS’s superiority is not merely due to accidental hyperparameter choices.

Specifically, we fix the encoder architecture of MAS and FullyNN to be identical (i.e., both models receive the same input features). Following the data generation procedure of the hawkes1 dataset, we independently sample 400 and 800 sequences within the time range \([0,100]\) for training. During training, we vary FullyNN’s number of fully connected layers from 1 to 3, its hidden widths between 16 and 32, and its activation functions between softplus and sigmoid (representing different levels of smoothness). FullyNN is trained for 100 epochs until convergence. We then compare its NLL and RMSE with those of MAS, as shown in ¿tbl:tab:32CrossVal?.

The results show that when the network depth exceeds two layers, the performance of FullyNN no longer improves significantly with additional depth. Moreover, regardless of whether softplus or sigmoid is used, MAS consistently outperforms FullyNN across all configurations. This cross-validation experiment further confirms the advantage of MAS over classical CCIF-based modeling methods.

11.3.2 Training Time↩︎

We align the parameter sizes of five models—THP, RMTPP, MAS, EMTPP, and FullyNN—and train them for 100 epochs on the Hawkes1 dataset. Experiments are rerun three times under the same conditions. We then measure the mean and the standard deviation of their running time (in minutes), and the results are summarized in 2. A visualized version of this comparison is shown in 4 (e).

12 Proof↩︎

The following mild assumptions are required to prove the theorems in the paper. Compared with the original version, we only make explicit the regularity conditions that were implicitly used in the derivations.

12.1 Proof of Theorem 1↩︎

12.2 Proof of Theorem 2↩︎

12.3 Proof of Theorem 3↩︎

12.4 Proof of Theorem 4↩︎

12.5 Proof of Theorem 5↩︎

When proving the generalization bound of the MAS model, we further assume that the horizontal interval of MAS interpolation is a fixed value \(\Delta\). That is, after each timestamp \(t_n\), there are \(M\) interpolation intervals, with a total length of \(L = M\Delta\). Equivalently, the interpolation support contains \(M+1\) knots; this off-by-one convention only changes constants. For the next timestamp \(t_{n+1}\), there is a probability that it falls into one of the \(M\) interpolation intervals, or into the extrapolation interval \([t_n + M\Delta, \infty)\). We denote the index of the interval into which \(t_{n+1}\) falls as \(r_{n+1}\), taking values in \(1, 2, \dots, M\).

We first prove 1:

Before proving 5, we introduce the following concepts.

Bracketing numbers quantify how many simple function pairs (brackets) are needed to approximate a complex function class within a given accuracy. The smaller the bracketing number, the simpler or more regular the function class is. Utilizing the following lemma in [47], we can acquire the bracketing number bound for MAS:

Utilizing 2 and 1, the bracketing number of MAS can be upper bounded:

An important feature of the bracketing number is that it bounds the complexity of a model. We first introduce the concept of Rademacher complexity [48].

Rademacher complexity can be understood as the degree to which a model fits noise. The higher the Rademacher complexity of a model, the stronger its ability to fit noise, and the more prone it is to overfitting. Since 3 bounds the negative log-likelihood class, we may normalize the loss functions into \([0,1]\); the scaling constant is absorbed into \(C_2\) and \(C_3\). We can derive the following lemma:

There exists the following relationship between the Rademacher complexity and the bracketing number of the model:

Now that we get: \[\begin{align} \vert \widehat{\mathcal{L}}-\mathbb{E}[\widehat{\mathcal{L}}]\vert &\le 2 \mathcal{R}_Z(\mathcal{F}) + \sqrt{\frac{\ln(1/\xi)}{2Z}} \\ &\le \frac{12}{\sqrt{Z}} \int_0^{\operatorname{diam}\Theta} \sqrt{\log N_{[\,]}(\varepsilon, \mathcal{F})} \, \text{d}\varepsilon+\sqrt{\frac{\ln(1/\xi)}{2Z}}\\ &\le \frac{12}{\sqrt{Z}} \int_0^{\operatorname{diam}\Theta} \sqrt{\log \left[ K \frac{B}{T}W^{\frac{1}{2}d} \left( \frac{ \operatorname{diam} \Theta}{\varepsilon} \right)^d\right]} \, \text{d}\varepsilon + \sqrt{\frac{\ln(1/\xi)}{2Z}}\\ &\le \frac{12}{\sqrt{Z}}\int_0^{\operatorname{diam}\Theta}{\sqrt{d}\sqrt{C_3 -\log \epsilon}}\,\text{d}\varepsilon +\sqrt{\frac{\ln \left( 1/\xi \right)}{2Z}}, \end{align}\] where \(W=\mathbb{E}\left(\sum_{n=1}^{N} r_n\right)^2\).

Note that the parameter dimension is \(N_0Mp\) (\(N_0\) is the maximum number of timestamps), \(M\) is the number of interpolation knots, and \(p\) is the number of parameters required for each interpolation interval. In 5, \(p\) is treated as a fixed per-knot degree of freedom and absorbed into \(C_2\); it is restored explicitly in the hyperparameter-selection corollary. By collecting each term in the above inequality, we acquire the final bound: \[\begin{align} \vert \widehat{\mathcal{L}}-\mathbb{E}[\widehat{\mathcal{L}}]\vert \le\frac{1}{\sqrt{Z}}\left(\frac{1}{2}\sqrt{{\log\frac{1}{\xi}}}+{C_2\sqrt{N_0M}}\int_0^{c}\sqrt{C_3-\log t}d t\right), \end{align}\] where \(c=\operatorname{diam}\Theta\). This finishes the proof.

12.6 Proof of Theorem 6↩︎

Suppose \(\mathcal{L}^*\) is the optimal negative log-likelihood, derived from the ground-truth intensity. Then, the generalized loss for MAS can be decomposed as: \[\begin{align} \mathbb{E}\widehat{\mathcal{L}}&\le \mathcal{L}^*+\left| \mathcal{L}^*-\widehat{\mathcal{L}} \right|+\left| \widehat{\mathcal{L}}-\mathbb{E}\widehat{\mathcal{L}} \right| \\ &\le \mathcal{L}^*+\frac{1}{ZT}\sum_{z=1}^Z{\left| \log p^*\left( S_z \right) -\log p_{\theta}\left( S_z \right) \right|}+\left| \widehat{\mathcal{L}}-\mathbb{E}\widehat{\mathcal{L}} \right|\\ &\le \mathcal{L}^* + \frac{1}{ZT}\sum_{z=1}^Z{\sum_{n=1}^N{\left| \log p^*\left( t^{\left( z \right)}_n \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|}}+\left| \widehat{\mathcal{L}}-\mathbb{E}\widehat{\mathcal{L}} \right|\\ &\le \underbrace{\frac{1}{ZT}\sum_{z,n=1}^{Z,N}{\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|}\mathcal{I}\left( t_{n}^{\left( z \right)}\ge t_{n-1}^{\left( z \right)}+M\Delta \right)}_{W_1} \\ &+\underbrace{\frac{1}{ZT}\sum_{z,n=1}^{Z,N}{\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|}\mathcal{I}\left( t_{n}^{\left( z \right)}<t_{n-1}^{\left( z \right)}+M\Delta \right)}_{W_2} +\mathcal{L}^* + \underbrace{\left| \widehat{\mathcal{L}}-\mathbb{E}\widehat{\mathcal{L}} \right|}_{W_3}. \end{align}\] \(W_1\) can be bounded by following inequality: \[\begin{align} \quad \;\;W_1&=\frac{1}{ZT}\sum_{z,n=1}^{Z,N}{\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|\mathcal{I}_n} \\ &\le \frac{B_0}{T}\mathcal{I}\left( t_{n}^{\left( z \right)}-t_{n-1}^{\left( z \right)}\ge M\Delta \right), \end{align}\] where \(\mathcal{I}_n=\mathcal{I}\left( t_{n}^{\left( z \right)}\ge t_{n-1}^{\left( z \right)}+M\Delta \right)\), and \(B_0\) is a bound on \(\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|\) from 3. The factor depending on the maximum sequence length is absorbed into the constant \(B\). Thus, we can take the expectation on both sides: \[\begin{align} &\quad \;\mathbb{E}W_1\le \frac{B_0}{T}\mathbb{P}\left( t_{n}^{\left( z \right)}-t_{n-1}^{\left( z \right)}\ge M\Delta \right) \\ &\le \frac{B_0}{T}\int_{t_{n}^{\left( z \right)}-t_{n-1}^{\left( z \right)}\ge M\Delta}{\frac{t_{n}^{\left( z \right)}-t_{n-1}^{\left( z \right)}}{M\Delta}}\text{d}\mu \le \frac{B_0}{T}\frac{\mathbb{E}t_i}{L}:=\frac{B}{TL}. \end{align}\] For \(W_2\), we have proven in 4 that there exists a MAS \(\Lambda^*_\theta\), s.t., \(\Vert\Lambda^*_\theta-\Lambda^*\Vert\le\frac{1}{2}C_0\Delta^2\), \(\Vert\Lambda^{*\prime}_\theta-\Lambda^{*\prime}\Vert\le C_0\Delta\). As a result, \(W_2\) can be bounded by: \[\begin{align} \;\;W_2&=\frac{1}{ZT}\sum_{z=1}^Z{\sum_{n=1}^N{\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|}}(1-\mathcal{I}_n) \\ &\le \frac{1}{ZT}\sum_{z=1}^Z{\sum_{n=1}^N{\left( \left| \log \Lambda _{\theta}^{*'}\left( t_{n}^{\left( z \right)} \right) -\log \Lambda ^{*'}\left( t_{n}^{\left( z \right)} \right) \right|\right.}}\\ &+\left.\left| \left( \Lambda _{\theta}^{*}\left( t_{n}^{\left( z \right)} \right) -\Lambda _{\theta}^{*}\left( t_{n-1}^{\left( z \right)} \right) \right) -\left( \Lambda ^{*}\left( t_{n}^{\left( z \right)} \right) -\Lambda ^{*}\left( t_{n-1}^{\left( z \right)} \right) \right) \right| \right)\\ &\le \frac{1}{ZT}\sum_{z=1}^Z{\sum_{n=1}^N{\left( \frac{1}{v_1}\lVert \Lambda _{\theta}^{*'}\left( t_{n}^{\left( z \right)} \right) -\Lambda ^{*'}\left( t_{n}^{\left( z \right)} \right) \rVert +C_0\Delta ^2 \right)}}\\ &\le \frac{1}{ZT}\sum_{z=1}^Z{\sum_{n=1}^N{\left( \frac{C_0}{v_1}\Delta +C_0\Delta ^2 \right)}} \le\frac{1}{T}{\sum_{n=1}^{N_0}{\left( \frac{C_0}{v_1}\Delta +C_0\Delta ^2 \right)}}:={C_0\left(\frac{\Delta}{C_4}+\frac{\Delta^2}{2}\right)}. \end{align}\] Here \(C_4\) absorbs the lower intensity bound \(v_1\) and the normalization constants involving \(T\) and \(N_0\). Finally, according to 5: \[\begin{align} W_3&=\vert \widehat{\mathcal{L}}-\mathbb{E}[\widehat{\mathcal{L}}]\vert \le \frac{1}{\sqrt{Z}}\left(\frac{1}{2}\sqrt{{\log\frac{1}{\xi}}}\right)+ \frac{1}{\sqrt{Z}}\left({C_2\sqrt{N_0M}}\int_0^{c}\sqrt{C_3-\log t}d t\right)\\ :&=\frac{1}{\sqrt{Z}}\left(\frac{1}{2}\sqrt{{\log\frac{1}{\xi}}}\right)+\frac{C_2}{\sqrt{Z}}\frac{\sqrt{L}}{\sqrt \Delta}R(L,\Delta). \end{align}\] By combining \(W_1,W_2\) and \(W_3\), we acquire: \[\label{equ:32bound3240app41} \mathbb{E}[\widehat{\mathcal{L}}]\le {\mathcal{L}^*}+ {\frac{1}{2\sqrt{Z}}\sqrt{{\log\frac{1}{\xi}}}}+{\frac{B}{T\cdot L}} +{C_0\left(\frac{\Delta}{C_4}+\frac{\Delta^2}{2}\right)}+ {\frac{C_2}{\sqrt{Z}}\frac{\sqrt{L}}{\sqrt{\Delta}}R(L,\Delta)}.\tag{13}\] This finishes the proof.

12.7 Proof of Theorem 7↩︎

12.8 Proof of Corollary 3↩︎

Finally, we estimate the bounds for two important parameters, \(\Delta\) and \(L\), by minimizing the generalization bound in 6. If we restore all constants in the bound 13 , we get:

\[\label{equ:32bound32detail} \mathbb{E}[\widehat{\mathcal{L}}]\le {\mathcal{L}^*} + {\frac{1}{2\sqrt{Z}}\sqrt{{\log\frac{1}{\xi}}}} + {\frac{N_0B_0}{p}\cdot\frac{\mathbb{E} t_i}{T\cdot L}} + {l_1\left(\frac{L}{pv_1}+\frac{L^2}{p^2}\right)} + {\frac{12}{\sqrt{Z}}{\sqrt{N_0pq}}\int_0^{\text{diam}\Theta}\sqrt{C_3-\log t}dt},\tag{14}\]

where \(\Delta=L/p\). We first focus on how to estimate \(L\). Typically, \(\Delta\ll1\), thus we could ignore the \(\Delta^2\) term in 13 and focus on the term \(\frac{l_1L}{v_1p}\). Notice that in 14 , only two terms \({N_0B_0\cdot\frac{\mathbb{E} t_i}{T\cdot L}}+\frac{l_1L}{v_1p}\) relate to \(L\). Then when 14 is minimized, \(L\) should satisfy: \[L\ge\sqrt{\frac{v_1N_0B_0\mathbb{E}[t_i]}{l_1T}}.\] As a bound between the ground truth \(\log p^*\left( t_{n}^{\left( z \right)} \right)\) and \(\log p_{\theta}\left( t_{n}^{\left( z \right)} \right)\), \(B_0=\sup_{t^{(z)}_n}\left| \log p^*\left( t_{n}^{\left( z \right)} \right) -\log p_{\theta}\left( t_{n}^{\left( z \right)} \right) \right|>1\) typically holds. Moreover, \(N_0\) is the maximum number of events that happen during the interval \([0,T]\). Thus, \(N_0B_0\ge T\), and the bound can be reduced to: \[L\ge\sqrt{\frac{v_1\mathbb{E}[t_i]}{l_1}},\] where \(l_1\) is the Lipschitz constant of the intensity during event intervals, and \(v_1\) is the base intensity during \([0,T]\).

We next focus on deriving the bound of \(p\). Because \(\Delta=L/p\), 14 can be rewritten as:

\[\begin{align} \label{equ:32bound32detail322} \mathbb{E}[\widehat{\mathcal{L}}] & \le {\mathcal{L}^*} + {\frac{1}{2\sqrt{Z}}\sqrt{{\log\frac{1}{\xi}}}} + {\frac{N_0B_0}{p}\cdot\frac{\mathbb{E} t_i}{T\cdot \Delta}} + {l_1\left(\frac{L}{pv_1}+\frac{L^2}{p^2}\right)} \notag \\ &+ {\frac{12}{\sqrt{Z}}{\sqrt{N_0pq}}\int_0^{\text{diam}\Theta}\sqrt{C_3-\log t}dt}, \notag \\ &\le {\mathcal{L}^*} + {\frac{1}{2\sqrt{Z}}\sqrt{{\log\frac{1}{\xi}}}} + {\frac{2N_0B_0}{p}\cdot\frac{\mathbb{E} t_i}{T\cdot \Delta}} + {\frac{12}{\sqrt{Z}}{\sqrt{N_0pq}}\int_0^{\text{diam}\Theta}\sqrt{C_3-\log t}dt}. \end{align}\tag{15}\]

The second step is due to the assumption that \(l_1\Delta^2\) is significantly smaller than \(\Delta\) and can be ignored. Let \(C_5=\int_0^{\text{diam}\Theta}\sqrt{C_3-\log t}dt\). For a fixed interpolation length \(L\), the leading \(p\)-dependent terms are \[\frac{l_1L}{v_1p} + \frac{12C_5}{\sqrt{Z}}\sqrt{N_0pq}.\] Minimizing the above expression equals the following convex optimization regarding \(p\): \(\frac{A}{p}+B\sqrt{p}\), \(A,B>0\), where \(A=\frac{l_1L}{v_1}\) and \(B=\frac{12C_5\sqrt{N_0q}}{\sqrt{Z}}\). \(f(p)=\frac{A}{p}+B\sqrt{p}\) is minimized when \(f'(p)=0\), which implies \(p=\left(\frac{2A}{B}\right)^{\frac{2}{3}}\). That is,

\[\begin{align} p^\star &= \left( \frac{l_1L\sqrt{Z}}{6v_1C_5\sqrt{N_0q}} \right)^{\frac{2}{3}}. \end{align}\]

Thus, the optimal number of interpolation intervals increases with the sample size \(Z\) and decreases with the per-knot freedom \(q\) and maximum sequence length \(N_0\). When \(L\) is set to the lower-bound scale derived above and \(v_1\) is identified with the baseline intensity \(\lambda_{\min}\), because \(\Delta=L/p^\star\), we acquire the corresponding interval scale:

\[\Delta =\frac{L}{p^\star} \le C_{\Delta} \frac{(N_0q)^{\frac{1}{3}}\lambda_{\text{min}}^{\frac{5}{6}}(\mathbb{E} [t_i])^{\frac{1}{6}}}{Z^{\frac{1}{3}}l_1^{\frac{5}{6}}},\] where \(C_{\Delta}\) absorbs the numerical and bracketing-integral constants such as \(C_5\). This finishes the proof.

References↩︎