2021

[1]Juan A. Ponciano

[1]Instituto de Investigación en Ciencias Físicas y Matemáticas, Universidad
San Carlos de Guatemala, Guatemala City, Guatemala

2]National Advanced Computing Laboratory, National High Technology Center, San
José, Costa Rica

3]Department of Biology, University of Florida, Gainesville, FL 32611,
Florida, USA

# COVID-19 Regional Waves and Spread Risk Assessment through the Analysis of
the Initial Outbreak in Guatemala

<EMAIL_ADDRESS>Juan D. Chang  Mariela Abdalah  Kevin Facey  José
M. Ponciano  *  [  [

###### Abstract

The initial surge of the COVID-19 pandemic hit Guatemala on March 2020. On a
country scale, the epidemic has undergone a fairly well-known and
distinguishable initial phase, reaching its peak on mid July 2020. However,
the detailed picture is more involved and reflects inter-regional variations
in the epidemic dynamics, presumably grounded on socio-demographic,
connectivity, and human mobility factors. Classifying the regional epidemic
curves and identifying the major hubs of regional COVID-19 spread can
contribute towards defining an evidence-based risk map for future outbreaks of
infectious diseases with similar transmissibility properties. In this work, we
make a regional wave decomposition of the initial epidemic phase registered in
Guatemala, and we use the Richards phenomenological model alongside
multivariate ordination techniques of its estimated model parameters to draw a
countrywide picture of the first epidemiological wave. By exploring
similarities in the model space parameters, we traced routes for the disease
spread across the country. We evaluated how well the proposed classification
can help to define a regional risk hierarchy comprising early stage focal
points, major hubs, and secondary regions of epidemic progression.

###### keywords:

COVID-19, Epidemic waves, Richards Model, Principal Component Analysis

##  1 Introduction

The novel coronavirus (SARS-CoV-2) outbreak has left a worldwide spatial and
temporal progression landscape of a pervasive communicable disease
whoDashboard . The scattered historical epidemic data represents a tremendous
volume of information about an explosive viral infection among individuals
who, at an earlier stage, were mostly uninformed and unaware about the
potential harm of the virus. Our focus is on the epidemic development in
countries with little to no resources to implement a strong epidemiological
response. Although in these countries the epidemic information has been
unevenly collected and presents several limitations, here we recognize its
potential value to explore the course of this epidemic and of future ones.

Shortly after the pandemic started, a large amount of scientific researchers
from many disciplines began to study the collected data and raise the
discussion about issues concerning the viral epidemiology and the evolution of
the pandemic to a thorough level of analysis yang2020modified ;
yan2020interpretable ; chatterjee2020147 ; Bertozzi16732 ; shakeel2021covid .
Despite the reach of detailed studies that make use of modern human mobility
information/databases, assessing the spread risk of COVID-19 in countries
where such information is not available remains a challenging goal. Even with
data from GPS technology from mobile phones, assessing spread risks has proven
to be puzzling given the high complexity of the structure of modern social
systems. Practical and theoretical challenges to characterize epidemic growth
and the fast adaptive evolution of the virus in the host population are all
factors that make of accurate risk assessments a difficult target. If these
factors are carefully studied and comprehended, a better understanding of the
infection behavior could be possible. For instance, chang2021mobility ;
acuna2020108370 ; lopez2020 ; jia2020population  show that an in-depth risk-
assessment of spatial spread modeling could be made feasible by solving the
problem of capturing human behavioral change. Furthermore, characterizing
human mobility in the initial stage of the epidemic could reveal the most
probable propagation routes of the infection in the absence of interventions
Li2021 . Any scientific progress along this line of work would be valuable to
practitioners attempting to manage future transmissible airborne epidemics.

To date, many studies have focused on mid and large-scale population flows to
predict geographical distribution of COVID-19 wu2020nowcasting ;
jia2020population . At a domestic scale, the GPS technology in mobile phones
offers a fast and accurate tracking of individuals movement. While this tool
is extremely valuable, refined nationwide data is not always available nor
easily reachable in many countries. The above is particularly true for the
Central American region. Hence, in this study, we focused on developing and
proposing an evidence-based approach to model the spread risk of an airborne
disease outbreak within one developing country. To do that, we characterize
the spread progression of COVID-19 only at an early stage. As a focus study we
chose the country of origin of three of us: the Latin American country
Guatemala. Our modeling strategy can be implemented to any epidemic outbreaks
with similar transmissibility to the initial phase of the COVID-19 epidemic,
with mild or no interventions.

Solving the puzzle to establish a time line for the spatial spread of COVID-19
initial phase in Guatemala begins with the identification of regional epidemic
waves and the proper characterization of their relevant phenomenological
features. On the other hand, a spatial risk assessment of this disease relies
on identifying key trends present in the COVID-19 regional wave landscape and
finding out plausible spatial associations.

In this study we develop a systematic procedure to identify the initial wave
from daily incidence curves and to search for suitable analytical
approximations from a phenomenological epidemic model. We show that the
Richards Model, a straightforward generalization of the classical Logistic
Model richards1959flexible , notably approximates a wide range of growth
curves configurations, thus provides a convenient tool for wave description.
Our approach to wave classification is rooted in unconstrained ordination
abdi2010principal  of the Richards curves parameter estimates for every
geographic locale to look for similarity patterns in the waves properties as
embodied in these parameters. We further take advantage of available estimates
about the country inter-regional population flow to combine wave patterns with
spatial mobility data for spread risk assessment.

Our approach on the spatial association of regional outbreaks may prove to be
useful to refine epidemic risk strategies, particularly in the context of
countries devoid of resources for scientific research and a robust epidemic
response infrastructure. Finally, we also illustrate how the same information
that typically is used to provide COVID-19 trackers at National levels without
any further insights can be used to generate actionable risk information for
this and future epidemics in these countries.

The rest of this paper is organized as follows: In Section 2 we explain the
bases of logistic models and Richards equations. A brief description of
related work and the results of our study are exposed and discussed in Section
4. The final conclusions of our analysis are presented in Section 5. Finally,
section 3 describes the followed methodology regarding data treatment, wave
identification and classification, parameter estimation, PCA analysis, and
spatial risk assessment.

##  2 Background

###  2.1 Logistic Models and Disease Outbreaks

There is a vast literature showing that logistic phenomenological models are
capable of capturing the main features of growth profiles in single infectious
disease outbreaks ma2014estimating ; chowell201666 ; Pell2018 ; Brauer2019 .
Uncontrolled propagation that leads to exponential curves is the most common
assumption at the very beginning of an emerging epidemic. Accordingly, the
increase of the number of infected cases $c(t)$ over time is properly
described by the exponential growth model, which states that the per capita
rate of change of the total number of cases equals a constant $r$
chowell201666 ; viboud2016generalized . Explicitly, the exponential model is

$\frac{1}{c(t)}\frac{dc(t)}{dt}=r$ (1)

The $r$ parameter sets the timescale of the epidemic progress and is known as
the intrinsic growth rate. At the initial stage of an epidemic, $r$ can be
used to obtain an estimate of the basic reproduction number $R_{0}$
heffernan2005 , the number of secondary infection cases that are derived from
one contagious individual.

The central hypothesis in the exponential model lies on the assumption that,
at all times, the per capita contribution to the growth of the epidemic is
unaffected by the total number of cases. However, a variety of factors acting
as “density-dependent processes” can affect the evolution of the epidemic,
causing the outbreak to slow down and to reach a phase with decaying
transmission rate chowell201666 . The exponential model fails to follow this
behavior, which is instead faithfully reproduced by the logistic approach.
These models can be conceived as variations of the standard logistic growth
model, formulated in terms of two parameters: the generalized growth rate $r$
and the carrying capacity $K$. This second parameter accounts for the total
number of cumulative infections over the whole outbreak. The standard logistic
model heretofore abbreviated as LM is defined as ma2014estimating .

$\frac{1}{c(t)}\frac{dc(t)}{dt}=r\left(1-\frac{c(t)}{K}\right)$ (2)

Logistic type models are commonly assumed to be of empirical nature and do not
rely on a specific basis of mechanistic assumptions. They are instead rooted
on general assumptions about the shape of incidence curves, making use of a
few free parameters. However, it can be demonstrated that they approximately
reproduce the underlying mechanistic dynamics of the standard SIR model where
susceptible individuals are being converted into infected individuals
ma2014estimating . In fact, a little known result from theoretical ecology is
the work of Dennis (1978) dennis1978analytical  who showed that anytime a
population of interest (here the cumulative number of COVID-19 cases) grows on
a supply of a single limiting nutrient in a open system (here the susceptible
individuals), the resulting growth is logistic like. Furthermore, in closed
systems, the growth of such population in these circumstances is exactly
logistic. Briefly, Dennis’ dennis1978analytical  mathematical argument starts
by modeling the process of conversion of nutrient into population growth using
a system of two differential equations (one for the nutrient, one for the
population of interest) that very much resembles mathematically the standard
SIR model. He then shows that reducing this system of equations to a single
equation is possible and importantly, the resulting equation is logistic like
and in a special case, exactly logistic. These results from theoretical
ecology bring more relevance to phenomenological models of the logistic type,
as they provide equations with less parameters that bear a close relation to
epidemic dynamics. Importantly, the work of Dennis dennis1978analytical  and
ma2014estimating  show that the logistic-like model parameters are themselves
the product of parameters carrying information regarding the process of
conversion of susceptible into infected. The link of these models with
mechanistic processes is also likely to result in stable statistical fits and
reliable estimates in data analyses wu2020generalized .

###  2.2 Richards Model

The Richards Model richards1959flexible , originally introduced for ecological
population growth as the “theta-logistic model” sibly2005regulation , provides
a useful generalization of the classical LM. This model includes an extra
parameter, resulting in a more flexible growth function and thus, accounting
for variations in the transition time from the accelerating to the
decelerating phase of the epidemic curve.

In the Richards Model (RM), the differential equation for the cumulative cases
$c(t)$ is given by

$\frac{dc(t)}{dt}=rc(t)\left[1-\left(\frac{c(t)}{K}\right)^{a}\right]$ (3)

When $a=1$, RM is reduced to the standard logistic model. On the other hand,
when $a\ll 1$, the density-dependent term is lost.

The exact solution to Equation (3) is explicitly given by

$c(t)=K\left(1+ae^{-ar(t-t_{c})}\right)^{-1/a}$ (4)

where $t_{c}$ is the time where the second derivative of $c(t)$ vanishes and
therefore, corresponds to the inflection point from increasing to decreasing
rate of cumulative cases.

Although some of the parameters in the Richards empirical model may not have a
straightforward connection with epidemiological concepts, contrary to the
widely used compartmental SIR models, Equation (4) describes the cumulative
case data in a surprisingly accurate manner. Besides, some terms in the
Richards equation have been found to contain intrinsic relations to the ones
of the SIR model (e.g., the exponential term). A drawback of using the
Richards model approach is the curve over fitting tendency, for which it is
necessary to add some constraints in the modeling process WANG201212 . In any
case, here we show that using this simple model in the face of incomplete
epidemic information is useful to assess the temporal and spatial epidemic
risk of COVID-19 and other similar communicable diseases.

##  3 Methods

###  3.1 Wave Identification and Maximum Likelihood Estimation

Prior to the fitting process, we estimated the inflection points of the
observed epidemic curve, by implementing a computational algorithm to
determine the zeros of the variation of daily reports. For each department, we
applied a digital filter to the running average data to obtain a smoothed
version of the curve behavior. This filter was performed through a Savitzky-
Golay filter Savitzky1964 . This filter was also employed to identify the
location in time of the inflection point from increasing to decreasing rate of
cumulative cases, within a wave interval.

Once the beginning point and the duration of each wave was approximated, we
truncated the smoothed version of the cumulative growth curve to a time
window, defined by the results of the wave identification in each location.
Then, we started the fitting procedure by adjusting the parameters of the
Richards curve.

We chose the set of parameters $r,a,K$, and $t_{c}$ of Equation (4) that
maximized the likelihood function $L(\bm{\theta})$, where $\bm{\theta}$ is the
vector of parameters. We used the customary observation error model assumption
that, for the $ith$ department, the observed number of cumulative cases
$C_{i}$ followed a normal distribution. The mean of this distribution was
given by a deterministic process predicted by the Richards model, that is

$C_{i}\sim\mathcal{N}(\mu=K_{i}\left(1+a_{i}e^{-a_{i}r(t-t_{c_{i}})}\right)^{-1/a_{i}},\sigma^{2}),\qquad
i=1,2,\cdots,22$ (5)

Our likelihood optimization routine followed an embedded algorithm process.
Initially, we optimized the $\log$-likelihood $\log(L(\bm{\theta}))$ via the
Nelder-Mead numerical optimization method but this strategy routinely got
trapped in local optima. To solve this problem, we embedded a Nelder-Mead
search into each Monte Carlo step of a Simulated Annealing Algorithm (SAA).
This strategy not only proved to do better explorations of the parameter space
in the $(K,t_{c})$ dimensions but also systematically improved on Nelder-Mead
searches that got trapped in local minima. To get reliable search results, we
implemented random changes in $K$ and $t_{c}$ values in the neighboring states
set of the Markov chain. We performed 300 Monte Carlo steps to generate sample
states, and an annealing schedule of 50 different temperatures. This strategy
was sufficient to warrant convergence at a prescribed tolerance (see R code in
GitHub: https://github.com/jmponciano/Ponciano-JA-etal-COVID-Guatemala).

###  3.2 Dimensional Reduction and Wave Classification

The wave classification analysis was carried out with Principal Component
Analysis ($PCA$) aimed at reducing the $4$-dimensional parameter space into a
lower dimensional space that revealed patterns of similarity in the fitted
curves. This reduction was achieved by defining an orthogonal projection of
the wave parameters space onto a lower dimensional linear space, such that the
variance of the projected points is maximized johnson2014applied . The reduced
space obtained via $PCA$ was generated by the eigenvectors of the parameters
covariance matrix corresponding to the two largest eigenvalues. We call this
reduced space the $PCA$ parameters subspace. The original data matrix
consisted of a $22\times 4$ data matrix $X$ corresponding respectively to 22
regional waves, each one described by 4 parameters. Because $PCA$ is known to
be unreliable when the scale of the ordination variables varies widely, we
first standardized all parameter estimates to the same scale in the matrix
$X$. We then computed the empirical covariance matrix $S=\rm{Cov}(X)$ of the
wave parameters. Next, we considered a 2D subspace, which is determined by the
$S$ eigenvectors with the two largest eigenvalues. By construction, the
orthogonal projection of the data points into this subspace maximizes the
variance of the projected points. We used our own algorithm for the eigen-
decomposition of $S$, as well as the one for computing the scores
$z_{i}^{(1)}$ and $z_{i}^{(2)}$, corresponding to the projections of the data
into the first principal components $PC1$ and $PC2$, respectively. By
considering individually the one-dimensional spaces defined by $PC1$ and
$PC2$, and computing the euclidean separations among the corresponding
principal components scores, i.e., $d_{ij}=|z_{i}^{\alpha}-z_{j}^{\alpha}|$
for $\alpha=1,2$, we established different classifications of the regional
waves. This was achieved by resorting to the hierarchical clustering technique
within the Cluster package in R.

###  3.3 Mobility Matrix and Spatial Association

Classifying the waves by region allowed us to establish a timeline for the
spatial spread of COVID-19 virus in Guatemala. This timeline allowed in turn
exploring the spatial association among different regions by constructing a
mobility matrix from internal migratory flows. The migratory statistics at
departmental level were retrieved from the National Institute of Statistics
INE2018 . In this study, we considered the data related to the amount of
workers and students in every municipality that traveled to another department
daily, as reported by the last population census carried out in the country in
2018 INE2018 . As the information was reported with a municipality granularity
level, we aggregated the data to determine the quantities per department.

Following Lai et al.Lai2020 , we define for each department a risk of virus
importation, $\rho$, as the number of travelers received divided by the total
number of persons leaving the departments that have a high risk of
experiencing international infection during one day. We used this information
to complement the spatial analysis association of growth curves.

All the algorithms that are mentioned in this section were programmed in R.

##  4 Results and Discussion

The Richards model has been widely employed for the general description of
several multiple-wave epidemic events worldwide Shanafelt2018 ;
sanna2018spatial ; Pell2018 ; wu2020generalized ; Zuhairoh2020 ; WANG201212 ;
Macedo2021 . Also, it has been useful to estimate quantities of
epidemiological relevance, such as the initial growth rate and the inflection
point of an outbreak in a given location ma2014estimating ; hsieh2010epidemic
.

Some authors have applied the Principal Component Analysis technique to select
the main COVID-19 severity predictors Markovic2021  or to evaluate the impact
of containment strategies by public polices classification Chen2021 . In terms
of model parameter estimation, Agarwal et al. Agarwal2021  utilized a PCA to
obtain the fundamental parameters of a complex BHRP network model of the
second wave of COVID-19. Mahmoudi et al. Mahmoudi2021  performed a PCA to
classify high-risk countries according to their respective numbers of patients
and deaths. Another study presents a similar analysis of French cities, but
employs the PCA technique to categorize departments in accordance with
multidimensional variables such as epidemic indicators and predominant age.
Among other findings, the authors reported a relationship between department
mortality and the time mobility restrictions were adopted in that specific
location Gaudart2021 . To the best of our knowledge, no studies about COVID-19
wave classification have been published.

Here, we employed the official surveillance data of daily reported COVID-19
cases in Guatemala to carry a wave classification analysis based on both, a
fit of the Richards model and a subsequent principal component analyses of the
parameter estimates for each region. The only information needed was retrieved
from the COVID-19 web platform, managed by the national Ministry of Public
Health (MSPAS) MSPAS  and consisted of the time series of the daily registered
cases, starting from March 13th, 2020 onwards. The infected cases were
registered at the onset time of the disease symptoms. In this research, we
limited the wave identification analysis to a departmental level (Guatemala is
organized into 22 different Departments equivalent to states in the USA, and
within each department there are multiple “municipalities” akin to counties in
the USA). In total, there are 340 municipalities in the country. As the data
granularity is based on municipalities (the second-level administrative
divisions of the country as mentioned above) we pooled the data from all
municipalities per department in the country to obtain the cumulative values
for the 22 departments in Guatemala.

For the wave identification process, we first calculated the moving average of
the daily cases with a 14 day window. For model fitting using the Richards
equation, we made use of the time series of cumulative cases of each
department (see Methods). The time windows used for the analyses started 65
days after the first reported case in the country, in order to reduce the
noisy regime with strong fluctuations of the very early outbreak. For each
department, the wave identification was performed by fitting the first phase
of the cumulative infected to the Richards Model via maximum likelihood (see
Methods).

The results of fitting the Richards Model to the first wave of registered
COVID-19 cases are displayed in Figure 1. Table 1 contains a summary for the
maximum likelihood model parameters estimates for each department. In general,
the inflection points of the regional first waves related to the Richards
model are found from 57 to 174 days, starting from day 65 after the first
reported case in the country.

Figure 1: Richards model fit for the observed COVID-19 daily infected cases in the 22 departments of Guatemala Table 1: Maximum likelihood parameter estimates of the Richards model for COVID-19 cases in the 22 departments of Guatemala Department | Population INE2018  | $a$ | Growth rate, $r$ |  Carrying capacity, $K$ |  Inflection point, $t_{c}$
---|---|---|---|---|---
Alta Verapaz | 1215038 | 1.028 | 0.040 | 1685 | 142
Baja Verapaz | 299476 | 1.016 | 0.023 | 1258 | 159
Chimaltenango | 615776 | 1.010 | 0.044 | 1233 | 96
Chiquimula | 415063 | 1.008 | 0.060 | 941 | 106
El Progreso | 176632 | 1.062 | 0.099 | 400 | 57
Escuintla | 733181 | 1.005 | 0.084 | 3195 | 82
Guatemala | 3015081 | 1.011 | 0.057 | 27230 | 85
Huehuetenango | 1170669 | 1.006 | 0.030 | 2509 | 171
Izabal | 408688 | 1.006 | 0.041 | 2374 | 129
Jalapa | 342923 | 1.009 | 0.033 | 618 | 146
Jutiapa | 488395 | 1.029 | 0.053 | 495 | 82
Petén | 545600 | 1.004 | 0.036 | 2415 | 156
Quezaltenango | 799101 | 1.004 | 0.049 | 3048 | 120
Quiché | 949262 | 1.004 | 0.030 | 1783 | 169
Retalhuleu | 326828 | 1.003 | 0.039 | 1448 | 125
Sacatepequez | 330469 | 1.010 | 0.045 | 3374 | 102
San Marcos | 1032277 | 1.022 | 0.075 | 1087 | 62
Santa Rosa | 396607 | 1.020 | 0.057 | 697 | 78
Sololá | 421583 | 1.089 | 0.038 | 1598 | 174
Suchitepequez | 554695 | 1.004 | 0.049 | 2142 | 111
Totonicapán | 418569 | 1.014 | 0.031 | 723 | 141
Zacapa | 245374 | 1.013 | 0.032 | 1518 | 129

The resulting waves can be represented by four parameter estimates from the
Richards model. The points in a $4-$dimensional space defined by these
estimates were taken as the basis for the PCA and calculation of the wave
parameter space onto the principal subspace.

We found that the first principal component, $PC1$, explains $48\%$ of the
overall variability, while $PC2$ accounts for $29\%$ of it. This means that
$77\%$ of the variability is related just to $PC1$ and $PC2$; hence, we
reduced the variability analysis to the space defined by these two components.
The projection of the wave parameter data onto the $2D$ principal subspace is
displayed in Figure 2(a) along with the projections of the vector basis given
by $r,a,t_{c}$ and $K$. Examining the correlation matrix between the model
parameters and the principal components displayed in table 2 as well as Figure
2(a) it can be seen that the first principal component scores mainly capture
similarities in the peak position ($t_{c}$) and hence the time span of the
regional waves. On the other hand, from Figure 2(a) we can see that the second
principal component dimension gives more information about the carrying
capacity $k$ and the $a$ exponent parameter.

(a)

(b)

Figure 2: Obtained principal component space (a) and relation of the obtained
principal component space with the inverse virus importation risk (b)

Figure 3 illustrates the results of the classification of regional waves by
using euclidean distances, according to the $PCI$ score associated to each
wave. We found six different groups of waves that were labeled from $G1$ to
$G6$. In this figure, daily infected cases for each department have been
normalized by their maximum value.

Table 2: Correlations of Richards model parameters with principal components | PC1 | PC2 | PC3 | PC4
---|---|---|---|---
$a$ | 0.28 | 0.76 | \- 0.59 | 0.06
$r$ | 0.96 | 0.06 | 0.11 | -0.25
$K$ | 0.24 | -0.74 | -0.63 | -0.01
$t_{c}$ | -0.93 | 0.14 | -0.23 | -0.24

Wave classification results using Principal Component scores allowed us to
propose associations in the spatial distribution of the infection through the
following scenario:

First, the group of regional waves with the lowest mean inflection point
$t_{c}$ seems closely related to hot spots of virus importation (groups
$G1-G2$). Second, departments characterized by mid $t_{c}$ values have
presumably high connectivity to those hot spots (groups $G3-G4$). Also,
departments with significant population mobility (group $G3$) can act as
secondary sources for regional epidemic spread, boosting in that way the
epidemic course of rural and more isolated regions that belong to groups $G5$
and $G6$. Finally, once the outbreak gets triggered, the interplay of mobility
patterns and population density shapes the main traits of the regional growth
curves.

The proposed scenario above strongly suggests that geographical connectivity
is encoded in the global picture of the regional COVID-19 wave landscape.
Following this rationale, the local growth curves characterization can help to
single out the main elements that are necessary to depict a spread risk map of
the disease. Our investigation of the spatial associations in the context of
the COVID-19 epidemic and the estimation of a risk map from these also suggest
that this idea can be extended to model the spread risk of any airborne
infectious disease with similar epidemiology.

We compared our findings with available information on areas that are most
exposed to virus importation, along with regional mobility indicators. In
Guatemala, the main international land migratory flows occurs across the
southeastern border with Mexico and the southwestern border with El Salvador.
Those areas belong to the departments of San Marcos and Jutiapa, and are
considered high-risk regions. In addition, the highly populated city of
Guatemala, located in the Guatemala department, receives the highest volume of
international commercial flights. It is thus the most vulnerable region to
virus importation from an international source. These three high-risk
departments are sorted in groups $G1$ and $G2$, according to the $PCI$ scores
in Figure 3. Also, Figure 4 shows a mobility map that represents risk in three
disease epicentres; namely, the departments of Guatemala, San Marcos, and
Jutiapa. The figure also displays the information about the cumulative cases
per department at two different moments of the initial outbreak: the beginning
and the peak of the first wave.

Figures 2(b) and 5 explore an hypothetical relationship between the principal
component scores and the importation risk index $\rho$ defined above. Figure
5(a) shows a probable approximate relation between $PCII$ and the inverse of
the importation risk index. In order to check this statement, figure 5(b)
plots the second principal component $PCII$ as a function of $\rho^{-1}$. The
figure suggests a positive correlation between $PCII$ and the inverse of the
risk of importation index $\rho$. In the same figure, we show in different
symbols and colors the groups $G1$ to $G6$ formed by the $PCI$ classification.

When we tested the correlation between $PCII$ scores and $\rho^{-1}$ using a
linear model we found that the results are consistent with the hypothesis that
the importation risk might be used as a covariable to explain the $PCII$
scores. Within the linear model, the estimated internal migratory flow data is
correlated (Adjusted R-squared 0.3717, p-value 0.004676, non-parametric
bootstrap Manly2018  95% confidence interval $(0.08,0.70)$), through the
importation risk index, to the wave classification provided by the PCA
technique. Although the above results show a clear signal for a positive
correlation, the relationship between $PCII$ and $\rho^{-1}$ is rather non-
linear as observed in figures 5(b) and 6.

One can clearly identify two departments, namely El Progreso ($PCII=1.85$) and
Sololá ($PCII=2.74$), which have exceptionally high $PCII$ scores compared to
the rest of the data set. Those departments have an importation risk within
the 1st quartile of the distribution of values. El Progreso is also the less
populated department in the country (176632 people) which can explain why its
$PCII$ score, negatively correlated to the carrying capacity $K$ (see table
2), is high. On the other hand, the department of Sololá seems to be a special
case. It has a population within the 2nd quartile (421583 people) and we
cannot attempt such a straight explanation.

By relaxing the condition for a linear response and withdrawing the two above
departments considered as outliers, we fit the resulting dataset to a
Generalized Additive Model (GAM). The resulting value corresponding to
deviance explained by the GAM, estimated against the null model is $47.3\%$.
As a standard practice, this value is identified as a pseudo-$R^{2}$ for the
model. Figure 6 shows the GAM fit to the data, together with the one-standard
deviation region estimated by the model.

Leaving aside the disease epicenters that we take as primary hot spots, Figure
5 demonstrates that the Escuintla department, another member of group $G1$,
can be considered as a high-risk department for virus spread, given its $PCII$
score. Indeed, Table 1 shows that the carrying capacity of Escuintla’s first
wave was relatively high, unlike departments such as El Progreso, which also
belongs to $G1$ but has a lower carrying capacity (high $PCII$ score). We used
this criterion to withdrawn El Progreso from the virus spread high-risk group.
The Santa Rosa and Chiquimula departments, belonging to the $G2$ group,
experienced an early peak within the initial COVID outbreak. However, they are
characterized by intermediate $PCII$ and $\rho$ values, therefore they
presumably do not contribute dominantly in shaping the spatial spread pattern.
Therefore we’ll consider them as belonging to a third category of risk.

Seemingly, a key role in epidemic spread is played by the $G3$ group, as its
members have all similar wave patterns and relatively high importation risk
index. Those localities have important urban centers which exhibited high
community-level transmission during the pandemic and, additionally, present
high connectivity with their nearby departments. According to the PCA scores,
these departments can be considered as potential high-risk secondary cities in
the epidemic outbreak pathway.

Finally, Figure 5(b) strengthens the hypothesis that the departments from
groups $G5$ and $G6$ have, on average, the lowest connectivity to the primary
hot spots and therefore their epidemic wave peaks will appear delayed with
respect to the first regional waves. We summarize the principal points of this
discussion in the spread risk map of Figure 7.

Heretofore, our findings strongly support employing the wave patterns
landscape analysis as a guide to evaluate future epidemic spread risks in
countries like Guatemala, which lack a strong epidemiological response data
management infrastructure.

Figure 3: Classification of COVID-19 waves by the first principal component
scores in the 22 departments of Guatemala

(a)

(b)

(c)

Figure 4: Cumulative cases and mobility map for the departments of (a)
Guatemala, (b) San Marcos, and (c) Jutiapa. For each locality, the left and
right maps correspond, respectively, to the beginning and the peak of the
initial wave.

(a)

(b)

Figure 5: Relation between the second principal component $PCII$ and the
importation risk (a), and the inverse virus importation risk (b). Figure 6:
Generalized Additive Model fit  Figure 7: Disease spread risk map of Guatemala
based on COVID-19 initial wave outbreak

##  5 Conclusion

In this work, we focused on the initial phase of the COVID-19 epidemic in
Guatemala and used the regional growth curves characterization for spread risk
assessment. Using the phenomenological Richards Model applied to the initial
COVID-19 infection outbreak, we proposed an evidence-based risk map that
highlights the regions which play a key role in the temporal progression of
the disease. The proposed risk map is grounded on a spatial association
analysis, mainly provided by the classification of regional waves via the PCA
technique and complemented by human mobility data. The main contribution of
the present work is the use of past epidemic growth profiles to reveal not
only spatial similarities among different geographical regions but also
previously unnoticed connections among these. As a result, better predictions
of an epidemiological process in spatial networks is obtained.

The utility of the proposed risk map goes beyond the development of the
COVID-19 disease in Guatemala, and might be considered as a tool to guide
intervention plans in future outbreaks of airborne communicable diseases with
similar epidemiology. This is particularly relevant, considering that the
early conditions of the COVID-19 epidemic progression experienced in
Guatemala, where lax public health measures were adopted, are likely to be
reproduced in new epidemic outbreaks. Further studies are required to better
understand the coupled dynamics on networks. We believe that the elements
presented here consist of an important insight and first step on top of which
a thorough mechanistic modeling approach can be built.

Acknowledgments

Funding for this study was provided by the International Development Research
Centre (IDRC) and SICA-CSUCA for J.A.P., J.D.C., M.A and K.F.

## References

  *   * (1) World Health Organization: WHO Coronavirus (COVID-19) Dashboard (2022). https://covid19.who.int/
  * (2) Yang, Z., Zeng, Z., Wang, K., Wong, S.S., Liang, W., Zanin, M., Liu, P., Cao, X., Gao, Z., Mai, Z., Liang, J., Liu, X., Li, S., Li, Y., Ye, F., Guan, W., Yang, Y., Li, F., Luo, S., Xie, Y., Liu, B., Wang, Z., Zhang, S., Wang, Y., Zhong, N., He, J.: Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions. Journal of Thoracic Disease 12(3), 165–174 (2020). https://doi.org/10.21037/jtd.2020.02.64
  * (3) Yan, L., Zhang, H.-T., Goncalves, J., Xiao, Y., Wang, M., Guo, Y., Sun, C., Tang, X., Jing, L., Zhang, M., Huang, X., Xiao, Y., Cao, H., Chen, Y., Ren, T., Wang, F., Xiao, Y., Huang, S., Tan, X., Huang, N., Jiao, B., Cheng, C., Zhang, Y., Luo, A., Mombaerts, L., Jin, J., Cao, Z., Li, S., Xu, H., Yuan, Y.: An interpretable mortality prediction model for COVID-19 patients. Nature Machine Intelligence 2(5), 283–288 (2020). https://doi.org/10.1038/s42256-020-0180-7
  * (4) Chatterjee, K., Chatterjee, K., Kumar, A., Shankar, S.: Healthcare impact of COVID-19 epidemic in India: A stochastic mathematical model. Medical Journal Armed Forces India 76(2), 147–155 (2020). https://doi.org/%****␣covidGT-arXiv.bbl␣Line␣150␣****10.1016/j.mjafi.2020.03.022
  * (5) Bertozzi, A.L., Franco, E., Mohler, G., Short, M.B., Sledge, D.: The challenges of modeling and forecasting the spread of COVID-19. Proceedings of the National Academy of Sciences of the United States of America 117(29), 16732–16738 (2020). https://doi.org/10.1073/pnas.2006520117
  * (6) Shakeel, S.M., Kumar, N.S., Madalli, P.P., Srinivasaiah, R., Swamy, D.R.: Covid-19 prediction models: A systematic literature review. Osong Public Health and Research Perspectives 12(4), 215–229 (2021). https://doi.org/10.24171/J.PHRP.2021.0100
  * (7) Chang, S., Pierson, E., Koh, P.W., Gerardin, J., Redbird, B., Grusky, D., Leskovec, J.: Mobility network models of COVID-19 explain inequities and inform reopening. Nature 589(7840), 82–87 (2021). https://doi.org/10.1038/s41586-020-2923-3
  * (8) Acuña-Zegarra, M.A., Santana-Cibrian, M., Velasco-Hernandez, J.X.: Modeling behavioral change and COVID-19 containment in Mexico: A trade-off between lockdown and compliance. Mathematical Biosciences 325, 108370 (2020). https://doi.org/10.1016/j.mbs.2020.108370
  * (9) López, L., Rodó, X.: The end of social confinement and COVID-19 re-emergence risk. Nature Human Behaviour 4(7), 746–755 (2020). https://doi.org/10.1038/s41562-020-0908-8
  * (10) Jia, J.S., Lu, X., Yuan, Y., Xu, G., Jia, J., Christakis, N.A.: Population flow drives spatio-temporal distribution of COVID-19 in China. Nature 582(7812), 389–394 (2020). https://doi.org/10.1038/s41586-020-2284-y
  * (11) Li, Z., Li, H., Zhang, X., Zhao, C.: Estimation of human mobility patterns for forecasting the early spread of disease. Healthcare 9(9) (2021). https://doi.org/10.3390/healthcare9091224
  * (12) Wu, J.T., Leung, K., Leung, G.M.: Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. The Lancet 395(10225), 689–697 (2020). https://doi.org/10.1016/S0140-6736(20)30260-9
  * (13) Richards, F.J.: A Flexible Growth Function for Empirical Use. Journal of Experimental Botany 10(29), 290–300 (1959)
  * (14) Abdi, H., Williams, L.J.: Principal component analysis. Wiley Interdisciplinary Reviews: Computational Statistics 2(4), 433–459 (2010). https://doi.org/10.1002/wics.101
  * (15) Ma, J., Dushoff, J., Bolker, B.M., Earn, D.J.D.: Estimating Initial Epidemic Growth Rates. Bulletin of Mathematical Biology 76(1), 245–260 (2014). https://doi.org/10.1007/s11538-013-9918-2
  * (16) Chowell, G., Sattenspiel, L., Bansal, S., Viboud, C.: Mathematical models to characterize early epidemic growth: A review. Physics of Life Reviews 18, 66–97 (2016). https://doi.org/10.1016/j.plrev.2016.07.005
  * (17) Pell, B., Kuang, Y., Viboud, C., Chowell, G.: Using phenomenological models for forecasting the 2015 Ebola challenge. Epidemics 22, 62–70 (2018). https://doi.org/10.1016/j.epidem.2016.11.002
  * (18) Brauer, F., Castillo-Chavez, C., Feng, Z.: Mathematical Models in Epidemiology, pp. 382–385. Springer, New York (2019). https://doi.org/10.1007/978-1-4939-9828-9_17
  * (19) Viboud, C., Simonsen, L., Chowell, G.: A generalized-growth model to characterize the early ascending phase of infectious disease outbreaks. Epidemics 15, 27–37 (2016). https://doi.org/10.1016/j.epidem.2016.01.002
  * (20) Heffernan, J.M., Smith, R.J., Wahl, L.M.: Perspectives on the basic reproductive ratio. Journal of the Royal Society Interface 2(4), 281–293 (2005). https://doi.org/10.1098/rsif.2005.0042
  * (21) Dennis, B.: Analytical solution to an open-system model of population growth. Mathematical Biosciences 40(1-2), 167–169 (1978)
  * (22) Wu, K., Darcet, D., Wang, Q., Sornette, D.: Generalized logistic growth modeling of the COVID-19 outbreak: comparing the dynamics in the 29 provinces in China and in the rest of the world. Nonlinear Dynamics 101(3), 1561–1581 (2020). https://doi.org/10.1007/s11071-020-05862-6
  * (23) Sibly, R.M., Barker, D., Denham, M.C., Hone, J., Pagel, M.: On the regulation of populations of mammals, birds, fish, and insects. Science 309(5734), 607–610 (2005)
  * (24) Wang, X.S., Wu, J., Yang, Y.: Richards model revisited: Validation by and application to infection dynamics. Journal of Theoretical Biology 313, 12–19 (2012). https://doi.org/10.1016/j.jtbi.2012.07.024
  * (25) Savitzky, A., Golay, M.J.E.: Smoothing and Differentiation of Data by Simplified Least Squares Procedure. Analytical Chemistry 36(8), 1627–1639 (1964). https://doi.org/10.1021/ac60214a047
  * (26) Johnson, R.A., Wichern, D.W., et al.: Applied Multivariate Statistical Analysis vol. 6. Pearson London, UK:, ??? (2014)
  * (27) National Institute of Statistics: Cruce de variables (2019). http://redatam.censopoblacion.gt/bingtm/RpWebEngine.exe/Portal?BASE=CPVGT2018
  * (28) Lai, S., Bogoch, I.I., Ruktanonchai, N.W., Watts, A., Lu, X., Yang, W., Yu, H., Khan, K., Tatem, A.J.: Assessing spread risk of Wuhan novel coronavirus within and beyond China, January-April 2020: a travel network-based modelling study (2020). https://doi.org/10.1101/2020.02.04.20020479
  * (29) Shanafelt, D.W., Jones, G., Lima, M., Perrings, C., Chowell, G.: Forecasting the 2001 Foot-and-Mouth Disease Epidemic in the UK. EcoHealth 15(2), 338–347 (2018). https://doi.org/10.1007/s10393-017-1293-2
  * (30) Sanna, M., Wu, J., Zhu, Y., Yang, Z., Lu, J., Hsieh, Y.H.: Spatial and Temporal Characteristics of 2014 Dengue Outbreak in Guangdong, China. Nature Scientific Reports 8 (2018). https://doi.org/10.1038/s41598-018-19168-6
  * (31) Zuhairoh, F., Rosadi, D.: Real-time forecasting of the COVID-19 epidemic using the richards model in South Sulawesi, Indonesia. Indonesian Journal of Science and Technology 5(3), 456–462 (2020). https://doi.org/10.17509/ijost.v5i3.26139
  * (32) Macêdo, A.M.S., Brum, A.A., Duarte-Filho, G.C., Almeida, F.A.G., Ospina, R., Vasconcelos, G.L.: A Comparative Analysis between a SIRD Compartmental Model and the Richards Growth Model. Trends in Computational and Applied Mathematics 22(4), 545–557 (2021). https://doi.org/10.5540/tcam.2021.022.04.00545
  * (33) Hsieh, Y.H., Fisman, D.N., Wu, J.: On epidemic modeling in real time: An application to the 2009 Novel A (H1N1) influenza outbreak in Canada. BMC Research Notes 3, 2–9 (2010). https://doi.org/10.1186/1756-0500-3-283
  * (34) Markovic, S., Rodic, A., Salom, I., Milicevic, O., Djordjevic, M., Djordjevic, M.: COVID-19 severity determinants inferred through ecological and epidemiological modeling. One Health 13(August) (2021). https://doi.org/10.1016/j.onehlt.2021.100355
  * (35) Chen, J., Gong, F., Xiang, S., Yu, T.: Application of principal component analysis in evaluation of epidemic situation policy implementation. Journal of Physics: Conference Series 1903 (2021). https://doi.org/10.1088/1742-6596/1903/1/012056
  * (36) Agarwal, D.P., Ch, M.A., Kharate, S., Fantin, E., Raj, I., Balamuralitharan, S.: Parameter estimation of covid-19 second wave b-h-r-p transmission model by using principle component analysis. Annals of R.S.C.B. 25(5), 446–457 (2021)
  * (37) Mahmoudi, M.R., Heydari, M.H., Qasem, S.N., Mosavi, A., Band, S.S.: Principal component analysis to study the relations between the spread rates of COVID-19 in high risks countries. Alexandria Engineering Journal 60, 457–464 (2021). https://doi.org/10.1016/j.aej.2020.09.013
  * (38) Gaudart, J., Landier, J., Huiart, L., Legendre, E., Lehot, L., Bendiane, M.K., Chiche, L., Petitjean, A., Mosnier, E., Kirakoya-Samadoulougou, F., Demongeot, J., Piarroux, R., Rebaudet, S.: Factors associated with the spatial heterogeneity of the first wave of COVID-19 in France: a nationwide geo-epidemiological study. The Lancet Public Health 6(4), 222–231 (2021). https://doi.org/10.1016/S2468-2667(21)00006-2
  * (39) Ministerio de Salud Pública y Asistencia Social: Situación de COVID-19 en Guatemala (2022). https://tablerocovid.mspas.gob.gt/
  * (40) Manly, B.F.J.: Randomization, Bootstrap and Monte Carlo Methods in Biology. Chapman and Hall/CRC, ??? (2018). https://doi.org/10.1201/9781315273075


# Constant-time connectivity tests

Philipp Klaus Krause

###### Abstract

We present implementations of constant-time algorithms for connectivity tests
and related problems. Some are implementations of slightly improved variants
of previously known algorithms; for other problems we present new algorithms
that have substantially better runtime than previously known algorithms
(estimates of the distance to and tolerant testers for connectivity, $2$-edge-
connectivity, $3$-edge-connectivity, eulerianity).

##  1 Introduction

Property testing is concerned with extremely fast (constant-time or other
sublinear) algorithms for approximate decision-making. While the runtime of
constant-time algorithms does not depend on the size of the input, the runtime
depending on parameters, such as the average degree of the input graph, or on
the maximum allowed error of the output are still the subject of research.

We implemented constant-time graph algorithms for testing connectivity,
$2$-edge connectivity and $3$-edge connectivity for sparse graphs. We
implemented an estimate of the number of connected components, estimates of
the distance to connectivity, $2$-edge-connectivity and eulerianity, tolerant
testers for connectivity, $2$-edge-connectivity and eulerianity. we also
present an algorithm for estimating the distance to $3$-edge-connectivity and
a tolerant tester for $3$-edge-connectivity that have not been implemented
yet. For connectivity eulerianity our approach has better worst-case runtime
than previous approaches (though the same expected runtime). For $2$-edge-
connectivity and $3$-edge-connectivity, our approaches have better runtime
than the known algorithm for general $k$-edge-connectivity. These advantages
carry over to the tolerant testers.

Full, compileable C source of our implementation can be found as free software
[18] at http://zshg.sourceforge.net.

##  2 Preliminaries

Let $G=(V,E)$ be a graph. For a set of nodes $U\subseteq V$, the degree $d(U)$
is the number of edges with exactly one endpoint in $U$. For a single node
$v\in V$ we also define the degree $d(v)\mathrel{\mathop{:}}=d(\\{v\\})$. The
number $\delta(G)\mathrel{\mathop{:}}=\min\\{d(v)\ |\ v\in V\\}$ is the
_minimum degree_ of $G$, $\Delta(G)\mathrel{\mathop{:}}=\max\\{d(v)\ |\ v\in
V\\}$ is the _maximum degree_ of $G$. The number

$d(G)\mathrel{\mathop{:}}=\frac{1}{|V|}\sum_{v\in V}d(v)=2\frac{|E|}{|V|}$

is the _average degree_ of $G$.

In the _bounded-degree model_ , a graph $G=(V,E),\Delta\geq\Delta(G)$ is
_$\epsilon$ -far_ from having a property $\mathcal{P}$, if it cannot be
transformed into a graph
$G^{\prime}\in\mathcal{P},\Delta(G^{\prime})\leq\Delta$ by at most
$\epsilon\Delta|V|$ edge modifications. Complexities are given as functions of
the maximum degree $\Delta$.

In the _sparse graph model_ (also called _unbounded-degree model_), a graph
$G=(V,E)$ is _$\epsilon$ -far_ from having a property $\mathcal{P}$, if it
cannot be transformed into a graph $G^{\prime}\in\mathcal{P}$ by at most
$\epsilon|E|$ edge modifications (i.e. edge insertions and edge deletions).
Complexities are usually given as functions of the average degree $d$. When
estimating distance to a property, often the distance is given in terms of
edge modifications relative to the number of nodes ($\delta$):
$\delta|V|=\epsilon|E|$.

The _incidence-lists model_ consists of the bounded-degree model and the
unbounded-degree model. It corresponds to an implementation of the graph data
structure using adjacency lists or incidence lists.

In the incidence lists model the available queries are _degree_ and _$i$ -th
neighbour_.

Other models include the _dense graph model_ and the _combined model_ (also
called the _general graph model_). The dense graph model corresponds to an
implementation of the graph data structure using an adjacency matrix. The
available query is _adjacency_. The combined model allows all queries from the
incidence lists and the dense graph models.

All models allow uniform random sampling of nodes.

A graph is _connected_ if it has at least one node and for any two nodes $u,v$
in the graph, there is a path from $u$ to $v$. A _connected component_ of a
graph is a maximal connected subgraph. A graph is called $k$_-edge-connected_
, if any graph obtained from it by removing less than $k$ edges is connected.
An edge in a graph is a _bridge_ , if all paths connecting its endpoints
contain the edge. A $k$-class in a graph $g$ is a $k$-edge connected subset of
nodes of a graph $G$ (note that for $k>2$ a $k$-class in $G$ is not
necessarily a $k$-edge-connected subgraph of $G$). A $k$-set in a graph is a
$k$-edge-connected subgraph that is connected to the rest of the graph by less
than $k$ edges. Every node in a graph is contained in at least one $k$-set.
The $1$-sets are the connected components. A minimal $k$-set is called a
$k$-class-leaf (the terminology is inspired by the use of the word “leaf” for
$2$-class-leaves [3] and by the creation of a auxiliary graph, that has one
node per $k$-class in the original graph. For $k=2$, this auxiliary graph is
the bridge tree. The $k$-class-leaves correspond to nodes of degree $1$ in the
auxiliary graph.

A set of nodes $U\subseteq V$ is _$\ell$ -extreme_, if $d(U)=\ell$ and each
subset of $U$ has degree larger than $\ell$. Every $\ell$-extreme set is an
$\ell^{\prime}$-class for some $\ell^{\prime}$ and any two extreme sets in a
graph are either disjoint or one is a subset of the other [13]. For
$U,W\subseteq V$ we use the notation $U\sqsubset W$ to denote that $U\subseteq
W$ and $U$ is extreme and there is no extreme set $U^{\prime}$, such that
$U\subseteq U^{\prime}\subseteq W$.

A graph is _eulerian_ if it contains a cycle that traverses each edge of the
graph exactly once.

##  3 Related Work

Property | Model | Complexity
---|---|---
Connectivity [8] | $\Delta$ | $O\left(\frac{\log^{2}\left(\frac{1}{\epsilon\Delta}\right)}{\epsilon}\right)$
Connectivity [15] | $d$ | $O\left(\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d}\right)$
Connectivity [14] | $d$ | $O\left(\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d^{2}}\right)$
Connectivity | $d$ | $O\left(-\log(1-p)\cdot\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d^{2}}\right)$
$2$-edge-connectivity [8] | $\Delta$ | $O\left(\frac{\log^{2}\left(\frac{1}{\epsilon\Delta}\right)}{\epsilon}\right)$
$2$-edge-connectivity | $d$ | $O\left(-\log(1-p)\cdot\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d^{2}}\right)$
$3$-edge-connectivity [8] | $\Delta$ | $O\left(\frac{\log\left(\frac{1}{\epsilon\Delta}\right)}{\epsilon^{2}}\right)$
$3$-edge-connectivity | $d$ | $O\left(-\log(1-p)\cdot\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{3}d^{3}}\right)$
$k$-edge-connectivity [8] | $\Delta$ | $O\left(\frac{k^{3}\log\left(\frac{1}{\epsilon\Delta}\right)}{\epsilon^{3-\frac{2}{k}}\Delta^{2-\frac{2}{k}}}\right)$
$k$-edge-connectivity [14] | $d$ | $\tilde{O}\left(\frac{k^{4}}{\epsilon^{4}d^{4}}\right)$
Eulerianity [8] | $\Delta$ | $O\left(\frac{\log^{2}\left(\frac{1}{\epsilon\Delta}\right)}{\epsilon}\right)$
Eulerianity [14] | $d$ | $O\left(\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d^{2}}\right)$
Eulerianity | $d$ | $O\left(-\log(1-p)\cdot\frac{\log\left(\frac{1}{\epsilon d}\right)}{\epsilon^{2}d^{2}}\right)$
Figure 1: Comparison of connectivity tests Property | Model | Complexity
---|---|---
Connectivity [5] | $d$ | $O\left(\frac{d\log\left(\frac{d}{\delta}\right)}{\delta^{2}}\right)$
Connectivity [12] | $d$ | $O\left(\frac{1}{\epsilon^{4}d^{4}}\right)=O\left(\frac{1}{\delta^{4}}\right)$
Connectivity [1] | $d$ |  $O\left(\frac{1}{\delta^{4}}\right)$, exp.: $O\left(\frac{\log\left(\frac{1}{\delta}\right)}{\delta^{2}}\right)$
Connectivity | $d$ | $O\left(\frac{\log\left(\frac{1}{\delta}\right)}{(1-p)\delta^{2}}\right)$
$2$-edge-connectivity | $d$ | $O\left(\frac{\log\left(\frac{1}{\delta}\right)}{(1-p)\delta^{2}}\right)$
$3$-edge-connectivity | $d$ | $O\left(\frac{\log\left(\frac{1}{\delta}\right)}{(1-p)\delta^{3}}\right)$
$k$-edge-connectivity [12] | $d$ | $O\left(\frac{k^{6}\log\left(\frac{k}{\epsilon d}\right)}{\epsilon^{6}d^{6}}\right)=O\left(\frac{k^{6}\log\left(\frac{k}{\delta}\right)}{\delta^{6}}\right)$
Eulerianity [12] | $d$ | $O\left(\frac{1}{\epsilon^{4}d^{4}}\right)=O\left(\frac{1}{\delta^{4}}\right)$
Eulerianity | $d$ | $O\left(\frac{\log\left(\frac{1}{\delta}\right)}{(1-p)\delta^{2}}\right)$
Figure 2: Comparison of estimates of distance to connectivity

For connectivity, there is a relatively simple algorithm with small constants
[8] in the bounded-degree model, that also works well in the sparse graph
model [15], where it has slightly higher complexity. With a more careful
analysis, one can achieve a slightly better complexity bound [14]. Based on
the existing approaches, we came up with a connectivity test that uses fewer
queries (but has the same asymptotic complexity).

The algorithm for connectivity testing can be generalized to $k$-edge-
connectivity testing [8]. These generalizations are more involved than the
algorithm for connectivity, but they still seem implementable, which also
holds for their sparse graph versions [14]. A graph is eulerian when it is
connected and has no nodes of odd degree [7]. Thus the algorithms for
connectivity can easily be modified into algorithms testing for Eulerianity.

Tolerant testers can be created from algorithms that estimate the distance to
having the property (as we do in Section 6). There are various approaches to
estimating the distance to connectivity. A recent one [1] has an expected
runtime of
$O\left(\frac{\log\left(\frac{1}{\delta}\right)}{\delta^{2}}\right)$, though
the worst-case complexity is substantially higher at
$O\left(\frac{1}{\delta^{4}}\right)$. An earlier approach [12] can estimate
the distance to $k$-edge-connectivity for any $k$. For $k=2$ and $k=3$ we can
estimate the distance with lower asymptotic complexity than the general
approach.

There was an attempt to construct tolerant testers for $k$-edge-connectivity
via local reconstructors [4]. The claimed complexity was substantially higher
than previous approaches and exponential in $k$; there are multiple flaws in
the proof, that cannot be fixed easily.

##  4 Connectivity Test

We implemented a connectivity (and $2$\- and $3$-edge connectivity) test for
the sparse graph model. It is very similar to existing connectivity test
algorithms [8, 15, 14] are very similar, which also carries over to the
proofs. Our approach (Figure 3) follows them closely, with a few minor changes
to reduce the query complexity by a constant factor. Instead of having a fixed
probability of $\frac{2}{3}$ for rejecting graphs that are $\epsilon$-far from
being connected, we allow an extra parameter $p$; compared to the approach of
running the algorithm repeatedly, this allows for reduced runtime and number
of queries when the acceptable probability $1-p$ of an unconnected graph that
is $\epsilon$-far from being connected being considered connected is not a
power of $\frac{1}{3}$. The implementations of the algorithms rand_index(),
which uniformly selects a random node, and k_set($s$, $2^{i}$), which return
the number of nodes in the minimal $k$-set containing $s$ connected to the
rest of the graph by less than $k$ edges are not shown.

For $k=1$, the algorithm is a connectivity test (k_set() can be implemented as
a depth-first search). For $k=2$ and $k=3$, k_set() can be implemented using
variants of algorithms for identifying 2-class-leaves and 3-class leaves [8].

⬇

bool zshg_k($\epsilon$, $p$, $d$, $n$)

{

$q=-\log_{e}(1-p)$;

if(!$n$)

return(false);

if($d<1$)

return($n$ == 1);

$\ell=\log_{2}\left(\frac{8\lceil\frac{k}{2}\rceil}{\epsilon d}\right)$;

if($n\leq\frac{8q\lceil\frac{k}{2}\rceil\ell}{\epsilon d}$)

return(k_set(0, $n$) == $n$);

for($i$ = 1; $i\leq\ell$; i++)

{

$m_{i}=\frac{8q\lceil\frac{k}{2}\rceil\ell}{2^{i}\epsilon d}$;

for($c$ = 0; $c\leq m_{i}$; c++)

{

$s$ = rand_index(n);

if(k_set($s$, $2^{i}$) < $2^{i}$)

return(false);

}

}

return(true);

}

Figure 3: $k$-edge-connectivity test algorithm in C-like pseudocode (the
implementation resides in zshg_impl() in zshg.c)

######  Lemma 1.

The algorithm in Figure 3 accepts all $k$-edge-connected graphs and rejects
graphs that are $\epsilon$-far from being $k$-edge-connected with probability
at least $p$.

###### Proof.

Clearly, the algorithm never rejects a $k$-edge-connected graph. If a graph is
$\epsilon$-far from being $k$-edge-connected, it has more than $\frac{\epsilon
m}{\lceil\frac{k}{2}\rceil}$ $k$-class-leaves (otherwise one could add
$\epsilon mk$ edges to connect all $k$-class-leaves by adding
$\lceil\frac{k}{2}\rceil$ cycles through them, resulting in a $k$-edge-
connected graph). Then the graph has at least $\frac{\epsilon
m}{2\lceil\frac{k}{2}\rceil}$ $k$-class-leaves consisting of less than
$\frac{4\lceil\frac{k}{2}\rceil}{\epsilon d}$ nodes each. Let $B_{i}$ be the
set of all such $k$-class-leaves of size at least $2^{i-1}$, but less than
$2^{i}$. Let
$\ell\mathrel{\mathop{:}}=\log_{2}(\frac{8\lceil\frac{k}{2}\rceil}{\epsilon
d})$. Since $\sum_{i=1}^{\ell}|B_{i}|\geq\frac{\epsilon
m}{2\lceil\frac{k}{2}\rceil}$, there exists an $i,1\leq i\leq\ell$, such that
$|B_{i}|\geq\frac{\epsilon m}{2\lceil\frac{k}{2}\rceil\ell}$. The number of
nodes residing in $k$-class-leaves of size at least $2^{i-1}$, but less than
$2^{i}$ is at least $2^{i-1}|B_{i}|\geq\frac{2^{i-2}\epsilon
m}{\lceil\frac{k}{2}\rceil\ell}$. The probability of selecting such a node by
uniform sampling is $\frac{2^{i-1}|B_{i}|}{n}\geq\frac{2^{i-2}\epsilon
m}{n\lceil\frac{k}{2}\rceil\ell}=\frac{2^{i-3}\epsilon
d}{\lceil\frac{k}{2}\rceil\ell}$. The probability of not selecting any such
node in iteration $i$ of the algorithm bounds the probability
$p_{\textrm{fail}}$ of not rejecting an unconnected graph that is
$\epsilon$-far from being connected. Using $1+x\leq e^{x}$, which follows from
Bernoulli’s inequality [17, 2], we get

$\displaystyle p_{\textrm{fail}}\leq\left(1-\frac{2^{i-3}\epsilon
d}{\left\lceil\frac{k}{2}\right\rceil\ell}\right)^{m_{i}}\leq
e^{-\left(\frac{2^{i-3}\epsilon
d}{\left\lceil\frac{k}{2}\right\rceil\ell}m_{i}\right)}=e^{-\left(\frac{2^{i-3}\epsilon
d}{k\ell}\frac{8q\left\lceil\frac{k}{2}\right\rceil\ell}{2^{i}\epsilon
d}\right)}=$ $\displaystyle=e^{-q}=e^{\log_{e}(1-p)}=1-p.$

∎

######  Lemma 2.

For $k=1,2$ the query complexity of the algorithm in Figure 3 is
$O\left(-\frac{\log(1-p)\log\left(\frac{1}{\epsilon
d}\right)}{\epsilon^{2}d^{2}}\right)$. For $k=3$ it is
$O\left(-\frac{\log(1-p)\log\left(\frac{1}{\epsilon
d}\right)}{\epsilon^{3}d^{3}}\right)$.

###### Proof.

In each iteration of the algorithm, $m_{i}$ nodes are selected. For $k=1$, for
each selected node, in k_set() a depth-first-search for up to $2^{i}$ nodes is
done (dfs() in zshg.c), resulting in less than $2^{i}\cdot 2^{i}$ queries for
each search. The query complexity is thus less than

$\displaystyle\sum_{i=1}^{\ell}m_{i}\cdot 2^{i}\cdot
2^{i}=\sum_{i=1}^{\ell}\frac{8q\ell}{2^{i}\epsilon d}\cdot 2^{i}\cdot
2^{i}\leq\frac{8q\ell}{\epsilon
d}2^{\ell+1}=\frac{8q\log_{2}\left(\frac{8}{\epsilon d}\right)}{\epsilon
d}2\frac{8}{\epsilon d}=$
$\displaystyle=128\frac{q\log_{2}\left(\frac{8}{\epsilon
d}\right)}{\epsilon^{2}d^{2}}=128\frac{-\log_{e}\left(1-p\right)\cdot\log_{2}\left(\frac{8}{\epsilon
d}\right)}{\epsilon^{2}d^{2}}\in
O\left(-\log\left(1-p\right)\cdot\frac{\log\left(\frac{1}{\epsilon
d}\right)}{\epsilon^{2}d^{2}}\right).$

For $k=2$, zshg uses a variant of an algorithm for 2-class-leaves [8] (comp2()
in zshg2.c). This algorithm essentially does two depth-first-searches for up
to $2^{i}$ nodes, with the second search somewhat restricted, for a total of
at most $(2^{i})^{2}$ queries. This results in a call to 2_set needing up to
$4$ times as many queries as 1_set, and thus a total factor of about $4$ in
the number of queries done by zshg_2 vs. zshg_1. For $k=3$, we use a variant
of an algorithm for 3-class-leaves [8]. This algorithm essentially does one
depth-first-search for up to $2^{i}$ nodes followed by an invocation of
2_set() on each node discovered. We get an additional factor of about
$\frac{16}{\epsilon d}$ in the number of queries by zshg_3 vs. zshg_1. ∎

From our connectivity test, we obtain a test for eulerianity based on the
following:

######  Lemma 3.

Let $G=(V,E)$ be a graph that is $\epsilon$-far from eulerianity. Then, it is
$\frac{\epsilon}{2}$-far from connectivity or it has more than $\epsilon dn$
nodes of odd degree.

###### Proof.

Assume that the graph is $\frac{\epsilon}{2}$-close to connectivity (i.e. it
has at most $\frac{\epsilon}{2}dn$ connected components) and at most
$\frac{\epsilon}{2}dn$ nodes of odd degree. We show that we can make the graph
Eulerian by adding $\epsilon dn$ edges, a contradiction. we first add
$\frac{\epsilon}{2}dn$ edges connecting nodes of odd degree to obtain a graph
in which all nodes have even degree. We then add $\frac{\epsilon}{2}dn$ edges,
choosing one node out of each connected component in the original graph and
adding a cycle through these. The resulting graph is connected and has only
nodes of even degree and is thus eulerian [7, 11, 9]. ∎

######  Theorem 1.

There are testers for $1$\- and $2$-edge-connectivity and for eulerianity with
complexity $O\left(-\frac{\log(1-p)\log\left(\frac{1}{\epsilon
d}\right)}{\epsilon^{2}d^{2}}\right)$, a tester for $3$-edge-connectivity with
complexity $O\left(-\frac{\log(1-p)\log\left(\frac{1}{\epsilon
d}\right)}{\epsilon^{3}d^{3}}\right)$..

###### Proof.

This follows from Lemmata 1, 2, 3. ∎

##  5 Distance to Connectivity Estimate

Our algorithm for estimating the number of connected components is based on an
earlier algorithm [1], which has a worst-case query complexity of
$O(\delta^{-4})$ and an expected query complexity of
$O(\delta^{-2}\log(\delta^{-1}))$. We made some modifications to reduce the
worst-case query complexity to $O(\delta^{-2}\log(\delta^{-1}))$ and allow the
probability $p$ of the error being outside the error bound to be specified
instead of always using a fixed probability of $\frac{3}{4}$.

rand_range($i$) returns a random integer in the range $1,\ldots,i$. Each
number $j$ is returned with probability $j^{-2}-(j+1)^{-2}$ for $j<i$, while
$i$ is returned with probability $i^{-2}$. The function 1_set() can be
implemented as above in Section 4; the additional parameter &$\ell$ is used to
limit the total number of queries made: $\ell$ is decremented each time a
query is made; when $\ell$ reaches $0$, no further queries are made.

⬇

float 1_sets($\delta$, $p$, $n$)

{

$r=\frac{2}{(1-p)\delta^{2}}$;

$a=0$;

$\ell=r\left(\log_{e}\left(\frac{2}{\delta}\right)+\frac{9}{2}\right)$;

for($i$ = 0; $i<r$; $i$++)

{

$s$ = rand_index($n$);

$x$ = rand_range($\frac{2}{\delta}$);

$b$ = 1_set($s$, $x+1$, &$\ell$);

if($b\leq x$)

$a$ += $b$;

}

return($\frac{an}{r}+\frac{\delta n}{4}$);

}

Figure 4: $1$-set counting algorithm in C-like pseudocode (the implementation
resides in zshg_components() in zshg_c.c)

######  Lemma 4.

For $\delta\leq 1$, with probability at least $p$, the return value of the
algorithm in Figure 4 is the number of connected components in the graph up to
an error of at most $\delta n$.

###### Proof.

Let $X_{i}$ be the value of $x$ in iteration $i$. Let $c$ be the number of
connected components. Let $c^{*}$ be the number of connected components of
size at most $\frac{2}{\delta}$. Let $C$ be the return value of the algorithm
in Figure 4. Let $\hat{c}$ be the return value when ignoring the query bound
$\ell$. Let $B_{i}$ be the value of $b$ in iteration $i$ if $b<X_{i}$ and $0$
otherwise. Let $C_{s}$ be the size of the connected components containing the
node $s$.

$\displaystyle E(B_{i})$ $\displaystyle=\frac{1}{n}\sum_{s\in
V}C_{s}P(X_{i}\geq C_{s})=\frac{1}{n}\sum_{\begin{subarray}{c}s\in V\\\
C_{s}\leq\frac{2}{\delta}\end{subarray}}\frac{1}{C_{s}}=\frac{c^{*}}{n}\leq
1.$ $\displaystyle E(B_{i}^{2})$ $\displaystyle=\frac{1}{n}\sum_{s\in
V}C_{s}^{2}P(X_{i}\geq C_{s})=\frac{1}{n}\sum_{\begin{subarray}{c}s\in V\\\
C_{s}\leq\frac{2}{\delta}\end{subarray}}1\leq 1.$ $\displaystyle Var(B_{i})$
$\displaystyle=E(B_{i}^{2})-E(B_{i})^{2}\leq E(B_{i}^{2})\leq 1.$
$\displaystyle\hat{c}$ $\displaystyle=\frac{a}{n}\sum_{i=0}^{r-1}+\frac{\delta
n}{4}.$ $\displaystyle E(\hat{c})$
$\displaystyle=\frac{rE(B_{i})n}{r}+\frac{\delta n}{4}=c^{*}+\frac{\delta
n}{4}.$ $\displaystyle Var(\hat{c})$
$\displaystyle=\frac{n^{2}}{r}Var(B_{i})\leq\frac{n^{2}}{r}.$

The number of connected components bigger than $\frac{2}{\delta}$ is at most
$\frac{\delta n}{2}$. This allows us to bound the probability of the estimate
being far off [19] (still ignoring the query bound $\ell$).

$\displaystyle c-\frac{\delta n}{4}<c^{*}+\frac{\delta
n}{4}=E(\hat{c})=c^{*}+\frac{\delta n}{4}\leq c+\frac{\delta n}{4}.$
$\displaystyle P\left(\left|\hat{c}-E(\hat{c})\right|>\frac{3\delta
n}{4}\right)<\frac{Var(\hat{c})}{\left(\frac{3\delta
n}{4}\right)^{2}}\leq\frac{8}{9}(1-p).$

Let $Q_{i}$ be the number of queries made by a call to 1_set($s$, $X_{i}$) for
a uniformly chosen node $s$ (i. e. we consider the behaviour of the algorithm
as if the limit $\ell$ wasn’t there). The procedure 1_set can be implemented
as a depth-first-search, resulting in at most $\frac{x(x+1)}{2}$ queries for
1_set($s$, $x$).

$\displaystyle E(Q_{i})$
$\displaystyle\leq\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{1}{j^{2}}-\frac{1}{(j+1)^{2}}\right)\frac{j(j+1)}{2}+\frac{1}{(\frac{2}{\delta})^{2}}\frac{\frac{2}{\delta}(\frac{2}{\delta}+1)}{2}=$
$\displaystyle=\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{2j+1}{2j(j+1)}\right)+\frac{2+\delta}{4}\leq\log_{e}\left(\frac{2}{\delta}\right)+\frac{1}{2}.$
$\displaystyle E(Q_{i}^{2})$
$\displaystyle\leq\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{1}{j^{2}}-\frac{1}{(j+1)^{2}}\right)\left(\frac{j(j+1)}{2}\right)^{2}+\frac{1}{(\frac{2}{\delta})^{2}}\left(\frac{\frac{2}{\delta}(\frac{2}{\delta}+1)}{2}\right)^{2}=$
$\displaystyle=\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{(j+1)^{2}}{4}-\frac{j^{2}}{4}\right)+\frac{(\frac{2}{\delta}+1)^{2}}{4}=$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{\frac{2}{\delta}-1}j+\frac{1}{4}\left(\frac{2}{\delta}-1\right)+\frac{(\frac{2}{\delta}+1)^{2}}{4}=\frac{2}{\delta^{2}}+\frac{1}{\delta}.$
$\displaystyle Var(Q_{i})$ $\displaystyle=E(Q_{i}^{2})-E(Q_{i})^{2}\leq
E(Q_{i}^{2})\leq\frac{2}{\delta^{2}}+\frac{1}{\delta}.$

This allows us to bound the probability of hitting the query bound $\ell$
[19].

$\displaystyle P\left(\sum_{i=0}^{r-1}Q_{i}>\ell\right)\leq
P\left(\sum_{i=0}^{r-1}Q_{i}>r\left(E(Q_{0})+4\right)\right)<\frac{Var(Q_{0})}{16r}<\frac{1}{9}(1-p).$

We can thus bound the probability of the algorithm returning a result outside
the error bounds.

$\displaystyle P\left(|C-c|>\delta n\right)\leq
P\left(\sum_{i=0}^{r-1}Q_{i}>\ell\textrm{ or
}\left|\hat{c}-E(\hat{c})\right|>\frac{3\delta n}{4}\right)\leq$
$\displaystyle\leq
P\left(\sum_{i=0}^{r-1}Q_{i}>\ell\right)+P\left(\left|\hat{c}-E(\hat{c})\right|>\frac{3\delta
n}{4}\right)<(1-p).$

∎

######  Lemma 5.

The algorithm in Figure 4 has query complexity
$O(\frac{1}{(1-p)\delta^{2}}\log_{e}(\frac{1}{\delta}))$.

###### Proof.

This follows directly from the choice of $\ell$. ∎

While our algorithm for estimating the distance to connectivity did directly
reuse the function 1_set from the connectivity test, this doesn’t work for
general $k$-edge-connectivity, since for $k\geq 2$, not every $k$-set is a
$k$-class-leaf. For $k=2$ we use:

######  Lemma 6 (Eswaran and Tarjan [6]).

Let $G=(V,E)$ be a graph. Let $c_{1}$ be the number of 2-class-leaves that are
1-class-leaves in $G$. Let $c_{2}$ be the number of 2-class-leaves that are
not 1-class-leaves in $G$. If $|V|>2$ and $c_{1}+c_{2}>1$, the minimal number
of edge modifications necessary to make $G$ $2$-edge-connected is
$\left\lceil\frac{c_{2}}{2}\right\rceil+c_{1}$.

Using classic terminology, $c_{1}$ is the number of connected components in
$G$ that are $2$-edge-connected and $c_{2}$ is the number of $2$-edge-
connected components in $G$ that are connected to the rest of the graph by a
single bridge. To estimate the distance to $2$-edge-connectivity we want to
estimate the numbers $c_{1}$ and $c_{2}$.

We use the algorithm in Figure 5, to find out if a node $s$ is in a $2$-edge-
connected component of size up to $\mathfrak{r}$, if this is also a connected
component, and its size. If $s$ is in a $2$-set of size up to $r+1$, that is
not a $1$-set, the first depth-first search will have a bridge connecting this
$2$-set to the rest of the graph in its search tree. Thus, after the second
search we know if $s$ is in a $2$-set of size up to $\mathfrak{r}$. If we are
in such a set, we do a third search to find out if this $2$-set is $2$-edge-
connected and thus a $2$-class-leaf. For every node it visits (except for the
last one), the first search will issue queries to find more nodes. At each
node, at most one query is made that finds a new node or tells us that there
are no more neighbours. At each node, we could also make queries that give us
edges to nodes that we already visited. The number of these queries is bounded
by the nodes already visited. Thus, in the first search there at most
$\sum_{i=0}^{r}(1+(i-1))=\frac{\mathfrak{r}(\mathfrak{r}+1)}{2}$ queries. The
second search is similar, but it will have to avoid certain edges, resulting
in some additional queries; the number of these edges is $\mathfrak{r}$, we
thus get a bound of $\frac{\mathfrak{r}(\mathfrak{r}+1)}{2}+\mathfrak{r}$
queries. The third search also needs to avoid edges to nodes not found by the
second search. But since the second search resulted in a $2$-set, there is
only one such edge. The number of queries in the second search is thus bounded
by
$\sum_{i=0}^{\mathfrak{r}-1}(1+(i-1))+(\mathfrak{r}-1)+1=\frac{\mathfrak{r}(\mathfrak{r}-1)}{2}+\mathfrak{r}$.
This gives us a bound of at most $\frac{3\mathfrak{r}(\mathfrak{r}+1)}{2}$
queries made by the algorithm in Figure 5.

⬇

float 1_2_leaf($s$, $\mathfrak{r}$)

{

Do first depth-first search from $s$ for up to $\mathfrak{r}+1$ nodes.

Let $n_{1}$ be the number of nodes found.

Do second depth-first search from $s$ for up to $\mathfrak{r}+1$ nodes,

never traversing an edge of the search tree of 1 in the same direction.

Let $n_{2}$ be the number of nodes found.

if ($n_{1}$ == $n_{2}$ == $\mathfrak{r}+1$) // Not in a 2-set of size at most
$\mathfrak{r}$

return $(0,0)$;

Do third depth-first search from $s$ for up to $\mathfrak{r}$ nodes,

only considering nodes found in 2,

never traversing an edge of the search tree of 2 in the same direction.

Let $n_{3}$ be the number of nodes found.

if ($n_{2}$ == $n_{3}$ == $n_{1}$) // In a 2-class-leaf of size $n_{3}$ that
is a 1-class-leaf

return $(n_{3},0)$;

else if ($n_{2}$ == $n_{3}$) // In a 2-class-leaf of size $n_{3}$ that is not
a 1-class-leaf

return $(0,n_{3})$;

else // In a $2$-set of size $n_{2}$, that is not a $2$-class-leaf

return $(0,0)$;

}

Figure 5: 1_2_leaf in pseudocode (the implementation resides in
zshg_component2() in zshg2_c.c)

⬇

float 1_2_leaves($\delta$, $p$, $n$)

{

$r=\frac{7}{2(1-p)\delta^{2}}$;

$a_{1}=a_{2}=0$;

$\ell=r(3\log_{e}(\frac{2}{\delta})+6)$;

for($i$ = 0; $i<r$; $i$++)

{

$s$ = rand_index($n$);

$x$ = rand_range($\frac{2}{\delta}$);

$(b_{1},b_{2})$ = 1_2_leaf($s$, $x$, &$\ell$);

if($b_{1}$)

$a_{1}$ += $b_{1}$;

if($b_{2}$)

$a_{2}$ += $b_{2}$;

}

return$\left(\frac{a_{1}n}{r},\frac{a_{2}n}{r}\right)$;

}

Figure 6: Counting algorithm for $2$-edge-connectivity in C-like pseudocode
(the implementation resides in zshg2_component() in zshg2_c.c)

######  Lemma 7.

For $\delta\leq 1$, with probability at least $p$, the return value of the
algorithm in Figure 6 is the minimum number of edge modifications necessary to
make the input graph $2$-edge-connected up to an error of at most
$\frac{1}{2}+\delta n$.

###### Proof.

Let $c_{1}^{*}$ be the number of 2-class-leaves that are 1-class-leaves in $G$
and of size at most $\frac{2}{\delta}$. Let $c_{2}^{*}$ be the number of
2-class-leaves that are not 1-class-leaves in $G$ and of size at most
$\frac{2}{\delta}$. The algorithm in Figure 6 gives us an estimate
$(\hat{C}_{1},\hat{C}_{2})$ for ($c_{1}^{*},c_{2}^{*})$ that we will later use
to estimate the distance to $2$-edge-connectivity. Let
$(\hat{c}_{1},\hat{c}_{2})$ be the return value of the algorithm in Figure 6
when ignoring the query bound $\ell$. Let $B_{1,i}$ be the value of $b_{1}$ in
iteration $i$, let $B_{2,i}$ be the value of $b_{2}$ in iteration $i$. Let
$C_{2,s}$ be the size of the $2$-class-leaf that contains $s$. Let
$\hat{c}\mathrel{\mathop{:}}=\frac{\hat{c}_{2}}{2}+\frac{1}{2}+\hat{c}_{1}+\frac{\delta
n}{2}$ be our estimate of the distance. Similarly to the calculation for the
distance to connectivity estimate, we get:

$\displaystyle E(B_{1,i})$ $\displaystyle=\frac{c_{1}^{*}}{n}\leq 1.$
$\displaystyle E(B_{2,i})$ $\displaystyle=\frac{c_{2}^{*}}{n}\leq 1.$
$\displaystyle Var(B_{1,i})$ $\displaystyle\leq 1.$ $\displaystyle
Var(B_{2,i})$ $\displaystyle\leq 1.$ $\displaystyle Cov(B_{1,i},B_{2,i})$
$\displaystyle\leq 1.$ $\displaystyle\hat{c_{1}}$
$\displaystyle=\frac{n}{r}\sum_{i=0}^{r-1}B_{1,i}.$ $\displaystyle\hat{c_{2}}$
$\displaystyle=\frac{n}{r}\sum_{i=0}^{r-1}B_{2,i}.$ $\displaystyle
E(\hat{c}_{1})$ $\displaystyle=r\frac{nE(B_{1,i})}{r}=c_{1}^{*}.$
$\displaystyle E(\hat{c}_{2})$
$\displaystyle=r\frac{nE(B_{2,i})}{r}=c_{2}^{*}.$ $\displaystyle
Var(\hat{c_{1}})$
$\displaystyle=\frac{n}{r^{2}}Var(B_{1,i})\leq\frac{n^{2}}{r}.$ $\displaystyle
Var(\hat{c}_{2})$
$\displaystyle=\frac{n}{r^{2}}Var(B_{2,i})\leq\frac{n^{2}}{r}.$ $\displaystyle
Cov(\hat{c}_{1},\hat{c}_{2})$
$\displaystyle\leq\sqrt{Var(\hat{c_{1}})}\sqrt{Var(\hat{c_{2}})}\leq\frac{n^{2}}{r}.$
$\displaystyle E(\hat{c})$
$\displaystyle=\frac{c_{2}^{*}}{2}+\frac{1}{2}+c_{1}^{*}+\frac{\delta n}{4}.$
$\displaystyle Var(\hat{c})$
$\displaystyle=\frac{1}{4}Var(\hat{c}_{2})+Var(\hat{c}_{1})+\frac{1}{2}Cov(\hat{c}_{1},\hat{c}_{2})\leq\frac{7}{4}\frac{n^{2}}{r}.$

The number of 2-class-leaves bigger than $\frac{2}{\delta}$ is at most
$\frac{\delta n}{2}$. Ignoring the query bound for now, we bound the
probability of the estimate being far off [19]:

$\displaystyle\left\lceil\frac{c_{2}}{2}\right\rceil+c_{1}-\frac{1}{2}-\frac{\delta
n}{4}\leq
E(\hat{c})\leq\left\lceil\frac{c_{2}}{2}\right\rceil+c_{1}+\frac{1}{2}+\frac{\delta
n}{4}.$ $\displaystyle P\left(\left|\hat{c}-E(\hat{c})\right|>\frac{3\delta
n}{4}\right)<\frac{Var(\hat{c})}{\left(\frac{3\delta
n}{4}\right)^{2}}\leq\frac{28}{9r\delta}=\frac{8}{9}(1-p).$

Let $Q_{i}$ be the number of queries made by a call to 1_2_leaves in the
algorithm in Figure 6.

$\displaystyle E(Q_{i})$
$\displaystyle\leq\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{1}{j}-\frac{1}{(j+1)^{2}}\right)\frac{3j(j+1)}{2}+\frac{1}{\left(\frac{2}{\delta}\right)^{2}}\frac{3\frac{2}{\delta}(\frac{2}{\delta}+1)}{2}=$
$\displaystyle=\sum_{j=1}^{\frac{2}{\delta}-1}\frac{6j+3}{3j(j+1)}+\frac{6+3\delta}{16}<3\log_{e}\left(\frac{2}{\delta}\right).$
$\displaystyle E(Q_{i}^{2})$
$\displaystyle\leq\sum_{j=1}^{\frac{2}{\delta}-1}\left(\frac{1}{j^{2}}-\frac{1}{(j+1)^{2}}\right)\left(\frac{3j(j+1)}{2}\right)^{2}+\frac{1}{\left(\frac{2}{\delta}\right)^{2}}\left(\frac{3\frac{2}{\delta}(\frac{2}{\delta}+1)}{2}\right)^{2}=$
$\displaystyle=\sum_{j=1}^{\delta{2}{\delta}-1}\left(\frac{9(j+1)^{2}}{4}-\frac{9j^{2}}{4}\right)+\frac{9(\frac{2}{\delta}+1)^{2}}{4}=\frac{9}{\delta^{2}}+\frac{9}{2\delta}.$
$\displaystyle Var(Q_{i})$
$\displaystyle=E(Q_{i}^{2})-E(Q_{i})^{2}\leq\frac{9}{\delta^{2}}+\frac{9}{2\delta}.$

Which we use to bound the probability of hitting the query bound [19].

$\displaystyle P\left(\sum_{i=0}^{r-1}Q_{i}>\ell\right)\leq
P\left(\sum_{i=0}^{r-1}Q_{i}>r(E(Q_{0})+6)\right)<\frac{Var(Q_{0})}{36r}<\frac{1}{9}(1-p).$
$\displaystyle P\left(\left|\hat{C}-E(\hat{c})\right|>\frac{3\delta
n}{4}\right)\leq P\left(\left|\hat{c}-E(\hat{c})\right|>\frac{3\delta
n}{4}\right)+P\left(\sum_{i=0}^{r-1}Q_{i}>\ell\right)<1-p.$

∎

######  Lemma 8.

The algorithm in Figure 6 has query complexity
$O(\frac{1}{(1-p)\delta^{2}}\log_{e}(\frac{1}{\delta}))$.

###### Proof.

This follows directly from the choice of $\ell$. ∎

For general $k$, we use a generalization of Lemma 6 above.

The _edge-demand_ of a set of nodes $U\subsetneq V$ is

$\Phi_{k}(U)\mathrel{\mathop{:}}=\max\left\\{0,k-d(U),\sum_{W\sqsubset
U}\Phi(W)\right\\},\\\ \Phi_{k}(V)\mathrel{\mathop{:}}=\sum_{W\sqsubset
V}\Phi(W).\\\ $

######  Lemma 9 (Naor et alii [13] for connected graphs, generalized by Marko
[12]).

Let $G=(V,E)$ be a graph and $k\in\mathbb{N}$. The distance of $G$ to
$k$-edge-connectivity is
$\frac{1}{|E|}\left\lceil\frac{\Phi_{k}(V)}{2}\right\rceil$.

To estimate $\Phi_{k}(V)$ we can estimate the number of extreme sets $U$, for
which $\Phi_{k}(U)\neq\sum_{W\sqsubset U}\Phi(W)$, and their $\Phi_{k}(U)$ and
then get the estimate of $\Phi_{k}(V)$ as the weighted sum.

For $k=3$, there are $4$ kinds of such $U$: 3-class-leaves that are connected
components (for these, $\Phi_{k}=3$), 3-class-leaves that are connected to the
rest of the graph by a single bridge (for these, $\Phi_{k}=2$), 3-class-leaves
that are connected to the rest of the graph by a exactly two edges and are not
inside the fourth kind (for these, $\Phi_{k}=1$), connected components that
are 2-edge-connected and contain exactly 2 3-edge-leaves (for these,
$\Phi_{k}=3$). Sets of the fourth kind contribute $3$ each to $\Phi_{3}$, but
for the algorithm it is simpler to make them contribute $1$ and instead ignore
the condition on sets of the third kind not being in the fourth kind. Let
$C_{0}$ be the number of of 3-class-leaves that are connected components, let
$C_{1}$ be the number of of 3-class-leaves that are connected to the rest of
the graph by a single bridge, let $C_{2}$ be the number of $3$-class-leaves
that are connected to the rest of the graph by a exactly two edges, let
$C_{3}$ be the number of connected components that are 2-edge-connected and
contain exactly 2 3-edge-leaves.

$\Phi_{3}(V)=3C_{0}+2C_{1}+1C_{2}+1C_{3}-k.$

To find out if a node is in any such set of size at most $\mathfrak{r}$, we
use the algorithm in Figure 7. To decide if a node is in a set of the fourth
kind, it uses the plant graph [16] of the connected component, which can be
computed in time quadratic in the number of nodes in the connected component
[20]. This results in total query complexity $O(\mathfrak{r}^{3})$

We can use this algorithm to construct an algorithm for estimating the
distance to 3-edge-connectivity similar to how we did so for 2-edge-
connectivity above.

⬇

float 1_2_3_leaf($s$, $\mathfrak{r}$)

{

Do first depth-first search from $s$ for up to $\mathfrak{r}+1$ nodes.

Let $n_{0}$ be the number of nodes found. Let $m_{0}$ be the set of nodes
found.

Do a second depth-first search from $s$ for up to $\mathfrak{r}+1$ nodes,

never traversing an edge of the search tree of 1 in the same direction.

Let $n_{1}$ be the number of nodes found.

if($n_{0}\leq\mathfrak{r}$ && $n_{0}$ == $n_{1}$) // $s$ is in a 2-connected
connected component of size at most $\mathfrak{r}$

{

Compute the plant graph $\Gamma$ of the subgraph induced by $m_{0}$.

if ($\Gamma$ has exactly 2 nodes of degree 1 and no cycles)

$n_{2}=n_{0}$;

else

$n_{2}$ = 0;

}

else

$n_{2}$ = 0;

For every edge traversed in the first depth-first search,

invoke 2\_set searching for up to $\mathfrak{r}+1$ nodes in graph with that
edge omitted.

Let $m$ be a smallest set found by 2\_set.

if ($|m|$ > $\mathfrak{r}$) // Not in a 3-set of size at most $\mathfrak{r}$.

return $(0,0,0,n_{2})$;

Do second depth-first search from $s$ for up to $\mathfrak{r}+1$ nodes,

never traversing an edge of the search tree of 1 in the same direction.

Let $n_{2}$ be the number of nodes found.

Do a first depth-first search from $s$ for up to $|m|$ nodes in the subgraph
induced by $m$.

Let $n_{1}$ be the number of nodes found.

For every edge traversed in the previous depth-first search,

invoke 2\_set searching for up to $|m|$ nodes in the subgraph induced by $m$
with that edge omitted.

Let $m^{\prime}$ be a smallest set found by 2\_set.

if ($|m|$ != $|m^{\prime}|$) // In a 3-set of size $|m|$ that is not a
3-class-leaf.

return $(0,0,0,n_{2})$;

if($d(m)$ == 0)

return $(|m|,0,0,n_{2})$;

else if ($d(m)$ == 1)

return $(0,|m|,0,n_{2})$;

else // $d(m)$ == 2

return $(0,0,|m|,n_{2})$;

}

Figure 7: 1_2_3_leaf in pseudocode

From our estimate of the number of components, we also obtain an estimate for
the distance to eulerianity (Figure 8). The algorithm 1_sets_even is a variant
of 1_sets that only counts components in which all nodes have even degree.

⬇

float euler_distance($\epsilon$, $p$, $n$)

{

$\delta=\frac{\epsilon}{d}$;

$c=$ 1_sets_even$\left(\frac{\delta}{2},\frac{1-p}{2}\right)$;

$m=\frac{2}{\delta^{2}(1+p)}$;

$a=0$;

for($i$ = 0; $i<m$; $i$++)

if($d(s)$ odd)

$a$++;

$u=\frac{an}{m}$;

$e=c+\frac{u}{2}$;

return$\left(\frac{e}{nd}\right)$;

}

Figure 8: Algorithm for estimating the distance to eulerianity in C-like
pseudocode (the implementation resides in zshg_euler_distance() in zshg_c.c)

######  Lemma 10.

With probability at least $p$, value returned by the algorithm in Figure 8 is
the distance of the graph to eulerianity with an error of at most $\epsilon$.

###### Proof.

A graph is eulerian if it is connected and all nodes have even degree [7, 11,
9].

Let $C$ be the number of connected components, in which all nodes have even
degree. Let $U$ be the number of nodes of odd degree. Then the minimum number
of edge modifications necessary to make the graph connected is
$C+\frac{U}{2}$: In every component, the number of nodes of odd degree is even
[10], so one can add a cycle that goes through all components (using
preexisting edges within components that contain nodes of odd degree),
requiring only that number of edges. It is not possible to do with fewer edge
modifications: At every component that contains only nodes of even degree, we
have to add at least 2 incident edges to a node, and at every node of odd
degree in the graph, we have to add or remove an edge. Let q be the
probability of the algorithm returning a value outside the error bound.

$\displaystyle
E(c)=C,E(a_{i})=\frac{U}{n},E(a_{i}^{2})=\frac{U}{n},Var(u)\leq\frac{U}{n}.$
$\displaystyle
E(u)=U,Var(u)=\frac{n^{2}}{m}Var(a_{i})\leq\frac{n}{m}U\leq\frac{n^{2}}{m}.$
$\displaystyle P(|u-U|>\epsilon
dn)<\frac{Var(u)}{\epsilon^{2}d^{2}n^{2}}\leq\frac{1}{m\epsilon^{2}d^{2}}=\frac{1-p}{2}.$
$\displaystyle q\leq P(|c-C|>\frac{\epsilon dn}{2})+P(|u-U|>\epsilon dn)\leq
1-\frac{1-p}{2}+1-\frac{1+p}{2}=1-p.$

∎

######  Lemma 11.

The algorithm in Figure 8 has query complexity
$O\left(\frac{1}{(1-p)\delta^{2}}\log_{e}(\frac{1}{\delta})\right)$.

###### Proof.

For 1_sets_even we get that complexity like in Lemma 5 and the number of
further queries is $m\in O\left(\frac{1}{\delta^{2}(1-p)}\right)\subseteq
O\left(\frac{1}{(1-p)\delta^{2}}\log_{e}(\frac{1}{\delta})\right)$. ∎

######  Theorem 2.

There are estimates for the distance to connectivity, $2$-edge-connectivity
and eulerianity with complexity
$O(\frac{1}{(1-p)\delta^{2}}\log_{e}(\frac{1}{\delta}))$. There is an estimate
for the distance to $3$-edge-connectivity with complexity
$O(\frac{1}{(1-p)\delta^{3}}\log_{e}(\frac{1}{\delta}))$.

###### Proof.

For the distance to connectivity this follows from Lemmata 4, 5. For the
distance to $2$-edge-connectivity, this follows from Lemmata 6, 7, 8. For the
distance to eulerianity, this follows from Lemmata 10, 11. For the distance to
$3$-edge-connectivity, this can be proven using Lemma 9 and the algorithm in
Figure 7. ∎

##  6 Tolerant Connectivity Test

⬇

bool zshg_tolerant($\epsilon_{1}$, $\epsilon_{2}$, $p$, $d$, $n$)

{

$\delta=\frac{\epsilon_{2}-\epsilon_{1}}{4}d$;

$e$ = 1_sets($\delta$, $p$, $n$) - 1;

return$\left(e\leq\frac{(\epsilon_{1}+\epsilon_{2})dn}{4}\right)$;

}

Figure 9: Tolerant connectivity testing algorithm in C-like pseudocode (the
implementation resides in zshg_tolerant() in zshg_c.c)

######  Lemma 12.

For any graph that is $\epsilon_{1}$-close to connectivity, the algorithm in
Figure 9 returns true with probability at least $p$. For any graph that is
$\epsilon_{2}$-far from connectivity, the algorithm in Figure 9 returns false
with probability at least $p$.

###### Proof.

Let $m=\frac{dn}{2}$ be the number of edges in the graph.

Case 1: The graph is $\epsilon_{1}$-close to connectivity, i. e. it consists
of at most $\epsilon_{1}m+1$ connected components.

$\displaystyle
P\left(e>\frac{(\epsilon_{1}+\epsilon_{2})dn}{4}\right)=P\left(e>\frac{\epsilon_{1}+\epsilon_{2}}{2}m\right)\leq$
$\displaystyle
P\left(\textrm{1\\_sets}(\delta,p,n)>\frac{\epsilon_{1}+\epsilon_{2}}{2}m+1\right)=P\left(\textrm{1\\_sets}(\delta,p,n)>(\epsilon_{1}m+1)+\delta
n\right)\leq$ $\displaystyle(1-p).$

Case 2: The graph is $\epsilon_{2}$-far from connectivity, i. e. it consists
of at least $\epsilon_{2}m+2$ connected components.

$\displaystyle
P\left(e\leq\frac{(\epsilon_{1}+\epsilon_{2})d}{4}\right)=P\left(e\leq\frac{\epsilon_{1}+\epsilon_{2}}{2}m\right)\leq$
$\displaystyle
P\left(\textrm{1\\_sets}(\delta,p,n)\leq\frac{\epsilon_{1}+\epsilon_{2}}{2}m+1\right)=P\left(\textrm{1\\_sets}(\delta,p,n)\leq(\epsilon_{2}m+1)-\delta
n\right)\leq$ $\displaystyle(1-p).$

Case 3: The graph is neither $\epsilon_{1}$-close to nor $\epsilon_{2}$-far
from connectivity. We don’t care. ∎

In a similar way, we can obtain tolerant testers for $2$-edge-connectivity,
$3$-edge-connectivity and eulerianity.

######  Theorem 3.

There are tolerant testers for connectivity, $2$-edge-connectivity and
eulerianity with complexity
$O\left(\frac{1}{(1-p)(\epsilon_{2}-\epsilon_{1})^{2}d^{2}}\log_{e}(\frac{1}{(\epsilon_{2}-\epsilon_{1})d})\right)$.
There is a tolerant tester for $3$-edge-connectivity with complexity
$O\left(\frac{1}{(1-p)(\epsilon_{2}-\epsilon_{1})^{3}d^{3}}\log_{e}(\frac{1}{(\epsilon_{2}-\epsilon_{1})d})\right)$.

###### Proof.

For connectivity this follows from Lemmata 4,5,12. The proofs for the others
are similar. ∎

##  7 Remarks

While the results on correctness also hold for multigraphs, the results on
query complexity don’t.

All of our algorithms can be parallelized easily and we did so in the
implementations. The parallel versions offer an advantage when multiple
pending queries can be answered more efficiently, e.g. in the case of large
(too big to fit into RAM) graphs stored on an SSD (SSDs typically achieve
maximum throughput for random reads at about 16 simultaneous pending reads) or
in the case of the queries being processed by a remote server on a network.

## References

  * [1] Petra Berenbrink, Bruce Krayenhoff, and Frederik Mallmann-Trenn.  Estimating the number of connected components in sublinear time.  Information Processing Letters, 114(11):639 – 642, 2014.
  * [2] Jacobo Bernoulli.  Positiones arithmeticæ de seriebus infinitis Earumque summa finita. Dissertation, 1689.
  * [3] Henry R. Brahana.  A proof of petersen’s theorem.  Annals of Mathematics, 19(1):59–63, 1917.
  * [4] Andrea Campagna, Alan Guo, and Ronitt Rubinfeld.  Local Reconstructors and Tolerant Testers for Connectivity and Diameter, pages 411–424.  Springer, Berlin, Heidelberg, 2013.
  * [5] Bernard Chazelle, Ronitt Rubinfeld, and Luca Trevisan.  Approximating the minimum spanning tree weight in sublinear time.  SIAM Journal on Computing, 34(6):1370–1379, 2005.
  * [6] Kapali P. Eswaran and Robert E. Tarjan.  Augmentation problems.  SIAM Journal on Computing, 5:653–665, 1976.
  * [7] Leonhard Eulero.  Solutio problematis ad geometriam situs pertinentis.  Commentarii academiae scientiarum Petropolitanae, 8:128–140, 1741\.
  * [8] Oded Goldreich and Dana Ron.  Property Testing in Bounded Degree Graphs.  Algorithmica, 32(2):302–343, Feb 2002.
  * [9] Carl Hierholzer.  Ueber die Möglichkeit, einen Linienzug ohne Wiederholung und ohne Unterbrechung zu umfahren.  Mathematische Annalen, 6:30–32, 1873.
  * [10] Dénes König.  Theorie der endlichen und unendlichen Graphen.  Chelsea, 1935.
  * [11] Johann B. Listing.  Vorstudien zur Topologie.  Vandenheock und Ruprecht, 1847.
  * [12] Sharon Marko.  Distance approximation in bounded-degree and general sparse graphs.  Master’s thesis, 2005.
  * [13] Dalit Naor, Dan Gusfield, and Charles Martel.  A fast algorithm for optimally increasing the edge connectivity.  SIAM Journal on Computing, 26:1139–1165, 08 1997.
  * [14] Michal Parnas and Dana Ron.  Testing the diameter of graphs.  Random Structures & Algorithms, 20(2):165–183, 2002.
  * [15] Dana Ron.  Algorithmic and Analysis Techniques in Property Testing.  Foundations and Trends in Theoretical Computer Science, 5(2):73–205, February 2010.
  * [16] Efim A. Диниц, Alexander В. Карзанов, and Michael В. Ломоносов.  О структуре системы минимальных реберных разрезов графа.  In А.А. Фридман, editor, Исследования по Дискретной Оптимизации, pages 290–306. Наука, Москва, 1976.
  * [17] Renatus F. Slusius.  Mesolabum.  Apud Guilielmum Henricum Streel, Serenissimæ Suæ Celsitudinis Typographum, 1668.
  * [18] Richard M. Stallman, Joshua Gay, and Lawrence Lessig.  Free Software, Free Society: Selected Essays of Richard M. Stallman.  GNU Press, 2002.
  * [19] Pafnuti L. Tchébychef.  Des valeurs moyennes.  Journal de mathématiques pures et appliquées $2^{e}$ série, 12:177–184, 1867.
  * [20] Alexander В. Карзанов and Евгений A. Тимофеев.  Эффективный алгоритм нахождения всех минимальных реберных разрезов неориентированного графа.  Кибернетика, (2):8–12, 1986\.