Member, IEEE Barbara Pascal
email: barbara.pascal@univ-lille.fr

Fellow, IEEE Patrice Abry

Member, IEEE Nelly Pustelnik

Stéphane Roux

Fellow, IEEE Rémi Gribonval

Patrick Flandrin

Life Fellow

Nonsmooth convex optimization to estimate the Covid-19 reproduction number space-time evolution with robustness against low quality data

Keywords: Covid-19, reproduction number, space-time evolution, nonsmooth convex optimization, outlier robustness

Daily pandemic surveillance, often achieved through the estimation of the reproduction number, constitutes a critical challenge for national health authorities to design countermeasures. In an earlier work, we proposed to formulate the estimation of the reproduction number as an optimization problem, combining data-model fidelity and space-time regularity constraints, solved by nonsmooth convex proximal minimizations. Though promising, that first formulation significantly lacks robustness against the Covid-19 data low quality (irrelevant or missing counts, pseudo-seasonalities,. . . ) stemming from the emergency and crisis context, which significantly impairs accurate pandemic evolution assessments. The present work aims to overcome these limitations by carefully crafting a functional permitting to estimate jointly, in a single step, the reproduction number and outliers defined to model low quality data. This functional also enforces epidemiology-driven regularity properties for the reproduction number estimates, while preserving convexity, thus permitting the design of efficient minimization algorithms, based on proximity operators that are derived analytically. The explicit convergence of the proposed algorithm is proven theoretically. Its relevance is quantified on real Covid-19 data, consisting of daily new infection counts for 200+ countries and for the 96 metropolitan France counties, publicly available at Johns Hopkins University and Santé-Publique-France. The procedure permits automated daily updates of these estimates, reported via animated and interactive maps. Open-source estimation procedures will be made publicly available.

I. INTRODUCTION

Context. The ongoing COVID-19 pandemic has produced an unprecedented health and economic crisis, urging for the construction of efficient monitoring procedures of its spreading across territories, a crucial step in designing sanitary, social and economic policies by national authorities [START_REF] Flahault | Covid-19 cacophony: is there any orchestra conductor?[END_REF]. It is however often not the value of the infection level per se that matters to design pandemic counter-measures (lockdown,. . . ), but its evolution along time and variations across territories. However, in the context of the Covid-19 pandemic outburst and of the sanitary crisis, collecting daily new infection counts, the basis material in any pandemic surveillance strategy, results in low-quality data (missing samples, outliers, seasonalities. . . ). Surprisingly, after 18 months of pandemic, the quality of the data collected and made available by most national health authorities in the world remains limited, which severely impairs an accurate and timely pandemic evolution assessment, the issue at the heart of the present work. Related works. Pandemic surveillance has been envisaged with several categories of tools from different fields of sciences (cf. [START_REF] Arino | Describing, modelling and forecasting the spatial and temporal spread of COVID-19-A short review[END_REF] for a review). Yet, it is classically performed using compartmental models, elaborating on the classical Susceptible-Infectious-Recovered reference. For realistic use and to match social realities (social groups, contact, . . . ), these models need to be refined (cf. e.g., [START_REF] Liu | Measurability of the epidemic reproduction number in data-driven contact networks[END_REF][START_REF] Brauer | Mathematical models in epidemiology[END_REF]), essentially by increasing the number of compartments, which implies increasing (in a quadratic way) the numbers of parameters to be estimated. Parameter estimation is usually achieved within Bayesian frameworks, maximizing the likelihood attached to the models, at the expense of heavy computational burdens. Such models are thus often used a posteriori (i.e., after the epidemic) with consolidated and accurate datasets. The low quality of the Covid-19 data collected by most national public health authorities across the world significantly impairs the use of such Bayesian schemes and thus the reliable estimations of these parameters. This thus precludes the use of such models for daily basis update of the pandemic evolution assessment.

Instead of compartmental models, epidemiologists massively use the so-called reproduction number, R, that measures how many new individuals are on average infected by an already contaminated person, a key marker of the pandemic strength (cf. e.g., [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Wallinga | Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF][START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF][START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF]). Estimated to around 3 during the outburst of the Covid-19 pandemic [START_REF] Di Domenico | Expected impact of lockdown in Ile-de-France and possible exit strategies[END_REF][START_REF] Salje | Estimating the burden of SARS-CoV-2 in France[END_REF], it is also used as a function of time R T to monitor the temporal evolution of the pandemic. The strength of such a pandemic monitoring is to be based on a single parameter to estimate R, while accounting for the basic and universal pandemic propagation mechanisms: The number of new infections, Z T , at day T , depends on R T and on the numbers of new infections measured at previous days, {Z T -1 , Z T -2 , Z T -3 , . . .} weighted by the so-called serial interval function Φ(t). The latter quantifies the probability that someone found Covid-positive today was actually contaminated several days ago [START_REF] Liu | Measurability of the epidemic reproduction number in data-driven contact networks[END_REF][START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF][START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF][START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF]. Following [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF], R T can be estimated using a Bayesian scheme 1 .

1 https://shiny.dide.imperial.ac.uk/epiestim/ Instead, in an earlier work [START_REF] Abry | Spatial and temporal regularization to estimate COVID-19 reproduction number R (t): Promoting piecewise smoothness via convex optimization[END_REF], we proposed an estimation procedure of the spatio-temporal evolution of the reproduction number, written as a nonsmooth convex optimization problem, combining data-model fidelity with time and space regularity constraints for the estimates of R to be pandemic realistic, and solved by proximal algorithms [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Cai | Image restoration: Total variation, wavelet frames, and beyond[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Parikh | Proximal algorithms[END_REF][START_REF] Pustelnik | Wavelet-based image deconvolution and reconstruction[END_REF]. Though promising, that proximal-optimization based estimation significantly lacks robustness against the low quality of the Covid-19 data, a generic and crucial issue addressed in the present contribution. Goals, contributions and outline. The overall goal of this work is to devise an inverse problem-type convex-optimization procedure to assess the spatio-temporal evolution of the reproduction number R, from daily new infection counts, measured simultaneously on several connected territories, such as the different counties of a same country. The proposed estimation procedure is designed to be: i) robust to the poor quality of Covid-19 data ; ii) efficiently applicable online to update estimates, on a daily basis or whenever new data are available, at moderate computational costs even for large datasets.

To that end, Section II recalls the pandemic model used here [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF], frames its classical estimation and related issues. Section III describes the core methodological contributions of the paper: i) It discusses the poor quality of the data and proposes an extension of the original pandemic modelling that handles poor quality data as outliers (cf. Section III-A) ; ii) It proposes an original and theoretical design of a functional whose minimization yields in a single step a joint estimation of both the reproduction number and the outliers, by combining the log-likelihood of the extended model as a data fidelity term with temporal and spatial regularity constraints for epidemiology-realistic estimates of R and sparsity constraints for the outliers (cf. Section III-B). This functional is carefully crafted to preserve convexity so that fast and efficient proximal algorithms can be devised for minimization ; iii) The theoretical properties of the functional and of its solution (existence, unicity) are studied theoretically in Section III-C ; iv) This permits to devise a proximal type iterative algorithm, with explicit analytical derivation of the proximal operators, whose convergence and stopping criterion are discussed theoretically and practically (cf. Section III-D).

The relevance of the proposed estimation procedure to assess the spatio-temporal evolution of the pandemic is illustrated on real Covid-19 data (daily counts of new infections) made available by national public health authorities (such as e.g., Santé-Publique-France, SPF) or collected by Johns Hopkins University, described in Section IV.

Estimation performance are reported in Section V: i) First, the role and choices of the hyperparameters balancing the regularization terms and the data-model fidelity term are discussed in detail (cf. Section V-A) ; ii) Second, the benefits of the proposed strategy to estimate the temporal evolution of R T at the level of a given country are illustrated by comparisons against different estimation strategies, such as a two-step procedure (outliers estimated and removed first, followed by estimation of R) (cf. Section V-B) ; iii) Third, the relevance of the space-time evolution of R estimated across territories that are connected (shared borders, massive population commutations,. . . ), here the 96 départements (counties) that organize administratively metropolitan France), is illustrated and discussed precisely (cf. Section V-C).

The proposed procedure is efficient enough to permit the automated daily updates of the estimates of R T for 200+ countries, for the 96 continental French départements of France and for the 50 US states. Links to the resulting animated and interactive maps of estimates are provided and open-source estimation procedures will be made publicly available.

II. PANDEMIC MODEL AND ESTIMATION

A. Model

The pandemic model used here, focused on the reproduction number R, was proposed in [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF] and relies on two key ideas: i) Conditionally to past counts Z 1:T -1 {Z 1 , . . . , Z T -1 }, the count of daily new infections, Z T , is modeled as a random variable drawn from a Poisson distribution Poiss(p T ), ii) whose parameter p T depends on the past counts Z 1:T -1 , on the causal serial interval function Φ T and on the current R T :

P (Z T |Z 1:T -1 ; R T ) ∼ Poiss p T R T × τ φ t=1 Φ t Z T -t (1)
The serial interval function Φ t is a key element of the model that accounts for the epidemiology evolution mechanisms, as it models the random delays between the onset of symptoms in a primary case and the onset of symptoms in secondary cases, [START_REF] Liu | Measurability of the epidemic reproduction number in data-driven contact networks[END_REF][START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF][START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF][START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF]. For the Covid-19 pandemic and for earlier pandemics of same types, it was shown that Φ can be approximated as a Gamma function, with shape and rate parameters 1.87 and 0.28, respectively, corresponding to mean and standard deviations of 6.6 and 3.5 days, indicating a high risk of infecting other persons from 3 to 10 days after the symptoms have appeared [START_REF] Ma | Epidemiological parameters of coronavirus disease 2019: A pooled analysis of publicly reported individual data of 1155 cases from seven countries[END_REF][START_REF] Riccardo | Epidemiological characteristics of COVID-19 cases in Italy and estimates of the reproductive numbers one month into the epidemic[END_REF][START_REF] Guzzetta | The impact of a nation-wide lockdown on COVID-19 transmissibility in Italy[END_REF]. It is assumed here that Φ is known and follows this Gamma approximation.

B. Estimation

Given the analytic expression of the Poisson distribution in Model (1), the log-likelihood ln

P (Z T |Z 1:T -1 ; R 1:T ) = Z T ln p T -p T -ln Z T !, permits to compute the Maximum- Likelihood Estimate, RMLE T = argmin(-ln P (Z T |Z 1:T -1 ; R 1:T )), (2)
in closed-form expression as:

RMLE T = Z T / t≥1 Φ t Z T -t .
It can be interpreted as a ratio of moving averages, whose size and shape are consistent with pandemic modeling (thus with Φ T ). The RMLE T computed from real Covid-19 data (cf. Section IV) are displayed in Figs. 2 and3 (middle plots) and display a far too large variability along time to be of any practical use in actual pandemic monitoring. Such variability mostly stems from the low quality of the data. The goal of the present work is thus to improve the estimation of R despite the low quality of the available Covid-19 data.

III. ROBUST ESTIMATION: METHODOLOGY A. Covid-19 data low quality and multivariate infection counts

Covid-19 daily new infection counts made available by public health authorities are for most countries severely corrupted, with missing samples, non meaningful negative counts followed by retrospective cumulated counts spread across the following days,. . . . They also show pseudo-seasonality effects, with significantly less counts on non working days. Seasonality are however uneasy to model, as non working days are not only week-ends but also additional day-offs that depend on countries. Also, the way the pandemic has been monitored along week-ends and non-working days has significantly changed along time, even within a given country. Therefore, rather than trying to devise necessarily ad-hoc and parametric models for noise in Covid-19 data, we have instead chosen to refer to these different types of data corruption under the generic term of outlier, denoted O T , and we propose to model observed daily new infection counts Z T still from a Poisson distribution conditionally to past counts Z 1:T -1 , yet with a Poisson parameter p T that depends both on the current reproduction number R T and current outliers O T .

Additionally, when new infection counts are monitored simultaneously for several connected territories, e.g., the counties or states of a given country, the observed daily counts consist of D time series {Z

(d) t , t ∈ [1, T ], d ∈ [1, D]},
(d) t = τΦ u=1 Φ u Z (d) t-u : P Z (d) t |Z (d) 1:t-1 ; R (d) t , O (d) t (3) ∼ Poiss p (d) t R (d) t × (ΦZ) (d) t + O (d) t .

B. Crafting the convex nonsmooth estimation functional

The goal is now to estimate jointly R = {R

(d) t , t ∈ [1, T ], d ∈ [1, D]} and O = {O (d) t , t ∈ [1, T ], d ∈ [1, D]} from {Z (d) t , t ∈ [1, T ], d ∈ [1, D]}.
The next original contribution is thus to frame the estimation into an inverse problem strategy, with a careful crafting of the functional to minimize combining the negative log-likelihood of extended Model (3) with positiveness and regularity constraints in time and space on R and structure on O.

1) Negative log-likelihood function: The negative loglikelihood associated to the extended Poisson Model (3) is defined in a separable manner as

F (R, O|Z) D d=1 T t=1 f (R (d) t , O (d) t |Z (d) t , (ΦZ) (d) t ), (4)
with f (r, o|z, (Φz))

ι {(0,0)} (r, o) if z = Φz = 0, d KL (z|r × (Φz) + o) otherwise.
The indicator function of the singleton {(0, 0)} permits to force the (MLE) estimates of R and O to vanish when the pandemic stops, i.e., when (ΦZ)

(d) t = Z (d) t = 0, ι {(0,0)} (r, o) = 0 if r = o = 0 ∞ otherwise, (5)
while the standard Kullback-Leibler divergence (KLD) is defined as

D KL (Z|P ) = T t=1 D d=1 d KL (Z (d) t |P (d) t ), with d KL (z|p)      z ln z p + p -z if z > 0, p > 0 p if z = 0, p ≥ 0 ∞ else. (6)
2) Regularity and positivity constraints on R: Unlike Model (1), from which a unique MLE estimate can be derived (Eq. ( 2)), for Model (3) the extended likelihood (Eq. ( 4)) is no longer strictly convex, thus does not have a unique maximum, and hence requires additional constraints on the estimates of R and O. i) Time regularity: By nature and to be of any potential use for actual pandemic monitoring, the estimated reproduction numbers for a given territory must vary only slowly and smoothly along time. Following [START_REF] Abry | Spatial and temporal regularization to estimate COVID-19 reproduction number R (t): Promoting piecewise smoothness via convex optimization[END_REF], it has been chosen to enforce piecewise linear time evolution of the estimate of the reproduction numbers. Such a property ensures that, besides an estimate of R at day t, one also gets an estimation of a local trend, assessing whether the pandemic is growing or decreasing. Following [START_REF] Colas | Nonlinear denoising for solid friction dynamics characterization[END_REF][START_REF] Debarre | Sparsest Continuous Piecewise-Linear Representation of Data[END_REF], to favor piecewise linear temporal evolution of the estimates, we penalize the 1 -norm of the multivariate time-domain Laplacian D 2 : R D×T → R D×(T -2) of the estimates of R:

D 2 R 1 D d=1 T -1 t=2 |(D 2 R) (d) t |, (7)
with

(D 2 R) (d) t 1 2 R (d) t-1 -R (d) t + 1 2 R (d)
t+1 . ii) Positivity: R is by nature positive. Yet, Eq. ( 4) indicates that, per se, the proposed extended negative log-likelihood ensures the positivity of p

(d) t = R (d) t × (ΦZ) (d) t + O (d) t but not of R (d) t
itself. Thus, to enforce the positivity of R (d) t , an indicator function is added in the constraints:

ι ≥0 (R) ∞ if ∃ t ∈ [1, T ], d ∈ [1, D], s.t. R (d) t < 0 0 otherwise. (8)
iii) Spatial regularity: When different territories d are connected, it is natural to expect that their pandemic dynamics are correlated and thus that the corresponding estimates of R show some spatial regularity. Following [START_REF] Abry | Spatial and temporal regularization to estimate COVID-19 reproduction number R (t): Promoting piecewise smoothness via convex optimization[END_REF], we use a graph, with D vertices corresponding to the different territories and E edges to connections between pairs of territories. We favor piecewise constant spatial estimates of R1:D t across territories by penalizing the Graph Total Variation, defined as:

GTV(R) = T t=1 d1∼d2 R (d1) t -R (d2) t , (9)
where d1∼d2 runs over vertices connected by an edge. The Graph Total Variation can be rewritten as a composition of a discrete difference operator and a sparsity-promoting 1norm, GTV(R) = GR 1 , defining the linear operator G:

G : R D×T → R E×T R → R (d1) t -R (d2) t , t ∈ [1, T ], d 1 ∼ d 2 .
(10)

3) Constraints on the outliers: The discussion in Section III-A suggests that outliers have a sparse structure: They occur with significant values only on specific days (sundays, day-offs,. . . ) and are essentially negligible elsewhere. Sparsity is classically enforced in inverse problem settings by constraining the 1 -norm. We thus further propose to impose the following constraint to the estimates of the outliers:

O 1 D d=1 T t=1 |O (d) t |. (11)
4) Estimation as an optimization problem: Combining the data-fidelity term of Eq. ( 4) and the penalizations of Eqs. ( 7), ( 8), [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF], and [START_REF] Salje | Estimating the burden of SARS-CoV-2 in France[END_REF], we propose to obtain the estimates of R and O by minimizing:

( R, Ô) ∈ Argmin R,O∈R D×T F (R, O|Z) + λ T D 2 R 1 + ι ≥0 (R) + λ S GR 1 + λ O O 1 , (12)
with λ T > 0, λ S > 0 and λ O > 0, regularization hyperparameters, balancing the strengths of the different constraints one against the others and against the likelihood.

5) Reformulation of the optimization problem:

To study the theoretical properties of Problem [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF], and to derive a minimization algorithm, it is useful to recast Eq. ( 12) into a generic formulation:

( R, Ô) = Argmin R,O∈R D×T F (R, O|Z) + H(L(R, O)), (13)
with L a linear operator and H a function enforcing constraints promoting certain sparsity / non-negativity, defined as:

L : R D×T 2 → R D×(T -2) × R D×T × R E×T × R D×T (R, O) → (λ T D 2 R, R, λ S GR, λ O O) , (14)
H(Q 1 , Q 2 , Q 3 , Q 4 ) Q 1 1 + ι ≥0 (Q 2 ) + Q 3 1 + Q 4 1 . (15)

C. Convexity, existence of a minimizer and uniqueness

Showing that Eq. ( 12) provides well-defined estimates of R and O requires to prove the existence of a minimizer of F (•) + H(L•). To that aim, preliminary results are necessary.

Lemma 1. The extended Kullback-Leibler data-fidelity term F (Eq. ( 4)) is jointly convex with respect to the pair (R, O).

Proof. F being fully separable in time and space, it is sufficient to demonstrate the joint convexity of f (r, o|z, (Φz)) with respect to the pair (r, o). If z = (Φz) = 0, the convexity of f (r, o|z, (Φz)) w.r.t. (r, o) reduces to the convexity of the indicator function of the convex set {(0, 0)}. Else, since the one-component KLD d KL (z|p) is convex with respect to its second argument [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF], and (r, o) → r × (Φz) + o being linear, (r, o) → f (r, o|z, (Φz)) is convex.

Corollary 1. The penalized Kullback-Leibler functional (Eq. ( 12)) is jointly convex with respect to the pair (R, O).

Proof. First, Lemma 1 ensures the convexity of the datafidelity term F . Then, the function H, appearing in Eq. ( 13), involves 1 -norms and an indicator function which are convex, thus it is convex and the composition of H with the linear operator L is convex as well. As the sum of two convex terms, the penalized Kullback-Leibler functional is convex.

Further, because the Kullback-Leibler data-fidelity term F (Eq. ( 4)), as well as the indicator function in Eq. ( 8), may take infinite values, the domain of the convex functional in Eq. ( 12) is not the whole space R D×T 2 , but restricts to a domain Ω of feasible variables, described in Lemma 2.

Lemma 2. The domain Ω of feasible points for Problem [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF]

consists of pairs (R, O) such that, ∀u ∈ [1, t], ∀d ∈ [1, D] (i) R (d) t (ΦZ) (d) t + O (d) t ≥ 0, (ii) if Z (d) t = (ΦZ) (d) t = 0 then R (d) t = O (d) t = 0, (iii) R (d) t ≥ 0.
It is convex and closed.

Proof. Conditions (i) and (ii) reflect the fact that the extended Kullback-Leibler data-fidelity term F (Eq. ( 4)) must be finite for feasible time series, while Condition (iii) stems from the positivity constraint enforced by the indicator function in Eq. ( 8). Further, each Condition (i) to (iii) defines a closed convex set of feasible time series, the finite intersection of closed convex sets being closed and convex as well, Ω is a closed convex set.

Then, the existence of a minimizer can be proven.

Theorem 1. The penalized Kullback-Leibler functional in Eq. ( 12) is lower bounded. Further, if (ΦZ)

(d) t > 0 for all t ∈ [1, T ], d ∈ [1, D],
then the minimization problem of Eq. (12) has, a least, one solution.

Proof. Since, for feasible points, the standard Kullback-Leibler divergence, the indicator functions and the 1 -norms are nonnegative, the functional is lower bounded by zero.

Let (R, O) 2 R 2 2 + O 2
2 be the canonical norm on the product space R D×T 2 . Let 1 (resp. 0) the multivariate time series with all entries equal to 1 (resp. 0),

(1, 0) ∈ Ω is a feasible point. Let µ F (1, 0|Z) + H(L(1, 0)) = F (1, 0|Z) ≥ 0 (16)
As shown below, the functional

F (•|Z) + H(L•) is coercive, thus there exists α > 0 such that (R, O) 2 > α ⇒ F (R, O|Z) + H(L(R, O)) > µ. ( 17
)
Let ∆ α be the closed ball of radius α in R D×T 2 for the canonical norm • 2 . Since the feasible set Ω is closed, the intersection Ω ∩ ∆ α is bounded and closed, hence compact as a subset of a finite-dimensional Hilbert space. Further, the functional F (•|Z) + H(L•) being continuous over Ω ∩ ∆ α it reaches its minimum over this set at some point ( R, Ô)

∈ Ω ∩ ∆ α . Then, since (1, 0) ∈ Ω∩∆ α , it follows that F ( R, Ô|Z)+ H(L( R, Ô)) ≤ µ. From Eq. ( 17
), if (R, O) ∈ Ω\∆ α , then one has F (R, O|Z) + H(L(R, O)) > µ,
showing that the relative minimum ( R, Ô) on Ω ∩ ∆ α is to be a global minimum on the whole domain of feasible points Ω. To conclude the proof, we need to compute α satisfying ( 17), which follows from:

i) If O 1 > C µ/λ O , then H(L(R, O)) ≥ λ O O 1 > µ. ii) Since ∀z ≥ 0, lim p→∞ d KL (z|p) = ∞ there is a con- stant c = c(Z) such that D KL (Z|Q) > µ as soon as Q ∞ > c. Therefore, with c min u,d (ΦZ) (d) t , as soon as O 1 ≤ C and R ∞ > (c + C)/c , we have R × (ΦZ) + O ∞ ≥ c R ∞ -O ∞ ≥ c R ∞ -O 1 > c,
and hence F (R, O|Z) ≥ µ. In finite dimension, the equivalence of the norms yields constants

C O and C R such that, O 2 ≥ C O ⇒ O 1 ≥ C and R 2 ≥ C R ⇒ R ∞ ≥ (c+C)/c . Finally, α C 2 O + C 2 R
satisfies [START_REF] Parikh | Proximal algorithms[END_REF].

The regularizing functional in Eq. ( 12) not showing strict convexity, investigating the uniqueness of its minimizer requires further analysis, beyond the scope of this work. Nevertheless, leveraging the strict convexity of the standard KLD in Eq. ( 4

x ∈ Argmin x∈R m g(Bx) + h(x). ( 18
)
If g is strictly convex, then there exists ŷ ∈ R n , such that for any minimizer x of (18), ŷ = Bx.

Proposition 1. Let ( R, Ô) denote a solution of Problem [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF]. Let P {p = 0: For these times and territories, the estimated instant Poisson parameter is forced to cancel by the indicator function involved in Eq. ( 4). Then, restricting to the indices for which Z and ΦZ are not vanishing simultaneously, the data-fidelity term F is strictly convex and the uniqueness of P follows by the same reasoning.

D. Minimization via a proximal algorithm

From Corollary 1, the penalized Kullback-Leibler functional is convex, and, from Theorem 1, it has at least one minimizer. The 1 -norm and the indicator functions in the objective function make the overall functional non-smooth, preventing from the use of standard gradient descent. Minimizing [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF] thus requires more advanced tools, handling nondifferentiable functions, such as proximal algorithms.

1) Adjoint and proximal operators: Let L a bounded linear operator L, its adjoint operator is denoted L * and its operator norm is L op = sup{ LX 2 , X| 2 ≤ 1}. The convex conjugate of the proper, lower-semicontinuous function H : X → R, is defined as H * (Y ) = sup X∈X Y, X -H(X) [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], with X a Hilbert space. For τ ∈ R + the proximal operator of τ H at Y ∈ X reads by definition: prox τ H (Y ) arg min X∈X 1/2 X -Y 2 + τ H(X).

2) Primal-dual algorithmic scheme: In this work, we chose to implement the Chambolle-Pock primal-dual minimization scheme [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Making use of the compact reformulation in Eq. ( 13), we perform here an explicit particularization of the Chambolle-Pock algorithm detailed in Algorithm 1. From [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], Algorithm 1 converges toward a minimizer of Problem ( 13), if the descent parameters τ , σ satisfy

τ σ L 2 op < 1. (19)
Therefore, the actual and practical use of Algorithm 1 and enforcement of Condition ( 19), requires to evaluate, at least, an upper bound of L op .

Proposition 2. Let L be defined in Eq. ( 14), then

L 2 op ≤ max{λ 2 T D 2 2 op + λ S G 2 op + 1, λ 2 O }. (20)
Proof. By definition of the operator norm • op , one has

L(R, O) 2 2 (21) = λ 2 T D 2 R 2 2 + R 2 2 + λ 2 S GR 2 2 + λ 2 O O 2 2 ≤ (λ 2 T D 2 2 op + 1 + λ 2 S G 2 op ) R 2 2 + λ 2 O O 2 2 ≤ max{λ 2 T D 2 2 op + 1 + λ 2 S G 2 op , λ 2 O } (R, O) 2 2 .
In practice, equal descent step sizes are chosen, so that they saturate Condition [START_REF] Ma | Epidemiological parameters of coronavirus disease 2019: A pooled analysis of publicly reported individual data of 1155 cases from seven countries[END_REF], replacing the unknown L op by its upper bound derived at Proposition 2, leading to

τ = σ = 0.99 max{λ 2 T D 2 2 op + λ S G 2 op + 1, λ 2 O } . (22)
3) Closed-form expression of the proximal and adjoint operators: For proximal algorithms, such as Algorithm 1, the key ingredient to achieve an actual and fast implementation is to derive a closed-form expression for each proximal operator involved. Proximal operator of the regularizing function: The dual variable update in Algorithm 1 involves the proximal operator of the convex conjugate of the penalty H defined in [START_REF] Cai | Image restoration: Total variation, wavelet frames, and beyond[END_REF]. In practice, first, the proximal operator of H is computed, and then, the proximal operator ofits convex conjugate H * is derived applying Moreau's identity, which states that, for σ > 0, prox σH * (X) = X -σprox H/σ (X/σ). Then, function H/σ being fully separable, both over its four entries and over the components of each entry, its proximal operator can be computed component-wise. These computations involve two standard closed-form proximal operators: i) the proximal operator of the absolute value, known as the soft-thresholding, prox σ|•| (q) = max {|q| -σ, 0} q/|q|, and ii) the proximal operator of the indicator function, which amounts to a projection: prox ι ≥0 (q) = max {0, q}.

Algorithm 1 Primal-dual minimization of the penalized Kullback-Leibler [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF] for the estimation of reproduction numbers and outliers. Require: Infection counts: Z ∈ R D×T and (ΦZ) ∈ R D×T Choose descent parameters: τ, σ > 0, precision: ε Max. iterations: k max , length of window:

k smooth Initialization R [0] = Z, O [0] = 0, Q [0] = L(R [0] , O [0] ), R [0] = R [0] , O [0] = O [0] Φ 0 = F (R [0] , O [0] |Z) + H(L(R [0] , O [0] ) Ψ 0 = ε + 1, Ψ 0 = ε + 1, k = 0 while Ψ k ≥ and k < k max do
Update the dual, primal and auxiliary variables

Q [k+1] = prox σH * (Q [k] + σL(R [k] , O [k] )) (R [k+1] , O [k+1] ) = prox τ F ((R [k+1] , O [k+1] ) -τ L * Q [k+1] ) (R [k+1] , O [k+1] ) = 2(R [k+1] , O [k+1] ) -(R [k] , O [k] )
Computation of the convergence criteria

Φ k+1 = F (R [k+1] , O [k+1] |Z) + H(L(R [k+1] ; O [k+1] )) Ψ k+1 = |Φ k+1 -Φ k |/Φ k Ψ k+1 = max {Ψ : |k -k smooth | + ≤ ≤ k + 1} (23) k ← k + 1 end while
Adjoint of the linear operator: Updating the primal variable requires first to compute

L * : R D×(T -2) × R D×T × R E×T × R D×T → R D×T 2 ,
which can be expressed in terms of the adjoint operators of D 2 and G, via

L * : (Q 1 , Q 2 , Q 3 , Q 4 ) → (λ T D * 2 Q 1 + Q 2 + λ S G * Q 3 , λ O Q 4 ). ( 24
)
Proximal operator of the data-fidelity term: The extended negative log-likelihood F of Eq. ( 4) being separable over the times t ∈ [1, T ] and territories d ∈ [1, D], its proximal operator can be computed in a component-wise manner. Each term in (4) corresponds to the composition of a linear map with the standard Kullback-Leibler divergence d KL (z|•), whose proximal operator writes explicitly [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]:

prox τ dKL(z|•) (p) = (p -τ + (p -τ ) 2 + 4τ z)/2. ( 25
)
Proposition 3. For β (Φz) 2 + 1, a : (r, o) → r × (Φz) + o, prox τ f (•,•|z,(Φz)) (r, o) (26) = (0, 0), if z = (Φz) = 0 (r, o) -a * β -1 I -prox τ βdKL(z|•) a(r, o), otherwise.
from which derives straightforwardly the proximal operator of the global data-fidelity function F .

Proof. Since aa * = β > 0, the formula in Eq. ( 26) is obtained applying Proposition 3.4 of [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF].

4) Stopping criterion: The final step in the actual of an iterative algorithm is to choose a stopping criterion. Classically, the convergence of Algorithm 1 is quantified by the decrease of the normalized increments of the objective functional.

Fig. 1 shows that such increments (dashed blue curves) display large and erractic fluctuations, thus precluding their use to construct a reliable stopping criterion. Instead, and to provide a robust convergence criterion, we propose to make use of a nonlinearly smoothed normalized increments (solid black curve in Fig. 1), obtained as a maximum over a sliding window over the past k smooth iterations, as detailed in Eq. ( 23). Hence, in our concrete experiments, we set k smooth = 500 and Algorithm 1 is stopped whenever the smoothed increment reaches a precision ε, set to 10 -7 , or when the number of iterations exceeds a maximum of k max iterations, set to 10 7 .

IV. COVID-19 DATA

To show the relevance and efficiency of the proposed procedure for the estimation of the evolution of the reproduction number of the pandemic, two real Covid-19 datasets made available by national health authorities are used, different in nature and from different international repositories.

1) Dataset1: Country population level daily new infection counts: Johns Hopkins University2 provides access to the cumulated counts of daily new infections (together with deceased and recovered persons), at the entire population level, on a per country basis, for 200+ countries worldwide, since January 1st, 2020, hence, essentially since the earliest stages of the Covid-19 pandemic3 . In the present work, use will be made of the daily new infection counts only, {Z 1:T }, from Feb., 15th, 2020 until until August, 31st, 2021, hence T = 586 days of pandemic.

2) Dataset2: French Counties hospital Level daily new infection counts: Santé-Publique-France4 , the French national health authorities, provides several different datasets related to the Covid pandemic in France. In the present work, use will be made of the daily hospital-recorded new infection counts, across France on a per département-basis 5 , départements (or counties) consisting of the usual French administrative granularity. These new infections counts are stacked into a Z ∈ R T ×D data matrix, for the 101 French départements, and estimation of the reproduction number will be conducted jointly for all départements. Such data are however available only after March 19th, 2020 and analyzed until August, 31st, 2021, hence T = 531 days of pandemic.

V. ROBUST ESTIMATION FROM COVID-19 DATA

A. Hyperparameter tuning

One of the issue at hand for a practical use of Eq. 12 to estimate jointly R and O lies in the selection of the hyperparameters, λ T , λ S , and λ O that balance the contributions of the regularizations terms one against the others as well as against the data-model fidelity term. Such choices have significant impacts on the achieved estimates. Accurate parameter tuning is thus critical yet may prove time consuming notably when this selection needs to be made for each of the 200+ countries studied here, with possibly very different population sizes or pandemic intensities.

Automated data driven selection of hyperparameters in functional minimization has been addressed theoretically in several different settings (cf. e.g. [START_REF] Pascal | Automated data-driven selection of the hyperparameters for total-variation-based texture segmentation[END_REF]) and references therein), often relying on a Stein Unbiased Risk Estimator (SURE). This is working well only for large size data, which is however not the case for the Covid-19 pandemic.

Instead, in the present work, hyperparameter selection was driven by the following dimensional analysis. Inspecting Eq. 12 shows that the functional to minimize is covariant under the change (Z, O; λ T , λ S , λ O ) → (αZ, αO; αλ T , αλ S , λ O ) for any α > 0, thus leading to identical minima. This led us to propose to use one same set of hyperparameters (λ T , λ S , λ O ) valid for all countries, irrespective of their population sizes or pandemic intensities, by applying minimizations in 12 to αZ instead of Z, where α are arbitrary multiplicative constants, that may depend on each country. We found empirically that setting systematically α to the empirical standard deviation of the observed Z constitutes an efficient way to account both for the population size and for the severity of the pandemic. Precisely, λ T = 3.5 is chosen as a universal (for all countries) time constant, after inspections of numerous different countries over the whole pandemic period, as a valid trade-off between too small -that would lead to too much variability in the estimation of R and would prevent the use of R to actually assess the strength of the pandemic -and too largethat would lead to too few variability hence preventing the possibility of detecting changes in the pandemic dynamics that may follow some sanitary policy decisions (lockdown, vaccination). Further, λ O is set to 0.025 to enable the outlier term in the functional to be versatile enough to account for pseudo-seasonal effects modeled as sparse irrelevant measures. Finally, for different D = 96 territories, the connectivity of the graph for metropolitan France contains 475 edges, indicating that the time and space regularizations will have comparable impacts when 96λ time 2 × 475λ space , which frames the selection of λ space . Inspection of the results lead to choose λ space = 0.002, thus corresponding to a mild spatial regularization.

B. Temporal evolution for each country independently

Because Public Health Policies were not or weakly coordinated across countries, even for closely intertwined spaces such as Western European countries, population level daily new infection counts (Dataset1, JHU) are first analyzed inde- pendently for country, that is, without spatial regularization, or setting λ S ≡ 0 in [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF].

Estimates of the reproduction numbers and the outliers, RO , ÔO are thus obtained for each country by applying [START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF] to daily infection counts Z 1:T , after normalisation by standard deviations computed across the whole pandemic period, using the same hyperparameters for all countries : (λ T , λ S , λ O ) = (3.5, 0, 0.025).

To assess the relevance of the proposed one-step procedure, RO , ÔO are compared against estimates obtained from a twostep procedure developed in [START_REF] Abry | Spatial and temporal regularization to estimate COVID-19 reproduction number R (t): Promoting piecewise smoothness via convex optimization[END_REF], RnoO : First, a slidingmedian over a 7-day window is applied to Z 1:T , and values that depart from the window median by ±2.5 local (in-window) standard deviation are replaced by the window median, yielding estimates of outliers referred to as ÔnoO ; Second, an estimate RnoO of R is obtained by applying Eq. 12 with (λ T , λ S , λ O ) = (3.5, 0, +∞), to the a priori denoised infection counts after normalization by standard deviation.

Fig. 2 compares the estimates RMLE , ( RO , ÔO ) and ( RnoO , ÔnoO ), for France across the whole pandemic period, and also focuses on a recent narrow period of time to ease illustration and analysis. Fig. 2 shows that, for the pseudoseasonal effect related to new infection count under-reporting due to non-working days, both procedures essentially correctly detect and estimate the corresponding outliers. For the subsequent working-day over-reporting performed by the French public health authority, Fig. 2 clearly shows that the proposed one-step procedure correctly accounts for such small outliers, thus yielding a smooth and outlier-robust estimation of R, while mis-reporting are missed by the two-step procedure, that thus produces estimates with numerous discontinuities irrelevant in terms of pandemic monitoring and clearly driven by these mis-reporting. The period of early April 2021 is of great interest, corresponding to a long week-end effect, with a normal week-end, followed by an alternation of working and non working days, resulting in anomalously low infection counts, followed by high counts (misreport correction), during intertwined and following working days. The two-step procedure yields for that period totally irrelevant estimates, while the proposed one-step procedure remains satisfactorily insensitive to such count misreports. These pseudo-seasonal or long week-end effects are not specific to France and can be observed across numerous other countries. Fig. 3 further compares the estimates RMLE , ( RO , ÔO ) and ( RnoO , ÔnoO ) for several different countries. Fig. 3 clearly shows that i) RMLE constitute a very irregular estimate of R along time that cannot be used at all in practice ; ii) compared to ( RnoO , ÔnoO ), ( RO , ÔO ) better accounts for both outliers and pseudo-seasonal effects and yields smooth and regular along time and outlier-robust estimations of R. RO is thus far more likely to represent the actual, less biased and realistic assessment of the pandemic time evolution in a given country. Estimates ( RO , ÔO ) are updated on a daily basis for 200+ countries across the world, and made publicly available via interactive and animated maps (cf. Section VI).

C. Joint space-time Evolution for metropolitan France

Let us now consider multivariate daily infections counts, Z ∈ R T ×D , for the D = 96 départments of metropolitan France, considered as the vertices of the graph G in Eq. 12, with edges between départments sharing a terrestrial boundary. Estimates ( RO,S , ÔO,S ) ∈ R T ×D are obtained by applying (not accounting for outliers but with spatial regularization, magenta), RO,noS (accounting for outliers but without spatial regularization, blue), RnoO,noS (not accounting for outliers and without spatial regularization, cyan).

Eq. 12 to Z ∈ R T ×D , after standardization independently per départment, with (λ T , λ S , λ O ) = (3.5, 0.002, 0.025) (cf. Section V-A).

To illustrate the relevance of ( RO,S , ÔO,S ), they are compared against estimates obtained without spatial regularization and without accounting for outliers ( RnoO,noS , ÔnoO,noS ) ((λ T , λ S , λ O ) = (3.5, 0, +∞)), without spatial regularization but accounting for outliers ( RO,noS , ÔO,noS ) ((λ T , λ S , λ O ) = (3.5, 0, 0.025)), with spatial regularization but without accounting for outliers ( RnoO,S , ÔnoO,S ) ((λ T , λ S , λ O ) = (3.5, 0.002, +∞)).

These estimates are first compared as functions of time for two départments in Fig. 4, which shows that estimates that do not account for outliers are far too irregular for being useful in pandemic monitoring, even with spatial regularization. To the converse, RO,S that results from both accounting for outliers and spatial regularization yields very regular estimates likely to reflect a relevant assessment of the pandemic. The proposed estimation procedure is of particular interest when applied to territories with granularity level as is the case for the French départments (with on average slightly less than one million inhabitants) or during low activity phases of the pandemic, such as early summer 2020 (days 95 to 155). Spatial regularization is of particular relevance for theses two chosen départments. Indeed, départment Paris (Id=75) corresponds to the city of Paris, thus a small territory yet with massive connection to the surrounding départments : estimates of R within Paris, can thus not notably differ from those of these so-called Paris-Crown départments. The same holds for départment Rhône (Id=69), of very small surface (for historical reasons) yet acting, with the city of Lyon as a massive transit hub, for a number of surrounding départments.

Further, the estimates of R are for continental France are reported in Fig. 5) for three different days. Estimates of RO,S with spatial regularization and accounting for outliers (right most plots) permits a far clearer assessment of the status of the pandemic across metropolitan France. Estimates RnoO,S and RnoO,noS show far too much spatial variability to be realistic estimates. Estimates RnoO,S improve spatial regularity, but still show significant variability induced by data low quality. Estimates RO,S , with spatial regularization

VI. CONCLUSION

The present work developed and assessed an inverse problem-type procedure to estimate the spatio-temporal evolution of the Covid-19 reproduction number, R, robust to the low quality of the pandemic data. The devised procedure relies first on a detailed analysis of the data quality that suggested that most count misreports, be they mild and seasonal (sundays, week-end, days-off) or large and random (report failure), can be efficiently accounted for as sparse outliers. Second, it is based on devising carefully a functional that balances a datamodel fidelity term (Poisson distribution and outliers), and regularity in time and space properties that R needs to fulfill to be of actual relevance in practical pandemic monitoring. Fast and efficient algorithms were devised to minimize this functional and their convergence was studied. Applied to real Covid-19, the procedure was shown to yield meaningful estimates of R, that thus provides a relevant assessment of the spatio-temporal evolution of the pandemic that can be involved as part of a decision strategy for designing Covid-19 counter measures. Indeed, at a current time, the procedure outputs estimates of R for the entire pandemic period, thus permitting to assess a posteriori if and how the implementation of a sanitary policy decision (lockdown, curfew,. . . ) impacted the evolution of the pandemic. Also and importantly, while the proposed procedure is not forecasting R, it is actually implementing a nowcasting approach: Not only estimates of R today are provided, but the imposed piecewise smoothness constraints provide epidemiologist with a local trend around at current time indicating whether the pandemic is progressing or regressing. Spatial regularization also permits to robustly assess homogeneity or heterogeneity of the pandemic across related territories, thus permitting to decide on local or global (nation-wide) measures.

Estimates are updated automatically on a daily basis for 200+ countries worldwide and for several related territories (Counties in France, States in the USA). They are made publicly available via downloadable text files or interactive and animated maps 6 . Procedures will be made publicly available.

Future investigations to improve estimation and monitoring include the construction of confidence intervals for the estimates and automated data-driven selection of the hyperparameters, using adaptive strategy that adjust the variations and phases of the pandemic.

We believe that the present work constitutes a significant step forward toward the practical use of the proposed procedure permits an actual real-time and on-the-fly monitoring of the Covid-19 pandemics. These tools are ready to be put at work immediately at the outbreak of potential future pandemics. Finally, the automated production of these daily up-dates is also intended as a scientific tool to favor transdisciplinary scientific work against Covid-19 impact on the society at large.

t

  , t ∈ {1, . . . T }, d ∈ {1, . . . D}}, then P is unique.Proof. In the trivial situation when the daily counts are uniformly vanishing, i.e., Z = ΦZ ≡ 0, the unique minimizer coincides the unique feasible point and consists in estimated reproduction numbers and outliers equal to zero at any time and in any territory. Otherwise, two situations are to be considered (i) If Z and ΦZ never vanish simultaneously: In this caseF (R, O|Z) = D KL (Z|B(R, O)), with B : (R, O) → R × (ΦZ) + O.Further, the standard Kullback-Leibler divergence g(P ) = D KL (Z|P ) is strictly convex. Hence, setting h = H(L•), Proposition 1 yields the result. (ii) If there exist some times u and territories d s.t. Z

Fig. 1 :

 1 Fig. 1: Stopping criterion. (Normalized) Objective function increments and non linearly smoothed increments as functions of the number of iterations in Algorithm 1 (computed from new infections in France).

Fig. 2 :

 2 Fig. 2: France (six first months in 2021). Top: Z and Z -ÔO , (black and red) ; Middle: RMLE , RO , RnoO (black, red and blue) ; Bottom: ÔO , ÔnoO (red and blue).

Fig. 3 :

 3 Fig. 3: Different countries. ,Top: Z and Z -ÔO (black and red) ; Middle: ( RMLE , RO , RnoO ) (black, red and blue) ; Bottom: ( ÔO , ÔnoO ) (red and blue). Estimates ( RO , ÔO ) are updated on a daily basis for 200+ countries across the world, and made publicly available via interactive and animated maps available (cf. Section VI).

Fig. 4 :

 4 Fig.4: Estimates of R t for départment Paris (Id=75) (right) and départment Rhône (Id=69) (left): RO,S (accounting for outliers and with spatial regularization, red), RnoO,S (not accounting for outliers but with spatial regularization, magenta), RO,noS (accounting for outliers but without spatial regularization, blue), RnoO,noS (not accounting for outliers and without spatial regularization, cyan).

Fig. 5 :

 5 Fig. 5: Compared estimates of R (d) t across metropolitan France for three different days along the pandemic.

  with [1, T ] the integers {1, 2, . . . , T }. The pandemic spread on each territory d is thus characterized by both a reproduction number time series and an outlier time series {(R

		(d) t , O t ), t ∈ (d)
	[1, T ]}, that must thus be estimated jointly for all territories.
	These considerations lead to a first key contribution of
	the present work: the modeling of corrupted daily infection
	counts across territories by extending Model (1) to multivariate
	daily new infection Poisson counts, each with an instantaneous
	Poisson parameter p	(d) t , with (ΦZ)

https://coronavirus.jhu.edu/

https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse covid 19 time series/

[START_REF] Brauer | Mathematical models in epidemiology[END_REF] https://www.santepubliquefrance.fr

/ 5 https://www.data.gouv.fr/fr/datasets/

Animated and interactive maps available at perso.ens-lyon.fr/patrice.abry/, http://barthes.enssib.fr/coronavirus/cartes/Rmonde/ and http://barthes.enssib. fr/coronavirus/cartes/RFrance/ are conceived and deployed in collaboration with Eric Guichard (ENSSIB, Triangle ENS de Lyon, France), who is here gratefully acknowledged.