# Sparklen: A Statistical Learning Toolkit for High-Dimensional Hawkes Processes in Python

Romain Edmond Lacoste

<figure><img src="image_1.png" /></figure>

## To cite this version:

Romain Edmond Lacoste. Sparklen: A Statistical Learning Toolkit for High-Dimensional Hawkes Processes in Python. 2025. hal-04966672v2

HAL Id: hal-04966672

https://hal.science/hal-04966672v2

Preprint submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sparklen: A Statistical Learning Toolkit for High-Dimensional Hawkes Processes in Python

Romain E. Lacoste

Université Gustave Eiffel

## Abstract

This paper introduces **Sparklen**, a statistical learning toolkit for Hawkes processes in Python, designed to bring together efficiency and ease of use. The purpose of this package is to provide the Python community with a complete suite of cutting-edge tools specifically tailored for the study of exponential Hawkes processes, with a particular focus on high-dimensional framework. It includes state-of-the-art estimation tools with built-in support for incorporating regularization techniques, and novel classification methods. To enhance computational performance, **Sparklen** leverages a high-performance C++ core for intensive tasks. This dual-language approach makes **Sparklen** a powerful solution for computationally demanding real-world applications. Here, we present its implementation framework and provide illustrative examples, demonstrating its capabilities and practical usage.

*Keywords:* Hawkes processes, high-dimension, exponential kernel, empirical risk minimization, maximum likelihood estimation, regularization, optimization, Python, C++.

## 1. Introduction: Hawkes process modeling in Python

This paper introduce the Python package **Sparklen** (see Lacoste (2025)), which implements a complete set of statistical learning methods for exponential Hawkes processes with an emphasize on high-dimension setting. Hawkes processes, introduced in Hawkes (1971), form a specific but rather versatile class of point processes. Such processes model time series in which the occurrence of one event temporarily increases the probability of other events occurring. This intrinsic ability to take into account self-exciting effects makes them particularly interesting for real data modeling. Historically applied in seismology (see Ogata (1988)), they have since been used in a wide variety of other fields, including neuroscience in Reynaud-Bouret, Rivoirard, and Tuleau-Malot (2013), finance in Bacry, Mastromatteo, and Muzy (2015), ecology in Denis, Dion-Blanc, Lacoste, Sansonnet, and Bas (2024). The multidimensional version, known as the Multivariate Hawkes Processes (MHP), captures additionally interactions among each univariate process within a network. This generalization enables the modeling of more intricate dynamics, significantly expanding the range of potential applications. For example, MHP has been applied to model action potentials within neural networks in Bonnet, Dion-Blanc, Gindraud, and Lemler (2022), or for trend detection in social networks in Pinto, Chahed, and Altman (2015).
**Software for MHP.** To support this virtuous trend of applications, considerable effort has gone into the development of statistical methods for MHP. Bridging the gap between theoretical advances and practical applications necessitates the development of software that provides comprehensive and user-friendly tools. Over the years, packages for MHP have emerged, aiming to fulfill this need. A key reference in this field is **tick** (Bacry, Bompaire, Deegan, Gaiffas, and Poulsen (2018)), a Python package that provides a list of tools for statistical learning with an emphasis on time-dependent models. Although **tick** offers a complete set of tools for MHP, it faces maintenance issues and is not installable in many recent Python configurations. Other packages, covering more specific uses and therefore being less comprehensive are also available. In the Bayesian framework, one can cite the **pyhawkes** (Linderman (2017)) Python package which proposes inference algorithms for discovering latent network structure given point process observations. In situations where count data are only observed in discrete time, the **hawkesbow** (Cheysson (2024)) R package provides an implementation of an estimation method for univariate Hawkes processes using a spectral approach derived from the Bartlett spectrum. Recently, the **FaDIn** package introduces a novel discretized inference procedure for MHP with parametric kernels of finite support (see Staerman, Allain, Gramfort, and Moreau (2023)). So, despite the growing popularity of MHP, the community lacks a comprehensive, accessible implementation environment that unifies the wide range of tools developed for these processes.

**Sparklen software.** In this context, we introduce the Python package **Sparklen** (see Lacoste (2025)). The aim of this package is to provide to the Python community a complete set of easy to use and efficient statistical tools for MHP, with an emphasis on high-dimension setting. Before delving into the package's features, let us first examine its implementation architecture. From a syntactic standpoint, the code rigorously adheres to object-oriented programming principles and conforms to the widely used **scikit-learn** API (see Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot, and Édouard Duchesnay (2011); Buitinck, Louppe, and Blondel (2013)). This ensures a familiar and efficient framework for users, simplifying and streamlining their workflow. The package's structure is thoughtfully organized, with each functionality carefully compartmentalized. A key strength lies in its use of core C++ code for computationally intensive components, ensuring both efficiency and performance. The binding between C++ and Python is handled through **SWIG** wrapper code (see Beazley (2003)). In essence, **SWIG** facilitates interoperability by generating Python code that wraps C++ classes through written interface files. In a nutshell, this dual-language approach combines the accessibility of Python with the computational efficiency of C++. This makes **Sparklen** highly user-friendly, while also ensuring it is well-suited for computationally demanding applications. Finally, this package is thoroughly documented, further enhancing its accessibility.

**Functionality and scope.** The **Sparklen** library offers a collection of modules for optimization with regularization, featuring classes for optimizers, learning rate schedulers, proximal operators, and calibration methods. These components can be flexibly combined and integrated to meet various needs. The core module of the package is dedicated to exponential MHP. The exponential kernel is a popular choice due to its computational efficiency and its structure, which facilitates straightforward interpretability. At this stage, we wish to highlight
that all the tools implemented in this package are equally applicable to univariate Hawkes processes, as these are naturally included as a special case of the MHP model. The **Sparklen** package integrates several tools for working with MHP, such as an efficient cluster-based simulation method for generating MHP events. The package's main features are categorized into two key areas: inference and multiclass classification methods for MHP, with a particular emphasis on the high-dimensional framework.

**Inference in Sparklen.** First, let us outline the inference tools provided in the **Sparklen** package. In the literature of MHP, particular attention has been paid to the design of efficient estimating procedures. For example, maximum likelihood approaches have been proposed, one can cite Ogata (1978). In addition, methods relying on least square contrast exist in the literature (Reynaud-Bouret, Rivoirard, Grammont, and Tuleau-Malot (2014); Bacry, Bompaire, Gaiffas, and Muzy (2020)). In this article, we consider an original setting, assuming that the data we have access to consist of repeated short-time paths. For this observational setup, one can cite the work of Lotz (2024) for likelihood approach and Dion-Blanc, Denis, Sansonnet, and Lacoste (2024) for least-squares. The **Sparklen** package implements both likelihood and least-squares approaches, handling the aforementioned setup. A key consideration is that the loss, gradient, and hessian of both functionals are computed exactly for precision and efficiently for speed. Furthermore, several regularizations are implemented, allowing users to incorporate prior knowledge about the nature of the problem. Notably, the tuning of the regularization constant, optimally selected internally based on various user-selectable criteria, is a very appealing contribution to this package. This makes the implemented tools highly versatile in a variety of situations. Specifically, the procedures are tailored for high-dimensional settings. In particular, the use of Lasso regularization combined with automatic tuning of the penalty constant can lead to successful recovery of the support of a high-dimensional MHP.

**Classification in Sparklen.** Second, let us describe the multiclass classification setting we consider and the tools implemented in the **Sparklen** package to address such problem. Time series classification is a fast-growing and flourishing field of research, the fruit of a vigorous research effort in recent years (see Ismail Fawaz, Forestier, Weber, Idoumghar, and Muller (2019)). Due to their complex nature, standard learning algorithms are generally not tailored for working on raw time series. On the other hand, the problem of classifying time series, which consists in labeling data given that each class is characterized by a specific underlying temporal dynamic, naturally arises in real life. Examples of applications are legion, for example for ECG pattern recognition in Chaovalitwongse, Prokopyev, and Pardalos (2006), for fault detection and cyber-security applications in Susto, Cenedese, and Terzi (2018). As a result, there is a strong need to develop cutting-edge algorithms to process such data. Recently, some work has focused on the classification of Hawkes processes paths. For example, work in natural language for rumour stance classification can be found in Lukasik, Srijith, Vu, Bontcheva, Zubiaga, and Cohn (2016). In our context, we assume that the learning sample consists of *i.i.d.* labeled data where the features are the jump times of a MHP observed on the fixed time interval. In such a setting, an application can be found in ecology for modeling echolocation calls of bats to classify their behavior in several sites throughout France in Denis et al. (2024). The used classification procedure follows the strategy presented in Denis, Dion-Blanc, and Sansonnet (2022) in the univariate case and have been recently extended to the multivariate case in Dion-Blanc et al. (2024). This novel procedure takes advantage of an
initial interaction recovery step, by class, followed by a refitting step based on a suitable classification criterion. The **Sparklen** package implements the two classification procedures. While both rely on empirical risk minimization, the second method, which includes a prior recovery step using a Lasso-type estimator, enables handling of high-dimensional settings.

**Outline of the paper.** The remainder of this paper is structured as follows. In Section 2, we provide an overview of the theoretical foundations underpinning **Sparklen**'s features. Section 3 presents the implementation and provides insightful examples that highlight **Sparklen**'s capabilities. Then in Section 4, the package is illustrated on MemeTracker data set. Finally, in Section 5, we conclude the paper and outline potential future developments for the package.

## 2. Hawkes processes

First, in Section 2.1, we introduce the Hawkes model. Then, Section 2.2 presents an efficient simulation algorithm that leverages its branching structure. Section 2.3 covers the inference tools available for estimating a MHP. Finally, in Section 2.4, we present the supervised classification framework and classification algorithms tailored to MHP paths. In the following, the set $\{1, \dots, n\}$ is denoted by $[n]$ for any integer $n \ge 1$.

### 2.1. Hawkes model

In many fields of application, data takes the form of points distributed in a space. For instance, it is usual to observe points distributed along the real line, such as the occurrence of events over time. A pivotal issue, for modeling and prediction purpose, is to provide a suitable probabilistic framework to describe the temporal distribution of these events. Point processes, with their associated formalism, constitute a powerful and flexible tool for modeling such distributions.

A point process, denoted $N$, is encoded by its events arrival times $(t_\ell)_{\ell \ge 1}$, given by the position of the points in $\mathbb{R}_+^*$. In addition, it can be described by its counting process, denoted $(N(t))_{t \ge 0}$ and defined by $N(t) = \sum_{\ell \ge 1} 1_{\{t_\ell \le t\}}$, which gives the number of events that have occurred before time $t$. In the paradigm of point processes, a key component is the intensity function, which characterizes the process by dictating the dynamics of event occurrences. Such function, for a given $t \ge 0$, takes the following form:

$$ \lambda(t) := \lim_{h \to 0^+} \frac{\mathbb{P}(N(t+h) - N(t) = 1 | \mathcal{F}_{t-})}{h} $$

where the filtration $\mathcal{F}_{t-}$ represents the information available strictly before time $t$. Heuristically, at a given time, it may be perceived as the infinitesimal probability of observing an event in the near future, conditionally on the past of the process.

A classical example of point processes is given by the homogeneous Poisson process, which models the occurrence of events at rate $\lambda(t) = \mu$. In this case, the intensity is not random and does not depend on $t$. In particular, successive events are independent of one another. Despite being a quite popular tool, this feature of independence is not always realistic from a modeling standpoint. Indeed, in many applications, the observed dynamic induces more tricky temporal dependencies that cannot be captured by such a model. On the other hand,
Hawkes processes are tailored to model time series for which the occurrence of events depends on past observations.

This work centers on the linear Hawkes process with an exponential kernel, a well-regarded choice due to its favorable properties. Notably, the recursive formulation of the likelihood enables efficient computation, while the model’s structure further supports straightforward interpretability, enhancing its usability in practice.

The multivariate version of these processes, called Multivariate Hawkes Processes (MHP), stands out as a natural extension enabling to describe more complex dynamics, thus broadening the spectrum of applications. In such a model, the evolution of a component’s dynamics is no longer only likely to be impacted by the presence of self-exciting effects, but also by exciting interactions of components connected to it within the network. A d-MHP, where $d \ge 1$ is the network’s dimension, is denoted as $N = (N_1, \dots, N_d)$. Each component $N_j$ is characterized by its intensity function, defined for each $t \ge 0$ as follows:

$$ \lambda_j(t) := \mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \int_{0}^{t} \beta e^{-\beta(t-s)} dN_{j'}(s) = \mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \sum_{\ell:t_{j',\ell}<t} \beta e^{-\beta(t-t_{j',\ell})} \quad (1) $$

where $\mu = (\mu_j)_{j \in [d]}$ is the vector of exogenous intensities, $\alpha = (\alpha_{j,j'})_{j,j' \in [d]}$ is the matrix of interactions, and $\beta$ is the decay. For each $j \in [d]$, the coefficient $\mu_j > 0$ expresses the arrival of spontaneous events for the j-th process. For each $j, j' \in [d]$, the coefficient $\alpha_{j,j'}$ is non-negative and reflects the positive influence of $N_{j'}$ on $N_j$. Finally, the scalar $\beta > 0$ dictates how quick these influences vanish over time.

## 2.2. Simulation and branching structure

An appealing property of Hawkes processes is their branching structure. In particular, they can be represented as Poisson cluster processes, which are age-dependent immigration-birth processes. Consequently, each cluster forms a continuous-time Galton-Watson tree with a Poisson offspring process. Introduced in Hawkes and Oakes (1974) for the univariate case, this bridge between branching theory and Hawkes processes enables to derive theoretical results of interest such as existence, asymptotic behavior as well as counting and interval properties of these processes. In addition to offering a rich theoretical framework, this cluster representation gives rise to an efficient simulation procedure (presented in Møller and Rasmussen (2005); Chen and Wang (2020)). Indeed, by exploiting the hierarchical structure of descendants, one can simulate the exact number of events without any rejection, in contrast to the classical thinning approach (see Ogata (1981)).

Consider a d-MHP and let $T > 0$ the (non-necessarily finite) horizon time. The simulation of the events on $[0, T]$ for this MHP is as follows. First, for each $j \in [d]$, we simulate the immigrants of the j-th component according to a homogeneous Poisson process with intensity $\lambda = \mu_j$ on $[0, T]$. These events are denoted $\{t_{j,\ell}^0\}_{\ell \in [m_j^0]}$ with $m_j^0 \sim P(\mu_j T)$. These immigrants become the ancestors and we must simulate their offspring. For each j and for each $\ell \in [m_j^0]$, consider the $\ell$-th ancestor of type j at generation 0. The latter, denoted $t_{j,\ell}^0$, generates descendants of generation 1 in each dimension $j' \in [d]$ according to a inhomogeneous Poisson process with intensity function $\lambda_{j,j'}(\cdot - t_{j,\ell}^0) = \alpha_{j,j'}\beta \exp(-\beta(\cdot - t_{j,\ell}^0))$ on $[t_{j,\ell}^0, T]$. For $j \in [d]$, the events are denoted $\{t_{j,\ell}^1\}_{\ell \in [m_j^1]}$ with $m_j^1 \sim P(\alpha_{j,j'}(1 - \exp(-\beta(T - t_{j,\ell}^0)))$. Then, in the
next step, these descendants become the ancestors and the descendants of generation 2 can be simulated, and so on. We continue iterating until there are no ancestors left in each dimension. This construction is precisely detailed in Algorithm 1.

<caption>**Algorithm 1:** Simulation by cluster representation</caption>

**Input:** $\mu \in \mathbb{R}_+^d$, $\alpha \in \mathbb{R}_+^{d \times d}$, $\beta > 0$ and $T > 0$
**Initialization:** $\mathcal{T} = \{\mathcal{T}_j | \mathcal{T}_j = \emptyset, j \in [d]\}$, $\mathcal{A} = \{\mathcal{A}_j | \mathcal{A}_j = \emptyset, j \in [d]\}$
**Immigrant simulation**

1. **for** $j = 1, \dots, d$ **do**
2. &nbsp;&nbsp;&nbsp;&nbsp;$k \sim \mathcal{P}(\mu_j T)$; &nbsp;&nbsp;&nbsp;&nbsp;▷ Number of immigrants of type $j$
3. &nbsp;&nbsp;&nbsp;&nbsp;$t_{j,1}, \dots, t_{j,k} \stackrel{\text{i.i.d.}}{\sim} \mathcal{U}([0, T])$;
4. &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{A}_j = \{t_{j,1}, \dots, t_{j,k}\}$;
5. &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{T}_j = \mathcal{T}_j \cup \mathcal{A}_j$;
6. **end**

**Offspring simulation**

7. **while** there exists at least one $j$ such that $\mathcal{A}_j \neq \emptyset$ **do**
8. &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{O} = \{\mathcal{O}_j | \mathcal{O}_j = \emptyset, j \in [d]\}$; &nbsp;&nbsp;&nbsp;&nbsp;▷ Offspring initialization
9. &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, \dots, d$ **do**
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $\mathcal{A}_j \neq \emptyset$ **then**
11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each $a_{j,\ell} \in \mathcal{A}_j$ **do**
12. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $j' = 1, \dots, d$ **do**
13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $\alpha_{j',j} > 0$ **then**
14. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \sim \mathcal{P}(\alpha_{j',j})$;
15. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$t_{j',1}, \dots, t_{j',k} \stackrel{\text{i.i.d.}}{\sim} \mathcal{E}(\beta)$;
16. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{O}_{j'} = \mathcal{O}_{j'} \cup \{a_{j,\ell} + t_{j',1}, \dots, a_{j,\ell} + t_{j',k}\}$;
17. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
18. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{T}_{j'} = \mathcal{T}_{j'} \cup \mathcal{O}_{j'}$;
19. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
20. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
21. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
22. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
23. &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{A} = \mathcal{O}$; &nbsp;&nbsp;&nbsp;&nbsp;▷ Offspring become ancestors

24. **end**

**Remove events beyond $T$ and sort**

25. **for** $j = 1, \dots, d$ **do**
26. &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{T}_j = \text{sort}(\{t_{j,\ell} | t_{j,\ell} < T\})$;

27. **end**

**Output:** $\mathcal{T} = \{\mathcal{T}_j, j \in [d]\}$

As outlined in Algorithm 1, the descendants are simulated according to an inhomogeneous Poisson process on $[a_{j,\ell}, +\infty[$, followed by the removal of events exceeding $T$. By doing so, the number of events simplifies into $\alpha_{j,j'}$ as can be seen in line 14 in Algorithm 1. This prevents from calculating the integral for every ancestor $a_{j,\ell}$ and each coefficient $\alpha_{j,j'}$. This trick helps accelerating the sampling procedure in practice.

## 2.3. Parametric inference methods
Before presenting the parametric estimation strategies, let us first define the statistical framework we consider. As it was implied in Equation (1), we suppose that the intensity function of the process N depends on an unknown parameter $\theta^*$. It is therefore through this unknown parameter that the global dynamic of the process is encoded. In contrast, the parameter $\beta$ is assumed to be known for optimization purposes (see below). The unknown parameter of interest, $\theta^* = (\mu^*, \alpha^*)$, is supposed to belong to the following parameter family:

$$ \Theta := \{\mu \in (\mathbb{R}_+^*)^d, \alpha \in \mathbb{R}_+^{d \times d}, \rho(\alpha) < 1\} \subseteq \mathbb{R}_+^{d \times d+1} $$

with $\rho(\alpha)$ the spectral radius of the matrix $\alpha$. This hypothesis, known as sub-critical stability hypothesis, is crucial to ensure the non-explosion of the process.

Let $T > 0$ be fixed. We denote $\mathcal{T}_T := \{\{t_{j,l}\}_{l \in [N_j(T)]}, j \in [d]\}$ the event times over $[0, T]$ of a MHP $N = (N_1, \dots, N_d)$ of intensity $(\lambda_{1,\theta^*}, \dots, \lambda_{d,\theta^*})$. We assume that we have access to a $n$-sample $\mathcal{D}_n := \{\mathcal{T}_T^{(1)}, \dots, \mathcal{T}_T^{(n)}\}$ which consists of independent copies of $\mathcal{T}_T$. In other words, we are dealing with short-time path repetitions. In particular, we do not assume that observed paths have reached the stationary regime. Given the observation of the sample $\mathcal{D}_n$, our strategy is to minimize the following optimization problem:

$$ \hat{\theta}_n(\kappa) \in \arg\min_{\theta \in \mathbb{R}^{d \times d+1}} \{F_{T,n}(\theta) + \kappa \Omega(\theta)\} \quad (2) $$

where $F_{T,n}(\cdot)$ is a loss function describing how well the model fits the data, $\Omega(\cdot)$ a regularization function that penalizes complex solutions and $\kappa > 0$ the regularization constant which controls the level of inductive bias.

The minimization problem given by Equation (2) is the cornerstone of this section. In what follows, we present an overview of the various options available for loss functions, regularization functions, penalty constants choice, and optimization methods. Each option is examined in light of current advancements in the field and its relevance in the context of the application concerned.

**Loss function.** We focus on two parametric estimation approaches, distinguished by the choice of the following two loss functions:

* Log-likelihood functional, $F_{T,n}(\theta) = L_{T,n}(\theta)$,

* Least-squares contrast, $F_{T,n}(\theta) = R_{T,n}(\theta)$.

The first, of more classical inspiration, is given by the maximum likelihood estimation method. Introduced in Ogata (1978), the Maximum Likelihood Estimator (MLE) is shown to be consistent, asymptotically normal and asymptotically efficient when the time horizon $T$ tends to infinity for a stationary Hawkes process. In the setting we consider, namely that of repeated observations with $T$ fixed, the MLE is shown to be consistent as the number of observations $n$ tends to infinity (see Lotz (2024)). From a practical point of view, we rather consider the negative log-likelihood in order to use first-order optimization algorithms to find its minimizer
(see below). The negative log-likelihood over the repetitions is given by:

$$L_{T,n}(\theta) := \sum_{i=1}^{n} \left( \sum_{j=1}^{d} \int_{0}^{T} \lambda_{j,\theta}^{(i)}(t) \, dt - \sum_{j=1}^{d} \sum_{\ell:t_{j,\ell}^{(i)} < T} \log \left( \lambda_{j,\theta}^{(i)}(t_{j,\ell}^{(i)}) \right) \right)$$

where $\lambda_{j,\theta}^{(i)}(\cdot)$ is defined from the observation $\mathcal{T}_T^{(i)}$. An alternative approach, receiving comparatively less attention, is founded on the least-squares contrast, drawing inspiration from the theory of empirical risk minimization (Reynaud-Bouret et al. (2014); Bacry et al. (2020)). Such a functional benefits from attractive computational and optimization-related properties (as discussed below). Recently, in Dion-Blanc et al. (2024), within our observational setup of repeated short-time paths, the authors propose a Lasso-type procedure based on this loss to handle high-dimensional issues for MHP. Notably, a result of uniform consistency of the Lasso estimator as well as the consistency of the estimated support is proven when $n$ tends to infinity. The least-squares contrast averaged over the repetitions takes the following form:

$$R_{T,n}(\theta) := \frac{1}{n} \sum_{i=1}^{n} \left( \frac{1}{T} \sum_{j=1}^{d} \int_{0}^{T} \lambda_{j,\theta}^{(i)2}(t) \, dt - \frac{2}{T} \sum_{j=1}^{d} \sum_{\ell:t_{j,\ell}^{(i)} < T} \lambda_{j,\theta}^{(i)}(t_{j,\ell}^{(i)}) \right).$$

For both functions, it is possible, by simplification, to isolate the terms that depend exclusively on the exponential function and the event arrival times. As these terms are independent of the parameter, they can be calculated once and reused thereafter. This calculation trick speeds up future evaluations of the loss, gradient or Hessian. The implementation adopts this approach for optimal efficiency. For more details, see Appendix A.

In addition, it is worth noting that the inclusion of $\beta$ among the parameters to be estimated would result in the loss of convexity for both functionals, thereby complicating the optimization process. As a result, the decision to fix $\beta$, as outlined previously, is a deliberate strategy aimed to circumvent such difficulties. In situations where $\beta$ is unknown, the practitioner may explore a grid of candidate values and select the most suitable $\beta$ based on an *ad hoc* criterion. In Table 1, a comparison of the two functions properties is provided.

<table><caption>Table 1: Comparison of the log-likelihood and least-squares functions in terms of algorithmic complexity and optimization properties. The total number of observed events is denoted by $m$. Note that for least-squares, once pre-calculations have been performed, the complexity of both the loss function and its gradient evaluation remains independent of $m$.</caption><thead><tr><th rowspan="2">Function</th><th colspan="3">Complexity</th><th colspan="2">Property</th></tr><tr><th>Pre-calculation</th><th>Loss</th><th>Gradient</th><th>Convex</th><th>L-Smooth</th></tr></thead><tbody><tr><td>Log-likelihood</td><td>O(md)</td><td>O(md)</td><td>O(md)</td><td>✓</td><td>✗</td></tr><tr><td>Least-squares</td><td>O(md)</td><td>O(d<sup>3</sup>)</td><td>O(d<sup>3</sup>)</td><td>✓</td><td>✓</td></tr></tbody></table>

**Regularization.** In machine learning, regularization encompasses a set of techniques that emerge naturally to prevent the effects of over-fitting and increase generalization ability (see
Bach (2024) for more details). The incorporation of a regularization term, inducing constrained optimization, reflects a compromise between minimizing the objective function and maintaining the structural simplicity of the model through the hyperparameter $κ$. The choice of the regularization function $\Omega(\cdot)$ may be dictated by a prior knowledge of the structure of the parameter $\theta^*$. In this paper, we consider three types of regularization:

* Lasso regularization by adding a $l_1$-penalty $\Omega(\cdot) = \|\cdot\|_1$ term (see Tibshirani (1996))

* Ridge (or Tikhonov) regularization by incorporating a $l_2$-penalty $\Omega(\cdot) = \|\cdot\|_2^2$ (see Hoerl and Kennard (1970))

* Elastic-Net regularization by considering a linear combination of $l_1$-penalty and $l_2$-penalty $\Omega(\cdot) = \zeta \|\cdot\|_1 + (1 - \zeta) \|\cdot\|_2^2$ (see Zou and Hastie (2005))

Each of these have its own typical use case. For instance, in the case of a large network with few effective interactions, one may use Lasso regularization for its intrinsic virtue of variable selection and improved interpretability.

**Calibration of $\kappa$.** The penalization constant $\kappa$ controls the power of the induced regularization. Small values lead to complex models that are likely to over fit, while on the other hand large values lead to simplistic models with poor prediction power. The question of choosing this value, which encapsulates the trade-off described before, is therefore decisive for the success of the learning phase. In this article, we examine two methodologies for making this choice:

* Cross validation (CV) (e.g. Arlot and Celisse (2010))

* Extended Bayesian Information Criterion (EBIC) (e.g. Chen and Chen (2008))

For both methods, the strategy involves exploring a grid $\Delta$ of $\kappa$ values and selecting the one that minimizes the criterion of interest. In the case of CV, this entails partitioning the dataset into folds and aiming to minimize the cross-validated risk. In the case of EBIC, given $\gamma \in [0, 1]$, the criterion takes the following form:

$$ \text{EBIC}_{\gamma}(\kappa) := -2L_{T,n}(\hat{\theta}_n(\kappa)) + |S_{\hat{\theta}_n(\kappa)}| \log(n) + 2\gamma \log\left(\left(\frac{M^2}{|S_{\hat{\theta}_n(\kappa)}|}\right)\right) $$

where $\kappa \in \Delta$ is the penalization constant under consideration; $\hat{\theta}_n(\kappa)$ is the minimizer of the problem given by Equation (2) (with any function $F_{T,n}(\cdot)$ and $\Omega(\cdot)$); $L_{T,n}$ is the log-likelihood of the model, $|S_{\hat{\theta}_n(\kappa)}|$ is the size of this support, namely the number of active coefficients of $\hat{\theta}_n(\kappa)$. Compared to a classical BIC criterion (namely $\gamma = 0$), an additional penalization is introduced. It that accounts for the number of possible active sets of the same size and is designed to better handle the high-dimensional setting.

**Optimization.** This paragraph is dedicated to the description of tools for solving the convex optimization problem given by Equation (2). The objective function is written as the sum of two functions. While the loss function $F_{T,n}(\cdot)$ is convex, differentiable and may even
be smooth (in the case of least-squares contrast), the regularization function $\Omega(\cdot)$ may not be differentiable (the $\ell_1$-norm is not differentiable in 0). To this extent, to carry out the minimization of such objective function, we consider first-order optimizations algorithm based on proximal methods. We focus on two optimization algorithms:

* gradient descent (GD) (see Bach (2024)),

* accelerated gradient descent (AGD) (see Nesterov (1983)).

The GD algorithm is often regarded as the workhorse of first-order optimization, known for its simplicity and effectiveness across a broad range of applications. While straightforward in implementation, it can be slow to converge, especially in poorly conditioned settings. AGD refines this classic approach by incorporating a momentum term that strategically combines information from the two most recent iterates. Although this refinement introduces a slight additional computational cost, it is marginal, as it requires only a single gradient evaluation per iteration, as in standard GD. In exchange, it achieves a significantly faster convergence rate, often leading to a considerable reduction in the number of iterations required to achieve the desired level of accuracy.

In our case, both of these optimizers are coupled with the proximal operator to incorporate regularization. For the considered regularization, the corresponding proximal operators admit closed-form solutions. In particular, in the case of Lasso regularization, the proximal operator reduces to the soft-thresholding function and theoretical properties of both optimizers, known in the literature as ISTA and FISTA, are well-studied and thoroughly documented (see Beck and Teboulle (2009)).

Finally, a critical consideration in Gradient Descent is the selection of an appropriate step size. When the loss function has a Lipschitz continuous gradient, a commonly recommended step size is $1/L$, where $L$ denotes the Lipschitz constant of the gradient. In cases where this assumption does not hold, a backtracking line-search method can be employed to dynamically determine an effective step size at each iteration, ensuring stable and efficient convergence of the descent process (see Armijo (1966)).

## 2.4. Multiclass classification

Before introducing classification algorithms tailored for MHP path classification, let us first define the supervised multiclass classification model.

**Multiclass classification model.** We consider the setting of $K$-multiclass classification with $K \ge 2$ the number of classes. Let $(\mathcal{T}_T, Y)$ be the couple of variables formed by the feature and its associated label. In this case, the feature $\mathcal{T}_T$ consists of time events of a MHP $N = (N_1, \dots, N_d)$ observed on $[0, T]$ whose dynamics of occurrence is governed by its label $Y \in [K]$. More precisely, conditional on $Y$, $N$ is a MHP of intensity $\lambda_{\theta_Y^*} = (\lambda_{\theta_Y^*,1}, \dots, \lambda_{\theta_Y^*,d})$. For each $j \in [d]$, the intensity function of the process $N_j$ depends on $Y$ and is given by:

$$ \lambda_{\theta_Y^*,j}(t) = \mu_{Y,j}^* + \sum_{j'=1}^{d} \alpha_{Y,j,j'}^* \int_{0}^{t} \beta e^{-\beta(t-s)} dN_{j'}(s). $$

So in this model, the classes are discriminated according to the exogenous intensity vectors $(\mu_Y^*)_{Y \in [K]}$ and the interaction matrices $(\alpha_Y^*)_{Y \in [K]}$. This implies that classes are characterized
by their underlying rate of spontaneous activity, as well as by the intrinsic structure of their interaction network. In the following, we assume that the parameter $\theta^* = (\mu_k^*, \alpha_k^*)_{k \in [K]}$ is unknown as well as the distribution of Y which is denoted by $p^* = (p_k^*)_{k \in [K]}$. On the other hand, the decay rate $\beta$ is assumed to be known and identical for each class. Such a modeling assumption implies that the speed at which influences between network components vanish over time cannot be used as a benchmark to distinguish between classes.

We follow the paradigm of supervised classification: although the law of the pair $(\mathcal{T}_T, Y)$ is unknown, we assume we have access to a learning sample $\mathcal{D}_n^L = \{(\mathcal{T}_T^{(1)}, Y^{(1)}), \dots, (\mathcal{T}_T^{(n)}, Y^{(n)})\}$ of size $n$ which consists of independent copies of $(\mathcal{T}_T, Y)$. From training data the goal is to learn a decision rule $\hat{g}$, i.e. a function of the data with values in $[K]$, which accurately predicts the class to which an observation belongs.

**ERM classifier.** As the distribution of the couple $(\mathcal{T}_T, Y)$ is unknown, it is not possible to calculate and therefore minimize the risk of misclassification. To mimic the Bayes classifier, which by definition minimizes this risk, a classic idea consists in minimizing the point error committed on each observation of the training sample over a chosen hypothesis space of classifier functions. This is the principle behind the empirical risk minimization (ERM) method, which has been covered abundantly in the literature (see Vapnik (1991); Massart and Nédélec (2006)).

Let us now present the ERM classification procedure tailored to classify Multivariate Hawkes Processes paths. This procedure was introduced in Dion-Blanc et al. (2024), see for more details on its principle and asymptotic properties. First, we estimate $p^*$ by its empirical counterpart $\hat{p}$. Then, rather than minimizing the empirical risk with the classic loss $0-1$, which leads to difficult optimization problems, we instead consider a convex surrogate (see Zhang (2004); Bartlett, Jordan, and McAuliffe (2006)).

Let $\hat{\Theta}$ be the set of interest, $\theta \in \hat{\Theta}$, and $f_{\hat{p},\theta}$ its associated score function (see Dion-Blanc et al. (2024) for more details). Consider the $L_2$-risk given by the following form:

$$ \hat{\mathcal{R}}_2(f_{\hat{p},\theta}) := \frac{1}{n} \sum_{i=1}^{n} \sum_{k=1}^{K} \left( Z_k^{(i)} - f_{\hat{p},\theta}^k (\mathcal{T}_T^{(i)}) \right)^2, \quad Z_k^{(i)} = 2 \mathbf{1}_{\{Y_i=k\}} - 1. \quad (3) $$

Then, we define the estimator of $\theta^*$ as the minimizer of the empirical $L_2$-risk:

$$ \hat{\theta}^R \in \arg \min_{\theta \in \hat{\Theta}} \hat{\mathcal{R}}_2(f_{\hat{p},\theta}). \quad (4) $$

From the estimator $\hat{\theta}^R$ of parameter $\theta^*$, we define the ERM classifier as follows

$$ \hat{g}(\mathcal{T}_T) \in \arg \max_{k \in \mathcal{Y}} f_{\hat{p}, \hat{\theta}^R}^k (\mathcal{T}_T) \quad (5) $$

Having outlined the procedure, we now delve deeper into the practical aspects of solving the optimization problem presented in Equation (4). To ensure the positivity of each coefficient, the minimization is conducted under inequality constraints, using a projected-type gradient descent algorithm. However, since this objective function is non-smooth and non-convex w.r.t. to the coefficients, its minimization demands careful consideration. In particular, tuning the step size in the descent process proves to be quite challenging. Furthermore, classical method
such as backtracking line-search with Armijo-Wolfe condition cannot be used due to the piece-wise constant nature of the projection operator (see Michael W. Ferry and Zhang (2023)).

On the other hand, adaptive gradient methods, such as ADAGRAD (see Duchi, Hazan, and Singer (2011)), have become widely adopted in large-scale optimization due to their ability to adjust the step size for each feature according to the geometry of the problem. In practice, ADAGRAD is recognized as an efficient method in non-convex settings, particularly for training deep neural networks (see Gupta, Bengio, and Weston (2014)). In addition, some theoretical guarantees for the convergence of ADAGRAD for non-convex functions have been established in the literature (see Ward, Wu, and Bottou (2020); Wang, Zhang, Ma, and Chen (2023)).

With this in mind, we use a parameter-free projected adaptive gradient descent method, in the likeness of ADAGRAD, called FREE ADAGRAD and introduced in Chzhen, Giraud, and Stoltz (2023). Compared with the classical algorithm, its key advantage lies in the fact that it is adaptive to the distance between the initialization and the optimum, and to the sum of the square norm of the gradients. The only input parameter of the algorithm is $\gamma_0$, which can be interpreted as an initial lower-bound guess for the distance between the starting point and the optimum. In practice, we set $\gamma_0 = 0.1$.

**High-dimension.** Consider the high-dimensional framework, meaning the dimension of the network may be large w.r.t. the number of observations. In such setting, the computational burden of the ERM procedure becomes substantial and the predictor $\hat{g}$ may not be consistent. Fortunately, high-dimensional data frequently exhibit a structure that is significantly more low-dimensional than their apparent complexity suggests. In this context, it is thus reasonable to postulate a sparse underlying structure and to endeavor to recover this lower-dimensional representation. Here, the sparsity hypothesis results in a sparse interaction matrix for each class.

In line with this, Dion-Blanc et al. (2024) introduced an ERMLR procedure whose methodology is the following. First, an initial interaction recovery step, by class, is performed using a Lasso-type estimator. Then, taking advantage of the estimation of the support of the adjacency matrix, by class, a refitting step based on an ERM criterion is carried out. In particular, the optimization of the $L_2$-risk is done on the following set of parameters:

$$ \hat{\Theta}^R := \{ \theta \in \hat{\Theta}, S_{\theta_k} = \hat{S}_k \} $$

where $\hat{S}_k$ is the estimated support of the $k$-th class given by the Lasso estimator. Hence, the minimization is performed on a set of parameters whose dimension is much smaller than $d^2$. The lasso step is significant as it selects effective interactions, thereby alleviating the complexity of the problem. On the other hand, an induced and undesired effect of $\ell_1$ penalization is the shrinkage of large coefficients. To bypass this issue, refitting strategies are commonly used and well covered in the literature, see for example Chzhen, Hebiri, and Salmon (2019). Performing a refitting stage is, therefore, a natural progression, particularly in this context, as we employ a criterion specifically tailored for classification. The advantages of this approach are threefold: first, the selection of interactions enhances interpretability; second, it improves the classifier's performance; and third, it reduces the overall computational cost of the procedure. The ERMLR procedure is delineated comprehensively in Algorithm 2.
<table><caption>**Algorithm 2:** ERMLR procedure</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>for</td><td>k = 1, ..., K do</td></tr><tr><td>2</td><td>|</td><td>p̂<sub>k</sub> = <sup>1</sup>&frasl;<sub>n</sub> &sum;<sub>i=1</sub><sup>n</sup> 1{Y<sup>i</sup>=k};</td></tr><tr><td>3</td><td>end</td><td></td></tr><tr><td>4</td><td>for</td><td>k = 1, ..., K do</td></tr><tr><td>5</td><td>|</td><td>Choose k̂ ∈ arg min<sub>κ∈Δ</sub> EBIC<sub>1</sub>(κ);</td></tr><tr><td>6</td><td>Solve</td><td>θ̂<sub>k</sub> ∈ arg min<sub>θ∈R<sup>d&times;d+1</sup></sub> {R<sub>T,n</sub>(θ) + k̂<sub>k</sub> ||θ||<sub>1</sub>} using FISTA;</td></tr><tr><td>7</td><td>Get</td><td>Ŝ<sub>k</sub> = {θ̂<sub>k,j,j'</sub> ≠ 0, (j,j') ∈ [d]};</td></tr><tr><td>8</td><td>end</td><td></td></tr><tr><td>9</td><td>ERM</td><td>reffiting</td></tr><tr><td>10</td><td>Solve</td><td>θ̂<sup>R</sup> ∈ arg min<sub>θ∈θ̂<sup>R</sup></sub> R̂<sub>2</sub>(f<sub>p̂,θ</sub>) using FREE ADAGRAD starting from (θ̂<sub>k</sub>)<sub>k∈[K]</sub>;</td></tr><tr><td></td><td>Get</td><td>ĝ(·) ∈ arg max<sub>k∈Y</sub> f<sup>k</sup><sub>θ̂<sup>R</sup></sub>(·);</td></tr><tr><td></td><td>Output:</td><td>ĝ(·)</td></tr></tbody></table>

# 3. Package presentation and functionality

This section begins with an overview of the global architecture and features of the package in Section 3.1. This is followed by a detailed explanation of the implementation of the statistical methods in Section 3.2, 3.3, 3.4, closely mirroring the structure of Section 2. In these sections, comprehensive and intuitive guides are provided on how to apply these methods in practice, complemented by insightful examples for better illustration.

## 3.1. Package overview

**Global architecture.** The directory structure of the package is thoughtfully organized to clearly compartmentalize its various components. At the core, the lib/ folder houses C++ header and source files, along with SWIG interface files. These elements are not intended for direct user interaction; rather, they serve the primary purpose of enabling rapid computation, allowing the package to perform complex calculations efficiently while abstracting away the underlying C++ implementation. All user-ended Python code resides exclusively within the sparklen/ directory, which encompasses the suite of modules and sub-modules offered. This directory serves as the primary interface for users, encapsulating the entire functionality of the package. It is meticulously organized to facilitate ease of use and maintainability, allowing users to access and possibly extend the package's capabilities if needed. The setup.py script plays a pivotal role in the package's architecture by orchestrating the build process. It ensures that SWIG compiles the interface files correctly into Python extension modules, making the high-performance C++ functionalities accessible from Python. This script also manages dependencies and installation requirements, streamlining the installation process for users. Additionally, a README.md file is included to offer documentation, detailing installation instructions, link to usage examples, and other resources to assist users in effectively navigating the package. A visual representation of this structure is given in Figure 1.
<figure><img src="image_3.png" /><figcaption>Figure 1: Tree structure of Sparklen library directory.</figcaption></figure>

**Features of the Sparklen package.** We will now delve into the features of the package, encapsulated within the `sparklen/` directory as a collection of modules and sub-modules. Each module is named after the package, with sub-modules separated from the parent module name by a dot, following Python's standard attribute syntax. This naming convention not only promotes clarity, but also makes it easier to navigate and use. There are essentially four modules: `sparklen.hawkes`, `sparklen.optim`, `sparklen.prox`, `sparklen.calibration`.

The flagship and most crucial module is `sparklen.hawkes`. Such a module integrates all functionalities associated with Hawkes processes, each organized within the following four sub-modules. The sub-module `sparklen.hawkes.simulation` is responsible for simulating Hawkes processes. The `sparklen.hawkes.model` manages the calculations related to the model, such as the evaluation of the loss or the calculation of the gradient for the log-likelihood and least-squares models. The `sparklen.hawkes.inference` contains all the tools for carrying out Hawkes processes inference. In particular, it is based on the previous sub-module and uses the `sparklen.optim` modules for the optimization phase, `sparklen.prox` to manage the integration of a regularization term and `sparklen.calibration` to choose the penalty constant. Finally in `sparklen.hawkes.classification` are implemented tools for classification. See the following sections for more details on these sub-modules. A visual representation of modules, sub-modules and potential dependencies is given in Figure 2.

**C++/Python binding.** An appealing feature of this package lies in its ability to leverage the strengths of both C++ and Python. By implementing the computationally intensive components in C++ (like gradient evaluation), we ensure that the most resource-heavy tasks are executed with optimal efficiency, taking full advantage of C++'s speed and performance. At the same time, the C++ code is seamlessly wrapped in Python using SWIG, which provides a user-friendly, Pythonic interface. This approach offers the best of both worlds: users benefit from the fast, low-level computations of C++, while interacting with high-level, intuitive, flexible Python objects that are easy to manipulate and integrate into their workflows. This dual-language strategy makes it possible to deal with resource-intensive applications (e.g. inference of high-dimensional Hawkes processes) while guaranteeing ease of use. A concise
<figure><img src="image_4.png" /><figcaption>Figure 2: The modules and sub-modules of **Sparklen** package are shown, in the form of a flowchart, along with potential dependencies specified by dashed arrows.</figcaption></figure>

yet comprehensive explanation of the mechanisms underpinning the interface between C++ and Python is provided in Appendix B.

In the following, we present the main features of the module `sparklen.hawkes`, namely simulation, inference and classification. This part complements Section 2 by shedding light on how these tools are implemented in practice, with examples provided for added clarity. In each example the workflow, illustrated with code snippets, is interspersed with explanatory comments to enhance explicitness.

## 3.2. Implementation of the simulation framework

**Description.** Everything related to the simulation of Hawkes processes is contained in the sub-module `sparklen.hawkes.simulation`. The centerpiece of this sub-module is the `SimuHawkesExp` class. Each object in this class is instantiated by giving the following parameters:

* `mu`: exogenous intensities (ndarray of shape $(d,)$),

* `alpha`: interaction matrix (ndarray of shape $(d, d)$),

* `beta`: common decay scalar (float),

* `end_time`: horizon time (float),

* `n Samples`: number of repeated paths (int),

* `random_state`: seed for random number generation (int or None).

Once this has been done, several methods of interest can be called. The main method is `simulate()`, which simulates the events of the considered Hawkes process. A brief list of the usable methods of this class is given in the following Table 2.
<table><caption>Table 2: Methods of the SimuHawkesExp class from the Sparklen package.</caption><thead><tr><td>Method</td><td>Description</td></tr></thead><tbody><tr><td>simulate()</td><td>Simulate process events.</td></tr><tr><td>compensator(t)</td><td>Compute the compensator evaluated in t.</td></tr><tr><td>spectral radius()</td><td>Return the spectral radius of the interaction matrix.</td></tr><tr><td>print_info()</td><td>Display key details about the object.</td></tr></tbody></table>

**Example 1: simulation of a MHP.** The following provides a comprehensive guide to using the SimuHawkesExp class to simulate the events of a MHP. In this example, we consider a 2-MHP whose parameters are given in Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Heatmap visualization of $\theta^* = (\mu^*, \alpha^*)$, where $\mu^*$ values are shown in red (gradient on the left colorbar) and $\alpha^*$ values are shown in blue (gradient on the right colorbar).</figcaption></figure>

It all starts with importation of the class SimuHawkesExp from the simulation sub-module.

```python
>>> from sparklen.hawkes.simulation import SimuHawkesExp
```

The next step is to instantiate an object of this class. At this stage, the parameters of the specified MHP, along with the desired time horizon and number of repetitions, must be provided. In this example, we set $\beta = 3.0$, $T = 5.0$ and, for simplicity, choose $n = 3$ repetitions. In addition, note that the seed is fixed, ensuring reproducible output across multiple runs of this example code.

```python
>>> hawkes = SimuHawkesExp(
...     mu=mu, alpha=alpha, beta=beta,
...     end_time=T, n Samples=n,
...     random_state=8)
```

With the object instantiated, we can invoke the simulate method to generate realizations.

```python
>>> hawkes.simulate()
```

Having done this, we can access the generated data via the timestamps property.

```python
>>> hawkes.timestamps
```
```python
[[array([0.53476798, 1.20676073, 1.86374452, 2.18940937, 2.39482727]),
   array([0.96932274, 1.155102 ])],
 [array([2.26133264, 3.08363972]),
   array([0.48054217, 0.70050504, 1.64644688, 2.26772207, 2.48926523])],
 [array([0.67427111, 2.29518797, 3.05520982]),
   array([1.86160981])]]
```

An illustration of these repeated paths is supplied in Figure 4.

<figure><img src="image_6.png" /><figcaption>Figure 4: Visualization of the simulated timestamps of the 2-MHP. On top, the repeated paths of the first component of the network; on the bottom, the repeated paths of the second component. The repetition of paths within each dimension is color-coded using a gradient of blue.</figcaption></figure>

### 3.3. Implementation of the inference framework

**Description.** The **Sparklen**'s inference tools are housed within the dedicated sub-module `sparklen.hawkes.inference`. This module relies on the `sparklen.optim`, `sparklen.prox`, and `sparklen.calibration` modules (see Figure 2) to effectively carry out the minimization of the problem given by Equation (2). The flagship `LearnerHawkesExp` class operates as a central unit, providing a unified interface that orchestrates the interplay between these three modules. This implementation directly mirrors the modular structure presented in Section 2.3, granting users access to the same extensive range of options. By abstracting the complexity of managing these interconnected modules, it ensures a smooth and seamless workflow. At the same time, the design maintains flexibility, allowing users to customize the estimation procedure through intuitive string-based arguments and configuration dictionaries. To sum up, the benefits of this approach are two-fold: it simplifies the user experience while ensuring the module remains adaptable to a wide variety of application needs. The following options are available:
* `loss`: 'least-squares' or 'log-likelihood';

* `penalty`: 'none', 'lasso', 'ridge' or 'elasticnet';

* `kappaChoice`: 'cv', 'bic' or 'ebic';

* `optimizer`: 'gd' or 'agd';

* `lr_scheduler`: 'lipschitz' or 'backtracking'.

Each object of this class is instantiated by specifying the configuration choices, along with additional arguments where applicable (e.g., the number of folds for cross-validation if 'cv' is chosen, or the hyperparameter $\gamma$ if the 'ebic' criterion is selected). After instantiation, the user must call the `fit()` method to fit the model to the data. Once the model has been fitted, several methods of interest can be invoked. A brief overview of the available methods for this class is provided in Table 3.

<table><caption>Table 3: Methods of the LearnerHawkesExp class from the Sparklen package</caption><thead><tr><td>Method</td><td>Description</td></tr></thead><tbody><tr><td>fit(data, end_time)</td><td>Fit Hawkes model to data.</td></tr><tr><td>score(data, end_time)</td><td>Compute score metric.</td></tr><tr><td>plot_estimated_values()</td><td>Plot estimated parameter values.</td></tr><tr><td>plot_estimated_support()</td><td>Plot estimated parameter support.</td></tr><tr><td>print_info()</td><td>Display key details about the object.</td></tr></tbody></table>

**Example 2: inference of a MHP.** The following provides a comprehensive guide to using the LearnerHawkesExp class to get an estimation of the true parameter $\theta^*$ of a MHP. In this example, we consider a 5-MHP whose parameter is given in Figure 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Heatmap visualization of $\theta^* = (\mu^*, \alpha^*)$, where $\mu^*$ values are shown in red (gradient on the left colorbar) and $\alpha^*$ values are shown in blue (gradient on the right colorbar).</figcaption></figure>

Before getting started, we need to import the SimuHawkesExp responsible for the simulation of synthetic data and LearnerHawkesExp to carry out the inference task.
>>> from sparklen.hawkes.simulation import SimuHawkesExp
>>> from sparklen.hawkes.inference import LearnerHawkesExp

The initial step involves simulating the repeated paths of the corresponding MHP, which will serve as the training set for subsequent inference. For this demonstration, we choose $\beta = 3.0$, we use a training sample of size $n = 1000$ and a time horizon of $T = 5.0$. Additional details on this step can be found in Example 3.2.

>>> hawkes = SimuHawkesExp(
...     mu=mu, alpha=alpha, beta=beta,
...     end_time=T, n Samples=n,
...     random_state=4)

>>> hawkes.simulate()

>>> data = hawkes.timestamps

The next step is to instantiate an object of this class. At this stage, we must give the decay parameter, and the chosen options of the learner. For readability, the `verbose_bar` argument is set to `'False'`. If set to `'True'`, a progress bar will be displayed to monitor the optimization progress. The `verbose` argument, on the other hand, controls whether recorded information during the optimization phase is printed in the end. If set to `'True'`, details will be displayed every `print_every` iterations.

>>> learner = LearnerHawkesExp(
...     decay=beta, loss="least-squares", penalty="none",
...     optimizer="agd", lr_scheduler="backtracking",
...     max_iter=200, tol=1e-5, verbose_bar=False,
...     verbose=True, print_every=10, record_every=10)

With the object instantiated, we can fit the model using the given training data. The output includes information recorded during optimization, the termination status, and the total duration of the procedure.

>>> learner.fit(data, T)

Optimization completed. Convergence achieved after 24 iterations.

Time elapsed: 0.40 seconds.

<table><thead><tr><th>Iteration</th><th>Loss</th><th>Tolerance</th></tr></thead><tbody><tr><td>0</td><td>-12.7803</td><td>0.0863339</td></tr><tr><td>10</td><td>-13.5298</td><td>0.000408739</td></tr><tr><td>20</td><td>-13.5624</td><td>2.35301e-05</td></tr></tbody></table>
Having done this, we can obtain the estimated parameters by calling the `estimated_params` property method.

>>> learner.estimated_params

array([0.56945931, 0.12918584, 0.09980257, 0.1054643, 0.10675509, 0.00803352], [0.55539019, 0.11320206, 0.09733585, 0.08572127, 0.10494816, 0.00640289], [0.63591031, 0.10224589, 0.08963574, 0.24939751, 0.24530404, 0.15960453], [0.56215877, 0.0958248, 0.08496727, 0.24182603, 0.24653997, 0.15607922], [0.59766823, 0.08496727, 0.13851235, 0.13871558, 0.1488716])

Additionally, for better visualization, we can plot the estimated values of the parameter by calling the method `plot_estimated_values()`, represented in Figure 6.

>>> learner.plot_estimated_values()

<figure><img src="image_8.png" /><figcaption>Figure 6: Heatmap visualization of $\hat{\theta} = (\hat{\mu}, \hat{\alpha})$, where $\mu^*$ values are shown in red (gradient on the left colorbar) and $\alpha^*$ values are shown in blue (gradient on the right colorbar).</figcaption></figure>

As shown in Figure 6, the estimation procedure provides a quality estimate of the true parameter $\theta^*$. The visualization demonstrates that the estimated parameter aligns well with the ground truth.

One can also get the score on training data by calling the `score()` method:

>>> learner.score(data, T)

**Example 3: inference of a sparse high-dimensional MHP.** Now, we shift our focus to a high-dimensional MHP where the interaction matrix exhibits a sparse structure. This example provides a detailed walkthrough for using the LearnerHawkesExp class, showcasing the specialized options designed for such cases. We consider a 25-dimensional MHP, with the parameter configuration shown in Figure 7.
<figure><img src="image_9.png" /><figcaption>Figure 7: Heatmap visualization of $\theta^* = (\mu^*, \alpha^*)$, where $\mu^*$ values are shown in red (gradient on the left colorbar) and $\alpha^*$ values are shown in blue (gradient on the right colorbar). The sparsity rate of the matrix $\alpha^*$, i.e. the proportion of zero coefficients, is equal to 0.88.</figcaption></figure>

First, as before, we need to import the classes and simulate the training dataset. For details on these steps, refer to the two examples above. We choose $\beta = 3.0$, $n = 250$ and $T = 5.0$. In this example, the high-dimensional ratio, given by $d/n$, is equal to 0.1

Next, we need to create an instance of the class `LearnerHawkesExp` and specify the options for this learner. In this context, we opt for Lasso regularization, and we compare the tuning of the penalization constant being done using EBIC with $\gamma = 1.0$ and using CV with a number of folds equal to 10.

```python
>>> learner_lasso_cv = LearnerHawkesExp(
...     decay=beta, loss="least-squares",
...     penalty="lasso", kappa.choice="cv"
...     optimizer="agd", lr_scheduler="backtracking",
...     max_iter=200, tol=1e-5,
...     verbose_bar=False, verbose=False,
...     cv=10)
```

```python
>>> learner_lasso_ebic = LearnerHawkesExp(
...     decay=beta, loss="least-squares",
...     penalty="lasso", kappa.choice="ebic"
...     optimizer="agd", lr_scheduler="backtracking",
...     max_iter=200, tol=1e-5,
...     verbose_bar=False, verbose=False,
...     gamma=1.0)
```

Then, we can fit the model using the provided training data.

```python
>>> learner_lasso_cv.fit(data, T)
```

```python
>>> learner_lasso_ebic.fit(data, T)
```
One can also get the score on training data by calling the `score()` method:

```python
>>> learner_lasso_cv.score(data, T)
```

```python
>>> learner_lasso_ebic.score(data, T)
```

-66.83

-66.29

In the context of sparse MHP in a high-dimensional setting, a key consideration is the ability to recover the true support of $\theta^*$. To demonstrate the efficacy of the Lasso estimator in interaction selection, we can visualize the estimated support by invoking the method `plot_estimated_support()`, represented in Figure 8.

```python
>>> learner_lasso_cv.plot_estimated_support()
```

```python
>>> learner_lasso_ebic.plot_estimated_support()
```

<figure><img src="image_10.png" /><figcaption>Figure 8: Estimated support with Lasso regularization.</figcaption></figure>

As it can be seen in Figure 8, the Lasso estimator, when using CV to select the penalization constant, tends to be overly conservative and insufficiently penalizes, failing to identify the true support. In contrast, tuning $\kappa$ using EBIC enables the Lasso estimator to successfully recover the underlying block structure of the interaction matrix.

## 3.4. Implementation of the classification framework

**Description.** The Sparklen's classification tools are housed within the dedicated sub-module `n.hawkes.classification`. The first element is the `makeClassification` function responsible for the generation of a random K-class classification problem. The arguments of this function are the following:
* bold Mu: exogenous intensities of each class (ndarray of shape (K, d, )),

* bold alpha: interaction matrix of each class (ndarray of shape (K, d, d)),

* beta: common decay scalar (float),

* end_time: horizon time (float),

* n Samples: the number of samples (int),

* n_classes: the number of classes (int),

* weights: proportions of samples assigned to each class (ndarray of shape (K, )),

* random_state: seed for random number generation (int or None).

The rest of the module consists of tools for performing classification. In particular, the two
classifier classes are the **ERMCLlassifier** and the **ERMLRCLlassifier** that implement the ERM
and ERMLR classifier presented in Section 2.4. The core implementation logic is shared
between the two classes and both are designed to follow the **scikit-learn** API (see Pedregosa
*et al.* (2011) and Buitinck *et al.* (2013)). A list of usable methods of these class is given in
Table 4.

<table><caption>Table 4: Usable methods shared by the ERMCLlassifier and ERMLRCLassifier class from the Sparklen package.</caption>
   <thead>
    <tr>
     <td>
      Method
     </td>
     <td>
      Description
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      fit(X, y, end_time)
     </td>
     <td>
      Build a classifier from the training set (X, y).
     </td>
    </tr>
    <tr>
     <td>
      predict(X, end_time)
     </td>
     <td>
      Predict class for X.
     </td>
    </tr>
    <tr>
     <td>
      predict_proba(X, end_time)
     </td>
     <td>
      Predict class probabilities for X.
     </td>
    </tr>
    <tr>
     <td>
      score(X, y, end_time)
     </td>
     <td>
      Return the mean accuracy on the given test set (X, y).
     </td>
    </tr>
    <tr>
     <td>
      print_info()
     </td>
     <td>
      Display key details about the object.
     </td>
    </tr>
   </tbody>
  </table>

**Example 4: classification of a MHP.** To illustrate the performance of both classifiers, we consider the following classification model with *K* = 3 of a 15-MHP. We represent the parameter associated to *k* = 2 in Figure 9. The parameters of the two other classes are obtained by interchanging the blocks of the interaction matrix. All three classes have the same exogenous intensity vector. Finally, we choose *β* = 3.0, *n* = 600, *T* = 5.0 and the seed is fixed.

To begin with, the function and class need to be imported.

>>> from sparklen.hawkes.classification import make_classification
>>> from sparklen.hawkes.classification import ERMCLassifier, ERMLRCLassifier

The first step consists of generating the training and testing datasets. This can be achieved
using the makeClassification() method, followed by splitting the data evenly, with 50%
allocated for training and the remaining 50% for testing.
<figure><img src="image_11.png" /><figcaption>Figure 9: Heatmap visualization of $\theta^* = (\mu^*, \alpha^*)$, where $\mu^*$ values are shown in red (gradient on the left colorbar) and $\alpha^*$ values are shown in blue (gradient on the right colorbar).</figcaption></figure>

```python
>>> X, y = make_classification(
...     bold Mu=bold Mu, bold alpha=bold alpha,
...     beta=beta, end_time=T,
...     n Samples=n, n_classes=K
...     random_state=4)
>>> X_train, X_test, y_train, y_test = train_test_split(
...     X, y, test_size=0.5, random_state=4)
```

Then, we need to instantiate the classifier-type objects and to specify the chosen arguments.

```python
>>> conf_erm = ERMLClassifier(
...     decay=beta, gamma0=0.1,
...     max_iter=500, tol=1e-6)
>>> conf_ermlr = ERMLRCLassifier(
...     decay=beta, gamma0=0.1,
...     max_iter=500, tol=1e-6)
```

Finally we can fit the model according to the given training set (X, y).

```python
>>> conf_erm.fit(X=X_train, y=y_train, end_time=T)
>>> conf_ermlr.fit(X=X_train, y=y_train, end_time=T)
```

Once this step is completed, the user gains access to a variety of useful methods. For instance, predictive probabilities can be obtained by calling the `predict_proba()` method.

```python
>>> conf_erm Predict_proba(X=X_test, end_time=T)
>>> conf_ermlr Predict_proba(X=X_test, end_time=T)
```
Additionally, one can also get the prediction of the classifier by calling the `predict()`.

>>> coef_erm Predicting (X=X_test, end_time=T)

>>> coef_ermlr Predicting (X=X_test, end_time=T)

Furthermore, the classifier's performance relative to the provided test labels can be evaluated by calling the `score()` method.

>>> coef_erm.score(X=X_test, y=y_test, end_time=T)

>>> coef_ermlr.score(X=X_test, y=y_test, end_time=T)

0.65

0.88

Notably, it can be seen that classifier ERMLR performs better than classifier ERM, validating the theoretical results.

Finally, the user can also plot the confusion matrix of the prediction with respect to the true labels by calling the `plot_score_cm()` method. The two confusion matrices associated with each classifier are presented in Figure 10.

>>> coef_erm.plot_score_cm(X=X_test, y=Y_test, end_time=T)

>>> coef_ermlr.plot_score_cm(X=X_test, y=Y_test, end_time=T)

<figure><img src="image_12.png" /><figcaption>Figure 10: Confusion matrices illustrating performance for the two classifiers. The diagonal shows the accuracy per class, while off-diagonal elements indicate misclassified instances.</figcaption></figure>
## 4. Illustration on MemeTracker dataset

In this section, we illustrate the **Sparklen** package using the MemeTracker dataset (see Leskovec, Backstrom, and Kleinberg (2009)). MemeTracker is a benchmark corpus containing blog posts published between August 2008 and April 2009. This dataset tracks frequently quoted phrases and memes, allowing for the analysis of information propagation across digital media. It is widely used in computational social science, natural language processing, and network analysis to study how narratives emerge, evolve, and spread across the internet.

**Pre-processing.** We use posts from August 2008 to December 2008, covering a period of 153 days. An event for a website is defined as the creation of a post on that website containing a hyperlink to any other website. Additionally, we consider that an edge exists between two websites if at least one hyperlink connects them in the training set. To focus on the most influential sources, we extract the top 100 media sites with the highest number of documents. For each of these 100 websites, we retrieve their corresponding events by using the code made available by Achab, Bacry, Gaïffas, Mastromatteo, and Muzy (2017). Finally, we segment the process by grouping events by day, resulting in 153 repetitions of events over the given period. Finally, it is important to note that an approximation of the interaction matrix for the network can be computed by leveraging hyperlinks to trace the flow of information (see Achab et al. (2017)). These causal relationships relies on the parenthood relations given by the clustering representation (see Bacry et al. (2015)).

**Benchmark.** For comparison purposes, we consider a benchmark approach, specifically ADM4 (see Zhou, Zha, and Song (2013)). Both our procedure and ADM4 operate with a given fixed input $\beta$. Since this hyperparameter is unknown, we explore a grid of decay values and select the one that maximizes the score, i.e., the loss value under consideration. Given the sparse nature of the ground truth interaction matrix, both procedures are coupled with Lasso regularization. For the **Sparklen** procedure, we use the least-squares loss for its ability to effectively handle a large number of events and the EBIC criterion to select the penalty constant. In the case of ADM4, no automatic tuning of the penalty constant is provided, so we explore a grid of regularization constants and select the one that maximizes the score.

**Metrics.** Let $\alpha^*$ be the ground truth interaction matrix and $\hat{\alpha}$ an estimated candidate. We consider the following three metrics. First, the Hamming Distance:

$$ \text{HammDist}(\alpha^*, \hat{\alpha}) := \frac{1}{d^2} \sum_{j,j'=1}^{d} \mathbb{1}_{\{\alpha_{j,j'}^* \neq \hat{\alpha}_{j,j'}\}} $$

which allows to assess the quality of the support recovery. Then, the Relative Error:

$$ \text{RelErr}(\alpha^*, \hat{\alpha}) := \frac{1}{d^2} \sum_{j,j'=1}^{d} \frac{|\alpha_{j,j'}^* - \hat{\alpha}_{j,j'}|}{|\alpha_{j,j'}^*|} \mathbb{1}_{\{\alpha_{j,j'}^* \neq 0\}} + |\hat{\alpha}_{j,j'}| \mathbb{1}_{\{\alpha_{j,j'}^* = 0\}} $$

which quantifies the averaged relative error between the true $\alpha^*$ and the estimated $\hat{\alpha}$ on non-zero entries of $\alpha^*$ and error on zeros entries. Finally, the mean Kendall rank Correlation:

$$ \text{RankCorr}(\alpha^*, \hat{\alpha}) := \frac{1}{d} \sum_{j=1}^{d} \tau(\alpha_{j,\cdot}^*, \hat{\alpha}_{j,\cdot}) $$
which corresponds to the averaged Kendall's rank correlation coefficient between each row of the true $\alpha^*$ and that of the estimated $\hat{\alpha}$. While relevant, we wish to highlight that these three metrics have limitations in terms of reliability when used to compare against a sparse matrix, as is the case for the ground truth here. For example, the null matrix has a Hamming Distance and a Relative Error of 0.06 compared to the ground truth matrix.

**Evaluation scheme.** We present the evaluation scheme that relies on Monte-Carlo repetitions. We independently repeat the following steps 5 times. Following the approach in Zhou et al. (2013), the dataset is split evenly, with 50% allocated for training and the remaining 50% for testing. On the training data, we run both procedures by iterating through a grid of beta values and retrieve the estimate $\hat{\alpha}$ associated with the optimal decay parameter. Once this is done, we compute the metrics of the estimate provided by the two procedures. Finally, on the test dataset, we compute the negative log-likelihood of the model based on the estimated parameters.

**Results.** The obtained results are given in Table 5.

<table><caption>Table 5: Averaged values of the metrics over the 5 repetitions for both procedures, with the standard deviation provided in parentheses.</caption><thead><tr><th>Metric</th><th>HammDist</th><th>RelErr</th><th>RankCorr</th><th>NegLikTest</th></tr></thead><tbody><tr><td>ADM4</td><td>0.80 (0.01)</td><td>75.56 (2.60)</td><td>0.136 (0.003)</td><td><strong>1.94e+07</strong> (4.32e+05)</td></tr><tr><td>Sparklen</td><td><strong>0.06</strong> (0.00)</td><td><strong>37.68</strong> (0.20)</td><td><strong>0.570</strong> (0.008)</td><td>2.20e+07 (5.42e+05)</td></tr></tbody></table>

As observed, our procedure yields a higher negative log-likelihood on the test sample compared to ADM4. This can be attributed to the fact that our method optimizes the least-squares loss, whereas ADM4 is based on log-likelihood. On the other hand, it can be seen that our procedure outperforms ADM4 across all three metrics. In particular, the estimated interaction matrix exhibits a very low Hamming Distance, suggesting that our method effectively recovers the support of the ground truth $\alpha^*$. This indicates that it captures the influence network more accurately.

# 5. Concluding remarks and future perspectives

**Conclusion.** Sparklen provides a collection of statistical learning tools designed to work with MHP in Python. In particular, it implements an efficient sampling procedure, a highly versatile and flexible framework for performing inference, and, finally, novel approaches to address the challenge of multiclass classification within the supervised learning framework. In this article, we have presented these tools and their use cases through instructive examples. Notably, we explored the potential and effectiveness of the proposed procedures using simulated data. Additionally, we demonstrated the inference framework by recovering a social network graph from the MemeTracker dataset. Together with its statistical contributions, its quality implementation, which ensures both accessibility and efficiency, makes it a highly
appealing software solution. We strongly believe that this package can greatly benefit users seeking to employ MHP as a modeling tool across various applications.

**Perspectives.** At the heart of this project is the ambition to provide the community with a cohesive environment of cutting-edge, user-friendly tools for manipulating MHP. While this package already establishes a solid foundation, plans for future developments are envisioned to further enhance and expand its scope of usability and accessibility. Let us discuss few of them.

A natural first step in this incremental goal would be to incorporate additional kernel func-
tions. While exponential kernels are widely used and cover a broad range of applications,
introducing more complex options, such as sums of exponential kernels (see Lemonnier and
Vayatis (2014)), would further expand the scope of applicability. From a implementation
standpoint, this addition would present minimal challenges, as the library is designed in a
modular fashion, enabling model classes to be seamlessly integrated into the inference machin-
ery. The only requirement would be to implement the actual model classes, thus encouraging
such an addition.

This article focuses on linear Hawkes processes, but incorporating tools to work with non-linear ones could be both interesting and highly motivating from an application standpoint. Allowing inhibition-type interactions would enhance modeling possibilities and extend applicability to a broader range of cases. For example, Bonnet, Martinez Herrera, and Sangnier (2023) introduces a method to estimate MHP processes with an exponential kernel in the presence of both excitation and inhibition effects.

Finally, incorporating a multi-threading option to leverage the computational capabilities of
modern processors could significantly enhance efficiency in computationally demanding tasks.

Computational details

The results in this paper were obtained using Python 3.11.10 with the **Sparklen** 0.1.0 package
on a CentOS Linux release 7.7.1908 operating system. The computations were run on a single
node with an Intel Xeon Gold 6230 CPU using one core and 16 GB of RAM.

Reproducibility and documentation

All results and figures can be reproduced by running the scripts available in the GitHub
repository https://github.com/romain-e-lacoste/sparklen. The repository includes all
necessary code, along with instructions for setting up the required environment and dependen-
cies. Additionally, to facilitate the use of **Sparklen**, the code is meticulously commented. Full
documentation, taking advantage of these comments with the Python tool **Sphinx** (**Sphinx de-
velopment team (2025)**), is currently under construction. The **Sparklen** package is distributed
under the OSI-approved 3-Clause BSD License.
Acknowledgments

This work has been supported by the Chaire "Modélisation Mathématique et Biodiversité" of Veolia-École polytechnique-Museum national d'Histoire naturelle-Fondation X, through a Ph.D. scholarship. The computations have been performed under the project "hawkes" on the HPC facility Cholesky operated by the IDCS / École polytechnique. This work is also part of the 2022 DAE 103 EMERGENCE(S) - PROCECO project supported by Ville de Paris. Finally, I am grateful to my advisors, Christophe Denis, Charlotte Dion-Blanc, and Laure Sansonnet, for their sound advice and their review work throughout the writing of this article.

References

Achab M, Bacry E, Gaïffas S, Mastromatteo I, Muzy JF (2017). "Uncovering Causality from Multivariate Hawkes Integrated Cumulants." In *Proceedings of the 34th International Conference on Machine Learning*, volume 70 of *Proceedings of Machine Learning Research*, pp. 1-10. PMLR. URL https://proceedings.mlr.press/v70/achab17a.html.

Arlot S, Celisse A (2010). "A survey of cross-validation procedures for model selection." *Statistics Surveys*, 4(none), 40-79. doi:10.1214/09-SS054.

Armijo L (1966). "Minimization of functions having Lipschitz continuous first partial derivatives." *Pacific Journal of mathematics*, **16**(1), 1-3.

Bach F (2024). *Learning theory from first principles*. MIT press. URL https://mitpress.mit.edu/9780262049443/learning-theory-from-first-principles/.

Bacry E, Bompaire M, Deegan P, Gaïffas S, Poulsen SV (2018). "tick: a Python library for statistical learning, with an emphasis on Hawkes processes and time-dependent models." *Journal of Machine Learning Research*, **18**(214), 1-5. URL https://www.jmlr.org/papers/v18/17-381.html.

Bacry E, Bompaire M, Gaïffas S, Muzy JF (2020). "Sparse and low-rank multivariate Hawkes processes." *Journal of Machine Learning Research*, **21**(50), 1-32. URL http://jmlr.org/papers/v21/15-114.html.

Bacry E, Mastromatteo I, Muzy J (2015). "Hawkes Processes in Finance." *Market Microstructure and Liquidity*, **01**(01), 1550005. doi:10.1142/S2382626615500057.

Bartlett PL, Jordan MI, McAuliffe JD (2006). "Convexity, Classification, and Risk Bounds." *Journal of the American Statistical Association*, **101**(473), 138-156. doi:10.1198/016214505000000907.

Beazley DM (2003). "Automated scientific software scripting with SWIG." *Future Generation Computer Systems*, **19**(5), 599-609.

Beck A, Teboulle M (2009). "A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems." *SIAM Journal on Imaging Sciences*, **2**(1), 183-202. doi:10.1137/080716542.
Bonnet A, Dion-Blanc C, Gindraud F, Lemler S (2022). “Neuronal network inference and membrane potential model using multivariate Hawkes processes.” *Journal of Neuroscience Methods*, **372**, 109550. ISSN 0165-0270. doi:https://doi.org/10.1016/j.jneumeth.2022.109550.

Bonnet A, Martinez Herrera M, Sangnier M (2023). “Inference of multivariate exponential Hawkes processes with inhibition and application to neuronal activity.” *Statistics and Computing*, **33**(4), 91.

Buitinck L, Louppe G, Blondel M (2013). “API design for machine learning software: experiences from the scikit-learn project.” In *ECML/PKDD 2013 Workshop: Languages for Data Mining and Machine Learning*.

Chaovalitwongse WA, Prokopyev OA, Pardalos PM (2006). “Electroencephalogram (EEG) time series classification: Applications in epilepsy.” *Annals of Operations Research*, **148**(1), 227–250.

Chen J, Chen Z (2008). “Extended Bayesian Information Criteria for Model Selection with Large Model Spaces.” *Biometrika*, **95**(3), 759–771. doi:10.1093/biomet/asn034.

Chen X, Wang X (2020). “Perfect Sampling of Multivariate Hawkes Processes.” In *2020 Winter Simulation Conference (WSC)*, pp. 469–480. doi:10.1109/WSC48552.2020.9384038.

Cheysson F (2024). *Estimation of Hawkes Processes from Binned Observations*. R package version 1.0.3, URL https://cran.r-project.org/web/packages/hawkesbow/hawkesbow.pdf.

Chzhen E, Giraud C, Stoltz G (2023). “Parameter-free projected gradient descent.” *arXiv preprint arXiv:2305.19605*.

Chzhen E, Hebiri M, Salmon J (2019). “On Lasso refitting strategies.” *Bernoulli*, **25**(4A), 3175–3200. doi:10.3150/18-BEJ1085.

Denis C, Dion-Blanc C, Lacoste RE, Sansonnet L, Bas Y (2024). “Bats Monitoring: A Classification Procedure of Bats Behaviors based on Hawkes Processes.” *Journal of the Royal Statistical Society Series C: Applied Statistics*, **73**(4), 1025–1041. doi:10.1093/jrsssc/qlae024.

Denis C, Dion-Blanc C, Sansonnet L (2022). “Multiclass classification for Hawkes processes.” In *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, volume 180, pp. 539–547. PMLR. URL https://proceedings.mlr.press/v180/denis22a.html.

Dion-Blanc C, Denis C, Sansonnet L, Lacoste RE (2024). “ERM-Lasso classification algorithm for Multivariate Hawkes Processes paths.” *arXiv preprint arXiv:2407.11455*.

Duchi J, Hazan E, Singer Y (2011). “Adaptive Subgradient Methods for Online Learning and Stochastic Optimization.” *Journal of Machine Learning Research*, **12**(61), 2121–2159. URL http://jmlr.org/papers/v12/duchi11a.html.
Gupta MR, Bengio S, Weston J (2014). "Training highly multiclass classifiers." *The Journal of Machine Learning Research*, **15**(1), 1461–1492. URL https://www.jmlr.org/papers/v15/gupta14a.html.

Hawkes AG (1971). "Spectra of some self-exciting and mutually exciting point processes." *Biometrika*, **58**(1), 83–90. doi:10.1093/biomet/58.1.83.

Hawkes AG, Oakes D (1974). "A cluster process representation of a self-exciting process." *Journal of Applied Probability*, **11**(3), 493–503. doi:10.2307/3212693.

Hoerl AE, Kennard RW (1970). "Ridge Regression: Biased Estimation for Nonorthogonal Problems." *Technometrics*, **12**(1), 55–67. doi:10.1080/00401706.1970.10488634.

Ismail Fawaz H, Forestier G, Weber J, Idoumghar L, Muller PA (2019). "Deep learning for time series classification: a review." *Data mining and knowledge discovery*, **33**(4), 917–963. doi:10.1007/s10618-019-00619-1.

Lacoste RE (2025). **Sparklen: A Statistical Learning Toolkit for High-Dimensional Hawkes Process in Python**. Python package version 1.0.0, URL https://github.com/romain-e-lacoste/sparklen.

Lemonnier R, Vayatis N (2014). "Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate hawkes processes." In *Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2014, Nancy, France, September 15-19, 2014. Proceedings, Part II* 14, pp. 161–176. Springer.

Leskovec J, Backstrom L, Kleinberg J (2009). "Meme-tracking and the dynamics of the news cycle." In *Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, p. 497–506. ACM. doi:10.1145/1557019.1557077.

Linderman SW (2017). **pyhawkes: Bayesian inference for network Hawkes processes**. Python package version 0.3.1, URL http://www.github.com/slinderman/pyhawkes.

Lotz A (2024). "A sparsity test for multivariate Hawkes processes." *arXiv preprint arXiv:2405.08640*.

Lukasik M, Srijith P, Vu D, Bontcheva K, Zubiaga A, Cohn T (2016). "Hawkes processes for continuous time sequence classification: an application to rumour stance classification in twitter." In *Proceedings of 54th Annual Meeting of the Association for Computational Linguistics*, pp. 393–398. Association for Computational Linguistics.

Massart P, Nédélec É (2006). "Risk bounds for statistical learning." *The Annals of Statistics*, **34**(5), 2326–2366. doi:10.1214/009053606000000786.

Michael W Ferry Philip E Gill EW, Zhang M (2023). "A class of projected-search methods for bound-constrained optimization." *Optimization Methods and Software*, **0**(0), 1–30. doi:10.1080/10556788.2023.2241769.

Møller J, Rasmussen JG (2005). "Perfect simulation of Hawkes processes." *Advances in Applied Probability*, **37**(3), 629–646. doi:10.1239/aap/1127483739.
Nesterov Y (1983). "A method for solving the convex programming problem with convergence rate O (1/k2)." In Doklady Akademii Nauk SSSR, volume 269, p. 543.

Ogata Y (1978). "The asymptotic behaviour of maximum likelihood estimators for stationary point processes." *Annals of the Institute of Statistical Mathematics*, **30**(2), 243-261. doi: **10.1007/BF02480216**.

Ogata Y (1981). "On Lewis' simulation method for point processes." *IEEE Transactions on Information Theory*, **27**(1), 23-31. doi: **10.1109/TIT.1981.1056305**.

Ogata Y (1988). "Statistical models for earthquake occurrences and residual analysis for point processes." *Journal of the American Statistical association*, **83**(401), 9-27.

Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Édouard Duchesnay (2011). "Scikit-learn: Machine Learning in Python." *Journal of Machine Learning Research*, **12**(85), 2825-2830. URL http://jmlr.org/papers/v12/pedregosa11a.html.

Pinto JCL, Chahed T, Altman E (2015). "Trend detection in social networks using Hawkes processes." In *Proceedings of the 2015 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining 2015*, p. 1441-1448. doi: **10.1145/2808797.2814178**.

**Sphinx** development team (2025). *The Sphinx documentation generator*. Python package version 8.2.1, URL https://github.com/sphinx-doc/sphinx.

Reynaud-Bouret P, Rivoirard V, Grammont F, Tuleau-Malot C (2014). "Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis." *The Journal of Mathematical Neuroscience*, **4**, 1-41. doi: **10.1186/2190-8567-4-3**.

Reynaud-Bouret P, Rivoirard V, Tuleau-Malot C (2013). "Inference of functional connectivity in Neurosciences via Hawkes processes." In *2013 IEEE Global Conference on Signal and Information Processing*, pp. 317-320. doi: **10.1109/GlobalSIP.2013.6736879**.

Staerman G, Allain C, Gramfort A, Moreau T (2023). "Fadin: Fast discretized inference for hawkes processes with general parametric kernels." In *International Conference on Machine Learning*, pp. 32575-32597. PMLR.

Susto GA, Cenedese A, Terzi M (2018). "Time-series classification methods: Review and applications to power systems data." *Big data application in power systems*, pp. 179-220.

Tibshirani R (1996). "Regression shrinkage and selection via the lasso." *Journal of the Royal Statistical Society Series B: Statistical Methodology*, **58**(1), 267-288. doi: **10.1111/j.2517-6161.1996.tb02080.x**.

Vapnik V (1991). "Principles of Risk Minimization for Learning Theory." In *Advances in Neural Information Processing Systems*, volume 4. Morgan-Kaufmann. URL https://proceedings.neurips.cc/paper_files/paper/1991/file/ff4d5fbbafdf976cfdc032e3bde78de5-Paper.pdf.
Wang B, Zhang H, Ma Z, Chen W (2023). “Convergence of adagrad for non-convex objectives: Simple proofs and relaxed assumptions.” In *The Thirty Sixth Annual Conference on Learning Theory*, pp. 161–190. PMLR. URL https://proceedings.mlr.press/v195/wang23a.html.

Ward R, Wu X, Bottou L (2020). “Adagrad stepsizes: Sharp convergence over nonconvex landscapes.” *Journal of Machine Learning Research*, **21**(219), 1–30.

Zhang T (2004). “Statistical behavior and consistency of classification methods based on convex risk minimization.” *The Annals of Statistics*, **32**(1), 56–85. doi:10.1214/aos/1079120130.

Zhou K, Zha H, Song L (2013). “Learning Social Infectivity in Sparse Low-rank Networks Using Multi-dimensional Hawkes Processes.” In *Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics*, volume 31 of *Proceedings of Machine Learning Research*, pp. 641–649. PMLR. URL https://proceedings.mlr.press/v31/zhou13a.html.

Zou H, Hastie T (2005). “Regularization and Variable Selection Via the Elastic Net.” *Journal of the Royal Statistical Society Series B: Statistical Methodology*, **67**(2), 301–320. doi:10.1111/j.1467-9868.2005.00503.x.
## A. Details on pre-calculation

To avoid overloading the notations, we consider the version of functional that is not averaged over multiple repetitions but rather calculated from a single trajectory.

**Log-likelihood.** The negative log-likelihood loss can be written as follows:

$$L_T(\theta) = \frac{1}{T} \sum_{j=1}^{d} L_{j,T}(\theta)$$

Developing, for any $j \in [d]$, we have:

$$L_{j,T}(\theta) = \int_{0}^{T} \left( \mu_j + \sum_{j'=1}^{d} \sum_{h:t_{j',h}<t} \alpha_{j,j'} \beta e^{-\beta(t-t_{j',h})} \right) dt \\ - \sum_{\ell:t_{j,\ell}<T} \log \left( \mu_j + \sum_{j'=1}^{d} \sum_{h:t_{j',h}<t_{j,\ell}} \alpha_{j,j'} \beta e^{-\beta(t_{j,\ell}-t_{j',h})} \right) \\ = \mu_j T + \sum_{j'=1}^{d} \alpha_{j,j'} \sum_{h:t_{j',h}<T} \left(1 - e^{-\beta(T-t_{j',h})}\right) \\ - \sum_{\ell:t_{j,\ell}<T} \log \left( \mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \sum_{h:t_{j',h}<t_{j,\ell}} \beta e^{-\beta(t_{j,\ell}-t_{j',h})} \right) \\ = \mu_j T + \sum_{j'=1}^{d} \mathcal{I}_{j'} - \sum_{\ell:t_{j,\ell}<T} \log \left( \mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \Psi_{j'}(t_{j,\ell}) \right)$$

For the gradient, it follows that for $j \in [d]$, we have:

$$\frac{\partial L_{j,T}(\theta)}{\partial \mu_j} = T - \sum_{\ell:t_{j,\ell}<T} \frac{1}{\mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \Psi_{j'}(t_{j,\ell})}.$$

In addition, for $j, j' \in [d]$, we have:

$$\frac{\partial L_{j,T}(\theta)}{\partial \alpha_{j,j'}} = \mathcal{I}_{j'} - \sum_{\ell:t_{j,\ell}<T} \frac{\Psi_{j'}(t_{j,\ell})}{\mu_j + \sum_{j'=1}^{d} \alpha_{j,j'} \Psi_{j'}(t_{j,\ell})}.$$

**Least-squares.** The least-squares loss can be written as follows:

$$R_T(\theta) = \frac{1}{T} \sum_{j=1}^{d} R_{j,T}(\theta)$$
Developing, for any $j \in [d]$, we have:

$$
\begin{align*}
R_{j,T}(\theta) &= \int_0^T \left( \mu_j + \sum_{j'=1}^d \sum_{h:t_{j',h}<t} \alpha_{j,j'} \beta e^{-\beta(t-t_{j',h})} \right)^2 dt \\
&\quad - 2 \sum_{\ell:t_{j,\ell}<T} \left( \mu_j + \sum_{j'=1}^d \sum_{h:t_{j',h}<t_{j,\ell}} \alpha_{j,j'} \beta e^{-\beta(t_{j,\ell}-t_{j',h})} \right) \\
&= \mu_j^2 + \sum_{j'=1}^d \alpha_{j,j'}^2 \sum_{h:t_{j',h}<T} \frac{\beta}{2} (1 - e^{-2\beta(T-t_{j',h})}) \\
&\quad + \sum_{j'=1}^d \sum_{j''=1}^d \alpha_{j,j'} \alpha_{j,j''} \sum_{h:t_{j',h}<T} (1 - e^{-2\beta(T-t_{j',h})}) \sum_{h':t_{j''},h'<t_{j',h}} \beta e^{-\beta(t_{j',h}-t_{j''},h')} \\
&\quad + 2\mu_j \sum_{j'=1}^d \alpha_{j,j'} \sum_{h:t_{j',h}<T} (1 - e^{-\beta(T-t_{j',h})}) \\
&\quad - 2\mu_j N_j(T) - 2 \sum_{j'=1}^d \alpha_{j,j'} \sum_{\ell:t_{j,\ell}<T} \sum_{h:t_{j',h}<t_{j,\ell}} \beta e^{-\beta(t_{j,\ell}-t_{j',h})} \\
&= \mu_j^2 + \sum_{j'=1}^d \alpha_{j,j'}^2 \mathcal{I}_{j'}^2 + \sum_{j'=1}^d \sum_{j''=1}^d \alpha_{j,j'} \alpha_{j,j''} \mathcal{W}_{j',j''} \\
&\quad + 2\mu_j \sum_{j'=1}^d \alpha_{j,j'} \mathcal{I}_{j'} - 2\mu_j N_j(T) - 2 \sum_{j'=1}^d \alpha_{j,j'} \mathcal{V}_{j,j'}
\end{align*}
$$

For the gradient, it follows that for $j \in [d]$, we have:

$$
\frac{\partial R_{j,T}(\theta)}{\partial \mu_j} = 2\mu_j T + 2 \sum_{j'=1}^{d} \alpha_{j,j'} H_{j'} - 2N_j(T)
$$

In addition, for $j, j' \in [d]$, we have:

$$
\frac{\partial R_{j,T}(\theta)}{\partial \alpha_{j,j'}} = 2\alpha_{j,j'} \mathcal{I}_{j'}^2 + \sum_{\substack{j''=1 \\ j'' \neq j'}}^{d} \alpha_{j,j''} \mathcal{W}_{j',j''} + 2\mu_j \mathcal{I}_{j'} - 2\mathcal{V}_{j,j'}
$$

Notably, such factorization unveils a remarkable computational property of the least-squares
functional. Indeed, pre-calculation performed once in advance removes the complexity's de-
pendency on the number of events for both loss and gradient evaluation. This makes the
loss function highly advantageous in studies involving networks with a very large number of
observed events.

B. Details on SWIG
**C++/Python binding architecture.** To make the C++ class usable from Python, it is wrapped in Python using **SWIG**. The proxy classes generated automatically by **SWIG** are placed in a build folder to indicate that, although they can be imported into Python, they are not intended for direct user interaction. Instead, a corresponding Python class is created to act as a shell around these **SWIG** proxy classes. This Python class contains an attribute referencing the proxy class, with each method being rewritten to include error management and to call the proxy method through the attribute for computation. This three-layer design allows error handling to be implemented and the classes to be given a fully Python interface leveraging **SWIG** generated code while benefiting from the high-level performance of C++.

**Internal interoperability mechanism.** The integration between C++ and Python in the **Sparklen** package is underpinned by a sophisticated mechanism designed to enable seamless data sharing between the two languages. At its heart lies the **SharedArray** class, a highly adaptable array container in C++ that leverages the `std::shared_ptr` smart pointer for efficient memory management. This implementation ensures automatic reference tracking, safeguarding against memory leaks.

Then, this interoperability is made possible by the utilization of **SWIG** typemaps, specifically `typemap(in)` and `typemap(out)`. These typemaps are crafted to facilitate data conversion without incurring additional memory allocation, while carefully preserving ownership semantics. For instance, when a Python `NumPy` array is passed to C++, the C++ `SharedArray` refrains from assuming ownership, instead accessing the memory directly via a pointer. This approach eliminates redundant copying and maximizes efficiency since no new memory allocation is done.

Consider the following example, which shows how a 2D `NumPy` array is converted into a `SharedArrayDouble2D` object:

```c
%typemap(in) (SharedArrayDouble2D &) (SharedArrayDouble2D res) {
    if (!PyArray_Check($input)) {
        PyErr_SetString(PyExc_TypeError,
                        "The input argument should be a NumPy Array");
        SWIG_fail;
    }
    if (PyArray_TYPE((PyArrayObject*)$input) != NPY_DOUBLE) {
        PyErr_SetString(PyExc_TypeError,
                        "The data type of the NumPy Array should be double");
        SWIG_fail;
    }
    if (PyArray_NDIM((PyArrayObject*)$input) != 2) {
        PyErr_SetString(PyExc_ValueError,
                        "The Expected Numpy Array should be 2-dimensional");
        SWIG_fail;
    }
    // Collect NumPy Array features
    npy_intp n_rows = PyArray DIM((PyArrayObject*)$input, 0);
    npy_intp n_COLS = PyArray DIM((PyArrayObject*)$input, 1);
    double *data = staticCast<double>(PyArray_DATA((PyArrayObject*)$input));
```
```c
// Instantiate SharedArray2D from NumPy array features
res = SharedArrayDouble2D(n_rows, n_cols, data);
res.setPythonOwner(true);
$1 = &res;
}
```

The NumPy array's dimensions and data pointer are extracted, enabling the creation of a SharedArrayDouble2D object directly linked to the same memory as the NumPy array. Ownership is carefully managed by flagging the Python ownership attribute, ensuring memory efficiency and conflict-free operations.

## Affiliation:

Romain E. Lacoste
Université Gustave Eiffel
5 Boulevard Descartes, 77420, Champs-sur-Marne, France
E-mail: romain.lacoste@polytechnique.edu
URL: https://romain-e-lacoste.github.io/