# One-Bit Total Variation Denoising over Networks with Applications to
Partially Observed Epidemics

Claire DONNAT<EMAIL_ADDRESS>Department of Statistics, University of
Chicago,
5747 S Ellis Ave, Chicago, Illinois, 60637, United States of America  Olga
KLOPP<EMAIL_ADDRESS>ESSEC Business School
3 Av. Bernard Hirsch, 95000 Cergy, France
CREST, ENSAE, Institut Polytechnique de Paris,
5 Av. Le Chatelier, 91120 Palaiseau, France  Nicolas VERZELEN
<EMAIL_ADDRESS>MISTEA, Université Montpellier, INRAE, Institut
Agro
2 Pl. Pierre Viala, 34000 Montpellier, France

###### Abstract

This paper introduces a novel approach for epidemic nowcasting and forecasting
over networks using total variation (TV) denoising, a method inspired by
classical signal processing techniques. Considering a network that models a
population as a set of $n$ nodes characterised by their infection statuses
$Y_{i}$ and that represents contacts as edges, we prove the consistency of
graph-TV denoising for estimating the underlying infection probabilities
$\\{p_{i}\\}_{i\in\\{1,\cdots,n\\}}$ in the presence of Bernoulli noise. Our
results provide an important extension of existing bounds derived in the
Gaussian case to the study of binary variables — an approach hereafter
referred to as one-bit total variation denoising. The methodology is further
extended to handle incomplete observations, thereby expanding its relevance to
various real-world situations where observations over the full graph may not
be accessible. Focusing on the context of epidemics, we establish that one-bit
total variation denoising enhances both nowcasting and forecasting accuracy in
networks, as further evidenced by comprehensive numerical experiments and two
real-world examples. The contributions of this paper lie in its theoretical
developments, particularly in addressing the incomplete data case, thereby
paving the way for more precise epidemic modelling and enhanced surveillance
strategies in practical settings.

###### keywords:

High-dimensional statistics; Graph total-variation; Graph-Trend Filtering;
Gamma sparsity; Epidemic Forecasting; Epidemic Nowcasting.

##  1 Introduction

We consider the problem of signal denoising over networks and apply our
results to the analysis of infectious disease transmission over contact
networks. The spread of epidemics over networks is an important area of
research with implications for public health (Ottaviano et al., 2018; Sanchez
et al., 2022), epidemiology (Danon et al., 2011; Keeling & Eames, 2005a;
Moreno et al., 2002; Shirley & Rushton, 2005; Spricer & Britton, 2019), and
network analysis (Kim et al., 2021; Wan et al., 2014). While the study of
disease propagation through interconnected populations had already gained
significant attention over the past decade (Keeling, 2005; Newman, 2002a;
Pastor-Satorras et al., 2015a; Pellis et al., 2015), the recent emergence of
new infectious diseases, such as COVID-19, further underscores the importance
of developing comprehensive models that can capture the intricate interplay
between epidemic transmission processes and the underlying network structures
through which they spread (Keeling & Eames, 2005a; Manzo, 2020).

Studies of epidemic processes typically rely on compartmental models (Kermack
& McKendrick, 1927; Tolles & Luong, 2020), which split a population of size
$n$ into different compartments (or groups) based on their disease status, and
specify how agents transition between disease states (e.g., Susceptible,
Infectious and Recovered). These approaches can be further refined into two
categories, depending on how they account for randomness and uncertainty in
the disease transmission process: deterministic, versus stochastic epidemic
models. Deterministic models (Bowman et al., 2005; Carcione et al., 2020;
Eikenberry et al., 2020; Kumar et al., 2021; Mubarak et al., 2021) assume that
the disease transmission process can be precisely described using mathematical
equations with fixed parameters. These models typically use differential
equations to describe how the numbers of infected, susceptible, and recovered
individuals change over time, and do not capture individual-level variability
in the disease transmission process. Although deterministic models have
provided valuable insights into the dynamics of epidemics, they often fail to
capture the inherent randomness and unpredictability associated with real-
world outbreaks. By contrast, the consideration of variability recognises that
epidemic processes are influenced by random events, chance encounters, and
unforeseen interactions among individuals. This acknowledgment prompts the
need for stochastic models that can explicitly incorporate randomness into the
modeling framework (Ball et al., 2019).

Instead of using differential equations, stochastic models often rely on
techniques such as agent-based modelling (Hackl & Dubernet, 2019; Perez &
Dragicevic, 2009; Siettos et al., 2015; Teweldemedhin et al., 2004) or
stochastic differential equations on compartmental models (Andersson &
Britton, 2012; Britton, 2010; Bu et al., 2022; Huang et al., 2022; Morsomme &
Xu, 2022) to simulate individual interactions and disease transmission events.
In this framework, each interaction or transmission event is subject to
chance. The fluxes between compartments therefore follow a probability
distribution, rather than a deterministic equation. Although these approaches
are designed to better capture the uncertainty inherent in epidemic spreads,
most of them only typically incorporate randomness in terms of transmission
events: given a contact between an infected and susceptible individual, the
transmission is subject to chance, and typically modelled as a Bernouilli
variable. Specifically, these models do not attempt to remedy another major
limitation of compartmental models: the uniform mixing assumption. This
assumption implies that every infected individual in the population can
potentially transmit the disease to any susceptible individual in the network,
which is not always realistic (Keeling & Eames, 2005a). Such a simplification
overlooks critical dynamics in epidemic spread on networks, particularly in
large populations where the uniform mixing assumption is known to fail. These
approximations can significantly affect the accuracy of models that attempt to
capture how diseases spread through complex networks.

A possible solution to address this limitation is to model transmissions
within a given contact network, thereby capturing more realistically possible
transmission events. The structure of contact networks can in fact known to
significantly impact disease spread (Ganesh et al., 2005). Indeed, while the
standard epidemic models described in the previous paragraph correspond to
fully connected networks, accounting for a more realistic structure of
interactions and connections between individuals may lead to more accurate
predictions and more effective control strategies (Chakrabarti et al., 2008;
Wang et al., 2013). For example, Pastor-Satorras & Vespignani (2001) have
shown that the structure of the contact network greatly influences the
epidemic threshold and disease prevalence. Research on temporal networks
(Holme & Saramäki, 2012) has shown that the timing of interactions plays a
crucial role in shaping epidemic dynamics by affecting disease
transmissibility. Epidemic models on graphs have also highlighted the
importance of super-spreaders, that is, individuals or nodes with a
disproportionately high potential for transmission. Finally, Keeling & Eames
(2005b) introduced the notion of “assortativity” in networks, where highly
connected nodes preferentially interact with each other. This assortative
mixing has significant implications for the effective targeting interventions
and controlling epidemics.

Another key limitation of traditional deterministic and stochastic
differential equation models is the assumption that the current state of the
epidemic is known: all agents’ infectious, susceptible, or recovered states
are observed. Knowledge of the current state of the epidemic is, in fact,
essential for making accurate predictions about its future course, as small
inaccuracies can lead to substantial changes in epidemic trajectories over
time (Siegenfeld et al., 2020). However, accurately assessing the current
situation of an epidemic, a process often referred to as nowcasting, is
notoriously challenging (Desai et al., 2019; McGough et al., 2020; Wu et al.,
2021), due to the availability of quality data and the delays in reporting of
current cases (Chakraborty et al., 2022; Rosenfeld & Tibshirani, 2021). In
many situations, complete information about individuals’ health status cannot
be disclosed due to privacy concerns. It thus becomes essential to develop
methods that are adapted to the practical constraints of real data and respect
individuals’ rights while still providing valuable insights into epidemic
spread.

In the present paper, we propose a new approach to nowcasting and forecasting
of epidemic spreads on a network inspired by total-variation signal denoising.
Signal denoising based on the total variation (TV) penalty is a widely used
technique in signal processing and image reconstruction that aims to remove
noise from an observed signal while preserving important spatial or structural
features. This approach uses the TV penalty, which measures the variation in
the signal over a network. Letting $p_{i},i\in\\{1,\cdots,n\\}$ be a signal
over a network $\mathcal{G}$ on $n$ nodes and $E$ edges, the total variation
penalty can be written as: $\sum_{(i,j)\in E}|p_{i}-p_{j}|.$ By penalising
differences in the signal over neighbouring vertices, this term acts as a
regularization penalty that enforces smoothness in the reconstructed signal.
TV regularization is typically formulated as an optimisation problem, where
the goal is to minimise the TV of the signal subject to a fidelity constraint
in the form of a data-fitting loss, ensuring the reconstructed signal remains
close to the original data. A more extensive review of Total-variation
denoising is provided in Section 7.

In this paper, we study TV denoising in the setting of Bernoulli observations
and show that this approach improves nowcasting estimates and prediction
accuracy of dynamic behavior of epidemic processes over simple compartmental
approaches. To address the challenges of epidemic nowcasting, we also extend
our approach to the case of incomplete observations which allows a better
assessment of the disease’s spread and can help implement more targeted
surveillance and control measures.

### Contributions

Our contributions are threefold and summarised as follows. First, we prove
that the total variation denoiser that we propose – thereafter referred to as
one-bit TV denoiser — is consistent in the case of Bernoulli noise (see
Theorem 2.4). These results extend existing bounds that have been obtained in
the literature in the case of Gaussian noise (Hütter & Rigollet, 2016). In
particular, we obtain new bounds on the risk of TV denoiser measured in terms
of the $\ell_{1}$-norm. Second, we extend our analysis to the challenging case
of missing observations and demonstrate good performance of the one-bit TV
denoiser (see Theorem 4.1). To the best of our knowledge, such results are new
in the literature. Finally, we show how the one-bit total variation denoiser
can be applied to improve the prediction of virus spread dynamics on networks
(see Section 3). Our results are supported by extensive numerical studies that
confirm the advantages of the proposed method, as well as two real-data
experiments.

### Notations

For any $k$, $\mathbf{1}_{k}$ denotes the $k$-dimensional vector of ones and
$\mathbf{0}_{k}$ denotes the $k$-dimensional vector of zeros. For an integer
$k\in\mathbb{N}$, we write $[k]=\\{1,\dots,k\\}$. We denote by $A^{\dagger}$
the Moore-Penrose pseudo-inverse of a matrix $A$ and by $\|A\|$ its operator
norm. Let $S\subset[n]$ be a subset of the nodes in $[n]$, we will use the
notation $d_{\max}(S)$ to denote the maximum degree of nodes $i$ in subset
$S$.

##  2 Signal denoising

Consider a network where each node $i\in V$ has an underlying state
$y_{i}\in\\{0,1\\}$ where $0$ and $1$ may respectively represent, for example,
healthy and infected states. The probability that person $i$ is in the
infected state (or the proportion of the infected subpopulation at node $i$)
is indicated by $p^{*}_{i}$, that is, $p^{*}_{i}={\text{pr}}(y_{i}=1)$. Let
$p^{*}=(p^{*}_{i})_{i=1,\dots n}$ denote the vector that captures the infected
states of an interconnected population consisting of $n$ nodes. The state of
the network can thus be represented by the following “signal + noise” model:

$Y=p^{*}+\xi$ (1)

where $Y=(y_{i})_{i\in[n]}$ are our observations and $\xi$ is a vector of
centred Bernoulli noise. Given one vector of observations $Y$, we consider the
problem of estimating the underlying $(p^{*}_{i})_{i\in[n]}$.

Let $\mathcal{G}=(V,E)$ be the underlying undirected connected graph with
vertex set $V$ and edge set $E$. We set the cardinalities of both sets as
$|V|=n$ and $|E|=m$. We denote by $A=(a_{ij})_{(i,j)\in[n]\times[n]}$ the
graph’s adjacency matrix. It will be convenient for us to represent a graph by
its edge-vertex incidence matrix $D\in\\{-1,0,1\\}^{m\times n}$. To each edge
$e=(i,j)\in E$ corresponds a row $D_{e}$ of $D$ where the $k$th entry
$D_{e,k}$ of $D_{e}$ is given by

$D_{e,k}=\begin{cases}1&\text{if }k=\min(i,j)\\\ -1&\text{if }k=\max(i,j)\\\
0&\text{otherwise.}\\\ \end{cases}$

Note that that $L=D^{T}D=diag(A\mathbf{1}_{n})-A$ is the unnormalized
Laplacian of the graph $\mathcal{G}$. Here $diag(A\mathbf{1}_{n})$ is the
diagonal matrix with $j$th diagonal element given by the degree of vertex $j$.

The TV denoiser $\widehat{p}$ associated to $G$ is any solution of the
following minimization problem:

$\widehat{p}\in\underset{p\in[0,1]^{n}}{\mathop{\rm
arg\,min}}\left\\{\dfrac{1}{n}\overset{n}{\underset{i=1}{\sum}}\left(y_{i}-p_{i}\right)^{2}+\lambda\|Dp\|_{1}\right\\},$
(2)

where $\lambda>0$ is a regularization parameter and the TV penalty
$\|Dp\|_{1}$ is the convex relaxation of $\|Dp\|_{0}$, which corresponds to
the number of times $p$ changes values along the edges of the graph $G$. Note
that (2) is a convex problem that may be solved efficiently (Arnold &
Tibshirani, 2016). We start by defining the compatibility factor and the
inverse scaling factor that will be crucial in characterising the performance
of the TV denoiser:

######  Definition 2.1 (Compatibility factor).

Let $D\in\\{-1,0,1\\}^{m\times n}$ be the incidence matrix. The compatibility
factor of a set $T\subset[m]$ is defined as

$\kappa_{T}=\kappa_{T}(D)=\underset{p\in[0,1]^{n}}{\inf}\dfrac{\sqrt{|T|}\,\|p\|_{2}}{\left\|\left(Dp\right)_{T}\right\|_{1}}\quad\text{for}\quad
T\not=\emptyset.$

######  Definition 2.2 (Inverse scaling factor).

Let $S=D^{\dagger}=[s_{1},\dots,s_{m}]$. The inverse scaling factor of $D$ is
defined as

$\rho=\rho(D)=\underset{j\in[m]}{\max}\|s_{j}\|_{2}.$

The following result shown in Hütter & Rigollet (2016) provides a bound on
both $\rho$ and $\kappa_{T}$:

######  Proposition 2.3.

Let $D$ be the incidence matrix of a connected graph $G$ with maximal degree
$d_{\max}$ and $\emptyset\not=T\subseteq E$. Let
$0=\lambda_{1}<\lambda_{2}\leq\dots\leq\lambda_{n}$ be the eigenvalues of the
Laplacian $L$. Then,

$\kappa_{T}\geq\dfrac{1}{2\min\\{\sqrt{d_{\max}},\sqrt{|T|}\\}}\quad\text{and}\quad\rho\leq\dfrac{\sqrt{2}}{\lambda_{2}}.$

Equipped with these notations, we establish the following bound on the risk of
TV denoiser (2):

######  Theorem 2.4 (Risk bound for one-bit TV denoising).

Fix $\delta\in[0,1]$, $T\subset[m]$ and assume that $G$ is connected. Define
the regularization parameter

$\lambda=\frac{\sqrt{2}\rho}{n}\log\left(\frac{4n^{2}}{\delta}\right)\enspace.$
(3)

Then, the TV denoiser $\widehat{p}$ defined in (2) satisfies

$\begin{split}{\|\widehat{p}-p^{*}\|^{2}_{2}}\leq\dfrac{16\rho^{2}|T|\log^{2}\left(\frac{4n^{2}}{\delta}\right)}{\kappa^{2}_{T}}+4\sqrt{2}\rho\left\|\left(Dp^{*}\right)_{T^{c}}\right\|_{1}\log\left(\frac{4n^{2}}{\delta}\right)+\frac{4\|p^{*}\|_{0}\log\left(\tfrac{4}{\delta}\right)}{n}\end{split}$
(4)

with probability at least $1-\delta$.

The proof of Theorem 2.4 is provided in Appendix A. As we now show, inequality
(4) allows trading off $|T|$ and
$\left\|\left(Dp^{*}\right)_{T^{c}}\right\|_{1}$. Taking $T={\rm
supp}(Dp^{*})$ , we obtain the following bound on the estimation error:

${\|\widehat{p}-p^{*}\|^{2}_{2}}\leq\dfrac{16\rho^{2}\|Dp^{*}\|_{0}\log^{2}\left(\frac{4n^{2}}{\delta}\right)}{\kappa^{2}_{T}}+\frac{4\|p^{*}\|_{0}\log\left(\tfrac{4}{\delta}\right)}{n}.$

In particular, in the context of epidemic modelling, we can assume that in the
first stages of the epidemic, the number of infected persons is relatively
small and that $p^{*}$ is sparse. Let $\mathcal{S}={\rm supp}(p^{*})$ and
$s=|\mathcal{S}|=\|p^{*}\|_{0}$. Then, using $\left|{\rm
supp}(Dp^{*})\right|\leq sd_{\max}(\mathcal{S})$, the previous inequality can
be simplified to:

${\|\widehat{p}-p^{*}\|^{2}_{2}}\leq\dfrac{16\rho^{2}sd_{\max}(\mathcal{S})\log^{2}\left(\frac{4n^{2}}{\delta}\right)}{\kappa^{2}_{T}}+\frac{4s\log\left(\tfrac{4}{\delta}\right)}{n}.$

On the other hand, taking $T$ to be the empty set in (4), we obtain the
alternative error bound:

${\|\widehat{p}-p^{*}\|^{2}_{2}}\leq
4\sqrt{2}\rho\left\|Dp^{*}\right\|_{1}\log\left(\frac{4n^{2}}{\delta}\right)+\frac{4s\log\left(\tfrac{4}{\delta}\right)}{n}.$

The next proposition provides a risk bound in $\ell_{1}$ norm.

######  Proposition 2.5.

Fix $\delta\in[0,1]$, $T\in[m]$ and assume that $G$ is connected. Define the
regularization parameter

$\lambda=\frac{\sqrt{2}\rho}{n}\log\left(\frac{4n^{2}}{\delta}\right)$

Then, the TV denoiser $\widehat{p}$ defined in (2) satisfies

$\begin{split}{\|\widehat{p}-p^{*}\|_{1}}&\leq\dfrac{4\rho\sqrt{s|T}|}{\kappa_{T}}\log\left(\frac{4n^{2}}{\delta}\right)+2s\sqrt{\frac{\log\left(\frac{4}{\delta}\right)}{n}}\\\
&+3\sqrt{\rho\,s\left\|\left(Dp^{*}\right)_{T^{c}}\right\|_{1}\log\left(\frac{4n^{2}}{\delta}\right)}+\sqrt{2}\rho
sd_{\max}(\mathcal{S})\log\left(\frac{4n^{2}}{\delta}\right)\enspace,\end{split}$
(5)

with probability at least $1-\delta$. Here, we recall that $s$ stands for
$\|p^{*}\|_{0}$.

The proof of Proposition 2.5 is provided in Appendix B. In Proposition 2.5, we
can again trade off $|T|$ and
$\left\|\left(Dp^{*}\right)_{T^{c}}\right\|_{1}$. Taking $T={\rm
supp}(Dp^{*})$ and using the fact that $|{\rm supp}(Dp^{*})|\leq
sd_{\max}(S)$, the previous bound becomes:

${\|\widehat{p}-p^{*}\|_{1}}\leq\rho
s\left(\dfrac{4\sqrt{d_{\max}(\mathcal{S})}}{\kappa_{T}}+2d_{\max}(\mathcal{S})\right)\log\left(\frac{4n^{2}}{\delta}\right)+2s\sqrt{\frac{\log\left(\frac{4}{\delta}\right)}{n}}.$

On the other hand, taking $T$ to be the empty set, we obtain:

${\|\widehat{p}-p^{*}\|_{1}}\leq
3\sqrt{\rho\,s\left\|\left(Dp^{*}\right)\right\|_{1}\log\left(\frac{4n^{2}}{\delta}\right)}+\sqrt{2}s\left(\rho\,d_{\max}(\mathcal{S})\log\left(\frac{4n^{2}}{\delta}\right)+\sqrt{\frac{2\log\left(\frac{4}{\delta}\right)}{n}}\right).$

The errors terms given by (4) and (5) crucially depend on the inverse scaling
factor $\rho$ and the compatibility factor $\kappa_{T}$. Hütter & Rigollet
(2016) provides upper bounds for these quantities for different graphs. Using
$\lesssim$ for an inequality where the left-hand side is bounded by the right-
hand side up to a numerical constant, we specialize the upper bound on the
risk of the TV denoiser to specific graphs. All the following bounds ( case 1
through 4) hold with probability at least $1-4/n$.

Case 1: $2D-$grids. For the 2D grid, $\rho\lesssim\sqrt{\log n}$ (see Hütter &
Rigollet (2016)). The $\ell_{2}$ and $\ell_{1}$ bounds on the risk of
Proposition 2.5 and Theorem 2.4 thus yield that:

${\|\widehat{p}-p^{*}}\|_{2}\lesssim\left\\{\sqrt{\left\|Dp^{*}\right\|_{1}+\frac{s}{n\log^{1/2}(n)}}\;\right\\}\log^{3/4}\left(n\right)\quad\text{and}$
${\|\widehat{p}-p^{*}\|_{1}}\lesssim\left\\{\sqrt{\frac{s\left(\left\|Dp^{*}\right\|_{1}\right)}{\log^{3/2}(n)}+s}\right\\}\;\log^{3/2}\left(n\right).$

Case 2: Complete graphs. For a complete graph, $\rho\lesssim 1/n$ and
$\kappa_{T}\gtrsim 1/\sqrt{n}$ (Hütter & Rigollet, 2016). This implies:

${\|\widehat{p}-p^{*}\|_{2}}\lesssim\sqrt{\dfrac{\left\|Dp^{*}\right\|_{1}+s}{n}}\;\log^{1/2}\left(n\right)\quad\text{and}$
${\|\widehat{p}-p^{*}\|_{1}}\lesssim\left(\sqrt{\dfrac{s\left\|Dp^{*}\right\|_{1}}{n\log(n)}}+s\right)\log\left(n\right).$

Case 3: Star graphs. For a star graph, $\rho\lesssim 1$ and $\kappa_{T}\gtrsim
1/\sqrt{|T|}$, so that:

${\|\widehat{p}-p^{*}\|_{2}}\lesssim\left\\{\sqrt{\left\|Dp^{*}\right\|_{1}+\frac{s}{n}}\right\\}\log^{1/2}\left(n\right)\quad\text{and}$
${\|\widehat{p}-p^{*}\|_{1}}\lesssim\left(\sqrt{\frac{s\left\|Dp^{*}\right\|_{1}}{\log(n)}}+sd_{\max}(S)\right)\log\left(n\right)$

Case 4: Graphs with spectral gaps. If the graph Laplacian has a spectral gap,
i.e., there exists $\nu>0$ such that the second smallest eigenvalue of the
Laplacian is larger than $\nu$, then $\rho\leq 2/\nu$ (see Proposition 14 in
Hütter & Rigollet (2016)). By Proposition 2.3, we also know that
$\frac{1}{\kappa_{T}^{2}}\leq 2{\min(|T|,d_{\max}})$. In particular, for a
random $d_{n}$-regular graph with $d_{n}=d_{0}(\log n)^{\beta}$ for some
$\beta>0$ and $d_{0}>1$ and for Erdös-Rényi graphs $\mathcal{G}(n,p)$ with
$p=d_{n}/n$, we obtain:

${\|\widehat{p}-p^{*}\|_{2}}\lesssim\sqrt{\frac{d_{n}\|Dp^{*}\|_{0}}{\nu^{2}}+\frac{s}{n\log^{1/2}(n)}}\log\left(n\right)\quad\text{and}$
${\|\widehat{p}-p^{*}\|_{1}}\lesssim\left(\sqrt{\frac{sd_{n}\|Dp^{*}\|_{0}}{\nu^{2}}}+\frac{sd_{n}}{\nu}\right)\log\left(n\right)$

Discussion. These bounds allow to shed more light on the interplay between the
size of the support ${\rm supp}(Dp^{*})$, the size of the epidemic
$\|p^{*}\|_{0}$, as well as the topology of the graph. At a high level — and
as confirmed by our experiments in Section 5—, we expect our total variation
denoiser to be better suited to situations where the cardinality of the
signal’s support is small compared to the size of the graph
($\|p^{*}\|_{0}=s<n$), and the maximal degree of the graph remains of constant
order, so that $|{\rm supp}(Dp^{*})|\leq s{d_{\max}}(V)$ is also controlled.
This means for instance that the denoiser is likely to provide a better
control of the $\ell_{1}$-error in the $2D$-grid (where the $\ell_{1}$ error
is of the order of $s\log^{3/2}(n)$ compared to star graph, where the error
could be in the order of $sn\log(n)$.

##  3  Virus spread dynamics on network

We now apply our one-bit total variation denoiser to the analysis of epidemic
spreads over networks. As highlighted in the introduction, numerous
compartmental models have been suggested to depict the spread of epidemic
processes over networks. We refer the reader to the recent surveys by Draief &
Massoulié (2009); Pastor-Satorras et al. (2015b) or Paré et al. (2020a) for a
more extensive review of these methods. We first propose deploying our model
to the simple Susceptible - Infectious - Susceptible (SIS), before showing how
our approach can be used in a wider variety of epidemic models.

### SIS Model

We start by considering the networked version of the susceptible - infectious
- susceptible (SIS) compartmental model, originally introduced in Kermack &
McKendrick (1927). This widely used model has been indeed adapted for network
settings in various subsequent studies (Ahn & Hassibi, 2013; Paré et al.,
2020b), as detailed in the subsequent paragraph. We take the vector
$(p_{i})_{i\in[n]}$ to represent the probabilities of each node being
infected, with each entry corresponding to an individual node’s infection
probability. In the networked setting (Paré et al., 2020a), the evolution of
the $p_{i}$s is governed by the following differential equation:

$\dot{p}_{i}(t)=(1-p_{i}(t))\beta_{i}\sum_{j=1}^{n}\omega_{ij}p_{j}(t)-\gamma_{i}p_{i}(t).$
(6)

Here $\omega_{ij}$ quantifies the strength of the connection from node $j$ to
node $i$, $\beta_{i}>0$ is the infection (or susceptibility) rate and
$\gamma_{i}$ is the healing rate. The infection rate $\beta_{i}$ quantifies
the rate at which a susceptible individual (node $i$) in the population is
infected through contacts with infected individuals (denoted here by the sum
$\sum_{j}w_{ij}p_{j}$). In many cases, the susceptibility to a virus may not
be uniform across the population. Susceptibility to a virus often varies
across a population due to factors such as age, pre-existing health
conditions, or immunity from previous exposure or vaccination. To account for
this heterogeneity in infection rates, one solution consists in identifying
relevant subgroups within the population and determining the infection rate
for each subgroup.

The healing rate $\gamma_{i}$ models the rate at which an infected individual
$i$ transitions between infected and recovered states. When $\gamma_{i}$ is
assumed to be identical across all individuals
($\gamma_{i}=\gamma,\quad\forall i=1,\cdots,n$), it can be understood as the
proportion of infected individuals who recover from the virus in a given time
window $\Delta t$. In this case, the healing rate can be estimated by tracking
the number of individuals who transition from infected to recovered status
over time, and calculated as the number of recoveries divided by the number of
confirmed cases during that time window $\Delta t$. Usually, the time window
for calculating the recovery rate is determined by the average duration of
infection. In general, to estimate the healing rate, we need information on
the typical progression of the disease, including the duration of the
infection, the severity of symptoms, and the factors that contribute to
recovery. Note that the healing rate can vary depending on factors such as
age, underlying health conditions, and the availability and effectiveness of
treatments.

The SIS model introduced in (6) can be derived from a subpopulation
perspective or through a mean-field approximation of a $2^{n}$ state Markov
chain model (Paré et al., 2020b). As the epidemics dynamic is usually sampled
at discrete time steps, we will consider the discrete version of model (6):

$p^{k+1}_{i}=p^{k}_{i}+(1-p^{k}_{i})\beta_{i}\sum_{j=1}^{n}\omega_{ij}p^{k}_{i}-\gamma_{i}p^{k}_{i}$
(7)

where $k>0$ is the time index. In matrix form, the model can be re-written as:

$p^{k+1}=p^{k}+\left[(1-P^{k})B\Omega-\Gamma\right]p^{k}$ (8)

where $P^{k}=\mathop{\rm diag}\nolimits{p^{k}}$, $B=\mathop{\rm
diag}\nolimits{\beta_{i}}$, $\Omega=(\omega_{ij})$ and $\Gamma=\mathop{\rm
diag}\nolimits{\gamma_{i}}$. We denote by

$\mathcal{O}_{k}=1+(1-P^{k})B\Omega-\Gamma.$

For this model to be well defined, we need the following assumption: For all
$i\in[n]$, we have $\gamma_{i}<1$ and $\sum_{j=1}^{n}\beta_{i}\omega_{ij}<1$.
In this case, we have the following result:

######  Lemma 3.1 ((Paré et al., 2020b)).

Consider the model (8) under Assumption 3. Suppose $p^{0}_{i}\in[0,1]$ for all
$i\in[n]$. Then, for all $k>0$ and $i\in[n]$, $p^{k}_{i}\in[0,1]$.

Using data on healthy / infected individuals at time point $k^{0}$, we can
estimate ${p}^{k_{0}}$ using the TV estimator given by (2). Let

$\widehat{\mathcal{O}}_{k}=1+(1-\widehat{P}^{k})B\widehat{\Omega}-\Gamma$ (9)

be the evolution operator, where $\widehat{P}^{k}=\mathop{\rm
diag}\nolimits(\widehat{p}^{k})$ and $\widehat{\Omega}$ is an estimator of
$\Omega$ from noisy observations of the network (e.g. using variational
maximum likelihood from some (partial) observations of the network (Gaucher &
Klopp, 2021), as discussed in Remark 3 below). Then, by recursively computing
$\hat{p}^{k}=\widehat{\mathcal{O}}_{k-1}\hat{p}^{k-1}$ for $k>k^{0}$, we can
predict the evolution of the epidemic. In particular, the expected number of
infected individuals at time $k$, $E(I^{k})=\sum p^{k}_{i}=\|p^{k}\|_{l_{1}}$
can be estimated using Proposition 2.5.

######  Proposition 3.2.

Consider the discrete time SIS given by (7). Assume that we observe
$Y^{k_{0}}\in\\{0,1\\}^{n}$ which encodes the data on the healthy / infected
status of a population of $n$ connected individuals at time point $k^{0}$. Fix
$\delta\in[0,1]$, $T\subset[m]$ and assume that the contact graph $G$ is
connected. Then, with probability at least $1-\delta$, we have

$\begin{split}&\left|\|\widehat{p}^{k}\|_{l_{1}}-E(I^{k})\right|\leq\left(1-\gamma_{\min}+\beta_{\max}\|\Omega\|_{1,\infty}\right)^{k-k_{0}}\left(\dfrac{4\rho\sqrt{s|T}|}{\kappa_{T}}\log\left(\frac{4n^{2}}{\delta}\right)+\right.\\\
&\left.+2s\sqrt{\frac{\log\left(\frac{4}{\delta}\right)}{n}}+\textasciicaron
3\sqrt{\rho\,s\left\|\left(Dp^{k_{0}}\right)_{T^{c}}\right\|_{1}\log\left(\frac{4n^{2}}{\delta}\right)}+\sqrt{2}\rho
sd_{\max}(S)\log\left(\frac{4n^{2}}{\delta}\right)\right)\enspace.\end{split}$

where $\gamma_{\min}=\underset{i}{\min}\;\gamma_{i}$ and
$\beta_{\max}=\underset{i}{\max}\;\beta_{i}$

######  Proof 3.3.

We have $E(I^{k})=\|p^{k}\|_{l_{1}}$ and

$\left|\|\widehat{p}^{k}\|_{l_{1}}-\|p^{k}\|_{l_{1}}\right|\leq\|\widehat{p}^{k}-p^{k}\|_{l_{1}}\leq\left(1-\gamma_{\min}+\beta_{\max}\|\Omega\|_{1,\infty}\right)^{k-k_{0}}\|\widehat{p}^{k_{0}}-p^{k_{0}}\|_{l_{1}}.$

The result follows by bounding $\|\widehat{p}^{k_{0}}-p^{k_{0}}\|_{l_{1}}$
using the results of Proposition 2.5.

Note that Assumption 3 implies $\beta\|\Omega\|_{1,\infty}\leq 1$.

######  Remark 3.4.

Beyond SIS model, one can also consider Susceptible - Infected - Recovered
(SIR) and Susceptible - Exposed - Infected - Recovered (SEIR) models (Paré et
al., 2020a). For a SIR process, for example, let $r^{k}_{i}$ denotes the
probability of node $i$ being recovered at time $k$. Then, the following
equations govern the evolution of $p^{k}$ and $r^{k}$ :

$\begin{split}p^{k+1}_{i}&=p^{k}_{i}+(1-p^{k}_{i}-r^{k}_{i})\beta_{i}\sum_{j=1}^{n}\omega_{ij}p^{k}_{i}-\gamma_{i}p^{k}_{i}\\\
r^{k+1}_{i}&=r^{k}_{i}-\gamma_{i}p^{k}_{i}.\end{split}$

Our denoiser can therefore easily be deployed in the context of SIR and SEIR
models by simply substituting the evolution operator of Equation 9 with a
appropriate version.

######  Remark 3.5 (Taking into account false positives).

In the case of epidemic models on networks, it is also necessary to account
for errors that may occur in our observations $Y^{k_{0}}$ due to false
positives. Such false positives might be due to errors due to imperfect test
specificity for instance: while gold standard medical tests can have
specificity up to 99%, reporting based on symptom data is usually less
accurate, and can thus lead to false positives. Here we assume that we know
the false positive rate $\alpha$, and that our observations $y^{k}_{i}$ follow
Bernoulli distribution with probability
$\rho^{k}_{i}=(1-\alpha)p^{k}_{i}+\alpha$. In this case, to estimate $p^{k}$
we add a post-processing step thresholding to zero all the coordinates of
$\widehat{\rho}^{k}$ obtained from (2) that are smaller or equal to $\alpha$.
We test the performance of this method on graphs with various topologies,
using the experimental setup described in Section 5. As shown in Figure 1 for
the case of the k-NN graph, our procedure can significantly improve the
accuracy of the nowcasted estimates $\hat{p}$ over the naive estimates
$\hat{p}^{\text{naive}}=y_{\text{observed}}$, particularly as the epidemic
size increases. Plots for other graph topologies (small-world, preferential
attachment, etc) can be found in Appendix D.

Figure 1: $\ell_{1}$ error for the k-NN graph, as a function of the noise
level (false positive rate) $\alpha\in\\{0,1\%,5\%\\}$, $k_{0}$ and $\beta$.
The healing rate is here fixed at $\gamma=0.1.$. The benchmark values (in red)
have been slightly shifted on the x-axis to improve legibility. The results
are averaged over 100 independent experiments. Error bars indicate
interquartile ranges.

######  Remark 3.6 (Estimating the strength of the connections in the
network).

One of the key quantities that determines the evolution of the epidemic on the
network in (7) is the matrix $\Omega=(\omega_{ij})_{(ij)\in[n]\times[n]}$,
which quantifies the strength of the connections between nodes in the network.
We can estimate $\Omega$ from the observation of a network using the
inhomogeneous random graph model assuming that for $i<j$, the elements
$A_{ij}$ of the adjacency matrix $A$ are independent Bernoulli random
variables with the success probability $\omega_{ij}\in[0,1]$. In this case,
each entry $\omega_{ij}$ can be interpreted as the probability that the edge
$(i,j)$ is present in the graph $\mathcal{G}$.

The problem of estimating $\Omega$ from a single observation of the network
has been considered in several papers (Gao et al., 2015; Klopp et al., 2017;
Klopp & Verzelen, 2019; Gaucher & Klopp, 2021) where strong theoretical
guarantees have been obtained for this problem. One of the most popular
methods is based on the variational approximation of the likelihood function
(Celisse et al., 2012; Bickel et al., 2013). Loosely speaking, this method
attaches to the latent variables a distribution with free parameters. These
parameters are fitted in a way to obtain a distribution close in Kullback-
Leibler divergence to the true posterior. The crux of this approach lies in
the realisation that the variational distribution is simpler than the true
posterior, which allows to solve the optimization problem approximately, see,
for example, Gaucher & Klopp (2021).

######  Remark 3.7 (Learning parameters $\gamma$ and $\beta$.).

Assume that we know $p_{k},p_{k+1},\dots p_{K}$. We can write the following
equation:

$\begin{bmatrix}p^{k+1}-p^{k}\\\ p^{k+2}-p^{k+1}\\\ \dots\\\
p^{K-1}-p^{K}\end{bmatrix}=\Phi\begin{bmatrix}\beta\\\ \gamma\end{bmatrix}$

where

$\Phi=\begin{bmatrix}(I-P^{k})\Omega p^{k}&-p^{k}\\\ \dots\\\
(I-P^{K-1})\Omega p^{K-1}&-p^{K-1}\\\ \end{bmatrix}$

Vrabac et al. (2020) propose to estimate the vector of parameters
$[\beta,\gamma]^{T}$ using the inverse (or pseudoinverse) of matrix $\Phi$.
However, in practice, we do not have access to $p^{j}$ but only to its noisy
counterpart $Y_{i}$, as given by (1). In such situations, one can replace
$p^{j}$ by its estimator $\widehat{p}^{j}$ given by (2) or (11) and compute
the inverse (or pseudoinverse) of matrix $\Phi$ using $\widehat{p}^{j}$. We
illustrate the efficiency of this approach in numerical studies in Section
5.3.

##  4 Partially observed epidemics

In (2), we assume that we observe the state $y_{i}\in\\{0,1\\}$ for any node
in the network. However, in many real-world situations, it is challenging to
obtain complete information about each individual’s disease status. People may
not be aware of their infection or may choose not to disclose it. Therefore,
methods that can work with partial information become indispensable. In this
section, we extend our analysis to the case of partial observations of
$Y^{k}$.

Let $\mathfrak{M}=(m_{i})$ be a mask vector. That is, $m_{i}=1$ if we observe
the information about node $i$, and $m_{i}=0$ if not. We assume that $m_{i}$
are independent Bernoulli random variables with parameter $\pi_{i}$:
$m_{i}\sim\mathcal{B}(\pi_{i})$ and $\pi_{i}>0$. In the simplest setting, each
coefficient is observed with the same probability, i.e., for every $i\in[n]$,
$\pi_{i}=\nu$. In the case of partial observations, we observe two vectors,
the mask vector $\mathfrak{M}$ and the vector $\widetilde{Y}=\mathfrak{M}\circ
Y$ where $\circ$ denotes the entry-wise (Hadamar) product. For any
$v=(v_{i})_{i\in[n]}\in\mathbb{R}^{n}$ we define the weighted $l_{2}$-norm of
$v$ as:

$\|v\|_{l_{2}(\pi)}^{2}=\sum_{i}\pi_{i}v^{2}_{i}.$

Let

$\kappa_{\pi}=\left(\underset{i}{\min}\;\pi_{i}\right)^{-1}$

and $\pi^{-1}=\left(\pi_{i}^{-1}\right)_{i\in[n]}$. It is easy to see that for
any $v\in[0,1]^{n}$

$\|v\|_{2}\leq\kappa_{\pi}\|v\|_{l_{2}(\pi)}.$ (10)

We use the partial observations $(\mathfrak{M},\widetilde{Y})$ to construct
the TV denoiser $\widehat{p}_{miss}$ associated to $G$ as any solution of the
following minimization problem:

$\widehat{p}_{miss}\in\underset{p\in[0,1]^{n}}{\mathop{\rm
arg\,min}}\left\\{\dfrac{1}{n}\overset{n}{\underset{i=1}{\sum}}m_{i}\left(\widetilde{y}_{i}-p_{i}\right)^{2}+\lambda\|Dp\|_{1}\right\\}.$
(11)

We prove the following bound on the risk of the TV denoiser estimated on
partially observed data:

######  Theorem 4.1 (Risk bound for TV denoising with partial observations).

Assume that $G$ is connected and define the regularization parameter

$\lambda=\frac{9\sqrt{2}\rho\,\log\left(n\right)}{n}.$ (12)

For $T={\rm supp}\left(Dp^{*}\right)$, the TV denoiser from partial
observations, $\widehat{p}_{miss}$, defined in (11) satisfies

$\begin{split}{\|\widehat{p}_{miss}-p^{*}\|^{2}_{l_{2}(\pi)}}\leq
C\left\\{\dfrac{\rho^{2}\,\kappa_{\pi}|T|}{\kappa^{2}_{T}}+\dfrac{\|\pi\|_{1}\|\pi^{-1}\|_{1}}{n^{2}}\right\\}\log^{2}(n)\end{split}$
(13)

with probability at least $1-c/n$ where $C,c$ are fixed numerical constants.

The proof of Theorem 4.1 is provided in Appendix C. These results therefore
allow us to provide theoretical guarantees on the recovery of $p^{*}$ despite
missing observations. In particular, we consider the setting where
$\pi_{i}=v$, in which case, $v$ captures the fraction of observed nodes within
the network. In this case, the dependency between the number of missing
observations and the quality of the bound becomes clearer: as expected, the
smaller the fraction $v$, the larger the error.

##  5 Numerical experiments

In this section, we validate our method through a series of synthetic and
semi-synthetic experiments. These controlled experiments allow us to test a
variety of data regimes by varying the topology of the contact graph - a
dependency highlighted by the presence of the scaling factor $\rho$ in
Proposition 3.2 \- as well as varying the parameters of the epidemic model
($\beta$ and $\gamma$).

To simulate an epidemic process, we first sample a graph on $n$ nodes, where
$n$ is fixed to 1,000 unless otherwise specified. In the results presented
here, we consider a variety of random graph models, including Erdös-Rényi
graphs (where each edge is independently sampled with probability $p$),
Stochastic Block Models, Small World Networks (parametrised by their average
degree $m$ and a rewiring parameter $p$), Power-Law graphs (parametrised by a
parameter, $m$, indicating the number of edges to attach from a new node to
existing nodes as the graph is generated), as well as $k$-Nearest Neighbour
graphs (generated by uniformly sampling 2D coordinates between 0 and 1, and
connecting nodes to their $k$-nearest neighbours in Euclidean norm). To
further improve the realism of our experiments and consider graphs reflective
of real-world contact network characteristics, we also generate an epidemic
process on the Berkeley graph (Rossi & Ahmed, 2015). The latter is a large
social friendship network extracted from Facebook consisting of users (nodes)
with edges representing friendship ties. This also allows us to test the
performance of our algorithm on a larger-size graph, as this social network
has indeed a total of 22,900 nodes and 852,419 edges.

For each experiment, we randomly sample a “patient 0”, and propagate the
epidemic on the graph $\mathcal{G}$ as per Equation (8) for $k_{0}$ steps,
using a specified value of the infection rate $\beta$ and healing rate
$\gamma$. The probability of infection for each node here is chosen to be
identical: $\beta_{i}=\beta$. We select the interaction strength between each
node $i,j$ as $w_{ij}=\frac{1}{d_{i}\vee d_{j}}$. This models a situation in
which high degree nodes (e.g. influencers) have looser connections to their
neighbours, while peripheral nodes with fewer connections might have tighter
bounds to their neighbours. The rationale for this choice of prior is
reinforced by the observation that the rate of epidemic transmission increases
over time. Individuals with a high number of connections, denoted by high
degrees, are here presumed to have, in average, brief interactions with each
contact. In contrast, individuals with fewer connections tend to represent
closer-knit groups, such as family units, where extended interactions are more
common. Note that Assumption 3 imposes additional conditions on the topology
of the graph. For instance, in Erdös-Rényi (ER) graphs, Assumption 3 implies
that $np_{ER}\beta\leq 1$ so that $p_{ER}=O(\frac{1}{n})$ — in other words,
the ER graph cannot be too dense. This is by no means a limitation of our
model, but rather, a way of establishing graph topologies in which our
sparsity assumption can be verified.

In our experiments, unless specified otherwise, $\gamma$ is fixed to 0.1. This
corresponds to an average recovery time of 9 days, assuming that the recovery
is distributed as a negative binomial with a probability $\gamma$ of recovery
at each time step. The result of the simulation is an underlying probability
$(p^{\star})_{i\in[n]}$ that each node is infected after $k_{0}$ epidemic
diffusion steps. In this setting, the lower the value of $k_{0}$, the more
localized the epidemic. Finally, we generate a vector of observations,
sampling each entry (i.e. the infection status of each node) as
$y_{i}\sim\text{Bernouilli}(p^{\star}_{i})$. An example of this diffusion
process is presented in Figure 2.

(a) k-Nearest-Neighbour Graph: Time 0

(b) k-Nearest-Neighbour Graph: Time 10

(c) k-Nearest-Neighbour Graph: Time 20

Figure 2: Visualisation of a SIS epidemic process on a 5-nearest neighbour
graph at time 0 (left: one initial infected node), and after $k_{0}=10$
(middle) and $k_{0}=20$ (right) diffusion steps, for $\beta=0.9$.

###  5.1 Algorithm.

To efficiently solve this problem, we use the semi smooth-Newton augmented
Lagrangian method (SSNAL) of Sun et al (Sun et al., 2021). SSNAL is a scalable
and efficient algorithm that was originally designed to solve the following
optimization problem:

$\min_{X\in\mathbb{R}^{d\times
n}}\frac{1}{2}\sum_{i=1}^{n}\|x_{i}-a_{i}\|^{2}+\lambda\sum_{i<j}W_{ij}\|x_{i}-x_{j}\|_{2}$
(14)

where $W_{ij}\in[0,1]$ is a matrix of weights such that $W_{ij}=0$ if nodes
$i$ and $j$ are disconnected in the graph, and reflects the proximity between
node $i$ and $j$ otherwise. This formulation appears in several settings,
particularly in convex clustering problems, where it is used to define a
hierarchical clustering of the data. Note that in our case, since our vector
$(p_{i})_{i\in[n]}$ is of dimension 1, so
$\|p_{i}-p_{j}\|_{2}=\sqrt{(p_{i}-p_{j})^{2}}=|p_{i}-p_{j}|$, and our problem
shares the same objective function as Equation (14). In SSNAL, the
optimisation problem (14) is solved using an augmented Lagrangian method on
equation (14), and solving the corresponding sub-problems using semi-smooth
Newton conjugate gradient updates. In their work, Sun et al. (2021) show that
as long as the dimension of the feature vector is small (i.e., in our case,
equal to 1), this method has the potential to considerably improve existing
solvers. We provide a comparison of this solver against alternative solvers
(e.g. ADMM (Boyd et al., 2011) and CVX (Diamond & Boyd, 2016)) in Appendix D.
In our experiments, an appropriate value for the parameter $\lambda$ is chosen
by cross-validation111The code for all of our simulations and analyses is
available at the following link: https://github.com/donnate/epidemics..

###  5.2 Results.

Figures 3 and 4 show the results for the k-NN graph, as the values of $k$,
$k_{0}$ and $\beta$ vary. In these plots, the $\ell_{1}$ and $\ell_{2}$ errors
(Figure 3 and 4, respectively) of our estimator are presented as a function of
the probability of transmission $\beta$. Each column represents a number of
diffusion of steps ($k_{0}$): as emphasized earlier, the higher the value, the
more advanced — and therefore difficult to denoise, from an inference
perspective — the epidemic spread. Each row corresponds to a different graph
topology: in this case, the higher the value of $k$, the denser the contact
network. Table 1 presents the results for the Berkeley graph, and plots
corresponding to other topologies are displayed in Figures 11 and 12, Appendix
D. Overall, we observe that for graphs where the epidemic has not propagated,
the denoising method exhibits only marginal improvement over the naive
estimator $\hat{p}^{\text{naive}}=y_{\text{observed}}$: in the k-NN and Power-
Law graphs for instance, when the number of steps is small, the average number
of infections is low (less than 10— see Figure 10), and the one-bit TV
denoising method only improve the raw observations by 7%
($\|\hat{p}-p^{*}\|_{1}=2.03$ for the 2-NN denoised estimate vs 2.17 for the
raw observations at $k_{0}=10$). However, as the size of the epidemic
increases, the effect of the denoising estimator grows considerably. For the
5-NN graph, with $\beta=0.5$ and $k_{0}=30$, our denoising leads to an almost
50% decrease of the $\ell_{1}$ error for the estimated current state of the
epidemic.

Figure 3: Effect of the denoising (with $\lambda$ chosen through cross
validation), as a function of epidemic time (columns) and graph topology
(row), measured by the difference in $\ell_{1}$-norm between the ground-truth.
Points denote the mean $\ell_{1}=\|p^{*}-\hat{p}\|_{1}$ error over 200
simulations for $k=2$, and 100 simulations for the other topologies. Error
bars indicate interquartile ranges. A plot showing the size of the
corresponding epidemic can be found in figure 10, Appendix D. Figure 4: Effect
of the denoising (with $\lambda$ chosen through cross validation), as a
function of epidemic time (columns) and graph topology (row), measured by the
difference in square $\ell_{2}$-norm between the ground truth and estimated
infection parameter $p$. Points denote the mean
$\ell_{2}=\|p^{*}-\hat{p}\|_{2}^{2}$ error over 200 simulations for $k=2$, and
100 simulations for the other topologies.

Figures 5 shows the $\ell_{2}$ error in the predicted state of the epidemic in
2 days. These estimates are obtained by entering the estimated value of
$\hat{p}$ in Equation (7), and rolling out the updates for two time steps.
This yields potentially substantial differences in the predictions of the
state of the epidemic, as depicted in Figure 5.

Steps | $\beta$ | $\|\hat{p}-p^{*}\|_{1}$ | $\|\hat{p}^{\text{naive}}-p^{*}\|_{1}$
---|---|---|---
5 | 0.7 | 0.00014 | 0.00014
5 | 0.9 | 0.00024 | 0.00026
10 | 0.3 | 0.00011 | 0.00012
10 | 0.5 | 0.00026 | 0.00032
10 | 0.7 | 0.00094 | 0.00124
20 | 0.3 | 0.00033 | 0.00042
Table 1: Accuracy of the estimation (in $\ell_{1}$ norm) for the Berkeley
Graph. Note that in both cases, the accuracy is high, since the epidemic in
this case remains very localised. Figure 5: $\ell_{2}$ error on the state of
the epidemic in 2 days, using the predicted $\hat{p}^{(2)}$. Results are
compared to the naive estimator that solely uses the observed values:
$\hat{p}^{\text{naive}}=y_{\text{observed}}$ to nowcast the epidemic and is
then propagated twice as per Equation (9). The benchmark values (in red) have
been slightly shifted on the x-axis to improve legibility. The results are
averaged over 100 independent experiments. Error bars indicate interquartile
ranges.

###  5.3 Estimating parameters $\gamma$ and $\beta$.

To illustrate the efficiency of 1-bit denoiser in the case of estimation of
parameters $\gamma$ and $\beta$ we use a kNN-nearest neighbour graph and
propagate the epidemic for 20 steps. We vary the values of $\beta$ and
$\gamma$, using $k_{0}=30$. We compare the values of the estimated values of
$\beta$ and $\gamma$ as reported by our estimator with those using the same
method, but fitting instead the observed values of the state of the epidemic
(the naive estimator) . The results for the estimates of $\beta$, $\gamma$ and
the reproductive number $R_{0}=\frac{\beta}{\gamma}$ \- an important
characteristic of the disease - , averaged over 100 simulations, are reported
in Table 2. In most cases, our TV denoiser improves the parameter estimates,
achieving almost perfect reconstruction for $\beta=0.8$ and $m=5$. In Figure
6, we evaluate the effect of the number of observed time steps $k_{0}$, the
graph parameter $k$, and the transmission probability $\beta$ on the
estimation error. Overall, we observe that, for this method to be accurate,
the epidemic must have a size that is neither too big nor too small: if the
epidemic is too big (e.g. as it tends to be the case when $\beta=0.9)$, the
sparsity assumption ceases to hold, while if the epidemic is too small (e.g.
as it tends to be the case when $\beta=0.2)$, there are too few observations
to reliably estimate the parameters.

k-NN | $\beta$ | $\gamma$ | $R_{0}$ | Method | $\hat{\beta}$ | $\hat{\gamma}$ | $\hat{R}_{0}$
---|---|---|---|---|---|---|---
|  |  |  | TV denoiser | 0.78 (0.58, 1) | 0.32 (0.23, 0.44) | 2.4 (2.27, 2.55)
| 0.35 | 0.1 | 3.5 | Naive | 0.99 (0.72, 1) | 0.39 (0.27, 0.52) | 2.47 (1.91, 2.7)
|  |  |  | TV denoiser | 0.72 (0.58, 0.84) | 0.17 (0.13, 0.22) | 4.1 (3.89, 4.38)
| 0.5 | 0.1 | 5 | Naive | 0.85 (0.69, 1) | 0.21 (0.15, 0.27) | 4.11 (3.76, 4.39)
|  |  |  | TV denoiser | 0.84 (0.72, 0.98) | 0.11 (0.09, 0.14) | 7.53 (7.07, 8.2)
5 | 0.8 | 0.1 | 8 | Naive | 0.94 (0.8, 1) | 0.12 (0.09, 0.14) | 7.38 (6.21, 8.02)
|  |  |  | TV denoiser | 0.62 (0.48, 0.78) | 0.23 (0.17, 0.31) | 2.64 (2.51, 2.81)
| 0.35 | 0.1 | 3.5 | Naive | 0.78 (0.6, 0.97) | 0.29 (0.21, 0.38) | 2.67 (2.53, 2.81)
|  |  |  | TV denoiser | 0.66 (0.55, 0.76) | 0.15 (0.12, 0.19) | 4.29 (4.09, 4.6)
| 0.5 | 0.1 | 5 | Naive | 0.77 (0.64, 0.9) | 0.18 (0.14, 0.22) | 4.26 (4.08, 4.53))
10 |  |  |  | TV denoiser | 0.83 (0.74, 0.93) | 0.11 (0.09, 0.13) | 7.75 (7.38, 8.11)
| 0.8 | 0.1 | 8 | Naive | 0.91 (0.82, 1) | 0.12 (0.1, 0.14) | 7.62 (7.11, 7.97)

Table 2: Results for the estimation of the parameters $\gamma$ and $\beta$, as
well as the reproductive number $R_{0}=\frac{\beta}{\gamma}.$ Intervals
denotes interquartile ranges. Results averaged over 500 simulations, for
$K_{0}=30$. Figure 6: Error in the estimation of $\beta$, as a function of
graph density for the kNN graph (parametrized by $k$, represented by each
column) and number of observation steps $K$. Results are aggregated over 500
experiments. Points denote the mean error in estimating $\beta$, and error
bars indicate the interquartile range.

###  5.4 Partially observed epidemic

We also study the performance of the 1-bit TV denoiser in the case of partial
observations, that is the goodness-of-fit term of our estimator is only
applied to nodes whose infectious statuses are observed. We vary the
proportion of missing values, and plot the corresponding estimation error. The
results are presented in Figure 7. We observe that our estimator can provide
substantial improvement over the naive estimator (where all masked entries are
filled in with 0) for even large proportions of missing values. This
improvement is particularly substantial in cases where the epidemic has
progressed enough so that the trivial estimation $p^{*}_{1}=0$ is no longer
accurate, but not enough to violate our sparsity assumption on the size of the
support $\|p^{*}\|_{0}.$

Figure 7: $\ell_{1}$ error of our one-bit TV denoiser compared to
$\hat{p}^{\text{naive}}=y_{\text{observed}}$ in the presence of missing data
as a function of the proportion of missing values for the 5-NN graph. Here we
vary $k_{0}\in\\{10,20,30\\}$ (one value of $k_{0}$ per row), while columns
corresponds to various values of $\beta.$

##  6 Real Data Analysis

In this section, we propose to analyse two real epidemics datasets using our
method.

###  6.1 Analysis of the COVID-19 outbreak in California

To demonstrate the applicability of our method for real-time forecasting
(nowcasting), we propose to re-analyse the initial COVID-19 outbreak in
California. Our analysis considers the graph formed by the 58 California
counties and their corresponding incidence data from January 31, 2020, to
April 1, 2020. Incidence data are known to suffer from considerable variations
throughout the week and underreporting (Graham, 2020; Lipsitch et al., 2015),
which can significantly skew the estimation of the reproductive number and the
predicted epidemic trajectory (Siegenfeld et al., 2020). It is, in fact,
standard practice to first smooth the raw incidence data through a 7-day
rolling average to mitigate irregularities due to weekend reporting (see, for
instance, case reporting by the Johns Hopkins Coronavirus Resource Center
(Johns Hopkins CRC, 2024), or by the World Health Organization (World Health
Organization, 2024)). However, while the importance of accounting for temporal
variations seems widely recognized, the extension of this practice to the
spatial domain remains unexplored. Yet, due to movement of people between
counties, it seems logical to assume that neighbouring counties should expect
similar prevalence: any deviations from this pattern could likely be
attributed to measurement bias, and should therefore be corrected. We propose
here to correct such biases by employing our method, which we believe leads to
a more precise estimation of local disease prevalence.

We determine the optimal regularization parameter $\lambda$ by splitting the
data between training and testing. For a given $\lambda$, we train a simple
epidemic model on the smoothed data from January 31 to March 14, 2020. This
first part of the data is used to estimate the basic reproduction number
($R_{0}$) and to forecast the following two weeks. Fitting and predictions are
performed here using the packages R0 and earlyR, which are designed for the
estimation and prediction of early-stage epidemics. These packages use the
incidence matrix to estimate $R_{0}$ and predict the epidemic trajectory using
simple growth models, which only require as input an estimate of the mean and
standard deviation of the serial interval. Following Du et al. (2020), we set
these values at $\mu=3.95$ and $\sigma=4.75$, respectively.

We then choose the regularization parameter ($\lambda$) that yields the most
accurate forecast on the next three days, from March 15 to March 18, 2020.
Having chosen an appropriate $\lambda$, we evaluated our model in the period
from March 19 to March 31, 2020. The primary reason for this partitioning is
that COVID reporting began only in early March for a subset of California
counties. However, on March 16, 2020, California implemented a lockdown order,
aiming to aggressively reduce $R_{0}$ and therefore changing the dynamics of
the epidemic. However, we anticipate that the incidence immediately after the
lockdown reflects the initial dynamics of the epidemic spread, as such changes
are not instantaneous given the incubation period of the virus of around 4 to
5 days.

Given the lack of individual-level data and having only access county-level
aggregates, we find it easier to model each county as a ’super node’ - or
equivalently, the entire population as a graph such that each inhabitant is
connected to all other inhabitants of the same county, as well as those in
neighbouring counties. Letting $n_{c}$ denote the number of inhabitants in
county $c$, this approach thus requires solving the optimisation problem
defined as follows:

$\mathop{\rm
arg\,min}_{p\in\mathbb{R}^{C}}\sum_{c\in[C]}I_{c}(1-p_{c})^{2}+(n_{c}-I_{c})p_{c}^{2}+\lambda\sum_{c,c^{\prime}}n_{c}n_{c}^{\prime}w_{c,c^{\prime}}|p_{c}-p_{c^{\prime}}|$
(15)

where $w_{c,c^{\prime}}$ is a set of appropriately chosen weights. In the
previous equation, $I_{c}$ denotes the total number of infected inhabitants in
county $c$, while $n_{c}-I_{c}$ represents the number of susceptible
individuals in county $c$. Here, consistently with our simulation experiments,
we pick $w_{c,c^{\prime}}=\frac{1}{\sum_{c^{\prime}\sim c}n_{c^{\prime}}}$ the
inverse of the degree of each node.

The results of this additional preprocessing procedure are presented in Table
3. Overall, for 85% of the counties on which the epidemic could be estimated
(meaning that these counties had at least one COVID case before March 20th),
our estimator provides a significant improvement over the raw data. This
improvement is particularly striking in well populated and connected areas,
such as the San Francisco Bay, where our estimator systematically improves the
naive estimator.

county |  | MSE Prediction
---
(Raw Data)
| MSE Prediction
---
(Smoothed Data)
Contra Costa | 60 | 27
Fresno | 0 | 347
Humboldt | 0 | 0
Los Angeles | 652 | 425
Marin | 285 | 17
Orange | 136 | 73
Placer | 0 | 1
Riverside | 32 | 25
Sacramento | 25 | 18
San Benito | 1,382 | 297
San Diego | 2,140 | 716
San Francisco | 63 | 23
San Joaquin | 2,157 | 351
San Mateo | 158 | 114
Santa Clara | 101 | 78
Solano | 2 | 0
Sonoma | 147 | 3
Stanislaus | 784 | 211
Tulare | 4,558 | 3
Ventura | 8,679 | 3,043
Yolo | 3 | 32
Table 3: Comparison of the results in terms of the MSE between predicted
epidemic trajectory and realized trajectory for the period from March 20th to
March 30th 2023, using the raw data or the spatially smoothed data from our
estimator.

###  6.2 The ExFLU Dataset of Aiello et al. (Aiello et al., 2016)

In this second dataset, we reinvestigate the epidemic dataset of Aiello et al.
(2016). This dataset records a flu epidemic outbreak over a college campus in
Winter 2013. The study’s initial aim was to examine the influence of social
interventions on the transmission of respiratory infections. A total of 590
university students were enrolled in the study, engaging in weekly surveys
regarding symptoms resembling influenza-like illness (ILI) and their social
interactions. This enables the generation of a weekly contact graph, which
helps in mapping potential epidemic propagation phenomena. More specifically
here, to evaluate the potential of our method to denoise contact tracing
networks, we consider the ego network of each participant declaring
experiencing ILI sympoms induced by the contact graph spanning the week of,
and the week after their symptom onset. We define as a contagion event any
symptomatic or asymptomatic individual having contracted an ILI illness based
on the day of or after the contact’s reported symptom onset. This creates a
total of 97 ego networks. We then evaluate the ability of our one-bit
denoising approach in correctly estimating the ILI statuses of the nodes
within this contact network the following week (see Figure 8).

To select the appropriate regularisation parameter, we split the ego networks
into two sets, and select lambda on the first, and evaluate on the second. We
randomize this procedure and average the results over $100$ iterations. As in
the previous example, we quantify the error using the $l_{2}$ norm. We obtain,
in average, an $l_{2}$ error of $0.0482$, compared to $0.0493$ when setting
the regularisation parameter to $0$. Despite the noisy nature of this dataset,
we are therefore able to improve upon the naive estimator $(\lambda=0).$ This
suggests that our TV denoiser might be effective for epidemic nowcasting using
contact tracing networks.

Figure 8: Examples of contact networks induced by taking the 2-hop induced
subgraph around participants exhibiting ILI-like symptoms. The contact network
here represents all reported contacts between participants on the week of and
the week following the symptom onset of the egonet source node. Colours denote
different illness statuses. Grey: healthy participant. Yellow: participant
sick with ILI symptoms on the week prior to the source node’s symptom onset.
Red: participant healthy before the source node’s symptom onset and developing
ILI symptoms in the week following the source node’s symptom onset.

##  7 Related Works

In this section we provide a more in-depth overview of related works.

Total Variation Denoising. Total Variation (TV) denoising (also known as Fused
Lasso) has experienced a considerable practical success in the area of image
denoising. First introduced in Rudin et al. (1992), it benefits from the
strong theoretical guarantee in the case of Gaussian noise. We refer the
readers to the works of Mammen & Geer (1997); Dalalyan et al. (2017) for chain
graphs and Wang et al. (2016); Hütter & Rigollet (2016); Padilla et al. (2018)
for more general graphs. Earlier works on TV denoising (Mammen & Geer, 1997)
marked an important step by obtaining the first statistical guarantees.
Another notable development in understanding the behaviour of the TV denoiser
was presented in Needell & Ward (2013). In their work, the authors focused on
cases where the noise possesses a small $l_{2}$ norm but is otherwise
arbitrary. This framework is commonly encountered in the literature on noisy
compressed sensing. An analysis of the statistical performance of total
variation image denoising was performed in Wang et al. (2016) for sub-Gaussian
noise and for Gaussian white noise model in Hütter & Rigollet (2016) where, in
particular, optimal fast rates for the 2D grid are obtained. By contrast, in
this work, we built on the rich literature on total variation denoising and
provide an analysis of its statistical guarantees in the case of binary
observations (see Section 2).

Epidemic spreads on networks. Works on epidemic spread on networks and, more
generally, on graph dynamical systems, can be found in various areas. Common
types of network epidemic models include SIR and SEIR models, Agent-Based
models, Percolation models and etc. Network-based stochastic models, such as
the Susceptible-Infectious-Recovered (SIR) and Susceptible-Exposed-Infectious-
Recovered (SEIR) models, are commonly used to simulate epidemic dynamics.
Transition rates between compartments are estimated based on network structure
and contact patterns (Newman, 2002b; Kenah & Robins, 2007). Researchers have
extended these models to incorporate heterogeneous mixing patterns and complex
network structures (Miller & Volz, 2013). Agent-based models (ABMs) have
gained popularity for simulating realistic human behavior within network
contexts. For example, Balcan et al. (2009) used ABMs to investigate the
impact of human mobility on epidemic spread, considering factors like travel
patterns and contact behaviour. These models provide insights into the
effectiveness of public health interventions and behavioral changes during
outbreaks.

Stochastic network models, such as the temporal exponential random graph model
(TERGM), capture dynamic changes in network structure over time. They are
employed to study how changes in contact patterns impact epidemic dynamics
(see, for example, Volz & Meyers (2007)). These models provide a statistical
framework for exploring network evolution during epidemics.

Discrete-time Markov Processes (DTMP) and Continuous-time Markov Processes
(CTMP) are commonly used in network epidemic models to represent the dynamics
of disease spread or information diffusion within a population (Gómez et al.,
2010; Zhang & Moura, 2014). In DTMP, time is divided into discrete intervals
or time steps. For example, in a SIR (Susceptible-Infectious-Recovered) model,
susceptible individuals may become infected with a certain probability during
each time step. In Continuous-time Markov Processes (CTMP), time is treated as
a continuous variable, allowing for modeling events that occur at any moment
in time. State transitions can occur at any time, not just at discrete time
steps. Events such as contacts, infections, recoveries, or interventions are
represented as stochastic processes. The rates at which these events occur
depend on factors such as contact frequencies, transmission probabilities, and
recovery rates.

Bayesian inference methods have gained popularity for parameter estimation in
epidemic models. Markov Chain Monte Carlo (MCMC) algorithms, such as the
Metropolis-Hastings algorithm, are used to sample from posterior distributions
of model parameters (Jewell et al., 2009; Hu et al., 2017). Bayesian
frameworks allow, in particular for uncertainty quantification and model
selection. However, to the best of our knowledge, no statistical guarantees
have yet been provided for the recovery of the parameters.

Sampling and contact tracing data are essential for estimating epidemic
parameters. Lloyd-Smith et al. (2005) used network-based contact tracing to
reconstruct transmission chains in emerging infectious diseases. This approach
provides insights into superspreading events and helps identify high-risk
nodes. Statistical techniques for inferring contact networks from observed
epidemic data have been also developed. For instance, Rocha et al. (2011)
introduced a likelihood-based method for reconstructing contact networks using
observed infection times. Such methods are crucial for understanding the
underlying network structure in real-world epidemics.

The study of infectious disease spread in the context of partially observed
epidemics is essential in real-world situations where complete information
about individuals’ health status is often unavailable or uncertain. Research
in this area has addressed various aspects of the problem, including modeling,
inference, and control strategies (for a review of the latter, see for
instance Britton & Giardina (2016)). Several models have been developed that
incorporate partial information. For example, the Susceptible-Infectious-
Unknown (SIU) model extends the classic Susceptible-Infectious-Recovered (SIR)
model to account for individuals with unknown health status.

Within the broader context of partially-observed Markov processes, several
techniques have been incorporated into the R package POMP (King et al., 2016).
Notably, one of these methods is maximum iterated filtering (MIF) (Ionides et
al., 2011). Other approaches include Gaussian approximation of the epidemic
density-dependent Markovian jump process (Narci et al., 2021). Additionally,
likelihood-free approaches such as approximate Bayesian computation based on
sequential Monte Carlo (Toni et al., 2009; Sisson et al., 2007) and particle
Markov chain Monte Carlo (Andrieu et al., 2010) have emerged. These approaches
aim to infer the hidden states of individuals (infected or healthy) based on
observed data, such as reported cases or symptom data.

## Acknowledgments

Claire Donnat gratefully acknowledges support from NSF Award Number
RI:2238616, as well as the resources provided by the University of Chicago’s
Research Computing Center. The work of Olga Klopp was funded by CY Initiative
(grant “Investissements d’Avenir” ANR-16-IDEX-0008) and Labex MME-DII
(ANR11-LBX-0023-01). This work was done while O. Klopp and C. Donnat were
visiting the Simons Institute for the Theory of Computing. The work of N.
Verzelen has been partially supported by ANR-21-CE23-0035 (ASCAI).

## Supplementary material

The Supplementary Material includes the proofs of the theorems and proposition
established in this paper, as well as further results stemming from our
synthetic experiments.

##  Appendix A Proof of Theorem 2.4

The proof of Theorem 2.4 is close in spirit to the proof of the Theorem 2 in
(Hütter & Rigollet, 2016) the main difference being the control of the noise
term.

Applying the chain rule, we have that the subdifferential of the $l_{1}$ term
is given by

$\partial\|Dp\|_{1}=D^{T}\mathop{\rm sign}(Dp)$

where

$\mathop{\rm sign}(x)_{i}=\begin{cases}1&\text{if }x_{i}>0,\\\ [-1,1]&\text{if
}x_{i}=0,\\\ -1&\text{if}x_{i}<0.\\\ \end{cases}$

Using first order optimality condition for the convex problem (2), we derive
that there exists $z\in\mathop{\rm sign}(D\widehat{p})$ such that for any
$\bar{p}\in\mathbb{R}^{n}$, we have

$\dfrac{2}{n}\langle\bar{p},Y-\widehat{p}\rangle=\lambda\langle\bar{p},D^{T}z\rangle=\lambda\langle
D\bar{p},z\rangle\enspace.$ (16)

This implies that

$\dfrac{2}{n}\langle\widehat{p},Y-\widehat{p}\rangle=\lambda\|D\widehat{p}\|_{1}\quad\text{and}$
$\dfrac{2}{n}\langle p^{*},Y-\widehat{p}\rangle\leq\lambda\|D{p^{*}}\|_{1}.$

Subtracting the two terms above, we get

$\dfrac{2}{n}\|p^{*}-\widehat{p}\|^{2}_{2}\leq\dfrac{2}{n}\langle\xi,\widehat{p}-p^{*}\rangle+\lambda\|Dp^{*}\|_{1}-\lambda\|D\widehat{p}\|_{1}.$

Next, we control the error term. Let $\Pi$ denote the projection on $\ker D$,
the kernel of $D$. Then, by the definition of $D^{\dagger}$,
$D^{\dagger}D=I-\Pi$ is the projection on $(\ker D)^{\perp}$. We have that
$\ker D=\ker(D^{T}D)=\ker L$. We assume that $G$ is connected, which implies
that zero is an eigenvalue of $L$ of multiplicity one and $\mathop{\rm
dim}\nolimits\ker D=1$. We can write

$\displaystyle\langle\xi,\widehat{p}-p^{*}\rangle$
$\displaystyle=\langle\Pi\xi,\widehat{p}-p^{*}\rangle+\langle(I-\Pi)\xi,\widehat{p}-p^{*}\rangle$
$\displaystyle=\langle\Pi\xi,\Pi(\widehat{p}-p^{*})\rangle+\langle
D^{\dagger}D\,\xi,\widehat{p}-p^{*}\rangle$
$\displaystyle=\langle\Pi\xi,\Pi(\widehat{p}-p^{*})\rangle+\langle\xi,D^{\dagger}D\left(\widehat{p}-p^{*}\right)\rangle$
$\displaystyle=\langle\Pi\xi,\Pi(\widehat{p}-p^{*})\rangle+\langle\left(D^{\dagger}\right)^{T}\xi,D\left(\widehat{p}-p^{*}\right)\rangle\enspace.$
(17)

To bound the right hand side in (A), we first use the Hölder inequality to get

$\displaystyle\langle\Pi\xi,\Pi(\widehat{p}-p^{*})\rangle\leq\|\Pi\xi\|_{2}\|\Pi(\widehat{p}-p^{*})\|_{2}\enspace.$

Note that $L\,\mathbf{1}_{n}=\mathbf{0}_{n}$ and, for any vector
$v\in\mathbb{R}^{n}$, $\Pi v=\frac{1}{n}\langle
v,\mathbf{1}_{n}\rangle\mathbf{1}_{n}$ which implies

$\|\Pi\xi\|_{2}=\tfrac{1}{\sqrt{n}}|\langle\xi,\mathbf{1}_{n}\rangle|.$

Now, using Hoeffding’s inequality we get that, with probability at least
$1-\delta/2$,

$|\langle\xi,\mathbf{1}_{n}\rangle|=\left|\sum_{i=1}^{n}\xi_{i}\right|\leq\sqrt{2\|p^{*}\|_{0}\log\left(\tfrac{4}{\delta}\right)}\,\quad\text{and}\quad\|\Pi\xi\|_{2}\leq\sqrt{\frac{2\|p^{*}\|_{0}\,\log\left(\tfrac{4}{\delta}\right)}{n}}\enspace.$
(18)

On the other hand, we get

$\|\Pi(\widehat{p}-p^{*})\|_{2}=\tfrac{1}{\sqrt{n}}|\langle\widehat{p}-p^{*},\mathbf{1}_{n}\rangle|\leq\tfrac{1}{\sqrt{n}}\|\widehat{p}-p^{*}\|_{1}$
(19)

and (18–19) imply

$\langle\Pi\xi,\Pi(\widehat{p}-p^{*})\rangle\leq\frac{\sqrt{2\|p^{*}\|_{0}\,\log\left(\tfrac{4}{\delta}\right)}}{n}\|\widehat{p}-p^{*}\|_{1}.$
(20)

For the second term in (A) we have that, with probability at least
$1-\delta/2$,

$\displaystyle\langle\left(D^{\dagger}\right)^{T}\xi,D\left(\widehat{p}-p^{*}\right)\rangle$
$\displaystyle\leq\|\left(D^{\dagger}\right)^{T}\xi\|_{\infty}\|D\left(\widehat{p}-p^{*}\right)\|_{1}$
$\displaystyle\leq\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}\enspace,$
(21)

where the last inequality follows from Lemma A.1. Putting together (20) and
(A) we get that, with probability at least $1-\delta$,

$\displaystyle\langle\xi,\widehat{p}-p^{*}\rangle\leq\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}+\frac{\sqrt{2\|p^{*}\|_{0}\,}}{n}\log\left(\frac{4}{\delta}\right)\|\widehat{p}-p^{*}\|_{1}\enspace,$

since $\|p^{*}\|_{1}\leq\|p^{*}\|_{0}$. This leads us to

$\displaystyle\|p^{*}-\widehat{p}\|^{2}_{2}\leq\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}$
$\displaystyle+n\lambda\|Dp^{*}\|_{1}-n\lambda\|D\widehat{p}\|_{1}$
$\displaystyle+\frac{\sqrt{2\|p^{*}\|_{0}\,\log\left(\tfrac{4}{\delta}\right)}}{n}\|\widehat{p}-p^{*}\|_{1}\enspace.$

Using $n\lambda=\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)$ and the
triangle inequality we get

$\displaystyle\|p^{*}-\widehat{p}\|^{2}_{2}\leq
2n\lambda\left(\|D\left(\widehat{p}-p^{*}\right)_{T}\|_{1}+\|\left(Dp^{*}\right)_{T^{c}}\|_{1}\right)+\frac{2\|p^{*}\|_{0}\log\left(\tfrac{4}{\delta}\right)}{n}+\frac{1}{4}\|\widehat{p}-p^{*}\|^{2}_{2}\enspace.$

Now, applying Definition 2.1 of compatibility factor we can write

$\displaystyle\frac{1}{2}\|p^{*}-\widehat{p}\|^{2}_{2}\leq\dfrac{8\rho^{2}|T|}{\kappa^{2}_{T}}\log^{2}\left(\frac{4n^{2}}{\delta}\right)+2\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|\left(Dp^{*}\right)_{T^{c}}\|_{1}+\frac{2\|p^{*}\|_{0}\log\left(\tfrac{4}{\delta}\right)}{n}$

which proves the statement of Theorem 2.4.

###  A.1 Lemma A.1

######  Lemma A.1.

Fix $\delta\in[0,1]$. For $i\in[n]$, let $z_{i}$ be independent centered
Bernoulli random variables with parameter $q_{i}$ and $z=(z_{i})_{i\in[n]}$.
Then, with probability at least $1-\delta$,

$\|\left(D^{\dagger}\right)^{T}z\|_{\infty}\leq\sqrt{2}\rho\log\left(\frac{2n^{2}}{\delta}\right).$

######  Proof A.2.

We start by applying Bernstein’s inequality to
$\sum_{i=1}^{n}D^{\dagger}_{ij}z_{i}$. Note that
$\left|D^{\dagger}_{ij}z_{i}\right|\leq\rho$ and
$\sum_{i=1}^{n}\left(D^{\dagger}_{ij}\right)^{2}E\,z^{2}_{i}\leq\sum_{i=1}^{n}\left(D^{\dagger}_{ij}\right)^{2}q_{i}\leq\rho^{2}$.
Then, the Bernstein’s inequality implies that, with probability at least
$1-\frac{1}{n^{2}\delta}$,

$\left|\sum_{i=1}^{n}D^{\dagger}_{ij}z_{i}\right|\leq\sqrt{2}\rho\log\left(\frac{2n^{2}}{\delta}\right)\
.$

Now, applying the union bound and using $m\leq n^{2}$ we get the result of
Lemma A.1.

##  Appendix B Proof of Proposition 2.5

Given a $p\in\mathbb{R}^{n}$, let
$\mathcal{P}(p)=\left(\mathcal{P}(p_{i})\right)^{n}_{i=1}\in[0,1]^{n}$ be
defined as follows

$\mathcal{P}(p_{i})=\begin{cases}p_{i}&\text{if }\;p_{i}\in[0,1],\\\
0&\text{if }\;p_{i}<0,\\\ 1&\text{if }\;p_{i}>0.\\\ \end{cases}$

Let
$\mathcal{F}=\dfrac{1}{n}\overset{n}{\underset{i=1}{\sum}}\left(y^{k_{0}}_{i}-p_{i}\right)^{2}+\lambda\|Dp\|_{1}$.
Note that, for any $p$, $\mathcal{F}(\mathcal{P}(p))\leq\mathcal{F}(p)$. So,
without loss of generality, we assume that $\widehat{p}_{i}\in[0,1]^{n}$. Now
take $\bar{p}=(\bar{p}_{i})^{n}_{i=1}$ such that

$\bar{p}_{i}=\begin{cases}-1&\text{if }\;i\in S^{c},\\\ 0&\text{if not.
}\end{cases}$

Then, the first order optimality condition (16) implies

$\langle\bar{p},\left(p^{*}-\widehat{p}\right)_{S}\rangle-\langle\bar{p},\widehat{p}_{S^{c}}\rangle+\langle\bar{p},\xi\rangle=\frac{n\lambda}{2}\langle
D\bar{p},z\rangle.$

Note that $\xi_{S^{c}}=0$ which implies

$\|\widehat{p}_{S^{c}}\|_{1}=\frac{n\lambda}{2}\langle
D\bar{p},z\rangle\leq\frac{n\lambda}{2}\|D\bar{p}\|_{1}.$

We have that for any $e=(i,j)\in E$, $(D\bar{p})_{e}=0$ if $(i,j)\in S\times
S$ or if $(i,j)\in S^{c}\times S^{c}$ which implies that $\|D\bar{p}\|_{1}\leq
2sd_{\max}(S)$. Then we compute

$\displaystyle\|p^{*}-\widehat{p}\|_{1}=\|\left(p^{*}-\widehat{p}\right)_{S}\|_{1}+\|\widehat{p}_{S^{c}}\|_{1}\leq
n\lambda sd_{\max}(S)+\sqrt{s}\|p^{*}-\widehat{p}\|_{2}.$

Combining this bound with Theorem 2.4 we get the statement of Proposition 2.5.

##  Appendix C Proof of Theorem 4.1

For simplicity, we drop the $miss$ index, that is we write
$\widehat{p}_{miss}=\widehat{p}$. By first order optimality condition for the
convex problem (11), there exists $z\in\mathop{\rm sign}(D\widehat{p})$ such
that for any $\bar{p}\in\mathbb{R}$, we have

$\tfrac{2}{n}\left\langle\bar{p},\mathfrak{M}\circ\left(\widetilde{Y}-\widehat{p}\right)\right\rangle=\lambda\langle\bar{p},D^{T}z\rangle=\lambda\langle
D\bar{p},z\rangle\ ,$ (22)

which implies

$\tfrac{2}{n}\left\langle\widehat{p},\mathfrak{M}\circ\left(\widetilde{Y}-\widehat{p}\right)\right\rangle=\lambda\|D\widehat{p}\|_{1}\quad\text{and}$
$\tfrac{2}{n}\left\langle
p^{*},\mathfrak{M}\circ\left(\widetilde{Y}-\widehat{p}\right)\right\rangle\leq\lambda\|D{p^{*}}\|_{1}.$

Subtracting the above two, we obtain

$\tfrac{2}{n}\left\|\mathfrak{M}\circ\left(p^{*}-\widehat{p}\right)\right\|^{2}_{2}\leq\tfrac{2}{n}\langle\mathfrak{M}\circ\xi,\widehat{p}-p^{*}\rangle+\lambda\|Dp^{*}\|_{1}-\lambda\|D\widehat{p}\|_{1}.$

The control of the error term is similar to the case with complete
observations and we only sketch it. Let
$\eta=\mathfrak{M}\circ\xi=(\eta_{i})_{i\in[n]}$. Note that the $\eta_{i}$’s
are independent Bernoulli random variables with parameter $\pi_{i}p_{i}^{*}$.
Then, following similar argument as in the proof of Theorem 2.4 and using
Bernstein’s inequality, we get with probability at least $1-\delta/2$,

$\displaystyle\langle\Pi\eta,\Pi\left(\widehat{p}-p^{*}\right)\rangle$
$\displaystyle\leq\|\Pi\eta\|_{2}\left\|\Pi\left(\widehat{p}-p^{*}\right)\right\|_{2}$
(23) $\displaystyle\leq\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,}\log\left(\tfrac{4}{\delta}\right)}{n}\;\|\widehat{p}-p^{*}\|_{1}$
$\displaystyle\leq\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(\tfrac{4}{\delta}\right)}{n}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}$

Similarly, we have that with probability $1-\delta/2$

$\displaystyle\left\langle\left(D^{\dagger}\right)^{T}\eta,D\left(\widehat{p}-p^{*}\right)\right\rangle$
$\displaystyle\leq\|\left(D^{\dagger}\right)^{T}\eta\|_{\infty}\|D\left(\widehat{p}-p^{*}\right)\|_{1}$
$\displaystyle\leq\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}$
(24)

where the last inequality follows from Lemma A.1. Putting together (23) and
(C) we get that, with probability at least $1-\delta$,

$\displaystyle\langle\mathfrak{M}\circ\xi,\widehat{p}-p^{*}\rangle\leq$
$\displaystyle\sqrt{2}\rho\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}$
$\displaystyle\hskip 14.22636pt+\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(\tfrac{4}{\delta}\right)}{n}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}$

and

$\displaystyle\left\|\mathfrak{M}\circ\left(p^{*}-\widehat{p}\right)\right\|^{2}_{2}\leq\sqrt{2}\rho$
$\displaystyle\log\left(\frac{4n^{2}}{\delta}\right)\|D\left(\widehat{p}-p^{*}\right)\|_{1}+\dfrac{n\lambda}{2}\|Dp^{*}\|_{1}-\dfrac{n\lambda}{2}\|D\widehat{p}\|_{1}$
$\displaystyle+\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(\tfrac{4}{\delta}\right)}{n}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}\enspace.$

For $T={\rm supp}\left(Dp^{*}\right)$ and $\delta=4/n$, using
$n\lambda=9\sqrt{2}\rho\log\left(n\right)$, we can write

$\displaystyle\left\|\mathfrak{M}\circ\left(p^{*}-\widehat{p}\right)\right\|^{2}_{2}+\frac{n\lambda}{6}\left\|D\,\widehat{p}_{T^{C}}\right\|_{1}\leq$
$\displaystyle\frac{5n\lambda}{6}\left\|D(p^{*}-\widehat{p})_{T}\right\|_{1}$
(25) $\displaystyle\hskip 14.22636pt+\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(n\right)}{n}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}$

and

$\displaystyle\left\|\left(D\,\widehat{p}\right)_{T^{C}}\right\|_{1}\leq
5\left\|\left(D(p^{*}-\widehat{p})\right)_{T}\right\|_{1}+\frac{2\sqrt{\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(n\right)}{3n\rho\log\left(n\right)}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}.$
(26)

Next, we need to provide a lower bound for
$\left\|\mathfrak{M}\circ\left(p^{*}-\widehat{p}\right)\right\|^{2}_{2}$. It
can be done on a set of vectors that satisfy (26) if
$\left\|p^{*}-\widehat{p}\right\|^{2}_{l_{2}(\pi)}$ is not too small. We start
by introducing the corresponding set of vectors. Let

$\Delta=\frac{\sqrt{2\|\pi\circ p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}}{3n\rho}.$ (27)

For a set of indices $T$ and
$\epsilon\geq\dfrac{\log(n)}{0.01\,\log\left(6/5\right)}$, we denote by
$\mathcal{C}(T,\epsilon)$ the following set of vectors:

$\mathcal{C}(T,\epsilon)=\left\\{v\in[0,1]^{n}\,:\,\left\|v\right\|_{l_{2}(\pi)}^{2}\geq\epsilon,\|\left(Dv\right)_{T^{C}}\|_{1}\leq
5\|\left(Dv\right)_{T}\|_{1}+\Delta\|v\|_{l_{2}(\pi)}\right\\}.$

Next result provides partial isometry for vectors in
$\mathcal{C}(T,\epsilon)$:

######  Lemma C.1.

With probability higher than $1-8/n$, we have for all
$v\in\mathcal{C}(T,\epsilon)$

$\left\|\mathfrak{M}\circ
v\right\|^{2}_{2}\geq\frac{\|v\|_{l_{2}(\pi)}^{2}}{2}-77\left\\{\dfrac{36\rho\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}+\dfrac{\sqrt{\|\pi^{-1}\|_{1}}}{n}\left(5\sqrt{\|\pi\|_{1}}\log(n)+1\right)\right\\}^{2}.$

We now consider two cases, depending on whether the vector $p^{*}-\widehat{p}$
belongs to the set $\mathcal{C}(T,\epsilon)$ or not.
Case 1: Suppose first that
$\left\|p^{*}-\widehat{p}\right\|_{l_{2}(\pi)}^{2}<\epsilon$, then the
statement of Theorem 4.1 is true.

Case 2: It remains to consider the case
$\left\|p^{*}-\widehat{p}\right\|_{l_{2}(\pi)}^{2}\geq\epsilon$. Let $T={\rm
supp}(Dp^{*})$. Then, (26) implies that
$p^{*}-\widehat{p}\in\mathcal{C}(T,\epsilon)$ with probability at least
$1-4/n$. Then, combining (25) and Lemma C.1 we get

$\displaystyle\frac{1}{2}\|p^{*}-\widehat{p}\|^{2}_{l_{2}(\pi)}\leq
77\left\\{\dfrac{36\rho\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}+\dfrac{\sqrt{\|\pi^{-1}\|_{1}}}{n}\left(5\sqrt{\|\pi\|_{1}}\log(n)+1\right)\right\\}^{2}$
$\displaystyle+\frac{12\rho\,\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}\left\|p^{*}-\widehat{p}\right\|_{l_{2}(\pi)}+\frac{\sqrt{2\|\pi\circ
p^{*}\|_{1}\,\|\pi^{-1}\|_{1}}\log\left(n\right)}{n}\;\|\widehat{p}-p^{*}\|_{l_{2}(\pi)}.$
(28)

Now the statement of Theorem 4.1 follows applying $2ab\leq
8a^{2}+\tfrac{1}{8}b^{2}$ to the last two terms in (C).

###  C.1 Proof of Lemma C.1

Set $\nu=\dfrac{\log(n)}{0.01\,\log\left(6/5\right)}$, $\alpha=\dfrac{6}{5}$
and

$\mathcal{E}=77\left\\{\dfrac{36\rho\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}+\dfrac{\sqrt{\|\pi^{-1}\|_{1}}}{n}\left(5\sqrt{\|\pi\|_{1}}\log(n)+1\right)\right\\}^{2}.$

The proof of Lemma C.1 is based on the peeling argument. Denote by
$\mathcal{B}$ the set that contains the complement of the event that we are
interested in:

$\mathcal{B}=\left\\{\exists\,v\in\mathcal{C}(T,\epsilon)\,\text{such
that}\,\left|\left\|\mathfrak{M}\circ
v\right\|_{2}^{2}-\|v\|_{l_{2}(\pi)}^{2}\right|>\tfrac{1}{2}\|v\|_{l_{2}(\pi)}^{2}+\mathcal{E}\right\\}.$

For $r>\nu$ we will consider the following set of vectors:

$\mathcal{C}(T,\epsilon,r)=\left\\{v\in\mathcal{C}(T,\epsilon)\,:\,\left\|v\right\|_{l_{2}(\pi)}^{2}\leq
r\right\\}$

and, for $l\in\mathbb{N}$, the following events

$\mathcal{B}_{l}=\left\\{\exists\,v\in\mathcal{C}(T,\epsilon,\alpha^{l}\nu)\,:\,\left|\left\|\mathfrak{M}\circ
v\right\|_{2}^{2}-\|v\|_{L_{2}(\Pi)}^{2}\right|>\tfrac{5}{12}\alpha^{l}\nu+\mathcal{E}\right\\},$
$S_{l}=\left\\{v\in\mathcal{C}(T,\epsilon)\,:\,\alpha^{l-1}\nu\leq\|v\|_{l_{2}(\pi)}^{2}\leq\alpha^{l}\nu\right\\}.$

If the event $\mathcal{B}$ holds for some vector
$v\in\mathcal{C}(T,\epsilon)$, then $v$ belongs to some $S_{l}$ and

$\begin{split}\left|\left\|\mathfrak{M}\circ
v\right\|_{2}^{2}-\|v\|_{l_{2}(\pi)}^{2}\right|&>\tfrac{1}{2}\|v\|_{l_{2}(\pi)}^{2}+\mathcal{E}\\\
&>\tfrac{1}{2}\alpha^{l-1}\nu+\mathcal{E}\\\
&=\tfrac{5}{12}\alpha^{l}\nu+\mathcal{E}\end{split}$ (29)

which implies $\mathcal{B}\subset\cup\,\mathcal{B}_{l}$. Next, we apply the
union bound using the following lemma:

######  Lemma C.2.

Let

$Z_{r}=\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\left\|\mathfrak{M}\circ
v\right\|_{2}^{2}-\|v\|_{l_{2}(\pi)}^{2}\right|.$

We have that

$pr\left(Z_{r}\geq\frac{5}{12}r+\mathcal{E}\right)\leq 4e^{-c_{1}\,r}$

with $c_{1}\geq 0.01$.

Lemma C.2 implies that $\mathbb{P}\left(\mathcal{B}_{l}\right)\leq
4\exp(-c_{1}\,\alpha^{l}\nu)$. Using $e^{x}\geq x$, we obtain

$\begin{split}\mathbb{P}\left(\mathcal{B}\right)&\leq\underset{l=1}{\overset{\infty}{\Sigma}}\mathbb{P}\left(\mathcal{B}_{l}\right)\leq
4\underset{l=1}{\overset{\infty}{\Sigma}}\exp(-c_{1}\,\alpha^{l}\nu)\\\ &\leq
4\underset{l=1}{\overset{\infty}{\Sigma}}\exp\left(-c_{1}\,\nu\,\log(\alpha)\,l\right)\leq\dfrac{4\exp\left(-c_{1}\,\nu\,\log(\alpha)\right)}{1-\exp\left(-c_{1}\,\nu\,\log(\alpha)\right)}\\\
&=\dfrac{4\exp\left(-\log(n)\right)}{1-\exp\left(-\log(n)\right)}.\end{split}$

This completes the proof of Lemma C.1.

###  C.2 Proof of Lemma C.2

We first provide an upper bound on $\mathbb{E}(Z_{r})$ and then show that
$Z_{r}$ concentrates around its expectation. By definition of $Z_{r}$ we have

$Z_{r}=\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\sum_{i}m_{i}v_{i}^{2}-E\left(\sum_{i}m_{i}v_{i}^{2}\right)\right|.$

Using a standard symmetrization argument (see e.g. (Ledoux, 2001)) we obtain

$\begin{split}E\left(Z_{r}\right)&=E\left(\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\sum_{i}m_{i}v_{i}^{2}-E\left(m_{i}v_{i}^{2}\right)\right|\right)\\\
&\leq
2E\left(\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\sum_{i}\zeta_{i}m_{i}v_{i}^{2}\right|\right)\end{split}$

where $\\{\zeta_{i}\\}$ is an i.i.d. Rademacher sequence. Then, using that
$v_{i}\leq 1$, the contraction inequality (see e.g. (Koltchinskii, 2011,
Theorem 2.2)) yields

$E\left(Z_{r}\right)\leq
8E\left(\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\sum_{i}\zeta_{i}m_{i}v_{i}\right|\right)=8\,E\left(\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left|\left\langle\psi,v\right\rangle\right|\right)$

where $\psi=(\zeta_{i}m_{i})_{i\in[n]}$. We can write

$\displaystyle\left\langle\psi,v\right\rangle$
$\displaystyle=\langle\Pi\psi,v\rangle+\langle(I-\Pi)\psi,v\rangle=\langle\Pi\psi,\Pi
v\rangle+\left\langle\left(D^{\dagger}\right)^{T}\psi,Dv\right\rangle$ (30)

where $\Pi$ denote the projection on the kernel of $D$. To bound the first
term on the right hand side in (30), we first use the Hölder inequality:

$\displaystyle\left|\langle\Pi\psi,\Pi v\rangle\right|\leq\|\Pi\psi\|_{2}\|\Pi
v\|_{2}.$

Using the same argument as in the proof of Theorem 2.4, we have that

$\|\Pi\psi\|_{2}=\tfrac{1}{\sqrt{n}}|\langle\psi,\mathbf{1}_{n}\rangle|\quad\text{and}\quad\|\Pi
v\|_{2}=\tfrac{1}{\sqrt{n}}|\langle v,\mathbf{1}_{n}\rangle|\enspace.$

For $v\in\mathcal{C}(T,\epsilon,r)$ we get

$\|\Pi
v\|_{2}\leq\tfrac{1}{\sqrt{n}}\|v\|_{1}\leq\sqrt{\frac{\|\pi^{-1}\|_{1}}{n}}\left\|v\right\|_{l_{2}(\pi)}\leq\sqrt{\frac{\|\pi^{-1}\|_{1}\,r}{n}}\enspace.$
(31)

On the other hand, Bernstein’s inequality implies that that, with probability
at least $1-\frac{1}{n}$, we have

$\left|\sum_{i=1}^{n}\psi_{i}\right|\leq\sqrt{2\|\pi\|_{1}}\log\left(2n\right)$

and

$E\left(\tfrac{1}{\sqrt{n}}\left|\langle\psi,\mathbf{1}_{n}\rangle\right|\right)\leq\sqrt{\frac{2\|\pi\|_{1}}{n}}\log\left(2n\right)+\dfrac{1}{\sqrt{n}}$
(32)

where we used $\left|\langle\psi,\mathbf{1}_{n}\rangle\right|\leq n$. Putting
together (31) and (32) we get that for all

$\displaystyle\mathbb{E}\left[\sup_{v\in\mathcal{C}(T,\epsilon,r)}\left|\langle\Pi\psi,\Pi
v\rangle\right|\right]\leq\frac{\sqrt{\|\pi^{-1}\|_{1}\,r}}{n}\left(\sqrt{2\|\pi\|_{1}}\log\left(2n\right)+1\right).$
(33)

For the second term in (30), using again the Hölder inequality, we get

$\displaystyle\left|\left\langle\left(D^{\dagger}\right)^{T}\psi,Dv\right\rangle\right|$
$\displaystyle\leq\|\left(D^{\dagger}\right)^{T}\psi\|_{\infty}\|Dv\|_{1}.$
(34)

Lemma A.1 implies that, with probability at least $1-4/n$, we have

$\left\|\left(D^{\dagger}\right)^{T}\psi\right\|_{\infty}\leq
3\sqrt{2}\rho\log\left(8n\right).$ (35)

On the other hand, we have that

$\|\left(D^{\dagger}\right)^{T}\psi\|_{\infty}=\underset{j}{\sup}\left|\sum_{i=1}^{n}D^{\dagger}_{ij}\psi_{i}\right|\leq\underset{j}{\sup}\sqrt{\sum_{i}\left(D^{\dagger}_{ij}\right)^{2}}\sqrt{\sum_{i}\psi_{i}^{2}}\leq\rho\sqrt{n}$

which together with (35) imply

$E\left(\left\|\left(D^{\dagger}\right)^{T}\psi\right\|_{\infty}\right)\leq
3\sqrt{2}\rho\log\left(n\right)+\dfrac{4\rho}{\sqrt{n}}\leq
6\rho\log\left(n\right)$ (36)

for $n\geq 4$. Note that, as $v\in\mathcal{C}(T,\epsilon,r)$, we also obtain

$\displaystyle\|Dv\|_{1}$
$\displaystyle\leq\|\left(Dv\right)_{T}\|_{1}+\|\left(Dv\right)_{T^{C}}\|_{1}\leq
6\|\left(Dv\right)_{T}\|_{1}+\Delta\|v\|_{l_{2}(\pi)}$
$\displaystyle\leq\dfrac{6\sqrt{|T|}}{\kappa_{T}}\|v\|_{2}+\Delta\sqrt{r}$
$\displaystyle\leq\dfrac{6\sqrt{\kappa_{\pi}|T|}}{\kappa_{T}}\|v\|_{l_{2}(\pi)}+\Delta\sqrt{r}\quad\text{using
\eqref{relation_norms}}$
$\displaystyle\leq\left(\dfrac{6\sqrt{\kappa_{\pi}|T|}}{\kappa_{T}}+\Delta\right)\sqrt{r}.$
(37)

Using (36) and (C.2) we get

$\displaystyle
E\left(\underset{v\in\mathcal{C}(T,\epsilon,r)}{\sup}\left\langle\left(D^{\dagger}\right)^{T}\psi,Dv\right\rangle\right)\leq
6\rho\sqrt{r}\log(n)\left(\dfrac{6\sqrt{\kappa_{\pi}|T|}}{\kappa_{T}}+\Delta\right).$
(38)

Putting together (33), (38) and using the definition (27) of $\Delta$, we get

$E\left(Z_{r}\right)\leq
8\sqrt{r}\left\\{\dfrac{36\rho\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}+\dfrac{\sqrt{\|\pi^{-1}\|_{1}}}{n}\left(5\sqrt{\|\pi\|_{1}}\log(n)+1\right)\right\\}$

where we used $\|\pi\circ p^{*}\|_{1}\leq\|\pi\|_{1}$. Now, $2ab\leq
a^{2}+b^{2}$ implies

$\displaystyle E\left(Z_{r}\right)$
$\displaystyle\leq\dfrac{5}{24}r+77\left\\{\dfrac{36\rho\sqrt{\kappa_{\pi}|T|}\log(n)}{\kappa_{T}}+\dfrac{\sqrt{\|\pi^{-1}\|_{1}}}{n}\left(5\sqrt{\|\pi\|_{1}}\log(n)+1\right)\right\\}^{2}$
$\displaystyle\hskip 14.22636pt=\dfrac{5}{24}r+\mathcal{E}$ (39)

Now we can show that $Z_{r}$ concentrates around its expectation using
Talagrand’s concentration inequality (Talagrand, 1996), which in the current
can be obtained by inverting the tail bound in Theorem 3.3.16 in (Giné &
Nickl, 2016).

######  Theorem C.3.

Let $(S,\mathcal{S})$ be a measurable space and let $n\in\mathbb{N}$. Let
$X_{k},~{}k=1,\dots,n$ be independent ${S}$-valued random variables and let
$\mathcal{F}$ be a countable set of functions
$f=(f_{1},...,f_{n}):{S}^{n}\rightarrow[-K,K]^{n}$ such that
$\mathbb{E}f_{k}(X_{k})=0$ for all $f\in\mathcal{F}$ and $k=1,...,n$. Set

$Z:=\sup_{f\in\mathcal{F}}\sum_{k=1}^{n}f_{k}(X_{k})\enspace.$

Define the variance proxy

$V_{n}:=2K\mathbb{E}Z+\sup_{f\in\mathcal{F}}\sum_{k=1}^{n}\mathbb{E}\left[(f_{k}(X_{k}))^{2}\right].$

Then, for all $t\geq 0$,

$pr\left(Z-\mathbb{E}Z\geq
t\right)\leq\exp\left(\frac{-t^{2}}{4V_{n}+(9/2)Kt}\right).$

Note that the functional $v\to\left\|\mathfrak{M}\circ
v\right\|_{2}^{2}-\|v\|_{l_{2}(\pi)}^{2}$ is continuous on the set of vectors
$\mathcal{C}(T,\epsilon,r)$. Hence it’s enough to consider a dense countable
subset of $\mathcal{C}(T,\epsilon,r)$ and we may apply Talagrand’s inequality
to $Z_{r}$. We have for our particular case, since
$\sup_{f\in\mathcal{F}}|f(X)|=\sup_{f\in\\{\mathcal{F}\bigcup\\{-\mathcal{F}\\}\\}}f(X)$,
that

$X_{i}=m_{i}-\mathbb{E}\left(m_{i}\right),~{}~{}~{}S=[-1,1].$

We can compute:

$\displaystyle\sum_{i}\mathbb{E}\left[(f_{k}(X_{k}))^{2}\right]=\sum_{i}\mathbb{E}\left[m_{i}-\mathbb{E}\left(m_{i}\right)\right]^{2}v^{4}_{i}\leq\sum_{i}\pi_{i}v^{2}_{i}=\|v\|_{l_{2}(\pi)}^{2}\leq
r$

which implies that $V_{n}\leq\tfrac{17}{12}r+2\mathcal{E}$. Hence, using (C.2)
and Talagrand’s inequality with $t=\dfrac{5}{24}r+0.5\,\mathcal{E}$, we obtain

$\displaystyle
pr\left(Z_{r}>\frac{5}{12}r+1.5\mathcal{E}\right)\leq\exp(-c_{1}r)$

with $c_{1}\geq 0.01$ which completes the proof of Lemma C.2.

##  Appendix D Further Results: Simulation study

In this section, we report additional results of our simulation study.

###  D.1 Algorithms and Compute time

In our simulation study, we use the semismooth Newton based augmented
Lagrangian method (SSNAL) proposed by Sun et al. (2021). SSNAL is a scalable
algorithm that was originally designed to perform convex clustering to solve
the following optimization problem:

$\min_{X\in\mathbb{R}^{d\times
n}}\frac{1}{2}\sum_{i=1}^{n}\|x_{i}-a_{i}\|^{2}+\lambda\sum_{i<j}W_{ij}\|x_{i}-x_{j}\|_{2}$

Note that in our case, since our vector $p$ is one dimensional
$\|p_{i}-p_{j}\|_{2}=\sqrt{(p_{i}-p_{j})^{2}}=|p_{i}-p_{j}|$, and our problem
shares the same objective function as convex clustering. In their work, Sun et
al. (2021) prove their algorithm to be theoretically efficient and practically
efficient and robust. More specifically, the authors show that, as long as the
dimension of the feature vector is small (i.e, in our case, equal to 1), this
method has the potential of considerably improving upon existing solvers.

Evaluation. To highlight the scalability of SSNAL, we propose comparing here
four different solvers on a set of synthetic experiments. More specifically,
we compare the performances of:

  1. 1.

An off-the-shelf convex optimization solver (CVXPY),

  2. 2.

The dual-coordinate descent algorithm proposed by Tibshirani and Taylor
Tibshirani & Taylor (2011). This method relies on solving the dual problem:

$\text{minimize}\frac{1}{2}\|y-\Gamma^{T}u\|^{2}\text{ subject to
}\|u\|_{\infty}\leq\lambda,$

yielding a simple algorithm.

  3. 3.

An ADMM-based approach: ADMM is often the method of choice to deal with large
scale constrained and has been suggested in multiple instance as the method of
choice for graph-based penalties on large graphs Hallac et al. (2017),

  4. 4.

The SSNAL algorithm of Sun et al. (2021), described in the previous paragraph.

In this first set of experiments, we propose comparing the performance of the
different algorithms when the underlying is an Erdos Renyi random graph. The
reason for considering this topology here is that this is a difficult setting,
since the graph is quite dense. Since the complexity of the penalty depends on
the number of edges, we plot the time requirements of the algorithm as the
number of nodes increases and the penalty varies. The results are displayed in
Figure 9 below.

Figure 9: Computation time for the different algorithms

###  D.2 Synthetic Experiments: Additional Results

Figure 10: Number of infections for the kNN graph as a function of $\beta$,
$n_{\text{steps}}$ and the parameter $k$ in the k-NN graph. The healing rate
is here fixed to $\gamma=0.1.$ Results are averaged over 100 independent
experiments. Figure 11: $\ell_{1}$ error for various graph topologies as a
function of $\beta$, $n_{\text{steps}}$. The healing rate is here fixed to
$\gamma=0.1.$ Results are averaged over 100 independent experiments, and error
bars denote interquartile ranges. Figure 12: $\ell_{2}$ error for various
graph topologies as a function of $\beta$, $n_{\text{steps}}$. The healing
rate is here fixed to $\gamma=0.1.$ Results are averaged over 100 independent
experiments, and error bars denote interquartile ranges. Figure 13: Number of
infections for graphs with various topologies as a function of $\beta$,
$n_{\text{steps}}$. The healing rate is here fixed to $\gamma=0.1.$ Results
are averaged over 100 independent experiments, and error bars denote
interquartile ranges.

## References

  * Ahn & Hassibi (2013) Ahn, H. J. & Hassibi, B. (2013).  Global dynamics of epidemic spread over complex networks.  In 52nd IEEE Conference on Decision and Control.
  * Aiello et al. (2016) Aiello, A. E., Simanek, A. M., Eisenberg, M. C., Walsh, A. R., Davis, B., Volz, E., Cheng, C., Rainey, J. J., Uzicanin, A., Gao, H. et al. (2016).  Design and methods of a social network isolation study for reducing respiratory infection transmission: The ex-flu cluster randomized trial.  Epidemics 15, 38–55.
  * Andersson & Britton (2012) Andersson, H. & Britton, T. (2012).  Stochastic epidemic models and their statistical analysis, vol. 151.  Springer Science & Business Media.
  * Andrieu et al. (2010) Andrieu, C., Doucet, A. & Holenstein, R. (2010).  Particle markov chain monte carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 269–342.
  * Arnold & Tibshirani (2016) Arnold, T. B. & Tibshirani, R. J. (2016).  Efficient Implementations of the Generalized Lasso Dual Path Algorithm.  Journal of Computational and Graphical Statistics 25, 1–27.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/10618600.2015.1008638.
  * Balcan et al. (2009) Balcan, D., Colizza, V., Gonçalves, B., Hu, H., Ramasco, J. J. & Vespignani, A. (2009).  Multiscale mobility networks and the spatial spreading of infectious diseases.  Proceedings of the National Academy of Sciences 106, 21484–21489.  Publisher: Proceedings of the National Academy of Sciences.
  * Ball et al. (2019) Ball, F., Britton, T., Pardoux, E., Larédo, C., Sirl, D. & Tran, V. (2019).  Stochastic Epidemic Models with Inference.  Lecture Notes in Mathematics. Springer International Publishing.
  * Bickel et al. (2013) Bickel, P., Choi, D., Chang, X. & Zhang, H. (2013).  Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels.  The Annals of Statistics 41, 1922–1943.  Publisher: Institute of Mathematical Statistics.
  * Bowman et al. (2005) Bowman, C., Gumel, A., Van den Driessche, P., Wu, J. & Zhu, H. (2005).  A mathematical model for assessing control strategies against west nile virus.  Bulletin of mathematical biology 67, 1107–1133.
  * Boyd et al. (2011) Boyd, S., Parikh, N. & Chu, E. (2011).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc.
  * Britton (2010) Britton, T. (2010).  Stochastic epidemic models: a survey.  Mathematical biosciences 225, 24–35.
  * Britton & Giardina (2016) Britton, T. & Giardina, F. (2016).  Introduction to statistical inference for infectious diseases.  Journal de la société française de statistique 157, 53–70.
  * Bu et al. (2022) Bu, F., Aiello, A. E., Xu, J. & Volfovsky, A. (2022).  Likelihood-based inference for partially observed epidemics on dynamic networks.  Journal of the American Statistical Association 117, 510–526.
  * Carcione et al. (2020) Carcione, J. M., Santos, J. E., Bagaini, C. & Ba, J. (2020).  A simulation of a covid-19 epidemic based on a deterministic seir model.  Frontiers in public health 8, 230.
  * Celisse et al. (2012) Celisse, A., Daudin, J.-J. & Pierre, L. (2012).  Consistency of maximum-likelihood and variational estimators in the stochastic block model.  Electronic Journal of Statistics 6, 1847–1899.  Publisher: Institute of Mathematical Statistics and Bernoulli Society.
  * Chakrabarti et al. (2008) Chakrabarti, D., Wang, Y., Wang, C., Leskovec, J. & Faloutsos, C. (2008).  Epidemic thresholds in real networks.  ACM Transactions on Information and System Security (TISSEC) 10, 1–26.
  * Chakraborty et al. (2022) Chakraborty, T., Ghosh, I., Mahajan, T. & Arora, T. (2022).  Nowcasting of covid-19 confirmed cases: Foundations, trends, and challenges.  Modeling, Control and Drug Development for COVID-19 Outbreak Prevention , 1023–1064.
  * Dalalyan et al. (2017) Dalalyan, A. S., Hebiri, M. & Lederer, J. (2017).  On the prediction performance of the Lasso.  Bernoulli 23, 552–581.  Publisher: [Bernoulli Society for Mathematical Statistics and Probability, International Statistical Institute (ISI)].
  * Danon et al. (2011) Danon, L., Ford, A. P., House, T., Jewell, C. P., Keeling, M. J., Roberts, G. O., Ross, J. V., Vernon, M. C. et al. (2011).  Networks and the epidemiology of infectious disease.  Interdisciplinary perspectives on infectious diseases 2011.
  * Desai et al. (2019) Desai, A. N., Kraemer, M. U., Bhatia, S., Cori, A., Nouvellet, P., Herringer, M., Cohn, E. L., Carrion, M., Brownstein, J. S., Madoff, L. C. et al. (2019).  Real-time epidemic forecasting: challenges and opportunities.  Health security 17, 268–275.
  * Diamond & Boyd (2016) Diamond, S. & Boyd, S. (2016).  Cvxpy: A python-embedded modeling language for convex optimization.  Journal of Machine Learning Research 17, 1–5.
  * Draief & Massoulié (2009) Draief, M. & Massoulié, L. (2009).  Epidemics and Rumours in Complex Networks.  London Mathematical Society Lecture Note Series. Cambridge University Press.
  * Du et al. (2020) Du, Z., Xu, X., Wu, Y., Wang, L., Cowling, B. J. & Meyers, L. A. (2020).  Serial interval of covid-19 among publicly reported confirmed cases.  Emerging infectious diseases 26, 1341.
  * Eikenberry et al. (2020) Eikenberry, S. E., Mancuso, M., Iboi, E., Phan, T., Eikenberry, K., Kuang, Y., Kostelich, E. & Gumel, A. B. (2020).  To mask or not to mask: Modeling the potential for face mask use by the general public to curtail the covid-19 pandemic.  Infectious disease modelling 5, 293–308.
  * Ganesh et al. (2005) Ganesh, A., Massoulié, L. & Towsley, D. (2005).  The effect of network topology on the spread of epidemics.  In Proceedings IEEE 24th Annual Joint Conference of the IEEE Computer and Communications Societies., vol. 2. IEEE.
  * Gao et al. (2015) Gao, C., Lu, Y. & Zhou, H. H. (2015).  Rate-Optimal Graphon Estimation.  The Annals of Statistics 43, 2624–2652.  Publisher: Institute of Mathematical Statistics.
  * Gaucher & Klopp (2021) Gaucher, S. & Klopp, O. (2021).  Maximum likelihood estimation of sparse networks with missing observations.  Journal of Statistical Planning and Inference 215, 299–329.
  * Giné & Nickl (2016) Giné, E. & Nickl, R. (2016).  Mathematical foundations of infinite-dimensional statistical methods.  Cambridge University Press .
  * Graham (2020) Graham, F. (2020).  Daily briefing: Why the coronavirus death toll in wuhan, china, has just jumped by 50.  Nature (Lond.) .
  * Gómez et al. (2010) Gómez, S., Arenas, A., Borge-Holthoefer, J., Meloni, S. & Moreno, Y. (2010).  Discrete-time Markov chain approach to contact-based disease spreading in complex networks.  Europhysics Letters 89, 38009.
  * Hackl & Dubernet (2019) Hackl, J. & Dubernet, T. (2019).  Epidemic spreading in urban areas using agent-based transportation models.  Future internet 11, 92.
  * Hallac et al. (2017) Hallac, D., Park, Y., Boyd, S. & Leskovec, J. (2017).  Network inference via the time-varying graphical lasso.  In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.
  * Holme & Saramäki (2012) Holme, P. & Saramäki, J. (2012).  Temporal networks.  Physics Reports 519, 97–125.
  * Hu et al. (2017) Hu, B., Gonzales, J. L. & Gubbins, S. (2017).  Bayesian inference of epidemiological parameters from transmission experiments.  Scientific Reports 7, 16774.
  * Huang et al. (2022) Huang, J., Morsomme, R., Dunson, D. & Xu, J. (2022).  Detecting changes in the transmission rate of a stochastic epidemic model.  arXiv preprint arXiv:2211.14691 .
  * Hütter & Rigollet (2016) Hütter, J.-C. & Rigollet, P. (2016).  Optimal rates for total variation denoising.  In Conference on Learning Theory. PMLR.  ISSN: 1938-7228.
  * Ionides et al. (2011) Ionides, E. L., Bhadra, A., Atchadé, Y. & King, A. (2011).  Iterated filtering.  The Annals of Statistics 39, 1776–1802.
  * Jewell et al. (2009) Jewell, C. P., Kypraios, T., Neal, P. & Roberts, G. O. (2009).  Bayesian analysis for emerging infectious diseases.  Bayesian Analysis 4, 465–496.  Publisher: International Society for Bayesian Analysis.
  * Johns Hopkins CRC (2024) Johns Hopkins CRC, . (2024).  New covid-19 cases worldwide.  Available online:https://coronavirus.jhu.edu/data/new-cases.  Accessed on 29 April 2024.
  * Keeling (2005) Keeling, M. (2005).  The implications of network structure for epidemic dynamics.  Theoretical population biology 67, 1–8.
  * Keeling & Eames (2005a) Keeling, M. J. & Eames, K. T. (2005a).  Networks and epidemic models.  Journal of the royal society interface 2, 295–307.
  * Keeling & Eames (2005b) Keeling, M. J. & Eames, K. T. (2005b).  Networks and epidemic models.  Journal of The Royal Society Interface 2, 295–307.  Publisher: Royal Society.
  * Kenah & Robins (2007) Kenah, E. & Robins, J. M. (2007).  Network-based analysis of stochastic SIR epidemic models with random and proportionate mixing.  Journal of theoretical biology 249, 706–722.
  * Kermack & McKendrick (1927) Kermack, W. O. & McKendrick, A. G. (1927).  A contribution to the mathematical theory of epidemics.  Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115, 700–721.
  * Kim et al. (2021) Kim, K., Yoo, S., Lee, S., Lee, D. & Lee, K.-H. (2021).  Network analysis to identify the risk of epidemic spreading.  Applied Sciences 11, 2997.
  * King et al. (2016) King, A. A., Nguyen, D. & Ionides, E. L. (2016).  Statistical inference for partially observed markov processes via the r package pomp.  Journal of Statistical Software 69, 1–43.
  * Klopp et al. (2017) Klopp, O., Tsybakov, A. B. & Verzelen, N. (2017).  Oracle Inequalities for Network Models and Sparse Graphon Estimation.  The Annals of Statistics 45, 316–354.  Publisher: Institute of Mathematical Statistics.
  * Klopp & Verzelen (2019) Klopp, O. & Verzelen, N. (2019).  Optimal graphon estimation in cut distance.  Probability Theory and Related Fields 174, 1033–1090.
  * Koltchinskii (2011) Koltchinskii, V. (2011).  Oracle inequalities in empirical risk minimization and sparse recovery problems, vol. 2033 of Lecture Notes in Mathematics.  Springer, Heidelberg.  Lectures from the 38th Probability Summer School held in Saint-Flour, 2008, École d’Été de Probabilités de Saint-Flour. [Saint-Flour Probability Summer School].
  * Kumar et al. (2021) Kumar, A., Choi, T.-M., Wamba, S. F., Gupta, S. & Tan, K. H. (2021).  Infection vulnerability stratification risk modelling of covid-19 data: a deterministic seir epidemic model analysis.  Annals of Operations Research , 1–27.
  * Ledoux (2001) Ledoux, M. (2001).  The concentration of measure phenomenon, vol. 89 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI.
  * Lipsitch et al. (2015) Lipsitch, M., Donnelly, C. A., Fraser, C., Blake, I. M., Cori, A., Dorigatti, I., Ferguson, N. M., Garske, T., Mills, H. L., Riley, S. et al. (2015).  Potential biases in estimating absolute and relative case-fatality risks during outbreaks.  PLoS neglected tropical diseases 9, e0003846.
  * Lloyd-Smith et al. (2005) Lloyd-Smith, J. O., Schreiber, S. J., Kopp, P. E. & Getz, W. M. (2005).  Superspreading and the effect of individual variation on disease emergence.  Nature 438, 355–359.  Number: 7066 Publisher: Nature Publishing Group.
  * Mammen & Geer (1997) Mammen, E. & Geer, S. v. d. (1997).  Locally adaptive regression splines.  The Annals of Statistics 25, 387–413.  Publisher: Institute of Mathematical Statistics.
  * Manzo (2020) Manzo, G. (2020).  Complex social networks are missing in the dominant covid-19 epidemic models.  Sociologica 14, 31–49.
  * McGough et al. (2020) McGough, S. F., Johansson, M. A., Lipsitch, M. & Menzies, N. A. (2020).  Nowcasting by bayesian smoothing: A flexible, generalizable model for real-time epidemic tracking.  PLoS computational biology 16, e1007735.
  * Miller & Volz (2013) Miller, J. C. & Volz, E. M. (2013).  Incorporating Disease and Population Structure into Models of SIR Disease in Contact Networks.  PLOS ONE 8, e69162.  Publisher: Public Library of Science.
  * Moreno et al. (2002) Moreno, Y., Pastor-Satorras, R. & Vespignani, A. (2002).  Epidemic outbreaks in complex heterogeneous networks.  The European Physical Journal B-Condensed Matter and Complex Systems 26, 521–529.
  * Morsomme & Xu (2022) Morsomme, R. & Xu, J. (2022).  Exact inference for stochastic epidemic models via uniformly ergodic block sampling.  arXiv preprint arXiv:2201.09722 .
  * Mubarak et al. (2021) Mubarak, M., Berneburg, J. & Nowzari, C. (2021).  Stochastic vs. deterministic modeling for the spread of covid-19 in small networks.  In 2021 American Control Conference (ACC). IEEE.
  * Narci et al. (2021) Narci, R., Delattre, M., Larédo, C. & Vergu, E. (2021).  Inference for partially observed epidemic dynamics guided by kalman filtering techniques.  Computational Statistics and Data Analysis 164, 107319.
  * Needell & Ward (2013) Needell, D. & Ward, R. (2013).  Near-Optimal Compressed Sensing Guarantees for Total Variation Minimization.  IEEE Transactions on Image Processing 22, 3941–3949.  Conference Name: IEEE Transactions on Image Processing.
  * Newman (2002a) Newman, M. E. (2002a).  Spread of epidemic disease on networks.  Physical review E 66, 016128.
  * Newman (2002b) Newman, M. E. J. (2002b).  Spread of epidemic disease on networks.  Physical Review E 66, 016128.  Publisher: American Physical Society.
  * Ottaviano et al. (2018) Ottaviano, S., De Pellegrini, F., Bonaccorsi, S. & Van Mieghem, P. (2018).  Optimal curing policy for epidemic spreading over a community network with heterogeneous population.  Journal of Complex Networks 6, 800–829.
  * Padilla et al. (2018) Padilla, O. H. M., Sharpnack, J., Scott, J. G. & Tibshirani, R. J. (2018).  The DFS Fused Lasso: Linear-Time Denoising over General Graphs.  Journal of Machine Learning Research 18, 1–36.
  * Paré et al. (2020a) Paré, P. E., Beck, C. L. & Başar, T. (2020a).  Modeling, estimation, and analysis of epidemics over networks: An overview.  Annual Reviews in Control 50, 345–360.
  * Paré et al. (2020b) Paré, P. E., Liu, J., Beck, C. L., Kirwan, B. E. & Başar, T. (2020b).  Analysis, Estimation, and Validation of Discrete-Time Epidemic Processes.  IEEE Transactions on Control Systems Technology 28, 79–93.  Conference Name: IEEE Transactions on Control Systems Technology.
  * Pastor-Satorras et al. (2015a) Pastor-Satorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. (2015a).  Epidemic processes in complex networks.  Reviews of modern physics 87, 925.
  * Pastor-Satorras et al. (2015b) Pastor-Satorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. (2015b).  Epidemic processes in complex networks.  Rev. Mod. Phys. 87, 925–979.
  * Pastor-Satorras & Vespignani (2001) Pastor-Satorras, R. & Vespignani, A. (2001).  Epidemic Spreading in Scale-Free Networks.  Physical Review Letters 86, 3200.  Publisher: American Physical Society.
  * Pellis et al. (2015) Pellis, L., Ball, F., Bansal, S., Eames, K., House, T., Isham, V. & Trapman, P. (2015).  Eight challenges for network epidemic models.  Epidemics 10, 58–62.
  * Perez & Dragicevic (2009) Perez, L. & Dragicevic, S. (2009).  An agent-based approach for modeling dynamics of contagious disease spread.  International journal of health geographics 8, 1–17.
  * Rocha et al. (2011) Rocha, L. E. C., Liljeros, F. & Holme, P. (2011).  Simulated Epidemics in an Empirical Spatiotemporal Network of 50,185 Sexual Contacts.  PLoS Computational Biology 7, e1001109.
  * Rosenfeld & Tibshirani (2021) Rosenfeld, R. & Tibshirani, R. J. (2021).  Epidemic tracking and forecasting: Lessons learned from a tumultuous year.  Proceedings of the National Academy of Sciences 118, e2111456118.
  * Rossi & Ahmed (2015) Rossi, R. A. & Ahmed, N. K. (2015).  The network data repository with interactive graph analytics and visualization.  In AAAI.
  * Rudin et al. (1992) Rudin, L. I., Osher, S. & Fatemi, E. (1992).  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena 60, 259–268.
  * Sanchez et al. (2022) Sanchez, F., Calvo, J. G., Mery, G., García, Y. E., Vásquez, P., Barboza, L. A., Pérez, M. D. & Rivas, T. (2022).  A multilayer network model of covid-19: Implications in public health policy in costa rica.  Epidemics 39, 100577.
  * Shirley & Rushton (2005) Shirley, M. D. & Rushton, S. P. (2005).  The impacts of network topology on disease spread.  Ecological Complexity 2, 287–299.
  * Siegenfeld et al. (2020) Siegenfeld, A. F., Taleb, N. N. & Bar-Yam, Y. (2020).  What models can and cannot tell us about covid-19.  Proceedings of the National Academy of Sciences 117, 16092–16095.
  * Siettos et al. (2015) Siettos, C., Anastassopoulou, C., Russo, L., Grigoras, C. & Mylonakis, E. (2015).  Modeling the 2014 ebola virus epidemic–agent-based simulations, temporal analysis and future predictions for liberia and sierra leone.  PLoS currents 7.
  * Sisson et al. (2007) Sisson, S. A., Fan, Y. & Tanaka, M. M. (2007).  Sequential monte carlo without likelihoods.  Proceedings of the National Academy of Sciences 104, 1760 – 1765.
  * Spricer & Britton (2019) Spricer, K. & Britton, T. (2019).  An sir epidemic on a weighted network.  Network Science 7, 556–580.
  * Sun et al. (2021) Sun, D., Toh, K.-C. & Yuan, Y. (2021).  Convex clustering: Model, theoretical guarantee and efficient algorithm.  The Journal of Machine Learning Research 22, 427–458.
  * Talagrand (1996) Talagrand, M. (1996).  New concentration inequalities in product spaces.  Inventiones Mathematicae 126, 505–563.
  * Teweldemedhin et al. (2004) Teweldemedhin, E., Marwala, T. & Mueller, C. (2004).  Agent-based modelling: a case study in hiv epidemic.  In Fourth International Conference on Hybrid Intelligent Systems (HIS’04). IEEE.
  * Tibshirani & Taylor (2011) Tibshirani, R. J. & Taylor, J. (2011).  The solution path of the generalized lasso.  The annals of statistics 39, 1335–1371.
  * Tolles & Luong (2020) Tolles, J. & Luong, T. (2020).  Modeling epidemics with compartmental models.  Jama 323, 2515–2516.
  * Toni et al. (2009) Toni, T., Welch, D., Strelkowa, N., Ipsen, A. & Stumpf, M. P. H. (2009).  Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems.  Journal of the Royal Society, Interface 6, 187—202.
  * Volz & Meyers (2007) Volz, E. & Meyers, L. A. (2007).  Susceptible–infected–recovered epidemics in dynamic contact networks.  Proceedings of the Royal Society B: Biological Sciences 274, 2925–2933.
  * Vrabac et al. (2020) Vrabac, D., Paré, P. E., Sandberg, H. & Johansson, K. H. (2020).  Overcoming Challenges for Estimating Virus Spread Dynamics from Data.  In 2020 54th Annual Conference on Information Sciences and Systems (CISS).
  * Wan et al. (2014) Wan, X., Liu, J., Cheung, W. K. & Tong, T. (2014).  Inferring epidemic network topology from surveillance data.  PLoS One 9, e100661.
  * Wang et al. (2013) Wang, H., Li, Q., D’Agostino, G., Havlin, S., Stanley, H. E. & Van Mieghem, P. (2013).  Effect of the interconnected network structure on the epidemic threshold.  Physical Review E 88, 022801.
  * Wang et al. (2016) Wang, Y.-X., Sharpnack, J., Smola, A. J. & Tibshirani, R. J. (2016).  Trend Filtering on Graphs.  Journal of Machine Learning Research 17, 1–41.
  * World Health Organization (2024) World Health Organization, . (2024).  Daily covid-19 cases and deaths.  Data processed by Our World in Data. Available online:https://ourworldindata.org/grapher/daily-covid-cases-deaths-7-day-ra.  Accessed: 29 April 2024.
  * Wu et al. (2021) Wu, J. T., Leung, K., Lam, T. T., Ni, M. Y., Wong, C. K., Peiris, J. M. & Leung, G. M. (2021).  Nowcasting epidemics of novel pathogens: lessons from covid-19.  Nature medicine 27, 388–395.
  * Zhang & Moura (2014) Zhang, J. & Moura, J. M. F. (2014).  Diffusion in Social Networks as SIS Epidemics: Beyond Full Mixing and Complete Graphs.  IEEE Journal of Selected Topics in Signal Processing 8, 537–551.  Conference Name: IEEE Journal of Selected Topics in Signal Processing.