# Joint reproduction number and spatial connectivity structure estimation via graph sparsity-promoting penalized functional

Etienne Lasalle, Barbara Pascal

► To cite this version:

Etienne Lasalle, Barbara Pascal. Joint reproduction number and spatial connectivity structure estimation via graph sparsity-promoting penalized functional. 2025. hal-05271255

HAL Id: hal-05271255

https://hal.science/hal-05271255v1

Preprint submitted on 21 Sep 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JOINT REPRODUCTION NUMBER AND SPATIAL CONNECTIVITY STRUCTURE ESTIMATION VIA GRAPH SPARSITY-PROMOTING PENALIZED FUNCTIONAL

E. Lasalle¹, B. Pascal¹

¹ Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France, firstname.lastname@ls2n.fr

## ABSTRACT

During an epidemic outbreak, decision makers crucially need accurate and robust tools to monitor the pathogen propagation. The effective reproduction number, defined as the expected number of secondary infections stemming from one contaminated individual, is a state-of-the-art indicator quantifying the epidemic intensity. Numerous estimators have been developed to precisely track the reproduction number temporal evolution. Yet, COVID-19 pandemic surveillance raised unprecedented challenges due to the poor quality of worldwide reported infection counts. When monitoring the epidemic in different territories simultaneously, leveraging the spatial structure of data significantly enhances both the accuracy and robustness of reproduction number estimates. However, this requires a good estimate of the spatial structure. To tackle this major limitation, the present work proposes a joint estimator of the reproduction number and connectivity structure. The procedure is assessed through intensive numerical simulations on carefully designed synthetic data and illustrated on real COVID-19 spatiotemporal infection counts.

***Index Terms***— Regularized estimation, Spatiotemporal signals, Graph learning, Alternating optimization, Epidemiology

## 1. INTRODUCTION

**Context.** Accurate and efficient monitoring tools are of key importance during epidemic outbreaks. The effective reproduction number, a widely used indicator, captures the daily evolution of the intensity of an epidemic by quantifying the secondary infections stemming from a typical infectious individual [1–5]. A major challenge in real-time pandemic surveillance lies in the potentially very low quality of reported infection counts, showing pseudo-seasonalities, underestimation on days-off, delays [6, 7]. State-of-the-art reproduction number estimators, although enforcing a mild temporal regularity of the estimate [3, 5], are not robust enough to handle such low data quality. This motivated the design of more advanced variational strategies using sparsity-based regularizations to yield more realistic and readable estimates [6, 7]. Furthermore, as pandemics are simultaneously monitored in a large number of territories (e.g., different countries, counties of a state, states of a federal country), extensions of these variational estimators took advantage of the spatial structure of available data to enhance the accuracy of the reproduction number estimates via a graph-based spatial regularization [6–8]. Numerical simulations on synthetic data demonstrated that taking into account the spatial structure significantly improves the estimate quality, provided that the connectivity structure between territories is known [8]. In practice, this structure is inaccessible and, while using the unweighted vicinity graph, in which territories sharing terrestrial borders are connected, showed promising results on French

départements [6, 7], it is too simplistic to analyze the pandemic at larger scales. To that aim, data-driven procedures inferring the connectivity structure are crucially needed.

**Related works.** Graph learning is a long-standing problem taking roots not only in graph processing but also in statistics, notably leading to the very popular graphical lasso estimator [9, 10], in machine learning, for semi-supervised learning [11], and network analysis [12]. However, it only recently emerged as a topic of interest in graph signal processing [13, 14]. In this framework, a popular strategy is to learn a graph in the topology of which the signal behaves smoothly, where the smoothness can be for example enforced via Laplacian-based total variation [15] which has some connections with graphical lasso [13]. During the past decade, joint estimation of the underlying graph and of the signal, in graph signal processing, or labels, in machine learning, received much attention [11, 13, 16]. A major advantage of joint strategies is that they result in a significantly smaller estimation error than two-step procedures, estimating first the graph from noisy data and then the regularized signal.

**Goals, contributions and outline.** Given infection count time series in several territories, the overarching goal of the present work is to jointly learn the spatial connectivity structure and estimate the reproduction number along time and across territories in an unsupervised manner. To that aim, elaborating on the graph representation framework developed in [16], the spatial structure is encoded in a weighted graph, represented through its combinatorial Laplacian matrix. A major originality of the present work is to perform Laplacian estimation under Poisson data-dependent noise, contrary to most of previous Laplacian estimation procedures developed under Gaussianity assumption [16, 17]. The second major originality is that, contrary to standard graph signal processing frameworks in which the signal is real-valued (e.g., year of graduation in a social network, transportation capacities at hub in a transportation network), each territory has its own infection count time series. Section 2 recasts the univariate temporal epidemiological model [3], into a multivariate spatiotemporal setting following [6, 7], leverages the scaled Poisson model [18] to enable variance modulation across territories, and adapts the graph-regularized reproduction number estimator from [6] to this novel setting. The proposed joint estimation of the spatiotemporal reproduction number and connectivity structure is developed in Section 3. Section 4 demonstrates on realistic synthetic data both the significant estimation accuracy gain obtained leveraging the spatial structure and the ability of the proposed method to retrieve very accurately the underlying connectivity structure and illustrates the method on real COVID-19 infection counts. Concluding remarks and perspectives are provided in Section 5.

**Notations.** Real (resp. non-negative real) numbers are denoted by $\mathbb{R}$ (resp. $\mathbb{R}_+$) and positive integers by $\mathbb{N}^*$. Scalars (resp. time series) are written in plain (resp. bold) characters and linear operators in upper case sans serif font. The component-wise product (resp. division) is denoted by $\odot$ (resp. $\cdot$/). Ground truth quantities are stared.
<figure><img src="image_2.png" /><figcaption>Fig. 1: COVID-19 infection counts reported in France, Italy, and the United-Kingdom from Sept. 1, 2020 to Oct. 1, 2021.<fnref n="1" /></figcaption></figure>

## 2. SPATIOTEMPORAL EPIDEMIOLOGICAL MODEL

**Scaled Poisson spatiotemporal epidemiological model.** The very small variance of the Poisson epidemiological model proposed in [3] poorly accounts for the large variability of COVID-19 data illustrated in Fig. 1 displaying infection counts in France, Italy and United-Kingdom collected from National Health Authorities and made available by the Johns Hopkins University.<fnref n="1" /> To better model COVID-19 propagation, a scaled Poisson model, with an extra scale parameter controlling the infection counts variance, has been proposed in [18]. The first contribution of this work is to extend this univariate description into a multivariate model of infection counts in different territories, hereafter referred to as countries. In country c at time t the infection count $Z_{c,t}$, conditionally to past counts $Z_{c,1}, \dots, Z_{c,t-1}$, is modeled as a scaled Poisson random variable

$$ \frac{Z_{c,t}}{\gamma_c} |Z_{c,1}, \dots, Z_{c,t-1}; R_{c,t} \sim \mathcal{P} \left( \frac{R_{c,t} \Phi_{c,t}^Z}{\gamma_c} \right). \quad (1) $$

The spatiotemporal Poisson parameter combines the effective reproduction number $R_{c,t} \ge 0$, the global infectiousness $\Phi_{c,t}^Z = \sum_{s=1}^{T_{\Phi}} \phi(s) Z_{c,t-s}$, which is a weighted sum of past counts where the serial interval function $\phi$ accounts for the random delay between primary and secondary infections, and a scale parameter $\gamma_c > 0$, constant over time, modulating count variance in country c. For COVID-19 pandemic, $\phi$ is modeled by a gamma distribution of mean 6.6 days and standard deviation 3.5 days truncated at $\tau_{\Phi} = 25$ days [19, 20]. The infection count time series over T days in the C different countries are represented as column vectors $Z_{c,:} \in \mathbb{R}^T$ for $c \in \{1, \dots, C\}$; the concatenation of the C time series is encoded in matrix form $Z \in \mathbb{R}_{+}^{C \times T}$. The present model assumes that the C random times series $Z_{c,:}$, $c \in \{1, \dots, C\}$ are independent and thus, that pandemic dynamic similarities between countries are entirely encoded in the resemblance of their reproduction number profiles.

**Regularized spatiotemporal variational estimator.** Due to the low quality of reported counts illustrated in Fig. 1, the reproduction number estimation requires advanced regularizing strategies. A first building block of the present work consists in designing a variational estimator of spatiotemporal reproduction numbers adapted to the multivariate scaled Poisson model of Eq. (1) enforcing a realistic, hence smooth, behavior both in time and across connected countries given a fixed connectivity structure. Following [6–8, 18], the data-fidelity term of [6, Eq. (7)] is replaced by a multivariate scaled Kullback-Leibler divergence, designed from the negative log-likelihood of the scaled Poisson model derived in [18, Eq. (20)]. Temporal regularity is obtained by enforcing a piecewise linear temporal behavior of the reproduction number, independently on each country, via the penalization of the $\ell_1$-norm of its time second derivative. To favor a smooth spatial behavior, the connectivity structure is encoded in a graph whose vertices are the countries and with edges linking countries expected to share similar epidemiological dynamics and a graph total variation penalization [15] is used to regularize the spatial variations of the reproduction number. While [6] used an unweighted graph, the present work proposes to consider a weighted

graph, enabling to account for different similarity levels and hence for more intricate, and thus more realistic, connectivity structures. Finally, to ease the design of a Laplacian learning strategy in Section 3, the $\ell_1$ spatial regularization is turned into a squared $\ell_2$ term. Altogether, this leads to the customized variational estimator

$$ \widehat{\mathbf{R}}^{\text{fix-L}} = \underset{\mathbf{R} \in \mathbb{R}_{+}^{C \times T}}{\operatorname{argmin}} \sum_{c=1}^{C} (\omega_c \mathrm{D}_{\mathrm{KL}}(\mathbf{Z}_{c,:} | \mathbf{R}_{c,:} \odot \mathbf{\Phi}_{c,:}^Z) + \lambda_T \|\mathrm{D}_2 \mathbf{R}_{c,:}\|_1) \\ + \lambda_S \sum_{t=1}^{T} \|\mathbf{B}\mathbf{R}_{c,t}\|_2^2 \quad (2) $$

where $\mathrm{D}_{\mathrm{KL}}(\mathbf{Z}_{c,:} | \mathbf{R}_{c,:} \odot \mathbf{\Phi}_{c,:}^Z) = \sum_{t=1}^{T} \mathrm{D}_{\mathrm{KL}}(\mathbf{Z}_{c,t} | \mathbf{R}_{c,t} \mathbf{\Phi}_{c,t}^Z)$ with $\mathrm{D}_{\mathrm{KL}}$ the Kullback-Leibler divergence;<fnref n="2" /> the weights $\omega_c > 0$ in the data-fidelity term would ideally be $\omega_c = 1/\gamma_c$ but could be replaced by heuristics when the scale parameters are unknown; the matrix $\mathrm{D}_2 \in \mathbb{R}^{T-2 \times T}$ denotes the second-order difference operator; B is the incidence matrix<fnref n="3" /> of the fixed graph encoding connectivity between countries; and $\lambda_T > 0$ (resp. $\lambda_S > 0$) are hyperparameters controlling the level of temporal (resp. spatial) regularization. The usual definition of the graph total variation involves the Laplacian matrix of the graph, denoted by $L = B^T B$. Remarking that $\|\mathbf{B}\mathbf{R}_{c,t}\|_2^2 = \mathbf{R}_{c,t}^T \mathbf{R}_{c,t}$ shows that the graph-based penalization in (2) is a rewriting of the graph total variation. Furthermore, replacing the incidence matrix $B \in \mathbb{R}^{E \times C}$ by any matrix $\tilde{B} \in \mathbb{R}^{C \times C}$ satisfying $\tilde{B}^T \tilde{B} = L$ lets the graph total variation unchanged, while significantly improving computational and memory efficiency when the number of edges $E$ is far larger than the number of countries $C$. The objective function in (2) is lower semi-continuous, proper, lower bounded by zero and coercive. Hence there exists a solution to (2). As the sum of convex functions composed with linear operators, the objective function is convex. Furthermore, the Kullback-Leibler divergence is strictly convex with respect to its second argument, then provided that $\mathbf{\Phi}_{c,t}^Z > 0$ for all countries c and times t, the solution of (2) is unique, and hence the fix-L estimator is well-defined.

**Algorithmic resolution.** Due to the presence of the $\ell_1$-norm, the objective function is nonsmooth, discarding the use of gradient-based minimization algorithms. Instead one has to resort to proximal algorithms [21]. The objective function in (2) can be rewritten as the sum of functions with closed-formed proximity operators [22] possibly composed with linear operators: $F(\mathbf{R}) + G(\mathbf{K}\mathbf{R})$ with $F(\mathbf{R}) = \sum_{c=1}^{C} \omega_c \mathrm{D}_{\mathrm{KL}}(\mathbf{Z}_{c,:} | \mathbf{R}_{c,:} \odot \mathbf{\Phi}_{c,:}^Z)$, K combining the linear operators involved in the regularization terms of (2) such that, $\mathbf{K}\mathbf{R} = (\mathrm{D}_2 \mathbf{R}_{1,:}, \dots, \mathrm{D}_2 \mathbf{R}_{C,:}, \mathbf{B}\mathbf{R}_{1,1}, \dots, \mathbf{B}\mathbf{R}_{T,1}) = (\mathbf{Q}_T, \mathbf{Q}_S) = \mathbf{Q}$ (3) where B is potentially equivalently replaced by $\tilde{B}$, and G defined on the range of K as $G(\mathbf{Q}_T, \mathbf{Q}_S) = \lambda_T \|\mathbf{Q}_T\|_1 + \lambda_S \|\mathbf{Q}_S\|_2^2$, where for $p \in \{1, 2\}$, $\|\cdot\|_p$ applies in a componentwise manner. The functions F, G being proper, lower-semicontinuous, convex, following [7, Algorithm 1], the minimization problem (2) can be solved taking advantage of the Chambolle-Pock primal-dual scheme [23]. Furthermore, the robust criterion proposed in [24, Eq. (8)] is leveraged to stop the iterations if a predefined precision $\varepsilon > 0$ is reached or if the maximum budget of iterations $k_{\max}$ has been exhausted.

<footnote n="2">²The Kullback-Leibler divergence is defined as $\mathrm{D}_{\mathrm{KL}}(Z|P) = Z \ln(Z/P) + P - Z$, if $Z > 0$, $P > 0$; P, if $Z = 0$, $P \ge 0$; $\infty$, otherwise.</footnote>

<footnote n="3">³Given an undirected graph with C vertices and E edges, the incidence matrix $B \in \mathbb{R}^{E \times C}$ is constructed so that if $e = \{c, c'\}$ is an arbitrarily oriented edge between vertices c and c' with $c < c'$, then $B_{e,c} = -B_{e,c'} = \sqrt{W_{c,c'}}$ and $B_{e,c''} = 0, \forall c'' \notin \{c, c'\}$, where $W_{c,c'}$ is the weight associated to edge $\{c, c'\}$ with $W_{c,c'} \in \{0, 1\}$ for unweighted graphs.</footnote>

<footnote n="1">¹coronavirus.jhu.edu</footnote>
**Algorithm 1 Alternating optimization**

**Input:** Z ∈ ℝ<sup>T</sup>, (λ<sub>T</sub>, λ<sub>S</sub>, λ<sub>L</sub>) ∈ ℝ<sup>3</sup>, n<sub>max</sub> ∈ N<sup>*</sup>, R<sup>(0)</sup>;

1: L<sup>(0)</sup>, M<sup>(0)</sup> ← minimizeL(R<sup>(0)</sup>, λ<sub>S</sub>, λ<sub>L</sub>)

2: ˜B<sup>(0)</sup> ← Cholesky(L<sup>(0)</sup>); define K from (3) using ˜B<sup>(0)</sup>

3: Q<sup>(0)</sup> ← KR<sup>(0)</sup>

**Alternating optimization:**

4: for n = 0 to n<sub>max</sub> - 1 do

5: R<sup>(n+1)</sup>, Q<sup>(n+1)</sup> ← minimizeR(Z, λ<sub>T</sub>, λ<sub>S</sub>, ˜B<sup>(n)</sup>, R<sup>(n)</sup>, Q<sup>(n)</sup>)

6: L<sup>(n+1)</sup>, M<sup>(n+1)</sup> ← minimizeL(R<sup>(n+1)</sup>, λ<sub>S</sub>, λ<sub>L</sub>, L<sup>(n)</sup>, M<sup>(n)</sup>)

7: ˜B<sup>(n+1)</sup> ← Cholesky(L<sup>(n+1)</sup>)

8: end for

9: Output: number R<sup>(n<sub>max</sub>)</sup>, L<sup>(n<sub>max</sub>)</sup>

### 3. JOINT ESTIMATION STRATEGY

Given a multivariate count time series Z ∈ ℝ<sup>C×T</sup>, the overarching goal of the present work is to estimate jointly the reproduction number R and the underlying connectivity structure. Among the numerous methods proposed in both graph signal processing and machine learning to address the simultaneous estimation of graph signals and connectivity structure, three are adapted to the variational framework: bi-level optimization [11], and joint [25] and alternate [16] minimization. Both bi-level and joint optimization require strong regularity and convexity properties on the respective objective functions involved to get convergence guarantees, thus, following [16], alternating minimization has been preferred. Two major contributions of the present work consist in designing a joint objective function and deriving a rapid alternating minimization scheme.

**Regularized and sparsity-inducing objective function.** Sticking to the variational framework of Section 2, the joint estimation of R and L is formulated as a constrained optimization problem:

$$ \hat{\mathbf{R}}^{\text{Joint}}, \hat{\mathbf{L}}^{\text{Joint}} \in \operatorname{Argmin}_{\mathbf{R} \in \mathcal{R}, \mathbf{L} \in \mathcal{L}} J(\mathbf{R}, \mathbf{L}). \quad (4) $$

The smart design of the objective function J and admissible sets R, L enables to enforce simultaneously spatial and temporal regularity of R, so that estimates are epidemiologically realistic, and sparsity of L, so that the graph only contains the most salient connections. Combining the regularized estimate defined in (2) and the graph learning strategy adapted to scalar Gaussian graph signals from [16] yields

$$ J(\mathbf{R}, \mathbf{L}) = \sum_{c=1}^{C} \omega_c D_{KL} \left( \mathbf{Z}_{c,:} \middle| \mathbf{R}_{c,:} \odot \mathbf{\Phi}_{c,:}^Z \right) \\ + \lambda_T \sum_{c=1}^{C} \| D_2 \mathbf{R}_{c,:} \|_1 + \lambda_S \sum_{t=1}^{T} \mathbf{R}_{:,t}^\top \mathbf{L} \mathbf{R}_{:,t} + \lambda_L \| \mathbf{L} \|_{\text{Fro}}^2, \quad (5) $$

where the graph total variation has been written using the Laplacian matrix instead of the incidence matrix as in (2), ‖·‖<sub>Fro</sub> denotes the Frobenius norm, ‖·‖<sub>2</sub><sup>2</sup><sub>Fro</sub> = Σ<sub>c,c'</sub> L<sup>2</sup><sub>c,c'</sub>, and λ<sub>L</sub> > 0 controls the estimated ̂L<sup>Joint</sup> sparsity level. The admissible sets are chosen as

$$ \begin{align} \mathcal{R} &= \{ \mathbf{R} \in \mathbb{R}^{C \times T}, \text{s.t. } \forall c, t \ \mathbf{R}_{c,t} \ge 0 \}, \nonumber \\ \mathcal{L} &= \left\{ \mathbf{L} \in \mathbb{R}^{C \times C}, \text{s.t. } \operatorname{Trace}(\mathbf{L}) = C, \right. \tag{6} \\ & \phantom{=} \left. \forall c \neq c', \mathbf{L}_{c,c'} = \mathbf{L}_{c',c} \le 0 \text{ and } \mathbf{L} \mathbf{1}_C = \mathbf{0}_C \right\}, \nonumber \end{align} $$

where 0<sub>C</sub> (resp. 1<sub>C</sub>) is the vector of size C with constant value 0 (resp. 1). To be consistent with (1), R enforces nonnegativity of the estimated reproduction numbers, while L ensures the estimation of a proper Laplacian matrix with a fixed total weight, excluding the trivial empty graph solution. The smaller the regularization parameter λ<sub>L</sub>, the less the total weight is able to spread, hence the sparser the estimated graph. The objective function J is continuous, lower

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Spatially structured synthetic data. Top left: Graph with C = 9 vertices split into I = 3 clusters; top right: I = 3 realistic reproduction numbers constructed following [8,24]; bottom: synthetic multivariate scaled Poisson counts in each of the I = 3 clusters.</figcaption></figure>

bounded by zero, coercive in R and defined on the compact admissible set L for L, thus the minimization problem (4) has at least one solution. Yet, the objective function, while convex in R and L separately, is not jointly convex. Minimizing (5) thus requires to resort to an alternating scheme, successively minimizing J with respect to R (resp. L) at fixed L (resp. R). The designed scheme is detailed in Algorithm 1. The inner optimization procedures are both solved efficiently with primal-dual schemes and systematically initialized using the previous primal and dual solutions to speed up convergence. **Minimization with respect to R at fixed L.** By equivalence of the incidence and Laplacian matrices-based formulations of graph total variation, the terms in J which depends on the reproduction number exactly match the objective function in (2). Thus, the inner step of Algorithm 1 minimizing over R can be performed leveraging the Chambolle-Pock scheme customized to (2) with connectivity structure encoded in a matrix ˜B link to the fixed Laplacian through L = ˜B<sup>⊤</sup> ˜B. For the sake of computational efficiency, ˜B is chosen of size C×C and obtained from Cholesky decomposition of L.<fnref n="4" /> At first call in Algorithm 1, the primal and dual variables of the Chambolle-Pock algorithm are initialized at R<sup>[0]</sup> = Z. / Φ<sup>Z</sup> and Q<sup>[0]</sup> = KR<sup>[0]</sup>. **Minimizing with respect to L at fixed R.** The terms of J depending on L, together with the constraints encoded in L, rewrites

$$ \hat{\mathbf{L}}^{\text{fix-R}} = \operatorname{argmin}_{\mathbf{L}} \lambda_{\mathbf{L}} \sum_{c,c'=1}^{C} \mathbf{L}_{c,c'}^2 + \lambda_{S} \sum_{c,c'=1}^{C} (\mathbf{R}_{c:,c}^{\top} \mathbf{R}_{c',:}) \mathbf{L}_{c,c'}, \quad (7) $$

s.t. Trace(L) = C, L<sub>c,c'</sub> = L<sub>c',c</sub> ≤ 0, ∀c ≠ c', L1<sub>C</sub> = 0<sub>C</sub>. Since λ<sub>L</sub> > 0, the objective function in (7), defined on a convex domain, is non-negative, continuous, coercive, and strictly convex, therefore it admits a unique minimizer. Furthermore, remarking that the constrained minimization problem (7) is a quadratic program subjected to linear constraints enables to solve it very efficiently using an interior point method [26,27]. In practice, the resolution of (7) is performed by the quadratic program solver<fnref n="5" /> of the CVXOPT package [28]. At first call in Algorithm 1, the primal and dual variables of the quadratic program are set to the solver default values.

### 4. NUMERICAL EXPERIMENTS

The proposed Joint estimator (4) is validated through intensive numerical simulations on carefully designed realistic synthetic data and then illustrated on real multivariate COVID-19 counts.

<footnote n="4">⁴ numpy.org/doc/stable/reference/generated/numpy.linalg.cholesky</footnote>

<footnote n="5">⁵ cvxopt.org/userguide/coneprog.html#quadratic-programming</footnote>
<table><caption>**Table 1: Reproduction number estimation accuracy.** Performance on synthetic data following Fig. 2 are quantified by $10^{-4} \times \text{mRSE}$ (8), averaged over 20 generations and accompanied by 95% Gaussian confidence intervals.</caption><thead><tr><th>ML</th><th>EpiEstim</th><th>fix-L<sub>0</sub></th><th>fix-L<sub>b</sub></th><th>fix-L<sup>*</sup></th><th>Joint (ours)</th></tr></thead><tbody><tr><td>113.2(±4.5)</td><td>27.7(±1.2)</td><td>9.0(±1.1)</td><td>5.4(±0.6)</td><td><b>2.6</b>(±0.4)</td><td><b>2.7</b>(±0.4)</td></tr></tbody></table>

**Synthetic multivariate scaled Poisson counts.** Another major contribution of the present work consists in the design of realistic synthetic scaled Poisson spatiotemporal counts following (1) accompanied by ground truth reproduction numbers and a prescribed connectivity structure. Indeed, the graph diffusion-based procedure to generate spatiotemporal synthetic data developed in [8] induces non negligible correlations between reproduction numbers in path-connected territories not directly connected by an edge, departing from the desired connectivity structure, and restricts to standard Poisson data. To strictly respect the prescribed spatial structure, only graphs composed of fully connected subgraphs disconnected from each other are considered. In practice, the C abstract countries are partitioned into I clusters $C_1, \dots, C_I$, constituting the prescribed spatial structure. In the following experiments, $C = 9$ abstract countries are considered and split into I = 3 clusters, as illustrated in Fig. 2 (top left). Then, for each cluster $C_i$, one realistic COVID-19 reproduction number time series $\mathbf{R}_{i, \cdot}^\dagger \in \mathbb{R}_+^I$ is generated following [8, 24] and used as ground truth for all countries of the cluster: $\forall c \in C_i, \mathbf{R}_{c, \cdot}^\dagger = \mathbf{R}_{i, \cdot}^\dagger$, hence exactly enforcing the desired spatial structure. Finally, initial counts $\mathbf{Z}_{i, 0}$ are fixed and scale parameters are chosen large enough to allow variability between counts inside each cluster via $\gamma_c = 0.01 \times Z_c^0$. Finally multivariate scaled Poisson synthetic counts $\mathbf{Z}_{c, \cdot}$ are independently sampled, following Eq. (1). Fig. 2 displays the realistic reproduction numbers (top right) and a generation of synthetic counts (bottom) for each of the three clusters.

**Experimental protocol.** Performance of the reproduction number estimators are measured through the mean Relative Squared Error:

$$ \text{mRSE}(\mathbf{R}) = \frac{1}{CT} \sum_{c=1}^{C} \sum_{t=1}^{T} \left( \frac{\hat{\mathbf{R}}_{c,t} - \mathbf{R}_{c,t}^{*}}{\mathbf{R}_{c,t}^{*}} \right)^2 . \quad (8) $$

The proposed Joint strategy is compared with five methods. The ML estimator consists in maximizing the likelihood associated to (1), leading to $\hat{\mathbf{R}}_{c,t}^{\text{ML}} = Z_{c,t}/\Phi_{c,t}^Z$. The EpiEstim method [3, 5] processes independently infection counts of all countries in a Bayesian framework with a temporal smoothing on a window of length $\tau$. To test the robustness of the fixed graph method fix-L (2) to errors in the prescribed spatial structure, three spatial structures encoded in the Laplacian matrix are used: an empty graph $L_0$, a blurred version $L_b$ of the true graph, with spurious edges of small weights, and the true graph $L^*$. Regularization parameters are optimized through grid searches and the best mRSE is reported. Performance are averaged over 20 independent synthetic counts generations and uncertainty is assessed through 95% Gaussian confidence intervals.

**Algorithmic settings.** To mimic real world data for which the scale parameters $\gamma_1, \dots, \gamma_C$ are unknown, the Kullback-Leibler weights in (2) are chosen following the heuristic $\forall c \in \{1, \dots, C\}, \omega_c = 1/\text{std}(\mathbf{Z}_{c, \cdot})$ inspired from [6, 7]. In the iterative procedures computing the fix-L and Joint estimates, the reproduction numbers are initialized at the countrywise EpiEstim estimates with a sliding window of $\tau = 7$ days. The number of iterations in Algorithm 1 is set to $n_{\text{max}} = 10$. Regularization parameters are optimized through exhaustive search across: all odd numbers between 1 and 29 for the sliding window length $\tau$ in EpiEstim method; logarithmic grids of $\lambda_T, \lambda_S$ and $\lambda_L$ ranging respectively from 1 to 100, 0.01 to 1000 and

<figure><img src="image_4.png" /><figcaption>**Fig. 3: Left:** Map of the countries from the real data experiments. Colors indicate clusters derived from the spatial structure inferred by Joint. Hashed countries are isolated ones. Right: Reproduction numbers estimated from the Joint method organized by clusters.</figcaption></figure>

0.001 to 100 for the penalized estimates, exploring 64 values of $\lambda_T$ for fix-L$_0$, a grid of $16 \times 16$ couples ($\lambda_T, \lambda_S$) for fix-L and a grid of $8 \times 8 \times 8$ triplets ($\lambda_T, \lambda_S, \lambda_L$) for the Joint estimator.

**Discussion.** Table 1 reports estimation accuracies achieved by the six considered methods. Both the naive ML (first column) and state-of-the-art EpiEstim [3, 5] (second column) lie far behind the penalized fix-L and Joint strategies (third to sixth columns) in terms of mRSE. The fix-L method is robust to mild mismatches between the prescribed spatial structure used in fix-L and the true connectivity graph (fourth and fifth columns). Though, the fix-L$_0$ estimator enforcing only temporal regularity (third column) is superseded by the fix-L$_b$, fix-L$^*$ and Joint (fourth to sixth column) strategies including spatial regularization. Finally, the Joint method (sixth column) reaches the same excellent accuracy as the fix-L one (fifth column) using the true graph. Indeed, the estimated connectivity structure is almost exact, achieving a relative squared error of $\|\hat{\mathbf{L}}^{\text{Joint}} - \mathbf{L}^*\|_{\text{Fro}}^2/\|\mathbf{L}^*\|_{\text{Fro}}^2 = 2.3(\pm 1.0) \times 10^{-13}$ and a perfect estimation of the support of $\mathbf{L}^*$ when thresholding weights at $10^{-14}$.

**COVID-19 data.** The Joint method is illustrated on a selection of real COVID-19 infection counts in 39 countries in Europe and Africa between Sept. 1, 2020 and Oct. 1, 2021 downloaded from Johns Hopkins University's repository.¹ The optimal hyperparameters estimated through grid search on synthetic data are manually fine-tuned to adjust the regularization level.<fnref n="6" /> The inferred graph shows four exactly disconnected clusters (in red, orange, green and purple) and a few isolated countries (dashed regions), displayed in Fig. 3 (left). The estimated reproduction numbers in each country are grouped by cluster in Fig. 3 (right).<fnref n="6" /> The Joint method captures the main representative epidemiological dynamics across a large geographical area, which would have not been possible through visual inspection of the low quality counts, examples of which are displayed in Fig. 1. Interestingly, focusing on epidemiological dynamics, the Joint method does not only group together geographically close countries, e.g., United-Kingdom, Tunisia and Czech Republic in the red cluster.

## 5. CONCLUSIONS AND PERSPECTIVES

A joint estimator of reproduction numbers and spatial connectivity structure has been proposed. Intensive numerical experiments demonstrated its remarkable ability to infer spatiotemporal epidemiological indicators in an accurate and robust manner, while providing precise insights on the underlying spatial pattern. The robustness of this Joint method to change in the pool of monitored territories and in the phase of the pandemic will be systematically explored on real data. For reproducibility, the code is publicly available on Github.<fnref n="7" />

<footnote n="6">⁶ A notebook provides precise experimental setup and extra results.</footnote>

<footnote n="7">⁷github.com/elasalle/EpiJointSpatiотempEstim</footnote>
## 6. REFERENCES

[1] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz, "On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations," *J. Math. Biol.*, vol. 28, pp. 365–382, 1990.

[2] J. Wallinga and P. Teunis, "Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures," *Am. J. of Epidemiol.*, vol. 160, no. 6, pp. 509–516, 09 2004.

[3] A. Cori, N. M. Ferguson, C. Fraser, and S. Cauchemez, "A new framework and software to estimate time-varying reproduction numbers during epidemics," *Am. J. Epidemiol.*, vol. 178, no. 9, pp. 1505–1512, 2013.

[4] Q.-H. Liu, M. Ajelli, A. Aleta, S. Merler, Y. Moreno, and A. Vespignani, "Measurability of the epidemic reproduction number in data-driven contact networks," *Proc. Natl. Acad. Sci. U.S.A.*, vol. 115, pp. 12680–12685, 2018.

[5] Robin N Thompson, Jake E Stockwin, Rolina D van Gaalen, Jonny A Polonsky, Zhian N Kamvar, P Alex Demarsh, Elisabeth Dahlqwist, Siyang Li, Eve Miguel, Thibaut Jombart, et al., "Improved inference of time-varying reproduction numbers during infectious disease outbreaks," *Epidemics*, vol. 29, pp. 100356, 2019.

[6] P. Abry, N. Pustelnik, S. Roux, P. Jensen, P. Flandrin, R. Gribonval, C.-G. Lucas, É. Guichard, P. Borgnat, and N. Garnier, "Spatial and temporal regularization to estimate COVID-19 reproduction number $R(t)$: Promoting piecewise smoothness via convex optimization," *PLOS One*, vol. 15, no. 8, pp. e0237901, 2020.

[7] B. Pascal, P. Abry, N. Pustelnik, S. Roux, R. Gribonval, and P. Flandrin, "Nonsmooth convex optimization to estimate the Covid-19 reproduction number space-time evolution with robustness against low quality data," *IEEE Trans. Signal Process.*, vol. 70, pp. 2859–2868, 2022.

[8] J. Du, B. Pascal, and P. Abry, "Synthetic Spatiotemporal Covid19 Infection Counts to Assess Graph-Regularized Estimation of Multivariate Reproduction Numbers," Lyon, France, August. 26-30 2024.

[9] A. P. Dempster, "Covariance selection," *Biometrics*, pp. 157–175, 1972.

[10] J. Friedman, T. Hastie, and R. Tibshirani, "Sparse inverse covariance estimation with the graphical lasso," *Biostatistics*, vol. 9, no. 3, pp. 432–441, 2008.

[11] H. Ghanem, S. Vaiter, and N. Keriven, "Gradient Scarcity in Graph Learning with Bilevel Optimization," *Trans. Mach. Learn. Res.*, 2024.

[12] E. D. Kolaczyk, *Statistical Analysis of Network Data: Methods and Models*, vol. New York, NY, USA, Springer-Verlag, 2009.

[13] X. Dong, D. Thanou, Michael Rabbat, and P. Frossard, "Learning graphs from data: A signal representation perspective," *IEEE Signal Proces. Mag.*, vol. 36, no. 3, pp. 44–63, 2019.

[14] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro, "Connecting the dots: Identifying network structure via graph signal processing," *IEEE Signal Process. Mag.*, vol. 36, no. 3, pp. 16–43, 2019.

[15] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," *IEEE Signal Process. Mag.*, vol. 30, no. 3, pp. 83–98, 2013.

[16] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, "Learning Laplacian Matrix in Smooth Graph Signal Representations," *IEEE Trans. Signal Process.*, vol. 64, no. 23, pp. 6160–6173, 2016.

[17] H. E. Egilmez, E. Pavez, and A. Ortega, "Graph Learning From Data Under Laplacian and Structural Constraints," *IEEE J.J. Sel. Top. Signal Process.*, vol. 11, no. 6, pp. 825–841, 2017.

[18] B. Pascal and S. Vaiter, "Risk Estimate under a Nonstationary Autoregressive Model for Data-Driven Reproduction Number Estimation," *Preprint arXiv:2409.14937*, 2024.

[19] D. Cereda, M. Tirani, F. Rovida, V. Demicheli, M. Ajelli, P. Poretti, F. Trentini, G. Guzzetta, V. Marziano, and A. Barone, "The early phase of the COVID-19 outbreak in Lombardy, Italy," *Preprint arXiv:2003.09320*, 2020.

[20] F. Riccardo, M. Ajelli, X. D. Andrianou, A. Bella, M. Del Manso, M. Fabiani, S. Bellino, S. Boros, A. M. Ur diales, and V. Marziano, "Epidemiological characteristics of COVID-19 cases and estimates of the reproductive numbers 1 month into the epidemic, Italy, 28 January to 31 March 2020," *Eurosurveillance*, vol. 25, no. 49, pp. 2000790, 2020.

[21] N. Parikh and S. Boyd, "Proximal algorithms," *Foundations and Trends® in Optimization*, vol. 1, no. 3, pp. 127–239, 2014.

[22] G. Chierchia, E. Chouzenoux, P.L. Combettes, and J.-C. Pesquet, "The proximity operator repository," http://proximity-operator.net/.

[23] A. Chambolle and T. Pock, "A first-order primal-dual algorithm for convex problems with applications to imaging," *J. Math. Imag. Vis.*, vol. 40, no. 1, pp. 120–145, 2011.

[24] J. Du, B. Pascal, and P. Abry, "Compared performance of Covid19 reproduction number estimators based on realistic synthetic data," in *GRETSI'23 XXIXème Colloque Francophone de Traitement du Signal et des Images*, Grenoble, France, Aug. 28 - Sept. 1 2023.

[25] B. Fatemi, L. El Asri, and S. M. Kazemi, "SLAPS: Self-supervision improves structure learning for graph neural networks," *Adv. Neural Inf. Process. Syst.*, vol. 34, pp. 22667–22681, 2021.

[26] M. C. Grant and S. P. Boyd, "Graph implementations for nonsmooth convex programs," in *Recent advances in learning and control*. Springer, 2008, pp. 95–110.

[27] M. Andersen, J. Dahl, Z. Liu, L. Vandenberghe, S. Sra, S. Nowozin, and S. Wright, "Interior-point methods for large-scale cone programming," in *Optimization for Machine Learning*, vol. 5583. MIT Press, 2011.

[28] M. S. Andersen, J. Dahl, and L. Vandenberghe, "CVXOPT: A Python package for convex optimization," *Available at cvxopt.org*, vol. 54, 2013.